Actual source code: ex309.c
1: static const char help[] = "Test MATPRODUCT_AB (MatMatMult) with MATDIAGONAL and MATCONSTANTDIAGONAL against any matrix type\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: int main(int argc, char **args)
22: {
23: Mat A, B, D, D2, CD, CD2, result, ref, Cr;
24: Vec dvec, rdiag, dg;
25: MatType atype;
26: PetscInt n = 6, m, ncols = 3, rstart, rend;
27: PetscScalar cval = 1.5, cval2 = 2.5;
28: PetscBool equal = PETSC_FALSE, issame = PETSC_FALSE;
30: PetscFunctionBeginUser;
31: PetscCall(PetscInitialize(&argc, &args, NULL, help));
33: // A genuinely non-diagonal sparse matrix with a parallel-safe layout; its type
34: // is taken from -mat_type (e.g. aij, aijkokkos, aijcusparse, aijhipsparse).
35: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
36: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
37: PetscCall(MatSetFromOptions(A));
38: PetscCall(MatSetUp(A));
39: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
40: for (PetscInt i = rstart; i < rend; i++) {
41: PetscInt cols[2] = {i, (i + 1) % n};
42: PetscScalar vals[2] = {(PetscScalar)(i + 2), 1.0};
43: PetscCall(MatSetValues(A, 1, &i, 2, cols, vals, INSERT_VALUES));
44: }
45: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
46: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
47: PetscCall(MatGetType(A, &atype));
48: PetscCall(MatGetLocalSize(A, &m, NULL));
50: // A dense matrix
51: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
52: PetscCall(MatSetType(B, MATDENSE));
53: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, ncols));
54: PetscCall(MatSetOptionsPrefix(B, "dense_"));
55: PetscCall(MatSetFromOptions(B));
56: PetscCall(MatSetUp(B));
57: PetscCall(MatSetRandom(B, NULL));
59: // Two MATDIAGONAL matrices from Vecs matching A's layout and VecType (so on a
60: // device build with a device A everything stays on device).
61: PetscCall(MatCreateVecs(A, &dvec, NULL));
62: PetscCall(VecGetOwnershipRange(dvec, &rstart, &rend));
63: for (PetscInt i = rstart; i < rend; i++) PetscCall(VecSetValue(dvec, i, (PetscScalar)(i + 3), INSERT_VALUES));
64: PetscCall(VecAssemblyBegin(dvec));
65: PetscCall(VecAssemblyEnd(dvec));
66: PetscCall(MatCreateDiagonal(dvec, &D));
67: PetscCall(VecScale(dvec, -0.5));
68: PetscCall(MatCreateDiagonal(dvec, &D2));
69: PetscCall(VecDestroy(&dvec));
71: // Two MATCONSTANTDIAGONAL matrices.
72: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, m, m, n, n, cval, &CD));
73: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, m, m, n, n, cval2, &CD2));
75: // MATDIAGONAL against a general (sparse) matrix, both orientations.
76: PetscCall(CheckAB(A, D, "A * D (aij * MATDIAGONAL)"));
77: PetscCall(CheckAB(D, A, "D * A (MATDIAGONAL * aij)"));
79: // MATCONSTANTDIAGONAL against a general (sparse) matrix, both orientations.
80: PetscCall(CheckAB(A, CD, "A * CD (aij * MATCONSTANTDIAGONAL)"));
81: PetscCall(CheckAB(CD, A, "CD * A (MATCONSTANTDIAGONAL * aij)"));
83: // MATDIAGONAL against a dense matrix
84: PetscCall(CheckAB(D, B, "D * B (MATDIAGONAL * dense)"));
86: // MATCONSTANTDIAGONAL against a dense matrix
87: PetscCall(CheckAB(CD, B, "CD * B (MATCONSTANTDIAGONAL * dense)"));
89: // MATDIAGONAL/MATCONSTANTDIAGONAL against each other.
90: PetscCall(CheckAB(D, D2, "D * D (MATDIAGONAL * MATDIAGONAL)"));
91: PetscCall(CheckAB(CD, CD2, "CD * CD (MATCONSTANTDIAGONAL * MATCONSTANTDIAGONAL)"));
92: PetscCall(CheckAB(CD, D, "CD * D (MATCONSTANTDIAGONAL * MATDIAGONAL)"));
93: PetscCall(CheckAB(D, CD, "D * CD (MATDIAGONAL * MATCONSTANTDIAGONAL)"));
95: // Explicit symbolic-only phase, then numeric (mirrors callers that build the
96: // product structure without an immediate numeric); verify the symbolic result
97: // inherits the general operand's type, then check the numeric values.
98: PetscCall(MatProductCreate(A, D, NULL, &result));
99: PetscCall(MatProductSetType(result, MATPRODUCT_AB));
100: PetscCall(MatProductSetFromOptions(result));
101: PetscCall(MatProductSymbolic(result));
102: PetscCall(PetscObjectTypeCompare((PetscObject)result, atype, &issame));
103: PetscCheck(issame, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Symbolic AB product has an unexpected type");
104: PetscCall(MatProductNumeric(result));
105: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &ref));
106: PetscCall(MatDiagonalGetDiagonal(D, &rdiag));
107: PetscCall(MatDiagonalScale(ref, NULL, rdiag));
108: PetscCall(MatDiagonalRestoreDiagonal(D, &rdiag));
109: PetscCall(MatEqual(result, ref, &equal));
110: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Symbolic+numeric AB product does not match column scaling");
111: PetscCall(MatDestroy(&ref));
112: PetscCall(MatDestroy(&result));
114: // MAT_REUSE_MATRIX: reuse skips the symbolic phase and re-runs only the
115: // numeric, so mutating an operand and reusing must recompute the product.
116: // First the anytype-output orientations (D * A, A * D, A * CD).
117: PetscCall(MatMatMult(D, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
118: PetscCall(MatMatMultEqual(D, A, Cr, 10, &equal));
119: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * A (initial) gives the wrong result");
120: PetscCall(MatDiagonalGetDiagonal(D, &dg));
121: PetscCall(VecScale(dg, -2.5));
122: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
123: PetscCall(MatMatMult(D, A, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
124: PetscCall(MatMatMultEqual(D, A, Cr, 10, &equal));
125: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * A (reuse, diagonal changed) gives the wrong result");
126: PetscCall(MatDestroy(&Cr));
128: PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
129: PetscCall(MatDiagonalGetDiagonal(D, &dg));
130: PetscCall(VecScale(dg, 0.7));
131: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
132: PetscCall(MatMatMult(A, D, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
133: PetscCall(MatMatMultEqual(A, D, Cr, 10, &equal));
134: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A * D (reuse, diagonal changed) gives the wrong result");
135: PetscCall(MatDestroy(&Cr));
137: PetscCall(MatMatMult(A, CD, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
138: PetscCall(MatScale(CD, 2.0));
139: PetscCall(MatMatMult(A, CD, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
140: PetscCall(MatMatMultEqual(A, CD, Cr, 10, &equal));
141: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A * CD (reuse, constant scaled) gives the wrong result");
142: PetscCall(MatDestroy(&Cr));
144: // Then the MATDIAGONAL/MATCONSTANTDIAGONAL-output combinations (D * D, CD * CD,
145: // CD * D, D * CD); their numeric routines likewise recompute from the operands.
146: PetscCall(MatMatMult(D, D2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
147: PetscCall(MatDiagonalGetDiagonal(D, &dg));
148: PetscCall(VecScale(dg, 1.3));
149: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
150: PetscCall(MatMatMult(D, D2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
151: PetscCall(MatMatMultEqual(D, D2, Cr, 10, &equal));
152: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * D (reuse, diagonal changed) gives the wrong result");
153: PetscCall(MatDestroy(&Cr));
155: PetscCall(MatMatMult(CD, CD2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
156: PetscCall(MatScale(CD, 1.5));
157: PetscCall(MatMatMult(CD, CD2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
158: PetscCall(MatMatMultEqual(CD, CD2, Cr, 10, &equal));
159: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "CD * CD (reuse, constant scaled) gives the wrong result");
160: PetscCall(MatDestroy(&Cr));
162: PetscCall(MatMatMult(CD, D, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
163: PetscCall(MatDiagonalGetDiagonal(D, &dg));
164: PetscCall(VecScale(dg, -0.9));
165: PetscCall(MatDiagonalRestoreDiagonal(D, &dg));
166: PetscCall(MatMatMult(CD, D, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
167: PetscCall(MatMatMultEqual(CD, D, Cr, 10, &equal));
168: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "CD * D (reuse, diagonal changed) gives the wrong result");
169: PetscCall(MatDestroy(&Cr));
171: PetscCall(MatMatMult(D, CD, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cr));
172: PetscCall(MatScale(CD, 0.5));
173: PetscCall(MatMatMult(D, CD, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cr));
174: PetscCall(MatMatMultEqual(D, CD, Cr, 10, &equal));
175: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D * CD (reuse, constant scaled) gives the wrong result");
176: PetscCall(MatDestroy(&Cr));
178: PetscCall(MatDestroy(&A));
179: PetscCall(MatDestroy(&B));
180: PetscCall(MatDestroy(&D));
181: PetscCall(MatDestroy(&D2));
182: PetscCall(MatDestroy(&CD));
183: PetscCall(MatDestroy(&CD2));
185: PetscCall(PetscFinalize());
186: return 0;
187: }
189: /*TEST
190: test:
191: suffix: cpu
192: nsize: {{1 2}}
193: output_file: output/empty.out
195: test:
196: requires: kokkos_kernels
197: suffix: kokkos
198: nsize: {{1 2}}
199: args: -mat_type aijkokkos
200: output_file: output/empty.out
202: test:
203: requires: cuda
204: suffix: cuda
205: nsize: {{1 2}}
206: args: -mat_type aijcusparse -dense_mat_type densecuda
207: output_file: output/empty.out
209: test:
210: requires: hip
211: suffix: hip
212: nsize: {{1 2}}
213: args: -mat_type aijhipsparse -dense_mat_type densehip
214: output_file: output/empty.out
215: TEST*/