File : rg_ops.ads


with Ada.Text_IO, Reps, Fouriers_Init, Lin;
use Ada.Text_IO, Reps, Fouriers_Init, Lin;
with Fouriers, Fouriers.Ops, Fouriers.IO;

pragma Elaborate_All (Ada.Text_IO, Fouriers.Ops, Fouriers.IO);

generic

  Lmax, Wmax, Pdeg: in Positive;
  Dho: in Natural;
  Sigma, Kappa: Rep;
  type Scalar is private;

  with function IsNumeric(S: Scalar) return Boolean is <>;
  with function Scal(K: Integer) return Scalar is <>;
  with function Scal(R: Rep) return Scalar is <>;
  with function Ball(R: Rep) return Scalar is <>;
  with function Ball(S: Scalar) return Scalar is <>;
  with function Err0(S: Scalar) return Boolean is <>;
  with function "<"(S1,S2: Scalar) return Boolean is <>;
  with function IntFloor(S: Scalar) return Integer is <>;
  with function IntCeiling(S: Scalar) return Integer is <>;
  with function Approx(S: Scalar) return Rep is <>;
  with function Inf(S: Scalar) return Rep is <>;
  with function Sup(S: Scalar) return Rep is <>;
  with function SupAbs(S: Scalar) return Rep is <>;
  with function "abs"(S: Scalar) return Scalar is <>;
  with function Min(S1,S2: Scalar) return Scalar is <>;
  with function Max(S1,S2: Scalar) return Scalar is <>;
  with function Cap(R: Radius; S: Scalar) return Scalar is <>;
  with function "-"(S: Scalar) return Scalar is <>;

  with function Up(R: Rep; Dummy: Scalar) return Rep is <>;
  with function Center(S: Scalar) return Scalar is <>;
  with function Widen(S: Scalar) return Scalar is <>;
  with procedure ResetCenter(S: in out Scalar) is <>;
  with function ResetCenter(S: Scalar) return Scalar is <>;
  with procedure Split(S: in Scalar; S0,SE: out Scalar) is <>;
  with procedure ErrMult(R: in Rep; S: in out Scalar) is <>;
  with function "+"(S,T: Scalar) return Scalar is <>;
  with function "-"(S,T: Scalar) return Scalar is <>;
  with function "*"(R: Rep; S: Scalar) return Scalar is <>;
  with function "/"(S: Scalar; R: Rep) return Scalar is <>;
  with function Sqr(S: Scalar) return Scalar is <>;
  with function "*"(S,T: Scalar) return Scalar is <>;
  with function Inv(S: Scalar) return Scalar is <>;
  with function "/"(S,T: Scalar) return Scalar is <>;
  with function "**"(S: Scalar; I: Integer) return Scalar is <>;

  with function Cosh(S: Scalar) return Scalar is <>;
  with function Sinh(S: Scalar) return Scalar is <>;
  with function Exp(S: Scalar) return Scalar is <>;
  with function Log(S: Scalar) return Scalar is <>;
  with function Lambda(M,N: Integer) return Scalar is <>;

  with procedure Show(S: in Scalar; NewLine: in Boolean := True) is <>;
  with procedure Get(F: in File_Type; S: out Scalar; Decimal: in Boolean := True) is <>;
  with procedure Put(F: in File_Type; S: in Scalar; Decimal: in Boolean := True) is <>;

package RG_Ops is

  package Scal_Fou is new Fouriers (Lmax => Lmax, Wmax => Wmax, Pdeg => Pdeg, Dho => Dho, Scalar => Scalar);
  package Scal_Fou_Ops is new Scal_Fou.Ops;
  package Scal_Fou_IO is new Scal_Fou.IO;

  subtype Hamilton is Scal_Fou.Fourier;

  type HBall is
    record
      R: Weights;
      E: Scalar;
    end record;

  Cpp:      constant Rep    := Half;
  ScalZero: constant Scalar := Scal(Zero);
  ScalOne:  constant Scalar := Scal(One);
  Theta:    constant Scalar := Lambda(1,-1);
  Dminus:   constant Degrees(Scal_Fou.Freq1,Scal_Fou.Freq2) := Scal_Fou.DminusDef(Sigma,Kappa);

  procedure ShowParam;                                                    -- Does what is sais
  function UP(R: Rep) return Rep;                                         -- Used only to increase domains

  procedure Inv(Err: in Radius; R: in Weights;
                H1: in Hamilton; H2: out Hamilton);                       -- H2 := 1/H1
  function Inv(Err: Radius; R: Weights; H: Hamilton) return Hamilton;     -- Return 1/H
  procedure FixComp3(Err: in Radius; R3: in Radius; H: in Hamilton;
                     Hz: out Hamilton);                                   -- Hz Solves -H(.,.,.+Hz)=Hz
  procedure HoToU(Psi,Hq,Hzp,H: in Hamilton;
                  Hu: out Hamilton; Full: in Boolean := True);            -- Hu := (Theta/Mu0)H o T_mu0 o U_Psi

  procedure Prod(A: in Scalar; W: in out Sparse;
                 H1: in Hamilton; H2: out Hamilton);                      -- H2 := (A+W)*H1

  procedure Normalize(H1: in Hamilton; H2: out Hamilton);                 -- H2 := (H1 normalized)
  procedure NormalizeM(W: in out Sparse; H: in Hamilton;
                       F: in Hball; H2: out Hamilton);                    -- Normalize H+(I-W)*F
  procedure DNormalize(H,F: in Hamilton; F2: out Hamilton);               -- Derivative of Normalize at H_normalized
  procedure UpsiAux(Err: in Radius; R: in Weights;
                        Psi: in Hamilton; Hq,Hzp: out Hamilton);          -- Compute Hamiltonians used for Upsi

  function NormPrime(R: Weights; H: Hamilton) return Scalar;              -- Return ||H||'_R
  procedure Nash(Err: in Radius; R: in Weights; H: in out Hamilton);      -- Eliminate Iminus(H)       (for testing)
  procedure RG(Err: in Radius; Rpp,R: in Weights;
               H1,Psi0: in Hamilton; H2: out Hamilton);                   -- Plain RG transformation   (for testing)

  procedure Nash1(Err: in Radius; R: in Weights;
                   H1,H: in Hamilton; Hr: out Hamilton);                  -- First order Iminus elimination
  procedure DNash1(Err: in Radius; R: in Weights;
                   Hdq,Hdp,F: in Hamilton; Fr: out Hamilton);             -- Derivative of Nash1
  procedure RG1(Err: in Radius; R: in Weights;
                Psi0,Hq,Hzp,H1,H: in Hamilton; Hr: out Hamilton);         -- RG with Nash1 instead of Nash
  procedure DRG1(Err: in Radius; R: in Weights;
                 Hq,Hzp,H1,Hdq,Hdp,Hdi: in Hamilton;
                 F: in out Hamilton; Fr: out Hamilton);                   -- Derivative of RG1 at H1,  using Hdi
  procedure DRG1Norms(Err: in Radius; R: in Weights;
                     Hq,Hzp,H1,Hdq,Hdp: in Hamilton; W: in Sparse);       -- Partial norm of MM1, complements MM1Norm

  procedure WCheck(W: in out Sparse);                                     -- Check that W is legal
  procedure MM1(Err: in Radius; R: in Weights; W: in out Sparse;
                Hq,Hzp,H1,Hdq,Hdp,Hdi,F: in Hamilton; Fr: out Hamilton);  -- Fr := (W+DRG1*(I-W))F  using Hdi
  procedure MM1Norm(Err: in Radius; R: in Weights; W: in out Sparse;
                    Hq,Hzp,H1,Hdq,Hdp,Hdi: in Hamilton;
                    J1,J2,Lf: in Positive; S: out Scalar);                -- Partial norm of MM1, complements DRG1Norms

  procedure PsiBounds(Hq,Hzp: in Hamilton; Rpp: in Weights;
                      B,Rpmin: out Weights);                              -- Bounds for Psi
  procedure TildeSoM(R: in Weights; W: in out Sparse;
                     H1: in Hamilton; H2: in HBall; H3: out Hamilton);    -- H3 := ~S(H1,M*H2),  M=I-W
  procedure TildeRoM(R: in Weights; W: in out Sparse;
                 Hq,Hzp,H1: in Hamilton; H2: in HBall; H3: out Hamilton); -- H3 := ~R(H1,M*H2),  M=I-W

  procedure NMInit(R: in Weights; H1: in Hamilton;
                F1: in HBall; A,B,C: out Scalar);                         -- A := ||H1||'_R, B := ||H1||_R, C := ||F1||_R
  procedure NMStep(N: in Positive; Rp,R: in Weights;
          A1,B1,C1: in Scalar; DA,DB,C2: out Scalar; OK: out Boolean);    -- Bounds for one Nash-Moser step
  procedure NMIter(Rp,R: in Weights;
           A1,B1,C1: in Scalar; DB: out Scalar; OK: out Boolean);         -- Bound for infinitely many Nash-Moser steps

  procedure NminusN1(R: in Weights;
                     H1: in Hamilton; F1: in HBall; F2: out HBall);       -- F2 contains N(H1+F1)-N1(H1+F1)
  procedure RminusR1oM(R: in Weights; W: in out Sparse;
               Psi0,Hq,Hzp,H1: in Hamilton; H2: in HBall; H3: out HBall); -- H3 contains R(H1+M*H2)-R1(H1+M*H2)

--- canonical transformations stuff

  type UBall is
    record
      R: Weights;
      Ex,Ez: Scalar;
    end record;

  procedure Comp(U1,U2: in UBall; U3: out UBall);                         -- U3 := U1 o U2
  procedure NormalizeM(W: in out Sparse; H: in Hamilton; F: in Hball;
                       H2: out Hamilton; C: out Scalar);                  -- also returns scaling C
  procedure NminusN1(R: in Weights; H1: in Hamilton; F1: in HBall;
                     F2: out HBall; U: out UBall);                        -- also computes U
  procedure RminusR1oM(R: in Weights; W: in out Sparse;
                       Psi0,Hq,Hzp,H1: in Hamilton; H2: in HBall;
                       H3: out HBall; C: out Scalar; U: out UBall);       -- also returns scaling C and computes U
  procedure Write(Name: in String; U: in UBall;
                  Decimal: in Boolean := True);                           -- write U to Name_x and Name_z

private

  type BooArray is array(Scal_Fou.Freq1,Scal_Fou.Freq2,Scal_Fou.Indx3) of Boolean;

  function DminusCount return Integer;                                    -- Number of Iminus coefficients
  procedure HoToUerr(Hq,H: in Hamilton; Hu: in out Hamilton);             -- used by HoToU
  procedure ColumnList(W: in Sparse; B: out BooArray);                    -- Record column indices of W
  procedure ColNorm(R: in Weights; Mj,Nj,Kj: in Integer; H: in Hamilton;
                    E: in out Scalar);                                    -- Use to compute partial norm of W
  procedure Scale(S: in Scalar; H: in out HBall);                         -- H := H(.,.,S.)/S
  procedure Kn(Rp,R: in Weights; A,B,C,X: in Scalar;
               K: out Scalar; OK: out Boolean);                           -- Constant K used by NMStep
  procedure Sn(N: in Positive; A,B,C,Dn: in Scalar; S: out Scalar);       -- Parameter S used by NMStep

  use Scal_Fou_Ops;

  procedure HoToU(Psi,Hq,Hzp: in Hamilton;
                  H,Hu: in FouVecPtr; Full: in Boolean := True);          -- HoToU for each component
  procedure DRG1(Err: in Radius; R: in Weights;
                 Hq,Hzp,H1,Hdq,Hdp,Hdi: in Hamilton; F,Fr: in FouVecPtr); -- DRG1 for each component
  procedure MM1(Err: in Radius; R: in Weights; W: in out Sparse;
                Hq,Hzp,H1,Hdq,Hdp,Hdi: in Hamilton;
                Lf: in Positive; F,Fr: in FouVecPtr);                     -- MM1 for each component
  procedure MaxNorm(R: in Weights; Lf: in Positive;
                    F: in FouVecPtr; E: in out Scalar);                   -- Maximum of E and norms of F(1..Lf)

  Trunc: constant Boolean := IsNumeric(ScalZero);

end RG_Ops;