Actual source code: multequal.c

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

  3: /*
  4:   n; try the MatMult variant n times
  5:   flg: return the boolean result, equal or not
  6:   t: 0 => no transpose; 1 => transpose; 2 => Hermitian transpose
  7:   add:  0 => no add (e.g., y = Ax);  1 => add third vector (e.g., z = Ax + y); 2 => add update (e.g., y = Ax + y)
  8: */
  9: static PetscErrorCode MatMultEqual_Private(Mat A, Mat B, PetscInt n, PetscBool *flg, PetscInt t, PetscInt add)
 10: {
 11:   Vec         Ax = NULL, Bx = NULL, s1 = NULL, s2 = NULL, Ay = NULL, By = NULL;
 12:   PetscRandom rctx;
 13:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
 14:   PetscInt    am, an, bm, bn, k;
 15:   PetscScalar none = -1.0;
 16: #if defined(PETSC_USE_INFO)
 17:   const char *sops[] = {"MatMult", "MatMultAdd", "MatMultAdd (update)", "MatMultTranspose", "MatMultTransposeAdd", "MatMultTransposeAdd (update)", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd", "MatMultHermitianTransposeAdd (update)"};
 18:   const char *sop;
 19: #endif

 21:   PetscFunctionBegin;
 24:   PetscCheckSameComm(A, 1, B, 2);
 26:   PetscAssertPointer(flg, 4);
 29:   PetscCall(MatGetLocalSize(A, &am, &an));
 30:   PetscCall(MatGetLocalSize(B, &bm, &bn));
 31:   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);
 32: #if defined(PETSC_USE_INFO)
 33:   sop = sops[add + 3 * t];
 34: #endif
 35:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
 36:   PetscCall(PetscRandomSetFromOptions(rctx));
 37:   if (t) {
 38:     PetscCall(MatCreateVecs(A, &s1, &Ax));
 39:     PetscCall(MatCreateVecs(B, &s2, &Bx));
 40:   } else {
 41:     PetscCall(MatCreateVecs(A, &Ax, &s1));
 42:     PetscCall(MatCreateVecs(B, &Bx, &s2));
 43:   }
 44:   if (add) {
 45:     PetscCall(VecDuplicate(s1, &Ay));
 46:     PetscCall(VecDuplicate(s2, &By));
 47:   }

 49:   *flg = PETSC_TRUE;
 50:   for (k = 0; k < n; k++) {
 51:     Vec Aadd = NULL, Badd = NULL;

 53:     PetscCall(VecSetRandom(Ax, rctx));
 54:     PetscCall(VecCopy(Ax, Bx));
 55:     if (add) {
 56:       PetscCall(VecSetRandom(Ay, rctx));
 57:       PetscCall(VecCopy(Ay, By));
 58:       Aadd = Ay;
 59:       Badd = By;
 60:       if (add == 2) {
 61:         PetscCall(VecCopy(Ay, s1));
 62:         PetscCall(VecCopy(By, s2));
 63:         Aadd = s1;
 64:         Badd = s2;
 65:       }
 66:     }
 67:     if (t == 1) {
 68:       if (add) {
 69:         PetscCall(MatMultTransposeAdd(A, Ax, Aadd, s1));
 70:         PetscCall(MatMultTransposeAdd(B, Bx, Badd, s2));
 71:       } else {
 72:         PetscCall(MatMultTranspose(A, Ax, s1));
 73:         PetscCall(MatMultTranspose(B, Bx, s2));
 74:       }
 75:     } else if (t == 2) {
 76:       if (add) {
 77:         PetscCall(MatMultHermitianTransposeAdd(A, Ax, Aadd, s1));
 78:         PetscCall(MatMultHermitianTransposeAdd(B, Bx, Badd, s2));
 79:       } else {
 80:         PetscCall(MatMultHermitianTranspose(A, Ax, s1));
 81:         PetscCall(MatMultHermitianTranspose(B, Bx, s2));
 82:       }
 83:     } else {
 84:       if (add) {
 85:         PetscCall(MatMultAdd(A, Ax, Aadd, s1));
 86:         PetscCall(MatMultAdd(B, Bx, Badd, s2));
 87:       } else {
 88:         PetscCall(MatMult(A, Ax, s1));
 89:         PetscCall(MatMult(B, Bx, s2));
 90:       }
 91:     }
 92:     PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
 93:     if (r2 < tol) {
 94:       PetscCall(VecNorm(s1, NORM_INFINITY, &r1));
 95:     } else {
 96:       PetscCall(VecAXPY(s2, none, s1));
 97:       PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
 98:       r1 /= r2;
 99:     }
100:     if (r1 > tol) {
101:       *flg = PETSC_FALSE;
102:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1));
103:       break;
104:     }
105:   }
106:   PetscCall(PetscRandomDestroy(&rctx));
107:   PetscCall(VecDestroy(&Ax));
108:   PetscCall(VecDestroy(&Bx));
109:   PetscCall(VecDestroy(&Ay));
110:   PetscCall(VecDestroy(&By));
111:   PetscCall(VecDestroy(&s1));
112:   PetscCall(VecDestroy(&s2));
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
117: {
118:   Vec         Ax, Bx, Cx, s1, s2, s3;
119:   PetscRandom rctx;
120:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
121:   PetscInt    am, an, bm, bn, cm, cn, k;
122:   PetscScalar none = -1.0;
123: #if defined(PETSC_USE_INFO)
124:   const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
125:   const char *sop;
126: #endif

128:   PetscFunctionBegin;
131:   PetscCheckSameComm(A, 1, B, 2);
133:   PetscCheckSameComm(A, 1, C, 3);
135:   PetscAssertPointer(flg, 5);
138:   PetscCall(MatGetLocalSize(A, &am, &an));
139:   PetscCall(MatGetLocalSize(B, &bm, &bn));
140:   PetscCall(MatGetLocalSize(C, &cm, &cn));
141:   if (At) {
142:     PetscInt tt = an;
143:     an          = am;
144:     am          = tt;
145:   }
146:   if (Bt) {
147:     PetscInt tt = bn;
148:     bn          = bm;
149:     bm          = tt;
150:   }
151:   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);

153: #if defined(PETSC_USE_INFO)
154:   sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
155: #endif
156:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
157:   PetscCall(PetscRandomSetFromOptions(rctx));
158:   if (Bt) {
159:     PetscCall(MatCreateVecs(B, &s1, &Bx));
160:   } else {
161:     PetscCall(MatCreateVecs(B, &Bx, &s1));
162:   }
163:   if (At) {
164:     PetscCall(MatCreateVecs(A, &s2, &Ax));
165:   } else {
166:     PetscCall(MatCreateVecs(A, &Ax, &s2));
167:   }
168:   PetscCall(MatCreateVecs(C, &Cx, &s3));

170:   *flg = PETSC_TRUE;
171:   for (k = 0; k < n; k++) {
172:     PetscCall(VecSetRandom(Bx, rctx));
173:     if (Bt) {
174:       PetscCall(MatMultTranspose(B, Bx, s1));
175:     } else {
176:       PetscCall(MatMult(B, Bx, s1));
177:     }
178:     PetscCall(VecCopy(s1, Ax));
179:     if (At) {
180:       PetscCall(MatMultTranspose(A, Ax, s2));
181:     } else {
182:       PetscCall(MatMult(A, Ax, s2));
183:     }
184:     PetscCall(VecCopy(Bx, Cx));
185:     PetscCall(MatMult(C, Cx, s3));

187:     PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
188:     if (r2 < tol) {
189:       PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
190:     } else {
191:       PetscCall(VecAXPY(s2, none, s3));
192:       PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
193:       r1 /= r2;
194:     }
195:     if (r1 > tol) {
196:       *flg = PETSC_FALSE;
197:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
198:       break;
199:     }
200:   }
201:   PetscCall(PetscRandomDestroy(&rctx));
202:   PetscCall(VecDestroy(&Ax));
203:   PetscCall(VecDestroy(&Bx));
204:   PetscCall(VecDestroy(&Cx));
205:   PetscCall(VecDestroy(&s1));
206:   PetscCall(VecDestroy(&s2));
207:   PetscCall(VecDestroy(&s3));
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

211: /*@
212:   MatMultEqual - Compares matrix-vector products of two matrices using `n` random vectors

214:   Collective

216:   Input Parameters:
217: + A - the first matrix
218: . B - the second matrix
219: - n - number of random vectors to be tested

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

224:   Level: intermediate

226:   Note:
227:   The tolerance for equality is a generous `PETSC_SQRT_MACHINE_EPSILON` in the norm of the difference of the two computed vectors to
228:   allow for differences in the numerical computations. Hence this routine may indicate equality even if there is a small systematic difference
229:   between the two matrices.

231: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`, `MatEqual()`
232: @*/
233: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
234: {
235:   PetscFunctionBegin;
236:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 0));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

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

243:   Collective

245:   Input Parameters:
246: + A - the first matrix
247: . B - the second matrix
248: - n - number of random vectors to be tested

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

253:   Level: intermediate

255: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
256: @*/
257: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
258: {
259:   PetscFunctionBegin;
260:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 1));
261:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 2));
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

265: /*@
266:   MatMultTransposeEqual - Compares matrix-vector products of two matrices.

268:   Collective

270:   Input Parameters:
271: + A - the first matrix
272: . B - the second matrix
273: - n - number of random vectors to be tested

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

278:   Level: intermediate

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

289: /*@
290:   MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

292:   Collective

294:   Input Parameters:
295: + A - the first matrix
296: . B - the second matrix
297: - n - number of random vectors to be tested

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

302:   Level: intermediate

304: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
305: @*/
306: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
307: {
308:   PetscFunctionBegin;
309:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 1));
310:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 2));
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: /*@
315:   MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.

317:   Collective

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

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

327:   Level: intermediate

329: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
330: @*/
331: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
332: {
333:   PetscFunctionBegin;
334:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 0));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: /*@
339:   MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.

341:   Collective

343:   Input Parameters:
344: + A - the first matrix
345: . B - the second matrix
346: - n - number of random vectors to be tested

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

351:   Level: intermediate

353: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
354: @*/
355: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
356: {
357:   PetscFunctionBegin;
358:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 1));
359:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 2));
360:   PetscFunctionReturn(PETSC_SUCCESS);
361: }

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

366:   Collective

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

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

377:   Level: intermediate

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

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

391:   Collective

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

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

402:   Level: intermediate

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

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

416:   Collective

418:   Input Parameters:
419: + A - the first matrix
420: . B - the second matrix
421: . C - the third matrix
422: - n - number of random vectors to be tested

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

427:   Level: intermediate

429: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
430: @*/
431: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
432: {
433:   PetscFunctionBegin;
434:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
439: {
440:   Vec         x, v1, v2, v3, v4, Cx, Bx;
441:   PetscReal   norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
442:   PetscInt    i, am, an, bm, bn, cm, cn;
443:   PetscRandom rdm;
444:   PetscScalar none = -1.0;

446:   PetscFunctionBegin;
447:   PetscCall(MatGetLocalSize(A, &am, &an));
448:   PetscCall(MatGetLocalSize(B, &bm, &bn));
449:   if (rart) {
450:     PetscInt t = bm;
451:     bm         = bn;
452:     bn         = t;
453:   }
454:   PetscCall(MatGetLocalSize(C, &cm, &cn));
455:   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);

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

460:   /* Create right vectors of B: x, v3, v4 */
461:   if (rart) {
462:     PetscCall(MatCreateVecs(B, &v1, &x));
463:   } else {
464:     PetscCall(MatCreateVecs(B, &x, &v1));
465:   }
466:   PetscCall(VecDuplicate(x, &v3));

468:   PetscCall(MatCreateVecs(C, &Cx, &v4));
469:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
470:   PetscCall(PetscRandomSetFromOptions(rdm));

472:   *flg = PETSC_TRUE;
473:   for (i = 0; i < n; i++) {
474:     PetscCall(VecSetRandom(x, rdm));
475:     PetscCall(VecCopy(x, Cx));
476:     PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x   */
477:     if (rart) {
478:       PetscCall(MatMultTranspose(B, x, v1));
479:     } else {
480:       PetscCall(MatMult(B, x, v1));
481:     }
482:     PetscCall(VecCopy(v1, Bx));
483:     PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
484:     PetscCall(VecCopy(v2, v1));
485:     if (rart) {
486:       PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
487:     } else {
488:       PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
489:     }
490:     PetscCall(VecNorm(v4, NORM_2, &norm_abs));
491:     PetscCall(VecAXPY(v4, none, v3));
492:     PetscCall(VecNorm(v4, NORM_2, &norm_rel));

494:     if (norm_abs > tol) norm_rel /= norm_abs;
495:     if (norm_rel > tol) {
496:       *flg = PETSC_FALSE;
497:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
498:       break;
499:     }
500:   }

502:   PetscCall(PetscRandomDestroy(&rdm));
503:   PetscCall(VecDestroy(&x));
504:   PetscCall(VecDestroy(&Bx));
505:   PetscCall(VecDestroy(&Cx));
506:   PetscCall(VecDestroy(&v1));
507:   PetscCall(VecDestroy(&v2));
508:   PetscCall(VecDestroy(&v3));
509:   PetscCall(VecDestroy(&v4));
510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }

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

516:   Collective

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

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

527:   Level: intermediate

529: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
530: @*/
531: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
532: {
533:   PetscFunctionBegin;
534:   PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }

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

541:   Collective

543:   Input Parameters:
544: + A - the first matrix
545: . B - the second matrix
546: . C - the third matrix
547: - n - number of random vectors to be tested

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

552:   Level: intermediate

554: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
555: @*/
556: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
557: {
558:   PetscFunctionBegin;
559:   PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
560:   PetscFunctionReturn(PETSC_SUCCESS);
561: }

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

566:   Collective

568:   Input Parameters:
569: + A - the shell matrix
570: - n - number of random vectors to be tested

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

575:   Level: intermediate

577: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
578: @*/
579: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
580: {
581:   Vec         x, y, s1, s2;
582:   PetscRandom rctx;
583:   PetscScalar a;
584:   PetscInt    k;
585:   PetscReal   norm, normA;
586:   MPI_Comm    comm;
587:   PetscMPIInt rank;

589:   PetscFunctionBegin;
591:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
592:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

594:   PetscCall(PetscRandomCreate(comm, &rctx));
595:   PetscCall(PetscRandomSetFromOptions(rctx));
596:   PetscCall(MatCreateVecs(A, &x, &s1));
597:   PetscCall(VecDuplicate(x, &y));
598:   PetscCall(VecDuplicate(s1, &s2));

600:   *flg = PETSC_TRUE;
601:   for (k = 0; k < n; k++) {
602:     PetscCall(VecSetRandom(x, rctx));
603:     PetscCall(VecSetRandom(y, rctx));
604:     if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
605:     PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));

607:     /* s2 = a*A*x + A*y */
608:     PetscCall(MatMult(A, y, s2));  /* s2 = A*y */
609:     PetscCall(MatMult(A, x, s1));  /* s1 = A*x */
610:     PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */

612:     /* s1 = A * (a x + y) */
613:     PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
614:     PetscCall(MatMult(A, y, s1));
615:     PetscCall(VecNorm(s1, NORM_INFINITY, &normA));

617:     PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
618:     PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
619:     if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
620:       *flg = PETSC_FALSE;
621:       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)));
622:       break;
623:     }
624:   }
625:   PetscCall(PetscRandomDestroy(&rctx));
626:   PetscCall(VecDestroy(&x));
627:   PetscCall(VecDestroy(&y));
628:   PetscCall(VecDestroy(&s1));
629:   PetscCall(VecDestroy(&s2));
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }