Actual source code: ex93.c
petsc-3.6.4 2016-04-12
1: static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n";
3: #include <petscmat.h>
5: extern PetscErrorCode testPTAPRectangular(void);
9: int main(int argc,char **argv)
10: {
11: Mat A,B,C,D;
12: PetscScalar a[] ={1.,1.,0.,0.,1.,1.,0.,0.,1.};
13: PetscInt ij[]={0,1,2};
14: PetscScalar none=-1.;
16: PetscReal fill=4;
17: PetscReal norm;
18: PetscMPIInt size,rank;
20: PetscInitialize(&argc,&argv,(char*)0,help);
21: MPI_Comm_size(PETSC_COMM_WORLD,&size);
22: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: MatCreate(PETSC_COMM_WORLD,&A);
25: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
26: MatSetType(A,MATAIJ);
27: MatSetUp(A);
28: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
29: if (!rank) {
30: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
31: }
32: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
33: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
34: MatSetOptionsPrefix(A,"A_");
35: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
36: PetscPrintf(PETSC_COMM_WORLD,"\n");
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: MatMatMultNumeric(B,A,C); /* recompute C=B*A */
42: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C); /* recompute C=B*A */
43: MatSetOptionsPrefix(C,"C_");
44: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
45: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
47: MatMatMultSymbolic(C,A,fill,&D);
48: MatMatMultNumeric(C,A,D); /* D = C*A = (A^T*A)*A */
49: MatSetOptionsPrefix(D,"D_");
50: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
51: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
53: /* Repeat the numeric product to test reuse of the previous symbolic product */
54: MatMatMultNumeric(C,A,D);
55: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
56: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
58: MatDestroy(&B);
59: MatDestroy(&C);
61: /* Test PtAP routine. */
62: MatDuplicate(A,MAT_COPY_VALUES,&B); /* B = A */
63: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
64: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
65: MatNorm(D,NORM_FROBENIUS,&norm);
66: if (norm > 1.e-15 && !rank) {
67: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);
68: }
69: MatDestroy(&C);
70: MatDestroy(&D);
72: /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
73: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
74: MatSetOptionsPrefix(C,"C=BtAB_");
75: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
76: PetscPrintf(PETSC_COMM_WORLD,"\n");
78: MatPtAPSymbolic(A,B,fill,&D);
79: MatPtAPNumeric(A,B,D);
80: MatSetOptionsPrefix(D,"D=BtAB_");
81: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
82: PetscPrintf(PETSC_COMM_WORLD,"\n");
84: /* Repeat numeric product to test reuse of the previous symbolic product */
85: MatPtAPNumeric(A,B,D);
86: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
87: MatNorm(D,NORM_FROBENIUS,&norm);
88: if (norm > 1.e-15) {
89: PetscPrintf(PETSC_COMM_WORLD,"Error in symbolic/numeric MatPtAP: %g\n",norm);
90: }
91: MatDestroy(&B);
92: MatDestroy(&C);
93: MatDestroy(&D);
95: if (size == 1) {
96: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
97: testPTAPRectangular();
99: /* test MatMatTransposeMult(): A*B^T */
100: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
101: MatSetOptionsPrefix(D,"D=A*A^T_");
102: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
103: PetscPrintf(PETSC_COMM_WORLD,"\n");
105: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
106: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C=A*B */
107: MatSetOptionsPrefix(C,"D=A*B=A*A^T_");
108: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
109: PetscPrintf(PETSC_COMM_WORLD,"\n");
110: }
112: MatDestroy(&A);
113: MatDestroy(&B);
114: MatDestroy(&C);
115: MatDestroy(&D);
116: PetscFinalize();
117: return(0);
118: }
120: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
121: #define PETSc_CHKERRQ CHKERRQ
124: PetscErrorCode testPTAPRectangular(void)
125: {
127: const int rows = 3;
128: const int cols = 5;
129: PetscErrorCode _ierr;
130: int i;
131: Mat A;
132: Mat P;
133: Mat C;
136: /* set up A */
137: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
138: PETSc_CHKERRQ(_ierr);
139: for (i=0; i<rows; i++) {
140: _MatSetValue(A, i, i, 1.0, INSERT_VALUES);
141: PETSc_CHKERRQ(_ierr);
142: }
143: _MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
144: PETSc_CHKERRQ(_ierr);
145: _MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
146: PETSc_CHKERRQ(_ierr);
148: /* set up P */
149: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
150: PETSc_CHKERRQ(_ierr);
151: _MatSetValue(P, 0, 0, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
152: _MatSetValue(P, 0, 1, 2.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
153: _MatSetValue(P, 0, 2, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
155: _MatSetValue(P, 0, 3, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
157: _MatSetValue(P, 1, 0, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
158: _MatSetValue(P, 1, 1, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
159: _MatSetValue(P, 1, 2, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
161: _MatSetValue(P, 2, 0, 3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
162: _MatSetValue(P, 2, 1, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
163: _MatSetValue(P, 2, 2, -3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
165: _MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
166: PETSc_CHKERRQ(_ierr);
167: _MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
168: PETSc_CHKERRQ(_ierr);
170: /* compute C */
171: _MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
172: PETSc_CHKERRQ(_ierr);
174: _MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
175: PETSc_CHKERRQ(_ierr);
176: _MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
177: PETSc_CHKERRQ(_ierr);
179: /* compare results */
180: /*
181: printf("C:\n");
182: _MatView(C,PETSC_VIEWER_STDOUT_WORLD);PETSc_CHKERRQ(_ierr);
184: blitz::Array<double,2> actualC(cols, cols);
185: actualC = 0.0;
186: for (int i=0; i<cols; i++) {
187: for (int j=0; j<cols; j++) {
188: _MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
189: PETSc_CHKERRQ(_ierr); ;
190: }
191: }
192: blitz::Array<double,2> expectedC(cols, cols);
193: expectedC = 0.0;
195: expectedC(0,0) = 10.0;
196: expectedC(0,1) = 2.0;
197: expectedC(0,2) = -9.0;
198: expectedC(0,3) = -1.0;
199: expectedC(1,0) = 2.0;
200: expectedC(1,1) = 5.0;
201: expectedC(1,2) = -1.0;
202: expectedC(1,3) = -2.0;
203: expectedC(2,0) = -9.0;
204: expectedC(2,1) = -1.0;
205: expectedC(2,2) = 10.0;
206: expectedC(2,3) = 0.0;
207: expectedC(3,0) = -1.0;
208: expectedC(3,1) = -2.0;
209: expectedC(3,2) = 0.0;
210: expectedC(3,3) = 1.0;
212: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
213: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
214: */
216: _MatDestroy(&A);
217: PETSc_CHKERRQ(_ierr);
218: _MatDestroy(&P);
219: PETSc_CHKERRQ(_ierr);
220: _MatDestroy(&C);
221: PETSc_CHKERRQ(_ierr);
222: return(0);
223: }