Actual source code: ex30.c

  1: static char help[] = "Tests ILU and ICC factorization with and without matrix ordering on seqaij format, and illustrates drawing of matrix sparsity structure with MatView().\n\
  2:   Input parameters are:\n\
  3:   -lf <level> : level of fill for ILU (default is 0)\n\
  4:   -lu : use full LU or Cholesky factorization\n\
  5:   -m <value>,-n <value> : grid dimensions\n\
  6: Note that most users should employ the KSP interface to the\n\
  7: linear solvers instead of using the factorization routines\n\
  8: directly.\n\n";

 10: #include <petscmat.h>

 12: int main(int argc, char **args)
 13: {
 14:   Mat           C, A;
 15:   PetscInt      i, j, m = 5, n = 5, Ii, J, lf = 0;
 16:   PetscBool     LU = PETSC_FALSE, CHOLESKY = PETSC_FALSE, TRIANGULAR = PETSC_FALSE, MATDSPL = PETSC_FALSE, flg, matordering, use_mkl_pardiso = PETSC_FALSE;
 17:   PetscScalar   v;
 18:   IS            row, col;
 19:   PetscViewer   viewer1, viewer2;
 20:   MatFactorInfo info;
 21:   Vec           x, y, b, ytmp;
 22:   PetscReal     norm2, norm2_inplace, tol = 100. * PETSC_MACHINE_EPSILON;
 23:   PetscRandom   rdm;
 24:   PetscMPIInt   size;
 25:   char          pack[PETSC_MAX_PATH_LEN];

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 29:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 30:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 31:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 33:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL));

 35:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, 0, 0, 0, 400, 400, &viewer1));
 36:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, 0, 400, 0, 400, 400, &viewer2));

 38:   PetscCall(MatCreate(PETSC_COMM_SELF, &C));
 39:   PetscCall(MatSetSizes(C, m * n, m * n, m * n, m * n));
 40:   PetscCall(MatSetFromOptions(C));
 41:   PetscCall(MatSetUp(C));

 43:   /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */
 44:   for (i = 0; i < m; i++) {
 45:     for (j = 0; j < n; j++) {
 46:       v  = -1.0;
 47:       Ii = j + n * i;
 48:       J  = Ii - n;
 49:       if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 50:       J = Ii + n;
 51:       if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 52:       J = Ii - 1;
 53:       if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 54:       J = Ii + 1;
 55:       if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 56:       v = 4.0;
 57:       PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
 58:     }
 59:   }
 60:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 61:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 63:   PetscCall(MatIsSymmetric(C, 0.0, &flg));
 64:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "C is non-symmetric");

 66:   PetscCall(MatSetOption(C, MAT_SPD, PETSC_TRUE));

 68:   /* Create vectors for error checking */
 69:   PetscCall(MatCreateVecs(C, &x, &b));
 70:   PetscCall(VecDuplicate(x, &y));
 71:   PetscCall(VecDuplicate(x, &ytmp));
 72:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
 73:   PetscCall(PetscRandomSetFromOptions(rdm));
 74:   PetscCall(VecSetRandom(x, rdm));
 75:   PetscCall(MatMult(C, x, b));

 77:   PetscCall(PetscOptionsHasName(NULL, NULL, "-mat_ordering", &matordering));
 78:   if (matordering) {
 79:     PetscCall(MatGetOrdering(C, MATORDERINGRCM, &row, &col));
 80:   } else {
 81:     PetscCall(MatGetOrdering(C, MATORDERINGNATURAL, &row, &col));
 82:   }

 84:   PetscCall(PetscOptionsHasName(NULL, NULL, "-display_matrices", &MATDSPL));
 85:   if (MATDSPL) {
 86:     printf("original matrix:\n");
 87:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INFO));
 88:     PetscCall(MatView(C, PETSC_VIEWER_STDOUT_SELF));
 89:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF));
 90:     PetscCall(MatView(C, PETSC_VIEWER_STDOUT_SELF));
 91:     PetscCall(MatView(C, viewer1));
 92:   }

 94:   /* Compute LU or ILU factor A */
 95:   PetscCall(MatFactorInfoInitialize(&info));

 97:   info.fill          = 1.0;
 98:   info.diagonal_fill = 0;
 99:   info.zeropivot     = 0.0;

101:   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
102: #if defined(PETSC_HAVE_MKL_PARDISO)
103:   PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &use_mkl_pardiso));
104: #endif

106:   PetscCall(PetscOptionsHasName(NULL, NULL, "-lu", &LU));
107:   if (LU) {
108:     printf("Test LU...\n");
109:     if (use_mkl_pardiso) PetscCall(MatGetFactor(C, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &A));
110:     else PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A));
111:     PetscCall(MatLUFactorSymbolic(A, C, row, col, &info));
112:   } else {
113:     printf("Test ILU...\n");
114:     info.levels = lf;

116:     PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ILU, &A));
117:     PetscCall(MatILUFactorSymbolic(A, C, row, col, &info));
118:   }
119:   PetscCall(MatLUFactorNumeric(A, C, &info));

121:   /* test MatForwardSolve() and MatBackwardSolve() with MKL Pardiso*/
122:   if (LU && use_mkl_pardiso) {
123:     PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR));
124:     if (TRIANGULAR) {
125:       printf("Test MatForwardSolve...\n");
126:       PetscCall(MatForwardSolve(A, b, ytmp));
127:       printf("Test MatBackwardSolve...\n");
128:       PetscCall(MatBackwardSolve(A, ytmp, y));
129:       PetscCall(VecAXPY(y, -1.0, x));
130:       PetscCall(VecNorm(y, NORM_2, &norm2));
131:       if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));
132:     }
133:   }

135:   /* Solve A*y = b, then check the error */
136:   PetscCall(MatSolve(A, b, y));
137:   PetscCall(VecAXPY(y, -1.0, x));
138:   PetscCall(VecNorm(y, NORM_2, &norm2));
139:   PetscCall(MatDestroy(&A));

141:   /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
142:   if (!LU && lf == 0) {
143:     PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &A));
144:     PetscCall(MatILUFactor(A, row, col, &info));
145:     /*
146:     printf("In-place factored matrix:\n");
147:     PetscCall(MatView(C,PETSC_VIEWER_STDOUT_SELF));
148:     */
149:     PetscCall(MatSolve(A, b, y));
150:     PetscCall(VecAXPY(y, -1.0, x));
151:     PetscCall(VecNorm(y, NORM_2, &norm2_inplace));
152:     PetscCheck(PetscAbs(norm2 - norm2_inplace) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ILU(0) %g and in-place ILU(0) %g give different residuals", (double)norm2, (double)norm2_inplace);
153:     PetscCall(MatDestroy(&A));
154:   }

156:   /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
157:   PetscCall(PetscOptionsHasName(NULL, NULL, "-chol", &CHOLESKY));
158:   if (CHOLESKY) {
159:     printf("Test Cholesky...\n");
160:     lf = -1;
161:     if (use_mkl_pardiso) PetscCall(MatGetFactor(C, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &A));
162:     else PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &A));
163:     PetscCall(MatCholeskyFactorSymbolic(A, C, row, &info));
164:   } else {
165:     printf("Test ICC...\n");
166:     info.levels        = lf;
167:     info.fill          = 1.0;
168:     info.diagonal_fill = 0;
169:     info.zeropivot     = 0.0;

171:     PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ICC, &A));
172:     PetscCall(MatICCFactorSymbolic(A, C, row, &info));
173:   }
174:   PetscCall(MatCholeskyFactorNumeric(A, C, &info));

176:   /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
177:   if (CHOLESKY) {
178:     PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR));
179:     if (TRIANGULAR) {
180:       printf("Test MatForwardSolve...\n");
181:       PetscCall(MatForwardSolve(A, b, ytmp));
182:       printf("Test MatBackwardSolve...\n");
183:       PetscCall(MatBackwardSolve(A, ytmp, y));
184:       PetscCall(VecAXPY(y, -1.0, x));
185:       PetscCall(VecNorm(y, NORM_2, &norm2));
186:       if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));
187:     }
188:   }

190:   PetscCall(MatSolve(A, b, y));
191:   PetscCall(MatDestroy(&A));
192:   PetscCall(VecAXPY(y, -1.0, x));
193:   PetscCall(VecNorm(y, NORM_2, &norm2));
194:   if (lf == -1 && norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ:   Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2));

196:   /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
197:   if (!CHOLESKY && lf == 0 && !matordering) {
198:     PetscCall(MatConvert(C, MATSBAIJ, MAT_INITIAL_MATRIX, &A));
199:     PetscCall(MatICCFactor(A, row, &info));
200:     /*
201:     printf("In-place factored matrix:\n");
202:     PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF));
203:     */
204:     PetscCall(MatSolve(A, b, y));
205:     PetscCall(VecAXPY(y, -1.0, x));
206:     PetscCall(VecNorm(y, NORM_2, &norm2_inplace));
207:     PetscCheck(PetscAbs(norm2 - norm2_inplace) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ICC(0) %g and in-place ICC(0) %g give different residuals", (double)norm2, (double)norm2_inplace);
208:     PetscCall(MatDestroy(&A));
209:   }

211:   /* Free data structures */
212:   PetscCall(ISDestroy(&row));
213:   PetscCall(ISDestroy(&col));
214:   PetscCall(MatDestroy(&C));
215:   PetscCall(PetscViewerDestroy(&viewer1));
216:   PetscCall(PetscViewerDestroy(&viewer2));
217:   PetscCall(PetscRandomDestroy(&rdm));
218:   PetscCall(VecDestroy(&x));
219:   PetscCall(VecDestroy(&y));
220:   PetscCall(VecDestroy(&ytmp));
221:   PetscCall(VecDestroy(&b));
222:   PetscCall(PetscFinalize());
223:   return 0;
224: }

226: /*TEST

228:    test:
229:       args: -mat_ordering -display_matrices -nox
230:       filter: grep -v " MPI process"

232:    test:
233:       suffix: 2
234:       args: -mat_ordering -display_matrices -nox -lu -chol

236:    test:
237:       suffix: 3
238:       args: -mat_ordering -lu -chol -triangular_solve

240:    test:
241:       suffix: 4

243:    test:
244:       suffix: 5
245:       args: -lu -chol

247:    test:
248:       suffix: 6
249:       args: -lu -chol -triangular_solve
250:       output_file: output/ex30_3.out

252:    test:
253:       suffix: 7
254:       requires: mkl_pardiso
255:       args: -lu -chol -mat_solver_type mkl_pardiso
256:       output_file: output/ex30_5.out

258:    test:
259:       suffix: 8
260:       requires: mkl_pardiso
261:       args: -lu -mat_solver_type mkl_pardiso -triangular_solve

263: TEST*/