1 static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n";
2
3 #include <petscmat.h>
4
5 #define M 5
6 #define N 6
7
main(int argc,char ** args)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;
16
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));
22
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));
36
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;
105
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));
126
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));
134
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));
142
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));
149
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]);
154
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));
161
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));
168
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));
175
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));
180
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);
184
185 PetscCall(MatScale(Adense, -1.0));
186 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMin for seqdense matrix\n"));
187 PetscCall(MatGetRowMin(Adense, min, imin));
188
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]);
193
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);
199
200 PetscCall(MatDestroy(&Adense));
201 PetscCall(VecDestroy(&max_d));
202 PetscCall(VecDestroy(&maxabs_d));
203 }
204
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];
212
213 PetscCall(PetscMalloc2(M, &imaxabsB, bs * M, &imaxabsB2));
214
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);
224
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));
227
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));
232
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));
239
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));
254
255 PetscCall(MatGetRowMaxAbs(B, maxabsB2, imaxabsB2));
256
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));
263
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 }
270
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 }
279
280 /*TEST
281
282 test:
283 output_file: output/ex114.out
284
285 test:
286 suffix: 2
287 args: -testcase 1
288 output_file: output/ex114.out
289
290 test:
291 suffix: 3
292 args: -testcase 2
293 output_file: output/ex114_3.out
294
295 test:
296 suffix: 4
297 args: -testcase 3
298 output_file: output/ex114_4.out
299
300 test:
301 suffix: 5
302 nsize: 3
303 args: -testcase 0
304 output_file: output/ex114_5.out
305
306 test:
307 suffix: 6
308 nsize: 3
309 args: -testcase 1
310 output_file: output/ex114_6.out
311
312 test:
313 suffix: 7
314 nsize: 3
315 args: -testcase 2
316 output_file: output/ex114_7.out
317
318 test:
319 suffix: 8
320 nsize: 3
321 args: -testcase 3
322 output_file: output/ex114_8.out
323
324 TEST*/
325