Actual source code: ex93.c
petsc-3.3-p6 2013-02-11
1: static char help[] = "Test sequential 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: Mat A,B,C,D;
11: PetscScalar a[]={1.,1.,0.,0.,1.,1.,0.,0.,1.};
12: PetscInt ij[]={0,1,2};
13: PetscScalar none=-1.;
15: PetscReal fill=4;
16: PetscReal norm;
18: PetscInitialize(&argc,&argv,(char *)0,help);
19: MatCreate(PETSC_COMM_SELF,&A);
20: MatSetSizes(A,3,3,3,3);
21: MatSetType(A,MATSEQAIJ);
22: MatSetUp(A);
23: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
24: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
25: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
26: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
27: MatSetOptionsPrefix(A,"A_");
28: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
29: PetscPrintf(PETSC_COMM_SELF,"\n");
31: /* Test MatMatMult() */
32: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
33: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A */
34: MatSetOptionsPrefix(C,"C=B*A=A^T*A_");
35: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
36: PetscPrintf(PETSC_COMM_SELF,"\n");
38: MatMatMultSymbolic(C,A,fill,&D);
39: MatMatMultNumeric(C,A,D); /* D = C*A = (A^T*A)*A */
40: MatSetOptionsPrefix(D,"D=C*A=(A^T*A)*A_");
41: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
42: PetscPrintf(PETSC_COMM_SELF,"\n");
44: /* Repeat the numeric product to test reuse of the previous symbolic product */
45: MatMatMultNumeric(C,A,D);
46: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
47: PetscPrintf(PETSC_COMM_SELF,"\n");
49: MatDestroy(&B);
50: MatDestroy(&C);
52: /* Test PtAP routine. */
53: MatDuplicate(A,MAT_COPY_VALUES,&B); /* B = A */
54: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
55: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
56: MatNorm(D,NORM_FROBENIUS,&norm);
57: if (norm > 1.e-15){
58: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);
59: }
60: MatDestroy(&C);
61: MatDestroy(&D);
63: /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
64: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
65: MatSetOptionsPrefix(C,"C=BtAB_");
66: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
67: PetscPrintf(PETSC_COMM_SELF,"\n");
69: MatPtAPSymbolic(A,B,fill,&D);
70: MatPtAPNumeric(A,B,D);
71: MatSetOptionsPrefix(D,"D=BtAB_");
72: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
73: PetscPrintf(PETSC_COMM_SELF,"\n");
75: /* Repeat numeric product to test reuse of the previous symbolic product */
76: MatPtAPNumeric(A,B,D);
77: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
78: MatNorm(D,NORM_FROBENIUS,&norm);
79: if (norm > 1.e-15){
80: PetscPrintf(PETSC_COMM_SELF,"Error in symbolic/numeric MatPtAP: %g\n",norm);
81: }
82: MatDestroy(&B);
83: MatDestroy(&C);
84: MatDestroy(&D);
86: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
87: testPTAPRectangular();
89: /* test MatMatTransposeMult(): A*B^T */
90: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
91: MatSetOptionsPrefix(D,"D=A*A^T_");
92: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
93: PetscPrintf(PETSC_COMM_SELF,"\n");
95: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
96: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C=A*B */
97: MatSetOptionsPrefix(C,"D=A*B=A*A^T_");
98: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
99: PetscPrintf(PETSC_COMM_SELF,"\n");
101: MatDestroy(&A);
102: MatDestroy(&B);
103: MatDestroy(&C);
104: MatDestroy(&D);
105: PetscFinalize();
106: return(0);
107: }
109: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
110: #define PETSc_CHKERRQ CHKERRQ
113: PetscErrorCode testPTAPRectangular(void)
114: {
116: const int rows = 3;
117: const int cols = 5;
118: PetscErrorCode _ierr;
119: int i;
120: Mat A;
121: Mat P;
122: Mat C;
125: /* set up A */
126: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,
127: 1, PETSC_NULL, &A);
128: PETSc_CHKERRQ(_ierr);
129: for (i=0; i<rows; i++) {
130: _MatSetValue(A, i, i, 1.0, INSERT_VALUES);
131: PETSc_CHKERRQ(_ierr);
132: }
133: _MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
134: PETSc_CHKERRQ(_ierr);
135: _MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
136: PETSc_CHKERRQ(_ierr);
138: /* set up P */
139: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,
140: 5, PETSC_NULL, &P);
141: PETSc_CHKERRQ(_ierr);
142: _MatSetValue(P, 0, 0, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
143: _MatSetValue(P, 0, 1, 2.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
144: _MatSetValue(P, 0, 2, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
146: _MatSetValue(P, 0, 3, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
148: _MatSetValue(P, 1, 0, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
149: _MatSetValue(P, 1, 1, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
150: _MatSetValue(P, 1, 2, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
152: _MatSetValue(P, 2, 0, 3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
153: _MatSetValue(P, 2, 1, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
154: _MatSetValue(P, 2, 2, -3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
155:
156: _MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
157: PETSc_CHKERRQ(_ierr);
158: _MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
159: PETSc_CHKERRQ(_ierr);
161: /* compute C */
162: _MatPtAP( A, P, MAT_INITIAL_MATRIX, 1.0, &C);
163: PETSc_CHKERRQ(_ierr);
165: _MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
166: PETSc_CHKERRQ(_ierr);
167: _MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
168: PETSc_CHKERRQ(_ierr);
170: /* compare results */
171: /*
172: printf("C:\n");
173: _MatView(C,PETSC_VIEWER_STDOUT_WORLD);PETSc_CHKERRQ(_ierr);
175: blitz::Array<double,2> actualC(cols, cols);
176: actualC = 0.0;
177: for (int i=0; i<cols; i++) {
178: for (int j=0; j<cols; j++) {
179: _MatGetValues(C, 1, &i, 1, &j, &actualC(i,j) );
180: PETSc_CHKERRQ(_ierr); ;
181: }
182: }
183: blitz::Array<double,2> expectedC(cols, cols);
184: expectedC = 0.0;
185:
186: expectedC(0,0) = 10.0;
187: expectedC(0,1) = 2.0;
188: expectedC(0,2) = -9.0;
189: expectedC(0,3) = -1.0;
190: expectedC(1,0) = 2.0;
191: expectedC(1,1) = 5.0;
192: expectedC(1,2) = -1.0;
193: expectedC(1,3) = -2.0;
194: expectedC(2,0) = -9.0;
195: expectedC(2,1) = -1.0;
196: expectedC(2,2) = 10.0;
197: expectedC(2,3) = 0.0;
198: expectedC(3,0) = -1.0;
199: expectedC(3,1) = -2.0;
200: expectedC(3,2) = 0.0;
201: expectedC(3,3) = 1.0;
202:
203: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
204: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
205: */
206:
207: _MatDestroy(&A);
208: PETSc_CHKERRQ(_ierr);
209: _MatDestroy(&P);
210: PETSc_CHKERRQ(_ierr);
211: _MatDestroy(&C);
212: PETSc_CHKERRQ(_ierr);
213: return(0);
214: }