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