1 static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n";
2
3 #include <petscmat.h>
4
5 extern PetscErrorCode testPTAPRectangular(void);
6
main(int argc,char ** argv)7 int main(int argc, char **argv)
8 {
9 Mat A, B, C, D;
10 PetscScalar a[] = {1., 1., 0., 0., 1., 1., 0., 0., 1.};
11 PetscInt ij[] = {0, 1, 2};
12 PetscReal fill = 4.0;
13 PetscMPIInt size, rank;
14 PetscBool isequal;
15 #if defined(PETSC_HAVE_HYPRE)
16 PetscBool test_hypre = PETSC_FALSE;
17 #endif
18
19 PetscFunctionBeginUser;
20 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
21 #if defined(PETSC_HAVE_HYPRE)
22 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_hypre", &test_hypre, NULL));
23 #endif
24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
25 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
26
27 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
28 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3));
29 PetscCall(MatSetType(A, MATAIJ));
30 PetscCall(MatSetFromOptions(A));
31 PetscCall(MatSetUp(A));
32 PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
33
34 if (rank == 0) PetscCall(MatSetValues(A, 3, ij, 3, ij, a, ADD_VALUES));
35 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
36 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
37
38 /* Test MatMatMult() */
39 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B)); /* B = A^T */
40 PetscCall(MatMatMult(B, A, MAT_INITIAL_MATRIX, fill, &C)); /* C = B*A */
41 PetscCall(MatMatMult(B, A, MAT_REUSE_MATRIX, fill, &C)); /* recompute C=B*A */
42 PetscCall(MatSetOptionsPrefix(C, "C_"));
43 PetscCall(MatMatMultEqual(B, A, C, 10, &isequal));
44 PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatMult: C != B*A");
45
46 PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D)); /* D = C*A = (A^T*A)*A */
47 PetscCall(MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D));
48 PetscCall(MatMatMultEqual(C, A, D, 10, &isequal));
49 PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatMult: D != C*A");
50
51 PetscCall(MatDestroy(&B));
52 PetscCall(MatDestroy(&C));
53 PetscCall(MatDestroy(&D));
54
55 /* Test MatPtAP */
56 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); /* B = A */
57 PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, fill, &C)); /* C = B^T*A*B */
58 PetscCall(MatPtAPMultEqual(A, B, C, 10, &isequal));
59 PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatPtAP: C != B^T*A*B");
60
61 /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */
62 PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, fill, &C));
63 PetscCall(MatPtAPMultEqual(A, B, C, 10, &isequal));
64 PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatPtAP(reuse): C != B^T*A*B");
65
66 PetscCall(MatDestroy(&C));
67
68 /* Test MatPtAP with A as a dense matrix */
69 {
70 Mat Adense;
71 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense));
72 PetscCall(MatPtAP(Adense, B, MAT_INITIAL_MATRIX, fill, &C));
73
74 PetscCall(MatPtAPMultEqual(Adense, B, C, 10, &isequal));
75 PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatPtAP(reuse): C != B^T*Adense*B");
76 PetscCall(MatDestroy(&Adense));
77 }
78
79 if (size == 1) {
80 /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
81 PetscCall(testPTAPRectangular());
82
83 /* test MatMatTransposeMult(): A*B^T */
84 PetscCall(MatMatTransposeMult(A, A, MAT_INITIAL_MATRIX, fill, &D)); /* D = A*A^T */
85 PetscCall(MatScale(A, 2.0));
86 PetscCall(MatMatTransposeMult(A, A, MAT_REUSE_MATRIX, fill, &D));
87 PetscCall(MatMatTransposeMultEqual(A, A, D, 10, &isequal));
88 PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatTranspose: D != A*A^T");
89 }
90
91 PetscCall(MatDestroy(&A));
92 PetscCall(MatDestroy(&B));
93 PetscCall(MatDestroy(&C));
94 PetscCall(MatDestroy(&D));
95 PetscCall(PetscFinalize());
96 return 0;
97 }
98
99 /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
testPTAPRectangular(void)100 PetscErrorCode testPTAPRectangular(void)
101 {
102 const int rows = 3, cols = 5;
103 int i;
104 Mat A, P, C;
105
106 PetscFunctionBegin;
107 /* set up A */
108 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows, 1, NULL, &A));
109 for (i = 0; i < rows; i++) PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES));
110 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
111 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
112
113 /* set up P */
114 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols, 5, NULL, &P));
115 PetscCall(MatSetValue(P, 0, 0, 1.0, INSERT_VALUES));
116 PetscCall(MatSetValue(P, 0, 1, 2.0, INSERT_VALUES));
117 PetscCall(MatSetValue(P, 0, 2, 0.0, INSERT_VALUES));
118
119 PetscCall(MatSetValue(P, 0, 3, -1.0, INSERT_VALUES));
120
121 PetscCall(MatSetValue(P, 1, 0, 0.0, INSERT_VALUES));
122 PetscCall(MatSetValue(P, 1, 1, -1.0, INSERT_VALUES));
123 PetscCall(MatSetValue(P, 1, 2, 1.0, INSERT_VALUES));
124
125 PetscCall(MatSetValue(P, 2, 0, 3.0, INSERT_VALUES));
126 PetscCall(MatSetValue(P, 2, 1, 0.0, INSERT_VALUES));
127 PetscCall(MatSetValue(P, 2, 2, -3.0, INSERT_VALUES));
128
129 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
130 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
131
132 /* compute C */
133 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C));
134
135 /* compare results */
136 /*
137 printf("C:\n");
138 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD));
139
140 blitz::Array<double,2> actualC(cols, cols);
141 actualC = 0.0;
142 for (int i=0; i<cols; i++) {
143 for (int j=0; j<cols; j++) PetscCall(MatGetValues(C, 1, &i, 1, &j, &actualC(i,j)));
144 }
145 blitz::Array<double,2> expectedC(cols, cols);
146 expectedC = 0.0;
147
148 expectedC(0,0) = 10.0;
149 expectedC(0,1) = 2.0;
150 expectedC(0,2) = -9.0;
151 expectedC(0,3) = -1.0;
152 expectedC(1,0) = 2.0;
153 expectedC(1,1) = 5.0;
154 expectedC(1,2) = -1.0;
155 expectedC(1,3) = -2.0;
156 expectedC(2,0) = -9.0;
157 expectedC(2,1) = -1.0;
158 expectedC(2,2) = 10.0;
159 expectedC(2,3) = 0.0;
160 expectedC(3,0) = -1.0;
161 expectedC(3,1) = -2.0;
162 expectedC(3,2) = 0.0;
163 expectedC(3,3) = 1.0;
164
165 int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
166 validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
167 */
168
169 PetscCall(MatDestroy(&A));
170 PetscCall(MatDestroy(&P));
171 PetscCall(MatDestroy(&C));
172 PetscFunctionReturn(PETSC_SUCCESS);
173 }
174
175 /*TEST
176
177 test:
178 output_file: output/empty.out
179
180 test:
181 suffix: 2
182 nsize: 2
183 args: -matmatmult_via nonscalable
184 output_file: output/empty.out
185
186 test:
187 suffix: 3
188 nsize: 2
189 output_file: output/empty.out
190
191 test:
192 suffix: 4
193 nsize: 2
194 args: -matptap_via scalable
195 output_file: output/empty.out
196
197 test:
198 suffix: btheap
199 args: -matmatmult_via btheap -matmattransmult_via color
200 output_file: output/empty.out
201
202 test:
203 suffix: heap
204 args: -matmatmult_via heap
205 output_file: output/empty.out
206
207 #HYPRE PtAP is broken for complex numbers
208 test:
209 suffix: hypre
210 nsize: 3
211 requires: hypre !complex
212 args: -matmatmult_via hypre -matptap_via hypre -test_hypre
213 output_file: output/empty.out
214
215 test:
216 suffix: llcondensed
217 args: -matmatmult_via llcondensed
218 output_file: output/empty.out
219
220 test:
221 suffix: scalable
222 args: -matmatmult_via scalable
223 output_file: output/empty.out
224
225 test:
226 suffix: scalable_fast
227 args: -matmatmult_via scalable_fast
228 output_file: output/empty.out
229
230 TEST*/
231