Actual source code: ex309.c
1: static const char help[] = "Test MATPRODUCT_AB (MatMatMult) with MATDIAGONAL and MATCONSTANTDIAGONAL against any matrix type, and MATPRODUCT_AtB (MatTransposeMatMult) of MATDIAGONAL and MATCONSTANTDIAGONAL with a dense matrix\n\n";
3: // Contributed by: Steven Dargaville
5: #include <petscmat.h>
7: // Compute result = X * Y (MATPRODUCT_AB) and verify it against the action of X and Y.
8: static PetscErrorCode CheckAB(Mat X, Mat Y, const char *what)
9: {
10: Mat result;
11: PetscBool equal = PETSC_FALSE;
13: PetscFunctionBeginUser;
14: PetscCall(MatMatMult(X, Y, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &result));
15: PetscCall(MatMatMultEqual(X, Y, result, 10, &equal));
16: PetscCheck(equal, PetscObjectComm((PetscObject)X), PETSC_ERR_PLIB, "MatMatMult %s gives the wrong result", what);
17: PetscCall(MatDestroy(&result));
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: // Compare C with the host reference Ref and destroy Ref
22: static PetscErrorCode CheckHostReference(Mat C, Mat *Ref, const char *what)
23: {
24: Mat Ch;
25: PetscReal nrm, err;
27: PetscFunctionBeginUser;
28: PetscCall(MatConvert(C, MATDENSE, MAT_INITIAL_MATRIX, &Ch));
29: PetscCall(MatNorm(*Ref, NORM_INFINITY, &nrm));
30: PetscCall(MatAXPY(Ch, -1.0, *Ref, SAME_NONZERO_PATTERN));
31: PetscCall(MatNorm(Ch, NORM_INFINITY, &err));
32: PetscCheck(err <= 100.0 * PETSC_SMALL * PetscMax(nrm, 1.0), PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "MatTransposeMatMult %s gives the wrong result, error %g relative to %g", what, (double)err, (double)nrm);
33: PetscCall(MatDestroy(&Ch));
34: PetscCall(MatDestroy(Ref));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: /* Check C = D^T * B = D * B for a MATDIAGONAL D and a dense B against host copies of D and B. Using host copies also
39: works for a VECKOKKOS diagonal with a CUPM dense matrix, where MatMatMultEqual() fails */
40: static PetscErrorCode CheckDiagonalAtB(Mat D, Mat B, Mat C, const char *what)
41: {
42: Mat Bh;
43: Vec dg, dh;
44: const PetscScalar *d;
45: PetscScalar *h;
46: PetscInt nloc;
48: PetscFunctionBeginUser;
49: PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bh));
50: PetscCall(MatCreateVecs(Bh, NULL, &dh));
51: PetscCall(MatDiagonalGetDiagonal(D, &dg));
52: PetscCall(VecGetLocalSize(dg, &nloc));
53: PetscCall(VecGetArrayRead(dg, &d));
54: PetscCall(VecGetArrayWrite(dh, &h));
55: PetscCall(PetscArraycpy(h, d, nloc));
56: PetscCall(VecRestoreArrayWrite(dh, &h));
57: PetscCall(VecRestoreArrayRead(dg, &d));
58: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
59: PetscCall(MatDiagonalScale(Bh, dh, NULL));
60: PetscCall(VecDestroy(&dh));
61: PetscCall(CheckHostReference(C, &Bh, what));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: // Check C = CD^T * B = alpha * B for a MATCONSTANTDIAGONAL CD and a dense B against a host copy of B
66: static PetscErrorCode CheckConstantDiagonalAtB(Mat CD, Mat B, Mat C, const char *what)
67: {
68: Mat Bh;
69: PetscScalar alpha;
71: PetscFunctionBeginUser;
72: PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bh));
73: PetscCall(MatConstantDiagonalGetConstant(CD, &alpha));
74: PetscCall(MatScale(Bh, alpha));
75: PetscCall(CheckHostReference(C, &Bh, what));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: int main(int argc, char **args)
80: {
81: Mat A, B, D, D2, CD, CD2, result, ref, Cr;
82: Vec dvec, rdiag, dg;
83: MatType atype;
84: PetscInt n = 6, m, ncols = 3, rstart, rend;
85: PetscScalar cval = 1.5, cval2 = 2.5;
86: PetscBool equal = PETSC_FALSE, issame = PETSC_FALSE, check_copies = PETSC_FALSE, dense_ab = PETSC_TRUE;
87: PetscLogEvent event;
89: PetscFunctionBeginUser;
90: PetscCall(PetscInitialize(&argc, &args, NULL, help));
91: PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_copies", &check_copies, NULL));
92: // -dense_ab 0 skips the D * B and CD * B checks, since MatMatMultEqual() fails for a VECKOKKOS diagonal with a CUPM dense matrix
93: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dense_ab", &dense_ab, NULL));
94: if (check_copies) PetscCall(PetscLogDefaultBegin());
95: PetscCall(PetscLogEventRegister("ProductCheck", MAT_CLASSID, &event));
97: // A genuinely non-diagonal sparse matrix with a parallel-safe layout; its type
98: // is taken from -mat_type (e.g. aij, aijkokkos, aijcusparse, aijhipsparse).
99: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
100: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
101: PetscCall(MatSetFromOptions(A));
102: PetscCall(MatSetUp(A));
103: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
104: for (PetscInt i = rstart; i < rend; i++) {
105: PetscInt cols[2] = {i, (i + 1) % n};
106: PetscScalar vals[2] = {(PetscScalar)(i + 2), 1.0};
107: PetscCall(MatSetValues(A, 1, &i, 2, cols, vals, INSERT_VALUES));
108: }
109: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
110: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
111: PetscCall(MatGetType(A, &atype));
112: PetscCall(MatGetLocalSize(A, &m, NULL));
114: // A dense matrix
115: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
116: PetscCall(MatSetType(B, MATDENSE));
117: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, ncols));
118: PetscCall(MatSetOptionsPrefix(B, "dense_"));
119: PetscCall(MatSetFromOptions(B));
120: PetscCall(MatSetUp(B));
121: PetscCall(MatSetRandom(B, NULL));
123: // Two MATDIAGONAL matrices from Vecs matching A's layout and VecType (so on a
124: // device build with a device A everything stays on device).
125: PetscCall(MatCreateVecs(A, &dvec, NULL));
126: PetscCall(VecGetOwnershipRange(dvec, &rstart, &rend));
127: for (PetscInt i = rstart; i < rend; i++) PetscCall(VecSetValue(dvec, i, (PetscScalar)(i + 3), INSERT_VALUES));
128: PetscCall(VecAssemblyBegin(dvec));
129: PetscCall(VecAssemblyEnd(dvec));
130: PetscCall(MatCreateDiagonal(dvec, &D));
131: PetscCall(VecScale(dvec, -0.5));
132: PetscCall(MatCreateDiagonal(dvec, &D2));
133: PetscCall(VecDestroy(&dvec));
135: // Two MATCONSTANTDIAGONAL matrices.
136: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, m, m, n, n, cval, &CD));
137: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, m, m, n, n, cval2, &CD2));
139: // MATDIAGONAL against a general (sparse) matrix, both orientations.
140: PetscCall(CheckAB(A, D, "A * D (aij * MATDIAGONAL)"));
141: PetscCall(CheckAB(D, A, "D * A (MATDIAGONAL * aij)"));
143: // MATCONSTANTDIAGONAL against a general (sparse) matrix, both orientations.
144: PetscCall(CheckAB(A, CD, "A * CD (aij * MATCONSTANTDIAGONAL)"));
145: PetscCall(CheckAB(CD, A, "CD * A (MATCONSTANTDIAGONAL * aij)"));
147: // MATDIAGONAL against a dense matrix
148: if (dense_ab) PetscCall(CheckAB(D, B, "D * B (MATDIAGONAL * dense)"));
149: PetscCall(MatTransposeMatMult(D, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
150: PetscCall(CheckDiagonalAtB(D, B, Cr, "D^T * B (MATDIAGONAL^T * dense)"));
151: PetscCall(MatDestroy(&Cr));
153: // MATCONSTANTDIAGONAL against a dense matrix
154: if (dense_ab) PetscCall(CheckAB(CD, B, "CD * B (MATCONSTANTDIAGONAL * dense)"));
155: PetscCall(MatTransposeMatMult(CD, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
156: PetscCall(CheckConstantDiagonalAtB(CD, B, Cr, "CD^T * B (MATCONSTANTDIAGONAL^T * dense)"));
157: PetscCall(MatDestroy(&Cr));
159: // MATDIAGONAL/MATCONSTANTDIAGONAL against each other.
160: PetscCall(CheckAB(D, D2, "D * D (MATDIAGONAL * MATDIAGONAL)"));
161: PetscCall(CheckAB(CD, CD2, "CD * CD (MATCONSTANTDIAGONAL * MATCONSTANTDIAGONAL)"));
162: PetscCall(CheckAB(CD, D, "CD * D (MATCONSTANTDIAGONAL * MATDIAGONAL)"));
163: PetscCall(CheckAB(D, CD, "D * CD (MATDIAGONAL * MATCONSTANTDIAGONAL)"));
165: // Explicit symbolic-only phase, then numeric (mirrors callers that build the
166: // product structure without an immediate numeric); verify the symbolic result
167: // inherits the general operand's type, then check the numeric values.
168: PetscCall(MatProductCreate(A, D, NULL, &result));
169: PetscCall(MatProductSetType(result, MATPRODUCT_AB));
170: PetscCall(MatProductSetFromOptions(result));
171: PetscCall(MatProductSymbolic(result));
172: PetscCall(PetscObjectTypeCompare((PetscObject)result, atype, &issame));
173: PetscCheck(issame, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Symbolic AB product has an unexpected type");
174: PetscCall(MatProductNumeric(result));
175: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &ref));
176: PetscCall(MatDiagonalGetDiagonal(D, &rdiag));
177: PetscCall(MatDiagonalScale(ref, NULL, rdiag));
178: PetscCall(MatDiagonalRestoreDiagonal(D, &rdiag));
179: PetscCall(MatEqual(result, ref, &equal));
180: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Symbolic+numeric AB product does not match column scaling");
181: PetscCall(MatDestroy(&ref));
182: PetscCall(MatDestroy(&result));
184: // MAT_REUSE_MATRIX: reuse skips the symbolic phase and re-runs only the
185: // numeric, so mutating an operand and reusing must recompute the product.
186: // First the anytype-output orientations (D * A, A * D, A * CD).
187: PetscCall(MatMatMult(D, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
188: PetscCall(MatMatMultEqual(D, A, Cr, 10, &equal));
189: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * A (initial) gives the wrong result");
190: PetscCall(MatDiagonalGetDiagonal(D, &dg));
191: PetscCall(VecScale(dg, -2.5));
192: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
193: PetscCall(MatMatMult(D, A, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
194: PetscCall(MatMatMultEqual(D, A, Cr, 10, &equal));
195: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * A (reuse, diagonal changed) gives the wrong result");
196: PetscCall(MatDestroy(&Cr));
198: PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
199: PetscCall(MatDiagonalGetDiagonal(D, &dg));
200: PetscCall(VecScale(dg, 0.7));
201: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
202: PetscCall(MatMatMult(A, D, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
203: PetscCall(MatMatMultEqual(A, D, Cr, 10, &equal));
204: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A * D (reuse, diagonal changed) gives the wrong result");
205: PetscCall(MatDestroy(&Cr));
207: PetscCall(MatMatMult(A, CD, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
208: PetscCall(MatScale(CD, 2.0));
209: PetscCall(MatMatMult(A, CD, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
210: PetscCall(MatMatMultEqual(A, CD, Cr, 10, &equal));
211: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A * CD (reuse, constant scaled) gives the wrong result");
212: PetscCall(MatDestroy(&Cr));
214: // Then the MATDIAGONAL/MATCONSTANTDIAGONAL-output combinations (D * D, CD * CD,
215: // CD * D, D * CD); their numeric routines likewise recompute from the operands.
216: PetscCall(MatMatMult(D, D2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
217: PetscCall(MatDiagonalGetDiagonal(D, &dg));
218: PetscCall(VecScale(dg, 1.3));
219: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
220: PetscCall(MatMatMult(D, D2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
221: PetscCall(MatMatMultEqual(D, D2, Cr, 10, &equal));
222: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * D (reuse, diagonal changed) gives the wrong result");
223: PetscCall(MatDestroy(&Cr));
225: PetscCall(MatMatMult(CD, CD2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
226: PetscCall(MatScale(CD, 1.5));
227: PetscCall(MatMatMult(CD, CD2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
228: PetscCall(MatMatMultEqual(CD, CD2, Cr, 10, &equal));
229: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "CD * CD (reuse, constant scaled) gives the wrong result");
230: PetscCall(MatDestroy(&Cr));
232: PetscCall(MatMatMult(CD, D, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
233: PetscCall(MatDiagonalGetDiagonal(D, &dg));
234: PetscCall(VecScale(dg, -0.9));
235: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
236: PetscCall(MatMatMult(CD, D, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
237: PetscCall(MatMatMultEqual(CD, D, Cr, 10, &equal));
238: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "CD * D (reuse, diagonal changed) gives the wrong result");
239: PetscCall(MatDestroy(&Cr));
241: PetscCall(MatMatMult(D, CD, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
242: PetscCall(MatScale(CD, 0.5));
243: PetscCall(MatMatMult(D, CD, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
244: PetscCall(MatMatMultEqual(D, CD, Cr, 10, &equal));
245: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * CD (reuse, constant scaled) gives the wrong result");
246: PetscCall(MatDestroy(&Cr));
248: // Reuse D^T * B and CD^T * B after the diagonals change and check them, then reuse them again with the copies logged
249: PetscCall(MatTransposeMatMult(D, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
250: PetscCall(MatTransposeMatMult(CD, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &result));
251: PetscCall(MatDiagonalGetDiagonal(D, &dg));
252: PetscCall(VecScale(dg, 1.7));
253: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
254: PetscCall(MatScale(CD, -3.0));
255: PetscCall(MatTransposeMatMult(D, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
256: PetscCall(CheckDiagonalAtB(D, B, Cr, "D^T * B (reuse, diagonal changed)"));
257: PetscCall(MatTransposeMatMult(CD, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &result));
258: PetscCall(CheckConstantDiagonalAtB(CD, B, result, "CD^T * B (reuse, constant scaled)"));
259: PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0));
260: PetscCall(MatTransposeMatMult(D, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
261: PetscCall(MatTransposeMatMult(CD, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &result));
262: PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0));
263: PetscCall(CheckDiagonalAtB(D, B, Cr, "D^T * B (reuse again)"));
264: PetscCall(CheckConstantDiagonalAtB(CD, B, result, "CD^T * B (reuse again)"));
265: PetscCall(MatDestroy(&Cr));
266: PetscCall(MatDestroy(&result));
267: #if PetscDefined(HAVE_DEVICE)
268: // With device matrices, the logged reuse must not copy between host and device
269: if (check_copies) {
270: PetscEventPerfInfo info;
272: PetscCall(PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &info));
273: PetscCheck(info.GpuToCpuCount == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%g unexpected GPU to CPU copies (%g bytes) in MatTransposeMatMult()", info.GpuToCpuCount, info.GpuToCpuSize);
274: PetscCheck(info.CpuToGpuCount == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%g unexpected CPU to GPU copies (%g bytes) in MatTransposeMatMult()", info.CpuToGpuCount, info.CpuToGpuSize);
275: }
276: #endif
278: PetscCall(MatDestroy(&A));
279: PetscCall(MatDestroy(&B));
280: PetscCall(MatDestroy(&D));
281: PetscCall(MatDestroy(&D2));
282: PetscCall(MatDestroy(&CD));
283: PetscCall(MatDestroy(&CD2));
285: PetscCall(PetscFinalize());
286: return 0;
287: }
289: /*TEST
290: test:
291: suffix: cpu
292: nsize: {{1 2}}
293: output_file: output/empty.out
295: test:
296: requires: kokkos_kernels
297: suffix: kokkos
298: nsize: {{1 2}}
299: args: -mat_type aijkokkos
300: output_file: output/empty.out
302: test:
303: requires: cuda
304: suffix: cuda
305: nsize: {{1 2}}
306: args: -mat_type aijcusparse -dense_mat_type densecuda
307: output_file: output/empty.out
309: test:
310: requires: hip
311: suffix: hip
312: nsize: {{1 2}}
313: args: -mat_type aijhipsparse -dense_mat_type densehip
314: output_file: output/empty.out
316: # the copy checks on 2 processes do not need GPU-aware MPI, since these products do not communicate
317: test:
318: requires: cuda defined(PETSC_USE_LOG)
319: suffix: cuda_copies
320: nsize: {{1 2}}
321: args: -mat_type aijcusparse -dense_mat_type densecuda -check_copies
322: output_file: output/empty.out
324: test:
325: requires: hip defined(PETSC_USE_LOG)
326: suffix: hip_copies
327: nsize: {{1 2}}
328: args: -mat_type aijhipsparse -dense_mat_type densehip -check_copies
329: output_file: output/empty.out
331: test:
332: requires: cuda kokkos_kernels defined(PETSC_USE_LOG)
333: suffix: kokkos_cuda
334: nsize: {{1 2}}
335: args: -mat_type aijkokkos -dense_mat_type densecuda -dense_mat_vec_type kokkos -check_copies -dense_ab 0
336: output_file: output/empty.out
338: test:
339: requires: hip kokkos_kernels defined(PETSC_USE_LOG)
340: suffix: kokkos_hip
341: nsize: {{1 2}}
342: args: -mat_type aijkokkos -dense_mat_type densehip -dense_mat_vec_type kokkos -check_copies -dense_ab 0
343: output_file: output/empty.out
344: TEST*/