Actual source code: ex7.c
2: static char help[] = "Tests matrix factorization. Note that most users should\n\
3: employ the KSP interface to the linear solvers instead of using the factorization\n\
4: routines directly.\n\n";
6: #include <petscmat.h>
8: int main(int argc, char **args)
9: {
10: Mat C, LU;
11: MatInfo info;
12: PetscInt i, j, m, n, Ii, J;
13: PetscScalar v, one = 1.0;
14: IS perm, iperm;
15: Vec x, u, b;
16: PetscReal norm, fill;
17: MatFactorInfo luinfo;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
22: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Mat test ex7 options", "Mat");
23: m = 3;
24: n = 3;
25: fill = 2.0;
26: PetscCall(PetscOptionsInt("-m", "Number of rows in grid", NULL, m, &m, NULL));
27: PetscCall(PetscOptionsInt("-n", "Number of columns in grid", NULL, n, &n, NULL));
28: PetscCall(PetscOptionsReal("-fill", "Expected fill ratio for factorization", NULL, fill, &fill, NULL));
30: PetscOptionsEnd();
32: /* Create the matrix for the five point stencil, YET AGAIN */
33: PetscCall(MatCreate(PETSC_COMM_SELF, &C));
34: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
35: PetscCall(MatSetFromOptions(C));
36: PetscCall(MatSetUp(C));
37: for (i = 0; i < m; i++) {
38: for (j = 0; j < n; j++) {
39: v = -1.0;
40: Ii = j + n * i;
41: if (i > 0) {
42: J = Ii - n;
43: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
44: }
45: if (i < m - 1) {
46: J = Ii + n;
47: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
48: }
49: if (j > 0) {
50: J = Ii - 1;
51: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
52: }
53: if (j < n - 1) {
54: J = Ii + 1;
55: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
56: }
57: v = 4.0;
58: PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
59: }
60: }
61: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
62: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
63: PetscCall(MatGetOrdering(C, MATORDERINGRCM, &perm, &iperm));
64: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
65: PetscCall(ISView(perm, PETSC_VIEWER_STDOUT_SELF));
67: PetscCall(MatFactorInfoInitialize(&luinfo));
69: luinfo.fill = fill;
70: luinfo.dtcol = 0.0;
71: luinfo.zeropivot = 1.e-14;
72: luinfo.pivotinblocks = 1.0;
74: PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &LU));
75: PetscCall(MatLUFactorSymbolic(LU, C, perm, iperm, &luinfo));
76: PetscCall(MatLUFactorNumeric(LU, C, &luinfo));
78: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m * n, &u));
79: PetscCall(VecSet(u, one));
80: PetscCall(VecDuplicate(u, &x));
81: PetscCall(VecDuplicate(u, &b));
83: PetscCall(MatMult(C, u, b));
84: PetscCall(MatSolve(LU, b, x));
86: PetscCall(VecView(b, PETSC_VIEWER_STDOUT_SELF));
87: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_SELF));
89: PetscCall(VecAXPY(x, -1.0, u));
90: PetscCall(VecNorm(x, NORM_2, &norm));
91: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Norm of error %g\n", (double)norm));
93: PetscCall(MatGetInfo(C, MAT_LOCAL, &info));
94: PetscCall(PetscPrintf(PETSC_COMM_SELF, "original matrix nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used));
95: PetscCall(MatGetInfo(LU, MAT_LOCAL, &info));
96: PetscCall(PetscPrintf(PETSC_COMM_SELF, "factored matrix nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used));
98: PetscCall(VecDestroy(&u));
99: PetscCall(VecDestroy(&b));
100: PetscCall(VecDestroy(&x));
101: PetscCall(ISDestroy(&perm));
102: PetscCall(ISDestroy(&iperm));
103: PetscCall(MatDestroy(&C));
104: PetscCall(MatDestroy(&LU));
106: PetscCall(PetscFinalize());
107: return 0;
108: }
110: /*TEST
112: test:
113: suffix: 1
114: filter: grep -v " MPI process"
116: test:
117: suffix: 2
118: args: -m 1 -n 1 -fill 0.49
119: filter: grep -v " MPI process"
121: TEST*/