Actual source code: ex77.c

  2: static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n";

 4:  #include petscmat.h

  8: int main(int argc,char **args)
  9: {
 10:   Vec            x,y,b,s1,s2;
 11:   Mat            A;           /* linear system matrix */
 12:   Mat            sA;         /* symmetric part of the matrices */
 13:   PetscInt       n,mbs=16,bs=1,nz=3,prob=2,i,j,col[3],row,Ii,J,n1;
 14:   const PetscInt *ip_ptr;
 15:   PetscScalar    neg_one = -1.0,value[3],alpha=0.1;
 16:   PetscMPIInt    size;
 18:   IS             ip, isrow, iscol;
 19:   PetscRandom    rdm;
 20:   PetscTruth     reorder=PETSC_FALSE;
 21:   MatInfo        minfo1,minfo2;
 22:   PetscReal      norm1,norm2,tol=1.e-10;

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

 30:   n = mbs*bs;
 31:   ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &A);
 32:   ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &sA);

 34:   /* Test MatGetOwnershipRange() */
 35:   MatGetOwnershipRange(A,&Ii,&J);
 36:   MatGetOwnershipRange(sA,&i,&j);
 37:   if (i-Ii || j-J){
 38:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
 39:   }

 41:   /* Assemble matrix */
 42:   if (bs == 1){
 43:     PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
 44:     if (prob == 1){ /* tridiagonal matrix */
 45:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
 46:       for (i=1; i<n-1; i++) {
 47:         col[0] = i-1; col[1] = i; col[2] = i+1;
 48:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 49:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 50:       }
 51:       i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
 52:       value[0]= 0.1; value[1]=-1; value[2]=2;
 53:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 54:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);

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

110:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
111:       value[0]=4.0; value[1] = -1.0;
112:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
113:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
114:     }
115: #endif
116:     /* off-diagonal blocks */
117:     value[0]=-1.0;
118:     for (i=0; i<(n/bs-1)*bs; i++){
119:       col[0]=i+bs;
120:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
121:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
122:       col[0]=i; row=i+bs;
123:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
124:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
125:     }
126:   }
127:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
128:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
129:   /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
130:   MatView(A, VIEWER_DRAW_WORLD);
131:   MatView(A, VIEWER_STDOUT_WORLD); */

133:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
134:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
135:   /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
136:   MatView(sA, VIEWER_DRAW_WORLD); 
137:   MatView(sA, VIEWER_STDOUT_WORLD); 
138:   */

140:   /* Test MatNorm() */
141:   MatNorm(A,NORM_FROBENIUS,&norm1);
142:   MatNorm(sA,NORM_FROBENIUS,&norm2);
143:   norm1 -= norm2;
144:   if (norm1<-tol || norm1>tol){
145:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
146:   }
147:   MatNorm(A,NORM_INFINITY,&norm1);
148:   MatNorm(sA,NORM_INFINITY,&norm2);
149:   norm1 -= norm2;
150:   if (norm1<-tol || norm1>tol){
151:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
152:   }
153: 
154:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
155:   MatGetInfo(A,MAT_LOCAL,&minfo1);
156:   MatGetInfo(sA,MAT_LOCAL,&minfo2);
157:   /*
158:   printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated); 
159:   printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated); 
160:   */
161:   i = (int) (minfo1.nz_used - minfo2.nz_used);
162:   j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
163:   if (i<0 || j<0) {
164:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
165:   }

167:   MatGetSize(A,&Ii,&J);
168:   MatGetSize(sA,&i,&j);
169:   if (i-Ii || j-J) {
170:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
171:   }
172: 
173:   MatGetBlockSize(A, &Ii);
174:   MatGetBlockSize(sA, &i);
175:   if (i-Ii){
176:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
177:   }

179:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
180:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
181:   PetscRandomSetFromOptions(rdm);
182:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
183:   VecDuplicate(x,&s1);
184:   VecDuplicate(x,&s2);
185:   VecDuplicate(x,&y);
186:   VecDuplicate(x,&b);

188:   VecSetRandom(x,rdm);

190:   MatDiagonalScale(A,x,x);
191:   MatDiagonalScale(sA,x,x);

193:   MatGetDiagonal(A,s1);
194:   MatGetDiagonal(sA,s2);
195:   VecNorm(s1,NORM_1,&norm1);
196:   VecNorm(s2,NORM_1,&norm2);
197:   norm1 -= norm2;
198:   if (norm1<-tol || norm1>tol) {
199:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
200:   }

202:   MatScale(A,alpha);
203:   MatScale(sA,alpha);

205:   /* Test MatMult(), MatMultAdd() */
206:   for (i=0; i<40; i++) {
207:     VecSetRandom(x,rdm);
208:     MatMult(A,x,s1);
209:     MatMult(sA,x,s2);
210:     VecNorm(s1,NORM_1,&norm1);
211:     VecNorm(s2,NORM_1,&norm2);
212:     norm1 -= norm2;
213:     if (norm1<-tol || norm1>tol) {
214:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n");
215:     }
216:   }

218:   for (i=0; i<40; i++) {
219:     VecSetRandom(x,rdm);
220:     VecSetRandom(y,rdm);
221:     MatMultAdd(A,x,y,s1);
222:     MatMultAdd(sA,x,y,s2);
223:     VecNorm(s1,NORM_1,&norm1);
224:     VecNorm(s2,NORM_1,&norm2);
225:     norm1 -= norm2;
226:     if (norm1<-tol || norm1>tol) {
227:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
228:     }
229:   }

231:   /* Test MatReordering() */
232:   MatGetOrdering(A,MATORDERING_NATURAL,&isrow,&iscol);
233:   ip = isrow;

235:   if (reorder){
236:     IS       nip;
237:     PetscInt *nip_ptr;
238:     PetscMalloc(mbs*sizeof(PetscInt),&nip_ptr);
239:     ISGetIndices(ip,&ip_ptr);
240:     PetscMemcpy(nip_ptr,ip_ptr,mbs*sizeof(PetscInt));
241:     i = nip_ptr[1]; nip_ptr[1] = nip_ptr[mbs-2]; nip_ptr[mbs-2] = i;
242:     i = nip_ptr[0]; nip_ptr[0] = nip_ptr[mbs-1]; nip_ptr[mbs-1] = i;
243:     ISRestoreIndices(ip,&ip_ptr);
244:     ISCreateGeneral(PETSC_COMM_SELF,mbs,nip_ptr,&nip);
245:     PetscFree(nip_ptr);

247:     MatReorderingSeqSBAIJ(sA, ip);
248:     ISDestroy(nip);
249:     /* ISView(ip, VIEWER_STDOUT_SELF); 
250:        MatView(sA,VIEWER_DRAW_SELF); */
251:   }
252: 
253:   ISDestroy(iscol);
254:   /* ISDestroy(isrow);*/

256:   ISDestroy(isrow);
257:   MatDestroy(A);
258:   MatDestroy(sA);
259:   VecDestroy(x);
260:   VecDestroy(y);
261:   VecDestroy(s1);
262:   VecDestroy(s2);
263:   VecDestroy(b);
264:   PetscRandomDestroy(rdm);

266:   PetscFinalize();
267:   return 0;
268: }