Actual source code: multequal.c

petsc-3.4.4 2014-03-13
  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=1.e-10;
 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);

 40: #if defined(PETSC_USE_REAL_SINGLE)
 41:   tol = 1.e-5;
 42: #endif
 43:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
 44:   PetscRandomSetFromOptions(rctx);
 45:   VecCreate(PetscObjectComm((PetscObject)A),&x);
 46:   VecSetSizes(x,an,PETSC_DECIDE);
 47:   VecSetFromOptions(x);

 49:   VecCreate(PetscObjectComm((PetscObject)A),&s1);
 50:   VecSetSizes(s1,am,PETSC_DECIDE);
 51:   VecSetFromOptions(s1);
 52:   VecDuplicate(s1,&s2);

 54:   *flg = PETSC_TRUE;
 55:   for (k=0; k<n; k++) {
 56:     VecSetRandom(x,rctx);
 57:     MatMult(A,x,s1);
 58:     MatMult(B,x,s2);
 59:     VecNorm(s2,NORM_INFINITY,&r2);
 60:     if (r2 < tol) {
 61:       VecNorm(s1,NORM_INFINITY,&r1);
 62:     } else {
 63:       VecAXPY(s2,none,s1);
 64:       VecNorm(s2,NORM_INFINITY,&r1);
 65:       r1  /= r2;
 66:     }
 67:     if (r1 > tol) {
 68:       *flg = PETSC_FALSE;
 69:       PetscInfo2(A,"Error: %D-th MatMult() %G\n",k,r1);
 70:       break;
 71:     }
 72:   }
 73:   PetscRandomDestroy(&rctx);
 74:   VecDestroy(&x);
 75:   VecDestroy(&s1);
 76:   VecDestroy(&s2);
 77:   return(0);
 78: }

 82: /*@
 83:    MatMultAddEqual - Compares matrix-vector products of two matrices.

 85:    Collective on Mat

 87:    Input Parameters:
 88: +  A - the first matrix
 89: -  B - the second matrix
 90: -  n - number of random vectors to be tested

 92:    Output Parameter:
 93: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

 95:    Level: intermediate

 97:    Concepts: matrices^equality between
 98: @*/
 99: PetscErrorCode  MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
100: {
102:   Vec            x,y,s1,s2;
103:   PetscRandom    rctx;
104:   PetscReal      r1,r2,tol=1.e-10;
105:   PetscInt       am,an,bm,bn,k;
106:   PetscScalar    none = -1.0;

109:   MatGetLocalSize(A,&am,&an);
110:   MatGetLocalSize(B,&bm,&bn);
111:   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);
113:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
114:   PetscRandomSetFromOptions(rctx);
115:   VecCreate(PetscObjectComm((PetscObject)A),&x);
116:   VecSetSizes(x,an,PETSC_DECIDE);
117:   VecSetFromOptions(x);

119:   VecCreate(PetscObjectComm((PetscObject)A),&s1);
120:   VecSetSizes(s1,am,PETSC_DECIDE);
121:   VecSetFromOptions(s1);
122:   VecDuplicate(s1,&s2);
123:   VecDuplicate(s1,&y);

125:   *flg = PETSC_TRUE;
126:   for (k=0; k<n; k++) {
127:     VecSetRandom(x,rctx);
128:     VecSetRandom(y,rctx);
129:     MatMultAdd(A,x,y,s1);
130:     MatMultAdd(B,x,y,s2);
131:     VecNorm(s2,NORM_INFINITY,&r2);
132:     if (r2 < tol) {
133:       VecNorm(s1,NORM_INFINITY,&r1);
134:     } else {
135:       VecAXPY(s2,none,s1);
136:       VecNorm(s2,NORM_INFINITY,&r1);
137:       r1  /= r2;
138:     }
139:     if (r1 > tol) {
140:       *flg = PETSC_FALSE;
141:       PetscInfo2(A,"Error: %d-th MatMultAdd() %G\n",k,r1);
142:       break;
143:     }
144:   }
145:   PetscRandomDestroy(&rctx);
146:   VecDestroy(&x);
147:   VecDestroy(&y);
148:   VecDestroy(&s1);
149:   VecDestroy(&s2);
150:   return(0);
151: }

155: /*@
156:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

158:    Collective on Mat

160:    Input Parameters:
161: +  A - the first matrix
162: -  B - the second matrix
163: -  n - number of random vectors to be tested

165:    Output Parameter:
166: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

168:    Level: intermediate

170:    Concepts: matrices^equality between
171: @*/
172: PetscErrorCode  MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
173: {
175:   Vec            x,s1,s2;
176:   PetscRandom    rctx;
177:   PetscReal      r1,r2,tol=1.e-10;
178:   PetscInt       am,an,bm,bn,k;
179:   PetscScalar    none = -1.0;

182:   MatGetLocalSize(A,&am,&an);
183:   MatGetLocalSize(B,&bm,&bn);
184:   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);
186:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
187:   PetscRandomSetFromOptions(rctx);
188:   VecCreate(PetscObjectComm((PetscObject)A),&x);
189:   VecSetSizes(x,am,PETSC_DECIDE);
190:   VecSetFromOptions(x);

192:   VecCreate(PetscObjectComm((PetscObject)A),&s1);
193:   VecSetSizes(s1,an,PETSC_DECIDE);
194:   VecSetFromOptions(s1);
195:   VecDuplicate(s1,&s2);

197:   *flg = PETSC_TRUE;
198:   for (k=0; k<n; k++) {
199:     VecSetRandom(x,rctx);
200:     MatMultTranspose(A,x,s1);
201:     MatMultTranspose(B,x,s2);
202:     VecNorm(s2,NORM_INFINITY,&r2);
203:     if (r2 < tol) {
204:       VecNorm(s1,NORM_INFINITY,&r1);
205:     } else {
206:       VecAXPY(s2,none,s1);
207:       VecNorm(s2,NORM_INFINITY,&r1);
208:       r1  /= r2;
209:     }
210:     if (r1 > tol) {
211:       *flg = PETSC_FALSE;
212:       PetscInfo2(A,"Error: %d-th MatMultTranspose() %G\n",k,r1);
213:       break;
214:     }
215:   }
216:   PetscRandomDestroy(&rctx);
217:   VecDestroy(&x);
218:   VecDestroy(&s1);
219:   VecDestroy(&s2);
220:   return(0);
221: }

225: /*@
226:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

228:    Collective on Mat

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

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

238:    Level: intermediate

240:    Concepts: matrices^equality between
241: @*/
242: PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
243: {
245:   Vec            x,y,s1,s2;
246:   PetscRandom    rctx;
247:   PetscReal      r1,r2,tol=1.e-10;
248:   PetscInt       am,an,bm,bn,k;
249:   PetscScalar    none = -1.0;

252:   MatGetLocalSize(A,&am,&an);
253:   MatGetLocalSize(B,&bm,&bn);
254:   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);
256:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
257:   PetscRandomSetFromOptions(rctx);
258:   VecCreate(PetscObjectComm((PetscObject)A),&x);
259:   VecSetSizes(x,am,PETSC_DECIDE);
260:   VecSetFromOptions(x);

262:   VecCreate(PetscObjectComm((PetscObject)A),&s1);
263:   VecSetSizes(s1,an,PETSC_DECIDE);
264:   VecSetFromOptions(s1);
265:   VecDuplicate(s1,&s2);
266:   VecDuplicate(s1,&y);

268:   *flg = PETSC_TRUE;
269:   for (k=0; k<n; k++) {
270:     VecSetRandom(x,rctx);
271:     VecSetRandom(y,rctx);
272:     MatMultTransposeAdd(A,x,y,s1);
273:     MatMultTransposeAdd(B,x,y,s2);
274:     VecNorm(s2,NORM_INFINITY,&r2);
275:     if (r2 < tol) {
276:       VecNorm(s1,NORM_INFINITY,&r1);
277:     } else {
278:       VecAXPY(s2,none,s1);
279:       VecNorm(s2,NORM_INFINITY,&r1);
280:       r1  /= r2;
281:     }
282:     if (r1 > tol) {
283:       *flg = PETSC_FALSE;
284:       PetscInfo2(A,"Error: %d-th MatMultTransposeAdd() %G\n",k,r1);
285:       break;
286:     }
287:   }
288:   PetscRandomDestroy(&rctx);
289:   VecDestroy(&x);
290:   VecDestroy(&y);
291:   VecDestroy(&s1);
292:   VecDestroy(&s2);
293:   return(0);
294: }