Actual source code: multequal.c


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

  4: static PetscErrorCode MatMultEqual_Private(Mat A,Mat B,PetscInt n,PetscBool *flg,PetscBool t,PetscBool add)
  5: {
  7:   Vec            Ax = NULL,Bx = NULL,s1 = NULL,s2 = NULL,Ay = NULL, By = NULL;
  8:   PetscRandom    rctx;
  9:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
 10:   PetscInt       am,an,bm,bn,k;
 11:   PetscScalar    none = -1.0;
 12:   const char*    sops[] = {"MatMult","MatMultAdd","MatMultTranspose","MatMultTranposeAdd"};
 13:   const char*    sop;

 23:   MatGetLocalSize(A,&am,&an);
 24:   MatGetLocalSize(B,&bm,&bn);
 25:   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);
 26:   sop  = sops[(add ? 1 : 0) + 2 * (t ? 1 : 0)];
 27:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
 28:   PetscRandomSetFromOptions(rctx);
 29:   if (t) {
 30:     MatCreateVecs(A,&s1,&Ax);
 31:     MatCreateVecs(B,&s2,&Bx);
 32:   } else {
 33:     MatCreateVecs(A,&Ax,&s1);
 34:     MatCreateVecs(B,&Bx,&s2);
 35:   }
 36:   if (add) {
 37:     VecDuplicate(s1,&Ay);
 38:     VecDuplicate(s2,&By);
 39:   }

 41:   *flg = PETSC_TRUE;
 42:   for (k=0; k<n; k++) {
 43:     VecSetRandom(Ax,rctx);
 44:     VecCopy(Ax,Bx);
 45:     if (add) {
 46:       VecSetRandom(Ay,rctx);
 47:       VecCopy(Ay,By);
 48:     }
 49:     if (t) {
 50:       if (add) {
 51:         MatMultTransposeAdd(A,Ax,Ay,s1);
 52:         MatMultTransposeAdd(B,Bx,By,s2);
 53:       } else {
 54:         MatMultTranspose(A,Ax,s1);
 55:         MatMultTranspose(B,Bx,s2);
 56:       }
 57:     } else {
 58:       if (add) {
 59:         MatMultAdd(A,Ax,Ay,s1);
 60:         MatMultAdd(B,Bx,By,s2);
 61:       } else {
 62:         MatMult(A,Ax,s1);
 63:         MatMult(B,Bx,s2);
 64:       }
 65:     }
 66:     VecNorm(s2,NORM_INFINITY,&r2);
 67:     if (r2 < tol) {
 68:       VecNorm(s1,NORM_INFINITY,&r1);
 69:     } else {
 70:       VecAXPY(s2,none,s1);
 71:       VecNorm(s2,NORM_INFINITY,&r1);
 72:       r1  /= r2;
 73:     }
 74:     if (r1 > tol) {
 75:       *flg = PETSC_FALSE;
 76:       PetscInfo3(A,"Error: %D-th %s() %g\n",k,sop,(double)r1);
 77:       break;
 78:     }
 79:   }
 80:   PetscRandomDestroy(&rctx);
 81:   VecDestroy(&Ax);
 82:   VecDestroy(&Bx);
 83:   VecDestroy(&Ay);
 84:   VecDestroy(&By);
 85:   VecDestroy(&s1);
 86:   VecDestroy(&s2);
 87:   return(0);
 88: }

 90: static PetscErrorCode MatMatMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg,PetscBool At,PetscBool Bt)
 91: {
 93:   Vec            Ax,Bx,Cx,s1,s2,s3;
 94:   PetscRandom    rctx;
 95:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
 96:   PetscInt       am,an,bm,bn,cm,cn,k;
 97:   PetscScalar    none = -1.0;
 98:   const char*    sops[] = {"MatMatMult","MatTransposeMatMult","MatMatTransposeMult","MatTransposeMatTranposeMult"};
 99:   const char*    sop;

111:   MatGetLocalSize(A,&am,&an);
112:   MatGetLocalSize(B,&bm,&bn);
113:   MatGetLocalSize(C,&cm,&cn);
114:   if (At) { PetscInt tt = an; an = am; am = tt; };
115:   if (Bt) { PetscInt tt = bn; bn = bm; bm = tt; };
116:   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);

118:   sop  = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
119:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
120:   PetscRandomSetFromOptions(rctx);
121:   if (Bt) {
122:     MatCreateVecs(B,&s1,&Bx);
123:   } else {
124:     MatCreateVecs(B,&Bx,&s1);
125:   }
126:   if (At) {
127:     MatCreateVecs(A,&s2,&Ax);
128:   } else {
129:     MatCreateVecs(A,&Ax,&s2);
130:   }
131:   MatCreateVecs(C,&Cx,&s3);

133:   *flg = PETSC_TRUE;
134:   for (k=0; k<n; k++) {
135:     VecSetRandom(Bx,rctx);
136:     if (Bt) {
137:       MatMultTranspose(B,Bx,s1);
138:     } else {
139:       MatMult(B,Bx,s1);
140:     }
141:     VecCopy(s1,Ax);
142:     if (At) {
143:       MatMultTranspose(A,Ax,s2);
144:     } else {
145:       MatMult(A,Ax,s2);
146:     }
147:     VecCopy(Bx,Cx);
148:     MatMult(C,Cx,s3);

150:     VecNorm(s2,NORM_INFINITY,&r2);
151:     if (r2 < tol) {
152:       VecNorm(s3,NORM_INFINITY,&r1);
153:     } else {
154:       VecAXPY(s2,none,s3);
155:       VecNorm(s2,NORM_INFINITY,&r1);
156:       r1  /= r2;
157:     }
158:     if (r1 > tol) {
159:       *flg = PETSC_FALSE;
160:       PetscInfo3(A,"Error: %D-th %s %g\n",k,sop,(double)r1);
161:       break;
162:     }
163:   }
164:   PetscRandomDestroy(&rctx);
165:   VecDestroy(&Ax);
166:   VecDestroy(&Bx);
167:   VecDestroy(&Cx);
168:   VecDestroy(&s1);
169:   VecDestroy(&s2);
170:   VecDestroy(&s3);
171:   return(0);
172: }

174: /*@
175:    MatMultEqual - Compares matrix-vector products of two matrices.

177:    Collective on Mat

179:    Input Parameters:
180: +  A - the first matrix
181: .  B - the second matrix
182: -  n - number of random vectors to be tested

184:    Output Parameter:
185: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

187:    Level: intermediate

189: @*/
190: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
191: {

195:   MatMultEqual_Private(A,B,n,flg,PETSC_FALSE,PETSC_FALSE);
196:   return(0);
197: }

199: /*@
200:    MatMultAddEqual - Compares matrix-vector products of two matrices.

202:    Collective on Mat

204:    Input Parameters:
205: +  A - the first matrix
206: .  B - the second matrix
207: -  n - number of random vectors to be tested

209:    Output Parameter:
210: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

212:    Level: intermediate

214: @*/
215: PetscErrorCode  MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
216: {

220:   MatMultEqual_Private(A,B,n,flg,PETSC_FALSE,PETSC_TRUE);
221:   return(0);
222: }

224: /*@
225:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

227:    Collective on Mat

229:    Input Parameters:
230: +  A - the first matrix
231: .  B - the second matrix
232: -  n - number of random vectors to be tested

234:    Output Parameter:
235: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

237:    Level: intermediate

239: @*/
240: PetscErrorCode  MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
241: {

245:   MatMultEqual_Private(A,B,n,flg,PETSC_TRUE,PETSC_FALSE);
246:   return(0);
247: }

249: /*@
250:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

252:    Collective on Mat

254:    Input Parameters:
255: +  A - the first matrix
256: .  B - the second matrix
257: -  n - number of random vectors to be tested

259:    Output Parameter:
260: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

262:    Level: intermediate

264: @*/
265: PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
266: {

270:   MatMultEqual_Private(A,B,n,flg,PETSC_TRUE,PETSC_TRUE);
271:   return(0);
272: }

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

277:    Collective on Mat

279:    Input Parameters:
280: +  A - the first matrix
281: .  B - the second matrix
282: .  C - the third matrix
283: -  n - number of random vectors to be tested

285:    Output Parameter:
286: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

288:    Level: intermediate

290: @*/
291: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
292: {

296:   MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_FALSE);
297:   return(0);
298: }

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

303:    Collective on Mat

305:    Input Parameters:
306: +  A - the first matrix
307: .  B - the second matrix
308: .  C - the third matrix
309: -  n - number of random vectors to be tested

311:    Output Parameter:
312: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

314:    Level: intermediate

316: @*/
317: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
318: {

322:   MatMatMultEqual_Private(A,B,C,n,flg,PETSC_TRUE,PETSC_FALSE);
323:   return(0);
324: }

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

329:    Collective on Mat

331:    Input Parameters:
332: +  A - the first matrix
333: .  B - the second matrix
334: .  C - the third matrix
335: -  n - number of random vectors to be tested

337:    Output Parameter:
338: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

340:    Level: intermediate

342: @*/
343: PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
344: {

348:   MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_TRUE);
349:   return(0);
350: }

352: static PetscErrorCode MatProjMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool rart,PetscBool *flg)
353: {
355:   Vec            x,v1,v2,v3,v4,Cx,Bx;
356:   PetscReal      norm_abs,norm_rel,tol=PETSC_SQRT_MACHINE_EPSILON;
357:   PetscInt       i,am,an,bm,bn,cm,cn;
358:   PetscRandom    rdm;
359:   PetscScalar    none = -1.0;

362:   MatGetLocalSize(A,&am,&an);
363:   MatGetLocalSize(B,&bm,&bn);
364:   if (rart) { PetscInt t = bm; bm = bn; bn = t; }
365:   MatGetLocalSize(C,&cm,&cn);
366:   if (an != bm || bn != cm || bn != cn) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %D %D %D %D %D %D",am,an,bm,bn,cm,cn);

368:   /* Create left vector of A: v2 */
369:   MatCreateVecs(A,&Bx,&v2);

371:   /* Create right vectors of B: x, v3, v4 */
372:   if (rart) {
373:     MatCreateVecs(B,&v1,&x);
374:   } else {
375:     MatCreateVecs(B,&x,&v1);
376:   }
377:   VecDuplicate(x,&v3);

379:   MatCreateVecs(C,&Cx,&v4);
380:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
381:   PetscRandomSetFromOptions(rdm);

383:   *flg = PETSC_TRUE;
384:   for (i=0; i<n; i++) {
385:     VecSetRandom(x,rdm);
386:     VecCopy(x,Cx);
387:     MatMult(C,Cx,v4);           /* v4 = C*x   */
388:     if (rart) {
389:       MatMultTranspose(B,x,v1);
390:     } else {
391:       MatMult(B,x,v1);
392:     }
393:     VecCopy(v1,Bx);
394:     MatMult(A,Bx,v2);          /* v2 = A*B*x */
395:     VecCopy(v2,v1);
396:     if (rart) {
397:       MatMult(B,v1,v3); /* v3 = R*A*R^t*x */
398:     } else {
399:       MatMultTranspose(B,v1,v3); /* v3 = Bt*A*B*x */
400:     }
401:     VecNorm(v4,NORM_2,&norm_abs);
402:     VecAXPY(v4,none,v3);
403:     VecNorm(v4,NORM_2,&norm_rel);

405:     if (norm_abs > tol) norm_rel /= norm_abs;
406:     if (norm_rel > tol) {
407:       *flg = PETSC_FALSE;
408:       PetscInfo3(A,"Error: %D-th Mat%sMult() %g\n",i,rart ? "RARt" : "PtAP",(double)norm_rel);
409:       break;
410:     }
411:   }

413:   PetscRandomDestroy(&rdm);
414:   VecDestroy(&x);
415:   VecDestroy(&Bx);
416:   VecDestroy(&Cx);
417:   VecDestroy(&v1);
418:   VecDestroy(&v2);
419:   VecDestroy(&v3);
420:   VecDestroy(&v4);
421:   return(0);
422: }

424: /*@
425:    MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B

427:    Collective on Mat

429:    Input Parameters:
430: +  A - the first matrix
431: .  B - the second matrix
432: .  C - the third matrix
433: -  n - number of random vectors to be tested

435:    Output Parameter:
436: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

438:    Level: intermediate

440: @*/
441: PetscErrorCode MatPtAPMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
442: {

446:   MatProjMultEqual_Private(A,B,C,n,PETSC_FALSE,flg);
447:   return(0);
448: }

450: /*@
451:    MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t

453:    Collective on Mat

455:    Input Parameters:
456: +  A - the first matrix
457: .  B - the second matrix
458: .  C - the third matrix
459: -  n - number of random vectors to be tested

461:    Output Parameter:
462: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

464:    Level: intermediate

466: @*/
467: PetscErrorCode MatRARtMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
468: {

472:   MatProjMultEqual_Private(A,B,C,n,PETSC_TRUE,flg);
473:   return(0);
474: }

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

479:    Collective on Mat

481:    Input Parameters:
482: +  A - the shell matrix
483: -  n - number of random vectors to be tested

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

488:    Level: intermediate
489: @*/
490: PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool  *flg)
491: {
493:   Vec            x,y,s1,s2;
494:   PetscRandom    rctx;
495:   PetscScalar    a;
496:   PetscInt       k;
497:   PetscReal      norm,normA;
498:   MPI_Comm       comm;
499:   PetscMPIInt    rank;

503:   PetscObjectGetComm((PetscObject)A,&comm);
504:   MPI_Comm_rank(comm,&rank);

506:   PetscRandomCreate(comm,&rctx);
507:   PetscRandomSetFromOptions(rctx);
508:   MatCreateVecs(A,&x,&s1);
509:   VecDuplicate(x,&y);
510:   VecDuplicate(s1,&s2);

512:   *flg = PETSC_TRUE;
513:   for (k=0; k<n; k++) {
514:     VecSetRandom(x,rctx);
515:     VecSetRandom(y,rctx);
516:     if (!rank) {
517:       PetscRandomGetValue(rctx,&a);
518:     }
519:     MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);

521:     /* s2 = a*A*x + A*y */
522:     MatMult(A,y,s2); /* s2 = A*y */
523:     MatMult(A,x,s1); /* s1 = A*x */
524:     VecAXPY(s2,a,s1); /* s2 = a s1 + s2 */

526:     /* s1 = A * (a x + y) */
527:     VecAXPY(y,a,x); /* y = a x + y */
528:     MatMult(A,y,s1);
529:     VecNorm(s1,NORM_INFINITY,&normA);

531:     VecAXPY(s2,-1.0,s1); /* s2 = - s1 + s2 */
532:     VecNorm(s2,NORM_INFINITY,&norm);
533:     if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
534:       *flg = PETSC_FALSE;
535:       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);
536:       break;
537:     }
538:   }
539:   PetscRandomDestroy(&rctx);
540:   VecDestroy(&x);
541:   VecDestroy(&y);
542:   VecDestroy(&s1);
543:   VecDestroy(&s2);
544:   return(0);
545: }