Actual source code: multequal.c
petsc-3.11.4 2019-09-28
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: Concepts: matrices^equality between
20: @*/
21: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
22: {
24: Vec x,s1,s2;
25: PetscRandom rctx;
26: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
27: PetscInt am,an,bm,bn,k;
28: PetscScalar none = -1.0;
33: MatGetLocalSize(A,&am,&an);
34: MatGetLocalSize(B,&bm,&bn);
35: 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);
37: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
38: PetscRandomSetFromOptions(rctx);
39: MatCreateVecs(A,&x,&s1);
40: VecDuplicate(s1,&s2);
42: *flg = PETSC_TRUE;
43: for (k=0; k<n; k++) {
44: VecSetRandom(x,rctx);
45: MatMult(A,x,s1);
46: MatMult(B,x,s2);
47: VecNorm(s2,NORM_INFINITY,&r2);
48: if (r2 < tol) {
49: VecNorm(s1,NORM_INFINITY,&r1);
50: } else {
51: VecAXPY(s2,none,s1);
52: VecNorm(s2,NORM_INFINITY,&r1);
53: r1 /= r2;
54: }
55: if (r1 > tol) {
56: *flg = PETSC_FALSE;
57: PetscInfo2(A,"Error: %D-th MatMult() %g\n",k,(double)r1);
58: break;
59: }
60: }
61: PetscRandomDestroy(&rctx);
62: VecDestroy(&x);
63: VecDestroy(&s1);
64: VecDestroy(&s2);
65: return(0);
66: }
68: /*@
69: MatMultAddEqual - Compares matrix-vector products of two matrices.
71: Collective on Mat
73: Input Parameters:
74: + A - the first matrix
75: - B - the second matrix
76: - n - number of random vectors to be tested
78: Output Parameter:
79: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
81: Level: intermediate
83: Concepts: matrices^equality between
84: @*/
85: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
86: {
88: Vec x,y,s1,s2;
89: PetscRandom rctx;
90: PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
91: PetscInt am,an,bm,bn,k;
92: PetscScalar none = -1.0;
95: MatGetLocalSize(A,&am,&an);
96: MatGetLocalSize(B,&bm,&bn);
97: 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);
99: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
100: PetscRandomSetFromOptions(rctx);
101: MatCreateVecs(A,&x,&s1);
102: VecDuplicate(s1,&s2);
103: VecDuplicate(s1,&y);
105: *flg = PETSC_TRUE;
106: for (k=0; k<n; k++) {
107: VecSetRandom(x,rctx);
108: VecSetRandom(y,rctx);
109: MatMultAdd(A,x,y,s1);
110: MatMultAdd(B,x,y,s2);
111: VecNorm(s2,NORM_INFINITY,&r2);
112: if (r2 < tol) {
113: VecNorm(s1,NORM_INFINITY,&r1);
114: } else {
115: VecAXPY(s2,none,s1);
116: VecNorm(s2,NORM_INFINITY,&r1);
117: r1 /= r2;
118: }
119: if (r1 > tol) {
120: *flg = PETSC_FALSE;
121: PetscInfo2(A,"Error: %d-th MatMultAdd() %g\n",k,(double)r1);
122: break;
123: }
124: }
125: PetscRandomDestroy(&rctx);
126: VecDestroy(&x);
127: VecDestroy(&y);
128: VecDestroy(&s1);
129: VecDestroy(&s2);
130: return(0);
131: }
133: /*@
134: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
136: Collective on Mat
138: Input Parameters:
139: + A - the first matrix
140: - B - the second matrix
141: - n - number of random vectors to be tested
143: Output Parameter:
144: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
146: Level: intermediate
148: Concepts: matrices^equality between
149: @*/
150: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
151: {
153: Vec x,s1,s2;
154: PetscRandom rctx;
155: PetscReal r1,r2,tol= PETSC_SQRT_MACHINE_EPSILON;
156: PetscInt am,an,bm,bn,k;
157: PetscScalar none = -1.0;
160: MatGetLocalSize(A,&am,&an);
161: MatGetLocalSize(B,&bm,&bn);
162: 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);
164: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
165: PetscRandomSetFromOptions(rctx);
166: MatCreateVecs(A,&s1,&x);
167: VecDuplicate(s1,&s2);
169: *flg = PETSC_TRUE;
170: for (k=0; k<n; k++) {
171: VecSetRandom(x,rctx);
172: MatMultTranspose(A,x,s1);
173: MatMultTranspose(B,x,s2);
174: VecNorm(s2,NORM_INFINITY,&r2);
175: if (r2 < tol) {
176: VecNorm(s1,NORM_INFINITY,&r1);
177: } else {
178: VecAXPY(s2,none,s1);
179: VecNorm(s2,NORM_INFINITY,&r1);
180: r1 /= r2;
181: }
182: if (r1 > tol) {
183: *flg = PETSC_FALSE;
184: PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);
185: break;
186: }
187: }
188: PetscRandomDestroy(&rctx);
189: VecDestroy(&x);
190: VecDestroy(&s1);
191: VecDestroy(&s2);
192: return(0);
193: }
195: /*@
196: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
198: Collective on Mat
200: Input Parameters:
201: + A - the first matrix
202: - B - the second matrix
203: - n - number of random vectors to be tested
205: Output Parameter:
206: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
208: Level: intermediate
210: Concepts: matrices^equality between
211: @*/
212: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
213: {
215: Vec x,y,s1,s2;
216: PetscRandom rctx;
217: PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
218: PetscInt am,an,bm,bn,k;
219: PetscScalar none = -1.0;
222: MatGetLocalSize(A,&am,&an);
223: MatGetLocalSize(B,&bm,&bn);
224: 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);
226: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
227: PetscRandomSetFromOptions(rctx);
228: MatCreateVecs(A,&s1,&x);
229: VecDuplicate(s1,&s2);
230: VecDuplicate(s1,&y);
232: *flg = PETSC_TRUE;
233: for (k=0; k<n; k++) {
234: VecSetRandom(x,rctx);
235: VecSetRandom(y,rctx);
236: MatMultTransposeAdd(A,x,y,s1);
237: MatMultTransposeAdd(B,x,y,s2);
238: VecNorm(s2,NORM_INFINITY,&r2);
239: if (r2 < tol) {
240: VecNorm(s1,NORM_INFINITY,&r1);
241: } else {
242: VecAXPY(s2,none,s1);
243: VecNorm(s2,NORM_INFINITY,&r1);
244: r1 /= r2;
245: }
246: if (r1 > tol) {
247: *flg = PETSC_FALSE;
248: PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);
249: break;
250: }
251: }
252: PetscRandomDestroy(&rctx);
253: VecDestroy(&x);
254: VecDestroy(&y);
255: VecDestroy(&s1);
256: VecDestroy(&s2);
257: return(0);
258: }
260: /*@
261: MatMatMultEqual - Test A*B*x = C*x for n random vector x
263: Collective on Mat
265: Input Parameters:
266: + A - the first matrix
267: - B - the second matrix
268: - C - the third matrix
269: - n - number of random vectors to be tested
271: Output Parameter:
272: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
274: Level: intermediate
276: Concepts: matrices^equality between
277: @*/
278: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
279: {
281: Vec x,s1,s2,s3;
282: PetscRandom rctx;
283: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
284: PetscInt am,an,bm,bn,cm,cn,k;
285: PetscScalar none = -1.0;
288: MatGetLocalSize(A,&am,&an);
289: MatGetLocalSize(B,&bm,&bn);
290: MatGetLocalSize(C,&cm,&cn);
291: 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);
293: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
294: PetscRandomSetFromOptions(rctx);
295: MatCreateVecs(B,&x,&s1);
296: MatCreateVecs(C,NULL,&s2);
297: VecDuplicate(s2,&s3);
299: *flg = PETSC_TRUE;
300: for (k=0; k<n; k++) {
301: VecSetRandom(x,rctx);
302: MatMult(B,x,s1);
303: MatMult(A,s1,s2);
304: MatMult(C,x,s3);
306: VecNorm(s2,NORM_INFINITY,&r2);
307: if (r2 < tol) {
308: VecNorm(s3,NORM_INFINITY,&r1);
309: } else {
310: VecAXPY(s2,none,s3);
311: VecNorm(s2,NORM_INFINITY,&r1);
312: r1 /= r2;
313: }
314: if (r1 > tol) {
315: *flg = PETSC_FALSE;
316: PetscInfo2(A,"Error: %D-th MatMatMult() %g\n",k,(double)r1);
317: break;
318: }
319: }
320: PetscRandomDestroy(&rctx);
321: VecDestroy(&x);
322: VecDestroy(&s1);
323: VecDestroy(&s2);
324: VecDestroy(&s3);
325: return(0);
326: }
328: /*@
329: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
331: Collective on Mat
333: Input Parameters:
334: + A - the first matrix
335: - B - the second matrix
336: - C - the third matrix
337: - n - number of random vectors to be tested
339: Output Parameter:
340: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
342: Level: intermediate
344: Concepts: matrices^equality between
345: @*/
346: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
347: {
349: Vec x,s1,s2,s3;
350: PetscRandom rctx;
351: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
352: PetscInt am,an,bm,bn,cm,cn,k;
353: PetscScalar none = -1.0;
356: MatGetLocalSize(A,&am,&an);
357: MatGetLocalSize(B,&bm,&bn);
358: MatGetLocalSize(C,&cm,&cn);
359: 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);
361: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
362: PetscRandomSetFromOptions(rctx);
363: MatCreateVecs(B,&x,&s1);
364: MatCreateVecs(C,NULL,&s2);
365: VecDuplicate(s2,&s3);
367: *flg = PETSC_TRUE;
368: for (k=0; k<n; k++) {
369: VecSetRandom(x,rctx);
370: MatMult(B,x,s1);
371: MatMultTranspose(A,s1,s2);
372: MatMult(C,x,s3);
374: VecNorm(s2,NORM_INFINITY,&r2);
375: if (r2 < tol) {
376: VecNorm(s3,NORM_INFINITY,&r1);
377: } else {
378: VecAXPY(s2,none,s3);
379: VecNorm(s2,NORM_INFINITY,&r1);
380: r1 /= r2;
381: }
382: if (r1 > tol) {
383: *flg = PETSC_FALSE;
384: PetscInfo2(A,"Error: %D-th MatTransposeMatMult() %g\n",k,(double)r1);
385: break;
386: }
387: }
388: PetscRandomDestroy(&rctx);
389: VecDestroy(&x);
390: VecDestroy(&s1);
391: VecDestroy(&s2);
392: VecDestroy(&s3);
393: return(0);
394: }
396: /*@
397: MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
399: Collective on Mat
401: Input Parameters:
402: + A - the first matrix
403: - B - the second matrix
404: - C - the third matrix
405: - n - number of random vectors to be tested
407: Output Parameter:
408: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
410: Level: intermediate
412: Concepts: matrices^equality between
413: @*/
414: PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
415: {
417: Vec x,s1,s2,s3;
418: PetscRandom rctx;
419: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
420: PetscInt am,an,bm,bn,cm,cn,k;
421: PetscScalar none = -1.0;
424: MatGetLocalSize(A,&am,&an);
425: MatGetLocalSize(B,&bm,&bn);
426: MatGetLocalSize(C,&cm,&cn);
427: 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);
429: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
430: PetscRandomSetFromOptions(rctx);
431: MatCreateVecs(B,&s1,&x);
432: MatCreateVecs(C,NULL,&s2);
433: VecDuplicate(s2,&s3);
435: *flg = PETSC_TRUE;
436: for (k=0; k<n; k++) {
437: VecSetRandom(x,rctx);
438: MatMultTranspose(B,x,s1);
439: MatMult(A,s1,s2);
440: MatMult(C,x,s3);
442: VecNorm(s2,NORM_INFINITY,&r2);
443: if (r2 < tol) {
444: VecNorm(s3,NORM_INFINITY,&r1);
445: } else {
446: VecAXPY(s2,none,s3);
447: VecNorm(s2,NORM_INFINITY,&r1);
448: r1 /= r2;
449: }
450: if (r1 > tol) {
451: *flg = PETSC_FALSE;
452: PetscInfo2(A,"Error: %D-th MatMatTransposeMult() %g\n",k,(double)r1);
453: break;
454: }
455: }
456: PetscRandomDestroy(&rctx);
457: VecDestroy(&x);
458: VecDestroy(&s1);
459: VecDestroy(&s2);
460: VecDestroy(&s3);
461: return(0);
462: }
464: /*@
465: MatIsLinear - Check if a shell matrix A is a linear operator.
467: Collective on Mat
469: Input Parameters:
470: + A - the shell matrix
471: - n - number of random vectors to be tested
473: Output Parameter:
474: . flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.
476: Level: intermediate
477: @*/
478: PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool *flg)
479: {
481: Vec x,y,s1,s2;
482: PetscRandom rctx;
483: PetscScalar a;
484: PetscInt k;
485: PetscReal norm,normA;
486: MPI_Comm comm;
487: PetscMPIInt rank;
491: PetscObjectGetComm((PetscObject)A,&comm);
492: MPI_Comm_rank(comm,&rank);
494: PetscRandomCreate(comm,&rctx);
495: PetscRandomSetFromOptions(rctx);
496: MatCreateVecs(A,&x,&s1);
497: VecDuplicate(x,&y);
498: VecDuplicate(s1,&s2);
500: *flg = PETSC_TRUE;
501: for (k=0; k<n; k++) {
502: VecSetRandom(x,rctx);
503: VecSetRandom(y,rctx);
504: if (!rank) {
505: PetscRandomGetValue(rctx,&a);
506: }
507: MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);
509: /* s2 = a*A*x + A*y */
510: MatMult(A,y,s2); /* s2 = A*y */
511: MatMult(A,x,s1); /* s1 = A*x */
512: VecAXPY(s2,a,s1); /* s2 = a s1 + s2 */
514: /* s1 = A * (a x + y) */
515: VecAXPY(y,a,x); /* y = a x + y */
516: MatMult(A,y,s1);
517: VecNorm(s1,NORM_INFINITY,&normA);
519: VecAXPY(s2,-1.0,s1); /* s2 = - s1 + s2 */
520: VecNorm(s2,NORM_INFINITY,&norm);
521: if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
522: *flg = PETSC_FALSE;
523: 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);
524: break;
525: }
526: }
527: PetscRandomDestroy(&rctx);
528: VecDestroy(&x);
529: VecDestroy(&y);
530: VecDestroy(&s1);
531: VecDestroy(&s2);
532: return(0);
533: }