1 2 static char help[] = "Tests Cholesky factorization for a SBAIJ matrix, (bs=2).\n"; 3 /* 4 This code is modified from the code contributed by JUNWANG@uwm.edu on Apr 13, 2007 5 */ 6 7 #include <petscmat.h> 8 9 int main(int argc, char **args) 10 { 11 Mat mat, fact, B; 12 PetscInt ind1[2], ind2[2]; 13 PetscScalar temp[4]; 14 PetscInt nnz[3]; 15 IS perm, colp; 16 MatFactorInfo info; 17 PetscMPIInt size; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &args, 0, help)); 21 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 22 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 23 24 nnz[0] = 2; 25 nnz[1] = 1; 26 nnz[2] = 1; 27 PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_SELF, 2, 6, 6, 0, nnz, &mat)); 28 29 ind1[0] = 0; 30 ind1[1] = 1; 31 temp[0] = 3; 32 temp[1] = 2; 33 temp[2] = 0; 34 temp[3] = 3; 35 PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES)); 36 ind2[0] = 4; 37 ind2[1] = 5; 38 temp[0] = 1; 39 temp[1] = 1; 40 temp[2] = 2; 41 temp[3] = 1; 42 PetscCall(MatSetValues(mat, 2, ind1, 2, ind2, temp, INSERT_VALUES)); 43 ind1[0] = 2; 44 ind1[1] = 3; 45 temp[0] = 4; 46 temp[1] = 1; 47 temp[2] = 1; 48 temp[3] = 5; 49 PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES)); 50 ind1[0] = 4; 51 ind1[1] = 5; 52 temp[0] = 5; 53 temp[1] = 1; 54 temp[2] = 1; 55 temp[3] = 6; 56 PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES)); 57 58 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 59 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 60 61 PetscCall(MatDuplicate(mat, MAT_SHARE_NONZERO_PATTERN, &B)); 62 ind1[0] = 0; 63 ind1[1] = 1; 64 temp[0] = 3; 65 temp[1] = 2; 66 temp[2] = 0; 67 temp[3] = 3; 68 PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES)); 69 ind2[0] = 4; 70 ind2[1] = 5; 71 temp[0] = 1; 72 temp[1] = 1; 73 temp[2] = 2; 74 temp[3] = 1; 75 PetscCall(MatSetValues(mat, 2, ind1, 2, ind2, temp, INSERT_VALUES)); 76 ind1[0] = 2; 77 ind1[1] = 3; 78 temp[0] = 4; 79 temp[1] = 1; 80 temp[2] = 1; 81 temp[3] = 5; 82 PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES)); 83 ind1[0] = 4; 84 ind1[1] = 5; 85 temp[0] = 5; 86 temp[1] = 1; 87 temp[2] = 1; 88 temp[3] = 6; 89 PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES)); 90 91 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 93 94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mat: \n")); 95 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_SELF)); 96 97 /* begin cholesky factorization */ 98 PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &perm, &colp)); 99 PetscCall(ISDestroy(&colp)); 100 101 info.fill = 1.0; 102 PetscCall(MatGetFactor(mat, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &fact)); 103 PetscCall(MatCholeskyFactorSymbolic(fact, mat, perm, &info)); 104 PetscCall(MatCholeskyFactorNumeric(fact, mat, &info)); 105 106 PetscCall(ISDestroy(&perm)); 107 PetscCall(MatDestroy(&mat)); 108 PetscCall(MatDestroy(&fact)); 109 PetscCall(MatDestroy(&B)); 110 PetscCall(PetscFinalize()); 111 return 0; 112 } 113