Actual source code: ex15.c
2: static char help[] = "Tests MatNorm(), MatLUFactor(), MatSolve() and MatSolveAdd().\n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat C;
9: PetscInt i, j, m = 3, n = 3, Ii, J;
10: PetscBool flg;
11: PetscScalar v;
12: IS perm, iperm;
13: Vec x, u, b, y;
14: PetscReal norm, tol = PETSC_SMALL;
15: MatFactorInfo info;
16: PetscMPIInt size;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, (char *)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!");
22: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
23: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
24: PetscCall(MatSetFromOptions(C));
25: PetscCall(MatSetUp(C));
26: PetscCall(PetscOptionsHasName(NULL, NULL, "-symmetric", &flg));
27: if (flg) { /* Treat matrix as symmetric only if we set this flag */
28: PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE));
29: PetscCall(MatSetOption(C, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
30: }
32: /* Create the matrix for the five point stencil, YET AGAIN */
33: for (i = 0; i < m; i++) {
34: for (j = 0; j < n; j++) {
35: v = -1.0;
36: Ii = j + n * i;
37: if (i > 0) {
38: J = Ii - n;
39: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
40: }
41: if (i < m - 1) {
42: J = Ii + n;
43: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
44: }
45: if (j > 0) {
46: J = Ii - 1;
47: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
48: }
49: if (j < n - 1) {
50: J = Ii + 1;
51: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
52: }
53: v = 4.0;
54: PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
55: }
56: }
57: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
58: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
59: PetscCall(MatGetOrdering(C, MATORDERINGRCM, &perm, &iperm));
60: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
61: PetscCall(ISView(perm, PETSC_VIEWER_STDOUT_SELF));
62: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m * n, &u));
63: PetscCall(VecSet(u, 1.0));
64: PetscCall(VecDuplicate(u, &x));
65: PetscCall(VecDuplicate(u, &b));
66: PetscCall(VecDuplicate(u, &y));
67: PetscCall(MatMult(C, u, b));
68: PetscCall(VecCopy(b, y));
69: PetscCall(VecScale(y, 2.0));
71: PetscCall(MatNorm(C, NORM_FROBENIUS, &norm));
72: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Frobenius norm of matrix %g\n", (double)norm));
73: PetscCall(MatNorm(C, NORM_1, &norm));
74: PetscCall(PetscPrintf(PETSC_COMM_SELF, "One norm of matrix %g\n", (double)norm));
75: PetscCall(MatNorm(C, NORM_INFINITY, &norm));
76: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Infinity norm of matrix %g\n", (double)norm));
78: PetscCall(MatFactorInfoInitialize(&info));
79: info.fill = 2.0;
80: info.dtcol = 0.0;
81: info.zeropivot = 1.e-14;
82: info.pivotinblocks = 1.0;
84: PetscCall(MatLUFactor(C, perm, iperm, &info));
86: /* Test MatSolve */
87: PetscCall(MatSolve(C, b, x));
88: PetscCall(VecView(b, PETSC_VIEWER_STDOUT_SELF));
89: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_SELF));
90: PetscCall(VecAXPY(x, -1.0, u));
91: PetscCall(VecNorm(x, NORM_2, &norm));
92: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve: Norm of error %g\n", (double)norm));
94: /* Test MatSolveAdd */
95: PetscCall(MatSolveAdd(C, b, y, x));
96: PetscCall(VecAXPY(x, -1.0, y));
97: PetscCall(VecAXPY(x, -1.0, u));
98: PetscCall(VecNorm(x, NORM_2, &norm));
99: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolveAdd(): Norm of error %g\n", (double)norm));
101: PetscCall(ISDestroy(&perm));
102: PetscCall(ISDestroy(&iperm));
103: PetscCall(VecDestroy(&u));
104: PetscCall(VecDestroy(&y));
105: PetscCall(VecDestroy(&b));
106: PetscCall(VecDestroy(&x));
107: PetscCall(MatDestroy(&C));
108: PetscCall(PetscFinalize());
109: return 0;
110: }
112: /*TEST
114: test:
116: TEST*/