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