#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 "utilities.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 f^T(u) - u = 0 where f^T(u) is the time-T
// CFD integration of the fluid velocity field u. Viswanath's algorithm 
// combines a Newton-hookstep trust-region minimization of ||f^T(u) - u||^2 = 0.
// with iterative GMRES solution of the Newton-step equation 
// Df^T/du du + Df^T/dT dT = -f^T(u).

// 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.</em> 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, 
       TimeStepParams& dtp, 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) 
void G(const FlowField& u, Real T, const FieldSymmetry& tau, FlowField& G_u, 
       Real Reynolds, DNSFlags& flags, TimeStepParams& dtp, int& fcount, Real& CFL);

void DG(const FlowField& u, const FlowField& du, Real T, Real dT, const FieldSymmetry& tau, 
	const FlowField& GuT, FlowField& DG_du, Real Reynolds, DNSFlags& flags, 
	TimeStepParams& dtp, Real eps, 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, 
	TimeStepParams& dtp, Real eps, int& fcount, Real& CFL);

void params2timestep(const TimeStepParams& dtp, Real T, int& N, TimeStep& dt);

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);


void putBetween(Real& delta, Real deltaMin, Real deltaMax);

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 Real epsSearch = args.getreal("-es",  "--epsSearch",   1e-14, "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 int  Nnewton   = args.getint( "-Nn", "--Nnewton",   10,    "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 bool fixedDelta= args.getflag("-fd",   "--fixedDelta",        "don't adjust trust region radius");
  const Real deltaMin  = args.getreal("-dmin", "--deltaMin",   1e-6,  "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 deltaRate = args.getreal("-dr",   "--deltaRate",  2.0,   "delta *= deltaRate if hookstep is accurate");
  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");

  /***/ Real T         = args.getreal("-T",   "--maptime",    10.0,  "(initial) integration time for DNS map f^T(u)");
  const Real Tsearch   = args.getbool("-Ts",  "--Tsearch",    true,  "let integration time vary T in the search");
  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 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 Rsearch   = args.getbool("-r",  "--relative", false,  "search for relative equilib or periodic orbit");
  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 dt0       = args.getreal("-dt","--timestep",  0.02, "(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 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");

  const BasisType nrm = l2basis ? L2Basis : AdHocBasis;
  const Real deltaShrinkMin = 0.1;
  const Real deltaShrinkMax = 0.5;
  const int w = coutprec + 7; 

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

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

  int fcount_newton = 0;
  int fcount_hookstep = 0;

  // Used to communicate timestep parameters to integration functions.
  cout << "dt0    == " << dt0 << endl;
  cout << "dtmin  == " << dtmin << endl;
  cout << "dtmax  == " << dtmax << endl;
  cout << "CFLmin == " << CFLmin << endl;
  cout << "CFLmax == " << CFLmax << endl;
  TimeStepParams dtparams(dt0, dtmin, dtmax, CFLmin, CFLmax, vardt);

  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;
  cout << "Teqn  == " << Teqn << endl;
  cout << "Tunk  == " << Tunk << endl;

  sN = ColumnVector(Nunk);

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

  FlowField      du(Nx, Ny, Nz, 3, Lx, Lz, ya, yb); // Newton step for 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 velcoity 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);
  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"
	 << setw(14) << "dT"
	 << setw(14) << "dax"
	 << setw(14) << "daz"
	 << 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) = tau f(u,T) - u" << endl; 
  G(u, T, tau, Gu, Reynolds, flags, dtparams, 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
    // 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 T*Tscale ~= TscaleRel*(size of coeffs in u);
    // 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;
    osconv.setf(ios::left); osconv << setw(14) << dT;
    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));
    if (Tsearch) 
      save(T, outdir + "Tnewt" + i2s(newtonStep));
    if (Rsearch) 
      tau.save(outdir + "taunewt" + i2s(newtonStep));

    //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("usoln");
      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 << ',' << setw(2) << 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, dtparams, epsDu,  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);

    // ==================================================================
    // 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 Newton/hookstep" << 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, dtparams, 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, dtparams, epsDu,  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);

      // 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;
      //Real drds2_global = 2*(rS - ruT - drds*snorm)/square(snorm);
      //Real drds2_local  = 2*(rQ - ruT - drds*snorm)/square(snorm);
      //cout << "drds2_global == " << drds2_global << endl;
      //cout << "drds2_local  == " << drds2_local << endl;

      // ===========================================================================
      cout << "Determining what to do with current Newton/hookstep and trust region" << endl;
      
      // 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

      // Special cases, which are bands within above spectrum
      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;

      // 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;

      // 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 << "Unacceptable improvement from Newton-hook step." << endl; 
	
	if (have_backup) {
	  cout << "But we have a backup step that was acceptable." << endl;
	  cout << "Revert to backup step, 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 and recompute hookstep." << endl;

	  // Try to minimize a quadratic model of the residual
	  // r(s+ds) = r(s) + r' |ds| + 1/2 r'' |ds|^2
	  // But require delta
	  Real lambda = -0.5*drds*snorm/(rS - ruT - drds*snorm);
	  lambda = Greater(deltaShrinkMin, lambda); // greatest lower bound
	  lambda =  lesser(deltaShrinkMax, lambda); // least upper bound
	  cout << "  old delta == " << delta << endl;
	  delta *= lambda;
	  cout << "  new delta == " << delta << endl;
	  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 << "Setting delta to half the length of the Newton step." << endl;
	    cout << "  old delta == " << delta << endl;
	    delta = 0.5*snorm;
	    cout << "  new delta == " << delta << endl;
	  }
	  if (delta < deltaMin) {
	    cout << "Reached minimum trust radius delta. Search has ground to a halt. Exiting." << endl;
	    exit(1);
	  }
	  delta_has_decreased = true;
	  deltaMaxLocal = delta;
	  recompute_hookstep = true;
	}
      }

      // Improvement is so good that we're likely to get significantly better result
      // by increasing trust region. Exceptions are 
      //   rS < rQ              -- the actual improvement is better quadratic prediction.
      //   delta_has_decreased  -- we got here from a larger, unacceptable trust region
      //   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)
	       && !delta_has_decreased
	       && !hookstep_equals_newtonstep
	       && delta < deltaMax
	       ) {

	if (improvement == Accurate)
	  cout << "Improvement is Accurate: close to prediction of quadratic model." << endl;
	else 
	  cout << "Improvement is NegaCurve: better than predicted by linear model, suggesting negative curvature." << endl;

	cout << "And we still have room to expand trust region." << 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;

	Real new_delta = adjustDelta(delta, deltaRate, deltaMin, deltaMax);
	if (new_delta < deltaMaxLocal) {
	  recompute_hookstep = true;
	  delta = new_delta;
	}
	else 
	  recompute_hookstep = false;
      }
	
      // Remaining cases: Improvement is acceptable: either Poor, Ok, Good, or 
      // {Accurate/NegaCurve and (delta_has_decreased || Hookstep==NewtonStep || delta>=deltaMax)}. 
      // 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 << endl;

	recompute_hookstep = false;

	// Backup step is better. Take it instead of current step.
	if (have_backup && backup_rS < rS) {
	  cout << "But the backup step is better." << endl;
	  cout << "Take the 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.
	else {

	  if (have_backup) {
	    cout << "Current step is better than backup step." << endl;
	    cout << "Take current step and keep current delta." << endl;
	    // Keep current delta because current step was best guess at a better delta,
	    // and it turned out to have a worse improvement in residual
	  }

	  // In all the following cases, we have no backup information, so delta adjustment
	  // is based entirely on quality of current residual. 
	  else if (improvement == Poor) {
	    cout << "Reduce delta and take step." << endl;
	    //cout << "  old delta == " << delta << endl;
	    //delta /= deltaRate;
	    //cout << "  new delta == " << delta << endl;
	    //putBetween(delta, deltaMin, deltaMax);
	    delta = adjustDelta(delta, 1/deltaRate, deltaMin, deltaMax);

	  }
	  else if (improvement == Ok) {
	    cout << "Leave delta unchanged and take step." << endl;
	    cout << "  delta == " << delta << endl;
	    //putBetween(delta, deltaMin, deltaMax);
	  }
	  // improvement == Good, Accurate, or NegaCurve
	  else {
	    if (improvement == Accurate || improvement == NegaCurve) {
	      cout << "Don't recomputing hookstep because ... " << endl;
	      cout << "  delta has decreased : " << (delta_has_decreased ? "true" : "false") << endl;
	      cout << "  hook equals newton  : " << (hookstep_equals_newtonstep ? "true" : "false") << endl;
	      cout << "  delta >= deltaMax   : " << (delta >= deltaMax ? "true" : "false") << endl;
	      cout << "And leave delta unchanged. " << endl;
	      cout << "  delta == " << delta << endl;
	    }
	    else {
	      cout << "Increase delta for next hookstep, if there's room." << endl;
	      //cout << "  old delta == " << delta << endl;
	      //delta *= deltaRate;
	      //cout << "  new delta == " << delta << endl;
	      //putBetween(delta, deltaMin, deltaMax);
	      delta = adjustDelta(delta, deltaRate, 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(-hookstep) iteration." <<endl;
    u += du;
    T += dT;
    tau *= dtau;
    dunorm = L2Norm(du);
    unorm  = L2Norm(u);

    Gu = GS;
    
    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;
  }
  fftw_cleanup();
}

void params2timestep(const TimeStepParams& p, Real T, int& N, TimeStep& dt) {

  // Set CFL adjustment time dT to an integer division of T near 1
  N = Greater(iround(T), 1);
  const Real dT = T/N;
  if (p.variable)
    dt = TimeStep(p.dt, p.dtmin, p.dtmax, dT, p.CFLmin, p.CFLmax);
  else 
    dt = TimeStep(p.dt, p.dt, p.dt, dT, p.CFLmin, p.CFLmax);
}


// f^T(u) = time-T DNS integration of u
void f(const FlowField& u, Real T, FlowField& f_u, Real Reynolds, DNSFlags& flags, 
       TimeStepParams& dtp, 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;
  }

  //cout << "f(u,T) : L2Norm(u) == " << L2Norm(u) << ",  T == " << T << flush;
  

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

  int N;
  TimeStep dt;
  params2timestep(dtp, T, N, dt);
  //cout << "TimeStep dt == " << dt << endl;

  f_u = u;

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

  //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<=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 (dtp.variable && dt.adjust(CFL))
      dns.reset_dt(dt);
  }

  ++fcount;
  dtp.dt = Real(dt);

  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(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b());
  f_u = u;
  DNS dns(f_u, nu, dt, flags);
  dns.advance(f_u, p, N);
  return;
}


void G(const FlowField& u, Real T, const FieldSymmetry& tau, FlowField& G_u,  
       Real Reynolds, DNSFlags& flags, TimeStepParams& dtp, int& fcount, Real& CFL) {
  f(u, T, G_u, Reynolds, flags, dtp, fcount, CFL);
  G_u *= tau;
  G_u -= u;
}

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, TimeStepParams& dtp, 
	Real epsDu, int& fcount, Real& CFL) {

  FlowField u_du;
  FlowField tmp;
  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;

  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, dtp, fcount, CFL);
  DG_du -= GuT;
  DG_du *= 1.0/eps;
}

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;
 }

void savewave(FlowField& u, const string& name) {

  int Nx = u.Nx();
  int Ny2 = (u.Ny()-1)/2;
  
  u.makePhysical();
  Vector f(Nx+1);
  for (int nx=0; nx<Nx; ++nx)
    f(nx) = u(nx,Ny2,0,2);
  f(Nx) = u(0,Ny2,0,2);
  f.save(name);
  u.makeSpectral();
}
  
void putBetween(Real& delta, Real deltaMin, Real deltaMax) {
  assert(deltaMin <= deltaMax);
  if (delta < deltaMin) {
    delta = deltaMin;
    cout << "delta bottomed out at deltaMin" << endl;
    cout << "  new delta == " << delta << endl;
  }
  else if (delta > deltaMax) {
    delta = deltaMax;
    cout << "delta topped out at deltaMax" << endl;
    cout << "  new delta == " << delta << endl; 
  }
}

Real adjustDelta(Real delta, Real deltaRate, Real deltaMin, Real deltaMax) {
  cout << "  old delta == " << delta << endl;
  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;
}

// 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.</em> 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) = (tau f^t(u) - u) for periodic orbit searches
//  r(u,t) = 1/2 L2Norm2(G(u,t)), 
//   
// 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) dt = -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