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,neg_one=-1.0;
 15:   PetscRandom rand;
 16:   IS          is1,is2;
 17:   PetscReal   s1norm,s2norm,rnorm,tol = 1.e-8;
 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:     VecAXPY(&neg_one,s1,s2);
123:     VecNorm(s2,NORM_1,&s2norm);
124:     if (s2norm >tol) {
125:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), ||s1-s2||=%g\n",s2norm);
126:     }
127:     /*
128:     VecNorm(s1,NORM_2,&s1norm);
129:     VecNorm(s2,NORM_2,&s2norm);
130:     rnorm = s2norm-s1norm;
131:     if (rnorm<-tol || rnorm>tol) { 
132:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d \n",s1norm,s2norm,bs);
133:     */
134:   }
135: 
136:   /* Test MatMultAdd() */
137:   for (i=0; i<40; i++) {
138:     VecSetRandom(rand,xx);
139:     VecSetRandom(rand,yy);
140:     MatMultAdd(A,xx,yy,s1);
141:     MatMultAdd(B,xx,yy,s2);
142:     VecAXPY(&neg_one,s1,s2);
143:     VecNorm(s2,NORM_1,&s2norm);
144:     if (s2norm >tol) {
145:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd(), ||s1-s2||=%g\n",s2norm);
146:     }
147:     /*
148:     VecNorm(s1,NORM_2,&s1norm);
149:     VecNorm(s2,NORM_2,&s2norm);
150:     rnorm = s2norm-s1norm;
151:     if (rnorm<-tol || rnorm>tol) { 
152:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
153:     } 
154:     */
155:   }
156: 
157:   /* Test MatMultTranspose() */
158:   for (i=0; i<40; i++) {
159:     VecSetRandom(rand,xx);
160:     MatMultTranspose(A,xx,s1);
161:     MatMultTranspose(B,xx,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: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
167:     }
168:   }
169:   /* Test MatMultTransposeAdd() */
170:   for (i=0; i<40; i++) {
171:     VecSetRandom(rand,xx);
172:     VecSetRandom(rand,yy);
173:     MatMultTransposeAdd(A,xx,yy,s1);
174:     MatMultTransposeAdd(B,xx,yy,s2);
175:     VecNorm(s1,NORM_2,&s1norm);
176:     VecNorm(s2,NORM_2,&s2norm);
177:     rnorm = s2norm-s1norm;
178:     if (rnorm<-tol || rnorm>tol) {
179:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
180:     }
181:   }
182: 
183: 
184:   /* Do LUFactor() on both the matrices */
185:   PetscMalloc(M*sizeof(int),&idx);
186:   for (i=0; i<M; i++) idx[i] = i;
187:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is1);
188:   ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is2);
189:   PetscFree(idx);
190:   ISSetPermutation(is1);
191:   ISSetPermutation(is2);

193:   info.fill = 2.0;
194:   info.dtcol = 0.0;
195:   info.damping = 0.0;
196:   info.zeropivot = 1.e-14;
197:   info.pivotinblocks = 1.0;
198:   MatLUFactor(B,is1,is2,&info);
199:   MatLUFactor(A,is1,is2,&info);
200: 
201:   /* Test MatSolveAdd() */
202:   for (i=0; i<40; i++) {
203:     VecSetRandom(rand,xx);
204:     VecSetRandom(rand,yy);
205:     MatSolveAdd(B,xx,yy,s2);
206:     MatSolveAdd(A,xx,yy,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:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
212:     }
213:   }
214: 
215:   /* Test MatSolveAdd() when x = A'b +x */
216:   for (i=0; i<40; i++) {
217:     VecSetRandom(rand,xx);
218:     VecSetRandom(rand,s1);
219:     VecCopy(s2,s1);
220:     MatSolveAdd(B,xx,s2,s2);
221:     MatSolveAdd(A,xx,s1,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:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
227:     }
228:   }
229: 
230:   /* Test MatSolve() */
231:   for (i=0; i<40; i++) {
232:     VecSetRandom(rand,xx);
233:     MatSolve(B,xx,s2);
234:     MatSolve(A,xx,s1);
235:     VecNorm(s1,NORM_2,&s1norm);
236:     VecNorm(s2,NORM_2,&s2norm);
237:     rnorm = s2norm-s1norm;
238:     if (rnorm<-tol || rnorm>tol) {
239:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
240:     }
241:   }
242: 
243:   /* Test MatSolveTranspose() */
244:   if (bs < 8) {
245:     for (i=0; i<40; i++) {
246:       VecSetRandom(rand,xx);
247:       MatSolveTranspose(B,xx,s2);
248:       MatSolveTranspose(A,xx,s1);
249:       VecNorm(s1,NORM_2,&s1norm);
250:       VecNorm(s2,NORM_2,&s2norm);
251:       rnorm = s2norm-s1norm;
252:       if (rnorm<-tol || rnorm>tol) {
253:         PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
254:       }
255:     }
256:   }

258:   MatDestroy(A);
259:   MatDestroy(B);
260:   VecDestroy(xx);
261:   VecDestroy(s1);
262:   VecDestroy(s2);
263:   VecDestroy(yy);
264:   ISDestroy(is1);
265:   ISDestroy(is2);
266:   PetscRandomDestroy(rand);
267:   PetscFinalize();
268:   return 0;
269: }