1: static char help[] = "Test sequential MatMatMult() and MatPtAP() for AIJ matrices.\n\n";

 3:  #include petscmat.h

  7: int main(int argc,char **argv) {
  8:   Mat            A,B,C,D;
  9:   PetscScalar    a[]={1.,1.,0.,0.,1.,1.,0.,0.,1.};
 10:   PetscInt       ij[]={0,1,2};
 11:   PetscScalar    none=-1.;
 13:   PetscReal      fill=4;

 15:   PetscInitialize(&argc,&argv,(char *)0,help);
 16:   MatCreate(PETSC_COMM_SELF,&A);
 17:   MatSetSizes(A,3,3,3,3);
 18:   MatSetType(A,MATSEQAIJ);
 19:   MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES);

 21:   MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
 22:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 23:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 25:   /* Form A^T*A*A to test PtAP routine. */
 26:   MatTranspose(A,&B);
 27:   MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C);
 28:   MatMatMultSymbolic(C,A,fill,&D);
 29:   MatMatMultNumeric(C,A,D);
 30:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);

 32:   /* Repeat the numeric product to test reuse of the previous symbolic product */
 33:   MatMatMultNumeric(C,A,D);
 34:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);

 36:   MatDestroy(B);
 37:   MatDestroy(C);

 39:   MatDuplicate(A,MAT_COPY_VALUES,&B);
 40:   MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);

 42:   MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
 43:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);

 45:   MatDestroy(C);
 46:   MatDestroy(D);

 48:   /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
 49:   MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
 50:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
 51:   MatPtAPSymbolic(A,B,fill,&D);
 52:   MatPtAPNumeric(A,B,D);
 53:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);

 55:   /* Repeat numeric product to test reuse of the previous symbolic product */
 56:   MatPtAPNumeric(A,B,D);
 57:   MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
 58:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);

 60:   MatDestroy(A);
 61:   MatDestroy(B);
 62:   MatDestroy(C);
 63:   MatDestroy(D);
 64:   PetscFinalize();
 65:   return(0);
 66: }