File : rg_ops.adb
with Messages, Reps.Ops, Lin.Ops, Numerics.Ops.Fou;
use Messages;
pragma Elaborate_All (Messages, Lin.Ops, Numerics.Ops.Fou);
package body RG_Ops is
package Num_Fou is new Numerics.Ops.Fou (Lmax => Lmax, Wmax => Wmax, Pdeg => Pdeg, Dho => Dho);
use Num_Fou, Scal_Fou, Scal_Fou_IO;
Fuzz: constant Rep := Half**23;
Exp1: constant Scalar := Exp(ScalOne);
InvTheta: constant Scalar := Lambda(0,-1);
function IMin(I1,I2: Integer) return Integer renames Integer'Min;
function IMax(I1,I2: Integer) return Integer renames Integer'Max;
pragma Inline (Imin,Imax);
function DminusCount return Integer is
D: Integer := 0;
begin
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
D := D+Imin(Dminus(M,N),Dmax(M,N))+1;
end loop;
end loop;
return D;
end DminusCount;
procedure ShowParam is
begin
if Trunc then
Show(" Scalar = Numeric");
else
Show(" Scalar = Interval");
end if;
Show(" Lmax =",Lmax);
Show(" Wmax =",Wmax);
Show(" Size1 =",Size1);
Show(" Size =",Size2);
Show(" SizeE =",SizeE);
Show(" Pdeg =",Pdeg);
Show(" Dho =",Dho);
Show(" DmaxCount =",DmaxCount);
Show(" ErrCount =",ErrCount);
Show(" DminusCount =",DminusCount);
Show("Hamilton'Size =",Hamilton'Size);
Show(" Cpp =",Cpp);
Show(" Sigma =",Sigma);
Show(" Kappa =",Kappa);
end ShowParam;
function Up(R: Rep) return Rep is
begin
return Up(R,ScalZero);
end Up;
function Approx(Hs: Hamilton) return Polynom is
Hn: Polynom;
begin
SetZero(Par(Hs),Rho(Hs),Hn);
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
for K in 0 .. Deg(M,N,Hs) loop
SetComponent(M,N,K,Approx(Component(M,N,K,Hs)),Hn);
end loop;
end loop;
end loop;
return Hn;
end Approx;
procedure GetNumeric(Hn: in Polynom; Hs: out Hamilton) is
begin
SetZero(Par(Hn),Rho(Hn),Hs);
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
for K in 0 .. Deg(M,N,Hn) loop
SetComponent(M,N,K,Scal(Component(M,N,K,Hn)),Hs);
end loop;
end loop;
end loop;
SetPoly(True,Hs);
end GetNumeric;
function GetNumeric(Hn: Polynom) return Hamilton is
Hs: Fourier;
begin
GetNumeric(Hn,Hs);
return Hs;
end GetNumeric;
procedure Inv(Err: in Radius; R: in Weights; H1: in Hamilton; H2: out Hamilton) is
begin
if Trunc then
NumInv(Err,H1,H2);
SetRho(R,H2);
return;
end if;
GetNumeric(NumInv(Err,Approx(H1)),H2);
Inv(R,H1,H2);
end Inv;
function Inv(Err: Radius; R: Weights; H: Hamilton) return Hamilton is
Hr: Fourier;
begin
Inv(Err,R,H,Hr);
return Hr;
end Inv;
procedure FixComp3(Err: in Radius; R3: in Radius; H: in Hamilton; Hz: out Hamilton) is
--- works only if H is a polynomial in z, see below
Four: constant Rep := Two*Two;
R: Weights;
begin
TraceEnter("FixComp3");
R := (Rho(H).Q,R3);
if Trunc then
NumFixComp3(Err,H,Hz);
SetRho(R,Hz);
TraceLeave;
return;
end if;
GetNumeric(NumFixComp3(Err,Approx(H)),Hz);
SetRho(R,Hz);
declare
Ea,Eb,Ec,Er,Sa1: Scalar;
H3: Fourier;
begin
Copy(Hz,H3);
AddComponent(0,0,1,ScalOne,H3); -- H3(.,.,z)=z+Hz(.,.,z)
Ec := Norm(R,Comp3(H,H3)+Hz); -- -H(.,.,H3) should be approx Hz
R.P := Sup(Norm(H3));
Ea := Norm(R,Der3(R.P,H));
Sa1 := ScalOne-Ea;
Er := Two*Ec/Sa1; -- checks also that a<1
AddBall(Er,H3);
R.P := Sup(Norm(H3));
if Rho(H).P<R.P then
Message("Fouriers.Ops.FixComp3",Domain_Violation);
end if;
Eb := Norm(R,DerDer3(H)); -- complains if H is not a polynomial in z
if Sup(Four*Eb*Ec)>Inf(Sqr(Sa1)) then
Message("Fouriers.Ops.FixComp3",Domain_Violation);
end if;
H3 := Comp3(Der3(R.P,H),H3); -- now H3 contains A(xi)
AddComponent(0,0,0,-ScalOne,H3);
R.P:=R3; -- reset domain to ...
AddBall(Ec*Norm(Inv(Err,R,H3)),Hz);
end;
TraceLeave;
end FixComp3;
procedure HoToUerr(Hq,H: in Hamilton; Hu: in out Hamilton) is
--- assumes (and uses) that H is resonant
use Components, Components_Ops;
S0: constant Scalar := Theta+ScalOne; -- same as Theta^2
S2: constant Scalar := Sqr(InvTheta);
R: constant Weights := Rho(Hq); -- corresponds to (R2",B3) in notes
Rp: constant Weights := Rho(H); -- corresponds to R' in notes
K0: Integer;
L0: Rep;
G,SigHat,UHat: Scalar;
Pe: Taylor;
begin
if Poly(H) then return; end if;
UHatEtc(Hq,Rp,R,SigHat,UHat); -- domains get checked here
for I in FreqE loop
if ErrDeg(I,H) >= 0 then -- since H is resonant ...
K0 := IntCeiling(Scal(I)/Kappa); -- either p_degree >= K0 ...
L0 := Rep'Min(Rep(I),Inf(Scal(I)/Sigma)); -- or else |Omega*Nu| >= L0
G := GammaMax(R.Q*InvTheta/Rp.Q,SigHat/Rp.Q,Rp.Q*Scal(L0));
ErrComp(I,H,Pe);
Scale(S2,Pe);
Mult(S0,Pe);
for K in 0 .. Imin(K0-1,Dho-1) loop -- p_degree < K0 implies |Omega*Nu| >= L0
MultComponent(K,G,Pe); -- which gives contraction factor G
end loop;
for K in Dho .. Imin(K0-1,Pdeg) loop
MultComponent(K,Max(G,UHat*S2**(K0-K)),Pe); -- for degrees >= K0 did not scale enough
end loop;
for K in K0 .. Pdeg loop
MultComponent(K,UHat,Pe);
end loop;
MultErrComp(Max(UHat*S2**K0,G),Pe); -- for degrees >= K0 did not scale enough
AddErrComp(0,Pe,Hu);
end if;
end loop;
exception
when Undefined => Message("RG_Ops.HoToUerr",Parameter_Error);
end HoToUerr;
procedure HoToU(Psi,Hq,Hzp,H: in Hamilton; Hu: out Hamilton; Full: in Boolean := True) is
--- assuming Iminus(H)=0
--- make sure Psi has at least domain of Hq
R: constant Weights := Rho(Hq); -- corresponds to (R2",B3) in notes
begin
TraceEnter("HoToU");
CompOmegaTmu0(R,H,Hq,Hu,False); -- False means ...
HoToUerr(Hq,H,Hu); -- estimating errors here
if Full then Sub(Psi,Hu); end if;
Comp3(Hzp,Hu);
TraceLeave;
end HoToU;
procedure HoToU(Psi,Hq,Hzp: in Hamilton; H,Hu: in FouVecPtr; Full: in Boolean := True) is
--- assuming Iminus(H)=0
--- make sure Psi has at least domain of Hq
Lf: constant Positive := Imin(H'Last,Hu'Last);
R: constant Weights := Rho(Hq); -- corresponds to (R2",B3) in notes
begin
TraceEnter("HoToU");
CompOmegaTmu0(R,H,Hq,Hu,False); -- False means ...
for L in 1 .. Lf loop
HoToUerr(Hq,H(L),Hu(L)); -- estimating errors here
if Full then Sub(Psi,Hu(L)); end if;
end loop;
Comp3(Hzp,Hu);
TraceLeave;
end HoToU;
procedure MProd(W: in out Sparse; H1: in Hball; H2: out Hball; C2: out Scalar) is
use Lin.Ops;
OO2: constant Integer := Join(0,0,2);
Mi,Ni,Ki,Mj,Nj,Kj: Integer;
S,E: Scalar;
Fz: array(Indx3) of Scalar;
Ed,Eo: array(Freq1,Freq2,Indx3) of Scalar;
begin
Fz(0) := ScalOne;
for K in 1 .. Pdeg loop
Fz(K) := H1.R.P*Fz(K-1);
end loop;
for M in Freq1 loop
for N in Freq2 loop
for K in Indx3 loop
Ed(M,N,K) := ScalOne;
Eo(M,N,K) := ScalZero;
end loop;
end loop;
end loop;
Cleanup(W);
C2 := ScalZero;
for L in 1 .. Length(W) loop
UnJoin(W(L).J,Mj,Nj,Kj);
if W(L).I=W(L).J then
S := ScalOne-Scal(W(L).R);
Ed(Mj,Nj,Kj) := Abs(S); -- diagonal terms
else
UnJoin(W(L).I,Mi,Ni,Ki);
S := W(L).R*NormFactor(H1.R.Q,Mi,Ni,True)*Fz(Ki);
Eo(Mj,Nj,Kj) := Eo(Mj,Nj,Kj)+Abs(S); -- offdiagonal terms
end if;
if W(L).I=OO2 then
if W(L).J=OO2 then
S := (ScalOne-Scal(W(L).R))/Fz(2);
else
S := Scal(W(L).R)/(NormFactor(H1.R.Q,Mj,Nj,True)*Fz(Kj));
end if;
C2 := Max(Abs(S),C2);
end if;
end loop;
E := ScalZero;
for M in Freq1 loop
for N in Freq2 loop
for K in Indx3 loop
S := Ed(M,N,K)+Eo(M,N,K)/(NormFactor(H1.R.Q,M,N,True)*Fz(K));
E := Max(S,E);
end loop;
end loop;
end loop;
Show("Norm of M ",Sup(E));
Show("C2 factor ",Sup(C2));
H2.R := H1.R;
H2.E := E*H1.E;
C2 := Ball(C2)*H1.E;
end MProd;
procedure Prod(A: in Scalar; W: in out Sparse; H1: in Hamilton; H2: out Hamilton) is
--- not optimal unless Coeff returns same as Component
Mi,Ni,Ki,Mj,Nj,Kj: Integer;
begin
SetZero(Par(H1),Rho(H1),H2);
for L in reverse 1 .. Length(W) loop
UnJoin(W(L).I,Mi,Ni,Ki);
UnJoin(W(L).J,Mj,Nj,Kj);
AddComponent(Mi,Ni,Ki,W(L).R*Coeff(Mj,Nj,Kj,H1),H2);
end loop;
SetPoly(True,H2);
MultAdd(A,H1,H2);
end Prod;
procedure Normalize(H1: in Hamilton; H2: out Hamilton) is
S: Scalar;
begin
TraceEnter("Normalize");
S := Coeff(0,0,2,H1)/Cpp;
Scale(Inv(S),H1,H2);
S := Component(0,0,2,H2); --- rest is cleanup
SetComponent(0,0,2,Scal(Cpp),H2);
if not (Dho>2 or else Poly(H1)) then
AddComponent(0,0,2,ResetCenter(S),H2);
SetPoly(Poly(H1),H2);
end if;
S := Component(0,0,0,H2);
SetComponent(0,0,0,ScalZero,H2);
if not (Dho>0 or else Poly(H1)) then
AddComponent(0,0,0,ResetCenter(S),H2);
SetPoly(Poly(H1),H2);
end if;
TraceLeave;
end Normalize;
procedure NormalizeM(W: in out Sparse; H: in Hamilton; F: in Hball; H2: out Hamilton) is
S: Scalar;
F2: HBall;
Ht: Fourier;
begin
TraceEnter("NormalizeM");
MProd(W,F,F2,S); -- F2 := (I-W)*F, S := Coeff(0,0,2,F2)
Copy(H,Ht);
SetRho(Min(Rho(H),F2.R),Ht);
AddBall(F2.E,Ht); -- Ht := H+F2
S := Coeff(0,0,2,H)+Ball(S); -- Coeff(0,0,2,Ht)
S := S/Cpp; -- continue as in Normalize
Scale(Inv(S),Ht,H2);
Show("Scaling ",Inf(S),Sup(S),17);
S := Component(0,0,2,H2);
SetComponent(0,0,2,Scal(Cpp),H2);
if Dho <= 2 then
AddComponent(0,0,2,ResetCenter(S),H2);
end if;
S := Component(0,0,0,H2);
SetComponent(0,0,0,ScalZero,H2);
if Dho=0 then
AddComponent(0,0,0,ResetCenter(S),H2);
end if;
TraceLeave;
end NormalizeM;
procedure DNormalize(H,F: in Hamilton; F2: out Hamilton) is
PolyF2: constant Boolean := Poly(H) and Poly(F);
C: constant Scalar := Coeff(0,0,2,F);
S: Scalar;
begin
TraceEnter("DNormalize");
S := Coeff(0,0,2,H)/Cpp;
Scale(Inv(S),F,F2);
if C /= ScalZero then -- avoid by using Hdi, since ZxDer3 ...
declare
Ht: Fourier;
begin
ZxDer3(H,Ht); -- is implemented only for polynomials
Sub(H,Ht);
Scale(Inv(S),Ht);
S := C/(Cpp*S); -- or Coeff(0,0,2,F2)/Cpp
Mult(-S,Ht); -- could set Ht.P to False
SetPoly(PolyF2,Ht);
Add(Ht,F2);
end;
end if;
S := Component(0,0,2,F2);
SetComponent(0,0,2,ScalZero,F2);
if not (Dho>2 or else PolyF2) then
AddComponent(0,0,2,ResetCenter(S),F2);
SetPoly(PolyF2,F2);
end if;
S := Component(0,0,0,F2);
SetComponent(0,0,0,ScalZero,F2);
if not (Dho>0 or else PolyF2) then
AddComponent(0,0,0,ResetCenter(S),F2);
SetPoly(PolyF2,F2);
end if;
TraceLeave;
end DNormalize;
procedure UpsiAux(Err: in Radius; R: in Weights; Psi: in Hamilton; Hq,Hzp: out Hamilton) is
RP: Radius;
begin
TraceEnter("UpsiAux");
if not Poly(Psi) then
Message("RG_Ops.UpsiAux",Not_Implemented);
end if;
Gig(Sigma,Psi,Hq);
SetRho((R.Q,One),Hq); -- reasonable for FixComp3
FixComp3(Err,R.P,Hq,Hzp);
AddComponent(0,0,1,ScalOne,Hzp); -- Hzp(q,z)=z'
Dig(Kappa,Psi,Hq);
RP := Sup(Norm(R,Hzp))+Fuzz;
SetRho((R.Q,RP),Hq); -- Hq is a polynomial
Show("Hzp.R=",Rho(Hzp).Q,Rho(Hzp).P,17);
Show(" Hq.R=",Rho(Hq).Q,Rho(Hq).P,17);
TraceLeave;
end UpsiAux;
function NormPrime(R: Weights; H: Hamilton) return Scalar is
E1,E2: Scalar;
Ht: Fourier;
begin
OmegaDq(R.Q,H,Ht);
E1 := Norm(R,Ht)/R.P/Kappa;
Der3(R.P,H,Ht);
E2 := Norm(R,Ht)/Sigma;
return E1+E2;
end NormPrime;
procedure Nash(Err: in Radius; R: in Weights; H: in out Hamilton) is
--- for testing NminusN1 and RG - not used in proof
A: Rep;
F1,F2: HBall;
H1,H2: Fourier;
begin
TraceEnter("Nash");
H1 := GetNumeric(Approx(H)); -- contains a single element of H
Diff(H,H1,H2);
F1.R := Rho(H);
F1.E := Norm(H2)+Norm(Iminus(Dminus,H1));
Iplus(Dminus,H1); -- H1+F1 contains H, Iminus(H1)=0
Show(" |(I-)H| =",Sup(F1.E));
NminusN1(R,H1,F1,F2); -- F2 contains N(H)-N1(H)
Show("|N(H)-N1(H)| =",Sup(F2.E));
A := Sup(NormPrime(R,H1));
if A>Half then
F1.E:= (Scal(A)/(ScalOne-Scal(A)))*F1.E;
end if; -- now F1.E bounds norm of DN1(F)=N1(H)-H1
Copy(H1,H);
SetRho(R,H);
AddBall(F1.E+F2.E,H); -- add N(H)-H1 to H1
TraceLeave;
end Nash;
procedure RG(Err: in Radius; Rpp,R: in Weights; H1,Psi0: in Hamilton; H2: out Hamilton) is
--- for testing and completeness - not used in proof
Hq,Hzp,Ht: Fourier;
begin
TraceEnter("RG");
Normalize(H1,Ht);
UpsiAux(Err,Rpp,Psi0,Hq,Hzp);
HoToU(Psi0,Hq,Hzp,Ht,H2);
ZeroConstant(H2);
Nash(Err,R,H2);
ZeroConstant(H2);
Show("|R(h)-h| = ",Sup(Norm(R,H2-H1)));
TraceLeave;
end RG;
procedure Nash1(Err: in Radius; R: in Weights; H1,H: in Hamilton; Hr: out Hamilton) is
Hdq,Hdp: Fourier;
begin
Gradient(R,H1,Hdq,Hdp); -- Omega*Grad_q(H1) and D_p(H1)
DNash1(Err,R,Hdq,Hdp,H-H1,Hr);
Add(H1,Hr);
end Nash1;
procedure DNash1(Err: in Radius; R: in Weights; Hdq,Hdp,F: in Hamilton; Fr: out Hamilton) is
--- Fr := (I+)F-(I+)A(Psi) with (I-)[A+1]Psi=(I-)F
--- A(Psi)=Hdp*Gig(Psi)-Hdq*Dig(Psi)
A: Rep;
E1,E2,E: Scalar;
Psi0: Fourier;
begin
TraceEnter("DNash1");
GetNumeric(NumInvAplus1(Err,Dminus,Approx(Hdq),Approx(Hdp),Approx(F)),Psi0);
SetRho(R,Psi0);
Copy(F,Fr);
SetRho(R,Fr);
MultAdd( ScalOne,Hdq,Dig(Kappa,Psi0),Fr);
Multadd(-ScalOne,Hdp,Gig(Sigma,Psi0),Fr); -- Fr := F-A(Psi0)
if not Trunc then
E1 := Norm(R,Hdq)/R.P/Kappa;
E2 := Norm(R,Hdp)/Sigma;
A := Sup(E1+E2);
if not (A<One or Min(Rho(F),R)=R) then
Message("RG_Ops.DNash1",Domain_Violation);
end if;
E := Norm(R,Iminus(Dminus,Fr-Psi0))*Scal(A)/(ScalOne-Scal(A));
AddBall(E,Fr);
end if;
Iplus(Dminus,Fr);
TraceLeave;
end DNash1;
procedure RG1(Err: in Radius; R: in Weights; Psi0,Hq,Hzp,H1,H: in Hamilton; Hr: out Hamilton) is
Ht: Fourier;
begin
TraceEnter("RG1");
Normalize(H,Hr);
HoToU(Psi0,Hq,Hzp,Hr,Ht);
ZeroConstant(Ht);
Show("|LS(H)-H| = ",Sup(Norm(Ht-H)));
Nash1(Err,R,H1,Ht,Hr);
ZeroConstant(Hr);
Show("|R1(H)-H| = ",Sup(Norm(R,Hr-H)));
TraceLeave;
end RG1;
procedure DRG1(Err: in Radius; R: in Weights; Hq,Hzp,H1,Hdq,Hdp,Hdi: in Hamilton; F: in out Hamilton; Fr: out Hamilton) is
Spp: constant Scalar := Component(0,0,2,F);
begin
SetComponent(0,0,2,ScalZero,F); -- remove z**2 term
DNormalize(H1,F,Fr);
HoToU(Hq,Hq,Hzp,Fr,F,False); -- 1st argument arbitrary
MultAdd(Spp,Hdi,F); -- add image of z**2 term
ZeroConstant(F);
DNash1(Err,R,Hdq,Hdp,F,Fr);
ZeroConstant(Fr);
end DRG1;
procedure DRG1(Err: in Radius; R: in Weights; Hq,Hzp,H1,Hdq,Hdp,Hdi: in Hamilton; F,Fr: in FouVecPtr) is
--- destroys input F.all
Lf: constant Positive := Imin(F'Last,Fr'Last);
Spp: Components.ScalVec(1 .. Lf);
begin
for L in 1 .. Lf loop
Spp(L) := Component(0,0,2,F(L));
SetComponent(0,0,2,ScalZero,F(L)); -- remove z**2 term
DNormalize(H1,F(L),Fr(L));
end loop;
HoToU(Hq,Hq,Hzp,Fr,F,False); -- 1st argument arbitrary
for L in 1 .. Lf loop
MultAdd(Spp(L),Hdi,F(L)); -- add image of z**2 term
ZeroConstant(F(L));
DNash1(Err,R,Hdq,Hdp,F(L),Fr(L));
ZeroConstant(Fr(L));
end loop;
end DRG1;
procedure ColumnList(W: in Sparse; B: out BooArray) is
Mj,Nj,Kj: Integer;
begin
for M in Freq1 loop
for N in Freq2 loop
for K in Indx3 loop
B(M,N,K) := False;
end loop;
end loop;
end loop;
for L in 1 .. Length(W) loop
UnJoin(W(L).J,Mj,Nj,Kj);
B(Mj,Nj,Kj) := True;
if Kj >= Dmax(Mj,Nj) then
Message("RG_Ops.ColumnList",Not_Implemented);
end if;
end loop;
end ColumnList;
procedure ColNorm(R: in Weights; Mj,Nj,Kj: in Integer; H: in Hamilton; E: in out Scalar) is
S: Scalar;
begin
Show("Col",Join(Mj,Nj,Kj),False);
S := Norm(R,H)/(NormFactor(R.Q,Mj,Nj)*Scal(R.P)**Kj);
Show(" * ",Sup(S));
E := Max(E,S);
end ColNorm;
procedure DRG1Norms(Err: in Radius; R: in Weights; Hq,Hzp,H1,Hdq,Hdp: in Hamilton; W: in Sparse) is
--- estimate terms for MM1Norm that are unaffected by W, except error terms H.E
use Components;
Kh: Integer;
M: Freq1;
E: Scalar := ScalZero;
C,S1,S2: Scalar;
P1: Taylor;
Rt: Weights;
HasColumn: BooArray;
Hw1,Hw2,Hc1,Hs1,Hc2,Hs2,Hmc,Hms,Hmnc,Hmns: Fourier;
Hp3: FouVecPtr := new FouVec(1 .. Pdeg);
begin
ColumnList(W,HasColumn); -- columns to be skipped
C := Coeff(0,0,2,H1);
S1 := Sqr(Theta+ScalOne)*C/Cpp; -- Theta/Mu (Theta+1=Theta^2)
S2 := (-Cpp)*Sqr(InvTheta)/C; -- -Theta*Mu
SetZero(P1);
SetComponent(0,ScalOne,P1); -- P1(z)=1
if SupAbs(S2*Norm(R,Hzp))>R.P then
Message("RG_Ops.DRG1Norms",Domain_Violation);
end if;
Powers(R,Hzp,Hp3);
Rt := Rho(Hq); -- from now on mostly CompOmegaTmu0, except ...
CosSin((-InvTheta)*Hq,Hw1,Hw2,Hc1,Hs1);
CosSin(Sqr(InvTheta)*Hq,Hw1,Hw2,Hc2,Hs2);
SetZero(0,Rt,Hmc);
SetComponent(0,0,0,ScalOne,Hmc);
SetZero(1,Rt,Hms);
for Mabs in 0 .. Size1 loop
if Mabs>0 then
IncCS(Rt,Hc1,Hs1,Hw1,Hw2,Hmc,Hms);
end if;
for I in 0 .. 1 loop
M := (2*I-1)*Mabs;
Hms := -Hms;
if M<0 or else I=1 then
Copy(Hmc,Hmnc);
Copy(Hms,Hmns);
for N in Nmin(M) .. Nmax(M) loop
if N > 0 then
IncCS(Rt,Hc2,Hs2,Hw1,Hw2,Hmnc,Hmns);
end if;
SetZero(0,Rt,Hw2);
MultAdd(S1,N,M+N,P1,Hmnc,Hw2,Dmax,True); -- includes Nontrivial scaling ...
MultAdd(-S1,N,M+N,P1,Hmns,Hw2,Dmax,True);
Comp3(Hw2,Hp3,Hw1); -- and composition with Hzp
for K in Dminus(M,N)+1 .. Dmax(M,N)-1 loop
if not (HasColumn(M,N,K) or else (K=0 and then M=0 and then N=0)) then
if K=0 then
DNash1(Err,R,Hdq,Hdp,Hw1,Hw2); -- part of composition with Hzp
else
DNash1(Err,R,Hdq,Hdp,S2**K*Hp3(K)*Hw1,Hw2); -- part of composition with Hzp
end if;
ZeroConstant(Hw2);
ColNorm(R,M,N,K,Hw2,E); -- norm for monomials
end if;
end loop;
Kh := Imax(Dminus(M,N)+1,Dmax(M,N));
if (Kh mod Pdeg)=0 then
Copy(Hw1,Hw2);
else
Prod(R,Hp3(Kh mod Pdeg),Hw1,Hw2);
end if;
for J in 1 .. Kh/Pdeg loop
Mult(Hp3(Pdeg),Hw2);
end loop;
DNash1(Err,R,Hdq,Hdp,Ball(S2**Kh)*Hw2,Hw2); -- part of composition with Hzp
ZeroConstant(Hw2);
ColNorm(R,M,N,Kh,Hw2,E); -- norm for higher order terms in p
end loop;
end if;
end loop;
end loop;
Show("DRG1Norms: partial norm =",Sup(E));
end DRG1Norms;
procedure WCheck(W: in out Sparse) is
Lsave: constant Integer := Length(W);
Mi,Ni,Ki,Mj,Nj,Kj: Integer;
begin
for L in 1 .. Lsave loop
UnJoin(W(L).I,Mi,Ni,Ki);
UnJoin(W(L).J,Mj,Nj,Kj);
if Ki<0 or else Kj<0
or else (not InRange(Mi,Ni,Ki)) or else (not InRange(Mj,Nj,Kj))
or else (Ki <= Dminus(Mi,Ni)) or else (Kj <= Dminus(Mj,Nj)) then
Message("RG_Ops.WCheck",Index_Error);
end if;
end loop;
Lin.Ops.Cleanup(W,True);
if Length(W) /= Lsave then
Message("RG_Ops.WCheck",Index_Error);
end if;
end WCheck;
procedure MM1(Err: in Radius; R: in Weights; W: in out Sparse; Hq,Hzp,H1,Hdq,Hdp,Hdi,F: in Hamilton; Fr: out Hamilton) is
--- Fr is W+DRG1*(I-W) applied to F
Ft: Fourier;
begin
Prod(-ScalOne,W,F,Fr);
Neg(Fr,Ft);
DRG1(Err,R,Hq,Hzp,H1,Hdq,Hdp,Hdi,Ft,Fr);
Prod(ScalZero,W,F,Ft);
Add(Ft,Fr);
end MM1;
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) is
Ft: FouVecPtr := new FouVec(1 .. Lf);
begin
for L in 1 .. Lf loop
Prod(-ScalOne,W,F(L),Fr(L));
Neg(Fr(L),Ft(L));
end loop;
DRG1(Err,R,Hq,Hzp,H1,Hdq,Hdp,Hdi,Ft,Fr);
for L in 1 .. Lf loop
Prod(ScalZero,W,F(L),Ft(L));
Add(Ft(L),Fr(L));
end loop;
end MM1;
procedure MaxNorm(R: in Weights; Lf: in Positive; F: in FouVecPtr; E: in out Scalar) is
S: Scalar;
begin
for L in 1 .. Lf loop
S := Norm(R,F(L));
E := Max(E,S);
Show("norm",Sup(S));
end loop;
end MaxNorm;
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) is
--- consider only column space of W, and higher order q-error terms
I: constant Integer := Imin(Lmax,Wmax);
J,L: Integer := 0;
E: Scalar;
HasColumn: BooArray;
F,Fr: FouVecPtr := new FouVec(1 .. Lf);
begin
ColumnList(W,HasColumn); --- checks also assumption mentioned below
S := ScalZero;
for K in 0 .. Dho loop --- loop over q-errors
J := J+1;
if J in J1 .. J2 then
L := L+1;
SetZero(0,R,F(L));
E := Inv(Scal(R.P)**K);
SetErrComp(I,K,Ball(E),F(L)); --- q-error, fill into F
Show("err",Join(I,0,K));
if L=Lf then
MM1(Err,R,W,Hq,Hzp,H1,Hdq,Hdp,Hdi,Lf,F,Fr);
MaxNorm(R,Lf,Fr,S);
L := 0;
end if;
end if;
end loop;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
for K in Dminus(M,N)+1 .. Dmax(M,N)-1 loop --- loop over basis vectors, assuming ...
if HasColumn(M,N,K) then --- False unless K<Dmax(M,N)
J := J+1;
if J in J1 .. J2 then --- only vectors J1 .. J2
L := L+1;
SetZero(0,R,F(L));
E := Inv(NormFactor(R.Q,M,N)*Scal(R.P)**K);
if K<Pdeg then
SetComponent(M,N,K,E,F(L)); --- monomial, fill into F
SetPoly(True,F(L));
else --- W should not be this large
SetComponent(M,N,K,Ball(E),F(L)); --- p-errors, fill into F
end if;
Show("col",Join(M,N,K));
if L=Lf then
MM1(Err,R,W,Hq,Hzp,H1,Hdq,Hdp,Hdi,Lf,F,Fr);
MaxNorm(R,Lf,Fr,S);
L := 0;
end if;
end if;
end if;
end loop;
end loop;
end loop;
if L > 0 then
MM1(Err,R,W,Hq,Hzp,H1,Hdq,Hdp,Hdi,L,F,Fr);
MaxNorm(R,L,Fr,S);
end if;
Show("partial norm:",Sup(S));
end MM1Norm;
procedure Scale(S: in Scalar; H: in out HBall) is
begin
H.R.P := Inf(Scal(H.R.P)/Abs(S));
H.E := H.E/Abs(S);
end Scale;
procedure TildeSoM(R: in Weights; W: in out Sparse; H1: in Hamilton; H2: in HBall; H3: out Hamilton) is
C1,C2,Eps,Eps0,E: Scalar;
Hb: HBall;
Ht: Hamilton;
begin
TraceEnter("TildeSoM");
if not Poly(H1) then
Message("RG_Ops.TildeSoM",Not_Implemented);
end if;
C1 := Coeff(0,0,2,H1);
MProd(W,H2,Hb,C2); -- Hb := (I-W)*H2, C2 := Coeff(0,0,2,Hb)
Eps := C2/C1;
Scale(Inv(C1/Cpp),Hb);
Eps0 := (Scal(Hb.R.P)-Scal(R.P))/Hb.R.P;
if Inf(Eps0-Eps) <= Zero then
Message("RG_Ops.TildeSoM",Domain_Violation);
end if;
Normalize(H1,Ht);
ScaleRem2(Eps,Ht,H3);
SetRho(R,H3);
E := Eps*(ScalOne+Inv(Eps0-Eps))*Hb.E; -- bound on ScaleRem1 of Ht
AddBall(E,H3);
TraceLeave;
end TildeSoM;
procedure PsiBounds(Hq,Hzp: in Hamilton; Rpp: in Weights; B,Rpmin: out Weights) is
begin
B.P := Sup(Norm(Rpp,Hzp));
B.Q := Sup(Norm((Rpp.Q,B.P),Hq));
Rpmin.Q := Sup(Rpp.Q*InvTheta+B.Q*InvTheta);
Rpmin.P := Sup(B.P*Sqr(InvTheta));
end PsiBounds;
procedure TildeRoM(R: in Weights; W: in out Sparse; Hq,Hzp,H1: in Hamilton; H2: in HBall; H3: out Hamilton) is
Err: constant Radius := Half**30;
B,Rp: Weights;
Ht,Hdq,Hdp: Hamilton;
begin
TraceEnter("TildeRoM");
PsiBounds(Hq,Hzp,R,B,Rp); -- Rp is minimal domain needed for HoToU
Rp := (Rp.Q+Fuzz,Rp.P+Fuzz);
TildeSoM(Rp,W,H1,H2,H3); -- Uses up as much domain as possible
HoToU(Hq,Hq,Hzp,H3,Ht,False); -- 1st argument arbitrary
ZeroConstant(Ht);
Gradient(R,H1,Hdq,Hdp);
DNash1(Err,R,Hdq,Hdp,Ht,H3);
ZeroConstant(H3);
TraceLeave;
end TildeRoM;
procedure NMInit(R: in Weights; H1: in Hamilton; F1: in HBall; A,B,C: out Scalar) is
begin
if Min(Min(Rho(H1),F1.R),R) /= R then
Message("RG_Ops.NMInit",Domain_Violation);
end if;
A := Scal(Up(Sup(NormPrime(R,H1))));
B := Scal(Up(Sup(Norm(R,H1))));
C := Scal(Up(Sup(F1.E)));
end NMInit;
procedure Kn(Rp,R: in Weights; A,B,C,X: in Scalar; K: out Scalar; OK: out Boolean) is
A0,Tn: Scalar;
begin
OK := False;
K := ScalZero;
A0 := ScalOne-A;
if Inf(A0)<=Zero then return; end if;
Tn := C/(Sigma*X*A0);
if Sup(Tn)>=Half then return; end if;
OK := True;
K := ((Scal(2)-A)*X+B/Sigma)/(Sigma*(ScalOne-Tn)*Sqr(A0));
K := K*GammaMax(Scal(R.Q)/Rp.Q,ScalZero,ScalZero);
end Kn;
procedure Sn(N: in Positive; A,B,C,Dn: in Scalar; S: out Scalar) is
R: Rep;
A0,K,S3: Scalar;
begin
S := Scal(Half);
A0 := ScalOne-A;
if Inf(A0)<=Zero then
Show("hmm, taking S=1/2");
return;
end if;
K := B/Sqr(Sigma*A0);
S3 := Two*K*A0*C/Sqr(Dn);
if N>1 then S3 := S3/A; end if;
R := Reps.Ops.Exp(Reps.Ops.Log(Approx(S3))/3.0);
if (R<=Zero) or (R>=Half) then
Show("hmm, taking S=1/2");
return;
end if;
S := Scal(R);
end Sn;
procedure NMStep(N: in Positive; Rp,R: in Weights; A1,B1,C1: in Scalar; DA,DB,C2: out Scalar; OK: out Boolean) is
A0,DrQ,DrP,Dn,K,Bt,S: Scalar;
begin
OK := False;
DrQ := Scal(Rp.Q)-Scal(R.Q);
DrP := Scal(Rp.P)-Scal(R.P);
if Inf(DrQ)<=Zero or Inf(DrP)<=Zero then return; end if;
Dn := Min(Kappa*(Rp.P*DrQ)/Sigma,DrP);
if Inf(A1)<Zero or Inf(B1)<Zero or Inf(C1)<Zero then return; end if;
A0 := ScalOne-A1;
if Inf(Rep(3)*A0)<=One or Sup(Rep(6)*B1)>=One then return; end if;
Kn(Rp,R,A1,B1,C1,Dn,K,OK);
if not OK then return; end if;
C2 := K*Sqr(C1/Dn);
Bt := C1/A0;
if N>1 then Bt := Bt*A1; end if;
DB := Bt+C2;
Sn(N,A1,B1,C1,Dn,S);
Kn(Rp,R,A1,B1,C1,S*Dn,K,OK);
if not OK then return; end if;
DA := (Bt+K*Sqr(C1/(S*Dn)))/(Sigma*(ScalOne-S)*Dn);
DA := DA*GammaMax(Scal(R.Q)/Rp.Q,ScalZero,ScalZero);
end NMStep;
procedure NMIter(Rp,R: in Weights; A1,B1,C1: in Scalar; DB: out Scalar; OK: out Boolean) is
DrQ,DrP,G,Dn,Gmax,Amax,Bmax,Cmax: Scalar;
begin
OK := False;
DrQ := Scal(Rp.Q)-Scal(R.Q);
DrP := Scal(Rp.P)-Scal(R.P);
if Inf(DrQ)<=Zero or Inf(DrP)<=Zero then return; end if;
if Inf(A1)<Zero or Inf(B1)<Zero or Inf(C1)<Zero then return; end if;
G := GammaMax(Scal(R.Q)/Rp.Q,ScalZero,ScalZero);
Dn := Min(Kappa*(R.P*DrQ)/Sigma,DrP);
Gmax := Scal(6)/Rep(5);
Amax := Scal(2)/Rep(3);
Bmax := ScalOne/Rep(6)-Two*DrP;
Cmax := Sqr(Dn)/Rep(24);
if not (G<Gmax and A1<Amax and B1<Bmax and C1<Cmax) then return; end if;
OK := True;
DB := Rep(13)*C1/Rep(3);
end NMIter;
procedure NminusN1(R: in Weights; H1: in Hamilton; F1: in HBall; F2: out HBall) is
Iter: constant Integer := 4;
Jmax: constant Integer := 6;
W0: constant array(1..Jmax) of Rep := (0.00, 0.005, 0.01, 0.02, 0.05, 0.10);
W1: constant array(1..Jmax) of Rep := (0.08, 0.160, 0.32, 0.68, 0.84, 0.92);
OK: Boolean;
Rp: Weights;
procedure NMSteps(I: in Integer; R1,R3: in Weights; A1,B1,C1: in Scalar; E: out Scalar) is
R2: Weights;
DA,DB,A2,B2,C2,E2,E3: Scalar;
begin
E := ScalOne;
for J in 1 .. Jmax loop -- try out different domains
if I=0 then
R2.Q := R1.Q-W0(J)*(R1.Q-R3.Q);
R2.P := R1.P-W0(J)*(R1.P-R3.P);
NMInit(R2,H1,F1,A2,B2,C2);
E2 := ScalZero;
else
R2.Q := R1.Q-W1(J)*(R1.Q-R3.Q);
R2.P := R1.P-W1(J)*(R1.P-R3.P);
NMStep(I,R1,R2,A1,B1,C1,DA,DB,C2,OK);
if not OK then goto End_Loop; end if;
A2 := A1+DA;
B2 := B1+DB;
if I=1 then
E2 := C2;
else
E2 := DB;
end if;
if Sup(E2) >= Sup(E) then goto End_Loop; end if;
end if;
if I<Iter then
NMSteps(I+1,R2,R3,A2,B2,C2,E3);
else
NMIter(R2,R3,A2,B2,C2,DB,OK);
if not OK then goto End_Loop; end if;
E3 := DB;
end if;
E := Scal(Rep'Min(Sup(E2+E3),Sup(E)));
<<End_Loop>> null;
end loop;
end NMSteps;
begin
TraceEnter("NminusN1");
Rp := Min(Rho(H1),F1.R);
NMSteps(0,Rp,R,ScalZero,ScalZero,ScalZero,F2.E);
if F2.E=ScalOne then
Message("RG_Ops.NminusN1",Domain_Violation);
end if;
F2.R := R;
TraceLeave;
end NminusN1;
procedure RminusR1oM(R: in Weights; W: in out Sparse; Psi0,Hq,Hzp,H1: in Hamilton; H2: in HBall; H3: out HBall) is
E: Scalar;
H,Ht: Hamilton;
begin
TraceEnter("RminusR1oM");
NormalizeM(W,H1,H2,Ht);
HoToU(Psi0,Hq,Hzp,Ht,H); -- Hq and Hzp define domain
Sub(H1,H); -- H := L(S(H1+M*H2))-H1
ZeroConstant(H);
E := Norm(H);
Show(" |H|=",Sup(E));
NminusN1(R,H1,(Rho(H),E),H3); -- H3 := (N-N1)(H1+H)
Show("H3.E=",Sup(H3.E));
TraceLeave;
end RminusR1oM;
--- canonical transformations stuff
procedure Comp(U1,U2: in UBall; U3: out UBall) is
R: Weights;
G: Scalar;
begin
R.Q := Sup(Scal(U2.R.Q)+U2.Ex);
R.P := Sup(Scal(U2.R.P)+U2.Ez);
if R.Q>U1.R.Q or R.P>U1.R.P then
Message("RG.Ops.Comp",Domain_Violation);
end if;
G := GammaMax(Scal(U2.R.Q)/U1.R.Q,ScalZero,ScalZero);
U3.R := U2.R;
U3.Ex := U2.Ex+U1.Ex*G;
U3.Ez := U2.Ez+U1.Ez*G;
end Comp;
procedure NormalizeM(W: in out Sparse; H: in Hamilton; F: in Hball; H2: out Hamilton; C: out Scalar) is
S: Scalar;
F2: HBall;
Ht: Fourier;
begin
TraceEnter("NormalizeM");
MProd(W,F,F2,S); -- F2 := (I-W)*F, S := Coeff(0,0,2,F2)
Copy(H,Ht);
SetRho(Min(Rho(H),F2.R),Ht);
AddBall(F2.E,Ht); -- Ht := H+F2
S := Coeff(0,0,2,H)+Ball(S); -- Coeff(0,0,2,Ht)
S := S/Cpp; -- continue as in Normalize
Scale(Inv(S),Ht,H2);
Show("Scaling ",Inf(S),Sup(S),17);
C := S; -- new
S := Component(0,0,2,H2);
SetComponent(0,0,2,Scal(Cpp),H2);
if Dho <= 2 then
AddComponent(0,0,2,ResetCenter(S),H2);
end if;
S := Component(0,0,0,H2);
SetComponent(0,0,0,ScalZero,H2);
if Dho=0 then
AddComponent(0,0,0,ResetCenter(S),H2);
end if;
TraceLeave;
end NormalizeM;
procedure NminusN1(R: in Weights; H1: in Hamilton; F1: in HBall; F2: out HBall; U: out UBall) is
--- modified version of NminusN1 above that computes U
Iter: constant Integer := 4;
Jmax: constant Integer := 6;
W0: constant array(1..Jmax) of Rep := (0.00, 0.005, 0.01, 0.02, 0.05, 0.10);
W1: constant array(1..Jmax) of Rep := (0.08, 0.160, 0.32, 0.68, 0.84, 0.92);
OK: Boolean;
Rp: Weights;
procedure NMSteps(I: in Integer; R1,R3: in Weights; A1,B1,C1: in Scalar; E: out Scalar; U13: out UBall) is
G: constant Scalar := Scal(0.2785); -- related to Gamma
R2: Weights;
DA,DB,A2,B2,C2,E2,E3,F: Scalar;
U12,U23: UBall;
begin
E := ScalOne;
for J in 1 .. Jmax loop -- try out different domains
if I=0 then
R2.Q := R1.Q-W0(J)*(R1.Q-R3.Q);
R2.P := R1.P-W0(J)*(R1.P-R3.P);
NMInit(R2,H1,F1,A2,B2,C2);
E2 := ScalZero;
else
R2.Q := R1.Q-W1(J)*(R1.Q-R3.Q);
R2.P := R1.P-W1(J)*(R1.P-R3.P);
NMStep(I,R1,R2,A1,B1,C1,DA,DB,C2,OK);
if not OK then goto End_Loop; end if;
A2 := A1+DA;
B2 := B1+DB;
if I=1 then
E2 := C2;
else
E2 := DB;
end if;
if Sup(E2) >= Sup(E) then goto End_Loop; end if;
end if;
U12.R := R2; -- new
U12.Ex := C1/(ScalOne-A1)/R2.P/Kappa; -- new
U12.Ez := C1/(ScalOne-A1)/Sigma; -- new
if I<Iter then
NMSteps(I+1,R2,R3,A2,B2,C2,E3,U23);
else
NMIter(R2,R3,A2,B2,C2,DB,OK);
if not OK then goto End_Loop; end if;
E3 := DB;
F := Exp(G*Log(Scal(R2.Q)/R3.Q))*Scal(16)/Scal(3); -- new
U23.R := R3; -- new
U23.Ex := F*C2/(ScalOne-A2)/R3.P/Kappa; -- new
U23.Ez := F*C2/(ScalOne-A2)/Sigma; -- new
end if;
if Sup(E2+E3)<Sup(E) then
E := Scal(Sup(E2+E3));
Comp(U12,U23,U13); -- new
end if;
<<End_Loop>> null;
end loop;
end NMSteps;
begin
TraceEnter("NminusN1");
Rp := Min(Rho(H1),F1.R);
NMSteps(0,Rp,R,ScalZero,ScalZero,ScalZero,F2.E,U);
if F2.E=ScalOne then
Message("RG_Ops.NminusN1",Domain_Violation);
end if;
F2.R := R;
TraceLeave;
end NminusN1;
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) is
--- Using new versions of NormalizeM and NminusN1
E: Scalar;
H,Ht: Hamilton;
begin
TraceEnter("RminusR1oM");
NormalizeM(W,H1,H2,Ht,C);
HoToU(Psi0,Hq,Hzp,Ht,H); -- Hq and Hzp define domain
Sub(H1,H); -- H := L(S(H1+M*H2))-H1
ZeroConstant(H);
E := Norm(H);
Show(" |H|=",Sup(E));
NminusN1(R,H1,(Rho(H),E),H3,U); -- H3 := (N-N1)(H1+H)
Show("H3.E=",Sup(H3.E));
TraceLeave;
end RminusR1oM;
procedure Write(Name: in String; U: in UBall; Decimal: in Boolean := True) is
H: Hamilton;
begin
SetZero(1,U.R,H);
AddBall(U.Ex,H);
Write(Name & "_x",H,Decimal);
SetZero(0,U.R,H);
AddBall(U.Ez,H);
Write(Name & "_z",H,Decimal);
end Write;
end RG_Ops;