Actual source code: ex263.c

  1: static char help[] = "Tests MatForwardSolve and MatBackwardSolve for LU and Cholesky decompositions using C/Pardiso.\n\n";

  3: #include <petscdmda.h>
  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   DM            da;
  9:   DMDALocalInfo info;
 10:   Mat           A, F;
 11:   Vec           x, y, b, ytmp;
 12:   IS            rowis, colis;
 13:   PetscInt      i, j, k, n = 5;
 14:   PetscBool     CHOL = PETSC_FALSE;
 15:   MatStencil    row, cols[5];
 16:   PetscScalar   vals[5];
 17:   PetscReal     norm2, tol = 100. * PETSC_MACHINE_EPSILON;
 18:   PetscRandom   rdm;
 19:   MatFactorInfo finfo;

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 23:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 25:   PetscCall(DMDACreate2d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, n, n, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
 26:   PetscCall(DMSetUp(da));

 28:   PetscCall(DMCreateMatrix(da, &A));
 29:   PetscCall(PetscOptionsHasName(NULL, NULL, "-chol", &CHOL));
 30:   if (CHOL) PetscCall(MatSetType(A, MATSBAIJ));
 31:   else PetscCall(MatSetType(A, MATBAIJ));
 32:   PetscCall(MatSetUp(A));

 34:   PetscCall(DMDAGetLocalInfo(da, &info));
 35:   for (j = info.ys; j < info.ys + info.ym; j++) {
 36:     for (i = info.xs; i < info.xs + info.xm; i++) {
 37:       row.j = j;
 38:       row.i = i;

 40:       k = 0;
 41:       if (j != 0) {
 42:         cols[k].j = j - 1;
 43:         cols[k].i = i;
 44:         vals[k]   = -1;
 45:         ++k;
 46:       }
 47:       if (i != 0) {
 48:         cols[k].j = j;
 49:         cols[k].i = i - 1;
 50:         vals[k]   = -1;
 51:         ++k;
 52:       }
 53:       cols[k].j = j;
 54:       cols[k].i = i;
 55:       vals[k]   = 4;
 56:       ++k;
 57:       if (j != info.my - 1) {
 58:         cols[k].j = j + 1;
 59:         cols[k].i = i;
 60:         vals[k]   = -1;
 61:         ++k;
 62:       }
 63:       if (i != info.mx - 1) {
 64:         cols[k].j = j;
 65:         cols[k].i = i + 1;
 66:         vals[k]   = -1;
 67:         ++k;
 68:       }
 69:       PetscCall(MatSetValuesStencil(A, 1, &row, k, cols, vals, INSERT_VALUES));
 70:     }
 71:   }
 72:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 73:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 75:   if (CHOL) PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));

 77:   /* Create vectors for error checking */
 78:   PetscCall(MatCreateVecs(A, &x, &b));
 79:   PetscCall(VecDuplicate(x, &y));
 80:   PetscCall(VecDuplicate(x, &ytmp));
 81:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
 82:   PetscCall(PetscRandomSetFromOptions(rdm));
 83:   PetscCall(VecSetRandom(x, rdm));
 84:   PetscCall(MatMult(A, x, b));

 86:   PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &rowis, &colis));
 87:   if (CHOL) {
 88:     PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test Cholesky...\n"));
 89:     PetscCall(MatGetFactor(A, MATSOLVERMKL_CPARDISO, MAT_FACTOR_CHOLESKY, &F));
 90:     PetscCall(MatCholeskyFactorSymbolic(F, A, rowis, &finfo));
 91:     PetscCall(MatCholeskyFactorNumeric(F, A, &finfo));
 92:   } else {
 93:     PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test LU...\n"));
 94:     PetscCall(MatGetFactor(A, MATSOLVERMKL_CPARDISO, MAT_FACTOR_LU, &F));
 95:     PetscCall(MatLUFactorSymbolic(F, A, rowis, colis, &finfo));
 96:     PetscCall(MatLUFactorNumeric(F, A, &finfo));
 97:   }

 99:   PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test MatForwardSolve...\n"));
100:   PetscCall(MatForwardSolve(F, b, ytmp));
101:   PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test MatBackwardSolve...\n"));
102:   PetscCall(MatBackwardSolve(F, ytmp, y));
103:   PetscCall(VecAXPY(y, -1.0, x));
104:   PetscCall(VecNorm(y, NORM_2, &norm2));
105:   if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));

107:   /* Free data structures */
108:   PetscCall(ISDestroy(&rowis));
109:   PetscCall(ISDestroy(&colis));
110:   PetscCall(MatDestroy(&F));
111:   PetscCall(MatDestroy(&A));
112:   PetscCall(DMDestroy(&da));
113:   PetscCall(PetscRandomDestroy(&rdm));
114:   PetscCall(VecDestroy(&x));
115:   PetscCall(VecDestroy(&y));
116:   PetscCall(VecDestroy(&ytmp));
117:   PetscCall(VecDestroy(&b));
118:   PetscCall(PetscFinalize());
119:   return 0;
120: }

122: /*TEST

124:    test:
125:       requires: mkl_cpardiso
126:       nsize: 4

128:    test:
129:       suffix: 2
130:       requires: !complex mkl_cpardiso
131:       nsize: 4
132:       args: -chol

134: TEST*/