This is the source code of DMRG for S=1/2
Heisenberg Chain with Next Nearest
Neighbor Interaction.

This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)

This program is really primitive.Use it just for tutorial purpouse.

*#RUN:L=LIB/ASL7 OPT=3
C*#RUN:L=LIB/ASL7 OPT=0 NIAP ARGCHK SUBCHK
C123456789012345678901234567890123456789012345678901234567890123456789012
C
C     Sannomiya Kyuukou Version: as fast as the Hankyu Train.
C     Version 1.1 1995/1/15 by Nishino Tomotoshi: Free Soft.
C
C     ***** Subroutines whose name begins with "DCS" belong to *****
C     ***** Center Library of SX/ACOS (ASL)                    *****
C
C     +----  Density Matrix (Numerical) Renormalization Group  -----
C     -----  for the J-J' Heisenberg Chain. (Ladder mo dekiru!)-----
C     -----              ________   ________                   -----
C     -----             /   Joint\ /        \                  -----
C     ----- [L-Block===E====J]====O====O====[J====E===R-Block] -----
C     -----            End   \________/ \________/             -----
C     -----                                                    -----
C     -----         === J(0), __/ J(1), Ising D(0), XY D(1)    -----
C     -----             {Internal Expression: J(2),    J(3)}   -----
C     -----                                                    -----
C     +---- Special #: 0-1 -- Left/Right Sub System            ----+
C     |                100 -- Max.# of Site (0 to 99)              |
C     |             -1~101 -- Max.# of Sub Space = 100+1           |
C     |              10000 -- Matrix Stack Space (10Mbyte/Operator)|
C     |              40000 -- Extended Matrix Stack Space          |
C     +------------ 100000 -- Max. Lanczos Vector Length ----------+
C
      INTEGER  LV,  NS,  NSB,  NL,  MX,   DX, LP, NP,  RT
      INTEGER       TZ,        NR,  MM(0:16), RP, IP, MRT
      INTEGER  SZ(0:100,0:1),  MP(0:100,0:5),     X(-1:101,0:3,0:100)
      REAL*8    E0(0:3,0:16),  BH(0:40000,0:1), SM(0:10000,0:5,0:100)
      REAL*8   T1,TIME(0:16),  JH(0:4), E(0:3), HM(0:10000,0:100)
      REAL*8  FMX,  FM(0:16),  ER, ERE,          V(0:100000,0:3)
C
C     +------------------------------------------------------------+
C     |      Parameter Input and Initialization of Variables       |
C     +------------------------------------------------------------+
C .......... Debug Level: 0-None, 1-Normal, 2-Detail, 3-Dump, 4-Debug
C                         5-FullDump-Debug, 6-Stop
       LV = 0
C .......... Desiable Accuracy of the Density Matrix (Desiable...)
       ER = 1.0E-15
C .......... Desiable Accuracy of the Ground State Energy
      ERE = 1.0E-7
C .......... Minimum Step Number in Lanczos Diagonalization
C            when we consider Excitations.
      LLS = 120
C .......... Exchange Parameters (N.N. J(0); N.N.N. J(1))
      JH(0)= 1.00D0
      JH(1)= 1.00D0
C .......... Ising/XY Strength  (Ising J(2); XY J(3))
      JH(2)= 1.00D0
      JH(3)= 1.00D0
C .......... Number of Sites; it should be even number. (Ns = 100)
      DO 90 NS = 16, 16,  2
C .......... Total Sz (= an Integer, since Ns should be even.)
      DO 90 TZ =  0,  1
C .......... Maximum Number of Iteration (Max RT)
      MRT= 10
C .......... Initial "MAXIMUM Number" of Basis Kept
      MX = 20
C .......... Increase Rate of MX during the Finite Size Iteration
      DX = 10
C .......... Parameter Display/ Parameter Check
      CALL PRMDSP( LV, NS, TZ, MX, ER, JH)
C .......... Make Initial Blocks: BL [O===O]---, and ---[O===O] BR
      CALL INIBLK( LV, JH, X(-1,0,   1), SM(0,0,   1), HM(0,   1))
      CALL INIBLK( LV, JH, X(-1,0,NS-2), SM(0,0,NS-2), HM(0,NS-2))
C
C     +------------------------------------------------------------+
C     |    Iteration of Diagonalization and Renormalizatioin       |
C     +------------------------------------------------------------+
C ...........Infinite Chain Iteration to NS-2: (Rotation=0)
      RT = 0
      CALL GETSEC( T1 )
C ........... Dummy .....
      IP = 0
      FMX= 0.0D0
C ...........Initial Block Size = 2 for both L/R Blocks.
      NL = 2
      NR = 2
   10 CONTINUE
      LP = NL - 1
      RP = NS - NR
      CALL MATSTR(LV, NL,NR,TZ,NSB, X(-1,0,LP), X(-1,0,RP),SZ, MP)
      CALL BLOCKH(LV,NSB,   SZ, MP, X(-1,0,LP), X(-1,0,RP),JH,
     &           HM(0,LP),HM(0,RP),SM( 0,0,LP),SM( 0,0,RP),BH    )
      IF(NL .EQ. NS/2-1 ) THEN
          LS = LLS
        ELSE
          LS = 10
        END IF
      CALL LANZOS(LV,NSB,LS,SZ, MP, X(-1,0,LP), X(-1,0,RP),JH,
     &                             SM( 0,0,LP),SM( 0,0,RP),BH,E,V)
      IF( LV .GE. 1 ) WRITE(6,100) NL, NR,E(0), IP, FMX
      NL = NL + 1
      NP = LP + 1
      CALL RENRML(LV,NL, NSB,SZ(0,0),MP,X(-1,0,LP),X(-1,0,NP),
     &            SM(0,0,LP),SM(0,0,NP),   BH(0,0),  HM(0,NP),
     &                                     MX,  IP,  ER,  FMX, V)
      NR = NR + 1
      NP = RP - 1
      CALL RENRMR(LV,NR, NSB,SZ(0,1),MP,X(-1,0,RP),X(-1,0,NP),
     &            SM(0,0,RP),SM(0,0,NP),   BH(0,1),  HM(0,NP),
     &                                     MX,  IP,  ER,  FMX, V)
          IF( LV.GE.  6 ) STOP
      IF( NL+NR .LT. NS ) GOTO 10
          IF( LV.GE.  5 ) STOP
C ............ Second Step: Finite Chain Iteration for LEFT Block ....
         NR = NS - NL - 2
C
      E0(0,RT)= E(0)
      E0(1,RT)= E(1)
      E0(2,RT)= E(2)
      E0(3,RT)= E(3)
        MM(RT)= IP
        FM(RT)= FMX
      CALL EXTIME( T1, TIME(RT) )
      WRITE(6,*)
      WRITE(6,110) RT, MM(RT), E0(0,RT), E0(1,RT), E0(2,RT),
     &                                     FM(RT), TIME(RT)
         RT = RT + 1
         MX = MX + DX
   20 CONTINUE
      LP = NL - 1
      RP = NS - NR
      IF( NR .EQ. 2 ) GOTO 30
      CALL MATSTR(LV, NL,NR,TZ,NSB, X(-1,0,LP), X(-1,0,RP),SZ, MP)
      CALL BLOCKH(LV,NSB,   SZ, MP, X(-1,0,LP), X(-1,0,RP),JH,
     &           HM(0,LP),HM(0,RP),SM( 0,0,LP),SM( 0,0,RP),BH    )
      IF(NL .EQ. NS/2-1 ) THEN
          LS = LLS
        ELSE
          LS = 10
        END IF
      CALL LANZOS(LV,NSB,LS,SZ, MP, X(-1,0,LP), X(-1,0,RP),JH,
     &                             SM( 0,0,LP),SM( 0,0,RP),BH,E,V)
      IF( LV .GE. 1 ) WRITE(6,100) NL, NR,E(0), IP, FMX
      IF( NL .EQ. NS/2-1 ) THEN
          E0(0,RT)= E(0)
          E0(1,RT)= E(1)
          E0(2,RT)= E(2)
          E0(3,RT)= E(3)
            MM(RT)= IP
            FM(RT)= FMX
          CALL EXTIME( T1, TIME(RT) )
           WRITE(6,111) RT, MM(RT), E0(0,RT), E0(1,RT), E0(2,RT),
     &                                          FM(RT), TIME(RT)
C           ( Please Do Not Increase MX here! )
          IF((ABS(E0(0,RT)-E0(0,RT-1)).LT.ERE .AND. MM(RT).GT.MM(RT-1))
     &                                         .OR. RT.GT.MRT) GOTO 40
          RT  = RT + 1
        END IF
      NL = NL + 1
      NP = LP + 1
      NR = NR - 1
      CALL RENRML(LV,NL, NSB,SZ(0,0),MP,X(-1,0,LP),X(-1,0,NP),
     &            SM(0,0,LP),SM(0,0,NP),   BH(0,0),  HM(0,NP),
     &                                     MX,  IP,  ER,  FMX, V)
      GOTO 20
C ............ Third Step: Finite Chain Iteration for LEFT Block ....
   30 CONTINUE
      LP = NL - 1
      RP = NS - NR
      IF( NL .EQ. 2 ) GOTO 20
      CALL MATSTR(LV, NL,NR,TZ,NSB, X(-1,0,LP), X(-1,0,RP),SZ, MP)
      CALL BLOCKH(LV,NSB,   SZ, MP, X(-1,0,LP), X(-1,0,RP),JH,
     &           HM(0,LP),HM(0,RP),SM( 0,0,LP),SM( 0,0,RP),BH    )
      IF(NL .EQ. NS/2-1 ) THEN
          LS = LLS
        ELSE
          LS = 10
        END IF
      CALL LANZOS(LV,NSB,LS,SZ, MP, X(-1,0,LP), X(-1,0,RP),JH,
     &                             SM( 0,0,LP),SM( 0,0,RP),BH,E,V)
      IF( LV .GE. 1 ) WRITE(6,100) NL, NR,E(0), IP, FMX
      IF( NL .EQ. NS/2 - 1 ) THEN
          E0(0,RT)= E(0)
          E0(1,RT)= E(1)
          E0(2,RT)= E(2)
          E0(3,RT)= E(3)
            MM(RT)= IP
            FM(RT)= FMX
          CALL EXTIME( T1, TIME(RT) )
           WRITE(6,111) RT, MM(RT), E0(0,RT), E0(1,RT), E0(2,RT),
     &                                          FM(RT), TIME(RT)
          IF((ABS(E0(0,RT)-E0(0,RT-1)).LT.ERE .AND. MM(RT).GT.MM(RT-1))
     &                                         .OR. RT.GT.MRT) GOTO 40
          RT = RT + 1
          MX = MX + DX
        END IF
      NR = NR + 1
      NP = RP - 1
      NL = NL - 1
      CALL RENRMR(LV,NR, NSB,SZ(0,1),MP,X(-1,0,RP),X(-1,0,NP),
     &            SM(0,0,RP),SM(0,0,NP),   BH(0,1),  HM(0,NP),
     &                                     MX,  IP,  ER,  FMX, V)
      GOTO 30
C ............ Observation Routines ..................................
C     Sweep Once More and Create SL/SR for Observation.
   40 CONTINUE
C
      WRITE(6,112)           ERE
C     ............ Parameter Do-End ..................................
   90 CONTINUE
  100 FORMAT(1H ,I3,'-o-o-',I3,' E0:',F21.14,' m:',I3,' Min[DM]:',E10.2)
  110 FORMAT(1H ,I2,' m:',I3,' E',3F16.11,' co.',E7.1,' ',F6.2,'s')
  111 FORMAT(1H ,I2,'   ',I3,'  ',3F16.11,' co.',E7.1,' ',F6.2,'s')
  112 FORMAT(1H ,    'Conv.acc->',3F16.11)
C
      Write(6,*) 'Observation Routines has not been coded yet.'
      CALL MEMDSP( LV, NL, NSB, MP, X(-1,0,LP) )
      Write(6,*) 'Bye-Bye!'
C
      STOP
      END
CCC
CC
C
      SUBROUTINE GETSEC( TIME )
      REAL*8             TIME
C                 -NEC-
            CALL  CLOCK( TIME )
      RETURN
      END
CCC
      SUBROUTINE EXTIME( T1, T2 )
      REAL*8             T1, T2, T3
            CALL GETSEC( T3 )
            T2 = T3 - T1
            T1 = T3
      RETURN
      END
CCC
      SUBROUTINE MEMDSP( LV, N, NSB, MP, X )
      INTEGER            LV, N, NSB, MP(0:100,0:5), X(-1:101,0:3)
      IF( LV .GE. 0 ) THEN
          WRITE(6,*)
          Write(6,*) 'Typical Memory Requirement for each Variable.'
          Write(6,*) 'V ------ Vector Length: ', MP(NSB,3)
          Write(6,*) 'HM,SM --- Matrix Stack: ',  X(N+1,2)
          Write(6,*) 'BH(L/R) - Matrix Stack: ',  X(N+1,2)*4
        END IF
      RETURN
      END
CCC
      SUBROUTINE PRMDSP( LV, NS, TZ, MX, ER, JH)
      INTEGER            LV, NS, TZ, MX
      REAL*8                             ER, JH(0:4)
C
      WRITE(6,*)
      IF( LV .GE. 0 ) WRITE(6,100)   NS, TZ, JH(0),  JH(1)
      IF( LV .GE. 1 ) WRITE(6,*)'initial-m:',   MX,' DM-cut-off:',ER
      IF( LV .GE. 1 ) WRITE(6,*)    'Ising:',JH(2),       '  XY:',JH(3)
  100 FORMAT(1H ,'NS=',I3,' SZ=',I2,' J1=',F5.2,' J2=',F5.2 )
      IF( NS .LT. 8 ) THEN
          Write(6,*) 'Lattice Size is TOO SMALL. Set NS >= 8.'
          CALL ERRSTP( 'PRMDSP' )
        END IF
      IF( TZ .GE. 6 ) THEN
          Write(6,*) 'Total Sz is Larger than (Initial System Size)/2.'
          Write(6,*) 'You have to modify MAIN-ROUTINE! (Sorry....) '
          CALL ERRSTP( 'PRMDSP' )
        END IF
      IF( NS .GT. 100 ) THEN
          Write(6,*) 'Lattice Size is TOO LARGE. Set NS = 100,',
     &                              ' or increase Memory Space.'
          CALL ERRSTP( 'PRMDSP' )
        END IF
      IF( MOD(NS,2) .EQ. 1 ) THEN
          Write(6,*) 'Ns should be an even number.'
          CALL ERRSTP( 'PRMDSP' )
        END IF
      IF( ABS(TZ) .GE. NS/2 ) THEN
          Write(6,*) 'Total Sz (Var.TZ) should be smaller than NU.'
          CALL ERRSTP( 'PRMDSP' )
        END IF
      RETURN
      END
CCC
      SUBROUTINE MATDSP( N, X, ME, CH1, CH2 )
      INTEGER            N, X(-1:101,0:3)
      REAL*8         ME(0:*)
      CHARACTER*3                  CH1, CH2
      CHARACTER*1     A(0:40,0:40),IMG
C
      IF( X(N+1,1) .GE. 40 ) THEN
          Write(6,*)'Size (',X(N+1,1),') is too large to display.'
          RETURN
        END IF
C
      DO 10 J = 0, X(N+1,1)-1
      DO 10 I = 0, X(N+1,1)-1
        A(I,J)= '.'
   10 CONTINUE
C
      IF( CH1 .EQ. 'BLD' ) THEN
          DO 20 K = 0, N
          DO 20 J = 0, X(K,0)-1
          DO 20 I = 0, X(K,0)-1
          A( X(K,1)+I, X(K,1)+J) = IMG( ME(X(K,2)+J*X(K,0)+I) )
   20     CONTINUE
        END IF
      IF( CH1 .EQ. 'BL+' ) THEN
          DO 30 K = 0, N-1
          DO 30 J = 0, X(K+1,0)-1
          DO 30 I = 0, X(K  ,0)-1
          A( X(K,1)+I, X(K+1,1)+J) = IMG( ME(X(K,3)+J*X(K,0)+I) )
   30     CONTINUE
        END IF
      IF( CH1 .EQ. 'BL-' ) THEN
          DO 40 K = 0,  N-1
          DO 40 J = 0,  X(K  ,0)-1
          DO 40 I = 0,  X(K+1,0)-1
          A( X(K+1,1)+I,X(K,1)+J) = IMG( ME(X(K,3)+J*X(K+1,0)+I) )
   40     CONTINUE
        END IF
C
      Write(6,*) '---',CH2,'---'
      DO 50 I = 0, X(N+1,1)-1
   50 Write(6,100) (A(I,J), J=0, X(N+1,1)-1)
  100 FORMAT(1H , 39A2 )
      RETURN
      END
CCC
      CHARACTER*1 FUNCTION IMG( F )
      REAL*8                    F
                           IMG = ' '
      IF( F .LT. -1.0E-8 ) IMG = '-'
      IF( F .LT. -1.0E-4 ) IMG = '='
      IF( F .GT.  1.0E-8 ) IMG = '+'
      IF( F .GT.  1.0E-4 ) IMG = '#'
      RETURN
      END
CCC
      SUBROUTINE MDUMP( N, X, H, S )
      INTEGER           N, X(-1:101,0:3)
      REAL*8       H(0:*), S(0:10000,0:5)
      CALL MATDSP( N, X,      H, 'BLD', '-H-' )
      CALL MATDSP( N, X, S(0,0), 'BLD', 'SzE' )
      CALL MATDSP( N, X, S(0,1), 'BL+', 'S+E' )
      CALL MATDSP( N, X, S(0,2), 'BL-', 'S-E' )
      CALL MATDSP( N, X, S(0,3), 'BLD', 'SzJ' )
      CALL MATDSP( N, X, S(0,4), 'BL+', 'S+J' )
      CALL MATDSP( N, X, S(0,5), 'BL-', 'S-J' )
      RETURN
      END
CCC
CC
C
      SUBROUTINE INIBLK( LV, JH, X, S, H )
      INTEGER            LV,     X(-1:101,0:3)
      REAL*8         H(0:*), JH(0:3),  S(0:10000,0:5)
C ....[Up,Up] -- [Up,Dw][Dw,Up] -- [Dw,Dw]: Set Sub-Block Size
      X(-1,0)= 0
      X( 0,0)= 1
      X( 1,0)= 2
      X( 2,0)= 1
      X( 3,0)= 0
      CALL POINTR( LV, 2, X )
C ..................................Set Initial Hamiltonian
      H(0)= JH(0)*JH(2)/4.0D0
      H(1)=-JH(0)*JH(2)/4.0D0
      H(2)= JH(0)*JH(3)/2.0D0
      H(3)= JH(0)*JH(3)/2.0D0
      H(4)=-JH(0)*JH(2)/4.0D0
      H(5)= JH(0)*JH(2)/4.0D0
C .............(Send^Z)..Set Spin Operator at the End/Joint point.
C     End Point (0,1,2), Joint Point (3,4,5): [E===J]===O===O
      S(0,0)= 1.0D0/2.0D0
      S(1,0)= 1.0D0/2.0D0
      S(2,0)= 0.0D0
      S(3,0)= 0.0D0
      S(4,0)=-1.0D0/2.0D0
      S(5,0)=-1.0D0/2.0D0
C .........................(Send^Plus/Minus)
      S(0,1)= 0.0D0
      S(1,1)= 1.0D0
      S(2,1)= 1.0D0
      S(3,1)= 0.0D0
      CALL HCONJU( 2, X, S(0,1), S(0,2) )
C .........................(Sjoint^Z)
      S(0,3)= 1.0D0/2.0D0
      S(1,3)=-1.0D0/2.0D0
      S(2,3)= 0.0D0
      S(3,3)= 0.0D0
      S(4,3)= 1.0D0/2.0D0
      S(5,3)=-1.0D0/2.0D0
C .........................(Sjoint^Plus/Minus)
      S(0,4)= 1.0D0
      S(1,4)= 0.0D0
      S(2,4)= 0.0D0
      S(3,4)= 1.0D0
      CALL HCONJU( 2, X, S(0,4), S(0,5) )
C
      IF( LV .GE. 4 ) CALL MDUMP( 2, X, H, S )
      RETURN
      END
CCC
      SUBROUTINE HCONJU( N, X, MP, MM )
      INTEGER            N, X(-1:101,0:3)
      REAL*8       MP(0:*), MM(0:*)
      DO 10 I = 0, N-1
      DO 10 J = 0, X(I+1,0)-1
      DO 10 K = 0, X(I  ,0)-1
      MM(X(I,3)+X(I+1,0)*K+J)= MP(X(I,3)+X(I,0)*J+K)
   10 CONTINUE
      RETURN
      END
CCC
CC
C
      SUBROUTINE POINTR( LV, N, X )
      INTEGER            LV, N, X(-1:101,0:3)
C
      IF( N .GE. 100 ) THEN
          Write(6,*) 'System Size Over! ',
     &               'Increase Buffer Size [100].'
          CALL ERRSTP( 'POINTR' )
        END IF
C .....................Set Pointer to the Linear Matrix Address
       X(-1,1)= 0
      DO 10 I = 0, N+1
       X( I,1)= X(I-1,1) + X(I-1,0)
   10 CONTINUE
C .....................Set Pointer to Block-Diagonal Matrix-Element
       X(-1,2)= 0
      DO 20 I = 0, N+1
       X( I,2)= X(I-1,2) + X(I-1,0)**2
   20 CONTINUE
C .....................Set Pointer to Block-Offdiagonal Matrix-Element
       X(-1,3)= 0
      DO 30 I = 0, N
       X( I,3)= X(I-1,3) + X(I-1,0)*X(I,0)
   30 CONTINUE
C .......................... Pointer Display/Check ...
      IF( LV .GE. 3 ) THEN
          Write(6,*) ' Pointers to Matrix Structure: L/R'
          DO 40 I = -1, N+1
          IF( X(I,0) .NE. 0 ) Write(6,100) I,X(I,0),X(I,1),X(I,2),X(I,3)
  100     FORMAT(1H ,5I8)
   40     CONTINUE
        END IF
      IF( LV .GE. 2 ) THEN
          Write(6,*)'Pointer Maximums:'
          Write(6,100) N+1, X(N+1,0), X(N+1,1), X(N+1,2), X(N,3)
        END IF
      IF( X(N+1,2) .GT. 10000) THEN
          Write(6,*) 'Memory Over in Operator Matrices.'
          CALL ERRSTP( 'POINTR' )
        END IF
      RETURN
      END
CCC
CC
C
      SUBROUTINE MATSTR(LV, NL,NR, TZ,NSB, XL,XR, SZ, MP)
      INTEGER       EZ, LV, NL,NR, TZ,NSB, XL(-1:101)
      INTEGER               MP(0:100,0:5), XR(-1:101), SZ(0:100,0:1)
C .......... Parameter Tansform (Fermion Reprisentation)
              EZ =(NL + NR + 2)/2 - TZ
C............ Ns = NL + NR + 2
      IF( EZ.LT.0  .OR.  EZ.GT.NL+NR+2 .OR. NL+NR+2.GT.100 ) THEN
          Write(6,*)'Ez,NL,NR = ',EZ,NL,NR,'  Something is wrong.'
          CALL ERRSTP( 'MATSTR' )
        END IF
C............................. Index: 0 --- Left, 1 --- Right.
      IF( EZ .LE. NR+1 ) THEN
          SZ(0,0)= 0
          SZ(0,1)= EZ
        ELSE
          SZ(0,0)= EZ-(NR+1)
          SZ(0,1)=     NR+1
        END IF
C............................. Sub-Space Scan
      NSB = 1
   10 CONTINUE
      IF( SZ(NSB-1,0) .LT. NL+1 .AND. SZ(NSB-1,1) .GT. 0 ) THEN
          SZ(NSB,0)= SZ(NSB-1,0) +1
          SZ(NSB,1)= SZ(NSB-1,1) -1
             NSB   =    NSB+1
         GOTO 10
        END IF
C     ........................ Matrix Structure Pointers ...
      DO 20   I = 0, NSB-1
         MP(I,0)= XL(SZ(I,0)-1) + XL(SZ(I,0))
         MP(I,1)= XR(SZ(I,1)-1) + XR(SZ(I,1))
         MP(I,2)=    MP(I,0)    *    MP(I,1)
   20 CONTINUE
         MP(0,3)= 0
         MP(0,4)= 0
         MP(0,5)= 0
      DO 30   I = 1, NSB
         MP(I,3)= MP(I-1,3) + MP(I-1,2)
         MP(I,4)= MP(I-1,4) + MP(I-1,0)**2
         MP(I,5)= MP(I-1,5) + MP(I-1,1)**2
   30 CONTINUE
C
      IF( LV .GE. 3 ) THEN
          Write(6,*) ' Szl, Szr, Dim[L], Dim[R], Dim[LR], Pointers:'
          DO 40 I = 0, NSB
          IF( MP(I,2) .NE. 0 )
     &    WRITE(6,100) SZ(I,0), SZ(I,1), MP(I,0), MP(I,1),
     &                 MP(I,2), MP(I,3), MP(I,4), MP(I,5)
   40     CONTINUE
  100     FORMAT(1H ,8I8)
        END IF
C
      IF( MP(NSB,3) .GT. 100000 ) THEN
          Write(6,*) 'Memory over in Lanczos Vectors'
          CALL ERRSTP( 'MATSTR' )
        END IF
      IF( MAX(MP(NSB,4),MP(NSB,5)) .GT. 40000 ) THEN
          Write(6,*) 'Memory over in Block Hamiltonian.'
          CALL ERRSTP( 'MATSTR' )
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE RENRML(LV, N,NSB,SZ,MP,XO,XN,SO,SN, BH, H,MX,IP,
     &                                                  ER, FMAX, V)
      INTEGER  BA,EP,BP,LV, N,NSB,MX, MP(0:100,0:5), XO(-1:101,0:3)
      INTEGER  IP, J,  FLG,SZ(0:100), Y(-1:101,0:3), XN(-1:101,0:3)
      REAL*8       ER, BH(0:40000), SO(0:10000,0:5),  V(0:100000,0:3)
      REAL*8     FMAX,  H(0:10000), SN(0:10000,0:5), TP(0:40000)
C
      CALL EXINDX( LV, N, NSB, Y, SZ, MP(0,0) )
      CALL VEC0( Y(N+1,2), V(0,2) )
      CALL VEC0( Y(N+1,1), V(0,3) )
      DO 10 I = 0, NSB-1
      CALL GETDML( MP(I,0),MP(I,1), V(Y(SZ(I),2),2), V(   MP(I,3),1))
      CALL DIAGON( LV,     MP(I,0), V(Y(SZ(I),2),2), V(Y(SZ(I),1),3))
   10 CONTINUE
      CALL SELECT( LV,FLG, MX, Y(N+1,1), ER, V, IP, FMAX)
      CALL NEWIDX( LV,  N,  Y, XN(-1,0),     V(0,3) )
C
      DO 20 I = 0, N
            J =    I-1
           EP = XN(I,2)
           BP =  Y(I,2) + Y(I,0)*(Y(I,0)-XN(I,0))
           BA =  Y(I,2)
      CALL ROTATD(FLG,Y(I,0),XN(I,0),BH(BA), H(EP), V(BP,2),V(0,3))
      CALL CRESZE(    Y(I,0),XO(J,0),TP,             SO(XO(J,2),3))
      CALL ROTATD(FLG,Y(I,0),XN(I,0),TP,  SN(EP,0), V(BP,2),V(0,3))
      CALL CRESZJ(    Y(I,0),XO(J,0),TP )
      CALL ROTATD(FLG,Y(I,0),XN(I,0),TP,  SN(EP,3), V(BP,2),V(0,3))
      IF( I .EQ. N ) GOTO 20
           EP = XN(I,3)
           BA =  Y(I+1,2) + Y(I+1,0)*(Y(I+1,0)-XN(I+1,0))
      CALL CRESPE(    Y(I,0),XO(J,0),TP,                  SO(XO(J,3),4))
      CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,1),V(BP,2),V(BA,2),V(0,3))
      CALL CRESPJ(    Y(I,0),XO(J,0),TP )
      CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,4),V(BP,2),V(BA,2),V(0,3))
   20 CONTINUE
C
      CALL HCONJU( N, XN, SN(0,1), SN(0,2) )
      CALL HCONJU( N, XN, SN(0,4), SN(0,5) )
      IF( LV .GE. 4 ) CALL MDUMP( N, XN(-1,0), H, SN )
      RETURN
      END
CCC
      SUBROUTINE RENRMR(LV, N,NSB,SZ,MP,XO,XN,SO,SN, BH, H, MX,IP,
     &                                                  ER, FMAX, V)
      INTEGER  BA,EP,BP,LV, N,NSB,MX, MP(0:100,0:5), XO(-1:101,0:3)
      INTEGER  IP, J,  FLG,SZ(0:100), Y(-1:101,0:3), XN(-1:101,0:3)
      REAL*8       ER, BH(0:40000), SO(0:10000,0:5),  V(0:100000,0:3)
      REAL*8     FMAX,  H(0:10000), SN(0:10000,0:5), TP(0:40000)
C
      CALL EXINDX( LV, N, NSB, Y, SZ, MP(0,1) )
      CALL VEC0( Y(N+1,2), V(0,2) )
      CALL VEC0( Y(N+1,1), V(0,3) )
      DO 10 I = 0, NSB-1
      CALL GETDMR( MP(I,0),MP(I,1), V(Y(SZ(I),2),2), V(   MP(I,3),1))
      CALL DIAGON( LV,     MP(I,1), V(Y(SZ(I),2),2), V(Y(SZ(I),1),3))
   10 CONTINUE
      CALL SELECT( LV,FLG, MX, Y(N+1,1), ER, V, IP, FMAX )
      CALL NEWIDX( LV,  N,  Y, XN(-1,0),     V(0,3) )
C
      DO 20 I = 0, N
            J =    I-1
           EP = XN(I,2)
           BP =  Y(I,2) + Y(I,0)*(Y(I,0)-XN(I,0))
           BA =  Y(N+1,2) - Y(I+1,2)
      CALL ROTATD(FLG,Y(I,0),XN(I,0),BH(BA), H(EP), V(BP,2),V(0,3))
      CALL CRESZE(    Y(I,0),XO(J,0),TP,             SO(XO(J,2),3))
      CALL ROTATD(FLG,Y(I,0),XN(I,0),TP,  SN(EP,0), V(BP,2),V(0,3))
      CALL CRESZJ(    Y(I,0),XO(J,0),TP )
      CALL ROTATD(FLG,Y(I,0),XN(I,0),TP,  SN(EP,3), V(BP,2),V(0,3))
      IF( I .EQ. N ) GOTO 20
           EP = XN(I,3)
           BA =  Y(I+1,2) + Y(I+1,0)*(Y(I+1,0)-XN(I+1,0))
      CALL CRESPE(    Y(I,0),XO(J,0),TP,                  SO(XO(J,3),4))
      CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,1),V(BP,2),V(BA,2),V(0,3))
      CALL CRESPJ(    Y(I,0),XO(J,0),TP )
      CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,4),V(BP,2),V(BA,2),V(0,3))
   20 CONTINUE
C
      CALL HCONJU( N, XN, SN(0,1), SN(0,2) )
      CALL HCONJU( N, XN, SN(0,4), SN(0,5) )
      IF( LV .GE. 4 ) CALL MDUMP( N, XN(-1,0), H, SN )
      RETURN
      END
CCC
      SUBROUTINE EXINDX( LV, N, NSB, Y, SZ, MP )
      INTEGER MP(0:100), LV, N, NSB, Y(-1:101,0:3), SZ(0:100)
      DO 10 I =-1, N+1
        Y(I,0)= 0
   10 CONTINUE
      DO 20 I = 0, NSB-1
       Y(SZ(I),0)= MP(I)
   20 CONTINUE
      CALL POINTR( LV, N, Y )
      IF( Y(N+1,2) .GT. 100000 ) THEN
          Write(6,*) 'Vector Overflow!'
          CALL ERRSTP( 'RENRML' )
        END IF
      RETURN
      END
CCC
      SUBROUTINE GETDML( WL, WR, V2, V1 )
      INTEGER            WL, WR
      REAL*8                     V2(0:*), V1(0:*)
      IF( WL .GE. WR ) THEN
          DO 10 J = 0, WR-1
          DO 10 K = 0, WL-1
          DO 10 L = 0, WL-1
          V2(WL*K+L)= V2(WL*K+L) + V1(WR*K+J) * V1(WR*L+J)
   10     CONTINUE
        ELSE
          DO 20 K = 0, WL-1
          DO 20 L = 0, WL-1
          DO 20 J = 0, WR-1
          V2(WL*K+L)= V2(WL*K+L) + V1(WR*K+J) * V1(WR*L+J)
   20     CONTINUE
        END IF
      RETURN
      END
CCC
      SUBROUTINE GETDMR( WL, WR, V2, V1 )
      INTEGER            WL, WR
      REAL*8                     V2(0:*), V1(0:*)
      IF( WL .GE. WR ) THEN
          DO 10 J = 0, WR-1
          DO 10 K = 0, WR-1
          DO 10 L = 0, WL-1
          V2(WR*J+K)= V2(WR*J+K) + V1(WR*L+K) * V1(WR*L+J)
   10     CONTINUE
        ELSE
          DO 20 L = 0, WL-1
          DO 20 J = 0, WR-1
          DO 20 K = 0, WR-1
          V2(WR*J+K)= V2(WR*J+K) + V1(WR*L+K) * V1(WR*L+J)
   20     CONTINUE
        END IF
      RETURN
      END
CCC
      SUBROUTINE DIAGON( LV, N, A, E )
      INTEGER            LV, N, IERR
      REAL*8        W(0:10000), A(0:*), E(0:*)
C
      IF( N .EQ. 0 ) THEN
          RETURN
        END IF
      IF( N .EQ. 1 ) THEN
          E(0)= A(0)
          A(0)= 1.0D0
          IF( LV .GE. 3 ) THEN
              Write(6,*) '-------------------- Eigen Value', E(0)
          END IF
          RETURN
        END IF
C
C     +----------------------------------------------+
C     |      Diagonalization of the matrix A:        |
C     |      ASL routine: Please read NEC mannual    |
C     +----------------------------------------------+
C          -NEC-
      CALL DCSMAA( A, N, N, E, W, IERR )
      IF( IERR .NE. 0 ) THEN
          Write(6,*) 'ERROR in DCSMAA, Flag= ', IERR
          CALL ERRSTP( 'DIAGON' )
        END IF
      IF( LV .GE. 3 ) THEN
          DO 10 I = 0, N-1
   10         Write(6,*) '-------------------- Eigen Value', E(I)
        END IF
      RETURN
      END
CCC
C     *========================== CAUTION =========================*
C     |     Density Matrix Cut-off Routine (Subrouitne SELECT)     |
C     |     is NOT perfect now: Do check your result when the      |
C     |     ground state is collection of local states!!           |
C     *============================================================*
      SUBROUTINE SELECT( LV, FLG, MX, N, ER, V, IP,    FMAX )
      INTEGER            LV, FLG, MX, N,        IP
      REAL*8         TP,COR,      ER, V(0:100000,0:3), FMAX
C
      IF( N .LE. MX ) THEN
          FLG = 0
          IF( LV .GE. 2 ) Write(6,*) 'All States Preserved'
           IP = N
         FMAX = 0.0D0
          RETURN
        END IF
C
          FLG = 1
      CALL  VMOVE( N, 0, 3, V )
      DO 30 I =   0, N-2
           IP = I
         FMAX = V(I,0)
      DO 20 J = I+1, N-1
      IF(FMAX.LT.V(J,0) ) THEN
         FMAX =  V(J,0)
           IP =    J
        END IF
   20 CONTINUE
      IF( IP .NE. I ) THEN
               TP = V( I,0)
           V( I,0)= V(IP,0)
           V(IP,0)= TP
        END IF
   30 CONTINUE
C                  ------- Cut Off Ratio, COR --------
       COR     = 1.05D0
C                  ------- (Dummy for Debug)  --------
        IP     = MIN( MX-1,N-1 )
      V(IP+1,0)= V(IP,0)
   40 CONTINUE
      IF( IP .EQ. 0 ) THEN
          IF( LV .GE. 3 ) Write(6,*) 'Only One State Is Selected.'
          IF( LV .GE. 2 ) THEN
              WRITE(6,*) 'Your result might be incorrect,'
              WRITE(6,*) 'because of minor error in the '
              WRITE(6,*) 'subroutine SELECT. Do check your result!'
            END IF
          GOTO 45
        END IF
      IF( V(IP,0) .LT. ER ) THEN
          IP = IP-1
          GOTO 40
        END IF
      IF( V(IP,0)/V(IP+1,0) .LT. COR ) THEN
          IP = IP-1
          GOTO 40
        END IF
   45 CONTINUE
C
         FMAX = MAX( V(IP,0), ER )
      DO 50 I = 0, N-1
      IF( V(I,3) .LT. FMAX ) V(I,3)= -2.0D0
   50 CONTINUE
      IP = IP + 1
      IF( LV .GE. 2 ) WRITE(6,100) IP, FMAX
  100 FORMAT(1H ,'                                           ',
     &           'Kept;', I5,'  Min[DM]', E15.7 )
C
      RETURN
      END
CCC
      SUBROUTINE NEWIDX( LV, N, Y, X, V )
      INTEGER            LV, N, Y(-1:101,0:3), X(-1:101,0:3)
      REAL*8               CNT,  V(0:100000)
      DO 20 I = -1, N+1
          CNT =    Y(I,0)
      DO 10 J = 0, Y(I,0)-1
      IF( V(Y(I,1)+J) .LT. -1.0D0 ) CNT = CNT -1
   10 CONTINUE
       X(I,0)= CNT
   20 CONTINUE
      CALL POINTR( LV, N, X )
      IF( LV .GE. 2 ) THEN
          WRITE(6,100) (X(I,0),I=0,N)
  100     FORMAT(1H ,15I5)
        END IF
      RETURN
      END
CCC
      SUBROUTINE ROTATD( FLG, WO, WN,  GM, SM, OM, TM )
      INTEGER            FLG, WO, WN
      REAL*8         SM(0:*), GM(0:*), OM(0:*), TM(0:*)
      IF( FLG .EQ. 0 ) THEN
          DO 10 J = 0, WO-1
          DO 10 I = 0, WO-1
          SM(WO*J+I)= GM(WO*J+I)
   10     CONTINUE
        ELSE
C                      ------ TM(I,K) ++ GM(I,J)*OM(J,K)
          CALL   VEC0( WO*WN, TM )
          DO 20 K = 0, WN-1
          DO 20 J = 0, WO-1
          DO 20 I = 0, WO-1
          TM(WO*K+I)= TM(WO*K+I) + GM(WO*J+I)*OM(WO*K+J)
   20     CONTINUE
C                      ------ SM(I,K) ++ OM(J,I)*TM(J,K)
          CALL   VEC0( WN*WN, SM )
          DO 30 K = 0, WN-1
          DO 30 I = 0, WN-1
          DO 30 J = 0, WO-1
          SM(WN*K+I)= SM(WN*K+I) + OM(WO*I+J)*TM(WO*K+J)
   30     CONTINUE
        END IF
      RETURN
      END
CCC
      SUBROUTINE ROTATP( FLG, Y,      X,     GM, SM, OM, QM, TM )
      INTEGER            FLG, Y(0:1), X(0:1)
      REAL*8         GM(0:*),SM(0:*),OM(0:*),QM(0:*),TM(0:*)
      IF( FLG .EQ. 0 ) THEN
          DO 10 J = 0, Y(1)-1
          DO 10 I = 0, Y(0)-1
          SM(Y(0)*J+I)= GM(Y(0)*J+I)
   10     CONTINUE
        ELSE
C                      ------ TM(I,K) ++ GM(I,J)*QM(J,K)
          CALL   VEC0( Y(0)*X(1), TM )
          DO 20 K = 0, X(1)-1
          DO 20 J = 0, Y(1)-1
          DO 20 I = 0, Y(0)-1
          TM(Y(0)*K+I)= TM(Y(0)*K+I) + GM(Y(0)*J+I)*QM(Y(1)*K+J)
   20     CONTINUE
C                      ------ SM(I,K) ++ OM(J,I)*TM(J,K)
          CALL   VEC0( X(0)*X(1), SM )
          DO 30 K = 0, X(1)-1
          DO 30 I = 0, X(0)-1
          DO 30 J = 0, Y(0)-1
          SM(X(0)*K+I)= SM(X(0)*K+I) + OM(Y(0)*I+J)*TM(Y(0)*K+J)
   30     CONTINUE
        END IF
      RETURN
      END
CCC
      SUBROUTINE CRESZE( W, X, BH,      S    )
      INTEGER    MP, EP, W, X(0:1)
      REAL*8                   BH(0:*), S(0:*)
      CALL VEC0( W*W, BH )
      DO 10 K = 0, 1
           EP = K*X(0)*X(0)
           MP = K*(W  *X(0)+X(0))
      DO 10 J = 0, X(K)-1
      DO 10 I = 0, X(K)-1
      BH(MP+J*W+I)= S(EP+J*X(K)+I)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE CRESZJ( W, X, BH    )
      INTEGER            W, X(0:1)
      REAL*8                   BH(0:*)
      CALL VEC0( W*W, BH )
      DO 10 I = 0, X(0)-1
      BH((W+1)*I)=-1.0D0/2
   10 CONTINUE
      DO 20 I = X(0), X(0)+X(1)-1
      BH((W+1)*I)= 1.0D0/2
   20 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE CRESPE( Y,      X,      BH,      S    )
      INTEGER    MP, EP, Y(0:1), X(0:2)
      REAL*8                             BH(0:*), S(0:*)
      CALL  VEC0(  Y(0)*Y(1), BH )
      DO 10 K = 0, 1
           MP = K*(Y(0)*X(1)+X(0))
           EP = K*      X(1)*X(0)
      DO 10 J = 0, X(K+1)-1
      DO 10 I = 0, X(K  )-1
      BH(MP+J*Y(0)+I)= S(EP+J*X(K)+I)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE CRESPJ( Y,      X,      BH )
      INTEGER            Y(0:1), X(0:2)
      REAL*8                             BH(0:*)
      CALL VEC0( Y(0)*Y(1), BH )
      DO 10 I = 0, X(1)-1
      BH(X(0)+(Y(0)+1)*I)= 1.0D0
   10 CONTINUE
      RETURN
      END
CCC
CC
C
      SUBROUTINE BLOCKH( LV,NSB, SZ, MP, XL, XR, JH,HL,HR, SL,SR, BH)
      INTEGER  OL,DL,ZL, LV,NSB, SZ(0:100, 0:1), XL(-1:101,0:3)
      INTEGER  OR,DR,ZR,         MP(0:100, 0:5), XR(-1:101,0:3)
      REAL*8            JH(0:4),SL(0:10000,0:5),  HL(0:10000)
      REAL*8    BH(0:40000,0:1),SR(0:10000,0:5),  HR(0:10000)
C
      IF(   MAX( MP(NSB,4),MP(NSB,5)) .LT. 40000 ) THEN
          CALL VEC0( MP(NSB,4), BH(0,0) )
          CALL VEC0( MP(NSB,5), BH(0,1) )
        ELSE
          Write(6,*) 'Memory Over in Block Hamiltonian Stack'
          CALL ERRSTP( 'BLOCKH' )
        END IF
C
      DO 20 I =  0,  NSB-1
           ZL = SZ( I,0)-1
           ZR = SZ( I,1)-1
           DL = XL(ZL,2)
           OL = XL(ZL,3)
           DR = XR(ZR,2)
           OR = XR(ZR,3)
      IF( MP(I,0) .EQ. 0 ) GOTO 10
      CALL BLKHH( XL(ZL,0), HL(DL  ),               BH(MP(I,4),0) )
      CALL BLKSD( XL(ZL,0), SL(DL,0), SL(DL,3), JH, BH(MP(I,4),0) )
      CALL BLKSP( XL(ZL,0), SL(OL,1), SL(OL,4), JH, BH(MP(I,4),0) )
      CALL BLKSM( XL(ZL,0), SL(OL,2), SL(OL,5), JH, BH(MP(I,4),0) )
C     IF( LV .GE. 4 ) CALL  SMTDSP(I,  MP(I,0),     BH(MP(I,4),0) )
   10 CONTINUE
      IF( MP(I,1) .EQ. 0 ) GOTO 20
      CALL BLKHH( XR(ZR,0), HR(DR  ),               BH(MP(I,5),1) )
      CALL BLKSD( XR(ZR,0), SR(DR,0), SR(DR,3), JH, BH(MP(I,5),1) )
      CALL BLKSP( XR(ZR,0), SR(OR,1), SR(OR,4), JH, BH(MP(I,5),1) )
      CALL BLKSM( XR(ZR,0), SR(OR,2), SR(OR,5), JH, BH(MP(I,5),1) )
C     IF( LV .GE. 4 ) CALL  SMTDSP(I,  MP(I,1),     BH(MP(I,5),0) )
   20 CONTINUE
C
      IF( LV .GE. 3 ) Write(6,*) 'Block Hamiltonian is created.'
      RETURN
      END
CCC
      SUBROUTINE SMTDSP( L, N, A )
      INTEGER            L, N
      REAL*8                   A(0:*)
      WRITE(6,*) '--- Width:', L
      IF( N .GT. 12 ) RETURN
      DO 10 I = 0, N-1
       WRITE(6,100) ( A(I+N*J), J=0, N-1 )
  100 FORMAT(1H ,12F6.2)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE BLKHH( X,       H, BH   )
      INTEGER       AP, X(0:1),  W, EP
      REAL*8         H(0:*),     BH(0:*)
            W =      X(0)+X(1)
      DO 10 I = 0, 1
           AP = I*(W*X(0)+X(0))
           EP = I*   X(0)*X(0)
      DO 10 J = 0, X(I)-1
      DO 10 K = 0, X(I)-1
      BH(AP+W*J+K)= BH(AP+W*J+K) + H(EP+X(I)*J+K)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE BLKSD( X, S0, S1, JH, BH )
      INTEGER      X(0:1), EP, AP,  W
      REAL*8 F0,F1,S0(0:*), S1(0:*), JH(0:3), BH(0:*)
C     ---------------- End Point:S0, Joint Point:S1
            W =  X(0)+X(1)
      DO 10 I = 0, 1
C     ----------------  N.N.Int.:F0,  N.N.N.Int.:F1
           F0 = JH(0)*JH(2)*(I-0.50D0)
           F1 = JH(1)*JH(2)*(I-0.50D0)
           AP = I*(W*X(0)+X(0))
           EP = I*   X(0)*X(0)
      DO 10 J = 0, X(I)-1
      DO 10 K = 0, X(I)-1
      BH(AP+W*J+K)= BH(AP+W*J+K) +F0*S1(EP+X(I)*J+K) +F1*S0(EP+X(I)*J+K)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE BLKSP(X, S0, S1, JH, BH )
      INTEGER     X(0:1),  W, AP
      REAL*8  S0(0:*), S1(0:*), JH(0:4), BH(0:*), F0, F1
            W = X(0)+X(1)
           F0 = JH(0)*JH(3)/2.0D0
           F1 = JH(1)*JH(3)/2.0D0
           AP =    X(0)*W
      DO 10 J = 0, X(1)-1
      DO 10 K = 0, X(0)-1
      BH(AP+W*J+K)= BH(AP+W*J+K) +F0*S1(X(0)*J+K) +F1*S0(X(0)*J+K)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE BLKSM(X, S0, S1, JH, BH )
      INTEGER     X(0:1),  W, AP
      REAL*8  S0(0:*), S1(0:*), JH(0:4), BH(0:*), F0, F1
            W = X(0)+X(1)
           F0 = JH(0)*JH(3)/2.0D0
           F1 = JH(1)*JH(3)/2.0D0
           AP =    X(0)
      DO 10 J = 0, X(0)-1
      DO 10 K = 0, X(1)-1
      BH(AP+W*J+K)= BH(AP+W*J+K) +F0*S1(X(1)*J+K) +F1*S0(X(1)*J+K)
   10 CONTINUE
      RETURN
      END
CCC
CC
C
      SUBROUTINE MULT(LV, ID, IS, NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V)
      INTEGER     BA, LV, ID, IS, NSB, SZ(0:100,0:1), XL(-1:101,0:3)
      INTEGER     YR, ZL, ZR, WL,  WR, MP(0:100,0:5), XR(-1:101,0:3)
      REAL*8             JH(0:3),    BH(0:40000,0:1),V(0:100000,0:3)
      REAL*8     SL(0:10000,0:5),    SR(0:10000,0:5)
C
      CALL VEC0( MP(NSB,3), V(0,ID) )
      DO 20 I = 0, NSB-1
      WL = MP(I,0)
      WR = MP(I,1)
      IF(  MP(I,2) .EQ. 0 ) GOTO 10
      BA = MP(I,3)
      ZL = SZ(I,0)-1
      ZR = SZ(I,1)-1
      CALL MULBLK(ID,IS,BA,WL,WR, BH(MP(I,4),0), BH(MP(I,5),1),     V )
      CALL MULMMD(ID,IS,BA,WR,XL(ZL,0),XR(ZR,0),                 JH,V )
      CALL MULMRD(ID,IS,BA,WR,XL(ZL,0),XR(ZR,0), SR(XR(ZR,2),3), JH,V )
      CALL MULLMD(ID,IS,BA,WR,XL(ZL,0),XR(ZR,0), SL(XL(ZL,2),3), JH,V )
      IF( I .EQ. NSB-1  ) GOTO 10
      YR = SZ(I,1)-2
      CALL MULMMO(ID,IS,BA,WR,XL(ZL,0),XR(YR,0),                 JH,V )
      CALL MULMRO(ID,IS,BA,WR,XL(ZL,0),XR(YR,0), SR(XR(YR,3),5), JH,V )
      CALL MULLMO(ID,IS,BA,WR,XL(ZL,0),XR(YR,0), SL(XL(ZL,3),4), JH,V )
   10 CONTINUE
   20 CONTINUE
C
      IF( LV .GE. 5 ) Write(6,*) 'End of Matrix Multiplication.'
      RETURN
      END
CCC
      SUBROUTINE MULBLK(ID,IS, BA, WL, WR,  BHL,   BHR, V )
      INTEGER           ID,IS, BA, WL, WR
      REAL*8      BHL(0:*),  BHR(0:*), V(0:100000,0:3)
C
      IF( WR .GT. WL ) THEN
          DO 10 I = 0, WL-1
          DO 10 J = 0, WL-1
          DO 10 K = 0, WR-1
          V(BA+J*WR+K,ID)= V(BA+J*WR+K,ID) +BHL(I*WL+J)*V(BA+I*WR+K,IS)
   10     CONTINUE
          DO 20 I = 0, WL-1
          DO 20 J = 0, WR-1
          DO 20 K = 0, WR-1
          V(BA+I*WR+K,ID)= V(BA+I*WR+K,ID) +BHR(J*WR+K)*V(BA+I*WR+J,IS)
   20     CONTINUE
        ELSE
          DO 30 K = 0, WR-1
          DO 30 I = 0, WL-1
          DO 30 J = 0, WL-1
          V(BA+J*WR+K,ID)= V(BA+J*WR+K,ID) +BHL(I*WL+J)*V(BA+I*WR+K,IS)
   30     CONTINUE
          DO 40 J = 0, WR-1
          DO 40 K = 0, WR-1
          DO 40 I = 0, WL-1
          V(BA+I*WR+K,ID)= V(BA+I*WR+K,ID) +BHR(J*WR+K)*V(BA+I*WR+J,IS)
   40     CONTINUE
        END IF
      RETURN
      END
CCC
      SUBROUTINE MULMMD(ID,IS, BA, WR, XL,  XR, JH, V )
      INTEGER       AP, ID,IS, BA, WR, XL(0:1), XR(0:1)
      REAL*8            FC,   JH(0:3),  V(0:100000,0:3)
      DO 10 I = 0, 1
      DO 10 J = 0, 1
           FC = JH(0)* JH(2)*(I-0.50D0)*(J-0.50D0)
           AP = I* WR* XL(0) + J*XR(0) + BA
      DO 10 K = 0, WR*(XL(I)-1), WR
      DO 10 L = 0,     XR(J)-1
      V(AP+K+L,ID)= V(AP+K+L,ID) + FC*V(AP+K+L,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULMRD(ID,IS, BA, WR, XL,  XR, SR, JH, V )
      INTEGER    AP, EP,ID,IS, BA, WR, XL(0:1), XR(0:1)
      REAL*8     FC,  JH(0:3), SR(0:*),  V(0:100000,0:3)
      DO 10 I = 0, 1
      DO 10 J = 0, 1
           AP = I *WR*XL(0) + J*XR(0) + BA
           EP = J    *XR(0)**2
      DO 10 K = 0,    XR(J)-1
      DO 10 L = 0,    XR(J)-1
           FC = JH(1)*JH(2)*(I-0.50D0) *SR(EP+L+K*XR(J))
      DO 10 M = 0, WR*(XL(I)-1), WR
      V(AP+M+L,ID)= V(AP+M+L,ID) + FC*V(AP+M+K,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULLMD(ID,IS, BA, WR, XL,      XR, SL, JH, V )
      INTEGER    AP, EP,ID,IS, BA, WR, XL(0:1), XR(0:1)
      REAL*8     FC,  JH(0:3), SL(0:*), V(0:100000,0:3)
      DO 10 I = 0, 1
      DO 10 J = 0, 1
           AP = I *WR*XL(0) + J*XR(0) + BA
           EP = I    *XL(0)**2
      DO 10 K = 0,    XL(I)-1
      DO 10 L = 0,    XL(I)-1
           FC = JH(1)*JH(2)*(J-0.50D0) *SL(EP+L+K*XL(I))
      DO 10 M = 0,    XR(J)-1
      V(AP+WR*L+M,ID)= V(AP+WR*L+M,ID) + FC*V(AP+WR*K+M,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULMMO( ID, IS, BA, WV, XL,  XR, JH, V )
      INTEGER   VP,HP,WH,ID, IS, BA, WV, XL(0:1), XR(-1:1)
      REAL*8       FC,  JH(0:3), V(0:100000,0:3)
           VP = BA + WV* XL(0)
           HP = BA + WV*(XL(0)+XL(1)) + XR(-1)
           WH =         XR(-1)+XR(0)
           FC =          JH(0)*JH(3)/2.0D0
      DO 10 I = 0, XL(1)-1
*VDIR NODEP
      DO 10 J = 0, XR(0)-1
      V(VP+WV*I+J,ID)= V(VP+WV*I+J,ID) + FC*V(HP+WH*I+J,IS)
      V(HP+WH*I+J,ID)= V(HP+WH*I+J,ID) + FC*V(VP+WV*I+J,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULMRO( ID, IS, BA, WV, XL,  XR, SR, JH, V )
      INTEGER   VP,HP,WH,ID, IS, BA, WV, XL(0:1), XR(-1:1)
      REAL*8       FC,  JH(0:3), SR(0:*),  V(0:100000,0:3)
           WH =          XR(-1)+XR(0)
      DO 10 I = 0, 1
           VP = BA + WV* XL(0)        + I*XR( 0)
           HP = BA + WV*(XL(0)+XL(1)) + I*XR(-1)
           EP =       I* XR(0)*XR(-1)
      DO 10 J = 0, XR(I-1)-1
      DO 10 K = 0, XR(I)-1
           FC = JH(1)*JH(3)*SR(EP+XR(I)*J+K)/2.0D0
*VDIR NODEP
      DO 10 L = 0, XL(1)-1
      V(VP+WV*L+K,ID)= V(VP+WV*L+K,ID) + FC*V(HP+WH*L+J,IS)
      V(HP+WH*L+J,ID)= V(HP+WH*L+J,ID) + FC*V(VP+WV*L+K,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULLMO( ID, IS, BA, WV, XL,  XR, SL, JH, V )
      INTEGER   VP,HP,WH,ID, IS, BA, WV, XL(0:2), XR(-1:1)
      REAL*8       FC,  JH(0:3), SL(0:*),  V(0:100000,0:3)
           WH =          XR(-1)+XR(0)
      DO 10 I = 0, 1
           VP = BA +                    I*WV*XL(0)
           HP = BA + WV*(XL(0)+XL(1)) + I*WH*XL(1) + XR(-1)
           EP =       I* XL(0)*XL(1)
      DO 10 J = 0, XL(I+1)-1
      DO 10 K = 0, XL(I)-1
           FC = JH(1)*JH(3)*SL(EP+XL(I)*J+K)/2.0D0
*VDIR NODEP
      DO 10 L = 0, XR(0)-1
      V(VP+WV*K+L,ID)= V(VP+WV*K+L,ID) + FC*V(HP+WH*J+L,IS)
      V(HP+WH*J+L,ID)= V(HP+WH*J+L,ID) + FC*V(VP+WV*K+L,IS)
   10 CONTINUE
      RETURN
      END
CCC
CC-------------------------------------------
C Subroutines For the Lanczos Diagonalization
C -------------------------------------------
      SUBROUTINE ERRSTP( A )
      CHARACTER*6        A
      Write(6,*) 'Software Termination in SUBROUTINE ', A, '.'
      Write(6,*) 'Do check Program/Parameters.'
      STOP
      END
CCC
      SUBROUTINE INPRO( NV, ID, IS, PRD, V)
      INTEGER           NV, ID, IS
      REAL*8                        PRD, V(0:100000,0:3)
          PRD = 0.0D0
      DO 10 I = 0, NV-1
          PRD = PRD + V(I,ID)*V(I,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE ADDER( NV, ID, IS1, IS2, A, B, V)
      INTEGER           NV, ID, IS1, IS2
      REAL*8                              A, B, V(0:100000,0:3)
      DO 10 I  = 0, NV-1
        V(I,ID)= A*V(I,IS1) + B*V(I,IS2)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE DIVER( NV, ID, IS, D, V)
      INTEGER           NV, ID, IS
      REAL*8                        D, V(0:100000,0:3)
      DO 10 I  = 0, NV-1
        V(I,ID)= V(I,IS) / D
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE VEC0( NV, V)
      INTEGER          NV
      REAL*8               V(0:100000)
      DO 10 J = 0, NV-1
          V(J)= 0.0D0
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE VMOVE( NV, ID, IS, V )
      INTEGER           NV, ID, IS
      REAL*8                        V(0:100000,0:3)
      DO 10 I  = 0, NV-1
        V(I,ID)= V(I,IS)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE VECINI( NV, V)
      INTEGER            NV
      REAL*8            PRD, V(0:100000,0:3)
      CALL  VEC0( NV, V(0,0) )
      CALL  VEC0( NV, V(0,1) )
      CALL  VEC0( NV, V(0,2) )
      CALL  VEC0( NV, V(0,3) )
      DO 10 I = 0,NV-1
        V(I,0)=( DBLE(IRAND()) + 00.1 )/SQRT( DBLE(NV) )
   10 CONTINUE
      CALL INPRO( NV, 0, 0, PRD, V )
          PRD = SQRT(PRD)
      CALL DIVER( NV, 0, 0, PRD, V )
      RETURN
      END
CCC
CC
C
      SUBROUTINE LANZOS(LV, NSB,LLS,SZ, MP, XL,XR, JH, SL,SR, BH, E, V)
      INTEGER LS,MLS,LK,LV, NSB, NV,    MP(0:100,0:5)
      INTEGER NQ,LLS,SZ(*), XL(*), XR(*)
      REAL*8      JH(*), BH(*), SL(*), SR(*), E(0:3)
      REAL*8   ENG, E2 , ALP(0:999), TE(0:999),V(0:100000,0:3)
      REAL*8   EPS, E2P, BET(0:999), ME(0:999,0:4), QE(0:2999)
C
      IF( LV.GE.3 ) WRITE(6,*) '-vvvv- Sub. LANcZOS -vvvvv-'
C     +----------------------------------------+
C     | LLS - Minimum Lanczos Step MLS       |
C     | MLS - Maximum Lanczos Step  999       |
C     | EPS - Eigenvalue Convergence Parameter |
C     +----------------------------------------+
       NV = MP(NSB,3)
      MLS = 190
      EPS = 1.0E-10
C                          ++++++ Initialization Steps ++++++
       E2 = 0.0D0
      CALL VECINI( NV,    V )
      CALL  VMOVE( NV,3,0,V )
      CALL   MULT(LV,1,0,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V)
      IF( LV .GE. 3 ) WRITE(6,*)
      IF( LV .GE. 3 ) WRITE(6,*) '-+- Eigenvalues -+-'
C              ++++++ LS = Lanczos Step, LK = 5*(LS/5) ++++++
       LK = 0
       LS = 0
   10 CONTINUE
      E2P = E2
C                       ++++++ Five Laanczos Steps ++++++
      DO 20 J = 0, 4
      CALL MACRO1(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,V,TE,ALP, BET,LK,0 )
   20 CONTINUE
C
C     ++++++ Diagonalization of the Tri-Diagonal matrix ++++++
      CALL MALANZ(  LK, LV, ALP,BET,ENG, E2, ME)
      LS = LS + 1
C                         ++++++ Convergence Check ++++++
      IF( LS*5 .LE. LLS                              ) GOTO 10
      IF( LS*5 .LE. MLS  .AND.  ABS(E2P-E2) .GE. EPS ) GOTO 10
      LS = LS - 1
      IF( LV .GE. 2 ) WRITE(6,100) ENG
  100 FORMAT(1H ,' Energy Upper Bound by Lanczos = ',E30.20)
C                     ++++++ Lanczos Vector Creation ++++++
      DO 30 I = 0,999
        TE(I) = ME(I,0)
   30 CONTINUE
C                         ++++++ Two-pass Method ++++++
      CALL VMOVE( NV, 0, 3,                    V )
      CALL ADDER( NV, 3, 3,  0,  TE(1), 0.0D0, V )
      CALL   MULT(LV,1,0,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V)
            LK = 0
C     +--------------------------------+
C     | NQ - Stacked-Eigenvalue Number |
C     | QE - Stacked-Eigenvalue        |
C     +--------------------------------+
            NQ = 0
      DO  50 I = 0, LS
      DO  40 J = 0,  4
      CALL MACRO1(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,V,TE,ALP, BET,LK,1 )
   40 CONTINUE
C
C     ++++++ Stack Low-Lying Eigenvalues ++++++
      CALL MALAN2(  LK, NQ, ALP, BET, QE(NQ) )
      IF( LV .GE. 5 ) WRITE(6,*) LK, ' -th Second Pass'
      IF( NQ .GE. 2999 ) THEN
          WRITE(6,*) 'Eigenvalue Stack Over Flow.'
          CALL ERRSTP( 'LANCZOS' )
        END IF
   50 CONTINUE
C
C     ++++++ Sort/Search of Low-Lying states ++++++
      CALL EVSORT( NQ, QE )
      CALL SEARCH( NQ, QE )
      DO 60 I = 0, 3
   60     E(I)= QE(I)
C
C     ++++++ Inprove the Lanczos Vector via Inverse Iteration. ++++++
      CALL INVERS(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,ENG,V)
C     ++++++ Inproved Ground State Energy ++++++
      E(0)= ENG
      IF( LV .GE. 2 ) WRITE(6,200) ENG
  200 FORMAT(1H ,' Ground State Energy by I.I. / C.G.=', F20.15)
C
      RETURN
      END
CCC
      SUBROUTINE MACRO1(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,
     &                                         BH,V,TE,ALP,BET,LK,ISW)
      INTEGER  NV,LK,ISW,LV,NSB,  MP(0:100,0:5)
      INTEGER  XL(*), XR(*), SZ(*)
      REAL*8  PROD,  ALP(0:999), BET(0:999), TE(0:999)
      REAL*8 SL(*), SR(*), BH(*), JH(*), V(0:100000,0:3)
C
      NV = MP(NSB,3)
C     ++++++ A Lanczos Step ++++++
      CALL INPRO( NV, 0, 1,           PROD,  V )
      ALP(LK) = PROD
      PROD    = -1.0 * PROD
      CALL ADDER( NV, 2, 1, 0, 1.0D0, PROD,  V )
      CALL INPRO( NV, 2, 2,           PROD,  V )
      PROD    = SQRT(PROD)
      BET(LK) =      PROD
      CALL DIVER( NV, 2, 2,           PROD,  V )
      CALL  MULT(LV,1,2,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V)
      PROD = -1.0 * PROD
      CALL ADDER( NV, 1, 1, 0, 1.0D0, PROD,  V )
      CALL VMOVE( NV, 0, 2,                  V )
C     ++++++ Second Path: Lanczos Vector Creation ++++++
      IF( ISW .EQ. 1 ) CALL ADDER( NV, 3, 3, 0, 1.0D0, TE(LK), V)
      LK = LK + 1
      RETURN
      END
CCC
      SUBROUTINE MALANZ( LK,  LV, ALP, BET, ENG,  E2,  VE )
      INTEGER      N, M, LK,  LV, LNV, ISW, IW1(0:999),IERR
      REAL*8  ENG,E2, ALP(0:999), VE(0:999,0:4)
      REAL*8  EPS,    BET(0:999),  E(0:999),   W1(0:5999)
C
        N = LK
      EPS = -1.0
        M = 5
      LNV = 1000
      ISW = -1
C
C     +----------------------------------------------------+
C     |      Diagonalization of the tri-diagonal matrix:   |
C     |      ASL routine: Please read NEC mannual.         |
C     +----------------------------------------------------+
C          -NEC-
      CALL DCSTSS( ALP,N,BET,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR )
      IF( IERR .NE. 0 ) THEN
          WRITE(6,100)  IERR
  100     FORMAT(1H ,' ERROR NO. in DCSTSS(ASL)',I10)
        END IF
      ENG = E(0)
      E2  = E(1)
      IF( LV .GE. 2 ) WRITE(6,200) (E(J),J = 0,2), LK
  200 FORMAT(1H ,'***',3F20.15,'*** ',I3,'- Steps')
      RETURN
      END
CCC
      SUBROUTINE MALAN2( LK, NQ, ALP,  BET,  QE )
      INTEGER            LK, NQ, IER
      REAL*8     TP, ALP(*), BET(*), QE(*), W1(0:5999)
C          -NEC-(Save Bet[n] because it is not conserved, actually.)-
      TP = BET(LK)
C
C     +----------------------------------------------------+
C     |      Diagonalization of the tri-diagonal matrix:   |
C     |      ASL routine: Please read NEC mannual.         |
C     +----------------------------------------------------+
C          -NEC-
      CALL DCSTSN( ALP, LK, BET, -1.0D0, QE, LK/5 + 4, -1, W1, IER )
      BET(LK) = TP
      NQ = NQ + LK/5 + 4
      RETURN
      END
CCC
      SUBROUTINE EVSORT( NQ, QE )
      INTEGER        IP, NQ
      REAL*8        MIN, TP, QE(0:*)
      DO 20 I = 0, NQ-2
           IP =    I
          MIN = QE(I)
      DO 10 J = I+1, NQ-1
      IF( MIN .GT. QE(J) ) THEN
          MIN  =   QE(J)
           IP  =      J
        END IF
   10 CONTINUE
      IF( I .NE. IP ) THEN
          TP     = QE(I)
          QE(I)  = QE(IP)
          QE(IP) = TP
        END IF
   20 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE SEARCH( NQ, QE )
      INTEGER       MIP, NQ, NE, SEP
      REAL*8        EPS, DL, QE(0:*), E(0:100), DQ(0:100000)
C
C     ................... Eigenvalue Pick-Up Levels .........
          EPS = 1.0E-10
           DL = 1.0E-5
C     ................... Eigenvalue differential .........
      DO 10 I = 0, NQ- 2
         DQ(I)= QE(I+1)- QE(I)
   10 CONTINUE
C
       NE   = 1
       E(0) = QE(0)
       SEP  = 1
   20 CONTINUE
      IF( DQ(SEP) .LT. DL  .AND. SEP .LT. NQ-2 ) THEN
          SEP = SEP + 1
          GOTO 20
        END IF
   30 CONTINUE
      IF( DQ(SEP) .GT. EPS .AND. SEP .LT. NQ-2 ) THEN
          SEP = SEP + 1
          GOTO 30
        END IF
      IF( SEP .EQ. NQ-2 ) GOTO 50
      MIP = SEP
   40 IF( DQ(SEP) .LT. DL  .AND. SEP .LT. NQ-2 ) THEN
          SEP = SEP + 1
          IF( DQ(SEP) .LT. DQ(MIP) ) THEN
              MIP = SEP
            END IF
          GOTO 40
        END IF
      IF( SEP .LT. NQ-2 ) THEN
          E(NE)=  QE(MIP)
            NE   = NE + 1
          IF( NE .LE. 100) GOTO 30
        END IF
   50 CONTINUE
C
      NQ = NE
      DO 60 I = 0, NE-1
         QE(I)= E(I)
   60 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE INVERS(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,LMD,V)
      INTEGER NV,MLR,LR,LV,NSB, MP(0:100,0:5), SZ(*)
      INTEGER XL(*), XR(*)
      REAL*8  BH(*), JH(*), SL(*), SR(*)
      REAL*8 DTR, LMD, DTP,EPS,  V(0:100000,0:3)
C
      IF( LV.GE.3 ) WRITE(6,*) '-vvvv- Sub. INVERS -vvvvv-'
C     +-----------------------------------------+
C     | MLR - Maximum Inverse Iteration #       |
C     | EPS - Eigenvector Convergence Parameter |
C     +-----------------------------------------+
       LR =  0
      MLR = 20
      EPS = 0.5E-11
       NV = MP(NSB,3)
C
      CALL VMOVE( NV,  0, 3,                          V )
      GOTO 20
   10 CONTINUE
      CALL CGMETH( LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,LMD,V )
   20 CALL  VMOVE( NV, 1, 0,                          V )
      CALL  INPRO( NV, 1, 1, DTP,                     V )
      DTP =  SQRT(DTP)
      CALL  DIVER( NV, 1, 1, DTP,                     V )
      CALL   MULT( LV, 2, 1,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V)
      CALL  INPRO( NV, 2, 2, DTP,                     V )
      DTP =  SQRT(DTP)
      CALL  DIVER( NV, 2, 2, DTP,                     V )
      IF( LMD .LT. 0 ) THEN
          CALL ADDER( NV, 3, 2, 1, 1.0D0,1.0D0,       V )
        ELSE
          CALL ADDER( NV, 3, 2, 1,-1.0D0,1.0D0,       V )
        END IF
      CALL  INPRO( NV, 3, 3, DTR,                     V )
      IF( LV .GE. 3 ) WRITE(6,100) LMD, DTP, SQRT(DTR)
  100 FORMAT(1H ,'Lmd=',E17.10,'  Lmd2=',E17.10,' Difference=',E17.10)
C      ++++++ Ground State Energy Inprovement ++++++
      IF( LMD .GT. 0 .AND. DTP .LT. LMD ) LMD = DTP
      IF( LMD .LT. 0 .AND. DTP .GT.-LMD ) LMD =-DTP
C
      IF( ABS(LMD) .LE. EPS ) THEN
          WRITE(6,*)'Eigen value ',LMD,' might be Accidentaly Zero.'
          WRITE(6,*)'Please change chemical potential to get normal',
     &              ' result.'
          IF( LV .GE. 3 ) THEN
              CALL ERRSTP( 'INVERS' )
            END IF
        END IF
C
      LR = LR + 1
      IF( LR .GE. MLR ) THEN
          WRITE(6,*)'Eigen vector is not obtained in ',MLR,' steps.'
          WRITE(6,*)'Eigen Value = ', LMD, ' in INVERS'
          IF( LV .GE. 1 ) CALL ERRSTP( 'INVERS' )
        END IF
C     ++++++ Convergence Check ++++++
      IF( SQRT(DTR) .GE. EPS ) GOTO 10
C
      RETURN
      END
CCC
      SUBROUTINE CGMETH(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,LMD,V)
      INTEGER   NV, LR, LV,NSB, MP(0:100,0:5), MLR
      INTEGER   XL(*), XR(*), SZ(*)
      REAL*8    BH(*), JH(*), SL(*), SR(*)
      REAL*8   DB, DR1, DR2, ALP, BET, EPS, DTP, LMD, V(0:100000,0:3)
C
C     +-----------------------------------------+
C     | MLR - Maximum Inverse Iteration #       |
C     | EPS - Eigenvector Convergence Parameter |
C     +-----------------------------------------+
       LR = 0
      MLR = 20
      EPS =1.0E-8
       NV = MP(NSB,3)
C (I)
      CALL INPRO( NV, 1, 1, DB, V)
      CALL VMOVE( NV, 2, 1,     V)
      DO 10 I = 0,NV-1
        V(I,0)= 0.0D0
   10 CONTINUE
C (II)
      CALL INPRO( NV, 1, 1, DR1,                    V)
   20 CALL  MULT( LV, 3, 2,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V)
      DO 30 I = 0, NV-1
        V(I,3)= V(I,3) - LMD*V(I,2)
   30 CONTINUE
      CALL INPRO( NV, 3, 2,DTP,               V)
      ALP=DR1/DTP
      CALL ADDER( NV, 0, 0, 2, 1.0D0, ALP,    V)
      CALL ADDER( NV, 1, 1, 3, 1.0D0,-ALP,    V)
      CALL INPRO( NV, 1, 1,      DR2,         V)
C
      IF( LV .GE. 5 ) WRITE(6,100) DR2 / DB, LR
  100 FORMAT(1H ,' CG Conv. Parameter = ', E18.10,' -',I3,'-Step')
      IF( DR2/DB .LE. EPS ) GOTO 40
      BET = DR2 / DR1
      DR1 = DR2
      LR  = LR + 1
      CALL ADDER( NV, 2, 1, 2,1.0D0, BET, V)
      IF( LR .GE. MLR ) THEN
          IF(LV.GE.4) WRITE(6,*)'Eigen vector is not obtained in ',
     &                MLR,' CG steps.'
          GOTO 50
        END IF
      GOTO 20
   40 CONTINUE
   50 IF( LV .GE. 5 ) WRITE(6,*) LR,' Rotations in CG'
C                
      RETURN
      END