SIAMFANLCh3s.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\calf}{{\\cal F}}\n",
    "\\newcommand{\\dnu}{d \\nu}\n",
    "\\newcommand{\\mf}{{\\bf F}}\n",
    "\\newcommand{\\md}{{\\bf D}}\n",
    "\\newcommand{\\mP}{{\\bf P}}\n",
    "\\newcommand{\\mU}{{\\bf U}}\n",
    "\\newcommand{\\vu}{{\\bf u}}\n",
    "\\newcommand{\\vx}{{\\bf x}}\n",
    "\\newcommand{\\vw}{{\\bf w}}\n",
    "\\newcommand{\\vy}{{\\bf y}}\n",
    "\\newcommand{\\vf}{{\\bf f}}\n",
    "\\newcommand{\\vs}{{\\bf s}}\n",
    "\\newcommand{\\ve}{{\\bf e}}\n",
    "\\newcommand{\\vd}{{\\bf d}}\n",
    "\\newcommand{\\vb}{{\\bf b}}\n",
    "\\newcommand{\\vz}{{\\bf z}}\n",
    "\\newcommand{\\mg}{{\\bf G}}\n",
    "\\newcommand{\\ml}{{\\bf L}}\n",
    "\\newcommand{\\mg}{{\\bf G}}\n",
    "\\newcommand{\\mv}{{\\bf V}}\n",
    "\\newcommand{\\ma}{{\\bf A}}\n",
    "\\newcommand{\\mi}{{\\bf I}}\n",
    "\\newcommand{\\mm}{{\\bf M}}\n",
    "\\newcommand{\\mb}{{\\bf B}}\n",
    "\\newcommand{\\ball}{{\\cal B}}\n",
    "\\newcommand{\\ptc}{{\\Psi TC}}\n",
    "\\newcommand{\\diag}{\\mbox{diag}}\n",
    "\\newcommand{\\begeq}{{\\begin{equation}}}\n",
    "\\newcommand{\\endeq}{{\\end{equation}}}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"fanote_init.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3.8: Solvers for Chapter 3\n",
    "\n",
    "Contents for Section 3.8\n",
    "\n",
    "[Overview](#Overview)\n",
    "\n",
    "[nsoli.jl](#Section-3.8.1:-nsoli.jl)\n",
    "\n",
    "- [Benchmarking the H-equation with nsoli.jl](#Section-3.8.2:-Benchmarking-the-H-equation-with-nsoli.jl)\n",
    "\n",
    "- [ Preconditioning the Convection-Diffusion Equation](#Section-3.8.3:-Preconditioning-the-Convection-Diffusion-Equation)\n",
    "\n",
    "- [ptcsoli.jl](#Section-3.8.4:-ptcsoli.jl)\n",
    "\n",
    "- [Benchmarking $\\ptc$ with the buckling beam problem](#Section-3.8.5:-Benchmarking-$\\ptc$-with-the-buckling-beam-problem)\n",
    "\n",
    "[Projects](#Section-3.9:-Projects)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will follow the pattern of the previous chapters and present two solvers, a Newton code and a $\\ptc$ code. Both codes are for systems of equations and use Krylov methods to compute the step. We have two Krylov solvers, GMRES and BiCGstab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 3.8.1: nsoli.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__nsoli.jl__ solves systems of nonlinear equations with Newton-Krylov methods. As usual, we begin with the docstrings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "search: \u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1mi\u001b[22m \u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22m \u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22msc \u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22mheq ptcsoli rsplit \u001b[0m\u001b[1mN\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1mi\u001b[22mPDE ishold\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "nsoli(F!, x0, FS, FPS, Jvec=dirder; rtol=1.e-6, atol=1.e-12,\n",
       "           maxit=20, lmaxit=-1, lsolver=\"gmres\", eta=.1,\n",
       "           fixedeta=true, Pvec=nothing, pside=\"right\",\n",
       "           armmax=10, dx = 1.e-7, armfix=false, pdata = nothing,\n",
       "           printerr = true, keepsolhist = false, \n",
       "           Krylov_Data = nothing, stagnationok=false)\n",
       "\\end{verbatim}\n",
       ")\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  Herewith: nsoli\n",
       "\n",
       "You must allocate storage for the function and the Krylov basis in advance –> in the calling program <– ie. in FS and FPS\n",
       "\n",
       "Inputs:\n",
       "\n",
       "\\begin{itemize}\n",
       "\\item F!: function evaluation, the ! indicates that F! overwrites FS, your   preallocated storage for the function.\n",
       "\n",
       "So FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)\n",
       "\n",
       "Your function MUST have –> return FS <– at the end.   See the examples in the docstrings\n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\n",
       "\\item x0: initial iterate\n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\n",
       "\\item FS: Preallocated storage for function. It is a vector of size N\n",
       "\n",
       "You should store it as (N) and design F! to use vectors of size (N). If you use (N,1) consistently instead, the solvers may work, but I make no guarantees.\n",
       "\n",
       "\n",
       "\\item FPS: preallocated storage for the Krylov basis. It is an N x m matrix where      you plan to take at most m-1 GMRES iterations before a restart. \n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\n",
       "\\item Jvec: Jacobian vector product, If you leave this out the   default is a finite difference directional derivative.\n",
       "\n",
       "So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v. \n",
       "\n",
       "(v, FS, x) or (v, FS, x, pdata) must be the argument list,    even if FP does not need FS.   One reason for this is that the finite-difference derivative   does and that is the default in the solver.\n",
       "\n",
       "\n",
       "\\item Precision: Lemme tell ya 'bout precision. I designed this code for    full precision functions and linear algebra in any precision you want.    You can declare FPS as Float64 or Float32 and nsoli    will do the right thing. Float16 support is there, but not working well.\n",
       "\n",
       "If the Jacobian is reasonably well conditioned, you can cut the cost   of orthogonalization and storage (for GMRES) in half with no loss.    There is no benefit if your linear solver is not GMRES or if    orthogonalization and storage of the Krylov vectors is only a   small part of the cost of the computation. So if your preconditioner   is good and you only need a few Krylovs/Newton, reduced precision won't   help you much.\n",
       "\n",
       "BiCGSTAB does not benefit from reduced precision. \n",
       "\n",
       "\\end{itemize}\n",
       "\\rule{\\textwidth}{1pt}\n",
       "Keyword Arguments (kwargs):\n",
       "\n",
       "rtol and atol: relative and absolute error tolerances\n",
       "\n",
       "maxit: limit on nonlinear iterations\n",
       "\n",
       "lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has m columns, and you need more than m-1 linear iterations, then GMRES  will restart. \n",
       "\n",
       "The default is -1 for GMRES. This means that you'll take m-1 iterations,  where size(V) = (n,m), and get no restarts. For BiCGSTAB the default is 10.\n",
       "\n",
       "lsolver: the linear solver, default = \"gmres\"\n",
       "\n",
       "Your choices will be \"gmres\" or \"bicgstab\". However, gmres is the only option for now.\n",
       "\n",
       "eta and fixed eta: eta > 0 or there's an error\n",
       "\n",
       "The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||\n",
       "\n",
       "where \n",
       "\n",
       "etag = eta if fixedeta=true\n",
       "\n",
       "etag = Eisenstat-Walker as implemented in book if fixedeta=false\n",
       "\n",
       "The default, which may change, is eta=.1, fixedeta=true\n",
       "\n",
       "Pvec: Preconditioner-vector product. The rules are similar to Jvec     So, Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where     P(x) is the preconditioner. You must use x as an input even     if your preconditioner does not depend on x\n",
       "\n",
       "pside: apply preconditioner on pside, default = \"right\". I do not       recommend \"left\". See Chapter 3 for the story on this.\n",
       "\n",
       "armmax: upper bound on step size reductions in line search\n",
       "\n",
       "dx: default = 1.e-7\n",
       "\n",
       "difference increment in finite-difference derivatives       h=dx*norm(x,Inf)+1.e-8\n",
       "\n",
       "armfix: default = false\n",
       "\n",
       "The default is a parabolic line search (ie false). Set to true and the step size will be fixed at .5. Don't do this unless you are doing experiments for research.\n",
       "\n",
       "pdata:\n",
       "\n",
       "precomputed data for the function, Jacobian-vector, and Preconditioner-vector products.  Things will go better if you use this rather than hide the data  in global variables within the module for your function/Jacobian\n",
       "\n",
       "If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of them.\n",
       "\n",
       "printerr: default = true\n",
       "\n",
       "I print a helpful message when the solver fails. To suppress that message set printerr to false.\n",
       "\n",
       "keepsolhist: default = false\n",
       "\n",
       "Set this to true to get the history of the iteration in the output tuple. This is on by default for scalar equations and off for systems. Only turn it on if you have use for the data, which can get REALLY LARGE.\n",
       "\n",
       "Krylov\\_Data: default = nothing\n",
       "\n",
       "This is a structure where I put the internal storage for the solvers. You can (but probably should not) preallocate this your self with the nkl\\_init function.\n",
       "\n",
       "Krylov\\_Data \n",
       "\n",
       "= nkl\\_init(n,lsolver)\n",
       "\n",
       "This is a dangerous thing to mess with and I only recommend it if  the allocations in nsoli become a problem in continuation or IVP integration. Krylov\\_Data is where I store the solution at the end of the iteration and if you reuse it without copying the solution to somewhere else, you'll lose it and it will be overwritten with the new solution. The continuation case study uses this and you  should look at that to see what I did.\n",
       "\n",
       "stagnationok: default = false\n",
       "\n",
       "Set this to true if you want to disable the line search and either observe divergence or stagnation. This is only useful for research or writing a book.\n",
       "\n",
       "Output:\n",
       "\n",
       "\\begin{itemize}\n",
       "\\item A named tuple (solution, functionval, history, stats, idid,              errcode, solhist)\n",
       "\n",
       "\\end{itemize}\n",
       "where\n",
       "\n",
       "– solution = converged result\n",
       "\n",
       "– functionval = F(solution)\n",
       "\n",
       "– history = the vector of residual norms (||F(x)||) for the iteration\n",
       "\n",
       "– stats = named tuple of the history of (ifun, ijac, iarm, ikfail), the  number of functions/Jacobian-vector prods/steplength reductions/linear solver failures at each iteration. Linear solver failures DO NOT mean that the nonlinear solver will fail. You should look at this stat if, for example, the line search fails. Increasing the size of FPS and/or lmaxit might solve the problem.\n",
       "\n",
       "I do not count the function values for a finite-difference derivative because they count toward a Jacobian-vector product.\n",
       "\n",
       "– idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "– errcode = 0 if the iteration succeeded\n",
       "\n",
       "\\begin{verbatim}\n",
       "    = -1 if the initial iterate satisfies the termination criteria\n",
       "\n",
       "    = 10 if no convergence after maxit iterations\n",
       "\n",
       "    = 1  if the line search failed\n",
       "\\end{verbatim}\n",
       "– solhist:\n",
       "\n",
       "\\begin{verbatim}\n",
       "  This is the entire history of the iteration if you've set\n",
       "  keepsolhist=true\n",
       "\\end{verbatim}\n",
       "solhist is an N x K array where N is the length of x and K is the number of iteration + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "\\rule{\\textwidth}{1pt}\n",
       "\\subsubsection{Example for nsoli}\n",
       "\\paragraph{Simple 2D problem.}\n",
       "You should get the same results as for nsol.jl because GMRES will solve the equation for the step exactly in two iterations. Finite difference Jacobians and analytic Jacobian-vector products for full precision and finite difference Jacobian-vector products for single precision.\n",
       "\n",
       "BiCGSTAB converges in 5 iterations and each nonlinear iteration costs two Jacobian-vector products. Note that the storage for the Krylov space in GMRES (jvs) is replace by a single vector (fpv) when BiCGSTAB is the linear solver.\n",
       "\n",
       "\\begin{verbatim}\n",
       "julia> function f!(fv,x)\n",
       "       fv[1]=x[1] + sin(x[2])\n",
       "       fv[2]=cos(x[1]+x[2])\n",
       "       return fv\n",
       "       end\n",
       "f! (generic function with 1 method)\n",
       "\n",
       "julia> function JVec(v, fv, x)\n",
       "       jvec=zeros(2);\n",
       "       p=-sin(x[1]+x[2])\n",
       "       jvec[1]=v[1]+cos(x[2])*v[2]\n",
       "       jvec[2]=p*(v[1]+v[2])\n",
       "       return jvec\n",
       "       end\n",
       "JVec (generic function with 1 method)\n",
       "\n",
       "julia> x0=ones(2); fv=zeros(2); jv=zeros(2,2); \n",
       "\n",
       "julia> jv32=zeros(Float32,2,2);\n",
       "\n",
       "julia> jvs=zeros(2,3); jvs32=zeros(Float32,2,3);\n",
       "\n",
       "julia> nout=nsol(f!,x0,fv,jv; sham=1);\n",
       "\n",
       "julia> kout=nsoli(f!,x0,fv,jvs,JVec; \n",
       "                  fixedeta=true, eta=.1, lmaxit=2);\n",
       "\n",
       "julia> kout32=nsoli(f!,x0,fv,jvs32; \n",
       "                    fixedeta=true, eta=.1, lmaxit=2);\n",
       "\n",
       "julia> [nout.history kout.history kout32.history]\n",
       "5×3 Array{Float64,2}:\n",
       " 1.88791e+00  1.88791e+00  1.88791e+00\n",
       " 2.43119e-01  2.43120e-01  2.43119e-01\n",
       " 1.19231e-02  1.19231e-02  1.19230e-02\n",
       " 1.03266e-05  1.03261e-05  1.03264e-05\n",
       " 1.46388e-11  1.40862e-11  1.39825e-11\n",
       "\n",
       "julia> fpv=zeros(2);\n",
       "\n",
       "julia> koutb=nsoli(f!,x0,fv,fpv,JVec; \n",
       "            fixedeta=true, eta=.1, lmaxit=2, lsolver=\"bicgstab\");\n",
       "\n",
       "julia> koutb.history\n",
       "6-element Vector{Float64}:\n",
       " 1.88791e+00\n",
       " 2.43120e-01\n",
       " 1.19231e-02\n",
       " 4.87500e-04\n",
       " 7.54236e-06\n",
       " 3.84646e-07\n",
       "\\end{verbatim}\n"
      ],
      "text/markdown": [
       "```\n",
       "nsoli(F!, x0, FS, FPS, Jvec=dirder; rtol=1.e-6, atol=1.e-12,\n",
       "           maxit=20, lmaxit=-1, lsolver=\"gmres\", eta=.1,\n",
       "           fixedeta=true, Pvec=nothing, pside=\"right\",\n",
       "           armmax=10, dx = 1.e-7, armfix=false, pdata = nothing,\n",
       "           printerr = true, keepsolhist = false, \n",
       "           Krylov_Data = nothing, stagnationok=false)\n",
       "```\n",
       "\n",
       ")\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  Herewith: nsoli\n",
       "\n",
       "You must allocate storage for the function and the Krylov basis in advance –> in the calling program <– ie. in FS and FPS\n",
       "\n",
       "Inputs:\n",
       "\n",
       "  * F!: function evaluation, the ! indicates that F! overwrites FS, your   preallocated storage for the function.\n",
       "\n",
       "    So FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)\n",
       "\n",
       "    Your function MUST have –> return FS <– at the end.   See the examples in the docstrings\n",
       "\n",
       "  * x0: initial iterate\n",
       "\n",
       "  * FS: Preallocated storage for function. It is a vector of size N\n",
       "\n",
       "    You should store it as (N) and design F! to use vectors of size (N). If you use (N,1) consistently instead, the solvers may work, but I make no guarantees.\n",
       "  * FPS: preallocated storage for the Krylov basis. It is an N x m matrix where      you plan to take at most m-1 GMRES iterations before a restart.\n",
       "\n",
       "  * Jvec: Jacobian vector product, If you leave this out the   default is a finite difference directional derivative.\n",
       "\n",
       "    So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v. \n",
       "\n",
       "    (v, FS, x) or (v, FS, x, pdata) must be the argument list,    even if FP does not need FS.   One reason for this is that the finite-difference derivative   does and that is the default in the solver.\n",
       "  * Precision: Lemme tell ya 'bout precision. I designed this code for    full precision functions and linear algebra in any precision you want.    You can declare FPS as Float64 or Float32 and nsoli    will do the right thing. Float16 support is there, but not working well.\n",
       "\n",
       "    If the Jacobian is reasonably well conditioned, you can cut the cost   of orthogonalization and storage (for GMRES) in half with no loss.    There is no benefit if your linear solver is not GMRES or if    orthogonalization and storage of the Krylov vectors is only a   small part of the cost of the computation. So if your preconditioner   is good and you only need a few Krylovs/Newton, reduced precision won't   help you much.\n",
       "\n",
       "    BiCGSTAB does not benefit from reduced precision.\n",
       "\n",
       "---\n",
       "\n",
       "Keyword Arguments (kwargs):\n",
       "\n",
       "rtol and atol: relative and absolute error tolerances\n",
       "\n",
       "maxit: limit on nonlinear iterations\n",
       "\n",
       "lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has m columns, and you need more than m-1 linear iterations, then GMRES  will restart. \n",
       "\n",
       "The default is -1 for GMRES. This means that you'll take m-1 iterations,  where size(V) = (n,m), and get no restarts. For BiCGSTAB the default is 10.\n",
       "\n",
       "lsolver: the linear solver, default = \"gmres\"\n",
       "\n",
       "Your choices will be \"gmres\" or \"bicgstab\". However, gmres is the only option for now.\n",
       "\n",
       "eta and fixed eta: eta > 0 or there's an error\n",
       "\n",
       "The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||\n",
       "\n",
       "where \n",
       "\n",
       "etag = eta if fixedeta=true\n",
       "\n",
       "etag = Eisenstat-Walker as implemented in book if fixedeta=false\n",
       "\n",
       "The default, which may change, is eta=.1, fixedeta=true\n",
       "\n",
       "Pvec: Preconditioner-vector product. The rules are similar to Jvec     So, Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where     P(x) is the preconditioner. You must use x as an input even     if your preconditioner does not depend on x\n",
       "\n",
       "pside: apply preconditioner on pside, default = \"right\". I do not       recommend \"left\". See Chapter 3 for the story on this.\n",
       "\n",
       "armmax: upper bound on step size reductions in line search\n",
       "\n",
       "dx: default = 1.e-7\n",
       "\n",
       "difference increment in finite-difference derivatives       h=dx*norm(x,Inf)+1.e-8\n",
       "\n",
       "armfix: default = false\n",
       "\n",
       "The default is a parabolic line search (ie false). Set to true and the step size will be fixed at .5. Don't do this unless you are doing experiments for research.\n",
       "\n",
       "pdata:\n",
       "\n",
       "precomputed data for the function, Jacobian-vector, and Preconditioner-vector products.  Things will go better if you use this rather than hide the data  in global variables within the module for your function/Jacobian\n",
       "\n",
       "If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of them.\n",
       "\n",
       "printerr: default = true\n",
       "\n",
       "I print a helpful message when the solver fails. To suppress that message set printerr to false.\n",
       "\n",
       "keepsolhist: default = false\n",
       "\n",
       "Set this to true to get the history of the iteration in the output tuple. This is on by default for scalar equations and off for systems. Only turn it on if you have use for the data, which can get REALLY LARGE.\n",
       "\n",
       "Krylov_Data: default = nothing\n",
       "\n",
       "This is a structure where I put the internal storage for the solvers. You can (but probably should not) preallocate this your self with the nkl_init function.\n",
       "\n",
       "Krylov_Data \n",
       "\n",
       "= nkl_init(n,lsolver)\n",
       "\n",
       "This is a dangerous thing to mess with and I only recommend it if  the allocations in nsoli become a problem in continuation or IVP integration. Krylov_Data is where I store the solution at the end of the iteration and if you reuse it without copying the solution to somewhere else, you'll lose it and it will be overwritten with the new solution. The continuation case study uses this and you  should look at that to see what I did.\n",
       "\n",
       "stagnationok: default = false\n",
       "\n",
       "Set this to true if you want to disable the line search and either observe divergence or stagnation. This is only useful for research or writing a book.\n",
       "\n",
       "Output:\n",
       "\n",
       "  * A named tuple (solution, functionval, history, stats, idid,              errcode, solhist)\n",
       "\n",
       "where\n",
       "\n",
       "– solution = converged result\n",
       "\n",
       "– functionval = F(solution)\n",
       "\n",
       "– history = the vector of residual norms (||F(x)||) for the iteration\n",
       "\n",
       "– stats = named tuple of the history of (ifun, ijac, iarm, ikfail), the  number of functions/Jacobian-vector prods/steplength reductions/linear solver failures at each iteration. Linear solver failures DO NOT mean that the nonlinear solver will fail. You should look at this stat if, for example, the line search fails. Increasing the size of FPS and/or lmaxit might solve the problem.\n",
       "\n",
       "I do not count the function values for a finite-difference derivative because they count toward a Jacobian-vector product.\n",
       "\n",
       "– idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "– errcode = 0 if the iteration succeeded\n",
       "\n",
       "```\n",
       "    = -1 if the initial iterate satisfies the termination criteria\n",
       "\n",
       "    = 10 if no convergence after maxit iterations\n",
       "\n",
       "    = 1  if the line search failed\n",
       "```\n",
       "\n",
       "– solhist:\n",
       "\n",
       "```\n",
       "  This is the entire history of the iteration if you've set\n",
       "  keepsolhist=true\n",
       "```\n",
       "\n",
       "solhist is an N x K array where N is the length of x and K is the number of iteration + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "---\n",
       "\n",
       "### Example for nsoli\n",
       "\n",
       "#### Simple 2D problem.\n",
       "\n",
       "You should get the same results as for nsol.jl because GMRES will solve the equation for the step exactly in two iterations. Finite difference Jacobians and analytic Jacobian-vector products for full precision and finite difference Jacobian-vector products for single precision.\n",
       "\n",
       "BiCGSTAB converges in 5 iterations and each nonlinear iteration costs two Jacobian-vector products. Note that the storage for the Krylov space in GMRES (jvs) is replace by a single vector (fpv) when BiCGSTAB is the linear solver.\n",
       "\n",
       "```jldoctest\n",
       "julia> function f!(fv,x)\n",
       "       fv[1]=x[1] + sin(x[2])\n",
       "       fv[2]=cos(x[1]+x[2])\n",
       "       return fv\n",
       "       end\n",
       "f! (generic function with 1 method)\n",
       "\n",
       "julia> function JVec(v, fv, x)\n",
       "       jvec=zeros(2);\n",
       "       p=-sin(x[1]+x[2])\n",
       "       jvec[1]=v[1]+cos(x[2])*v[2]\n",
       "       jvec[2]=p*(v[1]+v[2])\n",
       "       return jvec\n",
       "       end\n",
       "JVec (generic function with 1 method)\n",
       "\n",
       "julia> x0=ones(2); fv=zeros(2); jv=zeros(2,2); \n",
       "\n",
       "julia> jv32=zeros(Float32,2,2);\n",
       "\n",
       "julia> jvs=zeros(2,3); jvs32=zeros(Float32,2,3);\n",
       "\n",
       "julia> nout=nsol(f!,x0,fv,jv; sham=1);\n",
       "\n",
       "julia> kout=nsoli(f!,x0,fv,jvs,JVec; \n",
       "                  fixedeta=true, eta=.1, lmaxit=2);\n",
       "\n",
       "julia> kout32=nsoli(f!,x0,fv,jvs32; \n",
       "                    fixedeta=true, eta=.1, lmaxit=2);\n",
       "\n",
       "julia> [nout.history kout.history kout32.history]\n",
       "5×3 Array{Float64,2}:\n",
       " 1.88791e+00  1.88791e+00  1.88791e+00\n",
       " 2.43119e-01  2.43120e-01  2.43119e-01\n",
       " 1.19231e-02  1.19231e-02  1.19230e-02\n",
       " 1.03266e-05  1.03261e-05  1.03264e-05\n",
       " 1.46388e-11  1.40862e-11  1.39825e-11\n",
       "\n",
       "julia> fpv=zeros(2);\n",
       "\n",
       "julia> koutb=nsoli(f!,x0,fv,fpv,JVec; \n",
       "            fixedeta=true, eta=.1, lmaxit=2, lsolver=\"bicgstab\");\n",
       "\n",
       "julia> koutb.history\n",
       "6-element Vector{Float64}:\n",
       " 1.88791e+00\n",
       " 2.43120e-01\n",
       " 1.19231e-02\n",
       " 4.87500e-04\n",
       " 7.54236e-06\n",
       " 3.84646e-07\n",
       "```\n"
      ],
      "text/plain": [
       "\u001b[36m  nsoli(F!, x0, FS, FPS, Jvec=dirder; rtol=1.e-6, atol=1.e-12,\u001b[39m\n",
       "\u001b[36m             maxit=20, lmaxit=-1, lsolver=\"gmres\", eta=.1,\u001b[39m\n",
       "\u001b[36m             fixedeta=true, Pvec=nothing, pside=\"right\",\u001b[39m\n",
       "\u001b[36m             armmax=10, dx = 1.e-7, armfix=false, pdata = nothing,\u001b[39m\n",
       "\u001b[36m             printerr = true, keepsolhist = false, \u001b[39m\n",
       "\u001b[36m             Krylov_Data = nothing, stagnationok=false)\u001b[39m\n",
       "\n",
       "  )\n",
       "\n",
       "  C. T. Kelley, 2022\n",
       "\n",
       "  Julia versions of the nonlinear solvers from my SIAM books. Herewith: nsoli\n",
       "\n",
       "  You must allocate storage for the function and the Krylov basis in advance\n",
       "  –> in the calling program <– ie. in FS and FPS\n",
       "\n",
       "  Inputs:\n",
       "\n",
       "    •  F!: function evaluation, the ! indicates that F! overwrites FS,\n",
       "       your preallocated storage for the function.\n",
       "       So FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)\n",
       "       Your function MUST have –> return FS <– at the end. See the\n",
       "       examples in the docstrings\n",
       "\n",
       "    •  x0: initial iterate\n",
       "\n",
       "    •  FS: Preallocated storage for function. It is a vector of size N\n",
       "       You should store it as (N) and design F! to use vectors of size\n",
       "       (N). If you use (N,1) consistently instead, the solvers may work,\n",
       "       but I make no guarantees.\n",
       "\n",
       "    •  FPS: preallocated storage for the Krylov basis. It is an N x m\n",
       "       matrix where you plan to take at most m-1 GMRES iterations before\n",
       "       a restart.\n",
       "\n",
       "    •  Jvec: Jacobian vector product, If you leave this out the default\n",
       "       is a finite difference directional derivative.\n",
       "       So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v.\n",
       "       (v, FS, x) or (v, FS, x, pdata) must be the argument list, even if\n",
       "       FP does not need FS. One reason for this is that the\n",
       "       finite-difference derivative does and that is the default in the\n",
       "       solver.\n",
       "\n",
       "    •  Precision: Lemme tell ya 'bout precision. I designed this code for\n",
       "       full precision functions and linear algebra in any precision you\n",
       "       want. You can declare FPS as Float64 or Float32 and nsoli will do\n",
       "       the right thing. Float16 support is there, but not working well.\n",
       "       If the Jacobian is reasonably well conditioned, you can cut the\n",
       "       cost of orthogonalization and storage (for GMRES) in half with no\n",
       "       loss. There is no benefit if your linear solver is not GMRES or if\n",
       "       orthogonalization and storage of the Krylov vectors is only a\n",
       "       small part of the cost of the computation. So if your\n",
       "       preconditioner is good and you only need a few Krylovs/Newton,\n",
       "       reduced precision won't help you much.\n",
       "       BiCGSTAB does not benefit from reduced precision.\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "  Keyword Arguments (kwargs):\n",
       "\n",
       "  rtol and atol: relative and absolute error tolerances\n",
       "\n",
       "  maxit: limit on nonlinear iterations\n",
       "\n",
       "  lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has m\n",
       "  columns, and you need more than m-1 linear iterations, then GMRES will\n",
       "  restart.\n",
       "\n",
       "  The default is -1 for GMRES. This means that you'll take m-1 iterations,\n",
       "  where size(V) = (n,m), and get no restarts. For BiCGSTAB the default is 10.\n",
       "\n",
       "  lsolver: the linear solver, default = \"gmres\"\n",
       "\n",
       "  Your choices will be \"gmres\" or \"bicgstab\". However, gmres is the only\n",
       "  option for now.\n",
       "\n",
       "  eta and fixed eta: eta > 0 or there's an error\n",
       "\n",
       "  The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||\n",
       "\n",
       "  where\n",
       "\n",
       "  etag = eta if fixedeta=true\n",
       "\n",
       "  etag = Eisenstat-Walker as implemented in book if fixedeta=false\n",
       "\n",
       "  The default, which may change, is eta=.1, fixedeta=true\n",
       "\n",
       "  Pvec: Preconditioner-vector product. The rules are similar to Jvec So,\n",
       "  Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where P(x) is the\n",
       "  preconditioner. You must use x as an input even if your preconditioner does\n",
       "  not depend on x\n",
       "\n",
       "  pside: apply preconditioner on pside, default = \"right\". I do not recommend\n",
       "  \"left\". See Chapter 3 for the story on this.\n",
       "\n",
       "  armmax: upper bound on step size reductions in line search\n",
       "\n",
       "  dx: default = 1.e-7\n",
       "\n",
       "  difference increment in finite-difference derivatives h=dx*norm(x,Inf)+1.e-8\n",
       "\n",
       "  armfix: default = false\n",
       "\n",
       "  The default is a parabolic line search (ie false). Set to true and the step\n",
       "  size will be fixed at .5. Don't do this unless you are doing experiments for\n",
       "  research.\n",
       "\n",
       "  pdata:\n",
       "\n",
       "  precomputed data for the function, Jacobian-vector, and\n",
       "  Preconditioner-vector products. Things will go better if you use this rather\n",
       "  than hide the data in global variables within the module for your\n",
       "  function/Jacobian\n",
       "\n",
       "  If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of\n",
       "  them.\n",
       "\n",
       "  printerr: default = true\n",
       "\n",
       "  I print a helpful message when the solver fails. To suppress that message\n",
       "  set printerr to false.\n",
       "\n",
       "  keepsolhist: default = false\n",
       "\n",
       "  Set this to true to get the history of the iteration in the output tuple.\n",
       "  This is on by default for scalar equations and off for systems. Only turn it\n",
       "  on if you have use for the data, which can get REALLY LARGE.\n",
       "\n",
       "  Krylov_Data: default = nothing\n",
       "\n",
       "  This is a structure where I put the internal storage for the solvers. You\n",
       "  can (but probably should not) preallocate this your self with the nkl_init\n",
       "  function.\n",
       "\n",
       "  Krylov_Data\n",
       "\n",
       "  = nkl_init(n,lsolver)\n",
       "\n",
       "  This is a dangerous thing to mess with and I only recommend it if the\n",
       "  allocations in nsoli become a problem in continuation or IVP integration.\n",
       "  Krylov_Data is where I store the solution at the end of the iteration and if\n",
       "  you reuse it without copying the solution to somewhere else, you'll lose it\n",
       "  and it will be overwritten with the new solution. The continuation case\n",
       "  study uses this and you should look at that to see what I did.\n",
       "\n",
       "  stagnationok: default = false\n",
       "\n",
       "  Set this to true if you want to disable the line search and either observe\n",
       "  divergence or stagnation. This is only useful for research or writing a\n",
       "  book.\n",
       "\n",
       "  Output:\n",
       "\n",
       "    •  A named tuple (solution, functionval, history, stats, idid,\n",
       "       errcode, solhist)\n",
       "\n",
       "  where\n",
       "\n",
       "  – solution = converged result\n",
       "\n",
       "  – functionval = F(solution)\n",
       "\n",
       "  – history = the vector of residual norms (||F(x)||) for the iteration\n",
       "\n",
       "  – stats = named tuple of the history of (ifun, ijac, iarm, ikfail), the\n",
       "  number of functions/Jacobian-vector prods/steplength reductions/linear\n",
       "  solver failures at each iteration. Linear solver failures DO NOT mean that\n",
       "  the nonlinear solver will fail. You should look at this stat if, for\n",
       "  example, the line search fails. Increasing the size of FPS and/or lmaxit\n",
       "  might solve the problem.\n",
       "\n",
       "  I do not count the function values for a finite-difference derivative\n",
       "  because they count toward a Jacobian-vector product.\n",
       "\n",
       "  – idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "  – errcode = 0 if the iteration succeeded\n",
       "\n",
       "\u001b[36m      = -1 if the initial iterate satisfies the termination criteria\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m      = 10 if no convergence after maxit iterations\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m      = 1  if the line search failed\u001b[39m\n",
       "\n",
       "  – solhist:\n",
       "\n",
       "\u001b[36m    This is the entire history of the iteration if you've set\u001b[39m\n",
       "\u001b[36m    keepsolhist=true\u001b[39m\n",
       "\n",
       "  solhist is an N x K array where N is the length of x and K is the number of\n",
       "  iteration + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "\u001b[1m  Example for nsoli\u001b[22m\n",
       "\u001b[1m  –––––––––––––––––\u001b[22m\n",
       "\n",
       "\u001b[1m  Simple 2D problem.\u001b[22m\n",
       "\u001b[1m  ------------------\u001b[22m\n",
       "\n",
       "  You should get the same results as for nsol.jl because GMRES will solve the\n",
       "  equation for the step exactly in two iterations. Finite difference Jacobians\n",
       "  and analytic Jacobian-vector products for full precision and finite\n",
       "  difference Jacobian-vector products for single precision.\n",
       "\n",
       "  BiCGSTAB converges in 5 iterations and each nonlinear iteration costs two\n",
       "  Jacobian-vector products. Note that the storage for the Krylov space in\n",
       "  GMRES (jvs) is replace by a single vector (fpv) when BiCGSTAB is the linear\n",
       "  solver.\n",
       "\n",
       "\u001b[36m  julia> function f!(fv,x)\u001b[39m\n",
       "\u001b[36m         fv[1]=x[1] + sin(x[2])\u001b[39m\n",
       "\u001b[36m         fv[2]=cos(x[1]+x[2])\u001b[39m\n",
       "\u001b[36m         return fv\u001b[39m\n",
       "\u001b[36m         end\u001b[39m\n",
       "\u001b[36m  f! (generic function with 1 method)\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> function JVec(v, fv, x)\u001b[39m\n",
       "\u001b[36m         jvec=zeros(2);\u001b[39m\n",
       "\u001b[36m         p=-sin(x[1]+x[2])\u001b[39m\n",
       "\u001b[36m         jvec[1]=v[1]+cos(x[2])*v[2]\u001b[39m\n",
       "\u001b[36m         jvec[2]=p*(v[1]+v[2])\u001b[39m\n",
       "\u001b[36m         return jvec\u001b[39m\n",
       "\u001b[36m         end\u001b[39m\n",
       "\u001b[36m  JVec (generic function with 1 method)\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> x0=ones(2); fv=zeros(2); jv=zeros(2,2); \u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> jv32=zeros(Float32,2,2);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> jvs=zeros(2,3); jvs32=zeros(Float32,2,3);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> nout=nsol(f!,x0,fv,jv; sham=1);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> kout=nsoli(f!,x0,fv,jvs,JVec; \u001b[39m\n",
       "\u001b[36m                    fixedeta=true, eta=.1, lmaxit=2);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> kout32=nsoli(f!,x0,fv,jvs32; \u001b[39m\n",
       "\u001b[36m                      fixedeta=true, eta=.1, lmaxit=2);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> [nout.history kout.history kout32.history]\u001b[39m\n",
       "\u001b[36m  5×3 Array{Float64,2}:\u001b[39m\n",
       "\u001b[36m   1.88791e+00  1.88791e+00  1.88791e+00\u001b[39m\n",
       "\u001b[36m   2.43119e-01  2.43120e-01  2.43119e-01\u001b[39m\n",
       "\u001b[36m   1.19231e-02  1.19231e-02  1.19230e-02\u001b[39m\n",
       "\u001b[36m   1.03266e-05  1.03261e-05  1.03264e-05\u001b[39m\n",
       "\u001b[36m   1.46388e-11  1.40862e-11  1.39825e-11\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> fpv=zeros(2);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> koutb=nsoli(f!,x0,fv,fpv,JVec; \u001b[39m\n",
       "\u001b[36m              fixedeta=true, eta=.1, lmaxit=2, lsolver=\"bicgstab\");\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> koutb.history\u001b[39m\n",
       "\u001b[36m  6-element Vector{Float64}:\u001b[39m\n",
       "\u001b[36m   1.88791e+00\u001b[39m\n",
       "\u001b[36m   2.43120e-01\u001b[39m\n",
       "\u001b[36m   1.19231e-02\u001b[39m\n",
       "\u001b[36m   4.87500e-04\u001b[39m\n",
       "\u001b[36m   7.54236e-06\u001b[39m\n",
       "\u001b[36m   3.84646e-07\u001b[39m"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?nsoli"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 3.8.2: Benchmarking the H-equation with nsoli.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will begin by comparing the fastest solution from Chapter 2 with two variants of Newton-GMRES, one with fixed $\\eta = .1$ and one with the Eisenstat-Walker forcing term with $\\eta_{max}=.9$ and $\\gamma = .9$. I'll allocate 20 vectors for the Krylov basis in the array FPK. First we will compare the residual histories for a small version of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=512;\n",
    "FS=ones(n,); FPS=ones(n,n); FPS32=ones(Float32,n,n); \n",
    "    x0=ones(n,); c=.5; hdata = heqinit(x0, c);\n",
    "bargs=(atol = 1.e-10, rtol = 1.e-10, sham = 5, \n",
    "    resdec = .1, pdata=hdata);\n",
    "FPK=zeros(n,20);\n",
    "# Fixed eta = .1\n",
    "kbargs=(atol = 1.e-10, rtol = 1.e-10, eta=.1, \n",
    "    fixedeta=true, pdata=hdata);\n",
    "# Eisenstat-Walker\n",
    "kbargsew=(atol = 1.e-10, rtol = 1.e-10, eta=.9, \n",
    "    fixedeta=false, pdata=hdata);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll run the winner from Chapter 2 and __nsoli.jl__ with two choices\n",
    "for the forcing term.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "nout=nsol(heqf!, x0, FS, FPS32, heqJ!; bargs...);\n",
    "kout=nsoli(heqf!, x0, FS, FPK; kbargs...);\n",
    "koutew=nsoli(heqf!, x0, FS, FPK; kbargsew...);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They residual histories are essentially the same. The nonlinear residuals are a bit larger in the middle of the iteration for the Newton-Krylov iterations because of the tolerances for the linear iteration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×3 Matrix{Float64}:\n",
       " 3.49504e+00  3.49504e+00  3.49504e+00\n",
       " 1.79699e-02  4.98627e-02  4.98627e-02\n",
       " 1.55516e-04  1.84641e-03  1.84641e-03\n",
       " 1.33171e-06  1.82363e-04  1.82363e-04\n",
       " 1.13966e-08  2.34294e-06  2.34294e-06\n",
       " 9.75228e-11  2.42620e-11  2.42620e-11"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[nout.history kout.history koutew.history]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing the costs is harder. While a Jacobian-vector product for this problem has the same cost as a call to the function, the cost per iteration for nsol.jl is harder to evaluate in these terms. It's better to look at the benchmark results for a larger problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=4096;\n",
    "FS=ones(n); FPS=ones(n,n); FPS32=ones(Float32,n,n); \n",
    "    x0=ones(n); c=.5; hdata = heqinit(x0, c);\n",
    "bargs=(atol = 1.e-10, rtol = 1.e-10, sham = 5, \n",
    "    resdec = .1, pdata=hdata);\n",
    "FPK=zeros(n,20);\n",
    "kbargs=(atol = 1.e-10, rtol = 1.e-10, eta=.1, \n",
    "    fixedeta=true, pdata=hdata);\n",
    "kbargsew=(atol = 1.e-10, rtol = 1.e-10, eta=.9,\n",
    "    fixedeta=false, pdata=hdata);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shamanskii, n=5\n",
      "  94.134 ms (117 allocations: 739.11 KiB)\n",
      "Newton-GMRES, fixed eta\n",
      "  1.163 ms (409 allocations: 512.09 KiB)\n",
      "Newton-GMRES, Eisenstat-Walker\n",
      "  1.166 ms (409 allocations: 512.09 KiB)\n"
     ]
    }
   ],
   "source": [
    "println(\"Shamanskii, n=5\"); \n",
    "@btime nsol(heqf!, $x0, $FS, $FPS32, heqJ!; bargs...);\n",
    "println(\"Newton-GMRES, fixed eta\"); \n",
    "@btime nsoli(heqf!, $x0, $FS, $FPK; kbargs...);\n",
    "println(\"Newton-GMRES, Eisenstat-Walker\"); \n",
    "@btime nsoli(heqf!, $x0, $FS, $FPK; kbargsew...);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton-Krylov code is over 40 times faster. This is not unique to this problem. If your Jacobian is well-conditioned or you have a good preconditioner, as we do in the PDE example, Newton-Krylov should perform much better than any variation of Newton's method using direct linear solvers.\n",
    "\n",
    "The other interesting thing in this example is that the two forcing term choices performed equally well. \n",
    "\n",
    "Finally we will see if storing the Krylov basis in single precision improves matters. It's easy to do this by simply replacing the double precision storage for the Krylov basis ```FPK``` with the single precision one ```FPK32```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Newton-GMRES, fixed eta\n",
      "  1.190 ms (404 allocations: 504.16 KiB)\n",
      "Newton-GMRES, Eisenstat-Walker\n",
      "  1.185 ms (404 allocations: 504.16 KiB)\n"
     ]
    }
   ],
   "source": [
    "n=4096; FS=ones(n); \n",
    "FPK32=zeros(Float32,n,20)\n",
    "println(\"Newton-GMRES, fixed eta\"); \n",
    "@btime nsoli(heqf!, $x0, $FS, $FPK32; kbargs...);\n",
    "println(\"Newton-GMRES, Eisenstat-Walker\"); \n",
    "@btime nsoli(heqf!, $x0, $FS, $FPK32; kbargsew...);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is essentially no difference in performance between storing the basis in single and double precision. It is easy in hindsight to see why. Each function evaluation and forward difference Jacobian-vector product is $O(N \\log N)$ work. The cost of orthogonalization for $k$ GMRES iterations with classical Gram-Schmidt twice is $k^2 N$ (can you see why). So if we do $k$ Krylov iterations per Newton the cost of orthogonalization is $k^2 N$ and the cost of calls to the residual is $O(k N \\log N)$. The computation is dominated by the calls to the residual unless $k$ is very large. \n",
    "\n",
    "We will quantify this with a computation to look at the iteration statistics. The ```ijac``` field counts Jacobian-vector projects for each iteration.\n",
    "It is sufficient to look at the\n",
    "fixed $\\eta = .1$ case. The results for the Eisenstat-Walker forcing term are exactly the same.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 1, 1, 1, 1, 2]\n"
     ]
    }
   ],
   "source": [
    "fixedetaout = nsoli(heqf!, x0, FS, FPK; kbargs...);\n",
    "println(fixedetaout.stats.ijac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The statistics indicate that we converge after a single GMRES iteration and are taking a single Krylov per Newton for most of the iteration (remember that the initial iteration is $\\vs = 0$ when computing the Newton step). So the orthogonalization cost is $N$ and the function evaluation cost is $O(N \\log N)$. We would expect that storing the Krylov basis  in single precision would have very little benefit, and that is exactly what we see.\n",
    "\n",
    "We invite the reader to increase $c$ and the dimension of the problem to see if anything changes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 3.8.3: Preconditioning the Convection-Diffusion Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we will benchmark the Newton-GMRES iteration against the direct solvers from Chapter 2 and explore the differences between left and right preconditioning. We will begin by repeating the computation for the fastest version using __nsol.jl__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsol, sham=5\n",
      "  3.609 ms (587 allocations: 6.29 MiB)\n"
     ]
    }
   ],
   "source": [
    "n=31;\n",
    "# Get some room for the residual\n",
    "u0=zeros(n*n,);\n",
    "FV=copy(u0);\n",
    "# Get the precomputed data from pdeinit\n",
    "pdata=pdeinit(n)\n",
    "# Storage for the Jacobian, same sparsity pattern \n",
    "# as the discrete Laplacian\n",
    "J=copy(pdata.D2);\n",
    "# Iteration Parameters\n",
    "rtol=1.e-7\n",
    "atol=1.e-10\n",
    "println(\"nsol, sham=5\"); \n",
    "@btime nsol(pdeF!, $u0, $FV, $J, pdeJ!; \n",
    "    resdec=.5, rtol=rtol, atol=atol, pdata=pdata, sham=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll set up the problem for nsoli. We need to allocate storage for the Krylov basis. One case will be no preconditioning at all, so the Kryov basis will need more storage. The analytic Jacobian-vector product is __Jvec2d.jl__, which is in __TestProblems/EllipticPDE.jl__. The preconditioner is __Pvec2d.jl__ from __TestProblems/PDE_Tools.jl__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsoli, not preconditioned\n",
      "  3.877 ms (7291 allocations: 984.33 KiB)\n"
     ]
    }
   ],
   "source": [
    "# Storage for the Krylov basis\n",
    "    JV = zeros(n * n, 100)\n",
    "    eta=.1\n",
    "    fixedeta=false\n",
    "println(\"nsoli, not preconditioned\")\n",
    "@btime nsoli(pdeF!, $u0, $FV, $JV, Jvec2d; \n",
    "            rtol=rtol, atol=atol, Pvec=nothing, pdata=pdata, eta=eta,\n",
    "            fixedeta=fixedeta, pside=\"right\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even with no preconditioning, the iterative solver is almost twice as fast as __nsol.jl__ using the direct method. When you precondition, which we will do from the right for now, the difference is a factor of almost two over the solve without preconditioning. This difference would increase with a finer mesh. Try it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsoli, preconditioned, \n",
      "    Eisenstat-Walker forcing term\n",
      "  1.340 ms (1049 allocations: 649.41 KiB)\n"
     ]
    }
   ],
   "source": [
    "println(\"nsoli, preconditioned, \n",
    "    Eisenstat-Walker forcing term\")\n",
    "@btime nsoli(pdeF!, $u0, $FV, $JV, Jvec2d; \n",
    "             rtol=rtol, atol=atol, Pvec=Pvec2d, \n",
    "             pdata=pdata, eta=eta,\n",
    "             fixedeta=fixedeta, pside=\"right\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will benchmark with a fixed forcing term for our next example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsoli, preconditioned, fixed eta\n",
      "  1.682 ms (1289 allocations: 922.70 KiB)\n"
     ]
    }
   ],
   "source": [
    "fixedeta=true;\n",
    "println(\"nsoli, preconditioned, fixed eta\")\n",
    "@btime nsoli(pdeF!, $u0, $FV, $JV, Jvec2d;\n",
    "             rtol=rtol, atol=atol, Pvec=Pvec2d, \n",
    "             pdata=pdata, eta=eta,\n",
    "             fixedeta=fixedeta, pside=\"right\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we see that Eisenstat-Walker is a bit better. Finally, we return to Eisenstat-Walker with $\\eta_{max} = .9$. We see very little difference from $\\eta_{max}=.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsoli, preconditioned, \n",
      "    Eisenstat-Walker forcing term\n",
      "  1.329 ms (1064 allocations: 737.52 KiB)\n"
     ]
    }
   ],
   "source": [
    "eta=.9; fixedeta=false;\n",
    "println(\"nsoli, preconditioned, \n",
    "    Eisenstat-Walker forcing term\")\n",
    "@btime nsoli(pdeF!, $u0, FV, $JV, Jvec2d; \n",
    "             rtol=rtol, atol=atol, Pvec=Pvec2d, \n",
    "             pdata=pdata, eta=eta,\n",
    "             fixedeta=fixedeta, pside=\"right\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Left preconditioning? We'll see that even with $\\eta_{max}=.1$ it's a bit slower that right preconditioning. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsoli, left preconditioned, \n",
      "    Eisenstat-Walker forcing term\n",
      "  1.409 ms (1152 allocations: 741.39 KiB)\n"
     ]
    }
   ],
   "source": [
    "eta=.1\n",
    "fixedeta=false\n",
    "println(\"nsoli, left preconditioned, \n",
    "    Eisenstat-Walker forcing term\")\n",
    "@btime nsoli(pdeF!, $u0, $FV, $JV, Jvec2d; \n",
    "            rtol=rtol, atol=atol, Pvec=Pvec2d, \n",
    "            pdata=pdata, eta=eta,\n",
    "            fixedeta=fixedeta, pside=\"left\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we try left preconditioning with $\\eta_{max} = .9$. We plotted the results in Figure 3.3. While the number of nonlinear iterations is roughly double that of the right preconditioned version, the solver time is less than the number of nonlinear iterations would indicate. Can you figure out why that is?\n",
    "\n",
    "Note that we have to increase ```maxit``` to give the nonlinear solver enough iterations to overcome the mismatch between left preconditioning and the inexact Newton condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.360 ms (2028 allocations: 2.05 MiB)\n"
     ]
    }
   ],
   "source": [
    "eta=.9;\n",
    "@btime nsoli(pdeF!, $u0, $FV, $JV, Jvec2d; \n",
    "            rtol=rtol, atol=atol, Pvec=Pvec2d, \n",
    "            pdata=pdata, eta=eta, maxit=100,\n",
    "            fixedeta=fixedeta, pside=\"left\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 3.8.4: ptcsoli.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__ptcsoli.jl__ is our Newton-Krylov $\\ptc$ code. Herewith the docstrings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "search: \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mc\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1mi\u001b[22m \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mc\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22m \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mc\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22msc nsoli\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "function ptcsoli( F!, x0, FS, FPS, Jvec = dirder; rtol = 1.e-6, atol = 1.e-12,     maxit = 20, lmaxit = -1, lsolver = \"gmres\", eta = 0.1, fixedeta = true,     Pvec = nothing, PvecKnowsdelta = false, pside = \"right\", delta0 = 1.e-6,     dx = 1.e-7, pdata = nothing, printerr = true, keepsolhist = false, )\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  New for this book ==> ptcsoli\n",
       "\n",
       "PTC finds the steady-state solution of u' = -F(u), u(0) = u\\_0. The - sign is a convention.\n",
       "\n",
       "You must allocate storage for the function and Krylov basis in advance –> in the calling program <– ie. in FS and FPS\n",
       "\n",
       "Inputs:\n",
       "\n",
       "\\begin{itemize}\n",
       "\\item F!: function evaluation, the ! indicates that F! overwrites FS, your   preallocated storage for the function.\n",
       "\n",
       "So, FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)\n",
       "\n",
       "Your function MUST have –> return FS <– at the end.   See the example in TestProblems/Systems/FBeam!.jl\n",
       "\n",
       "\n",
       "\\item x0: initial iterate\n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\n",
       "\\item FS: Preallocated storage for function. It is a vector of size N\n",
       "\n",
       "You should store it as (N) and design F! to use vectors of size (N). If you use (N,1) consistently instead, the solvers may work, but I make no guarantees.\n",
       "\n",
       "\n",
       "\\item FPS: preallocated storage for the Krylov basis. It is an N x m matrix where      you plan to take at most m-1 GMRES iterations before a restart. \n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\n",
       "\\item Jvec: Jacobian vector product, If you leave this out the   default is a finite difference directional derivative.\n",
       "\n",
       "So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v. \n",
       "\n",
       "(v, FS, x) or (v, FS, x, pdata) must be the argument list,   even if FP does not need FS.   One reason for this is that the finite-difference derivative   does and that is the default in the solver.\n",
       "\n",
       "\n",
       "\\item Precision: Lemme tell ya 'bout precision. I designed this code for    full precision functions and linear algebra in any precision you want.    You can declare FPS as Float64 or Float32 and ptcsoli    will do the right thing. Float16 support is there, but not working well.\n",
       "\n",
       "If the Jacobian is reasonably well conditioned, you can cut the cost   of orthogonalization and storage (for GMRES) in half with no loss.   There is no benefit if your linear solver is not GMRES or if   orthogonalization and storage of the Krylov vectors is only a   small part of the cost of the computation. So if your preconditioner   is good and you only need a few Krylovs/Newton, reduced precision won't   help you much.\n",
       "\n",
       "BiCGSTAB does not benefit from reduced precision.\n",
       "\n",
       "\\end{itemize}\n",
       "\\rule{\\textwidth}{1pt}\n",
       "Keyword Arguments (kwargs):\n",
       "\n",
       "rtol and atol: relative and absolute error tolerances\n",
       "\n",
       "delta0: initial pseudo time step. The default value of 1.e-3 is a bit conservative and is one option you really should play with. Look at the example where I set it to 1.0!\n",
       "\n",
       "maxit: limit on nonlinear iterations, default=100. \n",
       "\n",
       "This is coupled to delta0. If your choice of delta0 is too small (conservative) then you'll need many iterations to converge and will need a larger value of maxit\n",
       "\n",
       "For PTC you'll need more iterations than for a straight-up nonlinear solve. This is part of the price for finding the  stable solution. \n",
       "\n",
       "lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has m columns, and you need more than m-1 linear iterations, then GMRES will restart.\n",
       "\n",
       "The default is -1. For GMRES this means that you'll take m-1 iterations, where size(V) = (n,m), and get no restarts. For BiCGSTAB you'll then get the default of 10 iterations.\n",
       "\n",
       "lsolver: the linear solver, default = \"gmres\"\n",
       "\n",
       "Your choices will be \"gmres\" or \"bicgstab\". However, gmres is the only option for now. \n",
       "\n",
       "eta and fixed eta: eta > 0 or there's an error.\n",
       "\n",
       "The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||\n",
       "\n",
       "where\n",
       "\n",
       "etag = eta if fixedeta=true\n",
       "\n",
       "etag = Eisenstat-Walker as implemented in book if fixedeta=false\n",
       "\n",
       "The default, which may change, is eta=.1, fixedeta=true \n",
       "\n",
       "Pvec: Preconditioner-vector product. The rules are similar to Jvec     So, Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where     P(x) is the preconditioner. You must use x as an input even     if your preconditioner does not depend on x.\n",
       "\n",
       "PvecKnowsdelta: If you want your preconditioner-vector product to depend on      the pseudo-timestep delta, put an array deltaval in your precomputed     data. Initialize it as     deltaval = zeros(1,)     and let ptcsoli know about it by setting the kwarg     PvecKnowsdelta = true     ptcsoli will update the value in deltaval with every change     to delta with pdata.deltaval[1]=delta     so your preconditioner-vector product can get to it.\n",
       "\n",
       "pside: apply preconditioner on pside, default = \"right\". I do not       recommend \"left\". The problem with \"left\" for ptcsoli is       that it can fail to satisfy the inexact Newton condition for        the unpreconditioned equation, especially early in the iteration       and lead to an incorrect result (unstable solution or wrong        branch of steady state).       See Chapter 3 for the story on this. \n",
       "\n",
       "dx: default = 1.e-7\n",
       "\n",
       "difference increment in finite-difference derivatives       h=dx*norm(x)+1.e-8 \n",
       "\n",
       "pdata:\n",
       "\n",
       "precomputed data for the function, Jacobian-vector, and Preconditioner-vector products.  Things will go better if you use this rather than hide the data in global variables within the module for your function/Jacobian\n",
       "\n",
       "If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of them. precomputed data for the function/Jacobian.  Things will go better if you use this rather than hide the data  in global variables within the module for your function/Jacobian. \n",
       "\n",
       "printerr: default = true\n",
       "\n",
       "I print a helpful message when the solver fails. To suppress that message set printerr to false. \n",
       "\n",
       "keepsolhist: default = false\n",
       "\n",
       "Set this to true to get the history of the iteration in the output tuple. This is on by default for scalar equations and off for systems. Only turn it on if you have use for the data, which can get REALLY LARGE.\n",
       "\n",
       "Output:\n",
       "\n",
       "A named tuple (solution, functionval, history, stats, idid,                errcode, solhist) where\n",
       "\n",
       "solution = converged result functionval = F(solution) history = the vector of residual norms (||F(x)||) for the iteration stats = named tuple of the history of (ifun, ijac, ikfail), the number of functions/jacobian-vector products/linear solver failures at each iteration.\n",
       "\n",
       "I do not count the function values for a finite-difference derivative because they count toward a Jacobian-vector product.\n",
       "\n",
       "Linear solver failures need not cause the nonlinear iteration to fail.  You get a warning and that is all. \n",
       "\n",
       "idid=true if the iteration succeeded and false if not. \n",
       "\n",
       "errcode = 0 if the iteration succeeded \n",
       "\n",
       "\\begin{verbatim}\n",
       "    = -1 if the initial iterate satisfies the termination criteria\n",
       "    = 10 if no convergence after maxit iterations\n",
       "\\end{verbatim}\n",
       "solhist:\n",
       "\n",
       "This is the entire history of the iteration if you've set keepsolhist=true\n",
       "\n",
       "solhist is an N x K array where N is the length of x and K is the number of iteration + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "\\subsubsection{Example for ptcsol}\n",
       "\\paragraph{The buckling beam problem.}\n",
       "You'll need to use TestProblems for this to work. The preconditioner is a solver for the high order term.\n",
       "\n",
       "\\begin{verbatim}\n",
       "julia> using SIAMFANLEquations.TestProblems\n",
       "\n",
       "julia> function PreCondBeam(v, x, bdata)\n",
       "          J = bdata.D2\n",
       "          ptv = J\\v\n",
       "       end\n",
       "PreCondBeam (generic function with 1 method)\n",
       "\n",
       "julia> n=63; maxit=1000; delta0 = 0.01; lambda = 20.0;\n",
       "\n",
       "julia> # Set up the precomputed data\n",
       "\n",
       "julia> bdata = beaminit(n, 0.0, lambda);\n",
       "\n",
       "julia> x = bdata.x; u0 = x .* (1.0 .- x) .* (2.0 .- x); \n",
       "\n",
       "julia> u0 .*= exp.(-10.0 * u0); FS = copy(u0); FPJV=zeros(n,20);\n",
       "\n",
       "julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; \n",
       "       delta0 = delta0, pdata = bdata, eta = 1.e-2, \n",
       "       rtol = 1.e-10, maxit = maxit, Pvec = PreCondBeam);\n",
       "\n",
       "julia> # It takes a few iterations to get there.\n",
       "       length(pout.history)\n",
       "25\n",
       "\n",
       "julia> [pout.history[1:5] pout.history[21:25]]\n",
       "5×2 Matrix{Float64}:\n",
       " 6.31230e+01  1.79574e+00\n",
       " 7.45927e+00  2.65956e-01\n",
       " 8.73595e+00  6.58220e-03\n",
       " 2.91937e+01  8.34114e-06\n",
       " 3.47970e+01  5.06847e-09\n",
       "\n",
       "julia> # We get the nonnegative stedy state.\n",
       "julia> maximum(pout.solution)\n",
       "2.19086e+00\n",
       "\n",
       "julia> # Now use BiCGSTAB for the linear solver\n",
       "\n",
       "julia> FPJV=zeros(n);\n",
       "\n",
       "julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; \n",
       "       delta0 = delta0, pdata = bdata,\n",
       "       eta = 1.e-2, rtol = 1.e-10, maxit = maxit, \n",
       "       Pvec = PreCondBeam, lsolver=\"bicgstab\");\n",
       "\n",
       "julia> # Same number of iterations as GMRES, but each one costs double \n",
       "\n",
       "julia> # the Jacobian-vector products and much less storage\n",
       "\n",
       "julia> length(pout.history)\n",
       "25\n",
       "\n",
       "julia> [pout.history[1:5] pout.history[21:25]]\n",
       "5×2 Matrix{Float64}:\n",
       " 6.31230e+01  1.68032e+00\n",
       " 7.47081e+00  2.35073e-01\n",
       " 8.62095e+00  5.18260e-03\n",
       " 2.96495e+01  3.23803e-06\n",
       "\\end{verbatim}\n"
      ],
      "text/markdown": [
       "function ptcsoli( F!, x0, FS, FPS, Jvec = dirder; rtol = 1.e-6, atol = 1.e-12,     maxit = 20, lmaxit = -1, lsolver = \"gmres\", eta = 0.1, fixedeta = true,     Pvec = nothing, PvecKnowsdelta = false, pside = \"right\", delta0 = 1.e-6,     dx = 1.e-7, pdata = nothing, printerr = true, keepsolhist = false, )\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  New for this book ==> ptcsoli\n",
       "\n",
       "PTC finds the steady-state solution of u' = -F(u), u(0) = u_0. The - sign is a convention.\n",
       "\n",
       "You must allocate storage for the function and Krylov basis in advance –> in the calling program <– ie. in FS and FPS\n",
       "\n",
       "Inputs:\n",
       "\n",
       "  * F!: function evaluation, the ! indicates that F! overwrites FS, your   preallocated storage for the function.\n",
       "\n",
       "    So, FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)\n",
       "\n",
       "    Your function MUST have –> return FS <– at the end.   See the example in TestProblems/Systems/FBeam!.jl\n",
       "  * x0: initial iterate\n",
       "\n",
       "  * FS: Preallocated storage for function. It is a vector of size N\n",
       "\n",
       "    You should store it as (N) and design F! to use vectors of size (N). If you use (N,1) consistently instead, the solvers may work, but I make no guarantees.\n",
       "  * FPS: preallocated storage for the Krylov basis. It is an N x m matrix where      you plan to take at most m-1 GMRES iterations before a restart.\n",
       "\n",
       "  * Jvec: Jacobian vector product, If you leave this out the   default is a finite difference directional derivative.\n",
       "\n",
       "    So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v. \n",
       "\n",
       "    (v, FS, x) or (v, FS, x, pdata) must be the argument list,   even if FP does not need FS.   One reason for this is that the finite-difference derivative   does and that is the default in the solver.\n",
       "  * Precision: Lemme tell ya 'bout precision. I designed this code for    full precision functions and linear algebra in any precision you want.    You can declare FPS as Float64 or Float32 and ptcsoli    will do the right thing. Float16 support is there, but not working well.\n",
       "\n",
       "    If the Jacobian is reasonably well conditioned, you can cut the cost   of orthogonalization and storage (for GMRES) in half with no loss.   There is no benefit if your linear solver is not GMRES or if   orthogonalization and storage of the Krylov vectors is only a   small part of the cost of the computation. So if your preconditioner   is good and you only need a few Krylovs/Newton, reduced precision won't   help you much.\n",
       "\n",
       "    BiCGSTAB does not benefit from reduced precision.\n",
       "\n",
       "---\n",
       "\n",
       "Keyword Arguments (kwargs):\n",
       "\n",
       "rtol and atol: relative and absolute error tolerances\n",
       "\n",
       "delta0: initial pseudo time step. The default value of 1.e-3 is a bit conservative and is one option you really should play with. Look at the example where I set it to 1.0!\n",
       "\n",
       "maxit: limit on nonlinear iterations, default=100. \n",
       "\n",
       "This is coupled to delta0. If your choice of delta0 is too small (conservative) then you'll need many iterations to converge and will need a larger value of maxit\n",
       "\n",
       "For PTC you'll need more iterations than for a straight-up nonlinear solve. This is part of the price for finding the  stable solution. \n",
       "\n",
       "lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has m columns, and you need more than m-1 linear iterations, then GMRES will restart.\n",
       "\n",
       "The default is -1. For GMRES this means that you'll take m-1 iterations, where size(V) = (n,m), and get no restarts. For BiCGSTAB you'll then get the default of 10 iterations.\n",
       "\n",
       "lsolver: the linear solver, default = \"gmres\"\n",
       "\n",
       "Your choices will be \"gmres\" or \"bicgstab\". However, gmres is the only option for now. \n",
       "\n",
       "eta and fixed eta: eta > 0 or there's an error.\n",
       "\n",
       "The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||\n",
       "\n",
       "where\n",
       "\n",
       "etag = eta if fixedeta=true\n",
       "\n",
       "etag = Eisenstat-Walker as implemented in book if fixedeta=false\n",
       "\n",
       "The default, which may change, is eta=.1, fixedeta=true \n",
       "\n",
       "Pvec: Preconditioner-vector product. The rules are similar to Jvec     So, Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where     P(x) is the preconditioner. You must use x as an input even     if your preconditioner does not depend on x.\n",
       "\n",
       "PvecKnowsdelta: If you want your preconditioner-vector product to depend on      the pseudo-timestep delta, put an array deltaval in your precomputed     data. Initialize it as     deltaval = zeros(1,)     and let ptcsoli know about it by setting the kwarg     PvecKnowsdelta = true     ptcsoli will update the value in deltaval with every change     to delta with pdata.deltaval[1]=delta     so your preconditioner-vector product can get to it.\n",
       "\n",
       "pside: apply preconditioner on pside, default = \"right\". I do not       recommend \"left\". The problem with \"left\" for ptcsoli is       that it can fail to satisfy the inexact Newton condition for        the unpreconditioned equation, especially early in the iteration       and lead to an incorrect result (unstable solution or wrong        branch of steady state).       See Chapter 3 for the story on this. \n",
       "\n",
       "dx: default = 1.e-7\n",
       "\n",
       "difference increment in finite-difference derivatives       h=dx*norm(x)+1.e-8 \n",
       "\n",
       "pdata:\n",
       "\n",
       "precomputed data for the function, Jacobian-vector, and Preconditioner-vector products.  Things will go better if you use this rather than hide the data in global variables within the module for your function/Jacobian\n",
       "\n",
       "If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of them. precomputed data for the function/Jacobian.  Things will go better if you use this rather than hide the data  in global variables within the module for your function/Jacobian. \n",
       "\n",
       "printerr: default = true\n",
       "\n",
       "I print a helpful message when the solver fails. To suppress that message set printerr to false. \n",
       "\n",
       "keepsolhist: default = false\n",
       "\n",
       "Set this to true to get the history of the iteration in the output tuple. This is on by default for scalar equations and off for systems. Only turn it on if you have use for the data, which can get REALLY LARGE.\n",
       "\n",
       "Output:\n",
       "\n",
       "A named tuple (solution, functionval, history, stats, idid,                errcode, solhist) where\n",
       "\n",
       "solution = converged result functionval = F(solution) history = the vector of residual norms (||F(x)||) for the iteration stats = named tuple of the history of (ifun, ijac, ikfail), the number of functions/jacobian-vector products/linear solver failures at each iteration.\n",
       "\n",
       "I do not count the function values for a finite-difference derivative because they count toward a Jacobian-vector product.\n",
       "\n",
       "Linear solver failures need not cause the nonlinear iteration to fail.  You get a warning and that is all. \n",
       "\n",
       "idid=true if the iteration succeeded and false if not. \n",
       "\n",
       "errcode = 0 if the iteration succeeded \n",
       "\n",
       "```\n",
       "    = -1 if the initial iterate satisfies the termination criteria\n",
       "    = 10 if no convergence after maxit iterations\n",
       "```\n",
       "\n",
       "solhist:\n",
       "\n",
       "This is the entire history of the iteration if you've set keepsolhist=true\n",
       "\n",
       "solhist is an N x K array where N is the length of x and K is the number of iteration + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "### Example for ptcsol\n",
       "\n",
       "#### The buckling beam problem.\n",
       "\n",
       "You'll need to use TestProblems for this to work. The preconditioner is a solver for the high order term.\n",
       "\n",
       "```jldoctest\n",
       "julia> using SIAMFANLEquations.TestProblems\n",
       "\n",
       "julia> function PreCondBeam(v, x, bdata)\n",
       "          J = bdata.D2\n",
       "          ptv = J\\v\n",
       "       end\n",
       "PreCondBeam (generic function with 1 method)\n",
       "\n",
       "julia> n=63; maxit=1000; delta0 = 0.01; lambda = 20.0;\n",
       "\n",
       "julia> # Set up the precomputed data\n",
       "\n",
       "julia> bdata = beaminit(n, 0.0, lambda);\n",
       "\n",
       "julia> x = bdata.x; u0 = x .* (1.0 .- x) .* (2.0 .- x); \n",
       "\n",
       "julia> u0 .*= exp.(-10.0 * u0); FS = copy(u0); FPJV=zeros(n,20);\n",
       "\n",
       "julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; \n",
       "       delta0 = delta0, pdata = bdata, eta = 1.e-2, \n",
       "       rtol = 1.e-10, maxit = maxit, Pvec = PreCondBeam);\n",
       "\n",
       "julia> # It takes a few iterations to get there.\n",
       "       length(pout.history)\n",
       "25\n",
       "\n",
       "julia> [pout.history[1:5] pout.history[21:25]]\n",
       "5×2 Matrix{Float64}:\n",
       " 6.31230e+01  1.79574e+00\n",
       " 7.45927e+00  2.65956e-01\n",
       " 8.73595e+00  6.58220e-03\n",
       " 2.91937e+01  8.34114e-06\n",
       " 3.47970e+01  5.06847e-09\n",
       "\n",
       "julia> # We get the nonnegative stedy state.\n",
       "julia> maximum(pout.solution)\n",
       "2.19086e+00\n",
       "\n",
       "julia> # Now use BiCGSTAB for the linear solver\n",
       "\n",
       "julia> FPJV=zeros(n);\n",
       "\n",
       "julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; \n",
       "       delta0 = delta0, pdata = bdata,\n",
       "       eta = 1.e-2, rtol = 1.e-10, maxit = maxit, \n",
       "       Pvec = PreCondBeam, lsolver=\"bicgstab\");\n",
       "\n",
       "julia> # Same number of iterations as GMRES, but each one costs double \n",
       "\n",
       "julia> # the Jacobian-vector products and much less storage\n",
       "\n",
       "julia> length(pout.history)\n",
       "25\n",
       "\n",
       "julia> [pout.history[1:5] pout.history[21:25]]\n",
       "5×2 Matrix{Float64}:\n",
       " 6.31230e+01  1.68032e+00\n",
       " 7.47081e+00  2.35073e-01\n",
       " 8.62095e+00  5.18260e-03\n",
       " 2.96495e+01  3.23803e-06\n",
       "```\n"
      ],
      "text/plain": [
       "  function ptcsoli( F!, x0, FS, FPS, Jvec = dirder; rtol = 1.e-6, atol =\n",
       "  1.e-12, maxit = 20, lmaxit = -1, lsolver = \"gmres\", eta = 0.1, fixedeta =\n",
       "  true, Pvec = nothing, PvecKnowsdelta = false, pside = \"right\", delta0 =\n",
       "  1.e-6, dx = 1.e-7, pdata = nothing, printerr = true, keepsolhist = false, )\n",
       "\n",
       "  C. T. Kelley, 2022\n",
       "\n",
       "  Julia versions of the nonlinear solvers from my SIAM books. New for this\n",
       "  book ==> ptcsoli\n",
       "\n",
       "  PTC finds the steady-state solution of u' = -F(u), u(0) = u_0. The - sign is\n",
       "  a convention.\n",
       "\n",
       "  You must allocate storage for the function and Krylov basis in advance –> in\n",
       "  the calling program <– ie. in FS and FPS\n",
       "\n",
       "  Inputs:\n",
       "\n",
       "    •  F!: function evaluation, the ! indicates that F! overwrites FS,\n",
       "       your preallocated storage for the function.\n",
       "       So, FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)\n",
       "       Your function MUST have –> return FS <– at the end. See the\n",
       "       example in TestProblems/Systems/FBeam!.jl\n",
       "\n",
       "    •  x0: initial iterate\n",
       "\n",
       "    •  FS: Preallocated storage for function. It is a vector of size N\n",
       "       You should store it as (N) and design F! to use vectors of size\n",
       "       (N). If you use (N,1) consistently instead, the solvers may work,\n",
       "       but I make no guarantees.\n",
       "\n",
       "    •  FPS: preallocated storage for the Krylov basis. It is an N x m\n",
       "       matrix where you plan to take at most m-1 GMRES iterations before\n",
       "       a restart.\n",
       "\n",
       "    •  Jvec: Jacobian vector product, If you leave this out the default\n",
       "       is a finite difference directional derivative.\n",
       "       So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v.\n",
       "       (v, FS, x) or (v, FS, x, pdata) must be the argument list, even if\n",
       "       FP does not need FS. One reason for this is that the\n",
       "       finite-difference derivative does and that is the default in the\n",
       "       solver.\n",
       "\n",
       "    •  Precision: Lemme tell ya 'bout precision. I designed this code for\n",
       "       full precision functions and linear algebra in any precision you\n",
       "       want. You can declare FPS as Float64 or Float32 and ptcsoli will\n",
       "       do the right thing. Float16 support is there, but not working\n",
       "       well.\n",
       "       If the Jacobian is reasonably well conditioned, you can cut the\n",
       "       cost of orthogonalization and storage (for GMRES) in half with no\n",
       "       loss. There is no benefit if your linear solver is not GMRES or if\n",
       "       orthogonalization and storage of the Krylov vectors is only a\n",
       "       small part of the cost of the computation. So if your\n",
       "       preconditioner is good and you only need a few Krylovs/Newton,\n",
       "       reduced precision won't help you much.\n",
       "       BiCGSTAB does not benefit from reduced precision.\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "  Keyword Arguments (kwargs):\n",
       "\n",
       "  rtol and atol: relative and absolute error tolerances\n",
       "\n",
       "  delta0: initial pseudo time step. The default value of 1.e-3 is a bit\n",
       "  conservative and is one option you really should play with. Look at the\n",
       "  example where I set it to 1.0!\n",
       "\n",
       "  maxit: limit on nonlinear iterations, default=100.\n",
       "\n",
       "  This is coupled to delta0. If your choice of delta0 is too small\n",
       "  (conservative) then you'll need many iterations to converge and will need a\n",
       "  larger value of maxit\n",
       "\n",
       "  For PTC you'll need more iterations than for a straight-up nonlinear solve.\n",
       "  This is part of the price for finding the stable solution.\n",
       "\n",
       "  lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has m\n",
       "  columns, and you need more than m-1 linear iterations, then GMRES will\n",
       "  restart.\n",
       "\n",
       "  The default is -1. For GMRES this means that you'll take m-1 iterations,\n",
       "  where size(V) = (n,m), and get no restarts. For BiCGSTAB you'll then get the\n",
       "  default of 10 iterations.\n",
       "\n",
       "  lsolver: the linear solver, default = \"gmres\"\n",
       "\n",
       "  Your choices will be \"gmres\" or \"bicgstab\". However, gmres is the only\n",
       "  option for now.\n",
       "\n",
       "  eta and fixed eta: eta > 0 or there's an error.\n",
       "\n",
       "  The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||\n",
       "\n",
       "  where\n",
       "\n",
       "  etag = eta if fixedeta=true\n",
       "\n",
       "  etag = Eisenstat-Walker as implemented in book if fixedeta=false\n",
       "\n",
       "  The default, which may change, is eta=.1, fixedeta=true\n",
       "\n",
       "  Pvec: Preconditioner-vector product. The rules are similar to Jvec So,\n",
       "  Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where P(x) is the\n",
       "  preconditioner. You must use x as an input even if your preconditioner does\n",
       "  not depend on x.\n",
       "\n",
       "  PvecKnowsdelta: If you want your preconditioner-vector product to depend on\n",
       "  the pseudo-timestep delta, put an array deltaval in your precomputed data.\n",
       "  Initialize it as deltaval = zeros(1,) and let ptcsoli know about it by\n",
       "  setting the kwarg PvecKnowsdelta = true ptcsoli will update the value in\n",
       "  deltaval with every change to delta with pdata.deltaval[1]=delta so your\n",
       "  preconditioner-vector product can get to it.\n",
       "\n",
       "  pside: apply preconditioner on pside, default = \"right\". I do not recommend\n",
       "  \"left\". The problem with \"left\" for ptcsoli is that it can fail to satisfy\n",
       "  the inexact Newton condition for the unpreconditioned equation, especially\n",
       "  early in the iteration and lead to an incorrect result (unstable solution or\n",
       "  wrong branch of steady state). See Chapter 3 for the story on this.\n",
       "\n",
       "  dx: default = 1.e-7\n",
       "\n",
       "  difference increment in finite-difference derivatives h=dx*norm(x)+1.e-8\n",
       "\n",
       "  pdata:\n",
       "\n",
       "  precomputed data for the function, Jacobian-vector, and\n",
       "  Preconditioner-vector products. Things will go better if you use this rather\n",
       "  than hide the data in global variables within the module for your\n",
       "  function/Jacobian\n",
       "\n",
       "  If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of\n",
       "  them. precomputed data for the function/Jacobian. Things will go better if\n",
       "  you use this rather than hide the data in global variables within the module\n",
       "  for your function/Jacobian.\n",
       "\n",
       "  printerr: default = true\n",
       "\n",
       "  I print a helpful message when the solver fails. To suppress that message\n",
       "  set printerr to false.\n",
       "\n",
       "  keepsolhist: default = false\n",
       "\n",
       "  Set this to true to get the history of the iteration in the output tuple.\n",
       "  This is on by default for scalar equations and off for systems. Only turn it\n",
       "  on if you have use for the data, which can get REALLY LARGE.\n",
       "\n",
       "  Output:\n",
       "\n",
       "  A named tuple (solution, functionval, history, stats, idid, errcode,\n",
       "  solhist) where\n",
       "\n",
       "  solution = converged result functionval = F(solution) history = the vector\n",
       "  of residual norms (||F(x)||) for the iteration stats = named tuple of the\n",
       "  history of (ifun, ijac, ikfail), the number of functions/jacobian-vector\n",
       "  products/linear solver failures at each iteration.\n",
       "\n",
       "  I do not count the function values for a finite-difference derivative\n",
       "  because they count toward a Jacobian-vector product.\n",
       "\n",
       "  Linear solver failures need not cause the nonlinear iteration to fail. You\n",
       "  get a warning and that is all.\n",
       "\n",
       "  idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "  errcode = 0 if the iteration succeeded\n",
       "\n",
       "\u001b[36m      = -1 if the initial iterate satisfies the termination criteria\u001b[39m\n",
       "\u001b[36m      = 10 if no convergence after maxit iterations\u001b[39m\n",
       "\n",
       "  solhist:\n",
       "\n",
       "  This is the entire history of the iteration if you've set keepsolhist=true\n",
       "\n",
       "  solhist is an N x K array where N is the length of x and K is the number of\n",
       "  iteration + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "\u001b[1m  Example for ptcsol\u001b[22m\n",
       "\u001b[1m  ––––––––––––––––––\u001b[22m\n",
       "\n",
       "\u001b[1m  The buckling beam problem.\u001b[22m\n",
       "\u001b[1m  --------------------------\u001b[22m\n",
       "\n",
       "  You'll need to use TestProblems for this to work. The preconditioner is a\n",
       "  solver for the high order term.\n",
       "\n",
       "\u001b[36m  julia> using SIAMFANLEquations.TestProblems\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> function PreCondBeam(v, x, bdata)\u001b[39m\n",
       "\u001b[36m            J = bdata.D2\u001b[39m\n",
       "\u001b[36m            ptv = J\\v\u001b[39m\n",
       "\u001b[36m         end\u001b[39m\n",
       "\u001b[36m  PreCondBeam (generic function with 1 method)\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> n=63; maxit=1000; delta0 = 0.01; lambda = 20.0;\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # Set up the precomputed data\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> bdata = beaminit(n, 0.0, lambda);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> x = bdata.x; u0 = x .* (1.0 .- x) .* (2.0 .- x); \u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> u0 .*= exp.(-10.0 * u0); FS = copy(u0); FPJV=zeros(n,20);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; \u001b[39m\n",
       "\u001b[36m         delta0 = delta0, pdata = bdata, eta = 1.e-2, \u001b[39m\n",
       "\u001b[36m         rtol = 1.e-10, maxit = maxit, Pvec = PreCondBeam);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # It takes a few iterations to get there.\u001b[39m\n",
       "\u001b[36m         length(pout.history)\u001b[39m\n",
       "\u001b[36m  25\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> [pout.history[1:5] pout.history[21:25]]\u001b[39m\n",
       "\u001b[36m  5×2 Matrix{Float64}:\u001b[39m\n",
       "\u001b[36m   6.31230e+01  1.79574e+00\u001b[39m\n",
       "\u001b[36m   7.45927e+00  2.65956e-01\u001b[39m\n",
       "\u001b[36m   8.73595e+00  6.58220e-03\u001b[39m\n",
       "\u001b[36m   2.91937e+01  8.34114e-06\u001b[39m\n",
       "\u001b[36m   3.47970e+01  5.06847e-09\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # We get the nonnegative stedy state.\u001b[39m\n",
       "\u001b[36m  julia> maximum(pout.solution)\u001b[39m\n",
       "\u001b[36m  2.19086e+00\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # Now use BiCGSTAB for the linear solver\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> FPJV=zeros(n);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; \u001b[39m\n",
       "\u001b[36m         delta0 = delta0, pdata = bdata,\u001b[39m\n",
       "\u001b[36m         eta = 1.e-2, rtol = 1.e-10, maxit = maxit, \u001b[39m\n",
       "\u001b[36m         Pvec = PreCondBeam, lsolver=\"bicgstab\");\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # Same number of iterations as GMRES, but each one costs double \u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # the Jacobian-vector products and much less storage\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> length(pout.history)\u001b[39m\n",
       "\u001b[36m  25\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> [pout.history[1:5] pout.history[21:25]]\u001b[39m\n",
       "\u001b[36m  5×2 Matrix{Float64}:\u001b[39m\n",
       "\u001b[36m   6.31230e+01  1.68032e+00\u001b[39m\n",
       "\u001b[36m   7.47081e+00  2.35073e-01\u001b[39m\n",
       "\u001b[36m   8.62095e+00  5.18260e-03\u001b[39m\n",
       "\u001b[36m   2.96495e+01  3.23803e-06\u001b[39m"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?ptcsoli"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3.8.5: Benchmarking $\\ptc$ with the buckling beam problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will set up the beam problem as we did before. Remember that __bdata.D2__ is the discrete Laplacian in one space dimension, which we compute within the initialization function __beaminit__. We will start with __ptcsol.jl__ to remind you what we did before and solve a larger problem to compare using a direct solver with GMRES."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=1023; lambda=20; delta=.01; maxit=1000; \n",
    "bdata = beaminit(n, 0.0, lambda); \n",
    "x = bdata.x; u0 = x .* (1.0 .- x) .* (2.0 .- x); \n",
    "u0 .*= exp.(-10.0 * u0);\n",
    "FS = copy(u0); FPS=copy(bdata.D2); FPJV = zeros(n, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll benchmark the solve. Remember that ```FBeam!``` and ```BeamJ!``` are defined in the TestProblems submodule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  482.208 μs (1103 allocations: 2.42 MiB)\n"
     ]
    }
   ],
   "source": [
    "@btime ptcsol(FBeam!, $u0, $FS, $FPS, BeamJ!; \n",
    "             rtol=1.e-10, pdata=bdata, \n",
    "             delta0=delta, maxit=maxit);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test ptcsoli we will use the $\\delta$-dependent preconditioner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ptvbeamdelta (generic function with 1 method)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ptvbeamdelta(v, x, bdata)\n",
    "    delta = bdata.deltaval[1]\n",
    "    J = bdata.D2 + (1.0 / delta) * I\n",
    "    ptv = J \\ v\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.543 ms (3841 allocations: 5.97 MiB)\n"
     ]
    }
   ],
   "source": [
    "@btime ptcsoli(FBeam!, $u0, $FS, $FPJV; \n",
    "    lsolver=\"gmres\", delta0=delta, pdata=bdata, \n",
    "    lmaxit=19, eta=1.e-2, Pvec=ptvbeamdelta,\n",
    "    pside=\"right\", PvecKnowsdelta=true, maxit=maxit);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the iterative linear solver costs over two times as much as the direct solver. This is no surprise as the application of the preconditioner requires a tridiagonal solve, which is the same cost as solving the equation for the Newton step with a direct method. The buckling beam problem is simply not hard enough to benefit from an iterative linear solver. The reader should try increasing $n$ to see if anything changes, but should keep in mind that one may need to reduce $\\delta_0$ as $n$ increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3.9: Projects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Low Storage Solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Benchmark the solves for the H-equation and the convection-diffusion equation using BiCGSTAB and GMRES(m) for the linear solvers. How do the runtimes and memory allocations compare to full GMRES? How do the runtimes and allocations depend on the dimension and $m$ for GMRES(m)? Do things change for $c=1$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh Independence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An iteration for a discretization of a differential or integral equation is _mesh-independent_ if the iteration statistics are independent of the grid. Nonlinear iterations are usually mesh-independent if the discretization is reasonably well-done <cite data-cite=\"allg\"><a href=\"siamfa.html#allg\">(ABPR86)</cite>. That is not the case, however, for the linear solves. One can only get mesh-independence for the linear solve if the preconditioner is so good that it essentially converts the problem into an integral equation. For the H-equation and the convection-diffusion (both preconditioned and not), vary the grid size and see how the iteration statistics change. Use both full GMRES and the low-storage solvers. You will want to make figures like the ones earlier in this chapter that plot residual norm against both the number of nonlinear iterations and the number of Jacobian-vector products."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Playing with the convection term"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vary the convection term $C$ in the convection-diffusion equation\n",
    "$$\n",
    "-\\nabla^2 u + C u ( u_x + u_y) = f\n",
    "$$\n",
    "where you use the boundary conditions and exact solution $u^*$ from [Chapter 2](SIAMFANLCh2.ipynb). Vary $C$ from $C=20$ ( the choice in our examples ) to $C = 1000$ or larger. Does the choice of forcing term affect the iteration? What happens to the linear and nonlinear iteration statistics? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next notebook = [Chapter 4: Fixed Point Problems  and Anderson Acceleration](SIAMFANLCh4.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.0",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}