Lanczos Diagonalization for Hubbard Ring
@

0010*#RUN: TIME=10 L=LIB/ASL7
0020C     ------- *  Dispersion Relation of the   * -------
0030C     ------- *  Extended Hubbard Hamiltonian * -------
0040C     ------- *                               * -------
0050C     ------- *  H = - 1*( CiCj + CjCi )      * -------
0060C     ------- *      + Ud* nn                 * -------
0070C     ------- *      + Un(  n+n  )*(  n+n  )  * -------
0080C     ------- *                               * -------
0090C     ------- *  10 Site 364364  >    63504   * -------
0100C     ------- *  12            >>>   853776   * -------
0110C     ------- *  14          >>>>> 11778624   * -------
0120C     ------- *                               * -------
0130C
0140      INTEGER NS, NK, NCF,   CF(0:12870         ), DG(0:65536,     0:1)
0150      INTEGER    NST(0:1),   ST(0:12870,0:20,0:1), WG(0:12870,     0:1)
0160      INTEGER     NE(0:1),   SG(0:12870,0:20,0:1),  S(0:12870,0:20,0:1)
0170      INTEGER               CON(0:12870,0:20,0:1)
0180      INTEGER LD, LU(0:20),  TD(0:99999,     0:2), TU(0:9999 ,0:20,0:2)
0190      REAL*8  UD, UN,   MU, TIME1,TIME2,  VEC(0:364364 ,0:3)
0200      REAL*8      BC(0:20),    BS(0:20),  NRM(0:364364), DIAG(0:364364)
0210C
0220      NS   =  4
0230C     NK   =  0
0240      NE(0)=  2 
0250      NE(1)=  2
0260C
0270      UD   =  12.0D0
0280      UN   =   6.0D0
0290      MU   =   2.0D0*UN + UD/2.0D0
0300C
0310      WRITE(6,346) UD, UN, NE(0), NE(1)
0320  346 FORMAT(1H ,' UD= ',F10.7,'  UN=',F10.7, ' Neu=',I3,' Ned=',I3 )
0330C
0335C     DO 90 I = 0   , NS/2
0340      DO 90 I = 0   , NS/2
0350           NK = I
0360C
0370      CALL   CLOCK( TIME1                                             )
0380      CALL   CLOCK( TIME2                                             )
0390      CALL  ORTSET( NS,  NK   , BC, BS                                )
0400      CALL  CRESTA( NS,  NE(0), NCF   , CF                            )
0410      CALL  REDUCE( NS, NCF,CF, NST(0), ST(0,0,0)                     )
0420      CALL  CRESTA( NS,  NE(1), NST(1), ST(0,0,1)                     )
0430      CALL   CREDG(             NST(0), ST(0,0,0),          DG(0,0)   )
0440      CALL   CREDG(             NST(1), ST(0,0,1),          DG(0,1)   )
0450      CALL  ROTATE( NS,  NE   , NST   , ST, SG,WG, S,  CON, DG        )
0460      CALL  NORMAL( NS,  BC   , NST   ,        WG, S,  NRM            )
0470      CALL  PRJOUT( NS,         NST   ,    CON,WG, S,  NRM            )
0480      CALL  DIAGON( NS,  NE   , NST   , ST, UD,UN,MU,  NRM, DIAG      )
0490      CALL  TRANSF( NS,  NE   , NST   , ST, SG,LD,LU,TD,TU, DG        )
0500      CALL  LANTZ1( NS,BS,BC,NK,NST, S,CON, SG,LD,LU,TD,TU,DIAG,NRM,VEC)
0510C-WR- CALL   NDIST( NS,  NK   , NST   , ST,                     VEC   )
0520      CALL  CDWORD( NS,  NK   , NST   , ST,                     VEC   )
0530   90 CONTINUE
0540      CALL   DTIME(      TIME1,  TIME2,                   ' IN ALL'   )
0550C
0560      STOP
0570      END
0580CCC
0590CC
0600C
0610      SUBROUTINE ORTSET( NS, NK, BC      , BS     )
0620      INTEGER            NS, NK
0630      REAL*8                 PI, BC(0:20), BS(0:20)
0640C
0650      IF( MOD(NS,2) .NE. 0 .OR. NK .GT. NS/2 ) THEN
0660          WRITE(6,*)' Wrong Site Number '
0670          STOP
0680        END IF
0690C
0700           PI = 4.0D0 * ATAN( 1.0D0 )
0710      DO 10 I = 0, NS -1
0720         BC(I)= COS( 2.0D0*I*NK*PI/NS )
0730         BS(I)= SIN( 2.0D0*I*NK*PI/NS )
0740C-WR- WRITE(6,100)  BC(I)        ,   BS(I)
0750C 100 FORMAT(1H , 'COS= ', F10.5, ' SIN=' F10.5 )
0760   10 CONTINUE
0770C
0780      RETURN
0790      END
0800CCC
0810CC
0820C
0830      SUBROUTINE CRESTA( NS, NE, NCONF, CONF )
0840      INTEGER            NS, NE, NCONF, CONF(0:12870)
0850C
0860        NCONF=0
0870      DO 20 I=0,2**NS-1
0880      IF( MOD(I       ,2) +MOD(I/ 2** 1,2) +MOD(I/ 2** 2,2)
0890     &   +MOD(I/ 2** 3,2) +MOD(I/ 2** 4,2) +MOD(I/ 2** 5,2)
0900     &   +MOD(I/ 2** 6,2) +MOD(I/ 2** 7,2) +MOD(I/ 2** 8,2)
0910     &   +MOD(I/ 2** 9,2) +MOD(I/ 2**10,2) +MOD(I/ 2**11,2)
0920     &   +MOD(I/ 2**12,2) +MOD(I/ 2**13,2) +MOD(I/ 2**14,2)
0930     &   +MOD(I/ 2**15,2) +MOD(I/ 2**16,2) +MOD(I/ 2**17,2)
0940     &   +MOD(I/ 2**18,2) +MOD(I/ 2**19,2)        .EQ.  NE ) THEN
0950      CONF(NCONF)= I
0960           NCONF = NCONF+1
0970        END IF
0980   20 CONTINUE
0990C-WR- WRITE(6,100) NS,NE,NCONF
1000C 100 FORMAT(1H ,' NS=',I2,' NE=',I3,' NCONF=',I5)
1010C
1020      RETURN
1030      END
1040CCC
1050CC
1060C
1070      SUBROUTINE REDUCE( NS, NCF, CF         , NST, ST        )
1080      INTEGER       FLG, NS, NCF, CF(0:12870), NST, ST(0:12870)
1090C
1100          NST = 0
1110      DO 20 I = 0, NCF-1
1120          FLG = 0
1130      DO 10 J = 1, NS-1
1140      IF( MOD(CF(I)*2**J,2**NS) + CF(I)/2**(NS-J) .LT. CF(I) ) FLG = 1
1150   10 CONTINUE
1160      IF( FLG .EQ. 0 ) THEN
1170          ST(NST)= CF(I)
1180             NST =  NST + 1
1190        END IF
1200   20 CONTINUE
1210C
1220C-WR- WRITE(6,100)NST
1230C 100 FORMAT(1H ,' Reduced Basis Set =', I5 )
1240C
1250      RETURN
1260      END
1270CCC
1280CC
1290C
1300      SUBROUTINE  CREDG( NST, ST         , DG        )
1310      INTEGER           NST, ST(0:12870), DG(0:65536)
1320C
1330      DO 10 I =  0, 65536
1340         DG(I)= -1
1350   10 CONTINUE
1360C
1370      DO 20 I  =  0, NST - 1
1380      DG(ST(I))=  I
1390   20 CONTINUE
1400C
1410      RETURN
1420      END
1430CCC
1440CC
1450C
1460      SUBROUTINE ROTATE( NS, NE, NST,  ST,  SG, WG, S, CON,  DG    )
1470      INTEGER        SP, NS, NE(0:1), NST(0:1), ST(0:12870,0:20,0:1)
1480      INTEGER       FCT,   S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1)
1490      INTEGER            CON(0:12870,0:20,0:1), WG(0:12870,     0:1)
1500      INTEGER                                   DG(0:65536,     0:1)
1510C
1520      DO 90 SP = 0, 1
1530           FCT = (-1)**( NE(SP)-1 )
1540C                                     --- Initialization (m=0) ---
1550      DO 10    I = 0,   NST(SP)-1
1560       SG(I,0,SP)= 1
1570       WG(I,  SP)= 1
1580        S(I,0,SP)= 1
1590      CON(I,0,SP)= I
1600   10 CONTINUE
1610C
1620      DO 40    J = 1,  NS    -1
1630      DO 20    I = 0, NST(SP)-1
1640       ST(I,J,SP)= MOD(ST(I,J-1,SP)*2,2**NS) + ST(I,J-1,SP)/2**(NS-1)
1650       SG(I,J,SP)=     SG(I,J-1,SP)*FCT**   (  ST(I,J-1,SP)/2**(NS-1) )
1660        S(I,J,SP)= 0
1670   20 CONTINUE
1680      DO 30    I = 0, NST(SP)-1
1690      CON(I,J,SP)= DG( ST(I,J,SP),SP )
1700      IF( ST(I,J,SP) .EQ. ST(I,0,SP) ) THEN
1710          WG(I,  SP)   =  WG(I,  SP)+1
1720           S(I,J,SP)   =  SG(I,J,SP)
1730        END IF
1740   30 CONTINUE
1750   40 CONTINUE
1760C
1770      DO 50  I = 0, NST(SP)-1
1780C-WR- WRITE(6,100) WG(I,SP),(SG(I,J,SP),J=0,NS-1),(S(I,J,SP),J=0,NS-1)
1790C 100 FORMAT(1H ,'W= ',I2, 32I2 )
1800   50 CONTINUE
1810C-WR- WRITE(6,*)
1820   90 CONTINUE
1830C
1840      RETURN
1850      END
1860CCC
1870CC
1880C
1890      SUBROUTINE NORMAL( NS, BC, NST, WG, S, NRM )
1900      INTEGER  NST(0:1), NS,     WG(0:12870,0:1), S(0:12870,0:20,0:1)
1910      REAL*8   TP, BC(0:20),    NRM(0:364364)
1920C
1930      DO 90 I = 0, NST(0)-1
1940      IF( WG(I,0) .EQ. 1 ) THEN
1950          DO 10 J = 0, NST(1)-1
1960          NRM( I*NST(1)+J ) = 1.0D0
1970   10     CONTINUE
1980       ELSE
1990          DO 40 J = 0, NST(1) -1
2000               TP = 0.0D0
2010          DO 30 K = 0, NS -1
2020               TP = TP + BC(K)*S(I,K,0)*S(J,K,1)
2030   30     CONTINUE
2040          IF( TP .GT. 1.0E-8 ) THEN
2050              NRM( I*NST(1)+J ) = 1.0D0 / SQRT(TP)
2060            ELSE
2070              NRM( I*NST(1)+J ) = 0.0D0
2080            END IF
2090   40     CONTINUE
2100        END IF
2110   90 CONTINUE
2120C
2130      RETURN
2140      END
2150CCC
2160CC
2170C
2180      SUBROUTINE PRJOUT( NS, NST, CON, WG, S, NRM )
2190      INTEGER    NS, NST(0:1), CON(0:12870,0:20,0:1),  WG(0:12870,0:1)
2200      INTEGER                    S(0:12870,0:20,0:1)
2210      REAL*8                   NRM(0:364364        )
2220C
2230      DO 90 I = 0, NST(0)-1
2240      IF( WG(I,0) .NE. 1 ) THEN
2250          DO 80 J = 0, NST(1)-1
2260          IF(   NRM( I*NST(1)+J ) .GT. 1.0E-8 ) THEN
2270              DO 70 K = 1, NS-1
2280              IF( S(I,K,0) .NE. 0 .AND. CON(J,K,1) .NE. J ) THEN
2290                  NRM( I*NST(1) + CON(J,K,1) )= 0.0D0
2300                END IF
2310   70         CONTINUE
2320            END IF
2330   80     CONTINUE
2340        END IF
2350   90 CONTINUE
2360C
2370C-WR- WRITE(6,100) ( NRM(I), I=0, NST(0)*NST(1)-1 )
2380C 100 FORMAT(1H ,'Normalization= ', F10.6 )
2390C
2400      RETURN
2410      END
2420CCC
2430CC
2440C
2450      SUBROUTINE DTIME(TIME1,TIME2,ASCII)
2460      REAL*8   TIME1,TIME2
2470      CHARACTER*6 ASCII
2480      IF(TIME1.GE.TIME2) THEN
2490          CALL CLOCK(TIME2)
2500          WRITE(6,100) TIME2-TIME1,ASCII
2510  100 FORMAT(1H ,F15.8,' SEC WAS PASSED IN ',A10)
2520        ELSE
2530          CALL CLOCK(TIME1)
2540          WRITE(6,100) TIME1-TIME2,ASCII
2550        END IF
2560      RETURN
2570      END
2580CCC
2590CC
2600C
2610      SUBROUTINE DIAGON( NS, NE,  NST, ST, UD, UN,  MU, NRM, DIAG )
2620      INTEGER   NS, NE(0:1), NST(0:1), ST(0:12870,0:20,0:1)
2630      INTEGER   T1, T2,  T3, CNT,  IS, ID
2640      REAL*8    UD, UN,  MU, NRM(0:364364), DIAG(0:364364)
2650C
2660      DO 20  I = 0, NST(0)-1
2670      DO 10  J = 0, NST(1)-1
2680            T1 =     ST(I,0,0)
2690            T2 =     ST(J,0,1)
2700           CNT =  I*NST(1)+J
2710      DIAG(CNT)= MU*( -NE(0) -NE(1) )
2720C
2730             T3= IAND(T1,T2)
2740      DIAG(CNT)= DIAG(CNT) + UD*(
2750     &         + MOD( T3        ,2)+ MOD( T3 / 2** 1,2)
2760     &         + MOD( T3 / 2** 2,2)+ MOD( T3 / 2** 3,2)
2770     &         + MOD( T3 / 2** 4,2)+ MOD( T3 / 2** 5,2)
2780     &         + MOD( T3 / 2** 6,2)+ MOD( T3 / 2** 7,2)
2790     &         + MOD( T3 / 2** 8,2)+ MOD( T3 / 2** 9,2)
2800     &         + MOD( T3 / 2**10,2)+ MOD( T3 / 2**11,2)
2810     &         + MOD( T3 / 2**12,2)+ MOD( T3 / 2**13,2)
2820     &         + MOD( T3 / 2**14,2)+ MOD( T3 / 2**15,2)
2830     &         + MOD( T3 / 2**16,2)+ MOD( T3 / 2**17,2)
2840     &         + MOD( T3 / 2**18,2)+ MOD( T3 / 2**19,2) )
2850C
2860   10 CONTINUE
2870   20 CONTINUE
2880C
2890      DO 50 IS = 0, NS-1
2900      IF(   IS .NE. NS-1 ) THEN
2910            ID = IS + 1
2920        ELSE
2930            ID = 0
2940        END IF
2950C
2960      DO 40  I = 0, NST(0)-1
2970      DO 30  J = 0, NST(1)-1
2980            T1 =     ST(I,0,0)
2990            T2 =     ST(J,0,1)
3000           CNT =  I*NST(1)+J
3010      DIAG(CNT)= DIAG(CNT)
3020     &         + UN*( MOD(T1/2**IS,2) + MOD(T2/2**IS,2) )
3030     &             *( MOD(T1/2**ID,2) + MOD(T2/2**ID,2) )
3040C
3050   30 CONTINUE
3060   40 CONTINUE
3070   50 CONTINUE
3080C
3090   DO 60 I = 0, NST(0)*NST(1)-1
3100      IF( NRM(I) .LT. 1.0E-8 ) DIAG(I)= 0.0D0
3110   60 CONTINUE
3120C
3130C-WR- WRITE(6,100) ( DIAG(I),I=0,NST(0)*NST(1)-1 )
3140C 100 FORMAT(1H ,'    Diagonal = ',F10.6 )
3150C
3160      RETURN
3170      END
3180CCC
3190CC
3200C
3210      SUBROUTINE TRANSF(NS,  NE, NST, ST,  SG, LD, LU, TD, TU, DG )
3220      INTEGER    ITP,   NS,  NE(0:1), NST(0:1),LD, LU(0:20)
3230      INTEGER   ST(0:12870,0:20,0:1),  SG(0:12870,0:20,0:1)
3240      INTEGER   TD(0:99999,0:2)     ,  TU(0: 9999,0:20,0:2)
3250      INTEGER   DG(0:65536,0:1), FL , MSK
3260C                                    ---- Down Spin Part ----
3270           LD = 0
3280      DO 20 I = 0, NS-1
3290      IF(   I .NE. NS-1 ) THEN
3300           FL =       1
3310          MSK = 3* 2**I
3320        ELSE
3330           FL = (-1)**( NE(1)-1 )
3340          MSK = 1+ 2**( NS   -1 )
3350        END IF
3360      DO 10 J = 0, NST(1)-1
3370          ITP = IEOR( MSK, ST(J,0,1) )
3380      IF( DG(ITP,1) .NE. -1 ) THEN
3390          TD(LD ,0) = J
3400          TD(LD ,1) = DG(ITP,1)
3410          TD(LD ,2) =-FL
3420             LD     = LD + 1
3430        END IF
3440   10 CONTINUE
3450   20 CONTINUE
3460C                                    ---- Up Spin Part ----
3470     DO 50 M = 0, NS-1
3480        LU(M)= 0
3490     DO 40 I = 0, NS-1
3500     IF(   I .NE. NS-1 ) THEN
3510          FL =       1
3520         MSK = 3* 2**I
3530       ELSE
3540          FL = (-1)**( NE(0)-1 )
3550         MSK = 1+ 2**( NS   -1 )
3560       END IF
3570     DO 30 J = 0, NST(0)-1
3580         ITP = IEOR( MSK, ST(J,M,0) )
3590     IF( DG(ITP,0) .NE. -1 ) THEN
3600         TU(LU(M),M,0) = J
3610         TU(LU(M),M,1) = DG(ITP,0)
3620         TU(LU(M),M,2) =-FL*SG(J,M,0)
3630            LU(M)      = LU(M) + 1
3640        END IF
3650   30 CONTINUE
3660   40 CONTINUE
3670   50 CONTINUE
3680C
3690C-WR- WRITE(6,100) LD
3700C-WR- WRITE(6,110)(LU(I),I=0,NS-1)
3710C 100 FORMAT(1H ,' Transfer Down=',  I6 )
3720C 110 FORMAT(1H ,' Transfer Up  =',10I6 )
3730      RETURN
3740      END
3750CCC
3760CC
3770C
3780      SUBROUTINE INPRO(  NVEC,  IDIS,  ISRC,  PROD, VEC )
3790      INTEGER            NVEC,  IDIS,  ISRC
3800      REAL*8                                  PROD, VEC(0:364364,0:3)
3810      PROD=0.0D0
3820      DO 10 I=0, NVEC-1
3830      PROD=PROD+VEC(I,IDIS)*VEC(I,ISRC)
3840   10 CONTINUE
3850      RETURN
3860      END
3870CCC
3880CC
3890C
3900      SUBROUTINE ADDER(  NVEC,  IDIS, ISRC1, ISRC2,    A,    B, VEC )
3910      INTEGER            NVEC,  IDIS, ISRC1, ISRC2
3920      REAL*8     VEC(0:364364,0:3),                    A,    B
3930      DO 10 I=0, NVEC-1
3940      VEC(I,IDIS)=A*VEC(I,ISRC1)+B*VEC(I,ISRC2)
3950   10 CONTINUE
3960      RETURN
3970      END
3980CCC
3990CC
4000C
4010      SUBROUTINE DIVER(  NVEC,  IDIS,  ISRC,     D, VEC)
4020      INTEGER            NVEC,  IDIS,  ISRC
4030      REAL*8                                     D, VEC(0:364364,0:3)
4040      DO 10 I=0, NVEC-1
4050      VEC(I,IDIS)=VEC(I,ISRC)/D
4060   10 CONTINUE
4070      RETURN
4080      END
4090CCC
4100CC
4110C
4120      SUBROUTINE VEC0(NVEC,VEC)
4130      INTEGER         NVEC
4140      REAL*8               VEC(0:364364,0:3)
4150      DO 10 I=0,3
4160      DO  5 J=0,NVEC -1
4170      VEC(J,I)=0.0D0
4180    5 CONTINUE
4190   10 CONTINUE
4200      RETURN
4210      END
4220CCC
4230CC
4240C
4250      SUBROUTINE VMOVE(  NVEC,  IDIS,  ISRC, VEC)
4260      INTEGER            NVEC,  IDIS,  ISRC
4270      REAL*8                                VEC(0:364364,0:3)
4280      DO 10 I=0, NVEC-1
4290      VEC(I,IDIS)=VEC(I,ISRC)
4300   10 CONTINUE
4310      RETURN
4320      END
4330CCC
4340CC
4350C
4360      SUBROUTINE VECINI( NS, NK, NST, NRM,  VEC )
4370      INTEGER            NS, NK, NST(0:1),   NV
4380      REAL*8      PROD, VEC(0:364364,0:3),  NRM(0:364364)
4390C
4400                NV = NST(0)*NST(1)
4410      IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN
4420          CALL VEC0(NV,VEC)
4430           DO 10 I =  0, NV - 1
4440          VEC(I,0) = (MOD(I,3)+MOD(I,5)+1)/SQRT(DBLE(NV))*NRM(I)
4450   10     CONTINUE
4460        ELSE
4470          CALL VEC0(NV*2,VEC)
4480           DO 20  I    =  0, NV-1
4490           VEC(   I,0) = (MOD(I,3)+MOD(I,5)+1)/SQRT(DBLE(NV))*NRM(I)
4500           VEC(NV+I,0) = (MOD(I,7)+MOD(I,5)+1)/SQRT(DBLE(NV))*NRM(I)
4510   20     CONTINUE
4520                    NV = NV * 2
4530        END IF
4540      CALL INPRO(NV,0,0,PROD,VEC)
4550      PROD=SQRT(PROD)
4560      CALL DIVER(NV,0,0,PROD,VEC)
4570      RETURN
4580      END
4590CCC
4600CC
4610C
4620      SUBROUTINE MULT( ID, IS, NS, NK, NST, LU, LD, TU, TD, S, SG, CON,
4630     &                                          BS, BC,NRM,  DIAG, VEC )
4640      INTEGER       ID, IS, NS,  NK, NST(0:1),LU(0:20), LD
4650      INTEGER   S(0:12870,0:20,0:1), CON(0:12870,0:20,0:1)
4660      INTEGER  SG(0:12870,0:20,0:1), TU(0:9999,0:20,0:2),TD(0:99999,0:2)
4670      REAL*8   BS(        0:20)    , BC(       0:20)
4680      REAL*8  VEC(0:364364,0:3)    ,NRM(0:364364)       , DIAG(0:364364)
4690C
4700      CALL  SCALE( ID, IS, NS, NK, NST,                       NRM, VEC )
4710      CALL MULTUP( ID, IS, NS, NK, NST, LU, TU,    CON, SG, BS,BC, VEC )
4720      CALL MULTDW( ID, IS, NS, NK, NST, LD, TD, S, CON, SG, BS,BC, VEC )
4730      CALL RESCAL( ID, IS, NS, NK, NST,                       NRM, VEC )
4740      CALL MULDIA( ID, IS, NS, NK, NST,                      DIAG, VEC )
4750C     WRITE(6,*)
4760C     WRITE(6,100)(VEC(I,IS),I=0,2*NST(0)*NST(1)-1)
4770C     WRITE(6,100)(VEC(I,ID),I=0,2*NST(0)*NST(1)-1)
4780C 100 FORMAT( 1H ,8F9.3 )
4790C
4800      RETURN
4810      END
4820CCC
4830CC
4840C
4850      SUBROUTINE SCALE( ID, IS, NS, NK, NST, NRM, VEC )
4860      INTEGER       NV, ID, IS, NS, NK, NST(0:1)
4870      REAL*8            NRM(0:364364) , VEC(0:364364, 0:3)
4880C
4890                NV = NST(0)*NST(1)
4900C
4910      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
4920          DO 10 I  = 0, NV - 1
4930          VEC(I,ID)= 0.0D0
4940          VEC(I,IS)= VEC(I,IS) * NRM(I)
4950   10     CONTINUE
4960        ELSE
4970          DO    20 I  = 0, NV - 1
4980          VEC(I   ,ID)= 0.0D0
4990          VEC(I+NV,ID)= 0.0D0
5000          VEC(I   ,IS)= VEC(I   ,IS) * NRM(I)
5010          VEC(I+NV,IS)= VEC(I+NV,IS) * NRM(I)
5020   20     CONTINUE
5030        END IF
5040C
5050      RETURN
5060      END
5070CCC
5080CC
5090C
5100      SUBROUTINE RESCAL( ID, IS, NS, NK, NST, NRM, VEC )
5110      INTEGER        NV, ID, IS, NS, NK, NST(0:1)
5120      REAL*8             NRM(0:364364) , VEC(0:364364, 0:3)
5130C
5140                NV = NST(0)*NST(1)
5150C
5160      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
5170          DO  10  I    = 0, NV -1
5180              VEC(I,ID)= VEC(I,ID) * NRM(I)
5190          IF( NRM(I) .GT. 1.0E-8 )  THEN
5200              VEC(I,IS)= VEC(I,IS) / NRM(I)
5210            END IF
5220   10     CONTINUE
5230        ELSE
5240          DO  20  I       = 0, NV -1
5250              VEC(I   ,ID)= VEC(I   ,ID) * NRM(I)
5260              VEC(I+NV,ID)= VEC(I+NV,ID) * NRM(I)
5270          IF( NRM(I) .GT. 1.0E-8 ) THEN
5280              VEC(I   ,IS)= VEC(I   ,IS) / NRM(I)
5290              VEC(I+NV,IS)= VEC(I+NV,IS) / NRM(I)
5300            END IF
5310   20     CONTINUE
5320        END IF
5330C
5340      RETURN
5350      END
5360CCC
5370CC
5380C
5390      SUBROUTINE MULDIA( ID, IS, NS, NK, NST, DIAG, VEC )
5400      INTEGER        NV, ID, IS, NS, NK, NST(0:1)
5410      REAL*8             DIAG(0:364364), VEC(0:364364,0:3)
5420C
5430              NV = NST(0) * NST(1)
5440C
5450      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
5460          DO 10 I    = 0, NV -1
5470            VEC(I,ID)= VEC(I,ID) + DIAG(I)*VEC(I,IS)
5480   10     CONTINUE
5490        ELSE
5500          DO 20 I       = 0, NV -1
5510            VEC(I   ,ID)= VEC(I   ,ID) + DIAG(I)*VEC(I   ,IS)
5520            VEC(I+NV,ID)= VEC(I+NV,ID) + DIAG(I)*VEC(I+NV,IS)
5530   20     CONTINUE
5540      END IF
5550C
5560      RETURN
5570      END
5580CCC
5590CC
5600C
5610      SUBROUTINE MULTUP(ID,IS, NS,NK, NST,    LU ,TU, CON,SG,BS,BC, VEC)
5620      INTEGER BX,BY, NV,ID,IS, NS,NK, NST(0:1)   ,TU(0:9999 , 0:20, 0:2)
5630      INTEGER  CON(0:12870,0:20,0:1),  SG(0:12870,0:20,0:1) ,   LU(0:20)
5640      REAL*8   VEC(0:364364,    0:3),  BS(0:20)  , BC(0:20) , FCT
5650C
5660                NV = NST(0)*NST(1)
5670      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
5680          DO 30  M = 0, NS    -1
5690          DO 20  I = 0, LU(M) -1
5700                BX =    TU(I,M,0) * NST(1)
5710                BY =    TU(I,M,1) * NST(1)
5720               FCT =    TU(I,M,2) *  BC(M)
5730*VDIR NODEP
5740          DO 10  L = 0, NST(1)-1
5750          VEC( BY+CON(L,M,1), ID )       = VEC( BY+CON(L,M,1), ID )
5760     &           + SG(L,M,1) *    FCT    * VEC( BX+L         , IS )
5770   10     CONTINUE
5780   20     CONTINUE
5790   30     CONTINUE
5800        ELSE
5810          DO 80  M = 0, NS    -1
5820          IF( ABS( BC(M) ) .GT. 1.0E-8 ) THEN
5830              DO 50  I = 0, LU(M) -1
5840                    BX =    TU(I,M,0) * NST(1)
5850                    BY =    TU(I,M,1) * NST(1)
5860                   FCT =    TU(I,M,2) *  BC(M)
5870*VDIR NODEP
5880              DO 40  L = 0, NST(1)-1
5890              VEC(   BY+CON(L,M,1),ID )     = VEC(   BY+CON(L,M,1),ID)
5900     &              +  SG(L,M,1) *    FCT   * VEC(   BX+L         ,IS)
5910              VEC(NV+BY+CON(L,M,1),ID )     = VEC(NV+BY+CON(L,M,1),ID)
5920     &              +  SG(L,M,1) *    FCT   * VEC(NV+BX+L         ,IS)
5930   40         CONTINUE
5940   50         CONTINUE
5950            END IF
5960          IF( ABS( BS(M) ) .GT. 1.0E-8 ) THEN
5970              DO 70  I = 0, LU(M) -1
5980                    BX =    TU(I,M,0) * NST(1)
5990                    BY =    TU(I,M,1) * NST(1)
6000                   FCT =    TU(I,M,2) *  BS(M)
6010*VDIR NODEP
6020              DO 60  L = 0, NST(1)-1
6030C
6040              VEC(   BY+CON(L,M,1),ID )     = VEC(   BY+CON(L,M,1),ID)
6050     &              +  SG(L,M,1) *    FCT   * VEC(NV+BX+L         ,IS)
6060              VEC(NV+BY+CON(L,M,1),ID )     = VEC(NV+BY+CON(L,M,1),ID)
6070     &              -  SG(L,M,1) *    FCT   * VEC(   BX+L         ,IS)
6080   60         CONTINUE
6090   70         CONTINUE
6100            END IF
6110   80     CONTINUE
6120        END IF
6130C
6140      RETURN
6150      END
6160CCC
6170CC
6180C
6190      SUBROUTINE MULTDW(ID,IS, NS,NK, NST, LD, TD, S, CON,SG,BS,BC, VEC)
6200      INTEGER BX,BY, NV,ID,IS, NS,NK, NST(0:1),LD,TD(0:99999, 0:2)
6210      INTEGER  CON(0:12870,0:20,0:1),  SG(0:12870,0:20,0:1)
6220      INTEGER                           S(0:12870,0:20,0:1)
6230      REAL*8   VEC(0:364364,    0:3),  BS(0:20)  , BC(0:20) , FCT , TP
6240C
6250      NV = NST(0)*NST(1)
6260      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
6270          DO 20 I = 0, LD-1
6280               BX =    TD(I,0)
6290               BY =    TD(I,1)
6300               TP =    TD(I,2)
6310*VDIR NODEP
6320          DO 10 J = 0, NST(0)-1
6330            VEC(J*NST(1)+BY,ID) = VEC(J*NST(1)+BY,ID)
6340     &                          + VEC(J*NST(1)+BX,IS) * TP
6350   10     CONTINUE
6360   20     CONTINUE
6370C
6380          DO 50 M = 1, NS-1
6390          DO 40 I = 0, NST(0)-1
6400          IF( S(I,M,0) .NE. 0 ) THEN
6410                  FCT = S(I,M,0) * BC(M)
6420*VDIR NOVECTOR
6430              DO 30 J = 0, LD-1
6440                   BX = CON( TD(J,0), NS-M, 1 )
6450                   BY =      TD(J,1)
6460                   TP =      TD(J,2)* FCT* SG(TD(J,0), NS-M, 1 )
6470              VEC(I*NST(1)+BY,ID) = VEC(I*NST(1)+BY,ID)
6480     &                            + VEC(I*NST(1)+BX,IS) * TP
6490   30         CONTINUE
6500            END IF
6510   40     CONTINUE
6520   50     CONTINUE
6530        ELSE
6540          DO 70 I = 0, LD-1
6550               BX =    TD(I,0)
6560               BY =    TD(I,1)
6570               TP =    TD(I,2)
6580*VDIR NODEP
6590          DO 60 J = 0, NST(0)-1
6600            VEC(   J*NST(1)+BY,ID) = VEC(   J*NST(1)+BY,ID)
6610     &                             + VEC(   J*NST(1)+BX,IS) * TP
6620            VEC(NV+J*NST(1)+BY,ID) = VEC(NV+J*NST(1)+BY,ID)
6630     &                             + VEC(NV+J*NST(1)+BX,IS) * TP
6640   60     CONTINUE
6650   70     CONTINUE
6660C
6670          DO 120 M = 1, NS-1
6680          IF( ABS( BC(M) ) .GT. 1.0E-8 ) THEN
6690              DO  90 I = 0, NST(0)-1
6700              IF( S(I,M,0) .NE. 0 ) THEN
6710                      FCT = BC(M) * S(I,M,0)
6720*VDIR NOVECTOR
6730                  DO 80 J = 0, LD-1
6740                       BX = CON( TD(J,0), NS-M, 1 )
6750                       BY =      TD(J,1)
6760                       TP =      TD(J,2)*   FCT  *SG(TD(J,0),NS-M,1)
6770                  VEC(   I*NST(1)+BY,ID)  = VEC(   I*NST(1)+BY,ID)
6780     &                           + TP      *VEC(   I*NST(1)+BX,IS)
6790                  VEC(NV+I*NST(1)+BY,ID)  = VEC(NV+I*NST(1)+BY,ID)
6800     &                           + TP      *VEC(NV+I*NST(1)+BX,IS)
6810   80             CONTINUE
6820                END IF
6830   90         CONTINUE
6840            END IF
6850          IF( ABS( BS(M) ) .GT. 1.0E-8 ) THEN
6860              DO 110 I = 0, NST(0)-1
6870              IF( S(I,M,0) .NE. 0 ) THEN
6880                      FCT = BS(M) * S(I,M,0)
6890*VDIR NOVECTOR
6900                 DO 100 J = 0, LD-1
6910                       BX = CON( TD(J,0), NS-M, 1 )
6920                       BY =      TD(J,1)
6930                       TP =      TD(J,2)*   FCT  *SG(TD(J,0),NS-M,1)
6940                  VEC(   I*NST(1)+BY,ID)  = VEC(   I*NST(1)+BY,ID)
6950     &                           + TP      *VEC(NV+I*NST(1)+BX,IS)
6960                  VEC(NV+I*NST(1)+BY,ID)  = VEC(NV+I*NST(1)+BY,ID)
6970     &                           - TP      *VEC(   I*NST(1)+BX,IS)
6980  100             CONTINUE
6990                END IF
7000  110         CONTINUE
7010            END IF
7020  120     CONTINUE
7030        END IF
7040C
7050      RETURN
7060      END
7070CCC
7080CC
7090C
7100      SUBROUTINE LANTZ1(NS,BS,BC,NK,NST,S,CON,SG,LD,LU, TD, TU,
7110     &                                            DIAG,NRM,VEC)
7120      INTEGER  NVEC,LK,LLANTZ,NS,NK,NST(0:1),    LD,LU(0:20)
7130      INTEGER    S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1)
7140      INTEGER  CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2)
7150      INTEGER                         TD(0:99999,     0:2)
7160      REAL*8    DIAG(0:364364), NRM(0:364364), VEC(0:364364,0:3)
7170      REAL*8 MEER(1000,10),TEER(1000),ENERGY,ALPHA(1000),BEETA(1000)
7180      REAL*8   BS(0:20)  ,  BC(0:20),E2,E2P,  ESP(1000)
7190C
7200        E2  = 0.0D0
7210C     LLANTZ=20
7220      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
7230          NVEC=NST(0)*NST(1)
7240        ELSE
7250          NVEC=NST(0)*NST(1)*2
7260        END IF
7270C
7280      IF( NVEC .GT. 364364 ) THEN
7290          WRITE(6,*) 'Too Large Matrix. Please Increase VEC Size '
7300          STOP
7310        END IF
7320C
7330      CALL   VEC0(  NVEC,             VEC )
7340      CALL VECINI( NS,NK,  NST,  NRM, VEC )
7350      CALL  VMOVE(  NVEC,  3  ,  0  , VEC )
7360C (I)
7370      CALL   MULT(1,0,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC)
7380C (II)
7390C     ********  LK IS AN INDEX FOR ALPHA AND BEETA  ********
7400      LK=1
7410      LLANTZ=0
7420   11 CONTINUE
7430C     ********  REPORT PER 5 ROTATION  ********
7440      E2P=E2
7450      DO 40 J=0,4
7460      CALL MACRO1(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC,
7470     &              TEER, ALPHA, BEETA,    LK,     0               )
7480   40 CONTINUE
7490      CALL MALANZ(LK,ALPHA,BEETA,ENERGY,E2,MEER,ESP)
7500      LLANTZ=LLANTZ+1
7510      IF(LLANTZ.LE.40.AND.ABS(E2P-E2).GE.1.0E-14) GOTO 11
7520      IF(LLANTZ.LE.10                           ) GOTO 11
7530      LLANTZ=LLANTZ-1
7540      WRITE(6,123) ( ESP(IK),IK=1,10 )
7550C-WR- WRITE(8,124) ( ESP(IK),IK=1,21 )
7560  123 FORMAT( 1H ,10F7.3 )
7570C 124 FORMAT( 1H ,7F10.5 )
7580C
7590      DO 88 I=1,1000
7600      TEER(I)=MEER(I,1)
7610   88 CONTINUE
7620      DO 99 I=0,2
7630      DO 98 J=0,NVEC -1
7640        VEC(J,I)=0.0D0
7650   98 CONTINUE
7660   99 CONTINUE
7670      CALL VMOVE(NVEC,  0  ,  3  , VEC)
7680      CALL ADDER(NVEC,  3  ,  3  ,  0  ,TEER(1), 0.0D0,VEC)
7690      CALL   MULT(1,0,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC)
7700      LK=1
7710      DO 150 I=0,LLANTZ
7720      DO 140 J=0,4
7730      CALL MACRO1(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC,
7740     &              TEER, ALPHA, BEETA,    LK,     1        )
7750  140 CONTINUE
7760  150 CONTINUE
7770      CALL INVERS(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,
7780     &                                             NRM,DIAG,ENERGY,VEC)
7790C-WR- WRITE(6,179) ENERGY
7800C 179 FORMAT(1H ,F20.15)
7810      RETURN
7820      END
7830CCC
7840CC
7850C
7860      SUBROUTINE MACRO1(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,
7870     &             DIAG, VEC   ,  TEER, ALPHA, BEETA,    LK,   ISW )
7880      INTEGER  NVEC,LK,ISW,NS,NK,NST(0:1),LD,LU(0:20)
7890      INTEGER   S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1)
7900      INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2)
7910      INTEGER                        TD(0:99999,     0:2)
7920      REAL*8  DIAG(0:364364),NRM(0:364364),VEC(0:364364,0:3)
7930      REAL*8  TEER(1000),ALPHA(1000),BEETA(1000),PROD,BS(0:20),BC(0:20)
7940C
7950      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
7960          NVEC=NST(0)*NST(1)
7970        ELSE
7980          NVEC=NST(0)*NST(1)*2
7990        END IF
8000C
8010      CALL INPRO(  NVEC,  0  ,  1  ,            PROD,VEC)
8020      ALPHA(LK)=PROD
8030      PROD=-1.0*PROD
8040      CALL ADDER(  NVEC,  2  ,  1  ,  0  ,1.0D0,PROD,VEC)
8050      CALL INPRO(  NVEC,  2  ,  2  ,            PROD,VEC)
8060      PROD=SQRT(PROD)
8070      BEETA(LK)=PROD
8080      CALL DIVER(  NVEC,  2  ,  2  ,            PROD,VEC)
8090      CALL  MULT(1,2,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC)
8100      PROD=-1.0*PROD
8110      CALL ADDER(  NVEC,  1  ,  1  ,  0  ,1.0D0,PROD,VEC)
8120      CALL VMOVE(  NVEC,  0  ,  2  ,                 VEC)
8130      LK=LK+1
8140      IF(ISW.EQ.1) CALL ADDER(  NVEC,3,3,0,1.0D0,TEER(LK),VEC)
8150      RETURN
8160      END
8170CCC
8180CC
8190C
8200      SUBROUTINE MALANZ(LK,ALPHA,BEETA,ENERGY,E2,VE,E)
8210      INTEGER LK,N,M,LNV,ISW,IW1(1000),IERR
8220      REAL*8  ALPHA(1000),BEETA(1000),ENERGY,E2
8230      REAL*8  EPS,E(1000),VE(1000,10),W1(6000)
8240        N = LK-1
8250      EPS = -1.0
8260        M = MIN(N,10)
8270      LNV = 1000
8280      ISW = -1
8290      CALL DCSTSS(ALPHA,N,BEETA,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR)
8300      IF(IERR.NE.0) THEN
8310          WRITE(6,100)IERR
8320  100     FORMAT(1H ,' ERROR NO. ',I10)
8330        END IF
8340      ENERGY=E(1)
8350      E2    =E(2)
8360C-WR- WRITE(6,150)(E(J),J=1,10)
8370C 150 FORMAT(1H ,'*',10F7.3 )
8380C-WR- IF(LK.GE. 65)WRITE(6,200)(VE(J,1),J=N-1,N)
8390C 200 FORMAT(2E30.20)
8400      RETURN
8410      END
8420CCC
8430CC
8440C
8450      SUBROUTINE  INVERS(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,
8460     &                                         DIAG,LAMBDA,VEC)
8470      INTEGER NVEC, LROT, NS,NK,NST(0:1),LU(0:20),LD
8480      INTEGER   S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1)
8490      INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2)
8500      INTEGER                        TD(0:99999,     0:2)
8510      REAL*8  DIAG(0:364364),NRM(0:364364),       DTP,DTR
8520      REAL*8  LAMBDA,VEC(0:364364,0:3), BS(0:20), BC(0:20)
8530C
8540      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
8550          NVEC=NST(0)*NST(1)
8560        ELSE
8570          NVEC=NST(0)*NST(1)*2
8580        END IF
8590C
8600      CALL VMOVE( NVEC,  0  ,  3  ,                         VEC)
8610      LROT=0
8620      GOTO 11
8630   10 CONTINUE
8640      CALL CGMETH(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,
8650     &                                          DIAG,LAMBDA,VEC)
8660   11 CALL VMOVE( NVEC,  1  ,  0  ,                         VEC)
8670      CALL INPRO( NVEC,  1  ,  1  ,DTP,                     VEC)
8680      DTP=SQRT(DTP)
8690      CALL DIVER( NVEC,  1  ,  1  ,DTP,                     VEC)
8700      CALL   MULT(2,1,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC)
8710      CALL INPRO( NVEC,  2  ,  2  ,DTP,                     VEC)
8720      DTP=SQRT(DTP)
8730      CALL DIVER( NVEC,  2  ,  2  ,DTP,                     VEC)
8740C
8750      IF(LAMBDA.LT.0) THEN
8760      CALL ADDER( NVEC,  3  ,  2  ,  1  ,1.0D0,1.0D0,       VEC)
8770        ELSE
8780      CALL ADDER( NVEC,  3  ,  2  ,  1 ,-1.0D0,1.0D0,       VEC)
8790        END IF
8800C
8810      CALL INPRO( NVEC,  3  ,  3  ,DTR,                     VEC)
8820C-WR- WRITE(6,100)LAMBDA,DTP,SQRT(DTR)
8830      LROT=LROT+1
8840C 100 FORMAT(1H ,'LAMBDA=',E20.10,'  ?  ',E20.10,'  D  ',E20.10)
8850      IF(LROT.GE.30) THEN
8860          WRITE(6,281)
8870  281     FORMAT(1H ,'DID NOT CONVERGE IN INVERS')
8880        STOP
8890        END IF
8900      IF(SQRT(DTR).GE.0.5E- 9) GOTO 10
8910C-WR- WRITE(6,345)( VEC(K,1),K=0,MIN(NVEC-1,35) )
8920C 345 FORMAT(1H , 5E15.7)
8930      RETURN
8940      END
8950CCC
8960CC
8970C
8980       SUBROUTINE  CGMETH(NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,
8990      &                                            DIAG,LAMBDA,VEC)
9000       INTEGER NVEC,LROT, NS,NK,NST(0:1),LU(0:20),LD
9010       INTEGER   S(0:12870,0:20,0:1), SG(0:12870,0:20,0:1)
9020      INTEGER CON(0:12870,0:20,0:1), TU(0:9999 ,0:20,0:2)
9030      INTEGER                        TD(0:99999,     0:2)
9040      REAL*8  DIAG(0:364364), NRM(0:364364)
9050      REAL*8  LAMBDA,VEC(0:364364,0:3), BS(0:20), BC(0:20)
9060      REAL*8  DB,DR1,DR2,ALPHA,BEETA,EPSRON,DTP
9070C
9080      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
9090          NVEC=NST(0)*NST(1)
9100        ELSE
9110          NVEC=NST(0)*NST(1)*2
9120        END IF
9130C (I)
9140      CALL INPRO( NVEC,  1  ,  1  , DB, VEC )
9150      CALL VMOVE( NVEC,  2  ,  1  ,     VEC )
9160      DO 10 I=0,NVEC-1
9170        VEC(I,0)=0.0
9180   10 CONTINUE
9190C (II)
9200      EPSRON=1.0E-8
9210      LROT=0
9220      CALL INPRO( NVEC,  1  ,  1  ,DR1,                    VEC)
9230    5 CALL   MULT(3,2,NS,NK,NST,LU,LD,TU,TD,S,SG,CON,BS,BC,NRM,DIAG,VEC)
9240      DO 30 I=0,NVEC-1
9250      VEC(I,3)=VEC(I,3)-LAMBDA*VEC(I,2)
9260   30 CONTINUE
9270      CALL INPRO( NVEC,  3  ,  2  ,DTP,                    VEC)
9280      ALPHA=DR1/DTP
9290      CALL ADDER( NVEC,  0  ,  0  ,  2  ,1.0D0, ALPHA,     VEC)
9300      CALL ADDER( NVEC,  1  ,  1  ,  3  ,1.0D0,-ALPHA,     VEC)
9310      CALL INPRO( NVEC,  1  ,  1  ,DR2,                    VEC)
9320C-WR- WRITE(6,999) DR2/DB
9330C 999 FORMAT(1H ,E30.20)
9340      IF(DR2/DB.LE.EPSRON) GOTO 200
9350      BEETA=DR2/DR1
9360      DR1=DR2
9370      LROT=LROT+1
9380      CALL ADDER( NVEC, 2 , 1 , 2 ,1.0D0, BEETA,VEC)
9390      IF(LROT.GE. 20) THEN
9400          WRITE(6,132)
9410  132     FORMAT(1H ,'DID NOT CONVERGE IN CG')
9420          GOTO 154
9430        END IF
9440      GOTO 5
9450  200 CONTINUE
9460  154 CONTINUE
9470      WRITE(6,678) LROT
9480  678 FORMAT(1H ,I10,'ROT IN CG')
9490      RETURN
9500      END
9510CCC
9520CC
9530C
9540      SUBROUTINE  NDIST( NS, NK, NST     , ST, VEC )
9550      INTEGER       ADD, NS, NK, NST(0:1), ST(0:12870,0:20,0:1), NV
9560      REAL*8    A(0:20),B(0:20), VEC(0:364364,0:3)
9570C
9580            NV= NST(0)*NST(1)
9590      DO 40 IS= 0,NS-1
9600      A(IS)   = 0
9610      B(IS)   = 0
9620      DO 30 I = 0,NST(0)-1
9630      DO 20 J = 0,NST(1)-1
9640      ADD     = I*NST(1)+J
9650      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
9660          A(IS)   = A(IS)+VEC(ADD,1)*VEC(ADD,1)*MOD(ST(I,0,0)/2**IS,2)
9670          B(IS)   = B(IS)+VEC(ADD,1)*VEC(ADD,1)*MOD(ST(J,0,1)/2**IS,2)
9680        ELSE
9690C
9700          A(IS)   = A(IS)+MOD(ST(I,0,0)/2**IS,2)
9710     &            *(VEC(ADD,1)*VEC(ADD,1)+VEC(NV+ADD,1)*VEC(NV+ADD,1))
9720          B(IS)   = B(IS)+MOD(ST(J,0,1)/2**IS,2)
9730     &            *(VEC(ADD,1)*VEC(ADD,1)+VEC(NV+ADD,1)*VEC(NV+ADD,1))
9740        END IF
9750C
9760   20 CONTINUE
9770   30 CONTINUE
9780   40 CONTINUE
9790      CALL DISP(NS, A,'UPHOLE')
9800      CALL DISP(NS, B,'DWHOLE')
9810      RETURN
9820      END
9830CCC
9840CC
9850C
9860      SUBROUTINE CDWORD( NS, NK, NST     , ST, VEC )
9870      INTEGER IS,ID,ADD, NS, NK, NST(0:1), ST(0:12870,0:20,0:1), NV
9880      REAL*8      S,CDW,         VEC(0:364364,0:3)
9890C
9900            NV= NST(0)*NST(1)
9910           CDW= 0.0D0
9920C
9930      DO 50 IS= 1,NS-1
9940      DO 40 ID= 0,IS-1
9950            S = (-1.0D0)**( IS+ID )
9960      DO 30 I = 0,NST(0)-1
9970      DO 20 J = 0,NST(1)-1
9980      ADD     = I*NST(1)+J
9990      IF( NK .EQ. 0  .OR.  NK .EQ. NS/2 ) THEN
10000          CDW = CDW + VEC(ADD,1)*VEC(ADD,1) * S
10010     &        *( MOD(ST(I,0,0)/2**IS,2) + MOD(ST(J,0,1)/2**IS,2) -1 )
10020     &        *( MOD(ST(I,0,0)/2**ID,2) + MOD(ST(J,0,1)/2**ID,2) -1 )
10030        ELSE
10040          CDW = CDW + (  VEC(   ADD,1)*VEC(   ADD,1)
10050     &                  +VEC(NV+ADD,1)*VEC(NV+ADD,1) ) * S
10060     &        *( MOD(ST(I,0,0)/2**IS,2) + MOD(ST(J,0,1)/2**IS,2) -1 )
10070     &        *( MOD(ST(I,0,0)/2**ID,2) + MOD(ST(J,0,1)/2**ID,2) -1 )
10080        END IF
10090   20 CONTINUE
10100   30 CONTINUE
10110   40 CONTINUE
10120   50 CONTINUE
10130C
10140      CDW = 2.0D0 * CDW / ( NS*(NS-1) )
10150      WRITE(6,100) CDW
10160  100 FORMAT(1H ,'CDW ORDER=',F12.5 )
10170C
10180      RETURN
10190      END
10200CCC
10210CC
10220C
10230      SUBROUTINE DISP(NS, A,CH)
10240      INTEGER         NS
10250      CHARACTER*6         CH
10260      REAL*8            A(0:20)
10270C
10280      WRITE(6, * ) CH
10290      WRITE(6,100)( A(I),I=0,NS-1 )
10300  100 FORMAT(1H ,8F9.6)
10310C
10320      RETURN
10330      END