Actual source code: ex74.c


  2: static char help[] = "Tests the various sequential routines in MATSEQSBAIJ format.\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   PetscMPIInt   size;
  9:   Vec           x, y, b, s1, s2;
 10:   Mat           A;                     /* linear system matrix */
 11:   Mat           sA, sB, sFactor, B, C; /* symmetric matrices */
 12:   PetscInt      n, mbs = 16, bs = 1, nz = 3, prob = 1, i, j, k1, k2, col[3], lf, block, row, Ii, J, n1, inc;
 13:   PetscReal     norm1, norm2, rnorm, tol = 10 * PETSC_SMALL;
 14:   PetscScalar   neg_one = -1.0, four = 4.0, value[3];
 15:   IS            perm, iscol;
 16:   PetscRandom   rdm;
 17:   PetscBool     doIcc = PETSC_TRUE, equal;
 18:   MatInfo       minfo1, minfo2;
 19:   MatFactorInfo factinfo;
 20:   MatType       type;

 23:   PetscInitialize(&argc, &args, (char *)0, help);
 24:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 26:   PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
 27:   PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL);

 29:   n = mbs * bs;
 30:   MatCreate(PETSC_COMM_SELF, &A);
 31:   MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
 32:   MatSetType(A, MATSEQBAIJ);
 33:   MatSetFromOptions(A);
 34:   MatSeqBAIJSetPreallocation(A, bs, nz, NULL);

 36:   MatCreate(PETSC_COMM_SELF, &sA);
 37:   MatSetSizes(sA, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
 38:   MatSetType(sA, MATSEQSBAIJ);
 39:   MatSetFromOptions(sA);
 40:   MatGetType(sA, &type);
 41:   PetscObjectTypeCompare((PetscObject)sA, MATSEQSBAIJ, &doIcc);
 42:   MatSeqSBAIJSetPreallocation(sA, bs, nz, NULL);
 43:   MatSetOption(sA, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);

 45:   /* Test MatGetOwnershipRange() */
 46:   MatGetOwnershipRange(A, &Ii, &J);
 47:   MatGetOwnershipRange(sA, &i, &j);
 48:   if (i - Ii || j - J) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n");

 50:   /* Assemble matrix */
 51:   if (bs == 1) {
 52:     PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL);
 53:     if (prob == 1) { /* tridiagonal matrix */
 54:       value[0] = -1.0;
 55:       value[1] = 2.0;
 56:       value[2] = -1.0;
 57:       for (i = 1; i < n - 1; i++) {
 58:         col[0] = i - 1;
 59:         col[1] = i;
 60:         col[2] = i + 1;
 61:         MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
 62:         MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);
 63:       }
 64:       i      = n - 1;
 65:       col[0] = 0;
 66:       col[1] = n - 2;
 67:       col[2] = n - 1;

 69:       value[0] = 0.1;
 70:       value[1] = -1;
 71:       value[2] = 2;

 73:       MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
 74:       MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);

 76:       i        = 0;
 77:       col[0]   = n - 1;
 78:       col[1]   = 1;
 79:       col[2]   = 0;
 80:       value[0] = 0.1;
 81:       value[1] = -1.0;
 82:       value[2] = 2;

 84:       MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
 85:       MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);

 87:     } else if (prob == 2) { /* matrix for the five point stencil */
 88:       n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
 90:       for (i = 0; i < n1; i++) {
 91:         for (j = 0; j < n1; j++) {
 92:           Ii = j + n1 * i;
 93:           if (i > 0) {
 94:             J = Ii - n1;
 95:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 96:             MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 97:           }
 98:           if (i < n1 - 1) {
 99:             J = Ii + n1;
100:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
101:             MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
102:           }
103:           if (j > 0) {
104:             J = Ii - 1;
105:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
106:             MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
107:           }
108:           if (j < n1 - 1) {
109:             J = Ii + 1;
110:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
111:             MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
112:           }
113:           MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES);
114:           MatSetValues(sA, 1, &Ii, 1, &Ii, &four, INSERT_VALUES);
115:         }
116:       }
117:     }

119:   } else { /* bs > 1 */
120:     for (block = 0; block < n / bs; block++) {
121:       /* diagonal blocks */
122:       value[0] = -1.0;
123:       value[1] = 4.0;
124:       value[2] = -1.0;
125:       for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
126:         col[0] = i - 1;
127:         col[1] = i;
128:         col[2] = i + 1;
129:         MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
130:         MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);
131:       }
132:       i      = bs - 1 + block * bs;
133:       col[0] = bs - 2 + block * bs;
134:       col[1] = bs - 1 + block * bs;

136:       value[0] = -1.0;
137:       value[1] = 4.0;

139:       MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
140:       MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES);

142:       i      = 0 + block * bs;
143:       col[0] = 0 + block * bs;
144:       col[1] = 1 + block * bs;

146:       value[0] = 4.0;
147:       value[1] = -1.0;

149:       MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
150:       MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES);
151:     }
152:     /* off-diagonal blocks */
153:     value[0] = -1.0;
154:     for (i = 0; i < (n / bs - 1) * bs; i++) {
155:       col[0] = i + bs;

157:       MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
158:       MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES);

160:       col[0] = i;
161:       row    = i + bs;

163:       MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
164:       MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES);
165:     }
166:   }
167:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
168:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

170:   MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY);
171:   MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY);

173:   /* Test MatGetInfo() of A and sA */
174:   MatGetInfo(A, MAT_LOCAL, &minfo1);
175:   MatGetInfo(sA, MAT_LOCAL, &minfo2);
176:   i  = (int)(minfo1.nz_used - minfo2.nz_used);
177:   j  = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
178:   k1 = (int)(minfo1.nz_allocated - minfo1.nz_used);
179:   k2 = (int)(minfo2.nz_allocated - minfo2.nz_used);
180:   if (i < 0 || j < 0 || k1 < 0 || k2 < 0) PetscPrintf(PETSC_COMM_SELF, "Error (compare A and sA): MatGetInfo()\n");

182:   /* Test MatDuplicate() */
183:   MatNorm(A, NORM_FROBENIUS, &norm1);
184:   MatDuplicate(sA, MAT_COPY_VALUES, &sB);
185:   MatEqual(sA, sB, &equal);

188:   /* Test MatNorm() */
189:   MatNorm(A, NORM_FROBENIUS, &norm1);
190:   MatNorm(sB, NORM_FROBENIUS, &norm2);
191:   rnorm = PetscAbsReal(norm1 - norm2) / norm2;
192:   if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2);
193:   MatNorm(A, NORM_INFINITY, &norm1);
194:   MatNorm(sB, NORM_INFINITY, &norm2);
195:   rnorm = PetscAbsReal(norm1 - norm2) / norm2;
196:   if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2);
197:   MatNorm(A, NORM_1, &norm1);
198:   MatNorm(sB, NORM_1, &norm2);
199:   rnorm = PetscAbsReal(norm1 - norm2) / norm2;
200:   if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2);

202:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
203:   MatGetInfo(A, MAT_LOCAL, &minfo1);
204:   MatGetInfo(sB, MAT_LOCAL, &minfo2);
205:   i  = (int)(minfo1.nz_used - minfo2.nz_used);
206:   j  = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
207:   k1 = (int)(minfo1.nz_allocated - minfo1.nz_used);
208:   k2 = (int)(minfo2.nz_allocated - minfo2.nz_used);
209:   if (i < 0 || j < 0 || k1 < 0 || k2 < 0) PetscPrintf(PETSC_COMM_SELF, "Error(compare A and sB): MatGetInfo()\n");

211:   MatGetSize(A, &Ii, &J);
212:   MatGetSize(sB, &i, &j);
213:   if (i - Ii || j - J) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n");

215:   MatGetBlockSize(A, &Ii);
216:   MatGetBlockSize(sB, &i);
217:   if (i - Ii) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n");

219:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
220:   PetscRandomSetFromOptions(rdm);
221:   VecCreateSeq(PETSC_COMM_SELF, n, &x);
222:   VecDuplicate(x, &s1);
223:   VecDuplicate(x, &s2);
224:   VecDuplicate(x, &y);
225:   VecDuplicate(x, &b);
226:   VecSetRandom(x, rdm);

228:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
229: #if !defined(PETSC_USE_COMPLEX)
230:   /* Scaling matrix with complex numbers results non-spd matrix,
231:      causing crash of MatForwardSolve() and MatBackwardSolve() */
232:   MatDiagonalScale(A, x, x);
233:   MatDiagonalScale(sB, x, x);
234:   MatMultEqual(A, sB, 10, &equal);

237:   MatGetDiagonal(A, s1);
238:   MatGetDiagonal(sB, s2);
239:   VecAXPY(s2, neg_one, s1);
240:   VecNorm(s2, NORM_1, &norm1);
241:   if (norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal(), ||s1-s2||=%g\n", (double)norm1);

243:   {
244:     PetscScalar alpha = 0.1;
245:     MatScale(A, alpha);
246:     MatScale(sB, alpha);
247:   }
248: #endif

250:   /* Test MatGetRowMaxAbs() */
251:   MatGetRowMaxAbs(A, s1, NULL);
252:   MatGetRowMaxAbs(sB, s2, NULL);
253:   VecNorm(s1, NORM_1, &norm1);
254:   VecNorm(s2, NORM_1, &norm2);
255:   norm1 -= norm2;
256:   if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatGetRowMaxAbs() \n");

258:   /* Test MatMult() */
259:   for (i = 0; i < 40; i++) {
260:     VecSetRandom(x, rdm);
261:     MatMult(A, x, s1);
262:     MatMult(sB, x, s2);
263:     VecNorm(s1, NORM_1, &norm1);
264:     VecNorm(s2, NORM_1, &norm2);
265:     norm1 -= norm2;
266:     if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), norm1-norm2: %g\n", (double)norm1);
267:   }

269:   /* MatMultAdd() */
270:   for (i = 0; i < 40; i++) {
271:     VecSetRandom(x, rdm);
272:     VecSetRandom(y, rdm);
273:     MatMultAdd(A, x, y, s1);
274:     MatMultAdd(sB, x, y, s2);
275:     VecNorm(s1, NORM_1, &norm1);
276:     VecNorm(s2, NORM_1, &norm2);
277:     norm1 -= norm2;
278:     if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), norm1-norm2: %g\n", (double)norm1);
279:   }

281:   /* Test MatMatMult() for sbaij and dense matrices */
282:   MatCreateSeqDense(PETSC_COMM_SELF, n, 5 * n, NULL, &B);
283:   MatSetRandom(B, rdm);
284:   MatMatMult(sA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
285:   MatMatMultEqual(sA, B, C, 5 * n, &equal);
287:   MatDestroy(&C);
288:   MatDestroy(&B);

290:   /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
291:   MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iscol);
292:   ISDestroy(&iscol);
293:   norm1 = tol;
294:   inc   = bs;

296:   /* initialize factinfo */
297:   PetscMemzero(&factinfo, sizeof(MatFactorInfo));

299:   for (lf = -1; lf < 10; lf += inc) {
300:     if (lf == -1) { /* Cholesky factor of sB (duplicate sA) */
301:       factinfo.fill = 5.0;

303:       MatGetFactor(sB, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sFactor);
304:       MatCholeskyFactorSymbolic(sFactor, sB, perm, &factinfo);
305:     } else if (!doIcc) break;
306:     else { /* incomplete Cholesky factor */ factinfo.fill = 5.0;
307:       factinfo.levels                                     = lf;

309:       MatGetFactor(sB, MATSOLVERPETSC, MAT_FACTOR_ICC, &sFactor);
310:       MatICCFactorSymbolic(sFactor, sB, perm, &factinfo);
311:     }
312:     MatCholeskyFactorNumeric(sFactor, sB, &factinfo);
313:     /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */

315:     /* test MatGetDiagonal on numeric factor */
316:     /*
317:     if (lf == -1) {
318:       MatGetDiagonal(sFactor,s1);
319:       printf(" in ex74.c, diag: \n");
320:       VecView(s1,PETSC_VIEWER_STDOUT_SELF);
321:     }
322:     */

324:     MatMult(sB, x, b);

326:     /* test MatForwardSolve() and MatBackwardSolve() */
327:     if (lf == -1) {
328:       MatForwardSolve(sFactor, b, s1);
329:       MatBackwardSolve(sFactor, s1, s2);
330:       VecAXPY(s2, neg_one, x);
331:       VecNorm(s2, NORM_2, &norm2);
332:       if (10 * norm1 < norm2) PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%" PetscInt_FMT "\n", (double)norm2, bs);
333:     }

335:     /* test MatSolve() */
336:     MatSolve(sFactor, b, y);
337:     MatDestroy(&sFactor);
338:     /* Check the error */
339:     VecAXPY(y, neg_one, x);
340:     VecNorm(y, NORM_2, &norm2);
341:     if (10 * norm1 < norm2 && lf - inc != -1) PetscPrintf(PETSC_COMM_SELF, "lf=%" PetscInt_FMT ", %" PetscInt_FMT ", Norm of error=%g, %g\n", lf - inc, lf, (double)norm1, (double)norm2);
342:     norm1 = norm2;
343:     if (norm2 < tol && lf != -1) break;
344:   }

346: #if defined(PETSC_HAVE_MUMPS)
347:   MatGetFactor(sA, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &sFactor);
348:   MatCholeskyFactorSymbolic(sFactor, sA, NULL, NULL);
349:   MatCholeskyFactorNumeric(sFactor, sA, NULL);
350:   for (i = 0; i < 10; i++) {
351:     VecSetRandom(b, rdm);
352:     MatSolve(sFactor, b, y);
353:     /* Check the error */
354:     MatMult(sA, y, x);
355:     VecAXPY(x, neg_one, b);
356:     VecNorm(x, NORM_2, &norm2);
357:     if (norm2 > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolve(), norm2: %g\n", (double)norm2);
358:   }
359:   MatDestroy(&sFactor);
360: #endif

362:   ISDestroy(&perm);

364:   MatDestroy(&A);
365:   MatDestroy(&sB);
366:   MatDestroy(&sA);
367:   VecDestroy(&x);
368:   VecDestroy(&y);
369:   VecDestroy(&s1);
370:   VecDestroy(&s2);
371:   VecDestroy(&b);
372:   PetscRandomDestroy(&rdm);

374:   PetscFinalize();
375:   return 0;
376: }

378: /*TEST

380:    test:
381:       args: -bs {{1 2 3 4 5 6 7 8}}

383: TEST*/