1 static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** argv)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;
14
15 PetscFunctionBeginUser;
16 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
20
21 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
22 PetscCall(PetscRandomSetFromOptions(r));
23
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));
32
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));
60
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));
71
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));
82
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));
86
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));
98
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));
106
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));
113
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));
120
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 }
130
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));
137
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 }
155
156 PetscCall(MatDestroy(&C));
157 PetscCall(MatDestroy(&B));
158 PetscCall(MatDestroy(&A));
159 PetscCall(PetscFinalize());
160 return 0;
161 }
162
163 /*TEST
164
165 test:
166 args: -M 10 -N 10
167 output_file: output/empty.out
168
169 test:
170 suffix: 2
171 nsize: 3
172 output_file: output/empty.out
173
174 test:
175 suffix: 3
176 nsize: 2
177 args: -matmattransmult_mpidense_mpidense_via cyclic
178 output_file: output/empty.out
179
180 test:
181 suffix: 4
182 args: -test_userAPI
183 output_file: output/empty.out
184
185 test:
186 suffix: 5
187 nsize: 3
188 args: -test_userAPI
189 output_file: output/empty.out
190
191 TEST*/
192