Actual source code: multequal.c
2: #include <petsc/private/matimpl.h>
4: static PetscErrorCode MatMultEqual_Private(Mat A,Mat B,PetscInt n,PetscBool *flg,PetscBool t,PetscBool add)
5: {
7: Vec Ax = NULL,Bx = NULL,s1 = NULL,s2 = NULL,Ay = NULL, By = NULL;
8: PetscRandom rctx;
9: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
10: PetscInt am,an,bm,bn,k;
11: PetscScalar none = -1.0;
12: const char* sops[] = {"MatMult","MatMultAdd","MatMultTranspose","MatMultTranposeAdd"};
13: const char* sop;
23: MatGetLocalSize(A,&am,&an);
24: MatGetLocalSize(B,&bm,&bn);
25: if (am != bm || an != bn) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
26: sop = sops[(add ? 1 : 0) + 2 * (t ? 1 : 0)];
27: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
28: PetscRandomSetFromOptions(rctx);
29: if (t) {
30: MatCreateVecs(A,&s1,&Ax);
31: MatCreateVecs(B,&s2,&Bx);
32: } else {
33: MatCreateVecs(A,&Ax,&s1);
34: MatCreateVecs(B,&Bx,&s2);
35: }
36: if (add) {
37: VecDuplicate(s1,&Ay);
38: VecDuplicate(s2,&By);
39: }
41: *flg = PETSC_TRUE;
42: for (k=0; k<n; k++) {
43: VecSetRandom(Ax,rctx);
44: VecCopy(Ax,Bx);
45: if (add) {
46: VecSetRandom(Ay,rctx);
47: VecCopy(Ay,By);
48: }
49: if (t) {
50: if (add) {
51: MatMultTransposeAdd(A,Ax,Ay,s1);
52: MatMultTransposeAdd(B,Bx,By,s2);
53: } else {
54: MatMultTranspose(A,Ax,s1);
55: MatMultTranspose(B,Bx,s2);
56: }
57: } else {
58: if (add) {
59: MatMultAdd(A,Ax,Ay,s1);
60: MatMultAdd(B,Bx,By,s2);
61: } else {
62: MatMult(A,Ax,s1);
63: MatMult(B,Bx,s2);
64: }
65: }
66: VecNorm(s2,NORM_INFINITY,&r2);
67: if (r2 < tol) {
68: VecNorm(s1,NORM_INFINITY,&r1);
69: } else {
70: VecAXPY(s2,none,s1);
71: VecNorm(s2,NORM_INFINITY,&r1);
72: r1 /= r2;
73: }
74: if (r1 > tol) {
75: *flg = PETSC_FALSE;
76: PetscInfo3(A,"Error: %D-th %s() %g\n",k,sop,(double)r1);
77: break;
78: }
79: }
80: PetscRandomDestroy(&rctx);
81: VecDestroy(&Ax);
82: VecDestroy(&Bx);
83: VecDestroy(&Ay);
84: VecDestroy(&By);
85: VecDestroy(&s1);
86: VecDestroy(&s2);
87: return(0);
88: }
90: static PetscErrorCode MatMatMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg,PetscBool At,PetscBool Bt)
91: {
93: Vec Ax,Bx,Cx,s1,s2,s3;
94: PetscRandom rctx;
95: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
96: PetscInt am,an,bm,bn,cm,cn,k;
97: PetscScalar none = -1.0;
98: const char* sops[] = {"MatMatMult","MatTransposeMatMult","MatMatTransposeMult","MatTransposeMatTranposeMult"};
99: const char* sop;
111: MatGetLocalSize(A,&am,&an);
112: MatGetLocalSize(B,&bm,&bn);
113: MatGetLocalSize(C,&cm,&cn);
114: if (At) { PetscInt tt = an; an = am; am = tt; };
115: if (Bt) { PetscInt tt = bn; bn = bm; bm = tt; };
116: if (an != bm || am != cm || bn != cn) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %D %D %D %D",am,an,bm,bn,cm,cn);
118: sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
119: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
120: PetscRandomSetFromOptions(rctx);
121: if (Bt) {
122: MatCreateVecs(B,&s1,&Bx);
123: } else {
124: MatCreateVecs(B,&Bx,&s1);
125: }
126: if (At) {
127: MatCreateVecs(A,&s2,&Ax);
128: } else {
129: MatCreateVecs(A,&Ax,&s2);
130: }
131: MatCreateVecs(C,&Cx,&s3);
133: *flg = PETSC_TRUE;
134: for (k=0; k<n; k++) {
135: VecSetRandom(Bx,rctx);
136: if (Bt) {
137: MatMultTranspose(B,Bx,s1);
138: } else {
139: MatMult(B,Bx,s1);
140: }
141: VecCopy(s1,Ax);
142: if (At) {
143: MatMultTranspose(A,Ax,s2);
144: } else {
145: MatMult(A,Ax,s2);
146: }
147: VecCopy(Bx,Cx);
148: MatMult(C,Cx,s3);
150: VecNorm(s2,NORM_INFINITY,&r2);
151: if (r2 < tol) {
152: VecNorm(s3,NORM_INFINITY,&r1);
153: } else {
154: VecAXPY(s2,none,s3);
155: VecNorm(s2,NORM_INFINITY,&r1);
156: r1 /= r2;
157: }
158: if (r1 > tol) {
159: *flg = PETSC_FALSE;
160: PetscInfo3(A,"Error: %D-th %s %g\n",k,sop,(double)r1);
161: break;
162: }
163: }
164: PetscRandomDestroy(&rctx);
165: VecDestroy(&Ax);
166: VecDestroy(&Bx);
167: VecDestroy(&Cx);
168: VecDestroy(&s1);
169: VecDestroy(&s2);
170: VecDestroy(&s3);
171: return(0);
172: }
174: /*@
175: MatMultEqual - Compares matrix-vector products of two matrices.
177: Collective on Mat
179: Input Parameters:
180: + A - the first matrix
181: . B - the second matrix
182: - n - number of random vectors to be tested
184: Output Parameter:
185: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
187: Level: intermediate
189: @*/
190: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
191: {
195: MatMultEqual_Private(A,B,n,flg,PETSC_FALSE,PETSC_FALSE);
196: return(0);
197: }
199: /*@
200: MatMultAddEqual - Compares matrix-vector products of two matrices.
202: Collective on Mat
204: Input Parameters:
205: + A - the first matrix
206: . B - the second matrix
207: - n - number of random vectors to be tested
209: Output Parameter:
210: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
212: Level: intermediate
214: @*/
215: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
216: {
220: MatMultEqual_Private(A,B,n,flg,PETSC_FALSE,PETSC_TRUE);
221: return(0);
222: }
224: /*@
225: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
227: Collective on Mat
229: Input Parameters:
230: + A - the first matrix
231: . B - the second matrix
232: - n - number of random vectors to be tested
234: Output Parameter:
235: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
237: Level: intermediate
239: @*/
240: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
241: {
245: MatMultEqual_Private(A,B,n,flg,PETSC_TRUE,PETSC_FALSE);
246: return(0);
247: }
249: /*@
250: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
252: Collective on Mat
254: Input Parameters:
255: + A - the first matrix
256: . B - the second matrix
257: - n - number of random vectors to be tested
259: Output Parameter:
260: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
262: Level: intermediate
264: @*/
265: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
266: {
270: MatMultEqual_Private(A,B,n,flg,PETSC_TRUE,PETSC_TRUE);
271: return(0);
272: }
274: /*@
275: MatMatMultEqual - Test A*B*x = C*x for n random vector x
277: Collective on Mat
279: Input Parameters:
280: + A - the first matrix
281: . B - the second matrix
282: . C - the third matrix
283: - n - number of random vectors to be tested
285: Output Parameter:
286: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
288: Level: intermediate
290: @*/
291: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
292: {
296: MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_FALSE);
297: return(0);
298: }
300: /*@
301: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
303: Collective on Mat
305: Input Parameters:
306: + A - the first matrix
307: . B - the second matrix
308: . C - the third matrix
309: - n - number of random vectors to be tested
311: Output Parameter:
312: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
314: Level: intermediate
316: @*/
317: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
318: {
322: MatMatMultEqual_Private(A,B,C,n,flg,PETSC_TRUE,PETSC_FALSE);
323: return(0);
324: }
326: /*@
327: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
329: Collective on Mat
331: Input Parameters:
332: + A - the first matrix
333: . B - the second matrix
334: . C - the third matrix
335: - n - number of random vectors to be tested
337: Output Parameter:
338: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
340: Level: intermediate
342: @*/
343: PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
344: {
348: MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_TRUE);
349: return(0);
350: }
352: static PetscErrorCode MatProjMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool rart,PetscBool *flg)
353: {
355: Vec x,v1,v2,v3,v4,Cx,Bx;
356: PetscReal norm_abs,norm_rel,tol=PETSC_SQRT_MACHINE_EPSILON;
357: PetscInt i,am,an,bm,bn,cm,cn;
358: PetscRandom rdm;
359: PetscScalar none = -1.0;
362: MatGetLocalSize(A,&am,&an);
363: MatGetLocalSize(B,&bm,&bn);
364: if (rart) { PetscInt t = bm; bm = bn; bn = t; }
365: MatGetLocalSize(C,&cm,&cn);
366: if (an != bm || bn != cm || bn != cn) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %D %D %D %D %D %D",am,an,bm,bn,cm,cn);
368: /* Create left vector of A: v2 */
369: MatCreateVecs(A,&Bx,&v2);
371: /* Create right vectors of B: x, v3, v4 */
372: if (rart) {
373: MatCreateVecs(B,&v1,&x);
374: } else {
375: MatCreateVecs(B,&x,&v1);
376: }
377: VecDuplicate(x,&v3);
379: MatCreateVecs(C,&Cx,&v4);
380: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
381: PetscRandomSetFromOptions(rdm);
383: *flg = PETSC_TRUE;
384: for (i=0; i<n; i++) {
385: VecSetRandom(x,rdm);
386: VecCopy(x,Cx);
387: MatMult(C,Cx,v4); /* v4 = C*x */
388: if (rart) {
389: MatMultTranspose(B,x,v1);
390: } else {
391: MatMult(B,x,v1);
392: }
393: VecCopy(v1,Bx);
394: MatMult(A,Bx,v2); /* v2 = A*B*x */
395: VecCopy(v2,v1);
396: if (rart) {
397: MatMult(B,v1,v3); /* v3 = R*A*R^t*x */
398: } else {
399: MatMultTranspose(B,v1,v3); /* v3 = Bt*A*B*x */
400: }
401: VecNorm(v4,NORM_2,&norm_abs);
402: VecAXPY(v4,none,v3);
403: VecNorm(v4,NORM_2,&norm_rel);
405: if (norm_abs > tol) norm_rel /= norm_abs;
406: if (norm_rel > tol) {
407: *flg = PETSC_FALSE;
408: PetscInfo3(A,"Error: %D-th Mat%sMult() %g\n",i,rart ? "RARt" : "PtAP",(double)norm_rel);
409: break;
410: }
411: }
413: PetscRandomDestroy(&rdm);
414: VecDestroy(&x);
415: VecDestroy(&Bx);
416: VecDestroy(&Cx);
417: VecDestroy(&v1);
418: VecDestroy(&v2);
419: VecDestroy(&v3);
420: VecDestroy(&v4);
421: return(0);
422: }
424: /*@
425: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
427: Collective on Mat
429: Input Parameters:
430: + A - the first matrix
431: . B - the second matrix
432: . C - the third matrix
433: - n - number of random vectors to be tested
435: Output Parameter:
436: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
438: Level: intermediate
440: @*/
441: PetscErrorCode MatPtAPMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
442: {
446: MatProjMultEqual_Private(A,B,C,n,PETSC_FALSE,flg);
447: return(0);
448: }
450: /*@
451: MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
453: Collective on Mat
455: Input Parameters:
456: + A - the first matrix
457: . B - the second matrix
458: . C - the third matrix
459: - n - number of random vectors to be tested
461: Output Parameter:
462: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
464: Level: intermediate
466: @*/
467: PetscErrorCode MatRARtMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
468: {
472: MatProjMultEqual_Private(A,B,C,n,PETSC_TRUE,flg);
473: return(0);
474: }
476: /*@
477: MatIsLinear - Check if a shell matrix A is a linear operator.
479: Collective on Mat
481: Input Parameters:
482: + A - the shell matrix
483: - n - number of random vectors to be tested
485: Output Parameter:
486: . flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.
488: Level: intermediate
489: @*/
490: PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool *flg)
491: {
493: Vec x,y,s1,s2;
494: PetscRandom rctx;
495: PetscScalar a;
496: PetscInt k;
497: PetscReal norm,normA;
498: MPI_Comm comm;
499: PetscMPIInt rank;
503: PetscObjectGetComm((PetscObject)A,&comm);
504: MPI_Comm_rank(comm,&rank);
506: PetscRandomCreate(comm,&rctx);
507: PetscRandomSetFromOptions(rctx);
508: MatCreateVecs(A,&x,&s1);
509: VecDuplicate(x,&y);
510: VecDuplicate(s1,&s2);
512: *flg = PETSC_TRUE;
513: for (k=0; k<n; k++) {
514: VecSetRandom(x,rctx);
515: VecSetRandom(y,rctx);
516: if (!rank) {
517: PetscRandomGetValue(rctx,&a);
518: }
519: MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);
521: /* s2 = a*A*x + A*y */
522: MatMult(A,y,s2); /* s2 = A*y */
523: MatMult(A,x,s1); /* s1 = A*x */
524: VecAXPY(s2,a,s1); /* s2 = a s1 + s2 */
526: /* s1 = A * (a x + y) */
527: VecAXPY(y,a,x); /* y = a x + y */
528: MatMult(A,y,s1);
529: VecNorm(s1,NORM_INFINITY,&normA);
531: VecAXPY(s2,-1.0,s1); /* s2 = - s1 + s2 */
532: VecNorm(s2,NORM_INFINITY,&norm);
533: if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
534: *flg = PETSC_FALSE;
535: PetscInfo3(A,"Error: %D-th |A*(ax+y) - (a*A*x+A*y)|/|A(ax+y)| %g > tol %g\n",k,(double)norm/normA,100.*PETSC_MACHINE_EPSILON);
536: break;
537: }
538: }
539: PetscRandomDestroy(&rctx);
540: VecDestroy(&x);
541: VecDestroy(&y);
542: VecDestroy(&s1);
543: VecDestroy(&s2);
544: return(0);
545: }