C-----LEMULT------------------------------------------------------------
C
C     CONTAINS SUBROUTINES  LANCZS, USPEC, AND CMATV
C     TO BE USED WITH THE REAL SYMMETRIC VERSION OF THE LANCZOS
C     EIGENVALUE/EIGENVECTOR PROCEDURES.
C     ALSO CONTAINS SUBROUTINES FOR POISSON TEST MATRICES THAT ALLOW
C     COMPUTATION OF TRUE ERRORS IN COMPUTED EIGENVALUES AND
C     IN CORRESPONDING EIGENVECTORS.
C
C     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
C         CMATV.
C     2.  IN THE SAMPLE USPEC AND CMATV FOR DIAGONAL TEST MATRICES:
C         FREE FORMAT (8,*) AND THE FORMAT (20A4).
C     3.  IN THE POISSON SUBROUTINES PROVIDED, THE DATA MACHEP
C         DEFINITION AND MANY OF THE INDICES FOR ARRAYS ARE NOT
C         IN A PORTABLE CONSTRUCTION.  THESE PROGRAMS SHOULD BE
C         REMOVED FROM THE LEMULT FILE IF THE USER IS NOT USING THEM.
C
C-----LANCZS-COMPUTE THE LANCZOS TRIDIAGONAL MATRICES-------------------
C
      SUBROUTINE LANCZS(MATVEC,ALPHA,BETA,V1,V2,G,KMAX,MOLD1,N,IIX)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  ALPHA(1),BETA(1),V1(1),V2(1),SUM,TEMP,ONE,ZERO
      REAL  G(1)
      DOUBLE PRECISION FINPRO, DSQRT
      EXTERNAL MATVEC
C-----------------------------------------------------------------------
C
      ZERO = 0.D0
      ONE = 1.D0
C
      IF(MOLD1.GT.1)GO TO 30
C
C    ALPHA/BETA GENERATION STARTS AT I = 1
C    MOLD1 = 1 SET V1 = 0. AND V2 = RANDOM UNIT VECTOR
      BETA(1) = ZERO
      IIL=IIX
C
C-----------------------------------------------------------------------
      CALL GENRAN(IIL,G,N)
C-----------------------------------------------------------------------
C
      DO 10 I = 1,N
   10 V2(I) = G(I)
C
C-----------------------------------------------------------------------
      SUM = FINPRO(N,V2(1),1,V2(1),1)
C-----------------------------------------------------------------------
C
      SUM = ONE/DSQRT(SUM)
      DO 20 I = 1,N
      V1(I) = ZERO
   20 V2(I) = V2(I)*SUM
C
C    ALPHA BETA GENERATION LOOP
   30 CONTINUE
C
      DO 60 I=MOLD1,KMAX
      SUM = BETA(I)
C    MATVEC(V2,V1,SUM) CALCULATES  V1 = A*V2 - SUM*V1
C
C-----------------------------------------------------------------------
      CALL MATVEC(V2,V1,SUM)
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      SUM = FINPRO(N,V1(1),1,V2(1),1)
C-----------------------------------------------------------------------
C
      ALPHA(I) = SUM
      DO 40 J=1,N
   40 V1(J) = V1(J)-SUM*V2(J)
C
C-----------------------------------------------------------------------
      SUM = FINPRO(N,V1(1),1,V1(1),1)
C-----------------------------------------------------------------------
C
      IN = I+1
      BETA(IN) = DSQRT(SUM)
      SUM = ONE/BETA(IN)
      DO 50 J=1,N
      TEMP = SUM*V1(J)
      V1(J) = V2(J)
   50 V2(J) = TEMP
   60 CONTINUE
C
C     END ALPHA, BETA GENERATION LOOP
C
      RETURN
C-----END OF LANCZS-----------------------------------------------------
      END
C
C-----USPEC (GENERAL SYMMETRIC SPARSE MATRICES)-------------------------
C
C     SUBROUTINE USPEC(N,MATNO)
      SUBROUTINE GUSPEC(N,MATNO)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  A(10000),AD(5010)
      INTEGER  IROW(10000),ICOL(5010)
C-----------------------------------------------------------------------
C    USPEC DIMENSIONS AND INITIALIZES THE ARRAYS NEEDED TO DEFINE
C    THE USER-SPECIFIED MATRIX AND THEN PASSES THE STORAGE LOCATIONS
C    OF THESE ARRAYS TO THE MULTIPLY SUBROUTINE CMATV.
C
C    MATRIX IS STORED IN FOLLOWING SPARSE MATRIX FORMAT:
C    N = ORDER OF A-MATRIX,
C    NZS = NUMBER OF NONZERO SUBDIAGONAL ENTRIES,
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 CORRESPONDING ROW INDEX FOR A(K).
C    AD(I), I=1,N CONTAINS DIAGONAL ENTRIES (INCLUDING ANY 0
C           DIAGONAL ENTRIES).
C    A(K), K=1,NZS CONTAINS NONZERO SUBDIAGONAL ENTRIES, BY COLUMN
C    FOR J > NZL THERE ARE NO NONZERO SUBDIAGONAL ELEMENTS IN COLUMN J.
C    ICOL(J) = 0 IS ALLOWED
C
C-----------------------------------------------------------------------
C     ARRAYS THAT DEFINE THE MATRIX 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 FOR
     1 MATRIX DISAGREE')
      GO TO 70
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
C     DIAGONAL IS READ FIRST, THEN NONZERO BELOW DIAGONAL ENTRIES
      READ(8,60) (AD(K), K=1,N)
      READ(8,60) (A(K), K=1,NZS)
   60 FORMAT(4E19.10)
C
C-----------------------------------------------------------------------
C     PASS STORAGE LOCATIONS OF ARRAYS THAT DEFINE THE MATRIX TO
C     THE MATRIX-VECTOR MULTIPLY SUBROUTINE CMATV
C
      CALL CMATVE(A,AD,ICOL,IROW,N,NZL)
C-----------------------------------------------------------------------
C
      RETURN
   70 STOP
C-----END OF USPEC------------------------------------------------------
      END
C
C-----MATRIX-VECTOR MULTIPLY FOR REAL SPARSE SYMMETRIC MATRICES---------
C
C     SUBROUTINE CMATV(W,U,SUM)
      SUBROUTINE GCMATV(W,U,SUM)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION   U(1),W(1),A(1),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
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) + 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
C
      ENTRY CMATVE(A,AD,ICOL,IROW,N,NZL)
C
C-----------------------------------------------------------------------
C
      RETURN
C-----END OF CMATV------------------------------------------------------
      END
C
C-----MATRIX-VECTOR MULTIPLY FOR DIAGONAL TEST MATRICES-----------------
C
C     SUBROUTINE CMATV(W,U,SUM)
      SUBROUTINE DCMATV(W,U,SUM)
C
C     CMATV COMPUTES  U = (DIAGONAL MATRIX) * W - SUM * U
C-----------------------------------------------------------------------
      DOUBLE PRECISION  W(1),U(1),SUM
      DOUBLE PRECISION  D(1)
C-----------------------------------------------------------------------
C
      DO 10 I=1,N
   10 U(I)= D(I)*W(I) - SUM*U(I)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY MVDIAE(D,N)
C
C-----------------------------------------------------------------------
C
      RETURN
C-----END OF DIAGONAL TEST MATRIX  MULTIPLY-----------------------------
      END
C
C
C-----START OF USPEC FOR DIAGONAL TEST MATRIX---------------------------
C
C     SUBROUTINE USPEC(N,MATNO)
      SUBROUTINE DUSPEC(N,MATNO)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  D(1000), SHIFT, SPACE
      DOUBLE PRECISION  DABS, DFLOAT
      REAL  EXPLAN(20)
C-----------------------------------------------------------------------
C
      READ(8,10) EXPLAN
   10 FORMAT(20A4)
      READ(8,*) NOLD,NUNIF,SPACE,D(1),SHIFT
      NNUNIF = NOLD - NUNIF
      WRITE(6,20) NOLD,SPACE,NNUNIF,D(1),SHIFT
   20 FORMAT(/' DIAGONAL TEST MATRIX, SIZE = ',I4/' MOST ENTRIES ARE ',
     1E10.3,' UNITS APART.',I3,' ENTRIES'/' ARE IRREGULARLY SPACED. FIRS
     1T ENTRY IS ',E10.3,' SHIFT = ',E10.3/)
C
      IF(N.NE.NOLD) GO TO 90
C     COMPUTE THE UNIFORM PORTION OF THE SPECTRUM
      DO 30 J=2,NUNIF
   30 D(J) = D(1) - DFLOAT(J-1)*SPACE
      NUNIF1=NUNIF + 1
      READ(8,10)  EXPLAN
      DO 40 J=NUNIF1,N
   40 READ(8,*) D(J)
      NB = NUNIF - 2
C
      IF(SHIFT.EQ.0.) GO TO 60
      DO 50 J=1,N
   50 D(J) = D(J) + SHIFT
C
C     PRINT OUT THE EIGENVALUES OF INTEREST
   60 WRITE(6,70) (D(I), I=1,10 )
      WRITE(6,80) (D(I), I = NB,N)
   70 FORMAT(/' REAL SYMMETRIC LANCZOS TEST, 1ST 10 ENTRIES OF DIAGONAL
     1TEST MATRIX = '/(3E22.14))
   80 FORMAT(/' MIDDLE UNIFORM PORTION OF MATRIX IS NOT PRINTED OUT'/
     1' END OF UNIFORM PLUS NONUNIFORM SECTION = '/(3E25.16))
C
C     DIAGONAL GENERATION COMPLETE
C
C-----------------------------------------------------------------------
C    CALL ENTRY TO MATRIX-VECTOR MULTIPLY SUBROUTINE  TO PASS
C    STORAGE LOCATION OF D-ARRAY AND ORDER OF A-MATRIX.
       CALL MVDIAE(D,N)
C-----------------------------------------------------------------------
C
      RETURN
   90 WRITE(6,100) NOLD,N
  100 FORMAT(' PROGRAM TERMINATES BECAUSE NOLD = ',I5,'DOES NOT EQUAL N
     1 =',I5)
C-----END OF USPEC SUBROUTINE FOR DIAGONAL TEST MATRICES----------------
      STOP
      END
C
C-----POISSON TEST MATRICES---------------------------------------------
C
C     CONTAINS SUBROUTINES USPEC, CMATV, EXEVG, EXERR AND  EXVEC
C
C-----START OF USPEC----------------------------------------------------
C
      SUBROUTINE USPEC(N,MATNO)
C     SUBROUTINE PUSPEC(N,MATNO)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  C0,C1,C2,HALF,ONE
      REAL  EXPLAN(20)
C--------------------------------------------------------------------
      HALF = 0.5D0
      ONE  = 1.0D0
C
C     READ USER-SPECIFIED PARAMETERS FROM INPUT FILE 8 (FREE FORMAT)
C
      READ(8,10) EXPLAN
      WRITE(6,10) EXPLAN
   10 FORMAT(20A4)
C
      READ(8,10) EXPLAN
      READ(8,*) KX,KY,C0
      N = KX*KY
      C1 = HALF-C0
      C2 = ONE
C
      WRITE(6,20) N,KX,KY,C2,C0,C1
   20 FORMAT(/5X,'N',4X,'KX',4X,'KY',7X,'DIAGONAL',3X,'X-CODIAGONAL',
     1 3X,'Y-CODIAGONAL'/3I6,3E15.8/)
C
C-----------------------------------------------------------------------
      CALL PMATVE(C0,C1,C2,KX,KY)
      CALL EXEVE(C0,C1,C2,KX,KY)
      CALL EXERRP(C0,C1,C2,KX,KY)
C     CALL EXVECP(C0,C1,C2,KX,KY)
C-----------------------------------------------------------------------
C
      RETURN
C-----END OF USPEC------------------------------------------------------
      END
C
C-----START OF CMATV----------------------------------------------------
C     CALCULATE U = A*W - SUM*U FOR REAL POISSON MATRICES
C
      SUBROUTINE CMATV(W,U,SUM)
C     SUBROUTINE PMATV(W,U,SUM)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  U(1),W(1)
      DOUBLE PRECISION  C0,C1,C2,CC0,CC1,SUM
C-----------------------------------------------------------------------
C
      N = KX*KY
      KX1 = KX-1
      KY1 = KY-1
C
      KK = 1
      U(KK) = (C2*W(KK)+C0*W(KK+1)+C1*W(KK+KX)) - SUM*U(KK)
      KK = KX
      U(KK) = (C2*W(KK)+C0*W(KK-1)+C1*W(KK+KX)) - SUM*U(KK)
      KK = N - KX + 1
      U(KK) = (C2*W(KK)+C0*W(KK+1)+C1*W(KK-KX)) - SUM*U(KK)
      KK = N
      U(KK) = (C2*W(KK)+C0*W(KK-1)+C1*W(KK-KX)) - SUM*U(KK)
C
      DO 10 J = 2,KX1
      KK = J
      U(KK) = (C2*W(KK)+C0*W(KK-1)+C0*W(KK+1)+C1*W(KK+KX)) - SUM*U(KK)
      KK = J+N-KX
      U(KK) = (C2*W(KK)+C0*W(KK-1)+C0*W(KK+1)+C1*W(KK-KX))-SUM*U(KK)
   10 CONTINUE
C
      DO 30 J = 2,KY1
      KK = (J-1)*KX + 1
      U(KK) = (C2*W(KK)+C0*W(KK+1)+C1*W(KK-KX)+C1*W(KK+KX)) - SUM*U(KK)
      DO 20 I = 2,KX1
      KK = KK + 1
      U(KK) = (C2*W(KK)+C0*W(KK-1)+C0*W(KK+1)+C1*W(KK-KX)
     1 +C1*W(KK+KX)) - SUM*U(KK)
   20 CONTINUE
      KK = KK + 1
      U(KK) = (C2*W(KK)+C0*W(KK-1)+C1*W(KK-KX)+C1*W(KK+KX)) - SUM*U(KK)
   30 CONTINUE
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY  PMATVE(CC0,CC1,C2,KX,KY)
C
C-----------------------------------------------------------------------
C
      C0 = -CC0
      C1 = -CC1
C
      RETURN
C-----END OF CMATV------------------------------- ----------------------
      END
C
C-----START OF EXEVG----------------------------------------------------
C
C     COMPUTES TRUE EIGENVALUES OF POISSON MATRIX, GAPS BETWEEN
C     TRUE EIGENVALUES, AND MULTIPLICITIES OF TRUE EIGENVALUES
C     AND STORE THESE VALUES, RESPECTIVELY, IN U, G, AND MP.
C     THESE QUANTITIES ARE WRITTEN OUT TO FILE 9
C
      SUBROUTINE EXEVG(U,G,MP)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  U(1)
      DOUBLE PRECISION MACHEP,EPSM,C0,C1,C2,T0,T1,PIK,PIL,ONE,TWO
      DOUBLE PRECISION 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 TRUE 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 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 MULTIPLICITIES FOR TRUE 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 9
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,NEXACT
  150 MP(K) = K
C
      DO 170 K = 2,NEXACT
      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 TRUEEV-----------------------------'//)
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY EXEVE(C0,C1,C2,KX,KY)
C
C-----------------------------------------------------------------------
C
      RETURN
C-----END OF EXEVG------------------------------------------------------
      END
C
C-----START OF EXERR----------------------------------------------------
C
C     FOR GIVEN COMPUTED EIGENVALUES, V(I), I=1,2,...,NG
C     COMPUTES THE CLOSEST TRUE EIGENVALUES AND THE ERROR IN THE
C     COMPUTED EIGENVALUES, AND STORES THESE RESPECTIVELY
C     IN U(I) AND IN G(MEV+I).  THESE QUANTITIES ARE WRITTEN
C     TO FILE 10.
C
      SUBROUTINE EXERR(V,U,G,MP,MEV,NG,NEXACT,IWRITE)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION U(1),V(1)
      DOUBLE PRECISION EV,EE,T0,T1,C0,C1,C2,PIK,PIL
      DOUBLE PRECISION ATOLN,EPSM,MACHEP,ZERO,ONE,TWO
      REAL G(1)
      INTEGER MP(1)
C-----------------------------------------------------------------------
      DATA MACHEP/Z3410000000000000/
      EPSM = 2.0D0*MACHEP
C-----------------------------------------------------------------------
C     ON ENTRY V CONTAINS NG GOOD EIGENVALUES OF T(1,MEV)
C     MP CONTAINS THE MULTIPLICITIES OF THESE EIGENVALUES.
C     U(I) = GAP TO NEAREST DISTINCT TMEV I=1,NG
C     U < 0. MEANS GAP IS DUE TO SPURIOUS EV
C
C     ON EXIT G(MEV+I) = ERROR FOR V(I) < 0. IF MULT EV > 1
C     K = MP(I) MEANS |V(I) - U(K)| = MIN
C     MP < 0 MEANS MORE THAN ONE I USES SAME K
C
C     T0 = C2 - 2*C1*COS(PIL*J)
C     U(KP) = T0 - 2*C0*COS(PIK*I)
C     KP = (J-1)*KX + I
C     C2 = ONE
C     C0 = X-CODIAGONAL = INPUT
C     C1 = Y-CODIAGONAL = HALF - C0
C-----------------------------------------------------------------------
      N = KX*KY
      ZERO = 0.0D0
      ONE  = 1.0D0
      TWO  = 2.0D0
C
C     SET G(I) = GAP FROM GOOD T(MEV) TO NEAREST DISTINCT TMEV I=1,NG
      DO 10 I = 1,NG
      G(I) = U(I)
   10 CONTINUE
C
C     REGENERATE A-EIGENVALUES
      T0 = DARCOS(-ONE)
      T1 = DFLOAT(KX+1)
      PIK = T0/T1
      T1 = DFLOAT(KY+1)
      PIL = T0/T1
      KP = 0
C
      DO 30 J = 1,KY
      T1 = PIL*DFLOAT(J)
      T0 = C2 - TWO*C1*DCOS(T1)
      DO 20 I = 1,KX
      KP = KP+1
      T1 = PIK*DFLOAT(I)
   20  U(KP) = T0 - TWO*C0*DCOS(T1)
   30 CONTINUE
C
C     ORDER U VECTOR BY INCREASING ALGEBRAIC SIZE
      DO 50 K = 2,N
      KM1 = K-1
      DO 40 L = 1,KM1
      JJ = K-L
      IF (U(JJ+1).GE.U(JJ)) GO TO 50
      T0 = U(JJ)
      U(JJ) = U(JJ+1)
   40 U(JJ+1) = T0
   50 CONTINUE
C
      ATOLN = DMAX1(DABS(U(1)),DABS(U(N)))*EPSM
C
C     DETERMINE MULTIPLICITIES FOR TRUE EIGENVALUES
      I = 1
      J = 1
      NEXACT = 0
   60 J = J+1
      IF (J.GT.N) GO TO 70
      EE = DABS(U(J)-U(I))
      IF (EE.GT.ATOLN) GO TO 70
      IDEX = IDEX+1
      GO TO 60
   70 NEXACT = NEXACT+1
      U(NEXACT) = U(I)
      I = J
      IF (I.GT.N) GO TO 80
      GO TO 60
   80 CONTINUE
C
C     NEXACT = NUMBER OF DISTINCT A-EIGENVALUES
C     U CONTAINS TRUE DISTINCT A-EV ORDERED BY INCREASING SIZE
C
      IF ( IWRITE.EQ.1) WRITE(6,90)MEV,NG,NEXACT
   90 FORMAT(/3I6,' = MEV, NG, NEXACT, POISZ CASE'/
     1 ' TRUE ERRORS FOR GOOD EIGENVALUES'/)
C
C      WRITE(6,61) (K,U(K), K=1,NEXACT)
C  61 FORMAT(4(I5,E15.8))
C
C     CALCULATION OF THE TRUE ERRORS.
      KL = 1
      DO 110 ITEV = 1,NG
      EV = V(ITEV)
      K = KL
      T1 = DABS(EV - U(KL))
C
      DO 100 KP = KL,NEXACT
      T0 = DABS(EV - U(KP))
      IF (T0.GE.T1) GO TO 100
      K = KP
      T1 = T0
  100 CONTINUE
C
      IF (K.EQ.KL.AND.ITEV.GT.1) T1 = -T1
      KL = K
      MP(ITEV) = K
      G(MEV+ITEV) = T1
  110 CONTINUE
C
C     TRUE ERRORS HAVE BEEN COMPUTED OUTPUT THEM TO FILE 10
C     FORM HEADER FOR ERREXACT FILE 10
      WRITE(10,120)N,KX,KY,C2,C0,C1
  120 FORMAT(' POISSONZ TRUE ERROR FOR GOOD EIGENVALUES'/
     1 5X,'N',4X,'NX',4X,'NY'/3I6//
     2 5X,'A-DIAGONAL',3X,'X-CODIAGONAL',3X,'Y-CODIAGONAL'/3E15.8//)
C
      WRITE(10,130)MEV,NG,NEXACT
  130 FORMAT(/3I6,' = MEV,NG,NEXACT'/1X,'T-EV NO',1X,'A-EV NO',
     1 10X,'GOOD EIGENVALUE',5X,'TRUEERR0R',7X,'TMINGAP')
C
      WRITE(10,140)(I,MP(I),V(I),G(MEV+I),G(I), I=1,NG)
  140 FORMAT(2I8,E25.16,2E14.3)
C
      WRITE(10,150)
  150 FORMAT(' ABOVE ARE THE TRUE ERRORS FOR POISSON GOODEV'/
     1 ' IF A-EV NO LT 0 THEN GOODEV HAS MULTIPLICITY GT 1'/
     1 ' IF TRUE ERROR LT 0 THEN MORE THAN ONE GOODEV APPROXIMATES'/
     1 ' THE SAME TRUE POISSON EIGENVALUE'/
     1 ' IF TMINGAP LT 0 THE MINGAP IS DUE TO SPURIOUS EIGENVALUE'//)
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY EXERRP(C0,C1,C2,KX,KY)
C
C-----------------------------------------------------------------------
C
C-----END OF EXERR------------------------------------------------------
      RETURN
      END
C
C-----START OF EXVEC----------------------------------------------------
C
C     (JVEC = 1): FOR A GIVEN RITZ VECTOR V AND EIGENVALUE X1, COMPUTES
C     THE CLOSEST EIGENVALUE Y1 AND CORRESPONDING TRUE EIGENVECTOR U,
C     AND THEN CALCULATES THE NORM OF THE DIFFERENCE BETWEEN
C     V AND U AND THE MAXIMAL DIFFERENCE BETWEEN THE COMPONENTS.
C     THESE QUANTITIES ARE WRITTEN TO FILE 6.
C
C     (JVEC = 2):  COMPUTES THE PROJECTION OF EACH
C     OF THE TRUE EIGENVECTORS ON THE LANCZOS STARTING VECTOR
C     USED BY THE LANCZS SUBROUTINE AND WRITES THEM TO FILE 12.
C
      SUBROUTINE EXVEC(U,V,X1,Y1,G,MP,IIX,JVEC,ICOUNT)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION U(1),V(1)
      DOUBLE PRECISION WI(110),WJ(110),WII(110)
      DOUBLE PRECISION X1,Y1,EV,EE,WS,PIK,PIL,SUM,PROJ,TEMP,S
      DOUBLE PRECISION ATOLN,EPSM,MACHEP,ZERO,HALF,ONE,TWO
      DOUBLE PRECISION C0,C1,C2,T0,T1,T2
      REAL G(1),GG
      INTEGER MP(1)
      DOUBLE PRECISION FINPRO
C-----------------------------------------------------------------------
C     THIS PROGRAM CALCULATES THE TRUE EIGENVALUES AND EIGENVECTORS
C     OF THE POISSON MATRIX A OF ORDER  N = KX*KY
C     A CONSISTS OF KY TRIDIAGONAL BLOCKS OF ORDER KX
C     KX = X-DIMENSION     KY = Y-DIMENSION.
C
C     IIX = SEED FOR RANDOM NUMBER GENERATOR USED TO CALCULATE
C           STARTING LANCZOS VECTOR IN LANCZS
C     V = RANDOM UNIT STARTING VECTOR FOR LANCZS
C     A*U = EV*U   ||U|| = ONE
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     U(K) IS A-EV ORDERED BY INCREASING SIZE, K = 1,N
C     LATER U IS USED TO STORE THE TRUE EIGENVECTOR
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     U(KV) = SIN(PIK*I*IK)*SIN(PIL*J*JK)
C     IK = 1,KX    JK = 1,KY    KV = (JK-1)*KX + IK
C     U IS UNSCALED EIGENVECTOR FOR EV(I,J) = Y1
C     WS = 1/||U||: ||U|| = .5*DSQRT(T2*T3)  T2 = KX+1  T3 = KY+1
C
C    JVEC  = (1,2) FLAGS COMPUTATIONS TO BE PERFORMED.
C
C          = (1) MEANS 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 U, A*U = Y1*U
C                     T2 = ||V-U||  T1 = MAX(|V(K)-U(K)|, K= 1,N)
C                     MAX OCCURS FIRST AT K = KK
C
C          = (2) MEANS CALCULATION OF THE PROJECTION OF THE STARTING
C                     LANCZOS VECTOR ON EACH EIGENVECTOR OF A.
C
C-----------------------------------------------------------------------
      DATA MACHEP/Z3410000000000000/
      EPSM = 2.0D0*MACHEP
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
      IF (ICOUNT.EQ.1) WRITE(6,60) N,KX,KY,JVEC,C2,C0,C1,ATOLN
   60 FORMAT(/' TRUE ERRORS FOR CONVERGED GOODEV'/
     1 4I6,' = N KX KY JVEC'//
     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
      IF (JVEC.EQ.1) GO TO 180
C
C     JVEC = 2 SO CALCULATE PROJECTIONS AND WRITE IN FILE 12
C
      WRITE(12,70)
   70 FORMAT(' PROJECTIONS OF LANCZOS STARTING VECTOR ON A-EIGENVECS')
C
      WRITE(12,80)N,KX,KY,IIX,C2,C0,C1,ATOLN
   80 FORMAT(1X,'A-SIZE',2X,'X-DIM',2X,'Y-DIM',6X,'SVSEED'/3I7,I12/
     1 5X,'A-DIAGONAL',3X,'X-CODIAGONAL',3X,'Y-CODIAGONAL',5X,'ATOLN'/
     2 3E15.8,E10.3)
C
      WRITE(12,90)
   90 FORMAT(5X,'PROJECTION',8X,'TRUE A-EIGENVALUE',1X,'EV NO'
     1,2X,'VEC NO')
C
C     GENERATE SAME RANDOM UNIT VECTOR USED IN THE LANCZS RECURSIONS.
      IIL=IIX
C
C-----------------------------------------------------------------------
      CALL GENRAN(IIL,G,N)
C-----------------------------------------------------------------------
C
      DO 100 I = 1,N
  100 V(I) = G(I)
C
C-----------------------------------------------------------------------
      SUM = FINPRO(N,V(1),1,V(1),1)
C-----------------------------------------------------------------------
C
      SUM = 1.D0/DSQRT(SUM)
C
      DO 110 I = 1,N
  110 V(I) = V(I)*SUM
C
C     DETERMINE UNIT EIGENVECTOR W ASSOCIATED WITH EACH EV(I,J) = Y1
C     AND CALCULATE THE PROJECTION G(K) OF U ON THE STARTING VECTOR V
C     A*U = EV*U       WS = 1/||WU||: WU = UNSCALED EIGENVECTOR
C
      DO 160 K =1,N
C     DETERMINE I J FROM K: MP(K) = KP = (J-1)*KX+I
      KP = MP(K)
      I = MOD(KP,KX)
      IF (I.EQ.0) I = KX
      T1 = WI(I)
      J = 1 + (KP-1)/KX
      T0 = WJ(J)
      T0 = WJ(J)
C
      Y1 = C2 - TWO*C1*DCOS(WJ(J)) - TWO*C0*DCOS(WI(I))
C     Y1 = EV(I,J)
C
      DO 120 II = 1,KX
      T2 = T1*DFLOAT(II)
  120 WII(II) = WS*DSIN(T2)
C
      KV = 0
      DO 140 JJ = 1,KY
      T2 = T0*DFLOAT(JJ)
      T2 = DSIN(T2)
C
      DO 130 II = 1,KX
      KV = KV + 1
  130 U(KV) = T2*WII(II)
C
  140 CONTINUE
C
C     U IS UNIT EIGENVECTOR OF A ASSOCIATED WITH EV(I,J) = Y1
C     G(K) IS THE PROJECTION OF U ON V  FOR Y1
C
C-----------------------------------------------------------------------
      PROJ = FINPRO(N,U(1),1,V(1),1)
C-----------------------------------------------------------------------
C
      TEMP = DABS(PROJ)
      G(K) = TEMP
C
C     DESIRED PROJECTION HAS BEEN COMPUTED OUTPUT IT TO FILE 12.
      WRITE(12,150) G(K),Y1,K,MP(K)
  150 FORMAT(E15.8,E25.16,I6,I8)
C
  160 CONTINUE
C
      WRITE(12,170)
  170 FORMAT(' ----- END OF FILE 12 PROJECT -----------------------'//)
C
      GO TO 310
C
C     JVEC = 1
C
C     X1 IS AN INPUT PARAMETER. WE CALCULATE TRUE
C     A-EIGENVALUE WHICH IS CLOSEST TO X1, LABEL IT Y1 AND CALCULATE
C     UNIT EIGENVECTOR OF A ASSOCIATED WITH Y1. A*U = Y1*U, ||U|| = 1.
C     Y1 = EV(I,J). EIGENVALUES OF A ARE ORDERED BY INCREASING SIZE.
C     V = RITZ VECTOR ASSOCIATED WITH GOODEV X1
C
  180 CONTINUE
      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 200
C
      DO 190 KVEC = 2,N
      IF (X1.GE.U(KVEC)) GO TO 190
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 200
  190 CONTINUE
C
  200 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 210
      EE = DMIN1(U(KX1+1)-U(KX1),U(KX1)-U(KX1-1))
  210 CONTINUE
C
      T0 = DABS(ONE - X1/Y1)
C
      WRITE(6,220) N,KX1,ICOUNT,Y1,X1,T0,EE
  220 FORMAT(3I8,' = N, A-EV NUMBER,RITZ NUMBER'//
     1 18X,' TRUEEV',19X,'GOODEV',4X,'RELERROR',4X,'A-MINGAP'/
     1 2E25.16,2E12.3/)
C
      IF (EE.GT.ATOLN) GO TO 240
C
      WRITE(6,230)
  230 FORMAT(' Y1 IS A MULTIPLE EIGENVALUE OF A SO WE EXIT'/)
C
      GO TO 310
C
C     Y1 IS TOEPLITZ EIGENVALUE CLOSEST TO X1.
C     CALCULATION OF EIGENVECTOR ASSOCIATED WITH EIGENVALUE Y1
C     A*U = Y1*U
C     DETERMINE I J FROM K: MP(K) = KP = (J-1)*KX+I
  240 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 250 II = 1,KX
      T0 = T1*DFLOAT(II)
  250 WII(II) = WS*DSIN(T0)
C
      KV = 0
      DO 270 JJ = 1,KY
      T0 = T2*DFLOAT(JJ)
      T0 = DSIN(T0)
C
      DO 260 II = 1,KX
      KV = KV + 1
  260 U(KV) = T0*WII(II)
C
  270 CONTINUE
C
C     U IS UNIT TRUE EIGENVECTOR OF A ASSOCIATED WITH Y1
C     V IS UNIT RITZVECTOR OF A ASSOCIATED WITH X1
C
      KK = 0
      S = ONE
      T1 = ZERO
C
      DO 280 K = 1,N
      IF (DABS(U(K)).LE.T1) GO TO 280
      T1 = DABS(U(K))
      KK = K
  280 CONTINUE
      IF (U(KK)*V(KK).LT.ZERO) S = - ONE
C
      KK = 0
      T1 = ZERO
      T2 = ZERO
      DO 290 K = 1,N
      TEMP = DABS(S*U(K) - V(K))
      T2 = T2 + TEMP**2
      IF (TEMP.LE.T1) GO TO 290
      KK = K
      T1 = TEMP
  290 CONTINUE
C
      T2 = DSQRT(T2)
      WRITE(6,300) KK,T1,T2
  300 FORMAT(' EIGENVECTOR ERROR. MAX ERROR AT COMPONENT = ',I6/
     1 ' MAX DABS(TRUEVEC(K)-RITZVEC(K)) = ',E12.5/
     1 ' NORM(TRUEVEC-RITZVEC)  = ',E12.5/)
C
  310 CONTINUE
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY EXVECP(C0,C1,C2,KX,KY)
C
C-----------------------------------------------------------------------
C
      RETURN
C-----END OF EXVEC------------------------------------------------------
      END

.