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:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 20:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 21:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 22:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-testcase", &testcase, NULL));

 24:   PetscCall(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:       PetscCall(MatSetSizes(A, M, N, PETSC_DECIDE, PETSC_DECIDE));
 28:     } else {
 29:       PetscCall(MatSetSizes(A, 0, 0, PETSC_DECIDE, PETSC_DECIDE));
 30:     }
 31:     testcase = 0;
 32:   } else {
 33:     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
 34:   }
 35:   PetscCall(MatSetFromOptions(A));
 36:   PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(MatSetValues(A, 1, &row, N, indices, values, INSERT_VALUES));
116:       row = 1;
117:       PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES));
118:       row = 3;
119:       PetscCall(MatSetValues(A, 1, &row, 1, indices + 4, values + 4, INSERT_VALUES));
120:       row = 4;
121:       PetscCall(MatSetValues(A, 1, &row, 2, indices + 4, values + 4, INSERT_VALUES));
122:     }
123:   }
124:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
125:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
126:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));

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

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

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

151:   PetscCall(VecWAXPY(e, 1.0, max, min)); /* e = max + min */
152:   PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
153:   PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "max+min > PETSC_MACHINE_EPSILON ");
154:   for (j = 0; j < n; j++) PetscCheck(imin[j] == imax[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT, j, imin[j], imax[j]);

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

163:   /* MatGetRowMinAbs() */
164:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMinAbs\n"));
165:   PetscCall(MatGetRowMinAbs(A, min, NULL));
166:   PetscCall(MatGetRowMinAbs(A, min, imin));
167:   PetscCall(VecView(min, PETSC_VIEWER_STDOUT_WORLD));
168:   PetscCall(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:     PetscCall(VecDuplicate(min, &max_d));
175:     PetscCall(VecDuplicate(min, &maxabs_d));

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

182:     PetscCall(VecWAXPY(e, -1.0, max, max_d)); /* e = -max + max_d */
183:     PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
184:     PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "norm(-max + max_d) %g > PETSC_MACHINE_EPSILON", (double)enorm);

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

190:     PetscCall(VecWAXPY(e, 1.0, max, min)); /* e = max + min */
191:     PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
192:     PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "max+min > PETSC_MACHINE_EPSILON ");
193:     for (j = 0; j < n; j++) PetscCheck(imin[j] == imax[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT " for seqdense matrix", j, imin[j], imax[j]);

195:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMaxAbs for seqdense matrix\n"));
196:     PetscCall(MatGetRowMaxAbs(Adense, maxabs_d, imaxabs));
197:     PetscCall(VecWAXPY(e, -1.0, maxabs, maxabs_d)); /* e = -maxabs + maxabs_d */
198:     PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
199:     PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON", (double)enorm);

201:     PetscCall(MatDestroy(&Adense));
202:     PetscCall(VecDestroy(&max_d));
203:     PetscCall(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:     PetscCall(PetscMalloc2(M, &imaxabsB, bs * M, &imaxabsB2));

216:     /* bs = 1 */
217:     PetscCall(MatConvert(A, MATMPIBAIJ, MAT_INITIAL_MATRIX, &B));
218:     PetscCall(VecDuplicate(min, &maxabsB));
219:     PetscCall(MatGetRowMaxAbs(B, maxabsB, NULL));
220:     PetscCall(MatGetRowMaxAbs(B, maxabsB, imaxabsB));
221:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMaxAbs for BAIJ matrix\n"));
222:     PetscCall(VecWAXPY(e, -1.0, maxabs, maxabsB)); /* e = -maxabs + maxabsB */
223:     PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
224:     PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON", (double)enorm);

226:     for (j = 0; j < n; j++) PetscCheck(imaxabs[j] == imaxabsB[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "imaxabs[%" PetscInt_FMT "] %" PetscInt_FMT " != imaxabsB %" PetscInt_FMT, j, imin[j], imax[j]);
227:     PetscCall(MatDestroy(&B));

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

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

241:     for (row = rstart; row < rend; row++) {
242:       PetscCall(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:         PetscCall(MatSetValues(B, bs, Brows, bs, Bcols, Bvals, INSERT_VALUES));
250:       }
251:       PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
252:     }
253:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
254:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

256:     PetscCall(MatGetRowMaxAbs(B, maxabsB2, imaxabsB2));

258:     /* Check maxabsB2 and imaxabsB2 */
259:     PetscCall(VecGetArrayRead(maxabsB, &vals));
260:     PetscCall(VecGetArrayRead(maxabsB2, &vals2));
261:     for (row = 0; row < m; row++) PetscCheck(PetscAbsScalar(vals[row] - vals2[bs * row]) <= PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_PLIB, "row %" PetscInt_FMT " maxabsB != maxabsB2", row);
262:     PetscCall(VecRestoreArrayRead(maxabsB, &vals));
263:     PetscCall(VecRestoreArrayRead(maxabsB2, &vals2));

265:     for (col = 0; col < n; col++) PetscCheck(imaxabsB[col] == imaxabsB2[bs * col] / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "col %" PetscInt_FMT " imaxabsB != imaxabsB2", col);
266:     PetscCall(VecDestroy(&maxabsB));
267:     PetscCall(MatDestroy(&B));
268:     PetscCall(VecDestroy(&maxabsB2));
269:     PetscCall(PetscFree2(imaxabsB, imaxabsB2));
270:   }

272:   PetscCall(VecDestroy(&min));
273:   PetscCall(VecDestroy(&max));
274:   PetscCall(VecDestroy(&maxabs));
275:   PetscCall(VecDestroy(&e));
276:   PetscCall(MatDestroy(&A));
277:   PetscCall(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*/