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