Actual source code: multequal.c
petsc-3.4.5 2014-06-29
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=1.e-10;
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);
40: #if defined(PETSC_USE_REAL_SINGLE)
41: tol = 1.e-5;
42: #endif
43: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
44: PetscRandomSetFromOptions(rctx);
45: VecCreate(PetscObjectComm((PetscObject)A),&x);
46: VecSetSizes(x,an,PETSC_DECIDE);
47: VecSetFromOptions(x);
49: VecCreate(PetscObjectComm((PetscObject)A),&s1);
50: VecSetSizes(s1,am,PETSC_DECIDE);
51: VecSetFromOptions(s1);
52: VecDuplicate(s1,&s2);
54: *flg = PETSC_TRUE;
55: for (k=0; k<n; k++) {
56: VecSetRandom(x,rctx);
57: MatMult(A,x,s1);
58: MatMult(B,x,s2);
59: VecNorm(s2,NORM_INFINITY,&r2);
60: if (r2 < tol) {
61: VecNorm(s1,NORM_INFINITY,&r1);
62: } else {
63: VecAXPY(s2,none,s1);
64: VecNorm(s2,NORM_INFINITY,&r1);
65: r1 /= r2;
66: }
67: if (r1 > tol) {
68: *flg = PETSC_FALSE;
69: PetscInfo2(A,"Error: %D-th MatMult() %G\n",k,r1);
70: break;
71: }
72: }
73: PetscRandomDestroy(&rctx);
74: VecDestroy(&x);
75: VecDestroy(&s1);
76: VecDestroy(&s2);
77: return(0);
78: }
82: /*@
83: MatMultAddEqual - Compares matrix-vector products of two matrices.
85: Collective on Mat
87: Input Parameters:
88: + A - the first matrix
89: - B - the second matrix
90: - n - number of random vectors to be tested
92: Output Parameter:
93: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
95: Level: intermediate
97: Concepts: matrices^equality between
98: @*/
99: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
100: {
102: Vec x,y,s1,s2;
103: PetscRandom rctx;
104: PetscReal r1,r2,tol=1.e-10;
105: PetscInt am,an,bm,bn,k;
106: PetscScalar none = -1.0;
109: MatGetLocalSize(A,&am,&an);
110: MatGetLocalSize(B,&bm,&bn);
111: 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);
113: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
114: PetscRandomSetFromOptions(rctx);
115: VecCreate(PetscObjectComm((PetscObject)A),&x);
116: VecSetSizes(x,an,PETSC_DECIDE);
117: VecSetFromOptions(x);
119: VecCreate(PetscObjectComm((PetscObject)A),&s1);
120: VecSetSizes(s1,am,PETSC_DECIDE);
121: VecSetFromOptions(s1);
122: VecDuplicate(s1,&s2);
123: VecDuplicate(s1,&y);
125: *flg = PETSC_TRUE;
126: for (k=0; k<n; k++) {
127: VecSetRandom(x,rctx);
128: VecSetRandom(y,rctx);
129: MatMultAdd(A,x,y,s1);
130: MatMultAdd(B,x,y,s2);
131: VecNorm(s2,NORM_INFINITY,&r2);
132: if (r2 < tol) {
133: VecNorm(s1,NORM_INFINITY,&r1);
134: } else {
135: VecAXPY(s2,none,s1);
136: VecNorm(s2,NORM_INFINITY,&r1);
137: r1 /= r2;
138: }
139: if (r1 > tol) {
140: *flg = PETSC_FALSE;
141: PetscInfo2(A,"Error: %d-th MatMultAdd() %G\n",k,r1);
142: break;
143: }
144: }
145: PetscRandomDestroy(&rctx);
146: VecDestroy(&x);
147: VecDestroy(&y);
148: VecDestroy(&s1);
149: VecDestroy(&s2);
150: return(0);
151: }
155: /*@
156: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
158: Collective on Mat
160: Input Parameters:
161: + A - the first matrix
162: - B - the second matrix
163: - n - number of random vectors to be tested
165: Output Parameter:
166: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
168: Level: intermediate
170: Concepts: matrices^equality between
171: @*/
172: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
173: {
175: Vec x,s1,s2;
176: PetscRandom rctx;
177: PetscReal r1,r2,tol=1.e-10;
178: PetscInt am,an,bm,bn,k;
179: PetscScalar none = -1.0;
182: MatGetLocalSize(A,&am,&an);
183: MatGetLocalSize(B,&bm,&bn);
184: 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);
186: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
187: PetscRandomSetFromOptions(rctx);
188: VecCreate(PetscObjectComm((PetscObject)A),&x);
189: VecSetSizes(x,am,PETSC_DECIDE);
190: VecSetFromOptions(x);
192: VecCreate(PetscObjectComm((PetscObject)A),&s1);
193: VecSetSizes(s1,an,PETSC_DECIDE);
194: VecSetFromOptions(s1);
195: VecDuplicate(s1,&s2);
197: *flg = PETSC_TRUE;
198: for (k=0; k<n; k++) {
199: VecSetRandom(x,rctx);
200: MatMultTranspose(A,x,s1);
201: MatMultTranspose(B,x,s2);
202: VecNorm(s2,NORM_INFINITY,&r2);
203: if (r2 < tol) {
204: VecNorm(s1,NORM_INFINITY,&r1);
205: } else {
206: VecAXPY(s2,none,s1);
207: VecNorm(s2,NORM_INFINITY,&r1);
208: r1 /= r2;
209: }
210: if (r1 > tol) {
211: *flg = PETSC_FALSE;
212: PetscInfo2(A,"Error: %d-th MatMultTranspose() %G\n",k,r1);
213: break;
214: }
215: }
216: PetscRandomDestroy(&rctx);
217: VecDestroy(&x);
218: VecDestroy(&s1);
219: VecDestroy(&s2);
220: return(0);
221: }
225: /*@
226: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
228: Collective on Mat
230: Input Parameters:
231: + A - the first matrix
232: - B - the second matrix
233: - n - number of random vectors to be tested
235: Output Parameter:
236: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
238: Level: intermediate
240: Concepts: matrices^equality between
241: @*/
242: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
243: {
245: Vec x,y,s1,s2;
246: PetscRandom rctx;
247: PetscReal r1,r2,tol=1.e-10;
248: PetscInt am,an,bm,bn,k;
249: PetscScalar none = -1.0;
252: MatGetLocalSize(A,&am,&an);
253: MatGetLocalSize(B,&bm,&bn);
254: 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);
256: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
257: PetscRandomSetFromOptions(rctx);
258: VecCreate(PetscObjectComm((PetscObject)A),&x);
259: VecSetSizes(x,am,PETSC_DECIDE);
260: VecSetFromOptions(x);
262: VecCreate(PetscObjectComm((PetscObject)A),&s1);
263: VecSetSizes(s1,an,PETSC_DECIDE);
264: VecSetFromOptions(s1);
265: VecDuplicate(s1,&s2);
266: VecDuplicate(s1,&y);
268: *flg = PETSC_TRUE;
269: for (k=0; k<n; k++) {
270: VecSetRandom(x,rctx);
271: VecSetRandom(y,rctx);
272: MatMultTransposeAdd(A,x,y,s1);
273: MatMultTransposeAdd(B,x,y,s2);
274: VecNorm(s2,NORM_INFINITY,&r2);
275: if (r2 < tol) {
276: VecNorm(s1,NORM_INFINITY,&r1);
277: } else {
278: VecAXPY(s2,none,s1);
279: VecNorm(s2,NORM_INFINITY,&r1);
280: r1 /= r2;
281: }
282: if (r1 > tol) {
283: *flg = PETSC_FALSE;
284: PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %G\n",k,r1);
285: break;
286: }
287: }
288: PetscRandomDestroy(&rctx);
289: VecDestroy(&x);
290: VecDestroy(&y);
291: VecDestroy(&s1);
292: VecDestroy(&s2);
293: return(0);
294: }