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: PetscFunctionBegin;
24: /* printf("MatShellMult...user should implement this routine without using a matrix\n"); */
25: PetscCall(MatShellGetContext(As, &P));
26: PetscCall(MatMult(P, x, y));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: int main(int argc, char **args)
31: {
32: Vec x, b, u; /* approx solution, RHS, exact solution */
33: Mat P, As; /* preconditioner matrix, linear system (matrix-free) */
34: KSP ksp; /* linear solver context */
35: PC pc; /* preconditioner context */
36: PetscReal norm; /* norm of solution error */
37: PetscInt i, n = 100, col[3], its;
38: PetscMPIInt size;
39: PetscScalar one = 1.0, value[3];
40: PetscBool flg;
42: PetscFunctionBeginUser;
43: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
44: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
45: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Compute the matrix and right-hand-side vector that define
49: the linear system, As x = b.
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: /* Create vectors */
52: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
53: PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
54: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
55: PetscCall(VecSetFromOptions(x));
56: PetscCall(VecDuplicate(x, &b));
57: PetscCall(VecDuplicate(x, &u));
59: /* Create matrix P, to be used as preconditioner */
60: PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
61: PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, n, n));
62: PetscCall(MatSetFromOptions(P));
63: PetscCall(MatSetUp(P));
65: value[0] = -1.0;
66: value[1] = 2.0;
67: value[2] = -1.0;
68: for (i = 1; i < n - 1; i++) {
69: col[0] = i - 1;
70: col[1] = i;
71: col[2] = i + 1;
72: PetscCall(MatSetValues(P, 1, &i, 3, col, value, INSERT_VALUES));
73: }
74: i = n - 1;
75: col[0] = n - 2;
76: col[1] = n - 1;
77: PetscCall(MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES));
78: i = 0;
79: col[0] = 0;
80: col[1] = 1;
81: value[0] = 2.0;
82: value[1] = -1.0;
83: PetscCall(MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES));
84: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
87: /* Set exact solution */
88: PetscCall(VecSet(u, one));
90: /* Create a matrix-free matrix As, P is used as a data context in MyMatShellMult() */
91: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, P, &As));
92: PetscCall(MatSetFromOptions(As));
93: PetscCall(MatShellSetOperation(As, MATOP_MULT, (void (*)(void))MyMatShellMult));
95: /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
96: PetscCall(MatIsLinear(As, 10, &flg));
97: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Shell matrix As is non-linear! Use '-info |grep MatIsLinear' to get detailed report");
99: /* Compute right-hand-side vector. */
100: PetscCall(MatMult(As, u, b));
102: PetscCall(MatSetOption(As, MAT_SYMMETRIC, PETSC_TRUE));
103: PetscCall(MatMultTranspose(As, u, x));
104: PetscCall(VecAXPY(x, -1.0, b));
105: PetscCall(VecNorm(x, NORM_INFINITY, &norm));
106: PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)As), PETSC_ERR_PLIB, "Error ||A x-A^T x||_\\infty: %1.6e", (double)norm);
107: PetscCall(MatSetOption(As, MAT_HERMITIAN, PETSC_TRUE));
108: PetscCall(MatMultHermitianTranspose(As, u, x));
109: PetscCall(VecAXPY(x, -1.0, b));
110: PetscCall(VecNorm(x, NORM_INFINITY, &norm));
111: PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)As), PETSC_ERR_PLIB, "Error ||A x-A^H x||_\\infty: %1.6e", (double)norm);
113: /* Create the linear solver and set various options */
114: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
115: PetscCall(KSPSetOperators(ksp, As, P));
117: /* Set linear solver defaults for this problem (optional). */
118: PetscCall(KSPGetPC(ksp, &pc));
119: PetscCall(PCSetType(pc, PCNONE));
120: PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
122: /* Set runtime options */
123: PetscCall(KSPSetFromOptions(ksp));
125: /* Solve linear system */
126: PetscCall(KSPSolve(ksp, b, x));
128: /* Check the error */
129: PetscCall(VecAXPY(x, -1.0, u));
130: PetscCall(VecNorm(x, NORM_2, &norm));
131: PetscCall(KSPGetIterationNumber(ksp, &its));
132: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
134: /* Free work space. */
135: PetscCall(VecDestroy(&x));
136: PetscCall(VecDestroy(&u));
137: PetscCall(VecDestroy(&b));
138: PetscCall(MatDestroy(&P));
139: PetscCall(MatDestroy(&As));
140: PetscCall(KSPDestroy(&ksp));
142: PetscCall(PetscFinalize());
143: return 0;
144: }
146: /*TEST
148: test:
149: args: -ksp_monitor_short -ksp_max_it 10
150: test:
151: suffix: 2
152: args: -ksp_monitor_short -ksp_max_it 10
154: TEST*/