Actual source code: ex48.c


  2: static char help[] = "Tests various routines in MatSeqBAIJ format.\n";

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            A,B,C,D,Fact;
  9:   Vec            xx,s1,s2,yy;
 10:   PetscInt       m=45,rows[2],cols[2],bs=1,i,row,col,*idx,M;
 11:   PetscScalar    rval,vals1[4],vals2[4];
 12:   PetscRandom    rdm;
 13:   IS             is1,is2;
 14:   PetscReal      s1norm,s2norm,rnorm,tol = 1.e-4;
 15:   PetscBool      flg;
 16:   MatFactorInfo  info;

 18:   PetscInitialize(&argc,&args,(char*)0,help);
 19:   /* Test MatSetValues() and MatGetValues() */
 20:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 21:   PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
 22:   M    = m*bs;
 23:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
 24:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 25:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
 26:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 27:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 28:   PetscRandomSetFromOptions(rdm);
 29:   VecCreateSeq(PETSC_COMM_SELF,M,&xx);
 30:   VecDuplicate(xx,&s1);
 31:   VecDuplicate(xx,&s2);
 32:   VecDuplicate(xx,&yy);

 34:   /* For each row add atleast 15 elements */
 35:   for (row=0; row<M; row++) {
 36:     for (i=0; i<25*bs; i++) {
 37:       PetscRandomGetValue(rdm,&rval);
 38:       col  = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
 39:       MatSetValues(A,1,&row,1,&col,&rval,INSERT_VALUES);
 40:       MatSetValues(B,1,&row,1,&col,&rval,INSERT_VALUES);
 41:     }
 42:   }

 44:   /* Now set blocks of values */
 45:   for (i=0; i<20*bs; i++) {
 46:     PetscRandomGetValue(rdm,&rval);
 47:     cols[0]  = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
 48:     vals1[0] = rval;
 49:     PetscRandomGetValue(rdm,&rval);
 50:     cols[1]  = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
 51:     vals1[1] = rval;
 52:     PetscRandomGetValue(rdm,&rval);
 53:     rows[0]  = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
 54:     vals1[2] = rval;
 55:     PetscRandomGetValue(rdm,&rval);
 56:     rows[1]  = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
 57:     vals1[3] = rval;
 58:     MatSetValues(A,2,rows,2,cols,vals1,INSERT_VALUES);
 59:     MatSetValues(B,2,rows,2,cols,vals1,INSERT_VALUES);
 60:   }

 62:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 63:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 64:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 65:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 67:   /* Test MatNorm() */
 68:   MatNorm(A,NORM_FROBENIUS,&s1norm);
 69:   MatNorm(B,NORM_FROBENIUS,&s2norm);
 70:   rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
 71:   if (rnorm>tol) {
 72:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
 73:   }
 74:   MatNorm(A,NORM_INFINITY,&s1norm);
 75:   MatNorm(B,NORM_INFINITY,&s2norm);
 76:   rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
 77:   if (rnorm>tol) {
 78:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
 79:   }
 80:   MatNorm(A,NORM_1,&s1norm);
 81:   MatNorm(B,NORM_1,&s2norm);
 82:   rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
 83:   if (rnorm>tol) {
 84:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
 85:   }

 87:   /* MatShift() */
 88:   rval = 10*s1norm;
 89:   MatShift(A,rval);
 90:   MatShift(B,rval);

 92:   /* Test MatTranspose() */
 93:   MatTranspose(A,MAT_INPLACE_MATRIX,&A);
 94:   MatTranspose(B,MAT_INPLACE_MATRIX,&B);

 96:   /* Now do MatGetValues()  */
 97:   for (i=0; i<30; i++) {
 98:     PetscRandomGetValue(rdm,&rval);
 99:     cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
100:     PetscRandomGetValue(rdm,&rval);
101:     cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
102:     PetscRandomGetValue(rdm,&rval);
103:     rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
104:     PetscRandomGetValue(rdm,&rval);
105:     rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
106:     MatGetValues(A,2,rows,2,cols,vals1);
107:     MatGetValues(B,2,rows,2,cols,vals2);
108:     PetscArraycmp(vals1,vals2,4,&flg);
109:     if (!flg) {
110:       PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %" PetscInt_FMT "\n",bs);
111:     }
112:   }

114:   /* Test MatMult(), MatMultAdd() */
115:   for (i=0; i<40; i++) {
116:     VecSetRandom(xx,rdm);
117:     VecSet(s2,0.0);
118:     MatMult(A,xx,s1);
119:     MatMultAdd(A,xx,s2,s2);
120:     VecNorm(s1,NORM_2,&s1norm);
121:     VecNorm(s2,NORM_2,&s2norm);
122:     rnorm = s2norm-s1norm;
123:     if (rnorm<-tol || rnorm>tol) {
124:       PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
125:     }
126:   }

128:   /* Test MatMult() */
129:   MatMultEqual(A,B,10,&flg);
130:   if (!flg) {
131:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n");
132:   }

134:   /* Test MatMultAdd() */
135:   MatMultAddEqual(A,B,10,&flg);
136:   if (!flg) {
137:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
138:   }

140:   /* Test MatMultTranspose() */
141:   MatMultTransposeEqual(A,B,10,&flg);
142:   if (!flg) {
143:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n");
144:   }

146:   /* Test MatMultTransposeAdd() */
147:   MatMultTransposeAddEqual(A,B,10,&flg);
148:   if (!flg) {
149:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n");
150:   }

152:   /* Test MatMatMult() */
153:   MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&C);
154:   MatSetRandom(C,rdm);
155:   MatMatMult(A,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
156:   MatMatMultEqual(A,C,D,40,&flg);
157:   if (!flg) {
158:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
159:   }
160:   MatDestroy(&D);
161:   MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&D);
162:   MatMatMult(A,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&D);
163:   MatMatMultEqual(A,C,D,40,&flg);
164:   if (!flg) {
165:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
166:   }

168:   /* Do LUFactor() on both the matrices */
169:   PetscMalloc1(M,&idx);
170:   for (i=0; i<M; i++) idx[i] = i;
171:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1);
172:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2);
173:   PetscFree(idx);
174:   ISSetPermutation(is1);
175:   ISSetPermutation(is2);

177:   MatFactorInfoInitialize(&info);

179:   info.fill          = 2.0;
180:   info.dtcol         = 0.0;
181:   info.zeropivot     = 1.e-14;
182:   info.pivotinblocks = 1.0;

184:   if (bs < 4) {
185:     MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact);
186:     MatLUFactorSymbolic(Fact,A,is1,is2,&info);
187:     MatLUFactorNumeric(Fact,A,&info);
188:     VecSetRandom(yy,rdm);
189:     MatForwardSolve(Fact,yy,xx);
190:     MatBackwardSolve(Fact,xx,s1);
191:     MatDestroy(&Fact);
192:     VecScale(s1,-1.0);
193:     MatMultAdd(A,s1,yy,yy);
194:     VecNorm(yy,NORM_2,&rnorm);
195:     if (rnorm<-tol || rnorm>tol) {
196:       PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n",(double)rnorm,bs);
197:     }
198:   }

200:   MatLUFactor(B,is1,is2,&info);
201:   MatLUFactor(A,is1,is2,&info);

203:   /* Test MatSolveAdd() */
204:   for (i=0; i<10; i++) {
205:     VecSetRandom(xx,rdm);
206:     VecSetRandom(yy,rdm);
207:     MatSolveAdd(B,xx,yy,s2);
208:     MatSolveAdd(A,xx,yy,s1);
209:     VecNorm(s1,NORM_2,&s1norm);
210:     VecNorm(s2,NORM_2,&s2norm);
211:     rnorm = s2norm-s1norm;
212:     if (rnorm<-tol || rnorm>tol) {
213:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
214:     }
215:   }

217:   /* Test MatSolveAdd() when x = A'b +x */
218:   for (i=0; i<10; i++) {
219:     VecSetRandom(xx,rdm);
220:     VecSetRandom(s1,rdm);
221:     VecCopy(s2,s1);
222:     MatSolveAdd(B,xx,s2,s2);
223:     MatSolveAdd(A,xx,s1,s1);
224:     VecNorm(s1,NORM_2,&s1norm);
225:     VecNorm(s2,NORM_2,&s2norm);
226:     rnorm = s2norm-s1norm;
227:     if (rnorm<-tol || rnorm>tol) {
228:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
229:     }
230:   }

232:   /* Test MatSolve() */
233:   for (i=0; i<10; i++) {
234:     VecSetRandom(xx,rdm);
235:     MatSolve(B,xx,s2);
236:     MatSolve(A,xx,s1);
237:     VecNorm(s1,NORM_2,&s1norm);
238:     VecNorm(s2,NORM_2,&s2norm);
239:     rnorm = s2norm-s1norm;
240:     if (rnorm<-tol || rnorm>tol) {
241:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
242:     }
243:   }

245:   /* Test MatSolveTranspose() */
246:   if (bs < 8) {
247:     for (i=0; i<10; i++) {
248:       VecSetRandom(xx,rdm);
249:       MatSolveTranspose(B,xx,s2);
250:       MatSolveTranspose(A,xx,s1);
251:       VecNorm(s1,NORM_2,&s1norm);
252:       VecNorm(s2,NORM_2,&s2norm);
253:       rnorm = s2norm-s1norm;
254:       if (rnorm<-tol || rnorm>tol) {
255:         PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
256:       }
257:     }
258:   }

260:   MatDestroy(&A);
261:   MatDestroy(&B);
262:   MatDestroy(&C);
263:   MatDestroy(&D);
264:   VecDestroy(&xx);
265:   VecDestroy(&s1);
266:   VecDestroy(&s2);
267:   VecDestroy(&yy);
268:   ISDestroy(&is1);
269:   ISDestroy(&is2);
270:   PetscRandomDestroy(&rdm);
271:   PetscFinalize();
272:   return 0;
273: }

275: /*TEST

277:    test:
278:       args: -mat_block_size {{1 2 3 4 5 6 7 8}}

280: TEST*/