Actual source code: ex76.c

  1: static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec           x, y, b;
  8:   Mat           A;      /* linear system matrix */
  9:   Mat           sA, sC; /* symmetric part of the matrices */
 10:   PetscInt      n, mbs = 16, bs = 1, nz = 3, prob = 1, i, j, col[3], block, row, Ii, J, n1, lvl;
 11:   PetscMPIInt   size;
 12:   PetscReal     norm2;
 13:   PetscScalar   neg_one = -1.0, four = 4.0, value[3];
 14:   IS            perm, cperm;
 15:   PetscRandom   rdm;
 16:   PetscBool     reorder = PETSC_FALSE, displ = PETSC_FALSE;
 17:   MatFactorInfo factinfo;
 18:   PetscBool     equal;
 19:   PetscBool     TestAIJ = PETSC_FALSE, TestBAIJ = PETSC_TRUE;
 20:   PetscInt      TestShift = 0;

 22:   PetscFunctionBeginUser;
 23:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 24:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 25:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
 28:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-reorder", &reorder, NULL));
 29:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testaij", &TestAIJ, NULL));
 30:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-testShift", &TestShift, NULL));
 31:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL));

 33:   n = mbs * bs;
 34:   if (TestAIJ) { /* A is in aij format */
 35:     PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, n, n, nz, NULL, &A));
 36:     TestBAIJ = PETSC_FALSE;
 37:   } else { /* A is in baij format */
 38:     PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A));
 39:     TestAIJ = PETSC_FALSE;
 40:   }

 42:   /* Assemble matrix */
 43:   if (bs == 1) {
 44:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL));
 45:     if (prob == 1) { /* tridiagonal matrix */
 46:       value[0] = -1.0;
 47:       value[1] = 2.0;
 48:       value[2] = -1.0;
 49:       for (i = 1; i < n - 1; i++) {
 50:         col[0] = i - 1;
 51:         col[1] = i;
 52:         col[2] = i + 1;
 53:         PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 54:       }
 55:       i      = n - 1;
 56:       col[0] = 0;
 57:       col[1] = n - 2;
 58:       col[2] = n - 1;

 60:       value[0] = 0.1;
 61:       value[1] = -1;
 62:       value[2] = 2;
 63:       PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));

 65:       i      = 0;
 66:       col[0] = 0;
 67:       col[1] = 1;
 68:       col[2] = n - 1;

 70:       value[0] = 2.0;
 71:       value[1] = -1.0;
 72:       value[2] = 0.1;
 73:       PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 74:     } else if (prob == 2) { /* matrix for the five point stencil */
 75:       n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
 76:       PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!");
 77:       for (i = 0; i < n1; i++) {
 78:         for (j = 0; j < n1; j++) {
 79:           Ii = j + n1 * i;
 80:           if (i > 0) {
 81:             J = Ii - n1;
 82:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 83:           }
 84:           if (i < n1 - 1) {
 85:             J = Ii + n1;
 86:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 87:           }
 88:           if (j > 0) {
 89:             J = Ii - 1;
 90:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 91:           }
 92:           if (j < n1 - 1) {
 93:             J = Ii + 1;
 94:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 95:           }
 96:           PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES));
 97:         }
 98:       }
 99:     }
100:   } else { /* bs > 1 */
101:     for (block = 0; block < n / bs; block++) {
102:       /* diagonal blocks */
103:       value[0] = -1.0;
104:       value[1] = 4.0;
105:       value[2] = -1.0;
106:       for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
107:         col[0] = i - 1;
108:         col[1] = i;
109:         col[2] = i + 1;
110:         PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
111:       }
112:       i      = bs - 1 + block * bs;
113:       col[0] = bs - 2 + block * bs;
114:       col[1] = bs - 1 + block * bs;

116:       value[0] = -1.0;
117:       value[1] = 4.0;
118:       PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));

120:       i      = 0 + block * bs;
121:       col[0] = 0 + block * bs;
122:       col[1] = 1 + block * bs;

124:       value[0] = 4.0;
125:       value[1] = -1.0;
126:       PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
127:     }
128:     /* off-diagonal blocks */
129:     value[0] = -1.0;
130:     for (i = 0; i < (n / bs - 1) * bs; i++) {
131:       col[0] = i + bs;
132:       PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
133:       col[0] = i;
134:       row    = i + bs;
135:       PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
136:     }
137:   }

139:   if (TestShift) {
140:     /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
141:     for (i = 0; i < bs; i++) {
142:       row      = i;
143:       col[0]   = i;
144:       value[0] = 0.0;
145:       PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
146:     }
147:   }

149:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
150:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

152:   /* Test MatConvert */
153:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
154:   PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA));
155:   PetscCall(MatMultEqual(A, sA, 20, &equal));
156:   PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_USER, "A != sA");

158:   /* Test MatGetOwnershipRange() */
159:   PetscCall(MatGetOwnershipRange(A, &Ii, &J));
160:   PetscCall(MatGetOwnershipRange(sA, &i, &j));
161:   PetscCheck(i == Ii && j == J, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetOwnershipRange() in MatSBAIJ format");

163:   /* Vectors */
164:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
165:   PetscCall(PetscRandomSetFromOptions(rdm));
166:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
167:   PetscCall(VecDuplicate(x, &b));
168:   PetscCall(VecDuplicate(x, &y));
169:   PetscCall(VecSetRandom(x, rdm));

171:   /* Test MatReordering() - not work on sbaij matrix */
172:   if (reorder) {
173:     PetscCall(MatGetOrdering(A, MATORDERINGRCM, &perm, &cperm));
174:   } else {
175:     PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &cperm));
176:   }
177:   PetscCall(ISDestroy(&cperm));

179:   /* initialize factinfo */
180:   PetscCall(MatFactorInfoInitialize(&factinfo));
181:   if (TestShift == 1) {
182:     factinfo.shifttype   = (PetscReal)MAT_SHIFT_NONZERO;
183:     factinfo.shiftamount = 0.1;
184:   } else if (TestShift == 2) {
185:     factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
186:   }

188:   /* Test MatCholeskyFactor(), MatICCFactor() */
189:   /*------------------------------------------*/
190:   /* Test aij matrix A */
191:   if (TestAIJ) {
192:     if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "AIJ: \n"));
193:     i = 0;
194:     for (lvl = -1; lvl < 10; lvl++) {
195:       if (lvl == -1) { /* Cholesky factor */
196:         factinfo.fill = 5.0;

198:         PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
199:         PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo));
200:       } else { /* incomplete Cholesky factor */
201:         factinfo.fill   = 5.0;
202:         factinfo.levels = lvl;

204:         PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
205:         PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo));
206:       }
207:       PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo));

209:       PetscCall(MatMult(A, x, b));
210:       PetscCall(MatSolve(sC, b, y));
211:       PetscCall(MatDestroy(&sC));

213:       /* Check the residual */
214:       PetscCall(VecAXPY(y, neg_one, x));
215:       PetscCall(VecNorm(y, NORM_2, &norm2));

217:       if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
218:     }
219:   }

221:   /* Test baij matrix A */
222:   if (TestBAIJ) {
223:     if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "BAIJ: \n"));
224:     i = 0;
225:     for (lvl = -1; lvl < 10; lvl++) {
226:       if (lvl == -1) { /* Cholesky factor */
227:         factinfo.fill = 5.0;

229:         PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
230:         PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo));
231:       } else { /* incomplete Cholesky factor */
232:         factinfo.fill   = 5.0;
233:         factinfo.levels = lvl;

235:         PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
236:         PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo));
237:       }
238:       PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo));

240:       PetscCall(MatMult(A, x, b));
241:       PetscCall(MatSolve(sC, b, y));
242:       PetscCall(MatDestroy(&sC));

244:       /* Check the residual */
245:       PetscCall(VecAXPY(y, neg_one, x));
246:       PetscCall(VecNorm(y, NORM_2, &norm2));
247:       if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
248:     }
249:   }

251:   /* Test sbaij matrix sA */
252:   if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SBAIJ: \n"));
253:   i = 0;
254:   for (lvl = -1; lvl < 10; lvl++) {
255:     if (lvl == -1) { /* Cholesky factor */
256:       factinfo.fill = 5.0;

258:       PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
259:       PetscCall(MatCholeskyFactorSymbolic(sC, sA, perm, &factinfo));
260:     } else { /* incomplete Cholesky factor */
261:       factinfo.fill   = 5.0;
262:       factinfo.levels = lvl;

264:       PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
265:       PetscCall(MatICCFactorSymbolic(sC, sA, perm, &factinfo));
266:     }
267:     PetscCall(MatCholeskyFactorNumeric(sC, sA, &factinfo));

269:     if (lvl == 0 && bs == 1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
270:       /*
271:         Mat B;
272:         PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&B));
273:         PetscCall(MatICCFactor(B,perm,&factinfo));
274:         PetscCall(MatEqual(sC,B,&equal));
275:         if (!equal) {
276:           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
277:         }
278:         PetscCall(MatDestroy(&B));
279:       */
280:     }

282:     PetscCall(MatMult(sA, x, b));
283:     PetscCall(MatSolve(sC, b, y));

285:     /* Test MatSolves() */
286:     if (bs == 1) {
287:       Vecs xx, bb;
288:       PetscCall(VecsCreateSeq(PETSC_COMM_SELF, n, 4, &xx));
289:       PetscCall(VecsDuplicate(xx, &bb));
290:       PetscCall(MatSolves(sC, bb, xx));
291:       PetscCall(VecsDestroy(xx));
292:       PetscCall(VecsDestroy(bb));
293:     }
294:     PetscCall(MatDestroy(&sC));

296:     /* Check the residual */
297:     PetscCall(VecAXPY(y, neg_one, x));
298:     PetscCall(VecNorm(y, NORM_2, &norm2));
299:     if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
300:   }

302:   PetscCall(ISDestroy(&perm));
303:   PetscCall(MatDestroy(&A));
304:   PetscCall(MatDestroy(&sA));
305:   PetscCall(VecDestroy(&x));
306:   PetscCall(VecDestroy(&y));
307:   PetscCall(VecDestroy(&b));
308:   PetscCall(PetscRandomDestroy(&rdm));

310:   PetscCall(PetscFinalize());
311:   return 0;
312: }

314: /*TEST

316:    test:
317:       args: -bs {{1 2 3 4 5 6 7 8}}

319:    test:
320:       suffix: 3
321:       args: -testaij
322:       output_file: output/ex76_1.out

324: TEST*/