Actual source code: ex272.c
1: static const char help[] = "Test MatPtAP with MATDIAGONAL and MATCONSTANTDIAGONAL\n";
3: #include <petscmat.h>
5: /* KOKKOS: Following two cases will fail, as MatDiagonalScale_{Seq,MPI}AIJKOKKOS does
6: not support CPU diagonal vector against AIJ KOKKOS.
7: -a_mat_type diagonal -p_mat_type aijkokkos -a_mat_vec_type standard
8: -a_mat_type aijkokkos -p_mat_type diagonal -p_mat_vec_type standard
9: Pairing a MATDIAGONAL holding a KOKKOS Vec with MATDENSECUDA/MATDENSEHIP also fails,
10: as the MatPtAPMultEqual() verification requires VecCopy() from a CUDA/HIP Vec into a
11: KOKKOS Vec, which is not supported.
12: -a_mat_type diagonal -p_mat_type densecuda -a_mat_vec_type kokkos
13: -a_mat_type densecuda -p_mat_type diagonal -p_mat_vec_type kokkos */
14: static PetscErrorCode CreateTestMatrix(MPI_Comm comm, const char prefix[], PetscInt m, PetscInt n, Mat *M)
15: {
16: PetscFunctionBeginUser;
17: PetscCall(MatCreate(comm, M));
18: PetscCall(MatSetSizes(*M, PETSC_DECIDE, PETSC_DECIDE, m, n));
19: PetscCall(MatSetType(*M, MATAIJ));
20: PetscCall(MatSetOptionsPrefix(*M, prefix));
21: PetscCall(MatSetFromOptions(*M));
22: PetscCall(MatSeqAIJSetPreallocation(*M, n, NULL));
23: PetscCall(MatMPIAIJSetPreallocation(*M, n, NULL, n, NULL));
24: PetscCall(MatSetUp(*M));
25: PetscCall(MatSetRandom(*M, NULL));
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: /* Check correctness of C = P^T A P via mat-vec
30: For AIJ-type P matrix, check if factorization is possible */
31: static PetscErrorCode VerifyPtAP(MPI_Comm comm, Mat A, Mat P, Mat C)
32: {
33: PetscBool flg, can_chol;
35: PetscFunctionBeginUser;
36: PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg));
37: PetscCheck(flg, comm, PETSC_ERR_PLIB, "MatPtAPMultEqual() failed");
39: PetscCall(MatGetFactorAvailable(C, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &can_chol));
40: if (can_chol) {
41: Mat Cshift, fact;
42: IS perm, cperm;
43: MatFactorInfo info;
44: PetscReal norm;
46: PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &Cshift));
47: PetscCall(MatNorm(Cshift, NORM_INFINITY, &norm));
48: PetscCall(MatShift(Cshift, 2.0 * norm + 1.0));
49: PetscCall(MatGetFactor(Cshift, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &fact));
50: PetscCall(MatGetOrdering(Cshift, MATORDERINGNATURAL, &perm, &cperm));
51: PetscCall(MatFactorInfoInitialize(&info));
52: info.fill = 5.0;
53: PetscCall(MatCholeskyFactorSymbolic(fact, Cshift, perm, &info));
54: PetscCall(MatCholeskyFactorNumeric(fact, Cshift, &info));
55: PetscCall(ISDestroy(&perm));
56: PetscCall(ISDestroy(&cperm));
57: PetscCall(MatDestroy(&fact));
58: PetscCall(MatDestroy(&Cshift));
59: }
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: int main(int argc, char **argv)
64: {
65: Mat A, P, C;
66: MPI_Comm comm;
67: PetscInt m = 10, n = 8;
69: PetscFunctionBeginUser;
70: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
71: comm = PETSC_COMM_WORLD;
73: PetscOptionsBegin(comm, "", help, "none");
74: PetscCall(PetscOptionsInt("-m", "m size", "", m, &m, NULL));
75: PetscCall(PetscOptionsInt("-n", "n size", "", n, &n, NULL));
76: PetscOptionsEnd();
78: PetscCall(CreateTestMatrix(comm, "a_", m, m, &A));
79: PetscCall(CreateTestMatrix(comm, "p_", m, n, &P));
81: /* Initial PtAP */
82: PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
83: PetscCall(VerifyPtAP(comm, A, P, C));
85: /* Reuse with modified A */
86: PetscCall(MatScale(A, 2.0));
87: PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
88: PetscCall(VerifyPtAP(comm, A, P, C));
90: /* Reuse with modified P */
91: PetscCall(MatScale(P, 0.5));
92: PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
93: PetscCall(VerifyPtAP(comm, A, P, C));
95: /* Reuse with modified A */
96: PetscCall(MatScale(A, 1.1));
97: PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
98: PetscCall(VerifyPtAP(comm, A, P, C));
100: /* Reuse with modified P */
101: PetscCall(MatScale(P, 3.7));
102: PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
103: PetscCall(VerifyPtAP(comm, A, P, C));
105: /* Modify both A and P */
106: PetscCall(MatScale(A, 0.23));
107: PetscCall(MatScale(P, 1.43));
108: PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
109: PetscCall(VerifyPtAP(comm, A, P, C));
111: PetscCall(MatDestroy(&C));
112: PetscCall(MatDestroy(&P));
113: PetscCall(MatDestroy(&A));
114: PetscCall(PetscFinalize());
115: return 0;
116: }
118: /*TEST
120: test:
121: suffix: a_diag_cpu
122: nsize: {{1 2}}
123: args: -a_mat_type diagonal -p_mat_type {{aij dense}}
124: output_file: output/empty.out
126: test:
127: suffix: p_diag_cpu
128: nsize: {{1 2}}
129: args: -a_mat_type {{diagonal aij dense}} -p_mat_type diagonal -n 10 -m 10
130: output_file: output/empty.out
132: test:
133: suffix: diag_diag_kokkos
134: nsize: {{1 2}}
135: requires: kokkos_kernels
136: args: -a_mat_type diagonal -p_mat_type diagonal -a_mat_vec_type kokkos -p_mat_vec_type {{kokkos standard}} -n 10 -m 10
137: output_file: output/empty.out
139: test:
140: suffix: diag_standard_diag_kokkos
141: nsize: {{1 2}}
142: requires: kokkos_kernels
143: args: -a_mat_type diagonal -p_mat_type diagonal -a_mat_vec_type standard -p_mat_vec_type kokkos -n 10 -m 10
144: output_file: output/empty.out
146: test:
147: suffix: diag_diag_cuda
148: nsize: {{1 2}}
149: requires: cuda
150: args: -a_mat_type diagonal -p_mat_type diagonal -a_mat_vec_type cuda -p_mat_vec_type {{cuda standard}} -n 10 -m 10
151: output_file: output/empty.out
153: test:
154: suffix: diag_standard_diag_cuda
155: nsize: {{1 2}}
156: requires: cuda
157: args: -a_mat_type diagonal -p_mat_type diagonal -a_mat_vec_type standard -p_mat_vec_type cuda -n 10 -m 10
158: output_file: output/empty.out
160: test:
161: suffix: diag_cuda_mat
162: nsize: {{1 2}}
163: requires: cuda
164: args: -a_mat_type diagonal -p_mat_type {{aijcusparse densecuda}} -a_mat_vec_type cuda
165: output_file: output/empty.out
167: test:
168: suffix: cuda_mat_diag
169: nsize: {{1 2}}
170: requires: cuda
171: args: -a_mat_type {{aijcusparse densecuda}} -p_mat_type diagonal -p_mat_vec_type cuda -n 10 -m 10
172: output_file: output/empty.out
174: test:
175: suffix: diag_standard_densecuda
176: nsize: {{1 2}}
177: requires: cuda
178: args: -a_mat_type diagonal -p_mat_type densecuda -a_mat_vec_type standard
179: output_file: output/empty.out
181: test:
182: suffix: densecuda_diag_standard
183: nsize: {{1 2}}
184: requires: cuda
185: args: -a_mat_type densecuda -p_mat_type diagonal -p_mat_vec_type standard -n 10 -m 10
186: output_file: output/empty.out
188: test:
189: suffix: diag_diag_hip
190: nsize: {{1 2}}
191: requires: hip
192: args: -a_mat_type diagonal -p_mat_type diagonal -a_mat_vec_type hip -p_mat_vec_type {{hip standard}} -n 10 -m 10
193: output_file: output/empty.out
195: test:
196: suffix: diag_standard_diag_hip
197: nsize: {{1 2}}
198: requires: hip
199: args: -a_mat_type diagonal -p_mat_type diagonal -a_mat_vec_type standard -p_mat_vec_type hip -n 10 -m 10
200: output_file: output/empty.out
202: test:
203: suffix: diag_hip_mat
204: nsize: {{1 2}}
205: requires: hip
206: args: -a_mat_type diagonal -p_mat_type {{aijhipsparse densehip}} -a_mat_vec_type hip
207: output_file: output/empty.out
209: test:
210: suffix: hip_mat_diag
211: nsize: {{1 2}}
212: requires: hip
213: args: -a_mat_type {{aijhipsparse densehip}} -p_mat_type diagonal -p_mat_vec_type hip -n 10 -m 10
214: output_file: output/empty.out
216: test:
217: suffix: diag_standard_densehip
218: nsize: {{1 2}}
219: requires: hip
220: args: -a_mat_type diagonal -p_mat_type densehip -a_mat_vec_type standard
221: output_file: output/empty.out
223: test:
224: suffix: densehip_diag_standard
225: nsize: {{1 2}}
226: requires: hip
227: args: -a_mat_type densehip -p_mat_type diagonal -p_mat_vec_type standard -n 10 -m 10
228: output_file: output/empty.out
230: test:
231: suffix: a_cdiag_p_dense_cpu
232: nsize: {{1 2}}
233: args: -a_mat_type constantdiagonal -p_mat_type {{aij dense}}
234: output_file: output/empty.out
236: test:
237: suffix: a_cdiag_p_diag_cpu
238: nsize: {{1 2}}
239: args: -a_mat_type constantdiagonal -p_mat_type diagonal -n 10 -m 10
240: output_file: output/empty.out
242: test:
243: suffix: p_cdiag_cpu
244: nsize: {{1 2}}
245: args: -a_mat_type {{constantdiagonal diagonal aij dense}} -p_mat_type constantdiagonal -n 10 -m 10
246: output_file: output/empty.out
248: test:
249: suffix: cdiag_diag_cuda
250: nsize: {{1 2}}
251: requires: cuda
252: args: -a_mat_type constantdiagonal -p_mat_type diagonal -p_mat_vec_type cuda -n 10 -m 10
253: output_file: output/empty.out
255: test:
256: suffix: diag_cuda_cdiag
257: nsize: {{1 2}}
258: requires: cuda
259: args: -a_mat_type diagonal -a_mat_vec_type cuda -p_mat_type constantdiagonal -n 10 -m 10
260: output_file: output/empty.out
262: test:
263: suffix: cdiag_cuda_mat
264: nsize: {{1 2}}
265: requires: cuda
266: args: -a_mat_type constantdiagonal -p_mat_type {{aijcusparse densecuda}}
267: output_file: output/empty.out
269: test:
270: suffix: cuda_mat_cdiag
271: nsize: {{1 2}}
272: requires: cuda
273: args: -a_mat_type {{aijcusparse densecuda}} -p_mat_type constantdiagonal -n 10 -m 10
274: output_file: output/empty.out
276: test:
277: suffix: cdiag_diag_kokkos
278: nsize: {{1 2}}
279: requires: kokkos_kernels
280: args: -a_mat_type constantdiagonal -p_mat_type diagonal -p_mat_vec_type kokkos -n 10 -m 10
281: output_file: output/empty.out
283: test:
284: suffix: diag_kokkos_cdiag
285: nsize: {{1 2}}
286: requires: kokkos_kernels
287: args: -a_mat_type diagonal -a_mat_vec_type kokkos -p_mat_type constantdiagonal -n 10 -m 10
288: output_file: output/empty.out
290: test:
291: suffix: cdiag_diag_hip
292: nsize: {{1 2}}
293: requires: hip
294: args: -a_mat_type constantdiagonal -p_mat_type diagonal -p_mat_vec_type hip -n 10 -m 10
295: output_file: output/empty.out
297: test:
298: suffix: diag_hip_cdiag
299: nsize: {{1 2}}
300: requires: hip
301: args: -a_mat_type diagonal -a_mat_vec_type hip -p_mat_type constantdiagonal -n 10 -m 10
302: output_file: output/empty.out
304: test:
305: suffix: cdiag_hip_mat
306: nsize: {{1 2}}
307: requires: hip
308: args: -a_mat_type constantdiagonal -p_mat_type {{aijhipsparse densehip}}
309: output_file: output/empty.out
311: test:
312: suffix: hip_mat_cdiag
313: nsize: {{1 2}}
314: requires: hip
315: args: -a_mat_type {{aijhipsparse densehip}} -p_mat_type constantdiagonal -n 10 -m 10
316: output_file: output/empty.out
318: TEST*/