Actual source code: ex48.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: static char help[] = "Tests the vatious routines in MatSeqBAIJ format.\n";

  4: #include <petscmat.h>

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

 21:   PetscInitialize(&argc,&args,(char*)0,help);

 23:   /* Test MatSetValues() and MatGetValues() */
 24:   PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
 25:   PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
 26:   M    = m*bs;
 27:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
 28:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 29:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
 30:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 31:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 32:   PetscRandomSetFromOptions(rdm);
 33:   VecCreateSeq(PETSC_COMM_SELF,M,&xx);
 34:   VecDuplicate(xx,&s1);
 35:   VecDuplicate(xx,&s2);
 36:   VecDuplicate(xx,&yy);

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

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

 66:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

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

 91:   /* MatShift() */
 92:   rval = 10*s1norm;
 93:   MatShift(A,rval);
 94:   MatShift(B,rval);

 96:   /* Test MatTranspose() */
 97:   MatTranspose(A,MAT_REUSE_MATRIX,&A);
 98:   MatTranspose(B,MAT_REUSE_MATRIX,&B);

100:   /* Now do MatGetValues()  */
101:   for (i=0; i<30; i++) {
102:     PetscRandomGetValue(rdm,&rval);
103:     cols[0] = (PetscInt)(PetscRealPart(rval)*M);
104:     PetscRandomGetValue(rdm,&rval);
105:     cols[1] = (PetscInt)(PetscRealPart(rval)*M);
106:     PetscRandomGetValue(rdm,&rval);
107:     rows[0] = (PetscInt)(PetscRealPart(rval)*M);
108:     PetscRandomGetValue(rdm,&rval);
109:     rows[1] = (PetscInt)(PetscRealPart(rval)*M);
110:     MatGetValues(A,2,rows,2,cols,vals1);
111:     MatGetValues(B,2,rows,2,cols,vals2);
112:     PetscMemcmp(vals1,vals2,4*sizeof(PetscScalar),&flg);
113:     if (!flg) {
114:       PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %D\n",bs);
115:     }
116:   }

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

132:   /* Test MatMult() */
133:   MatMultEqual(A,B,10,&flg);
134:   if (!flg) {
135:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n");
136:   }

138:   /* Test MatMultAdd() */
139:   MatMultAddEqual(A,B,10,&flg);
140:   if (!flg) {
141:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
142:   }

144:   /* Test MatMultTranspose() */
145:   MatMultTransposeEqual(A,B,10,&flg);
146:   if (!flg) {
147:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n");
148:   }

150:   /* Test MatMultTransposeAdd() */
151:   MatMultTransposeAddEqual(A,B,10,&flg);
152:   if (!flg) {
153:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n");
154:   }

156:   /* Do LUFactor() on both the matrices */
157:   PetscMalloc1(M,&idx);
158:   for (i=0; i<M; i++) idx[i] = i;
159:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1);
160:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2);
161:   PetscFree(idx);
162:   ISSetPermutation(is1);
163:   ISSetPermutation(is2);

165:   MatFactorInfoInitialize(&info);

167:   info.fill          = 2.0;
168:   info.dtcol         = 0.0;
169:   info.zeropivot     = 1.e-14;
170:   info.pivotinblocks = 1.0;

172:   if (bs < 4) {
173:     MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact);
174:     MatLUFactorSymbolic(Fact,A,is1,is2,&info);
175:     MatLUFactorNumeric(Fact,A,&info);
176:     VecSetRandom(yy,rdm);
177:     MatForwardSolve(Fact,yy,xx);
178:     MatBackwardSolve(Fact,xx,s1);
179:     MatDestroy(&Fact);
180:     VecScale(s1,-1.0);
181:     MatMultAdd(A,s1,yy,yy);
182:     VecNorm(yy,NORM_2,&rnorm);
183:     if (rnorm<-tol || rnorm>tol) {
184:       PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %D\n",rnorm,bs);
185:     }
186:   }

188:   MatLUFactor(B,is1,is2,&info);
189:   MatLUFactor(A,is1,is2,&info);

191:   /* Test MatSolveAdd() */
192:   for (i=0; i<10; i++) {
193:     VecSetRandom(xx,rdm);
194:     VecSetRandom(yy,rdm);
195:     MatSolveAdd(B,xx,yy,s2);
196:     MatSolveAdd(A,xx,yy,s1);
197:     VecNorm(s1,NORM_2,&s1norm);
198:     VecNorm(s2,NORM_2,&s2norm);
199:     rnorm = s2norm-s1norm;
200:     if (rnorm<-tol || rnorm>tol) {
201:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
202:     }
203:   }

205:   /* Test MatSolveAdd() when x = A'b +x */
206:   for (i=0; i<10; i++) {
207:     VecSetRandom(xx,rdm);
208:     VecSetRandom(s1,rdm);
209:     VecCopy(s2,s1);
210:     MatSolveAdd(B,xx,s2,s2);
211:     MatSolveAdd(A,xx,s1,s1);
212:     VecNorm(s1,NORM_2,&s1norm);
213:     VecNorm(s2,NORM_2,&s2norm);
214:     rnorm = s2norm-s1norm;
215:     if (rnorm<-tol || rnorm>tol) {
216:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
217:     }
218:   }

220:   /* Test MatSolve() */
221:   for (i=0; i<10; i++) {
222:     VecSetRandom(xx,rdm);
223:     MatSolve(B,xx,s2);
224:     MatSolve(A,xx,s1);
225:     VecNorm(s1,NORM_2,&s1norm);
226:     VecNorm(s2,NORM_2,&s2norm);
227:     rnorm = s2norm-s1norm;
228:     if (rnorm<-tol || rnorm>tol) {
229:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
230:     }
231:   }

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

248:   MatDestroy(&A);
249:   MatDestroy(&B);
250:   VecDestroy(&xx);
251:   VecDestroy(&s1);
252:   VecDestroy(&s2);
253:   VecDestroy(&yy);
254:   ISDestroy(&is1);
255:   ISDestroy(&is2);
256:   PetscRandomDestroy(&rdm);
257:   PetscFinalize();
258:   return 0;
259: }