1 static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";
2
3 #include <petscmat.h>
4
TransposeAXPY(Mat C,PetscScalar alpha,Mat mat,PetscErrorCode (* f)(Mat,Mat *))5 static PetscErrorCode TransposeAXPY(Mat C, PetscScalar alpha, Mat mat, PetscErrorCode (*f)(Mat, Mat *))
6 {
7 Mat D, E, F, G;
8 MatType mtype;
9
10 PetscFunctionBegin;
11 PetscCall(MatGetType(mat, &mtype));
12 if (f == MatCreateTranspose) {
13 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY: (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
14 } else {
15 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY: (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
16 }
17 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
18 PetscCall(f(C, &D));
19 PetscCall(f(D, &E));
20 PetscCall(MatAXPY(E, alpha, mat, SAME_NONZERO_PATTERN));
21 PetscCall(MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E));
22 PetscCall(MatView(E, PETSC_VIEWER_STDOUT_WORLD));
23 PetscCall(MatDestroy(&E));
24 PetscCall(MatDestroy(&D));
25 PetscCall(MatDestroy(&C));
26 if (f == MatCreateTranspose) {
27 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n"));
28 } else {
29 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n"));
30 }
31 if (f == MatCreateTranspose) {
32 PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &D));
33 } else {
34 PetscCall(MatHermitianTranspose(mat, MAT_INITIAL_MATRIX, &D));
35 }
36 PetscCall(f(D, &E));
37 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
38 PetscCall(MatAXPY(C, alpha, E, SAME_NONZERO_PATTERN));
39 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
40 PetscCall(MatDestroy(&E));
41 PetscCall(MatDestroy(&D));
42 PetscCall(MatDestroy(&C));
43 if (f == MatCreateTranspose) {
44 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n"));
45 } else {
46 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n"));
47 }
48 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
49 PetscCall(f(C, &D));
50 PetscCall(f(D, &E));
51 PetscCall(f(mat, &F));
52 PetscCall(f(F, &G));
53 PetscCall(MatAXPY(E, alpha, G, SAME_NONZERO_PATTERN));
54 PetscCall(MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E));
55 PetscCall(MatView(E, PETSC_VIEWER_STDOUT_WORLD));
56 PetscCall(MatDestroy(&G));
57 PetscCall(MatDestroy(&F));
58 PetscCall(MatDestroy(&E));
59 PetscCall(MatDestroy(&D));
60 PetscCall(MatDestroy(&C));
61
62 /* Call f on a matrix that does not implement the transposition */
63 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: Now without the transposition operation\n"));
64 PetscCall(MatConvert(mat, MATSHELL, MAT_INITIAL_MATRIX, &C));
65 PetscCall(f(C, &D));
66 PetscCall(f(D, &E));
67 /* XXX cannot use MAT_INPLACE_MATRIX, it leaks mat */
68 PetscCall(MatConvert(E, mtype, MAT_INITIAL_MATRIX, &F));
69 PetscCall(MatAXPY(F, alpha, mat, SAME_NONZERO_PATTERN));
70 PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
71 PetscCall(MatDestroy(&F));
72 PetscCall(MatDestroy(&E));
73 PetscCall(MatDestroy(&D));
74 PetscCall(MatDestroy(&C));
75 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
main(int argc,char ** argv)78 int main(int argc, char **argv)
79 {
80 Mat mat, tmat = 0;
81 PetscInt m = 7, n, i, j, rstart, rend, rect = 0;
82 PetscMPIInt size, rank;
83 PetscBool flg;
84 PetscScalar v, alpha;
85 PetscReal normf, normi, norm1;
86
87 PetscFunctionBeginUser;
88 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
89 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
90 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
91 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
92 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
93 n = m;
94 PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg));
95 if (flg) {
96 n += 2;
97 rect = 1;
98 }
99 PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg));
100 if (flg) {
101 n -= 2;
102 rect = 1;
103 }
104
105 /* ------- Assemble matrix --------- */
106 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
107 PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
108 PetscCall(MatSetFromOptions(mat));
109 PetscCall(MatSetUp(mat));
110 PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
111 for (i = rstart; i < rend; i++) {
112 for (j = 0; j < n; j++) {
113 v = 10.0 * i + j + 1.0;
114 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
115 }
116 }
117 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
118 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
119
120 /* ----------------- Test MatNorm() ----------------- */
121 PetscCall(MatNorm(mat, NORM_FROBENIUS, &normf));
122 PetscCall(MatNorm(mat, NORM_1, &norm1));
123 PetscCall(MatNorm(mat, NORM_INFINITY, &normi));
124 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original A: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
125 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
126
127 /* --------------- Test MatTranspose() -------------- */
128 PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
129 if (!rect && flg) {
130 PetscCall(MatTranspose(mat, MAT_REUSE_MATRIX, &mat)); /* in-place transpose */
131 tmat = mat;
132 mat = NULL;
133 } else { /* out-of-place transpose */
134 PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
135 }
136
137 /* ----------------- Test MatNorm() ----------------- */
138 /* Print info about transpose matrix */
139 PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
140 PetscCall(MatNorm(tmat, NORM_1, &norm1));
141 PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
142 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
143 PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
144
145 /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
146 if (mat && !rect) {
147 alpha = 1.0;
148 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
149 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A\n"));
150 PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
151 PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
152
153 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAYPX: B = alpha*B + A\n"));
154 PetscCall(MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
155 PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
156 }
157
158 {
159 Mat C;
160 alpha = 1.0;
161 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
162 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
163 PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
164 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
165 PetscCall(MatDestroy(&C));
166 PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
167 PetscCall(TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose));
168 }
169
170 {
171 Mat matB;
172 /* get matB that has nonzeros of mat in all even numbers of row and col */
173 PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
174 PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
175 PetscCall(MatSetFromOptions(matB));
176 PetscCall(MatSetUp(matB));
177 PetscCall(MatGetOwnershipRange(matB, &rstart, &rend));
178 if (rstart % 2 != 0) rstart++;
179 for (i = rstart; i < rend; i += 2) {
180 for (j = 0; j < n; j += 2) {
181 v = 10.0 * i + j + 1.0;
182 PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
183 }
184 }
185 PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
186 PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
187 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
188 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
189 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
190 PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
191 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
192 PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
193 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
194 PetscCall(MatDestroy(&matB));
195 }
196
197 /* Test MatZeroRows */
198 j = rstart - 1;
199 if (j < 0) j = m - 1;
200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
201 PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
202 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
203
204 /* Test MatShift */
205 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
206 PetscCall(MatShift(mat, -2.0));
207 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
208
209 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
210 /* Free data structures */
211 PetscCall(MatDestroy(&mat));
212 PetscCall(MatDestroy(&tmat));
213 PetscCall(PetscFinalize());
214 return 0;
215 }
216
217 /*TEST
218
219 test:
220 suffix: 11_A
221 args: -mat_type seqaij -rectA
222 filter: grep -v "Mat Object"
223
224 test:
225 suffix: 12_A
226 args: -mat_type seqdense -rectA
227 filter: grep -v type | grep -v "Mat Object"
228
229 test:
230 requires: cuda
231 suffix: 12_A_cuda
232 args: -mat_type seqdensecuda -rectA
233 output_file: output/ex2_12_A.out
234 filter: grep -v type | grep -v "Mat Object"
235
236 test:
237 requires: kokkos_kernels
238 suffix: 12_A_kokkos
239 args: -mat_type aijkokkos -rectA
240 output_file: output/ex2_12_A.out
241 filter: grep -v type | grep -v "Mat Object"
242
243 test:
244 suffix: 11_B
245 args: -mat_type seqaij -rectB
246 filter: grep -v "Mat Object"
247
248 test:
249 suffix: 12_B
250 args: -mat_type seqdense -rectB
251 filter: grep -v type | grep -v "Mat Object"
252
253 testset:
254 args: -rectB
255 output_file: output/ex2_12_B.out
256 filter: grep -v type | grep -v "Mat Object"
257
258 test:
259 requires: cuda
260 suffix: 12_B_cuda
261 args: -mat_type {{seqdensecuda seqaijcusparse}}
262
263 test:
264 requires: kokkos_kernels
265 suffix: 12_B_kokkos
266 args: -mat_type aijkokkos
267
268 test:
269 suffix: 12_B_aij
270 args: -mat_type aij
271 test:
272 suffix: 21
273 args: -mat_type mpiaij
274 filter: grep -v type | grep -v " MPI process"
275
276 test:
277 suffix: 22
278 args: -mat_type mpidense
279 filter: grep -v type | grep -v "Mat Object"
280
281 test:
282 requires: cuda
283 suffix: 22_cuda
284 output_file: output/ex2_22.out
285 args: -mat_type mpidensecuda
286 filter: grep -v type | grep -v "Mat Object"
287
288 test:
289 requires: kokkos_kernels
290 suffix: 22_kokkos
291 output_file: output/ex2_22.out
292 args: -mat_type aijkokkos
293 filter: grep -v type | grep -v "Mat Object"
294
295 test:
296 suffix: 23
297 nsize: 3
298 args: -mat_type mpiaij
299 filter: grep -v type | grep -v " MPI process"
300
301 test:
302 suffix: 24
303 nsize: 3
304 args: -mat_type mpidense
305 filter: grep -v type | grep -v "Mat Object"
306
307 test:
308 requires: cuda
309 suffix: 24_cuda
310 nsize: 3
311 output_file: output/ex2_24.out
312 args: -mat_type mpidensecuda
313 filter: grep -v type | grep -v "Mat Object"
314
315 test:
316 suffix: 2_aijcusparse_1
317 args: -mat_type mpiaijcusparse
318 output_file: output/ex2_21.out
319 requires: cuda
320 filter: grep -v type | grep -v " MPI process"
321
322 test:
323 suffix: 2_aijkokkos_1
324 args: -mat_type aijkokkos
325 output_file: output/ex2_21.out
326 requires: kokkos_kernels
327 filter: grep -v type | grep -v " MPI process"
328
329 test:
330 suffix: 2_aijcusparse_2
331 nsize: 3
332 args: -mat_type mpiaijcusparse
333 output_file: output/ex2_23.out
334 requires: cuda
335 filter: grep -v type | grep -v " MPI process"
336
337 test:
338 suffix: 2_aijkokkos_2
339 nsize: 3
340 args: -mat_type aijkokkos
341 output_file: output/ex2_23.out
342 # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
343 requires: !hip kokkos_kernels
344 filter: grep -v type | grep -v "MPI processes"
345
346 test:
347 suffix: 3
348 nsize: 2
349 args: -mat_type mpiaij -rectA
350
351 test:
352 suffix: 3_aijcusparse
353 nsize: 2
354 args: -mat_type mpiaijcusparse -rectA
355 requires: cuda
356
357 test:
358 suffix: 4
359 nsize: 2
360 args: -mat_type mpidense -rectA
361 filter: grep -v type | grep -v " MPI process"
362
363 test:
364 requires: cuda
365 suffix: 4_cuda
366 nsize: 2
367 output_file: output/ex2_4.out
368 args: -mat_type mpidensecuda -rectA
369 filter: grep -v type | grep -v " MPI process"
370
371 test:
372 suffix: aijcusparse_1
373 args: -mat_type seqaijcusparse -rectA
374 filter: grep -v "Mat Object"
375 output_file: output/ex2_11_A_aijcusparse.out
376 requires: cuda
377
378 TEST*/
379