SIAMFANLCh2s.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\n",
    "\\newcommand{\\calf}{{\\cal F}}\n",
    "\\newcommand{\\dnu}{d \\nu}\n",
    "\\newcommand{\\mf}{{\\bf F}}\n",
    "\\newcommand{\\vu}{{\\bf u}}\n",
    "\\newcommand{\\vx}{{\\bf x}}\n",
    "\\newcommand{\\vy}{{\\bf y}}\n",
    "\\newcommand{\\vs}{{\\bf s}}\n",
    "\\newcommand{\\ve}{{\\bf e}}\n",
    "\\newcommand{\\vd}{{\\bf d}}\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{\\ball}{{\\cal B}}\n",
    "\\newcommand{\\ptc}{{\\Psi TC}}\n",
    "\\newcommand{\\diag}{\\mbox{diag}}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "include(\"fanote_init.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.8 Solvers for Chapter 2\n",
    "\n",
    "Contents for Section 2.8\n",
    "\n",
    "[Overview](#Overview)\n",
    "\n",
    "- [nsol.jl](#Section-2.8.1:-nsol.jl)\n",
    "\n",
    "- [H-equation Revisited](#Section-2.8.2:-H-Equation-revisited)\n",
    "\n",
    "- [More on the Two-Point BVP](#Section-2.8.3:-More-on-the-Two-Point-BVP)\n",
    "\n",
    "- [Shamanskii for the Convection-Diffusion Problem](#Section-2.8.4:-Shamanskii-for-the-Convection-Diffusion-Problem)\n",
    "\n",
    "- [ptcsol.jl](#Section-2.8.5:-ptcsol.jl)\n",
    "\n",
    "- [More on the Buckling Beam](#More-on-the-Buckling-Beam)\n",
    "\n",
    "- [Section 2.9: Projects](#Section-2.9:-Projects)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will follow the pattern of Chapter 1 and present two solvers, a Newton code and a $\\ptc$ code. Both codes are for systems of equations and use direct methods to compute the step. We returned the solution history for the simple two dimensional example in Section 2.6, but will not do that again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 2.8.1: nsol.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__nsol.jl__ solves systems of nonlinear equations and computes the Newton step with direct linear solvers. Let's look at the docstrings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "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[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22mi \u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22msc aasol \u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1ml\u001b[22mheq htol \u001b[0m\u001b[1mn\u001b[22mor \u001b[0m\u001b[1mn\u001b[22mtoh ptcsol ishold cool Bool\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "nsol(F!, x0, FS, FPS, J!=diffjac!; rtol=1.e-6, atol=1.e-12,\n",
       "           maxit=20, solver=\"newton\", sham=5, armmax=10, resdec=.1,\n",
       "           dx = 1.e-7, armfix=false, \n",
       "           pdata = nothing, jfact = klfact,\n",
       "           printerr = true, keepsolhist = false, stagnationok=false)\n",
       "\\end{verbatim}\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  Herewith: nsol\n",
       "\n",
       "You must allocate storage for the function and Jacobian 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 and in TestProblems/Systems/simple.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 Jacobian. It is an N x N matrix\n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\n",
       "\\item J!: Jacobian evaluation, the ! indicates that J! overwrites FPS, your   preallocated storage for the Jacobian. If you leave this out the   default is a finite difference Jacobian.\n",
       "\n",
       "So, FP=J!(FP,FS,x) or FP=J!(FP,FS,x,pdata) returns FP=F'(x). \n",
       "\n",
       "(FP,FS, x) must be the argument list, even if FP does not need FS.   One reason for this is that the finite-difference Jacobian   does and that is the default in the solver.\n",
       "\n",
       "Your Jacobian function MUST have –> return FP <– at the end.   See the examples in the docstrings and in TestProblems/Systems/simple.jl\n",
       "\n",
       "\\end{itemize}\n",
       "\\begin{itemize}\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, Float32, or Float16 and nsol will do the right thing if   YOU do not destroy the declaration in your J! function. I'm amazed   that this works so easily. If the Jacobian is reasonably well    conditioned, you can cut the cost of Jacobian factorization and   storage in half with no loss. For large dense Jacobians and inexpensive   functions, this is a good deal.\n",
       "\n",
       "BUT ... There is very limited support for direct sparse solvers in   anything other than Float64. I recommend that you only use Float64   with direct sparse solvers unless you really know what you're doing. I   have a couple examples in the notebook, but watch out.\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",
       "solver: default = \"newton\"\n",
       "\n",
       "Your choices are \"newton\" or \"chord\". However, you have sham at your disposal only if you chose newton. \"chord\" will keep using the initial derivative until the iterate converges, uses the iteration budget, or the line search fails. It is not the same as sham=Inf, which is smarter.\n",
       "\n",
       "sham: default = 5 (ie Newton)\n",
       "\n",
       "This is the Shamanskii method. If sham=1, you have Newton. The iteration updates the derivative every sham iterations. The convergence rate has local q-order sham+1 if you only count iterations where you update the derivative. You need not provide your own derivative function to use this option. sham=Inf is chord only if chord is converging well.\n",
       "\n",
       "I made sham=1 the default for scalar equations. For systems I'm more aggressive and want to invest as little energy in linear algebra as possible. So the default is sham=5.\n",
       "\n",
       "armmax: upper bound on step size reductions in line search\n",
       "\n",
       "resdec: default = .1\n",
       "\n",
       "This is the target value for residual reduction. The default value is .1. In the old MATLAB codes it was .5. I only turn Shamanskii on if the residuals are decreasing rapidly, at least a factor of resdec, and the line search is quiescent. If you want to eliminate resdec from the method ( you don't ) then set resdec = 1.0 and you will never hear from it again.\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.  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 either of F! or J!, you must use in in the  calling sequence of both.\n",
       "\n",
       "jfact: default = klfact (tries to figure out best choice) \n",
       "\n",
       "If your Jacobian has any special structure, please set jfact to the correct choice for a factorization.\n",
       "\n",
       "I use jfact when I call PrepareJac! to evaluate the Jacobian (using your J!) and factor it. The default is to use klfact (an internal function) to do something reasonable. For general dense matrices, klfact picks lu! to compute an LU factorization and share storage with the Jacobian.  You may change LU to something else by, for example, setting jfact = cholesky! if your Jacobian is spd.\n",
       "\n",
       "klfact knows about banded matrices and picks qr. You should, however RTFM, allocate the extra two upper bands, and use jfact=qr! to override klfact.\n",
       "\n",
       "klfact uses lu for general sparse matrices.\n",
       "\n",
       "If you give me something that klfact does not know how to dispatch on, then nothing happens. I just return the original Jacobian matrix and  nsol will use backslash to compute the Newton step. I know that this is probably not optimal in your situation, so it is  good to pick something else, like jfact = lu.\n",
       "\n",
       "If you want to manage your own factorization within your Jacobian  evaluation function, then set\n",
       "\n",
       "jfact = nofact\n",
       "\n",
       "and nsol will not attempt to factor your Jacobian. That is also what happens when klfact does not know what to do. Your Jacobian is sent directly to Julia's {\\textbackslash}  operation\n",
       "\n",
       "Please do not mess with the line that calls PrepareJac!. \n",
       "\n",
       "\\begin{verbatim}\n",
       "    FPF = PrepareJac!(FPS, FS, x, ItRules)\n",
       "\\end{verbatim}\n",
       "FPF is not the same as FPS (the storage you allocate for the Jacobian) for a reason. FPF and FPS do not have the same type, even though they share storage. So, FPS=PrepareJac!(FPS, FS, ...) will break things.\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",
       "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), the number of functions/derivatives/steplength reductions at each iteration.\n",
       "\n",
       "I do not count the function values for a finite-difference derivative because they count toward a Jacobian evaluation. \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",
       "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 iterations + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "\\rule{\\textwidth}{1pt}\n",
       "\\subsubsection{Examples for nsol}\n",
       "\\paragraph{World's easiest problem example.}\n",
       "Test 64 and 32 bit Jacobians. No meaningful difference in the residual histories or the converged solutions.\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",
       "#\n",
       "# The return fv part is important even though f! overwrites fv.\n",
       "#\n",
       "       return fv\n",
       "       end\n",
       "f (generic function with 1 method)\n",
       "\n",
       "julia> x=ones(2); fv=zeros(2); jv=zeros(2,2); jv32=zeros(Float32,2,2);\n",
       "julia> nout=nsol(f!,x,fv,jv; sham=1);\n",
       "julia> nout32=nsol(f!,x,fv,jv32; sham=1);\n",
       "julia> [nout.history nout32.history]\n",
       "5×2 Matrix{Float64}:\n",
       " 1.88791e+00  1.88791e+00\n",
       " 2.43119e-01  2.43120e-01\n",
       " 1.19231e-02  1.19231e-02\n",
       " 1.03266e-05  1.03265e-05\n",
       " 1.46388e-11  1.45995e-11\n",
       "\n",
       "julia> [nout.solution nout.solution-nout32.solution]\n",
       "2×2 Array{Float64,2}:\n",
       " -7.39085e-01  -5.48450e-14\n",
       "  2.30988e+00  -2.26485e-14\n",
       "\\end{verbatim}\n",
       "\\paragraph{H-equation example.}\n",
       "I'm taking the sham=5 default here, so the convergence is not quadratic. The good news is that we evaluate the Jacobian only once.\n",
       "\n",
       "\\begin{verbatim}\n",
       "julia> n=16; x0=ones(n); FV=ones(n); JV=ones(n,n);\n",
       "julia> hdata=heqinit(x0, .5);\n",
       "julia> hout=nsol(heqf!,x0,FV,JV;pdata=hdata);\n",
       "julia> hout.history\n",
       "4-element Array{Float64,1}:\n",
       " 6.17376e-01\n",
       " 3.17810e-03\n",
       " 2.75227e-05\n",
       " 2.35817e-07\n",
       "\\end{verbatim}\n"
      ],
      "text/markdown": [
       "```\n",
       "nsol(F!, x0, FS, FPS, J!=diffjac!; rtol=1.e-6, atol=1.e-12,\n",
       "           maxit=20, solver=\"newton\", sham=5, armmax=10, resdec=.1,\n",
       "           dx = 1.e-7, armfix=false, \n",
       "           pdata = nothing, jfact = klfact,\n",
       "           printerr = true, keepsolhist = false, stagnationok=false)\n",
       "```\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  Herewith: nsol\n",
       "\n",
       "You must allocate storage for the function and Jacobian 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 and in TestProblems/Systems/simple.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 Jacobian. It is an N x N matrix\n",
       "\n",
       "  * J!: Jacobian evaluation, the ! indicates that J! overwrites FPS, your   preallocated storage for the Jacobian. If you leave this out the   default is a finite difference Jacobian.\n",
       "\n",
       "    So, FP=J!(FP,FS,x) or FP=J!(FP,FS,x,pdata) returns FP=F'(x). \n",
       "\n",
       "    (FP,FS, x) must be the argument list, even if FP does not need FS.   One reason for this is that the finite-difference Jacobian   does and that is the default in the solver.\n",
       "\n",
       "    Your Jacobian function MUST have –> return FP <– at the end.   See the examples in the docstrings and in TestProblems/Systems/simple.jl\n",
       "\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, Float32, or Float16 and nsol will do the right thing if   YOU do not destroy the declaration in your J! function. I'm amazed   that this works so easily. If the Jacobian is reasonably well    conditioned, you can cut the cost of Jacobian factorization and   storage in half with no loss. For large dense Jacobians and inexpensive   functions, this is a good deal.\n",
       "\n",
       "    BUT ... There is very limited support for direct sparse solvers in   anything other than Float64. I recommend that you only use Float64   with direct sparse solvers unless you really know what you're doing. I   have a couple examples in the notebook, but watch out.\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",
       "solver: default = \"newton\"\n",
       "\n",
       "Your choices are \"newton\" or \"chord\". However, you have sham at your disposal only if you chose newton. \"chord\" will keep using the initial derivative until the iterate converges, uses the iteration budget, or the line search fails. It is not the same as sham=Inf, which is smarter.\n",
       "\n",
       "sham: default = 5 (ie Newton)\n",
       "\n",
       "This is the Shamanskii method. If sham=1, you have Newton. The iteration updates the derivative every sham iterations. The convergence rate has local q-order sham+1 if you only count iterations where you update the derivative. You need not provide your own derivative function to use this option. sham=Inf is chord only if chord is converging well.\n",
       "\n",
       "I made sham=1 the default for scalar equations. For systems I'm more aggressive and want to invest as little energy in linear algebra as possible. So the default is sham=5.\n",
       "\n",
       "armmax: upper bound on step size reductions in line search\n",
       "\n",
       "resdec: default = .1\n",
       "\n",
       "This is the target value for residual reduction. The default value is .1. In the old MATLAB codes it was .5. I only turn Shamanskii on if the residuals are decreasing rapidly, at least a factor of resdec, and the line search is quiescent. If you want to eliminate resdec from the method ( you don't ) then set resdec = 1.0 and you will never hear from it again.\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.  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 either of F! or J!, you must use in in the  calling sequence of both.\n",
       "\n",
       "jfact: default = klfact (tries to figure out best choice) \n",
       "\n",
       "If your Jacobian has any special structure, please set jfact to the correct choice for a factorization.\n",
       "\n",
       "I use jfact when I call PrepareJac! to evaluate the Jacobian (using your J!) and factor it. The default is to use klfact (an internal function) to do something reasonable. For general dense matrices, klfact picks lu! to compute an LU factorization and share storage with the Jacobian.  You may change LU to something else by, for example, setting jfact = cholesky! if your Jacobian is spd.\n",
       "\n",
       "klfact knows about banded matrices and picks qr. You should, however RTFM, allocate the extra two upper bands, and use jfact=qr! to override klfact.\n",
       "\n",
       "klfact uses lu for general sparse matrices.\n",
       "\n",
       "If you give me something that klfact does not know how to dispatch on, then nothing happens. I just return the original Jacobian matrix and  nsol will use backslash to compute the Newton step. I know that this is probably not optimal in your situation, so it is  good to pick something else, like jfact = lu.\n",
       "\n",
       "If you want to manage your own factorization within your Jacobian  evaluation function, then set\n",
       "\n",
       "jfact = nofact\n",
       "\n",
       "and nsol will not attempt to factor your Jacobian. That is also what happens when klfact does not know what to do. Your Jacobian is sent directly to Julia's \\  operation\n",
       "\n",
       "Please do not mess with the line that calls PrepareJac!. \n",
       "\n",
       "```\n",
       "    FPF = PrepareJac!(FPS, FS, x, ItRules)\n",
       "```\n",
       "\n",
       "FPF is not the same as FPS (the storage you allocate for the Jacobian) for a reason. FPF and FPS do not have the same type, even though they share storage. So, FPS=PrepareJac!(FPS, FS, ...) will break things.\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",
       "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), the number of functions/derivatives/steplength reductions at each iteration.\n",
       "\n",
       "I do not count the function values for a finite-difference derivative because they count toward a Jacobian evaluation. \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",
       "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 iterations + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "---\n",
       "\n",
       "### Examples for nsol\n",
       "\n",
       "#### World's easiest problem example.\n",
       "\n",
       "Test 64 and 32 bit Jacobians. No meaningful difference in the residual histories or the converged solutions.\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",
       "#\n",
       "# The return fv part is important even though f! overwrites fv.\n",
       "#\n",
       "       return fv\n",
       "       end\n",
       "f (generic function with 1 method)\n",
       "\n",
       "julia> x=ones(2); fv=zeros(2); jv=zeros(2,2); jv32=zeros(Float32,2,2);\n",
       "julia> nout=nsol(f!,x,fv,jv; sham=1);\n",
       "julia> nout32=nsol(f!,x,fv,jv32; sham=1);\n",
       "julia> [nout.history nout32.history]\n",
       "5×2 Matrix{Float64}:\n",
       " 1.88791e+00  1.88791e+00\n",
       " 2.43119e-01  2.43120e-01\n",
       " 1.19231e-02  1.19231e-02\n",
       " 1.03266e-05  1.03265e-05\n",
       " 1.46388e-11  1.45995e-11\n",
       "\n",
       "julia> [nout.solution nout.solution-nout32.solution]\n",
       "2×2 Array{Float64,2}:\n",
       " -7.39085e-01  -5.48450e-14\n",
       "  2.30988e+00  -2.26485e-14\n",
       "```\n",
       "\n",
       "#### H-equation example.\n",
       "\n",
       "I'm taking the sham=5 default here, so the convergence is not quadratic. The good news is that we evaluate the Jacobian only once.\n",
       "\n",
       "```jldoctest\n",
       "julia> n=16; x0=ones(n); FV=ones(n); JV=ones(n,n);\n",
       "julia> hdata=heqinit(x0, .5);\n",
       "julia> hout=nsol(heqf!,x0,FV,JV;pdata=hdata);\n",
       "julia> hout.history\n",
       "4-element Array{Float64,1}:\n",
       " 6.17376e-01\n",
       " 3.17810e-03\n",
       " 2.75227e-05\n",
       " 2.35817e-07\n",
       "```\n"
      ],
      "text/plain": [
       "\u001b[36m  nsol(F!, x0, FS, FPS, J!=diffjac!; rtol=1.e-6, atol=1.e-12,\u001b[39m\n",
       "\u001b[36m             maxit=20, solver=\"newton\", sham=5, armmax=10, resdec=.1,\u001b[39m\n",
       "\u001b[36m             dx = 1.e-7, armfix=false, \u001b[39m\n",
       "\u001b[36m             pdata = nothing, jfact = klfact,\u001b[39m\n",
       "\u001b[36m             printerr = true, keepsolhist = false, stagnationok=false)\u001b[39m\n",
       "\n",
       "  C. T. Kelley, 2022\n",
       "\n",
       "  Julia versions of the nonlinear solvers from my SIAM books. Herewith: nsol\n",
       "\n",
       "  You must allocate storage for the function and Jacobian in advance –> in the\n",
       "  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 and in TestProblems/Systems/simple.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 Jacobian. It is an N x N matrix\n",
       "\n",
       "    •  J!: Jacobian evaluation, the ! indicates that J! overwrites FPS,\n",
       "       your preallocated storage for the Jacobian. If you leave this out\n",
       "       the default is a finite difference Jacobian.\n",
       "       So, FP=J!(FP,FS,x) or FP=J!(FP,FS,x,pdata) returns FP=F'(x).\n",
       "       (FP,FS, x) must be the argument list, even if FP does not need FS.\n",
       "       One reason for this is that the finite-difference Jacobian does\n",
       "       and that is the default in the solver.\n",
       "       Your Jacobian function MUST have –> return FP <– at the end. See\n",
       "       the examples in the docstrings and in\n",
       "       TestProblems/Systems/simple.jl\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, Float32, or Float16 and nsol\n",
       "       will do the right thing if YOU do not destroy the declaration in\n",
       "       your J! function. I'm amazed that this works so easily. If the\n",
       "       Jacobian is reasonably well conditioned, you can cut the cost of\n",
       "       Jacobian factorization and storage in half with no loss. For large\n",
       "       dense Jacobians and inexpensive functions, this is a good deal.\n",
       "       BUT ... There is very limited support for direct sparse solvers in\n",
       "       anything other than Float64. I recommend that you only use Float64\n",
       "       with direct sparse solvers unless you really know what you're\n",
       "       doing. I have a couple examples in the notebook, but watch out.\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",
       "  solver: default = \"newton\"\n",
       "\n",
       "  Your choices are \"newton\" or \"chord\". However, you have sham at your\n",
       "  disposal only if you chose newton. \"chord\" will keep using the initial\n",
       "  derivative until the iterate converges, uses the iteration budget, or the\n",
       "  line search fails. It is not the same as sham=Inf, which is smarter.\n",
       "\n",
       "  sham: default = 5 (ie Newton)\n",
       "\n",
       "  This is the Shamanskii method. If sham=1, you have Newton. The iteration\n",
       "  updates the derivative every sham iterations. The convergence rate has local\n",
       "  q-order sham+1 if you only count iterations where you update the derivative.\n",
       "  You need not provide your own derivative function to use this option.\n",
       "  sham=Inf is chord only if chord is converging well.\n",
       "\n",
       "  I made sham=1 the default for scalar equations. For systems I'm more\n",
       "  aggressive and want to invest as little energy in linear algebra as\n",
       "  possible. So the default is sham=5.\n",
       "\n",
       "  armmax: upper bound on step size reductions in line search\n",
       "\n",
       "  resdec: default = .1\n",
       "\n",
       "  This is the target value for residual reduction. The default value is .1. In\n",
       "  the old MATLAB codes it was .5. I only turn Shamanskii on if the residuals\n",
       "  are decreasing rapidly, at least a factor of resdec, and the line search is\n",
       "  quiescent. If you want to eliminate resdec from the method ( you don't )\n",
       "  then set resdec = 1.0 and you will never hear from it again.\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. Things will go better if you use\n",
       "  this rather than hide the data in global variables within the module for\n",
       "  your function/Jacobian\n",
       "\n",
       "  If you use pdata in either of F! or J!, you must use in in the calling\n",
       "  sequence of both.\n",
       "\n",
       "  jfact: default = klfact (tries to figure out best choice)\n",
       "\n",
       "  If your Jacobian has any special structure, please set jfact to the correct\n",
       "  choice for a factorization.\n",
       "\n",
       "  I use jfact when I call PrepareJac! to evaluate the Jacobian (using your J!)\n",
       "  and factor it. The default is to use klfact (an internal function) to do\n",
       "  something reasonable. For general dense matrices, klfact picks lu! to\n",
       "  compute an LU factorization and share storage with the Jacobian. You may\n",
       "  change LU to something else by, for example, setting jfact = cholesky! if\n",
       "  your Jacobian is spd.\n",
       "\n",
       "  klfact knows about banded matrices and picks qr. You should, however RTFM,\n",
       "  allocate the extra two upper bands, and use jfact=qr! to override klfact.\n",
       "\n",
       "  klfact uses lu for general sparse matrices.\n",
       "\n",
       "  If you give me something that klfact does not know how to dispatch on, then\n",
       "  nothing happens. I just return the original Jacobian matrix and nsol will\n",
       "  use backslash to compute the Newton step. I know that this is probably not\n",
       "  optimal in your situation, so it is good to pick something else, like jfact\n",
       "  = lu.\n",
       "\n",
       "  If you want to manage your own factorization within your Jacobian evaluation\n",
       "  function, then set\n",
       "\n",
       "  jfact = nofact\n",
       "\n",
       "  and nsol will not attempt to factor your Jacobian. That is also what happens\n",
       "  when klfact does not know what to do. Your Jacobian is sent directly to\n",
       "  Julia's \\ operation\n",
       "\n",
       "  Please do not mess with the line that calls PrepareJac!.\n",
       "\n",
       "\u001b[36m      FPF = PrepareJac!(FPS, FS, x, ItRules)\u001b[39m\n",
       "\n",
       "  FPF is not the same as FPS (the storage you allocate for the Jacobian) for a\n",
       "  reason. FPF and FPS do not have the same type, even though they share\n",
       "  storage. So, FPS=PrepareJac!(FPS, FS, ...) will break things.\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",
       "  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), the number of\n",
       "  functions/derivatives/steplength reductions at each iteration.\n",
       "\n",
       "  I do not count the function values for a finite-difference derivative\n",
       "  because they count toward a Jacobian evaluation.\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",
       "  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",
       "  iterations + 1. So, for scalar equations, it's a row vector.\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "\u001b[1m  Examples for nsol\u001b[22m\n",
       "\u001b[1m  –––––––––––––––––\u001b[22m\n",
       "\n",
       "\u001b[1m  World's easiest problem example.\u001b[22m\n",
       "\u001b[1m  --------------------------------\u001b[22m\n",
       "\n",
       "  Test 64 and 32 bit Jacobians. No meaningful difference in the residual\n",
       "  histories or the converged solutions.\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  #\u001b[39m\n",
       "\u001b[36m  # The return fv part is important even though f! overwrites fv.\u001b[39m\n",
       "\u001b[36m  #\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> x=ones(2); fv=zeros(2); jv=zeros(2,2); jv32=zeros(Float32,2,2);\u001b[39m\n",
       "\u001b[36m  julia> nout=nsol(f!,x,fv,jv; sham=1);\u001b[39m\n",
       "\u001b[36m  julia> nout32=nsol(f!,x,fv,jv32; sham=1);\u001b[39m\n",
       "\u001b[36m  julia> [nout.history nout32.history]\u001b[39m\n",
       "\u001b[36m  5×2 Matrix{Float64}:\u001b[39m\n",
       "\u001b[36m   1.88791e+00  1.88791e+00\u001b[39m\n",
       "\u001b[36m   2.43119e-01  2.43120e-01\u001b[39m\n",
       "\u001b[36m   1.19231e-02  1.19231e-02\u001b[39m\n",
       "\u001b[36m   1.03266e-05  1.03265e-05\u001b[39m\n",
       "\u001b[36m   1.46388e-11  1.45995e-11\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> [nout.solution nout.solution-nout32.solution]\u001b[39m\n",
       "\u001b[36m  2×2 Array{Float64,2}:\u001b[39m\n",
       "\u001b[36m   -7.39085e-01  -5.48450e-14\u001b[39m\n",
       "\u001b[36m    2.30988e+00  -2.26485e-14\u001b[39m\n",
       "\n",
       "\u001b[1m  H-equation example.\u001b[22m\n",
       "\u001b[1m  -------------------\u001b[22m\n",
       "\n",
       "  I'm taking the sham=5 default here, so the convergence is not quadratic. The\n",
       "  good news is that we evaluate the Jacobian only once.\n",
       "\n",
       "\u001b[36m  julia> n=16; x0=ones(n); FV=ones(n); JV=ones(n,n);\u001b[39m\n",
       "\u001b[36m  julia> hdata=heqinit(x0, .5);\u001b[39m\n",
       "\u001b[36m  julia> hout=nsol(heqf!,x0,FV,JV;pdata=hdata);\u001b[39m\n",
       "\u001b[36m  julia> hout.history\u001b[39m\n",
       "\u001b[36m  4-element Array{Float64,1}:\u001b[39m\n",
       "\u001b[36m   6.17376e-01\u001b[39m\n",
       "\u001b[36m   3.17810e-03\u001b[39m\n",
       "\u001b[36m   2.75227e-05\u001b[39m\n",
       "\u001b[36m   2.35817e-07\u001b[39m"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?nsol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Input\n",
    "\n",
    "The calling sequence for the Newton solvers in this book are similar, differing mostly in the management of the linear solver and memory allocation. The calling sequence for __nsol.jl__ is\n",
    "\n",
    "```julia\n",
    "function nsol(\n",
    "    F!,\n",
    "    x0,\n",
    "    FS,\n",
    "    FPS,\n",
    "    J! = diffjac!;\n",
    "    rtol = 1.e-6,\n",
    "    atol = 1.e-12,\n",
    "    maxit = 20,\n",
    "    solver = \"newton\",\n",
    "    sham = 1,\n",
    "    armmax = 10,\n",
    "    resdec = 0.1,\n",
    "    dx = 1.e-7,\n",
    "    armfix = false,\n",
    "    pdata = nothing,\n",
    "    jfact = lu!,\n",
    "    printerr = true,\n",
    "    keepsolhist = false,\n",
    "    stagnationok = false,\n",
    ")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we said earlier in the chapter, the calling sequence has some new things which are not in __nsolsc.jl__. The most significant are the arrays __FS__ and __FPS__, which preallocate\n",
    "storage for the function and Jacobian. As we have pointed out earlier, the farther upstream one allocates memory the better, so __nsol.jl__ insists that you allocate an vector __FS__ of the same size as the initial iterate and a matrix \n",
    "__FPS__ for the Jacobian. \n",
    "\n",
    "For problems in several variables, the keyword argument __pdata__ is very important. This is the data structure for you to store any precomputed or preallocated data your function evaluation needs. You will almost surely need __pdata__ for any but the most trivial problems. Most of the examples from in this chapter use __pdata__ in a serious manner.\n",
    "\n",
    "You should dimension __x0__ either as $(N)$ (or equivalently $(N,)$) and dimension __FS__ the same way. __nsol.jl__ expects vectors to be in double precision (Float64). Unless you are porting old code you should\n",
    "dimension $\\vx$ and $\\mf$ as vectors $(N)$ not as two dimensional arrays with one column $(N,1)$. The solvers \n",
    "should handle both, but I cannot guarantee that.\n",
    "\n",
    "The __!__ in the function evaluation __F!__ is to indicate that __nsold__ expects __F__  to overwrite its input. So,\n",
    "the way to call __F!__ is to preallocate the storage for the function value in an array __FS__ and then call the function as\n",
    "```Julia\n",
    "F!(FS,x)\n",
    "```\n",
    "or\n",
    "```Julia\n",
    "F!(FS,x,pdata)\n",
    "```\n",
    "__nsol.jl__ will figure out if you have populated __pdata__ or left it alone as the default value of __nothing__.\n",
    "\n",
    "\n",
    "And now for the Jacobian. __nsol.jl__ uses direct methods for linear algebra. If your matrix is dense, the default is to use Julia's __lu!__ function to do an LU factorization. If your matrix is symmetric or symmetric positive definite you can use the __jfact__ keyword to change __lu!__ to __ldlt!__ or __cholesky!__ for example. __nsold.jl__ assumes that the factorization you ask for will overwrite the matrix. Hence, the __factorize__ function in Julia is not what you want for this application.\n",
    "\n",
    "If the Jacobian is sparse, __nsol.jl__ will use __lu__ instead of __lu!__. The reason\n",
    "for this is that one cannot overwrite the Jacobian with the factorization\n",
    "in the general sparse case.\n",
    "\n",
    "You will also need to preallocate storage for the Jacobian in the array __FPS__. You may use any legal real precision for __FPS__. Float64 is the default. If you use Float32 you cut the storage for the matrix and the time for the factorization in half. We recommend that you do this if your Jacobian is dense. If you are using the __Sparsesuite__ sparse solvers, then you must store the Jacobian in double precision. __Sparsesuite__ does not support lower precision.\n",
    "\n",
    "Your Jacobian computation __J!__ must also overwrite it's input. The call looks like\n",
    "```julia\n",
    "J!(FV,FP,x)\n",
    "```\n",
    "or \n",
    "```julia\n",
    "J!(FV,FP,x,pdata)\n",
    "\n",
    "```\n",
    "returns FP=F'(x). The input FP=F(x), which __nsol.jl__ has already computed, has to be there. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.8.2: H-Equation revisited"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we do several experiments to illustrate the advantages of infrequent evaluation and factorization of the Jacobian and mixed precision computation. We will begin with a function that will support our testing. We want to investigate combinations of Newton's method/Shamanskii with ```sham=5``` (the default), storing and factoring the Jacobian in double/single precision, and analytic/forward difference Jacobians. To make this easy we write a function that solves the H-equation with __nsol.jl__ and lets us vary these cases. \n",
    "\n",
    "The functions for the residual __heqf!.jl__, the Jacobian __heqJ!.jl__, and the precomputed data \n",
    "__heqinit.jl__ are in the large file \n",
    "[src/TestProblems/Systems/Hequation.jl](https://github.com/ctkelley/SIAMFANLEquations.jl/blob/master/src/TestProblems/Systems/Hequation.jl).\n",
    "\n",
    "\n",
    "I'm passing the precomputed data to the function rather than computing it within. This keeps the cost of the precomputed data out of the benchmarking I'll do later.\n",
    "\n",
    "We use __splat__ in this example. We populate a named tuple ```bargs``` to keep the keyword arguments in a convenient place and then, when it's time to give it to __nsol.jl__, the call looks like ```bargs...```. The three dots are the __splat__ and tell __nsol__ to expand bargs and harvest the keyword arguments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "htest (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function htest(x0, FS, FPS, hdata; analytic=false, hsham=5)\n",
    "    n=length(FS)\n",
    "    #\n",
    "    # I've preallocaed x0, FS, and FPS. But they may have \n",
    "    # been changed by previous runs.\n",
    "    # The cost of resetting their entries to 1.0 is insignificant. \n",
    "    #\n",
    "    FS.=1.0\n",
    "    FPS.=1.0\n",
    "    bargs=(atol = 1.e-10, rtol = 1.e-10, sham = hsham, \n",
    "        resdec = .1, pdata=hdata)\n",
    "    if analytic\n",
    "        nout=nsol( heqf!, x0, FS, FPS, heqJ!; bargs...)\n",
    "    else\n",
    "        nout=nsol( heqf!, x0, FS, FPS; bargs...)\n",
    "    end\n",
    "    return nout\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin we will compare the iteration histories for four cases. We will consider analytic and forward difference Jacobians with the storage and factorization of the Jacobian done in double and single precision. __Theorem 1.2__ says the results should be almost indistinguishable. We will begin with Newton's method.\n",
    "\n",
    "All we need to do to store and factor Jacobians is to allocate the storage in single precision. That allocation is the line ```FPS32=ones(Float32,n,n);```. Note that we must reset ```FS``` and ```FPS``` before each new call to __nsol.jl__ because the solver uses the storage for residuals and Jacobians for the entire iteration. We do this with __broadcast__ after the initial allocation ```.=1.0``` instead of ```=ones(n,n)``` to avoid reallocation of the Jacobian.\n",
    "\n",
    "We will print all the residual histories in an array. The history vectors are the same length and are very hard to tell apart until the residual norm is one iteration from stagnation. This is just what the theory predicts. The theory \n",
    "(see <cite data-cite=\"ctk:sirev20\"><a href=\"siamfa.html#ctk:sirev20\">(Kel22)</cite> ) also predicts that there will be little difference between double precision linear algebra and single precision. We observe this as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " 1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00\n",
       " 5.14148e-03  5.14137e-03  5.14148e-03  5.14103e-03\n",
       " 1.00479e-07  1.01127e-07  1.00479e-07  9.72316e-08\n",
       " 2.00602e-15  8.51532e-14  2.00602e-15  1.67562e-13"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n=1024; FS=ones(n); FPS=ones(n,n); FPS32=ones(Float32,n,n); \n",
    "x0=ones(n); c=.5; hdata = heqinit(x0, c);\n",
    "nouta64=htest(x0, FS, FPS, hdata;\n",
    "    analytic=true, hsham=1);\n",
    "nouta32=htest(x0, FS, FPS32, hdata; \n",
    "    analytic=true, hsham=1);\n",
    "noutfd64=htest(x0, FS, FPS, hdata; \n",
    "    analytic=true, hsham=1);\n",
    "noutfd32=htest(x0, FS, FPS32,hdata; \n",
    "    analytic=false, hsham=1);\n",
    "hinit=nouta64.history[1]\n",
    "ha64=nouta64.history./hinit\n",
    "ha32=nouta32.history./hinit\n",
    "hd64=noutfd64.history./hinit\n",
    "hd32=noutfd32.history./hinit\n",
    "[ha64 ha32 hd64 hd32]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now do the same thing with the default setting of ```sham=5```. The theory correctly predicts that we will see more nonlinear iterations. We will be using __BenchmarkTools.jl__ to compare the costs later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Matrix{Float64}:\n",
       " 1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00\n",
       " 5.14148e-03  5.14137e-03  5.14112e-03  5.14103e-03\n",
       " 4.44954e-05  4.44942e-05  4.44899e-05  4.44895e-05\n",
       " 3.81019e-07  3.81011e-07  3.80954e-07  3.80950e-07\n",
       " 3.26071e-09  3.26063e-09  3.26001e-09  3.25997e-09\n",
       " 2.79022e-11  2.79020e-11  2.78955e-11  2.78950e-11"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nouta64=htest(x0, FS, FPS, hdata; analytic=true);\n",
    "noutfd64=htest(x0, FS, FPS, hdata; analytic=false);\n",
    "nouta32=htest(x0, FS, FPS32, hdata; analytic=true);\n",
    "noutfd32=htest(x0, FS, FPS32, hdata; analytic=false);\n",
    "hinit=nouta64.history[1]\n",
    "ha64=nouta64.history./hinit\n",
    "ha32=nouta32.history./hinit\n",
    "hd64=noutfd64.history./hinit\n",
    "hd32=noutfd32.history./hinit\n",
    "[ha64 ha32 hd64 hd32]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is interesting. We need more iterations, but we evaluate the Jacobian only once for Shamanskii. We can see this by looking at the ```stats``` field of the output tuple. They are all the same, so we will use ```nouta64```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(ifun = [1, 1, 1, 1, 1, 1], ijac = [0, 1, 0, 0, 0, 0], iarm = [0, 0, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nouta64.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpretation is that we do a function evaluation at all iterations and a single Jacobian evaluation to compute $\\vx_1$. We do no Jacobian work after that. The default in __nsol.jl__ is to reevaluate the Jacobian if the reduction in the residual norm larger than ```resdec = .1```. You can change ```resdec``` in the keyword arguments.\n",
    "The ```iarm``` field in the stats tells us that a line search was not necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) package to look at performance. \n",
    "The ```@btime``` command will show compute time and memory allocations for an average of several runs. The averaging will mitigate the effects of the compile time for the first run. \n",
    "\n",
    "To begin, we will compare the four versions of Newton's method. Note the ```$``` in in front of the array arguments. This is __interpolation__ and using it is important if you want to get accurate results from \n",
    "```@btime```.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytic, double\n",
      "  12.430 ms (93 allocations: 178.53 KiB)\n",
      "finite difference, double\n",
      "  63.888 ms (9300 allocations: 634.34 KiB)\n",
      "analytic, single\n",
      "  10.234 ms (111 allocations: 202.91 KiB)\n",
      "finite difference, single\n",
      "  61.139 ms (9318 allocations: 658.72 KiB)\n"
     ]
    }
   ],
   "source": [
    "println(\"analytic, double\"); \n",
    "@btime htest($x0, $FS, $FPS, hdata; \n",
    "    analytic=true, hsham=1);\n",
    "println(\"finite difference, double\"); \n",
    "@btime htest($x0, $FS, $FPS, hdata; \n",
    "    analytic=false, hsham=1);\n",
    "println(\"analytic, single\"); \n",
    "@btime htest($x0, $FS, $FPS32, hdata; \n",
    "    analytic=true, hsham=1);\n",
    "println(\"finite difference, single\"); \n",
    "@btime htest($x0, $FS, $FPS32, hdata; \n",
    "    analytic=false, hsham=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even for a problem of only modest size, the differences between the analytic Jacobian and the forward difference are significant. The differences between single and double precision linear algebra are, at least for the analytic Jacobian, roughly the factor of two we'd expect if the matrix factorization dominated the computation. For the forward difference Jacobian, we see that the cost of the Jacobian evaluation dominates everything else.\n",
    "\n",
    "Next, we look at the default ```sham=5``` from __nsol.jl__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytic, double\n",
      "  5.179 ms (83 allocations: 146.34 KiB)\n",
      "finite difference, double\n",
      "  22.209 ms (3152 allocations: 298.28 KiB)\n",
      "analytic, single\n",
      "  3.916 ms (113 allocations: 186.97 KiB)\n",
      "finite difference, single\n",
      "  20.926 ms (3182 allocations: 338.91 KiB)\n"
     ]
    }
   ],
   "source": [
    "println(\"analytic, double\"); \n",
    "@btime htest($x0, $FS, $FPS, hdata; analytic=true);\n",
    "println(\"finite difference, double\"); \n",
    "@btime htest($x0, $FS, $FPS, hdata; analytic=false);\n",
    "println(\"analytic, single\"); \n",
    "@btime htest($x0, $FS, $FPS32, hdata; analytic=true);\n",
    "println(\"finite difference, single\"); \n",
    "@btime htest($x0, $FS, $FPS32, hdata; analytic=false);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The moral here is pretty clear. We see that compute time is cut by a factor of at least two over Newton's method in all cases. The result for an analytic Jacobian with linear algebra in single precision and ```sham=5``` is 19 times faster than our slowest computation (Newton + finite difference Jacobian + double precision linear algebra). So, do less linear algebra and do it in single precision.\n",
    "\n",
    "Finally, we will increase the dimension. As we do that the computation becomes more burdensome, so we will only do four cases, all with an analytic Jacobian. The size of this example is large enough to clearly show the factor of two reduction in cost one gets from a single precision Jacobian.\n",
    "\n",
    "The reader with lots of free time should try these cases with a forward difference Jacobian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytic, double, Newton\n",
      "  435.297 ms (93 allocations: 706.53 KiB)\n",
      "analytic, single, Newton\n",
      "  251.679 ms (111 allocations: 802.91 KiB)\n",
      "analytic, double, sham=5\n",
      "  167.976 ms (83 allocations: 578.34 KiB)\n",
      "analytic, single, sham=5\n",
      "  94.866 ms (113 allocations: 738.97 KiB)\n"
     ]
    }
   ],
   "source": [
    "n=4096; FS=ones(n); FPS=ones(n,n); \n",
    "FPS32=ones(Float32,n,n); \n",
    "x0=ones(n); c=.5; hdata = heqinit(x0, c);\n",
    "println(\"analytic, double, Newton\"); \n",
    "@btime htest($x0, $FS, $FPS, hdata; \n",
    "    analytic=true, hsham=1);\n",
    "println(\"analytic, single, Newton\"); \n",
    "@btime htest($x0, $FS, $FPS32, hdata; \n",
    "    analytic=true, hsham=1);\n",
    "println(\"analytic, double, sham=5\"); \n",
    "@btime htest($x0, $FS, $FPS, hdata; \n",
    "    analytic=true);\n",
    "println(\"analytic, single, sham=5\"); \n",
    "@btime htest($x0, $FS, $FPS32, hdata; \n",
    "    analytic=true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For sufficiently large dimension, the linear algebra cost will dominate the computation. Time should increase by roughly a factor of 8 as the dimension doubles because our LU factorization takes $O(N^3)$ operations. We see that with this computation. Increasing the dimension from 1024 to 4096 did increase the runtimes by (roughly) a factor of 64 and we see that using single precision Jacobians cuts the run time in half. Note also that the single precision Shamanskii run is now five times faster than the double precision Newton computation. One of the projects at the end of this chapter challenges you to increase the dimension and compare the timings as you do that. Remember that we allocated storage for the Jacobian when we defined ```FPS```, so @btime is not measuring the allocation for the Jacobian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 2.8.3: More on the Two Point BVP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the Jacobian is sparse and you use the solvers from __SuiteSparse__ then you cannot use single precision. Even if the structure of the Jacobian allows you to use the LAPACK solvers or a special-purpose package, there is less benefit in using single precision for linear algebra than in the dense case. We will explore that for the boundary value problem, where we can use [BandedMatrices.jl](https://github.com/JuliaMatrices/BandedMatrices.jl) and __qr!__ for the linear solver. __qr!__ supports single precision linear algebra, so we will use that. __lu!__ does not, so the support is not consistent.\n",
    "\n",
    "The we will set up the problem for a very fine mesh, far finer than one needs to get a useful result, to illustrate the performance. The band solver takes $O(N)$ work, so we would expect the solve to be fast. To set things up we mimic __bvp_solve.jl__ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set it up\n",
    "    n=10^5;\n",
    "    bdata = bvpinit(n, Float64);\n",
    "#\n",
    "    U0 = zeros(2n);\n",
    "    FV = zeros(2n);\n",
    "# Banded matrix with the correct number of bands\n",
    "# Make double and single precision copies\n",
    "    FPV = BandedMatrix{Float64}(Zeros(2n, 2n), (2, 4));\n",
    "    FPV32 = BandedMatrix{Float32}(Zeros(2n, 2n), (2, 4));\n",
    "#\n",
    "# Build the initial iterate\n",
    "#\n",
    "    tv = bdata.tv;\n",
    "    sv = -.1 * tv .* tv;\n",
    "    view(U0,1:2:2n-1) .= exp.(-.1 .* tv .* tv);\n",
    "    view(U0,2:2:2n).= -.2 .* view(U0,1:2:2n-1) .* tv;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will begin by comparing the default ```sham=5``` with Newton's method. As we did with the H-equation, we will write a test function that uses the data we allocated above. We only use an analytic Jacobian for this and other examples with sparse Jacobians. The reader might want to look at the project in this chapter on [sparse differencing](#Sparse-Differencing)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "bvptest (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function bvptest(U0, FS, FPS, bdata; bsham=5, bfact=qr!)\n",
    "    FS .*= 0.0\n",
    "    FPS .*= 0.\n",
    "        bvpout = nsol(Fbvp!, U0, FS, FPS, Jbvp!; \n",
    "        atol=1.e-8, rtol = 1.e-8, sham=bsham,\n",
    "             pdata = bdata, jfact=bfact)\n",
    "    return bvpout\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will look at the convergence by running four cases with single and double precision and the \n",
    "kwarg sham set to 1 and 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGxCAYAAACa3EfLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0u0lEQVR4nO3dd1gUV/828Ht36b1JC4IYUVQUFCwRBWwgUWONJUWNGkuM3SexxR6NRo0tGk2ImFhiEnuJighiiwXBHiuoKIigUqXtzvuHP/Z1BRFw2YHl/lzXXo87c3bmOwtP5ubMmTkSQRAEEBEREVVxUrELICIiIlIHhhoiIiLSCgw1REREpBUYaoiIiEgrMNQQERGRVmCoISIiIq3AUENERERagaGGiIiItAJDDREREWkFhhqiEoSGhkIikRT7mjRpEuLj4yGRSBAaGip2qeVy5MgRDB48GO7u7jA2NsY777yDbt26ITo6utTbiIiIQMeOHWFrawsTExM0btwYK1asgFwuL9Xnp0+fDmdnZ+jo6MDCwqKcR1J1vO3vzPz587Fz584y7Wvx4sXl2tfbiIyMhEQiwd9//63xfVP1pSN2AURVwfr16+Hu7q6yzNHREXZ2djh16hTeffddkSp7O2vWrEFqairGjh2LBg0a4PHjx1iyZAlatmyJgwcPol27diV+/vDhwwgKCoKfnx9+/vlnGBsbY/fu3Rg7dixu376N5cuXl/j5Xbt24dtvv8W0adMQHBwMfX19dR5epeTg4PBWvzPz589H79690b17d/UWRqQFGGqISsHDwwM+Pj7FrmvZsqWGq3khOzsbRkZGb7WNH3/8Eba2tirLOnXqhDp16mD+/PlvDDWhoaHQ1dXF3r17YWxsDADo0KEDrl+/jtDQ0DeGmsuXLwMAxowZU6SOVz1//hyGhoZvOqRKSy6Xo6CgAPr6+qL9zhBpO15+InoLr7uUsGvXLjRu3Bj6+vqoXbs2li9fjlmzZkEikbzxswAgkUgwa9Ys5fvCz54/fx69e/eGpaWl8i99QRCwevVqeHl5wdDQEJaWlujduzfu3LnzxvqLCxImJiZo0KAB7t+//8bP6+rqQk9Pr0jYsLCwgIGBQYmfrVWrFqZPnw4AsLOzUznmWrVqoUuXLti+fTuaNGkCAwMDzJ49G8CLINStWzdYWlrCwMAAXl5e2LBhg8q2Cy99bN68GV9//TUcHBxgYmKCrl274tGjR8jIyMCwYcNgY2MDGxsbfPbZZ8jMzHzj8QYEBMDDwwPHjh1Dy5YtYWhoiHfeeQfffPONyuW2wp/tokWLMG/ePLi6ukJfXx8RERHF/twLf75XrlxB//79YW5uDjs7OwwePBhpaWnKdhKJBFlZWdiwYYPyMmhAQMAb61YoFPj222/h7OwMAwMD+Pj4IDw8vEi748ePo3379jA1NYWRkRFatWqFffv2FWlXmp9BcdLT0xEUFAQ7OzucOXPmje2Jyoo9NUSlUPhX9st0dIr/v8+BAwfQs2dP+Pn5YevWrSgoKMDixYvx6NGjt66jZ8+e6NevH0aMGIGsrCwAwPDhwxEaGooxY8Zg4cKFePLkCebMmYNWrVrhwoULsLOzK9M+0tLScP78+Tf20gDAiBEjsGXLFowZMwZTp06FkZER9uzZgx07dmDBggUlfnbHjh348ccfERISggMHDsDc3BxOTk7K9efPn8e1a9cwffp0uLq6wtjYGNevX0erVq1ga2uLFStWwNraGhs3bsSgQYPw6NEjfPXVVyr7mDp1Ktq2bYvQ0FDEx8dj0qRJ6N+/P3R0dODp6YktW7YgJiYGU6dOhampKVasWPHGY05KSkK/fv0wefJkzJkzB/v27cO8efPw9OlTrFq1SqXtihUrULduXSxevBhmZmZwc3Mrcdu9evVC3759MWTIEFy6dAlTpkwBAPz6668AgFOnTqFdu3Zo27YtvvnmGwCAmZnZG2tetWoVXFxcsGzZMigUCixatAjBwcE4evQo3nvvPQDA0aNH0bFjRzRu3BghISHQ19fH6tWr0bVrV2zZsgV9+/YFgDL/DAolJCTg/fffR15eHk6dOoXatWu/sW6iMhOI6LXWr18vACj2lZ+fL8TFxQkAhPXr1ys/06xZM6FmzZpCbm6ucllGRoZgbW0tvPx/ueI+WwiAMHPmTOX7mTNnCgCEGTNmqLQ7deqUAEBYsmSJyvL79+8LhoaGwldffVXmY/74448FHR0d4dy5c6Vqf+LECcHR0VH5vchkMmHRokWl+mzhcT1+/FhluYuLiyCTyYTr16+rLO/Xr5+gr68v3Lt3T2V5cHCwYGRkJDx79kwQBEGIiIgQAAhdu3ZVaTdu3DgBgDBmzBiV5d27dxesrKzeWK+/v78AQNi1a5fK8s8//1yQSqXC3bt3BUH4/z/bd999V8jLy1NpW9zPvfB7ePV7++KLLwQDAwNBoVAolxkbGwsDBw58Y60v78vR0VF4/vy5cnl6erpgZWUldOjQQbmsZcuWgq2trZCRkaFcVlBQIHh4eAhOTk7KGsr6M/jrr7+EmJgYwdHRUWjTpo2QmppaqtqJyoOXn4hK4bfffsPZs2dVXsX11GRlZeHcuXPo3r079PT0lMsLL328rV69eqm837t3LyQSCT755BMUFBQoX/b29vD09ERkZGSZtv/NN99g06ZN+OGHH+Dt7f3G9tHR0ejRowe8vb2xZ88eHDlyBFOmTMH06dMxd+7cMu37VY0bN0bdunVVlh05cgTt27dHzZo1VZYPGjQI2dnZOHXqlMryLl26qLyvX78+AKBz585Flj958qRUl6BMTU3xwQcfqCz76KOPoFAoEBUVpbL8gw8+gK6u7hu3+XL7lzVu3Bg5OTlITk4u9TaK07NnT5XLgaampujatSuioqIgl8uRlZWF06dPo3fv3jAxMVG2k8lk+PTTT5GQkIDr168DKPvP4ODBg2jTpg38/PwQFhYGKyurtzoWopLw8hNRKdSvX/+1A4Vf9vTpUwiCUOwln7JeBiqOg4ODyvtHjx69dn8AytTFP3v2bMybNw/ffvstvvzyy1J9ZtSoUbCzs8OOHTsgk8kAAG3btoVUKsWsWbPw8ccfl/syw6vHCgCpqanFLnd0dFSuf9mrJ9DCoPm65Tk5OSon9eIU913b29sXu//iai2JtbW1yvvCu8GeP39epu28qrC+V5fl5eUhMzMTGRkZEAShVN9tWX8GO3fuxPPnzzFy5MhqcXcbiYuhhkiNLC0tIZFIih0/k5SUpPK+8C/n3NxcleWvnhRe9vJAYwCwsbGBRCLBsWPHij1hlPYkMnv2bMyaNQuzZs3C1KlTS/UZAIiNjUX//v2VgaZQs2bNoFAocO3atXKHmlePFXhx0k9MTCyy/OHDhwBefB8VraSf7auhpLhjEMOrv3uFy/T09GBiYgIdHR1IpdJSfbdl/Rn88MMP2Lp1K4KDg7Fjxw4EBga+9fEQvQ4vPxGpkbGxMXx8fLBz507k5eUpl2dmZmLv3r0qbe3s7GBgYICLFy+qLN+1a1ep99elSxcIgoAHDx7Ax8enyKtRo0Zv3MbcuXMxa9YsTJ8+HTNnziz1voEXf52fO3euyIP2Ci9BvDzwVx3at2+PI0eOKE+ghX777TcYGRlp5FbpjIwM7N69W2XZ5s2bIZVK4efnV+H719fXL3PPzfbt25GTk6N8n5GRgT179qBNmzaQyWQwNjZGixYtsH37dpVtKxQKbNy4EU5OTspLgWX9GRgYGGD79u3o0qULPvjggzL9fhOVFXtqiNRszpw56Ny5M4KCgjB27FjI5XJ8//33MDExwZMnT5TtCsfC/Prrr3j33Xfh6emJM2fOYPPmzaXel6+vL4YNG4bPPvsM586dg5+fH4yNjZGYmIjjx4+jUaNGGDly5Gs/v2TJEsyYMQOdOnVC586d8e+//6qsf/kENWTIEGzYsAG3b9+Gi4sLAGD8+PEYM2YMunbtiuHDh8PIyAjh4eFYsmQJOnToAE9Pz1IfS2nMnDkTe/fuRdu2bTFjxgxYWVlh06ZN2LdvHxYtWgRzc3O17q841tbWGDlyJO7du4e6deti//79+PnnnzFy5Eg4OztX+P4bNWqEyMhI7NmzBw4ODjA1NUW9evVK/IxMJkPHjh0xYcIEKBQKLFy4EOnp6crb5AFgwYIF6NixI9q2bYtJkyZBT08Pq1evxuXLl7FlyxZlr1N5fga6urrYsmULhg4dit69e+O3335D//791fvFEIGhhkjtOnXqhG3btmHGjBno27cv7O3t8cUXX+Dhw4f4/fffVdouWbIEALBo0SJkZmaiXbt22Lt3L2rVqlXq/a1duxYtW7bE2rVrsXr1aigUCjg6OsLX1xfNmzcv8bN79uwB8OI29AMHDhRZLwiC8t9yuRxyuVxl2ejRo/HOO+/ghx9+wNChQ/H8+XPUqlULM2fOxPjx40t9DKVVr149nDx5ElOnTsWoUaPw/Plz1K9fH+vXr8egQYPUvr/i2Nvb48cff8SkSZNw6dIlWFlZYerUqSoBoSItX74co0aNQr9+/ZCdnQ1/f/83Dgj/8ssvkZOTgzFjxiA5ORkNGzbEvn374Ovrq2zj7++PI0eOYObMmRg0aBAUCgU8PT2xe/dulQHX5f0ZSKVShISEwNTUFJ988gmysrIwdOjQt/06iFRIhJf/C0VEFSI/Px9eXl545513cOjQIbHLoXIKCAhASkqK8knIRFS5sKeGqAIMGTIEHTt2hIODA5KSkvDTTz/h2rVrb5w2gIiIyo+hhqgCZGRkYNKkSXj8+DF0dXXRtGlT7N+/Hx06dBC7NCIircXLT0RERKQVeEs3ERERaQWGGiIiItIKDDVERESkFarNQGGFQoGHDx/C1NS00jy6nIiIiEomCAIyMjLg6OgIqbTkvphqE2oePnxYZFZZIiIiqhru37//xqlXqk2oMTU1BfDiSzEzMxO5GiIiIiqN9PR01KxZU3keL0m1CTWFl5zMzMwYaoiIiKqY0gwd4UBhIiIi0goMNURERKQVGGqIiIhIKzDUEBERkVZgqCEiIiKtwFBDREREWoGhhoiIiLQCQw0RERFpBYYaIiIi0goMNURERKQVGGrUICEhAREREUhISBC7FCIiomqLoeYthYSEwMXFBe3atYOLiwtCQkLELqlUGMSIiEjbSARBEMQuQhPS09Nhbm6OtLQ0tU1omZCQABcXFygUCpXltra20NXVhUQigVQqRY0aNWBmZgapVIrc3Fw8fvwYMpkM9erVg0wmg0wmQ2JiIrKzs9GiRQv4+PhAKpXi3r17CA8Ph46ODqRSqbKtTCaDjo4OatSoASsrK+jo6CA/Px/JycnQ0dGBp6cnDAwMoKOjg6SkJKSlpaFhw4bw9PSETCbDb7/9hgULFkAQBEgkEvz8888YMmSIWr4TIiIidSrL+Zuh5i1ERESgXbt2atmW2PT09GBubg4HBwf4+PigU6dOaNmyJZycnEo1MyoREVFFYKgphiZ7aqytrSGVSqFQKCAIAiwsLGBgYACFQoHs7Gw8ffoUgiDAysoKgiBAoVAgPT0d+fn5qFWrFlxdXaFQKJCQkIAbN26g8Ef08o9KEATo6elBJpNBEATI5XLk5eUBAKTS/39VsbAGPT09mJiYIDc3F1lZWaU+RolEAolEAgMDA3h6esLDwwONGjWCu7s7fHx8YGlp+TZfIRERUYkYaopREaEGeDGmZvjw4ZDL5ZDJZFi7dm2lvpTzuiDWsGFDpKSkIC0tDYIgIC8vD6X51TAyMoKvry9cXV1hbGyMixcvomnTpggMDIS/vz90dXUr6lCIiKgaYKgpRkWFGuBFULh16xbq1KkDJycntW67IpQmiAmCgPv37+PAgQPYu3cvbt++DblcjuTkZGRkZKCgoKDU+yvsJXJ3d0erVq3QpEkTtGrVCi4uLry0RUREJWKoKUZFhpqq6G2DWG5uLiIjI5GRkYGMjAzExcVh//79iI2NhVwuL/V2JBIJHBwc8M4776B27dqoUaMGGjdujO7du6NGjRpFar558ybc3NyqRHgEqmbNRESVCUNNMRhqNOvJkyfYv38/IiMjcfnyZdy/fx+CICA9PR3Z2dmlurTl4OAAV1dX6Ojo4PLly3jy5IlyXbNmzeDu7g6ZTAYXFxdYWlpCV1cXz58/x+PHj2FgYAAXFxfo6upCR0cHGRkZEAQB7u7usLe3h0wmQ2ZmJpKSkqCvrw8DAwPo6elBX18fhoaG0NfXh7GxMQwMDCCTyZR3sr08XulNQkJCMGzYMCgUCkilUqxbt65SX5osxCBGRJUJQ00xGGoql7i4OOzatQvnz59HRkYG4uPjkZSUhOTk5CLjfaqCwrvHCm/bf/bsWbHtpFIpjIyMYGxsDJlMBgBIS0uDTCZDjRo1YGBgAF1dXWRnZyM3Nxfu7u4ICAiAnp4eMjMzceDAAejp6cHAwEDlZWRkBDs7O9SqVQsGBgbQ19dHWloaTExM4OzsDAsLCxgbG8PQ0BBGRkbQ19cvtr6qGsSISHsx1BSDoabqkMvlePz4Me7fv4+4uDj88MMP+Pfff4u0KxyPo6OjA4lEoryTrPDyV2FoKFwOAPr6+tDX14dCocDz58/LdKlM2xgZGUFPTw9SqVSlF+xVNjY2cHFxgYWFBSwsLJCRkQFbW1v4+/vDw8MDzs7OMDExQX5+PszMzDg4nIjUiqGmGAw1VVdxd2zJZDLEx8e/1eURhUKBgoICyOVyyOVy5OfnIycnBzk5OcjNzVXeDp+fn4+srCw8efIEeXl5cHZ2Rn5+PgoKChAXF4cnT56gTp06cHV1hVwux4MHD7B7926sXbu2yD7btm0LGxsbGBgYIDc3V9lLlZ+fD2dnZ+Wt+YmJiUhLS0OtWrXg5eWFvLw8PHz4ECdOnIAgCEVewP+//f7V5Zqkq6uLJk2awNLSEubm5tizZ4+yZ8nExATm5uawsrKCra0tnJ2dUbduXdja2sLCwgLm5uYwMzODnZ3da3uSXoeXzIi0F0NNMRhqqraqdus8ULlqzsvLQ2ZmJlJSUpCamgodHR3Y2NggLy8Pz549Q3h4OBITE7Fq1aoin23QoAGMjIyQnp6OjIwM5SXCivpPh4mJCSwsLGBiYgK5XI6HDx+idu3a+Oabb2BhYQFLS0ssX74cVlZWcHJywpUrV/Dbb79BEAReMiPSQgw1xWCoqfqq2q3zQNWruaxBTKFQICUlBXfv3sWdO3eQm5sLCwsLPH36FAkJCdiyZYtycHhOTo6yh6six02poxePiCoPhppiMNQQlY4mg1hOTg6ePn2Kx48f48SJE0hMTMSjR4+QmpqK1NRUpKSk4MmTJzA0NISjoyOePn2K1NRUJCYmlrjdiIgIBAQEVGjtRKQZZTl/62ioJiKqIpycnDTWy2FgYAAHBwc4ODigcePGZfpsZmYm/v33XwQGBha5FPbHH3/A39+fD3ckqmZK/9ANIqJKxMTEBB06dMDPP/+svNOt0Nq1azFs2DBRBksTkXh4+YmIqrzCS2aGhobo0qULUlJSAADBwcHYs2dPkdBDRFVHWc7f7KkhoirPyckJAQEBaNGiBS5dugQHBwcAwD///IMePXqIXB0RaQpDDRFpFXt7e0RFRUFPTw8AsGfPHvzyyy8iV0VEmsBQQ0Rap06dOrh58yZMTEwAAJ9//jmWLVtWJafgIKLSY6ghIq3k7OyM27dvw8rKCgAwfvx41KpVC6mpqSJXRkQVhaGGiLSWra0tbt++jRo1agAA7t+/j27duolcFRFVFIYaItJqFhYWiI+Ph729PQDg5MmT2Lhxo8hVEVFFYKghIq1nZGSEu3fvomfPnhAEAQMGDMCaNWsQHx8vdmlEpEYMNURULejp6eGvv/7CqFGjIAgCvvjiC9SrVw/nz58XuzQiUpMqFWp69OgBS0tL9O7dW+xSiKgKkkqlWLlypXKSzry8PDRv3hxHjx4VuTIiUocqFWrGjBmD3377TewyiKgKk0gkWLduHfz9/QEAcrkc7dq1w/79+0WujIjeVpUKNW3btoWpqanYZRBRFSeVShEZGYkuXboAABQKBbp06YKtW7eKXBkRvQ21hZqoqCh07doVjo6OkEgk2LlzZ5E2q1evhqurKwwMDODt7Y1jx46pa/dERGW2Z88e9OvXDwAgCAL69euHn3/+WeSqiKi8dNS1oaysLHh6euKzzz5Dr169iqzfunUrxo0bh9WrV8PX1xdr165FcHAwrl69CmdnZwCAt7c3cnNzi3z20KFDcHR0LFM9ubm5KttKT08v4xERUXWwZcsWmJmZYd26dQCAYcOGIT09HRMnThS5MiIqK7WFmuDgYAQHB792/dKlSzFkyBAMHToUALBs2TIcPHgQa9aswYIFCwAA0dHR6ioHCxYswOzZs9W2PSLSXmvXroWZmRkWL14MAJg0aRKePXuGOXPmQCKRiFwdEZWWRsbU5OXlITo6GoGBgSrLAwMDcfLkyQrZ55QpU5CWlqZ83b9/v0L2Q0Ta4fvvv8e+ffugo/Pib7158+bhiy++gCAIIldGRKWlkVCTkpICuVwOOzs7leV2dnZISkoq9XaCgoLw4YcfYv/+/XBycsLZs2df21ZfXx9mZmYqLyKikrz//vs4fPiwcobvn376SXn7NxFVfhq9++nVblxBEMrUtXvw4EE8fvwY2dnZSEhIQLNmzdRdIhFVc/7+/vjjjz+U7zdv3ozr16+LWBERlZZGQo2NjQ1kMlmRXpnk5OQivTdERGLr0aOHMtjk5ubC09MTFy9eFLkqInoTjYQaPT09eHt7IywsTGV5WFgYWrVqpYkSiIjKpG/fvtizZw+kUilyc3Ph7e0Nf39/3klJVImpLdRkZmYiNjYWsbGxAIC4uDjExsbi3r17AIAJEybgl19+wa+//opr165h/PjxuHfvHkaMGKGuEoiI1KpLly4IDw+HTCZDQUEBoqKi0L9/f7HLIqLXkAhqGtofGRmJtm3bFlk+cOBAhIaGAnjx8L1FixYhMTERHh4e+OGHH+Dn56eO3b9Reno6zM3NkZaWxkHDRFQmp0+fRuvWrVFQUAAzMzOcPHkSDRs2FLssomqhLOdvtYWayo6hhojexuXLl9GlSxfcvXsX1tbWOHToEN59912Ym5uLXRqRVivL+btKzf1ERCQWDw8PnD9/Hs2aNUNqaipatWqFWrVq4cqVK2KXRkT/h6GGiKiUrKyscPjwYfj4+CA3NxfPnj1D8+bNS3xmFhFpDkMNEVEZmJmZYd++fbCwsAAAZGdno3Xr1jh69Ki4hRERQw0RUVnZ2tri7t27sLW1BfBiKpj27dtj//79IldGVL0x1BARlYOZmRni4uJQs2ZNAIBcLkeXLl2wdetWkSsjqr4YaoiIysnIyAi3bt2Cm5sbgBdTv/Tr1w+//PKLyJURVU8MNUREb0FPTw9Xr15F48aNlcs+//xzLFmyRMSqiKonhhoioreko6ODmJgYjB07Vrls0qRJ+Oabb1BNHgVGVCkw1BARqYFUKsWyZcswb9485bJ58+Zh/PjxIlZFVL0w1BARqdG0adMwffp05futW7eioKBAxIqIqg+GGiIiNZs7dy4GDBgAAEhKSkK9evWQk5MjclVE2o+hhoioAmzYsAGjRo0CANy5cwe1a9fGyJEjkZWVJXJlRNqLoYaIqIKsWrUKkydPBgAkJibip59+wkcffSRyVUTai6GGiKgCLViwQGXw8MWLF5GZmYmEhAREREQgISFBxOqItAtDDRFRBZs2bRpWrFgBXV1dxMfHw9PTEy4uLmjXrh1cXFwQEhIidolEWkEiVJOHKKSnp8Pc3BxpaWkwMzMTuxwiqobOnDmDwMBApKWlqSyXSqW4ePEiGjZsKFJlRJVXWc7f7KkhItKQ5s2bY/HixUWWKxQKNG7cGM2bN8eUKVMQHh6O58+fi1AhUdXGnhoiIg1KSEiAi4sLFApFie309fXh6+uLDh06oH379vD29oZMJtNQlUSVB3tqiIgqKScnJ6xbtw5S6ev/82tkZITc3FwcOXIEU6dORYsWLWBtbY01a9ZosFKiqoehhohIw4YMGYK7d+8iIiICO3fuhKenp3KdTCbDxIkTcfXqVaxcuRJt27aFTCZDWloa8vPzle2io6Px6aef4q+//hLjEIgqJV5+IiISmVwuR0hICKZOnYrU1FQAwIcffojFixcjNDQUM2fOVLatV68eunbtiidPnuDXX39F7969lcFGEATs2bMHrVu3hpWVlSjHQqRuZTl/M9QQEVUST548wYwZM7BmzRooFAoYGhriiy++gIODAw4cOICjR4+q9NYYGBhg0aJFGD16NADg+vXrcHd3h0QiQdOmTZXjcXx9fWFkZCTWYRG9FYaaYjDUEFFVceHCBYwePRrHjh0DANSuXRvLli2Dn58fDh06hD179mD//v0oKCjA48ePoaurCwCYOHEiNm3ahEePHqlsT09PD61atVKGHB8fH+jo6Gj8uIjKg6GmGAw1RFSVCIKALVu24H//+x8ePnwIAAgODsayZctQt25dyOVy3Lp1C/Xq1VO2f/fddxEXFwcAcHR0hJWVFZKTk5GcnKyybTMzM/j7+ytDDp+PQ5UZQ00xGGqIqCrKzMzEvHnzsHTpUuTn50NXVxcTJkzA9OnTYWJiomxXUFCAlStXYs+ePYiKioJcLleus7S0hI+PD0xNTREREYGnT58q17m5ueHGjRvK96mpqbC2ttbMwRGVAm/pJiLSEiYmJvjuu+9w+fJlBAcHIz8/HwsXLkS9evWwefNmFP5dqqOjg/Hjx+PIkSN4/PgxNm/ejP79+8PCwgJPnz5FnTp1sG3bNjx+/BinTp1C9+7d0aZNG3Tt2lW5r/z8fNSqVQtubm7K3iGiqoQ9NUREVYQgCNi7dy/GjRuHO3fuAAD8/PywYsUKldvCX5afn48TJ07A1tYWDRo0AAAcO3YMfn5+AIBGjRqha9eu+OCDD6Crq4tmzZrBysoKjx49Uj5LZ/bs2cjKykL79u3RunVrGBsba+BoiV7g5adiMNQQkbbIycnBkiVL8O233+L58+eQSqUYOXIk5syZU6pbuY8ePYpvvvkGJ06cUHmysZ2dHTp27IguXbqgb9++AF4EKScnJ2XPjZmZGf766y8EBgZWzMERvYKhphgMNUSkbe7du4dJkyYpn1NjbW2N+fPnY8iQIaWaUiE1NRX//PMP9uzZgwMHDiA9PR0AcPr0aTRv3hwAcOvWLfzzzz+IiYlBWFgYEhISYGBggH/++QcBAQEVdmxEhRhqisFQQ0Ta6siRIxgzZgyuXLkCAPD29sbKlSvx3nvvlXobeXl5OHbsGA4fPoxvv/1Weelp6NChCAkJQZMmTfD+++/j2LFjiIqKgrGxMQ4ePAhfX98KOSaiQgw1xWCoISJtlp+fj9WrV2PmzJlIS0sDAAwcOBDfffcd7O3ty73dnj17YufOnXj5VOHh4YHLly/D1NQUhw8fVvbqEFUE3v1ERFTN6OrqYuzYsbhx4wYGDx4MANiwYQPq1q2rvB28PLZv345Hjx4hNDQUnTp1AgDcuHED3t7eyMjIQFBQEGJiYtR2HERvgz01RERa6PTp0xg9ejTOnj0LAKhfvz5WrFiBDh06lHubCoUCPXv2xK5du2BnZwdnZ2ecPXsW1tbWiIiIQKNGjdRVPpESe2qIiKq5Fi1a4N9//8Uvv/yCGjVq4Nq1a+jYsSN69eqFu3fvlmubUqkUv//+Oxo2bAh/f3/s3r0bzZs3R2pqKjp06ID//vtPzUdBVDbsqSEi0nLPnj3DzJkz8eOPP0Iul8PAwACTJ0/GV199BUNDwzJvLzU1FVZWVpBIJHj69Cnat2+PmJgYNGjQABcvXizVnVdEpcWeGiIiUrKwsMDy5csRExODgIAA5OTkYNasWWjQoAF27NiBsv5ta21tDYlEotz2/Pnz0bZtW/z+++8MNCQqhhoiomqiUaNGOHLkCLZu3QonJyfEx8ejZ8+e6NSpU7kuHeXn56Nv377o3LkzvvrqKzRt2lS5rppcBKBKhqGGiKgakUgk6NOnD/777z9MmzYNenp6OHToEBo1aoT//e9/ygfwlYaOjg5MTU2hUCjQr18/5cSYJ0+eRPPmzZGYmFhRh0FULIYaIqJqyNjYGPPmzcOVK1fQpUsXFBQUYPHixahXrx5+//33YntaEhISEBERgYSEBAAvAtLq1avRqlUrpKWloVu3bnj69CmGDx+Oc+fO4ZtvvtH0YVE1x1BDRFSN1alTB3v27MG+fftQp04dJCUlYcCAAWjdurXK82dCQkLg4uKCdu3awcXFBSEhIQAAfX19bNu2DU5OTvjvv//w6aefYvv27RgwYABWrFgh1mFRNcW7n4iICACQm5uLH374AfPmzUNWVhYkEgmGDx+OL774Al5eXiqTX8pkMsTHx8PJyQkAEB0djdatWyMnJweTJ0/GggULVLadk5MDAwMDjR4PaQetvPspIyMDzZo1g5eXFxo1aoSff/5Z7JKIiLSKvr4+Jk+ejP/++w/9+/eHIAj46aef4OvrqxJoAEAul+PWrVvK997e3vj1118BAEuXLkVcXByAFwOGZ8yYgYCAgDKN1yEqjyrTUyOXy5GbmwsjIyNkZ2fDw8ND+STL0mBPDRFR2URFRWH06NG4ePFikXWv9tQU+u6779CmTRvlRJeJiYnw8PDAkydP0Lp1axw4cADGxsYaqZ+0g1b21MhkMhgZGQF40Y0pl8t5yyARUQXy8/NDdHQ0Vq1apfKQPqlUirVr1xYJNAAwefJklZm7HRwccOjQIZibm+P48eP44IMP8Pz5c43UT9WP2kJNVFQUunbtCkdHR0gkEuzcubNIm9WrV8PV1RUGBgbw9vbGsWPHyrSPZ8+ewdPTE05OTvjqq69gY2OjpuqJiKg4Ojo6GDVqFO7evYvAwEAAL553UzhpZkkuXbqEjz76CB4eHjhw4ABMTExw5MgR9OjRA7m5uRVdOlVDags1WVlZ8PT0xKpVq4pdv3XrVowbNw7Tpk1DTEwM2rRpg+DgYNy7d0/ZxtvbGx4eHkVeDx8+BPDiyZUXLlxAXFwcNm/ejEePHqmrfCIiKkGNGjXw+++/w9jYGBcuXMDu3btLbJ+bm4v3338fW7ZswRdffIEWLVpg//79MDIywsGDB/Hhhx8iLy9PQ9VTtSFUAADCjh07VJY1b95cGDFihMoyd3d3YfLkyeXax4gRI4Q///zztetzcnKEtLQ05ev+/fsCACEtLa1c+yMiIkGYOnWqAEBo1KiRIJfLS2x78OBBQSqVCgCE5cuXC4IgCOHh4YKBgYEAQOjVq5eQn5+vibKpCktLSyv1+VsjY2ry8vIQHR2t7LosFBgYiJMnT5ZqG48ePVKOnE9PT0dUVBTq1av32vYLFiyAubm58lWzZs3yHwAREQEAJk2aBHNzc1y6dAl//vlniW0DAwPx/fffAwAmTJiA8PBwtGvXDjt27ICenh62bduGAQMGQC6Xa6J0qgY0EmpSUlIgl8thZ2enstzOzg5JSUml2kZCQgL8/Pzg6emJ1q1b48svv0Tjxo1f237KlClIS0tTvu7fv/9Wx0BERIClpSUmTZoEAJg5cyYKCgpKbD9+/Hh8+umnkMvl+PDDD3H79m106tQJf//9N3R0dLBlyxYMHTq0yC3jROWh0bufCmd1LSQIQpFlr+Pt7Y3Y2FhcuHABFy9exMiRI0tsr6+vDzMzM5UXERG9vbFjx8LGxgY3btzA77//XmJbiUSCdevWoXnz5nj69Cm6deuGjIwMdO3aFX/88QdkMhlCQ0PxxRdf8I5WemsaCTU2NjaQyWRFemWSk5OL9N4QEVHlZmpqismTJwMAZs+e/cY7mQwMDLBjxw44ODjAwcFB2bvTq1cv/Pbbb8rgc/r06QqvnbSbRkKNnp4evL29ERYWprI8LCwMrVq10kQJRESkRl988QUcHR1x9+5d/PLLL29s7+joiOPHj+Off/6BpaWlcvlHH32EX3/9FRs3bkTLli0rsmSqBtQWajIzMxEbG4vY2FgAQFxcHGJjY5W3bE+YMAG//PILfv31V1y7dg3jx4/HvXv3MGLECHWVQEREGmJoaIjp06cDAObNm4fs7Ow3fqZ27drQ0dFRvr9z5w4AYNCgQfjoo4+Uy7OystRcLVUXags1586dQ5MmTdCkSRMAL0JMkyZNMGPGDABA3759sWzZMsyZMwdeXl6IiorC/v374eLioq4SiIhIg4YMGYJatWohKSkJq1evLvXn8vPzMWLECDRq1Ej5h3Chhw8fwsfHB99++62aq6XqoMrM/fS2OPcTEZH6hYaG4rPPPoO1tTXu3LlTqv++FhQUoHPnzjh06BCcnZ1x7tw51KhRAwDw888/Y9iwYahZsyYuXrwICwuLCj4Cquy0cu4nIiKqfD755BPUq1cPqampWL58eak+o6Ojgz/++ANubm64d+8eevfurXy68Oeff46VK1ciKiqKgYbKjKGGiIjKTUdHB7NnzwYALF68GE+ePCnV5ywtLbFr1y6YmpoiKioKY8eOVa778ssvUatWLeX7Bw8eqLVm0l4MNURE9FY+/PBDeHp6Ij09XfkE4dKoX78+tmzZAolEgp9++gk//fRTkTa7du3Cu+++iw0bNqizZNJSDDVERPRWpFIp5s6dCwBYsWJFmSYb7ty5M+bPnw8AGDdunHIC40JHjx5Fbm4uBg8ejC1btqivaNJKDDVERPTWunTpgubNmyM7OxsLFiwo02e//vprDB8+HNu3b4ejo6PKuiVLlmDYsGFQKBT49NNPsW3bNnWWTVqGdz8REZFaHD58GB07doSenh5u3bqltomEFQoFhgwZgtDQUOjo6GD79u3o2rWrWrZNlR/vfiIiIo1r3749AgICkJeXh3nz5pV7O3fu3MHkyZOVk1xKpVL88ssv6N+/PwoKCtC7d28cOnRIXWWTFmGoISIitZBIJMow8+uvv+L27dtl3kZ2djZat26NhQsXqgQjmUyG3377Db169UJeXh66deuGiIgItdVO2oGhhoiI1MbX1xfBwcEoKCjArFmzyvx5IyMjZZiZOXMmduzYoVyno6ODzZs3o2vXrsjJyUGXLl1w/PhxdZVOWoChhoiI1KrwTqhNmzbh6tWrZf784MGDMXr0aADAp59+ikuXLinX6enp4a+//kJQUBCys7Px/vvvc3ZvUmKoISIitfL29kbPnj0hCIJy/r+yWrJkCdq1a4esrCx069YNKSkpynX6+vrYsWMH2rZti4yMDAQFBeH8+fPqKp+qMIYaIiJSuzlz5kAikWDbtm3lChy6urr4888/Ubt2bcTFxaFPnz7Iz89Xrjc0NMSePXvQunVrpKWl4c8//1Rn+VRFMdQQEZHaNWzYEB9//DEA4JtvvinXNqytrbFr1y6YmJggMzMTaWlpKuuNjY2xb98+LF++vMzPxiHtxOfUEBFRhbh16xbc3d0hl8tx4sQJtGrVqlzbOXPmDBo1agRDQ8M3ts3Pz0dSUpLanpFD4uNzaoiISHR16tTB4MGDAQDTpk1Def+Gbt68uUqgefr0abHtcnNz0bt3b7Rq1QpxcXHl2hdVbQw1RERUYb755hvo6ekhMjISR44ceattKRQKTJ8+He7u7rh3716R9RkZGbhx4wZSUlLK9YwcqvoYaoiIqMLUrFkTI0aMAPB2vTUAkJOTg3379iE5ORndu3dHdna2ynobGxuEh4fj4MGD6NChw1vVTVUTQw0REVWoqVOnwsjICKdPn8bevXvLvR0jIyPs3LkTNWrUQExMDAYPHlwkJDk6OsLPz0/5/ubNm0hOTi73PqlqYaghIqIKZWdnhzFjxgB4cTmqcE6n8nBxccG2bdugo6ODrVu34rvvvntt26tXr8LPzw8dOnRAampqufdJVQdDDRERVbj//e9/MDMzw4ULF/D333+/1bbatGmDVatWAXhxSWvPnj3FttPV1YVEIsGlS5eKPOeGtBNDDRERVTgrKytMnDgRADBjxgwUFBS81faGDx+OkSNHQhAEfPzxx8VeYnJzc8PBgwdhbGyMI0eOKPdP2ouhhoiINGLcuHGwtrbG9evXsWnTprfe3vLlyxEUFITVq1fD1ta22DaNGjXC77//DgBYuXIlQkJC3nq/VHkx1BARkUaYmZnh66+/BgDMmjULeXl5b7U9XV1d/PPPP/jkk09KbNejRw/Mnj0bADBy5EicOHHirfZLlRdDDRERacyoUaNgb2+P+Ph4tfSaSCQS5b8fPXqE1atXF9tu+vTp6N27N/Lz89GzZ0/cv3//rfdNlQ9DDRERaYyRkRGmT58OAJg3bx6eP3+ulu1mZGSgefPmGDVqFEJDQ4usl0qlCA0Nhaen52ufc0NVH0MNERFp1NChQ+Hs7IyHDx9izZo1atmmqakpBg0aBODFIOJTp04VaWNsbIxdu3bBxsYG58+fL/Y5N1S1MdQQEZFG6evrY+bMmQCABQsWICMjQy3bnTlzJnr06IG8vDz07NkTDx48KNKmLM+5oaqHoYaIiDRuwIABcHNzQ0pKClasWKGWbUqlUvz222/w8PBAUlISevToUezlLT8/P5Xn3Bw+fFgt+yfxMdQQEZHG6ejoKO9I+v77718783ZZmZiYYNeuXbCyssLZs2cxbNiwYi8xFT7npmfPnmjZsqVa9k3iY6ghIiJR9O3bF40aNUJaWhoWL16stu3Wrl0bf/31F2QyGc6ePYsnT54U227FihX4888/YWJiorZ9k7gYaoiISBRSqRRz584F8OJBeuqceLJdu3bYsWMHTp8+DWtr62Lb6OjoQCp9cRoUBAGbN29+6ycdk7gYaoiISDQffPABmjVrhqysLLUP2u3atSvMzc2V73Nzc1/bdtiwYfj444/x1VdfqbUG0iyGGiIiEo1EIsG8efMAAKtXr0ZCQoLa9yEIApYvX47GjRu/drbuTp06QVdXF/Xq1VP7/klzGGqIiEhUHTt2hJ+fH3Jzc/Htt9+qffuZmZlYvnw5bty4gd69exc7PUOvXr1w+/ZtDB8+XO37J81hqCEiIlG93Fvzyy+/4M6dO2rdvqmpKXbv3g0TExNERkZi9OjRxd4RVbNmTeW/U1NT8fDhQ7XWQRWPoYaIiETXpk0bBAUFoaCgQHmrtzp5eHhgy5YtkEgkWLdunfI5NcW5du0aWrRogW7duqltGgfSDIYaIiKqFAp7azZu3Ihr166pfftdunTBokWLAADjxo3DoUOHim1nYGCAZ8+e4dy5cxg6dCinUqhCGGqIiKhS8PHxQffu3aFQKJTTKKjbxIkTMXDgQCgUCvTp0weJiYlF2ri6uiqfc7N582Z8//33FVILqZ9EqCYRND09Hebm5khLS4OZmZnY5RARUTEuX76Mxo0bQxAEnD9/Hk2aNFH7PnJzcxEYGIgePXpg7NixkEgkxbb78ccf8eWXX0IikWDv3r14//331V4LvVlZzt8MNUREVKl8/PHH2Lx5M7p06YI9e/ZUyD7kcjlkMlmJbQRBwPDhw/Hzzz/DzMwMp0+fhru7e4XUQ69XlvM3Lz8REVGlMmvWLMhkMuzduxf//vtvhezj5UCTlpaGtWvXFmkjkUiwatUqtG7dGunp6ejWrRuePXtWIfWQejDUEBFRpeLm5oZBgwYBeDGLdkXKzc2Fr68vRowYgTVr1hRZr6enh23btsHZ2Rk3btxA//79IZfLK7QmKr8qFWp0dHTg5eUFLy8vDB06VOxyiIiogsyYMQN6eno4cuQIjhw5UmH70dfXxyeffAIAGD16dLH7srW1xa5du2BoaIgDBw5g8uTJFVYPvZ0qNabGxsYGKSkp5fosx9QQEVUto0ePxqpVq/Dee+/hxIkTrx3Q+7YEQcCAAQOwceNGWFpa4vTp03BzcyvS7s8//0Tfvn0BABs2bMCAAQMqpB5SxTE1RERU5U2bNg2GhoY4deoU9u/fX2H7kUgk+Pnnn9GyZUs8ffoUXbt2LXbsTJ8+fZSXw4YNG4azZ89WWE1UPmoLNVFRUejatSscHR0hkUiwc+fOIm1Wr14NV1dXGBgYwNvbG8eOHSvTPtLT0+Ht7Y3WrVvj6NGjaqqciIgqI3t7e4wePRoA8M0330ChUFTYvgwMDLBjxw44OTnh+vXr6NevHwoKCoq0mzNnDrp16wZfX1/Url27wuqh8lFbqMnKyoKnp+drHz29detWjBs3DtOmTUNMTAzatGmD4OBg3Lt3T9nG29sbHh4eRV6F82/Ex8cjOjoaP/30EwYMGID09PTX1pObm4v09HSVFxERVS1fffUVTE1NERMTg+3bt1fovuzt7bF7924YGRnhypUrePDgQZE2UqkUmzZtwoEDB2BtbV2h9VA5CBUAgLBjxw6VZc2bNxdGjBihsszd3V2YPHlyufbRqVMn4ezZs69dP3PmTAFAkVdaWlq59kdEROIo/O95/fr1hYKCggrfX1hYmJCYmFjq9uHh4YJCoajAiqq3tLS0Up+/NTKmJi8vD9HR0QgMDFRZHhgYiJMnT5ZqG0+fPkVubi4AICEhAVevXi2x62/KlClIS0tTvu7fv1/+AyAiItGMHz8eVlZWuHbtGjZv3lymzyYkJCAiIgIJCQml/kyHDh1gb2+vfF/SpJbjx49H+/btsXTp0jLVRRVDI6EmJSUFcrkcdnZ2Ksvt7OyQlJRUqm1cu3YNPj4+8PT0RJcuXbB8+XJYWVm9tr2+vj7MzMxUXkREVPWYm5vjq6++AvDiwXz5+fml+lxISAhcXFzQrl07uLi4ICQkpMz73rRpE9zc3HDnzp1i17u6ugJ48cc7iU+jdz+9ejueIAilvkWvVatWuHTpEi5cuIDY2Fh07969AiokIqLK6Msvv4SdnR3u3LmDX3/99Y3tExISMGzYMOXgYoVCgeHDh5epx6agoADLly/HgwcP0LVr12LHZo4ePRrnz5/HlClTSn8wVGE0EmpsbGwgk8mK9MokJycX6b0hIiJ6lbGxsfJ26rlz5yInJ6fE9jdv3ixyt5RcLsetW7dKvU8dHR3s3LkTjo6OuHr1arFPE5ZIJCqTbmZlZfHGFBFpJNTo6enB29sbYWFhKsvDwsLQqlUrTZRARERV3LBhw1CzZk08ePAAP/30U4lt3dzcIJWqnuJkMhnq1KlTpn06Ojpi165dMDAwwP79+/H111+/tu3du3fh6+vLqRREpLZQk5mZidjYWMTGxgIA4uLiEBsbq7xle8KECfjll1/w66+/4tq1axg/fjzu3buHESNGqKsEIiLSYvr6+pgxYwYAYP78+cjMzHxtWycnJ6xbt045caVMJsPatWvh5ORU5v36+PggNDQUALBkyRKsX7++2HYpKSm4fv069u/fX+FzVtFrqOuWq4iIiGJvoR44cKCyzY8//ii4uLgIenp6QtOmTYWjR4+qa/dvVJZbwoiIqHLKy8sT6tSpIwAQ5s+f/8b29+/fFyIiIoT79++/9b5nzJghABB0dXWFS5cuFdtm8+bNyvPfpk2b3nqfVLbzd5Wa++ltcO4nIiLtsGnTJnzyySewsLBAXFwcLCwsNLJfhUKBfv36oW7dupgzZ06Ry1uFpkyZgu+++w4GBgY4duwYfHx8NFKftirL+ZuhhoiIqhS5XA5PT09cuXIF06dPx9y5czW2b4VC8dowU0gul6Nbt27Yt28f3nnnHZw9exYODg4aqlD7cEJLIiLSWjKZTBlkli1bhsePH2ts3y8HmpycHCxfvrzIXVYymQybNm2Cu7s7Hjx4gF69eikfHksVi6GGiIiqnO7du8Pb2xuZmZlYuHChxvcvCAI6deqknNPwVebm5ti9ezcsLCxw6tQpjBw5EtXkwoioGGqIiKjKkUgkmDdvHgDgxx9/VE58rMn9Dxs2DADw3XffYePGjUXauLm5YevWrZBKpVi/fj1WrFih0RqrI4YaIiKqkoKCgtC6dWvk5OTg22+/1fj+P/roI0ydOhUAMHToUPz7779F2gQGBmLx4sUAgIkTJ+Lw4cMarbG64UBhIiKqso4ePYqAgADo6urixo0bqFWrlkb3r1Ao0KtXL+zcuRN2dnY4c+YMnJ2dVdoIgoDPPvsMGzZsgKWlJc6cOVPmhwBWZxwoTERE1YK/vz86duyI/Px8zJ49W+P7l0ql+P333+Hp6YlHjx6hW7duyMrKUmkjkUjw008/oUWLFnB1dYW+vr7G66wu2FNDRERV2pkzZ9CiRQtIpVJcvXoV9erV03gN9+7dQ7NmzVBQUIDw8HB4eXkVafP48WMYGxvDyMhI4/VVZeypISKiaqN58+b44IMPoFAoMHPmTFFqcHZ2xt69e3H27NliAw0A1KhRQyXQXL16VUPVVR8MNUREVOXNnTsXEokEW7duxYULF0SpoVmzZqhdu7byfXZ2drHtBEHArFmz4OHhgT/++ENT5VULDDVERFTlNW7cGH379gUA5aSXYtq/fz9q166NM2fOFFknkUiQlZUFQRBEC2DaimNqiIhIK1y/fh0NGjSAQqHAv//+ixYtWohWS/fu3bFr1y44ODjg7NmzeOedd1TWy+VyHDp0CMHBwSJVWHVwTA0REVU79erVw8CBAwEA06dPF7WW3377DQ0bNkRiYiK6detW5FKUTCZTCTR5eXnIy8vTdJlah6GGiIi0xowZM6Crq4vDhw8jMjJStDrMzMywZ88e2NjYIDo6Gp999tlrp0l4/PgxOnbsiFGjRnEqhbfEUENERFqjVq1a+PzzzwG86K0RMyS4urpi27Zt0NXVxZ9//ok5c+YU2y42NhbHjh3DL7/8gh9//FHDVWoXhhoiItIq06ZNg4GBAU6cOIEDBw6IWoufnx/WrFkDAJg1a1ax0yR07NgRixYtAgCMGzcOR44c0WiN2oShhoiItIqjoyO+/PJLAOL31gDAkCFDMH78eAwfPhz+/v7Ftpk4cSI++eQTyOVyfPjhh7hz546Gq9QOvPuJiIi0TkpKClxdXZGZmYlt27ahZ8+eotajUCggkUggkUhe2+b58+fw9/fH2bNn4eHhgZMnT8LU1FSDVVZOvPuJiIiqNRsbG4wfPx4A8M0330Aul4taj1QqVQaagoICLFmyBM+fP1dpY2hoiB07dsDBwQGXL1/G4MGDRe9lqmoYaoiISCtNmDABlpaWuHr1Knbu3Cl2OUoDBgzApEmTMGTIkCKh5Z133sH27duhq6uLv//+G0uXLhWpyqqJoYaIiLSShYUFhgwZAgDYu3evyNX8f59//jl0dHSwZcsWzJ8/v8j6li1bYvny5QCAr7/+WtRb06sahhoiItJagYGBAIDDhw9Xmks5bdu2xapVqwC8GMi8Y8eOIm1GjBiBTz/9FHK5HH379sWDBw80XWaVxFBDRERaq3Xr1tDX10dCQgKuX78udjlKw4cPx+jRowEAn3zyCWJjY1XWSyQS/PTTT2jcuDGSk5Px4Ycf8onDpcBQQ0REWsvQ0BCtW7cGgGKfESOmpUuXomPHjsjOzsYHH3yApKQklfVGRkbYvn07zM3N8fjx4yLrqSiGGiIi0modO3YEAISFhYlciSodHR1s3boVdevWxZMnT3Dt2rUibd59910cOHAA586dg7OzswhVVi06YhdARERUkTp06AAAiIiIQEFBAXR0Ks+pz9LSEnv27MHz58/h6elZbJuWLVuqvM/Ly4Oenp4myqty2FNDRERarUmTJrCyskJGRgbOnDkjdjlF1K1bVyXQZGVlFdtOEASsWLECXl5eePbsmYaqq1oYaoiISKtJpVK0b98eQOUbV/OqkydPok6dOti9e3eRdRkZGVi6dCmuXbuG9evXi1Bd5cdQQ0REWq9wXE1lDzV//PEHkpKS8OmnnyI1NVVlnZmZGf7++28sW7YM48aNE6fASo5zPxERkdaLi4tD7dq1oaOjgydPnlTaOZXy8/PRu3dvfPbZZ+jevbvY5VQKnPuJiIjoJa6urqhduzYKCgoQFRUldjmvpauri127dpUq0KSnp2Po0KG4d+9exRdWRTDUEBFRtVBZb+0ur+HDhyMkJAS9evVCTk6O2OVUCgw1RERULRTe2l3Zx9WU1oIFC2BlZYVz585h7NixYpdTKTDUEBFRtdCuXTtIJBJcuXIFDx8+FLuct1arVi1s2rQJEokE69at4x1RYKghIqJqwsrKCt7e3gCA8PBwkatRj06dOmH27NkAgJEjR+L8+fMiVyQuhhoiIqo2tG1cDQBMmzYNnTt3Rm5uLnr16oUnT56IXZJoGGqIiKjaeHlcjbY80UQqleL3339H7dq1ER8fj08++QQKhULsskTBUENERNVGq1atYGhoiMTERFy9elXsctTG0tIS27Ztg4GBAf755x/MnTtX7JJEwVBDRETVhoGBAdq0aQNAe+6CKuTl5YW1a9cCAGbPno39+/eLXJHmMdQQEVG1UlWmTCiPAQMGYOTIkRAEAR9//DHu3LkjdkkaxVBDRETVSuG4msjISOTn54tcjfr98MMPaNGiBZ49e4Yvv/xS7HI0qsqEmuvXr8PLy0v5MjQ0xM6dO8Uui4iIqpjGjRujRo0ayMzMxOnTp8UuR+309fXx999/o0ePHggJCRG7HI2qMqGmXr16iI2NRWxsLI4fPw5jY2NlFyIREVFpSaVStG/fHoB23dr9MicnJ2zfvh0ODg5il6JRVSbUvGz37t1o3749jI2NxS6FiIiqIG2bMuFN/vjjD5w5c0bsMiqc2kJNVFQUunbtCkdHR0gkkmIvDa1evRqurq4wMDCAt7c3jh07Vq59/fnnn+jbt+9bVkxERNVVYU//6dOnkZaWJnI1FSs0NBT9+/dH7969kZqaKnY5FUptoSYrKwuenp5YtWpVseu3bt2KcePGYdq0aYiJiUGbNm0QHBysMmW6t7c3PDw8irxenqMjPT0dJ06cwPvvv6+u0omIqJpxdnaGm5sb5HI5jh49KnY5FapHjx6oV68ePvnkE1hYWIhdToXSUdeGgoODERwc/Nr1S5cuxZAhQzB06FAAwLJly3Dw4EGsWbMGCxYsAABER0e/cT+7du1CUFAQDAwMSmyXm5uL3Nxc5fv09PTSHAYREVUTHTt2xM2bNxEWFoYPPvhA7HIqjLm5OaKjo6vFkA2NjKnJy8tDdHQ0AgMDVZYHBgbi5MmTZdpWaS89LViwAObm5spXzZo1y7QfIiLSbtVpXM3LgSY3N1drJ77USKhJSUmBXC6HnZ2dynI7OzskJSWVejtpaWk4c+YMgoKC3th2ypQpSEtLU77u379f5rqJiEh7tW3bFlKpFP/99x8SEhLELkcjHj9+DD8/P7Rt2xY3b94Uuxy10+jdTxKJROW9IAhFlpXE3Nwcjx49gp6e3hvb6uvrw8zMTOVFRERUyMLCAs2aNQNQPXprgBfHrKuri/T0dPTs2RNZWVlil6RWGgk1NjY2kMlkRXplkpOTi/TeEBERaUrhXVDa+ryaV+nq6uLPP/+EnZ0dLl++jGHDhmnNbOWAhkKNnp4evL29i/zShIWFoVWrVpoogYiIqIiXx9Vo08m9JI6Ojvjzzz8hk8mwefPm1961XBWpLdRkZmYqn/gLAHFxcYiNjVXesj1hwgT88ssv+PXXX3Ht2jWMHz8e9+7dw4gRI9RVAhERUZm0bNkSRkZGSE5OxuXLl8UuR2P8/Pzw/fffA3hxfj5x4oTIFamH2kLNuXPn0KRJEzRp0gTAiy+pSZMmmDFjBgCgb9++WLZsGebMmQMvLy9ERUVh//79cHFxUVcJREREZaKvrw9/f38A1ecSVKFx48ahT58+KCgoQJ8+fcp0405lJRGqSX9beno6zM3NkZaWxkHDRESktHTpUkycOBHBwcHYv3+/2OVoVGZmJpo3b45r167Bz88Phw8fhq6urthlqSjL+btKzv1ERESkLoWDhY8ePary0NbqwMTEBNu3b4epqSmioqIwZcoUsUt6Kww1RERUrXl4eMDOzg7Z2dn4999/xS5H49zd3bF+/XoAwJIlS/DXX3+JXFH5MdQQEVG1JpFIlHdBVbdxNYV69eqF//3vfwCAzz77DNeuXRO5ovJhqCEiomqvOk2Z8Drz589HQEAAateuXenG1ZQWBwoTEVG1l5CQgJo1a0IqlSIlJQWWlpZilySK1NRUGBgYVKrJLzlQmIiIqAycnJzg7u4OhUKByMhIscsRjbW1tUqgqWpzYjHUEBERofJMmZCQkICIiAhRA4UgCFiwYAHeffddHD16VLQ6yoqhhoiICJVjXE1ISAhcXFzQrl07uLi4ICQkRLRarl69iry8POzbt0+0GsqKY2qIiIjw4jxhZWUFuVyO+Ph4jT/xPiEhAS4uLlAoFMplMpkM8fHxcHJy0mgtAJCdnY09e/agb9++Gt/3yzimhoiIqIzMzMzQokULAOL01ty8eVMl0ACAXC7HrVu3NF4LABgZGakEmqrQB8JQQ0RE9H/EvATl5uYGqVT1tCyTyVCnTh2N1/Kq1NRUBAcHY9OmTWKXUiKGGiIiov9TOFj48OHDRXpNKpqTkxPWrVsHmUwG4EWgWbt2rSiXnl4VEhKCgwcPYtiwYbh06ZLY5bwWx9QQERH9n/z8fFhZWSEzMxMxMTHw8vLSeA0JCQm4desW6tSpUykCDfDiMlhwcDDCwsJQp04dnDt3Dubm5hrZN8fUEBERlYOuri4CAgIAiHdrt5OTEwICAipNoAFe9Bpt3rwZzs7OuHXrFgYOHKjxnqzSYKghIiJ6SWW4tbsysrGxwbZt26Cnp4ddu3Zh4cKFYpdUBEMNERHRSwrH1URFRSEnJ0fkaioXHx8frFq1CgAwffr0Shf8GGqIiIheUr9+fTg6OiInJwcnT54Uu5xKZ+jQoRg8eDAUCgX69++Pe/fuiV2SEkMNERHRSyQSifISlNhTJlRGEokEq1atQtOmTZGSkoLevXsjNzdX7LIAMNQQEREVwXE1JTM0NMS2bdtgZWWFs2fPYuzYsWKXBIChhoiIqIjCUBMdHY3U1FSRq6mcatWqhU2bNkEikWDt2rUIDQ0VuySGGiIiolc5ODigYcOGEAQBERERYpdTaXXq1AmzZs0CACxcuBD5+fmi1sNQQ0REVAxegiqd6dOnY/bs2Th+/Dh0dXVFrYWhhoiIqBiFt3ZzsHDJpFIpZsyYAWtra7FLYaghIiIqjp+fH3R0dHDnzh3cuXNH7HKoFBhqiIiIimFqaor33nsPAC9BVRUMNURERK/BcTVVC0MNERHRaxSOqwkPD4dcLhe5GnoThhoiIqLXaNasGczMzPDkyRPExsaKXQ69AUMNERHRa+jo6KBt27YAeBdUVcBQQ0REVAKOq6k6GGqIiIhKUDiu5vjx43j+/LnI1VBJGGqIiIhKULduXTg5OSE3NxfHjx8XuxwqAUMNERFRCSQSCS9BVREMNURERG/AKROqBoYaIiKiN2jfvj0AICYmBikpKSJXQ6/DUENERPQGdnZ2aNy4MYAXD+KjyomhhoiIqBQ4rqbyY6ghIiIqhZfH1QiCIHI1VByGGiIiolJo06YN9PT0cPfuXdy+fVvscqgYDDVERESlYGxsjFatWgHgXVCVFUMNERFRKXFcTeXGUENERFRKheNqjhw5ArlcLnI19KoqFWoWL16Mhg0bwsPDAxs3bhS7HCIiqma8vb1hbm6OZ8+eITo6Wuxy6BVVJtRcunQJmzdvRnR0NM6dO4c1a9bg2bNnYpdFRETViEwmQ7t27QBwXE1lVGVCzbVr19CqVSsYGBjAwMAAXl5eOHDggNhlERFRNVN4CYrjaioftYWaqKgodO3aFY6OjpBIJNi5c2eRNqtXr4arqysMDAzg7e2NY8eOlXr7Hh4eiIiIwLNnz/Ds2TMcOXIEDx48UFf5REREpVI4WPjkyZPIysoSuRp6mY66NpSVlQVPT0989tln6NWrV5H1W7duxbhx47B69Wr4+vpi7dq1CA4OxtWrV+Hs7AzgxbXK3NzcIp89dOgQGjRogDFjxqBdu3YwNzdHs2bNoKPz+vJzc3NVtpWenq6GoyQiouquTp06cHFxwd27d3Hs2DF06tRJ7JLo/0iECngsokQiwY4dO9C9e3flshYtWqBp06ZYs2aNcln9+vXRvXt3LFiwoMz7GDp0KHr06IHOnTsXu37WrFmYPXt2keVpaWkwMzMr8/6IiIgKDR06FCEhIZg4cSIWL14sdjlaLT09Hebm5qU6f2tkTE1eXh6io6MRGBiosjwwMBAnT54s9XaSk5MBANevX8eZM2cQFBT02rZTpkxBWlqa8nX//v3yFU9ERPSKl6dMoMpDbZefSpKSkgK5XA47OzuV5XZ2dkhKSir1drp3745nz57B2NgY69evL/Hyk76+PvT19ctdMxER0esU3gF18eJFPHr0qMj5jcShkVBTSCKRqLwXBKHIspKUpVeHiIiootSoUQNNmjRBTEwMwsPD8dFHH4ldEkFDl59sbGwgk8mK9MokJycz3RIRUZXEKRMqH42EGj09PXh7exe59hgWFqacHIyIiKgqeXlcTQXcc0PloLbLT5mZmbh165byfVxcHGJjY2FlZQVnZ2dMmDABn376KXx8fPDee+9h3bp1uHfvHkaMGKGuEoiIiDSmdevW0NfXR0JCAm7cuIF69eqJXVK1p7ZQc+7cObRt21b5fsKECQCAgQMHIjQ0FH379kVqairmzJmDxMREeHh4YP/+/XBxcVFXCWohl8uRn58vdhlUDejp6UEqrTIP9SaiVxgaGsLX1xdHjhxBWFgYQ00lUCHPqamM3nSfuyAISEpK4nxSpDFSqRSurq7Q09MTuxQiKqfvvvsOU6ZMQbdu3Yp9kj69vbI8p4ah5v8kJibi2bNnsLW1hZGRUZnuyiIqK4VCgYcPH0JXVxfOzs78fSOqos6dO4dmzZrBzMwMqampJT5qhMqnLKGG3z5eXHIqDDTW1tZil0PVRI0aNfDw4UMUFBRAV1dX7HKIqByaNGkCKysrPHnyBGfPnsV7770ndknVGi/oA8oxNEZGRiJXQtVJ4WUnuVwuciVEVF4ymUz5ID7e2i0+hpqX8BIAaRJ/34i0A6dMqDwYaoiIiN5C4UP4Tp06hczMTJGrqd4YaoiIiN5C7dq1Ubt2bRQUFODo0aNil1OtMdQQERG9JU6ZUDkw1FCpxcfHQyKRIDY2VuxSiIgqFY6rqRwYaoiIiN5S27ZtIZFIcOXKFSQmJopdTrXFUPMagiAgKytL46+yPgsxICAAY8aMwVdffQUrKyvY29tj1qxZyvVpaWkYNmwYbG1tYWZmhnbt2uHChQvKdTKZDNHR0cpjtrKyQrNmzZSf37JlCxwcHAAArq6uAF48l0EikSAgIADAiwfJzZkzB05OTtDX14eXlxcOHDig3EZhD8/27dvRtm1bGBkZwdPTE6dOnSrzz4WIqDKytrZG06ZNAfASlJgYal4jOzsbJiYmGn9lZ2eXudYNGzbA2NgYp0+fxqJFizBnzhzlrLGdO3dGUlIS9u/fj+joaDRt2hTt27fHkydPYG5uDi8vL0RGRgIALl68qPzf9PR0AEBkZCT8/f0BAGfOnAHw4v+wiYmJ2L59OwBg+fLlWLJkCRYvXoyLFy8iKCgIH3zwAW7evKlS57Rp0zBp0iTExsaibt266N+/PwoKCsr18yEiqmwKL0Ex1IiHoUYLNG7cGDNnzoSbmxsGDBgAHx8fhIeHIyIiApcuXcJff/0FHx8fuLm5YfHixbCwsMDff/8N4EVPT2GoiYyMRPv27eHh4YHjx48rlxX2yNSoUQPAi79I7O3tYWVlBQBYvHgxvv76a/Tr1w/16tXDwoUL4eXlhWXLlqnUOWnSJHTu3Bl169bF7NmzcffuXZWZ3YmIqrLCwcKFf1SS5nGahNcwMjIS5XkD5XmqcePGjVXeOzg4IDk5GdHR0cjMzCwy9cPz589x+/ZtAC9CTUhICBQKBY4ePYr27dvD2dkZR48eRdOmTXHjxg1lT01x0tPT8fDhQ/j6+qos9/X1VV7mKq7OwktaycnJcHd3L/MxExFVNr6+vjAwMEBiYiKuXbuGBg0aiF1StcNQ8xoSiQTGxsZil1Eqr84bJJFIoFAooFAo4ODgoOyJeZmFhQUAwM/PDxkZGTh//jyOHTuGuXPnombNmpg/fz68vLxga2uL+vXrv7GGV5+OKwhCkWUv11m4TqFQlOYQiYgqPQMDA7Rp0wZhYWE4fPgwQ40IePlJizVt2hRJSUnQ0dFBnTp1VF42NjYAoBxXs2rVKkgkEjRo0ABt2rRBTEwM9u7dq9JLU9xcRWZmZnB0dFRerip08uTJUoUhIiJtwlu7xcVQo8U6dOiA9957D927d8fBgwcRHx+PkydPYvr06Th37pyyXUBAADZu3Ah/f39IJBJYWlqiQYMG2Lp1q3I8DQDY2trC0NAQBw4cwKNHj5CWlgYA+N///oeFCxdi69atuH79OiZPnozY2FiMHTtW04dMRCSqwnE1kZGRysmSSXMYarSYRCLB/v374efnh8GDB6Nu3bro168f4uPjYWdnp2zXtm1byOVylQDj7+8PuVyu0lOjo6ODFStWYO3atXB0dES3bt0AAGPGjMHEiRMxceJENGrUCAcOHMDu3bvh5uamsWMlIqoMPD09YWNjg8zMTJw+fVrscqodiVBNhminp6fD3NwcaWlpMDMzU1mXk5ODuLg4uLq6wsDAQKQKqbrh7x2RdurXrx+2bt2KmTNnqjw3jMqnpPP3q9hTQ0REpEYcVyMehhoiIiI1KhxXc/r0aeWDTEkzGGqIiIjUyMXFBXXq1IFcLi/2kRpUcRhqiIiI1IxTJoiDoYaIiEjNXp4ygTSHoYaIiEjN2rZtC6lUiv/++w8JCQlil1NtMNQQERGpmaWlJXx8fAAA4eHhIldTfTDUEBERVQDe2q15DDVEREQVoHBczeHDh1FNnnMrOoYaLTVo0CB0795d7DKIiKqt9957D0ZGRnj06BEuX74sdjnVAkMNiS4gIAASiUTl1a9fP7HLIiJ6K/r6+vDz8wPAW7s1haGGKoXPP/8ciYmJytfatWvFLomI6K1xXI1mMdS8QVZW1mtfOTk5pW77/PnzN7Ytj7///huNGjWCoaEhrK2t0aFDB5VtLV68GA4ODrC2tsaoUaOQn5+vXLdx40b4+PjA1NQU9vb2+Oijj5CcnKxcHxkZCYlEgoMHD6JJkyYwNDREu3btkJycjH/++Qf169eHmZkZ+vfvj+zs7HLVX8jIyAj29vbKl7m5+Vttj4ioMigcV3P06FHk5eWJXI32Y6h5AxMTk9e+evXqpdLW1tb2tW2Dg4NV2taqVatIm7JKTExE//79MXjwYFy7dg2RkZHo2bOnckBaREQEbt++jYiICGzYsAGhoaEIDQ1Vfj4vLw9z587FhQsXsHPnTsTFxWHQoEFF9jNr1iysWrUKJ0+exP3799GnTx8sW7YMmzdvxr59+xAWFoaVK1cq28+fP7/E783ExATHjh1T2cemTZtgY2ODhg0bYtKkScjIyCjz90FEVNl4eHjA1tYW2dnZOHXqlNjlaD0dsQug8ktMTERBQQF69uwJFxcXAECjRo2U6y0tLbFq1SrIZDK4u7ujc+fOCA8Px+effw4AGDx4sLJt7dq1sWLFCjRv3hyZmZkqIWvevHnw9fUFAAwZMgRTpkzB7du3Ubt2bQBA7969ERERga+//hoAMGLECPTp06fE2t955x3lvz/++GO4urrC3t4ely9fxpQpU3DhwgV21xJRlSeVStGhQwds3rwZhw8fhr+/v9glaTWGmjfIzMx87TqZTKby/uVLN6+SSlU7xeLj49+qLgDw9PRE+/bt0ahRIwQFBSEwMBC9e/eGpaUlAKBhw4YqNTo4OODSpUvK9zExMZg1axZiY2Px5MkTKBQKAMC9e/fQoEEDZbvGjRsr/21nZwcjIyNloClcdubMGeV7KysrWFlZlfo4CkMW8OKvGjc3N/j4+OD8+fNo2rRpqbdDRFQZFYaasLAwzJ07V+xytBovP72BsbHxa18GBgalbmtoaPjGtmUlk8kQFhaGf/75Bw0aNMDKlStRr149xMXFAQB0dXVV2kskEmVwycrKQmBgIExMTLBx40acPXsWO3bsAIAi131f3o5EIilxu0D5Lj+9rGnTptDV1cXNmzfL/J0QEVU2heNqzp49i2fPnolbjJZjT00VJ5FI4OvrC19fX8yYMQMuLi7KcFKS//77DykpKfjuu+9Qs2ZNAMC5c+fUUlNZLz+96sqVK8jPz4eDg4Na6iEiElPNmjVRr149XL9+HREREejRo4fYJWkthpoq7PTp0wgPD0dgYCBsbW1x+vRpPH78GPXr18fFixdL/KyzszP09PSwcuVKjBgxApcvX1Zbt2hZLj/dvn0bmzZtwvvvvw8bGxtcvXoVEydORJMmTZTjeIiIqrqOHTvi+vXrOHz4MENNBeLlpyrMzMwMUVFReP/991G3bl1Mnz4dS5YsKXKnVXFq1KiB0NBQ/PXXX2jQoAG+++47LF68WANVq9LT00N4eDiCgoJQr149jBkzBoGBgTh8+HCRMUtERFXVy1MmUMWRCNVkQor09HSYm5sjLS0NZmZmKutycnIQFxcHV1fXIuNkiCoKf++Iqo+0tDRYW1tDLpfj7t27cHZ2FrukKqOk8/er2FNDRERUwczNzdG8eXMA7K2pSAw1REREGsApEypepQw1PXr0gKWlJXr37l2mdURERJVV4bia8PBwlcdgkPpUylAzZswY/Pbbb2VeR0REVFm1aNECxsbGePz48RvvUKXyqZShpm3btjA1NS3zOiIiospKT08PAQEBADiupqKUOdRERUWha9eucHR0hEQiwc6dO4u0Wb16tfKODm9v7xKfHktERFRdFF6C4riailHmUJOVlQVPT0+sWrWq2PVbt27FuHHjMG3aNMTExKBNmzYIDg7GvXv3lG28vb3h4eFR5PXw4cPyHwkREVElVzhY+NixY8jJyRG5Gu1T5icKBwcHl/hwt6VLl2LIkCEYOnQoAGDZsmU4ePAg1qxZgwULFgAAoqOjy1lu6eXm5iI3N1f5Pj09vcL3SUREVJIGDRrAwcEBiYmJOHnyJNq1ayd2SVpFrWNq8vLyEB0djcDAQJXlgYGBOHnypDp39UYLFiyAubm58lU4vxFpp1mzZsHLy0vt233dJVYiovKQSCR8unAFUmuoSUlJgVwuh52dncpyOzs7JCUllXo7QUFB+PDDD7F//344OTnh7NmzpVr3silTpiAtLU35un//fvkOipTi4+MhkUgQGxsrdilFTJo0CeHh4WKXQUT0Rgw1FadCJrSUSCQq7wVBKLKsJAcPHizXupfp6+tDX1+/1Pukqs3ExAQmJiZil0FE9EaFoebcuXN48uRJqScApjdTa0+NjY0NZDJZkV6Z5OTkIr032iohIQERERFISEio8H0FBARgzJgx+Oqrr2BlZQV7e3vMmjVLpU1aWhqGDRsGW1tbmJmZoV27drhw4YJynUwmU45xEgQBVlZWaNasmfLzW7ZsgYODAwDA1dUVANCkSRNIJBLlrYkKhQJz5syBk5MT9PX14eXlhQMHDii3UdjDs337drRt2xZGRkbw9PTEqVOnynS8kZGRaN68OYyNjWFhYQFfX1/cvXsXQNHLT4MGDUL37t2xePFiODg4wNraGqNGjUJ+fr6yTWJiIjp37gxDQ0O4urpi8+bNqFWrFpYtW/baGh48eIC+ffvC0tIS1tbW6NatG+Lj48t0HERUvTk6OqJBgwYQBAERERFil6NV1Bpq9PT04O3tXeRWtbCwMLRq1Uqdu6qUQkJC4OLignbt2sHFxQUhISEVvs8NGzbA2NgYp0+fxqJFizBnzhzl9y8IAjp37oykpCTs378f0dHRaNq0Kdq3b48nT57A3NwcXl5eiIyMBADlw6AuXryoHFgdGRkJf39/AMCZM2cAvOgyTUxMxPbt2wEAy5cvx5IlS7B48WJcvHgRQUFB+OCDD3Dz5k2VWqdNm4ZJkyYhNjYWdevWRf/+/VFQUFCq4ywoKED37t3h7++Pixcv4tSpUxg2bFiJPYARERG4ffs2IiIisGHDBoSGhiI0NFS5fsCAAXj48CEiIyOxbds2rFu3DsnJya/dXnZ2Ntq2bQsTExNERUXh+PHjMDExQadOnZCXl1eq4yAiAjhlQoURyigjI0OIiYkRYmJiBADC0qVLhZiYGOHu3buCIAjCH3/8Iejq6gohISHC1atXhXHjxgnGxsZCfHx8WXelVmlpaQIAIS0trci658+fC1evXhWeP39e7u3fv39fkEqlAgDlSyaTCffv33+bskvk7+8vtG7dWmVZs2bNhK+//loQBEEIDw8XzMzMhJycHJU27777rrB27VpBEARhwoQJQpcuXQRBEIRly5YJvXv3Fpo2bSrs27dPEARBqFu3rrBmzRpBEAQhLi5OACDExMSobM/R0VH49ttvi9TxxRdfqHzul19+Ua6/cuWKAEC4du1aqY41NTVVACBERkYWu37mzJmCp6en8v3AgQMFFxcXoaCgQLnsww8/FPr27SsIgiBcu3ZNACCcPXtWuf7mzZsCAOGHH35QLgMg7NixQxAEQQgJCRHq1asnKBQK5frc3FzB0NBQOHjwYKmO42Xq+L0joqppz549AgDh3XffFbuUSq+k8/eryjym5ty5c2jbtq3y/YQJEwAAAwcORGhoKPr27YvU1FTMmTMHiYmJ8PDwwP79++Hi4vI22avSu3nzZpG5PORyOW7dugUnJ6cK22/jxo1V3js4OCh7G6Kjo5GZmQlra2uVNs+fP8ft27cBvLiEFRISAoVCgaNHj6J9+/ZwdnbG0aNH0bRpU9y4cUPZU1Oc9PR0PHz4EL6+virLfX19lZe5iqu18JJWcnIy3N3d33icVlZWGDRoEIKCgtCxY0d06NABffr0UW6nOA0bNoRMJlPZ56VLlwAA169fh46ODpo2bapcX6dOHVhaWr52e9HR0bh161aRJ1rn5OQov08iotLw9/eHjo4Obt++jbi4OOXlfXo7ZQ41AQEBEAShxDZffPEFvvjii3IXVRW5ublBKpWqBBuZTIY6depU6H51dXVV3kskEmUNCoUCDg4OystLL7OwsAAA+Pn5ISMjA+fPn8exY8cwd+5c1KxZE/Pnz4eXlxdsbW1Rv379N9ZRmsHhL9dauK4sk7qtX78eY8aMwYEDB7B161ZMnz4dYWFhaNmyZbHtS/puXvc7XNLvtkKhgLe3NzZt2lRkXY0aNUp7GEREMDU1RcuWLXH8+HEcPnwYn3/+udglaYVKOfdTVeTk5IR169YpewZkMhnWrl1bob00b9K0aVMkJSVBR0cHderUUXnZ2NgAgHJczapVqyCRSNCgQQO0adMGMTEx2Lt3r0ovjZ6eHoAXPVCFzMzM4OjoiOPHj6vs++TJk6UKQ2XVpEkTTJkyBSdPnoSHhwc2b95cru24u7ujoKAAMTExymW3bt3Cs2fPXvuZpk2b4ubNm7C1tS3yfZqbm5erDiKqvjhlgvox1KjRkCFDEB8fj4iICMTHx2PIkCGi1tOhQwe899576N69Ow4ePIj4+HicPHkS06dPx7lz55TtAgICsHHjRvj7+0MikcDS0hINGjTA1q1blXc4AYCtrS0MDQ1x4MABPHr0CGlpaQCA//3vf1i4cCG2bt2K69evY/LkyYiNjcXYsWPVdixxcXGYMmUKTp06hbt37+LQoUO4ceNGuYOTu7s7OnTogGHDhuHMmTOIiYnBsGHDYGho+NrBxx9//DFsbGzQrVs3HDt2DHFxcTh69CjGjh2rkbvdiEi7FA4WDg8PL1OvNb0eQ42aOTk5ISAgQNQemkISiQT79++Hn58fBg8ejLp166Jfv36Ij49XucW+bdu2kMvlKgHG398fcrlcpadGR0cHK1aswNq1a+Ho6Ihu3boBAMaMGYOJEydi4sSJaNSoEQ4cOIDdu3fDzc2tzPW+fHfSy4yMjPDff/+hV69eqFu3LoYNG4Yvv/wSw4cPL9M+Xvbbb7/Bzs4Ofn5+6NGjBz7//HOYmprCwMDgtTVERUXB2dkZPXv2RP369TF48GA8f/4cZmZm5a6DiKqnZs2awdTUFE+ePFHpNabykwhvGiCjJdLT02Fubo60tLQiJ6CcnBzlQK3XndCoYsXHx8PNzQ1Xr14tcxhSl4SEBNSsWROHDx9G+/btK3x//L0jom7dumH37t347rvv8PXXX4tdTqVU0vn7VeypoUrhwIEDGDZsmEYDzZEjR7B7927ExcXh5MmT6NevH2rVqgU/Pz+N1UBE1RvH1ahXhUyTQFRWI0aM0Pg+8/PzMXXqVNy5cwempqZo1aoVNm3aVOSuKSKiilI4rub48eN4/vw5DA0NRa6oamOooWorKCgIQUFBYpdBRNVYvXr18M477+DBgwc4ceKEsueGyoeXn4iIiEQikUg4ZYIaMdQQERGJqLB35vDhwyJXUvUx1BAREYmoMNTExMQgJSVF5GqqNoYaIiIiEdnZ2aFRo0YQBAFHjhwRu5wqjaGGiIhIZLy1Wz0YaoiIiET28mDhavJM3ArBUKOlBg0ahO7du4tdRqUQGRkJiURS4mSV5REQEIBx48apdZtEVD35+flBV1cXd+/exe3bt8Uup8piqCHRBQQEQCKRqLz69euntu23atUKiYmJnEmbiCotY2NjtGrVCgDvgnobDDVUKXz++edITExUvtauXau2bevp6cHe3v61s28TEVUGHFfz9hhqqrC///4bjRo1gqGhIaytrdGhQwdkZWWptFm8eDEcHBxgbW2NUaNGIT8/X7lu48aN8PHxgampKezt7fHRRx8hOTlZub7wss3BgwfRpEkTGBoaol27dkhOTsY///yD+vXrw8zMDP3790d2dvZbHYuRkRHs7e2Vr7L2qty9exddu3aFpaUljI2N0bBhQ+zfv1/lOAovP4WGhsLCwgIHDx5E/fr1YWJigk6dOiExMVG5vYKCAowZMwYWFhawtrbG119/jYEDB5Z4SS8vLw9fffUV3nnnHRgbG6NFixaIjIws61dBRNVU4biaI0eOQC6Xi1xN1cRQ8wZZWVllfhUUFCg/X1BQgKysLDx//vyN2y2LxMRE9O/fH4MHD8a1a9cQGRmJnj17qgwwi4iIwO3btxEREYENGzYgNDQUoaGhyvV5eXmYO3cuLly4gJ07dyIuLg6DBg0qsq9Zs2Zh1apVOHnyJO7fv48+ffpg2bJl2Lx5M/bt24ewsDCsXLlS2X7+/PkwMTEp8XXs2DGVfWzatAk2NjZo2LAhJk2ahIyMjDJ9H6NGjUJubi6ioqJw6dIlLFy4ECYmJq9tn52djcWLF+P3339HVFQU7t27h0mTJinXL1y4EJs2bcL69etx4sQJpKenY+fOnSXW8Nlnn+HEiRP4448/cPHiRXz44Yfo1KkTbt68WaZjIaLqydvbG+bm5nj27BnOnz8vdjlVk1BNpKWlCQCEtLS0IuueP38uXL16VXj+/HmRdQDK/Przzz+Vn//zzz8FAIK/v7/Kdm1sbIp8riyio6MFAEJ8fHyx6wcOHCi4uLgIBQUFymUffvih0Ldv39du88yZMwIAISMjQxAEQYiIiBAACIcPH1a2WbBggQBAuH37tnLZ8OHDhaCgIOX71NRU4ebNmyW+srOzle3XrVsnhIWFCZcuXRK2bNki1KpVS+jQoUOZvo9GjRoJs2bNKnZd4XE8ffpUEARBWL9+vQBAuHXrlrLNjz/+KNjZ2Snf29nZCd9//73yfUFBgeDs7Cx069ZNuczf318YO3asIAiCcOvWLUEikQgPHjxQ2Xf79u2FKVOmFFtXSb93RFQ99ejRQwAgfPvtt2KXUmmUdP5+FSe0rKI8PT3Rvn17NGrUCEFBQQgMDETv3r1haWmpbNOwYUPIZDLlewcHB1y6dEn5PiYmBrNmzUJsbCyePHkChUIBALh37x4aNGigbNe4cWPlv+3s7GBkZITatWurLDtz5ozyvZWVFaysrEp9LJ9//rny3x4eHnBzc4OPjw/Onz+Ppk2blmobY8aMwciRI3Ho0CF06NABvXr1Uqn7VUZGRnj33XeV7x0cHJSX3tLS0vDo0SM0b95cuV4mk8Hb21v5Hb3q/PnzEAQBdevWVVmem5sLa2vrUh0DEVGHDh2wY8cOHD58GFOnThW7nCqHoeYNMjMzy/wZfX195b979OiBzMxMSKWqV/ri4+Pfqi6ZTIawsDCcPHkShw4dwsqVKzFt2jScPn0arq6uAABdXV2Vz0gkEuVJOSsrC4GBgQgMDMTGjRtRo0YN3Lt3D0FBQcjLy1P53MvbkUgkJW4XeHH5af78+SXW/88//6BNmzbFrmvatCl0dXVx8+bNUoeaoUOHIigoCPv27cOhQ4ewYMECLFmyBKNHjy62fXHHILzybIhXBxa/uv5lCoUCMpkM0dHRKkESQImXwYiIXlY4rubEiRPIzs6GkZGRyBVVLQw1b2BsbPxWn9fR0YGOTtGv+W23C7w46fr6+sLX1xczZsyAi4sLduzYgQkTJrzxs//99x9SUlLw3XffoWbNmgCAc+fOvXVNADBixAj06dOnxDbvvPPOa9dduXIF+fn5cHBwKNN+a9asiREjRmDEiBGYMmUKfv7559eGmpKYm5sre58Kg5dcLkdMTAy8vLyK/UyTJk0gl8uRnJz82rBGRPQmderUgbOzM+7du4djx44hKChI7JKqFIaaKur06dMIDw9HYGAgbG1tcfr0aTx+/Bj169cv1eednZ2hp6eHlStXYsSIEbh8+TLmzp2rltrKcvnp9u3b2LRpE95//33Y2Njg6tWrmDhxIpo0aQJfX99S73PcuHEIDg5G3bp18fTpUxw5cqTU30VxRo8ejQULFqBOnTpwd3fHypUr8fTp09feFl63bl18/PHHGDBgAJYsWYImTZogJSUFR44cQaNGjfD++++XuxYiqj4kEgk6dOiAX3/9FWFhYQw1ZcRQU0WZmZkhKioKy5YtQ3p6OlxcXLBkyRIEBweX6vM1atRAaGgopk6dihUrVqBp06ZYvHgxPvjggwquXJWenh7Cw8OxfPlyZGZmombNmujcuTNmzpypchknICAAtWrVUrl762VyuRyjRo1CQkICzMzM0KlTJ/zwww/lruvrr79GUlISBgwYAJlMhmHDhiEoKKjIpaWXrV+/HvPmzcPEiRPx4MEDWFtb47333mOgIaIy6du3LywtLdGzZ0+xS6lyJEJJAwW0SHp6OszNzZGWlgYzMzOVdTk5OYiLi4OrqysMDAxEqpBKUqtWLcyaNavYW841QaFQoH79+ujTp4/aerT4e0dE9GYlnb9fxZ4aqvT+++8/mJqaYsCAARrb5927d3Ho0CH4+/sjNzcXq1atQlxcHD766CON1UBERGXDh+9Rpefu7o5Lly4VuYOsIkmlUoSGhqJZs2bw9fXFpUuXcPjw4bcap0NERBWLPTVExahZsyZOnDghdhlERFQG7KkhIiIircBQQ0RERFqBoeYlr3sEPlFFqCY3HhIRaQzH1ODFs1KkUikePnyIGjVqQE9P77UPWSNSB0EQ8Pjx42KnnSAiovJhqMGLO11cXV2RmJiIhw8fil0OVRMSiQROTk4lPtCPiIhKj6Hm/+jp6cHZ2RkFBQWQy+Vil0PVgK6uLgMNEZEaMdS8pPBSAC8HEBERVT0cKExERERagaGGiIiItAJDDREREWmFajOmpvCZIOnp6SJXQkRERKVVeN4uzbO9qk2oycjIAPBiTh8iIiKqWjIyMmBubl5iG4lQTR5rqlAo8PDhQ5iamqr9wXrp6emoWbMm7t+/DzMzM7Vum/4/fs+awe9ZM/g9awa/Z82pqO9aEARkZGTA0dERUmnJo2aqTU+NVCqFk5NThe7DzMyM/6fRAH7PmsHvWTP4PWsGv2fNqYjv+k09NIU4UJiIiIi0AkMNERERaQWGGjXQ19fHzJkzoa+vL3YpWo3fs2bwe9YMfs+awe9ZcyrDd11tBgoTERGRdmNPDREREWkFhhoiIiLSCgw1REREpBUYaoiIiEgrMNQQERGRVmCoeUurV6+Gq6srDAwM4O3tjWPHjoldklZZsGABmjVrBlNTU9ja2qJ79+64fv262GVpvQULFkAikWDcuHFil6KVHjx4gE8++QTW1tYwMjKCl5cXoqOjxS5LqxQUFGD69OlwdXWFoaEhateujTlz5kChUIhdWpUWFRWFrl27wtHRERKJBDt37lRZLwgCZs2aBUdHRxgaGiIgIABXrlzRWH0MNW9h69atGDduHKZNm4aYmBi0adMGwcHBuHfvntilaY2jR49i1KhR+PfffxEWFoaCggIEBgYiKytL7NK01tmzZ7Fu3To0btxY7FK00tOnT+Hr6wtdXV38888/uHr1KpYsWQILCwuxS9MqCxcuxE8//YRVq1bh2rVrWLRoEb7//nusXLlS7NKqtKysLHh6emLVqlXFrl+0aBGWLl2KVatW4ezZs7C3t0fHjh2Vk0pXOIHKrXnz5sKIESNUlrm7uwuTJ08WqSLtl5ycLAAQjh49KnYpWikjI0Nwc3MTwsLCBH9/f2Hs2LFil6R1vv76a6F169Zil6H1OnfuLAwePFhlWc+ePYVPPvlEpIq0DwBhx44dyvcKhUKwt7cXvvvuO+WynJwcwdzcXPjpp580UhN7asopLy8P0dHRCAwMVFkeGBiIkydPilSV9ktLSwMAWFlZiVyJdho1ahQ6d+6MDh06iF2K1tq9ezd8fHzw4YcfwtbWFk2aNMHPP/8sdllap3Xr1ggPD8eNGzcAABcuXMDx48fx/vvvi1yZ9oqLi0NSUpLKeVFfXx/+/v4aOy9Wm1m61S0lJQVyuRx2dnYqy+3s7JCUlCRSVdpNEARMmDABrVu3hoeHh9jlaJ0//vgD58+fx9mzZ8UuRavduXMHa9aswYQJEzB16lScOXMGY8aMgb6+PgYMGCB2eVrj66+/RlpaGtzd3SGTySCXy/Htt9+if//+YpemtQrPfcWdF+/evauRGhhq3pJEIlF5LwhCkWWkHl9++SUuXryI48ePi12K1rl//z7Gjh2LQ4cOwcDAQOxytJpCoYCPjw/mz58PAGjSpAmuXLmCNWvWMNSo0datW7Fx40Zs3rwZDRs2RGxsLMaNGwdHR0cMHDhQ7PK0mpjnRYaacrKxsYFMJivSK5OcnFwkpdLbGz16NHbv3o2oqCg4OTmJXY7WiY6ORnJyMry9vZXL5HI5oqKisGrVKuTm5kImk4lYofZwcHBAgwYNVJbVr18f27ZtE6ki7fS///0PkydPRr9+/QAAjRo1wt27d7FgwQKGmgpib28P4EWPjYODg3K5Js+LHFNTTnp6evD29kZYWJjK8rCwMLRq1UqkqrSPIAj48ssvsX37dhw5cgSurq5il6SV2rdvj0uXLiE2Nlb58vHxwccff4zY2FgGGjXy9fUt8liCGzduwMXFRaSKtFN2djakUtVTnEwm4y3dFcjV1RX29vYq58W8vDwcPXpUY+dF9tS8hQkTJuDTTz+Fj48P3nvvPaxbtw737t3DiBEjxC5Na4waNQqbN2/Grl27YGpqquwZMzc3h6GhocjVaQ9TU9Mi45SMjY1hbW3N8UtqNn78eLRq1Qrz589Hnz59cObMGaxbtw7r1q0TuzSt0rVrV3z77bdwdnZGw4YNERMTg6VLl2Lw4MFil1alZWZm4tatW8r3cXFxiI2NhZWVFZydnTFu3DjMnz8fbm5ucHNzw/z582FkZISPPvpIMwVq5B4rLfbjjz8KLi4ugp6entC0aVPeaqxmAIp9rV+/XuzStB5v6a44e/bsETw8PAR9fX3B3d1dWLdundglaZ309HRh7NixgrOzs2BgYCDUrl1bmDZtmpCbmyt2aVVaREREsf9NHjhwoCAIL27rnjlzpmBvby/o6+sLfn5+wqVLlzRWn0QQBEEz8YmIiIio4nBMDREREWkFhhoiIiLSCgw1REREpBUYaoiIiEgrMNQQERGRVmCoISIiIq3AUENERERagaGGiIiItAJDDREREWkFhhoiIiLSCgw1REREpBX+HxFS2LQNJDXeAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bvpoutn64=bvptest(U0, FV, FPV, bdata;  bsham=1);\n",
    "bvpouts64=bvptest(U0, FV, FPV, bdata; bsham=5);\n",
    "bvpoutn32=bvptest(U0, FV, FPV32, bdata;  bsham=1);\n",
    "bvpouts32=bvptest(U0, FV, FPV32, bdata; bsham=5);\n",
    "newtonhist=bvpoutn64.history./bvpoutn64.history[1];\n",
    "shamhist=bvpouts64.history./bvpoutn64.history[1];\n",
    "newtonhist32=bvpoutn32.history./bvpoutn64.history[1];\n",
    "shamhist32=bvpouts32.history./bvpoutn64.history[1];\n",
    "nn=length(newtonhist);\n",
    "ns=length(shamhist);\n",
    "nn32=length(newtonhist32);\n",
    "ns32=length(shamhist32);\n",
    "semilogy(0:nn-1,newtonhist,\"k-\",0:ns-1,shamhist,\"k--\",\n",
    "    0:nn32-1,newtonhist32,\"k.\",0:ns32-1,shamhist32,\"k-.\")\n",
    "legend([\"newton\",\"sham=5\",\"newton, single\",\n",
    "        \"sham=5, single\"]);\n",
    "title(\"Figure 2.8 from print book\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now query the iteration statistics to see what the line search did for\n",
    "Newton's method in double precision. We see that the line search\n",
    "was active in the middle of the iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×10 adjoint(::Vector{Int64}) with eltype Int64:\n",
       " 0  0  0  2  2  1  0  0  0  0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bvpoutn64.stats.iarm'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Updating the Jacobian less frequently has a smaller impact than in the case of the H-equation. You should expect this because the linear solver costs $O(N)$ work rather than $O(N^3)$ and the Jacobian is updated for the early iterations because the residual is not decreasing fast enough.\n",
    "\n",
    "The other interesting feature of this computation is that the single precision results differ in a meaningful way from the double precision results. The reason for this is that the Jacobian is ill-conditioned enough to affect the quality of the single precision solver. The condition number of the Jacobian is $O(N)$, so for $N=10^5$ one should expect different results from the linear solve in single precision, and that's what you get.\n",
    "\n",
    "Moreover, the benchmark results say that there is no benefit from doing the linear algebra in single precision. This is also no surprise since the factorization is $O(N)$ work.\n",
    "\n",
    "The reader should try this problem with $N=10^6$ and watch the line search fail for the single precision linear algebra computations. That is also no surprise with an inaccurate Newton direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytic, double, Newton\n",
      "  100.058 ms (147 allocations: 47.31 MiB)\n",
      "analytic, double, sham=5\n",
      "  87.431 ms (139 allocations: 44.26 MiB)\n",
      "analytic, single, Newton\n",
      "  92.764 ms (201 allocations: 54.18 MiB)\n",
      "analytic, single, sham=5\n",
      "  85.293 ms (206 allocations: 57.23 MiB)\n"
     ]
    }
   ],
   "source": [
    "println(\"analytic, double, Newton\"); \n",
    "@btime bvptest($U0, $FV, $FPV, $bdata;  bsham=1);\n",
    "println(\"analytic, double, sham=5\"); \n",
    "@btime bvptest($U0, $FV, $FPV, $bdata;  bsham=5);\n",
    "println(\"analytic, single, Newton\"); \n",
    "@btime bvptest($U0, $FV, $FPV32, $bdata;  bsham=1);\n",
    "println(\"analytic, single, sham=5\"); \n",
    "@btime bvptest($U0, $FV, $FPV32, $bdata;  bsham=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 2.8.4: Shamanskii for the Convection-Diffusion Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dominant cost in the solution of the PDE problem is the computation and factorization of the sparse Jacobian. Even for the small problem with $31^2 = 961$ unknowns, one can see the effects. We will compare a Newton iteration with the default strategy using ```@btime```. The code that generated the results in section 2.7.5 contained the initialization. It's important to separate that from the solver if you're doing benchmarking. So we will set up the problem first and then benchmark the solver phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Newton\n",
      "  5.326 ms (823 allocations: 9.33 MiB)\n",
      "sham=5\n",
      "  3.671 ms (586 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, \n",
    "# same sparsity pattern as the discrete Laplacian\n",
    "J=copy(pdata.D2);\n",
    "println(\"Newton\"); \n",
    "@btime nsol(pdeF!, u0, FV, J, pdeJ!; \n",
    "    resdec=.5, rtol=1.e-7, atol=1.e-10, \n",
    "    pdata=pdata, sham=1);\n",
    "println(\"sham=5\"); \n",
    "@btime nsol(pdeF!, u0, FV, J, pdeJ!; \n",
    "    resdec=.5, rtol=1.e-7, atol=1.e-10, \n",
    "    pdata=pdata, sham=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Shamanskii iteration, which forms and factors the Jacobian only twice, saves about 40% in time. This is a larger savings than we saw with the two-point boundary value problem, but far from dramatic. The __SuiteSparse__ solvers are very good. The reader should try the same test with a finer mesh to see if the results change. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 2.8.5: ptcsol.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__ptcsol.jl__ is our $\\ptc$ solve. As usual, we begin with the docstrings ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "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[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[22mi \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 htol aasol nsol cool\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "ptcsol(F!, x0, FS, FPS, J! = diffjac!; rtol=1.e-6, atol=1.e-12,              maxit=20, delta0=1.e-6, dx=1.e-7, pdata = nothing, jfact = klfact,                printerr = true, keepsolhist = false, jknowsdt = false)\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  Herewith: some new stuff ==> ptcsol\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 Jacobian 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 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 Jacobian. It is an N x N matrix\n",
       "\n",
       "If FPS is sparse, you \\textbf{must} allocate storage for the diagonal so I will have room to put 1/dt in there.\n",
       "\n",
       "\n",
       "\\item J!: Jacobian evaluation, the ! indicates that J! overwrites FPS, your   preallocated storage for the Jacobian. If you leave this out the   default is a finite difference Jacobian.\n",
       "\n",
       "So, FP=J!(FP,FS,x) or FP=J!(FP,FS,x,pdata) returns FP=F'(x);   (FP,FS, x) must be the argument list, even if FP does not need FS.   One reason for this is that the finite-difference Jacobian   does and that is the default in the solver.\n",
       "\n",
       "Your Jacobian function MUST have –> return FP <– at the end.    See the examples in the TestProblems/Systems/FBeam!.jl\n",
       "\n",
       "You may have a better way to add (1/dt) I to your Jacobian. If you   want to do this yourself then your Jacobian function should be   FP=J!(FP,FS,x,dt) or FP=J!(FP,FS,x,dt,pdata) and return   F'(x) + (1.0/dt)*I. \n",
       "\n",
       "You will also have to set the kwarg \\textbf{jknowsdt} to true.\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, Float32, or Float16 and ptcsol will do the right thing if   YOU do not destroy the declaration in your J! function. I'm amazed   that this works so easily. If the Jacobian is reasonably well   conditioned, you can cut the cost of Jacobian factorization and   storage in half with no loss. For large dense Jacobians and inexpensive   functions, this is a good deal.\n",
       "\n",
       "BUT ... There is very limited support for direct sparse solvers in   anything other than Float64. I recommend that you only use Float64   with direct sparse solvers unless you really know what you're doing. I   have a couple examples in the notebook, but watch out.\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",
       "dx: default = 1.e-7\n",
       "\n",
       "difference increment in finite-difference derivatives       h=dx*norm(x)+1.e-6\n",
       "\n",
       "pdata:\n",
       "\n",
       "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",
       "jfact: default = klfact (tries to figure out best choice) \n",
       "\n",
       "If your Jacobian has any special structure, please set jfact to the correct choice for a factorization.\n",
       "\n",
       "I use jfact when I call PTCUpdate to evaluate the Jacobian (using your J!) and factor it. The default is to use klfact (an internal function) to do something reasonable. For general dense matrices, klfact picks lu! to compute an LU factorization and share storage with the Jacobian.  You may change LU to something else by, for example, setting jfact = cholseky! if your Jacobian is spd.\n",
       "\n",
       "klfact knows about banded matrices and picks qr. You should, however RTFM, allocate the extra two upper bands, and use jfact=qr! to override klfact.\n",
       "\n",
       "klfact uses lu for general sparse matrices.\n",
       "\n",
       "If you give me something that klfact does not know how to dispatch on, then nothing happens. I just return the original Jacobian matrix and  ptcsol will use backslash to compute the Newton step.\n",
       "\n",
       "I know that this is probably not optimal in your situation, so it is  good to pick something else, like jfact = lu.\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",
       "jknowsdt: default = false\n",
       "\n",
       "Set this to true if your Jacobian evaluation function returns F'(x) + (1/dt) I. You'll also need to follow the rules above for the Jacobian evaluation function. I do not recommend this and if your Jacobian is anything other than a matrix I can't promise anything. I've tested this for matrix outputs only.\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\n",
       "\n",
       "Unlike nsol, nsoli, or even ptcsoli, ptcsol has a fixed cost per  iteration of one function, one Jacobian, and one Factorization. Hence iteration statistics are not interesting and not in the output. \n",
       "\n",
       "idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "errcode = 0 if the iteration succeeded         = -1 if the initial iterate satisfies the termination criteria         = 10 if no convergence after maxit iterations\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",
       "\\subsubsection{Example for ptcsol}\n",
       "\\paragraph{The buckling beam problem.}\n",
       "You'll need to use TestProblems for this to work.\n",
       "\n",
       "\\begin{verbatim}\n",
       "julia> using SIAMFANLEquations.TestProblems\n",
       "\n",
       "julia> n=63; maxit=1000; delta = 0.01; lambda = 20.0;\n",
       "\n",
       "julia> bdata = beaminit(n, 0.0, lambda); x = bdata.x;\n",
       "\n",
       "julia> u0 = x .* (1.0 .- x) .* (2.0 .- x);\n",
       "\n",
       "julia> u0 .*= exp.(-10.0 * u0);\n",
       "\n",
       "julia> FS = copy(u0); FPS = copy(bdata.D2);\n",
       "\n",
       "julia> pout = ptcsol( FBeam!, u0, FS, FPS, BeamJ!; \n",
       " rtol = 1.e-10, pdata = bdata, delta0 = delta, maxit = maxit);\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 Array{Float64,2}:\n",
       " 6.31230e+01  9.75412e-01\n",
       " 7.52624e+00  8.35295e-02\n",
       " 8.31545e+00  6.58797e-04\n",
       " 3.15455e+01  4.12697e-08\n",
       " 3.66566e+01  6.29295e-12\n",
       "\n",
       "julia> # We get the nonnegative steady state.\n",
       "       maximum(pout.solution)\n",
       "2.19086e+00\n",
       "\\end{verbatim}\n"
      ],
      "text/markdown": [
       "ptcsol(F!, x0, FS, FPS, J! = diffjac!; rtol=1.e-6, atol=1.e-12,              maxit=20, delta0=1.e-6, dx=1.e-7, pdata = nothing, jfact = klfact,                printerr = true, keepsolhist = false, jknowsdt = false)\n",
       "\n",
       "C. T. Kelley, 2022\n",
       "\n",
       "Julia versions of the nonlinear solvers from my SIAM books.  Herewith: some new stuff ==> ptcsol\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 Jacobian 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 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 Jacobian. It is an N x N matrix\n",
       "\n",
       "    If FPS is sparse, you **must** allocate storage for the diagonal so I will have room to put 1/dt in there.\n",
       "  * J!: Jacobian evaluation, the ! indicates that J! overwrites FPS, your   preallocated storage for the Jacobian. If you leave this out the   default is a finite difference Jacobian.\n",
       "\n",
       "    So, FP=J!(FP,FS,x) or FP=J!(FP,FS,x,pdata) returns FP=F'(x);   (FP,FS, x) must be the argument list, even if FP does not need FS.   One reason for this is that the finite-difference Jacobian   does and that is the default in the solver.\n",
       "\n",
       "    Your Jacobian function MUST have –> return FP <– at the end.    See the examples in the TestProblems/Systems/FBeam!.jl\n",
       "\n",
       "    You may have a better way to add (1/dt) I to your Jacobian. If you   want to do this yourself then your Jacobian function should be   FP=J!(FP,FS,x,dt) or FP=J!(FP,FS,x,dt,pdata) and return   F'(x) + (1.0/dt)*I. \n",
       "\n",
       "    You will also have to set the kwarg **jknowsdt** to true.\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, Float32, or Float16 and ptcsol will do the right thing if   YOU do not destroy the declaration in your J! function. I'm amazed   that this works so easily. If the Jacobian is reasonably well   conditioned, you can cut the cost of Jacobian factorization and   storage in half with no loss. For large dense Jacobians and inexpensive   functions, this is a good deal.\n",
       "\n",
       "    BUT ... There is very limited support for direct sparse solvers in   anything other than Float64. I recommend that you only use Float64   with direct sparse solvers unless you really know what you're doing. I   have a couple examples in the notebook, but watch out.\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",
       "dx: default = 1.e-7\n",
       "\n",
       "difference increment in finite-difference derivatives       h=dx*norm(x)+1.e-6\n",
       "\n",
       "pdata:\n",
       "\n",
       "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",
       "jfact: default = klfact (tries to figure out best choice) \n",
       "\n",
       "If your Jacobian has any special structure, please set jfact to the correct choice for a factorization.\n",
       "\n",
       "I use jfact when I call PTCUpdate to evaluate the Jacobian (using your J!) and factor it. The default is to use klfact (an internal function) to do something reasonable. For general dense matrices, klfact picks lu! to compute an LU factorization and share storage with the Jacobian.  You may change LU to something else by, for example, setting jfact = cholseky! if your Jacobian is spd.\n",
       "\n",
       "klfact knows about banded matrices and picks qr. You should, however RTFM, allocate the extra two upper bands, and use jfact=qr! to override klfact.\n",
       "\n",
       "klfact uses lu for general sparse matrices.\n",
       "\n",
       "If you give me something that klfact does not know how to dispatch on, then nothing happens. I just return the original Jacobian matrix and  ptcsol will use backslash to compute the Newton step.\n",
       "\n",
       "I know that this is probably not optimal in your situation, so it is  good to pick something else, like jfact = lu.\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",
       "jknowsdt: default = false\n",
       "\n",
       "Set this to true if your Jacobian evaluation function returns F'(x) + (1/dt) I. You'll also need to follow the rules above for the Jacobian evaluation function. I do not recommend this and if your Jacobian is anything other than a matrix I can't promise anything. I've tested this for matrix outputs only.\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\n",
       "\n",
       "Unlike nsol, nsoli, or even ptcsoli, ptcsol has a fixed cost per  iteration of one function, one Jacobian, and one Factorization. Hence iteration statistics are not interesting and not in the output. \n",
       "\n",
       "idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "errcode = 0 if the iteration succeeded         = -1 if the initial iterate satisfies the termination criteria         = 10 if no convergence after maxit iterations\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.\n",
       "\n",
       "```jldoctest\n",
       "julia> using SIAMFANLEquations.TestProblems\n",
       "\n",
       "julia> n=63; maxit=1000; delta = 0.01; lambda = 20.0;\n",
       "\n",
       "julia> bdata = beaminit(n, 0.0, lambda); x = bdata.x;\n",
       "\n",
       "julia> u0 = x .* (1.0 .- x) .* (2.0 .- x);\n",
       "\n",
       "julia> u0 .*= exp.(-10.0 * u0);\n",
       "\n",
       "julia> FS = copy(u0); FPS = copy(bdata.D2);\n",
       "\n",
       "julia> pout = ptcsol( FBeam!, u0, FS, FPS, BeamJ!; \n",
       " rtol = 1.e-10, pdata = bdata, delta0 = delta, maxit = maxit);\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 Array{Float64,2}:\n",
       " 6.31230e+01  9.75412e-01\n",
       " 7.52624e+00  8.35295e-02\n",
       " 8.31545e+00  6.58797e-04\n",
       " 3.15455e+01  4.12697e-08\n",
       " 3.66566e+01  6.29295e-12\n",
       "\n",
       "julia> # We get the nonnegative steady state.\n",
       "       maximum(pout.solution)\n",
       "2.19086e+00\n",
       "```\n"
      ],
      "text/plain": [
       "  ptcsol(F!, x0, FS, FPS, J! = diffjac!; rtol=1.e-6, atol=1.e-12, maxit=20,\n",
       "  delta0=1.e-6, dx=1.e-7, pdata = nothing, jfact = klfact, printerr = true,\n",
       "  keepsolhist = false, jknowsdt = false)\n",
       "\n",
       "  C. T. Kelley, 2022\n",
       "\n",
       "  Julia versions of the nonlinear solvers from my SIAM books. Herewith: some\n",
       "  new stuff ==> ptcsol\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 Jacobian in advance –> in the\n",
       "  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 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 Jacobian. It is an N x N matrix\n",
       "       If FPS is sparse, you \u001b[1mmust\u001b[22m allocate storage for the diagonal so I\n",
       "       will have room to put 1/dt in there.\n",
       "\n",
       "    •  J!: Jacobian evaluation, the ! indicates that J! overwrites FPS,\n",
       "       your preallocated storage for the Jacobian. If you leave this out\n",
       "       the default is a finite difference Jacobian.\n",
       "       So, FP=J!(FP,FS,x) or FP=J!(FP,FS,x,pdata) returns FP=F'(x);\n",
       "       (FP,FS, x) must be the argument list, even if FP does not need FS.\n",
       "       One reason for this is that the finite-difference Jacobian does\n",
       "       and that is the default in the solver.\n",
       "       Your Jacobian function MUST have –> return FP <– at the end. See\n",
       "       the examples in the TestProblems/Systems/FBeam!.jl\n",
       "       You may have a better way to add (1/dt) I to your Jacobian. If you\n",
       "       want to do this yourself then your Jacobian function should be\n",
       "       FP=J!(FP,FS,x,dt) or FP=J!(FP,FS,x,dt,pdata) and return F'(x) +\n",
       "       (1.0/dt)*I.\n",
       "       You will also have to set the kwarg \u001b[1mjknowsdt\u001b[22m to true.\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, Float32, or Float16 and\n",
       "       ptcsol will do the right thing if YOU do not destroy the\n",
       "       declaration in your J! function. I'm amazed that this works so\n",
       "       easily. If the Jacobian is reasonably well conditioned, you can\n",
       "       cut the cost of Jacobian factorization and storage in half with no\n",
       "       loss. For large dense Jacobians and inexpensive functions, this is\n",
       "       a good deal.\n",
       "       BUT ... There is very limited support for direct sparse solvers in\n",
       "       anything other than Float64. I recommend that you only use Float64\n",
       "       with direct sparse solvers unless you really know what you're\n",
       "       doing. I have a couple examples in the notebook, but watch out.\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",
       "  dx: default = 1.e-7\n",
       "\n",
       "  difference increment in finite-difference derivatives h=dx*norm(x)+1.e-6\n",
       "\n",
       "  pdata:\n",
       "\n",
       "  precomputed data for the function/Jacobian. Things will go better if you use\n",
       "  this rather than hide the data in global variables within the module for\n",
       "  your function/Jacobian\n",
       "\n",
       "  jfact: default = klfact (tries to figure out best choice)\n",
       "\n",
       "  If your Jacobian has any special structure, please set jfact to the correct\n",
       "  choice for a factorization.\n",
       "\n",
       "  I use jfact when I call PTCUpdate to evaluate the Jacobian (using your J!)\n",
       "  and factor it. The default is to use klfact (an internal function) to do\n",
       "  something reasonable. For general dense matrices, klfact picks lu! to\n",
       "  compute an LU factorization and share storage with the Jacobian. You may\n",
       "  change LU to something else by, for example, setting jfact = cholseky! if\n",
       "  your Jacobian is spd.\n",
       "\n",
       "  klfact knows about banded matrices and picks qr. You should, however RTFM,\n",
       "  allocate the extra two upper bands, and use jfact=qr! to override klfact.\n",
       "\n",
       "  klfact uses lu for general sparse matrices.\n",
       "\n",
       "  If you give me something that klfact does not know how to dispatch on, then\n",
       "  nothing happens. I just return the original Jacobian matrix and ptcsol will\n",
       "  use backslash to compute the Newton step.\n",
       "\n",
       "  I know that this is probably not optimal in your situation, so it is good to\n",
       "  pick something else, like jfact = lu.\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",
       "  jknowsdt: default = false\n",
       "\n",
       "  Set this to true if your Jacobian evaluation function returns F'(x) + (1/dt)\n",
       "  I. You'll also need to follow the rules above for the Jacobian evaluation\n",
       "  function. I do not recommend this and if your Jacobian is anything other\n",
       "  than a matrix I can't promise anything. I've tested this for matrix outputs\n",
       "  only.\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\n",
       "\n",
       "  Unlike nsol, nsoli, or even ptcsoli, ptcsol has a fixed cost per iteration\n",
       "  of one function, one Jacobian, and one Factorization. Hence iteration\n",
       "  statistics are not interesting and not in the output.\n",
       "\n",
       "  idid=true if the iteration succeeded and false if not.\n",
       "\n",
       "  errcode = 0 if the iteration succeeded = -1 if the initial iterate satisfies\n",
       "  the termination criteria = 10 if no convergence after maxit iterations\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.\n",
       "\n",
       "\u001b[36m  julia> using SIAMFANLEquations.TestProblems\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> n=63; maxit=1000; delta = 0.01; lambda = 20.0;\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> bdata = beaminit(n, 0.0, lambda); x = bdata.x;\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> u0 = x .* (1.0 .- x) .* (2.0 .- x);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> u0 .*= exp.(-10.0 * u0);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> FS = copy(u0); FPS = copy(bdata.D2);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> pout = ptcsol( FBeam!, u0, FS, FPS, BeamJ!; \u001b[39m\n",
       "\u001b[36m   rtol = 1.e-10, pdata = bdata, delta0 = delta, maxit = maxit);\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 Array{Float64,2}:\u001b[39m\n",
       "\u001b[36m   6.31230e+01  9.75412e-01\u001b[39m\n",
       "\u001b[36m   7.52624e+00  8.35295e-02\u001b[39m\n",
       "\u001b[36m   8.31545e+00  6.58797e-04\u001b[39m\n",
       "\u001b[36m   3.15455e+01  4.12697e-08\u001b[39m\n",
       "\u001b[36m   3.66566e+01  6.29295e-12\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> # We get the nonnegative steady state.\u001b[39m\n",
       "\u001b[36m         maximum(pout.solution)\u001b[39m\n",
       "\u001b[36m  2.19086e+00\u001b[39m"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?ptcsol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### More on the Buckling Beam\n",
    "\n",
    "The residual function for the beam uses precomputed data to store the discrete second derivative and the bifurcation parameter $\\lambda$. These data are smaller and far less complex that the data for the H-equation and it's worthwhile to look at the functions. These are part of the __TestProblems__ submodule for __SIAMFANLEquations.jl__. We will not cover all the details here, but will give you enough to see how precomputed data is useful. The file is __TestProblems/Systems/FBeam!.jl__.\n",
    "\n",
    "First the function itself\n",
    "```Julia\n",
    "\"\"\"\n",
    "FBeam!(FV, U, bdata)\n",
    "Function evaluation for PTC example.\n",
    "F(u) = -u'' - lambda sin(u)\n",
    "\"\"\"\n",
    "function FBeam!(FV, U, bdata)\n",
    "    D2 = bdata.D2\n",
    "    lambda = bdata.lambda\n",
    "    su = lambda * sin.(U)\n",
    "    FV .= (D2 * U - su)\n",
    "end\n",
    "```\n",
    "Note how the function ```FBeam``` harvests $\\lambda$ and the discrete second derivative ```D2``` from the precomputed data ```bdata```. In this case ```bdata``` is a __named tuple__ which we create with the\n",
    "__beaminit__ function.\n",
    "\n",
    "```Julia\n",
    "\"\"\"\n",
    "beaminit(n,dt,lambda=20.0)\n",
    "\n",
    "Set up the beam problem with n interior grid points.\n",
    "\"\"\"\n",
    "\n",
    "function beaminit(n, dt, lambda = 20.0)\n",
    "    D2 = Lap1d(n)\n",
    "    dx = 1.0 / (n + 1)\n",
    "    x = collect(dx:dx:1.0-dx)\n",
    "    UN = zeros(size(x))\n",
    "    bdata = (D2 = D2, x = x, dx = dx, dt = dt, \n",
    "             lambda = lambda, UN = UN)\n",
    "end\n",
    "```\n",
    "The bdata structure has more that just ```D2``` and $\\lambda$. There are also data for the solver and the construction of the Jacobian.\n",
    "\n",
    "Finally, the function __Lap1D__ computes ```D2``` as a tridiagonal matrix. This is part of the LinearAlgebra package that is part of Julia.Base.\n",
    "\n",
    "```Julia\n",
    "\"\"\"\n",
    "Lap1d(n)\n",
    "\n",
    "returns -d^2/dx^2 on [0,1] zero BC\n",
    "\"\"\"\n",
    "function Lap1d(n)\n",
    "    dx = 1 / (n + 1)\n",
    "    d = 2.0 * ones(n)\n",
    "    sup = -ones(n - 1)\n",
    "    slo = -ones(n - 1)\n",
    "    D2 = Tridiagonal(slo, d, sup)\n",
    "    D2 = D2 / (dx * dx)\n",
    "    return D2\n",
    "end\n",
    "\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.9: Projects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More on the H-equation\n",
    "\n",
    "- Run the H-equation tests with $c=.99$ and larger dimensions ($2^k$ for $k=12, 13, 14$). Are the values of ```sham``` and ```resdec``` optimal for this case? How, for example, is ```resdec = Inf``` (your author's favorite value). To see why I like ```sham=Inf``` have a look at\n",
    "<cite data-cite=\"brent\"><a href=\"siamfa.html#brent\">(Bre73)</cite>.\n",
    "\n",
    "- What happens with $c=1$? Is the chord method a good idea? \n",
    "Don't look at\n",
    "<cite data-cite=\"ctk:sirev20\"><a href=\"siamfa.html#ctk:sirev20\">(Kel22)</cite>,\n",
    "<cite data-cite=\"ctk:n1\"><a href=\"siamfa.html#ctk:n1\">(DK80)</cite>, or\n",
    "<cite data-cite=\"ctk:chord\"><a href=\"siamfa.html#ctk:chord\">(DK83)</cite> before trying to figure things out on your own.\n",
    "\n",
    "- Increase the dimension as much as you can and compare single and double precision linear algebra.\n",
    "    \n",
    "- Duplicate the table on page 125 of <cite data-cite=\"chand\"><a href=\"siamfa.html#chand\">(Cha60)</cite>.    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__(Not trivial!)__ Try using __automatic differentiation__ to compute the Jacobian for the examples in this chapter and use the results to write a Jacobian evaluation function to pass to __nsol.jl__.\n",
    "You will first have to decide how you want to do that. Two mature packages are [ForwardDiff.jl](#https://github.com/JuliaDiff/ForwardDiff.jl)\n",
    "<cite data-cite=\"forwarddiff\"><a href=\"siamfa.html#forwarddiff\">(RPL16)</cite>\n",
    "and [Zygote.jl](#https://github.com/FluxML/Zygote.jl) \n",
    "<cite data-cite=\"zygote\"><a href=\"siamfa.html#zygote\">(CoRR)</cite>.\n",
    "Books like\n",
    "<cite data-cite=\"grautodiff\"><a href=\"siamfa.html#grautodiff\">(Gri00)</cite> explain the\n",
    "algorithmic differences between these packages.\n",
    "    \n",
    "How does the performance compare to an analytic or forward difference Jacobian? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nested Iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving a differential or integral equation\n",
    "by __nested iteration__ or __grid sequencing__ means resolving\n",
    "the rough features of the solution of a differential or\n",
    "integral equation on a coarse mesh, interpolating the\n",
    "solution to a finer mesh, resolving on the finer mesh, and\n",
    "then repeating the process until you have a solution on a target\n",
    "mesh.\n",
    "\n",
    "Apply this idea to some of the examples in the text, using piecewise\n",
    "linear interpolation to move from coarse to fine meshes. If the\n",
    "discretization is second-order accurate and you halve the mesh\n",
    "width at each level, how should you terminate the solver at\n",
    "each level? What kind of iteration statistics would tell you that\n",
    "you've done a satisfactory job?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse Differencing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write Jacobian evaluation functions using sparse differencing for the boundary value problem and the buckling beam examples. You can try, for example, the methods from \n",
    "<cite data-cite=\"curtispr\"><a href=\"siamfa.html#curtisor\">(CPR74)</cite> or\n",
    "<cite data-cite=\"colmore\"><a href=\"siamfa.html#colmore\">(CM83)</cite>.  <cite data-cite=\"ctk:newton\"><a href=\"siamfa.html#ctk:newton\">(Kel03)</cite> has a Matlab code for banded matrices that, while convertible to Julia, needs some work to explicitly store the Jacobian as a BandedMatrix so\n",
    "[BandedMatrices.jl](https://github.com/JuliaMatrices/BandedMatrices.jl) will solve the linear system efficiently.\n",
    "Another alternative is to use __sparsefdiff__ \n",
    "<cite data-cite=\"sparsediff\"><a href=\"siamfa.html#sparsediff\">(RMGH20)</cite> from\n",
    "[SparseDiffTools.jl](#https://github.com/JuliaDiff/SparseDiffTools.jl). \n",
    "    \n",
    "No matter what you do, you'll need to think about fill-in, symbolic factorization in the general case, and storage. \n",
    "    \n",
    "Have fun!    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  ```lu!``` for Sparse Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This project is for Julia 1.5 and higher. The option ```jfact=nofact``` in __nsol.jl__ is there so you can build and manage your own factorization inside your Jacobian evaluation function. Use this feature to store a factorization of the Jacobian for the PDE problem and reuse that storage rather than reallocating with every call to lu. What are the savings in time and allocations? You'll need to figure out how to use ```lu!``` for sparse matrices to do this project."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does $\\delta_0$ depend on $\\delta_x$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vary $\\delta_x$ and $\\delta_0$ in the buckling beam example. For \n",
    "the fixed $\\delta_x = 1/64$ in the example, explore the effects of\n",
    "increasing $\\delta_0$. Can you converge to the ``wrong'' (\\ie negative)\n",
    "solution in this way? What happens if you fix $\\delta_0$ and reduce\n",
    "$\\delta_x$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next notebook = [Chapter 3: Newton-Krylov Methods](SIAMFANLCh3.ipynb)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": []
  }
 ],
 "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
}