Actual source code: ex114.c

petsc-3.14.0 2020-09-29
Report Typos and Errors

  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];
 16:   PetscMPIInt    size,rank;
 17:   PetscReal      enorm;

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

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

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

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

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

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

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

117:   VecWAXPY(e,1.0,max,min); /* e = max + min */
118:   VecNorm(e,NORM_INFINITY,&enorm);
119:   if (enorm > PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON ");
120:   for (j = 0; j < n; j++) {
121:     if (imin[j] != imax[j]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%D] %D != imax %D",j,imin[j],imax[j]);
122:   }

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

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

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

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

150:     VecWAXPY(e,-1.0,max,max_d); /* e = -max + max_d */
151:     VecNorm(e,NORM_INFINITY,&enorm);
152:     if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-max + max_d) %g > PETSC_MACHINE_EPSILON",enorm);

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

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

167:     PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n");
168:     MatGetRowMaxAbs(Adense,maxabs_d,imaxabs);
169:     VecWAXPY(e,-1.0,maxabs,maxabs_d); /* e = -maxabs + maxabs_d */
170:     VecNorm(e,NORM_INFINITY,&enorm);
171:     if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",enorm);

173:     MatDestroy(&Adense);
174:     VecDestroy(&max_d);
175:     VecDestroy(&maxabs_d);
176:   } else {
177:     /* BAIJ */
178:     Mat      B;
179:     Vec      maxabsB;
180:     PetscInt imaxabsB[M];
181:     MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B);
182:     VecDuplicate(min,&maxabsB);
183:     MatGetRowMaxAbs(A,maxabsB,NULL);
184:     MatGetRowMaxAbs(A,maxabsB,imaxabsB);
185:     PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for MPIBAIJ matrix\n");
186:     VecWAXPY(e,-1.0,maxabs,maxabsB); /* e = -maxabs + maxabsB */
187:     VecNorm(e,NORM_INFINITY,&enorm);
188:     if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",enorm);

190:     for (j = 0; j < n; j++) {
191:       if (imaxabs[j] != imaxabsB[j]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imaxabs[%D] %D != imaxabsB %D",j,imin[j],imax[j]);
192:     }
193:     MatDestroy(&B);
194:     VecDestroy(&maxabsB);
195:   }

197:   VecDestroy(&min);
198:   VecDestroy(&max);
199:   VecDestroy(&maxabs);
200:   VecDestroy(&e);
201:   MatDestroy(&A);
202:   PetscFinalize();
203:   return ierr;
204: }

206: /*TEST

208:    test:
209:       output_file: output/ex114.out

211:    test:
212:       suffix: 2
213:       args: -testcase 1
214:       output_file: output/ex114.out

216:    test:
217:       suffix: 3
218:       args: -testcase 2
219:       output_file: output/ex114_3.out

221:    test:
222:       suffix: 4
223:       args: -testcase 3
224:       output_file: output/ex114_4.out

226:    test:
227:       suffix: 5
228:       nsize: 3
229:       args: -testcase 0
230:       output_file: output/ex114_5.out

232:    test:
233:       suffix: 6
234:       nsize: 3
235:       args: -testcase 1
236:       output_file: output/ex114_6.out

238:    test:
239:       suffix: 7
240:       nsize: 3
241:       args: -testcase 2
242:       output_file: output/ex114_7.out

244:    test:
245:       suffix: 8
246:       nsize: 3
247:       args: -testcase 3
248:       output_file: output/ex114_8.out

250: TEST*/