Actual source code: ex318.c
1: static const char help[] = "Tests MatDenseGetColumnVec() and friends on dense matrices created from a VecType, and on the Mats a MatProduct and MatDenseGetSubMatrix() derive from them\n\n";
3: #include <petscmat.h>
5: /* The Mat a product creates only keeps the VecType when it is the same kind of Mat as the one it is built from */
6: static PetscErrorCode CheckProductVecType(Mat C, Mat A, const char *what)
7: {
8: MatType amtype, cmtype;
9: VecType avtype, cvtype;
10: PetscBool same;
12: PetscFunctionBeginUser;
13: PetscCall(MatGetType(A, &amtype));
14: PetscCall(MatGetType(C, &cmtype));
15: PetscCall(PetscStrcmp(amtype, cmtype, &same));
16: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
17: PetscCall(MatGetVecType(A, &avtype));
18: PetscCall(MatGetVecType(C, &cvtype));
19: PetscCall(PetscStrcmp(avtype, cvtype, &same));
20: PetscCheck(same, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "%s has VecType %s, expected %s", what, cvtype, avtype);
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: /* A dense Mat of the MatType of A, which has the default VecType of that MatType */
25: static PetscErrorCode CreateDenseDefaultVecType(Mat A, PetscInt M, PetscInt N, Mat *B)
26: {
27: MatType mtype;
29: PetscFunctionBeginUser;
30: PetscCall(MatGetType(A, &mtype));
31: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
32: PetscCall(MatSetSizes(*B, PETSC_DECIDE, PETSC_DECIDE, M, N));
33: PetscCall(MatSetType(*B, mtype));
34: PetscCall(MatSetUp(*B));
35: PetscCall(MatZeroEntries(*B));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: int main(int argc, char **argv)
40: {
41: Mat A, C, P, S, S2, C2, C3, C4, C5, C6, D, E, F;
42: Vec v, w;
43: char vtype[64] = VECSTANDARD;
44: VecType avtype, pvtype;
45: PetscBool same;
46: PetscInt M = 9, N = 3, lda, rstart, rend, i, j;
47: PetscReal norm;
48: PetscScalar sum;
49: const PetscScalar *array;
51: PetscFunctionBeginUser;
52: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
53: PetscCall(PetscOptionsGetString(NULL, NULL, "-vec_type", vtype, sizeof(vtype), NULL));
54: /* A VECKOKKOS type gives a MATDENSECUDA/MATDENSEHIP with VECKOKKOS vectors when the Kokkos backend runs on that device */
55: PetscCall(MatCreateDenseFromVecType(PETSC_COMM_WORLD, vtype, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, &A));
56: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
58: /* Write column j through its column vector, then scale it through a read/write column vector: A(:, j) = 2 (j + 1) */
59: for (j = 0; j < N; j++) {
60: PetscCall(MatDenseGetColumnVecWrite(A, j, &v));
61: PetscCall(VecSet(v, (PetscScalar)(j + 1)));
62: PetscCall(MatDenseRestoreColumnVecWrite(A, j, &v));
63: }
64: for (j = 0; j < N; j++) {
65: PetscCall(MatDenseGetColumnVec(A, j, &v));
66: PetscCall(VecScale(v, 2.0));
67: PetscCall(MatDenseRestoreColumnVec(A, j, &v));
68: }
70: /* Check through read-only column vectors, and independently through the matrix itself */
71: for (j = 0; j < N; j++) {
72: PetscCall(MatDenseGetColumnVecRead(A, j, &v));
73: PetscCall(VecSum(v, &sum));
74: PetscCheck(PetscAbsScalar(sum - 2.0 * (j + 1) * M) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Column %" PetscInt_FMT " read back through MatDenseGetColumnVecRead() sums to %g, expected %g", j, (double)PetscRealPart(sum), 2.0 * (j + 1) * M);
75: PetscCall(MatDenseRestoreColumnVecRead(A, j, &v));
76: }
77: PetscCall(MatNorm(A, NORM_INFINITY, &norm));
78: PetscCheck(PetscAbsReal(norm - N * (N + 1)) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatNorm() gives %g, expected %g", (double)norm, (double)(N * (N + 1)));
79: PetscCall(MatDenseGetLDA(A, &lda));
80: PetscCall(MatDenseGetArrayRead(A, &array));
81: for (j = 0; j < N; j++) {
82: for (i = 0; i < rend - rstart; i++)
83: PetscCheck(array[i + j * lda] == 2.0 * (j + 1), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Entry (%" PetscInt_FMT ", %" PetscInt_FMT ") is %g, expected %g", rstart + i, j, (double)PetscRealPart(array[i + j * lda]), 2.0 * (j + 1));
84: }
85: PetscCall(MatDenseRestoreArrayRead(A, &array));
87: /* The Mat a MatProduct creates must have the VecType of the Mat it is built from */
88: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, rend - rstart, rend - rstart, M, M, 3.0, &S));
89: PetscCall(MatMatMult(S, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
90: PetscCall(MatGetVecType(A, &avtype));
91: PetscCall(CheckProductVecType(C, A, "The Mat the product created"));
93: /* C is 3 A, so the columns cancel, which needs both column Vecs to be of the same type */
94: for (j = 0; j < N; j++) {
95: PetscCall(MatDenseGetColumnVec(C, j, &v));
96: PetscCall(MatDenseGetColumnVecRead(A, j, &w));
97: PetscCall(VecAXPY(v, -3.0, w));
98: PetscCall(VecNorm(v, NORM_INFINITY, &norm));
99: PetscCall(MatDenseRestoreColumnVecRead(A, j, &w));
100: PetscCall(MatDenseRestoreColumnVec(C, j, &v));
101: PetscCheck(norm < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Column %" PetscInt_FMT " of the product differs from three times the column it was built from by %g", j, (double)norm);
102: }
104: /* A submatrix must have the VecType of its parent Mat */
105: PetscCall(MatDenseGetSubMatrix(C, PETSC_DECIDE, PETSC_DECIDE, 1, N, &P));
106: PetscCall(MatGetVecType(P, &pvtype));
107: PetscCall(PetscStrcmp(avtype, pvtype, &same));
108: PetscCheck(same, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The submatrix has VecType %s, expected %s", pvtype, avtype);
110: /* The columns of C are zero after the loop above, so adding back a column of A recovers it */
111: for (j = 0; j < N - 1; j++) {
112: PetscCall(MatDenseGetColumnVec(P, j, &v));
113: PetscCall(MatDenseGetColumnVecRead(A, j + 1, &w));
114: PetscCall(VecAXPY(v, 1.0, w));
115: PetscCall(VecNorm(v, NORM_INFINITY, &norm));
116: PetscCall(MatDenseRestoreColumnVecRead(A, j + 1, &w));
117: PetscCall(MatDenseRestoreColumnVec(P, j, &v));
118: PetscCheck(PetscAbsReal(norm - 2.0 * (j + 2)) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Column %" PetscInt_FMT " of the submatrix gives %g, expected %g", j, (double)norm, 2.0 * (j + 2));
119: }
120: PetscCall(MatDenseRestoreSubMatrix(C, &P));
122: /* An AIJ times dense product takes a different symbolic route; it only keeps the VecType when the Mat it
123: creates is the same kind of Mat as the block, which is not the case for a host AIJ and a device block */
124: PetscCall(MatCreate(PETSC_COMM_WORLD, &S2));
125: PetscCall(MatSetSizes(S2, rend - rstart, rend - rstart, M, M));
126: PetscCall(MatSetType(S2, MATAIJ));
127: PetscCall(MatSetUp(S2));
128: for (i = rstart; i < rend; i++) PetscCall(MatSetValue(S2, i, i, 3.0, INSERT_VALUES));
129: PetscCall(MatAssemblyBegin(S2, MAT_FINAL_ASSEMBLY));
130: PetscCall(MatAssemblyEnd(S2, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatMatMult(S2, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C2));
132: PetscCall(CheckProductVecType(C2, A, "The Mat the AIJ product created"));
134: /* and the two dense times dense products, which take a further symbolic route */
135: PetscCall(MatCreateDenseFromVecType(PETSC_COMM_WORLD, vtype, PETSC_DECIDE, PETSC_DECIDE, N, N, PETSC_DECIDE, NULL, &D));
136: PetscCall(MatZeroEntries(D));
137: PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C3));
138: PetscCall(CheckProductVecType(C3, A, "The Mat the dense product created"));
139: PetscCall(MatMatTransposeMult(A, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C4));
140: PetscCall(CheckProductVecType(C4, A, "The Mat the dense transpose product created"));
142: /* and a dense times AIJ product, which is not available for a sequential MATDENSECUDA or MATDENSEHIP */
143: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &same, MATSEQDENSE, MATMPIDENSE, ""));
144: if (same) {
145: PetscCall(MatMatMult(C4, S2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C5));
146: PetscCall(CheckProductVecType(C5, A, "The Mat the dense times AIJ product created"));
147: PetscCall(MatDestroy(&C5));
148: }
150: /* with dense operands of different VecType the Mat a product creates takes that of A, whatever the number of ranks */
151: PetscCall(CreateDenseDefaultVecType(A, N, N, &E));
152: PetscCall(CreateDenseDefaultVecType(A, M, N, &F));
153: PetscCall(MatMatMult(A, E, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C6));
154: PetscCall(CheckProductVecType(C6, A, "The Mat the dense product with operands of different VecType created"));
155: PetscCall(MatDestroy(&C6));
156: PetscCall(MatMatTransposeMult(A, F, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C6));
157: PetscCall(CheckProductVecType(C6, A, "The Mat the dense transpose product with operands of different VecType created"));
158: PetscCall(MatDestroy(&C6));
159: PetscCall(MatTransposeMatMult(A, F, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C6));
160: PetscCall(CheckProductVecType(C6, A, "The Mat the transpose dense product with operands of different VecType created"));
161: PetscCall(MatDestroy(&C6));
163: PetscCall(MatDestroy(&F));
164: PetscCall(MatDestroy(&E));
165: PetscCall(MatDestroy(&C4));
166: PetscCall(MatDestroy(&C3));
167: PetscCall(MatDestroy(&D));
168: PetscCall(MatDestroy(&C2));
169: PetscCall(MatDestroy(&S2));
170: PetscCall(MatDestroy(&C));
171: PetscCall(MatDestroy(&S));
172: PetscCall(MatDestroy(&A));
173: PetscCall(PetscFinalize());
174: return 0;
175: }
177: /*TEST
179: testset:
180: output_file: output/empty.out
181: nsize: {{1 2}}
182: test:
183: suffix: standard
184: args: -vec_type standard
185: test:
186: suffix: cuda
187: requires: cuda
188: args: -vec_type cuda
189: test:
190: suffix: hip
191: requires: hip
192: args: -vec_type hip
193: test:
194: suffix: kokkos
195: requires: kokkos_kernels !sycl
196: args: -vec_type kokkos
198: TEST*/