1C   VERSION 2    DOES NOT USE EISPACK
2C
3C ------------------------------------------------------------------
4C
5      SUBROUTINE DNLASO(OP, IOVECT, N, NVAL, NFIG, NPERM,
6     *   NMVAL, VAL, NMVEC, VEC, NBLOCK, MAXOP, MAXJ, WORK,
7     *   IND, IERR)
8C
9      INTEGER N, NVAL, NFIG, NPERM, NMVAL, NMVEC, NBLOCK,
10     *   MAXOP, MAXJ, IND(1), IERR
11      DOUBLE PRECISION VEC(NMVEC,1), VAL(NMVAL,4), WORK(1)
12      EXTERNAL OP, IOVECT
13C
14C AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT
15C
16C COMPUTER SCIENCES DEPARTMENT
17C UNIVERSITY OF TEXAS AT AUSTIN
18C AUSTIN, TX 78712
19C
20C VERSION 2 ORIGINATED APRIL 1982
21C
22C CURRENT VERSION  JUNE 1983
23C
24C DNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF
25C THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX.  THE SUBROUTINE
26C DNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE DNWLA WHICH IMPLEMENTS
27C THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND
28C SUBROUTINE DNPPLA WHICH POST PROCESSES THE OUTPUT OF DNWLA.
29C HOWEVER DNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING
30C PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS.
31C DNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES.  IF
32C THE LARGEST EIGENVALUES ARE DESIRED DNLASO IMPLICITLY USES THE
33C NEGATIVE OF THE MATRIX.
34C
35C
36C ON INPUT
37C
38C
39C   OP   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE
40C     OP(N,M,P,Q).  P AND Q ARE N X M MATRICES AND Q IS
41C     RETURNED AS THE MATRIX TIMES P.
42C
43C   IOVECT   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE
44C     IOVECT(N,M,Q,J,K).  Q IS AN N X M MATRIX.  IF K = 0
45C     THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH
46C     THE JTH LANCZOS VECTORS.  IF K = 1 THEN Q IS RETURNED
47C     AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS.  SEE
48C     DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES.
49C
50C   N   THE ORDER OF THE MATRIX.
51C
52C   NVAL   NVAL SPECIFIES THE EIGENVALUES TO BE FOUND.
53C     DABS(NVAL)  IS THE NUMBER OF EIGENVALUES DESIRED.
54C     IF NVAL < 0 THE ALGEBRAICALLY SMALLEST (LEFTMOST)
55C     EIGENVALUES ARE FOUND.  IF NVAL > 0 THE ALGEBRAICALLY
56C     LARGEST (RIGHTMOST) EIGENVALUES ARE FOUND.  NVAL MUST NOT
57C     BE ZERO.  DABS(NVAL) MUST BE LESS THAN  MAXJ/2.
58C
59C   NFIG   THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE
60C     EIGENVALUES.  NFIG MUST BE GREATER THAN OR EQUAL TO 1.
61C
62C   NPERM   AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER
63C     SUPPLIED EIGENPAIRS.  IN MOST CASES NPERM WILL BE ZERO.  SEE
64C     DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER
65C     THAN ZERO.  NPERM MUST NOT BE LESS THAN ZERO.
66C
67C   NMVAL   THE ROW DIMENSION OF THE ARRAY VAL.  NMVAL MUST BE GREATER
68C     THAN OR EQUAL TO DABS(NVAL).
69C
70C   VAL   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW
71C     DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4.  IF NPERM
72C     IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED
73C     IN VAL.  SEE DOCUMENTATION FOR DETAILS.
74C
75C   NMVEC   THE ROW DIMENSION OF THE ARRAY VEC.  NMVEC MUST BE GREATER
76C     THAN OR EQUAL TO N.
77C
78C   VEC   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW
79C     DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST DABS(NVAL).  IF
80C     NPERM > 0 THEN THE FIRST NPERM COLUMNS OF VEC MUST
81C     CONTAIN THE USER SUPPLIED EIGENVECTORS.
82C
83C   NBLOCK   THE BLOCK SIZE.  SEE DOCUMENTATION FOR CHOOSING
84C     AN APPROPRIATE VALUE FOR NBLOCK.  NBLOCK MUST BE GREATER
85C     THAN ZERO AND LESS THAN  MAXJ/6.
86C
87C   MAXOP   AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE
88C     OP.  DNLASO TERMINATES WHEN MAXOP IS EXCEEDED.  SEE
89C     DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP.
90C
91C   MAXJ   AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND
92C     DOCUMENTATION ON IOVECT).  FOR THE FASTEST CONVERGENCE MAXJ
93C     SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE
94C     MAXJ LARGER THAN MAXOP*NBLOCK.
95C
96C   WORK   A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS
97C     LARGE AS
98C
99C         2*N*NBLOCK + MAXJ*(NBLOCK+NV+2) + 2*NBLOCK*NBLOCK + 3*NV
100C
101C            + THE MAXIMUM OF
102C                 N*NBLOCK
103C                   AND
104C         MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1)
105C
106C     WHERE NV = DABS(NVAL)
107C
108C     THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED
109C     STARTING VECTORS.  SEE DOCUMENTATION FOR GUIDELINES IN
110C     CHOOSING STARTING VECTORS.
111C
112C   IND   AN INTEGER ARRAY OF DIMENSION AT LEAST DABS(NVAL).
113C
114C   IERR   AN INTEGER VARIABLE.
115C
116C
117C ON OUTPUT
118C
119C
120C   NPERM   THE NUMBER OF EIGENPAIRS NOW KNOWN.
121C
122C   VEC   THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS.
123C
124C   VAL   THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING
125C     EIGENVALUES.  THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF
126C     THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN-
127C     VALUES.  THE THIRD COLUMN CONTAINS MORE DOUBLE PRECISIONISTIC ESTIMATES
128C     OF THE ACCURACY OF THE EIGENVALUES.  THE FOURTH COLUMN CONTAINS
129C     ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS.  SEE
130C     DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES.
131C
132C   WORK   IF WORK IS TERMINATED BEFORE COMPLETION (IERR = -2)
133C     THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS
134C     FOR RESTARTING THE ALGORITHM AND DNLASO CAN BE IMMEDIATELY
135C     RECALLED TO CONTINUE WORKING ON THE PROBLEM.
136C
137C   IND   IND(1)  CONTAINS THE ACTUAL NUMBER OF CALLS TO OP.  ON SOME
138C     OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER
139C     THAN MAXOP.
140C
141C   IERR   AN ERROR COMPLETION CODE.  THE NORMAL COMPLETION CODE IS
142C     ZERO.  SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO
143C     COMPLETION CODES.
144C
145C
146C INTERNAL VARIABLES.
147C
148C
149      INTEGER I, I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11,
150     *  I12, I13, M, NBAND, NOP, NV, IABS, MAX0
151      LOGICAL RARITZ, SMALL
152      DOUBLE PRECISION DELTA, EPS, TEMP, DNRM2, DABS, TARR(1)
153      EXTERNAL DNPPLA, DNWLA, DORTQR, DCOPY, DNRM2, DVSORT
154C
155C NOP   RETURNED FROM DNWLA AS THE NUMBER OF CALLS TO THE
156C   SUBROUTINE OP.
157C
158C NV   SET EQUAL TO DABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED,
159C   AND PASSED TO DNWLA.
160C
161C SMALL   SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED.
162C
163C RARITZ   RETURNED FROM DNWLA AND PASSED TO DNPPLA.  RARITZ IS .TRUE.
164C   IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED.
165C
166C DELTA   RETURNED FROM DNWLA AS THE EIGENVALUE OF THE MATRIX
167C   WHICH IS CLOSEST TO THE DESIRED EIGENVALUES.
168C
169C DNPPLA   A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED
170C   BY DNWLA.
171C
172C DNWLA   A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM
173C   WITH SELECTIVE ORTHOGONALIZATION.
174C
175C DMVPC   A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND
176C   ORTHOGONALITY COEFFICIENT OF GIVEN RITZ VECTORS.
177C
178C DORTQR   A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS
179C   USING HOUSEHOLDER REFLECTIONS.
180C
181C DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP   A SUBSET OF THE BASIC LINEAR
182C   ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION.
183C
184C DLARAN   A SUBROUTINE TO GENERATE RANDOM VECTORS
185C
186C DLAEIG, DLAGER, DLABCM, DLABFC   SUBROUTINES FOR BAND EIGENVALUE
187C   CALCULATIONS.
188C
189C ------------------------------------------------------------------
190C
191C THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS.
192C
193      NV = IABS(NVAL)
194      IND(1) = 0
195      IERR = 0
196      IF (N.LT.6*NBLOCK) IERR = 1
197      IF (NFIG.LE.0) IERR = IERR + 2
198      IF (NMVEC.LT.N) IERR = IERR + 4
199      IF (NPERM.LT.0) IERR = IERR + 8
200      IF (MAXJ.LT.6*NBLOCK) IERR = IERR + 16
201      IF (NV.LT.MAX0(1,NPERM)) IERR = IERR + 32
202      IF (NV.GT.NMVAL) IERR = IERR + 64
203      IF (NV.GT.MAXOP) IERR = IERR + 128
204      IF (NV.GE.MAXJ/2) IERR = IERR + 256
205      IF (NBLOCK.LT.1) IERR = IERR + 512
206      IF (IERR.NE.0) RETURN
207C
208      SMALL = NVAL.LT.0
209C
210C ------------------------------------------------------------------
211C
212C THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS.
213C IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION
214C OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO  -1
215C AND DNLASO TERMINATES.
216C
217      IF (NPERM.EQ.0) GO TO 110
218C
219C THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST
220C EIGENVALUES ARE DESIRED, SINCE DNWLA WILL IMPLICITLY USE THE
221C NEGATIVE OF THE MATRIX.
222C
223      IF (SMALL) GO TO 20
224      DO 10 I=1,NPERM
225         VAL(I,1) = -VAL(I,1)
226   10 CONTINUE
227C
228C THIS SORTS THE USER SUPPLIED VALUES AND VECTORS.
229C
230   20 CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TARR, NMVEC, N, VEC)
231C
232C THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON.
233C IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE.
234C
235      DO 60 I=1,NPERM
236         VAL(I,2) = DABS(VAL(I,2))
237         VAL(I,3) = DNRM2(N,VEC(1,I),1)
238   60 CONTINUE
239C
240C THIS PERFORMS THE ORTHONORMALIZATION.
241C
242      M = N*NBLOCK + 1
243      CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M))
244      M = N*NBLOCK - NPERM
245      DO 70 I = 1, NPERM
246         M = M + NPERM + 1
247         IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 70
248         IERR = -1
249         RETURN
250C
251   70 CONTINUE
252C
253C THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN
254C THE ARRAY WORK FOR LATER REFERENCE IN DNWLA.
255C
256      M = 2*N*NBLOCK + 1
257      CALL DCOPY(NPERM, VAL(1,2), 1, WORK(M), 1)
258C
259C THIS SETS EPS TO AN APPROXIMATION OF THE RELATIVE MACHINE
260C PRECISION
261C
262C ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT
263C ***IN A PRODUCTION CODE
264C
265  110 EPS = 1.0D0
266      DO 120 I = 1,1000
267         EPS = 0.5D0*EPS
268         TEMP = 1.0D0 + EPS
269         IF(TEMP.EQ.1.0D0) GO TO 130
270  120 CONTINUE
271C
272C ------------------------------------------------------------------
273C
274C THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM
275C WITH SELECTIVE ORTHOGONALIZATION.
276C
277  130 NBAND = NBLOCK + 1
278      I1 = 1 + N*NBLOCK
279      I2 = I1 + N*NBLOCK
280      I3 = I2 + NV
281      I4 = I3 + NV
282      I5 = I4 + NV
283      I6 = I5 + MAXJ*NBAND
284      I7 = I6 + NBLOCK*NBLOCK
285      I8 = I7 + NBLOCK*NBLOCK
286      I9 = I8 + MAXJ*(NV+1)
287      I10 = I9 + NBLOCK
288      I11 = I10 + 2*NV + 6
289      I12 = I11 + MAXJ*(2*NBLOCK+1)
290      I13 = I12 + MAXJ
291      CALL DNWLA(OP, IOVECT, N, NBAND, NV, NFIG, NPERM, VAL, NMVEC,
292     *   VEC, NBLOCK, MAXOP, MAXJ, NOP, WORK(1), WORK(I1),
293     *   WORK(I2), WORK(I3), WORK(I4), WORK(I5), WORK(I6),
294     *   WORK(I7), WORK(I8), WORK(I9), WORK(I10), WORK(I11),
295     *   WORK(I12), WORK(I13), IND, SMALL, RARITZ, DELTA, EPS, IERR)
296C
297C ------------------------------------------------------------------
298C
299C THIS SECTION CALLS DNPPLA (THE POST PROCESSOR).
300C
301      IF (NPERM.EQ.0) GO TO 140
302      I1 = N*NBLOCK + 1
303      I2 = I1 + NPERM*NPERM
304      I3 = I2 + NPERM*NPERM
305      I4 = I3 + MAX0(N*NBLOCK,2*NPERM*NPERM)
306      I5 = I4 + N*NBLOCK
307      I6 = I5 + 2*NPERM + 4
308      CALL DNPPLA(OP, IOVECT, N, NPERM, NOP, NMVAL, VAL, NMVEC,
309     *  VEC, NBLOCK, WORK(I1), WORK(I2), WORK(I3), WORK(I4),
310     *  WORK(I5), WORK(I6), DELTA, SMALL, RARITZ, EPS)
311C
312  140 IND(1) = NOP
313      RETURN
314      END
315C
316C ------------------------------------------------------------------
317C
318      SUBROUTINE DNWLA(OP, IOVECT, N, NBAND, NVAL, NFIG, NPERM, VAL,
319     *   NMVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP, P1, P0,
320     *   RES, TAU, OTAU, T, ALP, BET, S, P2, BOUND, ATEMP, VTEMP, D,
321     *   IND, SMALL, RARITZ, DELTA, EPS, IERR)
322C
323      INTEGER N, NBAND, NVAL, NFIG, NPERM, NMVEC, NBLOCK, MAXOP, MAXJ,
324     *   NOP, IND(1), IERR
325      LOGICAL RARITZ, SMALL
326      DOUBLE PRECISION VAL(1), VEC(NMVEC,1), P0(N,1), P1(N,1),
327     *   P2(N,1), RES(1), TAU(1), OTAU(1), T(NBAND,1),
328     *   ALP(NBLOCK,1), BET(NBLOCK,1), BOUND(1), ATEMP(1),
329     *   VTEMP(1), D(1), S(MAXJ,1), DELTA, EPS
330      EXTERNAL OP, IOVECT
331C
332C DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE
333C ORTHOGONALIZATION.
334C
335C   NBAND  NBLOCK + 1  THE BAND WIDTH OF T.
336C
337C   NVAL   THE NUMBER OF DESIRED EIGENVALUES.
338C
339C   NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS
340C     INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE
341C     ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL
342C     TO THE CURRENT KRYLOV SUBSPACE.
343C
344C   NOP   THE NUMBER OF CALLS TO OP.
345C
346C   P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS.
347C
348C   RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS.
349C
350C   TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION.
351C
352C   T   THE BAND MATRIX.
353C
354C   ALP   THE CURRENT DIAGONAL BLOCK.
355C
356C   BET   THE CURRENT OFF DIAGONAL BLOCK.
357C
358C   BOUND, ATEMP, VTEMP, D  TEMPORARY STORAGE USED BY THE BAND
359C      EIGENVALUE SOLVER DLAEIG.
360C
361C   S   EIGENVECTORS OF T.
362C
363C   SMALL   .TRUE.  IF THE SMALL EIGENVALUES ARE DESIRED.
364C
365C   RARITZ  RETURNED AS  .TRUE.  IF A FINAL RAYLEIGH-RITZ PROCEDURE
366C     IS TO BE DONE.
367C
368C   DELTA   RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE
369C     OF THE MATRIX.  USED IN ESTIMATING THE ACCURACY OF THE
370C     COMPUTED EIGENVALUES.
371C
372C
373C INTERNAL VARIABLES USED
374C
375      INTEGER I, I1, II, J, K, L, M, NG, NGOOD,
376     *   NLEFT, NSTART, NTHETA, NUMBER, MIN0, MTEMP
377      LOGICAL ENOUGH, TEST
378      DOUBLE PRECISION ALPMAX, ALPMIN, ANORM, BETMAX, BETMIN,
379     *   EPSRT, PNORM, RNORM, TEMP,
380     *   TMAX, TMIN, TOLA, TOLG, UTOL, DABS,
381     *   DMAX1, DMIN1, DSQRT, DDOT, DNRM2, TARR(1), DZERO(1)
382      EXTERNAL DMVPC, DORTQR, DAXPY, DCOPY, DDOT,
383     *   DNRM2, DSCAL, DLAEIG, DLAGER, DLARAN, DVSORT
384C
385C J   THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT
386C    KRYLOV SUBSPACE.
387C
388C NGOOD   THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS
389C    LIE IN THE CURRENT KRYLOV SUBSPACE).
390C
391C NLEFT   THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED,
392C    I.E.  NLEFT = NVAL - NPERM.
393C
394C NUMBER = NPERM + NGOOD.
395C
396C ANORM   AN ESTIMATE OF THE NORM OF THE MATRIX.
397C
398C EPS   THE RELATIVE MACHINE PRECISION.
399C
400C UTOL   THE USER TOLERANCE.
401C
402C TARR  AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN CALLS TO
403C       DLAEIG
404C
405C DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE
406C       CONSISTENCY IN CALLS TO DCOPY
407C
408      DZERO(1) = 0.0D0
409      RNORM = 0.0D0
410      IF (NPERM.NE.0) RNORM = DMAX1(-VAL(1),VAL(NPERM))
411      PNORM = RNORM
412      DELTA = 10.D30
413      EPSRT = DSQRT(EPS)
414      NLEFT = NVAL - NPERM
415      NOP = 0
416      NUMBER = NPERM
417      RARITZ = .FALSE.
418      UTOL = DMAX1(DBLE(FLOAT(N))*EPS,10.0D0**DBLE((-FLOAT(NFIG))))
419      J = MAXJ
420C
421C ------------------------------------------------------------------
422C
423C ANY ITERATION OF THE ALGORITHM BEGINS HERE.
424C
425   30 DO 50 I=1,NBLOCK
426         TEMP = DNRM2(N,P1(1,I),1)
427         IF (TEMP.EQ.0D0) CALL DLARAN(N, P1(1,I))
428   50 CONTINUE
429      IF (NPERM.EQ.0) GO TO 70
430      DO 60 I=1,NPERM
431         TAU(I) = 1.0D0
432         OTAU(I) = 0.0D0
433   60 CONTINUE
434   70 CALL DCOPY(N*NBLOCK, DZERO, 0, P0, 1)
435      CALL DCOPY(NBLOCK*NBLOCK, DZERO, 0, BET, 1)
436      CALL DCOPY(J*NBAND, DZERO, 0, T, 1)
437      MTEMP = NVAL + 1
438      DO 75 I = 1, MTEMP
439         CALL DCOPY(J, DZERO, 0, S(1,I), 1)
440   75 CONTINUE
441      NGOOD = 0
442      TMIN = 1.0D30
443      TMAX = -1.0D30
444      TEST = .TRUE.
445      ENOUGH = .FALSE.
446      BETMAX = 0.0D0
447      J = 0
448C
449C ------------------------------------------------------------------
450C
451C THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP.
452C
453   80 J = J + NBLOCK
454C
455C THIS IS THE SELECTIVE ORTHOGONALIZATION.
456C
457      IF (NUMBER.EQ.0) GO TO 110
458      DO 100 I=1,NUMBER
459         IF (TAU(I).LT.EPSRT) GO TO 100
460         TEST = .TRUE.
461         TAU(I) = 0.0D0
462         IF (OTAU(I).NE.0.0D0) OTAU(I) = 1.0D0
463         DO 90 K=1,NBLOCK
464            TEMP = -DDOT(N,VEC(1,I),1,P1(1,K),1)
465            CALL DAXPY(N, TEMP, VEC(1,I), 1, P1(1,K), 1)
466C
467C THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A
468C NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS
469C TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST.
470C
471            IF (DABS(TEMP*BET(K,K)).GT.DBLE(FLOAT(N))*EPSRT*
472     *        ANORM .AND. I.GT.NPERM) GO TO 380
473   90    CONTINUE
474  100 CONTINUE
475C
476C IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET.
477C
478  110 IF(.NOT. TEST) GO TO 160
479      CALL DORTQR(N, N, NBLOCK, P1, ALP)
480      TEST = .FALSE.
481      IF(J .EQ. NBLOCK) GO TO 160
482      DO 130 I = 1,NBLOCK
483         IF(ALP(I,I) .GT. 0.0D0) GO TO 130
484         M = J - 2*NBLOCK + I
485         L = NBLOCK + 1
486         DO 120 K = I,NBLOCK
487            BET(I,K) = -BET(I,K)
488            T(L,M) = -T(L,M)
489            L = L - 1
490            M = M + 1
491  120    CONTINUE
492  130 CONTINUE
493C
494C THIS IS THE LANCZOS STEP.
495C
496  160 CALL OP(N, NBLOCK, P1, P2)
497      NOP = NOP + 1
498      CALL IOVECT(N, NBLOCK, P1, J, 0)
499C
500C THIS COMPUTES P2=P2-P0*BET(TRANSPOSE)
501C
502      DO 180 I=1,NBLOCK
503         DO 170 K=I,NBLOCK
504            CALL DAXPY(N, -BET(I,K), P0(1,K), 1, P2(1,I), 1)
505  170    CONTINUE
506  180 CONTINUE
507C
508C THIS COMPUTES ALP AND P2=P2-P1*ALP.
509C
510      DO 200 I=1,NBLOCK
511         DO 190 K=1,I
512            II = I - K + 1
513            ALP(II,K) = DDOT(N,P1(1,I),1,P2(1,K),1)
514            CALL DAXPY(N, -ALP(II,K), P1(1,I), 1, P2(1,K), 1)
515            IF (K.NE.I) CALL DAXPY(N, -ALP(II,K), P1(1,K),
516     *        1, P2(1,I), 1)
517  190   CONTINUE
518  200 CONTINUE
519C
520C  REORTHOGONALIZATION OF THE SECOND BLOCK
521C
522      IF(J .NE. NBLOCK) GO TO 220
523      DO 215 I=1,NBLOCK
524         DO 210 K=1,I
525            TEMP = DDOT(N,P1(1,I),1,P2(1,K),1)
526            CALL DAXPY(N, -TEMP, P1(1,I), 1, P2(1,K), 1)
527            IF (K.NE.I) CALL DAXPY(N, -TEMP, P1(1,K),
528     *        1, P2(1,I), 1)
529            II = I - K + 1
530            ALP(II,K) = ALP(II,K) + TEMP
531  210   CONTINUE
532  215 CONTINUE
533C
534C THIS ORTHONORMALIZES THE NEXT BLOCK
535C
536  220 CALL DORTQR(N, N, NBLOCK, P2, BET)
537C
538C THIS STORES ALP AND BET IN T.
539C
540      DO 250 I=1,NBLOCK
541         M = J - NBLOCK + I
542         DO 230 K=I,NBLOCK
543            L = K - I + 1
544            T(L,M) = ALP(L,I)
545  230    CONTINUE
546         DO 240 K=1,I
547            L = NBLOCK - I + K + 1
548            T(L,M) = BET(K,I)
549  240    CONTINUE
550  250 CONTINUE
551C
552C THIS NEGATES T IF SMALL IS FALSE.
553C
554      IF (SMALL) GO TO 280
555      M = J - NBLOCK + 1
556      DO 270 I=M,J
557         DO 260 K=1,L
558            T(K,I) = -T(K,I)
559  260    CONTINUE
560  270 CONTINUE
561C
562C THIS SHIFTS THE LANCZOS VECTORS
563C
564  280 CALL DCOPY(NBLOCK*N, P1, 1, P0, 1)
565      CALL DCOPY(NBLOCK*N, P2, 1, P1, 1)
566      CALL DLAGER(J, NBAND, J-NBLOCK+1, T, TMIN, TMAX)
567      ANORM = DMAX1(RNORM, TMAX, -TMIN)
568      IF (NUMBER.EQ.0) GO TO 305
569C
570C THIS COMPUTES THE EXTREME EIGENVALUES OF ALP.
571C
572      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
573      CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK,
574     1   P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
575      ALPMIN = TARR(1)
576      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
577      CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR,
578     1  NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
579      ALPMAX = TARR(1)
580C
581C THIS COMPUTES ALP = BET(TRANSPOSE)*BET.
582C
583  305 DO 310 I = 1, NBLOCK
584         DO 300 K = 1, I
585            L = I - K + 1
586            ALP(L,K) = DDOT(NBLOCK-I+1, BET(I,I), NBLOCK, BET(K,I),
587     1        NBLOCK)
588  300    CONTINUE
589  310 CONTINUE
590      IF(NUMBER .EQ. 0) GO TO 330
591C
592C THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET.
593C
594      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
595      CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK,
596     1  P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, ANORM*ANORM)
597      BETMIN = DSQRT(TARR(1))
598C
599C THIS UPDATES TAU AND OTAU.
600C
601      DO 320 I=1,NUMBER
602         TEMP = (TAU(I)*DMAX1(ALPMAX-VAL(I),VAL(I)-ALPMIN)
603     *     +OTAU(I)*BETMAX+EPS*ANORM)/BETMIN
604         IF (I.LE.NPERM) TEMP = TEMP + RES(I)/BETMIN
605         OTAU(I) = TAU(I)
606         TAU(I) = TEMP
607  320 CONTINUE
608C
609C THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET.
610C
611  330 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
612      CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR,
613     1  NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0,
614     2  ANORM*ANORM)
615      BETMAX = DSQRT(TARR(1))
616      IF (J.LE.2*NBLOCK) GO TO 80
617C
618C ------------------------------------------------------------------
619C
620C THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND
621C LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK
622C IS JUSTIFIED.
623C
624      TOLG = EPSRT*ANORM
625      TOLA = UTOL*RNORM
626      IF(MAXJ-J .LT. NBLOCK .OR. (NOP .GE. MAXOP .AND.
627     1  NLEFT .NE. 0)) GO TO 390
628      GO TO 400
629
630C
631C ------------------------------------------------------------------
632C
633C THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO
634C SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN
635C ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY
636C IS AT 400.
637C
638  380 J = J - NBLOCK
639      IERR = -8
640  390 IF (NLEFT.EQ.0) RETURN
641      TEST = .TRUE.
642  400 NTHETA = MIN0(J/2,NLEFT+1)
643      CALL DLAEIG(J, NBAND, 1, NTHETA, T, VAL(NUMBER+1),
644     1  MAXJ, S, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
645      CALL DMVPC(NBLOCK, BET, MAXJ, J, S, NTHETA, ATEMP, VTEMP, D)
646C
647C THIS CHECKS FOR TERMINATION OF A CHECK RUN
648C
649      IF(NLEFT .NE. 0 .OR. J .LT. 6*NBLOCK) GO TO 410
650      IF(VAL(NUMBER+1)-ATEMP(1) .GT. VAL(NPERM) - TOLA) GO TO 790
651C
652C THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T
653C TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED.
654C
655 410  IF (NTHETA.LE.NLEFT) GO TO 470
656      IF (NPERM.EQ.0) GO TO 430
657      M = NUMBER + NLEFT + 1
658      IF (VAL(M).GE.VAL(NPERM)) GO TO 430
659      NPERM = NPERM - 1
660      NGOOD = 0
661      NUMBER = NPERM
662      NLEFT = NLEFT + 1
663      GO TO 400
664C
665C THIS UPDATES DELTA.
666C
667  430 M = NUMBER + NLEFT + 1
668      DELTA = DMIN1(DELTA,VAL(M))
669      ENOUGH = .TRUE.
670      IF(NLEFT .EQ. 0) GO TO 80
671      NTHETA = NLEFT
672      VTEMP(NTHETA+1) = 1
673C
674C ------------------------------------------------------------------
675C
676C THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL.
677C
678C THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES.
679C
680      IF (.NOT.(TEST .OR. ENOUGH)) GO TO 470
681      DELTA = DMIN1(DELTA,ANORM)
682      PNORM = DMAX1(RNORM,DMAX1(-VAL(NUMBER+1),DELTA))
683      TOLA = UTOL*PNORM
684      NSTART = 0
685      DO 460 I=1,NTHETA
686         M = NUMBER + I
687         IF (DMIN1(ATEMP(I)*ATEMP(I)/(DELTA-VAL(M)),ATEMP(I))
688     *     .GT.TOLA) GO TO 450
689         IND(I) = -1
690         GO TO 460
691C
692  450    ENOUGH = .FALSE.
693         IF (.NOT.TEST) GO TO 470
694         IND(I) = 1
695         NSTART = NSTART + 1
696  460 CONTINUE
697C
698C  COPY VALUES OF IND INTO VTEMP
699C
700      DO 465 I = 1,NTHETA
701         VTEMP(I) = DBLE(FLOAT(IND(I)))
702  465 CONTINUE
703      GO TO 500
704C
705C THIS CHECKS FOR NEW GOOD VECTORS.
706C
707  470 NG = 0
708      DO 490 I=1,NTHETA
709         IF (VTEMP(I).GT.TOLG) GO TO 480
710         NG = NG + 1
711         VTEMP(I) = -1
712         GO TO 490
713C
714  480    VTEMP(I) = 1
715  490 CONTINUE
716C
717      IF (NG.LE.NGOOD) GO TO 80
718      NSTART = NTHETA - NG
719C
720C ------------------------------------------------------------------
721C
722C THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS.
723C IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED.
724C
725  500 TEST = TEST .AND. .NOT.ENOUGH
726      NGOOD = NTHETA - NSTART
727      NSTART = NSTART + 1
728      NTHETA = NTHETA + 1
729C
730C THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND
731C EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR
732C FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE
733C EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS
734C PAUSE.
735C
736      CALL DCOPY(NTHETA, VAL(NUMBER+1), 1, VAL(NPERM+1), 1)
737      IF (NSTART.EQ.0) GO TO 580
738      IF (NSTART.EQ.NTHETA) GO TO 530
739      CALL DVSORT(NTHETA, VTEMP, ATEMP, 1, VAL(NPERM+1), MAXJ,
740     *  J, S)
741C
742C THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING
743C VECTORS.
744C
745  530 IF (.NOT.TEST) NSTART = 0
746      IF (.NOT.TEST) GO TO 580
747C
748C  FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW
749C
750      TEMP = ATEMP(1)
751      DO 535 I = 1, NSTART
752         TEMP = DMIN1(TEMP, ATEMP(I))
753  535 CONTINUE
754      M = NGOOD + 1
755      L = NGOOD + MIN0(NSTART,NBLOCK)
756      DO 540 I=M,L
757         CALL DSCAL(J, TEMP/ATEMP(I), S(1,I), 1)
758  540 CONTINUE
759      M = (NSTART-1)/NBLOCK
760      IF (M.EQ.0) GO TO 570
761      L = NGOOD + NBLOCK
762      DO 560 I=1,M
763         DO 550 K=1,NBLOCK
764            L = L + 1
765            IF (L.GT.NTHETA) GO TO 570
766            I1 = NGOOD + K
767            CALL DAXPY(J, TEMP/ATEMP(L), S(1,L), 1, S(1,I1), 1)
768  550    CONTINUE
769  560 CONTINUE
770  570 NSTART = MIN0(NSTART,NBLOCK)
771C
772C THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS.
773C
774  580 IF (NGOOD.EQ.0 .OR. .NOT.(TEST .OR. ENOUGH)) GO TO 600
775      DO 590 I=1,NGOOD
776         M = NPERM + I
777         RES(M) = ATEMP(I)
778  590 CONTINUE
779C
780C THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE
781C LANCZOS VECTORS.
782C
783  600 NUMBER = NPERM + NGOOD
784      IF (TEST .OR. ENOUGH) CALL DCOPY(N*NBLOCK, DZERO, 0, P1, 1)
785      IF (NGOOD.EQ.0) GO TO 620
786      M = NPERM + 1
787      DO 610 I=M,NUMBER
788         CALL DCOPY(N, DZERO, 0, VEC(1,I), 1)
789  610 CONTINUE
790  620 DO 670 I=NBLOCK,J,NBLOCK
791         CALL IOVECT(N, NBLOCK, P2, I, 1)
792         DO 660 K=1,NBLOCK
793            M = I - NBLOCK + K
794            IF (NSTART.EQ.0) GO TO 640
795            DO 630 L=1,NSTART
796               I1 = NGOOD + L
797               CALL DAXPY(N, S(M,I1), P2(1,K), 1, P1(1,L), 1)
798  630       CONTINUE
799  640       IF (NGOOD.EQ.0) GO TO 660
800            DO 650 L=1,NGOOD
801               I1 = L + NPERM
802               CALL DAXPY(N, S(M,L), P2(1,K), 1, VEC(1,I1), 1)
803  650       CONTINUE
804  660    CONTINUE
805  670 CONTINUE
806      IF (TEST .OR. ENOUGH) GO TO 690
807C
808C THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE
809C TAU RECURRENCE.
810C
811      M = NPERM + 1
812      DO 680 I=M,NUMBER
813         TEMP = 1.0D0/DNRM2(N,VEC(1,I),1)
814         CALL DSCAL(N, TEMP, VEC(1,I), 1)
815         TAU(I) = 1.0D0
816         OTAU(I) = 1.0D0
817  680 CONTINUE
818C
819C  SHIFT S VECTORS TO ALIGN FOR LATER CALL TO DLAEIG
820C
821      CALL DCOPY(NTHETA, VAL(NPERM+1), 1, VTEMP, 1)
822      CALL DVSORT(NTHETA, VTEMP, ATEMP, 0, TARR, MAXJ, J, S)
823      GO TO 80
824C
825C ------------------------------------------------------------------
826C
827C THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE
828C PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING
829C THE PERMANENT VECTORS.
830C
831  690 IF (NGOOD.EQ.0 .AND. NOP.GE.MAXOP) GO TO 810
832      IF (NGOOD.EQ.0) GO TO 30
833C
834C THIS ORTHONORMALIZES THE VECTORS
835C
836      CALL DORTQR(NMVEC, N, NPERM+NGOOD, VEC, S)
837C
838C THIS SORTS THE VALUES AND VECTORS.
839C
840      IF(NPERM .NE. 0) CALL DVSORT(NPERM+NGOOD, VAL, RES, 0, TEMP,
841     *   NMVEC, N, VEC)
842      NPERM = NPERM + NGOOD
843      NLEFT = NLEFT - NGOOD
844      RNORM = DMAX1(-VAL(1),VAL(NPERM))
845C
846C THIS DECIDES WHERE TO GO NEXT.
847C
848      IF (NOP.GE.MAXOP .AND. NLEFT.NE.0) GO TO 810
849      IF (NLEFT.NE.0) GO TO 30
850      IF (VAL(NVAL)-VAL(1).LT.TOLA) GO TO 790
851C
852C THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED
853C TO LOOK FOR UNDISCLOSED MULTIPLICITIES.
854C
855      M = NPERM - NBLOCK + 1
856      IF (M.LE.0) RETURN
857      DO 780 I=1,M
858         L = I + NBLOCK - 1
859         IF (VAL(L)-VAL(I).LT.TOLA) GO TO 30
860  780 CONTINUE
861C
862C THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ
863C PROCEDURE IS NEEDED.
864C
865  790 M = NPERM - NBLOCK
866      IF (M.LE.0) RETURN
867      DO 800 I=1,M
868         L = I + NBLOCK
869         IF (VAL(L)-VAL(I).GE.TOLA) GO TO 800
870         RARITZ = .TRUE.
871         RETURN
872  800 CONTINUE
873C
874      RETURN
875C
876C ------------------------------------------------------------------
877C
878C THIS REPORTS THAT MAXOP WAS EXCEEDED.
879C
880  810 IERR = -2
881      GO TO 790
882C
883      END
884