Actual source code: multequal.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
6: /*@
7: MatMultEqual - Compares matrix-vector products of two matrices.
9: Collective on Mat
11: Input Parameters:
12: + A - the first matrix
13: - B - the second matrix
14: - n - number of random vectors to be tested
16: Output Parameter:
17: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
19: Level: intermediate
21: Concepts: matrices^equality between
22: @*/
23: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
24: {
26: Vec x,s1,s2;
27: PetscRandom rctx;
28: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
29: PetscInt am,an,bm,bn,k;
30: PetscScalar none = -1.0;
35: MatGetLocalSize(A,&am,&an);
36: MatGetLocalSize(B,&bm,&bn);
37: 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);
39: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
40: PetscRandomSetFromOptions(rctx);
41: MatCreateVecs(A,&x,&s1);
42: VecDuplicate(s1,&s2);
44: *flg = PETSC_TRUE;
45: for (k=0; k<n; k++) {
46: VecSetRandom(x,rctx);
47: MatMult(A,x,s1);
48: MatMult(B,x,s2);
49: VecNorm(s2,NORM_INFINITY,&r2);
50: if (r2 < tol) {
51: VecNorm(s1,NORM_INFINITY,&r1);
52: } else {
53: VecAXPY(s2,none,s1);
54: VecNorm(s2,NORM_INFINITY,&r1);
55: r1 /= r2;
56: }
57: if (r1 > tol) {
58: *flg = PETSC_FALSE;
59: PetscInfo2(A,"Error: %D-th MatMult() %g\n",k,(double)r1);
60: break;
61: }
62: }
63: PetscRandomDestroy(&rctx);
64: VecDestroy(&x);
65: VecDestroy(&s1);
66: VecDestroy(&s2);
67: return(0);
68: }
72: /*@
73: MatMultAddEqual - Compares matrix-vector products of two matrices.
75: Collective on Mat
77: Input Parameters:
78: + A - the first matrix
79: - B - the second matrix
80: - n - number of random vectors to be tested
82: Output Parameter:
83: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
85: Level: intermediate
87: Concepts: matrices^equality between
88: @*/
89: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
90: {
92: Vec x,y,s1,s2;
93: PetscRandom rctx;
94: PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
95: PetscInt am,an,bm,bn,k;
96: PetscScalar none = -1.0;
99: MatGetLocalSize(A,&am,&an);
100: MatGetLocalSize(B,&bm,&bn);
101: 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);
103: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
104: PetscRandomSetFromOptions(rctx);
105: MatCreateVecs(A,&x,&s1);
106: VecDuplicate(s1,&s2);
107: VecDuplicate(s1,&y);
109: *flg = PETSC_TRUE;
110: for (k=0; k<n; k++) {
111: VecSetRandom(x,rctx);
112: VecSetRandom(y,rctx);
113: MatMultAdd(A,x,y,s1);
114: MatMultAdd(B,x,y,s2);
115: VecNorm(s2,NORM_INFINITY,&r2);
116: if (r2 < tol) {
117: VecNorm(s1,NORM_INFINITY,&r1);
118: } else {
119: VecAXPY(s2,none,s1);
120: VecNorm(s2,NORM_INFINITY,&r1);
121: r1 /= r2;
122: }
123: if (r1 > tol) {
124: *flg = PETSC_FALSE;
125: PetscInfo2(A,"Error: %d-th MatMultAdd() %g\n",k,(double)r1);
126: break;
127: }
128: }
129: PetscRandomDestroy(&rctx);
130: VecDestroy(&x);
131: VecDestroy(&y);
132: VecDestroy(&s1);
133: VecDestroy(&s2);
134: return(0);
135: }
139: /*@
140: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
142: Collective on Mat
144: Input Parameters:
145: + A - the first matrix
146: - B - the second matrix
147: - n - number of random vectors to be tested
149: Output Parameter:
150: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
152: Level: intermediate
154: Concepts: matrices^equality between
155: @*/
156: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
157: {
159: Vec x,s1,s2;
160: PetscRandom rctx;
161: PetscReal r1,r2,tol= PETSC_SQRT_MACHINE_EPSILON;
162: PetscInt am,an,bm,bn,k;
163: PetscScalar none = -1.0;
166: MatGetLocalSize(A,&am,&an);
167: MatGetLocalSize(B,&bm,&bn);
168: 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);
170: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
171: PetscRandomSetFromOptions(rctx);
172: MatCreateVecs(A,&s1,&x);
173: VecDuplicate(s1,&s2);
175: *flg = PETSC_TRUE;
176: for (k=0; k<n; k++) {
177: VecSetRandom(x,rctx);
178: MatMultTranspose(A,x,s1);
179: MatMultTranspose(B,x,s2);
180: VecNorm(s2,NORM_INFINITY,&r2);
181: if (r2 < tol) {
182: VecNorm(s1,NORM_INFINITY,&r1);
183: } else {
184: VecAXPY(s2,none,s1);
185: VecNorm(s2,NORM_INFINITY,&r1);
186: r1 /= r2;
187: }
188: if (r1 > tol) {
189: *flg = PETSC_FALSE;
190: PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);
191: break;
192: }
193: }
194: PetscRandomDestroy(&rctx);
195: VecDestroy(&x);
196: VecDestroy(&s1);
197: VecDestroy(&s2);
198: return(0);
199: }
203: /*@
204: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
206: Collective on Mat
208: Input Parameters:
209: + A - the first matrix
210: - B - the second matrix
211: - n - number of random vectors to be tested
213: Output Parameter:
214: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
216: Level: intermediate
218: Concepts: matrices^equality between
219: @*/
220: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
221: {
223: Vec x,y,s1,s2;
224: PetscRandom rctx;
225: PetscReal r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
226: PetscInt am,an,bm,bn,k;
227: PetscScalar none = -1.0;
230: MatGetLocalSize(A,&am,&an);
231: MatGetLocalSize(B,&bm,&bn);
232: 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);
234: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
235: PetscRandomSetFromOptions(rctx);
236: MatCreateVecs(A,&s1,&x);
237: VecDuplicate(s1,&s2);
238: VecDuplicate(s1,&y);
240: *flg = PETSC_TRUE;
241: for (k=0; k<n; k++) {
242: VecSetRandom(x,rctx);
243: VecSetRandom(y,rctx);
244: MatMultTransposeAdd(A,x,y,s1);
245: MatMultTransposeAdd(B,x,y,s2);
246: VecNorm(s2,NORM_INFINITY,&r2);
247: if (r2 < tol) {
248: VecNorm(s1,NORM_INFINITY,&r1);
249: } else {
250: VecAXPY(s2,none,s1);
251: VecNorm(s2,NORM_INFINITY,&r1);
252: r1 /= r2;
253: }
254: if (r1 > tol) {
255: *flg = PETSC_FALSE;
256: PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);
257: break;
258: }
259: }
260: PetscRandomDestroy(&rctx);
261: VecDestroy(&x);
262: VecDestroy(&y);
263: VecDestroy(&s1);
264: VecDestroy(&s2);
265: return(0);
266: }
270: /*@
271: MatMatMultEqual - Test A*B*x = C*x for n random vector x
273: Collective on Mat
275: Input Parameters:
276: + A - the first matrix
277: - B - the second matrix
278: - C - the third matrix
279: - n - number of random vectors to be tested
281: Output Parameter:
282: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
284: Level: intermediate
286: Concepts: matrices^equality between
287: @*/
288: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
289: {
291: Vec x,s1,s2,s3;
292: PetscRandom rctx;
293: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
294: PetscInt am,an,bm,bn,cm,cn,k;
295: PetscScalar none = -1.0;
298: MatGetLocalSize(A,&am,&an);
299: MatGetLocalSize(B,&bm,&bn);
300: MatGetLocalSize(C,&cm,&cn);
301: 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);
303: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
304: PetscRandomSetFromOptions(rctx);
305: MatCreateVecs(B,&x,&s1);
306: MatCreateVecs(C,NULL,&s2);
307: VecDuplicate(s2,&s3);
309: *flg = PETSC_TRUE;
310: for (k=0; k<n; k++) {
311: VecSetRandom(x,rctx);
312: MatMult(B,x,s1);
313: MatMult(A,s1,s2);
314: MatMult(C,x,s3);
316: VecNorm(s2,NORM_INFINITY,&r2);
317: if (r2 < tol) {
318: VecNorm(s3,NORM_INFINITY,&r1);
319: } else {
320: VecAXPY(s2,none,s3);
321: VecNorm(s2,NORM_INFINITY,&r1);
322: r1 /= r2;
323: }
324: if (r1 > tol) {
325: *flg = PETSC_FALSE;
326: PetscInfo2(A,"Error: %D-th MatMatMult() %g\n",k,(double)r1);
327: break;
328: }
329: }
330: PetscRandomDestroy(&rctx);
331: VecDestroy(&x);
332: VecDestroy(&s1);
333: VecDestroy(&s2);
334: VecDestroy(&s3);
335: return(0);
336: }
340: /*@
341: MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
343: Collective on Mat
345: Input Parameters:
346: + A - the first matrix
347: - B - the second matrix
348: - C - the third matrix
349: - n - number of random vectors to be tested
351: Output Parameter:
352: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
354: Level: intermediate
356: Concepts: matrices^equality between
357: @*/
358: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
359: {
361: Vec x,s1,s2,s3;
362: PetscRandom rctx;
363: PetscReal r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
364: PetscInt am,an,bm,bn,cm,cn,k;
365: PetscScalar none = -1.0;
368: MatGetLocalSize(A,&am,&an);
369: MatGetLocalSize(B,&bm,&bn);
370: MatGetLocalSize(C,&cm,&cn);
371: 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);
373: PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
374: PetscRandomSetFromOptions(rctx);
375: MatCreateVecs(B,&x,&s1);
376: MatCreateVecs(C,NULL,&s2);
377: VecDuplicate(s2,&s3);
379: *flg = PETSC_TRUE;
380: for (k=0; k<n; k++) {
381: VecSetRandom(x,rctx);
382: MatMult(B,x,s1);
383: MatMultTranspose(A,s1,s2);
384: MatMult(C,x,s3);
386: VecNorm(s2,NORM_INFINITY,&r2);
387: if (r2 < tol) {
388: VecNorm(s3,NORM_INFINITY,&r1);
389: } else {
390: VecAXPY(s2,none,s3);
391: VecNorm(s2,NORM_INFINITY,&r1);
392: r1 /= r2;
393: }
394: if (r1 > tol) {
395: *flg = PETSC_FALSE;
396: PetscInfo2(A,"Error: %D-th MatTransposeMatMult() %g\n",k,(double)r1);
397: break;
398: }
399: }
400: PetscRandomDestroy(&rctx);
401: VecDestroy(&x);
402: VecDestroy(&s1);
403: VecDestroy(&s2);
404: VecDestroy(&s3);
405: return(0);
406: }