#include <sys/stat.h>
#include <fstream>
#include <iomanip>
#include <string>
#include <iostream>
#include <octave/oct.h>
#include <octave/CColVector.h>
#include <octave/dColVector.h>
#include <octave/dMatrix.h>
#include <octave/EIG.h>
#include <octave/dbleSVD.h>
#include "channelflow/flowfield.h"
#include "channelflow/vector.h"
#include "channelflow/chebyshev.h"
#include "channelflow/tausolver.h"
#include "channelflow/dns.h"
#include "channelflow/symmetry.h"
#include "channelflow/periodicfunc.h"
#include "octutils.h"

// This program is an implementation of Divakar Viswanath's Newton-hookstep
// GMRES algorithm for computing invariant solutions of Navier-Stokes. The
// basic idea is to find solutions of sigma f^T(u) - u = 0 where f^T(u) is the
// time-T CFD integration of the fluid velocity field u and sigma is a symmetry
// operation. T and sigma can be held fixed or varied. Viswanath's algorithm
// combines a Newton-hookstep trust-region minimization of ||sigma f^T(u) - u||^2
// with iterative GMRES solution of the Newton-step equations.

// In this program, the details of GMRES algorithm are modeled on Lecture 35
// of Trefethen and Bau's "Numerical Linear Algebra". The details of the hookstep
// algorithm are based on section 6.4.2 of Dennis and Schnabel's "Numerical Methods
// for Unconstrained Optimization and Nonlinear Equations."
//
// See end of file for a detailed description of the algorithm and notation.
//
// References:
//
// Divakar Viswanath, "Recurrent motions within plane Couette turbulence",
// J. Fluid Mech. 580 (2007), http://arxiv.org/abs/physics/0604062
//
// Lloyd N. Trefethen and David Bau, "Numerical Linear Algebra", SIAM,
// Philadelphia, 1997.
//
// Dennis and Schnabel, "Numerical Methods for Unconstrained Optimization
// and Nonlinear Equations", Prentice-Hall Series in Computational Mathematics,
// Englewood Cliffs, New Jersey, 1983.
//
// John Gibson, Wed Oct 10 16:10:46 EDT 2007

using namespace std;
using namespace channelflow;

enum BasisType {L2Basis, AdHocBasis};

Real V2Norm2(const FlowField& u); // L2Norm2(field2vector(u));
Real V2Norm(const FlowField& u);
Real V2IP(const FlowField& u, const FlowField& v);

Real Norm2(const FlowField& u, BasisType norm);
Real Norm(const FlowField& u, BasisType norm);
Real IP(const FlowField& u, const FlowField& v, BasisType norm);

// f^T(u) = time-T DNS integration of u
void f(const FlowField& u, Real T, FlowField& f_u, Real Reynolds, DNSFlags& flags,
       TimeStep& dt, int& fcount, Real& CFL);

// f^(n.dt)(u) == time-(n.dt) DNS integration of u
void f(const FlowField& u, int n, Real dt, FlowField& f_u, Real Reynolds,
       DNSFlags& flags);

// G(u,T) = (tau f^T(u) - u) for orbits and (tau f^T(u) - u)/T for eqbs
void G(const FlowField& u, Real T, const FieldSymmetry& tau, FlowField& G_u,
       Real Reynolds, DNSFlags& flags, TimeStep& dt, bool Tnormalize, int& fcount, Real& CFL);

void DG(const FlowField& u, const FlowField& du, Real T, Real dT,
	const FieldSymmetry& tau,  const FieldSymmetry& dtau,
	const FlowField& GuT, FlowField& DG_du, Real Reynolds, DNSFlags& flags,
	TimeStep& dt, bool Tnormalize, Real eps, bool centerdiff, int& fcount, Real& CFL);

//enum HookstepPhase {ConstantDelta, ReducingDelta, IncreasingDelta, Finished};
//ostream& operator<<(ostream& os, HookstepPhase p);

//enum ResidualImprovement {Unacceptable, Poor, Ok, Good, Accurate, NegaCurve};
//ostream& operator<<(ostream& os, ResidualImprovement i);

//Real adjustLambda(Real lambda, Real lambdaMin, Real lambdaMax);
//Real adjustDelta(Real delta, Real rescale, Real deltaMin, Real deltaMax);

int main(int argc, char* argv[]) {

  string purpose("Newton-Krylov-hookstep search for (relative) periodic orbit in U_S");
  ArgList args(argc, argv, purpose);

  const bool eqb       = args.getflag("-eqb",  "--equilibrium",   "search for equilibrium or relative equilibrium (trav wave)");
  const bool orb       = args.getflag("-orb",  "--periodicorbit", "search for periodic orbit or relative periodic orbit");
  const bool relative  = args.getflag("-rel",  "--relative",      "search for relative periodic orbit or relative equilibrium");
  /***/ Real T         = args.getreal("-T",   "--maptime",    10.0,  "initial guess for orbit period or time of eqb/reqb map f^T(u)");

  const Real epsSearch = args.getreal("-es",  "--epsSearch",   1e-13, "stop search if L2Norm(s f^T(u) - u) < epsEQB");
  const Real epsKrylov = args.getreal("-ek",  "--epsKrylov",   1e-14, "min. condition # of Krylov vectors");
  const Real epsDu     = args.getreal("-edu", "--epsDuLinear", 1e-7,  "relative size of du to u in linearization");
  const Real epsDt     = args.getreal("-edt", "--epsDtLinear", 1e-7,  "size of dT in linearization of f^T about T");
  const Real epsGMRES  = args.getreal("-eg",  "--epsGMRES",    1e-3,  "stop GMRES iteration when Ax=b residual is < this");
  const Real epsGMRESf = args.getreal("-egf",  "--epsGMRESfinal",  0.05,  "accept final GMRES iterate if residual is < this");
  const bool centdiff = args.getflag("-cd", "--centerdiff", "centered differencing to estimate differentials");

  const int  Nnewton   = args.getint( "-Nn", "--Nnewton",   20,    "max number of Newton steps ");
  const int  Ngmres    = args.getint( "-Ng", "--Ngmres",    40,    "max number of GMRES iterations per restart");
  //const int  Nrestart  = args.getint( "-Nr", "--Nrestart",  5,     "max number of GMRES restarts per Newton step");
  const int  Nhook     = args.getint( "-Nh", "--Nhook",     20,    "max number of hookstep iterations per Newton");

  /***/ Real delta     = args.getreal("-d",    "--delta",      0.01,  "initial radius of trust region");
  const Real deltaMin  = args.getreal("-dmin", "--deltaMin",   1e-12, "stop if radius of trust region gets this small");
  const Real deltaMax  = args.getreal("-dmax", "--deltaMax",   0.1,   "maximum radius of trust region");
  const Real deltaFuzz = args.getreal("-df",   "--deltaFuzz",  0.01,  "accept steps within (1+/-deltaFuzz)*delta");
  const Real lambdaMin = args.getreal("-lmin",   "--lambdaMin", 0.2,  "minimum delta shrink rate");
  const Real lambdaMax = args.getreal("-lmax",   "--lambdaMax", 1.5,  "maximum delta expansion rate");
  const Real lambdaRequiredReduction = 0.5; // when reducing delta, reduce by at least this factor.
  const Real improvReq = args.getreal("-irq",  "--improveReq", 1e-3,  "reduce delta and recompute hookstep if improvement "
				      " is worse than this fraction of what we'd expect from gradient");
  const Real improvOk  = args.getreal("-iok",  "--improveOk",   0.10, "accept step and keep same delta if improvement "
				      "is better than this fraction of quadratic model");
  const Real improvGood= args.getreal("-igd",  "--improveGood", 0.75, "accept step and increase delta if improvement "
				      "is better than this fraction of quadratic model");
  const Real improvAcc = args.getreal("-iac",  "--improveAcc",  0.10, "recompute hookstep with larger trust region if "
				      "improvement is within this fraction quadratic prediction.");
  //const bool pureNewton= args.getflag("-pn",   "--pureNewton",        "do pure Newton search instead of hookstep");

  const Real TscaleRel = args.getreal("-Tsc", "--Tscale",     20,    "scale time in hookstep: |T| = Ts*|u|/L2Norm(du/dt)");
  const bool dudtortho = args.getbool("-tc",  "--dudtConstr", true,  "require orthogonality of step to dudt");
  const Real orthoscale= args.getreal("-os",  "--orthoScale", 1,     "rescale orthogonality constraint");

  const string symname = args.getstr ("-sf", "--symmfile", "", "file containing T,su,sx,sy,sz,ax,az");
  const bool sx        = args.getflag("-sx", "--x-sign", "change u,x sign");
  const bool sy        = args.getflag("-sy", "--y-sign", "change v,y sign");
  const bool sz        = args.getflag("-sz", "--z-sign", "change w,z sign");
  /***/ Real ax        = args.getreal("-ax", "--axshift", 0.0, "translate by ax*Lx");
  /***/ Real az        = args.getreal("-az", "--azshift", 0.0, "translate by az*Lz");
  const bool anti      = args.getflag("-a",  "--anti",         "antisymmetry instead of symmetry");
  const bool Rscale    = args.getreal("-rs", "--relativeScale", 1, "scale relative-search variables by this factor");
  const bool dudxortho = args.getbool("-xc", "--dudxConstraint", true,"require orthogonality of step to dudx and dudz");

  const bool l2basis   = args.getflag("-b",  "--l2basis", "use an explicit L2 basis");
  const bool s1        = args.getflag("-s1", "--s1", "restrict to s1-symmetric subspace");
  const bool s2        = args.getflag("-s2", "--s2", "restrict to s2-symmetric subspace");
  const bool s3        = args.getflag("-s3", "--s3", "restrict to s3-symmetric subspace");

  const Real Reynolds  = args.getreal("-R", "--Reynolds", 400, "Reynolds number");
  const bool vardt     = args.getflag("-vdt","--variableDt", "adjust dt to keep CFLmin<=CFL<CFLmax");
  const Real dtarg     = args.getreal("-dt","--timestep",  0.03125, "(initial) integration timestep");
  const Real dtmin     = args.getreal("-dtmin", "--dtmin", 0.001, "minimum timestep");
  const Real dtmax     = args.getreal("-dtmax", "--dtmax", 0.04,  "maximum timestep");
  const Real dTCFL     = args.getreal("-dTCFL", "--dTCFL", 1.00,  "check CFL # at this interval");
  const Real CFLmin    = args.getreal("-CFLmin", "--CFLmin", 0.40, "minimum CFL number");
  const Real CFLmax    = args.getreal("-CFLmax", "--CFLmax", 0.60, "maximum CFL number");
  const string stepstr = args.getstr("-ts", "--timestepping", "SBDF3", "timestep algorithm: See README.");
  const string nonlstr = args.getstr("-nl", "--nonlinearity", "Rotational", "nonlinearity method: See README.");

  const int  coutprec  = args.getint("-cp", "--coutprec", 6, "precison for std output");
  const bool pausing   = args.getflag("-p", "--pausing", "pause between Newton steps");
  const string outdir  = args.getpath("-o", "--outdir", "./", "output directory");
  const string uname   = args.getstr(1, "<flowfield>", "initial guess for Newton search");

  args.check();
  if (!(eqb^orb)) {
    cerr << "Please choose either -eqb or -orb option to search for (relative) equilibrium or (relative) periodic orbit" << endl;
    exit(1);
  }
  args.save(outdir);

  const BasisType nrm = l2basis ? L2Basis : AdHocBasis;
  const int w = coutprec + 7;
  const bool Tsearch    = orb ? true : false;
  const bool Tnormalize = orb ? false : true;
  const bool Rsearch  = relative ? true : false;

  cout << "Working directory == " << pwd() << endl;
  cout << "Command-line args == ";
  for (int i=0; i<argc; ++i) cout << argv[i] << ' '; cout << endl;

  DNSFlags flags;
  flags.baseflow     = PlaneCouette;
  flags.timestepping = s2stepmethod(stepstr);
  flags.dealiasing   = DealiasXZ;
  flags.nonlinearity = s2nonlmethod(nonlstr);
  flags.constraint   = PressureGradient;
  flags.dPdx         = 0.0;
  flags.verbosity    = Silent;
  cout << "DNSFlags == " << flags << endl << endl;

  int fcount_newton = 0;
  int fcount_hookstep = 0;

  TimeStep dt(dtarg, dtmin, dtmax, dTCFL, CFLmin, CFLmax, vardt);
  cout << "dt    == " << dt << endl;

  FlowField u(uname);

  const int Nx = u.Nx();
  const int Ny = u.Ny();
  const int Nz = u.Nz();
  const Real Lx = u.Lx();
  const Real ya = u.a();
  const Real yb = u.b();
  const Real Lz = u.Lz();
  const int kxmax = u.kxmaxDealiased();
  const int kzmax = u.kzmaxDealiased();

  vector<RealProfileNG> basis;
  if (l2basis || s1 || s2 || s3) {
    cout << "constructing basis set (takes a long time)..." << endl;
    basis = realBasisNG( Ny, kxmax, kzmax, Lx, Lz, ya, yb);

    vector<FieldSymmetry> symms;
    if (s1) symms.push_back(FieldSymmetry( 1, 1,-1, 0.5, 0.0));
    if (s2) symms.push_back(FieldSymmetry(-1,-1, 1, 0.5, 0.5));
    if (s3) symms.push_back(FieldSymmetry(-1,-1,-1, 0.0, 0.5));
    if (s1 || s2 || s3)
      selectSymmetries(basis, symms);
    orthonormalize(basis);
  }
  ColumnVector sN;      // unconstrained vector representation of eqb u
  field2vector(u, sN, basis);

  // Mcon == number of constraints
  // Neqn == number of equations (total)
  // Nunk == number of unknowns
  const int uunk = sN.length();                                  // # of variables for u unknonwn
  const int Tunk = Tsearch ?  uunk : -1;                         // index for T unknown
  const int xunk = Rsearch ? (uunk + (Tsearch ? 1 : 0)) : -1;    // index for x-shift unknown
  const int zunk = Rsearch ? (xunk + 1) : -1;                    // index for z-shift unknown
  const int Nunk = uunk + (Tsearch ? 1 : 0) + (Rsearch ? 2 : 0);
  const int xeqn = xunk;
  const int zeqn = zunk;
  const int Teqn = Tunk;
  const int Neqn = Nunk;
  const int Ngmres2 = (Nunk < Ngmres) ? Nunk : Ngmres;

  cout << Neqn << " equations" << endl;
  cout << Nunk << " unknowns" << endl;
  if (Tsearch) {
    cout << "Teqn  == " << Teqn << endl;
    cout << "Tunk  == " << Tunk << endl;
  }
  if (Rsearch) {
    cout << "xeqn  == " << xeqn << endl;
    cout << "xunk  == " << xunk << endl;
    cout << "zeqn  == " << zeqn << endl;
    cout << "zunk  == " << zunk << endl;
  }

  sN = ColumnVector(Nunk);

  // ==============================================================
  // Initialize Newton iteration

  FlowField      du(Nx, Ny, Nz, 3, Lx, Lz, ya, yb); // temporary for increments in orbit velocity field
  FlowField      Gu(Nx, Ny, Nz, 3, Lx, Lz, ya, yb); // G(u,T) = f^T(u) - u
  FlowField   DG_du(Nx, Ny, Nz, 3, Lx, Lz, ya, yb); // "DG*du" = D[u] G(u,T) du = 1/e (G(u+e.du,T)
  FlowField       p(Nx, Ny, Nz, 1, Lx, Lz, ya, yb); // pressure field for time integration
  FlowField      uS(Nx, Ny, Nz, 3, Lx, Lz, ya, yb); // stores uN and uH, the Newton and Hookstep velocity field
  FlowField      GS(Nx, Ny, Nz, 3, Lx, Lz, ya, yb); // stores GN and GH, GS = G(uG,tG)

  Real dT  = 0.0;                                    // Newton step for orbit period T
  Real tS = 0.0;
  Real dax = 0.0;                                    // Newton step for x phase shift
  Real daz = 0.0;                                    // Newton step for z phase shift

  FieldSymmetry tau(sx,sy,sz,ax,az,anti);
  if (symname.length() != 0)
    load(symname, T, tau);
  cout << "Initial integration time T == " << T << endl;
  cout << "Initial symmetry tau       == " << tau << endl;

  FieldSymmetry dtau; // identity, for now
  FieldSymmetry tauS(tau);

  ChebyCoeff U(Ny, ya, yb, Spectral);
  U[1] = 1;

  Real ruT = 0.0;        // Dennis & Schnabel residual r(u,T) = 1/2 ||G(u,T)||^2 == 1/2 VNorm2(G(u,T))
  Real guT = 0.0;        // g(u,T) = L2Norm(G(u,T))
  Real init_ruT = 0.0;
  Real init_guT = 0.0;
  Real prev_ruT = 0.0;
  Real prev_guT = 0.0;
  Real gmres_residual = 0.0;
  Real snorm   = 0.0;   // norm of hookstep
  Real sNnorm  = 0.0;   // norm of newton step
  Real dunorm  = 0.0;
  Real unorm   = 0.0;
  Real CFL = 0.0;
  bool stuck = false;

  // Determine which of four half-cell translations tau minimizes 1/2 || tau f^T(u) - u ||^2
  //cout << "Determining which of four half-cell translations tau minimizes residual" << endl;
  cout << "Search for solutions of tau f^T(u) - u == 0 with " << endl;
  cout << "  tau == " << tau << endl;

  ofstream osconv((outdir + "convergence.asc").c_str());
  osconv.setf(ios::left);
  osconv << setw(14) << "% L2Norm(G)" << setw(14) << "r(u,T)" << setw(14) << "delta";
  if (Tsearch)
    osconv << setw(14) << "dT";
  if (Rsearch) {
    osconv << setw(14) << "dax";
    osconv << setw(14) << "daz";
  }
  osconv << setw(14) << "L2Norm(du)" << setw(14) << "L2Norm(sN)" << setw(14) << "L2Norm(sH)"
	 << setw(14) << "|Ax-b|/|b|" << setw(10) << "ftotal" << setw(10) << "fnewt"
	 << setw(10) << "fhook" << endl;

  cout << setprecision(coutprec);

  cout << "Computing G(u,T) = " << ((Tnormalize) ? "1/T (tau f(u,T) - u)" : "tau f(u,T) - u") << endl;
  G(u, T, tau, Gu, Reynolds, flags, dt, Tnormalize, fcount_hookstep, CFL);
  cout << endl;

  for (int newtonStep = 0; newtonStep<=Nnewton; ++newtonStep) {

    cout << "========================================================" << endl;
    cout << "Newton iteration number " << newtonStep << endl;

    // Compute quantities that will be used multiple times in GMRES loop
    // Gu     = (tau f^T(u) - u)  or  1/T (tau f^T(u) - u)
    // dfuTdT = 1/(dt) (f^(T+dt)(u) - f^(T)(u))
    // dudt   = 1/(dt) (f^(dt)(u) - u)
    //        = tangent vector to flow
    // edudx  = tangent vector of x translation
    // edudz  = tangent vector of z translation

    dunorm = L2Norm(du);
    unorm = L2Norm(u);
    ruT = 0.5*Norm2(Gu, nrm);
    guT = L2Norm(Gu);

    cout << "Computing normalized dudt, edudx, edudz" << endl;
    // Calculate three unit vectors to which du must be orthogonal
    FlowField edudx = xdiff(u);
    FlowField edudz = zdiff(u);
    FlowField edudt(u);
    f(u, 1, epsDt, edudt, Reynolds, flags);   // using dt=epsDt timestep
    edudt -= u;
    edudt *= 1.0/epsDt;

    // Set Tscale so that |dT|*Tscale = TscaleRel |du|
    // Tscale = TscaleRel |du|/|dt|
    // O(t)*Tscale = O(u) => Tscale = u/t
    //const Real Tscale = TscaleRel*V2Norm(edudt);
    const Real Tscale = TscaleRel;
    cout << "delta  == " << delta << endl;
    cout << "Tscale == " << Tscale << endl;
    edudt *= 1.0/L2Norm(edudt);
    edudx *= 1.0/L2Norm(edudx);
    edudz *= 1.0/L2Norm(edudz);

    ruT = 0.5*Norm2(Gu, nrm);
    guT = L2Norm(Gu);

    if (newtonStep == 0) {
      init_guT = guT;
      prev_guT = guT;
      init_ruT = ruT;
      prev_ruT = ruT;
    }

    FlowField uU(u);
    uU += U;

    cout << "Current state of Newton iteration:" << endl;
    //cout << "dissip (u)     == " << dissipation(u) << endl;
    //cout << "forcing(u)     == " << forcing(u) << endl;
    //cout << "energy (u)     == " << 0.5*square(unorm) << endl;
    //cout << "dissip (u+U)   == " << dissipation(uU) << endl;
    //cout << "forcing(u+U)   == " << forcing(uU) << endl;
    //cout << "energy (u+U)   == " << 0.5*L2Norm2(uU) << endl;
    cout << "fcount_newton  == " << fcount_newton << endl;
    cout << "fcount_hookstep== " << fcount_hookstep << endl;
    cout << "    L2Norm(u)  == " << unorm << endl;
    cout << "   L2Norm(du)  == " << dunorm <<endl;
    cout << "   L2Norm(sH)  == " << snorm <<endl;
    cout << "           T   == " << T << endl;
    cout << "          dT   == " << dT << endl;
    cout << "         tau   == " << tau << endl;
    cout << "        dtau   == " << dtau << endl;
    cout << "L2Dist(u,u0)   == " << L2Dist(u, FlowField(uname)) << endl;
    cout << "guT == L2Norm(G(u,T)) : " << endl;
    cout << "  initial  guT == " << init_guT << endl;
    cout << "  previous guT == " << prev_guT << endl;
    cout << "  current  guT == " << guT << endl;
    cout << "ruT == 1/2 Norm2(G(u,T)) : " << endl;
    cout << "  initial  ruT == " << init_ruT << endl;
    cout << "  previous ruT == " << prev_ruT << endl;
    cout << "  current  ruT == " << ruT << endl;
    cout << "1/2 L2Norm2(G,u,T) : " << endl;
    cout << "               == " << 0.5*square(guT) << endl;
    cout << "CFL == " << CFL << endl;

    osconv.setf(ios::left); osconv << setw(14) << guT;
    osconv.setf(ios::left); osconv << setw(14) << ruT;
    osconv.setf(ios::left); osconv << setw(14) << delta;
    if (Tsearch) {
      osconv.setf(ios::left); osconv << setw(14) << dT;
    }
    if (Rsearch) {
      osconv.setf(ios::left); osconv << setw(14) << dtau.ax();
      osconv.setf(ios::left); osconv << setw(14) << dtau.az();
    }
    osconv.setf(ios::left); osconv << setw(14) << dunorm;
    osconv.setf(ios::left); osconv << setw(14) << sNnorm;
    osconv.setf(ios::left); osconv << setw(14) << snorm;
    osconv.setf(ios::left); osconv << setw(14) << gmres_residual;
    osconv.setf(ios::left); osconv << setw(10) << fcount_newton + fcount_hookstep;
    osconv.setf(ios::left); osconv << setw(10) << fcount_newton;
    osconv.setf(ios::left); osconv << setw(10) << fcount_hookstep;
    osconv << endl;

    u.binarySave(outdir + "unewt" + i2s(newtonStep));
    u.binarySave(outdir + "ubest");
    if (Tsearch || Rsearch) {
      save(outdir + "taunewt" + i2s(newtonStep), T, tau);
      save(outdir + "taubest", T, tau);
    }

    //for (int p=0; p<P; ++p)
    //stepos << L2IP(u, e[p]) << ' ';
    //stepos << endl;

    if (guT < epsSearch) {
      cout << "Newton search converged. Breaking." << endl;
      cout << "L2Norm(G(u,T))) == " << guT << " < "
	   << epsSearch << " == " << epsSearch << endl;
      //u.binarySave(outdir + "usoln");
      //if (Tsearch || Rsearch)
      //save(outdir + "tausoln", T, tau);
      break;
    }
    else if (stuck) {
      cout << "Newton search is stuck. Breaking." << endl;
      break;
    }
    else if (newtonStep == Nnewton) {
      cout << "Reached maximum number of Newton steps. Breaking." << endl;
      break;
    }
    if (pausing)
      cfpause();


    prev_guT = guT;
    prev_ruT = ruT;

    // Set up RHS vector b = (-field2vector(G(u)), 0, 0, 0);
    // (0, 0, 0) is the RHS of eqn enforcing orthog. to (edudt, edudx, edudz)
    ColumnVector b(Neqn);
    setToZero(b);
    field2vector(Gu, b, basis);
    if (Tsearch)
      b(Teqn) = 0;
    if (Rsearch) {
      b(xeqn) = 0;
      b(zeqn) = 0;
    }
    b *= -1.0;

    GMRES gmres(b, Ngmres2, epsKrylov);

    // ===============================================================
    // GMRES iteration to solve D[u]G(u,t) du + D[t]G(u,t) dt = -G(u) for du,dt
    for (int n=0; n<Ngmres2; ++n) {

      cout << "Newt,GMRES == " << newtonStep << ',' << n << ", " << flush;

      // Compute v = Ab in Arnoldi iteration terms, where b is Q.column(n)
      // In Navier-Stokes terms, the main quantity to compute is
      //   L v = D[u] G(u,T) du + D[T] G(u,T) dT

      // Use approximations
      // D[u]G(u,T) du = 1/e (G(u + e du, T) - G(u)) for e << 1
      // D[T]G(u,T) dT = 1/e (G(u, T + e dT) - G(u)) for e << 1

      ColumnVector q = gmres.testVector();
      vector2field(q, du, basis);
      dT  = Tsearch ? q(Tunk)*Tscale : 0.0;
      dax = Rsearch ? q(xunk)*Rscale : 0.0;
      daz = Rsearch ? q(zunk)*Rscale : 0.0;
      dtau = FieldSymmetry(1,1,1,dax,daz);

      dunorm = L2Norm(du);
      //cout << "L2Norm(du) == " << L2Norm(du) << endl;
      //cout << "      dax  == " << dax << endl;
      //cout << "      daz  == " << daz << endl;
      //cout << "     dtau  == " << dtau << endl;
      //cout << "       dT  == " << dT << endl;

      // Compute D[u]G(u,T) du = 1/e (G(u+e.du,T) - G(u,T))
      //                       = 1/e (tau f^T(u+e.du) - u+e.du - (tau f^T(u) - u)
      //                       = 1/e (e.du + tau(f^T(u) - f^T(u+e.du))
      // Hmmm. could be better numerically to use this last form, but I'll use first.
      //cout << "Computing D[u]G(u,T) du..." << flush;

      DG(u, du, T, dT, tau, dtau, Gu, DG_du, Reynolds, flags, dt, Tnormalize, epsDu,  centdiff, fcount_newton, CFL);

      ColumnVector Lq(Neqn);
      field2vector(DG_du, Lq, basis);
      if (Tsearch)
	Lq(Teqn) = dudtortho ? orthoscale*L2InnerProduct(du, edudt) : 0.0;
      if (Rsearch) {
	Lq(xeqn) = dudxortho ? orthoscale*L2InnerProduct(du, edudx) : 0.0;
	Lq(zeqn) = dudxortho ? orthoscale*L2InnerProduct(du, edudz) : 0.0;
      }

      gmres.iterate(Lq);
      gmres_residual = gmres.residual();

      cout << " res == " << gmres_residual << endl;

      if (gmres_residual < epsGMRES) {
        cout << "GMRES converged. Breaking." <<endl;
	sN = gmres.solution();
	break;
      }
      else if (n==Ngmres-1 && gmres_residual < epsGMRESf ) {
        cout << "GMRES has not converged, but the final iterate is acceptable. Breaking." <<endl;
	sN = gmres.solution();
      }
      else if (n==Ngmres-1) {
	cout << "GMRES failed to converge. Exiting." << endl;
	exit(1);
      }
    } // end GMRES iteration

    sNnorm = L2Norm(sN);

    // Diagnostic
    //vector2field(sN, du, basis);
    //du.binarySave("duN" + i2s(newtonStep));

    // ==================================================================
    // Hookstep algorithm

    cout << "------------------------------------------------" << endl;
    cout << "Beginning hookstep calculations." << endl;

    int Nk = gmres.n(); // Krylov dimension
    Matrix Hn = gmres.Hn();
    Matrix Qn = gmres.Qn();
    Matrix Qn1t = gmres.Qn1().transpose(); // wasteful of memory
    SVD svd(Hn, SVD::std);

    Matrix V = svd.right_singular_matrix();
    Matrix Ut = svd.left_singular_matrix().transpose();
    DiagMatrix D = svd.singular_values();

    ColumnVector btmp = Qn1t*b;
    ColumnVector bh = Ut*btmp; // b hat

    int hookcount = 0;

    bool have_backup  = false;
    Real backup_snorm = 0.0;
    Real backup_delta = 0.0;
    Real backup_rS    = 0.0;
    Real backup_dT    = 0;
    FieldSymmetry backup_dtau;

    FlowField backup_du;
    FlowField backup_GS;

    ColumnVector sH;

    //bool delta_has_decreased = false;
    Real deltaMaxLocal = deltaMax;

    // Find a good trust region and the optimal hookstep in it.
    // Quality of trust region
    while (true) {

      cout << "-------------------------------------------" << endl;
      cout << "Newton, hookstep number " << newtonStep << ", "
	   << hookcount++ << endl;
      cout << "delta == " << delta << endl;

      bool hookstep_equals_newtonstep;

      if (L2Norm(sN) <= delta) {
	cout << "Newton step is within trust region: " << endl;
	cout << "L2Norm(sN) == " << L2Norm(sN) << " <= " << delta << " == delta" << endl;
	hookstep_equals_newtonstep = true;
	sH = sN;
      }
      else {
	cout << "Newton step is outside trust region: " << endl;
	cout << "L2Norm(sN) == " << L2Norm(sN) << " > " << delta << " == delta" << endl;
	cout << "Calculate hookstep sH with radius L2Norm(sH) == delta" << endl;

	hookstep_equals_newtonstep = false;

	// This for-loop determines the hookstep. Search for value of mu that
	// produces ||shmu|| == delta. That provides optimal reduction of
	// quadratic model of residual in trust region norm(shmu) <= delta.
	Real mu = 0;
	Real nshmu = 0;
	ColumnVector shmu(Nk); // (s hat)(mu)]

	for (int hookSearch=0; hookSearch<Nhook; ++hookSearch) {
	  for (int i=0; i<Nk; ++i)
	    shmu(i) = bh(i)/(D(i,i) + mu); // research notes
	  nshmu = L2Norm(shmu);
	  Real Phi = nshmu*nshmu - delta*delta;
	  cout << "mu, L2Norm(sH(mu)) == " << mu << ", " << nshmu << endl;

	  // Found satisfactory value of mu and thus shmu and sH s.t. |sH| < delta
	  if (nshmu < delta || (nshmu > (1-deltaFuzz)*delta && nshmu < (1+deltaFuzz)*delta))
	    break;

	  // Update value of mu and try again. Update rule is a Newton search for
	  // Phi(mu)==0 based on a model of form a/(b+mu) for norm(sh(mu)). See
	  // Dennis & Schnabel.
	  else if (hookSearch<Nhook-1) {
	    Real PhiPrime = 0.0;
	    for (int i=0; i<Nk; ++i) {
	      Real di_mu = D(i,i) + mu;
	      Real bi = bh(i);
	      PhiPrime -= 2*bi*bi/(di_mu*di_mu*di_mu);
	    }
	    mu -= (nshmu/delta) * (Phi/PhiPrime);
	  }
	  else {
	    cout << "Couldn't find solution of hookstep optimization eqn Phi(mu)==0" << endl;
	    cout << "This shouldn't happen. It indicates an error in the algorithm." << endl;
	    cout << "Exiting.\n";
	    exit(1);
	  }
	} // search over mu for norm(s) == delta

	cout << "Found hookstep of proper radius" << endl;
	cout << "nshmu == " << L2Norm(shmu) << " ~= " << delta << " == delta " << endl;
	sH = Qn*(V*shmu);
      } // else {Newton step outside trust region, so

      // Compute
      // (1) actual residual from evaluation of
      //     rS == r(s+ds) act  == 1/2 (G(u+du,T+dT),G(u+du,T+dT))
      //
      // (2) predicted residual from quadratic model of r(s)
      //     rP == r(s+ds) pred  == 1/2 (G(u,T),G(u,T)) + (G(u,T), D[u]G du + D[t]G dT)
      //
      // (3) slope of r(s)
      //     dr/ds == (r(s + ds) - r(s))/|ds|
      //           == (G(u,T), D[u]G du + D[t]G dT)/|ds|
      // where ( , ) is FlowField inner product V2IP, which equals L2 norm of vector rep

      // (1) Compute actual residual of step, rS
      cout << "Computing residual of hookstep sH" << endl;
      vector2field(sH, du, basis);
      dunorm = L2Norm(du);
      uS  = u;
      uS += du;
      dT  = Tsearch ? sH(Tunk)*Tscale : 0;
      dax = Rsearch ? sH(xunk)*Rscale : 0.0;
      daz = Rsearch ? sH(zunk)*Rscale : 0.0;
      dtau = FieldSymmetry(1,1,1,dax,daz);

      cout << "L2Norm(du) == " << L2Norm(du) << endl;
      cout << "       dT  == " << dT << endl;
      cout << "     dtau  == " << dtau << endl;
      tS  = T + dT;
      tauS = dtau*tau;
      G(uS, tS, tauS, GS, Reynolds, flags, dt, Tnormalize, fcount_hookstep, CFL);
      cout << endl;

      Real rS = 0.5*Norm2(GS, nrm); // actual residual of hookstep
      Real Delta_rS = rS - ruT;          // improvement in residual
      //cout << "r(s), r(s+ds) == " << ruT << ", " << rS << endl;

      // (2) and (3) Compute quadratic model and slope
      cout << "Computing local quadratic model of residual" << endl;
      DG(u, du, T, dT, tau, dtau, Gu, DG_du, Reynolds, flags, dt, Tnormalize, epsDu,  centdiff, fcount_newton, CFL);
      cout << endl;
      // Local quadratic and linear models of residual, based on Taylor exp at current position
      // Quadratic model of r(s+ds): rQ == 1/2 |G(u,T) + D[u]G du + D[t]G dT)|^2
      // Linear  model   of r(s+ds): rL == 1/2 (G(u,T),G(u,T)) + (G(u,T), D[u]G du + D[t]G dT)

      Real Delta_rL = IP(Gu,DG_du,nrm); // rL - 1/2 (G(u,T),G(u,T)) == (G(u,T), D[u]G du + D[t]G dT)
      Real rL = ruT + Delta_rL;         // rL == 1/2 (G(u,T),G(u,T)) + (G(u,T), D[u]G du + D[t]G dT)

      if (rL >= ruT) {
	cout << "error : local linear model of residual is increasing, indicating\n"
	     << "        that the solution to the Newton equations is inaccurate\n"
	     << "exiting." << endl;
	exit(1);
      }

      DG_du += Gu;
      Real rQ = 0.5*Norm2(DG_du,nrm);   // rQ == 1/2 |G(u,T) + D[u]G du + D[t]G dT)|^2
      Real Delta_rQ = rQ - ruT;         // rQ == 1/2 (G(u,T),G(u,T)) - 1/2 |G(u,T) + D[u]G du + D[t]G dT)|^2

      snorm = L2Norm(sH);

      // ===========================================================================
      cout << "Determining what to do with current Newton/hookstep and trust region" << endl;

      // Coefficients for non-local quadratic model of residual, based on r(s), r'(s) and r(d+ds)
      // dr/ds == (G(u,T), D[u]G du + D[t]G dT)/|ds|
      Real drds  = Delta_rL/snorm;

      // Try to minimize a quadratic model of the residual
      // r(s+ds) = r(s) + r' |ds| + 1/2 r'' |ds|^2
      Real lambda = -0.5*drds*snorm/(rS - ruT - drds*snorm);

      // Compare the actual reduction in residual to the quadratic model of residual.
      // How well the model matches the actual reduction will determine, later on, how and
      // when the radius of the trust region should be changed.
      Real Delta_rS_req  = improvReq*drds*snorm;    // the minimum acceptable change in residual
      Real Delta_rS_ok   = improvOk*Delta_rQ;       // acceptable, but reduce trust region for next newton step
      Real Delta_rS_good = improvGood*Delta_rQ;     // acceptable, keep same trust region in next newton step
      Real Delta_rQ_acc  = abs(improvAcc*Delta_rQ); // for accurate models, increase trust region and recompute hookstep
      Real Delta_rS_accP = Delta_rQ + Delta_rQ_acc; //   upper bound for accurate predictions
      Real Delta_rS_accM = Delta_rQ - Delta_rQ_acc; //   lower bound for accurate predictions

      // Characterise change in residual Delta_rS
      ResidualImprovement improvement;  // fine-grained characterization

      // Place improvement in contiguous spectrum: Unacceptable > Poor > Ok > Good.
      if (Delta_rS > Delta_rS_req)
	improvement = Unacceptable;     // not even a tiny fraction of linear rediction
      else if (Delta_rS > Delta_rS_ok)
	improvement = Poor;             // worse than small fraction of quadratic prediction
      else if (Delta_rS < Delta_rS_ok && Delta_rS > Delta_rS_good)
	improvement = Ok;
      else {
	improvement = Good;             // not much worse or better than large fraction of prediction
	if (Delta_rS_accM <= Delta_rS &&  Delta_rS <= Delta_rS_accP)
	  improvement = Accurate;         // close to quadratic prediction
	else if (Delta_rS < Delta_rL)
	  improvement = NegaCurve;        // negative curvature in r(|s|) => try bigger step
      }

      cout << "ruT      == " << setw(w) << ruT           << " residual at current position" << endl;
      cout << " rS      == " << setw(w) << rS            << " residual of newton/hookstep" << endl;
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " actual improvement in residual from newton/hookstep" << endl;
      cout << endl;

      if (improvement==Unacceptable)
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " ------> Unacceptable <------ " << endl;
      cout << "            " << setw(w) << Delta_rS_req  << " lower bound for acceptable improvement" << endl;
      if (improvement==Poor)
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " ------> Poor <------" << endl;
      cout << "            " << setw(w) << Delta_rS_ok   << " upper bound for ok improvement." << endl;
      if (improvement==Ok)
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " ------> Ok <------" << endl;
      cout << "            " << setw(w) << Delta_rS_good << " upper bound for good improvement." << endl;
      if (improvement==Good)
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " ------> Good <------" << endl;
      cout << "            " << setw(w) << Delta_rS_accP << " upper bound for accurate prediction." << endl;
      if (improvement==Accurate)
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " ------> Accurate <------" << endl;
      cout << "            " << setw(w) << Delta_rS_accM << " lower bound for accurate prediction." << endl;
      cout << "            " << setw(w) << Delta_rL      << " local linear model of improvement" << endl;
      if (improvement==NegaCurve)
      cout << "Delta_rS == " << setw(w) << Delta_rS      << " ------> NegativeCurvature <------" << endl;

      bool recompute_hookstep = false;

      cout << "lambda       == " << lambda << " is the reduction/increase factor for delta suggested by quadratic model" << endl;
      cout << "lambda*delta == " << lambda*delta << " is the delta suggested by quadratic model" << endl;

      // See comments at end for outline of this control structure

      // Following Dennis and Schnable, if the Newton step lies within the trust region,
      // reset trust region radius to the length of the Newton step, and then adjust it
      // further according to the quality of the residual.
      //if (hookstep_equals_newtonstep)
      //delta = snorm;

      // CASE 1: UNACCEPTABLE IMPROVEMENT
      // Improvement is so bad (relative to gradient at current position) that we'll
      // in all likelihood get better results in a smaller trust region.
      if (improvement == Unacceptable) {
	cout << "Improvement is unacceptable." << endl;

	if (have_backup) {
	  cout << "But we have a backup step that was acceptable." << endl;
	  cout << "Revert to backup step and backup delta, take the step, and go to next Newton-GMRES iteration." << endl;
	  dtau = backup_dtau;
	  dT = backup_dT;
	  du = backup_du;
	  GS = backup_GS;
	  rS = backup_rS;
	  snorm = backup_snorm;
	  delta = backup_delta;

	  recompute_hookstep = false;
	}
	else {
	  cout << "No backup step is available." << endl;
	  cout << "Reduce trust region by minimizing local quadratic model and recompute hookstep." << endl;
	  deltaMaxLocal = delta;
	  lambda = adjustLambda(lambda, lambdaMin, lambdaRequiredReduction);
	  delta  = adjustDelta(delta, lambda, deltaMin, deltaMax);

	  if (delta > snorm) {
	    cout << "That delta is still bigger than the Newton step." << endl;
	    //cout << "This shouldn't happen. Control structure needs review." << endl;
	    cout << "Reducing delta to half the length of the Newton step." << endl;
	    cout << "  old delta == " << delta << endl;
	    delta = 0.5*snorm;
	    cout << "  new delta == " << delta << endl;
	  }
	  //delta_has_decreased = true;
	  recompute_hookstep = true;
	}
      }

      // CASE 2: EXCELLENT IMPROVEMENT AND ROOM TO GROW
      // Improvement == Accurate or Negacurve means we're likely to get a significant improvement
      // by increasing trust region. Increase trust region by a fixed factor (rather than quadratic
      // model) so that trust-region search is monotonic.
      // Exceptions are
      //   have_backup && backup_rS < rS -- residual is getting wrose as we increase delta
      //   hookstep==newtonstep          -- increasing trust region won't change answer
      //   delta>=deltamax               -- we're already at the largest allowable trust region
      else if ((improvement == NegaCurve || improvement == Accurate)
	       && !(have_backup && backup_rS < rS)
	       && !hookstep_equals_newtonstep
	       && !(delta >= deltaMax)) {

	cout << "Improvement is " << improvement << endl;

	//lambda = adjustLambda(lambda, lambdaMin, lambdaMax);
	Real new_delta = adjustDelta(delta, lambdaMax, deltaMin, deltaMax);

	if (new_delta < deltaMaxLocal) {
	  cout << "Continue adjusting trust region radius delta because" << endl;
	  cout << "Suggested delta has room: new_delta < deltaMaxLocal " << new_delta << " < " << deltaMaxLocal << endl;
	  if (have_backup)
	    cout << "And residual is improving:     rS < backup_sS     " << rS << " < " << backup_rS << endl;

	  cout << "Increase delta and recompute hookstep." << endl;
	  have_backup = true;
	  backup_dtau = dtau;
	  backup_dT = dT;
	  backup_du = du;
	  backup_GS = GS;
	  backup_rS = rS;
	  backup_snorm = snorm;
	  backup_delta = delta;

	  recompute_hookstep = true;
	  cout << " old delta == " << delta << endl;
	  delta = new_delta;
	  cout << " new delta == " << delta << endl;
	}
	else {
	  cout << "Stop adjusting trust region radius delta and take step because the new delta" << endl;
	  cout << "reached a local limit:  new_delta >= deltaMaxLocal " << new_delta << " >= " << deltaMaxLocal << endl;
	  cout << "Reset delta to local limit and go to next Newton iteration" << endl;
	  delta = deltaMaxLocal;
	  cout << "  delta == " << delta << endl;
	  recompute_hookstep = false;
	}
      }

      // CASE 3: MODERATE IMPROVEMENT, NO ROOM TO GROW, OR BACKUP IS BETTER
      // Remaining cases: Improvement is acceptable: either Poor, Ok, Good, or
      // {Accurate/NegaCurve and (backup is superior || Hookstep==NewtonStep || delta>=deltaMaxLocal)}.
      // In all these cases take the current step or backup if better.
      // Adjust delta according to accuracy of quadratic prediction of residual (Poor, Ok, etc).
      else {
	cout << "Improvement is " << improvement <<  " (some form of acceptable)." << endl;
	cout << "Stop adjusting trust region and take a step because" << endl;
	cout << "Improvement is merely Poor, Ok, or Good : " << (improvement == Poor || improvement == Ok || improvement == Good) << endl;
	cout << "Newton step is within trust region      : " << hookstep_equals_newtonstep << endl;
	cout << "Backup step is better                   : " << (have_backup && backup_rS < rS) << endl;
	cout << "Delta has reached local limit           : " << (delta >= deltaMax) << endl;

	recompute_hookstep = false;

	// Backup step is better. Take it instead of current step.
	if (have_backup && backup_rS < rS) {
	  cout << "Take backup step and set delta to backup value." << endl;
	  dtau = backup_dtau;
	  dT = backup_dT;
	  du = backup_du;
	  GS = backup_GS;
	  rS = backup_rS;

	  cout << "  old delta == " << delta << endl;
	  snorm = backup_snorm;
	  delta = backup_delta;
	  cout << "  new delta == " << delta << endl;
	}

	// Current step is best available. Take it. Adjust delta according to Poor, Ok, etc.
	// Will start new monotonic trust-region search in Newton-hookstep, so be more
	// flexible about adjusting delta and use quadratic model.
	else {
	  if (have_backup) {
	    cout << "Take current step and keep current delta, since it's produced best result." << endl;
	    cout << "delta == " << delta << endl;
	    // Keep current delta because current step was arrived at through a sequence of
	    // adjustments in delta, and this value yielded best results.
	  }
	  else if (improvement == Poor) {
	    cout << "Take current step and reduce delta, because improvement is poor." << endl;
	    lambda = adjustLambda(lambda, lambdaMin, 1.0);
	    delta  = adjustDelta(delta, lambda, deltaMin, deltaMax);
	  }
	  else if (improvement == Ok || hookstep_equals_newtonstep || (delta >= deltaMax)) {
	    cout << "Take current step and leave delta unchanged, for the following reasons:" << endl;
	    cout << "improvement            : " << improvement << endl;
	    cout << "|newtonstep| <= delta  : " << (hookstep_equals_newtonstep ? "true" : "false") << endl;
	    cout << "delta >= deltaMax      : " << (delta >= deltaMax    ? "true" : "false") << endl;
	    //cout << "delta has decreased    : " << (delta_has_decreased ? "true" : "false") << endl;
	    cout << "delta == " << delta << endl;
	  }
	  else {  // improvement == Good, Accurate, or NegaCurve and no restriction on increasing apply
	    cout << "Take step and increase delta, if there's room." << endl;
	    lambda = adjustLambda(lambda, 1.0, lambdaMax);
	    delta  = adjustDelta(delta, lambda, deltaMin, deltaMax);
	  }
	}
      }
      //cout << "Recompute hookstep ? " << (recompute_hookstep ? "yes" : "no") << endl;

      if (pausing)
	cfpause();

      if (recompute_hookstep)
	continue;
      else
	break;
    }

    cout << "Taking best step and continuing Newton iteration." <<endl;
    u += du;
    T += dT;
    tau *= dtau;
    dunorm = L2Norm(du);
    unorm  = L2Norm(u);

    Gu = GS;

    cout << "L2Norm(G(u)) == " << L2Norm(Gu) << endl;
    cout << "ruT          == " << ruT << endl;
    cout << "L2Norm(sN)   == " << L2Norm(sN) << endl;
    cout << "L2Norm(sH)   == " << L2Norm(sH) << endl;
    cout << "L2Norm(du)   == " << dunorm << endl;
    cout << "       dT    == " << dT << endl;
    cout << "L2IP(du,edudx)/L2Norm(du) == " << L2IP(du,edudx)/dunorm << endl;
    cout << "L2IP(du,edudz)/L2Norm(du) == " << L2IP(du,edudz)/dunorm << endl;
    cout << "L2IP(du,edudt)/L2Norm(du) == " << L2IP(du,edudt)/dunorm << endl;
    //du.binarySave("duS" + i2s(newtonStep));
  }
  fftw_cleanup();
}

// f^T(u) = time-T DNS integration of u
void f(const FlowField& u, Real T, FlowField& f_u, Real Reynolds,
       DNSFlags& flags, TimeStep& dt, int& fcount, Real& CFL) {

  if (T<0) {
    cout << "f: negative integration time T == " << T << endl;
    cout << "returning f(u,T) == (1+abs(T))*u" << endl;
    f_u = u;
    f_u *= 1+abs(T);
    return;
  }
  if (!isfinite(L2Norm(u))) {
    cout << "error in f: u is not finite. exiting." << endl;
    exit(1);
  }

  // Special case #1: no time steps
  if (T == 0) {
    cout << "f: T==0, no integration, returning" << endl;
    f_u = u;
    return;
  }

  const Real nu = 1.0/Reynolds;
  FlowField p(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b());
  f_u = u;
  dt.adjust_for_T(T, false);

  // Adjust dt for CFL if necessary
  DNS dns(f_u, nu, dt, flags);
  if (dt.variable()) {
    dns.advance(f_u, p, 1);
    dt.adjust(dns.CFL(), false);
    dns.reset_dt(dt);
    f_u = u;
  }

  //cout << "\ndt == " << dt << endl;
  //cout << ", CFL == " << dns.CFL() << endl;
  //cout << "t == " << flush;

  //  t == current time in integration
  //  T == total integration time
  // dT == CFL-check interval
  // dt == DNS time-step
  //  N == T/dT,  n == dT/dt;
  //  T == N dT, dT == n dt
  //  t == s dT (s is loop index)
  cout << "f : 0" << flush;
  for (int s=1; s<=dt.N(); ++s) {
    dns.advance(f_u, p, dt.n());

    Real t = s*dt.dT();
    if (s % 10 == 0)
      cout << iround(t) << flush;
    else if (s % 2 == 0)
      cout << '.' << flush;

    CFL = dns.CFL();
    if (CFL < dt.CFLmin())
      cout << '<' << flush;
    else if (CFL > dt.CFLmax())
      cout << '>' << flush;
    if (dt.adjust(CFL, false))
      dns.reset_dt(dt);
  }

  ++fcount;

  if (!isfinite(L2Norm(f_u))) {
    cout << "error in f: f(u,t) is not finite. exiting." << endl;
    exit(1);
  }
}

/// f^{N dt}(u) = time-(N dt) DNS integration of u
void f(const FlowField& u, int N, Real dt, FlowField& f_u, Real Reynolds,
        DNSFlags& flags) {

  cout << "f(u, N, dt, f_u, Reynolds, flags, dt) : " << flush;
  const Real nu = 1.0/Reynolds;
  FlowField p;
  f_u = u;
  DNS dns(f_u, nu, dt, flags);
  dns.advance(f_u, p, N);
  return;
}

// G(u) = (tau f^T(u) - u) for orbits and (tau f^T(u) - u)/T for eqbs
void G(const FlowField& u, Real T, const FieldSymmetry& tau, FlowField& G_u,
       Real Reynolds, DNSFlags& flags, TimeStep& dt, bool Tnormalize, int& fcount, Real& CFL) {
  f(u, T, G_u, Reynolds, flags, dt, fcount, CFL);
  G_u *= tau;
  G_u -= u;
  if (Tnormalize)
    G_u *= 1.0/T;
}

void DG(const FlowField& u, const FlowField& du, Real T, Real dT,
	const FieldSymmetry& tau,  const FieldSymmetry& dtauArg,
	const FlowField& GuT, FlowField& DG_du,
	Real Reynolds, DNSFlags& flags, TimeStep& dt, bool Tnormalize,
	Real epsDu, bool centdiff, int& fcount, Real& CFL) {

  FlowField u_du;
  FieldSymmetry dtau;
  FieldSymmetry tau_dtau;
  Real T_dT;

  Real step_magn = sqrt(L2Norm2(du) + square(dT/T) + square(dtau.ax()) + square(dtau.az()));
  Real eps = step_magn < epsDu ? 1 : epsDu/step_magn;

  if (centdiff) {
    Real eps2 = 0.5*eps;
    u_du  = du;
    u_du *= eps2;
    u_du += u;
    T_dT = T + eps2*dT;
    dtau = dtauArg;
    dtau *= eps2;
    tau_dtau = dtau*tau;

    G(u_du, T_dT, tau_dtau, DG_du, Reynolds, flags, dt, Tnormalize, fcount, CFL);
    u_du  = du;
    u_du *= -eps2;
    u_du += u;
    T_dT = T - eps2*dT;
    dtau = dtauArg;
    dtau *= -eps2;
    tau_dtau = dtau*tau;
    FlowField tmp;
    G(u_du, T_dT, tau_dtau, tmp, Reynolds, flags, dt, Tnormalize, fcount, CFL);

    DG_du -= tmp;
    DG_du *= 1.0/eps;
  }
  else {
    u_du  = du;
    u_du *= eps;
    u_du += u;
    T_dT = T + eps*dT;
    dtau = dtauArg;
    dtau *= eps;
    tau_dtau = dtau*tau;

    G(u_du, T_dT, tau_dtau, DG_du, Reynolds, flags, dt, Tnormalize, fcount, CFL);
    DG_du -= GuT;
    DG_du *= 1.0/eps;
  }
  //if (Tnormalize)
  //DG_du *= 1.0/T;
}

Real V2Norm2(const FlowField& u) {
  ColumnVector v;
  field2vector(u,v);
  return L2Norm2(v);
}

Real V2Norm(const FlowField& u) {
  return sqrt(V2Norm2(u));
}

Real V2IP(const FlowField& u, const FlowField& v) {
  ColumnVector uV;
  ColumnVector vV;
  field2vector(u,uV);
  field2vector(v,vV);
  return (uV.transpose())*vV;
}

Real Norm2(const FlowField& u, BasisType norm) {
  if (norm == L2Basis)
    return L2Norm2(u);
  else
    return V2Norm2(u);
}

Real Norm(const FlowField& u, BasisType norm) {
  if (norm == L2Basis)
    return L2Norm(u);
  else
    return V2Norm(u);
}

Real IP(const FlowField& u, const FlowField& v, BasisType norm) {
  if (norm == L2Basis)
    return L2IP(u,v);
  else
    return V2IP(u,v);
}

/*********************************
ostream& operator<<(ostream& os, HookstepPhase p) {
  if      (p == ReducingDelta)   os << "ReducingDelta";
  else if (p == IncreasingDelta) os << "IncreasingDelta";
  else                           os << "Finished";
  return os;
}
ostream& operator<<(ostream& os, ResidualImprovement i) {
  if      (i == Unacceptable) os << "Unacceptable";
  else if (i == Poor)         os << "Poor";
  else if (i == Ok)           os << "Ok";
  else if (i == Good)         os << "Good";
  else if (i == Accurate)     os << "Accurate";
  else                        os << "NegativeCurvature";
  return os;
 }

Real adjustDelta(Real delta, Real deltaRate, Real deltaMin, Real deltaMax) {
  cout << "  old delta == " << delta << endl;

  // Check to see if we're already running with delta == deltaMin.
  // If so, and if deltaRate calls for furhter reduction, search has
  // ground to a halt. Error and exit.
  if (delta <= 1.00001*deltaMin && deltaRate < 1.0) {
    cout << "    delta == " << delta << endl;
    cout << "delta new == " << delta*deltaRate << endl;
    cout << "delta min == " << deltaMin << endl;
    cout << "Sorry, can't go below deltaMin. Exiting." << endl;
    exit(1);
  }

  delta *= deltaRate;

  if (delta < deltaMin) {
    delta = deltaMin;
    cout << "delta bottomed out at deltaMin" << endl;
  }
  else if (delta > deltaMax) {
    delta = deltaMax;
    cout << "delta topped out at deltaMax" << endl;
  }
  cout << "  new delta == " << delta << endl;
  return delta;
}

Real adjustLambda(Real lambda, Real lambdaMin, Real lambdaMax) {
  if (lambda < lambdaMin) {
    cout << "lambda == " << lambda << " is too small. Resetting to " << endl;
    lambda = lambdaMin;
    cout << "lambda == " << lambda << endl;
  }
  else if (lambda < lambdaMin) {
    cout << "lambda == " << lambda << " is too large. Resetting to " << endl;
    lambda = lambdaMax;
    cout << "lambda == " << lambda << endl;
  }
  return lambda;
}
*************************************/

// The following is meant as a guide to the structure and notation
// of the above code. See Divakar Viswanath, "Recurrent motions
// within plane Couette turbulence", J. Fluid Mech. 580 (2007),
// http://arxiv.org/abs/physics/0604062 for a higher-level presentation
// of the algorithm.
//
// I've biased the notation towards conventions from fluids and from
// Predrag Cvitanivic's ChaosBook, and away from Dennis & Schnabel and
// Trefethen. That is, f stands for the finite-time map of plane Couette
// dynamics, not the residual being minimized.
//
// Note: in code, dt stands for the integration timestep and
//                dT stands for Newton-step increment of period T
// But for beauty's sake I'm using t and dt for the period and
// its Newton-step increment in the follwoing comments.
//
// Let
//  f^t(u) = time-t map of Navier-Stokes computed by DNS.
//  G(u,t) = (sigma f^t(u) - u) = 0 is the eqn to be solved
//  r(u,t) = 1/2 L2Norm2(G(u,t)) residual, to be minimized
//
// We want to find zeros of G(u,t) or equiv minimum of r(u,t).
// Note that the residual is reported to the user in the form
//   residual(u,t) == L2Norm(G(u,T)) == sqrt(2*r(u,t))
//
// The Newton algorithm for solving G(u,t) = 0 is
//
// Let u* be the solution, i.e. G(u*,t*) = 0
// Suppose we start with an initial guess u near u*.
// Let u* = u + du and t* = t + dt.
// Then G(u*,t*) = G(u+du, t+dt)
//             0 = G(u,t) + D[u]G(u,t) du + D[t]G(u,t) dt + O(|du|^2, dt^2)
// Drop h.o.t and solve
//   D[u]G(u,t) du + D[t]G(u,t) = -G(u)
//                  (du, du/dt) = 0
// Then start with new guess u' = u+du t'=t+dt and iterate.
//
// Use GMRES to solve the Newton-step equation. GMRES requires multiple
// calculations of the LHS DG(u,t) du + dG(u,t)/dt dt for trial values
// du, dt. These can be approximated with finite differencing:
//   D[u]G(u,t) du = 1/e (G(u + e du, t) - G(u,t)) where e |du| << 1.
//   D[t]G(u,t) dt = 1/e (G(u, t + e dt) - G(u,t)) where e |dt| << 1.
//
// To use GMRES, we need to cast velocity fields such as du and DG(u) du
// into vectors and vice versa. Do this with a real-valued, orthonormal,
// div-free, dirichlet basis {Phi_n, n=0,...,N-1}
//   a_n = (u, Phi_n), n=0,..,N-1
//     u = a_n Phi_n,  n=0,..,N-1
//     t = a_N

// Then the Newton step equations
//  D[u]G(u,t) du + D[t]G(u,t) dt = -G(u)
//                  (du, edudt) = 0
// can be translated to a linear algebra problem
//      A s = b
// via
//      s_n = (du, Phi_n)    n=0,..,N-1
//      b_n = (-G(u,t), Phi_n) n=0,..,N-1
//  (A s)_n = (D[u]G(u,t) du + D[t]G(u,t) dt, Phi_n)
//      s_N = dt
//      b_N = 0
//  (A s)_N = (du, du/dt)

// We then solve this system approximately with GMRES.
//
// Hookstep algorithm: If the Newton step du = s_n Phi_n does not decrease
// the residual L2Norm(G(u+du, u+dt)) sufficiently, we obtain the a smaller step
// s by minimizing || A s - b ||^2 subject to ||s||^2 <= delta^2, rather
// solving A s = b directly.

// But since A is M x M, with M = O(10^5), we minimize the norm of the
// projection of A s - b onto the Krylov subspace from GMRES. GMRES gives
//   A Qn = Qn1 Hn, where
// where Qn is M x n, Qn1 is M x (n+1), and Hn is (n+1) x n.
// Projection of A s - b on (n+1)th Krylov subspace is Qn1^T (A s - b)
// Confine s to nth Krylov subspace: s = Qn sn, where sn is n x 1.
// Then minimize
//   || Qn1^T (A s - b) || = || Qn1^T A Qn sn - Qn1^T b ||
//                         = || Qn1^T Qn1 Hn sn - Qn1^T b ||
//                         = || Hn sn - Qn1^T b ||
// subject to ||s||^2 = ||sn||^2 < delta^2
//
// Do minimization via SVD of Hn = U D V^T.
//   || Qn1^T (A s - b) || = || Hn sn - Qn1^T b ||
//                         = || U D V^T sn - Qn1^T b ||
//                         = || D V^T sn - U^T Qn1^T b ||
//                         = || D sh - bh ||
// where sh = V^T sn and bh = U^T Qn1^T b.
//
// From Divakar Viswanath (personal communication):
//
// Since D is diagonal D_{i,i} = d_i, the constrained minimization
// problem can be solved easily with lagrange multipliers:
//   sh_i = (bh_i d_i)/(d^2_i + mu)
// for the mu such that ||sh||^2 = delta^2. The value of mu can be
// found with a 1d Newton search for zeros of
//   Phi(mu) = ||sh(mu)||^2 - delta^2
// A straight Newton search a la mu -= Phi(mu)/Phi'(mu) is suboptimal
// since Phi'(mu) -> 0 as mu -> infty with Phi''(mu) > 0 everywhere, so
// we use a slightly modified update rule for the Newton search. Please
// refer to Dennis and Schnabel regarding that. Then, given the solution
// sh(mu), we compute the hookstep solution s from
//   s = Qn sn = Qn V sh

// D[t]G(u,t) dt == D[t] (tau f^(t+dt)(u) - u) - (tau f^(t)(u) - u)
//               == tau (f^(t+dt)(u) - f^(t)(u))
//               == tau (f^dt(f^t(u)) - f^t(u))
//               == tau (f^dt(f^t(u)) - f^t(u))
//               == tau (df^T(u)/dT) * dT