xref: /petsc/src/mat/utils/multequal.c (revision 4fcfbe5ff1195060b3dcb0a53ce01d395fc1078b) 
1 #define PETSCMAT_DLL
2 
3 #include "private/matimpl.h"  /*I   "petscmat.h"  I*/
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "MatMultEqual"
7 /*@
8    MatMultEqual - Compares matrix-vector products of two matrices.
9 
10    Collective on Mat
11 
12    Input Parameters:
13 +  A - the first matrix
14 -  B - the second matrix
15 -  n - number of random vectors to be tested
16 
17    Output Parameter:
18 .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
19 
20    Level: intermediate
21 
22    Concepts: matrices^equality between
23 @*/
24 PetscErrorCode PETSCMAT_DLLEXPORT MatMultEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
25 {
26   PetscErrorCode ierr;
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;
32 
33   PetscFunctionBegin;
34   PetscValidHeaderSpecific(A,MAT_COOKIE,1);
35   PetscValidHeaderSpecific(B,MAT_COOKIE,2);
36   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
37   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
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);
39   PetscCheckSameComm(A,1,B,2);
40 
41 #if defined(PETSC_USE_SCALAR_SINGLE)
42   tol  = 1.e-5;
43 #endif
44   ierr = PetscRandomCreate(((PetscObject)A)->comm,&rctx);CHKERRQ(ierr);
45   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
46   ierr = VecCreate(((PetscObject)A)->comm,&x);CHKERRQ(ierr);
47   ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr);
48   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
49 
50   ierr = VecCreate(((PetscObject)A)->comm,&s1);CHKERRQ(ierr);
51   ierr = VecSetSizes(s1,am,PETSC_DECIDE);CHKERRQ(ierr);
52   ierr = VecSetFromOptions(s1);CHKERRQ(ierr);
53   ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
54 
55   *flg = PETSC_TRUE;
56   for (k=0; k<n; k++) {
57     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
58     ierr = MatMult(A,x,s1);CHKERRQ(ierr);
59     ierr = MatMult(B,x,s2);CHKERRQ(ierr);
60     ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
61     if (r2 < tol){
62       ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr);
63     } else {
64       ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr);
65       ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
66       r1 /= r2;
67     }
68     if (r1 > tol) {
69       *flg = PETSC_FALSE;
70       ierr = PetscInfo2(A,"Error: %D-th MatMult() %G\n",k,r1);CHKERRQ(ierr);
71       break;
72     }
73   }
74   ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
75   ierr = VecDestroy(x);CHKERRQ(ierr);
76   ierr = VecDestroy(s1);CHKERRQ(ierr);
77   ierr = VecDestroy(s2);CHKERRQ(ierr);
78   PetscFunctionReturn(0);
79 }
80 
81 #undef __FUNCT__
82 #define __FUNCT__ "MatMultAddEqual"
83 /*@
84    MatMultAddEqual - Compares matrix-vector products of two matrices.
85 
86    Collective on Mat
87 
88    Input Parameters:
89 +  A - the first matrix
90 -  B - the second matrix
91 -  n - number of random vectors to be tested
92 
93    Output Parameter:
94 .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
95 
96    Level: intermediate
97 
98    Concepts: matrices^equality between
99 @*/
100 PetscErrorCode PETSCMAT_DLLEXPORT MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
101 {
102   PetscErrorCode ierr;
103   Vec            x,y,s1,s2;
104   PetscRandom    rctx;
105   PetscReal      r1,r2,tol=1.e-10;
106   PetscInt       am,an,bm,bn,k;
107   PetscScalar    none = -1.0;
108 
109   PetscFunctionBegin;
110   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
111   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
112   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
113   PetscCheckSameComm(A,1,B,2);
114   ierr = PetscRandomCreate(((PetscObject)A)->comm,&rctx);CHKERRQ(ierr);
115   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
116   ierr = VecCreate(((PetscObject)A)->comm,&x);CHKERRQ(ierr);
117   ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr);
118   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
119 
120   ierr = VecCreate(((PetscObject)A)->comm,&s1);CHKERRQ(ierr);
121   ierr = VecSetSizes(s1,am,PETSC_DECIDE);CHKERRQ(ierr);
122   ierr = VecSetFromOptions(s1);CHKERRQ(ierr);
123   ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
124   ierr = VecDuplicate(s1,&y);CHKERRQ(ierr);
125 
126   *flg = PETSC_TRUE;
127   for (k=0; k<n; k++) {
128     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
129     ierr = VecSetRandom(y,rctx);CHKERRQ(ierr);
130     ierr = MatMultAdd(A,x,y,s1);CHKERRQ(ierr);
131     ierr = MatMultAdd(B,x,y,s2);CHKERRQ(ierr);
132     ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
133     if (r2 < tol){
134       ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr);
135     } else {
136       ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr);
137       ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
138       r1 /= r2;
139     }
140     if (r1 > tol) {
141       *flg = PETSC_FALSE;
142       ierr = PetscInfo2(A,"Error: %d-th MatMultAdd() %G\n",k,r1);CHKERRQ(ierr);
143       break;
144     }
145   }
146   ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
147   ierr = VecDestroy(x);CHKERRQ(ierr);
148   ierr = VecDestroy(y);CHKERRQ(ierr);
149   ierr = VecDestroy(s1);CHKERRQ(ierr);
150   ierr = VecDestroy(s2);CHKERRQ(ierr);
151   PetscFunctionReturn(0);
152 }
153 
154 #undef __FUNCT__
155 #define __FUNCT__ "MatMultTransposeEqual"
156 /*@
157    MatMultTransposeEqual - Compares matrix-vector products of two matrices.
158 
159    Collective on Mat
160 
161    Input Parameters:
162 +  A - the first matrix
163 -  B - the second matrix
164 -  n - number of random vectors to be tested
165 
166    Output Parameter:
167 .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
168 
169    Level: intermediate
170 
171    Concepts: matrices^equality between
172 @*/
173 PetscErrorCode PETSCMAT_DLLEXPORT MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
174 {
175   PetscErrorCode ierr;
176   Vec            x,s1,s2;
177   PetscRandom    rctx;
178   PetscReal      r1,r2,tol=1.e-10;
179   PetscInt       am,an,bm,bn,k;
180   PetscScalar    none = -1.0;
181 
182   PetscFunctionBegin;
183   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
184   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
185   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
186   PetscCheckSameComm(A,1,B,2);
187   ierr = PetscRandomCreate(((PetscObject)A)->comm,&rctx);CHKERRQ(ierr);
188   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
189   ierr = VecCreate(((PetscObject)A)->comm,&x);CHKERRQ(ierr);
190   ierr = VecSetSizes(x,am,PETSC_DECIDE);CHKERRQ(ierr);
191   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
192 
193   ierr = VecCreate(((PetscObject)A)->comm,&s1);CHKERRQ(ierr);
194   ierr = VecSetSizes(s1,an,PETSC_DECIDE);CHKERRQ(ierr);
195   ierr = VecSetFromOptions(s1);CHKERRQ(ierr);
196   ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
197 
198   *flg = PETSC_TRUE;
199   for (k=0; k<n; k++) {
200     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
201     ierr = MatMultTranspose(A,x,s1);CHKERRQ(ierr);
202     ierr = MatMultTranspose(B,x,s2);CHKERRQ(ierr);
203     ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
204     if (r2 < tol){
205       ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr);
206     } else {
207       ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr);
208       ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
209       r1 /= r2;
210     }
211     if (r1 > tol) {
212       *flg = PETSC_FALSE;
213       ierr = PetscInfo2(A,"Error: %d-th MatMultTranspose() %G\n",k,r1);CHKERRQ(ierr);
214       break;
215     }
216   }
217   ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
218   ierr = VecDestroy(x);CHKERRQ(ierr);
219   ierr = VecDestroy(s1);CHKERRQ(ierr);
220   ierr = VecDestroy(s2);CHKERRQ(ierr);
221   PetscFunctionReturn(0);
222 }
223 
224 #undef __FUNCT__
225 #define __FUNCT__ "MatMultTransposeAddEqual"
226 /*@
227    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
228 
229    Collective on Mat
230 
231    Input Parameters:
232 +  A - the first matrix
233 -  B - the second matrix
234 -  n - number of random vectors to be tested
235 
236    Output Parameter:
237 .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
238 
239    Level: intermediate
240 
241    Concepts: matrices^equality between
242 @*/
243 PetscErrorCode PETSCMAT_DLLEXPORT MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
244 {
245   PetscErrorCode ierr;
246   Vec            x,y,s1,s2;
247   PetscRandom    rctx;
248   PetscReal      r1,r2,tol=1.e-10;
249   PetscInt       am,an,bm,bn,k;
250   PetscScalar    none = -1.0;
251 
252   PetscFunctionBegin;
253   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
254   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
255   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
256   PetscCheckSameComm(A,1,B,2);
257   ierr = PetscRandomCreate(((PetscObject)A)->comm,&rctx);CHKERRQ(ierr);
258   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
259   ierr = VecCreate(((PetscObject)A)->comm,&x);CHKERRQ(ierr);
260   ierr = VecSetSizes(x,am,PETSC_DECIDE);CHKERRQ(ierr);
261   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
262 
263   ierr = VecCreate(((PetscObject)A)->comm,&s1);CHKERRQ(ierr);
264   ierr = VecSetSizes(s1,an,PETSC_DECIDE);CHKERRQ(ierr);
265   ierr = VecSetFromOptions(s1);CHKERRQ(ierr);
266   ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
267   ierr = VecDuplicate(s1,&y);CHKERRQ(ierr);
268 
269   *flg = PETSC_TRUE;
270   for (k=0; k<n; k++) {
271     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
272     ierr = VecSetRandom(y,rctx);CHKERRQ(ierr);
273     ierr = MatMultTransposeAdd(A,x,y,s1);CHKERRQ(ierr);
274     ierr = MatMultTransposeAdd(B,x,y,s2);CHKERRQ(ierr);
275     ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
276     if (r2 < tol){
277       ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr);
278     } else {
279       ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr);
280       ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
281       r1 /= r2;
282     }
283     if (r1 > tol) {
284       *flg = PETSC_FALSE;
285       ierr = PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %G\n",k,r1);CHKERRQ(ierr);
286       break;
287     }
288   }
289   ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
290   ierr = VecDestroy(x);CHKERRQ(ierr);
291   ierr = VecDestroy(y);CHKERRQ(ierr);
292   ierr = VecDestroy(s1);CHKERRQ(ierr);
293   ierr = VecDestroy(s2);CHKERRQ(ierr);
294   PetscFunctionReturn(0);
295 }
296