Actual source code: multequal.c
petsc-3.5.4 2015-05-23
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: MatGetVecs(A,&x,&s1);
46: VecDuplicate(s1,&s2);
48: *flg = PETSC_TRUE;
49: for (k=0; k<n; k++) {
50: VecSetRandom(x,rctx);
51: MatMult(A,x,s1);
52: MatMult(B,x,s2);
53: VecNorm(s2,NORM_INFINITY,&r2);
54: if (r2 < tol) {
55: VecNorm(s1,NORM_INFINITY,&r1);
56: } else {
57: VecAXPY(s2,none,s1);
58: VecNorm(s2,NORM_INFINITY,&r1);
59: r1 /= r2;
60: }
61: if (r1 > tol) {
62: *flg = PETSC_FALSE;
63: PetscInfo2(A,"Error: %D-th MatMult() %g\n",k,(double)r1);
64: break;
65: }
66: }
67: PetscRandomDestroy(&rctx);
68: VecDestroy(&x);
69: VecDestroy(&s1);
70: VecDestroy(&s2);
71: return(0);
72: }
76: /*@
77: MatMultAddEqual - Compares matrix-vector products of two matrices.
79: Collective on Mat
81: Input Parameters:
82: + A - the first matrix
83: - B - the second matrix
84: - n - number of random vectors to be tested
86: Output Parameter:
87: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
89: Level: intermediate
91: Concepts: matrices^equality between
92: @*/
93: PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
94: {
96: Vec x,y,s1,s2;
97: PetscRandom rctx;
98: PetscReal r1,r2,tol=1.e-10;
99: PetscInt am,an,bm,bn,k;
100: PetscScalar none = -1.0;
103: MatGetLocalSize(A,&am,&an);
104: MatGetLocalSize(B,&bm,&bn);
105: 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);
107: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
108: PetscRandomSetFromOptions(rctx);
109: MatGetVecs(A,&x,&s1);
110: VecDuplicate(s1,&s2);
111: VecDuplicate(s1,&y);
113: *flg = PETSC_TRUE;
114: for (k=0; k<n; k++) {
115: VecSetRandom(x,rctx);
116: VecSetRandom(y,rctx);
117: MatMultAdd(A,x,y,s1);
118: MatMultAdd(B,x,y,s2);
119: VecNorm(s2,NORM_INFINITY,&r2);
120: if (r2 < tol) {
121: VecNorm(s1,NORM_INFINITY,&r1);
122: } else {
123: VecAXPY(s2,none,s1);
124: VecNorm(s2,NORM_INFINITY,&r1);
125: r1 /= r2;
126: }
127: if (r1 > tol) {
128: *flg = PETSC_FALSE;
129: PetscInfo2(A,"Error: %d-th MatMultAdd() %g\n",k,(double)r1);
130: break;
131: }
132: }
133: PetscRandomDestroy(&rctx);
134: VecDestroy(&x);
135: VecDestroy(&y);
136: VecDestroy(&s1);
137: VecDestroy(&s2);
138: return(0);
139: }
143: /*@
144: MatMultTransposeEqual - Compares matrix-vector products of two matrices.
146: Collective on Mat
148: Input Parameters:
149: + A - the first matrix
150: - B - the second matrix
151: - n - number of random vectors to be tested
153: Output Parameter:
154: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
156: Level: intermediate
158: Concepts: matrices^equality between
159: @*/
160: PetscErrorCode MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
161: {
163: Vec x,s1,s2;
164: PetscRandom rctx;
165: PetscReal r1,r2,tol=1.e-10;
166: PetscInt am,an,bm,bn,k;
167: PetscScalar none = -1.0;
170: MatGetLocalSize(A,&am,&an);
171: MatGetLocalSize(B,&bm,&bn);
172: 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);
174: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
175: PetscRandomSetFromOptions(rctx);
176: MatGetVecs(A,&s1,&x);
177: VecDuplicate(s1,&s2);
179: *flg = PETSC_TRUE;
180: for (k=0; k<n; k++) {
181: VecSetRandom(x,rctx);
182: MatMultTranspose(A,x,s1);
183: MatMultTranspose(B,x,s2);
184: VecNorm(s2,NORM_INFINITY,&r2);
185: if (r2 < tol) {
186: VecNorm(s1,NORM_INFINITY,&r1);
187: } else {
188: VecAXPY(s2,none,s1);
189: VecNorm(s2,NORM_INFINITY,&r1);
190: r1 /= r2;
191: }
192: if (r1 > tol) {
193: *flg = PETSC_FALSE;
194: PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);
195: break;
196: }
197: }
198: PetscRandomDestroy(&rctx);
199: VecDestroy(&x);
200: VecDestroy(&s1);
201: VecDestroy(&s2);
202: return(0);
203: }
207: /*@
208: MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
210: Collective on Mat
212: Input Parameters:
213: + A - the first matrix
214: - B - the second matrix
215: - n - number of random vectors to be tested
217: Output Parameter:
218: . flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
220: Level: intermediate
222: Concepts: matrices^equality between
223: @*/
224: PetscErrorCode MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
225: {
227: Vec x,y,s1,s2;
228: PetscRandom rctx;
229: PetscReal r1,r2,tol=1.e-10;
230: PetscInt am,an,bm,bn,k;
231: PetscScalar none = -1.0;
234: MatGetLocalSize(A,&am,&an);
235: MatGetLocalSize(B,&bm,&bn);
236: 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);
238: PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
239: PetscRandomSetFromOptions(rctx);
240: MatGetVecs(A,&s1,&x);
241: VecDuplicate(s1,&s2);
242: VecDuplicate(s1,&y);
244: *flg = PETSC_TRUE;
245: for (k=0; k<n; k++) {
246: VecSetRandom(x,rctx);
247: VecSetRandom(y,rctx);
248: MatMultTransposeAdd(A,x,y,s1);
249: MatMultTransposeAdd(B,x,y,s2);
250: VecNorm(s2,NORM_INFINITY,&r2);
251: if (r2 < tol) {
252: VecNorm(s1,NORM_INFINITY,&r1);
253: } else {
254: VecAXPY(s2,none,s1);
255: VecNorm(s2,NORM_INFINITY,&r1);
256: r1 /= r2;
257: }
258: if (r1 > tol) {
259: *flg = PETSC_FALSE;
260: PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);
261: break;
262: }
263: }
264: PetscRandomDestroy(&rctx);
265: VecDestroy(&x);
266: VecDestroy(&y);
267: VecDestroy(&s1);
268: VecDestroy(&s2);
269: return(0);
270: }