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.
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: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`
227: @*/
228: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
229: {
230: PetscFunctionBegin;
231: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 0));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@
236: MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.
238: Collective
240: Input Parameters:
241: + A - the first matrix
242: . B - the second matrix
243: - n - number of random vectors to be tested
245: Output Parameter:
246: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
248: Level: intermediate
250: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
251: @*/
252: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
253: {
254: PetscFunctionBegin;
255: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 1));
256: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 2));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@
261: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
263: Collective
265: Input Parameters:
266: + A - the first matrix
267: . B - the second matrix
268: - n - number of random vectors to be tested
270: Output Parameter:
271: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
273: Level: intermediate
275: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
276: @*/
277: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
278: {
279: PetscFunctionBegin;
280: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 0));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@
285: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
287: Collective
289: Input Parameters:
290: + A - the first matrix
291: . B - the second matrix
292: - n - number of random vectors to be tested
294: Output Parameter:
295: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
297: Level: intermediate
299: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
300: @*/
301: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
302: {
303: PetscFunctionBegin;
304: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 1));
305: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 2));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
312: Collective
314: Input Parameters:
315: + A - the first matrix
316: . B - the second matrix
317: - n - number of random vectors to be tested
319: Output Parameter:
320: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
322: Level: intermediate
324: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
325: @*/
326: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
327: {
328: PetscFunctionBegin;
329: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 0));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /*@
334: MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
336: Collective
338: Input Parameters:
339: + A - the first matrix
340: . B - the second matrix
341: - n - number of random vectors to be tested
343: Output Parameter:
344: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
346: Level: intermediate
348: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
349: @*/
350: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
351: {
352: PetscFunctionBegin;
353: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 1));
354: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 2));
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@
359: MatMatMultEqual - Test A*B*x = C*x for n random vector x
361: Collective
363: Input Parameters:
364: + A - the first matrix
365: . B - the second matrix
366: . C - the third matrix
367: - n - number of random vectors to be tested
369: Output Parameter:
370: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
372: Level: intermediate
374: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
375: @*/
376: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
377: {
378: PetscFunctionBegin;
379: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: /*@
384: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
386: Collective
388: Input Parameters:
389: + A - the first matrix
390: . B - the second matrix
391: . C - the third matrix
392: - n - number of random vectors to be tested
394: Output Parameter:
395: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
397: Level: intermediate
399: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
400: @*/
401: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
402: {
403: PetscFunctionBegin;
404: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
411: Collective
413: Input Parameters:
414: + A - the first matrix
415: . B - the second matrix
416: . C - the third matrix
417: - n - number of random vectors to be tested
419: Output Parameter:
420: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
422: Level: intermediate
424: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
425: @*/
426: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
427: {
428: PetscFunctionBegin;
429: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
434: {
435: Vec x, v1, v2, v3, v4, Cx, Bx;
436: PetscReal norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
437: PetscInt i, am, an, bm, bn, cm, cn;
438: PetscRandom rdm;
439: PetscScalar none = -1.0;
441: PetscFunctionBegin;
442: PetscCall(MatGetLocalSize(A, &am, &an));
443: PetscCall(MatGetLocalSize(B, &bm, &bn));
444: if (rart) {
445: PetscInt t = bm;
446: bm = bn;
447: bn = t;
448: }
449: PetscCall(MatGetLocalSize(C, &cm, &cn));
450: 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);
452: /* Create left vector of A: v2 */
453: PetscCall(MatCreateVecs(A, &Bx, &v2));
455: /* Create right vectors of B: x, v3, v4 */
456: if (rart) {
457: PetscCall(MatCreateVecs(B, &v1, &x));
458: } else {
459: PetscCall(MatCreateVecs(B, &x, &v1));
460: }
461: PetscCall(VecDuplicate(x, &v3));
463: PetscCall(MatCreateVecs(C, &Cx, &v4));
464: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
465: PetscCall(PetscRandomSetFromOptions(rdm));
467: *flg = PETSC_TRUE;
468: for (i = 0; i < n; i++) {
469: PetscCall(VecSetRandom(x, rdm));
470: PetscCall(VecCopy(x, Cx));
471: PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x */
472: if (rart) {
473: PetscCall(MatMultTranspose(B, x, v1));
474: } else {
475: PetscCall(MatMult(B, x, v1));
476: }
477: PetscCall(VecCopy(v1, Bx));
478: PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
479: PetscCall(VecCopy(v2, v1));
480: if (rart) {
481: PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
482: } else {
483: PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
484: }
485: PetscCall(VecNorm(v4, NORM_2, &norm_abs));
486: PetscCall(VecAXPY(v4, none, v3));
487: PetscCall(VecNorm(v4, NORM_2, &norm_rel));
489: if (norm_abs > tol) norm_rel /= norm_abs;
490: if (norm_rel > tol) {
491: *flg = PETSC_FALSE;
492: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
493: break;
494: }
495: }
497: PetscCall(PetscRandomDestroy(&rdm));
498: PetscCall(VecDestroy(&x));
499: PetscCall(VecDestroy(&Bx));
500: PetscCall(VecDestroy(&Cx));
501: PetscCall(VecDestroy(&v1));
502: PetscCall(VecDestroy(&v2));
503: PetscCall(VecDestroy(&v3));
504: PetscCall(VecDestroy(&v4));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: /*@
509: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
511: Collective
513: Input Parameters:
514: + A - the first matrix
515: . B - the second matrix
516: . C - the third matrix
517: - n - number of random vectors to be tested
519: Output Parameter:
520: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
522: Level: intermediate
524: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
525: @*/
526: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
527: {
528: PetscFunctionBegin;
529: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@
534: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
536: Collective
538: Input Parameters:
539: + A - the first matrix
540: . B - the second matrix
541: . C - the third matrix
542: - n - number of random vectors to be tested
544: Output Parameter:
545: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
547: Level: intermediate
549: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
550: @*/
551: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
552: {
553: PetscFunctionBegin;
554: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: MatIsLinear - Check if a shell matrix `A` is a linear operator.
561: Collective
563: Input Parameters:
564: + A - the shell matrix
565: - n - number of random vectors to be tested
567: Output Parameter:
568: . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.
570: Level: intermediate
572: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
573: @*/
574: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
575: {
576: Vec x, y, s1, s2;
577: PetscRandom rctx;
578: PetscScalar a;
579: PetscInt k;
580: PetscReal norm, normA;
581: MPI_Comm comm;
582: PetscMPIInt rank;
584: PetscFunctionBegin;
586: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
587: PetscCallMPI(MPI_Comm_rank(comm, &rank));
589: PetscCall(PetscRandomCreate(comm, &rctx));
590: PetscCall(PetscRandomSetFromOptions(rctx));
591: PetscCall(MatCreateVecs(A, &x, &s1));
592: PetscCall(VecDuplicate(x, &y));
593: PetscCall(VecDuplicate(s1, &s2));
595: *flg = PETSC_TRUE;
596: for (k = 0; k < n; k++) {
597: PetscCall(VecSetRandom(x, rctx));
598: PetscCall(VecSetRandom(y, rctx));
599: if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
600: PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));
602: /* s2 = a*A*x + A*y */
603: PetscCall(MatMult(A, y, s2)); /* s2 = A*y */
604: PetscCall(MatMult(A, x, s1)); /* s1 = A*x */
605: PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */
607: /* s1 = A * (a x + y) */
608: PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
609: PetscCall(MatMult(A, y, s1));
610: PetscCall(VecNorm(s1, NORM_INFINITY, &normA));
612: PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
613: PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
614: if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
615: *flg = PETSC_FALSE;
616: 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)));
617: break;
618: }
619: }
620: PetscCall(PetscRandomDestroy(&rctx));
621: PetscCall(VecDestroy(&x));
622: PetscCall(VecDestroy(&y));
623: PetscCall(VecDestroy(&s1));
624: PetscCall(VecDestroy(&s2));
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }