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;

 19:   PetscInitialize(&argc, &args, (char *)0, help);
 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 == 0) {
 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 == 0) { /* 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;
 46:       values[1] = -2.0;
 47:       values[2] = -2.0;
 48:       values[3] = -4.0;
 49:       values[4] = 1.0;
 50:       values[5] = 1.0;
 51:       MatSetValues(A, 1, &row, N, indices, values, INSERT_VALUES);
 52:       row        = 2;
 53:       indices[0] = 0;
 54:       indices[1] = 3;
 55:       indices[2] = 5;
 56:       values[0]  = -2.0;
 57:       values[1]  = -2.0;
 58:       values[2]  = -2.0;
 59:       MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES);
 60:       row        = 3;
 61:       indices[0] = 0;
 62:       indices[1] = 1;
 63:       indices[2] = 4;
 64:       values[0]  = -2.0;
 65:       values[1]  = -2.0;
 66:       values[2]  = -2.0;
 67:       MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES);
 68:       row        = 4;
 69:       indices[0] = 0;
 70:       indices[1] = 1;
 71:       indices[2] = 2;
 72:       values[0]  = -2.0;
 73:       values[1]  = -2.0;
 74:       values[2]  = -2.0;
 75:       MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES);
 76:       break;
 77:     case 3:
 78:       row       = 0;
 79:       values[0] = -2.0;
 80:       values[1] = -2.0;
 81:       values[2] = -2.0;
 82:       MatSetValues(A, 1, &row, 3, indices + 1, values, INSERT_VALUES);
 83:       row       = 1;
 84:       values[0] = -2.0;
 85:       values[1] = -2.0;
 86:       values[2] = -2.0;
 87:       MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES);
 88:       row       = 2;
 89:       values[0] = -2.0;
 90:       values[1] = -2.0;
 91:       values[2] = -2.0;
 92:       MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES);
 93:       row       = 3;
 94:       values[0] = -2.0;
 95:       values[1] = -2.0;
 96:       values[2] = -2.0;
 97:       values[3] = -1.0;
 98:       MatSetValues(A, 1, &row, 4, indices, values, INSERT_VALUES);
 99:       row       = 4;
100:       values[0] = -2.0;
101:       values[1] = -2.0;
102:       values[2] = -2.0;
103:       values[3] = -1.0;
104:       MatSetValues(A, 1, &row, 4, indices, values, INSERT_VALUES);
105:       break;

107:     default:
108:       row       = 0;
109:       values[0] = -1.0;
110:       values[1] = 0.0;
111:       values[2] = 1.0;
112:       values[3] = 3.0;
113:       values[4] = 4.0;
114:       values[5] = -5.0;
115:       MatSetValues(A, 1, &row, N, indices, values, INSERT_VALUES);
116:       row = 1;
117:       MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES);
118:       row = 3;
119:       MatSetValues(A, 1, &row, 1, indices + 4, values + 4, INSERT_VALUES);
120:       row = 4;
121:       MatSetValues(A, 1, &row, 2, indices + 4, values + 4, INSERT_VALUES);
122:     }
123:   }
124:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
125:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
126:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);

128:   MatGetLocalSize(A, &m, &n);
129:   VecCreate(PETSC_COMM_WORLD, &min);
130:   VecSetSizes(min, m, PETSC_DECIDE);
131:   VecSetFromOptions(min);
132:   VecDuplicate(min, &max);
133:   VecDuplicate(min, &maxabs);
134:   VecDuplicate(min, &e);

136:   /* MatGetRowMax() */
137:   PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMax\n");
138:   MatGetRowMax(A, max, NULL);
139:   MatGetRowMax(A, max, imax);
140:   VecView(max, PETSC_VIEWER_STDOUT_WORLD);
141:   VecGetLocalSize(max, &n);
142:   PetscIntView(n, imax, PETSC_VIEWER_STDOUT_WORLD);

144:   /* MatGetRowMin() */
145:   MatScale(A, -1.0);
146:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
147:   PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMin\n");
148:   MatGetRowMin(A, min, NULL);
149:   MatGetRowMin(A, min, imin);

151:   VecWAXPY(e, 1.0, max, min); /* e = max + min */
152:   VecNorm(e, NORM_INFINITY, &enorm);

156:   /* MatGetRowMaxAbs() */
157:   PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMaxAbs\n");
158:   MatGetRowMaxAbs(A, maxabs, NULL);
159:   MatGetRowMaxAbs(A, maxabs, imaxabs);
160:   VecView(maxabs, PETSC_VIEWER_STDOUT_WORLD);
161:   PetscIntView(n, imaxabs, PETSC_VIEWER_STDOUT_WORLD);

163:   /* MatGetRowMinAbs() */
164:   PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMinAbs\n");
165:   MatGetRowMinAbs(A, min, NULL);
166:   MatGetRowMinAbs(A, min, imin);
167:   VecView(min, PETSC_VIEWER_STDOUT_WORLD);
168:   PetscIntView(n, imin, PETSC_VIEWER_STDOUT_WORLD);

170:   if (size == 1) {
171:     /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */
172:     Mat Adense;
173:     Vec max_d, maxabs_d;
174:     VecDuplicate(min, &max_d);
175:     VecDuplicate(min, &maxabs_d);

177:     MatScale(A, -1.0);
178:     MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense);
179:     PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMax for seqdense matrix\n");
180:     MatGetRowMax(Adense, max_d, imax);

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

186:     MatScale(Adense, -1.0);
187:     PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMin for seqdense matrix\n");
188:     MatGetRowMin(Adense, min, imin);

190:     VecWAXPY(e, 1.0, max, min); /* e = max + min */
191:     VecNorm(e, NORM_INFINITY, &enorm);

195:     PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMaxAbs for seqdense matrix\n");
196:     MatGetRowMaxAbs(Adense, maxabs_d, imaxabs);
197:     VecWAXPY(e, -1.0, maxabs, maxabs_d); /* e = -maxabs + maxabs_d */
198:     VecNorm(e, NORM_INFINITY, &enorm);

201:     MatDestroy(&Adense);
202:     VecDestroy(&max_d);
203:     VecDestroy(&maxabs_d);
204:   }

206:   { /* BAIJ matrix */
207:     Mat                B;
208:     Vec                maxabsB, maxabsB2;
209:     PetscInt           bs = 2, *imaxabsB, *imaxabsB2, rstart, rend, cstart, cend, ncols, col, Brows[2], Bcols[2];
210:     const PetscInt    *cols;
211:     const PetscScalar *vals, *vals2;
212:     PetscScalar        Bvals[4];

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

216:     /* bs = 1 */
217:     MatConvert(A, MATMPIBAIJ, MAT_INITIAL_MATRIX, &B);
218:     VecDuplicate(min, &maxabsB);
219:     MatGetRowMaxAbs(B, maxabsB, NULL);
220:     MatGetRowMaxAbs(B, maxabsB, imaxabsB);
221:     PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMaxAbs for BAIJ matrix\n");
222:     VecWAXPY(e, -1.0, maxabs, maxabsB); /* e = -maxabs + maxabsB */
223:     VecNorm(e, NORM_INFINITY, &enorm);

227:     MatDestroy(&B);

229:     /* Test bs = 2: Create B with bs*bs block structure of A */
230:     VecCreate(PETSC_COMM_WORLD, &maxabsB2);
231:     VecSetSizes(maxabsB2, bs * m, PETSC_DECIDE);
232:     VecSetFromOptions(maxabsB2);

234:     MatGetOwnershipRange(A, &rstart, &rend);
235:     MatGetOwnershipRangeColumn(A, &cstart, &cend);
236:     MatCreate(PETSC_COMM_WORLD, &B);
237:     MatSetSizes(B, bs * (rend - rstart), bs * (cend - cstart), PETSC_DECIDE, PETSC_DECIDE);
238:     MatSetFromOptions(B);
239:     MatSetUp(B);

241:     for (row = rstart; row < rend; row++) {
242:       MatGetRow(A, row, &ncols, &cols, &vals);
243:       for (col = 0; col < ncols; col++) {
244:         for (j = 0; j < bs; j++) {
245:           Brows[j] = bs * row + j;
246:           Bcols[j] = bs * cols[col] + j;
247:         }
248:         for (j = 0; j < bs * bs; j++) Bvals[j] = vals[col];
249:         MatSetValues(B, bs, Brows, bs, Bcols, Bvals, INSERT_VALUES);
250:       }
251:       MatRestoreRow(A, row, &ncols, &cols, &vals);
252:     }
253:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
254:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

256:     MatGetRowMaxAbs(B, maxabsB2, imaxabsB2);

258:     /* Check maxabsB2 and imaxabsB2 */
259:     VecGetArrayRead(maxabsB, &vals);
260:     VecGetArrayRead(maxabsB2, &vals2);
262:     VecRestoreArrayRead(maxabsB, &vals);
263:     VecRestoreArrayRead(maxabsB2, &vals2);

266:     VecDestroy(&maxabsB);
267:     MatDestroy(&B);
268:     VecDestroy(&maxabsB2);
269:     PetscFree2(imaxabsB, imaxabsB2);
270:   }

272:   VecDestroy(&min);
273:   VecDestroy(&max);
274:   VecDestroy(&maxabs);
275:   VecDestroy(&e);
276:   MatDestroy(&A);
277:   PetscFinalize();
278:   return 0;
279: }

281: /*TEST

283:    test:
284:       output_file: output/ex114.out

286:    test:
287:       suffix: 2
288:       args: -testcase 1
289:       output_file: output/ex114.out

291:    test:
292:       suffix: 3
293:       args: -testcase 2
294:       output_file: output/ex114_3.out

296:    test:
297:       suffix: 4
298:       args: -testcase 3
299:       output_file: output/ex114_4.out

301:    test:
302:       suffix: 5
303:       nsize: 3
304:       args: -testcase 0
305:       output_file: output/ex114_5.out

307:    test:
308:       suffix: 6
309:       nsize: 3
310:       args: -testcase 1
311:       output_file: output/ex114_6.out

313:    test:
314:       suffix: 7
315:       nsize: 3
316:       args: -testcase 2
317:       output_file: output/ex114_7.out

319:    test:
320:       suffix: 8
321:       nsize: 3
322:       args: -testcase 3
323:       output_file: output/ex114_8.out

325: TEST*/