Actual source code: multequal.c
petsc-3.13.6 2020-09-29
2: #include <petsc/private/matimpl.h>
4: /*@
5: MatMultEqual - Compares matrix-vector products of two matrices.
7: Collective on Mat
9: Input Parameters:
10: + A - the first matrix
11: . B - the second matrix
12: - n - number of random vectors to be tested
14: Output Parameter:
15: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
17: Level: intermediate
19: @*/
20: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
21: {
23: Vec x,s1,s2;
24: PetscRandom rctx;
25: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
26: PetscInt am,an,bm,bn,k;
27: PetscScalar none = -1.0;
32: MatGetLocalSize(A,&am,&an);
33: MatGetLocalSize(B,&bm,&bn);
34: 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);
36: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
37: PetscRandomSetFromOptions(rctx);
38: MatCreateVecs(A,&x,&s1);
39: VecDuplicate(s1,&s2);
41: *flg = PETSC_TRUE;
42: for (k=0; k<n; k++) {
43: VecSetRandom(x,rctx);
44: MatMult(A,x,s1);
45: MatMult(B,x,s2);
46: VecNorm(s2,NORM_INFINITY,&r2);
47: if (r2 < tol) {
48: VecNorm(s1,NORM_INFINITY,&r1);
49: } else {
50: VecAXPY(s2,none,s1);
51: VecNorm(s2,NORM_INFINITY,&r1);
52: r1 /= r2;
53: }
54: if (r1 > tol) {
55: *flg = PETSC_FALSE;
56: PetscInfo2(A,"Error: %D-th MatMult() %g\n",k,(double)r1);
57: break;
58: }
59: }
60: PetscRandomDestroy(&rctx);
61: VecDestroy(&x);
62: VecDestroy(&s1);
63: VecDestroy(&s2);
64: return(0);
65: }
67: /*@
68: MatMultAddEqual - Compares matrix-vector products of two matrices.
70: Collective on Mat
72: Input Parameters:
73: + A - the first matrix
74: . B - the second matrix
75: - n - number of random vectors to be tested
77: Output Parameter:
78: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
80: Level: intermediate
82: @*/
83: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
84: {
86: Vec x,y,s1,s2;
87: PetscRandom rctx;
88: PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
89: PetscInt am,an,bm,bn,k;
90: PetscScalar none = -1.0;
93: MatGetLocalSize(A,&am,&an);
94: MatGetLocalSize(B,&bm,&bn);
95: 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);
97: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
98: PetscRandomSetFromOptions(rctx);
99: MatCreateVecs(A,&x,&s1);
100: VecDuplicate(s1,&s2);
101: VecDuplicate(s1,&y);
103: *flg = PETSC_TRUE;
104: for (k=0; k<n; k++) {
105: VecSetRandom(x,rctx);
106: VecSetRandom(y,rctx);
107: MatMultAdd(A,x,y,s1);
108: MatMultAdd(B,x,y,s2);
109: VecNorm(s2,NORM_INFINITY,&r2);
110: if (r2 < tol) {
111: VecNorm(s1,NORM_INFINITY,&r1);
112: } else {
113: VecAXPY(s2,none,s1);
114: VecNorm(s2,NORM_INFINITY,&r1);
115: r1 /= r2;
116: }
117: if (r1 > tol) {
118: *flg = PETSC_FALSE;
119: PetscInfo2(A,"Error: %d-th MatMultAdd() %g\n",k,(double)r1);
120: break;
121: }
122: }
123: PetscRandomDestroy(&rctx);
124: VecDestroy(&x);
125: VecDestroy(&y);
126: VecDestroy(&s1);
127: VecDestroy(&s2);
128: return(0);
129: }
131: /*@
132: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
134: Collective on Mat
136: Input Parameters:
137: + A - the first matrix
138: . B - the second matrix
139: - n - number of random vectors to be tested
141: Output Parameter:
142: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
144: Level: intermediate
146: @*/
147: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
148: {
150: Vec x,s1,s2;
151: PetscRandom rctx;
152: PetscReal r1,r2,tol= PETSC_SQRT_MACHINE_EPSILON;
153: PetscInt am,an,bm,bn,k;
154: PetscScalar none = -1.0;
157: MatGetLocalSize(A,&am,&an);
158: MatGetLocalSize(B,&bm,&bn);
159: 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);
161: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
162: PetscRandomSetFromOptions(rctx);
163: MatCreateVecs(A,&s1,&x);
164: VecDuplicate(s1,&s2);
166: *flg = PETSC_TRUE;
167: for (k=0; k<n; k++) {
168: VecSetRandom(x,rctx);
169: MatMultTranspose(A,x,s1);
170: MatMultTranspose(B,x,s2);
171: VecNorm(s2,NORM_INFINITY,&r2);
172: if (r2 < tol) {
173: VecNorm(s1,NORM_INFINITY,&r1);
174: } else {
175: VecAXPY(s2,none,s1);
176: VecNorm(s2,NORM_INFINITY,&r1);
177: r1 /= r2;
178: }
179: if (r1 > tol) {
180: *flg = PETSC_FALSE;
181: PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);
182: break;
183: }
184: }
185: PetscRandomDestroy(&rctx);
186: VecDestroy(&x);
187: VecDestroy(&s1);
188: VecDestroy(&s2);
189: return(0);
190: }
192: /*@
193: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
195: Collective on Mat
197: Input Parameters:
198: + A - the first matrix
199: . B - the second matrix
200: - n - number of random vectors to be tested
202: Output Parameter:
203: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
205: Level: intermediate
207: @*/
208: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
209: {
211: Vec x,y,s1,s2;
212: PetscRandom rctx;
213: PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
214: PetscInt am,an,bm,bn,k;
215: PetscScalar none = -1.0;
218: MatGetLocalSize(A,&am,&an);
219: MatGetLocalSize(B,&bm,&bn);
220: 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);
222: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
223: PetscRandomSetFromOptions(rctx);
224: MatCreateVecs(A,&s1,&x);
225: VecDuplicate(s1,&s2);
226: VecDuplicate(s1,&y);
228: *flg = PETSC_TRUE;
229: for (k=0; k<n; k++) {
230: VecSetRandom(x,rctx);
231: VecSetRandom(y,rctx);
232: MatMultTransposeAdd(A,x,y,s1);
233: MatMultTransposeAdd(B,x,y,s2);
234: VecNorm(s2,NORM_INFINITY,&r2);
235: if (r2 < tol) {
236: VecNorm(s1,NORM_INFINITY,&r1);
237: } else {
238: VecAXPY(s2,none,s1);
239: VecNorm(s2,NORM_INFINITY,&r1);
240: r1 /= r2;
241: }
242: if (r1 > tol) {
243: *flg = PETSC_FALSE;
244: PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);
245: break;
246: }
247: }
248: PetscRandomDestroy(&rctx);
249: VecDestroy(&x);
250: VecDestroy(&y);
251: VecDestroy(&s1);
252: VecDestroy(&s2);
253: return(0);
254: }
256: /*@
257: MatMatMultEqual - Test A*B*x = C*x for n random vector x
259: Collective on Mat
261: Input Parameters:
262: + A - the first matrix
263: . B - the second matrix
264: . C - the third matrix
265: - n - number of random vectors to be tested
267: Output Parameter:
268: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
270: Level: intermediate
272: @*/
273: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
274: {
276: Vec x,s1,s2,s3;
277: PetscRandom rctx;
278: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
279: PetscInt am,an,bm,bn,cm,cn,k;
280: PetscScalar none = -1.0;
283: MatGetLocalSize(A,&am,&an);
284: MatGetLocalSize(B,&bm,&bn);
285: MatGetLocalSize(C,&cm,&cn);
286: 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);
288: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
289: PetscRandomSetFromOptions(rctx);
290: MatCreateVecs(B,&x,&s1);
291: MatCreateVecs(C,NULL,&s2);
292: VecDuplicate(s2,&s3);
294: *flg = PETSC_TRUE;
295: for (k=0; k<n; k++) {
296: VecSetRandom(x,rctx);
297: MatMult(B,x,s1);
298: MatMult(A,s1,s2);
299: MatMult(C,x,s3);
301: VecNorm(s2,NORM_INFINITY,&r2);
302: if (r2 < tol) {
303: VecNorm(s3,NORM_INFINITY,&r1);
304: } else {
305: VecAXPY(s2,none,s3);
306: VecNorm(s2,NORM_INFINITY,&r1);
307: r1 /= r2;
308: }
309: if (r1 > tol) {
310: *flg = PETSC_FALSE;
311: PetscInfo2(A,"Error: %D-th MatMatMult() %g\n",k,(double)r1);
312: break;
313: }
314: }
315: PetscRandomDestroy(&rctx);
316: VecDestroy(&x);
317: VecDestroy(&s1);
318: VecDestroy(&s2);
319: VecDestroy(&s3);
320: return(0);
321: }
323: /*@
324: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
326: Collective on Mat
328: Input Parameters:
329: + A - the first matrix
330: . B - the second matrix
331: . C - the third matrix
332: - n - number of random vectors to be tested
334: Output Parameter:
335: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
337: Level: intermediate
339: @*/
340: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
341: {
343: Vec x,s1,s2,s3;
344: PetscRandom rctx;
345: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
346: PetscInt am,an,bm,bn,cm,cn,k;
347: PetscScalar none = -1.0;
350: MatGetLocalSize(A,&am,&an);
351: MatGetLocalSize(B,&bm,&bn);
352: MatGetLocalSize(C,&cm,&cn);
353: if (am != bm || an != 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);
355: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
356: PetscRandomSetFromOptions(rctx);
357: MatCreateVecs(B,&x,&s1);
358: MatCreateVecs(C,NULL,&s2);
359: VecDuplicate(s2,&s3);
361: *flg = PETSC_TRUE;
362: for (k=0; k<n; k++) {
363: VecSetRandom(x,rctx);
364: MatMult(B,x,s1);
365: MatMultTranspose(A,s1,s2);
366: MatMult(C,x,s3);
368: VecNorm(s2,NORM_INFINITY,&r2);
369: if (r2 < tol) {
370: VecNorm(s3,NORM_INFINITY,&r1);
371: } else {
372: VecAXPY(s2,none,s3);
373: VecNorm(s2,NORM_INFINITY,&r1);
374: r1 /= r2;
375: }
376: if (r1 > tol) {
377: *flg = PETSC_FALSE;
378: PetscInfo2(A,"Error: %D-th MatTransposeMatMult() %g\n",k,(double)r1);
379: break;
380: }
381: }
382: PetscRandomDestroy(&rctx);
383: VecDestroy(&x);
384: VecDestroy(&s1);
385: VecDestroy(&s2);
386: VecDestroy(&s3);
387: return(0);
388: }
390: /*@
391: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
393: Collective on Mat
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: @*/
407: PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
408: {
410: Vec x,s1,s2,s3;
411: PetscRandom rctx;
412: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
413: PetscInt am,an,bm,bn,cm,cn,k;
414: PetscScalar none = -1.0;
417: MatGetLocalSize(A,&am,&an);
418: MatGetLocalSize(B,&bm,&bn);
419: MatGetLocalSize(C,&cm,&cn);
420: if (an != bn || am != cm || bm != 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);
422: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
423: PetscRandomSetFromOptions(rctx);
424: MatCreateVecs(B,&s1,&x);
425: MatCreateVecs(C,NULL,&s2);
426: VecDuplicate(s2,&s3);
428: *flg = PETSC_TRUE;
429: for (k=0; k<n; k++) {
430: VecSetRandom(x,rctx);
431: MatMultTranspose(B,x,s1);
432: MatMult(A,s1,s2);
433: MatMult(C,x,s3);
435: VecNorm(s2,NORM_INFINITY,&r2);
436: if (r2 < tol) {
437: VecNorm(s3,NORM_INFINITY,&r1);
438: } else {
439: VecAXPY(s2,none,s3);
440: VecNorm(s2,NORM_INFINITY,&r1);
441: r1 /= r2;
442: }
443: if (r1 > tol) {
444: *flg = PETSC_FALSE;
445: PetscInfo2(A,"Error: %D-th MatMatTransposeMult() %g\n",k,(double)r1);
446: break;
447: }
448: }
449: PetscRandomDestroy(&rctx);
450: VecDestroy(&x);
451: VecDestroy(&s1);
452: VecDestroy(&s2);
453: VecDestroy(&s3);
454: return(0);
455: }
457: /*@
458: MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
460: Collective on Mat
462: Input Parameters:
463: + A - the first matrix
464: . B - the second matrix
465: . C - the third matrix
466: - n - number of random vectors to be tested
468: Output Parameter:
469: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
471: Level: intermediate
473: @*/
474: PetscErrorCode MatPtAPMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
475: {
477: Vec x,v1,v2,v3,v4;
478: PetscReal norm_abs,norm_rel,tol=PETSC_SQRT_MACHINE_EPSILON;
479: PetscInt i,am,an,bm,bn,cm,cn;
480: PetscRandom rdm;
481: PetscScalar none = -1.0;
484: MatGetLocalSize(A,&am,&an);
485: MatGetLocalSize(B,&bm,&bn);
486: MatGetLocalSize(C,&cm,&cn);
487: 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);
489: /* Create left vector of A: v2 */
490: MatCreateVecs(A,NULL,&v2);
492: /* Create right vectors of B: x, v3, v4 */
493: MatCreateVecs(B,&x,&v1);
494: VecDuplicate(x,&v3);
495: VecDuplicate(x,&v4);
497: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
498: PetscRandomSetFromOptions(rdm);
500: *flg = PETSC_TRUE;
501: for (i=0; i<n; i++) {
502: VecSetRandom(x,rdm);
503: MatMult(B,x,v1);
504: MatMult(A,v1,v2); /* v2 = A*B*x */
506: MatMultTranspose(B,v2,v3); /* v3 = Bt*A*B*x */
507: MatMult(C,x,v4); /* v4 = C*x */
508: VecNorm(v4,NORM_2,&norm_abs);
509: VecAXPY(v4,none,v3);
510: VecNorm(v4,NORM_2,&norm_rel);
512: if (norm_abs > tol) norm_rel /= norm_abs;
513: if (norm_rel > tol) {
514: *flg = PETSC_FALSE;
515: PetscInfo2(A,"Error: %D-th MatPtAPMult() %g\n",i,(double)norm_rel);
516: break;
517: }
518: }
520: PetscRandomDestroy(&rdm);
521: VecDestroy(&x);
522: VecDestroy(&v1);
523: VecDestroy(&v2);
524: VecDestroy(&v3);
525: VecDestroy(&v4);
526: return(0);
527: }
529: /*@
530: MatIsLinear - Check if a shell matrix A is a linear operator.
532: Collective on Mat
534: Input Parameters:
535: + A - the shell matrix
536: - n - number of random vectors to be tested
538: Output Parameter:
539: . flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.
541: Level: intermediate
542: @*/
543: PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool *flg)
544: {
546: Vec x,y,s1,s2;
547: PetscRandom rctx;
548: PetscScalar a;
549: PetscInt k;
550: PetscReal norm,normA;
551: MPI_Comm comm;
552: PetscMPIInt rank;
556: PetscObjectGetComm((PetscObject)A,&comm);
557: MPI_Comm_rank(comm,&rank);
559: PetscRandomCreate(comm,&rctx);
560: PetscRandomSetFromOptions(rctx);
561: MatCreateVecs(A,&x,&s1);
562: VecDuplicate(x,&y);
563: VecDuplicate(s1,&s2);
565: *flg = PETSC_TRUE;
566: for (k=0; k<n; k++) {
567: VecSetRandom(x,rctx);
568: VecSetRandom(y,rctx);
569: if (!rank) {
570: PetscRandomGetValue(rctx,&a);
571: }
572: MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);
574: /* s2 = a*A*x + A*y */
575: MatMult(A,y,s2); /* s2 = A*y */
576: MatMult(A,x,s1); /* s1 = A*x */
577: VecAXPY(s2,a,s1); /* s2 = a s1 + s2 */
579: /* s1 = A * (a x + y) */
580: VecAXPY(y,a,x); /* y = a x + y */
581: MatMult(A,y,s1);
582: VecNorm(s1,NORM_INFINITY,&normA);
584: VecAXPY(s2,-1.0,s1); /* s2 = - s1 + s2 */
585: VecNorm(s2,NORM_INFINITY,&norm);
586: if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
587: *flg = PETSC_FALSE;
588: 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);
589: break;
590: }
591: }
592: PetscRandomDestroy(&rctx);
593: VecDestroy(&x);
594: VecDestroy(&y);
595: VecDestroy(&s1);
596: VecDestroy(&s2);
597: return(0);
598: }