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