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: #if defined(PETSC_USE_INFO)
 12:   const char *sops[] = {"MatMult", "MatMultAdd", "MatMultTranspose", "MatMultTransposeAdd", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd"};
 13:   const char *sop;
 14: #endif

 16:   PetscFunctionBegin;
 19:   PetscCheckSameComm(A, 1, B, 2);
 24:   PetscCall(MatGetLocalSize(A, &am, &an));
 25:   PetscCall(MatGetLocalSize(B, &bm, &bn));
 26:   PetscCheck(am == bm && an == bn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat A,Mat B: local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, am, bm, an, bn);
 27: #if defined(PETSC_USE_INFO)
 28:   sop = sops[(add ? 1 : 0) + 2 * t]; /* t = 0 => no transpose, t = 1 => transpose, t = 2 => Hermitian transpose */
 29: #endif
 30:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
 31:   PetscCall(PetscRandomSetFromOptions(rctx));
 32:   if (t) {
 33:     PetscCall(MatCreateVecs(A, &s1, &Ax));
 34:     PetscCall(MatCreateVecs(B, &s2, &Bx));
 35:   } else {
 36:     PetscCall(MatCreateVecs(A, &Ax, &s1));
 37:     PetscCall(MatCreateVecs(B, &Bx, &s2));
 38:   }
 39:   if (add) {
 40:     PetscCall(VecDuplicate(s1, &Ay));
 41:     PetscCall(VecDuplicate(s2, &By));
 42:   }

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

101: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
102: {
103:   Vec         Ax, Bx, Cx, s1, s2, s3;
104:   PetscRandom rctx;
105:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
106:   PetscInt    am, an, bm, bn, cm, cn, k;
107:   PetscScalar none = -1.0;
108: #if defined(PETSC_USE_INFO)
109:   const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
110:   const char *sop;
111: #endif

113:   PetscFunctionBegin;
116:   PetscCheckSameComm(A, 1, B, 2);
118:   PetscCheckSameComm(A, 1, C, 3);
123:   PetscCall(MatGetLocalSize(A, &am, &an));
124:   PetscCall(MatGetLocalSize(B, &bm, &bn));
125:   PetscCall(MatGetLocalSize(C, &cm, &cn));
126:   if (At) {
127:     PetscInt tt = an;
128:     an          = am;
129:     am          = tt;
130:   };
131:   if (Bt) {
132:     PetscInt tt = bn;
133:     bn          = bm;
134:     bm          = tt;
135:   };
136:   PetscCheck(an == bm && am == cm && bn == cn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat A, B, C local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, am, an, bm, bn, cm, cn);

138: #if defined(PETSC_USE_INFO)
139:   sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
140: #endif
141:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
142:   PetscCall(PetscRandomSetFromOptions(rctx));
143:   if (Bt) {
144:     PetscCall(MatCreateVecs(B, &s1, &Bx));
145:   } else {
146:     PetscCall(MatCreateVecs(B, &Bx, &s1));
147:   }
148:   if (At) {
149:     PetscCall(MatCreateVecs(A, &s2, &Ax));
150:   } else {
151:     PetscCall(MatCreateVecs(A, &Ax, &s2));
152:   }
153:   PetscCall(MatCreateVecs(C, &Cx, &s3));

155:   *flg = PETSC_TRUE;
156:   for (k = 0; k < n; k++) {
157:     PetscCall(VecSetRandom(Bx, rctx));
158:     if (Bt) {
159:       PetscCall(MatMultTranspose(B, Bx, s1));
160:     } else {
161:       PetscCall(MatMult(B, Bx, s1));
162:     }
163:     PetscCall(VecCopy(s1, Ax));
164:     if (At) {
165:       PetscCall(MatMultTranspose(A, Ax, s2));
166:     } else {
167:       PetscCall(MatMult(A, Ax, s2));
168:     }
169:     PetscCall(VecCopy(Bx, Cx));
170:     PetscCall(MatMult(C, Cx, s3));

172:     PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
173:     if (r2 < tol) {
174:       PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
175:     } else {
176:       PetscCall(VecAXPY(s2, none, s3));
177:       PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
178:       r1 /= r2;
179:     }
180:     if (r1 > tol) {
181:       *flg = PETSC_FALSE;
182:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
183:       break;
184:     }
185:   }
186:   PetscCall(PetscRandomDestroy(&rctx));
187:   PetscCall(VecDestroy(&Ax));
188:   PetscCall(VecDestroy(&Bx));
189:   PetscCall(VecDestroy(&Cx));
190:   PetscCall(VecDestroy(&s1));
191:   PetscCall(VecDestroy(&s2));
192:   PetscCall(VecDestroy(&s3));
193:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

196: /*@
197:    MatMultEqual - Compares matrix-vector products of two matrices.

199:    Collective

201:    Input Parameters:
202: +  A - the first matrix
203: .  B - the second matrix
204: -  n - number of random vectors to be tested

206:    Output Parameter:
207: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

209:    Level: intermediate

211: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`
212: @*/
213: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
214: {
215:   PetscFunctionBegin;
216:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, PETSC_FALSE));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: /*@
221:    MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.

223:    Collective

225:    Input Parameters:
226: +  A - the first matrix
227: .  B - the second matrix
228: -  n - number of random vectors to be tested

230:    Output Parameter:
231: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

233:    Level: intermediate

235: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
236: @*/
237: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
238: {
239:   PetscFunctionBegin;
240:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, PETSC_TRUE));
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: /*@
245:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

247:    Collective

249:    Input Parameters:
250: +  A - the first matrix
251: .  B - the second matrix
252: -  n - number of random vectors to be tested

254:    Output Parameter:
255: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

257:    Level: intermediate

259: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
260: @*/
261: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
262: {
263:   PetscFunctionBegin;
264:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, PETSC_FALSE));
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: /*@
269:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

271:    Collective

273:    Input Parameters:
274: +  A - the first matrix
275: .  B - the second matrix
276: -  n - number of random vectors to be tested

278:    Output Parameter:
279: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

281:    Level: intermediate

283: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
284: @*/
285: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
286: {
287:   PetscFunctionBegin;
288:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, PETSC_TRUE));
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: /*@
293:    MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.

295:    Collective

297:    Input Parameters:
298: +  A - the first matrix
299: .  B - the second matrix
300: -  n - number of random vectors to be tested

302:    Output Parameter:
303: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

305:    Level: intermediate

307: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
308: @*/
309: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
310: {
311:   PetscFunctionBegin;
312:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, PETSC_FALSE));
313:   PetscFunctionReturn(PETSC_SUCCESS);
314: }

316: /*@
317:    MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.

319:    Collective

321:    Input Parameters:
322: +  A - the first matrix
323: .  B - the second matrix
324: -  n - number of random vectors to be tested

326:    Output Parameter:
327: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

329:    Level: intermediate

331: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
332: @*/
333: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
334: {
335:   PetscFunctionBegin;
336:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, PETSC_TRUE));
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

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

343:    Collective

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

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

354:    Level: intermediate

356: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
357: @*/
358: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
359: {
360:   PetscFunctionBegin;
361:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }

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

368:    Collective

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

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

379:    Level: intermediate

381: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
382: @*/
383: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
384: {
385:   PetscFunctionBegin;
386:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

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

393:    Collective

395:    Input Parameters:
396: +  A - the first matrix
397: .  B - the second matrix
398: .  C - the third matrix
399: -  n - number of random vectors to be tested

401:    Output Parameter:
402: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

404:    Level: intermediate

406: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
407: @*/
408: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
409: {
410:   PetscFunctionBegin;
411:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
416: {
417:   Vec         x, v1, v2, v3, v4, Cx, Bx;
418:   PetscReal   norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
419:   PetscInt    i, am, an, bm, bn, cm, cn;
420:   PetscRandom rdm;
421:   PetscScalar none = -1.0;

423:   PetscFunctionBegin;
424:   PetscCall(MatGetLocalSize(A, &am, &an));
425:   PetscCall(MatGetLocalSize(B, &bm, &bn));
426:   if (rart) {
427:     PetscInt t = bm;
428:     bm         = bn;
429:     bn         = t;
430:   }
431:   PetscCall(MatGetLocalSize(C, &cm, &cn));
432:   PetscCheck(an == bm && bn == cm && bn == cn, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat A, B, C local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, am, an, bm, bn, cm, cn);

434:   /* Create left vector of A: v2 */
435:   PetscCall(MatCreateVecs(A, &Bx, &v2));

437:   /* Create right vectors of B: x, v3, v4 */
438:   if (rart) {
439:     PetscCall(MatCreateVecs(B, &v1, &x));
440:   } else {
441:     PetscCall(MatCreateVecs(B, &x, &v1));
442:   }
443:   PetscCall(VecDuplicate(x, &v3));

445:   PetscCall(MatCreateVecs(C, &Cx, &v4));
446:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
447:   PetscCall(PetscRandomSetFromOptions(rdm));

449:   *flg = PETSC_TRUE;
450:   for (i = 0; i < n; i++) {
451:     PetscCall(VecSetRandom(x, rdm));
452:     PetscCall(VecCopy(x, Cx));
453:     PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x   */
454:     if (rart) {
455:       PetscCall(MatMultTranspose(B, x, v1));
456:     } else {
457:       PetscCall(MatMult(B, x, v1));
458:     }
459:     PetscCall(VecCopy(v1, Bx));
460:     PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
461:     PetscCall(VecCopy(v2, v1));
462:     if (rart) {
463:       PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
464:     } else {
465:       PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
466:     }
467:     PetscCall(VecNorm(v4, NORM_2, &norm_abs));
468:     PetscCall(VecAXPY(v4, none, v3));
469:     PetscCall(VecNorm(v4, NORM_2, &norm_rel));

471:     if (norm_abs > tol) norm_rel /= norm_abs;
472:     if (norm_rel > tol) {
473:       *flg = PETSC_FALSE;
474:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
475:       break;
476:     }
477:   }

479:   PetscCall(PetscRandomDestroy(&rdm));
480:   PetscCall(VecDestroy(&x));
481:   PetscCall(VecDestroy(&Bx));
482:   PetscCall(VecDestroy(&Cx));
483:   PetscCall(VecDestroy(&v1));
484:   PetscCall(VecDestroy(&v2));
485:   PetscCall(VecDestroy(&v3));
486:   PetscCall(VecDestroy(&v4));
487:   PetscFunctionReturn(PETSC_SUCCESS);
488: }

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

493:    Collective

495:    Input Parameters:
496: +  A - the first matrix
497: .  B - the second matrix
498: .  C - the third matrix
499: -  n - number of random vectors to be tested

501:    Output Parameter:
502: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

504:    Level: intermediate

506: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
507: @*/
508: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
509: {
510:   PetscFunctionBegin;
511:   PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

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

518:    Collective

520:    Input Parameters:
521: +  A - the first matrix
522: .  B - the second matrix
523: .  C - the third matrix
524: -  n - number of random vectors to be tested

526:    Output Parameter:
527: .  flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

529:    Level: intermediate

531: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
532: @*/
533: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
534: {
535:   PetscFunctionBegin;
536:   PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
537:   PetscFunctionReturn(PETSC_SUCCESS);
538: }

540: /*@
541:    MatIsLinear - Check if a shell matrix `A` is a linear operator.

543:    Collective

545:    Input Parameters:
546: +  A - the shell matrix
547: -  n - number of random vectors to be tested

549:    Output Parameter:
550: .  flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.

552:    Level: intermediate

554: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
555: @*/
556: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
557: {
558:   Vec         x, y, s1, s2;
559:   PetscRandom rctx;
560:   PetscScalar a;
561:   PetscInt    k;
562:   PetscReal   norm, normA;
563:   MPI_Comm    comm;
564:   PetscMPIInt rank;

566:   PetscFunctionBegin;
568:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
569:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

571:   PetscCall(PetscRandomCreate(comm, &rctx));
572:   PetscCall(PetscRandomSetFromOptions(rctx));
573:   PetscCall(MatCreateVecs(A, &x, &s1));
574:   PetscCall(VecDuplicate(x, &y));
575:   PetscCall(VecDuplicate(s1, &s2));

577:   *flg = PETSC_TRUE;
578:   for (k = 0; k < n; k++) {
579:     PetscCall(VecSetRandom(x, rctx));
580:     PetscCall(VecSetRandom(y, rctx));
581:     if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
582:     PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));

584:     /* s2 = a*A*x + A*y */
585:     PetscCall(MatMult(A, y, s2));  /* s2 = A*y */
586:     PetscCall(MatMult(A, x, s1));  /* s1 = A*x */
587:     PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */

589:     /* s1 = A * (a x + y) */
590:     PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
591:     PetscCall(MatMult(A, y, s1));
592:     PetscCall(VecNorm(s1, NORM_INFINITY, &normA));

594:     PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
595:     PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
596:     if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
597:       *flg = PETSC_FALSE;
598:       PetscCall(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)));
599:       break;
600:     }
601:   }
602:   PetscCall(PetscRandomDestroy(&rctx));
603:   PetscCall(VecDestroy(&x));
604:   PetscCall(VecDestroy(&y));
605:   PetscCall(VecDestroy(&s1));
606:   PetscCall(VecDestroy(&s2));
607:   PetscFunctionReturn(PETSC_SUCCESS);
608: }