Actual source code: ex74.c

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

  3: #include <petscmat.h>

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

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 23:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 24:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

169:   PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY));
170:   PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY));

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

181:   /* Test MatDuplicate() */
182:   PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1));
183:   PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB));
184:   PetscCall(MatEqual(sA, sB, &equal));
185:   PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDuplicate()");

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

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

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

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

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

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

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

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

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

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

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

280:   /* Test MatMatMult() for sbaij and dense matrices */
281:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, 5 * n, NULL, &B));
282:   PetscCall(MatSetRandom(B, rdm));
283:   PetscCall(MatMatMult(sA, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
284:   PetscCall(MatMatMultEqual(sA, B, C, 5 * n, &equal));
285:   PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()");
286:   PetscCall(MatDestroy(&C));
287:   PetscCall(MatDestroy(&B));

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

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

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

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

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

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

323:     PetscCall(MatMult(sB, x, b));

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

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

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

361:   PetscCall(ISDestroy(&perm));

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

373:   PetscCall(PetscFinalize());
374:   return 0;
375: }

377: /*TEST

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

382: TEST*/