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;

 24:   PetscInitialize(&argc, &args, (char *)0, help);
 25:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 27:   PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
 28:   PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL);
 29:   PetscOptionsGetBool(NULL, NULL, "-reorder", &reorder, NULL);
 30:   PetscOptionsGetBool(NULL, NULL, "-testaij", &TestAIJ, NULL);
 31:   PetscOptionsGetInt(NULL, NULL, "-testShift", &TestShift, NULL);
 32:   PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL);

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

 43:   /* Assemble matrix */
 44:   if (bs == 1) {
 45:     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:         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:       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:       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);
 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:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 84:           }
 85:           if (i < n1 - 1) {
 86:             J = Ii + n1;
 87:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 88:           }
 89:           if (j > 0) {
 90:             J = Ii - 1;
 91:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 92:           }
 93:           if (j < n1 - 1) {
 94:             J = Ii + 1;
 95:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 96:           }
 97:           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:         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:       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:       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:       MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
134:       col[0] = i;
135:       row    = i + bs;
136:       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:       MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
147:     }
148:   }

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

153:   /* Test MatConvert */
154:   MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
155:   MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA);
156:   MatMultEqual(A, sA, 20, &equal);

159:   /* Test MatGetOwnershipRange() */
160:   MatGetOwnershipRange(A, &Ii, &J);
161:   MatGetOwnershipRange(sA, &i, &j);

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

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

180:   /* initialize factinfo */
181:   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) 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:         MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC);
200:         MatCholeskyFactorSymbolic(sC, A, perm, &factinfo);
201:       } else { /* incomplete Cholesky factor */
202:         factinfo.fill   = 5.0;
203:         factinfo.levels = lvl;

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

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

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

218:       if (displ) 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) 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:         MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC);
231:         MatCholeskyFactorSymbolic(sC, A, perm, &factinfo);
232:       } else { /* incomplete Cholesky factor */
233:         factinfo.fill   = 5.0;
234:         factinfo.levels = lvl;

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

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

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

252:   /* Test sbaij matrix sA */
253:   if (displ) 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:       MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC);
260:       MatCholeskyFactorSymbolic(sC, sA, perm, &factinfo);
261:     } else { /* incomplete Cholesky factor */
262:       factinfo.fill   = 5.0;
263:       factinfo.levels = lvl;

265:       MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC);
266:       MatICCFactorSymbolic(sC, sA, perm, &factinfo);
267:     }
268:     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:         MatDuplicate(sA,MAT_COPY_VALUES,&B);
274:         MatICCFactor(B,perm,&factinfo);
275:         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:         MatDestroy(&B);
280:       */
281:     }

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

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

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

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

311:   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*/