PROGRAM DIAG
C                                                  Version date: 28-Apr-2002
C
C     The program calculates the lowest 2 (3) eigenvalues of the Hamiltonian 
C     truncated to N states. If N>0 it also calculates Temple's inequality.  
C     dE0=(<H^2>-<H>^2)/(E1-E0). Table for b+i*bs with i=0 to nb-1. 
C     
C     Reduced (angular) matrix elements made with Mathematica are read from:
C     Htab2r: 'I N M X', with N<=M and I=0,2,4 (X=Hv[N,M]), I=1 (X=Hf[N,M]), 
C     I=3 (X=2L(L+7)=(l-3)(l+4)) and I=9 (END). The Hred^2 matrix elements
C     are read taken from HHtab2r: with I=10,12,14,16,18 (X=Hv^2[N,M]), 
C     I=11,13,15 X=(Hv*Hf)[N,M], I=20,22 (X=(Hf^2)[N,M]) and I=19 (END). 
C     Radial matrix elements <q,l|r^I|k,p> made with Mathematica are read 
C     from Rtab2 and RRtab2: 'I q p l k X' for k>2, k-l=0,1,2,3,4,5,6,8 and
C     q,p=0,1,....,pmax (but for k=l q<=p), with I=1,4 (X*Hf for dl=1 and 
C     X*Hv for dl=0,2,4), I=2,5,8 (X*Hf^2 for dl=0,2, X*Hv*Hf for dl=1,3,5 
C     and X*Hv^2 for dl=0,2,4,6,8) and I=0 (END).
C
C     IMPORTANT: entries in (H)Htab2r, (R)Rtab2 and Bjtab should occur only
C     once. Reset NM=nmax (420), LM=lmax (42) and PM=pmax (20) when required. 
C     plist contains the number of radial modes for each of the spherical
C     harmonics. Stops reading at "nmax p" -- nmax should match NM exactly.
C
C              Full F=2 Hamiltonian by combining angular and radial parts: 
C     Formats: N<0 b bs nb dr=0 Table for N-truncated H to 'tab-out' as
C              {b,  0.  ,E0,E1,0.} (requiring 'Htab2r', 'Rtab2' and 'plist')
C     Formats: N>0 b bs nb dr=0 for N-truncated H to 'tab-out' as
C              {b,E0-dE0,E0,E1,0.} (requiring in addition 'HHtab2r', 'RRtab2')
C
C              Full Hamiltonian as constructed in mathematica programmes,
C              Htab0 made with H0t.ma for F=0 or H2h.ma for F=2 (harmonic):
C     Formats: N b=w bs nb dr=-1 for N-truncated H (F=1-sg(w)) to 'tab-out' as
C              {w,  0.  ,E0,E1,0.} (requiring 'Htab0/2h', no 'HHtab0', N<=NM)
C     Formats: N b bs nb dr<-1 for N-truncated H (F=2) to 'tab-out' as
C              {b,E0-dE0,E0,E1,0.} (requiring 'Htab2' and 'HHtab2', N<=NM)
C
C              If 0<|dr|<1 diagonalize the reduced Hamiltonian at fixed r 
C              and calculate the norm of the derivative of the wavefunction, 
C              dZ/dr=(Z(r-dr)-Z(r+dr))/(2dr), Veff=E0+6/r^2+(dZ/dr)^2/2. 
C     Formats: N b=r bs nb 0>dr>-1 for N-truncated Hred (F=0) to 'tab-out' as
C              {r,E0-dE0,E0,E1,Veff} (requiring 'Htab0r' and 'HHtab0r')
C     Formats: N b=r bs nb dr>0 for N-truncated Hred (F=2) to 'tab-out' as
C              {r,E0-dE0,E0,E1,Veff} (requiring 'Htab2r' and 'HHtab2r')
C
C     For nb=1 wavefunction written to the file 'psi' and for dr=0, fn*fn and
C              fn are computed (diagonalizes Hred for each tabulated r values). 
C     For nb=0 avoids writing psi, and when dr=0 compute fn*fn (but put fn=0).
C     For nb>1 when bs=0 and N<0 (else nb->1) table of first nb eigenvalues.
C     For nb<0 (->nb=1) pruning is applied with level TP (|bs| when nonzero), 
C              until relative error below TE (|dr| when nonzero), or number 
C              of basis vectors decreases (but always restart after first run). 
C              It is useful to test by editing plist, P(I)=0->1, and restart.
C
C     Batch use: diag <tab-in >>tout; format of tab-in: N b bs nb dr
C     writes result to tab-out. Without '>>tout' ALSO on the screen!
C
C     Minimum requirement array declarations (for Htab2h may reset NM to NMT):
C     DOUBLE PRECISION AM(N(NC+1)),WORK(MAX(N(3NC+6),N**2)),Z(3N),HD(MAX(N,NM))
C     DOUBLE PRECISION D(NMT),H(NMT,NMT),HH(NM,NM),FF(NM,NM)
C     INTEGER L(NM),P(NM),NF(NM),NS(LM),NE(LM)
C
C     If all P(I)=P, NC=INT(4SQRT(NM+(P-1)/P)-3)P, e.g. for P=9 we get
C     N=3780,NC=711, and AM(2691360),WORK(14288400),Z(11340),HD(3780)
C     With pruning not all states are used, and pmax>9 is allowed when
C     N<=NNM and NC<=NCM, but NC>NCM need not be fatal as long as N<NNM.
C     Warn if NC>NCM or N>NNM and when PM=pmax is reset to (R)Rtab2 value. 
C     STOP if LM=lmax, NM=nmax do not match (R)Rtab2 and (H)tab2r values.
C
      INTEGER PM,PN,PT,L(420),P(420),NF(420),NS(42),NE(42)
      DOUBLE PRECISION X,T,B,BS,D1,D2,DD,TP,TE,DR
C
C     In case of limited memory, adjust the following defaults!
C
      DOUBLE PRECISION D(3080),H(3080,3080),HH(420,420),FF(420,420)
      DOUBLE PRECISION AM(2691360),WORK(14288400),Z(11340),HD(3780)
      NNM=3780
      NMT=3080
      NCM=711
      PM=20
C
C     Set remaining defaults. 
C
      TE=1.D-2
      TP=2.D-5
      NM=420
      LM=42
      II=0
      IQ=0
      IP=0
      MM=2
      IJOB=-2
C
C     Read parameters
C
      PRINT *,"Enter: N b bs nb dr (0 0 0 0 0 for usage)"
      READ *,N,B,BS,NB,DR
C
C     When N=0 give info (NB=0) on input format and ask if user wants
C     'plist' initialized; skip info for NB non-zero and set P(I)=|NB|.
C
      IF(N.NE.0)GOTO 3
      IF(NB.NE.0)GOTO 1
      PRINT *,"Tabulate from b, with nb steps of bs; N<0 no Temple"
      PRINT *,"dr>=1: {dr:E0-1}->{dr-10:E1-2}; dr<-1 uses (H)Htab2"
      PRINT *,"dr=-1, b=w|-w uses Htab0/2h; 0<|dr|<1, b=r, Veff(r)"
      PRINT *,"from chi0'=dchi0/dr: dr<0, (F=0) uses (H)Htab0r and"
      PRINT *,"dr>0 (F=2) uses (H)Htab2r, dr>=11 uses psi (created"
      PRINT *,"when dr<11, nb=1) & plist, resets dr=0 and computes"
      PRINT *,"fn*fn (nb=0) and fn (nb=1). dr=0 uses plist, Bjtab,"
      PRINT *,"(R)Rtab2, (H)Htab2r. nb<0 prunes: if |Z(n)|<tp=|bs|"
      PRINT *,"remove radial mode, else P->P+1, till error<te=|dr|"
      PRINT *,"or N stabalizes; tp=2.D-5 or te=1.D-2 when bs|dr=0."
      PRINT *,"Without Temple: bs=0, nb>2|3, gives nb eigenvalues."
      PRINT *,"N=0 initializes plist with P(I)=|nb|, if 0 it shows"
      PRINT *,"info & asks for number of radial modes at each Y_n."
      PRINT *,""
      PRINT *,"Initialize plist (no=0, yes=#radial modes<pmax)?"
      READ *,NB
      IF(NB.EQ.0)STOP
    1 OPEN(1,FILE='plist')
      PRINT *,"Initializing plist with P(I)=",ABS(NB)
      DO 2 I=1,NM
         WRITE (1,*) I,ABS(NB)
    2 CONTINUE
      WRITE (1,*) NB*NM,B,TE,TP
      CLOSE(1)
      STOP
C
C     When dr>1, reset dr=dr-10 & {E0,E1,Veff0}->(E1,E2,Veff1} (we do NOT 
C     allow for pruning in this case by nb->|nb|). When dr>=11, reset dr=0 
C     and compute fn from psi (IQ=1) (we do NOT allow for pruning by 
C     nb->1, when nb is non-zero). Reset defaults according to input 
C
    3 IF(DR.LT.11)GOTO 4
      IF(NB.NE.0)NB=1
      DR=0.
      IQ=1
    4 IF(DR.LT.1)GOTO 5
      IF(NB.LT.0)NB=-NB
      DR=DR-10.
      II=1
    5 K=1
      KK=1
      IF(N.LT.0)K=-1
      NN=K*N
      IF(II.EQ.1)MM=3
      IF(NB.LT.0)IP=1
      IF((BS.EQ.0).AND.(NB.GT.1).AND.(K.NE.1))MM=NB
      IF((BS.EQ.0).AND.(NB.NE.0))NB=1
      IF((K.NE.1).AND.(IP.EQ.0).AND.(NB.NE.0))IJOB=-1
      OPEN(3,FILE='tab-out')
      IF(IP.EQ.0)GOTO 6
      IF(DR.NE.0)TE=DABS(DR)
      IF(BS.NE.0)TP=DABS(BS)
      NB=1
      DR=0.
C
C     To make absolutely sure everything we assume to be zero, is zero!
C     Return to this point after pruning or making a table for dr non-zero.
C
    6 DO 7 I=1,NM
         DO 7 J=1,NM
            H(I,J)=0.
            HH(I,J)=0.
            FF(I,J)=0.
    7 CONTINUE
C
C     Read matrix elements of reduced Hamiltonian, I=dl
C     I=0,2,4 -> H=Hv/2, I=1 -> H=Hf/2, I=3 -> (l-3)(l+4).
C     T=b for Htab2, (-)w for Htab0(2h), r for Htab2(0)r, dr=0->T=1.
C     Htab2,Htab0,Htab0r,Htab2r for dr<-1, dr=-1, -1<dr<0, dr>=0.
C     Reset NM to NMT when using all of Htab2h (w<0) [or Htab2 (b<0)]
C
      T=1.
      IF(DR.NE.0)T=B
      IF(DR.GE.0)OPEN(1,FILE='Htab2r')
      IF(DR.LT.-1)OPEN(1,FILE='Htab2')
      IF((DR.EQ.-1).AND.(B.GE.0))OPEN(1,FILE='Htab0')
      IF((DR.EQ.-1).AND.(B.LT.0))T=-B
      IF((DR.LE.-1).AND.(B.LT.0))NM=NMT
      IF((DR.EQ.-1).AND.(B.LT.0))OPEN(1,FILE='Htab2h')
      IF((DR.LT.0).AND.(DR.GT.-1))OPEN(1,FILE='Htab0r')
   10 READ (1,*) I,N,M,X
      IF((I.EQ.9).OR.(M.GT.NM))GOTO 11
      IF(I.EQ.1)H(N,M)=T*X
      IF((I.EQ.5).AND.(DR.EQ.-1))H(N,M)=X/DSQRT(T)
      IF((I.EQ.0).OR.(I.EQ.2).OR.(I.EQ.4))H(N,M)=T**4*X/2
      IF((I.EQ.3).AND.(DR.EQ.0))L(N)=INT(0.5*SQRT(49.+4*X)-0.49)
      IF((I.EQ.3).AND.(DR.NE.0).AND.(N.EQ.M))HD(N)=X/T/T/2
C
C     Only when reading Htab0 add to existing H!
C
      IF((I.EQ.3).AND.(DR.NE.0).AND.(N.NE.M))H(N,M)=H(N,M)+X/T/T/2
      H(M,N)=H(N,M)
      GOTO 10
   11 CLOSE(1)
      IF(N.EQ.NM)GOTO 13
      IF(DR.LT.0)GOTO 12
      PRINT *,"Reset nmax=",NM," to Htab2r value ",N
      STOP
   12 IF((NN.LE.NM).OR.(M.LE.NM))GOTO 13
      PRINT *,"Warning: only",NM," states fit in H()"
      NN=NM
C
C     Read matrix elements of Hred*Hred if K=1 (N>0), I=dl+10, 
C     I=10,12,14,16,18 -> HH=Hvv/4, I=11,13,15 -> HH=Hvf; 
C     I=dl+20, I=20,22 -> FF=Hff (FF, to avoid mix with HH!)
C     HHtab(2,[0|2h,]0r,2r) for dr<-1, [dr=-1,] -1<dr<0, dr>=0 
C     If HHtab0 and HHtab2h unavailable, K is reset to -1, 
C     otherwise reading HHtab needs small modifications.
C
   13 IF(DR.EQ.-1)K=-1
      IF((DR.LE.-1).AND.(B.LT.0))K=-1
      IF(K.NE.1) GOTO 22
      IF(DR.GE.0)OPEN(1,FILE='HHtab2r')
      IF(DR.LT.-1)OPEN(1,FILE='HHtab2')
      IF((DR.EQ.-1).AND.(B.GE.0))OPEN(1,FILE='HHtab0')
      IF((DR.EQ.-1).AND.(B.LT.0))OPEN(1,FILE='HHtab2h')
      IF((DR.LT.0).AND.(DR.GT.-1))OPEN(1,FILE='HHtab0r')
   20 READ (1,*) I,N,M,X
      IF((I.EQ.19).OR.(M.GT.NM))GOTO 21
      IF((I.EQ.20).OR.(I.EQ.22))FF(N,M)=T**2*X
      IF((I.EQ.11).OR.(I.EQ.13).OR.(I.EQ.15))HH(N,M)=T**5*X
      IF((I.EQ.10).OR.(I.EQ.12).OR.(I.EQ.14).OR.
     >                (I.EQ.16).OR.(I.EQ.18))HH(N,M)=T**8*X/4
      HH(M,N)=HH(N,M)
      FF(M,N)=FF(N,M)
      GOTO 20
   21 CLOSE(1)
      IF((N.EQ.NM).OR.(DR.LT.0))GOTO 22
      PRINT *,"Reset nmax=",NM," to HHtab2r value ",N
      STOP
   22 IF(DR.NE.0)GOTO 43
C
C     Find the N values with l given, i.e. invert l(N) (needed if dr=0)
C
      I=1
      DO 31 M=3,LM
         NS(M)=I
         DO 30 I=NS(M),NM-1
            If(L(I).GT.M)GOTO 31
            NE(M)=I
   30    CONTINUE
   31 CONTINUE
      NE(LM)=NM
C
C     Read #radials for Y(n) (used as variational set), stop at "nmax p"
C     Check if # states stays below NN. Reset the last p(I) to match NN
C     If P(I)>pmax, reset P(I) to pmax=PM. Change te and tp with pruning.
C     When computing fsq and fn from psi (dr>=11), take P(I),b,te,tp from 
C     PROPER plist and read psi. Sets N to N(psi). If N(plist)<N(psi) or
C     b,te,tp do not match STOP. With IQ=1 we do not allow II=1.
C
   34 OPEN(1,FILE='plist')
      M=0
   35 READ (1,*) I,J
      P(I)=J
      IF(J.LT.0)P(I)=1
      IF(J.GT.PM)P(I)=PM
      IF((M+J).GT.NN)P(I)=NN-M
      M=M+P(I)
      IF(I.EQ.NM)GOTO 36
      GOTO 35
   36 IF(IP.EQ.0)READ (1,*) J,X,TE,TP
      CLOSE(1)
      IF(IQ.EQ.0)GOTO 41
      B=X
      OPEN(1,FILE='psi')
      READ (1,*) IF,N,X,I,DD,D1,D2
      IF((IF.EQ.2).AND.(J.GE.N).AND.(X.EQ.B)
     >           .AND.(I.EQ.II).AND.(D1.EQ.TE).AND.(D2.EQ.TP))GOTO 37
      PRINT *,"Wrong F,N,B,II,TE,TP for psi (B,TE,TP set by plist)"
      STOP
   37 IF(NN.EQ.N)GOTO 38
      CLOSE(1)
      NN=N
      GOTO 34
   38 IF(J.GT.N)PRINT *,"Warning: N(psi)<N(plist)"
      PRINT *,"Computing fsq and fn by reading psi"
      IF(II.EQ.0)WRITE (3,*) "(* F=2: N,b,E0,te,tp=",NN,B,DD,TE,TP," *)"
      IF(II.EQ.1)WRITE (3,*) "(* F=2: N,b,E1,te,tp=",NN,B,DD,TE,TP," *)"
   39 READ (1,*) N,X 
      Z(N+II*NN)=X
      IF(N.LT.NN)GOTO 39
      CLOSE(1)
C
C     Enumerate basisvectors (N,P) -> NF(N)+P, P=0,1,...,P(N)-1
C
   41 NF(1)=1
      DO 42 N=2,NM
         NF(N)=NF(N-1)+P(N-1)
   42 CONTINUE
C
C     When too big, reset NN to total number of available states!
C
   43 IF(M.LT.NN)NN=M
      IF(K.NE.1)MM=MIN(MM,NN)
      IF(IQ.EQ.1)GOTO 99
C
C     Print parameters, to screen and `tab-out' for a table.
C
      IF(KK.NE.1)GOTO 46
      IF(II.EQ.1)GOTO 44
      IF(DR.EQ.0)PRINT *,"F=2: N,b,bs,nb,te,tp=",
     >            NN,B,BS,NB,TE,TP,", {b,E0-dE0,E0,E1,0.}:"
      IF(DR.EQ.0)WRITE (3,*) "(* F=2: N,b,bs,nb,te,tp=",
     >            NN,B,BS,NB,TE,TP,", {b,E0-dE0,E0,E1,0.}: *)"
      IF(DR.LT.-1)PRINT *,"F=2: N,b,bs,nb=",
     >            NN,B,BS,NB,", {b,E0-dE0,E0,E1,0.}:"
      IF(DR.LT.-1)WRITE (3,*) "(* F=2: N,b,bs,nb=",
     >            NN,B,BS,NB,", {b,E0-dE0,E0,E1,0.}: *)"
      IF(DR.GT.0)PRINT *,"F=2: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {r,E0-dE0,E0,E1,Veff0}:"
      IF(DR.GT.0)WRITE (3,*) "(* F=2: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {r,E0-dE0,E0,E1,Veff0}: *)"
      IF((DR.LT.0).AND.(DR.GT.-1))PRINT *,"F=0: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {r,E0-dE0,E0,E1,Veff0}:"
      IF((DR.LT.0).AND.(DR.GT.-1))WRITE (3,*) "(* F=0: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {b,E0-dE0,E0,E1,Veff0}: *)"
      IF(DR.EQ.-1)PRINT *,"F=",INT(1.1-B/DABS(B)),": N,w,ws,nw=",
     >            NN,B,BS,NB,", {w,E0-dE0,E0,E1,0.}:"
      IF(DR.EQ.-1)WRITE (3,*) "(* F=",INT(1.1-B/DABS(B)),": N,w,ws,nw=",
     >            NN,B,BS,NB,", {w,E0-dE0,E0,E1,0.}: *)"
      GOTO 45
   44 IF(DR.EQ.0)PRINT *,"F=2: N,b,bs,nb,te,tp=",
     >            NN,B,BS,NB,TE,TP,", {b,E1-dE1,E1,E2,0.}:"
      IF(DR.EQ.0)WRITE (3,*) "(* F=2: N,b,bs,nb,te,tp=",
     >            NN,B,BS,NB,TE,TP,", {b,E1-dE1,E1,E2,0.}: *)"
      IF(DR.LT.-1)PRINT *,"F=2: N,b,bs,nb=",
     >            NN,B,BS,NB,", {b,E1-dE1,E1,E2,0.}:"
      IF(DR.LT.-1)WRITE (3,*) "(* F=2: N,b,bs,nb=",
     >            NN,B,BS,NB,", {b,E1-dE1,E1,E2,0.}: *)"
      IF(DR.GT.0)PRINT *,"F=2: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {r,E1-dE1,E1,E2,Veff1}:"
      IF(DR.GT.0)WRITE (3,*) "(* F=2: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {r,E1-dE1,E1,E2,Veff1}: *)"
      IF((DR.LT.0).AND.(DR.GT.-1))PRINT *,"F=0: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {r,E1-dE1,E1,E2,Veff1}:"
      IF((DR.LT.0).AND.(DR.GT.-1))WRITE (3,*) "(* F=0: N,r,rs,nr,dr=",
     >            NN,B,BS,NB,DR,", {b,E1-dE1,E1,E2,Veff1}: *)"
      IF(DR.EQ.-1)PRINT *,"F=",INT(1.1-B/DABS(B)),": N,w,ws,nw=",
     >            NN,B,BS,NB,", {w,E1-dE1,E1,E2,0.}:"
      IF(DR.EQ.-1)WRITE (3,*) "(* F=",INT(1.1-B/DABS(B)),": N,w,ws,nw=",
     >            NN,B,BS,NB,", {w,E1-dE1,E1,E2,0.}: *)"
   45 WRITE (3,*) "Es={"
C
C     While reading radial matrix elements (dr=0), build H (in WORK) for 
C     each b, diagonalize, build H*H, reset WORK=AM=0 and move to b+bs.
C     Retrun to this point when making a table for dr=0.
C
   46 DO 47 I=1,NN*NN
         WORK(I)=0.
   47 CONTINUE
C
C     For dr=0 read Rtab2 to build H, else use Htab values.
C
      IF(DR.EQ.0)GOTO 49
      DO 48 IM=1,NN
         WORK((IM-1)*NN+IM)=H(IM,IM)
         DO 48 IN=IM+1,NN
            WORK((IN-1)*NN+IM)=H(IN,IM)
            WORK((IM-1)*NN+IN)=H(IM,IN)
   48 CONTINUE
      GOTO 55
C
C     I labels the power of B to be multiplied with for each matrix element,
C     and identifies its type: I=4 -> Hv, I=1 -> Hf (no mixing so both in H),
C     I=-2 -> Z(l(N),P)^2/2, the diagonal kinetic term stored in HD.
C
   49 OPEN(1,FILE='Rtab2')
   50 READ (1,*) I,NP,MP,NL,ML,X
      IF(I.EQ.0)GOTO 53
      DO 52 IM=NS(ML),NE(ML)
         IF(MP.GE.P(IM))GOTO 52
         JM=NF(IM)+MP
         IF(I.EQ.-2)HD(JM)=X*(B**I)
         IF(I.EQ.-2)GOTO 52
         DO 51 IN=NS(NL),NE(NL)
            IF(NP.GE.P(IN))GOTO 51
            JN=NF(IN)+NP
C
C           When NL=ML (equal momentum l), we listed in Rtab2 for NP<=MP.
C           If IN=IM (same Y[n]) this is ok, else we also need to exchange
C           NP and MP, when unequal. Easiest is just to symmetrize.
C
            WORK((JN-1)*NN+JM)=H(IN,IM)*X*(B**I)
            WORK((JM-1)*NN+JN)=H(IN,IM)*X*(B**I)
   51    CONTINUE         
   52 CONTINUE         
      GOTO 50
   53 CLOSE(1)
      IF(NL.EQ.LM)GOTO 54
      PRINT *,"Reset lmax=",LM," to Rtab2 value ",NL
      STOP
   54 IF(PM.LE.(NP+1))GOTO 55
      PRINT *,"Warning: pmax=",PM," was reset to Rtab2 value ",NP+1
      PM=NP+1
   55 IER=0
C
C     First find the number of codiagonals. 
C
      NC=0
      DO 56 I=1,NN
         DO 56 J=I+NC,NN
            IF(WORK((I-1)*NN+J).NE.0)NC=J-I
   56 CONTINUE
C
C     Warn if NN>NNM or NC>CM and set up the matrix AM for H
C
      IF((NN.GT.NNM).OR.(NC.GT.NCM))PRINT *,"Warning: NC=",NC,
     >                         " or N=",NN," might be too big"
      DO 57 I=0,NC
        DO 57 J=1,NN
	   IR=I+J-NC 
           IF(IR.GT.0)AM(I*NN+J)=WORK((IR-1)*NN+J)
           IF(IR.EQ.J)AM(I*NN+J)=AM(I*NN+J)+HD(J)
   57 CONTINUE
      CALL EIGBS (AM,NN,NN,IJOB,NC,MM,D,Z,NN,WORK,IER)
C
C     Writing the eigenvector to psi if NB=1 (when dr=0, keep with plist!)
C
      IF((NB.NE.1).OR.(IJOB.EQ.-1))GOTO 67
      PRINT *,"Writing psi"
      IF(DR.EQ.0)PRINT *,"To be kept with plist!"
      OPEN(2,FILE='psi')
      IF(DR.EQ.0) WRITE (2,*) 2,NN,B,II,D(1+II),TE,TP
      IF(DR.GT.0) WRITE (2,*) 2,NN,B,II,D(1+II),0.,0.
      IF(DR.LT.-1)WRITE (2,*) 2,NN,B,II,D(1+II),0.,0.
      IF((DR.EQ.-1).AND.(B.GE.0))WRITE (2,*) 0,NN,B,II,D(1+II),0.,0.
      IF((DR.EQ.-1).AND.(B.LT.0))WRITE (2,*) 2,NN,B,II,D(1+II),0.,0.
      IF((DR.LT.0).AND.(DR.GT.-1))WRITE (2,*) 0,NN,B,II,D(1+II),0.,0.
      DO 66,I=1,NN
         WRITE (2,*) I,Z(I+II*NN)
   66 CONTINUE
      CLOSE(2)
   67 IF(K.NE.1) GOTO 77
C
C     Temple's inequality -- for dr=0 read RRtab2, else use HHtab values
C     
      T=0.
      IF(DR.EQ.0)GOTO 69
      DO 68 IM=1,NN
         T=T+Z(IM+II*NN)*(HH(IM,IM)+FF(IM,IM))*Z(IM+II*NN)
         DO 68 IN=IM+1,NN
            T=T+2*Z(IN+II*NN)*(HH(IN,IM)+FF(IN,IM))*Z(IM+II*NN)
   68 CONTINUE
      GOTO 75
C
C     Do Temple's inequality while reading RRtab2. I labels the power 
C     of B multiplied with for each matrix element, and identifies its 
C     type: I=8 -> Hvv, I=5 -> Hvf (HH, no mixing), I=2 -> Hff (FF).
C
   69 OPEN(1,FILE='RRtab2')
   70 READ (1,*) I,NP,MP,NL,ML,X
      IF(I.EQ.0)GOTO 73
      DO 72 IM=NS(ML),NE(ML)
         IF(MP.GE.P(IM))GOTO 72
         JM=NF(IM)+MP
         DO 71 IN=NS(NL),NE(NL)
            IF(NP.GE.P(IN))GOTO 71
            IF(I.EQ.2)DD=FF(IN,IM)*X*(B**I)
            IF(I.NE.2)DD=HH(IN,IM)*X*(B**I)
            JN=NF(IN)+NP
C
C           When NL=ML (equal momentum l), we listed in RRtab2 for NP<=MP.
C           If IN=IM (same Y[n]) this is ok, else we also need to exchange
C           NP and MP, when unequal. Summing over all IN,IM with NL=ML
C           double counts if NP=MP, so insist (NP!=MP)||(IN<IM)
C
            IF((NP.NE.MP).OR.(IN.LE.IM))T=T+2*Z(JN+II*NN)*DD*Z(JM+II*NN)
            IF(JN.EQ.JM)T=T-Z(JN+II*NN)*DD*Z(JM+II*NN)
   71    CONTINUE         
   72 CONTINUE         
      GOTO 70
   73 CLOSE(1)
      IF(NL.EQ.LM)GOTO 74
      PRINT *,"Reset lmax=",LM," to RRtab2 value ",NL
      STOP 
   74 IF(PM.LE.(NP+1))GOTO 75
      PRINT *,"Warning: pmax=",PM," was reset to RRtab2 value ",NP+1
      PM=NP+1
   75 DO 76 J=1,NN
         T=T-((D(1+II)-HD(J))**2)*(Z(J+II*NN)**2)
   76 CONTINUE
   77 D1=D(1+II)
      D2=D(2+II)
      DD=D1-T/(D2-D1)
C
C     Compute derivative of the wavefunction, by making two more passes with
C     r-dr and r+dr. Requires re-reading Htab0r (-1<dr<0) or Htab2r (dr>0).
C     Htab2r: NC=INT(4Sqrt(NN)-3), Htab0r: NC=INT(1.3Sqrt(NN)), but NO need
C     to compute the number of codiagonals again! Here we use HD to store 
C     the wavefunction after the first pass. WORK now contains all of H.
C
      T=0.
      IF((DR.EQ.0).OR.(DR.LE.-1))GOTO 90
      IJOB=-2
      DO 80 I=1,NN*NN
         WORK(I)=0.
   80 CONTINUE
      IF(DR.GE.0)OPEN(1,FILE='Htab2r')
      IF(DR.LT.0)OPEN(1,FILE='Htab0r')
   81 READ (1,*) I,N,M,X
      IF(I.EQ.9)GOTO 82
      IF(I.EQ.1)WORK((N-1)*NN+M)=(B-DR)*X
      IF((I.EQ.0).OR.(I.EQ.2).OR.(I.EQ.4))WORK((N-1)*NN+M)=(B-DR)**4*X/2
      IF(I.EQ.3)WORK((N-1)*NN+M)=WORK((N-1)*NN+M)+X/(B-DR)**2/2
      WORK((M-1)*NN+N)=WORK((N-1)*NN+M)
      GOTO 81
   82 CLOSE(1)
      DO 83 I=0,NC
        DO 83 J=1,NN
           IR=I+J-NC
           IF(IR.GT.0)AM(I*NN+J)=WORK((IR-1)*NN+J)
           IF(IR.EQ.J)AM(I*NN+J)=AM(I*NN+J)
   83 CONTINUE
      CALL EIGBS (AM,NN,NN,IJOB,NC,MM,D,Z,NN,WORK,IER)
      DO 84 N=1,NN
         HD(N)=Z(N+II*NN)
   84 CONTINUE
      DO 85 I=1,NN*NN
         WORK(I)=0.
   85 CONTINUE
      IF(DR.GE.0)OPEN(1,FILE='Htab2r')
      IF(DR.LT.0)OPEN(1,FILE='Htab0r')
   86 READ (1,*) I,N,M,X
      IF(I.EQ.9)GOTO 87
      IF(I.EQ.1)WORK((N-1)*NN+M)=(B+DR)*X
      IF((I.EQ.0).OR.(I.EQ.2).OR.(I.EQ.4))WORK((N-1)*NN+M)=(B+DR)**4*X/2
      IF(I.EQ.3)WORK((N-1)*NN+M)=WORK((N-1)*NN+M)+X/(B+DR)**2/2
      WORK((M-1)*NN+N)=WORK((N-1)*NN+M)
      GOTO 86
   87 CLOSE(1)
      DO 88 I=0,NC
        DO 88 J=1,NN
           IR=I+J-NC
           IF(IR.GT.0)AM(I*NN+J)=WORK((IR-1)*NN+J)
           IF(IR.EQ.J)AM(I*NN+J)=AM(I*NN+J)
   88 CONTINUE
      CALL EIGBS (AM,NN,NN,IJOB,NC,MM,D,Z,NN,WORK,IER)
C
C     Z(1+II*NN) assumed non-zero. Correct for overall sign ambiguity.
C
      X=1.
      IF((HD(1)/Z(1+II*NN)).LT.0)X=-1.
      DO 89 N=1,NN
         T=T+(Z(N+II*NN)-X*HD(N))**2
   89 CONTINUE
C
C     (dX/dr)^2=T/DR/DR/4+O(DR^2), used to compute Veff.
C
      T=D1+T/DR/DR/8+6/B/B
C
C     Print the results -- DD=0. without Temple, T=0. without Veff.
C
   90 IF(K.NE.1)DD=0
      IF(KK.LT.NB)WRITE (3,*) "    {",B,",",DD,",",D1,",",D2,",",T,"},"
      IF(KK.GE.NB)WRITE (3,*) "    {",B,",",DD,",",D1,",",D2,",",T,"}}"
      PRINT *,"{",B,",",DD,",",D1,",",D2,",",T,"}"
C
C     If bs=0 and nb>2/3 we use nb to set the number of eigenvalues to 
C     be computed and write a table. Switched off when we compute H^2.
C
      IF(MM.LE.2+II) GOTO 92
      WRITE (3,*) "Etab={"
      DO 91 I=1,MM-1
          PRINT *,"{",I,",",D(I),"},"
          WRITE (3,*) "{",I,",",D(I),"},"
   91 CONTINUE
      PRINT *,"{",MM,",",D(MM),"}"
      WRITE (3,*) "{",MM,",",D(MM),"}}"
C
C     Making a table (nb>1), for dr=0 it suffices to re-read (R)Rtab2.
C     Otherwize we need to re-read everything so as to reset b (w or r).
C
   92 B=B+BS
      KK=KK+1
      IF((DR.EQ.0).AND.(KK.LE.NB))GOTO 46
      IF((DR.NE.0).AND.(KK.LE.NB))GOTO 6
      IF((DR.NE.0).OR.(IJOB.EQ.-1).OR.(NB.GT.1))STOP
      B=B-BS
      IF(IP.EQ.0)GOTO 99
C
C     Prune when dr=0 and nb<0 -> IP=1, nb=1. We prune basis vectors with
C     |Z[n]|<T. Reset P(I) accordingly, add 1 otherwise to improve accuracy,
C     writing to plist. In plist the last line gives the sum of the P(I), 
C     and the pruning level and minimum relative error (TP and TE) in use. 
C
      IF(DABS(1.0-DD/D1).LT.TE)STOP
      PRINT *,"Replacing old plist with new one"
      OPEN(2,FILE='plist')
      KK=0
      DO 94 I=1,NM 
         J=0
         DO 93 N=0,P(I)-1
            IF(DABS(Z(NF(I)+N)).GE.TP)J=J+1
   93    CONTINUE
C
C        If no eigenvector survives the pruning for given I, keep it to zero.
C        This requires some care -- replace after done 0's by 1's and see if
C        these are set to 0 again (one iteration, since KK decreases).
C
         IF(J.EQ.0)J=-1
         WRITE (2,*) I,MIN(J+1,PM)
         KK=KK+MIN(J+1,PM)
   94 CONTINUE
      WRITE (2,*) KK,B,TE,TP
      CLOSE(2)
C
C     We can continue this process until N no longer changes (pmax saturated).
C     The first pruning tends to remove many states, so always restart (IP=1).
C     Last pruning removes as many states as it adds in plist. Set IP=0 and 
C     restart. This properly pairs plist to psi and also makes fsq and fn. 
C
      IF((KK.LE.NN).AND.(IP.GT.1))IP=-1
      IP=IP+1
      NN=KK
      KK=1
      GOTO 6
C
C     If dr=0, nb=1, N>0, compute fn*fn at r=I*b/IT for I=1 to IT using 
C     Bjtab (made with Bjtab.ma), which contains tabulated values of 
C     Bj[l,Z[l,p]I/IT]/(Sqrt[Nr[Z[l,p],l]]Bj[l,Z[l,p]]) with Z[l,p-1]
C     the p-th root of the derivative of the spherical Bessel function 
C     Bj[l,z]. In Bjtab we first increase p, then l and finally I. We
C     diagonalize H at each r, store chi0 in HD (Z already used for Psi)  
C     and compute f0=<Psi|chi0>. NOTE: NO factor 4pi is extracted and
C     \int_0^b dr\sum_n(r*f_n(r))^2=1 (instead of 1/4pi). 
C
   99 T=0.
      D1=0.
      IR=0
      PRINT *,"{r,fsq,f0}="
      WRITE (3,*) "fsqfn={"
      OPEN(1,FILE='Bjtab')
      READ (1,*) IT,LT,PT,X
      IF((LT.EQ.LM).AND.(PT.GE.PM))GOTO 100
      PRINT *,"Range of l in Bjtab,",LT,", does not match lmax=",LM 
      PRINT *,"or range of p,",PT,", does not extend to pmax=",PM 
      STOP
  100 IR=IR+1
      IF(NB.EQ.0)GOTO 105
      DO 101 I=1,NM*NM
         WORK(I)=0.
  101 CONTINUE
      OPEN(2,FILE='Htab2r')
  102 READ (2,*) I,N,M,X
      IF(I.EQ.9)GOTO 103
      IF(I.EQ.1)WORK((N-1)*NM+M)=IR*B*X/IT
      IF((I.EQ.0).OR.(I.EQ.2).OR.(I.EQ.4))
     >                           WORK((N-1)*NM+M)=(IR*B/IT)**4*X/2
      IF(I.EQ.3)WORK((N-1)*NM+M)=WORK((N-1)*NM+M)+X/(IR*B/IT)**2/2
      WORK((M-1)*NM+N)=WORK((N-1)*NM+M)
      GOTO 102
  103 CLOSE(2)
      NC=INT(4*SQRT(NM*1.)-3)
      DO 104 I=0,NC
        DO 104 J=1,NM
           N=I+J-NC
           IF(N.GT.0)AM(I*NM+J)=WORK((N-1)*NM+J)
           IF(N.EQ.J)AM(I*NM+J)=AM(I*NM+J)
  104 CONTINUE
      CALL EIGBS (AM,NM,NM,-2,NC,1,D,HD,NM,WORK,IER)
  105 READ (1,*) IR,LN,PN,X
      WORK(PN+1)=X
      IF(PN.NE.(PT-1))GOTO 105
      DO 107 N=NS(LN),NE(LN)
         X=0.
         DO 106 I=0,P(N)-1
            X=X+Z(NF(N)+I+II*NN)*WORK(I+1)
  106    CONTINUE
         T=T+X*X
         IF(NB.EQ.0)GOTO 107
         IF(HD(1).LT.0)D1=D1-X*HD(N)
         IF(HD(1).GE.0)D1=D1+X*HD(N)
  107 CONTINUE
      IF(LN.LT.LM)GOTO 105
      T=T/B**3
      D1=D1/B/DSQRT(B)
      IF(Z(1+II*NN).LT.0)D1=-D1
      PRINT *,"{",IR*B/IT,",",T,",",D1,"}"
      IF(IR.LT.IT)WRITE (3,*)
     >            "       {",IR*B/IT,",",T,",",D1,"},"
      IF(IR.EQ.IT)WRITE (3,*) 
     >            "       {",IR*B/IT,",",T,",",D1,"}}"
      D1=0.
      T=0.
      IF(IR.LT.IT)GOTO 100
      CLOSE(1)
      STOP
      END