Actual source code: ex117.c
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: */
7: #include <petscmat.h>
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;
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!");
24: nnz[0] = 2;
25: nnz[1] = 1;
26: nnz[2] = 1;
27: PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_SELF, 2, 6, 6, 0, nnz, &mat));
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));
58: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
59: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
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));
91: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
92: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
94: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mat: \n"));
95: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_SELF));
97: /* begin cholesky factorization */
98: PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &perm, &colp));
99: PetscCall(ISDestroy(&colp));
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));
106: PetscCall(ISDestroy(&perm));
107: PetscCall(MatDestroy(&mat));
108: PetscCall(MatDestroy(&fact));
109: PetscCall(MatDestroy(&B));
110: PetscCall(PetscFinalize());
111: return 0;
112: }