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;
445: PetscFunctionBegin;
446: PetscCall(MatGetLocalSize(A, &am, &an));
447: PetscCall(MatGetLocalSize(B, &bm, &bn));
448: if (rart) {
449: PetscInt t = bm;
450: bm = bn;
451: bn = t;
452: }
453: PetscCall(MatGetLocalSize(C, &cm, &cn));
454: 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);
456: /* Create left vector of A: v2 */
457: PetscCall(MatCreateVecs(A, &Bx, &v2));
459: /* Create right vectors of B: x, v3, v4 */
460: if (rart) {
461: PetscCall(MatCreateVecs(B, &v1, &x));
462: } else {
463: PetscCall(MatCreateVecs(B, &x, &v1));
464: }
465: PetscCall(VecDuplicate(x, &v3));
467: PetscCall(MatCreateVecs(C, &Cx, &v4));
468: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
469: PetscCall(PetscRandomSetFromOptions(rdm));
471: *flg = PETSC_TRUE;
472: for (i = 0; i < n; i++) {
473: PetscCall(VecSetRandom(x, rdm));
474: PetscCall(VecCopy(x, Cx));
475: PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x */
476: if (rart) {
477: PetscCall(MatMultTranspose(B, x, v1));
478: } else {
479: PetscCall(MatMult(B, x, v1));
480: }
481: PetscCall(VecCopy(v1, Bx));
482: PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
483: PetscCall(VecCopy(v2, v1));
484: if (rart) {
485: PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
486: } else {
487: PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
488: }
489: PetscCall(VecNorm(v4, NORM_2, &norm_abs));
490: PetscCall(VecAXPY(v4, -1.0, v3));
491: PetscCall(VecNorm(v4, NORM_2, &norm_rel));
493: if (norm_abs > tol) norm_rel /= norm_abs;
494: if (norm_rel > tol || PetscIsInfOrNanReal(norm_rel)) {
495: *flg = PETSC_FALSE;
496: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
497: break;
498: }
499: }
501: PetscCall(PetscRandomDestroy(&rdm));
502: PetscCall(VecDestroy(&x));
503: PetscCall(VecDestroy(&Bx));
504: PetscCall(VecDestroy(&Cx));
505: PetscCall(VecDestroy(&v1));
506: PetscCall(VecDestroy(&v2));
507: PetscCall(VecDestroy(&v3));
508: PetscCall(VecDestroy(&v4));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: /*@
513: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
515: Collective
517: Input Parameters:
518: + A - the first matrix
519: . B - the second matrix
520: . C - the third matrix
521: - n - number of random vectors to be tested
523: Output Parameter:
524: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
526: Level: intermediate
528: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
529: @*/
530: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
531: {
532: PetscFunctionBegin;
533: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
540: Collective
542: Input Parameters:
543: + A - the first matrix
544: . B - the second matrix
545: . C - the third matrix
546: - n - number of random vectors to be tested
548: Output Parameter:
549: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
551: Level: intermediate
553: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
554: @*/
555: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
556: {
557: PetscFunctionBegin;
558: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: MatIsLinear - Check if a shell matrix `A` is a linear operator.
565: Collective
567: Input Parameters:
568: + A - the shell matrix
569: - n - number of random vectors to be tested
571: Output Parameter:
572: . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.
574: Level: intermediate
576: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
577: @*/
578: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
579: {
580: Vec x, y, s1, s2;
581: PetscRandom rctx;
582: PetscScalar a;
583: PetscInt k;
584: PetscReal norm, normA;
585: MPI_Comm comm;
586: PetscMPIInt rank;
588: PetscFunctionBegin;
590: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
591: PetscCallMPI(MPI_Comm_rank(comm, &rank));
593: PetscCall(PetscRandomCreate(comm, &rctx));
594: PetscCall(PetscRandomSetFromOptions(rctx));
595: PetscCall(MatCreateVecs(A, &x, &s1));
596: PetscCall(VecDuplicate(x, &y));
597: PetscCall(VecDuplicate(s1, &s2));
599: *flg = PETSC_TRUE;
600: for (k = 0; k < n; k++) {
601: PetscCall(VecSetRandom(x, rctx));
602: PetscCall(VecSetRandom(y, rctx));
603: if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
604: PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));
606: /* s2 = a*A*x + A*y */
607: PetscCall(MatMult(A, y, s2)); /* s2 = A*y */
608: PetscCall(MatMult(A, x, s1)); /* s1 = A*x */
609: PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */
611: /* s1 = A * (a x + y) */
612: PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
613: PetscCall(MatMult(A, y, s1));
614: PetscCall(VecNorm(s1, NORM_INFINITY, &normA));
616: PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
617: PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
618: if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
619: *flg = PETSC_FALSE;
620: 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)));
621: break;
622: }
623: }
624: PetscCall(PetscRandomDestroy(&rctx));
625: PetscCall(VecDestroy(&x));
626: PetscCall(VecDestroy(&y));
627: PetscCall(VecDestroy(&s1));
628: PetscCall(VecDestroy(&s2));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }