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: PetscInitialize(&argc,&argv,(char*)0,help);
20: #if defined(PETSC_HAVE_HYPRE)
21: PetscOptionsGetBool(NULL,NULL,"-test_hypre",&test_hypre,NULL);
22: #endif
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: MatCreate(PETSC_COMM_WORLD,&A);
27: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
28: MatSetType(A,MATAIJ);
29: MatSetFromOptions(A);
30: MatSetUp(A);
31: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
33: if (rank == 0) {
34: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
35: }
36: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
37: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
39: /* Test MatMatMult() */
40: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
41: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A */
42: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C); /* recompute C=B*A */
43: MatSetOptionsPrefix(C,"C_");
44: MatMatMultEqual(B,A,C,10,&isequal);
47: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D); /* D = C*A = (A^T*A)*A */
48: MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D);
49: MatMatMultEqual(C,A,D,10,&isequal);
52: MatDestroy(&B);
53: MatDestroy(&C);
54: MatDestroy(&D);
56: /* Test MatPtAP */
57: MatDuplicate(A,MAT_COPY_VALUES,&B); /* B = A */
58: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
59: MatPtAPMultEqual(A,B,C,10,&isequal);
62: /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */
63: MatPtAP(A,B,MAT_REUSE_MATRIX,fill,&C);
64: MatPtAPMultEqual(A,B,C,10,&isequal);
67: MatDestroy(&C);
69: /* Test MatPtAP with A as a dense matrix */
70: {
71: Mat Adense;
72: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
73: MatPtAP(Adense,B,MAT_INITIAL_MATRIX,fill,&C);
75: MatPtAPMultEqual(Adense,B,C,10,&isequal);
77: MatDestroy(&Adense);
78: }
80: if (size == 1) {
81: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
82: testPTAPRectangular();
84: /* test MatMatTransposeMult(): A*B^T */
85: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
86: MatScale(A,2.0);
87: MatMatTransposeMult(A,A,MAT_REUSE_MATRIX,fill,&D);
88: MatMatTransposeMultEqual(A,A,D,10,&isequal);
90: }
92: MatDestroy(&A);
93: MatDestroy(&B);
94: MatDestroy(&C);
95: MatDestroy(&D);
96: PetscFinalize();
97: return 0;
98: }
100: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
101: PetscErrorCode testPTAPRectangular(void)
102: {
103: const int rows = 3,cols = 5;
104: int i;
105: Mat A,P,C;
107: /* set up A */
108: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
109: for (i=0; i<rows; i++) {
110: MatSetValue(A, i, i, 1.0, INSERT_VALUES);
111: }
112: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
115: /* set up P */
116: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
117: MatSetValue(P, 0, 0, 1.0, INSERT_VALUES);
118: MatSetValue(P, 0, 1, 2.0, INSERT_VALUES);
119: MatSetValue(P, 0, 2, 0.0, INSERT_VALUES);
121: MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);
123: MatSetValue(P, 1, 0, 0.0, INSERT_VALUES);
124: MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);
125: MatSetValue(P, 1, 2, 1.0, INSERT_VALUES);
127: MatSetValue(P, 2, 0, 3.0, INSERT_VALUES);
128: MatSetValue(P, 2, 1, 0.0, INSERT_VALUES);
129: MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);
131: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
134: /* compute C */
135: MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
137: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
140: /* compare results */
141: /*
142: printf("C:\n");
143: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
145: blitz::Array<double,2> actualC(cols, cols);
146: actualC = 0.0;
147: for (int i=0; i<cols; i++) {
148: for (int j=0; j<cols; j++) {
149: MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
150: }
151: }
152: blitz::Array<double,2> expectedC(cols, cols);
153: expectedC = 0.0;
155: expectedC(0,0) = 10.0;
156: expectedC(0,1) = 2.0;
157: expectedC(0,2) = -9.0;
158: expectedC(0,3) = -1.0;
159: expectedC(1,0) = 2.0;
160: expectedC(1,1) = 5.0;
161: expectedC(1,2) = -1.0;
162: expectedC(1,3) = -2.0;
163: expectedC(2,0) = -9.0;
164: expectedC(2,1) = -1.0;
165: expectedC(2,2) = 10.0;
166: expectedC(2,3) = 0.0;
167: expectedC(3,0) = -1.0;
168: expectedC(3,1) = -2.0;
169: expectedC(3,2) = 0.0;
170: expectedC(3,3) = 1.0;
172: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
173: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
174: */
176: MatDestroy(&A);
177: MatDestroy(&P);
178: MatDestroy(&C);
179: return 0;
180: }
182: /*TEST
184: test:
186: test:
187: suffix: 2
188: nsize: 2
189: args: -matmatmult_via nonscalable
190: output_file: output/ex93_1.out
192: test:
193: suffix: 3
194: nsize: 2
195: output_file: output/ex93_1.out
197: test:
198: suffix: 4
199: nsize: 2
200: args: -matptap_via scalable
201: output_file: output/ex93_1.out
203: test:
204: suffix: btheap
205: args: -matmatmult_via btheap -matmattransmult_via color
206: output_file: output/ex93_1.out
208: test:
209: suffix: heap
210: args: -matmatmult_via heap
211: output_file: output/ex93_1.out
213: #HYPRE PtAP is broken for complex numbers
214: test:
215: suffix: hypre
216: nsize: 3
217: requires: hypre !complex
218: args: -matmatmult_via hypre -matptap_via hypre -test_hypre
219: output_file: output/ex93_hypre.out
221: test:
222: suffix: llcondensed
223: args: -matmatmult_via llcondensed
224: output_file: output/ex93_1.out
226: test:
227: suffix: scalable
228: args: -matmatmult_via scalable
229: output_file: output/ex93_1.out
231: test:
232: suffix: scalable_fast
233: args: -matmatmult_via scalable_fast
234: output_file: output/ex93_1.out
236: TEST*/