Actual source code: ex93.c
1: static char help[] = "Test sequential MatMatMult() and MatPtAP() for AIJ matrices.\n\n";
3: #include petscmat.h
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;
17: PetscInitialize(&argc,&argv,(char *)0,help);
18: MatCreate(PETSC_COMM_SELF,&A);
19: MatSetSizes(A,3,3,3,3);
20: MatSetType(A,MATSEQAIJ);
21: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
23: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
24: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
25: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
27: /* Form A^T*A*A to test PtAP routine. */
28: MatTranspose(A,MAT_INITIAL_MATRIX,&B);
29: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C);
30: MatMatMultSymbolic(C,A,fill,&D);
31: MatMatMultNumeric(C,A,D);
32: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
34: /* Repeat the numeric product to test reuse of the previous symbolic product */
35: MatMatMultNumeric(C,A,D);
36: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
38: MatDestroy(B);
39: MatDestroy(C);
41: MatDuplicate(A,MAT_COPY_VALUES,&B);
42: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
44: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
45: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
47: MatDestroy(C);
48: MatDestroy(D);
50: /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
51: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
52: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
53: MatPtAPSymbolic(A,B,fill,&D);
54: MatPtAPNumeric(A,B,D);
55: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
57: /* Repeat numeric product to test reuse of the previous symbolic product */
58: MatPtAPNumeric(A,B,D);
59: MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
60: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
62: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
63: testPTAPRectangular();
65: MatDestroy(A);
66: MatDestroy(B);
67: MatDestroy(C);
68: MatDestroy(D);
69: PetscFinalize();
70: return(0);
71: }
73: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
74: #define PETSc_CHKERRQ CHKERRQ
77: PetscErrorCode testPTAPRectangular(void)
78: {
80: const int rows = 3;
81: const int cols = 5;
82: PetscErrorCode _ierr;
83: int i;
84: Mat A;
85: Mat P;
86: Mat C;
89: /* set up A */
90: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,
91: 1, PETSC_NULL, &A);
92: PETSc_CHKERRQ(_ierr);
93: for (i=0; i<rows; i++) {
94: _MatSetValue(A, i, i, 1.0, INSERT_VALUES);
95: PETSc_CHKERRQ(_ierr);
96: }
97: _MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
98: PETSc_CHKERRQ(_ierr);
99: _MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
100: PETSc_CHKERRQ(_ierr);
102: /* set up P */
103: _MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,
104: 5, PETSC_NULL, &P);
105: PETSc_CHKERRQ(_ierr);
106: _MatSetValue(P, 0, 0, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
107: _MatSetValue(P, 0, 1, 2.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
108: _MatSetValue(P, 0, 2, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
110: _MatSetValue(P, 0, 3, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
112: _MatSetValue(P, 1, 0, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
113: _MatSetValue(P, 1, 1, -1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
114: _MatSetValue(P, 1, 2, 1.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
116: _MatSetValue(P, 2, 0, 3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
117: _MatSetValue(P, 2, 1, 0.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
118: _MatSetValue(P, 2, 2, -3.0, INSERT_VALUES); PETSc_CHKERRQ(_ierr);
119:
120: _MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
121: PETSc_CHKERRQ(_ierr);
122: _MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
123: PETSc_CHKERRQ(_ierr);
125: /* compute C */
126: _MatPtAP( A, P, MAT_INITIAL_MATRIX, 1.0, &C);
127: PETSc_CHKERRQ(_ierr);
129: _MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
130: PETSc_CHKERRQ(_ierr);
131: _MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
132: PETSc_CHKERRQ(_ierr);
134: /* compare results */
135: /*
136: printf("C:\n");
137: _MatView(C,PETSC_VIEWER_STDOUT_WORLD);PETSc_CHKERRQ(_ierr);
139: blitz::Array<double,2> actualC(cols, cols);
140: actualC = 0.0;
141: for (int i=0; i<cols; i++) {
142: for (int j=0; j<cols; j++) {
143: _MatGetValues(C, 1, &i, 1, &j, &actualC(i,j) );
144: PETSc_CHKERRQ(_ierr); ;
145: }
146: }
147: blitz::Array<double,2> expectedC(cols, cols);
148: expectedC = 0.0;
149:
150: expectedC(0,0) = 10.0;
151: expectedC(0,1) = 2.0;
152: expectedC(0,2) = -9.0;
153: expectedC(0,3) = -1.0;
154: expectedC(1,0) = 2.0;
155: expectedC(1,1) = 5.0;
156: expectedC(1,2) = -1.0;
157: expectedC(1,3) = -2.0;
158: expectedC(2,0) = -9.0;
159: expectedC(2,1) = -1.0;
160: expectedC(2,2) = 10.0;
161: expectedC(2,3) = 0.0;
162: expectedC(3,0) = -1.0;
163: expectedC(3,1) = -2.0;
164: expectedC(3,2) = 0.0;
165: expectedC(3,3) = 1.0;
166:
167: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
168: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
169: */
170:
171: _MatDestroy(A);
172: PETSc_CHKERRQ(_ierr);
173: _MatDestroy(P);
174: PETSc_CHKERRQ(_ierr);
175: _MatDestroy(C);
176: PETSc_CHKERRQ(_ierr);
177: return(0);
178: }