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