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;
19: PetscInitialize(&argc, &args, (char *)0, help);
20: MPI_Comm_size(PETSC_COMM_WORLD, &size);
22: MatCreate(PETSC_COMM_WORLD, &C);
23: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
24: MatSetFromOptions(C);
25: MatSetUp(C);
26: PetscOptionsHasName(NULL, NULL, "-symmetric", &flg);
27: if (flg) { /* Treat matrix as symmetric only if we set this flag */
28: MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE);
29: 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: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
40: }
41: if (i < m - 1) {
42: J = Ii + n;
43: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
44: }
45: if (j > 0) {
46: J = Ii - 1;
47: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
48: }
49: if (j < n - 1) {
50: J = Ii + 1;
51: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
52: }
53: v = 4.0;
54: MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
55: }
56: }
57: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
59: MatGetOrdering(C, MATORDERINGRCM, &perm, &iperm);
60: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
61: ISView(perm, PETSC_VIEWER_STDOUT_SELF);
62: VecCreateSeq(PETSC_COMM_SELF, m * n, &u);
63: VecSet(u, 1.0);
64: VecDuplicate(u, &x);
65: VecDuplicate(u, &b);
66: VecDuplicate(u, &y);
67: MatMult(C, u, b);
68: VecCopy(b, y);
69: VecScale(y, 2.0);
71: MatNorm(C, NORM_FROBENIUS, &norm);
72: PetscPrintf(PETSC_COMM_SELF, "Frobenius norm of matrix %g\n", (double)norm);
73: MatNorm(C, NORM_1, &norm);
74: PetscPrintf(PETSC_COMM_SELF, "One norm of matrix %g\n", (double)norm);
75: MatNorm(C, NORM_INFINITY, &norm);
76: PetscPrintf(PETSC_COMM_SELF, "Infinity norm of matrix %g\n", (double)norm);
78: MatFactorInfoInitialize(&info);
79: info.fill = 2.0;
80: info.dtcol = 0.0;
81: info.zeropivot = 1.e-14;
82: info.pivotinblocks = 1.0;
84: MatLUFactor(C, perm, iperm, &info);
86: /* Test MatSolve */
87: MatSolve(C, b, x);
88: VecView(b, PETSC_VIEWER_STDOUT_SELF);
89: VecView(x, PETSC_VIEWER_STDOUT_SELF);
90: VecAXPY(x, -1.0, u);
91: VecNorm(x, NORM_2, &norm);
92: if (norm > tol) PetscPrintf(PETSC_COMM_SELF, "MatSolve: Norm of error %g\n", (double)norm);
94: /* Test MatSolveAdd */
95: MatSolveAdd(C, b, y, x);
96: VecAXPY(x, -1.0, y);
97: VecAXPY(x, -1.0, u);
98: VecNorm(x, NORM_2, &norm);
99: if (norm > tol) PetscPrintf(PETSC_COMM_SELF, "MatSolveAdd(): Norm of error %g\n", (double)norm);
101: ISDestroy(&perm);
102: ISDestroy(&iperm);
103: VecDestroy(&u);
104: VecDestroy(&y);
105: VecDestroy(&b);
106: VecDestroy(&x);
107: MatDestroy(&C);
108: PetscFinalize();
109: return 0;
110: }
112: /*TEST
114: test:
116: TEST*/