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