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
main(int argc,char ** args)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