1: static char help[] = "Testing PtAP for SeqMAIJ matrix, P, with SeqAIJ matrix, A.\n\n";

 3:  #include petscmat.h

  7: int main(int argc,char **argv) {
  8:   Mat            pA,P,aijP;
  9:   PetscScalar    pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.};
 10:   PetscInt       pij[]={0,1,2};
 11:   PetscInt       aij[3][3]={{0,1,2},{3,4,5},{6,7,8}};
 12:   Mat            A,mC,C;
 13:   PetscScalar    one=1.;

 16:   PetscInitialize(&argc,&argv,(char *)0,help);

 18:   /* Create MAIJ matrix, P */
 19:   MatCreate(PETSC_COMM_SELF,&pA);
 20:   MatSetSizes(pA,3,3,3,3);
 21:   MatSetType(pA,MATSEQAIJ);
 22:   MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES);
 23:   MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES);
 24:   MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY);
 25:   MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY);
 26:   MatCreateMAIJ(pA,3,&P);
 27:   MatDestroy(pA);

 29:   /* Create AIJ equivalent matrix, aijP, for comparison testing */
 30:   MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP);

 32:   /* Create AIJ matrix, A */
 33:   MatCreate(PETSC_COMM_SELF,&A);
 34:   MatSetSizes(A,9,9,9,9);
 35:   MatSetType(A,MATSEQAIJ);
 36:   MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES);
 37:   MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES);
 38:   MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES);
 39:   MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES);
 40:   {int i;
 41:    for (i=0;i<9;i++) {
 42:     MatSetValue(A,i,i,one,ADD_VALUES);
 43:    }
 44:   }
 45:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 46:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 47: 
 48:   /* Perform SeqAIJ_SeqMAIJ PtAP */
 49:   MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
 50:   MatView(mC,PETSC_VIEWER_STDOUT_SELF);

 52:   /* Perform SeqAIJ_SeqAIJ PtAP for comparison testing */
 53:   MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);
 54:   MatView(C,PETSC_VIEWER_STDOUT_SELF);

 56:   /* Perform diff of two matrices */
 57:   MatAXPY(C,-1.0,mC,DIFFERENT_NONZERO_PATTERN);
 58:   /* Note: We should be able to use SAME_NONZERO_PATTERN on the line above, */
 59:   /*       but don't because this flag doesn't assist testing. */
 60:   MatView(C,PETSC_VIEWER_STDOUT_SELF);

 62:   /* Cleanup */
 63:   MatDestroy(P);
 64:   MatDestroy(aijP);
 65:   MatDestroy(A);
 66:   MatDestroy(C);
 67:   MatDestroy(mC);
 68:   PetscFinalize();
 69:   return(0);
 70: }