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
20: PetscInitialize(&argc, &argv, (char *)0, help);
21: #if defined(PETSC_HAVE_HYPRE)
22: PetscOptionsGetBool(NULL, NULL, "-test_hypre", &test_hypre, NULL);
23: #endif
24: MPI_Comm_size(PETSC_COMM_WORLD, &size);
25: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
27: MatCreate(PETSC_COMM_WORLD, &A);
28: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3);
29: MatSetType(A, MATAIJ);
30: MatSetFromOptions(A);
31: MatSetUp(A);
32: MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
34: if (rank == 0) MatSetValues(A, 3, ij, 3, ij, a, ADD_VALUES);
35: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
36: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
38: /* Test MatMatMult() */
39: MatTranspose(A, MAT_INITIAL_MATRIX, &B); /* B = A^T */
40: MatMatMult(B, A, MAT_INITIAL_MATRIX, fill, &C); /* C = B*A */
41: MatMatMult(B, A, MAT_REUSE_MATRIX, fill, &C); /* recompute C=B*A */
42: MatSetOptionsPrefix(C, "C_");
43: MatMatMultEqual(B, A, C, 10, &isequal);
46: MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D); /* D = C*A = (A^T*A)*A */
47: MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D);
48: MatMatMultEqual(C, A, D, 10, &isequal);
51: MatDestroy(&B);
52: MatDestroy(&C);
53: MatDestroy(&D);
55: /* Test MatPtAP */
56: MatDuplicate(A, MAT_COPY_VALUES, &B); /* B = A */
57: MatPtAP(A, B, MAT_INITIAL_MATRIX, fill, &C); /* C = B^T*A*B */
58: MatPtAPMultEqual(A, B, C, 10, &isequal);
61: /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */
62: MatPtAP(A, B, MAT_REUSE_MATRIX, fill, &C);
63: MatPtAPMultEqual(A, B, C, 10, &isequal);
66: MatDestroy(&C);
68: /* Test MatPtAP with A as a dense matrix */
69: {
70: Mat Adense;
71: MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense);
72: MatPtAP(Adense, B, MAT_INITIAL_MATRIX, fill, &C);
74: MatPtAPMultEqual(Adense, B, C, 10, &isequal);
76: MatDestroy(&Adense);
77: }
79: if (size == 1) {
80: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
81: testPTAPRectangular();
83: /* test MatMatTransposeMult(): A*B^T */
84: MatMatTransposeMult(A, A, MAT_INITIAL_MATRIX, fill, &D); /* D = A*A^T */
85: MatScale(A, 2.0);
86: MatMatTransposeMult(A, A, MAT_REUSE_MATRIX, fill, &D);
87: MatMatTransposeMultEqual(A, A, D, 10, &isequal);
89: }
91: MatDestroy(&A);
92: MatDestroy(&B);
93: MatDestroy(&C);
94: MatDestroy(&D);
95: 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: /* set up A */
107: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows, 1, NULL, &A);
108: for (i = 0; i < rows; i++) MatSetValue(A, i, i, 1.0, INSERT_VALUES);
109: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
110: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
112: /* set up P */
113: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols, 5, NULL, &P);
114: MatSetValue(P, 0, 0, 1.0, INSERT_VALUES);
115: MatSetValue(P, 0, 1, 2.0, INSERT_VALUES);
116: MatSetValue(P, 0, 2, 0.0, INSERT_VALUES);
118: MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);
120: MatSetValue(P, 1, 0, 0.0, INSERT_VALUES);
121: MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);
122: MatSetValue(P, 1, 2, 1.0, INSERT_VALUES);
124: MatSetValue(P, 2, 0, 3.0, INSERT_VALUES);
125: MatSetValue(P, 2, 1, 0.0, INSERT_VALUES);
126: MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);
128: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
131: /* compute C */
132: MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
134: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
137: /* compare results */
138: /*
139: printf("C:\n");
140: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
142: blitz::Array<double,2> actualC(cols, cols);
143: actualC = 0.0;
144: for (int i=0; i<cols; i++) {
145: for (int j=0; j<cols; j++) {
146: MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
147: }
148: }
149: blitz::Array<double,2> expectedC(cols, cols);
150: expectedC = 0.0;
152: expectedC(0,0) = 10.0;
153: expectedC(0,1) = 2.0;
154: expectedC(0,2) = -9.0;
155: expectedC(0,3) = -1.0;
156: expectedC(1,0) = 2.0;
157: expectedC(1,1) = 5.0;
158: expectedC(1,2) = -1.0;
159: expectedC(1,3) = -2.0;
160: expectedC(2,0) = -9.0;
161: expectedC(2,1) = -1.0;
162: expectedC(2,2) = 10.0;
163: expectedC(2,3) = 0.0;
164: expectedC(3,0) = -1.0;
165: expectedC(3,1) = -2.0;
166: expectedC(3,2) = 0.0;
167: expectedC(3,3) = 1.0;
169: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
170: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
171: */
173: MatDestroy(&A);
174: MatDestroy(&P);
175: MatDestroy(&C);
176: return 0;
177: }
179: /*TEST
181: test:
183: test:
184: suffix: 2
185: nsize: 2
186: args: -matmatmult_via nonscalable
187: output_file: output/ex93_1.out
189: test:
190: suffix: 3
191: nsize: 2
192: output_file: output/ex93_1.out
194: test:
195: suffix: 4
196: nsize: 2
197: args: -matptap_via scalable
198: output_file: output/ex93_1.out
200: test:
201: suffix: btheap
202: args: -matmatmult_via btheap -matmattransmult_via color
203: output_file: output/ex93_1.out
205: test:
206: suffix: heap
207: args: -matmatmult_via heap
208: output_file: output/ex93_1.out
210: #HYPRE PtAP is broken for complex numbers
211: test:
212: suffix: hypre
213: nsize: 3
214: requires: hypre !complex
215: args: -matmatmult_via hypre -matptap_via hypre -test_hypre
216: output_file: output/ex93_hypre.out
218: test:
219: suffix: llcondensed
220: args: -matmatmult_via llcondensed
221: output_file: output/ex93_1.out
223: test:
224: suffix: scalable
225: args: -matmatmult_via scalable
226: output_file: output/ex93_1.out
228: test:
229: suffix: scalable_fast
230: args: -matmatmult_via scalable_fast
231: output_file: output/ex93_1.out
233: TEST*/