Actual source code: ex93.c
1: static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n";
3: #include <petscmat.h>
5: extern PetscErrorCode testPTAPRectangular(void);
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
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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));
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));
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));
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");
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");
51: PetscCall(MatDestroy(&B));
52: PetscCall(MatDestroy(&C));
53: PetscCall(MatDestroy(&D));
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");
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");
66: PetscCall(MatDestroy(&C));
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));
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: }
79: if (size == 1) {
80: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
81: PetscCall(testPTAPRectangular());
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: }
91: PetscCall(MatDestroy(&A));
92: PetscCall(MatDestroy(&B));
93: PetscCall(MatDestroy(&C));
94: PetscCall(MatDestroy(&D));
95: PetscCall(PetscFinalize());
96: return 0;
97: }
99: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
100: PetscErrorCode testPTAPRectangular(void)
101: {
102: const int rows = 3, cols = 5;
103: int i;
104: Mat A, P, C;
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));
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));
119: PetscCall(MatSetValue(P, 0, 3, -1.0, INSERT_VALUES));
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));
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));
129: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
130: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
132: /* compute C */
133: PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C));
135: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
136: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
138: /* compare results */
139: /*
140: printf("C:\n");
141: PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD));
143: blitz::Array<double,2> actualC(cols, cols);
144: actualC = 0.0;
145: for (int i=0; i<cols; i++) {
146: for (int j=0; j<cols; j++) {
147: PetscCall(MatGetValues(C, 1, &i, 1, &j, &actualC(i,j)));
148: }
149: }
150: blitz::Array<double,2> expectedC(cols, cols);
151: expectedC = 0.0;
153: expectedC(0,0) = 10.0;
154: expectedC(0,1) = 2.0;
155: expectedC(0,2) = -9.0;
156: expectedC(0,3) = -1.0;
157: expectedC(1,0) = 2.0;
158: expectedC(1,1) = 5.0;
159: expectedC(1,2) = -1.0;
160: expectedC(1,3) = -2.0;
161: expectedC(2,0) = -9.0;
162: expectedC(2,1) = -1.0;
163: expectedC(2,2) = 10.0;
164: expectedC(2,3) = 0.0;
165: expectedC(3,0) = -1.0;
166: expectedC(3,1) = -2.0;
167: expectedC(3,2) = 0.0;
168: expectedC(3,3) = 1.0;
170: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
171: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
172: */
174: PetscCall(MatDestroy(&A));
175: PetscCall(MatDestroy(&P));
176: PetscCall(MatDestroy(&C));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*TEST
182: test:
184: test:
185: suffix: 2
186: nsize: 2
187: args: -matmatmult_via nonscalable
188: output_file: output/ex93_1.out
190: test:
191: suffix: 3
192: nsize: 2
193: output_file: output/ex93_1.out
195: test:
196: suffix: 4
197: nsize: 2
198: args: -matptap_via scalable
199: output_file: output/ex93_1.out
201: test:
202: suffix: btheap
203: args: -matmatmult_via btheap -matmattransmult_via color
204: output_file: output/ex93_1.out
206: test:
207: suffix: heap
208: args: -matmatmult_via heap
209: output_file: output/ex93_1.out
211: #HYPRE PtAP is broken for complex numbers
212: test:
213: suffix: hypre
214: nsize: 3
215: requires: hypre !complex
216: args: -matmatmult_via hypre -matptap_via hypre -test_hypre
217: output_file: output/ex93_hypre.out
219: test:
220: suffix: llcondensed
221: args: -matmatmult_via llcondensed
222: output_file: output/ex93_1.out
224: test:
225: suffix: scalable
226: args: -matmatmult_via scalable
227: output_file: output/ex93_1.out
229: test:
230: suffix: scalable_fast
231: args: -matmatmult_via scalable_fast
232: output_file: output/ex93_1.out
234: TEST*/