Actual source code: ex30.c


  2: 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\
  3:   Input parameters are:\n\
  4:   -lf <level> : level of fill for ILU (default is 0)\n\
  5:   -lu : use full LU or Cholesky factorization\n\
  6:   -m <value>,-n <value> : grid dimensions\n\
  7: Note that most users should employ the KSP interface to the\n\
  8: linear solvers instead of using the factorization routines\n\
  9: directly.\n\n";

 11: #include <petscmat.h>

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

 28:   PetscInitialize(&argc, &args, (char *)0, help);
 29:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 31:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 32:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 33:   PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL);

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

 38:   MatCreate(PETSC_COMM_SELF, &C);
 39:   MatSetSizes(C, m * n, m * n, m * n, m * n);
 40:   MatSetFromOptions(C);
 41:   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) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 50:       J = Ii + n;
 51:       if (J < m * n) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 52:       J = Ii - 1;
 53:       if (J >= 0) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 54:       J = Ii + 1;
 55:       if (J < m * n) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 56:       v = 4.0;
 57:       MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
 58:     }
 59:   }
 60:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
 61:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

 63:   MatIsSymmetric(C, 0.0, &flg);

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

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

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

 92:   /* Compute LU or ILU factor A */
 93:   MatFactorInfoInitialize(&info);

 95:   info.fill          = 1.0;
 96:   info.diagonal_fill = 0;
 97:   info.zeropivot     = 0.0;

 99:   PetscOptionsHasName(NULL, NULL, "-lu", &LU);
100:   if (LU) {
101:     printf("Test LU...\n");
102:     MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A);
103:     MatLUFactorSymbolic(A, C, row, col, &info);
104:   } else {
105:     printf("Test ILU...\n");
106:     info.levels = lf;

108:     MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ILU, &A);
109:     MatILUFactorSymbolic(A, C, row, col, &info);
110:   }
111:   MatLUFactorNumeric(A, C, &info);

113:   /* Solve A*y = b, then check the error */
114:   MatSolve(A, b, y);
115:   VecAXPY(y, -1.0, x);
116:   VecNorm(y, NORM_2, &norm2);
117:   MatDestroy(&A);

119:   /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
120:   if (!LU && lf == 0) {
121:     MatDuplicate(C, MAT_COPY_VALUES, &A);
122:     MatILUFactor(A, row, col, &info);
123:     /*
124:     printf("In-place factored matrix:\n");
125:     MatView(C,PETSC_VIEWER_STDOUT_SELF);
126:     */
127:     MatSolve(A, b, y);
128:     VecAXPY(y, -1.0, x);
129:     VecNorm(y, NORM_2, &norm2_inplace);
131:     MatDestroy(&A);
132:   }

134:   /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
135:   CHOLESKY = LU;
136:   if (CHOLESKY) {
137:     printf("Test Cholesky...\n");
138:     lf = -1;
139:     MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &A);
140:     MatCholeskyFactorSymbolic(A, C, row, &info);
141:   } else {
142:     printf("Test ICC...\n");
143:     info.levels        = lf;
144:     info.fill          = 1.0;
145:     info.diagonal_fill = 0;
146:     info.zeropivot     = 0.0;

148:     MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ICC, &A);
149:     MatICCFactorSymbolic(A, C, row, &info);
150:   }
151:   MatCholeskyFactorNumeric(A, C, &info);

153:   /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
154:   if (lf == -1) {
155:     PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR);
156:     if (TRIANGULAR) {
157:       printf("Test MatForwardSolve...\n");
158:       MatForwardSolve(A, b, ytmp);
159:       printf("Test MatBackwardSolve...\n");
160:       MatBackwardSolve(A, ytmp, y);
161:       VecAXPY(y, -1.0, x);
162:       VecNorm(y, NORM_2, &norm2);
163:       if (norm2 > tol) PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2);
164:     }
165:   }

167:   MatSolve(A, b, y);
168:   MatDestroy(&A);
169:   VecAXPY(y, -1.0, x);
170:   VecNorm(y, NORM_2, &norm2);
171:   if (lf == -1 && norm2 > tol) PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ:   Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2);

173:   /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
174:   if (!CHOLESKY && lf == 0 && !matordering) {
175:     MatConvert(C, MATSBAIJ, MAT_INITIAL_MATRIX, &A);
176:     MatICCFactor(A, row, &info);
177:     /*
178:     printf("In-place factored matrix:\n");
179:     MatView(A,PETSC_VIEWER_STDOUT_SELF);
180:     */
181:     MatSolve(A, b, y);
182:     VecAXPY(y, -1.0, x);
183:     VecNorm(y, NORM_2, &norm2_inplace);
185:     MatDestroy(&A);
186:   }

188:   /* Free data structures */
189:   ISDestroy(&row);
190:   ISDestroy(&col);
191:   MatDestroy(&C);
192:   PetscViewerDestroy(&viewer1);
193:   PetscViewerDestroy(&viewer2);
194:   PetscRandomDestroy(&rdm);
195:   VecDestroy(&x);
196:   VecDestroy(&y);
197:   VecDestroy(&ytmp);
198:   VecDestroy(&b);
199:   PetscFinalize();
200:   return 0;
201: }

203: /*TEST

205:    test:
206:       args: -mat_ordering -display_matrices -nox
207:       filter: grep -v " MPI process"

209:    test:
210:       suffix: 2
211:       args: -mat_ordering -display_matrices -nox -lu

213:    test:
214:       suffix: 3
215:       args: -mat_ordering -lu -triangular_solve

217:    test:
218:       suffix: 4

220:    test:
221:       suffix: 5
222:       args: -lu

224:    test:
225:       suffix: 6
226:       args: -lu -triangular_solve
227:       output_file: output/ex30_3.out

229: TEST*/