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