1 static char help[] = "KSP linear solver on an operator with a null space.\n\n";
2
3 #include <petscksp.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Vec x, b, u; /* approx solution, RHS, exact solution */
8 Mat A; /* linear system matrix */
9 KSP ksp; /* KSP context */
10 PetscInt i, n = 10, col[3], its, i1, i2;
11 PetscScalar none = -1.0, value[3], avalue;
12 PetscReal norm;
13 PC pc;
14
15 PetscFunctionBeginUser;
16 PetscCall(PetscInitialize(&argc, &args, NULL, help));
17 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
18
19 /* Create vectors */
20 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
21 PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
22 PetscCall(VecSetFromOptions(x));
23 PetscCall(VecDuplicate(x, &b));
24 PetscCall(VecDuplicate(x, &u));
25
26 /* create a solution that is orthogonal to the constants */
27 PetscCall(VecGetOwnershipRange(u, &i1, &i2));
28 for (i = i1; i < i2; i++) {
29 avalue = i;
30 VecSetValues(u, 1, &i, &avalue, INSERT_VALUES);
31 }
32 PetscCall(VecAssemblyBegin(u));
33 PetscCall(VecAssemblyEnd(u));
34 PetscCall(VecSum(u, &avalue));
35 avalue = -avalue / (PetscReal)n;
36 PetscCall(VecShift(u, avalue));
37
38 /* Create and assemble matrix */
39 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
40 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
41 PetscCall(MatSetFromOptions(A));
42 value[0] = -1.0;
43 value[1] = 2.0;
44 value[2] = -1.0;
45 for (i = 1; i < n - 1; i++) {
46 col[0] = i - 1;
47 col[1] = i;
48 col[2] = i + 1;
49 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
50 }
51 i = n - 1;
52 col[0] = n - 2;
53 col[1] = n - 1;
54 value[1] = 1.0;
55 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
56 i = 0;
57 col[0] = 0;
58 col[1] = 1;
59 value[0] = 1.0;
60 value[1] = -1.0;
61 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
62 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
63 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
64 PetscCall(MatMult(A, u, b));
65
66 /* Create KSP context; set operators and options; solve linear system */
67 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
68 PetscCall(KSPSetOperators(ksp, A, A));
69
70 /* Insure that preconditioner has same null space as matrix */
71 /* currently does not do anything */
72 PetscCall(KSPGetPC(ksp, &pc));
73
74 PetscCall(KSPSetFromOptions(ksp));
75 PetscCall(KSPSolve(ksp, b, x));
76 /* PetscCall(KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD)); */
77
78 /* Check error */
79 PetscCall(VecAXPY(x, none, u));
80 PetscCall(VecNorm(x, NORM_2, &norm));
81 PetscCall(KSPGetIterationNumber(ksp, &its));
82 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
83
84 /* Free work space */
85 PetscCall(VecDestroy(&x));
86 PetscCall(VecDestroy(&u));
87 PetscCall(VecDestroy(&b));
88 PetscCall(MatDestroy(&A));
89 PetscCall(KSPDestroy(&ksp));
90 PetscCall(PetscFinalize());
91 return 0;
92 }
93