Actual source code: multequal.c
2: #include <petsc/private/matimpl.h>
4: static PetscErrorCode MatMultEqual_Private(Mat A, Mat B, PetscInt n, PetscBool *flg, PetscInt t, PetscBool add)
5: {
6: Vec Ax = NULL, Bx = NULL, s1 = NULL, s2 = NULL, Ay = NULL, By = NULL;
7: PetscRandom rctx;
8: PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
9: PetscInt am, an, bm, bn, k;
10: PetscScalar none = -1.0;
11: #if defined(PETSC_USE_INFO)
12: const char *sops[] = {"MatMult", "MatMultAdd", "MatMultTranspose", "MatMultTransposeAdd", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd"};
13: const char *sop;
14: #endif
16: PetscFunctionBegin;
19: PetscCheckSameComm(A, 1, B, 2);
24: PetscCall(MatGetLocalSize(A, &am, &an));
25: PetscCall(MatGetLocalSize(B, &bm, &bn));
26: 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);
27: #if defined(PETSC_USE_INFO)
28: sop = sops[(add ? 1 : 0) + 2 * t]; /* t = 0 => no transpose, t = 1 => transpose, t = 2 => Hermitian transpose */
29: #endif
30: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
31: PetscCall(PetscRandomSetFromOptions(rctx));
32: if (t) {
33: PetscCall(MatCreateVecs(A, &s1, &Ax));
34: PetscCall(MatCreateVecs(B, &s2, &Bx));
35: } else {
36: PetscCall(MatCreateVecs(A, &Ax, &s1));
37: PetscCall(MatCreateVecs(B, &Bx, &s2));
38: }
39: if (add) {
40: PetscCall(VecDuplicate(s1, &Ay));
41: PetscCall(VecDuplicate(s2, &By));
42: }
44: *flg = PETSC_TRUE;
45: for (k = 0; k < n; k++) {
46: PetscCall(VecSetRandom(Ax, rctx));
47: PetscCall(VecCopy(Ax, Bx));
48: if (add) {
49: PetscCall(VecSetRandom(Ay, rctx));
50: PetscCall(VecCopy(Ay, By));
51: }
52: if (t == 1) {
53: if (add) {
54: PetscCall(MatMultTransposeAdd(A, Ax, Ay, s1));
55: PetscCall(MatMultTransposeAdd(B, Bx, By, s2));
56: } else {
57: PetscCall(MatMultTranspose(A, Ax, s1));
58: PetscCall(MatMultTranspose(B, Bx, s2));
59: }
60: } else if (t == 2) {
61: if (add) {
62: PetscCall(MatMultHermitianTransposeAdd(A, Ax, Ay, s1));
63: PetscCall(MatMultHermitianTransposeAdd(B, Bx, By, s2));
64: } else {
65: PetscCall(MatMultHermitianTranspose(A, Ax, s1));
66: PetscCall(MatMultHermitianTranspose(B, Bx, s2));
67: }
68: } else {
69: if (add) {
70: PetscCall(MatMultAdd(A, Ax, Ay, s1));
71: PetscCall(MatMultAdd(B, Bx, By, s2));
72: } else {
73: PetscCall(MatMult(A, Ax, s1));
74: PetscCall(MatMult(B, Bx, s2));
75: }
76: }
77: PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
78: if (r2 < tol) {
79: PetscCall(VecNorm(s1, NORM_INFINITY, &r1));
80: } else {
81: PetscCall(VecAXPY(s2, none, s1));
82: PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
83: r1 /= r2;
84: }
85: if (r1 > tol) {
86: *flg = PETSC_FALSE;
87: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1));
88: break;
89: }
90: }
91: PetscCall(PetscRandomDestroy(&rctx));
92: PetscCall(VecDestroy(&Ax));
93: PetscCall(VecDestroy(&Bx));
94: PetscCall(VecDestroy(&Ay));
95: PetscCall(VecDestroy(&By));
96: PetscCall(VecDestroy(&s1));
97: PetscCall(VecDestroy(&s2));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
102: {
103: Vec Ax, Bx, Cx, s1, s2, s3;
104: PetscRandom rctx;
105: PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
106: PetscInt am, an, bm, bn, cm, cn, k;
107: PetscScalar none = -1.0;
108: #if defined(PETSC_USE_INFO)
109: const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
110: const char *sop;
111: #endif
113: PetscFunctionBegin;
116: PetscCheckSameComm(A, 1, B, 2);
118: PetscCheckSameComm(A, 1, C, 3);
123: PetscCall(MatGetLocalSize(A, &am, &an));
124: PetscCall(MatGetLocalSize(B, &bm, &bn));
125: PetscCall(MatGetLocalSize(C, &cm, &cn));
126: if (At) {
127: PetscInt tt = an;
128: an = am;
129: am = tt;
130: };
131: if (Bt) {
132: PetscInt tt = bn;
133: bn = bm;
134: bm = tt;
135: };
136: 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);
138: #if defined(PETSC_USE_INFO)
139: sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
140: #endif
141: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
142: PetscCall(PetscRandomSetFromOptions(rctx));
143: if (Bt) {
144: PetscCall(MatCreateVecs(B, &s1, &Bx));
145: } else {
146: PetscCall(MatCreateVecs(B, &Bx, &s1));
147: }
148: if (At) {
149: PetscCall(MatCreateVecs(A, &s2, &Ax));
150: } else {
151: PetscCall(MatCreateVecs(A, &Ax, &s2));
152: }
153: PetscCall(MatCreateVecs(C, &Cx, &s3));
155: *flg = PETSC_TRUE;
156: for (k = 0; k < n; k++) {
157: PetscCall(VecSetRandom(Bx, rctx));
158: if (Bt) {
159: PetscCall(MatMultTranspose(B, Bx, s1));
160: } else {
161: PetscCall(MatMult(B, Bx, s1));
162: }
163: PetscCall(VecCopy(s1, Ax));
164: if (At) {
165: PetscCall(MatMultTranspose(A, Ax, s2));
166: } else {
167: PetscCall(MatMult(A, Ax, s2));
168: }
169: PetscCall(VecCopy(Bx, Cx));
170: PetscCall(MatMult(C, Cx, s3));
172: PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
173: if (r2 < tol) {
174: PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
175: } else {
176: PetscCall(VecAXPY(s2, none, s3));
177: PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
178: r1 /= r2;
179: }
180: if (r1 > tol) {
181: *flg = PETSC_FALSE;
182: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
183: break;
184: }
185: }
186: PetscCall(PetscRandomDestroy(&rctx));
187: PetscCall(VecDestroy(&Ax));
188: PetscCall(VecDestroy(&Bx));
189: PetscCall(VecDestroy(&Cx));
190: PetscCall(VecDestroy(&s1));
191: PetscCall(VecDestroy(&s2));
192: PetscCall(VecDestroy(&s3));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: /*@
197: MatMultEqual - Compares matrix-vector products of two matrices.
199: Collective
201: Input Parameters:
202: + A - the first matrix
203: . B - the second matrix
204: - n - number of random vectors to be tested
206: Output Parameter:
207: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
209: Level: intermediate
211: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`
212: @*/
213: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
214: {
215: PetscFunctionBegin;
216: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, PETSC_FALSE));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@
221: MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.
223: Collective
225: Input Parameters:
226: + A - the first matrix
227: . B - the second matrix
228: - n - number of random vectors to be tested
230: Output Parameter:
231: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
233: Level: intermediate
235: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
236: @*/
237: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
238: {
239: PetscFunctionBegin;
240: PetscCall(MatMultEqual_Private(A, B, n, flg, 0, PETSC_TRUE));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*@
245: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
247: Collective
249: Input Parameters:
250: + A - the first matrix
251: . B - the second matrix
252: - n - number of random vectors to be tested
254: Output Parameter:
255: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
257: Level: intermediate
259: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
260: @*/
261: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
262: {
263: PetscFunctionBegin;
264: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, PETSC_FALSE));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
271: Collective
273: Input Parameters:
274: + A - the first matrix
275: . B - the second matrix
276: - n - number of random vectors to be tested
278: Output Parameter:
279: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
281: Level: intermediate
283: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
284: @*/
285: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
286: {
287: PetscFunctionBegin;
288: PetscCall(MatMultEqual_Private(A, B, n, flg, 1, PETSC_TRUE));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /*@
293: MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
295: Collective
297: Input Parameters:
298: + A - the first matrix
299: . B - the second matrix
300: - n - number of random vectors to be tested
302: Output Parameter:
303: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
305: Level: intermediate
307: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
308: @*/
309: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
310: {
311: PetscFunctionBegin;
312: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, PETSC_FALSE));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*@
317: MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
319: Collective
321: Input Parameters:
322: + A - the first matrix
323: . B - the second matrix
324: - n - number of random vectors to be tested
326: Output Parameter:
327: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
329: Level: intermediate
331: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
332: @*/
333: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
334: {
335: PetscFunctionBegin;
336: PetscCall(MatMultEqual_Private(A, B, n, flg, 2, PETSC_TRUE));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: MatMatMultEqual - Test A*B*x = C*x for n random vector x
343: Collective
345: Input Parameters:
346: + A - the first matrix
347: . B - the second matrix
348: . C - the third matrix
349: - n - number of random vectors to be tested
351: Output Parameter:
352: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
354: Level: intermediate
356: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
357: @*/
358: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
359: {
360: PetscFunctionBegin;
361: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
368: Collective
370: Input Parameters:
371: + A - the first matrix
372: . B - the second matrix
373: . C - the third matrix
374: - n - number of random vectors to be tested
376: Output Parameter:
377: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
379: Level: intermediate
381: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
382: @*/
383: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
384: {
385: PetscFunctionBegin;
386: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@
391: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
393: Collective
395: Input Parameters:
396: + A - the first matrix
397: . B - the second matrix
398: . C - the third matrix
399: - n - number of random vectors to be tested
401: Output Parameter:
402: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
404: Level: intermediate
406: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
407: @*/
408: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
409: {
410: PetscFunctionBegin;
411: PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
416: {
417: Vec x, v1, v2, v3, v4, Cx, Bx;
418: PetscReal norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
419: PetscInt i, am, an, bm, bn, cm, cn;
420: PetscRandom rdm;
421: PetscScalar none = -1.0;
423: PetscFunctionBegin;
424: PetscCall(MatGetLocalSize(A, &am, &an));
425: PetscCall(MatGetLocalSize(B, &bm, &bn));
426: if (rart) {
427: PetscInt t = bm;
428: bm = bn;
429: bn = t;
430: }
431: PetscCall(MatGetLocalSize(C, &cm, &cn));
432: 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);
434: /* Create left vector of A: v2 */
435: PetscCall(MatCreateVecs(A, &Bx, &v2));
437: /* Create right vectors of B: x, v3, v4 */
438: if (rart) {
439: PetscCall(MatCreateVecs(B, &v1, &x));
440: } else {
441: PetscCall(MatCreateVecs(B, &x, &v1));
442: }
443: PetscCall(VecDuplicate(x, &v3));
445: PetscCall(MatCreateVecs(C, &Cx, &v4));
446: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
447: PetscCall(PetscRandomSetFromOptions(rdm));
449: *flg = PETSC_TRUE;
450: for (i = 0; i < n; i++) {
451: PetscCall(VecSetRandom(x, rdm));
452: PetscCall(VecCopy(x, Cx));
453: PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x */
454: if (rart) {
455: PetscCall(MatMultTranspose(B, x, v1));
456: } else {
457: PetscCall(MatMult(B, x, v1));
458: }
459: PetscCall(VecCopy(v1, Bx));
460: PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
461: PetscCall(VecCopy(v2, v1));
462: if (rart) {
463: PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
464: } else {
465: PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
466: }
467: PetscCall(VecNorm(v4, NORM_2, &norm_abs));
468: PetscCall(VecAXPY(v4, none, v3));
469: PetscCall(VecNorm(v4, NORM_2, &norm_rel));
471: if (norm_abs > tol) norm_rel /= norm_abs;
472: if (norm_rel > tol) {
473: *flg = PETSC_FALSE;
474: PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
475: break;
476: }
477: }
479: PetscCall(PetscRandomDestroy(&rdm));
480: PetscCall(VecDestroy(&x));
481: PetscCall(VecDestroy(&Bx));
482: PetscCall(VecDestroy(&Cx));
483: PetscCall(VecDestroy(&v1));
484: PetscCall(VecDestroy(&v2));
485: PetscCall(VecDestroy(&v3));
486: PetscCall(VecDestroy(&v4));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
493: Collective
495: Input Parameters:
496: + A - the first matrix
497: . B - the second matrix
498: . C - the third matrix
499: - n - number of random vectors to be tested
501: Output Parameter:
502: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
504: Level: intermediate
506: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
507: @*/
508: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
509: {
510: PetscFunctionBegin;
511: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: /*@
516: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
518: Collective
520: Input Parameters:
521: + A - the first matrix
522: . B - the second matrix
523: . C - the third matrix
524: - n - number of random vectors to be tested
526: Output Parameter:
527: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.
529: Level: intermediate
531: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
532: @*/
533: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
534: {
535: PetscFunctionBegin;
536: PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: /*@
541: MatIsLinear - Check if a shell matrix `A` is a linear operator.
543: Collective
545: Input Parameters:
546: + A - the shell matrix
547: - n - number of random vectors to be tested
549: Output Parameter:
550: . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.
552: Level: intermediate
554: .seealso: `Mat`, `MatMatMultAddEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
555: @*/
556: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
557: {
558: Vec x, y, s1, s2;
559: PetscRandom rctx;
560: PetscScalar a;
561: PetscInt k;
562: PetscReal norm, normA;
563: MPI_Comm comm;
564: PetscMPIInt rank;
566: PetscFunctionBegin;
568: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
569: PetscCallMPI(MPI_Comm_rank(comm, &rank));
571: PetscCall(PetscRandomCreate(comm, &rctx));
572: PetscCall(PetscRandomSetFromOptions(rctx));
573: PetscCall(MatCreateVecs(A, &x, &s1));
574: PetscCall(VecDuplicate(x, &y));
575: PetscCall(VecDuplicate(s1, &s2));
577: *flg = PETSC_TRUE;
578: for (k = 0; k < n; k++) {
579: PetscCall(VecSetRandom(x, rctx));
580: PetscCall(VecSetRandom(y, rctx));
581: if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
582: PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));
584: /* s2 = a*A*x + A*y */
585: PetscCall(MatMult(A, y, s2)); /* s2 = A*y */
586: PetscCall(MatMult(A, x, s1)); /* s1 = A*x */
587: PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */
589: /* s1 = A * (a x + y) */
590: PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
591: PetscCall(MatMult(A, y, s1));
592: PetscCall(VecNorm(s1, NORM_INFINITY, &normA));
594: PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
595: PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
596: if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
597: *flg = PETSC_FALSE;
598: 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)));
599: break;
600: }
601: }
602: PetscCall(PetscRandomDestroy(&rctx));
603: PetscCall(VecDestroy(&x));
604: PetscCall(VecDestroy(&y));
605: PetscCall(VecDestroy(&s1));
606: PetscCall(VecDestroy(&s2));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }