Actual source code: ex74.c

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

 4:  #include petscmat.h

  8: int main(int argc,char **args)
  9: {
 10:   PetscMPIInt        size;
 11:   PetscErrorCode     ierr;
 12:   Vec                x,y,b,s1,s2;
 13:   Mat                A;             /* linear system matrix */
 14:   Mat                sA,sB,sC;         /* symmetric part of the matrices */
 15:   PetscInt           n,mbs=16,bs=1,nz=3,prob=1,i,j,col[3],lf,block, row,I,J,n1,inc;
 16:   PetscReal          norm1,norm2,rnorm,tol=1.e-10;
 17:   PetscScalar        neg_one = -1.0,four=4.0,value[3],alpha=0.1;
 18:   IS                 perm, iscol;
 19:   PetscRandom        rdm;
 20:   PetscTruth         doIcc=PETSC_TRUE,equal;
 21:   MatInfo            minfo1,minfo2;
 22:   MatFactorInfo      factinfo;
 23:   MatType            type;

 25:   PetscInitialize(&argc,&args,(char *)0,help);
 26:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 27:   if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");
 28:   PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
 29:   PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);

 31:   n = mbs*bs;
 32:   ierr=MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,n,n,nz,PETSC_NULL, &A);
 33:   MatCreate(PETSC_COMM_SELF,&sA);
 34:   MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
 35:   MatSetType(sA,MATSEQSBAIJ);
 36:   /* -mat_type <seqsbaij_derived type>, e.g., seqsbaijspooles, sbaijmumps */
 37:   MatSetFromOptions(sA);
 38:   MatGetType(sA,&type);
 39:   PetscTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
 40:   MatSeqSBAIJSetPreallocation(sA,bs,nz,PETSC_NULL);
 41:   MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR);

 43:   /* Test MatGetOwnershipRange() */
 44:   MatGetOwnershipRange(A,&I,&J);
 45:   MatGetOwnershipRange(sA,&i,&j);
 46:   if (i-I || j-J){
 47:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
 48:   }

 50:   /* Assemble matrix */
 51:   if (bs == 1){
 52:     PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
 53:     if (prob == 1){ /* tridiagonal matrix */
 54:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
 55:       for (i=1; i<n-1; i++) {
 56:         col[0] = i-1; col[1] = i; col[2] = i+1;
 57:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 58:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 59:       }
 60:       i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
 61:       value[0]= 0.1; value[1]=-1; value[2]=2;
 62:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 63:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);

 65:       i = 0;
 66:       col[0] = n-1;  col[1] = 1; col[2]=0;
 67:       value[0] = 0.1; value[1] = -1.0; value[2]=2;
 68:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 69:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 70:     }
 71:     else if (prob ==2){ /* matrix for the five point stencil */
 72:       n1 = (int) (sqrt((PetscReal)n) + 0.001);
 73:       if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
 74:       for (i=0; i<n1; i++) {
 75:         for (j=0; j<n1; j++) {
 76:           I = j + n1*i;
 77:           if (i>0)   {
 78:             J = I - n1;
 79:             MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
 80:             MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
 81:           }
 82:           if (i<n1-1) {
 83:             J = I + n1;
 84:             MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
 85:             MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
 86:           }
 87:           if (j>0)   {
 88:             J = I - 1;
 89:             MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
 90:             MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
 91:           }
 92:           if (j<n1-1) {
 93:             J = I + 1;
 94:             MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
 95:             MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
 96:           }
 97:           MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
 98:           MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
 99:         }
100:       }
101:     }
102:   }
103:   else { /* bs > 1 */
104:     for (block=0; block<n/bs; block++){
105:       /* diagonal blocks */
106:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
107:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
108:         col[0] = i-1; col[1] = i; col[2] = i+1;
109:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
110:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
111:       }
112:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
113:       value[0]=-1.0; value[1]=4.0;
114:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
115:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

117:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
118:       value[0]=4.0; value[1] = -1.0;
119:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
120:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
121:     }
122:     /* off-diagonal blocks */
123:     value[0]=-1.0;
124:     for (i=0; i<(n/bs-1)*bs; i++){
125:       col[0]=i+bs;
126:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
127:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
128:       col[0]=i; row=i+bs;
129:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
130:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
131:     }
132:   }
133:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
134:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
135:   /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
136:   MatView(A, PETSC_VIEWER_DRAW_WORLD);
137:   MatView(A, PETSC_VIEWER_STDOUT_WORLD); */

139:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
140:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
141:   /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
142:   MatView(sA, PETSC_VIEWER_DRAW_WORLD); 
143:   MatView(sA, PETSC_VIEWER_STDOUT_WORLD); 
144:   */

146:   /* Test MatDuplicate() */
147:   MatNorm(A,NORM_FROBENIUS,&norm1);
148:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
149:   MatEqual(sA,sB,&equal);
150:   if (!equal) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()");

152:   /* Test MatNorm() */
153:   MatNorm(A,NORM_FROBENIUS,&norm1);
154:   MatNorm(sB,NORM_FROBENIUS,&norm2);
155:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
156:   if (rnorm > tol){
157:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
158:   }
159:   MatNorm(A,NORM_INFINITY,&norm1);
160:   MatNorm(sB,NORM_INFINITY,&norm2);
161:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
162:   if (rnorm > tol){
163:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
164:   }
165:   MatNorm(A,NORM_1,&norm1);
166:   MatNorm(sB,NORM_1,&norm2);
167:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
168:   if (rnorm > tol){
169:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
170:   }

172:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
173:   MatGetInfo(A,MAT_LOCAL,&minfo1);
174:   MatGetInfo(sB,MAT_LOCAL,&minfo2);
175:   /*
176:   printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated); 
177:   printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated); 
178:   */
179:   i = (int) (minfo1.nz_used - minfo2.nz_used);
180:   j = (int) (minfo2.nz_allocated - minfo2.nz_used);
181:   if (i<0 || j<0) {
182:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
183:   }

185:   MatGetSize(A,&I,&J);
186:   MatGetSize(sB,&i,&j);
187:   if (i-I || j-J) {
188:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
189:   }
190: 
191:   MatGetBlockSize(A, &I);
192:   MatGetBlockSize(sB, &i);
193:   if (i-I){
194:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
195:   }

197:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
198:   PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rdm);
199:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
200:   VecDuplicate(x,&s1);
201:   VecDuplicate(x,&s2);
202:   VecDuplicate(x,&y);
203:   VecDuplicate(x,&b);
204:   VecSetRandom(x,rdm);

206:   MatDiagonalScale(A,x,x);
207:   MatDiagonalScale(sB,x,x);
208:   MatMultEqual(A,sB,10,&equal);
209:   if (!equal) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");

211:   MatGetDiagonal(A,s1);
212:   MatGetDiagonal(sB,s2);
213: 
214:   VecAXPY(s2,neg_one,s1);
215:   VecNorm(s2,NORM_1,&norm1);
216:   if ( norm1>tol) {
217:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%G\n",norm1);
218:   }
219:   /*
220:   VecNorm(s1,NORM_1,&norm1);
221:   VecNorm(s2,NORM_1,&norm2);
222:   norm1 -= norm2;
223:   if (norm1<-tol || norm1>tol) { 
224:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
225:   } 
226:   */

228:   MatScale(A,alpha);
229:   MatScale(sB,alpha);

231:   /* Test MatGetRowMax() */
232:   MatGetRowMax(A,s1);
233:   MatGetRowMax(sB,s2);
234:   VecNorm(s1,NORM_1,&norm1);
235:   VecNorm(s2,NORM_1,&norm2);
236:   norm1 -= norm2;
237:   if (norm1<-tol || norm1>tol) {
238:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMax() \n");
239:   }

241:   /* Test MatMult() */
242:   for (i=0; i<40; i++) {
243:     VecSetRandom(x,rdm);
244:     MatMult(A,x,s1);
245:     MatMult(sB,x,s2);
246:     VecNorm(s1,NORM_1,&norm1);
247:     VecNorm(s2,NORM_1,&norm2);
248:     norm1 -= norm2;
249:     if (norm1<-tol || norm1>tol) {
250:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %G\n",norm1);
251:     }
252:   }

254:   /* MatMultAdd() */
255:   for (i=0; i<40; i++) {
256:     VecSetRandom(x,rdm);
257:     VecSetRandom(y,rdm);
258:     MatMultAdd(A,x,y,s1);
259:     MatMultAdd(sB,x,y,s2);
260:     VecNorm(s1,NORM_1,&norm1);
261:     VecNorm(s2,NORM_1,&norm2);
262:     norm1 -= norm2;
263:     if (norm1<-tol || norm1>tol) {
264:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(),  norm1-norm2: %G\n",norm1);
265:     }
266:   }

268:   /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
269:   MatGetOrdering(A,MATORDERING_NATURAL,&perm,&iscol);
270:   ISDestroy(iscol);
271:   norm1 = tol;
272:   inc   = bs;

274:   /* initialize factinfo */
275:   PetscMemzero(&factinfo,sizeof(MatFactorInfo));

277:   for (lf=-1; lf<10; lf += inc){
278:     if (lf==-1) {  /* Cholesky factor of sB (duplicate sA) */
279:       factinfo.fill = 5.0;
280:       MatCholeskyFactorSymbolic(sB,perm,&factinfo,&sC);
281:     } else if (!doIcc){
282:       break;
283:     } else {       /* incomplete Cholesky factor */
284:       factinfo.fill   = 5.0;
285:       factinfo.levels = lf;
286:       MatICCFactorSymbolic(sB,perm,&factinfo,&sC);
287:     }
288:     MatCholeskyFactorNumeric(sB,&factinfo,&sC);
289:     /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */

291:     /* test MatGetDiagonal on numeric factor */
292:     /*
293:     if (lf == -1) {
294:       MatGetDiagonal(sC,s1);  
295:       printf(" in ex74.c, diag: \n");
296:       VecView(s1,PETSC_VIEWER_STDOUT_SELF);
297:     }
298:     */

300:     MatMult(sB,x,b);
301:     MatSolve(sC,b,y);
302:     MatDestroy(sC);
303: 
304:     /* Check the error */
305:     VecAXPY(y,neg_one,x);
306:     VecNorm(y,NORM_2,&norm2);
307:     /* printf("lf: %d, error: %G\n", lf,norm2); */
308:     if (10*norm1 < norm2 && lf-inc != -1){
309:       PetscPrintf(PETSC_COMM_SELF,"lf=%D, %D, Norm of error=%G, %G\n",lf-inc,lf,norm1,norm2);
310:     }
311:     norm1 = norm2;
312:     if (norm2 < tol && lf != -1) break;
313:   }

315:   ISDestroy(perm);

317:   MatDestroy(A);
318:   MatDestroy(sB);
319:   MatDestroy(sA);
320:   VecDestroy(x);
321:   VecDestroy(y);
322:   VecDestroy(s1);
323:   VecDestroy(s2);
324:   VecDestroy(b);
325:   PetscRandomDestroy(rdm);

327:   PetscFinalize();
328:   return 0;
329: }