Actual source code: ex48.c

  1: /*$Id: ex48.c,v 1.25 2001/08/07 21:30:08 bsmith Exp $*/

  3: static char help[] = "Tests the vatious routines in MatBAIJ format.\n";

 5:  #include petscmat.h

  9: int main(int argc,char **args)
 10: {
 11:   Mat         A,B;
 12:   Vec         xx,s1,s2,yy;
 13:   int         m=45,ierr,rows[2],cols[2],bs=1,i,row,col,*idx,M;
 14:   PetscScalar rval,vals1[4],vals2[4],zero=0.0;
 15:   PetscRandom rand;
 16:   IS          is1,is2;
 17:   PetscReal   s1norm,s2norm,rnorm,tol = 1.e-10;
 18:   PetscTruth  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:   MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,PETSC_NULL,&B);
 29:   PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);
 30:   VecCreateSeq(PETSC_COMM_SELF,M,&xx);
 31:   VecDuplicate(xx,&s1);
 32:   VecDuplicate(xx,&s2);
 33:   VecDuplicate(xx,&yy);
 34: 
 35:   /* For each row add atleast 15 elements */
 36:   for (row=0; row<M; row++) {
 37:     for (i=0; i<25*bs; i++) {
 38:       PetscRandomGetValue(rand,&rval);
 39:       col  = (int)(PetscRealPart(rval)*M);
 40:       MatSetValues(A,1,&row,1,&col,&rval,ADD_VALUES);
 41:       MatSetValues(B,1,&row,1,&col,&rval,ADD_VALUES);
 42:     }
 43:   }
 44: 
 45:   /* Now set blocks of values */
 46:   for (i=0; i<20*bs; i++) {
 47:     PetscRandomGetValue(rand,&rval);
 48:     cols[0] = (int)(PetscRealPart(rval)*M);
 49:     vals1[0] = rval;
 50:     PetscRandomGetValue(rand,&rval);
 51:     cols[1] = (int)(PetscRealPart(rval)*M);
 52:     vals1[1] = rval;
 53:     PetscRandomGetValue(rand,&rval);
 54:     rows[0] = (int)(PetscRealPart(rval)*M);
 55:     vals1[2] = rval;
 56:     PetscRandomGetValue(rand,&rval);
 57:     rows[1] = (int)(PetscRealPart(rval)*M);
 58:     vals1[3] = rval;
 59:     MatSetValues(A,2,rows,2,cols,vals1,ADD_VALUES);
 60:     MatSetValues(B,2,rows,2,cols,vals1,ADD_VALUES);
 61:   }
 62: 
 63:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 64:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 65:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 67: 
 68:   /* Test MatNorm() */
 69:   MatNorm(A,NORM_FROBENIUS,&s1norm);
 70:   MatNorm(B,NORM_FROBENIUS,&s2norm);
 71:   rnorm = s2norm-s1norm;
 72:   if (rnorm<-tol || rnorm>tol) {
 73:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm()- Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
 74:   }
 75:   /* MatScale() */
 76:   rval = 10*s1norm;
 77:   MatShift(&rval,A);
 78:   MatShift(&rval,B);
 79: 
 80:   /* Test MatTranspose() */
 81:   MatTranspose(A,PETSC_NULL);
 82:   MatTranspose(B,PETSC_NULL);
 83: 
 84: 
 85:   /* Now do MatGetValues()  */
 86:   for (i=0; i<30; i++) {
 87:     PetscRandomGetValue(rand,&rval);
 88:     cols[0] = (int)(PetscRealPart(rval)*M);
 89:     PetscRandomGetValue(rand,&rval);
 90:     cols[1] = (int)(PetscRealPart(rval)*M);
 91:     PetscRandomGetValue(rand,&rval);
 92:     rows[0] = (int)(PetscRealPart(rval)*M);
 93:     PetscRandomGetValue(rand,&rval);
 94:     rows[1] = (int)(PetscRealPart(rval)*M);
 95:     MatGetValues(A,2,rows,2,cols,vals1);
 96:     MatGetValues(B,2,rows,2,cols,vals2);
 97:     PetscMemcmp(vals1,vals2,4*sizeof(PetscScalar),&flg);
 98:     if (!flg) {
 99:       PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %d\n",bs);
100:     }
101:   }
102: 
103:   /* Test MatMult(), MatMultAdd() */
104:   for (i=0; i<40; i++) {
105:     VecSetRandom(rand,xx);
106:     VecSet(&zero,s2);
107:     MatMult(A,xx,s1);
108:     MatMultAdd(A,xx,s2,s2);
109:     VecNorm(s1,NORM_2,&s1norm);
110:     VecNorm(s2,NORM_2,&s2norm);
111:     rnorm = s2norm-s1norm;
112:     if (rnorm<-tol || rnorm>tol) {
113:       PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %d\n",s1norm,s2norm,bs);
114:     }
115:   }

117:   /* Test MatMult() */
118:   for (i=0; i<40; i++) {
119:     VecSetRandom(rand,xx);
120:     MatMult(A,xx,s1);
121:     MatMult(B,xx,s2);
122:     VecNorm(s1,NORM_2,&s1norm);
123:     VecNorm(s2,NORM_2,&s2norm);
124:     rnorm = s2norm-s1norm;
125:     if (rnorm<-tol || rnorm>tol) {
126:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d \n",s1norm,s2norm,bs);
127:     }
128:   }
129: 
130:   /* Test MatMultAdd() */
131:   for (i=0; i<40; i++) {
132:     VecSetRandom(rand,xx);
133:     VecSetRandom(rand,yy);
134:     MatMultAdd(A,xx,yy,s1);
135:     MatMultAdd(B,xx,yy,s2);
136:     VecNorm(s1,NORM_2,&s1norm);
137:     VecNorm(s2,NORM_2,&s2norm);
138:     rnorm = s2norm-s1norm;
139:     if (rnorm<-tol || rnorm>tol) {
140:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
141:     }
142:   }
143: 
144:   /* Test MatMultTranspose() */
145:   for (i=0; i<40; i++) {
146:     VecSetRandom(rand,xx);
147:     MatMultTranspose(A,xx,s1);
148:     MatMultTranspose(B,xx,s2);
149:     VecNorm(s1,NORM_2,&s1norm);
150:     VecNorm(s2,NORM_2,&s2norm);
151:     rnorm = s2norm-s1norm;
152:     if (rnorm<-tol || rnorm>tol) {
153:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
154:     }
155:   }
156:   /* Test MatMultTransposeAdd() */
157:   for (i=0; i<40; i++) {
158:     VecSetRandom(rand,xx);
159:     VecSetRandom(rand,yy);
160:     MatMultTransposeAdd(A,xx,yy,s1);
161:     MatMultTransposeAdd(B,xx,yy,s2);
162:     VecNorm(s1,NORM_2,&s1norm);
163:     VecNorm(s2,NORM_2,&s2norm);
164:     rnorm = s2norm-s1norm;
165:     if (rnorm<-tol || rnorm>tol) {
166:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
167:     }
168:   }
169: 
170: 
171:   /* Do LUFactor() on both the matrices */
172:   PetscMalloc(M*sizeof(int),&idx);
173:   for (i=0; i<M; i++) idx[i] = i;
174:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is1);
175:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is2);
176:   PetscFree(idx);
177:   ISSetPermutation(is1);
178:   ISSetPermutation(is2);

180:   info.fill = 2.0;
181:   info.dtcol = 0.0;
182:   info.damping = 0.0;
183:   info.zeropivot = 1.e-14;
184:   info.pivotinblocks = 1.0;
185:   MatLUFactor(B,is1,is2,&info);
186:   MatLUFactor(A,is1,is2,&info);
187: 
188:   /* Test MatSolveAdd() */
189:   for (i=0; i<40; i++) {
190:     VecSetRandom(rand,xx);
191:     VecSetRandom(rand,yy);
192:     MatSolveAdd(B,xx,yy,s2);
193:     MatSolveAdd(A,xx,yy,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 - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
199:     }
200:   }
201: 
202:   /* Test MatSolveAdd() when x = A'b +x */
203:   for (i=0; i<40; i++) {
204:     VecSetRandom(rand,xx);
205:     VecSetRandom(rand,s1);
206:     VecCopy(s2,s1);
207:     MatSolveAdd(B,xx,s2,s2);
208:     MatSolveAdd(A,xx,s1,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(same) - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
214:     }
215:   }
216: 
217:   /* Test MatSolve() */
218:   for (i=0; i<40; i++) {
219:     VecSetRandom(rand,xx);
220:     MatSolve(B,xx,s2);
221:     MatSolve(A,xx,s1);
222:     VecNorm(s1,NORM_2,&s1norm);
223:     VecNorm(s2,NORM_2,&s2norm);
224:     rnorm = s2norm-s1norm;
225:     if (rnorm<-tol || rnorm>tol) {
226:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
227:     }
228:   }
229: 
230:   /* Test MatSolveTranspose() */
231:   if (bs < 8) {
232:     for (i=0; i<40; i++) {
233:       VecSetRandom(rand,xx);
234:       MatSolveTranspose(B,xx,s2);
235:       MatSolveTranspose(A,xx,s1);
236:       VecNorm(s1,NORM_2,&s1norm);
237:       VecNorm(s2,NORM_2,&s2norm);
238:       rnorm = s2norm-s1norm;
239:       if (rnorm<-tol || rnorm>tol) {
240:         PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
241:       }
242:     }
243:   }

245:   MatDestroy(A);
246:   MatDestroy(B);
247:   VecDestroy(xx);
248:   VecDestroy(s1);
249:   VecDestroy(s2);
250:   VecDestroy(yy);
251:   ISDestroy(is1);
252:   ISDestroy(is2);
253:   PetscRandomDestroy(rand);
254:   PetscFinalize();
255:   return 0;
256: }