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*/