Actual source code: multequal.c

petsc-3.5.1 2014-08-06
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=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:   MatGetVecs(A,&x,&s1);
 46:   VecDuplicate(s1,&s2);

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

 76: /*@
 77:    MatMultAddEqual - Compares matrix-vector products of two matrices.

 79:    Collective on Mat

 81:    Input Parameters:
 82: +  A - the first matrix
 83: -  B - the second matrix
 84: -  n - number of random vectors to be tested

 86:    Output Parameter:
 87: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

 89:    Level: intermediate

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

103:   MatGetLocalSize(A,&am,&an);
104:   MatGetLocalSize(B,&bm,&bn);
105:   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);
107:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
108:   PetscRandomSetFromOptions(rctx);
109:   MatGetVecs(A,&x,&s1);
110:   VecDuplicate(s1,&s2);
111:   VecDuplicate(s1,&y);

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

143: /*@
144:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

146:    Collective on Mat

148:    Input Parameters:
149: +  A - the first matrix
150: -  B - the second matrix
151: -  n - number of random vectors to be tested

153:    Output Parameter:
154: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

156:    Level: intermediate

158:    Concepts: matrices^equality between
159: @*/
160: PetscErrorCode  MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
161: {
163:   Vec            x,s1,s2;
164:   PetscRandom    rctx;
165:   PetscReal      r1,r2,tol=1.e-10;
166:   PetscInt       am,an,bm,bn,k;
167:   PetscScalar    none = -1.0;

170:   MatGetLocalSize(A,&am,&an);
171:   MatGetLocalSize(B,&bm,&bn);
172:   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);
174:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
175:   PetscRandomSetFromOptions(rctx);
176:   MatGetVecs(A,&s1,&x);
177:   VecDuplicate(s1,&s2);

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

207: /*@
208:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

210:    Collective on Mat

212:    Input Parameters:
213: +  A - the first matrix
214: -  B - the second matrix
215: -  n - number of random vectors to be tested

217:    Output Parameter:
218: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

220:    Level: intermediate

222:    Concepts: matrices^equality between
223: @*/
224: PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
225: {
227:   Vec            x,y,s1,s2;
228:   PetscRandom    rctx;
229:   PetscReal      r1,r2,tol=1.e-10;
230:   PetscInt       am,an,bm,bn,k;
231:   PetscScalar    none = -1.0;

234:   MatGetLocalSize(A,&am,&an);
235:   MatGetLocalSize(B,&bm,&bn);
236:   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);
238:   PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);
239:   PetscRandomSetFromOptions(rctx);
240:   MatGetVecs(A,&s1,&x);
241:   VecDuplicate(s1,&s2);
242:   VecDuplicate(s1,&y);

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