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;
20: 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: PetscOptionsInt("-m", "Number of rows in grid", NULL, m, &m, NULL);
27: PetscOptionsInt("-n", "Number of columns in grid", NULL, n, &n, NULL);
28: 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: MatCreate(PETSC_COMM_SELF, &C);
34: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
35: MatSetFromOptions(C);
36: 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: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
44: }
45: if (i < m - 1) {
46: J = Ii + n;
47: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
48: }
49: if (j > 0) {
50: J = Ii - 1;
51: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
52: }
53: if (j < n - 1) {
54: J = Ii + 1;
55: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
56: }
57: v = 4.0;
58: MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
59: }
60: }
61: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
63: MatGetOrdering(C, MATORDERINGRCM, &perm, &iperm);
64: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
65: ISView(perm, PETSC_VIEWER_STDOUT_SELF);
67: MatFactorInfoInitialize(&luinfo);
69: luinfo.fill = fill;
70: luinfo.dtcol = 0.0;
71: luinfo.zeropivot = 1.e-14;
72: luinfo.pivotinblocks = 1.0;
74: MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &LU);
75: MatLUFactorSymbolic(LU, C, perm, iperm, &luinfo);
76: MatLUFactorNumeric(LU, C, &luinfo);
78: VecCreateSeq(PETSC_COMM_SELF, m * n, &u);
79: VecSet(u, one);
80: VecDuplicate(u, &x);
81: VecDuplicate(u, &b);
83: MatMult(C, u, b);
84: MatSolve(LU, b, x);
86: VecView(b, PETSC_VIEWER_STDOUT_SELF);
87: VecView(x, PETSC_VIEWER_STDOUT_SELF);
89: VecAXPY(x, -1.0, u);
90: VecNorm(x, NORM_2, &norm);
91: PetscPrintf(PETSC_COMM_SELF, "Norm of error %g\n", (double)norm);
93: MatGetInfo(C, MAT_LOCAL, &info);
94: PetscPrintf(PETSC_COMM_SELF, "original matrix nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used);
95: MatGetInfo(LU, MAT_LOCAL, &info);
96: PetscPrintf(PETSC_COMM_SELF, "factored matrix nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used);
98: VecDestroy(&u);
99: VecDestroy(&b);
100: VecDestroy(&x);
101: ISDestroy(&perm);
102: ISDestroy(&iperm);
103: MatDestroy(&C);
104: MatDestroy(&LU);
106: 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*/