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: }