Actual source code: ex93.c
petsc-3.7.3 2016-08-01
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: MatSetFromOptions(A);
28: MatSetUp(A);
29: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
31: if (!rank) {
32: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
33: }
34: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
35: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
36: MatSetOptionsPrefix(A,"A_");
37: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
38: PetscPrintf(PETSC_COMM_WORLD,"\n");
40: /* Test MatMatMult() */
41: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
42: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A */
43: MatMatMultNumeric(B,A,C); /* recompute C=B*A */
44: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C); /* recompute C=B*A */
45: MatSetOptionsPrefix(C,"C_");
46: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
47: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
49: MatMatMultSymbolic(C,A,fill,&D);
50: MatMatMultNumeric(C,A,D); /* D = C*A = (A^T*A)*A */
51: MatSetOptionsPrefix(D,"D_");
52: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
53: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
55: /* Repeat the numeric product to test reuse of the previous symbolic product */
56: MatMatMultNumeric(C,A,D);
57: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
58: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}
60: MatDestroy(&B);
61: MatDestroy(&C);
63: /* Test PtAP routine. */
64: MatDuplicate(A,MAT_COPY_VALUES,&B); /* B = A */
65: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
66: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
67: MatNorm(D,NORM_FROBENIUS,&norm);
68: if (norm > 1.e-15 && !rank) {
69: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);
70: }
71: MatDestroy(&C);
72: MatDestroy(&D);
74: /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
75: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
76: MatSetOptionsPrefix(C,"C=BtAB_");
77: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
78: PetscPrintf(PETSC_COMM_WORLD,"\n");
80: MatPtAPSymbolic(A,B,fill,&D);
81: MatPtAPNumeric(A,B,D);
82: MatSetOptionsPrefix(D,"D=BtAB_");
83: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
84: PetscPrintf(PETSC_COMM_WORLD,"\n");
86: /* Repeat numeric product to test reuse of the previous symbolic product */
87: MatPtAPNumeric(A,B,D);
88: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
89: MatNorm(D,NORM_FROBENIUS,&norm);
90: if (norm > 1.e-15) {
91: PetscPrintf(PETSC_COMM_WORLD,"Error in symbolic/numeric MatPtAP: %g\n",norm);
92: }
93: MatDestroy(&B);
94: MatDestroy(&C);
95: MatDestroy(&D);
97: if (size == 1) {
98: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
99: testPTAPRectangular();
101: /* test MatMatTransposeMult(): A*B^T */
102: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
103: MatSetOptionsPrefix(D,"D=A*A^T_");
104: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
105: PetscPrintf(PETSC_COMM_WORLD,"\n");
107: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
108: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C=A*B */
109: MatSetOptionsPrefix(C,"D=A*B=A*A^T_");
110: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
111: PetscPrintf(PETSC_COMM_WORLD,"\n");
112: }
114: MatDestroy(&A);
115: MatDestroy(&B);
116: MatDestroy(&C);
117: MatDestroy(&D);
118: PetscFinalize();
119: return(0);
120: }
122: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
123: #define PETSc_CHKERRQ CHKERRQ
126: PetscErrorCode testPTAPRectangular(void)
127: {
129: const int rows = 3;
130: const int cols = 5;
131: PetscErrorCode _ierr;
132: int i;
133: Mat A;
134: Mat P;
135: Mat C;
138: /* set up A */
139: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
140: PETSc_CHKERRQ(_ierr);
141: for (i=0; i<rows; i++) {
142: _MatSetValue(A, i, i, 1.0, INSERT_VALUES);
143: PETSc_CHKERRQ(_ierr);
144: }
145: _MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
146: PETSc_CHKERRQ(_ierr);
147: _MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
148: PETSc_CHKERRQ(_ierr);
150: /* set up P */
151: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
152: PETSc_CHKERRQ(_ierr);
153: _MatSetValue(P, 0, 0, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
154: _MatSetValue(P, 0, 1, 2.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
155: _MatSetValue(P, 0, 2, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
157: _MatSetValue(P, 0, 3, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
159: _MatSetValue(P, 1, 0, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
160: _MatSetValue(P, 1, 1, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
161: _MatSetValue(P, 1, 2, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
163: _MatSetValue(P, 2, 0, 3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
164: _MatSetValue(P, 2, 1, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
165: _MatSetValue(P, 2, 2, -3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
167: _MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
168: PETSc_CHKERRQ(_ierr);
169: _MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
170: PETSc_CHKERRQ(_ierr);
172: /* compute C */
173: _MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
174: PETSc_CHKERRQ(_ierr);
176: _MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
177: PETSc_CHKERRQ(_ierr);
178: _MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
179: PETSc_CHKERRQ(_ierr);
181: /* compare results */
182: /*
183: printf("C:\n");
184: _MatView(C,PETSC_VIEWER_STDOUT_WORLD);PETSc_CHKERRQ(_ierr);
186: blitz::Array<double,2> actualC(cols, cols);
187: actualC = 0.0;
188: for (int i=0; i<cols; i++) {
189: for (int j=0; j<cols; j++) {
190: _MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
191: PETSc_CHKERRQ(_ierr); ;
192: }
193: }
194: blitz::Array<double,2> expectedC(cols, cols);
195: expectedC = 0.0;
197: expectedC(0,0) = 10.0;
198: expectedC(0,1) = 2.0;
199: expectedC(0,2) = -9.0;
200: expectedC(0,3) = -1.0;
201: expectedC(1,0) = 2.0;
202: expectedC(1,1) = 5.0;
203: expectedC(1,2) = -1.0;
204: expectedC(1,3) = -2.0;
205: expectedC(2,0) = -9.0;
206: expectedC(2,1) = -1.0;
207: expectedC(2,2) = 10.0;
208: expectedC(2,3) = 0.0;
209: expectedC(3,0) = -1.0;
210: expectedC(3,1) = -2.0;
211: expectedC(3,2) = 0.0;
212: expectedC(3,3) = 1.0;
214: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
215: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
216: */
218: _MatDestroy(&A);
219: PETSc_CHKERRQ(_ierr);
220: _MatDestroy(&P);
221: PETSc_CHKERRQ(_ierr);
222: _MatDestroy(&C);
223: PETSc_CHKERRQ(_ierr);
224: return(0);
225: }