SIMD GLM schemes

Design features for SIMD aware GLM schemes.

Needs to be aware and adjustable for vector length (Subscript[m, 1]) and total core count/shared memory core limit (Subscript[m, 2]).  For Gossamer CUDA (Subscript[m, 1])=32  and (Subscript[m, 2])=15*32=480.  We will use the multiplier (Subscript[m, 3])=15 rather than the total count (Subscript[m, 2])=480.

Needs to always use (Subscript[m, 1]) function evaluations at a time.

Need to operate on and update either a large set of internals or have a large multi-value vector.   These are much larger than have been considered to date.  I will call this design variable (Subscript[m, 4]).

Most of the arithmetic needs to be designed for (Subscript[m, 1]) and needs to be under a simple control stucture.  I believe that the hyper-power algorithms are very usefull for this.

The plan is to first write a collection of WR algorithms as a stationary GLM.

The plan is then to try to move the WR algorithm by a variable number of sub steps.

In order to change step size we need to be able to interpolate points in computation.

I have a sub-step size of h and a super step size of H=(Subscript[m, 4])h

Fixed Step-Size Stationary WR

For y'=f[y] (Subscript[y, 0])=y0 on the interval 0≤t≤H with h=(H/((Subscript[m, 4])-1)) WR is the vector iteration for a 
	({{(Subscript[WR, 0]), =, Init[y0,h,f]}, {(Subscript[WR, n+1]), =, y0+Quad[(Subscript[WR, n]),h,f]}})
for any choice of initialization "Init" operator and any choice of quadrature opertaor "Quad".   For clarity 
	Init:{(ℝ^dim),ℝ,(ℝ^dim)-> (ℝ^dim)}->(ℝ^((Subscript[m, 4])×dim)) 
	Quad:{(ℝ^((Subscript[m, 4])×dim)),(ℝ^dim)-> (ℝ^dim)}->(ℝ^((Subscript[m, 4])×dim)) 
The original WR algorithm is simply to replicate the initial condition
	Init[(Subscript[y, 0]),h,f]:=(Underscript[{(Subscript[y, 0]),…,(Subscript[y, 0])}, (Underscript[⏟, (Subscript[m, 4])])])
and (using zero based indexing to match the math) approximate the quadrature in the integral equation
	y[t]=(Subscript[y, 0])+(∫_(τ=0)^t)f[y[τ]dτ
usning the trapezoid rule which gives
	 F=f[(Subscript[W, n])]
Quad[(Subscript[W, n]),h,f]:=(h/2){0,(Subscript[F, 0])+(Subscript[F, 1]),(Subscript[F, 0])+2(Subscript[F, 1])+(Subscript[F, 2]),(Subscript[F, 0])+2(Subscript[F, 1])+2(Subscript[F, 2])+(Subscript[F, 3]),…,(Subscript[F, 0])+2(∑_(i=0)^((Subscript[m, 4])-1))((Subscript[F, i])+(Subscript[F, i+1]))+(Subscript[F, (Subscript[m, 4])])}
which can be written concisely as 
	Quad[(Subscript[W, n]),h,f]=h (Subscript[QTrap, (Subscript[m, 4])]).f[(Subscript[W, n])]
where (Subscript[QTrap, m]):=(1/2)(({{0, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 0, 0}, {1, 2, 1, 0, 0, 0}, {1, 2, 2, 1, 0, 0}, {⋮, ⋮,  , ⋱, ⋱, 0}, {1, 2, ⋯, ⋯, 2, 1}}))
Now there is no point in writing a stationary algorithm and recomputing the 0th point as (Subscript[y, 0]).  We might as well compute
	Quad2[(Subscript[W, n]),h,f]:=(h/2){(Subscript[F, 0])+(Subscript[F, 1]),(Subscript[F, 0])+2(Subscript[F, 1])+(Subscript[F, 2]),(Subscript[F, 0])+2(Subscript[F, 1])+2(Subscript[F, 2])+(Subscript[F, 3]),…,(Subscript[F, 0])+2(∑_(i=0)^((Subscript[m, 4])-1))((Subscript[F, i])+(Subscript[F, i+1]))+(Subscript[F, (Subscript[m, 4])]),(OverVector[α]).F + (OverVector[β]).(Subscript[W, n])}
where I am deliberately using the notation from LMM schemes for the predictor at the front.  This moves the scheme forward through one step and would be very easy to code.  This is essentially an explicit LMM scheme as a predictor with (Subscript[m, 4])-1 very simple correction steps using the WR algorithm. As long as each step involves the same number of FMADs it should be possible to pipeline the operations.  If the constants in the FMADs are the same you can even save some memory. There is no reason that you can not decide wether to advance or not based on the convergence of one or more entries. There is no reason you can not try to advance more than one point although to do that you need additional LMM like steps.

The order of the WR scheme at the end is estimated using the integral equation.  This is a little different from the order estimate for an advancing scheme.  WR accuracy is expressed by moving everything accross to the other side and substituting the exact answer into the numerical scheme at the stationary point.  An advancing WR has the same defn it just has to have a weasle word in it that it does not leave untill the point IS stationary. So for clarity I assume 
	y[t]-(Subscript[y, 0])-(∫_(τ=0)^t)f[y[τ]dτ=0
and I have a "vector" Y≈{y[0],y[h],y[2h],y[3h],…y[((Subscript[m, 4])-1)h]}
	Y-((Subscript[y, 0])-Quad.F[Y])=0
So what I do is I substitute the correct values in to the discrete equations and expand in h as h->0 to match the advancing asymtotic order definition.  To match the WR relaxation definition I let (Subscript[m, 4])->∞ and hold H=((Subscript[m, 4])-1)/h constant. As always they should be different by 1 power.

Tests and Trials and …

Setup tests

GravityP29[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=1/81.45, d1,d2},
d1=((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(-3/2));
d2=((((y2^2)+(y3^2)+(((y1+μ))^2)))^(-3/2));
{y4,y5,y6,
2y5 + y1 - μ(y1+μ-1)d1 -(1-μ)(y1+μ)d2,
-2y4 + y2 - μ y2 d1 -(1-μ)y2 d2,
- μ y3 d1 -(1-μ)y3 d2
}
]
IC1 = {0.994,0,0,0,-2.001585106370825224,0};
IC=IC1;
TMax=17.06521656015796; 
ySol[t_]=y[t]/.NDSolve[{
y'[t]==GravityP29[y[t]],
y[0]==IC
},
y[t],
{t,0, TMax}][[1]];
ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax},
Epilog->{PointSize[0.02],Red, Point[ySol[TMax/3][[{1,2}]]]}
]

Trapezoid Integration

Trying our WR Trapezoid algorithm on a small chunk.  This is a low accuracy scheme and appears to be unable to track the initial corner accurately.  So I moved it.

Clear[ATrap,h]
TSmallMax=TMax/5;
TStart=TMax/3;
{nWR,nIter}={64,128};
WR0 = Table[ySol[TStart],{nWR}];
ATrap=ConstantArray[0,{nWR,nWR}];
F=Array[f,{nWR}];
Do[
ATrap[[i,i]]=ATrap[[i,1]]=1;
Do[ATrap[[i,j]]=2,{j,2,i-1}],
{i,2,nWR}];
h=TSmallMax/(nWR-1);
WR=WR0;
F=Map[GravityP29,WR];
Data=Table[
F=Map[GravityP29,WR];
WR =WR0+(h/2)ATrap.F,
{nIter}];
Show[
ListPlot[WR[[All,{1,2}]]],
ParametricPlot[ySol[t][[{1,2}]],{t,TStart,TStart + TSmallMax},
PlotStyle->Red],
FrameLabel->{"y1","y2"},
PlotLabel->StringForm[
"nIter=``, nWR=`` in Trap",
nIter, nWR] 
]
ConvData=Table[
Map[Log[10,#]&,Map[Norm,Data[[i]]-WR]],
{i, 2,Length[Data]}];
ListContourPlot[ConvData,
(*DataRange->{{0, TSmallMax},{1, nIter}},*)
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),
PlotLabel->StringForm["Algorithm=Trap with (Subscript[Log, 10]) contours"],
FrameLabel->{"n","t"},
Contours->Table[n,{n,-16,2,1}]
]

ErrData=Table[Norm[ySol[TStart+(i-1) h]-WR[[i]]],{i, 2, Length[WR]}];
FitFun[hh_]=Fit[Take[ErrData,5],{hh}, hh]
Show[
ListLogLogPlot[ 
ErrData,
PlotLabel->StringForm["Error last iterate: h=``",h],
PlotRange->All
],
LogLogPlot[
FitFun[hh],
{hh, 0.001, Length[ErrData]}]
]

The error looks linear in step count and the final error looks quadratic in h. 
{
{0.11,0.050},
{0.05, 0.010},
{0.026, 0.0034}
}

Data={{0.11,0.050},{0.05,0.010},{0.026,0.0034}};
ListPlot[
Table[
{h,err}=Data[[i]];
(err/(h^2)),
{i, Length[Data]}]
]

h

ListLogPlot[ Hmm=Table[Norm[ySol[(i-1) h]-WR[[i]]],{i,  Length[WR]}],
PlotLabel->"Error last iterate",
PlotRange->All]

Take[
Hmm,3]

Conclusions are that the algorithm is converging fine but is not sufficiently accurate.

Mid Point Integration

The mid point formula for a quadratic going through the middle of 
	{(Subscript[t, n]),{(Subscript[y, n]),(Subscript[f, n])}}-{(Subscript[t, n+1]),{(Subscript[y, n+1]),(Subscript[f, n+1])}}
at (Subscript[t, n+1/2]) is 
	(Subscript[y, n+(1/2)])=(1/2)((Subscript[y, n+1])+(Subscript[y, n]))+(h/8)((Subscript[f, n+1])-(Subscript[f, n]))

Clear[AMid]
TSmallMax=TMax/25;
{nWR,nIter}={512,128};
WR0 = Table[IC,{nWR}];
AMid=ConstantArray[0,{nWR,nWR-1}];
Do[
Do[AMid[[i,j]]=1,{j,1,i-1}],
{i,2,nWR}]
h=TSmallMax/(nWR-1);
WR=WR0;
F=Map[GravityP29,WR];
Data=Table[
FTemp=Map[Gravityp29, WR];
yMid=Table[(1/2)(WR[[i]]+WR[[i+1]])+(h/8)(FTemp[[i+1]]-FTemp[[i]]),{i,1,nWR-1}];
FMid=Map[GravityP29,yMid];
WR =WR0+h AMid.FMid,
{nIter}];
Show[
ListPlot[WR[[All,{1,2}]]],
ParametricPlot[ySol[t][[{1,2}]],{t,0,TSmallMax},
PlotStyle->Red],
FrameLabel->{"y1","y2"},
PlotLabel->StringForm[
"nIter=``, nWR=`` in Mid",
nIter, nWR] 
]
ConvData=Table[
Map[Log[10,#]&,Map[Norm,Data[[i]]-WR]],
{i, Length[Data]}];
ListContourPlot[ConvData,
DataRange->{{0, TSmallMax},{1, nIter}},
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),
PlotLabel->StringForm["Algorithm=Mid with (Subscript[Log, 10]) contours"],
FrameLabel->{"n","t"}]

ySol[0]
WR[[1]]

ListLogPlot[ Table[Norm[ySol[(i-1) h]-WR[[i]]],{i,  Length[WR]}]]

Dimensions[WR]

Conclusions are that the algorithm is converging fine but is not sufficiently accurate.

Simp Integration

OK I am going to make the Simpsons rule by combining the Mid and Trap rules.

Clear[AMid]
TSmallMax=TMax/25;
{nWR,nIter}={512,128};
WR0 = Table[IC,{nWR}];
AMid=ConstantArray[0,{nWR,nWR-1}];
Do[
Do[AMid[[i,j]]=1,{j,1,i-1}],
{i,2,nWR}]
h=TSmallMax/(nWR-1);
WR=WR0;
F=Map[GravityP29,WR];
Data=Table[
FMid=Map[GravityP29,Table[((WR[[i]]+WR[[i+1]])/2),{i,1,nWR-1}]];
WR =WR0+h AMid.FMid,
{nIter}];
Show[
ListPlot[WR[[All,{1,2}]]],
ParametricPlot[ySol[t][[{1,2}]],{t,0,TSmallMax},
PlotStyle->Red],
FrameLabel->{"y1","y2"},
PlotLabel->StringForm[
"nIter=``, nWR=`` in Mid",
nIter, nWR] 
]
ConvData=Table[
Map[Log[10,#]&,Map[Norm,Data[[i]]-WR]],
{i, Length[Data]}];
ListContourPlot[ConvData,
DataRange->{{0, TSmallMax},{1, nIter}},
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),
PlotLabel->StringForm["Algorithm=Mid with (Subscript[Log, 10]) contours"],
FrameLabel->{"n","t"}]

Conclusions are that the algorithm is converging fine but is not sufficiently accurate.

Interpolation Stuff.

I need formulae for a quadratic through (Subscript[y, n]) with slope (Subscript[f, n]) and (Subscript[y, n+1])with slope (Subscript[f, n+1]). I then need to 
evaluate that at the mid point.

Clear[f,n,h]
Expand[
Simplify[InterpolatingPolynomial[ {
{(Subscript[t, n]),{(Subscript[y, n]),(Subscript[f, n])}},
{(Subscript[t, n])+h,{(Subscript[y, n+1]),(Subscript[f, n+1])}}
},
tt]/.tt->(Subscript[t, n])+h/2
]
]

the magic formula is (Subscript[y, n+(1/2)])=(1/2)((Subscript[y, n+1])+(Subscript[y, n]))+(h/8)((Subscript[f, n+1])-(Subscript[f, n]))

Organization

Closed and Open Newton Cotes

Closed Formula

It is easy to generate the appropriate interpolating polynomials and integrate them to produce schemes based on consistent interpolation schemes.

Clear[x0,h]
StencilSpread=1;
f[x_]=InterpolatingPolynomial[
Table[
{x0+i h,(Subscript[y, i]),(Subscript[f, i]),(((Subscript[f, i+1])-(Subscript[f, i-1]))/(2h))},
{i, -StencilSpread,StencilSpread}],
x] ;
Expand[Integrate[ f[x],{x,x0-StencilSpread h,x0+StencilSpread h}]]

GaussianMatrix

Needs["FunctionApproximations`"]
f[x_]=(E^(- 10x));
{xMin, xMax}={0,3}
g[x_]=RationalInterpolation[f[x],{x, 3,3},Range[0,3,0.5]]
Plot[ {f[x],g[x]},{x,xMin,xMax}]

LogPlot[ Abs[f[x]-g[x]],{x, xMin, xMax},
PlotRange->All]

g[0]
f[0]

f[x_]=a0 + a1 x + a2 (x^2);
{xMin, xMax}={1,3}
g[x_]=RationalInterpolation[f[x],{x, 3,3},{x,xMin,xMax}]

Clear[x0,h]
StencilSpread=1;
f[x_]=InterpolatingPolynomial[
Table[
{x0+i h,(Subscript[y, i]),(Subscript[f, i])},
{i, -StencilSpread,StencilSpread}],
x] ;
Expand[Integrate[ f[x],{x,x0-StencilSpread h,x0+StencilSpread h}]]

Clear[x0,h]
StencilSpread=1;
f[x_]=InterpolatingPolynomial[
Table[
{x0+i h,(Subscript[y, i]),(Subscript[f, i]),θ ((((Subscript[f, i+1])-(Subscript[f, i-1])))/(2h))+ (1-θ) ((((Subscript[y, i-1])-2 (Subscript[y, i])+(Subscript[y, i+1])))/(h^2))},
{i, -StencilSpread,StencilSpread}],
x];
Collect[Expand[Integrate[ f[x],{x,x0-StencilSpread h,x0+StencilSpread h}]],θ]

Visualization

I should check That I get what I think I get;

Clear[x0,h,g,x]
MaxOrder=12;
h=0.1;
g[x_]= Sum[RandomReal[{-12,12}] Cos[i 2 π x],{i, 0,MaxOrder}];
x0=RandomReal[{-1,1}];
θ=RandomReal[{0,1}];
Do[ (Subscript[y, i])=g[x0+i h]; (Subscript[f, i])=g'[x0+ i h],
{i,-MaxOrder, MaxOrder}];
StencilSpread=1;
f[x_]=InterpolatingPolynomial[
Table[
{x0+i h,(Subscript[y, i]),(Subscript[f, i]),θ(( (Subscript[f, i+1])-(Subscript[f, i-1]))/(2 h))+(1-θ)(( (Subscript[y, i+1])-2 (Subscript[y, i])+(Subscript[y, i-1]))/(h^2))},
{i, -StencilSpread,StencilSpread}],
x] ;
Plot[ {f[x],g[x]},{x,x0-(StencilSpread+1) h,x0+(StencilSpread+1) h},
PlotStyle->{Red, Blue},
Epilog->Map[Point, Table[{x0+i h,(Subscript[y, i])},{i, -StencilSpread,StencilSpread}]]
]
Plot[ {f'[x],g'[x]},{x,x0-(StencilSpread+1) h,x0+(StencilSpread+1) h},
PlotStyle->{Red, Blue},
Epilog->Map[Point, Table[{x0+i h,(Subscript[f, i])},{i, -StencilSpread,StencilSpread}]]
]
Plot[ {f''[x],g''[x]},{x,x0-(StencilSpread+1) h,x0+(StencilSpread+1) h},
PlotStyle->{Red, Blue},
Epilog->Map[Point, Table[{x0+i h,f''[x0+i h]},{i, -StencilSpread,StencilSpread}]]
]
fInt=Integrate[ f[x],{x,x0-StencilSpread h,x0+StencilSpread h}]
gInt=Integrate[ g[x],{x,x0-StencilSpread h,x0+StencilSpread h}]
100(fInt-gInt)/fInt

Length[ Collect[ f[x]/.{θ->0},x]]

Simplify[D[g[x],{x,8}]]
Simplify[D[f[x],{x,8}]]

Polynomials

olve

Clear[x0,h]
StencilSpread=1;
f[x_]=InterpolatingPolynomial[
Table[
{x0+i h,(Subscript[y, i]),(Subscript[f, i])},
{i, -StencilSpread,StencilSpread}],
x] ;
Expand[Integrate[ f[x],{x,x0-StencilSpread h,x0+StencilSpread h}]]

Clear[x0,h]
StencilSpread=1;
f[x_]=InterpolatingPolynomial[
Table[
{x0+i h,(Subscript[y, i]),(Subscript[f, i]),θ ((((Subscript[f, i+1])-(Subscript[f, i-1])))/(2h))+ (1-θ) ((((Subscript[y, i-1])-2 (Subscript[y, i])+(Subscript[y, i+1])))/(h^2))},
{i, -StencilSpread,StencilSpread}],
x];
Collect[Expand[Integrate[ f[x],{x,x0-StencilSpread h,x0+StencilSpread h}]],θ]