Actual source code: multequal.c

petsc-master 2019-09-20
Report Typos and Errors

  2:  #include <petsc/private/matimpl.h>

  4: /*@
  5:    MatMultEqual - Compares matrix-vector products of two matrices.

  7:    Collective on Mat

  9:    Input Parameters:
 10: +  A - the first matrix
 11: -  B - the second matrix
 12: -  n - number of random vectors to be tested

 14:    Output Parameter:
 15: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

 17:    Level: intermediate

 19: @*/
 20: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
 21: {
 23:   Vec            x,s1,s2;
 24:   PetscRandom    rctx;
 25:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
 26:   PetscInt       am,an,bm,bn,k;
 27:   PetscScalar    none = -1.0;

 32:   MatGetLocalSize(A,&am,&an);
 33:   MatGetLocalSize(B,&bm,&bn);
 34:   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);
 36:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
 37:   PetscRandomSetFromOptions(rctx);
 38:   MatCreateVecs(A,&x,&s1);
 39:   VecDuplicate(s1,&s2);

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

 67: /*@
 68:    MatMultAddEqual - Compares matrix-vector products of two matrices.

 70:    Collective on Mat

 72:    Input Parameters:
 73: +  A - the first matrix
 74: -  B - the second matrix
 75: -  n - number of random vectors to be tested

 77:    Output Parameter:
 78: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

 80:    Level: intermediate

 82: @*/
 83: PetscErrorCode  MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
 84: {
 86:   Vec            x,y,s1,s2;
 87:   PetscRandom    rctx;
 88:   PetscReal      r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
 89:   PetscInt       am,an,bm,bn,k;
 90:   PetscScalar    none = -1.0;

 93:   MatGetLocalSize(A,&am,&an);
 94:   MatGetLocalSize(B,&bm,&bn);
 95:   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);
 97:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
 98:   PetscRandomSetFromOptions(rctx);
 99:   MatCreateVecs(A,&x,&s1);
100:   VecDuplicate(s1,&s2);
101:   VecDuplicate(s1,&y);

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

131: /*@
132:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

134:    Collective on Mat

136:    Input Parameters:
137: +  A - the first matrix
138: -  B - the second matrix
139: -  n - number of random vectors to be tested

141:    Output Parameter:
142: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

144:    Level: intermediate

146: @*/
147: PetscErrorCode  MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
148: {
150:   Vec            x,s1,s2;
151:   PetscRandom    rctx;
152:   PetscReal      r1,r2,tol= PETSC_SQRT_MACHINE_EPSILON;
153:   PetscInt       am,an,bm,bn,k;
154:   PetscScalar    none = -1.0;

157:   MatGetLocalSize(A,&am,&an);
158:   MatGetLocalSize(B,&bm,&bn);
159:   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);
161:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
162:   PetscRandomSetFromOptions(rctx);
163:   MatCreateVecs(A,&s1,&x);
164:   VecDuplicate(s1,&s2);

166:   *flg = PETSC_TRUE;
167:   for (k=0; k<n; k++) {
168:     VecSetRandom(x,rctx);
169:     MatMultTranspose(A,x,s1);
170:     MatMultTranspose(B,x,s2);
171:     VecNorm(s2,NORM_INFINITY,&r2);
172:     if (r2 < tol) {
173:       VecNorm(s1,NORM_INFINITY,&r1);
174:     } else {
175:       VecAXPY(s2,none,s1);
176:       VecNorm(s2,NORM_INFINITY,&r1);
177:       r1  /= r2;
178:     }
179:     if (r1 > tol) {
180:       *flg = PETSC_FALSE;
181:       PetscInfo2(A,"Error: %d-th MatMultTranspose() %g\n",k,(double)r1);
182:       break;
183:     }
184:   }
185:   PetscRandomDestroy(&rctx);
186:   VecDestroy(&x);
187:   VecDestroy(&s1);
188:   VecDestroy(&s2);
189:   return(0);
190: }

192: /*@
193:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

195:    Collective on Mat

197:    Input Parameters:
198: +  A - the first matrix
199: -  B - the second matrix
200: -  n - number of random vectors to be tested

202:    Output Parameter:
203: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

205:    Level: intermediate

207: @*/
208: PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
209: {
211:   Vec            x,y,s1,s2;
212:   PetscRandom    rctx;
213:   PetscReal      r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
214:   PetscInt       am,an,bm,bn,k;
215:   PetscScalar    none = -1.0;

218:   MatGetLocalSize(A,&am,&an);
219:   MatGetLocalSize(B,&bm,&bn);
220:   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);
222:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
223:   PetscRandomSetFromOptions(rctx);
224:   MatCreateVecs(A,&s1,&x);
225:   VecDuplicate(s1,&s2);
226:   VecDuplicate(s1,&y);

228:   *flg = PETSC_TRUE;
229:   for (k=0; k<n; k++) {
230:     VecSetRandom(x,rctx);
231:     VecSetRandom(y,rctx);
232:     MatMultTransposeAdd(A,x,y,s1);
233:     MatMultTransposeAdd(B,x,y,s2);
234:     VecNorm(s2,NORM_INFINITY,&r2);
235:     if (r2 < tol) {
236:       VecNorm(s1,NORM_INFINITY,&r1);
237:     } else {
238:       VecAXPY(s2,none,s1);
239:       VecNorm(s2,NORM_INFINITY,&r1);
240:       r1  /= r2;
241:     }
242:     if (r1 > tol) {
243:       *flg = PETSC_FALSE;
244:       PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %g\n",k,(double)r1);
245:       break;
246:     }
247:   }
248:   PetscRandomDestroy(&rctx);
249:   VecDestroy(&x);
250:   VecDestroy(&y);
251:   VecDestroy(&s1);
252:   VecDestroy(&s2);
253:   return(0);
254: }

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

259:    Collective on Mat

261:    Input Parameters:
262: +  A - the first matrix
263: -  B - the second matrix
264: -  C - the third matrix
265: -  n - number of random vectors to be tested

267:    Output Parameter:
268: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

270:    Level: intermediate

272: @*/
273: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
274: {
276:   Vec            x,s1,s2,s3;
277:   PetscRandom    rctx;
278:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
279:   PetscInt       am,an,bm,bn,cm,cn,k;
280:   PetscScalar    none = -1.0;

283:   MatGetLocalSize(A,&am,&an);
284:   MatGetLocalSize(B,&bm,&bn);
285:   MatGetLocalSize(C,&cm,&cn);
286:   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);

288:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
289:   PetscRandomSetFromOptions(rctx);
290:   MatCreateVecs(B,&x,&s1);
291:   MatCreateVecs(C,NULL,&s2);
292:   VecDuplicate(s2,&s3);

294:   *flg = PETSC_TRUE;
295:   for (k=0; k<n; k++) {
296:     VecSetRandom(x,rctx);
297:     MatMult(B,x,s1);
298:     MatMult(A,s1,s2);
299:     MatMult(C,x,s3);

301:     VecNorm(s2,NORM_INFINITY,&r2);
302:     if (r2 < tol) {
303:       VecNorm(s3,NORM_INFINITY,&r1);
304:     } else {
305:       VecAXPY(s2,none,s3);
306:       VecNorm(s2,NORM_INFINITY,&r1);
307:       r1  /= r2;
308:     }
309:     if (r1 > tol) {
310:       *flg = PETSC_FALSE;
311:       PetscInfo2(A,"Error: %D-th MatMatMult() %g\n",k,(double)r1);
312:       break;
313:     }
314:   }
315:   PetscRandomDestroy(&rctx);
316:   VecDestroy(&x);
317:   VecDestroy(&s1);
318:   VecDestroy(&s2);
319:   VecDestroy(&s3);
320:   return(0);
321: }

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

326:    Collective on Mat

328:    Input Parameters:
329: +  A - the first matrix
330: -  B - the second matrix
331: -  C - the third matrix
332: -  n - number of random vectors to be tested

334:    Output Parameter:
335: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

337:    Level: intermediate

339: @*/
340: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
341: {
343:   Vec            x,s1,s2,s3;
344:   PetscRandom    rctx;
345:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
346:   PetscInt       am,an,bm,bn,cm,cn,k;
347:   PetscScalar    none = -1.0;

350:   MatGetLocalSize(A,&am,&an);
351:   MatGetLocalSize(B,&bm,&bn);
352:   MatGetLocalSize(C,&cm,&cn);
353:   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);

355:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
356:   PetscRandomSetFromOptions(rctx);
357:   MatCreateVecs(B,&x,&s1);
358:   MatCreateVecs(C,NULL,&s2);
359:   VecDuplicate(s2,&s3);

361:   *flg = PETSC_TRUE;
362:   for (k=0; k<n; k++) {
363:     VecSetRandom(x,rctx);
364:     MatMult(B,x,s1);
365:     MatMultTranspose(A,s1,s2);
366:     MatMult(C,x,s3);

368:     VecNorm(s2,NORM_INFINITY,&r2);
369:     if (r2 < tol) {
370:       VecNorm(s3,NORM_INFINITY,&r1);
371:     } else {
372:       VecAXPY(s2,none,s3);
373:       VecNorm(s2,NORM_INFINITY,&r1);
374:       r1  /= r2;
375:     }
376:     if (r1 > tol) {
377:       *flg = PETSC_FALSE;
378:       PetscInfo2(A,"Error: %D-th MatTransposeMatMult() %g\n",k,(double)r1);
379:       break;
380:     }
381:   }
382:   PetscRandomDestroy(&rctx);
383:   VecDestroy(&x);
384:   VecDestroy(&s1);
385:   VecDestroy(&s2);
386:   VecDestroy(&s3);
387:   return(0);
388: }

390: /*@
391:    MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x

393:    Collective on Mat

395:    Input Parameters:
396: +  A - the first matrix
397: -  B - the second matrix
398: -  C - the third matrix
399: -  n - number of random vectors to be tested

401:    Output Parameter:
402: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

404:    Level: intermediate

406: @*/
407: PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
408: {
410:   Vec            x,s1,s2,s3;
411:   PetscRandom    rctx;
412:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
413:   PetscInt       am,an,bm,bn,cm,cn,k;
414:   PetscScalar    none = -1.0;

417:   MatGetLocalSize(A,&am,&an);
418:   MatGetLocalSize(B,&bm,&bn);
419:   MatGetLocalSize(C,&cm,&cn);
420:   if (an != bn || am != cm || bm != 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);

422:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
423:   PetscRandomSetFromOptions(rctx);
424:   MatCreateVecs(B,&s1,&x);
425:   MatCreateVecs(C,NULL,&s2);
426:   VecDuplicate(s2,&s3);

428:   *flg = PETSC_TRUE;
429:   for (k=0; k<n; k++) {
430:     VecSetRandom(x,rctx);
431:     MatMultTranspose(B,x,s1);
432:     MatMult(A,s1,s2);
433:     MatMult(C,x,s3);

435:     VecNorm(s2,NORM_INFINITY,&r2);
436:     if (r2 < tol) {
437:       VecNorm(s3,NORM_INFINITY,&r1);
438:     } else {
439:       VecAXPY(s2,none,s3);
440:       VecNorm(s2,NORM_INFINITY,&r1);
441:       r1  /= r2;
442:     }
443:     if (r1 > tol) {
444:       *flg = PETSC_FALSE;
445:       PetscInfo2(A,"Error: %D-th MatMatTransposeMult() %g\n",k,(double)r1);
446:       break;
447:     }
448:   }
449:   PetscRandomDestroy(&rctx);
450:   VecDestroy(&x);
451:   VecDestroy(&s1);
452:   VecDestroy(&s2);
453:   VecDestroy(&s3);
454:   return(0);
455: }

457: /*@
458:    MatIsLinear - Check if a shell matrix A is a linear operator.

460:    Collective on Mat

462:    Input Parameters:
463: +  A - the shell matrix
464: -  n - number of random vectors to be tested

466:    Output Parameter:
467: .  flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.

469:    Level: intermediate
470: @*/
471: PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool  *flg)
472: {
474:   Vec            x,y,s1,s2;
475:   PetscRandom    rctx;
476:   PetscScalar    a;
477:   PetscInt       k;
478:   PetscReal      norm,normA;
479:   MPI_Comm       comm;
480:   PetscMPIInt    rank;

484:   PetscObjectGetComm((PetscObject)A,&comm);
485:   MPI_Comm_rank(comm,&rank);

487:   PetscRandomCreate(comm,&rctx);
488:   PetscRandomSetFromOptions(rctx);
489:   MatCreateVecs(A,&x,&s1);
490:   VecDuplicate(x,&y);
491:   VecDuplicate(s1,&s2);

493:   *flg = PETSC_TRUE;
494:   for (k=0; k<n; k++) {
495:     VecSetRandom(x,rctx);
496:     VecSetRandom(y,rctx);
497:     if (!rank) {
498:       PetscRandomGetValue(rctx,&a);
499:     }
500:     MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);

502:     /* s2 = a*A*x + A*y */
503:     MatMult(A,y,s2); /* s2 = A*y */
504:     MatMult(A,x,s1); /* s1 = A*x */
505:     VecAXPY(s2,a,s1); /* s2 = a s1 + s2 */

507:     /* s1 = A * (a x + y) */
508:     VecAXPY(y,a,x); /* y = a x + y */
509:     MatMult(A,y,s1);
510:     VecNorm(s1,NORM_INFINITY,&normA);

512:     VecAXPY(s2,-1.0,s1); /* s2 = - s1 + s2 */
513:     VecNorm(s2,NORM_INFINITY,&norm);
514:     if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
515:       *flg = PETSC_FALSE;
516:       PetscInfo3(A,"Error: %D-th |A*(ax+y) - (a*A*x+A*y)|/|A(ax+y)| %g > tol %g\n",k,(double)norm/normA,100.*PETSC_MACHINE_EPSILON);
517:       break;
518:     }
519:   }
520:   PetscRandomDestroy(&rctx);
521:   VecDestroy(&x);
522:   VecDestroy(&y);
523:   VecDestroy(&s1);
524:   VecDestroy(&s2);
525:   return(0);
526: }