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: #if PetscDefined(USE_INFO)
16: const char *sops[] = {"MatMult", "MatMultAdd", "MatMultAdd (update)", "MatMultTranspose", "MatMultTransposeAdd", "MatMultTransposeAdd (update)", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd", "MatMultHermitianTransposeAdd (update)"};
17: const char *sop;
18: #endif
20: PetscFunctionBegin;
23: PetscCheckSameComm(A, 1, B, 2);
25: PetscAssertPointer(flg, 4);
28: PetscCall(MatGetLocalSize(A, &am, &an));
29: PetscCall(MatGetLocalSize(B, &bm, &bn));
30: 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);
31: #if PetscDefined(USE_INFO)
32: sop = sops[add + 3 * t];
33: #endif
34: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
35: PetscCall(PetscRandomSetFromOptions(rctx));
36: if (t) {
37: PetscCall(MatCreateVecs(A, &s1, &Ax));
38: PetscCall(MatCreateVecs(B, &s2, &Bx));
39: } else {
40: PetscCall(MatCreateVecs(A, &Ax, &s1));
41: PetscCall(MatCreateVecs(B, &Bx, &s2));
42: }
43: if (add) {
44: PetscCall(VecDuplicate(s1, &Ay));
45: PetscCall(VecDuplicate(s2, &By));
46: }
48: *flg = PETSC_TRUE;
49: for (k = 0; k < n; k++) {
50: Vec Aadd = NULL, Badd = NULL;
52: PetscCall(VecSetRandom(Ax, rctx));
53: PetscCall(VecCopy(Ax, Bx));
54: if (add) {
55: PetscCall(VecSetRandom(Ay, rctx));
56: PetscCall(VecCopy(Ay, By));
57: Aadd = Ay;
58: Badd = By;
59: if (add == 2) {
60: PetscCall(VecCopy(Ay, s1));
61: PetscCall(VecCopy(By, s2));
62: Aadd = s1;
63: Badd = s2;
64: }
65: }
66: if (t == 1) {
67: if (add) {
68: PetscCall(MatMultTransposeAdd(A, Ax, Aadd, s1));
69: PetscCall(MatMultTransposeAdd(B, Bx, Badd, s2));
70: } else {
71: PetscCall(MatMultTranspose(A, Ax, s1));
72: PetscCall(MatMultTranspose(B, Bx, s2));
73: }
74: } else if (t == 2) {
75: if (add) {
76: PetscCall(MatMultHermitianTransposeAdd(A, Ax, Aadd, s1));
77: PetscCall(MatMultHermitianTransposeAdd(B, Bx, Badd, s2));
78: } else {
79: PetscCall(MatMultHermitianTranspose(A, Ax, s1));
80: PetscCall(MatMultHermitianTranspose(B, Bx, s2));
81: }
82: } else {
83: if (add) {
84: PetscCall(MatMultAdd(A, Ax, Aadd, s1));
85: PetscCall(MatMultAdd(B, Bx, Badd, s2));
86: } else {
87: PetscCall(MatMult(A, Ax, s1));
88: PetscCall(MatMult(B, Bx, s2));
89: }
90: }
91: PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
92: if (r2 < tol) {
93: PetscCall(VecNorm(s1, NORM_INFINITY, &r1));
94: } else {
95: PetscCall(VecAXPY(s2, -1.0, s1));
96: PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
97: r1 /= r2;
98: }
99: if (r1 > tol) {
100: *flg = PETSC_FALSE;
101: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1));
102: break;
103: }
104: }
105: PetscCall(PetscRandomDestroy(&rctx));
106: PetscCall(VecDestroy(&Ax));
107: PetscCall(VecDestroy(&Bx));
108: PetscCall(VecDestroy(&Ay));
109: PetscCall(VecDestroy(&By));
110: PetscCall(VecDestroy(&s1));
111: PetscCall(VecDestroy(&s2));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
116: {
117: Vec Ax, Bx, Cx, s1, s2, s3;
118: PetscRandom rctx;
119: PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
120: PetscInt am, an, bm, bn, cm, cn, k;
121: #if PetscDefined(USE_INFO)
122: const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
123: const char *sop;
124: #endif
126: PetscFunctionBegin;
129: PetscCheckSameComm(A, 1, B, 2);
131: PetscCheckSameComm(A, 1, C, 3);
133: PetscAssertPointer(flg, 5);
136: PetscCall(MatGetLocalSize(A, &am, &an));
137: PetscCall(MatGetLocalSize(B, &bm, &bn));
138: PetscCall(MatGetLocalSize(C, &cm, &cn));
139: if (At) {
140: PetscInt tt = an;
141: an = am;
142: am = tt;
143: }
144: if (Bt) {
145: PetscInt tt = bn;
146: bn = bm;
147: bm = tt;
148: }
149: 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);
151: #if PetscDefined(USE_INFO)
152: sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
153: #endif
154: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
155: PetscCall(PetscRandomSetFromOptions(rctx));
156: if (Bt) {
157: PetscCall(MatCreateVecs(B, &s1, &Bx));
158: } else {
159: PetscCall(MatCreateVecs(B, &Bx, &s1));
160: }
161: if (At) {
162: PetscCall(MatCreateVecs(A, &s2, &Ax));
163: } else {
164: PetscCall(MatCreateVecs(A, &Ax, &s2));
165: }
166: PetscCall(MatCreateVecs(C, &Cx, &s3));
168: *flg = PETSC_TRUE;
169: for (k = 0; k < n; k++) {
170: PetscCall(VecSetRandom(Bx, rctx));
171: if (Bt) {
172: PetscCall(MatMultTranspose(B, Bx, s1));
173: } else {
174: PetscCall(MatMult(B, Bx, s1));
175: }
176: PetscCall(VecCopy(s1, Ax));
177: if (At) {
178: PetscCall(MatMultTranspose(A, Ax, s2));
179: } else {
180: PetscCall(MatMult(A, Ax, s2));
181: }
182: PetscCall(VecCopy(Bx, Cx));
183: PetscCall(MatMult(C, Cx, s3));
185: PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
186: if (r2 < tol) {
187: PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
188: } else {
189: PetscCall(VecAXPY(s2, -1.0, s3));
190: PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
191: r1 /= r2;
192: }
193: if (r1 > tol) {
194: *flg = PETSC_FALSE;
195: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
196: break;
197: }
198: }
199: PetscCall(PetscRandomDestroy(&rctx));
200: PetscCall(VecDestroy(&Ax));
201: PetscCall(VecDestroy(&Bx));
202: PetscCall(VecDestroy(&Cx));
203: PetscCall(VecDestroy(&s1));
204: PetscCall(VecDestroy(&s2));
205: PetscCall(VecDestroy(&s3));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*@
210: MatMultEqual - Compares matrix-vector products of two matrices using `n` random vectors
212: Collective
214: Input Parameters:
215: + A - the first matrix
216: . B - the second matrix
217: - n - number of random vectors to be tested
219: Output Parameter:
220: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
222: Level: intermediate
224: Note:
225: The tolerance for equality is a generous `PETSC_SQRT_MACHINE_EPSILON` in the norm of the difference of the two computed vectors to
226: allow for differences in the numerical computations. Hence this routine may indicate equality even if there is a small systematic difference
227: between the two matrices.
229: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`, `MatEqual()`
230: @*/
231: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
232: {
233: PetscFunctionBegin;
234: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 0));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.
241: Collective
243: Input Parameters:
244: + A - the first matrix
245: . B - the second matrix
246: - n - number of random vectors to be tested
248: Output Parameter:
249: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
251: Level: intermediate
253: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
254: @*/
255: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
256: {
257: PetscFunctionBegin;
258: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 1));
259: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 2));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
266: Collective
268: Input Parameters:
269: + A - the first matrix
270: . B - the second matrix
271: - n - number of random vectors to be tested
273: Output Parameter:
274: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
276: Level: intermediate
278: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
279: @*/
280: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
281: {
282: PetscFunctionBegin;
283: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 0));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@
288: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
290: Collective
292: Input Parameters:
293: + A - the first matrix
294: . B - the second matrix
295: - n - number of random vectors to be tested
297: Output Parameter:
298: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
300: Level: intermediate
302: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
303: @*/
304: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
305: {
306: PetscFunctionBegin;
307: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 1));
308: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 2));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@
313: MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
315: Collective
317: Input Parameters:
318: + A - the first matrix
319: . B - the second matrix
320: - n - number of random vectors to be tested
322: Output Parameter:
323: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
325: Level: intermediate
327: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
328: @*/
329: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
330: {
331: PetscFunctionBegin;
332: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 0));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@
337: MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
339: Collective
341: Input Parameters:
342: + A - the first matrix
343: . B - the second matrix
344: - n - number of random vectors to be tested
346: Output Parameter:
347: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
349: Level: intermediate
351: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
352: @*/
353: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
354: {
355: PetscFunctionBegin;
356: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 1));
357: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 2));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*@
362: MatMatMultEqual - Test A*B*x = C*x for n random vector x
364: Collective
366: Input Parameters:
367: + A - the first matrix
368: . B - the second matrix
369: . C - the third matrix
370: - n - number of random vectors to be tested
372: Output Parameter:
373: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
375: Level: intermediate
377: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
378: @*/
379: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
380: {
381: PetscFunctionBegin;
382: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
389: Collective
391: Input Parameters:
392: + A - the first matrix
393: . B - the second matrix
394: . C - the third matrix
395: - n - number of random vectors to be tested
397: Output Parameter:
398: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
400: Level: intermediate
402: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
403: @*/
404: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
405: {
406: PetscFunctionBegin;
407: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: /*@
412: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
414: Collective
416: Input Parameters:
417: + A - the first matrix
418: . B - the second matrix
419: . C - the third matrix
420: - n - number of random vectors to be tested
422: Output Parameter:
423: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
425: Level: intermediate
427: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
428: @*/
429: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
430: {
431: PetscFunctionBegin;
432: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
437: {
438: Vec x, v1, v2, v3, v4, Cx, Bx;
439: PetscReal norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
440: PetscInt i, am, an, bm, bn, cm, cn;
441: PetscRandom rdm;
443: PetscFunctionBegin;
444: PetscCall(MatGetLocalSize(A, &am, &an));
445: PetscCall(MatGetLocalSize(B, &bm, &bn));
446: if (rart) {
447: PetscInt t = bm;
448: bm = bn;
449: bn = t;
450: }
451: PetscCall(MatGetLocalSize(C, &cm, &cn));
452: 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);
454: /* Create left vector of A: v2 */
455: PetscCall(MatCreateVecs(A, &Bx, &v2));
457: /* Create right vectors of B: x, v3, v4 */
458: if (rart) {
459: PetscCall(MatCreateVecs(B, &v1, &x));
460: } else {
461: PetscCall(MatCreateVecs(B, &x, &v1));
462: }
463: PetscCall(VecDuplicate(x, &v3));
465: PetscCall(MatCreateVecs(C, &Cx, &v4));
466: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
467: PetscCall(PetscRandomSetFromOptions(rdm));
469: *flg = PETSC_TRUE;
470: for (i = 0; i < n; i++) {
471: PetscCall(VecSetRandom(x, rdm));
472: PetscCall(VecCopy(x, Cx));
473: PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x */
474: if (rart) {
475: PetscCall(MatMultTranspose(B, x, v1));
476: } else {
477: PetscCall(MatMult(B, x, v1));
478: }
479: PetscCall(VecCopy(v1, Bx));
480: PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
481: PetscCall(VecCopy(v2, v1));
482: if (rart) {
483: PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
484: } else {
485: PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
486: }
487: PetscCall(VecNorm(v4, NORM_2, &norm_abs));
488: PetscCall(VecAXPY(v4, -1.0, v3));
489: PetscCall(VecNorm(v4, NORM_2, &norm_rel));
491: if (norm_abs > tol) norm_rel /= norm_abs;
492: if (norm_rel > tol || PetscIsInfOrNanReal(norm_rel)) {
493: *flg = PETSC_FALSE;
494: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
495: break;
496: }
497: }
499: PetscCall(PetscRandomDestroy(&rdm));
500: PetscCall(VecDestroy(&x));
501: PetscCall(VecDestroy(&Bx));
502: PetscCall(VecDestroy(&Cx));
503: PetscCall(VecDestroy(&v1));
504: PetscCall(VecDestroy(&v2));
505: PetscCall(VecDestroy(&v3));
506: PetscCall(VecDestroy(&v4));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: /*@
511: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
513: Collective
515: Input Parameters:
516: + A - the first matrix
517: . B - the second matrix
518: . C - the third matrix
519: - n - number of random vectors to be tested
521: Output Parameter:
522: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
524: Level: intermediate
526: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
527: @*/
528: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
529: {
530: PetscFunctionBegin;
531: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*@
536: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
538: Collective
540: Input Parameters:
541: + A - the first matrix
542: . B - the second matrix
543: . C - the third matrix
544: - n - number of random vectors to be tested
546: Output Parameter:
547: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
549: Level: intermediate
551: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
552: @*/
553: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
554: {
555: PetscFunctionBegin;
556: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: /*@
561: MatIsLinear - Check if a shell matrix `A` is a linear operator.
563: Collective
565: Input Parameters:
566: + A - the shell matrix
567: - n - number of random vectors to be tested
569: Output Parameter:
570: . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.
572: Level: intermediate
574: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
575: @*/
576: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
577: {
578: Vec x, y, s1, s2;
579: PetscRandom rctx;
580: PetscScalar a;
581: PetscReal norm, normA;
582: MPI_Comm comm;
583: PetscMPIInt rank;
585: PetscFunctionBegin;
587: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
588: PetscCallMPI(MPI_Comm_rank(comm, &rank));
590: PetscCall(PetscRandomCreate(comm, &rctx));
591: PetscCall(PetscRandomSetFromOptions(rctx));
592: PetscCall(MatCreateVecs(A, &x, &s1));
593: PetscCall(VecDuplicate(x, &y));
594: PetscCall(VecDuplicate(s1, &s2));
596: *flg = PETSC_TRUE;
597: for (PetscInt k = 0; k < n; k++) {
598: PetscCall(VecSetRandom(x, rctx));
599: PetscCall(VecSetRandom(y, rctx));
600: if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
601: PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));
603: /* s2 = a*A*x + A*y */
604: PetscCall(MatMult(A, y, s2)); /* s2 = A*y */
605: PetscCall(MatMult(A, x, s1)); /* s1 = A*x */
606: PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */
608: /* s1 = A * (a x + y) */
609: PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
610: PetscCall(MatMult(A, y, s1));
611: PetscCall(VecNorm(s1, NORM_INFINITY, &normA));
613: PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
614: PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
615: if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
616: *flg = PETSC_FALSE;
617: 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)));
618: break;
619: }
620: }
621: PetscCall(PetscRandomDestroy(&rctx));
622: PetscCall(VecDestroy(&x));
623: PetscCall(VecDestroy(&y));
624: PetscCall(VecDestroy(&s1));
625: PetscCall(VecDestroy(&s2));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }