Actual source code: multequal.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  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=PETSC_SQRT_MACHINE_EPSILON;
 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);
 39:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
 40:   PetscRandomSetFromOptions(rctx);
 41:   MatCreateVecs(A,&x,&s1);
 42:   VecDuplicate(s1,&s2);

 44:   *flg = PETSC_TRUE;
 45:   for (k=0; k<n; k++) {
 46:     VecSetRandom(x,rctx);
 47:     MatMult(A,x,s1);
 48:     MatMult(B,x,s2);
 49:     VecNorm(s2,NORM_INFINITY,&r2);
 50:     if (r2 < tol) {
 51:       VecNorm(s1,NORM_INFINITY,&r1);
 52:     } else {
 53:       VecAXPY(s2,none,s1);
 54:       VecNorm(s2,NORM_INFINITY,&r1);
 55:       r1  /= r2;
 56:     }
 57:     if (r1 > tol) {
 58:       *flg = PETSC_FALSE;
 59:       PetscInfo2(A,"Error: %D-th MatMult() %g\n",k,(double)r1);
 60:       break;
 61:     }
 62:   }
 63:   PetscRandomDestroy(&rctx);
 64:   VecDestroy(&x);
 65:   VecDestroy(&s1);
 66:   VecDestroy(&s2);
 67:   return(0);
 68: }

 72: /*@
 73:    MatMultAddEqual - Compares matrix-vector products of two matrices.

 75:    Collective on Mat

 77:    Input Parameters:
 78: +  A - the first matrix
 79: -  B - the second matrix
 80: -  n - number of random vectors to be tested

 82:    Output Parameter:
 83: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

 85:    Level: intermediate

 87:    Concepts: matrices^equality between
 88: @*/
 89: PetscErrorCode  MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
 90: {
 92:   Vec            x,y,s1,s2;
 93:   PetscRandom    rctx;
 94:   PetscReal      r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
 95:   PetscInt       am,an,bm,bn,k;
 96:   PetscScalar    none = -1.0;

 99:   MatGetLocalSize(A,&am,&an);
100:   MatGetLocalSize(B,&bm,&bn);
101:   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);
103:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
104:   PetscRandomSetFromOptions(rctx);
105:   MatCreateVecs(A,&x,&s1);
106:   VecDuplicate(s1,&s2);
107:   VecDuplicate(s1,&y);

109:   *flg = PETSC_TRUE;
110:   for (k=0; k<n; k++) {
111:     VecSetRandom(x,rctx);
112:     VecSetRandom(y,rctx);
113:     MatMultAdd(A,x,y,s1);
114:     MatMultAdd(B,x,y,s2);
115:     VecNorm(s2,NORM_INFINITY,&r2);
116:     if (r2 < tol) {
117:       VecNorm(s1,NORM_INFINITY,&r1);
118:     } else {
119:       VecAXPY(s2,none,s1);
120:       VecNorm(s2,NORM_INFINITY,&r1);
121:       r1  /= r2;
122:     }
123:     if (r1 > tol) {
124:       *flg = PETSC_FALSE;
125:       PetscInfo2(A,"Error: %d-th MatMultAdd() %g\n",k,(double)r1);
126:       break;
127:     }
128:   }
129:   PetscRandomDestroy(&rctx);
130:   VecDestroy(&x);
131:   VecDestroy(&y);
132:   VecDestroy(&s1);
133:   VecDestroy(&s2);
134:   return(0);
135: }

139: /*@
140:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

142:    Collective on Mat

144:    Input Parameters:
145: +  A - the first matrix
146: -  B - the second matrix
147: -  n - number of random vectors to be tested

149:    Output Parameter:
150: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

152:    Level: intermediate

154:    Concepts: matrices^equality between
155: @*/
156: PetscErrorCode  MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
157: {
159:   Vec            x,s1,s2;
160:   PetscRandom    rctx;
161:   PetscReal      r1,r2,tol= PETSC_SQRT_MACHINE_EPSILON;
162:   PetscInt       am,an,bm,bn,k;
163:   PetscScalar    none = -1.0;

166:   MatGetLocalSize(A,&am,&an);
167:   MatGetLocalSize(B,&bm,&bn);
168:   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);
170:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
171:   PetscRandomSetFromOptions(rctx);
172:   MatCreateVecs(A,&s1,&x);
173:   VecDuplicate(s1,&s2);

175:   *flg = PETSC_TRUE;
176:   for (k=0; k<n; k++) {
177:     VecSetRandom(x,rctx);
178:     MatMultTranspose(A,x,s1);
179:     MatMultTranspose(B,x,s2);
180:     VecNorm(s2,NORM_INFINITY,&r2);
181:     if (r2 < tol) {
182:       VecNorm(s1,NORM_INFINITY,&r1);
183:     } else {
184:       VecAXPY(s2,none,s1);
185:       VecNorm(s2,NORM_INFINITY,&r1);
186:       r1  /= r2;
187:     }
188:     if (r1 > tol) {
189:       *flg = PETSC_FALSE;
190:       PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);
191:       break;
192:     }
193:   }
194:   PetscRandomDestroy(&rctx);
195:   VecDestroy(&x);
196:   VecDestroy(&s1);
197:   VecDestroy(&s2);
198:   return(0);
199: }

203: /*@
204:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

206:    Collective on Mat

208:    Input Parameters:
209: +  A - the first matrix
210: -  B - the second matrix
211: -  n - number of random vectors to be tested

213:    Output Parameter:
214: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

216:    Level: intermediate

218:    Concepts: matrices^equality between
219: @*/
220: PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
221: {
223:   Vec            x,y,s1,s2;
224:   PetscRandom    rctx;
225:   PetscReal      r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
226:   PetscInt       am,an,bm,bn,k;
227:   PetscScalar    none = -1.0;

230:   MatGetLocalSize(A,&am,&an);
231:   MatGetLocalSize(B,&bm,&bn);
232:   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);
234:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
235:   PetscRandomSetFromOptions(rctx);
236:   MatCreateVecs(A,&s1,&x);
237:   VecDuplicate(s1,&s2);
238:   VecDuplicate(s1,&y);

240:   *flg = PETSC_TRUE;
241:   for (k=0; k<n; k++) {
242:     VecSetRandom(x,rctx);
243:     VecSetRandom(y,rctx);
244:     MatMultTransposeAdd(A,x,y,s1);
245:     MatMultTransposeAdd(B,x,y,s2);
246:     VecNorm(s2,NORM_INFINITY,&r2);
247:     if (r2 < tol) {
248:       VecNorm(s1,NORM_INFINITY,&r1);
249:     } else {
250:       VecAXPY(s2,none,s1);
251:       VecNorm(s2,NORM_INFINITY,&r1);
252:       r1  /= r2;
253:     }
254:     if (r1 > tol) {
255:       *flg = PETSC_FALSE;
256:       PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);
257:       break;
258:     }
259:   }
260:   PetscRandomDestroy(&rctx);
261:   VecDestroy(&x);
262:   VecDestroy(&y);
263:   VecDestroy(&s1);
264:   VecDestroy(&s2);
265:   return(0);
266: }

270: /*@
271:    MatMatMultEqual - Test A*B*x = C*x for n random vector x 

273:    Collective on Mat

275:    Input Parameters:
276: +  A - the first matrix
277: -  B - the second matrix
278: -  C - the third matrix
279: -  n - number of random vectors to be tested

281:    Output Parameter:
282: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

284:    Level: intermediate

286:    Concepts: matrices^equality between
287: @*/
288: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
289: {
291:   Vec            x,s1,s2,s3;
292:   PetscRandom    rctx;
293:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
294:   PetscInt       am,an,bm,bn,cm,cn,k;
295:   PetscScalar    none = -1.0;

298:   MatGetLocalSize(A,&am,&an);
299:   MatGetLocalSize(B,&bm,&bn);
300:   MatGetLocalSize(C,&cm,&cn);
301:   if (an != bm || am != cm || bn != cn) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %D %D %D %D",am,an,bm,bn,cm, cn);

303:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
304:   PetscRandomSetFromOptions(rctx);
305:   MatCreateVecs(B,&x,&s1);
306:   MatCreateVecs(C,NULL,&s2);
307:   VecDuplicate(s2,&s3);

309:   *flg = PETSC_TRUE;
310:   for (k=0; k<n; k++) {
311:     VecSetRandom(x,rctx);
312:     MatMult(B,x,s1);
313:     MatMult(A,s1,s2);
314:     MatMult(C,x,s3);

316:     VecNorm(s2,NORM_INFINITY,&r2);
317:     if (r2 < tol) {
318:       VecNorm(s3,NORM_INFINITY,&r1);
319:     } else {
320:       VecAXPY(s2,none,s3);
321:       VecNorm(s2,NORM_INFINITY,&r1);
322:       r1  /= r2;
323:     }
324:     if (r1 > tol) {
325:       *flg = PETSC_FALSE;
326:       PetscInfo2(A,"Error: %D-th MatMatMult() %g\n",k,(double)r1);
327:       break;
328:     }
329:   }
330:   PetscRandomDestroy(&rctx);
331:   VecDestroy(&x);
332:   VecDestroy(&s1);
333:   VecDestroy(&s2);
334:   VecDestroy(&s3);
335:   return(0);
336: }

340: /*@
341:    MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x 

343:    Collective on Mat

345:    Input Parameters:
346: +  A - the first matrix
347: -  B - the second matrix
348: -  C - the third matrix
349: -  n - number of random vectors to be tested

351:    Output Parameter:
352: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

354:    Level: intermediate

356:    Concepts: matrices^equality between
357: @*/
358: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
359: {
361:   Vec            x,s1,s2,s3;
362:   PetscRandom    rctx;
363:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
364:   PetscInt       am,an,bm,bn,cm,cn,k;
365:   PetscScalar    none = -1.0;

368:   MatGetLocalSize(A,&am,&an);
369:   MatGetLocalSize(B,&bm,&bn);
370:   MatGetLocalSize(C,&cm,&cn);
371:   if (am != bm || an != cm || bn != cn) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %D %D %D %D",am,an,bm,bn,cm, cn);

373:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
374:   PetscRandomSetFromOptions(rctx);
375:   MatCreateVecs(B,&x,&s1);
376:   MatCreateVecs(C,NULL,&s2);
377:   VecDuplicate(s2,&s3);

379:   *flg = PETSC_TRUE;
380:   for (k=0; k<n; k++) {
381:     VecSetRandom(x,rctx);
382:     MatMult(B,x,s1);
383:     MatMultTranspose(A,s1,s2);
384:     MatMult(C,x,s3);

386:     VecNorm(s2,NORM_INFINITY,&r2);
387:     if (r2 < tol) {
388:       VecNorm(s3,NORM_INFINITY,&r1);
389:     } else {
390:       VecAXPY(s2,none,s3);
391:       VecNorm(s2,NORM_INFINITY,&r1);
392:       r1  /= r2;
393:     }
394:     if (r1 > tol) {
395:       *flg = PETSC_FALSE;
396:       PetscInfo2(A,"Error: %D-th MatTransposeMatMult() %g\n",k,(double)r1);
397:       break;
398:     }
399:   }
400:   PetscRandomDestroy(&rctx);
401:   VecDestroy(&x);
402:   VecDestroy(&s1);
403:   VecDestroy(&s2);
404:   VecDestroy(&s3);
405:   return(0);
406: }