Actual source code: ex109.c
1: static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **argv)
6: {
7: Mat A, B, C, D, AT;
8: PetscInt i, M, N, Istart, Iend, n = 7, j, J, Ii, m = 8, am, an;
9: PetscScalar v;
10: PetscRandom r;
11: PetscBool equal = PETSC_FALSE, flg;
12: PetscReal fill = 1.0, norm;
13: PetscMPIInt size;
15: PetscFunctionBeginUser;
16: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
17: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
19: PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL));
21: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
22: PetscCall(PetscRandomSetFromOptions(r));
24: /* Create a aij matrix A */
25: M = N = m * n;
26: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
27: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
28: PetscCall(MatSetType(A, MATAIJ));
29: PetscCall(MatSetFromOptions(A));
30: PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
31: PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
33: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
34: am = Iend - Istart;
35: for (Ii = Istart; Ii < Iend; Ii++) {
36: v = -1.0;
37: i = Ii / n;
38: j = Ii - i * n;
39: if (i > 0) {
40: J = Ii - n;
41: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
42: }
43: if (i < m - 1) {
44: J = Ii + n;
45: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
46: }
47: if (j > 0) {
48: J = Ii - 1;
49: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
50: }
51: if (j < n - 1) {
52: J = Ii + 1;
53: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
54: }
55: v = 4.0;
56: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
57: }
58: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
59: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
61: /* Create a dense matrix B */
62: PetscCall(MatGetLocalSize(A, &am, &an));
63: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
64: PetscCall(MatSetSizes(B, an, PETSC_DECIDE, PETSC_DECIDE, M));
65: PetscCall(MatSetType(B, MATDENSE));
66: PetscCall(MatSeqDenseSetPreallocation(B, NULL));
67: PetscCall(MatMPIDenseSetPreallocation(B, NULL));
68: PetscCall(MatSetFromOptions(B));
69: PetscCall(MatSetRandom(B, r));
70: PetscCall(PetscRandomDestroy(&r));
72: /* Test reuse of user-provided dense C (unassembled) -- not recommended usage */
73: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
74: PetscCall(MatSetType(C, MATDENSE));
75: PetscCall(MatSetSizes(C, am, PETSC_DECIDE, PETSC_DECIDE, N));
76: PetscCall(MatSetFromOptions(C));
77: PetscCall(MatSetUp(C));
78: PetscCall(MatZeroEntries(C));
79: PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C));
80: PetscCall(MatNorm(C, NORM_INFINITY, &norm));
81: PetscCall(MatDestroy(&C));
83: /* Test C = A*B (aij*dense) */
84: PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C));
85: PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C));
87: /* Test developer API */
88: PetscCall(MatProductCreate(A, B, NULL, &D));
89: PetscCall(MatProductSetType(D, MATPRODUCT_AB));
90: PetscCall(MatProductSetAlgorithm(D, "default"));
91: PetscCall(MatProductSetFill(D, fill));
92: PetscCall(MatProductSetFromOptions(D));
93: PetscCall(MatProductSymbolic(D));
94: for (i = 0; i < 2; i++) PetscCall(MatProductNumeric(D));
95: PetscCall(MatEqual(C, D, &equal));
96: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "C != D");
97: PetscCall(MatDestroy(&D));
99: /* Test D = AT*B (transpose(aij)*dense) */
100: PetscCall(MatCreateTranspose(A, &AT));
101: PetscCall(MatMatMult(AT, B, MAT_INITIAL_MATRIX, fill, &D));
102: PetscCall(MatMatMultEqual(AT, B, D, 10, &equal));
103: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != AT*B (transpose(aij)*dense)");
104: PetscCall(MatDestroy(&D));
105: PetscCall(MatDestroy(&AT));
107: /* Test D = C*A (dense*aij) */
108: PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D));
109: PetscCall(MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D));
110: PetscCall(MatMatMultEqual(C, A, D, 10, &equal));
111: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != C*A (dense*aij)");
112: PetscCall(MatDestroy(&D));
114: /* Test D = A*C (aij*dense) */
115: PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, fill, &D));
116: PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &D));
117: PetscCall(MatMatMultEqual(A, C, D, 10, &equal));
118: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != A*C (aij*dense)");
119: PetscCall(MatDestroy(&D));
121: /* Test D = B*C (dense*dense) */
122: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
123: if (size == 1) {
124: PetscCall(MatMatMult(B, C, MAT_INITIAL_MATRIX, fill, &D));
125: PetscCall(MatMatMult(B, C, MAT_REUSE_MATRIX, fill, &D));
126: PetscCall(MatMatMultEqual(B, C, D, 10, &equal));
127: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != B*C (dense*dense)");
128: PetscCall(MatDestroy(&D));
129: }
131: /* Test D = B*C^T (dense*dense) */
132: PetscCall(MatMatTransposeMult(B, C, MAT_INITIAL_MATRIX, fill, &D));
133: PetscCall(MatMatTransposeMult(B, C, MAT_REUSE_MATRIX, fill, &D));
134: PetscCall(MatMatTransposeMultEqual(B, C, D, 10, &equal));
135: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != B*C^T (dense*dense)");
136: PetscCall(MatDestroy(&D));
138: /* Test MatProductCreateWithMat() and reuse C and B for B = A*C */
139: flg = PETSC_FALSE;
140: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_userAPI", &flg));
141: if (flg) {
142: /* user driver */
143: PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &B));
144: } else {
145: /* clear internal data structures related with previous products to avoid circular references */
146: PetscCall(MatProductClear(A));
147: PetscCall(MatProductClear(B));
148: PetscCall(MatProductClear(C));
149: PetscCall(MatProductCreateWithMat(A, C, NULL, B));
150: PetscCall(MatProductSetType(B, MATPRODUCT_AB));
151: PetscCall(MatProductSetFromOptions(B));
152: PetscCall(MatProductSymbolic(B));
153: PetscCall(MatProductNumeric(B));
154: }
156: PetscCall(MatDestroy(&C));
157: PetscCall(MatDestroy(&B));
158: PetscCall(MatDestroy(&A));
159: PetscCall(PetscFinalize());
160: return 0;
161: }
163: /*TEST
165: test:
166: args: -M 10 -N 10
167: output_file: output/ex109.out
169: test:
170: suffix: 2
171: nsize: 3
172: output_file: output/ex109.out
174: test:
175: suffix: 3
176: nsize: 2
177: args: -matmattransmult_mpidense_mpidense_via cyclic
178: output_file: output/ex109.out
180: test:
181: suffix: 4
182: args: -test_userAPI
183: output_file: output/ex109.out
185: test:
186: suffix: 5
187: nsize: 3
188: args: -test_userAPI
189: output_file: output/ex109.out
191: TEST*/