Actual source code: ex93.c
petsc-3.10.5 2019-03-28
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: PetscScalar none=-1.;
14: PetscReal fill=4;
15: PetscReal norm;
16: PetscMPIInt size,rank;
17: PetscBool test_hypre=PETSC_FALSE;
19: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
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) {
34: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
35: }
36: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
37: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
38: MatSetOptionsPrefix(A,"A_");
39: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
40: PetscPrintf(PETSC_COMM_WORLD,"\n");
42: /* Test MatMatMult() */
43: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
44: MatSetOptionsPrefix(B,"B_");
45: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A */
46: MatMatMultNumeric(B,A,C); /* recompute C=B*A */
47: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C); /* recompute C=B*A */
48: MatSetOptionsPrefix(C,"C_");
49: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
50: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
52: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D);
53: MatMatMultNumeric(C,A,D); /* D = C*A = (A^T*A)*A */
54: MatSetOptionsPrefix(D,"D_");
55: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
56: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
58: /* Repeat the numeric product to test reuse of the previous symbolic product */
59: MatMatMultNumeric(C,A,D);
60: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
61: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
63: MatDestroy(&B);
64: MatDestroy(&C);
66: /* Test PtAP routine. */
67: MatDuplicate(A,MAT_COPY_VALUES,&B); /* B = A */
68: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
69: if (!test_hypre) {
70: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
71: MatNorm(D,NORM_FROBENIUS,&norm);
72: if (norm > 1.e-15 && !rank) {
73: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);
74: }
75: }
76: MatDestroy(&C);
77: MatDestroy(&D);
79: /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
80: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
81: MatSetOptionsPrefix(C,"C=BtAB_");
82: if (!test_hypre) {
83: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
84: PetscPrintf(PETSC_COMM_WORLD,"\n");
85: }
87: if (!test_hypre) {
88: MatPtAPSymbolic(A,B,fill,&D);
89: MatPtAPNumeric(A,B,D);
90: MatSetOptionsPrefix(D,"D=BtAB_");
91: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
92: PetscPrintf(PETSC_COMM_WORLD,"\n");
94: /* Repeat numeric product to test reuse of the previous symbolic product */
95: MatPtAPNumeric(A,B,D);
96: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
97: MatNorm(D,NORM_FROBENIUS,&norm);
98: if (norm > 1.e-15) {
99: PetscPrintf(PETSC_COMM_WORLD,"Error in symbolic/numeric MatPtAP: %g\n",norm);
100: }
101: MatDestroy(&D);
102: }
103: MatDestroy(&C);
104: MatDestroy(&B);
106: if (size == 1) {
107: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
108: testPTAPRectangular();
110: /* test MatMatTransposeMult(): A*B^T */
111: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
112: MatSetOptionsPrefix(D,"D=A*A^T_");
113: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
114: PetscPrintf(PETSC_COMM_WORLD,"\n");
116: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
117: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C=A*B */
118: MatSetOptionsPrefix(C,"D=A*B=A*A^T_");
119: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
120: PetscPrintf(PETSC_COMM_WORLD,"\n");
121: }
123: MatDestroy(&A);
124: MatDestroy(&B);
125: MatDestroy(&C);
126: MatDestroy(&D);
127: PetscFinalize();
128: return ierr;
129: }
131: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
132: PetscErrorCode testPTAPRectangular(void)
133: {
135: const int rows = 3;
136: const int cols = 5;
138: int i;
139: Mat A;
140: Mat P;
141: Mat C;
144: /* set up A */
145: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
146: for (i=0; i<rows; i++) {
147: MatSetValue(A, i, i, 1.0, INSERT_VALUES);
148: }
149: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
150: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
152: /* set up P */
153: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
154: MatSetValue(P, 0, 0, 1.0, INSERT_VALUES);
155: MatSetValue(P, 0, 1, 2.0, INSERT_VALUES);
156: MatSetValue(P, 0, 2, 0.0, INSERT_VALUES);
158: MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);
160: MatSetValue(P, 1, 0, 0.0, INSERT_VALUES);
161: MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);
162: MatSetValue(P, 1, 2, 1.0, INSERT_VALUES);
164: MatSetValue(P, 2, 0, 3.0, INSERT_VALUES);
165: MatSetValue(P, 2, 1, 0.0, INSERT_VALUES);
166: MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);
168: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
169: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
171: /* compute C */
172: MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
174: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
175: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
177: /* compare results */
178: /*
179: printf("C:\n");
180: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
182: blitz::Array<double,2> actualC(cols, cols);
183: actualC = 0.0;
184: for (int i=0; i<cols; i++) {
185: for (int j=0; j<cols; j++) {
186: MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
187: ;
188: }
189: }
190: blitz::Array<double,2> expectedC(cols, cols);
191: expectedC = 0.0;
193: expectedC(0,0) = 10.0;
194: expectedC(0,1) = 2.0;
195: expectedC(0,2) = -9.0;
196: expectedC(0,3) = -1.0;
197: expectedC(1,0) = 2.0;
198: expectedC(1,1) = 5.0;
199: expectedC(1,2) = -1.0;
200: expectedC(1,3) = -2.0;
201: expectedC(2,0) = -9.0;
202: expectedC(2,1) = -1.0;
203: expectedC(2,2) = 10.0;
204: expectedC(2,3) = 0.0;
205: expectedC(3,0) = -1.0;
206: expectedC(3,1) = -2.0;
207: expectedC(3,2) = 0.0;
208: expectedC(3,3) = 1.0;
210: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
211: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
212: */
214: MatDestroy(&A);
215: MatDestroy(&P);
216: MatDestroy(&C);
217: return(0);
218: }
220: /*TEST
222: test:
224: test:
225: suffix: 2
226: nsize: 2
227: args: -B_matmatmult_via nonscalable
229: test:
230: suffix: 3
231: nsize: 2
232: output_file: output/ex93_2.out
234: test:
235: suffix: 4
236: nsize: 2
237: args: -A_matptap_via scalable
238: output_file: output/ex93_2.out
240: test:
241: suffix: btheap
242: args: -B_matmatmult_via btheap
243: output_file: output/ex93_1.out
245: test:
246: suffix: heap
247: args: -B_matmatmult_via heap
248: output_file: output/ex93_1.out
250: test:
251: suffix: hypre
252: nsize: 3
253: requires: hypre
254: args: -B_matmatmult_via hypre -A_matptap_via hypre -test_hypre
256: test:
257: suffix: llcondensed
258: args: -B_matmatmult_via llcondensed
259: output_file: output/ex93_1.out
261: test:
262: suffix: scalable
263: args: -B_matmatmult_via scalable
264: output_file: output/ex93_1.out
266: test:
267: suffix: scalable_fast
268: args: -B_matmatmult_via scalable_fast
269: output_file: output/ex93_1.out
271: TEST*/