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: {
  9:   Mat            pA,P,aijP;
 10:   PetscScalar    pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.};
 11:   PetscInt       i,pij[]={0,1,2};
 12:   PetscInt       aij[3][3]={{0,1,2},{3,4,5},{6,7,8}};
 13:   Mat            A,mC,C;
 14:   PetscScalar    one=1.;
 16:   PetscMPIInt    size;

 18:   PetscInitialize(&argc,&argv,(char*)0,help);
 19:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 20:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

 22:   /* Create MAIJ matrix, P */
 23:   MatCreate(PETSC_COMM_SELF,&pA);
 24:   MatSetSizes(pA,3,3,3,3);
 25:   MatSetType(pA,MATSEQAIJ);
 26:   MatSetUp(pA);
 27:   MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
 28:   MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES);
 29:   MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY);
 30:   MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY);
 31:   MatCreateMAIJ(pA,3,&P);
 32:   MatDestroy(&pA);

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

 37:   /* Create AIJ matrix, A */
 38:   MatCreate(PETSC_COMM_SELF,&A);
 39:   MatSetSizes(A,9,9,9,9);
 40:   MatSetType(A,MATSEQAIJ);
 41:   MatSetUp(A);
 42:   MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
 43:   MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES);
 44:   MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES);
 45:   MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES);
 46:   for (i=0; i<9; i++) {
 47:     MatSetValue(A,i,i,one,ADD_VALUES);
 48:   }
 49:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 50:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 52:   /* Perform PtAP_SeqAIJ_SeqMAIJ */
 53:   MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
 54:   MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC);
 55:   MatView(mC,PETSC_VIEWER_STDOUT_SELF);

 57:   /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */
 58:   MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);
 59:   MatView(C,PETSC_VIEWER_STDOUT_SELF);

 61:   /* Check mC = C */
 62:   MatAXPY(C,-1.0,mC,DIFFERENT_NONZERO_PATTERN);
 63:   /* Note: We should be able to use SAME_NONZERO_PATTERN on the line above, */
 64:   /*       but don't because this flag doesn't assist testing. */
 65:   MatView(C,PETSC_VIEWER_STDOUT_SELF);

 67:   /* Cleanup */
 68:   MatDestroy(&P);
 69:   MatDestroy(&aijP);
 70:   MatDestroy(&A);
 71:   MatDestroy(&C);
 72:   MatDestroy(&mC);
 73:   PetscFinalize();
 74:   return(0);
 75: }