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