Actual source code: multequal.c


  2: #include <petsc/private/matimpl.h>

  4: static PetscErrorCode MatMultEqual_Private(Mat A, Mat B, PetscInt n, PetscBool *flg, PetscInt t, PetscBool add)
  5: {
  6:   Vec         Ax = NULL, Bx = NULL, s1 = NULL, s2 = NULL, Ay = NULL, By = NULL;
  7:   PetscRandom rctx;
  8:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
  9:   PetscInt    am, an, bm, bn, k;
 10:   PetscScalar none   = -1.0;
 11:   const char *sops[] = {"MatMult", "MatMultAdd", "MatMultTranspose", "MatMultTransposeAdd", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd"};
 12:   const char *sop;

 21:   MatGetLocalSize(A, &am, &an);
 22:   MatGetLocalSize(B, &bm, &bn);
 24:   sop = sops[(add ? 1 : 0) + 2 * t]; /* t = 0 => no transpose, t = 1 => transpose, t = 2 => Hermitian transpose */
 25:   PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx);
 26:   PetscRandomSetFromOptions(rctx);
 27:   if (t) {
 28:     MatCreateVecs(A, &s1, &Ax);
 29:     MatCreateVecs(B, &s2, &Bx);
 30:   } else {
 31:     MatCreateVecs(A, &Ax, &s1);
 32:     MatCreateVecs(B, &Bx, &s2);
 33:   }
 34:   if (add) {
 35:     VecDuplicate(s1, &Ay);
 36:     VecDuplicate(s2, &By);
 37:   }

 39:   *flg = PETSC_TRUE;
 40:   for (k = 0; k < n; k++) {
 41:     VecSetRandom(Ax, rctx);
 42:     VecCopy(Ax, Bx);
 43:     if (add) {
 44:       VecSetRandom(Ay, rctx);
 45:       VecCopy(Ay, By);
 46:     }
 47:     if (t == 1) {
 48:       if (add) {
 49:         MatMultTransposeAdd(A, Ax, Ay, s1);
 50:         MatMultTransposeAdd(B, Bx, By, s2);
 51:       } else {
 52:         MatMultTranspose(A, Ax, s1);
 53:         MatMultTranspose(B, Bx, s2);
 54:       }
 55:     } else if (t == 2) {
 56:       if (add) {
 57:         MatMultHermitianTransposeAdd(A, Ax, Ay, s1);
 58:         MatMultHermitianTransposeAdd(B, Bx, By, s2);
 59:       } else {
 60:         MatMultHermitianTranspose(A, Ax, s1);
 61:         MatMultHermitianTranspose(B, Bx, s2);
 62:       }
 63:     } else {
 64:       if (add) {
 65:         MatMultAdd(A, Ax, Ay, s1);
 66:         MatMultAdd(B, Bx, By, s2);
 67:       } else {
 68:         MatMult(A, Ax, s1);
 69:         MatMult(B, Bx, s2);
 70:       }
 71:     }
 72:     VecNorm(s2, NORM_INFINITY, &r2);
 73:     if (r2 < tol) {
 74:       VecNorm(s1, NORM_INFINITY, &r1);
 75:     } else {
 76:       VecAXPY(s2, none, s1);
 77:       VecNorm(s2, NORM_INFINITY, &r1);
 78:       r1 /= r2;
 79:     }
 80:     if (r1 > tol) {
 81:       *flg = PETSC_FALSE;
 82:       PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1);
 83:       break;
 84:     }
 85:   }
 86:   PetscRandomDestroy(&rctx);
 87:   VecDestroy(&Ax);
 88:   VecDestroy(&Bx);
 89:   VecDestroy(&Ay);
 90:   VecDestroy(&By);
 91:   VecDestroy(&s1);
 92:   VecDestroy(&s2);
 93:   return 0;
 94: }

 96: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
 97: {
 98:   Vec         Ax, Bx, Cx, s1, s2, s3;
 99:   PetscRandom rctx;
100:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
101:   PetscInt    am, an, bm, bn, cm, cn, k;
102:   PetscScalar none   = -1.0;
103:   const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
104:   const char *sop;

115:   MatGetLocalSize(A, &am, &an);
116:   MatGetLocalSize(B, &bm, &bn);
117:   MatGetLocalSize(C, &cm, &cn);
118:   if (At) {
119:     PetscInt tt = an;
120:     an          = am;
121:     am          = tt;
122:   };
123:   if (Bt) {
124:     PetscInt tt = bn;
125:     bn          = bm;
126:     bm          = tt;
127:   };

130:   sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
131:   PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx);
132:   PetscRandomSetFromOptions(rctx);
133:   if (Bt) {
134:     MatCreateVecs(B, &s1, &Bx);
135:   } else {
136:     MatCreateVecs(B, &Bx, &s1);
137:   }
138:   if (At) {
139:     MatCreateVecs(A, &s2, &Ax);
140:   } else {
141:     MatCreateVecs(A, &Ax, &s2);
142:   }
143:   MatCreateVecs(C, &Cx, &s3);

145:   *flg = PETSC_TRUE;
146:   for (k = 0; k < n; k++) {
147:     VecSetRandom(Bx, rctx);
148:     if (Bt) {
149:       MatMultTranspose(B, Bx, s1);
150:     } else {
151:       MatMult(B, Bx, s1);
152:     }
153:     VecCopy(s1, Ax);
154:     if (At) {
155:       MatMultTranspose(A, Ax, s2);
156:     } else {
157:       MatMult(A, Ax, s2);
158:     }
159:     VecCopy(Bx, Cx);
160:     MatMult(C, Cx, s3);

162:     VecNorm(s2, NORM_INFINITY, &r2);
163:     if (r2 < tol) {
164:       VecNorm(s3, NORM_INFINITY, &r1);
165:     } else {
166:       VecAXPY(s2, none, s3);
167:       VecNorm(s2, NORM_INFINITY, &r1);
168:       r1 /= r2;
169:     }
170:     if (r1 > tol) {
171:       *flg = PETSC_FALSE;
172:       PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1);
173:       break;
174:     }
175:   }
176:   PetscRandomDestroy(&rctx);
177:   VecDestroy(&Ax);
178:   VecDestroy(&Bx);
179:   VecDestroy(&Cx);
180:   VecDestroy(&s1);
181:   VecDestroy(&s2);
182:   VecDestroy(&s3);
183:   return 0;
184: }

186: /*@
187:    MatMultEqual - Compares matrix-vector products of two matrices.

189:    Collective

191:    Input Parameters:
192: +  A - the first matrix
193: .  B - the second matrix
194: -  n - number of random vectors to be tested

196:    Output Parameter:
197: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

199:    Level: intermediate

201: @*/
202: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
203: {
204:   MatMultEqual_Private(A, B, n, flg, 0, PETSC_FALSE);
205:   return 0;
206: }

208: /*@
209:    MatMultAddEqual - Compares matrix-vector products of two matrices.

211:    Collective

213:    Input Parameters:
214: +  A - the first matrix
215: .  B - the second matrix
216: -  n - number of random vectors to be tested

218:    Output Parameter:
219: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

221:    Level: intermediate

223: @*/
224: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
225: {
226:   MatMultEqual_Private(A, B, n, flg, 0, PETSC_TRUE);
227:   return 0;
228: }

230: /*@
231:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

233:    Collective

235:    Input Parameters:
236: +  A - the first matrix
237: .  B - the second matrix
238: -  n - number of random vectors to be tested

240:    Output Parameter:
241: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

243:    Level: intermediate

245: @*/
246: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
247: {
248:   MatMultEqual_Private(A, B, n, flg, 1, PETSC_FALSE);
249:   return 0;
250: }

252: /*@
253:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

255:    Collective

257:    Input Parameters:
258: +  A - the first matrix
259: .  B - the second matrix
260: -  n - number of random vectors to be tested

262:    Output Parameter:
263: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

265:    Level: intermediate

267: @*/
268: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
269: {
270:   MatMultEqual_Private(A, B, n, flg, 1, PETSC_TRUE);
271:   return 0;
272: }

274: /*@
275:    MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.

277:    Collective

279:    Input Parameters:
280: +  A - the first matrix
281: .  B - the second matrix
282: -  n - number of random vectors to be tested

284:    Output Parameter:
285: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

287:    Level: intermediate

289: @*/
290: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
291: {
292:   MatMultEqual_Private(A, B, n, flg, 2, PETSC_FALSE);
293:   return 0;
294: }

296: /*@
297:    MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.

299:    Collective

301:    Input Parameters:
302: +  A - the first matrix
303: .  B - the second matrix
304: -  n - number of random vectors to be tested

306:    Output Parameter:
307: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

309:    Level: intermediate

311: @*/
312: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
313: {
314:   MatMultEqual_Private(A, B, n, flg, 2, PETSC_TRUE);
315:   return 0;
316: }

318: /*@
319:    MatMatMultEqual - Test A*B*x = C*x for n random vector x

321:    Collective

323:    Input Parameters:
324: +  A - the first matrix
325: .  B - the second matrix
326: .  C - the third matrix
327: -  n - number of random vectors to be tested

329:    Output Parameter:
330: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

332:    Level: intermediate

334: @*/
335: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
336: {
337:   MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE);
338:   return 0;
339: }

341: /*@
342:    MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x

344:    Collective

346:    Input Parameters:
347: +  A - the first matrix
348: .  B - the second matrix
349: .  C - the third matrix
350: -  n - number of random vectors to be tested

352:    Output Parameter:
353: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

355:    Level: intermediate

357: @*/
358: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
359: {
360:   MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE);
361:   return 0;
362: }

364: /*@
365:    MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x

367:    Collective

369:    Input Parameters:
370: +  A - the first matrix
371: .  B - the second matrix
372: .  C - the third matrix
373: -  n - number of random vectors to be tested

375:    Output Parameter:
376: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

378:    Level: intermediate

380: @*/
381: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
382: {
383:   MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE);
384:   return 0;
385: }

387: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
388: {
389:   Vec         x, v1, v2, v3, v4, Cx, Bx;
390:   PetscReal   norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
391:   PetscInt    i, am, an, bm, bn, cm, cn;
392:   PetscRandom rdm;
393:   PetscScalar none = -1.0;

395:   MatGetLocalSize(A, &am, &an);
396:   MatGetLocalSize(B, &bm, &bn);
397:   if (rart) {
398:     PetscInt t = bm;
399:     bm         = bn;
400:     bn         = t;
401:   }
402:   MatGetLocalSize(C, &cm, &cn);

405:   /* Create left vector of A: v2 */
406:   MatCreateVecs(A, &Bx, &v2);

408:   /* Create right vectors of B: x, v3, v4 */
409:   if (rart) {
410:     MatCreateVecs(B, &v1, &x);
411:   } else {
412:     MatCreateVecs(B, &x, &v1);
413:   }
414:   VecDuplicate(x, &v3);

416:   MatCreateVecs(C, &Cx, &v4);
417:   PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
418:   PetscRandomSetFromOptions(rdm);

420:   *flg = PETSC_TRUE;
421:   for (i = 0; i < n; i++) {
422:     VecSetRandom(x, rdm);
423:     VecCopy(x, Cx);
424:     MatMult(C, Cx, v4); /* v4 = C*x   */
425:     if (rart) {
426:       MatMultTranspose(B, x, v1);
427:     } else {
428:       MatMult(B, x, v1);
429:     }
430:     VecCopy(v1, Bx);
431:     MatMult(A, Bx, v2); /* v2 = A*B*x */
432:     VecCopy(v2, v1);
433:     if (rart) {
434:       MatMult(B, v1, v3); /* v3 = R*A*R^t*x */
435:     } else {
436:       MatMultTranspose(B, v1, v3); /* v3 = Bt*A*B*x */
437:     }
438:     VecNorm(v4, NORM_2, &norm_abs);
439:     VecAXPY(v4, none, v3);
440:     VecNorm(v4, NORM_2, &norm_rel);

442:     if (norm_abs > tol) norm_rel /= norm_abs;
443:     if (norm_rel > tol) {
444:       *flg = PETSC_FALSE;
445:       PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel);
446:       break;
447:     }
448:   }

450:   PetscRandomDestroy(&rdm);
451:   VecDestroy(&x);
452:   VecDestroy(&Bx);
453:   VecDestroy(&Cx);
454:   VecDestroy(&v1);
455:   VecDestroy(&v2);
456:   VecDestroy(&v3);
457:   VecDestroy(&v4);
458:   return 0;
459: }

461: /*@
462:    MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B

464:    Collective

466:    Input Parameters:
467: +  A - the first matrix
468: .  B - the second matrix
469: .  C - the third matrix
470: -  n - number of random vectors to be tested

472:    Output Parameter:
473: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

475:    Level: intermediate

477: @*/
478: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
479: {
480:   MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg);
481:   return 0;
482: }

484: /*@
485:    MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t

487:    Collective

489:    Input Parameters:
490: +  A - the first matrix
491: .  B - the second matrix
492: .  C - the third matrix
493: -  n - number of random vectors to be tested

495:    Output Parameter:
496: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

498:    Level: intermediate

500: @*/
501: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
502: {
503:   MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg);
504:   return 0;
505: }

507: /*@
508:    MatIsLinear - Check if a shell matrix A is a linear operator.

510:    Collective

512:    Input Parameters:
513: +  A - the shell matrix
514: -  n - number of random vectors to be tested

516:    Output Parameter:
517: .  flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.

519:    Level: intermediate
520: @*/
521: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
522: {
523:   Vec         x, y, s1, s2;
524:   PetscRandom rctx;
525:   PetscScalar a;
526:   PetscInt    k;
527:   PetscReal   norm, normA;
528:   MPI_Comm    comm;
529:   PetscMPIInt rank;

532:   PetscObjectGetComm((PetscObject)A, &comm);
533:   MPI_Comm_rank(comm, &rank);

535:   PetscRandomCreate(comm, &rctx);
536:   PetscRandomSetFromOptions(rctx);
537:   MatCreateVecs(A, &x, &s1);
538:   VecDuplicate(x, &y);
539:   VecDuplicate(s1, &s2);

541:   *flg = PETSC_TRUE;
542:   for (k = 0; k < n; k++) {
543:     VecSetRandom(x, rctx);
544:     VecSetRandom(y, rctx);
545:     if (rank == 0) PetscRandomGetValue(rctx, &a);
546:     MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);

548:     /* s2 = a*A*x + A*y */
549:     MatMult(A, y, s2);  /* s2 = A*y */
550:     MatMult(A, x, s1);  /* s1 = A*x */
551:     VecAXPY(s2, a, s1); /* s2 = a s1 + s2 */

553:     /* s1 = A * (a x + y) */
554:     VecAXPY(y, a, x); /* y = a x + y */
555:     MatMult(A, y, s1);
556:     VecNorm(s1, NORM_INFINITY, &normA);

558:     VecAXPY(s2, -1.0, s1); /* s2 = - s1 + s2 */
559:     VecNorm(s2, NORM_INFINITY, &norm);
560:     if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
561:       *flg = PETSC_FALSE;
562:       PetscInfo(A, "Error: %" PetscInt_FMT "-th |A*(ax+y) - (a*A*x+A*y)|/|A(ax+y)| %g > tol %g\n", k, (double)(norm / normA), (double)(100. * PETSC_MACHINE_EPSILON));
563:       break;
564:     }
565:   }
566:   PetscRandomDestroy(&rctx);
567:   VecDestroy(&x);
568:   VecDestroy(&y);
569:   VecDestroy(&s1);
570:   VecDestroy(&s2);
571:   return 0;
572: }