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