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