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