Lanczos Diagonalization for S=1 Heisenberg Spin Ring
@

C123456789012345678901234567890123456789012345678901234567890123456789012
C
C     +-----------------------------------------------------+
C     |                                                     |
C     |   Diagonalization of S = 1 Heisenberg Model         |
C     |                                                     |
C     |   H = S*S + a(S*S)^2, Generalized S=1 chain...      |
C     |                                                     |
C     |   Bisection Search Method is used.                  |
C     |   Translational Invariance and Parity are used.     |
C     |   Periodic Boundary Condition is assumed.           |
C     |                                                     |
C     | * This program assumes "Long Integer Mode"; some *  |
C     | * integers should be 64-bit integer.             *  |
C     |                                                     |
C     +-----------------------------------------------------+
C
C    Site#       Config. #       Reduced Config.
C                (Sz=0)          (*About*)
C      3:               7                3
C      4:              19                7
C      5:              51               11
C      6:             141               23*  Dump mode
C      7:             393               56*
C      8:           1,107              138*
C      9:           3,139              348*
C     10:           8,953              895*
C     11:          25,653            2,332*
C     12:          73,789            6,149*  
C     13:         212,941           16,380*
C     14:         616,227           44,016*
C     15:       1,787,607          119,173*  Debug Mode
C     16:       5,196,627          324,789*  Standard
C     17:      15,134,931          890,290*
C     18:      44,152,809        2,452,933*  (3^ 9 =  19683)
C     19:     128,996,853        6,789,308*
C     20:     377,379,369       18,868,968*  (3^10 =  59049)
C     21:   1,105,350,729       52,635,749*  (3^11 = 177147)
C
C     +---------------------------------------------------+
C     |   Limits:    NS =< 40.                            |
C     |   Patch : if NS  > 20, modify 59049 ( = 3^10).    |
C     |           if NS  > 20, modify  9953 in Subroutine |
C     |           "CRESTA" (>Largest Subspace for NS=10)  |
C     |           if NS  > 15, modify 122,222 and 244,444 |
C     |                            (=Vector Dimensions)   |
C     |           if NS  > 32  modify Subroutine "CRENUM" |
C     +---------------------------------------------------+
C
      INTEGER      LV, NS,  NST,  BH(0:59049,0:2),   ST(0:122222)
      INTEGER      NK, NX,   TZ,  BL(0:59049,0:2),  EXS(0:999)
      REAL*8        A, CS(0:40),  SN(0:40),         NRM(0:999)
      REAL*8  VEC(0:244444,0:3), DIA(0:122222)
C
C     +---------------------------------------------------------------+
C     | Debug Level: 0(None) 1(Normal) 2(Detail) 3(Dump) 4(Core Dump) |
C     +---------------------------------------------------------------+
         LV   =  1
C     +------------------------+
C     |   System Size "NS"     |
C     |   NS Should be even!!  |
C     |   Total  Sz   "Tz"     |
C     +------------------------+
            NS = 8
      DO 95  J = NS-2, 0, -1
            TZ = J
C
C     *** Parity is not used yet ***
C     +---------------------------------+
C     | Coefficient of the (S*S)^2 term |
C     | A=1: Permutation Model          |
C     | A=?: AKLT Model                 |
C     +---------------------------------+
            A = 0.0D0
C
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) '*** Parameters:'
      WRITE(6,*) '*** Ns, Total Sz, Alpha =', NS, TZ, A
C
C     +----------------------------+
C     | Total Momentum = pi*(I/NS) |
C     +----------------------------+
      DO 90 I = 0, NS/2
C-DB- DO 90 I = 0, 0
           NK = I
      WRITE(6,*) 'Total Momentum = ', 2*NK, '/', NS, '*pi'
      CALL CRENUM( NS,                    BH,BL               )
      CALL CRESTA( NS,LV,TZ,   NST,ST,    BH,BL               )
C
      IF( NST .EQ. 1 ) THEN
          WRITE(6,*) 'You Can Get Analytic Formula'
          GOTO 80
        END IF
C
      CALL ORTSET( NS,LV,NK,              SN,CS               )
      CALL NORMAL( NS,LV,   NX,NST,ST,EXS,   CS,NRM,      VEC )
      CALL CREDIA( NS,LV,      NST,ST,              A,DIA     )
      CALL LANZOS( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC )
      CALL VECDMP( NS,LV,NK,   NST,ST,                    VEC )
      CALL SZDIST( NS,   NK,   NST,ST,    BH,BL,          VEC )
   80 CONTINUE
   90 CONTINUE
   95 CONTINUE
C
      STOP
      END
CCC
CC
C
      SUBROUTINE CRENUM(   NS, BH,              BL )
      INTEGER  M, NS, T(0:20), BH(0:59049,0:2), BL(0:59049,0:2)
      INTEGER         B0(0:2), B1(0:2),         B2(0:2)
      INTEGER  A0, A1, A2, A3, A4, A5, A6, A7,  A8,  A9
      INTEGER  AA, AB, AC, AD, AE, AF
C
C            Bh                  Bl
C     [NS-1 .... NS/2]  [NS/2-1 .... 0]
C
            M = NS / 2
      DO 10 I =  0, 20
      IF( I .LT. M ) THEN
          T(I)= 1
        ELSE
          T(I)= 0
        END IF
   10 CONTINUE
C
      DO 20 I = 0, 2
        B0(I) = MOD(I+1,2) - I/2
        B1(I) = MOD(I  ,2)
        B2(I) =              I/2
   20 CONTINUE
C
      DO 30 I = 0, 3**M - 1
           A0 = MOD(I      , 3)
           A1 = MOD(I/3    , 3)
           A2 = MOD(I/3** 2, 3)
           A3 = MOD(I/3** 3, 3)
           A4 = MOD(I/3** 4, 3)
           A5 = MOD(I/3** 5, 3)
           A6 = MOD(I/3** 6, 3)
           A7 = MOD(I/3** 7, 3)
           A8 = MOD(I/3** 8, 3)
           A9 = MOD(I/3** 9, 3)
           AA = MOD(I/3**10, 3)
           AB = MOD(I/3**11, 3)
           AC = MOD(I/3**12, 3)
           AD = MOD(I/3**13, 3)
           AE = MOD(I/3**14, 3)
           AF = MOD(I/3**15, 3)
      BH(I,0) = T( 0)*B0(A0) +T( 1)*B0(A1) +T( 2)*B0(A2) +T( 3)*B0(A3)
     &        + T( 4)*B0(A4) +T( 5)*B0(A5) +T( 6)*B0(A6) +T( 7)*B0(A7)
     &        + T( 8)*B0(A8) +T( 9)*B0(A9) +T(10)*B0(AA) +T(11)*B0(AB)
     &        + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF)
      BH(I,1) = T( 0)*B1(A0) +T( 1)*B1(A1) +T( 2)*B1(A2) +T( 3)*B1(A3)
     &        + T( 4)*B1(A4) +T( 5)*B1(A5) +T( 6)*B1(A6) +T( 7)*B1(A7)
     &        + T( 8)*B1(A8) +T( 9)*B1(A9) +T(10)*B1(AA) +T(11)*B1(AB)
     &        + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF)
      BH(I,2) = T( 0)*B2(A0) +T( 1)*B2(A1) +T( 2)*B2(A2) +T( 3)*B2(A3)
     &        + T( 4)*B2(A4) +T( 5)*B2(A5) +T( 6)*B2(A6) +T( 7)*B2(A7)
     &        + T( 8)*B2(A8) +T( 9)*B2(A9) +T(10)*B2(AA) +T(11)*B2(AB)
     &        + T(12)*B2(AC) +T(13)*B2(AD) +T(14)*B2(AE) +T(15)*B2(AF)
   30 CONTINUE
C
            M = NS - M
      DO 40 I =  0, 20
      IF( I .LT. M ) THEN
          T(I)= 1
        ELSE
          T(I)= 0
        END IF
   40 CONTINUE
C
      DO 50 I = 0, 3**M - 1
           A0 = MOD(I      , 3)
           A1 = MOD(I/3    , 3)
           A2 = MOD(I/3** 2, 3)
           A3 = MOD(I/3** 3, 3)
           A4 = MOD(I/3** 4, 3)
           A5 = MOD(I/3** 5, 3)
           A6 = MOD(I/3** 6, 3)
           A7 = MOD(I/3** 7, 3)
           A8 = MOD(I/3** 8, 3)
           A9 = MOD(I/3** 9, 3)
           AA = MOD(I/3**10, 3)
           AB = MOD(I/3**11, 3)
           AC = MOD(I/3**12, 3)
           AD = MOD(I/3**13, 3)
           AE = MOD(I/3**14, 3)
           AF = MOD(I/3**15, 3)
      BL(I,0) = T( 0)*B0(A0) +T( 1)*B0(A1) +T( 2)*B0(A2) +T( 3)*B0(A3)
     &        + T( 4)*B0(A4) +T( 5)*B0(A5) +T( 6)*B0(A6) +T( 7)*B0(A7)
     &        + T( 8)*B0(A8) +T( 9)*B0(A9) +T(10)*B0(AA) +T(11)*B0(AB)
     &        + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF)
      BL(I,1) = T( 0)*B1(A0) +T( 1)*B1(A1) +T( 2)*B1(A2) +T( 3)*B1(A3)
     &        + T( 4)*B1(A4) +T( 5)*B1(A5) +T( 6)*B1(A6) +T( 7)*B1(A7)
     &        + T( 8)*B1(A8) +T( 9)*B1(A9) +T(10)*B1(AA) +T(11)*B1(AB)
     &        + T(12)*B1(AC) +T(13)*B1(AD) +T(14)*B1(AE) +T(15)*B1(AF)
      BL(I,2) = T( 0)*B2(A0) +T( 1)*B2(A1) +T( 2)*B2(A2) +T( 3)*B2(A3)
     &        + T( 4)*B2(A4) +T( 5)*B2(A5) +T( 6)*B2(A6) +T( 7)*B2(A7)
     &        + T( 8)*B2(A8) +T( 9)*B2(A9) +T(10)*B2(AA) +T(11)*B2(AB)
     &        + T(12)*B2(AC) +T(13)*B2(AD) +T(14)*B2(AE) +T(15)*B2(AF)
   50 CONTINUE
C
      RETURN
      END 
CCC
CC
C
      SUBROUTINE CRESTA( NS, LV, TZ, NST,      ST, BH, BL)
      INTEGER  NL,LZ, M, NS, LV, TZ, NST,      ST(0:122222)
      INTEGER   BH(0:59049,0:2), PT(0:131441), FG(0:131441)
      INTEGER   BL(0:59049,0:2), TP(0:131441),  NLS(-10:10)
      INTEGER                             LS(0:9953,-10:10)
C
            M = NS / 2
C                                   --- Make Table ---
      DO 20 I = -M, M
           NL =  0
      DO 10 J =  0, 3**M -1
           LZ = BL(J,2) - BL(J,0)
      IF(  LZ .EQ. I ) THEN
        LS(NL,I) = J
           NL    = NL + 1
        END IF
   10 CONTINUE
        NLS(I)= NL
   20 CONTINUE
C
          NST = 0
      DO 80 I = 0, 3**(NS-M) -1
           LZ = TZ - BH(I,2) + BH(I,0)
C                                   --- Select Polalization ---
      IF( ABS(LZ) .LE. M ) THEN
          DO 30 J = 0, NLS(LZ) -1
             PT(J)= I*3**M + LS(J,LZ)
             FG(J)=   0
   30     CONTINUE
C                                   --- Copy States ---
          DO 40 K = 0, NLS(LZ) -1
             TP(K)= PT(K)
   40     CONTINUE
C                                   --- Mark States ---
          DO 60 K = 1, NS-1
          DO 50 J = 0, NLS(LZ) -1
             TP(J)= MOD(TP(J)*3,3**NS) + TP(J)/3**(NS-1)
          IF( PT(J) .GT. TP(J) ) FG(J) = 1
   50     CONTINUE
   60     CONTINUE
C                                   --- Select States ---
          DO 70 J = 0, NLS(LZ) -1
          IF( FG(J) .EQ. 0 ) THEN
              ST(NST) = PT(J)
                 NST  = NST + 1
            END IF
   70     CONTINUE
        END IF
   80 CONTINUE
C
      IF( LV .GE. 1 ) WRITE(6,*) 'Total # of basis=', NST
      IF( LV .GE. 3 ) THEN
          DO 90 I = 0, NST -1
          CALL STDISP( NS, ST(I), 0.0D0 )
   90     CONTINUE
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE STDISP( NS, ST, F )
      INTEGER        TP, NS, ST, LW(0:40)
      REAL*8                     F
           TP = ST
      DO 10 I = NS -1, 0, -1
         LW(I)= MOD(TP,3)
           TP =     TP/3
   10 CONTINUE
      WRITE(6,100)ST, F, (LW(I),I=0,NS-1)
  100 FORMAT(1H ,I14,' ---', F20.10,' ', 40I1)
      RETURN
      END
CCC
CC
C
      SUBROUTINE ORTSET( NS, LV, NK, SN      , CS     )
      INTEGER            NS, LV, NK
      REAL*8                         SN(0:40), CS(0:40)
      REAL*8            PI
C
      IF( LV .GE. 2 ) WRITE(6,*) '--- Subroutine ORTSET ---'
      PI = 4.0D0*ATAN(1.0D0)
C
      DO 10 I = 0, NS-1
         SN(I)= SIN( 2.0D0*I*NK*PI/NS )
         CS(I)= COS( 2.0D0*I*NK*PI/NS )
      IF( LV .GE. 2 ) WRITE(6,100) SN(I), CS(I)
  100 FORMAT(1H ,' SIN=', F10.5, '   COS=', F10.5 )
   10 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE NORMAL( NS, LV, NX, NST, ST, EXS, CS, NRM, VEC )
      INTEGER        RT, NS, LV, NX, NST, ST(0:122222),EXS(0:999)
      INTEGER                               TP(0:122222)
      REAL*8  F, CS(0:40), NRM(0:999), VEC(0:244444,0:3)
C
      IF( LV .GE. 2 ) WRITE(6,*) '--- Subroutine NORMAL ---'
C
      DO 10 I = 0, NST -1
      VEC(I,0)= 0.0D0
       TP(I  )= ST(I)
   10 CONTINUE
C                             --- Mark Self Similar States ---
      DO 30 I = 1, NS  -1
      DO 20 J = 0, NST -1
         TP(J)= MOD(TP(J)*3,3**NS) + TP(J)/3**(NS-1)
      IF( ST(J) .EQ. TP(J) ) THEN
          VEC(J,0) = 1.0D0
        END IF
   20 CONTINUE
   30 CONTINUE
C                           --- Pick up Self Similar States ---
           NX = 0
      DO 40 I = 0, NST -1
      IF( VEC(I,0) .GT. 1.0E-8 ) THEN
          EXS(NX)= I
              NX = NX + 1
        END IF
   40 CONTINUE
C                           --- Normalize ---
      DO 60 I = 0, NX -1
           RT = ST(EXS(I))
            F = 1.0D0
      DO 50 J = 1, NS -1
           RT = RT/3**(NS-1) + MOD(RT*3,3**NS)
      IF(  RT .EQ. ST(EXS(I))       )  F = F + CS(J)
   50 CONTINUE
      IF( F .GT. 1.0E-8 ) THEN
          NRM(I) = 1.0D0 / SQRT( F )
        ELSE
          NRM(I) = 0.0D0
        END IF
   60 CONTINUE
C
      IF( NX .GT. 999) THEN
          WRITE(6,*)'Degenerate Basis # > ',999
          do 99 i = 0, 100
          call STDISP(NS,ST(EXS(i)),0.0D0)
   99     continue
          WRITE(6,*)'...Increase 999 in the source list'
          STOP
        END IF
C
      IF( LV .GE. 1 ) WRITE(6,200) NX
  200 FORMAT(1H ,' # of Auto Correlated Basis=', I10 )
C
      IF( LV .GE. 2 ) THEN
          DO 70 I = 0, NX -1
          CALL STDISP( NS, ST(EXS(I)), NRM(I) )
   70     CONTINUE
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE CREDIA( NS, LV, NST, ST, A, DIA )
      INTEGER            NS, LV, NST, ST(0:122222)
      REAL*8            AA(0:4),   A,DIA(0:122222)
C
      IF( LV .GE. 2 ) WRITE(6,*)'SUBROUTINE CREDIA'
C
      AA(0)= A
      AA(1)= A
      AA(2)= A*2
      AA(3)= A
      AA(4)= A
C
C     --- Boundary Part ---
      DO 10 J = 0, NST -1
        DIA(J)=   ( MOD(ST(J),3)-1 )*( ST(J)/3**(NS-1)-1 )
&             + AA( MOD(ST(J),3)    +  ST(J)/3**(NS-1)   )
   10 CONTINUE
C
      DO 30 I = 0, NS  -2
      DO 20 J = 0, NST -1
        DIA(J)= DIA(J)
     &     +   ( MOD(ST(J)/3**I,3)-1 )*( MOD(ST(J)/3**(I+1),3)-1 )
     &     + AA( MOD(ST(J)/3**I,3)    +  MOD(ST(J)/3**(I+1),3)   )
   20 CONTINUE
   30 CONTINUE
C
      IF( LV .GE. 3 ) THEN
          DO 40 I = 0, NST-1
          WRITE(6,100) I, DIA(I)
   40     CONTINUE
  100     FORMAT(1H , I4, F15.7 )
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE INPRO( NV, ID, IS, PROD, VEC)
      INTEGER           NV, ID, IS
      REAL*8                        PROD, VEC(0:244444,0:3)
         PROD = 0.0D0
      DO 10 I = 0, NV-1
         PROD = PROD + VEC(I,ID) * VEC(I,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE ADDER( NV, ID, IS1, IS2, A, B, VEC)
      INTEGER           NV, ID, IS1, IS2
      REAL*8                              A, B, VEC(0:244444,0:3)
      DO 10  I   = 0, NV-1
       VEC(I,ID) = A*VEC(I,IS1) + B*VEC(I,IS2)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE DIVER( NV, ID, IS, D,VEC)
      INTEGER           NV, ID, IS
      REAL*8                        D,VEC(0:244444,0:3)
      DO 10 I   = 0, NV-1
      VEC(I,ID) = VEC(I,IS) / D
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE VEC0( NV, VEC)
      INTEGER          NV
      REAL*8               VEC(0:244444,0:3)
      DO 20 I = 0, 3
      DO 10 J = 0, NV -1
      VEC(J,I)= 0.0D0
   10 CONTINUE
   20 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE VMOVE( NV, ID, IS, VEC)
      INTEGER           NV, ID, IS
      REAL*8                        VEC(0:244444,0:3)
      DO 10 I  = 0, NV-1
      VEC(I,ID)= VEC(I,IS)
   10 CONTINUE
      RETURN
      END
CCC
CC
C
      SUBROUTINE DISPLN( NS, VAL )
      INTEGER            NS
      REAL*8                 VAL(0:40)
      WRITE(6,100) ( VAL(I),I=0,NS-1 )
  100 FORMAT(1H ,8F9.6)
      RETURN
      END
CCC
CC
C
      SUBROUTINE VECINI( NS, NK, NV, NX, EXS, NRM, VEC)
      INTEGER            NS, NK, NV, NX, EXS(0:999)
      REAL*8    PROD,VEC(0:244444,0:3), NRM(0:999)
C
      IF( NK.EQ.0 .OR. NS.EQ.NK*2 ) THEN
          CALL VEC0(   NV, VEC )
          DO 10 I = 0, NV-1
          VEC(I,0)= (MOD(I,13)+MOD(I,19)+00.1) / SQRT(DBLE(NV))
   10     CONTINUE
C
*VDIR NODEP
          DO 20 I = 0, NX-1
          VEC(EXS(I),0)= VEC(EXS(I),0) * NRM(I)
   20     CONTINUE
C
          CALL INPRO(NV,0,0,PROD,VEC)
                PROD = SQRT(PROD)
          CALL DIVER(NV,0,0,PROD,VEC)
        ELSE
          CALL VEC0( 2*NV, VEC )
          DO 30 I = 0, NV-1
          VEC(I   ,0)= (MOD(I,13)+MOD(I,19)+00.1) / SQRT(DBLE(NV))
          VEC(I+NV,0)= (MOD(I,17)+MOD(I,23)+00.1) / SQRT(DBLE(NV))
   30     CONTINUE
C
*VDIR NODEP
          DO 40 I = 0, NX-1
          VEC(EXS(I)   ,0)= VEC(EXS(I)   ,0) * NRM(I)
          VEC(EXS(I)+NV,0)= VEC(EXS(I)+NV,0) * NRM(I)
   40     CONTINUE
C
          CALL INPRO(NV*2,0,0,PROD,VEC)
                PROD  =  SQRT(PROD)
          CALL DIVER(NV*2,0,0,PROD,VEC)
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE  SCALE( ID, IS, NS, LV, NK,NST, NX,EXS,NRM,VEC)
      INTEGER            ID, IS, NS, LV, NK,NST, NX,EXS(0:999)
      REAL*8           NRM(0:999),     VEC(0:244444,0:3)
C
      IF( LV .GE. 3 ) WRITE(6,*) '--- Subroutine SCALE ---'
C
      IF( NK.EQ.0 .OR. NS.EQ.NK*2 ) THEN
          DO 10 I  = 0, NST -1
          VEC(I,ID)= 0.0D0
   10     CONTINUE
*VDIR NODEP
          DO 20 I  = 0, NX -1
          VEC(EXS(I),IS)= VEC(EXS(I),IS) * NRM(I)
   20     CONTINUE
        ELSE
          DO    30  I  = 0, NST -1
          VEC(I    ,ID)= 0.0D0
          VEC(I+NST,ID)= 0.0D0
   30     CONTINUE
*VDIR NODEP
          DO    40 I = 0,NX -1
          VEC(EXS(I)    ,IS)= VEC(EXS(I)    ,IS) * NRM(I)
          VEC(EXS(I)+NST,IS)= VEC(EXS(I)+NST,IS) * NRM(I)
   40     CONTINUE
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE RESCAL( ID, IS, NS, LV, NK,NST, NX,EXS,NRM,VEC)
      INTEGER            ID, IS, NS, LV, NK,NST, NX,EXS(0:999)
      REAL*8       NRM(0:999),     VEC(0:244444,0:3)
C
      IF( LV .GE. 3 ) WRITE(6,*) '--- Subroutine RESCAL ---'
C
      IF( NK.EQ.0 .OR. NS.EQ.NK*2 ) THEN
*VDIR NODEP
          DO 10 I  = 0, NX -1
              VEC(EXS(I),ID)= VEC(EXS(I),ID) * NRM(I)
          IF( NRM(I) .GT. 1.0E-8 ) THEN
              VEC(EXS(I),IS)= VEC(EXS(I),IS) / NRM(I)
            END IF
   10     CONTINUE
        ELSE
*VDIR NODEP
          DO    20 I  = 0, NX -1
              VEC(EXS(I)    ,ID)= VEC(EXS(I)    ,ID) * NRM(I)
              VEC(EXS(I)+NST,ID)= VEC(EXS(I)+NST,ID) * NRM(I)
          IF( NRM(I) .GT. 1.0E-8 ) THEN
              VEC(EXS(I)    ,IS)= VEC(EXS(I)    ,IS) / NRM(I)
              VEC(EXS(I)+NST,IS)= VEC(EXS(I)+NST,IS) / NRM(I)
            END IF
   20     CONTINUE
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE MULT(ID,IS,NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC)
      INTEGER IP,  NW,ID,IS,NS,LV,NK,NX,NST,ST(0:122222)
      INTEGER     VLE,   IX,EXS(0:999)
      INTEGER  V0(0:141072), V2(0:141072), V4(0:141072), V5(0:141072)
      REAL*8   F0(0:141072), F1(0:141072), SN(0:40),  A, CS(0:40)
      REAL*8  DIA(0:122222), VEC(0:244444,0:3), NRM(0:999)
C
      IF( LV .GE. 3 ) WRITE(6,*)  '--- Subroutine MULT ---'
C                                --- Scale with VEC clear ---
      CALL   SCALE( ID, IS, NS, LV, NK,NST, NX, EXS, NRM, VEC)
C
      DO 90 I = 0, NST-1, 141072
           IX = I
          VLE = MIN(141072,NST)
C                                     --- State Copy ---
          DO 10 J  = 0, VLE -1
             V0(J) = ST(IX+J)
   10     CONTINUE
C                                  --- Loop: (Shift)^N ---
      DO 80 N = 0, NS-1
C                             --- Hamiltonian Decomposition ---
      IF( LV .GE. 3 ) WRITE(6,*) 'Rotation Level =', N
      DO 60 K = 0, NS-1
           IP = K
C                                      --- S+S- Part ---
      CALL NNSPSM(  NS,        IP, VLE,       V0,     V2,  A, F0 )
      CALL COMPRS(  ST(NST-1), IX, VLE,  NW,  V2, F0, V4, V5, F1 )
      CALL BISECT(  ST,NST,              NW,  V2, V4             )
C
      IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN
*VDIR NODEP
          DO 15 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
            END IF
   15     CONTINUE
        ELSE
*VDIR NODEP
          DO 20 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
     &                          - VEC(V2(L)+NST,IS)*F1(L)*SN(N)
              VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*SN(N)
     &                          + VEC(V2(L)+NST,IS)*F1(L)*CS(N)
            END IF
   20     CONTINUE
        END IF
C                                      --- S++S-- Part ---
      CALL NNPPMM(  NS,        IP, VLE,       V0,     V2,  A, F0 )
      CALL COMPRS(  ST(NST-1), IX, VLE,  NW,  V2, F0, V4, V5, F1 )
      CALL BISECT(  ST,NST,              NW,  V2, V4             )
C
      IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN
*VDIR NODEP
          DO 25 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
            END IF
   25     CONTINUE
        ELSE
*VDIR NODEP
          DO 30 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
     &                          - VEC(V2(L)+NST,IS)*F1(L)*SN(N)
              VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*SN(N)
     &                          + VEC(V2(L)+NST,IS)*F1(L)*CS(N)
            END IF
   30     CONTINUE
        END IF
C                                    --- S-S+ Part ---
      CALL NNSMSP(  NS,        IP, VLE,       V0,     V2,  A, F0 )
      CALL COMPRS(  ST(NST-1), IX, VLE,  NW,  V2, F0, V4, V5, F1 )
      CALL BISECT(  ST,NST,              NW,  V2, V4             )
C
      IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN
*VDIR NODEP
          DO 35 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
            END IF
   35     CONTINUE
        ELSE
*VDIR NODEP
          DO 40 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
     &                          - VEC(V2(L)+NST,IS)*F1(L)*SN(N)
              VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*SN(N)
     &                          + VEC(V2(L)+NST,IS)*F1(L)*CS(N)
            END IF
   40     CONTINUE
        END IF
C                                    --- S--S++ Part ---
      CALL NNMMPP(  NS,        IP, VLE,       V0,     V2,  A, F0 )
      CALL COMPRS(  ST(NST-1), IX, VLE,  NW,  V2, F0, V4, V5, F1 )
      CALL BISECT(  ST,NST,              NW,  V2, V4             )
C
      IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN
*VDIR NODEP
          DO 45 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
            END IF
   45     CONTINUE
        ELSE
*VDIR NODEP
          DO 50 L = 0, NW -1
          IF( V4(L) .EQ. ST(V2(L)) ) THEN
              VEC(V5(L)    ,ID) = VEC(V5(L)    ,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*CS(N)
     &                          - VEC(V2(L)+NST,IS)*F1(L)*SN(N)
              VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
     &                          + VEC(V2(L)    ,IS)*F1(L)*SN(N)
     &                          + VEC(V2(L)+NST,IS)*F1(L)*CS(N)
            END IF
   50     CONTINUE
        END IF
   60 CONTINUE
C                                --- State Translation ---
      IF( N .EQ. NS-1 ) GOTO 80
      DO 70 J = 0, VLE -1
          V0(J)= MOD(V0(J)*3,3**NS) + V0(J)/3**(NS-1)
   70 CONTINUE
   80 CONTINUE
   90 CONTINUE
      CALL RESCAL( ID, IS, NS, LV, NK,NST, NX,EXS,NRM,VEC)
C
C                                   --- Diagonal Part ---
      IF( NK.EQ.0 .OR. NK.EQ.NS/2 ) THEN 
          DO 95     J    = 0, NST -1
            VEC(    J,ID)= VEC(    J,ID) + DIA(J)*VEC(    J,IS)
   95     CONTINUE
        ELSE
          DO 96     J    = 0, NST -1
            VEC(    J,ID)= VEC(    J,ID) + DIA(J)*VEC(    J,IS)
            VEC(NST+J,ID)= VEC(NST+J,ID) + DIA(J)*VEC(NST+J,IS)
   96     CONTINUE
        END IF
C
      IF( LV .GE. 3 ) THEN
          WRITE(6,*) 'Stopped by Debug Switch LV=', LV
          STOP
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE NNSPSM( NS, IP, VLE, V0, V2, A, F0 )
      INTEGER            NS, IP, VLE, V0(0:141072), V2(0:141072)
      INTEGER    P0, P1, T0, T1
      REAL*8                       A, F0(0:141072)
C
C     --- S+[T1] S-[T0] Part ---
C
      IF( IP .NE. NS-1 ) THEN
               P0 = 3** IP
               P1 = 3**(IP+1)
          DO 10 I = 0, VLE -1
               T0 = MOD(V0(I)/P0,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.NE.2 .AND. T0.NE.0 ) THEN
             V2(I)= V0(I) + P1 - P0 
             F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) )
            ELSE
             F0(I)= 0.0D0
            END IF
   10     CONTINUE
        ELSE
               P0 = 1
               P1 = 3**(NS-1)
          DO 20 I = 0, VLE -1
               T0 = MOD(V0(I)   ,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.NE.2 .AND. T0.NE.0 ) THEN
             V2(I)= V0(I) + P1 - P0 
             F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) )
            ELSE
             F0(I)= 0.0D0
            END IF
   20     CONTINUE
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE NNPPMM( NS, IP, VLE, V0, V2, A, F0 )
      INTEGER            NS, IP, VLE, V0(0:141072), V2(0:141072)
      INTEGER    P0, P1, T0, T1
      REAL*8                       A, F0(0:141072)
C
C     --- S++[T1] S--[T0] Part ---
C
      IF( IP .NE. NS-1 ) THEN
               P0 = 3** IP
               P1 = 3**(IP+1)
          DO 10 I = 0, VLE -1
               T0 = MOD(V0(I)/P0,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.EQ.0 .AND. T0.EQ.2 ) THEN
             V2(I)= V0(I) + P1 + P1 - P0  - P0
             F0(I)=  A
            ELSE
             F0(I)= 0.0D0
            END IF
   10     CONTINUE
        ELSE
               P0 = 1
               P1 = 3**(NS-1)
          DO 20 I = 0, VLE -1
               T0 = MOD(V0(I)   ,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.EQ.0 .AND. T0.EQ.2 ) THEN
             V2(I)= V0(I) + P1 + P1 - P0  - P0
             F0(I)=  A
            ELSE
             F0(I)= 0.0D0
            END IF
   20     CONTINUE
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE NNSMSP( NS, IP, VLE, V0, V2, A, F0 )
      INTEGER            NS, IP, VLE, V0(0:141072), V2(0:141072)
      INTEGER    P0, P1, T0, T1
      REAL*8                  A,      F0(0:141072)
C
C     --- S-[T1] S+[T0] Part ---
C
      IF( IP .NE. NS-1 ) THEN
               P0 = 3** IP
               P1 = 3**(IP+1)
          DO 10 I = 0, VLE -1
               T0 = MOD(V0(I)/P0,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.NE.0 .AND. T0.NE.2 ) THEN
             V2(I)= V0(I) - P1 + P0 
             F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) )
            ELSE
             F0(I)= 0.0D0
            END IF
   10     CONTINUE
        ELSE
               P0 = 1
               P1 = 3**(NS-1)
          DO 20 I = 0, VLE -1
               T0 = MOD(V0(I)   ,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.NE.0 .AND. T0.NE.2 ) THEN
             V2(I)= V0(I) - P1 + P0 
             F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) )
            ELSE
             F0(I)= 0.0D0
            END IF
   20     CONTINUE
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE NNMMPP( NS, IP, VLE, V0, V2, A, F0 )
      INTEGER            NS, IP, VLE, V0(0:141072), V2(0:141072)
      INTEGER    P0, P1, T0, T1
      REAL*8                  A,      F0(0:141072)
C
C     --- S--[T1] S++[T0] Part ---
C
      IF( IP .NE. NS-1 ) THEN
               P0 = 3** IP
               P1 = 3**(IP+1)
          DO 10 I = 0, VLE -1
               T0 = MOD(V0(I)/P0,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.EQ.2 .AND. T0.EQ.0 ) THEN
             V2(I)= V0(I) - P1 - P1 + P0 + P0
             F0(I)=  A
            ELSE
             F0(I)= 0.0D0
            END IF
   10     CONTINUE
        ELSE
               P0 = 1
               P1 = 3**(NS-1)
          DO 20 I = 0, VLE -1
               T0 = MOD(V0(I)   ,3)
               T1 = MOD(V0(I)/P1,3)
           IF( T1.EQ.2 .AND. T0.EQ.0 ) THEN
             V2(I)= V0(I) - P1 - P1 + P0 + P0
             F0(I)=  A
            ELSE
             F0(I)= 0.0D0
            END IF
   20     CONTINUE
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE COMPRS( MST, IX, VLE, NW, V2, F0, V4, V5, F1 )
      INTEGER            MST, IX, VLE, NW, V2(0:141072)
      INTEGER                V4(0:141072), V5(0:141072)
      REAL*8                 F0(0:141072), F1(0:141072)
C
           NW = 0
      DO 10 I = 0, VLE -1
      IF( ABS(F0(I)) .GT. 1.0E-16 .AND. V2(I) .LE. MST ) THEN
          V4(NW) = V2(I)
          F1(NW) = F0(I)
          V5(NW) = IX + I
             NW  = NW + 1
        END IF
   10 CONTINUE
C
      RETURN
      END
CCC
      SUBROUTINE BISECT(     ST, NST, NW, V2, V4 )
      INTEGER DW,  ST(0:122222), NST, NW, V2(0:141072), V4(0:141072)
C
           DW = NST - NST/2
      DO 10 I = 0, NW -1
      IF( V4(I) .GE. ST(DW) ) THEN
          V2(I) = DW
        ELSE
          V2(I) = 0
        END IF
   10 CONTINUE
C
   20 CONTINUE
           DW = DW - DW/2
      DO 30 I = 0, NW -1
      IF( V4(I) .GE. ST( MIN(NST-1,V2(I)+DW) )  ) THEN
          V2(I) =        MIN(NST-1,V2(I)+DW)
        END IF
   30 CONTINUE
      IF( DW .NE. 1 ) GOTO 20
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE LANZOS(NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC )
      INTEGER LK,LLANCZ,NS,LV,NK,NX,NST,NV,EXS(0:999),   ST(0:122222)
      REAL*8  E2,   E2P,  ENERGY,          VEC(0:244444,0:3)
      REAL*8   SN(0:40),CS(0:40),          NRM(0:999),A,DIA(0:122222)
      REAL*8         MER(1000,5),TER(1000), ALP(1000),  BET(1000)
C
      E2  = 0.0D0
C
      IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN
          NV =   NST
        ELSE
          NV = 2*NST
        END IF
C
      CALL   VEC0(          NV,               VEC)
      CALL VECINI( NS, NK, NST, NX, EXS, NRM, VEC)
      CALL  VMOVE(          NV,  3,   0,      VEC)
C (I)
      CALL   MULT( 1,0, NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC)
C (II)
C     ********  LK is an Index for ALP and BET  ********
           LK = 1
       LLANCZ = 0
   10 CONTINUE
C     ********  Report every 5 Steps  ********
          E2P = E2
      DO 20 J = 0, 4
      CALL MACRO1( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC,
     &                                   TER,ALP,BET, LK,  0 )
   20 CONTINUE
      CALL MALANZ( LK, LV, ALP,BET, ENERGY,E2, MER)
          LLANCZ = LLANCZ + 1
      IF( LLANCZ.LE.40 .AND. ABS(E2P-E2).GE.1.0E-16) GOTO 10
          LLANCZ = LLANCZ - 1
       WRITE(6,100)ENERGY
  100 FORMAT(1H , 'ENERGY= ',E30.20)
C
      IF( ABS(ENERGY) .LT. 1.0E-6 ) THEN
          WRITE(6,*) 'Energy is accidentally Zero'
          WRITE(6,*) 'Do not believe the following data'
          RETURN
        END IF
C
      DO 30 I = 1,1000
        TER(I)= MER(I,1)
   30 CONTINUE
C
      CALL VMOVE(  NV, 0, 3,                                VEC)
      CALL ADDER(  NV, 3, 3, 0,TER(1) ,  0.0D0,             VEC)
      CALL  MULT( 1,0, NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC)
      LK=1
      DO 70 I = 0, LLANCZ
      DO 60 J = 0, 4
      CALL MACRO1( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC,
     &                                  TER,ALP,BET, LK,  1   )
   60 CONTINUE
   70 CONTINUE
 
      CALL INVERS( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,ENERGY,VEC)
       WRITE(6,110) ENERGY
  110 FORMAT(1H ,'Energy = ', F20.15)
C
      RETURN
      END
CCC
      SUBROUTINE MACRO1( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC,
     &                                       TER,ALP,BET, LK,ISW )
      INTEGER    ISW, LK,NS,LV,NK,NX,NST,NV,EXS(0:999)
      INTEGER                                ST(0:122222)
      REAL*8          PROD,                 VEC(0:244444,0:3)
      REAL*8      SN(0:40), CS(0:40),  NRM(0:999),A,DIA(0:122222)
      REAL*8     TER(1000), ALP(1000), BET(1000)
C
      IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN
          NV =   NST
        ELSE
          NV = 2*NST
        END IF
C
      CALL INPRO(  NV, 0, 1,         PROD,VEC)
      ALP(LK) =                      PROD
      PROD    =                     -PROD
      CALL ADDER(  NV, 2, 1, 0,1.0D0,PROD,VEC)
      CALL INPRO(  NV, 2, 2,         PROD,VEC)
      PROD    =                 SQRT(PROD)
      BET(LK) =                      PROD
      CALL DIVER(  NV, 2, 2,         PROD,VEC)
      CALL  MULT(      1, 2,NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC)
      PROD    =                     -PROD
      CALL ADDER(  NV, 1, 1, 0,1.0D0,PROD,VEC)
      CALL VMOVE(  NV, 0, 2,              VEC)
          LK  = LK + 1
      IF(ISW.EQ.1) CALL ADDER(  NV, 3, 3, 0, 1.0D0, TER(LK), VEC)
C
      RETURN
      END
CCC
      SUBROUTINE MALANZ( LK, LV,  ALP, BET, ENERGY, E2, VE)
      INTEGER  N,M, LNV, LK, LV,  ISW,IERR,         IW1(1000)
      REAL*8     ENERGY, E2,ALP(1000), BET(1000)
      REAL*8  TP,   EPS,      E(1000),  W1(6000),   VE(1000,5)
C
        N = LK -1
      EPS =-1.0
        M = 5
      LNV = 1000
      ISW =-1
       TP = BET(N)
      CALL BCSTSS(ALP,N,BET,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR)
      BET(N)= TP
      IF(IERR.NE.0) THEN
          WRITE(6,100) IERR
  100     FORMAT(1H ,' Error No. ',I10,' in RCSTSS')
        END IF
      ENERGY= E(1)
      E2    = E(2)
      IF( LV .GE. 2 ) WRITE(6,110) (E(J),J=1,3)
  110 FORMAT(1H ,'*****',3F20.15,'*****')
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE INVERS(NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,LMD,VEC)
      INTEGER      LROT,NS,LV,NK,NX,NST,NV, EXS(0:999)
      INTEGER                                ST(0:122222)
      REAL*8            LMD,A,DIA(0:122222),VEC(0:244444,0:3)
      REAL*8         SN(0:40),CS(0:40),DTP,  DTR,  NRM(0:999)
C
      IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN
          NV =   NST
        ELSE
          NV = 2*NST
        END IF
C
      CALL VMOVE(  NV, 0, 3, VEC)
      LROT = 0
      GOTO 20
C
   10 CONTINUE
      CALL CGMETH(NS,LV,NK,   NX,NST,ST,EXS,SN,CS,NRM,A,DIA,LMD,VEC)
   20 CALL VMOVE(               NV, 1,  0,                   VEC)
      CALL INPRO(               NV, 1,  1,        DTP,       VEC)
      DTP = SQRT( DTP )
      CALL DIVER(               NV, 1,  1,        DTP,       VEC)
      CALL  MULT( 2,1, NS, LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC)
      CALL INPRO(               NV, 2,  2,        DTP,       VEC)
      DTP = SQRT( DTP )
      CALL DIVER(               NV, 2,  2,        DTP,       VEC)
C
      IF( LMD .LT. 0 ) THEN
      CALL ADDER(               NV, 3,  2,  1,  1.0D0, 1.0D0,VEC)
        ELSE
      CALL ADDER(               NV, 3,  2,  1, -1.0D0, 1.0D0,VEC)
        END IF
C
      CALL INPRO(               NV, 3,  3,        DTR,       VEC)
      WRITE(6,100)  LMD, DTP, SQRT(DTR)
      LROT = LROT + 1
  100 FORMAT(1H , 'LMD=',E20.10, '  ?  ', E20.10, '  D  ',E20.10)
      IF( LROT .GE. 20 ) THEN
          WRITE(6,110)
  110     FORMAT(1H ,'Did not converge in INVERS')
          STOP
        END IF
      IF( SQRT(DTR) .GE. 0.5E-11 ) GOTO 10
C
      RETURN
      END
CCC
      SUBROUTINE CGMETH(NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,LMD,VEC)
      INTEGER NV,  LROT,NS,LV,NK,NX,NST,EXS(0:999), ST(0:122222)
      REAL*8   LMD,DB,A,ALP,BET, VEC(0:244444,0:3),DIA(0:122222)
      REAL*8   EPS, DTP,SN(0:40), CS(0:40), DR1,DR2, NRM(0:999)
C
      IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN
          NV =   NST
        ELSE
          NV = 2*NST
        END IF
C (I)
      CALL INPRO(  NV, 1, 1, DB, VEC)
      CALL VMOVE(  NV, 2, 1,     VEC)
      DO 10 I = 0, NV -1
      VEC(I,0)= 0.0D0
   10 CONTINUE
C (II)
       EPS = 1.0E-8
      LROT = 0
      CALL INPRO(               NV, 1,  1,        DR1,       VEC)
   20 CALL  MULT( 3,2, NS, LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC)
      DO 30 I = 0, NV -1
      VEC(I,3)= VEC(I,3) - LMD*VEC(I,2)
   30 CONTINUE
      CALL INPRO(               NV, 3,  2,        DTP,       VEC)
      ALP =  DR1 / DTP
      CALL ADDER(               NV, 0,  0,  2,  1.0D0,   ALP,VEC)
      CALL ADDER(               NV, 1,  1,  3,  1.0D0,  -ALP,VEC)
      CALL INPRO(               NV, 1,  1,        DR2,       VEC)
C
C-WR- WRITE(6,100) DR2 / DB
C 100 FORMAT(1H ,E30.20)
      IF( DR2/DB .LE. EPS ) GOTO 40
C
       BET = DR2 / DR1
       DR1 = DR2
      LROT = LROT + 1
      CALL ADDER(               NV, 2,  1,  2,  1.0D0,   BET,VEC)
      IF( LROT .GE. 20 ) THEN
          WRITE(6,110)
  110     FORMAT(1H ,'Did not converge in CG')
          GOTO 50
        END IF
      GOTO 20
   40 CONTINUE
   50 WRITE(6,120) LROT
  120 FORMAT(1H ,I10,'After Rotation in CG')
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE VECDMP( NS, LV, NK, NST, ST, VEC )
      INTEGER        IP, NS, LV, NK, NST, ST(0:122222)
      REAL*8                    FTP, VEC(0:244444,0:3)
C
      IF( LV  .EQ.   0 ) RETURN
C
      IF( NST .GT. 500 ) THEN
          WRITE(6,*) ' Too Many Basis '
          RETURN
        END IF
C
      DO 10 I = 0, NST -1
      IF( NK.EQ.0 .OR. NK.EQ.NS/2 ) THEN
          VEC(I,3) = VEC(I,1)*VEC(I,1)
        ELSE
          VEC(I,3) = VEC(I,1)*VEC(I,1) + VEC(NST+I,1)*VEC(NST+I,1)
        END IF
   10 CONTINUE
C
      DO 30 I = 0, MIN(25,NST) -1
          FTP = 0.0D0
           IP = 0
      DO 20 J = 0, NST -1
      IF( FTP .LT. VEC(J,3) ) THEN
          FTP =    VEC(J,3)
           IP = J
        END IF
   20 CONTINUE
      IF( FTP .GT. 1.0E-8 ) CALL STDISP( NS, ST(IP), VEC(IP,3) )
      VEC(IP,3) = 0.0D0
   30 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE SZDIST( NS, NK, NST, ST, BH, BL, VEC )
      INTEGER     TP, M, NS, NK, NST, ST(0:122222)
      INTEGER  TH,TL,BH(0:59049,0:2), BL(0:59049,0:2)
      REAL*8         FC, F2, F1,  F0,VEC(0:244444,0:3)
C
       M = NS /2
      TP = 3** M
      F0 = 0.0D0
      F1 = 0.0D0
      F2 = 0.0D0
C
      IF( NK.EQ.0 .OR. NK.EQ.NS/2 ) THEN 
          DO 10 I = 0, NST -1
               TH =     ST(I)/TP
               TL = MOD(ST(I),TP)
               FC = VEC(I,1) * VEC(I,1)
               F0 = F0 + FC*( BH(TH,0) + BL(TL,0) )
               F1 = F1 + FC*( BH(TH,1) + BL(TL,1) )
               F2 = F2 + FC*( BH(TH,2) + BL(TL,2) )
   10     CONTINUE
        ELSE
          DO 20 I = 0, NST -1
               TH =     ST(I)/TP
               TL = MOD(ST(I),TP)
               FC = VEC(I,1) *VEC(I,1) + VEC(NST+I,1) *VEC(NST+I,1)
               F0 = F0 + FC*( BH(TH,0) + BL(TL,0) )
               F1 = F1 + FC*( BH(TH,1) + BL(TL,1) )
               F2 = F2 + FC*( BH(TH,2) + BL(TL,2) )
   20     CONTINUE
        END IF
C
      WRITE(6,100) F0/NS, F1/NS, F2/NS
  100 FORMAT(1H ,'v 0 ^ Ratio = ', 3F15.7)
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE BCSTSS(A,N,B,EPS,E,M,VE,LNV,ISW,WN,W1,IERR)
      INTEGER MWN,NC,N,M,SP,BP,ISW,LNV,IERR,WN(0:1000)
      REAL*8 EPS,A(0:1000),B(0:1000),W1(0:6000),E(0:1000)
      REAL*8 MIN,MDL,BB(0:1000),VE(0:1000,0:2)
C -DUMMY-
      IERR = ISW
      IERR = LNV
      IERR = 0
        BP = 0
      MIN = ABS(A(0)) + ABS(B(0))
      DO 10 I = 0,N -3
        MIN = MAX(MIN,ABS(B(I))+ABS(A(I+1))+ABS(B(I+1)))
        BB(I)= B(I)*B(I)
   10 CONTINUE
      MIN = MAX(MIN,ABS(B(N-2))+ABS(A(N-1)))
      BB(N-2)= B(N-2)*B(N-2)
C
      CALL COUNT(N,NC,-MIN,A,BB)
      MIN = -MIN
      MWN = NC
      CALL COUNT(N,NC,-MIN,A,BB)
      SP = 0
      W1(SP)= -MIN
      WN(SP)= NC
C
   20 CONTINUE
      MDL =(MIN+W1(SP)) / 2
      IF(W1(SP)-MIN.LT.EPS.OR.W1(SP)-MDL.LT.ABS(MDL)*1.0E-17.OR.MDL
     & -MIN.LT.ABS(MDL)*1.0E-17)GOTO 25
C
      CALL COUNT(N,NC,MDL,A,BB)
C -WR- WRITE(6,198) SP,MWN-WN(SP), W1(SP)-MDL, MDL-MIN
C  198 FORMAT(1H , 'SP,MWN-WN(SP), W1(SP)-MDL, MDL-MIN', 2I3, 2E20.10 )
      IF(MWN-NC.EQ.0)THEN
        MIN = MDL
        MWN = NC
        GOTO 20
      END IF
C
      SP = SP + 1
      W1(SP)= MDL
      WN(SP)= NC
      GOTO 20
C
   25 CONTINUE
C
      DO 30 I = 0,MWN - WN(SP) - 1
        E(BP+I)= MDL
   30 CONTINUE
      BP = BP + MWN - WN(SP)
C
      IF(BP.LT.M) THEN
   40   CONTINUE
        IF(WN(SP).EQ.WN(SP-1)) THEN
          SP = SP-1
          GOTO 40
        END IF
        MIN = W1(SP)
        MWN = WN(SP)
        SP = SP - 1
        GOTO 20
      END IF
C
      CALL INVERT(N,E(0),A,B,W1,VE)
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE COUNT(N,NC,LAM,A,BB)
      INTEGER I,N,NC
      REAL*8 G(0:1000),GRN,EPS,LAM,A(0:1000),BB(0:1000)
C
      NC = 0
      EPS = 1.0E-18
      GRN = 1.0E+18
      G(0) = LAM - A(0)
      DO 10 I = 1,N- 1
        IF(ABS(G(I-1)).GT.EPS) THEN
          G(I) = LAM - A(I) - BB(I-1)/G(I-1)
        ELSE
          G(I) = GRN
        END IF
   10 CONTINUE
C
      DO 20 I = 0,N-1
        IF(G(I).LE.0) NC = NC + 1
   20 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE INVERT(N,E,A,B,W1,VE)
      INTEGER N
      REAL*8 PROD,MROD
      REAL*8 E,A(0:1000),B(0:1000),W1(0:1000,0:5),VE(0:1000,0:2)
C
      VE(0,0)= 1.0D0
      DO 10 I = 1,N-1
        VE(I,0)= 0.0D0
   10 CONTINUE
C
      CALL CHRSKY(N,E,A,B,W1)
C
   20 CONTINUE
      DO 30 I = 0,N-1
        VE(I,1)= VE(I,0)
   30 CONTINUE
      CALL VSOLVE(N,W1,VE)
C
      PROD = 0.0D0
      DO 40 I = 0,N-1
        PROD = PROD + VE(I,0)*VE(I,0)
   40 CONTINUE
C
      PROD = SQRT(PROD)
      DO 50 I = 0,N-1
        VE(I,0)= VE(I,0) / PROD
   50 CONTINUE
C
      PROD = 0.0D0
      MROD = 0.0D0
      DO 60 I = 0,N-1
        PROD = PROD +(VE(I,1)-VE(I,0))*(VE(I,1)-VE(I,0))
        MROD = MROD +(VE(I,1)+VE(I,0))*(VE(I,1)+VE(I,0))
   60 CONTINUE
C -WR- WRITE(6, * )' MROD, PROD ', MROD, PROD
      IF(PROD.GT.1.0E-16.AND.MROD.GT.1.0E-16) GOTO 20
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE CHRSKY(N,E,A,B,W1)
      INTEGER N
      REAL*8 E,A(0:1000),B(0:1000),W1(0:1000,0:5)
C
      W1(0,0)= A(0)- E
      DO 10 I = 0,N - 2
        W1(I,1)= B(I)/W1(I,0)
        W1(I+1,0)= A(I+1)- E - B(I)*W1(I,1)
   10 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE VSOLVE(N,W1,VE)
      INTEGER N
      REAL*8 VE(0:1000,0:2),W1(0:1000,0:5)
C
      DO 20 I = 1,N-1
        VE(I,0)= VE(I,0) - W1(I-1,1)*VE(I-1,0)
   20 CONTINUE
C
      DO 30 I = 0,N-1
        VE(I,0)= VE(I,0)/W1(I,0)
   30 CONTINUE
C
      DO 40 I = N -2,0,-1
        VE(I,0)= VE(I,0) - W1(I,1)*VE(I+1,0)
   40 CONTINUE
C
      RETURN
      END