Actual source code: multequal.c
1: #include <petsc/private/matimpl.h>
3: static PetscErrorCode MatMultEqual_Private(Mat A, Mat B, PetscInt n, PetscBool *flg, PetscInt t, PetscInt add)
4: {
5: Vec Ax = NULL, Bx = NULL, s1 = NULL, s2 = NULL, Ay = NULL, By = NULL;
6: PetscRandom rctx;
7: PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
8: PetscInt am, an, bm, bn, k;
9: PetscScalar none = -1.0;
10: #if defined(PETSC_USE_INFO)
11: const char *sops[] = {"MatMult", "MatMultAdd", "MatMultAdd (update)", "MatMultTranspose", "MatMultTransposeAdd", "MatMultTransposeAdd (update)", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd", "MatMultHermitianTransposeAdd (update)"};
12: const char *sop;
13: #endif
15: PetscFunctionBegin;
18: PetscCheckSameComm(A, 1, B, 2);
20: PetscAssertPointer(flg, 4);
23: PetscCall(MatGetLocalSize(A, &am, &an));
24: PetscCall(MatGetLocalSize(B, &bm, &bn));
25: 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);
26: #if defined(PETSC_USE_INFO)
27: sop = sops[add + 3 * t]; /* add = 0 => no add, add = 1 => add third vector, add = 2 => add update, t = 0 => no transpose, t = 1 => transpose, t = 2 => Hermitian transpose */
28: #endif
29: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
30: PetscCall(PetscRandomSetFromOptions(rctx));
31: if (t) {
32: PetscCall(MatCreateVecs(A, &s1, &Ax));
33: PetscCall(MatCreateVecs(B, &s2, &Bx));
34: } else {
35: PetscCall(MatCreateVecs(A, &Ax, &s1));
36: PetscCall(MatCreateVecs(B, &Bx, &s2));
37: }
38: if (add) {
39: PetscCall(VecDuplicate(s1, &Ay));
40: PetscCall(VecDuplicate(s2, &By));
41: }
43: *flg = PETSC_TRUE;
44: for (k = 0; k < n; k++) {
45: Vec Aadd = NULL, Badd = NULL;
47: PetscCall(VecSetRandom(Ax, rctx));
48: PetscCall(VecCopy(Ax, Bx));
49: if (add) {
50: PetscCall(VecSetRandom(Ay, rctx));
51: PetscCall(VecCopy(Ay, By));
52: Aadd = Ay;
53: Badd = By;
54: if (add == 2) {
55: PetscCall(VecCopy(Ay, s1));
56: PetscCall(VecCopy(By, s2));
57: Aadd = s1;
58: Badd = s2;
59: }
60: }
61: if (t == 1) {
62: if (add) {
63: PetscCall(MatMultTransposeAdd(A, Ax, Aadd, s1));
64: PetscCall(MatMultTransposeAdd(B, Bx, Badd, s2));
65: } else {
66: PetscCall(MatMultTranspose(A, Ax, s1));
67: PetscCall(MatMultTranspose(B, Bx, s2));
68: }
69: } else if (t == 2) {
70: if (add) {
71: PetscCall(MatMultHermitianTransposeAdd(A, Ax, Aadd, s1));
72: PetscCall(MatMultHermitianTransposeAdd(B, Bx, Badd, s2));
73: } else {
74: PetscCall(MatMultHermitianTranspose(A, Ax, s1));
75: PetscCall(MatMultHermitianTranspose(B, Bx, s2));
76: }
77: } else {
78: if (add) {
79: PetscCall(MatMultAdd(A, Ax, Aadd, s1));
80: PetscCall(MatMultAdd(B, Bx, Badd, s2));
81: } else {
82: PetscCall(MatMult(A, Ax, s1));
83: PetscCall(MatMult(B, Bx, s2));
84: }
85: }
86: PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
87: if (r2 < tol) {
88: PetscCall(VecNorm(s1, NORM_INFINITY, &r1));
89: } else {
90: PetscCall(VecAXPY(s2, none, s1));
91: PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
92: r1 /= r2;
93: }
94: if (r1 > tol) {
95: *flg = PETSC_FALSE;
96: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1));
97: break;
98: }
99: }
100: PetscCall(PetscRandomDestroy(&rctx));
101: PetscCall(VecDestroy(&Ax));
102: PetscCall(VecDestroy(&Bx));
103: PetscCall(VecDestroy(&Ay));
104: PetscCall(VecDestroy(&By));
105: PetscCall(VecDestroy(&s1));
106: PetscCall(VecDestroy(&s2));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
111: {
112: Vec Ax, Bx, Cx, s1, s2, s3;
113: PetscRandom rctx;
114: PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
115: PetscInt am, an, bm, bn, cm, cn, k;
116: PetscScalar none = -1.0;
117: #if defined(PETSC_USE_INFO)
118: const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
119: const char *sop;
120: #endif
122: PetscFunctionBegin;
125: PetscCheckSameComm(A, 1, B, 2);
127: PetscCheckSameComm(A, 1, C, 3);
129: PetscAssertPointer(flg, 5);
132: PetscCall(MatGetLocalSize(A, &am, &an));
133: PetscCall(MatGetLocalSize(B, &bm, &bn));
134: PetscCall(MatGetLocalSize(C, &cm, &cn));
135: if (At) {
136: PetscInt tt = an;
137: an = am;
138: am = tt;
139: }
140: if (Bt) {
141: PetscInt tt = bn;
142: bn = bm;
143: bm = tt;
144: }
145: 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);
147: #if defined(PETSC_USE_INFO)
148: sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
149: #endif
150: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
151: PetscCall(PetscRandomSetFromOptions(rctx));
152: if (Bt) {
153: PetscCall(MatCreateVecs(B, &s1, &Bx));
154: } else {
155: PetscCall(MatCreateVecs(B, &Bx, &s1));
156: }
157: if (At) {
158: PetscCall(MatCreateVecs(A, &s2, &Ax));
159: } else {
160: PetscCall(MatCreateVecs(A, &Ax, &s2));
161: }
162: PetscCall(MatCreateVecs(C, &Cx, &s3));
164: *flg = PETSC_TRUE;
165: for (k = 0; k < n; k++) {
166: PetscCall(VecSetRandom(Bx, rctx));
167: if (Bt) {
168: PetscCall(MatMultTranspose(B, Bx, s1));
169: } else {
170: PetscCall(MatMult(B, Bx, s1));
171: }
172: PetscCall(VecCopy(s1, Ax));
173: if (At) {
174: PetscCall(MatMultTranspose(A, Ax, s2));
175: } else {
176: PetscCall(MatMult(A, Ax, s2));
177: }
178: PetscCall(VecCopy(Bx, Cx));
179: PetscCall(MatMult(C, Cx, s3));
181: PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
182: if (r2 < tol) {
183: PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
184: } else {
185: PetscCall(VecAXPY(s2, none, s3));
186: PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
187: r1 /= r2;
188: }
189: if (r1 > tol) {
190: *flg = PETSC_FALSE;
191: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
192: break;
193: }
194: }
195: PetscCall(PetscRandomDestroy(&rctx));
196: PetscCall(VecDestroy(&Ax));
197: PetscCall(VecDestroy(&Bx));
198: PetscCall(VecDestroy(&Cx));
199: PetscCall(VecDestroy(&s1));
200: PetscCall(VecDestroy(&s2));
201: PetscCall(VecDestroy(&s3));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /*@
206: MatMultEqual - Compares matrix-vector products of two matrices.
208: Collective
210: Input Parameters:
211: + A - the first matrix
212: . B - the second matrix
213: - n - number of random vectors to be tested
215: Output Parameter:
216: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
218: Level: intermediate
220: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`
221: @*/
222: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
223: {
224: PetscFunctionBegin;
225: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 0));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@
230: MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.
232: Collective
234: Input Parameters:
235: + A - the first matrix
236: . B - the second matrix
237: - n - number of random vectors to be tested
239: Output Parameter:
240: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
242: Level: intermediate
244: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
245: @*/
246: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
247: {
248: PetscFunctionBegin;
249: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 1));
250: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 2));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*@
255: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
257: Collective
259: Input Parameters:
260: + A - the first matrix
261: . B - the second matrix
262: - n - number of random vectors to be tested
264: Output Parameter:
265: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
267: Level: intermediate
269: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
270: @*/
271: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
272: {
273: PetscFunctionBegin;
274: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 0));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: /*@
279: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
281: Collective
283: Input Parameters:
284: + A - the first matrix
285: . B - the second matrix
286: - n - number of random vectors to be tested
288: Output Parameter:
289: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
291: Level: intermediate
293: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
294: @*/
295: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
296: {
297: PetscFunctionBegin;
298: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 1));
299: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 2));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /*@
304: MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
306: Collective
308: Input Parameters:
309: + A - the first matrix
310: . B - the second matrix
311: - n - number of random vectors to be tested
313: Output Parameter:
314: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
316: Level: intermediate
318: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
319: @*/
320: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
321: {
322: PetscFunctionBegin;
323: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 0));
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
330: Collective
332: Input Parameters:
333: + A - the first matrix
334: . B - the second matrix
335: - n - number of random vectors to be tested
337: Output Parameter:
338: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
340: Level: intermediate
342: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
343: @*/
344: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
345: {
346: PetscFunctionBegin;
347: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 1));
348: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 2));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@
353: MatMatMultEqual - Test A*B*x = C*x for n random vector x
355: Collective
357: Input Parameters:
358: + A - the first matrix
359: . B - the second matrix
360: . C - the third matrix
361: - n - number of random vectors to be tested
363: Output Parameter:
364: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
366: Level: intermediate
368: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
369: @*/
370: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
371: {
372: PetscFunctionBegin;
373: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: /*@
378: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
380: Collective
382: Input Parameters:
383: + A - the first matrix
384: . B - the second matrix
385: . C - the third matrix
386: - n - number of random vectors to be tested
388: Output Parameter:
389: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
391: Level: intermediate
393: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
394: @*/
395: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
396: {
397: PetscFunctionBegin;
398: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*@
403: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
405: Collective
407: Input Parameters:
408: + A - the first matrix
409: . B - the second matrix
410: . C - the third matrix
411: - n - number of random vectors to be tested
413: Output Parameter:
414: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
416: Level: intermediate
418: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
419: @*/
420: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
421: {
422: PetscFunctionBegin;
423: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
428: {
429: Vec x, v1, v2, v3, v4, Cx, Bx;
430: PetscReal norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
431: PetscInt i, am, an, bm, bn, cm, cn;
432: PetscRandom rdm;
433: PetscScalar none = -1.0;
435: PetscFunctionBegin;
436: PetscCall(MatGetLocalSize(A, &am, &an));
437: PetscCall(MatGetLocalSize(B, &bm, &bn));
438: if (rart) {
439: PetscInt t = bm;
440: bm = bn;
441: bn = t;
442: }
443: PetscCall(MatGetLocalSize(C, &cm, &cn));
444: 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);
446: /* Create left vector of A: v2 */
447: PetscCall(MatCreateVecs(A, &Bx, &v2));
449: /* Create right vectors of B: x, v3, v4 */
450: if (rart) {
451: PetscCall(MatCreateVecs(B, &v1, &x));
452: } else {
453: PetscCall(MatCreateVecs(B, &x, &v1));
454: }
455: PetscCall(VecDuplicate(x, &v3));
457: PetscCall(MatCreateVecs(C, &Cx, &v4));
458: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
459: PetscCall(PetscRandomSetFromOptions(rdm));
461: *flg = PETSC_TRUE;
462: for (i = 0; i < n; i++) {
463: PetscCall(VecSetRandom(x, rdm));
464: PetscCall(VecCopy(x, Cx));
465: PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x */
466: if (rart) {
467: PetscCall(MatMultTranspose(B, x, v1));
468: } else {
469: PetscCall(MatMult(B, x, v1));
470: }
471: PetscCall(VecCopy(v1, Bx));
472: PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
473: PetscCall(VecCopy(v2, v1));
474: if (rart) {
475: PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
476: } else {
477: PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
478: }
479: PetscCall(VecNorm(v4, NORM_2, &norm_abs));
480: PetscCall(VecAXPY(v4, none, v3));
481: PetscCall(VecNorm(v4, NORM_2, &norm_rel));
483: if (norm_abs > tol) norm_rel /= norm_abs;
484: if (norm_rel > tol) {
485: *flg = PETSC_FALSE;
486: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
487: break;
488: }
489: }
491: PetscCall(PetscRandomDestroy(&rdm));
492: PetscCall(VecDestroy(&x));
493: PetscCall(VecDestroy(&Bx));
494: PetscCall(VecDestroy(&Cx));
495: PetscCall(VecDestroy(&v1));
496: PetscCall(VecDestroy(&v2));
497: PetscCall(VecDestroy(&v3));
498: PetscCall(VecDestroy(&v4));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /*@
503: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
505: Collective
507: Input Parameters:
508: + A - the first matrix
509: . B - the second matrix
510: . C - the third matrix
511: - n - number of random vectors to be tested
513: Output Parameter:
514: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
516: Level: intermediate
518: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
519: @*/
520: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
521: {
522: PetscFunctionBegin;
523: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
530: Collective
532: Input Parameters:
533: + A - the first matrix
534: . B - the second matrix
535: . C - the third matrix
536: - n - number of random vectors to be tested
538: Output Parameter:
539: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
541: Level: intermediate
543: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
544: @*/
545: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
546: {
547: PetscFunctionBegin;
548: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: /*@
553: MatIsLinear - Check if a shell matrix `A` is a linear operator.
555: Collective
557: Input Parameters:
558: + A - the shell matrix
559: - n - number of random vectors to be tested
561: Output Parameter:
562: . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.
564: Level: intermediate
566: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
567: @*/
568: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
569: {
570: Vec x, y, s1, s2;
571: PetscRandom rctx;
572: PetscScalar a;
573: PetscInt k;
574: PetscReal norm, normA;
575: MPI_Comm comm;
576: PetscMPIInt rank;
578: PetscFunctionBegin;
580: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
581: PetscCallMPI(MPI_Comm_rank(comm, &rank));
583: PetscCall(PetscRandomCreate(comm, &rctx));
584: PetscCall(PetscRandomSetFromOptions(rctx));
585: PetscCall(MatCreateVecs(A, &x, &s1));
586: PetscCall(VecDuplicate(x, &y));
587: PetscCall(VecDuplicate(s1, &s2));
589: *flg = PETSC_TRUE;
590: for (k = 0; k < n; k++) {
591: PetscCall(VecSetRandom(x, rctx));
592: PetscCall(VecSetRandom(y, rctx));
593: if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
594: PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));
596: /* s2 = a*A*x + A*y */
597: PetscCall(MatMult(A, y, s2)); /* s2 = A*y */
598: PetscCall(MatMult(A, x, s1)); /* s1 = A*x */
599: PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */
601: /* s1 = A * (a x + y) */
602: PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
603: PetscCall(MatMult(A, y, s1));
604: PetscCall(VecNorm(s1, NORM_INFINITY, &normA));
606: PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
607: PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
608: if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
609: *flg = PETSC_FALSE;
610: 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)));
611: break;
612: }
613: }
614: PetscCall(PetscRandomDestroy(&rctx));
615: PetscCall(VecDestroy(&x));
616: PetscCall(VecDestroy(&y));
617: PetscCall(VecDestroy(&s1));
618: PetscCall(VecDestroy(&s2));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }