Actual source code: ex75.c


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

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Vec            x,y,u,s1,s2;
  9:   Mat            A,sA,sB;
 10:   PetscRandom    rctx;
 11:   PetscReal      r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
 12:   PetscScalar    one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
 13:   PetscInt       n,col[3],n1,block,row,i,j,i2,j2,Ii,J,rstart,rend,bs=1,mbs=16,d_nz=3,o_nz=3,prob=1;
 14:   PetscMPIInt    size,rank;
 15:   PetscBool      flg;

 17:   PetscInitialize(&argc,&args,(char*)0,help);
 18:   PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
 19:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 20:   PetscOptionsGetInt(NULL,NULL,"-prob",&prob,NULL);

 22:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 25:   /* Create a BAIJ matrix A */
 26:   n = mbs*bs;
 27:   MatCreate(PETSC_COMM_WORLD,&A);
 28:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 29:   MatSetType(A,MATBAIJ);
 30:   MatSetFromOptions(A);
 31:   MatMPIBAIJSetPreallocation(A,bs,d_nz,NULL,o_nz,NULL);
 32:   MatSeqBAIJSetPreallocation(A,bs,d_nz,NULL);
 33:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 35:   if (bs == 1) {
 36:     if (prob == 1) { /* tridiagonal matrix */
 37:       value[0] = -1.0; value[1] = 0.0; value[2] = -1.0;
 38:       for (i=1; i<n-1; i++) {
 39:         col[0] = i-1; col[1] = i; col[2] = i+1;
 40:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 41:       }
 42:       i       = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
 43:       value[0]= 0.1; value[1]=-1.0; value[2]=0.0;
 44:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);

 46:       i        = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
 47:       value[0] = 0.0; value[1] = -1.0; value[2]=0.1;
 48:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 49:     } else if (prob ==2) { /* matrix for the five point stencil */
 50:       n1 = (int) PetscSqrtReal((PetscReal)n);
 51:       for (i=0; i<n1; i++) {
 52:         for (j=0; j<n1; j++) {
 53:           Ii = j + n1*i;
 54:           if (i>0)    {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 55:           if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 56:           if (j>0)    {J = Ii - 1;  MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 57:           if (j<n1-1) {J = Ii + 1;  MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 58:           MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
 59:         }
 60:       }
 61:     }
 62:     /* end of if (bs == 1) */
 63:   } else {  /* bs > 1 */
 64:     for (block=0; block<n/bs; block++) {
 65:       /* diagonal blocks */
 66:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 67:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 68:         col[0] = i-1; col[1] = i; col[2] = i+1;
 69:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 70:       }
 71:       i       = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
 72:       value[0]=-1.0; value[1]=4.0;
 73:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);

 75:       i       = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
 76:       value[0]=4.0; value[1] = -1.0;
 77:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
 78:     }
 79:     /* off-diagonal blocks */
 80:     value[0]=-1.0;
 81:     for (i=0; i<(n/bs-1)*bs; i++) {
 82:       col[0]=i+bs;
 83:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
 84:       col[0]=i; row=i+bs;
 85:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
 86:     }
 87:   }
 88:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 89:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 90:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

 92:   /* Get SBAIJ matrix sA from A */
 93:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);

 95:   /* Test MatGetSize(), MatGetLocalSize() */
 96:   MatGetSize(sA, &i,&j);
 97:   MatGetSize(A, &i2,&j2);
 98:   i   -= i2; j -= j2;
 99:   if (i || j) {
100:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
101:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
102:   }

104:   MatGetLocalSize(sA, &i,&j);
105:   MatGetLocalSize(A, &i2,&j2);
106:   i2  -= i; j2 -= j;
107:   if (i2 || j2) {
108:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
109:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
110:   }

112:   /* vectors */
113:   /*--------------------*/
114:   /* i is obtained from MatGetLocalSize() */
115:   VecCreate(PETSC_COMM_WORLD,&x);
116:   VecSetSizes(x,i,PETSC_DECIDE);
117:   VecSetFromOptions(x);
118:   VecDuplicate(x,&y);
119:   VecDuplicate(x,&u);
120:   VecDuplicate(x,&s1);
121:   VecDuplicate(x,&s2);

123:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
124:   PetscRandomSetFromOptions(rctx);
125:   VecSetRandom(x,rctx);
126:   VecSet(u,one);

128:   /* Test MatNorm() */
129:   MatNorm(A,NORM_FROBENIUS,&r1);
130:   MatNorm(sA,NORM_FROBENIUS,&r2);
131:   rnorm = PetscAbsReal(r1-r2)/r2;
132:   if (rnorm > tol && rank == 0) {
133:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs);
134:   }
135:   MatNorm(A,NORM_INFINITY,&r1);
136:   MatNorm(sA,NORM_INFINITY,&r2);
137:   rnorm = PetscAbsReal(r1-r2)/r2;
138:   if (rnorm > tol && rank == 0) {
139:     PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs);
140:   }
141:   MatNorm(A,NORM_1,&r1);
142:   MatNorm(sA,NORM_1,&r2);
143:   rnorm = PetscAbsReal(r1-r2)/r2;
144:   if (rnorm > tol && rank == 0) {
145:     PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs);
146:   }

148:   /* Test MatGetOwnershipRange() */
149:   MatGetOwnershipRange(sA,&rstart,&rend);
150:   MatGetOwnershipRange(A,&i2,&j2);
151:   i2  -= rstart; j2 -= rend;
152:   if (i2 || j2) {
153:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
154:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
155:   }

157:   /* Test MatDiagonalScale() */
158:   MatDiagonalScale(A,x,x);
159:   MatDiagonalScale(sA,x,x);
160:   MatMultEqual(A,sA,10,&flg);

163:   /* Test MatGetDiagonal(), MatScale() */
164:   MatGetDiagonal(A,s1);
165:   MatGetDiagonal(sA,s2);
166:   VecNorm(s1,NORM_1,&r1);
167:   VecNorm(s2,NORM_1,&r2);
168:   r1  -= r2;
169:   if (r1<-tol || r1>tol) {
170:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);
171:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
172:   }

174:   MatScale(A,alpha);
175:   MatScale(sA,alpha);

177:   /* Test MatGetRowMaxAbs() */
178:   MatGetRowMaxAbs(A,s1,NULL);
179:   MatGetRowMaxAbs(sA,s2,NULL);

181:   VecNorm(s1,NORM_1,&r1);
182:   VecNorm(s2,NORM_1,&r2);
183:   r1  -= r2;
184:   if (r1<-tol || r1>tol) {
185:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
186:   }

188:   /* Test MatMult(), MatMultAdd() */
189:   MatMultEqual(A,sA,10,&flg);
190:   if (!flg) {
191:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
192:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
193:   }

195:   MatMultAddEqual(A,sA,10,&flg);
196:   if (!flg) {
197:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
198:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
199:   }

201:   /* Test MatMultTranspose(), MatMultTransposeAdd() */
202:   for (i=0; i<10; i++) {
203:     VecSetRandom(x,rctx);
204:     MatMultTranspose(A,x,s1);
205:     MatMultTranspose(sA,x,s2);
206:     VecNorm(s1,NORM_1,&r1);
207:     VecNorm(s2,NORM_1,&r2);
208:     r1  -= r2;
209:     if (r1<-tol || r1>tol) {
210:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);
211:       PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
212:     }
213:   }
214:   for (i=0; i<10; i++) {
215:     VecSetRandom(x,rctx);
216:     VecSetRandom(y,rctx);
217:     MatMultTransposeAdd(A,x,y,s1);
218:     MatMultTransposeAdd(sA,x,y,s2);
219:     VecNorm(s1,NORM_1,&r1);
220:     VecNorm(s2,NORM_1,&r2);
221:     r1  -= r2;
222:     if (r1<-tol || r1>tol) {
223:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);
224:       PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
225:     }
226:   }

228:   /* Test MatDuplicate() */
229:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
230:   MatEqual(sA,sB,&flg);
231:   if (!flg) {
232:     PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
233:   }
234:   MatMultEqual(sA,sB,5,&flg);
235:   if (!flg) {
236:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
237:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
238:   }
239:   MatMultAddEqual(sA,sB,5,&flg);
240:   if (!flg) {
241:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
242:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
243:   }
244:   MatDestroy(&sB);
245:   VecDestroy(&u);
246:   VecDestroy(&x);
247:   VecDestroy(&y);
248:   VecDestroy(&s1);
249:   VecDestroy(&s2);
250:   MatDestroy(&sA);
251:   MatDestroy(&A);
252:   PetscRandomDestroy(&rctx);
253:   PetscFinalize();
254:   return 0;
255: }

257: /*TEST

259:    test:
260:       nsize: {{1 3}}
261:       args: -bs {{1 2 3  5  7 8}} -mat_ignore_lower_triangular -prob {{1 2}}

263: TEST*/