Actual source code: multequal.c

petsc-3.8.3 2017-12-09
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:    Concepts: matrices^equality between
 20: @*/
 21: PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
 22: {
 24:   Vec            x,s1,s2;
 25:   PetscRandom    rctx;
 26:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
 27:   PetscInt       am,an,bm,bn,k;
 28:   PetscScalar    none = -1.0;

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

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

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

 71:    Collective on Mat

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

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

 81:    Level: intermediate

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

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

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

133: /*@
134:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

136:    Collective on Mat

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

143:    Output Parameter:
144: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

146:    Level: intermediate

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

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

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

195: /*@
196:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

198:    Collective on Mat

200:    Input Parameters:
201: +  A - the first matrix
202: -  B - the second matrix
203: -  n - number of random vectors to be tested

205:    Output Parameter:
206: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

208:    Level: intermediate

210:    Concepts: matrices^equality between
211: @*/
212: PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
213: {
215:   Vec            x,y,s1,s2;
216:   PetscRandom    rctx;
217:   PetscReal      r1,r2,tol = PETSC_SQRT_MACHINE_EPSILON;
218:   PetscInt       am,an,bm,bn,k;
219:   PetscScalar    none = -1.0;

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

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

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

263:    Collective on Mat

265:    Input Parameters:
266: +  A - the first matrix
267: -  B - the second matrix
268: -  C - the third matrix
269: -  n - number of random vectors to be tested

271:    Output Parameter:
272: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

274:    Level: intermediate

276:    Concepts: matrices^equality between
277: @*/
278: PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
279: {
281:   Vec            x,s1,s2,s3;
282:   PetscRandom    rctx;
283:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
284:   PetscInt       am,an,bm,bn,cm,cn,k;
285:   PetscScalar    none = -1.0;

288:   MatGetLocalSize(A,&am,&an);
289:   MatGetLocalSize(B,&bm,&bn);
290:   MatGetLocalSize(C,&cm,&cn);
291:   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);

293:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
294:   PetscRandomSetFromOptions(rctx);
295:   MatCreateVecs(B,&x,&s1);
296:   MatCreateVecs(C,NULL,&s2);
297:   VecDuplicate(s2,&s3);

299:   *flg = PETSC_TRUE;
300:   for (k=0; k<n; k++) {
301:     VecSetRandom(x,rctx);
302:     MatMult(B,x,s1);
303:     MatMult(A,s1,s2);
304:     MatMult(C,x,s3);

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

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

331:    Collective on Mat

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

339:    Output Parameter:
340: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

342:    Level: intermediate

344:    Concepts: matrices^equality between
345: @*/
346: PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
347: {
349:   Vec            x,s1,s2,s3;
350:   PetscRandom    rctx;
351:   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
352:   PetscInt       am,an,bm,bn,cm,cn,k;
353:   PetscScalar    none = -1.0;

356:   MatGetLocalSize(A,&am,&an);
357:   MatGetLocalSize(B,&bm,&bn);
358:   MatGetLocalSize(C,&cm,&cn);
359:   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);

361:   PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);
362:   PetscRandomSetFromOptions(rctx);
363:   MatCreateVecs(B,&x,&s1);
364:   MatCreateVecs(C,NULL,&s2);
365:   VecDuplicate(s2,&s3);

367:   *flg = PETSC_TRUE;
368:   for (k=0; k<n; k++) {
369:     VecSetRandom(x,rctx);
370:     MatMult(B,x,s1);
371:     MatMultTranspose(A,s1,s2);
372:     MatMult(C,x,s3);

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