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";
4
5 #include <petscmat.h>
6
main(int argc,char ** args)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;
17
18 PetscFunctionBeginUser;
19 PetscCall(PetscInitialize(&argc, &args, NULL, help));
20
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));
28
29 PetscOptionsEnd();
30
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));
65
66 PetscCall(MatFactorInfoInitialize(&luinfo));
67
68 luinfo.fill = fill;
69 luinfo.dtcol = 0.0;
70 luinfo.zeropivot = 1.e-14;
71 luinfo.pivotinblocks = 1.0;
72
73 PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &LU));
74 PetscCall(MatLUFactorSymbolic(LU, C, perm, iperm, &luinfo));
75 PetscCall(MatLUFactorNumeric(LU, C, &luinfo));
76
77 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m * n, &u));
78 PetscCall(VecSet(u, one));
79 PetscCall(VecDuplicate(u, &x));
80 PetscCall(VecDuplicate(u, &b));
81
82 PetscCall(MatMult(C, u, b));
83 PetscCall(MatSolve(LU, b, x));
84
85 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_SELF));
86 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_SELF));
87
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));
91
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));
96
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));
104
105 PetscCall(PetscFinalize());
106 return 0;
107 }
108
109 /*TEST
110
111 test:
112 suffix: 1
113 filter: grep -v " MPI process"
114
115 test:
116 suffix: 2
117 args: -m 1 -n 1 -fill 0.49
118 filter: grep -v " MPI process"
119
120 TEST*/
121