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) { PetscInt tt = an; an = am; am = tt; };
119: if (Bt) { PetscInt tt = bn; bn = bm; bm = tt; };
122: sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
123: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
124: PetscRandomSetFromOptions(rctx);
125: if (Bt) {
126: MatCreateVecs(B,&s1,&Bx);
127: } else {
128: MatCreateVecs(B,&Bx,&s1);
129: }
130: if (At) {
131: MatCreateVecs(A,&s2,&Ax);
132: } else {
133: MatCreateVecs(A,&Ax,&s2);
134: }
135: MatCreateVecs(C,&Cx,&s3);
137: *flg = PETSC_TRUE;
138: for (k=0; k<n; k++) {
139: VecSetRandom(Bx,rctx);
140: if (Bt) {
141: MatMultTranspose(B,Bx,s1);
142: } else {
143: MatMult(B,Bx,s1);
144: }
145: VecCopy(s1,Ax);
146: if (At) {
147: MatMultTranspose(A,Ax,s2);
148: } else {
149: MatMult(A,Ax,s2);
150: }
151: VecCopy(Bx,Cx);
152: MatMult(C,Cx,s3);
154: VecNorm(s2,NORM_INFINITY,&r2);
155: if (r2 < tol) {
156: VecNorm(s3,NORM_INFINITY,&r1);
157: } else {
158: VecAXPY(s2,none,s3);
159: VecNorm(s2,NORM_INFINITY,&r1);
160: r1 /= r2;
161: }
162: if (r1 > tol) {
163: *flg = PETSC_FALSE;
164: PetscInfo(A,"Error: %" PetscInt_FMT "-th %s %g\n",k,sop,(double)r1);
165: break;
166: }
167: }
168: PetscRandomDestroy(&rctx);
169: VecDestroy(&Ax);
170: VecDestroy(&Bx);
171: VecDestroy(&Cx);
172: VecDestroy(&s1);
173: VecDestroy(&s2);
174: VecDestroy(&s3);
175: return 0;
176: }
178: /*@
179: MatMultEqual - Compares matrix-vector products of two matrices.
181: Collective on Mat
183: Input Parameters:
184: + A - the first matrix
185: . B - the second matrix
186: - n - number of random vectors to be tested
188: Output Parameter:
189: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
191: Level: intermediate
193: @*/
194: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
195: {
196: MatMultEqual_Private(A,B,n,flg,0,PETSC_FALSE);
197: return 0;
198: }
200: /*@
201: MatMultAddEqual - Compares matrix-vector products of two matrices.
203: Collective on Mat
205: Input Parameters:
206: + A - the first matrix
207: . B - the second matrix
208: - n - number of random vectors to be tested
210: Output Parameter:
211: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
213: Level: intermediate
215: @*/
216: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
217: {
218: MatMultEqual_Private(A,B,n,flg,0,PETSC_TRUE);
219: return 0;
220: }
222: /*@
223: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
225: Collective on Mat
227: Input Parameters:
228: + A - the first matrix
229: . B - the second matrix
230: - n - number of random vectors to be tested
232: Output Parameter:
233: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
235: Level: intermediate
237: @*/
238: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
239: {
240: MatMultEqual_Private(A,B,n,flg,1,PETSC_FALSE);
241: return 0;
242: }
244: /*@
245: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
247: Collective on Mat
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: @*/
260: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
261: {
262: MatMultEqual_Private(A,B,n,flg,1,PETSC_TRUE);
263: return 0;
264: }
266: /*@
267: MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.
269: Collective on Mat
271: Input Parameters:
272: + A - the first matrix
273: . B - the second matrix
274: - n - number of random vectors to be tested
276: Output Parameter:
277: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
279: Level: intermediate
281: @*/
282: PetscErrorCode MatMultHermitianTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
283: {
284: MatMultEqual_Private(A,B,n,flg,2,PETSC_FALSE);
285: return 0;
286: }
288: /*@
289: MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.
291: Collective on Mat
293: Input Parameters:
294: + A - the first matrix
295: . B - the second matrix
296: - n - number of random vectors to be tested
298: Output Parameter:
299: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
301: Level: intermediate
303: @*/
304: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
305: {
306: MatMultEqual_Private(A,B,n,flg,2,PETSC_TRUE);
307: return 0;
308: }
310: /*@
311: MatMatMultEqual - Test A*B*x = C*x for n random vector x
313: Collective on Mat
315: Input Parameters:
316: + A - the first matrix
317: . B - the second matrix
318: . C - the third matrix
319: - n - number of random vectors to be tested
321: Output Parameter:
322: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
324: Level: intermediate
326: @*/
327: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
328: {
329: MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_FALSE);
330: return 0;
331: }
333: /*@
334: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
336: Collective on Mat
338: Input Parameters:
339: + A - the first matrix
340: . B - the second matrix
341: . C - the third matrix
342: - n - number of random vectors to be tested
344: Output Parameter:
345: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
347: Level: intermediate
349: @*/
350: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
351: {
352: MatMatMultEqual_Private(A,B,C,n,flg,PETSC_TRUE,PETSC_FALSE);
353: return 0;
354: }
356: /*@
357: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
359: Collective on Mat
361: Input Parameters:
362: + A - the first matrix
363: . B - the second matrix
364: . C - the third matrix
365: - n - number of random vectors to be tested
367: Output Parameter:
368: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
370: Level: intermediate
372: @*/
373: PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
374: {
375: MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_TRUE);
376: return 0;
377: }
379: static PetscErrorCode MatProjMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool rart,PetscBool *flg)
380: {
381: Vec x,v1,v2,v3,v4,Cx,Bx;
382: PetscReal norm_abs,norm_rel,tol=PETSC_SQRT_MACHINE_EPSILON;
383: PetscInt i,am,an,bm,bn,cm,cn;
384: PetscRandom rdm;
385: PetscScalar none = -1.0;
387: MatGetLocalSize(A,&am,&an);
388: MatGetLocalSize(B,&bm,&bn);
389: if (rart) { PetscInt t = bm; bm = bn; bn = t; }
390: MatGetLocalSize(C,&cm,&cn);
393: /* Create left vector of A: v2 */
394: MatCreateVecs(A,&Bx,&v2);
396: /* Create right vectors of B: x, v3, v4 */
397: if (rart) {
398: MatCreateVecs(B,&v1,&x);
399: } else {
400: MatCreateVecs(B,&x,&v1);
401: }
402: VecDuplicate(x,&v3);
404: MatCreateVecs(C,&Cx,&v4);
405: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
406: PetscRandomSetFromOptions(rdm);
408: *flg = PETSC_TRUE;
409: for (i=0; i<n; i++) {
410: VecSetRandom(x,rdm);
411: VecCopy(x,Cx);
412: MatMult(C,Cx,v4); /* v4 = C*x */
413: if (rart) {
414: MatMultTranspose(B,x,v1);
415: } else {
416: MatMult(B,x,v1);
417: }
418: VecCopy(v1,Bx);
419: MatMult(A,Bx,v2); /* v2 = A*B*x */
420: VecCopy(v2,v1);
421: if (rart) {
422: MatMult(B,v1,v3); /* v3 = R*A*R^t*x */
423: } else {
424: MatMultTranspose(B,v1,v3); /* v3 = Bt*A*B*x */
425: }
426: VecNorm(v4,NORM_2,&norm_abs);
427: VecAXPY(v4,none,v3);
428: VecNorm(v4,NORM_2,&norm_rel);
430: if (norm_abs > tol) norm_rel /= norm_abs;
431: if (norm_rel > tol) {
432: *flg = PETSC_FALSE;
433: PetscInfo(A,"Error: %" PetscInt_FMT "-th Mat%sMult() %g\n",i,rart ? "RARt" : "PtAP",(double)norm_rel);
434: break;
435: }
436: }
438: PetscRandomDestroy(&rdm);
439: VecDestroy(&x);
440: VecDestroy(&Bx);
441: VecDestroy(&Cx);
442: VecDestroy(&v1);
443: VecDestroy(&v2);
444: VecDestroy(&v3);
445: VecDestroy(&v4);
446: return 0;
447: }
449: /*@
450: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
452: Collective on Mat
454: Input Parameters:
455: + A - the first matrix
456: . B - the second matrix
457: . C - the third matrix
458: - n - number of random vectors to be tested
460: Output Parameter:
461: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
463: Level: intermediate
465: @*/
466: PetscErrorCode MatPtAPMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
467: {
468: MatProjMultEqual_Private(A,B,C,n,PETSC_FALSE,flg);
469: return 0;
470: }
472: /*@
473: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
475: Collective on Mat
477: Input Parameters:
478: + A - the first matrix
479: . B - the second matrix
480: . C - the third matrix
481: - n - number of random vectors to be tested
483: Output Parameter:
484: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
486: Level: intermediate
488: @*/
489: PetscErrorCode MatRARtMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
490: {
491: MatProjMultEqual_Private(A,B,C,n,PETSC_TRUE,flg);
492: return 0;
493: }
495: /*@
496: MatIsLinear - Check if a shell matrix A is a linear operator.
498: Collective on Mat
500: Input Parameters:
501: + A - the shell matrix
502: - n - number of random vectors to be tested
504: Output Parameter:
505: . flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.
507: Level: intermediate
508: @*/
509: PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool *flg)
510: {
511: Vec x,y,s1,s2;
512: PetscRandom rctx;
513: PetscScalar a;
514: PetscInt k;
515: PetscReal norm,normA;
516: MPI_Comm comm;
517: PetscMPIInt rank;
520: PetscObjectGetComm((PetscObject)A,&comm);
521: MPI_Comm_rank(comm,&rank);
523: PetscRandomCreate(comm,&rctx);
524: PetscRandomSetFromOptions(rctx);
525: MatCreateVecs(A,&x,&s1);
526: VecDuplicate(x,&y);
527: VecDuplicate(s1,&s2);
529: *flg = PETSC_TRUE;
530: for (k=0; k<n; k++) {
531: VecSetRandom(x,rctx);
532: VecSetRandom(y,rctx);
533: if (rank == 0) {
534: PetscRandomGetValue(rctx,&a);
535: }
536: MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);
538: /* s2 = a*A*x + A*y */
539: MatMult(A,y,s2); /* s2 = A*y */
540: MatMult(A,x,s1); /* s1 = A*x */
541: VecAXPY(s2,a,s1); /* s2 = a s1 + s2 */
543: /* s1 = A * (a x + y) */
544: VecAXPY(y,a,x); /* y = a x + y */
545: MatMult(A,y,s1);
546: VecNorm(s1,NORM_INFINITY,&normA);
548: VecAXPY(s2,-1.0,s1); /* s2 = - s1 + s2 */
549: VecNorm(s2,NORM_INFINITY,&norm);
550: if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
551: *flg = PETSC_FALSE;
552: 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));
553: break;
554: }
555: }
556: PetscRandomDestroy(&rctx);
557: VecDestroy(&x);
558: VecDestroy(&y);
559: VecDestroy(&s1);
560: VecDestroy(&s2);
561: return 0;
562: }