Actual source code: ex7.c
1: static char help[] = "Illustrate how to solves a matrix-free linear system with KSP.\n\n";
3: /*
4: Note: modified from ~src/ksp/ksp/tutorials/ex1.c
5: */
6: #include <petscksp.h>
8: /*
9: MatShellMult - Computes the matrix-vector product, y = As x.
11: Input Parameters:
12: As - the matrix-free matrix
13: x - vector
15: Output Parameter:
16: y - vector
17: */
18: PetscErrorCode MyMatShellMult(Mat As, Vec x, Vec y)
19: {
20: Mat P;
22: PetscFunctionBegin;
23: /* printf("MatShellMult...user should implement this routine without using a matrix\n"); */
24: PetscCall(MatShellGetContext(As, &P));
25: PetscCall(MatMult(P, x, y));
26: PetscFunctionReturn(PETSC_SUCCESS);
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;
41: PetscFunctionBeginUser;
42: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
43: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
44: PetscCall(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: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
52: PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
53: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
54: PetscCall(VecSetFromOptions(x));
55: PetscCall(VecDuplicate(x, &b));
56: PetscCall(VecDuplicate(x, &u));
58: /* Create matrix P, to be used as preconditioner */
59: PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
60: PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, n, n));
61: PetscCall(MatSetFromOptions(P));
62: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES));
83: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
84: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
86: /* Set exact solution */
87: PetscCall(VecSet(u, one));
89: /* Create a matrix-free matrix As, P is used as a data context in MyMatShellMult() */
90: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, P, &As));
91: PetscCall(MatSetFromOptions(As));
92: PetscCall(MatShellSetOperation(As, MATOP_MULT, (void (*)(void))MyMatShellMult));
94: /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
95: PetscCall(MatIsLinear(As, 10, &flg));
96: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Shell matrix As is non-linear! Use '-info |grep MatIsLinear' to get detailed report");
98: /* Compute right-hand-side vector. */
99: PetscCall(MatMult(As, u, b));
101: PetscCall(MatSetOption(As, MAT_SYMMETRIC, PETSC_TRUE));
102: PetscCall(MatMultTranspose(As, u, x));
103: PetscCall(VecAXPY(x, -1.0, b));
104: PetscCall(VecNorm(x, NORM_INFINITY, &norm));
105: PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)As), PETSC_ERR_PLIB, "Error ||A x-A^T x||_\\infty: %1.6e", (double)norm);
106: PetscCall(MatSetOption(As, MAT_HERMITIAN, PETSC_TRUE));
107: PetscCall(MatMultHermitianTranspose(As, u, x));
108: PetscCall(VecAXPY(x, -1.0, b));
109: PetscCall(VecNorm(x, NORM_INFINITY, &norm));
110: PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)As), PETSC_ERR_PLIB, "Error ||A x-A^H x||_\\infty: %1.6e", (double)norm);
112: /* Create the linear solver and set various options */
113: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
114: PetscCall(KSPSetOperators(ksp, As, P));
116: /* Set linear solver defaults for this problem (optional). */
117: PetscCall(KSPGetPC(ksp, &pc));
118: PetscCall(PCSetType(pc, PCNONE));
119: PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
121: /* Set runtime options */
122: PetscCall(KSPSetFromOptions(ksp));
124: /* Solve linear system */
125: PetscCall(KSPSolve(ksp, b, x));
127: /* Check the error */
128: PetscCall(VecAXPY(x, -1.0, u));
129: PetscCall(VecNorm(x, NORM_2, &norm));
130: PetscCall(KSPGetIterationNumber(ksp, &its));
131: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
133: /* Free work space. */
134: PetscCall(VecDestroy(&x));
135: PetscCall(VecDestroy(&u));
136: PetscCall(VecDestroy(&b));
137: PetscCall(MatDestroy(&P));
138: PetscCall(MatDestroy(&As));
139: PetscCall(KSPDestroy(&ksp));
141: PetscCall(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*/