Actual source code: ex7.c
2: static char help[] = "Illustrate how to solves a matrix-free linear system with KSP.\n\n";
4: /*
5: Note: modified from ~src/ksp/ksp/tutorials/ex1.c
6: */
7: #include <petscksp.h>
9: /*
10: MatShellMult - Computes the matrix-vector product, y = As x.
12: Input Parameters:
13: As - the matrix-free matrix
14: x - vector
16: Output Parameter:
17: y - vector
18: */
19: PetscErrorCode MyMatShellMult(Mat As, Vec x, Vec y)
20: {
21: Mat P;
23: /* printf("MatShellMult...user should implement this routine without using a matrix\n"); */
24: MatShellGetContext(As, &P);
25: MatMult(P, x, y);
26: return 0;
27: }
29: int main(int argc, char **args)
30: {
31: Vec x, b, u; /* approx solution, RHS, exact solution */
32: Mat P, As; /* preconditioner matrix, linear system (matrix-free) */
33: KSP ksp; /* linear solver context */
34: PC pc; /* preconditioner context */
35: PetscReal norm; /* norm of solution error */
36: PetscInt i, n = 100, col[3], its;
37: PetscMPIInt size;
38: PetscScalar one = 1.0, value[3];
39: PetscBool flg;
42: PetscInitialize(&argc, &args, (char *)0, help);
43: MPI_Comm_size(PETSC_COMM_WORLD, &size);
44: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Compute the matrix and right-hand-side vector that define
48: the linear system, As x = b.
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: /* Create vectors */
51: VecCreate(PETSC_COMM_WORLD, &x);
52: PetscObjectSetName((PetscObject)x, "Solution");
53: VecSetSizes(x, PETSC_DECIDE, n);
54: VecSetFromOptions(x);
55: VecDuplicate(x, &b);
56: VecDuplicate(x, &u);
58: /* Create matrix P, to be used as preconditioner */
59: MatCreate(PETSC_COMM_WORLD, &P);
60: MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, n, n);
61: MatSetFromOptions(P);
62: MatSetUp(P);
64: value[0] = -1.0;
65: value[1] = 2.0;
66: value[2] = -1.0;
67: for (i = 1; i < n - 1; i++) {
68: col[0] = i - 1;
69: col[1] = i;
70: col[2] = i + 1;
71: MatSetValues(P, 1, &i, 3, col, value, INSERT_VALUES);
72: }
73: i = n - 1;
74: col[0] = n - 2;
75: col[1] = n - 1;
76: MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES);
77: i = 0;
78: col[0] = 0;
79: col[1] = 1;
80: value[0] = 2.0;
81: value[1] = -1.0;
82: MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES);
83: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
86: /* Set exact solution */
87: VecSet(u, one);
89: /* Create a matrix-free matrix As, P is used as a data context in MyMatShellMult() */
90: MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, P, &As);
91: MatSetFromOptions(As);
92: MatShellSetOperation(As, MATOP_MULT, (void (*)(void))MyMatShellMult);
94: /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
95: MatIsLinear(As, 10, &flg);
98: /* Compute right-hand-side vector. */
99: MatMult(As, u, b);
101: MatSetOption(As, MAT_SYMMETRIC, PETSC_TRUE);
102: MatMultTranspose(As, u, x);
103: VecAXPY(x, -1.0, b);
104: VecNorm(x, NORM_INFINITY, &norm);
106: MatSetOption(As, MAT_HERMITIAN, PETSC_TRUE);
107: MatMultHermitianTranspose(As, u, x);
108: VecAXPY(x, -1.0, b);
109: VecNorm(x, NORM_INFINITY, &norm);
112: /* Create the linear solver and set various options */
113: KSPCreate(PETSC_COMM_WORLD, &ksp);
114: KSPSetOperators(ksp, As, P);
116: /* Set linear solver defaults for this problem (optional). */
117: KSPGetPC(ksp, &pc);
118: PCSetType(pc, PCNONE);
119: KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
121: /* Set runtime options */
122: KSPSetFromOptions(ksp);
124: /* Solve linear system */
125: KSPSolve(ksp, b, x);
127: /* Check the error */
128: VecAXPY(x, -1.0, u);
129: VecNorm(x, NORM_2, &norm);
130: KSPGetIterationNumber(ksp, &its);
131: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
133: /* Free work space. */
134: VecDestroy(&x);
135: VecDestroy(&u);
136: VecDestroy(&b);
137: MatDestroy(&P);
138: MatDestroy(&As);
139: KSPDestroy(&ksp);
141: PetscFinalize();
142: return 0;
143: }
145: /*TEST
147: test:
148: args: -ksp_monitor_short -ksp_max_it 10
149: test:
150: suffix: 2
151: args: -ksp_monitor_short -ksp_max_it 10
153: TEST*/