Actual source code: ex48.c

petsc-3.3-p6 2013-02-11
  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;
 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;
 20: 
 21:   PetscInitialize(&argc,&args,(char *)0,help);
 22: 
 23:   /* Test MatSetValues() and MatGetValues() */
 24:   PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
 25:   PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
 26:   M    = m*bs;
 27:   MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
 28:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 29:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,PETSC_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);
 37: 
 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:   }
 47: 
 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:   }
 65: 
 66:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 70: 
 71:   /* Test MatNorm() */
 72:   MatNorm(A,NORM_FROBENIUS,&s1norm);
 73:   MatNorm(B,NORM_FROBENIUS,&s2norm);
 74:   rnorm = PetscAbsScalar(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 = PetscAbsScalar(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 = PetscAbsScalar(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:   }
117: 
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:   }
137: 
138:   /* Test MatMultAdd() */
139:   MatMultAddEqual(A,B,10,&flg);
140:   if (!flg){
141:     PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
142:   }
143: 
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:   PetscMalloc(M*sizeof(PetscInt),&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);
166:   info.fill      = 2.0;
167:   info.dtcol     = 0.0;
168:   info.zeropivot = 1.e-14;
169:   info.pivotinblocks = 1.0;
170:   MatLUFactor(B,is1,is2,&info);
171:   MatLUFactor(A,is1,is2,&info);
172: 
173:   /* Test MatSolveAdd() */
174:   for (i=0; i<10; i++) {
175:     VecSetRandom(xx,rdm);
176:     VecSetRandom(yy,rdm);
177:     MatSolveAdd(B,xx,yy,s2);
178:     MatSolveAdd(A,xx,yy,s1);
179:     VecNorm(s1,NORM_2,&s1norm);
180:     VecNorm(s2,NORM_2,&s2norm);
181:     rnorm = s2norm-s1norm;
182:     if (rnorm<-tol || rnorm>tol) {
183:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
184:     }
185:   }
186: 
187:   /* Test MatSolveAdd() when x = A'b +x */
188:   for (i=0; i<10; i++) {
189:     VecSetRandom(xx,rdm);
190:     VecSetRandom(s1,rdm);
191:     VecCopy(s2,s1);
192:     MatSolveAdd(B,xx,s2,s2);
193:     MatSolveAdd(A,xx,s1,s1);
194:     VecNorm(s1,NORM_2,&s1norm);
195:     VecNorm(s2,NORM_2,&s2norm);
196:     rnorm = s2norm-s1norm;
197:     if (rnorm<-tol || rnorm>tol) {
198:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
199:     }
200:   }
201: 
202:   /* Test MatSolve() */
203:   for (i=0; i<10; i++) {
204:     VecSetRandom(xx,rdm);
205:     MatSolve(B,xx,s2);
206:     MatSolve(A,xx,s1);
207:     VecNorm(s1,NORM_2,&s1norm);
208:     VecNorm(s2,NORM_2,&s2norm);
209:     rnorm = s2norm-s1norm;
210:     if (rnorm<-tol || rnorm>tol) {
211:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
212:     }
213:   }
214: 
215:   /* Test MatSolveTranspose() */
216:   if (bs < 8) {
217:     for (i=0; i<10; i++) {
218:       VecSetRandom(xx,rdm);
219:       MatSolveTranspose(B,xx,s2);
220:       MatSolveTranspose(A,xx,s1);
221:       VecNorm(s1,NORM_2,&s1norm);
222:       VecNorm(s2,NORM_2,&s2norm);
223:       rnorm = s2norm-s1norm;
224:       if (rnorm<-tol || rnorm>tol) {
225:         PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
226:       }
227:     }
228:   }

230:   MatDestroy(&A);
231:   MatDestroy(&B);
232:   VecDestroy(&xx);
233:   VecDestroy(&s1);
234:   VecDestroy(&s2);
235:   VecDestroy(&yy);
236:   ISDestroy(&is1);
237:   ISDestroy(&is2);
238:   PetscRandomDestroy(&rdm);
239:   PetscFinalize();
240:   return 0;
241: }