1 /*******************************************************
 2  *  This file is part of the FreeFEM tutorial made by  *
 3  *      Pierre Jolivet <pierre@joliv.et>               *
 4  *                                                     *
 5  *  See https://joliv.et/FreeFem-tutorial for more     *
 6  *                                                     *
 7  *   Description: matrix-free solvers                  *
 8  *******************************************************/
  
10 load "msh3"
11 load "PETSc"
12 include "getARGV.idp"
13 int n = max(10, getARGV("-n", 10));
14 meshL Th = Sline(n, [x, y, z]);
 
16 func real[int] op(real[int]& in) {
17     real[int] out(in.n);
18     out = 0;
19     for(int i = 0; i < in.n; ++i) {
20         if(i != 0 && i != in.n - 1) {
21             out[i] = (2 * in[i] - in[i - 1] - in[i + 1])*n^2;
22         }
23         else if(i == 0)
24             out[0] = (2 * in[0] - in[1])*n^2;
25         else
26             out[in.n - 1] = (2 * in[in.n - 1] - in[in.n - 2])*n^2;
27     }
28     return out;
29 }
30 func real[int] prec(real[int]& in) {
31     real[int, int] A(n + 1, n + 1);
32     A = 0;
33     for(int i = 0; i < n + 1; ++i)
34         A(i, i) = 2 * n^2;
35     for(int i = 1; i < n + 1; ++i) {
36         A(i - 1, i) = -n^2;
37         A(i, i - 1) = -n^2;
38     }
39     matrix sparse = A;
40     set(sparse, solver = LU);
41     real[int] out(in.n);
42     out = sparse^-1 * in;
43     return out;
44 }
 
46 fespace Vh(Th, P1);
47 Vh sol;
48 real[int] rhs(n + 1);
49 rhs = 1;
50 KSPSolve(op, rhs, sol[], precon = prec, sparams = "-ksp_type gmres -ksp_view_singularvalues -ksp_converged_reason -ksp_gmres_restart 1000 -pc_type none");
51 plot(sol, fill = 1, value = 1);