Actual source code: multequal.c
1: #define PETSCMAT_DLL
3: #include include/private/matimpl.h
7: /*@
8: MatMultEqual - Compares matrix-vector products of two matrices.
10: Collective on Mat
12: Input Parameters:
13: + A - the first matrix
14: - B - the second matrix
15: - n - number of random vectors to be tested
17: Output Parameter:
18: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
20: Level: intermediate
22: Concepts: matrices^equality between
23: @*/
24: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
25: {
27: Vec x,s1,s2;
28: PetscRandom rctx;
29: PetscReal r1,r2,tol=1.e-10;
30: PetscInt am,an,bm,bn,k;
31: PetscScalar none = -1.0;
36: MatGetLocalSize(A,&am,&an);
37: MatGetLocalSize(B,&bm,&bn);
38: if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
40: PetscRandomCreate(A->comm,&rctx);
41: PetscRandomSetFromOptions(rctx);
42: VecCreate(A->comm,&x);
43: VecSetSizes(x,an,PETSC_DECIDE);
44: VecSetFromOptions(x);
45:
46: VecCreate(A->comm,&s1);
47: VecSetSizes(s1,am,PETSC_DECIDE);
48: VecSetFromOptions(s1);
49: VecDuplicate(s1,&s2);
50:
51: *flg = PETSC_TRUE;
52: for (k=0; k<n; k++) {
53: VecSetRandom(x,rctx);
54: MatMult(A,x,s1);
55: MatMult(B,x,s2);
56: VecNorm(s2,NORM_INFINITY,&r2);
57: if (r2 < tol){
58: VecNorm(s1,NORM_INFINITY,&r1);
59: } else {
60: VecAXPY(s2,none,s1);
61: VecNorm(s2,NORM_INFINITY,&r1);
62: r1 /= r2;
63: }
64: if (r1 > tol) {
65: *flg = PETSC_FALSE;
66: PetscInfo2(0,"Error: %D-th MatMult() %G\n",k,r1);
67: break;
68: }
69: }
70: PetscRandomDestroy(rctx);
71: VecDestroy(x);
72: VecDestroy(s1);
73: VecDestroy(s2);
74: return(0);
75: }
79: /*@
80: MatMultAddEqual - Compares matrix-vector products of two matrices.
82: Collective on Mat
84: Input Parameters:
85: + A - the first matrix
86: - B - the second matrix
87: - n - number of random vectors to be tested
89: Output Parameter:
90: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
92: Level: intermediate
94: Concepts: matrices^equality between
95: @*/
96: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
97: {
99: Vec x,y,s1,s2;
100: PetscRandom rctx;
101: PetscReal r1,r2,tol=1.e-10;
102: PetscInt am,an,bm,bn,k;
103: PetscScalar none = -1.0;
106: MatGetLocalSize(A,&am,&an);
107: MatGetLocalSize(B,&bm,&bn);
108: if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
110: PetscRandomCreate(A->comm,&rctx);
111: PetscRandomSetFromOptions(rctx);
112: VecCreate(A->comm,&x);
113: VecSetSizes(x,an,PETSC_DECIDE);
114: VecSetFromOptions(x);
116: VecCreate(A->comm,&s1);
117: VecSetSizes(s1,am,PETSC_DECIDE);
118: VecSetFromOptions(s1);
119: VecDuplicate(s1,&s2);
120: VecDuplicate(s1,&y);
121:
122: *flg = PETSC_TRUE;
123: for (k=0; k<n; k++) {
124: VecSetRandom(x,rctx);
125: VecSetRandom(y,rctx);
126: MatMultAdd(A,x,y,s1);
127: MatMultAdd(B,x,y,s2);
128: VecNorm(s2,NORM_INFINITY,&r2);
129: if (r2 < tol){
130: VecNorm(s1,NORM_INFINITY,&r1);
131: } else {
132: VecAXPY(s2,none,s1);
133: VecNorm(s2,NORM_INFINITY,&r1);
134: r1 /= r2;
135: }
136: if (r1 > tol) {
137: *flg = PETSC_FALSE;
138: PetscInfo2(0,"Error: %d-th MatMultAdd() %G\n",k,r1);
139: break;
140: }
141: }
142: PetscRandomDestroy(rctx);
143: VecDestroy(x);
144: VecDestroy(y);
145: VecDestroy(s1);
146: VecDestroy(s2);
147: return(0);
148: }
152: /*@
153: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
155: Collective on Mat
157: Input Parameters:
158: + A - the first matrix
159: - B - the second matrix
160: - n - number of random vectors to be tested
162: Output Parameter:
163: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
165: Level: intermediate
167: Concepts: matrices^equality between
168: @*/
169: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
170: {
172: Vec x,s1,s2;
173: PetscRandom rctx;
174: PetscReal r1,r2,tol=1.e-10;
175: PetscInt am,an,bm,bn,k;
176: PetscScalar none = -1.0;
179: MatGetLocalSize(A,&am,&an);
180: MatGetLocalSize(B,&bm,&bn);
181: if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
183: PetscRandomCreate(A->comm,&rctx);
184: PetscRandomSetFromOptions(rctx);
185: VecCreate(A->comm,&x);
186: VecSetSizes(x,am,PETSC_DECIDE);
187: VecSetFromOptions(x);
188:
189: VecCreate(A->comm,&s1);
190: VecSetSizes(s1,an,PETSC_DECIDE);
191: VecSetFromOptions(s1);
192: VecDuplicate(s1,&s2);
193:
194: *flg = PETSC_TRUE;
195: for (k=0; k<n; k++) {
196: VecSetRandom(x,rctx);
197: MatMultTranspose(A,x,s1);
198: MatMultTranspose(B,x,s2);
199: VecNorm(s2,NORM_INFINITY,&r2);
200: if (r2 < tol){
201: VecNorm(s1,NORM_INFINITY,&r1);
202: } else {
203: VecAXPY(s2,none,s1);
204: VecNorm(s2,NORM_INFINITY,&r1);
205: r1 /= r2;
206: }
207: if (r1 > tol) {
208: *flg = PETSC_FALSE;
209: PetscInfo2(0,"Error: %d-th MatMultTranspose() %G\n",k,r1);
210: break;
211: }
212: }
213: PetscRandomDestroy(rctx);
214: VecDestroy(x);
215: VecDestroy(s1);
216: VecDestroy(s2);
217: return(0);
218: }
222: /*@
223: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
225: Collective on Mat
227: Input Parameters:
228: + A - the first matrix
229: - B - the second matrix
230: - n - number of random vectors to be tested
232: Output Parameter:
233: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
235: Level: intermediate
237: Concepts: matrices^equality between
238: @*/
239: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
240: {
242: Vec x,y,s1,s2;
243: PetscRandom rctx;
244: PetscReal r1,r2,tol=1.e-10;
245: PetscInt am,an,bm,bn,k;
246: PetscScalar none = -1.0;
249: MatGetLocalSize(A,&am,&an);
250: MatGetLocalSize(B,&bm,&bn);
251: if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
253: PetscRandomCreate(A->comm,&rctx);
254: PetscRandomSetFromOptions(rctx);
255: VecCreate(A->comm,&x);
256: VecSetSizes(x,am,PETSC_DECIDE);
257: VecSetFromOptions(x);
259: VecCreate(A->comm,&s1);
260: VecSetSizes(s1,an,PETSC_DECIDE);
261: VecSetFromOptions(s1);
262: VecDuplicate(s1,&s2);
263: VecDuplicate(s1,&y);
264:
265: *flg = PETSC_TRUE;
266: for (k=0; k<n; k++) {
267: VecSetRandom(x,rctx);
268: VecSetRandom(y,rctx);
269: MatMultTransposeAdd(A,x,y,s1);
270: MatMultTransposeAdd(B,x,y,s2);
271: VecNorm(s2,NORM_INFINITY,&r2);
272: if (r2 < tol){
273: VecNorm(s1,NORM_INFINITY,&r1);
274: } else {
275: VecAXPY(s2,none,s1);
276: VecNorm(s2,NORM_INFINITY,&r1);
277: r1 /= r2;
278: }
279: if (r1 > tol) {
280: *flg = PETSC_FALSE;
281: PetscInfo2(0,"Error: %d-th MatMultTransposeAdd() %G\n",k,r1);
282: break;
283: }
284: }
285: PetscRandomDestroy(rctx);
286: VecDestroy(x);
287: VecDestroy(y);
288: VecDestroy(s1);
289: VecDestroy(s2);
290: return(0);
291: }