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: }