Actual source code: ex93.c
petsc-3.8.4 2018-03-24
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: MatMatMultSymbolic(C,A,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: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
70: MatNorm(D,NORM_FROBENIUS,&norm);
71: if (norm > 1.e-15 && !rank) {
72: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);
73: }
74: MatDestroy(&C);
75: MatDestroy(&D);
77: /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
78: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
79: MatSetOptionsPrefix(C,"C=BtAB_");
80: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
81: PetscPrintf(PETSC_COMM_WORLD,"\n");
83: if (!test_hypre) {
84: MatPtAPSymbolic(A,B,fill,&D);
85: MatPtAPNumeric(A,B,D);
86: MatSetOptionsPrefix(D,"D=BtAB_");
87: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
88: PetscPrintf(PETSC_COMM_WORLD,"\n");
90: /* Repeat numeric product to test reuse of the previous symbolic product */
91: MatPtAPNumeric(A,B,D);
92: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
93: MatNorm(D,NORM_FROBENIUS,&norm);
94: if (norm > 1.e-15) {
95: PetscPrintf(PETSC_COMM_WORLD,"Error in symbolic/numeric MatPtAP: %g\n",norm);
96: }
97: MatDestroy(&B);
98: MatDestroy(&C);
99: MatDestroy(&D);
100: }
102: if (size == 1) {
103: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
104: testPTAPRectangular();
106: /* test MatMatTransposeMult(): A*B^T */
107: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
108: MatSetOptionsPrefix(D,"D=A*A^T_");
109: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
110: PetscPrintf(PETSC_COMM_WORLD,"\n");
112: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
113: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C=A*B */
114: MatSetOptionsPrefix(C,"D=A*B=A*A^T_");
115: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
116: PetscPrintf(PETSC_COMM_WORLD,"\n");
117: }
119: MatDestroy(&A);
120: MatDestroy(&B);
121: MatDestroy(&C);
122: MatDestroy(&D);
123: PetscFinalize();
124: return ierr;
125: }
127: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
128: #define PETSc_CHKERRQ CHKERRQ
129: PetscErrorCode testPTAPRectangular(void)
130: {
132: const int rows = 3;
133: const int cols = 5;
134: PetscErrorCode _ierr;
135: int i;
136: Mat A;
137: Mat P;
138: Mat C;
141: /* set up A */
142: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
143: PETSc_CHKERRQ(_ierr);
144: for (i=0; i<rows; i++) {
145: _MatSetValue(A, i, i, 1.0, INSERT_VALUES);
146: PETSc_CHKERRQ(_ierr);
147: }
148: _MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
149: PETSc_CHKERRQ(_ierr);
150: _MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
151: PETSc_CHKERRQ(_ierr);
153: /* set up P */
154: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
155: PETSc_CHKERRQ(_ierr);
156: _MatSetValue(P, 0, 0, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
157: _MatSetValue(P, 0, 1, 2.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
158: _MatSetValue(P, 0, 2, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
160: _MatSetValue(P, 0, 3, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
162: _MatSetValue(P, 1, 0, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
163: _MatSetValue(P, 1, 1, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
164: _MatSetValue(P, 1, 2, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
166: _MatSetValue(P, 2, 0, 3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
167: _MatSetValue(P, 2, 1, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
168: _MatSetValue(P, 2, 2, -3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
170: _MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
171: PETSc_CHKERRQ(_ierr);
172: _MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
173: PETSc_CHKERRQ(_ierr);
175: /* compute C */
176: _MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
177: PETSc_CHKERRQ(_ierr);
179: _MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
180: PETSc_CHKERRQ(_ierr);
181: _MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
182: PETSc_CHKERRQ(_ierr);
184: /* compare results */
185: /*
186: printf("C:\n");
187: _MatView(C,PETSC_VIEWER_STDOUT_WORLD);PETSc_CHKERRQ(_ierr);
189: blitz::Array<double,2> actualC(cols, cols);
190: actualC = 0.0;
191: for (int i=0; i<cols; i++) {
192: for (int j=0; j<cols; j++) {
193: _MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
194: PETSc_CHKERRQ(_ierr); ;
195: }
196: }
197: blitz::Array<double,2> expectedC(cols, cols);
198: expectedC = 0.0;
200: expectedC(0,0) = 10.0;
201: expectedC(0,1) = 2.0;
202: expectedC(0,2) = -9.0;
203: expectedC(0,3) = -1.0;
204: expectedC(1,0) = 2.0;
205: expectedC(1,1) = 5.0;
206: expectedC(1,2) = -1.0;
207: expectedC(1,3) = -2.0;
208: expectedC(2,0) = -9.0;
209: expectedC(2,1) = -1.0;
210: expectedC(2,2) = 10.0;
211: expectedC(2,3) = 0.0;
212: expectedC(3,0) = -1.0;
213: expectedC(3,1) = -2.0;
214: expectedC(3,2) = 0.0;
215: expectedC(3,3) = 1.0;
217: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
218: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
219: */
221: _MatDestroy(&A);
222: PETSc_CHKERRQ(_ierr);
223: _MatDestroy(&P);
224: PETSc_CHKERRQ(_ierr);
225: _MatDestroy(&C);
226: PETSc_CHKERRQ(_ierr);
227: return(0);
228: }