Actual source code: ex76.c


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

  4: #include <petscmat.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

315: /*TEST

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

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

325: TEST*/