(* Calculating of the energy (including resonances) for helium isoelectronic 
    series by finding a stationary point of the variational functional *)
(* With advanced differentiation over a, b, c *)
(* With optimized performance *)
(* Using algebraic extrapolation for guessing *)
(* As a function of Z, without re-scaling *)

(********** INPUT STARTS **********)

{nmin,nmax}={3,3};
choice=0; (* 0 for i+j+k<=N, 1 for i+j^2+k^2<=N *)
ndig=96; (* Precision of calculations *)

(**********  INPUT ENDS  **********)

ndig1=Min[ndig-10,22]; (* AccuracyGoal in FindRoot *)
ndim=3; (* dimensionality - 3 or 5 or 7 ...*)
miter=100; (* max n. of iterations *)
zmin=0;
zmax=1;
zstep=0.005;
ipts=22; (* Number of intermediate points for alg. interpolation *)
maxdegree=6; (* Max. degree of algebraic interpolation *)
x=" ";
x3="   ";
eps=1/20 (1+I); (* for guessing *)
epseps=10^-11 (1+I); (* for guessing *)
{$MinPrecision,$MaxPrecision}={0,5000}; (* To avoid errors at running the following line *)
$MinPrecision=$MaxPrecision=ndig;
Off[Precision::"precsm"];
Off[Solve::"precsm"];
Clear[t];
chop[t_]=If[t==0,0,t];
{zmin,zmax,zstep}=Rationalize[{zmin,zmax,zstep}];
(* Guesses at Z = 1 *)
 eigstart=-0.527751; (* pre-estimated eigenvalue *)
 Clear[nm];
 abcstart=Which[ nm==0, {0.483743, 1.075018,-0.146564},
                 nm==1, {0.457994, 0.937309, 0.052705},
                 nm==2, {0.533515, 1.092917,-0.043595},
                 nm==3, {0.382676, 1.021425, 0.126637},
                 nm==4, {0.594582, 1.101342,-0.059193},
                 nm>=5, {0.339739, 1.107264, 0.197078}
 ]; (* Estimated for choice=0 (i+j+k<=N) *)
eigstart=SetPrecision[eigstart,ndig];
abcstart=SetPrecision[abcstart,ndig];
Print["{Nrange,choice,precision} = ",{{nmin,nmax},If[choice==0," i+j+k<=N "," i+j^2+k^2<=N "],ndig}];

p=ToFileName["..","packages"];
If[!MemberQ[$Path,p],$Path=Append[$Path,p]];
<<hylleraa`;
<<alginter`;
<<printing` (* no ";" *)

Clear[a,b,c,i,j,k,i1,j1,k1,i2,j2,k2,int1,int2,ints,int1a,int2a,intsa,int1b,int2b,intsb,int1c,int2c,intsc];
{T,V0,V1,S}=hmatrix[a,b,c,i1,j1,k1,i2,j2,k2,ndim];
(* Differentiating the matrix *)
{T,V0,V1,S}={T,V0,V1,S}/.{
   hylleraas`private`int1[i_,j_,k_]->int1[i,j,k][a,b,c],
   hylleraas`private`int2[i_,j_,k_]->int2[i,j,k][a,b,c],
   hylleraas`private`ints[i_,j_,k_]->ints[i,j,k][a,b,c]
		};
{Ta,V0a,V1a,Sa}=D[{T,V0,V1,S},a]/.{
		Derivative[1,0,0][int1[i_,j_,k_]][a,b,c]->int1a[i,j,k],
		Derivative[1,0,0][int2[i_,j_,k_]][a,b,c]->int2a[i,j,k],
		Derivative[1,0,0][ints[i_,j_,k_]][a,b,c]->intsa[i,j,k],
		int1[i_,j_,k_][a,b,c]->int1[i,j,k],
		int2[i_,j_,k_][a,b,c]->int2[i,j,k],
		ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
{Tb,V0b,V1b,Sb}=D[{T,V0,V1,S},b]/.{
		Derivative[0,1,0][int1[i_,j_,k_]][a,b,c]->int1b[i,j,k],
		Derivative[0,1,0][int2[i_,j_,k_]][a,b,c]->int2b[i,j,k],
		Derivative[0,1,0][ints[i_,j_,k_]][a,b,c]->intsb[i,j,k],
		int1[i_,j_,k_][a,b,c]->int1[i,j,k],
		int2[i_,j_,k_][a,b,c]->int2[i,j,k],
		ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
{Tc,V0c,V1c,Sc}=D[{T,V0,V1,S},c]/.{
		Derivative[0,0,1][int1[i_,j_,k_]][a,b,c]->int1c[i,j,k],
		Derivative[0,0,1][int2[i_,j_,k_]][a,b,c]->int2c[i,j,k],
		Derivative[0,0,1][ints[i_,j_,k_]][a,b,c]->intsc[i,j,k],
		int1[i_,j_,k_][a,b,c]->int1[i,j,k],
		int2[i_,j_,k_][a,b,c]->int2[i,j,k],
		ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
{T,V0,V1,S}={T,V0,V1,S}/.{
			int1[i_,j_,k_][a,b,c]->int1[i,j,k],
			int2[i_,j_,k_][a,b,c]->int2[i,j,k],
			ints[i_,j_,k_][a,b,c]->ints[i,j,k]
		};
Do[ (* nm = nmin, ..., nmax *)
Print[nm//Definition];
outf="t"<>ToString[nm]<>".dat";
If[FileType[outf]===File, (* If the data file already exists *)
 inpf=outf;
 data=ReadList[inpf,Number,RecordSeparators->"\n",RecordLists->True];
 data=ReadList[inpf,Number,RecordSeparators->"\n",RecordLists->True];
(* One line does not work correctly for unknown reason.
   It is related somehow to setting $MinPrecision and $MaxPrecision. *)
 {nmdata,choicedata,ndigdata}=data//First;
 data=Drop[data,1];
 dt=Transpose[data];
 zdt=dt[[1]]//Rationalize; (* Rounding Z to exact number *)
 dt=SetPrecision[dt,ndig];
 Print["Using the previous output file ",outf];
 Print["Z-previous: ",printF[zdt//First,8,4],"  - ",printF[zdt//Last,8,4],
    "   (",Length[zdt]," data points)"];
 {apar,bpar,cpar,eguess}=Transpose/@{{zdt,dt[[2]]+I chop/@dt[[3]]},
    {zdt,dt[[4]]+I chop/@dt[[5]]},
    {zdt,dt[[6]]+I chop/@dt[[7]]},
    {zdt,dt[[8]]+I chop/@dt[[9]]} };
 zlast=Last[zdt];
 zstart=If[zmax<zlast,zmax,zlast-zstep]
,  (* If the data file does not exist *)
 {apar,bpar,cpar,eguess}={{},{},{},{}};
 zstart=zmax;
 os=OpenWrite[outf,PageWidth->Infinity,FormatType->OutputForm];
 Write[os,nm,x3,choice,x3,ndig];
 Close[os];
];
Clear[ind,mmax];
mm=hmapping[nm,choice,ind,mmax];
Print[Definition[mm]];
L=M=Table[Null,{mm},{mm}];
(* *)
Clear[a,b,c];
time0=Timing[
hintegrals[a,b,c,ndim,nm];
(* Diffirentiating integrals *)
nmax=2 nm;
nmadd=2 (ndim-3);
nmax3=nmax+nmadd+3;
Do[
  Do[namax=n-nc;namax1=Floor[namax/2];
    Do[nb=namax-na;
       integ1[na-1,nb-1,nc-1]=integ2[nb-1,na-1,nc-1]=hylleraas`private`int1[na-1,nb-1,nc-1];
       integ2[na-1,nb-1,nc-1]=integ1[nb-1,na-1,nc-1]=hylleraas`private`int2[na-1,nb-1,nc-1];
       integs[na-1,nb-1,nc-1]=integs[nb-1,na-1,nc-1]=hylleraas`private`ints[na-1,nb-1,nc-1];
       integ1a[na-1,nb-1,nc-1]=integ2a[nb-1,na-1,nc-1]=
              D[hylleraas`private`int1[na-1,nb-1,nc-1],a]//Together;
       integ2a[na-1,nb-1,nc-1]=integ1a[nb-1,na-1,nc-1]=
              D[hylleraas`private`int2[na-1,nb-1,nc-1],a]//Together;
       integsa[na-1,nb-1,nc-1]=integsa[nb-1,na-1,nc-1]=
              D[hylleraas`private`ints[na-1,nb-1,nc-1],a]//Together;
       integ1b[na-1,nb-1,nc-1]=integ2b[nb-1,na-1,nc-1]=
              D[hylleraas`private`int1[na-1,nb-1,nc-1],b]//Together;
       integ2b[na-1,nb-1,nc-1]=integ1b[nb-1,na-1,nc-1]=
              D[hylleraas`private`int2[na-1,nb-1,nc-1],b]//Together;
       integsb[na-1,nb-1,nc-1]=integsb[nb-1,na-1,nc-1]=
              D[hylleraas`private`ints[na-1,nb-1,nc-1],b]//Together;
       integ1c[na-1,nb-1,nc-1]=integ2c[nb-1,na-1,nc-1]=
              D[hylleraas`private`int1[na-1,nb-1,nc-1],c]//Together;
       integ2c[na-1,nb-1,nc-1]=integ1c[nb-1,na-1,nc-1]=
              D[hylleraas`private`int2[na-1,nb-1,nc-1],c]//Together;
       integsc[na-1,nb-1,nc-1]=integsc[nb-1,na-1,nc-1]=
              D[hylleraas`private`ints[na-1,nb-1,nc-1],c]//Together;
      ,{na,0,namax1}],{nc,0,n}],{n,0,nmax3}];
]//First;Print["time-matrix-general form = ",time0];
(* *)
Do[ (* Z = Zstart, Zstart-Zstep, ..., Zmin *)
(* Guessing *)
mdata=Length[eguess];
dabc={eps,eps,eps};
ddabc={epseps,epseps,epseps};
Clear[n];
nn=Max[Floor[n/.Solve[(n+1)(n+2)/2-1==mdata,n]]];
nn=Min[nn,maxdegree]; (* up to max. degree *)
Print["Z = ",z//InputForm,"     degree of alg. f. = ",nn];
eig=If[mdata==0,eigstart,alginterpolation[eguess,nn,z,intermediatePoints->ipts] ];
abc1=If[mdata==0,abcstart,
      {alginterpolation[apar,nn,z,intermediatePoints->ipts],
       alginterpolation[bpar,nn,z,intermediatePoints->ipts],
       alginterpolation[cpar,nn,z,intermediatePoints->ipts]} ];
abc2=If[mdata<=1,abcstart+dabc,
      {alginterpolation[apar,nn-1,z,intermediatePoints->ipts],
       alginterpolation[bpar,nn-1,z,intermediatePoints->ipts],
       alginterpolation[cpar,nn-1,z,intermediatePoints->ipts]}+ddabc ];
guess={{atrial,{abc1[[1]],abc2[[1]]}},
       {btrial,{abc1[[2]],abc2[[2]]}},
       {ctrial,{abc1[[3]],abc2[[3]]}}};
(******************************* START find-root ******************************)
dampf=1;
iter=0;
abciter={};
rt=FindRoot[
iter++;
{a,b,c}=
   SetPrecision[{atrial,btrial,ctrial},ndig]//N; (* increase of precision *)
(* finding the eigenvalue *)
time1=Timing[
Clear[int1,int2,ints,int1a,int2a,intsa,int1b,int2b,intsb,int1c,int2c,intsc,na,nb,nc];
int1[na_,nb_,nc_]:=int1[na,nb,nc]=int2[nb,na,nc]=integ1[na,nb,nc];
int2[na_,nb_,nc_]:=int2[na,nb,nc]=int1[nb,na,nc]=integ2[na,nb,nc];
ints[na_,nb_,nc_]:=ints[na,nb,nc]=ints[nb,na,nc]=integs[na,nb,nc];
int1a[na_,nb_,nc_]:=int1a[na,nb,nc]=int2a[nb,na,nc]=integ1a[na,nb,nc];
int2a[na_,nb_,nc_]:=int2a[na,nb,nc]=int1a[nb,na,nc]=integ2a[na,nb,nc];
intsa[na_,nb_,nc_]:=intsa[na,nb,nc]=intsa[nb,na,nc]=integsa[na,nb,nc];
int1b[na_,nb_,nc_]:=int1b[na,nb,nc]=int2b[nb,na,nc]=integ1b[na,nb,nc];
int2b[na_,nb_,nc_]:=int2b[na,nb,nc]=int1b[nb,na,nc]=integ2b[na,nb,nc];
intsb[na_,nb_,nc_]:=intsb[na,nb,nc]=intsb[nb,na,nc]=integsb[na,nb,nc];
int1c[na_,nb_,nc_]:=int1c[na,nb,nc]=int2c[nb,na,nc]=integ1c[na,nb,nc];
int2c[na_,nb_,nc_]:=int2c[na,nb,nc]=int1c[nb,na,nc]=integ2c[na,nb,nc];
intsc[na_,nb_,nc_]:=intsc[na,nb,nc]=intsc[nb,na,nc]=integsc[na,nb,nc];
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=T-z V0+V1;
M[[m1,m2]]=M[[m2,m1]]=S;
,{m1,m2,mm}],{m2,mm}];
A=Inverse[M].L;
			]//First;(*Print["time-matrix = ",time1];*)
time2=Timing[
{val0,vec0}=
              Sort[Transpose[Eigensystem[A]],
                  Abs[#1[[1]]-eig]<Abs[#2[[1]]-eig]&][[1]];
			]//First;(*Print["time-eigenvalues = ",time2];*)
(* finding derivative of the eigenvalue d/da *)
time3=Timing[
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=Ta-z V0a+V1a;
M[[m1,m2]]=M[[m2,m1]]=Sa;
,{m1,m2,mm}],{m2,mm}];
vala=(vec0.L.vec0)/(vec0.M.vec0)-val0;
(* finding derivative of the eigenvalue d/db *)
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=Tb-z V0b+V1b;
M[[m1,m2]]=M[[m2,m1]]=Sb;
,{m1,m2,mm}],{m2,mm}];
valb=(vec0.L.vec0)/(vec0.M.vec0)-val0;
(* finding derivative of the eigenvalue d/dc *)
Do[{i2,j2,k2}=ind[m2];
Do[{i1,j1,k1}=ind[m1];
L[[m1,m2]]=L[[m2,m1]]=Tc-z V0c+V1c;
M[[m1,m2]]=M[[m2,m1]]=Sc;
,{m1,m2,mm}],{m2,mm}];
valc=(vec0.L.vec0)/(vec0.M.vec0)-val0;
			]//First;
delt={vala,valb,valc};
deltabs=Plus@@Abs/@delt;
deltpr=PaddedForm[deltabs//FortranForm,{3,1},ExponentFunction->Identity];
abciter=Append[abciter,{deltabs,{atrial,btrial,ctrial},val0}];
If[iter==1&&z==zstart,Print["Time (one step): matrix, eigenv., derivatives:  ",time1,", ",time2,", ",time3]];
Print[iter,"  --- ",deltpr,"  ",{printF[atrial,9,6],printF[btrial,9,6],printF[ctrial,9,6]}
     ];delt=={0,0,0},
      Sequence@@guess//Release,
WorkingPrecision->ndig,AccuracyGoal->ndig1,MaxIterations->miter,DampingFactor->dampf];
(******************************* END find-root ******************************)
{delt,abc0,eig0}=Sort[abciter,#1[[1]]<#2[[1]]&][[1]];
zpr=printF[z,8,4];
{apr,bpr,cpr}={printF[abc0[[1]]+0. I,ndig1+2,ndig1-2],
 printF[abc0[[2]]+0. I,ndig1+2,ndig1-2],
 printF[abc0[[3]]+0. I,ndig1+2,ndig1-2]};
enpr=printF[eig0+0. I,ndig1+4,ndig1];
deltpr=PaddedForm[delt//FortranForm,{3,1},ExponentFunction->Identity];
Print["===========   Z = ",zpr];
Print["Non-linear parameters:  ",{apr,bpr,cpr}];
Print["Energy:  ",enpr];
Print["L.h.s. of the equations (should be zero), n. of iter.:  ",deltpr,", ",iter];
If[iter>=miter,
 Print["Number of iterations exceeded the limit of ",miter," iterations."];
 Break[] ];
(*
If[iter>12&&nn==maxdegree,maxdegree=Min[12,maxdegree+1]]; (* Increase degree if too many iterations *)
If[iter<6,maxdegree=Max[5,maxdegree-1]]; (* Decrease degree *)
*)
os=OpenAppend[outf,PageWidth->Infinity,FormatType->OutputForm];
Write[os,zpr,x3,apr,x,bpr,x,cpr,x3,enpr,x3,deltpr];
Close[os];
{apar,bpar,cpar,eguess}=
  {Append[apar,{z,abc0[[1]]}],Append[bpar,{z,abc0[[2]]}],Append[cpar,{z,abc0[[3]]}],
   Append[eguess,{z,eig0}]};
,{z,zstart,zmin,-zstep}]; 
,{nm,nmin,nmax}];
If[$MachineType!="PC",Exit[]];