C  VERSION 2   DOES NOT USE EISPACK
C
C ------------------------------------------------------------------
C
      SUBROUTINE DILASO(OP, IOVECT, N, XL, XR, NFIG, NPERM,
     *   NMVAL, VAL, NMVEC, MAXVEC, VEC, NBLOCK, MAXOP, MAXJ,
     *   WORK, IND, IERR)
C
      INTEGER N, NFIG, NPERM, NMVAL, NMVEC, MAXVEC, NBLOCK,
     *   MAXOP, MAXJ, IND(1), IERR
      DOUBLE PRECISION XL, XR, VEC(NMVEC,1), VAL(NMVAL,1), WORK(1)
      EXTERNAL OP, IOVECT
C
C AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT
C
C COMPUTER SCIENCES DEPARTMENT
C UNIVERSITY OF TEXAS AT AUSTIN
C AUSTIN, TX 78712
C
C VERSION 2 ORIGINATED APRIL 1982
C
C CURRENT VERSION  JUNE 1983
C
C DILASO FINDS ALL THE EIGENVALUES AND EIGENVECTORS OF A LARGE
C SPARSE SYMMETRIC MATRIX OUTSIDE A USER DEFINED EXCLUDED INTERVAL.
C THE SUBROUTINE DILASO IS PRIMARILY A DRIVER FOR SUBROUTINE DIWLA
C WHICH IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE
C ORTHOGONALIZATION AND SUBROUTINE DIPPLA WHICH POST PROCESSES
C THE OUTPUT OF DIWLA.  HOWEVER DILASO DOES CHECK FOR INCONSISTENCIES
C IN THE CALLING PARAMETERS AND DOES PREPROCESS AND USER SUPPLIED
C EIGENPAIRS.
C
C
C ON INPUT
C
C
C   OP   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE
C     OP(N,M,P,Q).  P AND Q ARE N X M MATRICES AND Q IS
C     RETURNED AS THE MATRIX TIMES P.
C
C   IOVECT   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE
C     IOVECT(N,M,Q,J,K).  Q IS AN N X M MATRIX.  IF K = 0
C     THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH
C     THE JTH LANCZOS VECTORS.  IF K = 1 THEN Q IS RETURNED
C     AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS.  SEE
C     DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES.
C
C   N   THE ORDER OF THE MATRIX.
C
C   XL   THE LEFT ENDPOINT OF THE EXCLUDED INTERVAL.
C
C   XR   THE RIGHT ENDPOINT OF THE EXCLUDED INTERVAL.
C
C   NFIG   THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE
C     EIGENVALUES.  NFIG MUST BE GREATER THAN OR EQUAL TO 1.
C
C   NPERM   AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER
C     SUPPLIED EIGENPAIRS.  IN MOST CASES NPERM WILL BE ZERO.  SEE
C     DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER
C     THAN ZERO.  NPERM MUST NOT BE LESS THAN ZERO.
C
C   NMVAL   THE ROW DIMENSION OF THE ARRAY VAL.  NMVAL MUST BE GREATER
C     THAN OR EQUAL TO MAXVEC.
C
C   VAL   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW
C     DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4.  IF NPERM
C     IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED
C     IN VAL.  SEE DOCUMENTATION FOR DETAILS.
C
C   NMVEC   THE ROW DIMENSION OF THE ARRAY VEC.  NMVEC MUST BE GREATER
C     THAN OR EQUAL TO N.
C
C   MAXVEC   THE MAXIMUM NUMBER OF EIGENPAIRS TO BE DETERMINED.
C     MAXVEC MUST NOT EXCEED THE COLUMN DIMENSION OF THE ARRAY VEC.
C
C   VEC   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW
C     DIMENSION NMVEC AND COLUMN DIMENSION MAXVEC.  IF NPERM > 0
C     THEN THE FIRST NPERM COLUMNS OF VEC MUST CONTAIN
C     THE USER SUPPLIED EIGENVECTORS.
C
C   NBLOCK   THE BLOCK SIZE.  SEE DOCUMENTATION FOR CHOOSING
C     AN APPROPRIATE VALUE FOR NBLOCK.  NBLOCK MUST BE GREATER
C     THAN ZERO .AND. LESS THAN  MAXJ/6.
C
C   MAXOP   AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE
C     OP.  DILASO TERMINATES WHEN MAXOP IS EXCEEDED.  SEE
C     DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP.
C
C   MAXJ   AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND
C     DOCUMENTATION ON IOVECT).  FOR THE FASTEST CONVERGENCE MAXJ
C     SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE
C     MAXJ LARGER THAN MAXOP*NBLOCK.
C
C   WORK   A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS
C     LARGE AS
C
C     2*N*NBLOCK + MAXJ*(NBLOCK+1+MAXVEC+2) + 2*NBLOCK*NBLOCK + 3*MAXVEC
C
C     + THE MAXIMUM OF
C          N*NBLOCK
C             AND
C          MAXJ*(2*NBLOCK+2) + 2*MAXVEC + 8 + (2*NBLOCK+2)*(NBLOCK+1)
C
C     THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED
C     STARTING VECTORS.  SEE DOCUMENTATION FOR GUIDELINES IN
C     CHOOSING STARTING VECTORS.
C
C   IND   AN INTEGER ARRAY OF DIMENSION AT LEAST 4*MAXVEC.
C
C   IERR   AN INTEGER VARIABLE.
C
C
C ON OUTPUT
C
C
C   NPERM   THE NUMBER OF EIGENPAIRS NOW KNOWN.
C
C   VEC   THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS.
C
C   VAL   THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING
C     EIGENVALUES.  THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF
C     THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN-
C     VALUES.  THE THIRD COLUMN CONTAINS MORE DOUBLE PRECISIONISTIC ESTIMATES
C     OF THE ACCURACY OF THE EIGENVALUES.  THE FOURTH COLUMN CONTAIN
C     ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS.  SEE
C     DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES.
C
C   WORK   IF EXECUTION IS TERMINATED BEFORE ALL THE DESIRED
C     EIGENVALUES ARE FOUND (IERR = -2, -5, OR -6 MOD(10) ) THEN
C     THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS
C     FOR RESTARTING THE ALGORITHM AND DILASO CAN BE IMMEDIATELY
C     RECALLED TO CONTINUE WORKING ON THE PROBLEM.
C
C   IND   IND(1) CONTAINS THE ACTUAL NUMBER OF CALLS TO OP.  ON SOME
C     OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER
C     THAN MAXOP.
C
C   IERR   AN ERROR COMPLETION CODE.  THE NORMAL COMPLETION CODE IS
C     ZERO.  SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO
C     COMPLETION CODES.
C
C
C INTERNAL VARIABLES.
C
C
      INTEGER I, I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11,
     *   I12, I13, M, NBAND, NOP, NP
      LOGICAL RARITZ
      DOUBLE PRECISION DELTAL, DELTAR, EPS, TEMP, DDOT, DNRM2, DABS
      EXTERNAL DIPPLA, DIWLA, DMVPC, DORTQR, DAXPY,
     *   DCOPY, DDOT, DNRM2, DSCAL, DSWAP, DLAEIG, DLABCM, DLABFC,
     *   DLAGER, DLARAN, DVSORT
C
C NOP   RETURNED FROM DIWLA AS THE NUMBER OF CALLS TO THE
C   SUBROUTINE OP.
C
C RARITZ   RETURNED FROM DIWLA AS .TRUE.  IF A FINAL RAYLEIGH-RITZ
C   PROCEDURE IS NEEDED.
C
C DELTAL   RETURNED FROM DIWLA AS THE EXCLUDED EIGENVALUE CLOSEST
C   TO XL.
C
C DELTAR   RETURNED FROM DIWLA AS THE EXCLUDED EIGENVALUE CLOSEST
C   TO XR.
C
C DIPPLA   A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED
C   BY DIWLA.
C
C DIWLA   A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM
C   WITH SELECTIVE ORTHOGONALIZATION.
C
C DMVPC   A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND
C   ORTHOGONALITY COEFFICIENT OF A GIVEN RITZ VECTOR.
C
C DORTQR   A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS
C   USING HOUSEHOLDER REFLECTIONS.
C
C DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP   A SUBSET OF THE BASIC LINEAR
C   ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION.
C
C DLAEIG, DLABCM, DLABFC, DLAGER  USED FOR BAND EIGENVALUE PROBLEMS
C
C DLARAN  SUBROUTINE WHICH RETURNS RANDOM VECTORS
C
C ------------------------------------------------------------------
C
C THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS.
C
      IND(1) = 0
      IERR = 0
C
      IF (N.LT.6*NBLOCK) IERR = 1
      IF (NFIG.LE.0) IERR = IERR + 2
      IF (NMVEC.LT.N) IERR = IERR + 4
      IF (NPERM.LT.0) IERR = IERR + 8
      IF (MAXJ.LT.6*NBLOCK) IERR = IERR + 16
      IF (MAXVEC.LT.NPERM) IERR = IERR + 32
      IF (MAXVEC.GT.NMVAL) IERR = IERR + 64
      IF (MAXVEC.GT.MAXOP) IERR = IERR + 128
      IF (XL.GT.XR) IERR = IERR + 256
      IF (NBLOCK.LT.1) IERR = IERR + 512
      IF (IERR.NE.0) RETURN
C
C ------------------------------------------------------------------
C
C THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS.
C IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION
C OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO  -1
C AND DILASO TERMINATES.  IF A USER SUPPLIED EIGENVALUE LIES INSIDE
C THE EXCLUDED INTERVAL THEN IERR IS SET TO  -4 AND DILASO TERMINATES.
C
      IF (NPERM.EQ.0) GO TO 100
C
C THIS SORTS THE USER SUPPLIED VALUES AND VECTORS.
C
      CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TEMP, NMVEC, N, VEC)
C
C THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON.
C IT ALSO INSURES THAT THE REIDUAL NORMS ARE POSITIVE.
C
      DO 40 I=1,NPERM
         VAL(I,2) = DABS(VAL(I,2))
         VAL(I,3) = DNRM2(N, VEC(1,I), 1)
   40 CONTINUE
C
C THIS CHECKS FOR VALUES OUT OF RANGE
C
      DO 50 I=1,NPERM
         IF (VAL(I,1).LE.XL .OR. VAL(I,1).GE.XR) GO TO 50
         IERR = -4
         RETURN
C
   50 CONTINUE
C
C THIS ORTHONORMALIZES THE VECTORS
C
      M = N*NBLOCK + 1
      CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M))
      M = N*NBLOCK - NPERM
      DO 60 I = 1, NPERM
         M = M + NPERM + 1
         IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 60
         IERR = -1
         RETURN
C
   60 CONTINUE
C
C THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN
C THE ARRAY WORK FOR LATER REFERENCE IN DIWLA.
C
      M = 2*N*NBLOCK + 1
      CALL DCOPY(NPERM, VAL(1,2), 1, WORK(M), 1)
C
C THIS SETS EPS TO AN APPROXIMATION OF THE MACHINE PRECISION
C
C ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT
C ***IN A PRODUCTION CODE
C
  100 EPS = 1.0D0
      DO 110 I =1,1000
         EPS = 0.5D0*EPS
         TEMP = 1.0D0 + EPS
         IF(TEMP.EQ.1.0D0) GO TO 120
  110 CONTINUE
C
C ------------------------------------------------------------------
C
C THIS SECTION CALLS DIWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM
C WITH SELECTIVE ORTHOGONALIZATION.
C
  120 NBAND = NBLOCK + 1
      I1 = 1 + N*NBLOCK
      I2 = I1 + N*NBLOCK
      I3 = I2 + MAXVEC
      I4 = I3 + MAXVEC
      I5 = I4 + MAXVEC
      I6 = I5 + MAXJ*NBAND
      I7 = I6 + NBLOCK*NBLOCK
      I8 = I7 + NBLOCK*NBLOCK
      I9 = I8 + MAXJ*(MAXVEC+2)
      I10 = I9 + NBLOCK
      I11 = I10 + 2*MAXVEC + 8
      I12 = I11 + MAXJ*(2*NBLOCK+1)
      I13 = I12 + MAXJ
      CALL DIWLA(OP, IOVECT, N, NBAND, XL, XR, NFIG, NPERM, VAL,
     *   NMVEC, MAXVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP,
     *   WORK(1), WORK(I1), WORK(I2), WORK(I3), WORK(I4),
     *   WORK(I5), WORK(I6), WORK(I7), WORK(I8), WORK(I9),
     *   WORK(I10), WORK(I11), WORK(I12), WORK(I13),
     *   IND, RARITZ, DELTAL, DELTAR, EPS, IERR)
C
C ------------------------------------------------------------------
C
C THIS SECTION CALLS DIPPLA (THE POST PROCESSOR).  NP IS USED
C TO DIMENSION AN ARRAY IN DIPPLA.  IT IS NEEDED TO LEGALIZE
C CHANGING THE VALUE OF NPERM, WHICH MAY OCCUR IN DIPPLA.
C
      IF (NPERM.EQ.0) GO TO 130
      I1 = N*NBLOCK + 1
      I2 = I1 + NPERM*NPERM
      I3 = I2 + NPERM*NPERM
      I4 = I3 + MAX0(N*NBLOCK, 2*NPERM*NPERM)
      I5 = I4 + N*NBLOCK
      I6 = I5 + 2*NPERM + 4
      NP = NPERM
      CALL DIPPLA(OP, IOVECT, N, XL, XR, NP, NPERM, NOP, NMVAL,
     *   VAL, NMVEC, VEC, NBLOCK, WORK(I1), WORK(I2),
     *   WORK(I3), WORK(I4), WORK(I5), WORK(I6), DELTAL,
     *   DELTAR, RARITZ, EPS, IERR)
C
  130 IND(1) = NOP
      RETURN
      END