C-----HLEMULT----HERMITIAN MATRICES-------------------------------------
C
C     CONTAINS SUBROUTINE LANCZS AND SAMPLE USPEC, CMATV
C     USED BY THE HERMITIAN VERSION OF THE LANCZOS ALGORITHMS
C
C     PORTABILITY:
C     THESE PROGRAMS ARE NOT PORTABLE DUE TO THE USE OF COMPLEX*16
C     VARIABLES.  MOREOVER, THE PFORT VERIFIER IDENTIFIED THE
C     FOLLOWING ADDITIONAL NONPORTABLE CONSTRUCTIONS:
C     1.  THE ENTRY MECHANISM USED TO PASS THE STORAGE
C         LOCATIONS OF THE USER-SPECIFIED MATRIX FROM THE
C         SUBROUTINE USPEC TO THE MATRIX-VECTOR SUBROUTINE CMATV.
C     2.  IN THE PROGRAMS PROVIDED FOR 'HERMITIAN POISSON' TEST MATRICES
C         USPEC CONTAINS FREE FORMAT (8,*), AND FORMAT (20A4); AND
C         EXACT ERROR SUBROUTINE CONTAINS DATA/MACHEP DEFINITION.
C
C
C-----LANCZS-COMPUTE THE LANCZOS TRIDIAGONAL MATRICES-------------------
C
C     GRAM-SCHMIDT ORTHOGONALIZATION WITHOUT MODIFICATION
C     REQUIRES EXTRA VECTOR VS IN LANCZS.  MODIFICATION IS NOT
C     PERMISSIBLE IN THE HERMITIAN CASE BECAUSE COMPLEX PORTION
C     OF THE MODIFICATION COULD NOT BE INCORPORATED.
C
      SUBROUTINE LANCZS(MATVEC,V1,V2,VS,ALPHA,BETA,GR,GC,G,KMAX,MOLD1,N,
     1 IIX)
C
C-----------------------------------------------------------------------
      COMPLEX*16 V1(1), V2(1), VS(1), ZEROC, TEMP
      DOUBLE PRECISION  ALPHA(1), BETA(1), BATA, SUM, ONE, ZERO
      DOUBLE PRECISION GR(1),GC(1)
      REAL  G(1)
      EXTERNAL MATVEC
      DOUBLE PRECISION  DSQRT
C-----------------------------------------------------------------------
C
      ZERO = 0.D0
      ONE = 1.D0
      ZEROC = DCMPLX(ZERO,ZERO)
C
      IF(MOLD1.GT.1)GO TO 50
C
C     ALPHA/BETA GENERATION STARTS AT I = 1
C     MOLD1 = 1 SET V1 = 0. AND V2 = RANDOM UNIT VECTOR
      IIL=IIX
C
C-----------------------------------------------------------------------
      CALL GENRAN(IIL,G,N)
C-----------------------------------------------------------------------
C
      DO 10 I = 1,N
   10 GR(I) = G(I)
C
C-----------------------------------------------------------------------
      CALL GENRAN(IIL,G,N)
C-----------------------------------------------------------------------
C
      DO 20 I = 1,N
   20 GC(I) = G(I)
C
      DO 30 I = 1,N
   30 V2(I) = DCMPLX(GR(I),GC(I))
C
C-----------------------------------------------------------------------
      CALL CINPRD(V2,V2,SUM,N)
C-----------------------------------------------------------------------
C
      SUM = ONE/DSQRT(SUM)
      DO 40 I = 1,N
      V1(I) = ZEROC
   40 V2(I) = V2(I)*SUM
      BETA(1) = ZERO
C
C     ALPHA BETA GENERATION LOOP
   50 CONTINUE
C
      DO 80 I=MOLD1,KMAX
      SUM = ZERO
C
C-----------------------------------------------------------------------
C     MATVEC(V2,VS,SUM) CALCULATES  VS = A*V2 - SUM*VS
      CALL MATVEC(V2,VS,SUM)
      CALL CINPRD(V2,VS,SUM,N)
C-----------------------------------------------------------------------
C
      ALPHA(I) = SUM
      BATA = BETA(I)
      DO 60 J=1,N
   60 V1(J) = (VS(J)-BATA*V1(J)) - SUM*V2(J)
C
C-----------------------------------------------------------------------
      CALL CINPRD(V1,V1,SUM,N)
C-----------------------------------------------------------------------
C
      IN = I+1
      BETA(IN) = DSQRT(SUM)
      SUM = ONE/BETA(IN)
      DO 70 J=1,N
      TEMP = SUM*V1(J)
      V1(J) = V2(J)
   70 V2(J) = TEMP
   80 CONTINUE
C     END ALPHA, BETA GENERATION LOOP
C
C-----END OF LANCZS-----------------------------------------------------
C
      RETURN
      END
C
C-----USPEC-GENERAL SPARSE, HERMITIAN MATRIX----------------------------
C
C     SUBROUTINE USPEC(N,MATNO)
      SUBROUTINE GUSPEC(N,MATNO)
C
C-----------------------------------------------------------------------
      COMPLEX*16 A(3000)
      DOUBLE PRECISION  AD(1000)
      INTEGER  IROW(3000),ICOL(1000)
C-----------------------------------------------------------------------
C    DIMENSION ARRAYS NEEDED TO DEFINE MATRIX, READ IN VALUES FOR
C    ARRAYS AND THEN PASS THE STORAGE LOCATIONS OF THESE ARRAYS TO
C    THE MATRIX-VECTOR MULTIPLY SUBROUTINE CMATV.
C
C    USER-SUPPLIED MATRIX IS STORED IN FOLLOWING SPARSE FORMAT:
C    N = ORDER OF A-MATRIX
C    NZS = NUMBER OF NONZERO SUBDIAGONAL ENTRIES IN A
C    NZL = INDEX OF LAST COLUMN CONTAINING NONZERO SUBDIAGONAL ENTRIES
C    ICOL(J), J=1,NZL IS THE NUMBER OF NONZERO SUBDIAGONAL ELEMENTS
C             IN COLUMN J.
C    IROW(K), K = 1,NZS, IS THE ROW INDEX FOR CORRESPONDING A(K).
C    AD(I), I=1,N ARE DIAGONAL ENTRIES (INCLUDING ANY 0 DIAGONAL
C           ENTRIES)
C    A(K), K=1,NZS ARE NONZERO SUBDIAGONAL ENTRIES, LISTED BY COLUMN.
C    FOR J > NZL THERE ARE NO NONZERO SUBDIAGONAL ELEMENTS IN COLUMN J.
C    ICOL(J) = 0 IS ALLOWED
C
C-----------------------------------------------------------------------
C     IN THIS SAMPLE SUBROUTINE THE ARRAYS ARE READ IN FROM FILE 8
C
      READ(8,10) NZS,NOLD,NZL,MATOLD
   10 FORMAT(I10,2I6,I8)
C
      WRITE(6,20) NZS,NOLD,NZL,MATOLD
   20 FORMAT(I10,2I6,I8,' = NZS,NOLD,NZL,MATOLD'/)
C
C     TEST OF PARAMETER CORRECTNESS
      ITEMP = (NOLD-N)**2 + (MATNO-MATOLD)**2
C
      IF(ITEMP.EQ.0) GO TO 40
C
      WRITE(6,30)
   30 FORMAT(/' PROGRAM TERMINATES BECAUSE EITHER ORDERS OF OR LABELS FO
     1R MATRIX DISAGREE'/)
      GO TO 80
C
   40 CONTINUE
C
C     NUMBER OF NONZERO SUBDIAGONAL ENTRIES IN EACH COLUMN IS READ
C     THEN THE CORRESPONDING ROW INDEX FOR EACH SUCH ENTRY IS READ
      READ(8,50) (ICOL(K), K=1,NZL)
      READ(8,50) (IROW(K), K=1,NZS)
   50 FORMAT(13I6)
C     DIAGONAL IS READ FIRST, THEN NONZERO BELOW DIAGONAL ENTRIES
      READ(8,60) (AD(K), K=1,N)
   60 FORMAT(4E20.12)
      READ(8,70) (A(K), K=1,NZS)
C  50 FORMAT(4Z20)
   70 FORMAT(4E20.12)
C
C-----------------------------------------------------------------------
C     PASS STORAGE LOCATIONS OF ARRAYS THAT DEFINE THE MATRIX TO
C     THE MATRIX-VECTOR MULTIPLY SUBROUTINE CMATV
      CALL CMATVE(A,AD,ICOL,IROW,N,NZL)
C-----------------------------------------------------------------------
C
      RETURN
   80 STOP
C
C-----END OF USPEC FOR GENERAL, SPARSE HERMITIAN MATRICES---------------
      END
C
C-----START OF MATRIX-VECTOR MULTIPLY-GENERAL SPARSE HERMITIAN---------
C
C     SUBROUTINE CMATV(W,U,SUM)
      SUBROUTINE GCMATV(W,U,SUM)
C
C-----------------------------------------------------------------------
      COMPLEX*16  U(1),W(1),A(1)
      DOUBLE PRECISION  AD(1),SUM
      INTEGER  IROW(1),ICOL(1)
C-----------------------------------------------------------------------
C    SPARSE MATRIX-VECTOR MULTIPLY FOR LANCZS  U = A*W - SUM*U
C    SEE USPEC SUBROUTINE FOR DESCRIPTION OF THE ARRAYS THAT DEFINE
C    THE MATRIX
C
C     COMPUTE THE DIAGONAL TERMS
      DO 10 I = 1,N
   10 U(I) = AD(I)*W(I)-SUM*U(I)
C
C     COMPUTE BY COLUMN
      LLAST = 0
      DO 30 J = 1,NZL
C
      IF (ICOL(J).EQ.0) GO TO 30
      LFIRST = LLAST + 1
      LLAST = LLAST + ICOL(J)
C
      DO 20 L = LFIRST,LLAST
      I = IROW(L)
C
      U(I) = U(I) + A(L)*W(J)
      U(J) = U(J) + DCONJG(A(L))*W(I)
C
   20 CONTINUE
C
   30 CONTINUE
C
      RETURN
C
C-----------------------------------------------------------------------
C     STORAGE LOCATIONS OF ARRAYS ARE PASSED TO CMATV FROM USPEC
      ENTRY CMATVE(A,AD,ICOL,IROW,N,NZL)
C-----------------------------------------------------------------------
C
C-----END OF CMATV-GENERAL, SPARSE, HERMITIAN MATRICES -----------------
      RETURN
      END
C
C-----USPEC, CMATV, EXEVG, AND HEXVEC FOR HERMITIAN 'POISSON' MATRICES--
C
C-----USPEC (HERMITIAN POISSON MATRICES)--------------------------------
C
C     SUBROUTINE HUSPEC(N,MATNO)
      SUBROUTINE USPEC(N,MATNO)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  C0,C1,C2,HALF,ONE,SCR,SCI,ANGLE,TEMP
      COMPLEX*16 SC,TC,CL0,CL1,CL3,CL4
      REAL EXPLAN(20)
      DOUBLE PRECISION EIGVAL(1000)
      REAL GAPS(1000)
      INTEGER MULTS(1000)
C-----------------------------------------------------------------------
      HALF = 0.5D0
      ONE  = 1.0D0
C
C     READ IN PARAMETERS TO DEFINE MATRIX
C     MATRIX IS COMPLEX DIAGONAL SIMILITARY TRANSFORM OF REAL SYMMETRIC
C     POISSON MATRIX WHICH HAS SYMMETRIC TOEPLITZ BLOCKS ALONG
C     THE DIAGONAL, EACH ONE OF WHICH HAS THE PARAMETER C2 ALONG THE
C     DIAGONAL AND -C0 ABOVE AND BELOW THE DIAGONAL, AND OFF-DIAGONAL
C     BLOCKS THAT ARE DIAGONAL WITH DIAGONAL ENTRIES -C1.  EACH BLOCK
C     IS KX*KX AND THERE ARE KY BLOCKS.  THE HERMITIAN VERSION IS
C     OBTAINED BY APPLYING A DIAGONAL SIMILARITY TRANSFORM TO THE
C     REAL MATRIX WHERE THIS TRANSFORMATION IS SUCH THAT ITS
C     DIAGONAL ENTRIES ARE (SC)**(K-1), K = 1,...,N, WHERE SC
C     HAS MODULUS 1.
C
      READ(8,10) EXPLAN
      WRITE(6,10) EXPLAN
      READ(8,10) EXPLAN
   10 FORMAT(20A4)
C     IF MTYPE = 0 WE HAVE ZERO BOUNDARY CONDITIONS
C     IF MTYPE = 1 WE HAVE NORMAL DERIVATIVE BOUNDARY CONDITIONS
C     NOTE THAT SUBROUTINES EXEVG AND HEXVEC ARE VALID ONLY FOR
C     MTYPE = 0.
      READ(8,*) NOLD,MATOLD,IVEC,MTYPE
      WRITE(6,20) NOLD,MATOLD
   20 FORMAT(' ORDER OF MATRIX READ FROM FILE =',I6/' MATRIX NUMBER =',
     1I8/)
      IF(MTYPE.EQ.0) WRITE(6,30)
   30 FORMAT(/' HERMITIAN POISSON CORRESPONDING TO ZERO BOUNDARY CONDITI
     1ONS'/)
      IF(MTYPE.EQ.1) WRITE(6,40)
   40 FORMAT(/' HERMITIAN POISSON CORRESPONDING TO NORMAL DERIVATIVE BOU
     1NDARY CONDITIONS'/)
      IF(IVEC.NE.0.AND.MTYPE.EQ.0) WRITE(6,50)
   50 FORMAT(' COMPUTE THE TRUE EIGENVALUES AND PUT IN FN TRUEEVAL'/)
C
C     TEST OF PARAMETER CORRECTNESS
      ITEMP = (NOLD-N)**2 + (MATNO-MATOLD)**2
C
      IF(ITEMP.EQ.0) GO TO 70
C
      WRITE(6,60)
   60 FORMAT(' PROGRAM TERMINATES BECAUSE EITHER ORDERS OF OR LABELS FOR
     1 MATRIX DISAGREE')
      GO TO 150
C
   70 CONTINUE
C
      READ(8,10) EXPLAN
      READ(8,*)  C0,KX,KY
      IF (KX.GT.4.AND.KY.GT.4) GO TO 90
      WRITE(6,80) KX,KY
   80 FORMAT(2I6,' = KX KY  ONE OR BOTH OF KX KY TOO SMALL SO STOP'/)
      GO TO 150
   90 CONTINUE
      READ(8,10) EXPLAN
C     BELOW SC = COS(ANGLE) + I SIN(ANGLE)
C     READ IN DESIRED COSINE, COMPUTE ANGLE, THEN SINE
      READ(8,*) SCR
      ANGLE = DARCOS(SCR)
      SCI = DSIN(ANGLE)
      SC = DCMPLX(SCR,SCI)
      WRITE(6,100) SC
C     IF (IVEC.NE.0.AND.MTYPE.EQ.0) WRITE(9,7) SC
  100 FORMAT(' GENERATOR OF DIAGONAL TRANSFORMATION ='/2E20.12)
C
      TC = SC
      DO 110 J=2,KX
  110 TC = SC*TC
      WRITE(6,120) TC
  120 FORMAT(' TC = ',2E20.12)
C
      N = KX*KY
      C2 = ONE
      C1 = HALF-C0
      TEMP = DSQRT(2.0D0)
      IF (MTYPE.EQ.0) TEMP = ONE
      CL0 = -SC*C0
      CL1 = -TC*C1
      CL3 = -SC*C0*TEMP
      CL4 = -TC*C1*TEMP
C
      WRITE(6,130) N,MTYPE,KX,KY,C2,C0,C1
  130 FORMAT(/5X,'N',1X,'MTYPE',4X,'KX',4X,'KY',7X,'DIAGONAL',
     1 3X,'X-CODIAGONAL',3X,'Y-CODIAGONAL'/4I6,3E15.8/)
C
C-----------------------------------------------------------------------
       CALL HMATVE(C2,CL0,CL1,CL3,CL4,KX,KY)
C-----------------------------------------------------------------------
C
      IF(IVEC.EQ.0.OR.MTYPE.NE.0) GO TO 140
C
C     COMPUTE THE EXACT EIGENVALUES
C
C-----------------------------------------------------------------------
      CALL EXEVG(EIGVAL,C0,C1,C2,GAPS,MULTS,KX,KY)
C-----------------------------------------------------------------------
C
      IF(IVEC.LT.0) GO TO 150
C
  140 CONTINUE
      RETURN
C
C-----END OF USPEC------------------------------------------------------
  150 STOP
      END
C
C-----START OF CMATV FOR HERMITIAN POISSON MATRICES---------------------
C
C     SUBROUTINE HMATV(W,U,SUM)
      SUBROUTINE CMATV(W,U,SUM)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  C2,SUM
      COMPLEX*16 U(1),W(1)
      COMPLEX*16 CL0,CL1,CL3,CL4,CR0,CR1,CR3,CR4
C-----------------------------------------------------------------------
C     CALCULATES U = A*W - SUM*U
C
      N = KK*LL
      CR0 = DCONJG(CL0)
      CR1 = DCONJG(CL1)
      CR3 = DCONJG(CL3)
      CR4 = DCONJG(CL4)
C
C-----------------------------------------------------------------------
C    FIRST AND LAST BLOCKS
      J = 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CR1*W(J+KK)) - SUM*U(J)
      J = 2
      U(J)=(C2*W(J)+CL3*W(J-1)+CR0*W(J+1)+CR1*W(J+KK))-SUM*U(J)
      J = KK
      U(J)=(C2*W(J)+CL3*W(J-1)+CR1*W(J+KK))-SUM*U(J)
      J = KK - 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL0*W(J-1)+CR1*W(J+KK))-SUM*U(J)
      J = N - KK + 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL4*W(J-KK))-SUM*U(J)
      J = N - KK + 2
      U(J)=(C2*W(J)+CL3*W(J-1)+CR0*W(J+1)+CL4*W(J-KK))-SUM*U(J)
      J = N
      U(J)=(C2*W(J)+CL3*W(J-1)+CL4*W(J-KK))-SUM*U(J)
      J = N - 1
      U(J)=(C2*W(J)+CL0*W(J-1)+CR3*W(J+1)+CL4*W(J-KK))-SUM*U(J)
C
      KK2 = KK - 2
      DO 10 JJ = 3,KK2
      J = JJ
      U(J)=(C2*W(J)+CL0*W(J-1)+CR0*W(J+1)+CR1*W(J+KK))-SUM*U(J)
      J = N - KK + JJ
   10 U(J)=(C2*W(J)+CL0*W(J-1)+CR0*W(J+1)+CL4*W(J-KK))-SUM*U(J)
C
C    START BLOCKS 2 AND LL-1
      J = KK + 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL1*W(J-KK)+CR1*W(J+KK))-SUM*U(J)
      J = KK + 2
      U(J)=(C2*W(J)+CL3*W(J-1)+CR0*W(J+1)+CL1*W(J-KK)+CR1*W(J+KK))
     1 -SUM*U(J)
      J = KK + KK
      U(J)=(C2*W(J)+CL3*W(J-1)+CL1*W(J-KK)+CR1*W(J+KK))-SUM*U(J)
      J = KK + KK - 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL0*W(J-1)+CL1*W(J-KK)+CR1*W(J+KK))
     1  -SUM*U(J)
      J = N - 2*KK + 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CR4*W(J+KK)+CL1*W(J-KK))
     1 -SUM*U(J)
      J = N - 2*KK + 2
      U(J)=(C2*W(J)+CL3*W(J-1)+CR0*W(J+1)+CR4*W(J+KK)+CL1*W(J-KK))
     1 -SUM*U(J)
      J = N - KK
      U(J)=(C2*W(J)+CL3*W(J-1)+CR4*W(J+KK)+CL1*W(J-KK))-SUM*U(J)
      J = N - KK - 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL0*W(J-1)+CR4*W(J+KK)+CL1*W(J-KK))
     1 -SUM*U(J)
C
      DO 20 JJ = 3,KK2
      J = KK + JJ
      U(J)=(C2*W(J)+CL0*W(J-1)+CR0*W(J+1)+CL1*W(J-KK)+CR1*W(J+KK))
     1 -SUM*U(J)
      J = N - 2*KK + JJ
      U(J)=(C2*W(J)+CL0*W(J-1)+CR0*W(J+1)+CR4*W(J+KK)+CL1*W(J-KK))
     1 -SUM*U(J)
   20 CONTINUE
C
C    MIDDLE BLOCKS
      LL2 = LL - 2
      JP = KK
      DO 40 JJ = 3,LL2
      JP = JP + KK
C     JP = (JJ-1)*KK
      J = JP + 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL1*W(J-KK)+CR1*W(J+KK))-SUM*U(J)
      J = J + 1
      U(J)=(C2*W(J)+CL3*W(J-1)+CR0*W(J+1)+CL1*W(J-KK)+
     1 CR1*W(J+KK))-SUM*U(J)
      J = J + KK - 2
      U(J) = (C2*W(J)+CL3*W(J-1)+CL1*W(J-KK)+CR1*W(J+KK))-SUM*U(J)
      J = J - 1
      U(J)=(C2*W(J)+CR3*W(J+1)+CL0*W(J-1)+CL1*W(J-KK)+
     1 CR1*W(J+KK))-SUM*U(J)
C
      DO 30 II = 3,KK2
      J = JP + II
      U(J)=(C2*W(J)+CL0*W(J-1)+CR0*W(J+1)+CL1*W(J-KK)+CR1*W(J+KK))
     1 -SUM*U(J)
   30 CONTINUE
C
   40 CONTINUE
C
      RETURN
C
      ENTRY HMATVE(C2,CL0,CL1,CL3,CL4,KK,LL)
C
C-----END OF HMATV------------------------------------------------------
      RETURN
      END
C
C-----START OF EXEVG----------------------------------------------------
C
C     FOR MTYPE = 0, ZERO BOUNDARY CONDITIONS:
C     COMPUTES EXACT EIGENVALUES OF HERMITIAN POISSON MATRIX,
C     THEIR MULTIPLICITIES, AND THE GAPS BETWEEN THE EIGENVALUES AND
C     PUTS THEM RESPECTIVELY INTO VECTORS U, MP, AND G.  THESE
C     QUANTITIES ARE ALL WRITTEN TO FILE 9.
C
      SUBROUTINE EXEVG(U,C0,C1,C2,G,MP,KX,KY)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  U(1),MACHEP
      DOUBLE PRECISION  EPSM,C0,C1,C2,T0,T1,PIK,PIL,ONE,TWO,ATOLN,EE
      REAL G(1)
      INTEGER MP(1)
C-----------------------------------------------------------------------
      DATA MACHEP/Z3410000000000000/
      EPSM = 2.0D0*MACHEP
C-----------------------------------------------------------------------
      N = KX*KY
      ONE  = 1.0D0
      TWO  = 2.0D0
      T0 = DARCOS(-ONE)
      T1 = DFLOAT(KX+1)
      PIK = T0/T1
      T1 = DFLOAT(KY+1)
      PIL = T0/T1
C     GENERATE EXACT EIGENVALUES
      KP = 0
      DO 20 J = 1,KY
      T1 = PIL*DFLOAT(J)
      T0 = C2 - TWO*C1*DCOS(T1)
      DO 10 I = 1,KX
      KP = KP+1
      T1 = PIK*DFLOAT(I)
   10 U(KP) = T0 - TWO*C0*DCOS(T1)
   20 CONTINUE
C
C     ORDER U VECTOR BY INCREASING ALGEBRAIC SIZE
      DO 40 K = 2,N
      KM1 = K-1
      DO 30 L = 1,KM1
      JJ = K-L
      IF (U(JJ+1).GE.U(JJ)) GO TO 40
      T0 = U(JJ)
      U(JJ) = U(JJ+1)
   30 U(JJ+1) = T0
   40 CONTINUE
      ATOLN = DMAX1(DABS(U(1)),DABS(U(N)))*EPSM
C
      WRITE(9,50)
   50 FORMAT(' TRUE EIGENVALUES FOR HERMITIAN POISSON')
C
      WRITE(9,60)N,KX,KY,C2,C0,C1,ATOLN
      WRITE(6,60) N,KX,KY,C2,C0,C1,ATOLN
   60 FORMAT(1X,'A-SIZE',2X,'X-DIM',2X,'Y-DIM'/3I7/
     1 5X,'A-DIAGONAL',3X,'X-CODIAGONAL',3X,'Y-CODIAGONAL',10X,'ATOLN'/
     2 4E15.8)
C
C     DETERMINE TRUE MULTIPLICITIES FOR EXACT EIGENVALUES
      I = 1
      IDEX = 1
      J = 1
      NEXACT = 0
   70 J = J+1
      IF (J.GT.N) GO TO 80
      EE = DABS(U(J)-U(I))
      IF (EE.GT.ATOLN) GO TO 80
      IDEX = IDEX+1
      GO TO 70
   80 NEXACT = NEXACT+1
      U(NEXACT) = U(I)
      MP(NEXACT) = IDEX
C     MP(K) = MULTIPLICITY OF KTH EIGENVALUE CLUSTER FOR A
      IDEX = 1
      I = J
      IF (I.GT.N) GO TO 90
      GO TO 70
   90 CONTINUE
C
C     MULTIPLICITIES HAVE BEEN DETERMINED
C     NEXACT = NUMBER OF DISTINCT A-EIGENVALUES
C
      WRITE(9,100)NEXACT
      WRITE(6,100)NEXACT
  100 FORMAT(I6,' = NUMBER OF TRUE A-EIGENVALUES WHICH ARE DISTINCT'/)
C
C     MINGAP CALCULATION FOR DISTINCT A-EIGENVALUES
      NM1 = NEXACT - 1
      G(NEXACT) = U(NM1)-U(NEXACT)
      G(1) = U(2)-U(1)
C
      DO 110 J = 2,NM1
      T0 = U(J)-U(J-1)
      T1 = U(J+1)-U(J)
      G(J) = T1
      IF (T0.LT.T1) G(J) = -T0
  110 CONTINUE
C
C     NEXACT DISTINCT A-EIGENVALUES ARE IN U IN ASCENDING ORDER
C     MP = MULTIPLICITIES OF THE DISTINCT EIGENVALUES OF A
C     G = TRUE MINIMUM GAP IN A FOR EACH OF THESE EIGENVALUES
C     G < 0 INDICATES THE LEFT-HAND GAP WAS MINIMAL.
C     OUTPUT MULTIPLICITIES, DISTINCT EVS, AND MINGAPS TO FILE 11
C
      WRITE(9,120)
  120 FORMAT(5X,'I',1X,'AMULT',5X,'TRUE A-EIGENVALUE(I)',
     1 3X,'A-MINGAP(I)')
C
      WRITE(9,130)(J,MP(J),U(J),G(J), J=1,NEXACT)
  130 FORMAT(2I6,E25.16,E14.3)
C
      WRITE(9,140)
  140 FORMAT(' NEXACT DISTINCT A-EIGENVALUES ARE IN ASCENDING ORDER'/
     1 ' AMULT = MULTIPLICITIES OF THE DISTINCT EIGENVALUES OF A.'/
     2 ' A-MINGAP(I) = TRUE MINIMUM GAP IN A FOR EACH EIGENVALUE.'/
     3 ' A-MINGAP(I) LT 0 INDICATES THE LEFT-HAND GAP WAS MINIMAL.'//)
C
C     WE ORDER U VECTOR BY INCREASING SIZE OF THE GAPS
C
      DO 150 K = 1,N
  150 MP(K) = K
C
      DO 170 K = 2,N
      KM1 = K-1
C
      DO 160 L = 1,KM1
      JJ = K - L
      IF (ABS(G(JJ+1)).GE.ABS(G(JJ))) GO TO 170
      EE = U(JJ)
      U(JJ) = U(JJ+1)
      U(JJ+1) = EE
      GG = G(JJ)
      G(JJ) = G(JJ+1)
      G(JJ+1) = GG
      IEE = MP(JJ)
      MP(JJ) = MP(JJ+1)
  160 MP(JJ+1) = IEE
C
  170 CONTINUE
C
      WRITE(9,180)
  180 FORMAT(5X,'K',6X,'A-MINGAP',5X,'TRUE A-EIGENVALUE(I)',2X,'A-EVNO')
C
      WRITE(9,190)(J,G(J),U(J),MP(J), J=1,NEXACT)
  190 FORMAT(I6,E14.3,E25.16,I8)
C
      WRITE(9,200)
  200 FORMAT(' NEXACT DISTINCT A-EIGENVALUES. GAPS IN ASCENDING ORDER'/
     2 ' A-MINGAP(I) = TRUE MINIMUM GAP IN A FOR EACH EIGENVALUE.'/
     3 ' A-MINGAP(I) LT 0 INDICATES THE LEFT-HAND GAP WAS MINIMAL.'/
     3 ' A-MATRIX IS BLOCK TRIDIAGONAL AND EACH DIAGONAL BLOCK IS OF ORD
     3ER NX.'/
     4 ' NX = NUMBER OF POINTS ON EACH X-LINE. THERE ARE NY DIAGONAL BLO
     4CKS.'/
     5 ' NY = NUMBER OF POINTS ON EACH Y-LINE.'/
     5 ' A-DIAGONAL   = A(K,K)'/
     6 ' X-CODIAGONAL = A(I,I+1)'/
     7 ' Y-CODIAGONAL = A(I,I+NX)'/
     8 ' ----- END OF FILE 9 EXACTEV----------------------------'//)
C
C-----END OF EXEVG------------------------------------------------------
C
      RETURN
      END
C
C-----START OF HEXVEC---------------------------------------------------
C
C     FOR THE HERMITIAN POISSON TEST CASES WITH MTYPE = 0 ONLY:
C     FOR A GIVEN RITZ VECTOR V AND EIGENVALUE X1, COMPUTES
C     THE CLOSEST TRUE EIGENVALUE Y1 AND CORRESPONDING TRUE
C     EIGENVECTOR Z, CALCULATES THE NORM OF V-Z AND THE MAXIMAL
C     DIFFERENCE OF THE COMPONENTS.  USER WOULD HAVE TO
C     INCORPORATE ENTRY AND CALL TO THIS SUBROUTINE INTO
C     HLEVEC PROGRAM IF THESE QUANTITIES ARE DESIRED.
C     U CONTAINS THE COMPUTED TRUE EIGENVALUES.
C     W CONTAINS THE TRUE EIGENVECTOR FOR THE REAL POISSON MATRIX
C
      SUBROUTINE HEXVEC(Z,V,U,W,X1,Y1,MP,JNUM)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  U(1),W(1)
      DOUBLE PRECISION WI(110),WJ(110),WII(110)
      DOUBLE PRECISION X1,Y1,EV,EE,WS,PIK,PIL,SUM,TEMP
      DOUBLE PRECISION ATOLN,EPSM,ZERO,HALF,ONE,TWO,MACHEP
      DOUBLE PRECISION C0,C1,C2,T0,T1,T2
      COMPLEX*16 CONE,S,SB,STEMP,V(1),Z(1)
      INTEGER MP(1)
C-----------------------------------------------------------------------
      DATA MACHEP/Z3410000000000000/
      EPSM = 2.0D0*MACHEP
C-----------------------------------------------------------------------
C     THIS PROGRAM CALCULATES THE EXACT EIGENVALUES AND EIGENVECTORS
C     OF THE HERMITIAN POISSON MATRIX A OF ORDER  N = KX BY KY
C     A CONSISTS OF KY TRIDIAGONAL BLOCKS OF ORDER KX
C     KX = X-DIMENSION     KY = Y-DIMENSION.
C
C     C2 = DIAGONAL OF KX BY KX MATRIX
C    -C0 = CO-DIAGONAL OF THE KX BY KX MATRIX.
C    -C1 = Y-CODIAGONAL.
C
C     NOTE THAT THE VECTORS WI,WJ,WII ARE DIMENSIONED INTERNALLY
C     THEY ARE USED JUST TO KEEP FROM REGENERATING INFORMATION.
C     WI,WII = REAL*8 ARRAYS OF DIMENSION AT LEAST KX
C     WJ     = REAL*8 ARRAY  OF DIMENSION AT LEAST KY.
C
C     NOTATION USED IN PROGRAM
C
C     PIK = ARCOS(-1)/(KX+1)    PIL = ARCOS(-1)/(KY+1)
C     WI(I) = PIK*I        WJ(J) = PIL*J
C
C     T0 = C2 - 2*C1*COS(PIL*J)    EV(I,J) = T0 - 2*C0*COS(PIK*I)
C     I = 1,KX    J = 1,KY    KP = (J-1)*KX + I
C
C     W(KV) = SIN(PIK*I*IK)*SIN(PIL*J*JK)
C     IK = 1,KX    JK = 1,KY    KV = (JK-1)*KX + IK
C     W IS UNSCALED EIGENVECTOR FOR EV(I,J)
C     WS = 1/||W||: ||W|| = .5*DSQRT(T2*T3)  T2 = KX+1  T3 = KY+1
C     U(K) IS A-EV ORDERED BY INCREASING SIZE, K = 1,N
C
C     GIVEN X1 FIND Y1 AND KVEC SUCH THAT
C           Y1 = EV(KVEC) AND |X1-Y1| = MIN
C           ALSO GIVEN UNIT RITZ VECTOR ASSOCIATED WITH X1
C           CALCULATE UNIT EIGENVECTOR W, A*W = Y1*W
C           T2 = ||V-W||  T1 = MAX(|V(K)-W(K)|, K= 1,N)
C           MAX OCCURS FIRST AT K = KK
C
C-----------------------------------------------------------------------
C     C2 = A(K,K)
C     C0 = A(K,K+1) = A(K+1,K)
C     C1 = A(K,K+KX) = A(K+KX,K)
C     C0 + C1 = HALF
C
C     SPECIFY PARAMETERS
      N = KX*KY
      ZERO = 0.0D0
      HALF = 0.5D0
      ONE  = 1.0D0
      TWO  = 2.0D0
      T0 = DARCOS(-ONE)
      T1 = DFLOAT(KX+1)
      PIK = T0/T1
      T2 = DFLOAT(KY+1)
      PIL = T0/T2
      WS = TWO/DSQRT(T1*T2)
C
C     GENERATE WI WJ VECTORS
      KP = 0
      DO 20 J = 1,KY
      T1 = PIL*DFLOAT(J)
      WJ(J) = T1
      T0 = C2 - TWO*C1*DCOS(T1)
      DO 10 I = 1,KX
      KP = KP+1
      T1 = PIK*DFLOAT(I)
      WI(I) = T1
   10 U(KP) = T0 - TWO*C0*DCOS(T1)
   20 CONTINUE
C     U(KP) = EV(I,J) = C2 - 2*C1*COS(PIL*J) - 2*C0*COS(PIK*I)
C
C     INITIALIZE MP VECTOR
      DO 30 K = 1,N
   30 MP(K) = K
C
C     WE ORDER U VECTOR BY INCREASING SIZE OF THE EVS
      DO 50 K = 2,N
      KM1 = K-1
C
      DO 40 L = 1,KM1
      JJ = K - L
      IF (U(JJ+1).GE.U(JJ)) GO TO 50
      EE = U(JJ)
      U(JJ) = U(JJ+1)
      U(JJ+1) = EE
      IEE = MP(JJ)
      MP(JJ) = MP(JJ+1)
   40 MP(JJ+1) = IEE
C
   50 CONTINUE
C
      ATOLN = DMAX1(DABS(U(1)),DABS(U(N)))*EPSM
C
      WRITE(6,60) N,KX,KY,C2,C0,C1,ATOLN
   60 FORMAT(/' EXACT ERRORS FOR CONVERGED GOODEV'/
     1 4I6,' = N KX KY'//
     1 4E12.5,' = C2 C0 C1 ATOLN'//)
C
C     KP = MP(K) MEANS EIGENVALUE U(K) CORRESPONDS TO EIGENVECTOR W(KP)
C     COMPUTE TOLERANCE USED IN COMPUTING TRUE MULTIPLICITIES
C
C     X1 IS AN INPUT PARAMETER. WE CALCULATE EXACT
C     A-EIGENVALUE WHICH IS CLOSEST TO X1, LABEL IT Y1 AND CALCULATE
C     UNIT EIGENVECTOR OF A ASSOCIATED WITH Y1. A*W = Y1*W, ||W|| = 1.
C     Y1 = U(KEV). EIGENVALUES OF A ARE ORDERED BY INCREASING SIZE.
C     V = COMPLEX RITZ VECTOR ASSOCIATED WITH GOODEV X1
C     WE SHOULD HAVE V = D*W WHERE D = DIAG(D(1),D(2),..,D(N))
C     D(1) = ONE, D(K+1)/D(K) = SB, |SB| = ONE
C
      KX1 = 0
      IF (X1.LE.U(1)) KX1 = 1
      IF (X1.GE.U(N)) KX1 = N
      NM1 = N-1
      IF (KX1.NE.0) GO TO 80
C
      DO 70 KVEC = 2,N
      IF (X1.GE.U(KVEC)) GO TO 70
C     U(KVEC-1).LE.X1.LT.U(KVEC)
      T1 = X1 - U(KVEC-1)
      T2 = U(KVEC) - X1
      KX1 = KVEC - 1
      IF (T1.GT.T2) KX1 = KVEC
      GO TO 80
   70 CONTINUE
C
   80 Y1 = U(KX1)
C
      IF (KX1.EQ.1) EE = U(2) - U(1)
      IF (KX1.EQ.N) EE = U(N) - U(NM1)
      IF (KX1.EQ.1.OR.KX1.EQ.N) GO TO 90
      EE = DMIN1(U(KX1+1)-U(KX1),U(KX1)-U(KX1-1))
   90 CONTINUE
C
      T0 = DABS(ONE - X1/Y1)
C
      WRITE(6,100) N,KX1,JNUM,Y1,X1,T0,EE
  100 FORMAT(3I8,' = N, A-EV NUMBER,GOODEV NO'//
     1 18X,'EXACTEV',19X,'GOODEV',4X,'RELERROR',4X,'A-MINGAP'/
     1 2E25.16,2E12.3/)
C
      IF (EE.GT.ATOLN) GO TO 120
C
      WRITE(6,110)
  110 FORMAT(' Y1 IS A MULTIPLE EIGENVALUE OF A SO WE EXIT'/)
C
      GO TO 200
C
C     Y1 IS TOEPLITZ EIGENVALUE CLOSEST TO X1.
C     CALCULATION OF EIGENVECTOR ASSOCIATED WITH EIGENVALUE Y1
C     A*W = Y1*W
C
C     DETERMINE I J FROM K: MP(K) = KP = (J-1)*KX+I
  120 CONTINUE
      K = KX1
      KP = MP(K)
      I = MOD(KP,KX)
      IF (I.EQ.0) I = KX
      T1 = WI(I)
      J = 1 + (KP-1)/KX
      T2 = WJ(J)
C
      DO 130 II = 1,KX
      T0 = T1*DFLOAT(II)
  130 WII(II) = WS*DSIN(T0)
C
      KV = 0
      DO 150 JJ = 1,KY
      T0 = T2*DFLOAT(JJ)
      T0 = DSIN(T0)
C
      DO 140 II = 1,KX
      KV = KV + 1
  140 W(KV) = T0*WII(II)
C
  150 CONTINUE
C
C     W IS UNIT EXACT EIGENVECTOR OF A ASSOCIATED WITH Y1
C     V IS UNIT COMPLEX RITZVECTOR OF B ASSOCIATED WITH X1
C
      CONE = DCMPLX(ONE,ZERO)
      STEMP = CONE
      DO 160 K = 1,N
      Z(K) = STEMP*W(K)
  160 STEMP = STEMP*SB
C
      T1 = ZERO
      S = ONE
      KK = 0
      DO 170 K = 1,N
      IF (CDABS(Z(K)).LE.T1) GO TO 170
      T1 = CDABS(Z(K))
      KK = K
  170 CONTINUE
C
      S = V(KK)/Z(KK)
C
      KK = 0
      T1 = ZERO
      T2 = ZERO
      DO 180 K = 1,N
      TEMP = CDABS(S*Z(K) - V(K))
      T2 = T2 + TEMP**2
      IF (TEMP.LE.T1) GO TO 180
      KK = K
      T1 = TEMP
  180 CONTINUE
C
      T2 = DSQRT(T2)
      WRITE(6,190) KK,T1,T2
  190 FORMAT(' EIGENVECTOR ERROR. MAX ERROR AT COMPONENT = ',I6/
     1 ' MAX CDABS(EXACTVEC(K)-RITZVEC(K)) = ',E12.5/
     1 ' NORM(EXACTVEC-RITZVEC) = ',E12.5/)
C
  200 CONTINUE
C
      RETURN
C
C-----------------------------------------------------------------------
      ENTRY EXVECP(SB,C0,C1,C2,KX,KY)
C-----------------------------------------------------------------------
C
C-----END OF HEXVEC-----------------------------------------------------
      RETURN
      END
C
C-----USPEC (TRIDIAGONAL HERMITIAN MATRICES)----------------------------
C
C     SUBROUTINE USPEC(N,MATNO)
      SUBROUTINE TSPEC(N,MATNO)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  D(100), DAR(100),DAI(100), PI, EIGVAL(100)
      DOUBLE PRECISION SPACE
      COMPLEX*16   DA(100),DB(100)
      REAL EXPLAN(20)
C-----------------------------------------------------------------------
C    DIMENSION ARRAYS NEEDED TO DEFINE MATRIX.  THEN
C    PASS THE STORAGE LOCATIONS OF THESE ARRAYS TO THE MATRIX-VECTOR
C    MULTIPLY SUBROUTINE CMATV.
C
C    DIAGONAL ENTRY = D, ABOVE DIAGONAL ENTRY = DA, BELOW DIAGONAL = DB.
C
      READ(8,10) EXPLAN
   10 FORMAT(20A4)
      READ(8,*) NOLD,MATOLD
C
      WRITE(6,20) N,MATOLD
   20 FORMAT(I10,2I6,I8,' = N,MATOLD'/)
C
C     TEST OF PARAMETER CORRECTNESS
      ITEMP = (NOLD-N)**2 + (MATNO-MATOLD)**2
C
      IF(ITEMP.EQ.0) GO TO 40
C
      WRITE(6,30)
   30 FORMAT(' PROGRAM TERMINATES BECAUSE EITHER ORDERS OF OR LABELS FOR
     1 MATRIX DISAGREE')
      GO TO 250
C
   40 CONTINUE
C
C     IF ITOEP = 1 THEN MATRIX IS TOEPLITZ AND WE PRINT OUT TRUE
C     EIGENVALUES
      READ(8,10) EXPLAN
      READ(8,*) ITOEP
      READ(8,10) EXPLAN
C
      IF(ITOEP.EQ.1) WRITE(6,50)
   50 FORMAT(/' TEST MATRIX IS HERMITIAN TOEPLITZ'/)
      IF(ITOEP.NE.1) GO TO 110
C
      READ(8,*) DAR(1),DAI(1),D(1)
      DA(1) = DCMPLX(DAR(1),DAI(1))
      DB(1) = DCONJG(DA(1))
      DO 60 J=2,N
      D(J) = D(1)
      DA(J) = DA(1)
   60 DB(J) = DB(1)
      WRITE(6,70) DB(1),D(1),DA(1)
      WRITE(9,70) DB(1),D(1),DA(1)
   70 FORMAT(' HERMITIAN TOEPLITZ MATRIX IS USED.'/' BELOW DIAGONAL ENTR
     1Y = ',2E12.3/' DIAGONAL ENTRY = ',E12.3/' ABOVE DIAGONAL ENTRY ='
     1,2E12.3)
C
C     COMPUTE THE TRUE EIGENVALUES.  FORMULA IS CORRECT ONLY FOR THOSE
C     MATRICES WHOSE DIAGONAL = 2., ABOVE DIAGONAL = A, BELOW DIAGONAL
C     = A-CONJUGATE, AND A HAS NORM 1.
C
      PI = DARCOS(-1.D0)
      DO 80 J=1,N
   80 EIGVAL(J) = 2.D0 * (1.D0 -DCOS(PI*DFLOAT(J)/DFLOAT(N+1)))
      WRITE(9,90) N
   90 FORMAT(I6, ' = ORDER OF MATRIX'/' TRUE EIGENVALUES ARE'/)
      WRITE(9,100) (J, EIGVAL(J), J=1,N)
  100 FORMAT(I5,4X,E25.16,6X,I5,4X,E25.16)
      GO TO 240
C
C     NONTOEPLITZ HERMITIAN.  DIAGONAL ENTRIES ARE EQUALLY-SPACED.
C     ABOVE DIAGONAL ENTRIES ARE GENERATED BY GENERATING EQUALLY-SPACED
C     REAL PARTS, AND EQUALLY-SPACED IMAGINARY PARTS.  THE BELOW
C     DIAGONAL ENTRIES ARE THEN OBTAINED BY TAKING THE COMPLEX CONJUGATE
C     OF THE ABOVE DIAGONAL ENTRIES
C
  110 READ(8,*) D(1), SPACE
      WRITE(6,120) D(1),SPACE
  120 FORMAT(' 1ST DIAGONAL ENTRY =',E20.12,' SPACING =',E20.12)
      DO 130 J=2,N
  130 D(J) = D(J-1) + SPACE
      WRITE(6,140) (D(J), J=1,3)
  140 FORMAT(' 1ST THREE DIAGONAL ENTRIES ='/(2E20.12))
      READ(8,10) EXPLAN
      READ(8,*) DAR(1), SPACE
      WRITE(6,150) DAR(1),SPACE
  150 FORMAT(' REAL PART OF 1ST ABOVE DIAGONAL ENTRY =',E20.12,/
     1' SPACING = ',E20.12)
      DO 160 J=2,N
  160 DAR(J) = DAR(J-1) + SPACE
      WRITE(6,170) (DAR(J), J=1,3)
  170 FORMAT(' REAL PARTS OF 1ST THREE ABOVE DIAGONAL ENTRIES ='/
     1(2E20.12))
      READ(8,10) EXPLAN
      READ(8,*) DAI(1), SPACE
      WRITE(6,180) DAI(1),SPACE
  180 FORMAT(' IMAGINARY PART OF 1ST ABOVE =',E20.12,/' SPACING =',
     1 E20.12)
      DO 190 J=2,N
  190 DAI(J) = DAI(J-1) + SPACE
      WRITE(6,200) (DAI(J), J = 1,3)
  200 FORMAT(' IMAGINARY PARTS OF 1ST THREE ABOVE DIAGONAL ENTRIES ='/
     1 (2E20.12))
      DO 210 J=1,N
      DA(J) = DCMPLX(DAR(J),DAI(J))
  210 DB(J) = DCONJG(DA(J))
C
      WRITE(9,220) (D(J), J=1,N)
  220 FORMAT(' DIAGONAL ENTRIES ='/(4E20.12))
      WRITE(9,230) (DA(J), J=1,N)
  230 FORMAT(' ABOVE DIAGONAL ENTRIES'/(4E20.12))
C
C     PASS STORAGE LOCATIONS OF ARRAYS THAT DEFINE THE MATRIX TO
C     THE MATRIX-VECTOR MULTIPLY SUBROUTINE CMATV
C
  240 CONTINUE
C
C-----------------------------------------------------------------------
      CALL TMATVE(DA,DB,D,N)
C-----------------------------------------------------------------------
C
      RETURN
  250 STOP
C
C-----END OF USPEC------------------------------------------------------
      END
C
C-----START OF MATRIX-VECTOR MULTIPLY (HERMITIAN TRIDIAGONAL)-----------
C
C     SUBROUTINE CMATV(W,U,SUM)
      SUBROUTINE TMATV(W,U,SUM)
C
C-----------------------------------------------------------------------
      COMPLEX*16   U(1),W(1),DA(1),DB(1)
      DOUBLE PRECISION  D(1),SUM
C-----------------------------------------------------------------------
C    HERMITIAN MATRIX-VECTOR MULTIPLY FOR LANCZS  U = A*W - SUM*U
C    MATRIX IS TRIDIAGONAL HERMITIAN TOEPLITZ
C-----------------------------------------------------------------------
C
C     COMPUTE A*W - SUM*U
      U(1) = D(1)*W(1) + DA(1)*W(2) - SUM*U(1)
      N1 = N-1
      DO 10 I = 2,N1
   10 U(I) = DB(I-1)*W(I-1)+D(I)*W(I) + DA(I)*W(I+1) -SUM*U(I)
      U(N) = DB(N-1)*W(N-1) + D(N)*W(N) - SUM*U(N)
C
      RETURN
C
C    STORAGE LOCATIONS ARE PASSED TO CMATV FROM USPEC
C
C-----------------------------------------------------------------------
      ENTRY TMATVE(DA,DB,D,N)
C-----------------------------------------------------------------------
C
C-----END OF CMATV------------------------------------------------------
      RETURN
      END


.