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>

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

 21:   PetscInitialize(&argc,&args,(char*)0,help);
 22:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 24:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 25:   PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);

 27:   n   = mbs*bs;
 28:   MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &A);
 29:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 30:   MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &sA);
 31:   MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

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

 40:   /* Assemble matrix */
 41:   if (bs == 1) {
 42:     PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
 43:     if (prob == 1) { /* tridiagonal matrix */
 44:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
 45:       for (i=1; i<n-1; i++) {
 46:         col[0] = i-1; col[1] = i; col[2] = i+1;
 47:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 48:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 49:       }
 50:       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;

 58:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
 59:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 60:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 61:     } else if (prob ==2) { /* matrix for the five point stencil */
 62:       n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
 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:       }
 89:     }
 90:   } else { /* bs > 1 */
 91: #if defined(DIAGB)
 92:     for (block=0; block<n/bs; block++) {
 93:       /* diagonal blocks */
 94:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 95:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 96:         col[0] = i-1; col[1] = i; col[2] = i+1;
 97:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 98:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 99:       }
100:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;

102:       value[0]=-1.0; value[1]=4.0;
103:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
104:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

106:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;

108:       value[0]=4.0; value[1] = -1.0;
109:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
110:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
111:     }
112: #endif
113:     /* off-diagonal blocks */
114:     value[0]=-1.0;
115:     for (i=0; i<(n/bs-1)*bs; i++) {
116:       col[0]=i+bs;
117:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
118:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
119:       col[0]=i; row=i+bs;
120:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
121:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
122:     }
123:   }
124:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
125:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

127:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
128:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);

130:   /* Test MatNorm() */
131:   MatNorm(A,NORM_FROBENIUS,&norm1);
132:   MatNorm(sA,NORM_FROBENIUS,&norm2);
133:   norm1 -= norm2;
134:   if (norm1<-tol || norm1>tol) {
135:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",(double)norm1);
136:   }
137:   MatNorm(A,NORM_INFINITY,&norm1);
138:   MatNorm(sA,NORM_INFINITY,&norm2);
139:   norm1 -= norm2;
140:   if (norm1<-tol || norm1>tol) {
141:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",(double)norm1);
142:   }

144:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
145:   MatGetInfo(A,MAT_LOCAL,&minfo1);
146:   MatGetInfo(sA,MAT_LOCAL,&minfo2);
147:   i = (int) (minfo1.nz_used - minfo2.nz_used);
148:   j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
149:   if (i<0 || j<0) {
150:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
151:   }

153:   MatGetSize(A,&Ii,&J);
154:   MatGetSize(sA,&i,&j);
155:   if (i-Ii || j-J) {
156:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
157:   }

159:   MatGetBlockSize(A, &Ii);
160:   MatGetBlockSize(sA, &i);
161:   if (i-Ii) {
162:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
163:   }

165:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
166:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
167:   PetscRandomSetFromOptions(rdm);
168:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
169:   VecDuplicate(x,&s1);
170:   VecDuplicate(x,&s2);
171:   VecDuplicate(x,&y);
172:   VecDuplicate(x,&b);

174:   VecSetRandom(x,rdm);

176:   MatDiagonalScale(A,x,x);
177:   MatDiagonalScale(sA,x,x);

179:   MatGetDiagonal(A,s1);
180:   MatGetDiagonal(sA,s2);
181:   VecNorm(s1,NORM_1,&norm1);
182:   VecNorm(s2,NORM_1,&norm2);
183:   norm1 -= norm2;
184:   if (norm1<-tol || norm1>tol) {
185:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
186:   }

188:   MatScale(A,alpha);
189:   MatScale(sA,alpha);

191:   /* Test MatMult(), MatMultAdd() */
192:   for (i=0; i<40; i++) {
193:     VecSetRandom(x,rdm);
194:     MatMult(A,x,s1);
195:     MatMult(sA,x,s2);
196:     VecNorm(s1,NORM_1,&norm1);
197:     VecNorm(s2,NORM_1,&norm2);
198:     norm1 -= norm2;
199:     if (norm1<-tol || norm1>tol) {
200:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n");
201:     }
202:   }

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

217:   /* Test MatReordering() */
218:   MatGetOrdering(A,MATORDERINGNATURAL,&isrow,&iscol);
219:   ip   = isrow;

221:   if (reorder) {
222:     IS       nip;
223:     PetscInt *nip_ptr;
224:     PetscMalloc1(mbs,&nip_ptr);
225:     ISGetIndices(ip,&ip_ptr);
226:     PetscArraycpy(nip_ptr,ip_ptr,mbs);
227:     i    = nip_ptr[1]; nip_ptr[1] = nip_ptr[mbs-2]; nip_ptr[mbs-2] = i;
228:     i    = nip_ptr[0]; nip_ptr[0] = nip_ptr[mbs-1]; nip_ptr[mbs-1] = i;
229:     ISRestoreIndices(ip,&ip_ptr);
230:     ISCreateGeneral(PETSC_COMM_SELF,mbs,nip_ptr,PETSC_COPY_VALUES,&nip);
231:     PetscFree(nip_ptr);

233:     MatReorderingSeqSBAIJ(sA, ip);
234:     ISDestroy(&nip);
235:   }

237:   ISDestroy(&iscol);
238:   ISDestroy(&isrow);
239:   MatDestroy(&A);
240:   MatDestroy(&sA);
241:   VecDestroy(&x);
242:   VecDestroy(&y);
243:   VecDestroy(&s1);
244:   VecDestroy(&s2);
245:   VecDestroy(&b);
246:   PetscRandomDestroy(&rdm);

248:   PetscFinalize();
249:   return 0;
250: }

252: /*TEST

254:    test:
255:       args: -bs {{1 2 3 4 5 6 7 8}}

257: TEST*/