1c4762a1bSJed Brown static char help[] = "Testing PtAP for SeqMAIJ matrix, P, with SeqAIJ matrix, A.\n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown #include <petscmat.h> 4c4762a1bSJed Brown 5c4762a1bSJed Brown int main(int argc,char **argv) 6c4762a1bSJed Brown { 7c4762a1bSJed Brown Mat pA,P,aijP; 8c4762a1bSJed Brown PetscScalar pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.}; 9c4762a1bSJed Brown PetscInt i,pij[]={0,1,2}; 10c4762a1bSJed Brown PetscInt aij[3][3]={{0,1,2},{3,4,5},{6,7,8}}; 11c4762a1bSJed Brown Mat A,mC,C; 12c4762a1bSJed Brown PetscScalar one=1.; 13c4762a1bSJed Brown PetscErrorCode ierr; 14c4762a1bSJed Brown PetscMPIInt size; 15c20d7725SJed Brown PetscBool flg; 16c4762a1bSJed Brown 17c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 18*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 192c71b3e2SJacob Faibussowitsch PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 20c4762a1bSJed Brown 21c4762a1bSJed Brown /* Create MAIJ matrix, P */ 22*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_SELF,&pA)); 23*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(pA,3,3,3,3)); 24*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetType(pA,MATSEQAIJ)); 25*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(pA)); 26*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 27*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES)); 28*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY)); 29*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY)); 30*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateMAIJ(pA,3,&P)); 31*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&pA)); 32c4762a1bSJed Brown 33c4762a1bSJed Brown /* Create AIJ equivalent matrix, aijP, for comparison testing */ 34*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP)); 35c4762a1bSJed Brown 36c20d7725SJed Brown /* Create AIJ matrix A */ 37*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_SELF,&A)); 38*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(A,9,9,9,9)); 39*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetType(A,MATSEQAIJ)); 40*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(A)); 41*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 42*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES)); 43*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES)); 44*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES)); 45c4762a1bSJed Brown for (i=0; i<9; i++) { 46*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValue(A,i,i,one,ADD_VALUES)); 47c4762a1bSJed Brown } 48*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 49*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 50c4762a1bSJed Brown 51c4762a1bSJed Brown /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */ 52*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C)); 53c20d7725SJed Brown 54c20d7725SJed Brown /* Perform PtAP_SeqAIJ_SeqMAIJ */ 55c20d7725SJed Brown /* Developer API */ 56*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductCreate(A,P,NULL,&mC)); 57*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetType(mC,MATPRODUCT_PtAP)); 58*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetAlgorithm(mC,"default")); 59*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetFill(mC,PETSC_DEFAULT)); 60*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetFromOptions(mC)); 61*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSymbolic(mC)); 62*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductNumeric(mC)); 63*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductNumeric(mC)); 64c4762a1bSJed Brown 65c4762a1bSJed Brown /* Check mC = C */ 66*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatEqual(C,mC,&flg)); 672c71b3e2SJacob Faibussowitsch PetscCheckFalse(!flg,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"MatProduct C != mC"); 68*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&mC)); 69c20d7725SJed Brown 70c20d7725SJed Brown /* User API */ 71*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC)); 72*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC)); 73c20d7725SJed Brown 74c20d7725SJed Brown /* Check mC = C */ 75*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatEqual(C,mC,&flg)); 762c71b3e2SJacob Faibussowitsch PetscCheckFalse(!flg,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"MatPtAP C != mC"); 77*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&mC)); 78c4762a1bSJed Brown 79c4762a1bSJed Brown /* Cleanup */ 80*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&P)); 81*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&aijP)); 82*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 83*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&C)); 84c4762a1bSJed Brown ierr = PetscFinalize(); 85c4762a1bSJed Brown return ierr; 86c4762a1bSJed Brown } 87c4762a1bSJed Brown 88c4762a1bSJed Brown /*TEST 89c4762a1bSJed Brown 90c4762a1bSJed Brown test: 91c4762a1bSJed Brown output_file: output/ex101.out 92c4762a1bSJed Brown 93c4762a1bSJed Brown TEST*/ 94