Actual source code: ex114.c

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

  3: #include <petscmat.h>

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

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

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

 23:   PetscCall(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:       PetscCall(MatSetSizes(A, M, N, PETSC_DECIDE, PETSC_DECIDE));
 27:     } else {
 28:       PetscCall(MatSetSizes(A, 0, 0, PETSC_DECIDE, PETSC_DECIDE));
 29:     }
 30:     testcase = 0;
 31:   } else {
 32:     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
 33:   }
 34:   PetscCall(MatSetFromOptions(A));
 35:   PetscCall(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;
 45:       values[1] = -2.0;
 46:       values[2] = -2.0;
 47:       values[3] = -4.0;
 48:       values[4] = 1.0;
 49:       values[5] = 1.0;
 50:       PetscCall(MatSetValues(A, 1, &row, N, indices, values, INSERT_VALUES));
 51:       row        = 2;
 52:       indices[0] = 0;
 53:       indices[1] = 3;
 54:       indices[2] = 5;
 55:       values[0]  = -2.0;
 56:       values[1]  = -2.0;
 57:       values[2]  = -2.0;
 58:       PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES));
 59:       row        = 3;
 60:       indices[0] = 0;
 61:       indices[1] = 1;
 62:       indices[2] = 4;
 63:       values[0]  = -2.0;
 64:       values[1]  = -2.0;
 65:       values[2]  = -2.0;
 66:       PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES));
 67:       row        = 4;
 68:       indices[0] = 0;
 69:       indices[1] = 1;
 70:       indices[2] = 2;
 71:       values[0]  = -2.0;
 72:       values[1]  = -2.0;
 73:       values[2]  = -2.0;
 74:       PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES));
 75:       break;
 76:     case 3:
 77:       row       = 0;
 78:       values[0] = -2.0;
 79:       values[1] = -2.0;
 80:       values[2] = -2.0;
 81:       PetscCall(MatSetValues(A, 1, &row, 3, indices + 1, values, INSERT_VALUES));
 82:       row       = 1;
 83:       values[0] = -2.0;
 84:       values[1] = -2.0;
 85:       values[2] = -2.0;
 86:       PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES));
 87:       row       = 2;
 88:       values[0] = -2.0;
 89:       values[1] = -2.0;
 90:       values[2] = -2.0;
 91:       PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES));
 92:       row       = 3;
 93:       values[0] = -2.0;
 94:       values[1] = -2.0;
 95:       values[2] = -2.0;
 96:       values[3] = -1.0;
 97:       PetscCall(MatSetValues(A, 1, &row, 4, indices, values, INSERT_VALUES));
 98:       row       = 4;
 99:       values[0] = -2.0;
100:       values[1] = -2.0;
101:       values[2] = -2.0;
102:       values[3] = -1.0;
103:       PetscCall(MatSetValues(A, 1, &row, 4, indices, values, INSERT_VALUES));
104:       break;

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

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

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

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

150:   PetscCall(VecWAXPY(e, 1.0, max, min)); /* e = max + min */
151:   PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
152:   PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "max+min > PETSC_MACHINE_EPSILON ");
153:   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]);

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

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

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

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

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

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

189:     PetscCall(VecWAXPY(e, 1.0, max, min)); /* e = max + min */
190:     PetscCall(VecNorm(e, NORM_INFINITY, &enorm));
191:     PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "max+min > PETSC_MACHINE_EPSILON ");
192:     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]);

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

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

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

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

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

225:     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]);
226:     PetscCall(MatDestroy(&B));

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

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

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

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

257:     /* Check maxabsB2 and imaxabsB2 */
258:     PetscCall(VecGetArrayRead(maxabsB, &vals));
259:     PetscCall(VecGetArrayRead(maxabsB2, &vals2));
260:     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);
261:     PetscCall(VecRestoreArrayRead(maxabsB, &vals));
262:     PetscCall(VecRestoreArrayRead(maxabsB2, &vals2));

264:     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);
265:     PetscCall(VecDestroy(&maxabsB));
266:     PetscCall(MatDestroy(&B));
267:     PetscCall(VecDestroy(&maxabsB2));
268:     PetscCall(PetscFree2(imaxabsB, imaxabsB2));
269:   }

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

280: /*TEST

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

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

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

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

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

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

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

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

324: TEST*/