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: const char *sops[] = {"MatMult", "MatMultAdd", "MatMultTranspose", "MatMultTransposeAdd", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd"};
12: const char *sop;
21: MatGetLocalSize(A, &am, &an);
22: MatGetLocalSize(B, &bm, &bn);
24: sop = sops[(add ? 1 : 0) + 2 * t]; /* t = 0 => no transpose, t = 1 => transpose, t = 2 => Hermitian transpose */
25: PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx);
26: PetscRandomSetFromOptions(rctx);
27: if (t) {
28: MatCreateVecs(A, &s1, &Ax);
29: MatCreateVecs(B, &s2, &Bx);
30: } else {
31: MatCreateVecs(A, &Ax, &s1);
32: MatCreateVecs(B, &Bx, &s2);
33: }
34: if (add) {
35: VecDuplicate(s1, &Ay);
36: VecDuplicate(s2, &By);
37: }
39: *flg = PETSC_TRUE;
40: for (k = 0; k < n; k++) {
41: VecSetRandom(Ax, rctx);
42: VecCopy(Ax, Bx);
43: if (add) {
44: VecSetRandom(Ay, rctx);
45: VecCopy(Ay, By);
46: }
47: if (t == 1) {
48: if (add) {
49: MatMultTransposeAdd(A, Ax, Ay, s1);
50: MatMultTransposeAdd(B, Bx, By, s2);
51: } else {
52: MatMultTranspose(A, Ax, s1);
53: MatMultTranspose(B, Bx, s2);
54: }
55: } else if (t == 2) {
56: if (add) {
57: MatMultHermitianTransposeAdd(A, Ax, Ay, s1);
58: MatMultHermitianTransposeAdd(B, Bx, By, s2);
59: } else {
60: MatMultHermitianTranspose(A, Ax, s1);
61: MatMultHermitianTranspose(B, Bx, s2);
62: }
63: } else {
64: if (add) {
65: MatMultAdd(A, Ax, Ay, s1);
66: MatMultAdd(B, Bx, By, s2);
67: } else {
68: MatMult(A, Ax, s1);
69: MatMult(B, Bx, s2);
70: }
71: }
72: VecNorm(s2, NORM_INFINITY, &r2);
73: if (r2 < tol) {
74: VecNorm(s1, NORM_INFINITY, &r1);
75: } else {
76: VecAXPY(s2, none, s1);
77: VecNorm(s2, NORM_INFINITY, &r1);
78: r1 /= r2;
79: }
80: if (r1 > tol) {
81: *flg = PETSC_FALSE;
82: PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1);
83: break;
84: }
85: }
86: PetscRandomDestroy(&rctx);
87: VecDestroy(&Ax);
88: VecDestroy(&Bx);
89: VecDestroy(&Ay);
90: VecDestroy(&By);
91: VecDestroy(&s1);
92: VecDestroy(&s2);
93: return 0;
94: }
96: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
97: {
98: Vec Ax, Bx, Cx, s1, s2, s3;
99: PetscRandom rctx;
100: PetscReal r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
101: PetscInt am, an, bm, bn, cm, cn, k;
102: PetscScalar none = -1.0;
103: const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
104: const char *sop;
115: MatGetLocalSize(A, &am, &an);
116: MatGetLocalSize(B, &bm, &bn);
117: MatGetLocalSize(C, &cm, &cn);
118: if (At) {
119: PetscInt tt = an;
120: an = am;
121: am = tt;
122: };
123: if (Bt) {
124: PetscInt tt = bn;
125: bn = bm;
126: bm = tt;
127: };
130: sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
131: PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx);
132: PetscRandomSetFromOptions(rctx);
133: if (Bt) {
134: MatCreateVecs(B, &s1, &Bx);
135: } else {
136: MatCreateVecs(B, &Bx, &s1);
137: }
138: if (At) {
139: MatCreateVecs(A, &s2, &Ax);
140: } else {
141: MatCreateVecs(A, &Ax, &s2);
142: }
143: MatCreateVecs(C, &Cx, &s3);
145: *flg = PETSC_TRUE;
146: for (k = 0; k < n; k++) {
147: VecSetRandom(Bx, rctx);
148: if (Bt) {
149: MatMultTranspose(B, Bx, s1);
150: } else {
151: MatMult(B, Bx, s1);
152: }
153: VecCopy(s1, Ax);
154: if (At) {
155: MatMultTranspose(A, Ax, s2);
156: } else {
157: MatMult(A, Ax, s2);
158: }
159: VecCopy(Bx, Cx);
160: MatMult(C, Cx, s3);
162: VecNorm(s2, NORM_INFINITY, &r2);
163: if (r2 < tol) {
164: VecNorm(s3, NORM_INFINITY, &r1);
165: } else {
166: VecAXPY(s2, none, s3);
167: VecNorm(s2, NORM_INFINITY, &r1);
168: r1 /= r2;
169: }
170: if (r1 > tol) {
171: *flg = PETSC_FALSE;
172: PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1);
173: break;
174: }
175: }
176: PetscRandomDestroy(&rctx);
177: VecDestroy(&Ax);
178: VecDestroy(&Bx);
179: VecDestroy(&Cx);
180: VecDestroy(&s1);
181: VecDestroy(&s2);
182: VecDestroy(&s3);
183: return 0;
184: }
186: /*@
187: MatMultEqual - Compares matrix-vector products of two matrices.
189: Collective
191: Input Parameters:
192: + A - the first matrix
193: . B - the second matrix
194: - n - number of random vectors to be tested
196: Output Parameter:
197: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
199: Level: intermediate
201: @*/
202: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
203: {
204: MatMultEqual_Private(A, B, n, flg, 0, PETSC_FALSE);
205: return 0;
206: }
208: /*@
209: MatMultAddEqual - Compares matrix-vector products of two matrices.
211: Collective
213: Input Parameters:
214: + A - the first matrix
215: . B - the second matrix
216: - n - number of random vectors to be tested
218: Output Parameter:
219: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
221: Level: intermediate
223: @*/
224: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
225: {
226: MatMultEqual_Private(A, B, n, flg, 0, PETSC_TRUE);
227: return 0;
228: }
230: /*@
231: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
233: Collective
235: Input Parameters:
236: + A - the first matrix
237: . B - the second matrix
238: - n - number of random vectors to be tested
240: Output Parameter:
241: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
243: Level: intermediate
245: @*/
246: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
247: {
248: MatMultEqual_Private(A, B, n, flg, 1, PETSC_FALSE);
249: return 0;
250: }
252: /*@
253: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
255: Collective
257: Input Parameters:
258: + A - the first matrix
259: . B - the second matrix
260: - n - number of random vectors to be tested
262: Output Parameter:
263: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
265: Level: intermediate
267: @*/
268: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
269: {
270: MatMultEqual_Private(A, B, n, flg, 1, PETSC_TRUE);
271: return 0;
272: }
274: /*@
275: MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
277: Collective
279: Input Parameters:
280: + A - the first matrix
281: . B - the second matrix
282: - n - number of random vectors to be tested
284: Output Parameter:
285: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
287: Level: intermediate
289: @*/
290: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
291: {
292: MatMultEqual_Private(A, B, n, flg, 2, PETSC_FALSE);
293: return 0;
294: }
296: /*@
297: MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
299: Collective
301: Input Parameters:
302: + A - the first matrix
303: . B - the second matrix
304: - n - number of random vectors to be tested
306: Output Parameter:
307: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
309: Level: intermediate
311: @*/
312: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
313: {
314: MatMultEqual_Private(A, B, n, flg, 2, PETSC_TRUE);
315: return 0;
316: }
318: /*@
319: MatMatMultEqual - Test A*B*x = C*x for n random vector x
321: Collective
323: Input Parameters:
324: + A - the first matrix
325: . B - the second matrix
326: . C - the third matrix
327: - n - number of random vectors to be tested
329: Output Parameter:
330: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
332: Level: intermediate
334: @*/
335: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
336: {
337: MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE);
338: return 0;
339: }
341: /*@
342: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
344: Collective
346: Input Parameters:
347: + A - the first matrix
348: . B - the second matrix
349: . C - the third matrix
350: - n - number of random vectors to be tested
352: Output Parameter:
353: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
355: Level: intermediate
357: @*/
358: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
359: {
360: MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE);
361: return 0;
362: }
364: /*@
365: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
367: Collective
369: Input Parameters:
370: + A - the first matrix
371: . B - the second matrix
372: . C - the third matrix
373: - n - number of random vectors to be tested
375: Output Parameter:
376: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
378: Level: intermediate
380: @*/
381: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
382: {
383: MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE);
384: return 0;
385: }
387: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
388: {
389: Vec x, v1, v2, v3, v4, Cx, Bx;
390: PetscReal norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
391: PetscInt i, am, an, bm, bn, cm, cn;
392: PetscRandom rdm;
393: PetscScalar none = -1.0;
395: MatGetLocalSize(A, &am, &an);
396: MatGetLocalSize(B, &bm, &bn);
397: if (rart) {
398: PetscInt t = bm;
399: bm = bn;
400: bn = t;
401: }
402: MatGetLocalSize(C, &cm, &cn);
405: /* Create left vector of A: v2 */
406: MatCreateVecs(A, &Bx, &v2);
408: /* Create right vectors of B: x, v3, v4 */
409: if (rart) {
410: MatCreateVecs(B, &v1, &x);
411: } else {
412: MatCreateVecs(B, &x, &v1);
413: }
414: VecDuplicate(x, &v3);
416: MatCreateVecs(C, &Cx, &v4);
417: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
418: PetscRandomSetFromOptions(rdm);
420: *flg = PETSC_TRUE;
421: for (i = 0; i < n; i++) {
422: VecSetRandom(x, rdm);
423: VecCopy(x, Cx);
424: MatMult(C, Cx, v4); /* v4 = C*x */
425: if (rart) {
426: MatMultTranspose(B, x, v1);
427: } else {
428: MatMult(B, x, v1);
429: }
430: VecCopy(v1, Bx);
431: MatMult(A, Bx, v2); /* v2 = A*B*x */
432: VecCopy(v2, v1);
433: if (rart) {
434: MatMult(B, v1, v3); /* v3 = R*A*R^t*x */
435: } else {
436: MatMultTranspose(B, v1, v3); /* v3 = Bt*A*B*x */
437: }
438: VecNorm(v4, NORM_2, &norm_abs);
439: VecAXPY(v4, none, v3);
440: VecNorm(v4, NORM_2, &norm_rel);
442: if (norm_abs > tol) norm_rel /= norm_abs;
443: if (norm_rel > tol) {
444: *flg = PETSC_FALSE;
445: PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel);
446: break;
447: }
448: }
450: PetscRandomDestroy(&rdm);
451: VecDestroy(&x);
452: VecDestroy(&Bx);
453: VecDestroy(&Cx);
454: VecDestroy(&v1);
455: VecDestroy(&v2);
456: VecDestroy(&v3);
457: VecDestroy(&v4);
458: return 0;
459: }
461: /*@
462: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
464: Collective
466: Input Parameters:
467: + A - the first matrix
468: . B - the second matrix
469: . C - the third matrix
470: - n - number of random vectors to be tested
472: Output Parameter:
473: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
475: Level: intermediate
477: @*/
478: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
479: {
480: MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg);
481: return 0;
482: }
484: /*@
485: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
487: Collective
489: Input Parameters:
490: + A - the first matrix
491: . B - the second matrix
492: . C - the third matrix
493: - n - number of random vectors to be tested
495: Output Parameter:
496: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
498: Level: intermediate
500: @*/
501: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
502: {
503: MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg);
504: return 0;
505: }
507: /*@
508: MatIsLinear - Check if a shell matrix A is a linear operator.
510: Collective
512: Input Parameters:
513: + A - the shell matrix
514: - n - number of random vectors to be tested
516: Output Parameter:
517: . flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.
519: Level: intermediate
520: @*/
521: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
522: {
523: Vec x, y, s1, s2;
524: PetscRandom rctx;
525: PetscScalar a;
526: PetscInt k;
527: PetscReal norm, normA;
528: MPI_Comm comm;
529: PetscMPIInt rank;
532: PetscObjectGetComm((PetscObject)A, &comm);
533: MPI_Comm_rank(comm, &rank);
535: PetscRandomCreate(comm, &rctx);
536: PetscRandomSetFromOptions(rctx);
537: MatCreateVecs(A, &x, &s1);
538: VecDuplicate(x, &y);
539: VecDuplicate(s1, &s2);
541: *flg = PETSC_TRUE;
542: for (k = 0; k < n; k++) {
543: VecSetRandom(x, rctx);
544: VecSetRandom(y, rctx);
545: if (rank == 0) PetscRandomGetValue(rctx, &a);
546: MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);
548: /* s2 = a*A*x + A*y */
549: MatMult(A, y, s2); /* s2 = A*y */
550: MatMult(A, x, s1); /* s1 = A*x */
551: VecAXPY(s2, a, s1); /* s2 = a s1 + s2 */
553: /* s1 = A * (a x + y) */
554: VecAXPY(y, a, x); /* y = a x + y */
555: MatMult(A, y, s1);
556: VecNorm(s1, NORM_INFINITY, &normA);
558: VecAXPY(s2, -1.0, s1); /* s2 = - s1 + s2 */
559: VecNorm(s2, NORM_INFINITY, &norm);
560: if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
561: *flg = PETSC_FALSE;
562: 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));
563: break;
564: }
565: }
566: PetscRandomDestroy(&rctx);
567: VecDestroy(&x);
568: VecDestroy(&y);
569: VecDestroy(&s1);
570: VecDestroy(&s2);
571: return 0;
572: }