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

  3:  #include <petscmat.h>

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

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

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

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

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

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

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

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

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