Actual source code: ex114.c


  2: static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n";

  4: #include <petscmat.h>

  6: #define M 5
  7: #define N 6

  9: int main(int argc,char **args)
 10: {
 11:   Mat            A;
 12:   Vec            min,max,maxabs,e;
 13:   PetscInt       m,n,j,imin[M],imax[M],imaxabs[M],indices[N],row,testcase=0;
 14:   PetscScalar    values[N];
 15:   PetscMPIInt    size,rank;
 16:   PetscReal      enorm;

 18:   PetscInitialize(&argc,&args,(char*)0,help);
 19:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 20:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 21:   PetscOptionsGetInt(NULL,NULL,"-testcase",&testcase,NULL);

 23:   MatCreate(PETSC_COMM_WORLD,&A);
 24:   if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */
 25:     if (rank == 0) {
 26:       MatSetSizes(A,M,N,PETSC_DECIDE,PETSC_DECIDE);
 27:     } else {
 28:       MatSetSizes(A,0,0,PETSC_DECIDE,PETSC_DECIDE);
 29:     }
 30:     testcase = 0;
 31:   } else {
 32:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
 33:   }
 34:   MatSetFromOptions(A);
 35:   MatSetUp(A);

 37:   if (rank == 0) { /* proc[0] sets matrix A */
 38:     for (j=0; j<N; j++) indices[j] = j;
 39:     switch (testcase) {
 40:     case 1: /* see testcast 0 */
 41:       break;
 42:     case 2:
 43:       row = 0;
 44:       values[0]  = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -4.0; values[4] = 1.0; values[5] = 1.0;
 45:       MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);
 46:       row = 2;
 47:       indices[0] = 0;    indices[1] = 3;    indices[2] = 5;
 48:       values[0]  = -2.0; values[1]  = -2.0; values[2]  = -2.0;
 49:       MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
 50:       row = 3;
 51:       indices[0] = 0;    indices[1] = 1;    indices[2] = 4;
 52:       values[0]  = -2.0; values[1]  = -2.0; values[2]  = -2.0;
 53:       MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
 54:       row = 4;
 55:       indices[0] = 0;    indices[1] = 1;    indices[2] = 2;
 56:       values[0]  = -2.0; values[1]  = -2.0; values[2]  = -2.0;
 57:       MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
 58:       break;
 59:     case 3:
 60:       row = 0;
 61:       values[0]  = -2.0; values[1] = -2.0; values[2] = -2.0;
 62:       MatSetValues(A,1,&row,3,indices+1,values,INSERT_VALUES);
 63:       row = 1;
 64:       values[0]  = -2.0; values[1] = -2.0; values[2] = -2.0;
 65:       MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
 66:       row = 2;
 67:       values[0]  = -2.0; values[1] = -2.0; values[2]  = -2.0;
 68:       MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
 69:       row = 3;
 70:       values[0]  = -2.0; values[1] = -2.0; values[2]  = -2.0; values[3] = -1.0;
 71:       MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);
 72:       row = 4;
 73:       values[0]  = -2.0; values[1] = -2.0; values[2]  = -2.0; values[3] = -1.0;
 74:       MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);
 75:       break;

 77:     default:
 78:       row  = 0;
 79:       values[0]  = -1.0; values[1] = 0.0; values[2] = 1.0; values[3] = 3.0; values[4] = 4.0; values[5] = -5.0;
 80:       MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);
 81:       row  = 1;
 82:       MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
 83:       row  = 3;
 84:       MatSetValues(A,1,&row,1,indices+4,values+4,INSERT_VALUES);
 85:       row  = 4;
 86:       MatSetValues(A,1,&row,2,indices+4,values+4,INSERT_VALUES);
 87:     }
 88:   }
 89:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 90:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 91:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

 93:   MatGetLocalSize(A, &m,&n);
 94:   VecCreate(PETSC_COMM_WORLD,&min);
 95:   VecSetSizes(min,m,PETSC_DECIDE);
 96:   VecSetFromOptions(min);
 97:   VecDuplicate(min,&max);
 98:   VecDuplicate(min,&maxabs);
 99:   VecDuplicate(min,&e);

101:   /* MatGetRowMax() */
102:   PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMax\n");
103:   MatGetRowMax(A,max,NULL);
104:   MatGetRowMax(A,max,imax);
105:   VecView(max,PETSC_VIEWER_STDOUT_WORLD);
106:   VecGetLocalSize(max,&n);
107:   PetscIntView(n,imax,PETSC_VIEWER_STDOUT_WORLD);

109:   /* MatGetRowMin() */
110:   MatScale(A,-1.0);
111:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);
112:   PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMin\n");
113:   MatGetRowMin(A,min,NULL);
114:   MatGetRowMin(A,min,imin);

116:   VecWAXPY(e,1.0,max,min); /* e = max + min */
117:   VecNorm(e,NORM_INFINITY,&enorm);
119:   for (j = 0; j < n; j++) {
121:   }

123:   /* MatGetRowMaxAbs() */
124:   PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs\n");
125:   MatGetRowMaxAbs(A,maxabs,NULL);
126:   MatGetRowMaxAbs(A,maxabs,imaxabs);
127:   VecView(maxabs,PETSC_VIEWER_STDOUT_WORLD);
128:   PetscIntView(n,imaxabs,PETSC_VIEWER_STDOUT_WORLD);

130:   /* MatGetRowMinAbs() */
131:   PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMinAbs\n");
132:   MatGetRowMinAbs(A,min,NULL);
133:   MatGetRowMinAbs(A,min,imin);
134:   VecView(min,PETSC_VIEWER_STDOUT_WORLD);
135:   PetscIntView(n,imin,PETSC_VIEWER_STDOUT_WORLD);

137:   if (size == 1) {
138:     /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */
139:     Mat Adense;
140:     Vec max_d,maxabs_d;
141:     VecDuplicate(min,&max_d);
142:     VecDuplicate(min,&maxabs_d);

144:     MatScale(A,-1.0);
145:     MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
146:     PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMax for seqdense matrix\n");
147:     MatGetRowMax(Adense,max_d,imax);

149:     VecWAXPY(e,-1.0,max,max_d); /* e = -max + max_d */
150:     VecNorm(e,NORM_INFINITY,&enorm);

153:     MatScale(Adense,-1.0);
154:     PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMin for seqdense matrix\n");
155:     MatGetRowMin(Adense,min,imin);

157:     VecWAXPY(e,1.0,max,min); /* e = max + min */
158:     VecNorm(e,NORM_INFINITY,&enorm);
160:     for (j = 0; j < n; j++) {
161:       if (imin[j] != imax[j]) {
162:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT " for seqdense matrix",j,imin[j],imax[j]);
163:       }
164:     }

166:     PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n");
167:     MatGetRowMaxAbs(Adense,maxabs_d,imaxabs);
168:     VecWAXPY(e,-1.0,maxabs,maxabs_d); /* e = -maxabs + maxabs_d */
169:     VecNorm(e,NORM_INFINITY,&enorm);

172:     MatDestroy(&Adense);
173:     VecDestroy(&max_d);
174:     VecDestroy(&maxabs_d);
175:   }

177:   { /* BAIJ matrix */
178:     Mat               B;
179:     Vec               maxabsB,maxabsB2;
180:     PetscInt          bs=2,*imaxabsB,*imaxabsB2,rstart,rend,cstart,cend,ncols,col,Brows[2],Bcols[2];
181:     const PetscInt    *cols;
182:     const PetscScalar *vals,*vals2;
183:     PetscScalar       Bvals[4];

185:     PetscMalloc2(M,&imaxabsB,bs*M,&imaxabsB2);

187:     /* bs = 1 */
188:     MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B);
189:     VecDuplicate(min,&maxabsB);
190:     MatGetRowMaxAbs(B,maxabsB,NULL);
191:     MatGetRowMaxAbs(B,maxabsB,imaxabsB);
192:     PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for BAIJ matrix\n");
193:     VecWAXPY(e,-1.0,maxabs,maxabsB); /* e = -maxabs + maxabsB */
194:     VecNorm(e,NORM_INFINITY,&enorm);

197:     for (j = 0; j < n; j++) {
199:     }
200:     MatDestroy(&B);

202:     /* Test bs = 2: Create B with bs*bs block structure of A */
203:     VecCreate(PETSC_COMM_WORLD,&maxabsB2);
204:     VecSetSizes(maxabsB2,bs*m,PETSC_DECIDE);
205:     VecSetFromOptions(maxabsB2);

207:     MatGetOwnershipRange(A,&rstart,&rend);
208:     MatGetOwnershipRangeColumn(A,&cstart,&cend);
209:     MatCreate(PETSC_COMM_WORLD,&B);
210:     MatSetSizes(B,bs*(rend-rstart),bs*(cend-cstart),PETSC_DECIDE,PETSC_DECIDE);
211:     MatSetFromOptions(B);
212:     MatSetUp(B);

214:     for (row=rstart; row<rend; row++) {
215:       MatGetRow(A,row,&ncols,&cols,&vals);
216:       for (col=0; col<ncols; col++) {
217:         for (j=0; j<bs; j++) {
218:           Brows[j] = bs*row + j;
219:           Bcols[j] = bs*cols[col]+j;
220:         }
221:         for (j=0; j<bs*bs; j++) Bvals[j] = vals[col];
222:         MatSetValues(B,bs,Brows,bs,Bcols,Bvals,INSERT_VALUES);
223:       }
224:       MatRestoreRow(A,row,&ncols,&cols,&vals);
225:     }
226:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
227:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

229:     MatGetRowMaxAbs(B,maxabsB2,imaxabsB2);

231:     /* Check maxabsB2 and imaxabsB2 */
232:     VecGetArrayRead(maxabsB,&vals);
233:     VecGetArrayRead(maxabsB2,&vals2);
234:     for (row=0; row<m; row++) {
235:       if (PetscAbsScalar(vals[row] - vals2[bs*row]) > PETSC_MACHINE_EPSILON)
236:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"row %" PetscInt_FMT " maxabsB != maxabsB2",row);
237:     }
238:     VecRestoreArrayRead(maxabsB,&vals);
239:     VecRestoreArrayRead(maxabsB2,&vals2);

241:     for (col=0; col<n; col++) {
242:       if (imaxabsB[col] != imaxabsB2[bs*col]/bs)
243:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"col %" PetscInt_FMT " imaxabsB != imaxabsB2",col);
244:     }
245:     VecDestroy(&maxabsB);
246:     MatDestroy(&B);
247:     VecDestroy(&maxabsB2);
248:     PetscFree2(imaxabsB,imaxabsB2);
249:   }

251:   VecDestroy(&min);
252:   VecDestroy(&max);
253:   VecDestroy(&maxabs);
254:   VecDestroy(&e);
255:   MatDestroy(&A);
256:   PetscFinalize();
257:   return 0;
258: }

260: /*TEST

262:    test:
263:       output_file: output/ex114.out

265:    test:
266:       suffix: 2
267:       args: -testcase 1
268:       output_file: output/ex114.out

270:    test:
271:       suffix: 3
272:       args: -testcase 2
273:       output_file: output/ex114_3.out

275:    test:
276:       suffix: 4
277:       args: -testcase 3
278:       output_file: output/ex114_4.out

280:    test:
281:       suffix: 5
282:       nsize: 3
283:       args: -testcase 0
284:       output_file: output/ex114_5.out

286:    test:
287:       suffix: 6
288:       nsize: 3
289:       args: -testcase 1
290:       output_file: output/ex114_6.out

292:    test:
293:       suffix: 7
294:       nsize: 3
295:       args: -testcase 2
296:       output_file: output/ex114_7.out

298:    test:
299:       suffix: 8
300:       nsize: 3
301:       args: -testcase 3
302:       output_file: output/ex114_8.out

304: TEST*/