SUBROUTINE BT(N,X,F,G,COMPFG,FM,EPS,MAXCOM,MAXIT,RESET, * IPRINT,IFAIL,IWORK,WORK,LWORK) C C B(undle)T(rust)-Algorithm (trying to minimize a convex function) C C Double Precision version C C N : dimension of the problem C X : vector, dimension N C Input : starting point x0 C Output: solution point C F : Input : f(x0) C Output: f(solution) C G : Input : subgradient g(x0) C COMPFG : SUBROUTINE, computes f(x) and a subgradient g(x), C Parameter: N,X,F,G C FM : estimate of the minimum of f C EPS : stopping criterion, C eucl. norm of alpha-subgradient =< EPS and alpha =< eps C MAXCOM : maximal number of calls to COMPFG C MAXIT : maximal number of iterations C RESET : maximal number of stored subgradients, RESET > 3 C IPRINT : to control printout C IPRINT=0: no printout C IPRINT=1: final printout C IPRINT=2: in iteration k C niter= number of the iteration C ncomp= number of calls to compfg C f = f(xk) C gn = norm of an alpha-subgradient at xk C ( gives information about optimality: C f(y) >= f(xk)-gn*norm(y-xk)-alpha for all y ) C IPRINT>2: debug printout C IFAIL : information about the result of BT C IFAIL=0: everything is okay C IFAIL=1: number of calls of compfg > MAXCOM C IFAIL=2: number of iterations > MAXIT C IFAIL=3: 'T too small' C IFAIL=4: not enough working space C IFAIL=5: failure in quadratic program C WORK : working array of dimension LWORK C LWORK : dimension of WORK, C LWORK >= (RESET+N+9)*RESET+4*N+10 C C*********************************************** C Version 1.0 of Bundle Trust, March 2, 1989 * C * C Copyright by: * C * C Helga Schramm * C Mathematisches Institut * C Universitaet Bayreuth * C Postfach 101251 * C D-8580 Bayreuth * C*********************************************** C INTEGER N,MAXIT,MAXCOM,IPRINT,IFAIL,LWORK,RESET DOUBLE PRECISION F,FM,EPS INTEGER IWORK(RESET) DOUBLE PRECISION X(N),G(N),WORK(LWORK) EXTERNAL COMPFG C C local variables C INTEGER LALPHA,LXN,LA,LHESS,LY,LYA,LD,LDL,LDU, * LWORK1,LWORK2,LWORKT C C LALPHA= 1 LXN = LALPHA+RESET LA = LXN+N LHESS = LA+RESET*N LY = LHESS+RESET*(RESET+1)/2 LYA = LY+RESET LD = LYA+RESET LDL = LD+N LDU = LDL+N LWORK1= LDU+N LWORK2= (RESET+11)*RESET/2+10 LWORKT= LWORK1+LWORK2-1 IFAIL = -1 IF (LWORK .LT. LWORKT) IFAIL= 4 IF (IPRINT .GT. 0) THEN WRITE(*,*) WRITE(*,*) 'BT-Algorithm' WRITE(*,*) '============' WRITE(*,*) WRITE(*,*) 'workspace provided is WORK(',LWORK,')' WRITE(*,*) 'to solve problem we need WORK(',LWORKT,')' WRITE(*,*) ENDIF IF (IFAIL .GT. 0) THEN IF (IPRINT .GT. 5) THEN WRITE(*,*) 'Not enough workspace to solve problem' WRITE(*,*) 'IFAIL= ', IFAIL ENDIF RETURN ENDIF CALL BT1(N,X,F,G,COMPFG,WORK(LALPHA),FM,EPS,MAXCOM,MAXIT, * RESET,IPRINT,IFAIL,WORK(LXN),WORK(LA),WORK(LHESS), * IWORK,WORK(LWORK1),LWORK2,WORK(LY),WORK(LYA), * WORK(LD),WORK(LDL),WORK(LDU)) END SUBROUTINE BT1(N,X,F,G,COMPFG,ALPHA,FM,EPS,MAXCOM,MAXIT,RESET, * IPRINT,IFAIL,XN,A,HESS,IWORK,WORK,LWORK,Y,YA, * D,DL,DU) INTEGER N,MAXIT,MAXCOM,IPRINT,IFAIL,LWORK,RESET DOUBLE PRECISION F,FM,EPS INTEGER IWORK(RESET) DOUBLE PRECISION X(N),G(N),ALPHA(RESET),XN(N),A(RESET,N), * HESS(RESET*(RESET+1)/2),WORK(LWORK), * Y(RESET),YA(RESET),D(N),DL(N),DU(N) EXTERNAL COMPFG,MLIS4 C DOUBLE PRECISION MU,ZERO PARAMETER (MU = 1.0D-1, * ZERO= 1.0D-10) C C local variables C INTEGER I,J,NITER,NITER1,NN,NCOMP,INFORM,LALT,K,L,MODE, * LOGIC,IS,LIN,KK,K1,LINF,NAC DOUBLE PRECISION T,FN,V,VDA,VSS,ADA,VGXAL,W,VK,GD,TL,TU,SU, * ALN,ALM,ZM,ZP,ALP,GG,TP,TMIN,FPN,AM,D2,TPS,TNC, * TOL,NU,VV,TA,TLL,TUL,ALL,ALU LOGICAL KRIT,LLIN INTRINSIC DABS,DSQRT c c Graph Partitioning additions: c temporary storage of eigenvalues real eigval(5), xbst c istop tells BT to abort if we get stuck in Lanczos common / stop / istop, mulply c xbst is current best feasible solution common / eigvls / eigval, xbst c Only write out the nperm reliable eigenvalues common / iconst / matmul, nfig, nblock, nval, nperm C C start C c GPP change: DL should be initialised else may be used before assigned DO 338 I=1,N DL(I)= 0.0d0 338 CONTINUE KRIT = .TRUE. LLIN = .FALSE. LINF = 0 NCOMP = 0 NITER = 0 NITER1= 0 IFAIL = -1 INFORM= 0 IS= 0 NN= 0 ALPHA(1)= 0.0D0 DO 1 I= 1,RESET YA(I)= -1.0D0 Y(I) = 1.0D0 1 CONTINUE SU = DMAX1(100.0D0*(F-FM),1.0D0) T = 0.0D0 TP = 0.0D0 TLL= 0.0D0 TUL= 0.0D0 GG = 0.0D0 DO 2 I=1,N GG= GG+G(I)*G(I) 2 CONTINUE ZP = DSQRT(GG) ALP= 0.0D0 ALM= 0.0D0 ALN= 0.0D0 ALL= 0.0D0 ALU= 0.0D0 IF (ZP .LE. EPS) THEN IFAIL= 0 IF (IPRINT .GT. 1) WRITE(*,*) 'x already optimal' RETURN ENDIF c gpp addition: print eigval and matmul IF (IPRINT .GT. 1) WRITE(*,'(3X,A,1X,A,7X,A,10X,A,9X,A)') * 'niter','ncomp','ubd','bfs','gap' xgap = f / xbst - 1.0 IF (IPRINT .GT. 1) * WRITE(*,'(A,I4,A,I4,A,f15.2,A,f15.2,A,f7.3,a,f12.4,a,f12.4)') * ' ',NITER1+1,' ',NCOMP+1,' ',F,' ',xbst,' ',xgap, * ' ', zp,' ',alp 998 format(1x, 5f13.5) C C iteration C 100 CONTINUE TA= T IF (KRIT) THEN IF (NITER .GE. RESET .OR. NN .EQ. RESET+100) THEN C C reset C IF (NN .EQ. RESET+100) THEN DO 3 I=1,N A(1,I)= -D(I)/T 3 CONTINUE ALPHA(1)= ALP DO 5 J=1,NITER IF (ABS(ALPHA(J)) .LT. ZERO) THEN ALPHA(2)= DMAX1(ALPHA(J),0.0D0) DO 4 I=1,N A(2,I)= A(J,I) 4 CONTINUE ENDIF 5 CONTINUE NN= 2 NITER= 2 GOTO 111 ENDIF NN= 0 DO 7 J=1,NITER IF ((Y(J).GE.ZERO) .OR. (ABS(ALPHA(J)) .LT. ZERO)) THEN NN= NN+1 ALPHA(NN)= DMAX1(ALPHA(J),0.0D0) DO 6 I=1,N A(NN,I)= A(J,I) 6 CONTINUE ENDIF 7 CONTINUE IF (NN .GE. RESET) THEN NN= 2 ALPHA(1)= ALP DO 8 I=1,N A(1,I)= -D(I)/T 8 CONTINUE DO 10 J=1,NITER IF (ABS(ALPHA(J)) .LT. ZERO) THEN ALPHA(2)= DMAX1(ALPHA(J),0.0D0) DO 9 I=1,N A(2,I)= A(J,I) 9 CONTINUE ENDIF 10 CONTINUE ENDIF 111 NITER= NN DO 11 I=1,RESET Y(I)= 0.0D0 YA(I)= -1.0D0 11 CONTINUE DO 12 I=NITER+1,RESET ALPHA(I)= 0.0D0 12 CONTINUE C C hess update C DO 15 K=1,NITER LALT= K*(K-1)/2 DO 14 J=1,K K1= LALT+J HESS(K1)= 0.0D0 DO 13 I=1,N HESS(K1)= HESS(K1)+A(J,I)*A(K,I) 13 CONTINUE 14 CONTINUE 15 CONTINUE KK= 1 Y(1)= 1.0D0 MODE= 1 IF (IPRINT .GT. 1) THEN WRITE(*,*) WRITE(*,*) 'Reset, nn=',NN WRITE(*,*) ENDIF ENDIF C C end of reset C NITER1= NITER1+1 IF (IFAIL .LT. 0 .AND. NITER1 .GE. MAXIT) THEN IFAIL= 2 IF (IPRINT .GT. 1) WRITE(*,*) 'max. number of iterations' MAXCOM= NCOMP ENDIF C C data for quad C NITER= NITER+1 ALPHA(NITER)= ALN DO 16 I=1,N A(NITER,I)= G(I) 16 CONTINUE LALT= NITER*(NITER-1)/2 DO 18 J=1,NITER K =LALT+J HESS(K)= 0.0D0 DO 17 I=1,N HESS(K)= HESS(K)+A(J,I)*A(NITER,I) 17 CONTINUE 18 CONTINUE IF (NITER .EQ. 1) Y(1)= 1.0D0 TL = 0.0D0 TU = SU ZM = ZP ALM= ALP ENDIF C C end krit C IF (NITER .EQ. 1 .AND. KRIT) THEN C C compute first t via line search C DO 19 I=1,N D(I) = -G(I) XN(I)= X(I) 19 CONTINUE FN= F T= 2.0D0*(F-FM)/ZP T= DMAX1(T,1.0D0/ZP) TMIN= 1.0D-20 FPN= -GG D2= GG AM= 0.0D0 TOL= DMAX1((F-FM),EPS) CALL MLIS4 (COMPFG,N,XN,FN,FPN,T,TMIN,SU,D,D2,G,0.2D0,0.1D0, * LOGIC,NCOMP,MAXCOM,X,TOL,AM,TPS,TNC) IF (IPRINT. GE. 10) WRITE(*,'(5X,A,I4,3X,E16.8)') * 'line search:',LOGIC,T V= -T*GG IF (IPRINT. GE. 5) THEN WRITE(*,'(48X,A,E16.8)') 't =',T WRITE(*,'(48X,A,E16.8)') 'rho=',0.5*T*T*GG ENDIF IF (LOGIC .EQ. 3) THEN C C null step C ALN= DMAX1(TPS,0.0D0) DO 20 I=1,N X(I)= XN(I) 20 CONTINUE KRIT= .TRUE. GOTO 100 ENDIF IF (LOGIC .GE. 0 .AND. LOGIC. LE. 2) THEN C C serious step C DO 22 I=1,NITER GD= 0.0D0 DO 21 J=1,N GD= GD+T*D(J)*A(I,J) 21 CONTINUE ALPHA(I)= ALPHA(I)+FN-F-GD ALPHA(I)= DMAX1(ALPHA(I),0.0D0) 22 CONTINUE ALN= 0.0D0 DO 23 I=1,N X(I)= XN(I) 23 CONTINUE F= FN KRIT= .TRUE. GOTO 100 ENDIF T= 0.1D0*(F-FM)/ZP T= DMAX1(T,1.0D0/ZP) DO 210 I=1,N X(I)= XN(I) 210 CONTINUE ENDIF C C compute d(t) C MODE = 2 LLIN= .FALSE. IF (KRIT .OR. NITER .LE. 2) THEN MODE = 1 LIN = 0 ENDIF IF (MODE .EQ. 2) THEN Y(NITER) = 1.0D0 IWORK(NITER)= 1 ENDIF IF (LIN .GT. 2 .AND. T .GT. TLL .AND. T .LT. TUL) THEN C C piecewise linear C LLIN= .TRUE. NU= (T-TLL)/(TUL-TLL) DO 24 I=1,N D(I)= DL(I)+NU*(DU(I)-DL(I)) 24 CONTINUE ALP= ALL+NU*(ALU-ALL) V= -ALPHA(1)+A(1,1)*D(1) DO 25 I=2,N V= V+A(1,I)*D(I) 25 CONTINUE DO 27 J= 2,NITER VV= -ALPHA(J)+A(J,1)*D(1) DO 26 I=2,N VV= VV+A(J,I)*D(I) 26 CONTINUE V= DMAX1(V,VV) 27 CONTINUE ELSE C C quadratic program C CALL QUAD01(NITER,N,HESS,T,ALPHA,MODE,KK,IWORK,Y,ADA,VK, * VDA,VSS,W,L,VGXAL,WORK,LWORK,INFORM) C V= T*VDA C DO 29 I=1,N D(I)= 0.0D0 DO 28 J=1,NITER D(I)= D(I)-T*Y(J)*A(J,I) 28 CONTINUE 29 CONTINUE C IF (IPRINT .GT. 10) THEN NAC= 0 WRITE(*,*) WRITE(*,'(/A,I2,A,E23.15)') 'inform',INFORM,' w',W WRITE(*,*) ' active subgradients:' DO 30 I=1,NITER IF (Y(I) .GE. ZERO) THEN WRITE(*,*) I,Y(I),ALPHA(I) NAC= NAC+1 ENDIF 30 CONTINUE WRITE(*,*) 'NAC=',NAC ENDIF ENDIF C LIN= LIN+1 DO 31 I=1,NITER IF (YA(I)*Y(I) .LT. ZERO .AND. (DABS(YA(I))+Y(I)) .GT. ZERO) * LIN= 0 YA(I)= Y(I) 31 CONTINUE C C next trialpoint C DO 32 I=1,N XN(I)= X(I)+D(I) 32 CONTINUE C C compute zp,alpha C ZP= 0.0D0 DO 33 I=1,N ZP= ZP+D(I)*D(I) 33 CONTINUE ZP= ZP/(T**2) IF (.NOT. LLIN) ALP= ADA C IF (IPRINT .GT. 20) WRITE(*,*) 'zp2,alp :',ZP,ALP IF (ZP .LE. -ZERO .OR. ALP .LE. -ZERO * .OR. INFORM .GT. 0 .AND. LINF .LE. 3) THEN LINF= LINF+1 KRIT= .TRUE. NN = RESET+100 IF (IPRINT .GT. 5) WRITE(*,*) 'linf= ',LINF IF (LINF .LE. 3) GOTO 100 ENDIF IF (ZP .LE. -ZERO .OR. ALP .LE. -ZERO * .OR. INFORM .GT. 0 .OR. LINF .GT. 3) THEN IF (IPRINT .GT. 1) WRITE(*,*) 'failure in qpdf2s' IF (IPRINT .GT. 1) WRITE(*,*) 'inform= ',INFORM IFAIL = 5 MAXIT = NITER1 MAXCOM= NCOMP ENDIF C ZP= DSQRT(DABS(ZP)) C c GPP addition: write out the eigenvalues xgap = f / xbst - 1.0 IF ((IPRINT .GT. 1 .AND. KRIT) .OR. IPRINT .GE. 5) then WRITE(*,'(A,I4,A,I4,A,f15.2,A,f15.2,A,f7.3,a,f12.4,a,f12.4)') * ' ',NITER1,' ',NCOMP+1,' ',F,' ',xbst,' ',xgap, * ' ', zp,' ',alp endif IF (IPRINT. GE. 5) THEN WRITE(*,'(48X,A,E16.8)') 't =',T WRITE(*,'(48X,A,E16.8)') 'rho=',0.5D0*(T*ZP)**2 ENDIF C C compute f(xn),g(xn) C CALL COMPFG(N,XN,FN,G) NCOMP= NCOMP+1 C IF (IFAIL .LT. 0 .AND. NCOMP .GE. MAXCOM) THEN IFAIL= 1 IF (IPRINT .GT. 1) WRITE(*,*) 'max. number of calls to compfg' MAXIT= NITER1 ENDIF c GPP change: abort if we have become stuck in Lanczos if (istop .eq. 1) then ifail = 1 maxit = niter1 IF (IPRINT .GT. 1) * WRITE(*,*) ' Lanczos does not converge under given setup' return endif C C stopping criterion C IF ((DABS(ZP) .LE. EPS) .AND. (ALP .LE. EPS)) THEN IF (IPRINT .GE. 1) THEN WRITE(*,*) WRITE(*,*) 'convergence' WRITE(*,*) ENDIF IFAIL = 0 MAXCOM= NCOMP MAXIT = NITER1 ELSE C C compute t C KRIT= .FALSE. C C test serious criterion C IF (FN-F .LT. 0.1D0*V) THEN TL= T GD= 0.0D0 DO 34 I=1,N GD= GD+G(I)*D(I) 34 CONTINUE IF ((GD .GE. 0.2D0*V) .OR. (TU-T .LE. MU)) THEN C C serious iteration C KRIT= .TRUE. IS= MAX(1,IS+1) IF ((FN-F .LE. 0.7D0*V) .OR. (IS .GT. 2)) THEN IF (ABS(V-FN+F).GE.ZERO) * TP= -0.5D0*V*T/ABS(V-FN+F) T= DMAX1(TP,5.0D0*T) T= DMIN1(T,SU) IS= 0 ENDIF C C update alpha C DO 36 I=1,NITER GD= 0.0D0 DO 35 J=1,N GD= GD+A(I,J)*D(J) 35 CONTINUE ALPHA(I)= ALPHA(I)+FN-F-GD ALPHA(I)= DMAX1(ALPHA(I),0.0D0) 36 CONTINUE ALN= 0.0D0 DO 37 I=1,N X(I)= XN(I) 37 CONTINUE F= FN ELSE IF (TU .EQ. SU) THEN IF (ABS(V-FN+F).GE.ZERO) TP= -0.5D0*V*T/ABS(V-FN+F) T= DMAX1(TP,5.0D0*T) T= DMIN1(T,SU) ELSE IF (ABS(V-FN+F).GE.ZERO) TP= -0.5D0*V*T/ABS(V-FN+F) T= DMAX1(TP,TL+0.25D0*(TU-TL)) T= DMIN1(T,TL+0.5D0*(TU-TL)) ENDIF IF (LIN .EQ. 1 .OR. LIN .EQ. 2) THEN DO 38 I=1,N DL(I)= D(I) 38 CONTINUE TLL= TL ALL= ALP ENDIF ENDIF ELSE TU= T IF (TL .EQ. 0.0D0) THEN C C test null criterion C ALN= F-FN DO 39 I=1,N ALN= ALN+G(I)*D(I) 39 CONTINUE C IF ((ALN. LE. DMAX1(ALP,EPS)) .OR. * (DABS(F-FN). LE. 1.0D0*(ZM+ALM))) THEN C C null iteration C KRIT= .TRUE. ALN= DMAX1(ALN,0.0D0) IS= MIN(-1,IS-1) IF (IS .LT. -N) THEN T= 0.5D0*T IS= 0 ENDIF ELSE T= TL+0.1D0*(TU-TL) ENDIF ELSE T= TL+0.5D0*(TU-TL) IF (LIN .EQ. 1 .OR. LIN .EQ. 2) THEN DO 40 I=1,N DU(I)= D(I) 40 CONTINUE TUL= TU ALU= ALP ENDIF ENDIF ENDIF ENDIF IF (IFAIL .LT. 0 .AND. DABS(T-TA) .LT. ZERO*1.0D-3 * .AND. .NOT. KRIT) THEN IFAIL = 3 IF (IPRINT .GT. 1) WRITE(*,*) 'T is wrong' MAXIT = NITER1 MAXCOM= NCOMP ENDIF IF (IFAIL .LT. 0) GOTO 100 END SUBROUTINE QUAD01(M,N,G,T,ALPHA,MODE,K,JS,Y,ADA,V,VDA, * VSS,W,L,VGXAL,WORK,LWORK,INFORM) INTEGER M,N,MODE,K,L,LWORK,INFORM DOUBLE PRECISION T,ADA,V,VDA,VSS,W,VGXAL INTEGER JS(M) DOUBLE PRECISION G(M*(M+1)/2),ALPHA(M),Y(M),WORK(LWORK) C C local variables C INTEGER ITMAX,IDELMX,ITER DOUBLE PRECISION EPSTOP,EPSCON PARAMETER ( ITMAX = 500, * EPSTOP = 1.0D-12, * EPSCON = 1.0D-12, * IDELMX = 500 ) C CALL QUAD02(ITMAX,M,N,M,G,M*(M+1)/2,1.0D0/T,ALPHA,EPSTOP, * EPSCON,IDELMX,MODE,K,JS,Y,ADA,V,VDA,VSS, * W,L,VGXAL,WORK,LWORK,ITER,INFORM) END SUBROUTINE QUAD02(ITMAX,M,N,KMAX,G,LENG,U,A,EPSTOP,EPSCON, * IDELMX,MODE,K,JS,X,ADD,V,VDD,VSS,W, * L,VGXAL,WORK,LWORK,ITER,INFORM) C C INTEGER ITMAX,M,N,KMAX,LENG,IDELMX,MODE,L,LWORK, * ITER,INFORM DOUBLE PRECISION U,EPSTOP,EPSCON,ADD,V,VDD,VSS,W,VGXAL INTEGER JS(KMAX) DOUBLE PRECISION G(LENG),A(M),X(M),WORK(LWORK) C INTEGER LWORK1,K,KMAX1,K1,IDELET,KMAXA DOUBLE PRECISION WOLD C C IF (MODE .LE. 1) THEN LWORK1=10+KMAX*(KMAX+11)/2 IF (M .LT. 1 .OR. N .LT. 1 .OR. KMAX .LT. MIN0(M,N+1).OR. * LENG .LT. M*(M+1)/2 .OR. * EPSTOP .LT. 0.0D0 .OR. EPSCON .LE. 0.0D0 .OR. * MODE .LT. 1 .OR. LWORK .LT. LWORK1) THEN INFORM = 5 RETURN ENDIF KMAX1= KMAX WORK(1)= LWORK1 WORK(2)= KMAX1 ELSE LWORK1= WORK(1)+0.1D0 IF (LWORK .LT. LWORK1) THEN INFORM= 5 RETURN ENDIF KMAX1= WORK(2)+0.1D0 K1 = WORK(3)+0.1D0 IF (KMAX1 .LT. MIN0(M,N+1) .OR. EPSTOP .LT. 0.0D0 .OR. * EPSCON .LE. 0.0D0 .OR. K .NE. K1) THEN INFORM= 5 RETURN ENDIF WOLD = WORK(4) IDELET= WORK(5)+0.1D0 ENDIF KMAXA = MIN0(M,N+1) CALL QUAD03(ITMAX,M,N,KMAXA,G,LENG,U,A,EPSTOP,EPSCON,IDELMX,MODE, * K,JS,X,ADD,V,VDD,VSS,W,L,VGXAL,ITER,INFORM, * WORK(10),WORK(KMAX1+10),WORK(2*KMAX1+10), * WORK(3*KMAX1+10),WORK(4*KMAX1+10),WORK(5*KMAX1+10), * KMAXA*(KMAXA+1)/2,WOLD,IDELET) WORK(3)= K WORK(4)= WOLD WORK(5)= IDELET END SUBROUTINE QUAD03(ITMAX,M,N,KMAX,G,LENG,U,A,EPSTOP,EPSCON,IDELMX, * MODE,K,JS,X,ADD,V,VDD,VSS,W,L,VGXAL,ITER, * INFORM,Y,YHAT,YSS,YA,YE,R,MDIMR,WOLD,IDELET) C INTEGER ITMAX,M,N,KMAX,LENG,IDELMX,MODE,K,L,ITER, * INFORM,MDIMR,IDELET DOUBLE PRECISION U,EPSTOP,EPSCON,ADD,V,VDD,VSS,W,VGXAL, * WOLD INTEGER JS(KMAX) DOUBLE PRECISION G(LENG),A(M),X(M),Y(KMAX),YHAT(KMAX), * YSS(KMAX),YA(KMAX),YE(KMAX),R(MDIMR) C INTEGER NINCRW,INCRW,INCRWM,KDIMR,ISETYA,ISETYE,NAUG, * NEXC,NDEL,NREF,IWOLD,JDEL,JADD,I,J,JA,LG,KG, * IA,JR,JG,IG,KR,JROT,JROTAT,ISTEP,KMIN1,KA,J1, * KREFAC DOUBLE PRECISION T,ERRKT,ERRX,EPS1,EPS2,HALF,ONE,TWO,S,DELTA, * DELTAP,DELTAW,DMU,ERRORV,R1,R2,RHOSQ,VIOLAT, * WMIN,YBARL,YANEW,YENEW,DTAU,DGAMMA,DSIGMA,DNU, * DALPHA,S1,YAYE,YESQR,DABSS C COMMON /QPLOG/ ERRKT,ERRX,NAUG,NEXC,NDEL DATA EPS1, EPS2, HALF, ONE, TWO * /0.5D0,1.0D-2,0.5D0,1.0D0,2.0D0/ C C NINCRW= 3 INCRW =-1 INCRWM=-1 IF (MODE .GT. 1) KDIMR= K*(K+1)/2 ISETYA= 0 IF (MODE .EQ. 2) ISETYA=1 ISETYE= 0 NAUG = 0 NEXC = 0 NDEL = 0 NREF = 0 ITER =-1 IWOLD = 0 JDEL = 0 JADD = 0 T = 0.0D0 GOTO (10,500,100,151,500,500,500,513),MODE 10 CONTINUE 50 W= HALF*G(1)+U*A(1) L= 1 KG=0 DO 60 J=1,M KG= KG+J S= HALF*G(KG)+U*A(J) IF (W .GT. S) THEN W=S L=J ENDIF 60 CONTINUE 70 CONTINUE K = 1 KDIMR = 1 JS(1) = L YHAT(1)= ONE LG = L*(L+1)/2 V =-(G(LG)+U*A(L)) IDELET = 0 R(1) = DSQRT(ONE+G(LG)) YA(1) = -A(L)/R(1) YE(1) = ONE/R(1) JADD = L 100 CONTINUE DO 110 J=1,M X(J)=0.0D0 110 CONTINUE ERRX= 0.0D0 ADD= 0.0D0 DO 120 J=1,K JA = JS(J) X(JA)= YHAT(J) ERRX = ERRX+YHAT(J) ADD = ADD+A(JA)*YHAT(J) 120 CONTINUE ERRX = ONE -ERRX VDD = 0.0D0 VGXAL= 0.0D0 ERRKT= 0.0D0 DO 150 I=1,M S=0.0D0 DO 140 J=1,K JA= JS(J) IG= MIN0(I,JA) JG= MAX0(I,JA) KG= (JG-1)*JG/2+IG S= S+G(KG)*YHAT(J) 140 CONTINUE C S=SUM OF G(I,J)*X(J) FOR J=1 THROUGH M. VDD= VDD+X(I)*S S = S+U*A(I) IF (I .EQ. 1) VSS= -S VSS= DMAX1(VSS,-S) S= S+V DABSS= DABS(S) IF (X(I) .GT. 0.0D0) ERRKT=DMAX1(ERRKT,DABSS) IF (S .GE. VGXAL .AND. I .GT. 1) GOTO 150 VGXAL=S L=I 150 CONTINUE W = HALF*VDD+U*ADD VDD=-(VDD+U*ADD) 151 CONTINUE LG= L*(L+1)/2 VIOLAT= VGXAL+EPSTOP*(ONE+G(LG)) ERRORV= DMAX1(DABS(V-VDD),DABS(V-VSS),DABS(VDD-VSS)) * /(ONE+DABS(V)) IF (IWOLD .EQ. 0) WOLD= W IWOLD = 1 DELTAW= WOLD-W WOLD = W IF (DELTAW .LE. 0.0D0) INCRW= INCRW+1 IF(VIOLAT.GE.0.0D0) GOTO 800 DO 160 J=1,K IF (L .EQ. JS(J)) GOTO 8160 160 CONTINUE IF (ITER .GE. ITMAX) GOTO 8161 IF (INCRW .EQ. NINCRW) GOTO 8162 IF(INCRWM .EQ. -1) WMIN= W IF(W .EQ. WMIN) INCRWM= INCRWM+1 IF(W .GE. WMIN) GOTO 170 WMIN = W INCRWM= 0 170 CONTINUE IF (INCRWM .EQ. NINCRW) GOTO 8170 DO 212 J=1,K JA = JS(J) IG = MIN0(L,JA) JG = MAX0(L,JA) KG = (JG-1)*JG/2+IG X(J)= G(KG) YSS(J)= ONE+G(KG) 212 CONTINUE CALL QUAD04(R,YSS,KDIMR,K,1) DO 215 J=1,K Y(J)= YSS(J) 215 CONTINUE RHOSQ= 0.0D0 DO 220 J=1,K RHOSQ=RHOSQ+YSS(J)**2 220 CONTINUE RHOSQ= (ONE+G(LG))-RHOSQ IF (RHOSQ .GT. EPSCON*(ONE+G(LG)) .AND. K .LT. KMAX) GOTO 270 CALL QUAD04(R,YSS,KDIMR,K,0) YBARL= 0.0D0 DO 230 J=1,K YBARL= YBARL+YSS(J) 230 CONTINUE DELTA = YBARL-ONE DELTAP= 0.0D0 DO 235 I=1,K IA= JS(I) S = 0.0D0 DO 232 J=1,K JA= JS(J) IG= MIN0(IA,JA) JG= MAX0(IA,JA) KG= (JG-1)*JG/2+IG S=S+G(KG)*YSS(J) 232 CONTINUE 235 DELTAP= DELTAP+YSS(I)*S S= 0.0D0 DO 240 J=1,K S=S+YSS(J)*X(J) 240 CONTINUE DELTAW= (DELTAP+G(LG)*YBARL**2)-TWO*YBARL*S DELTAP= DABS((DELTAP+G(LG))-TWO*S) S = RHOSQ-(DELTA**2+DELTAP) IF (DELTA .LT. -EPS1 .AND. K .LT. KMAX) GOTO 260 T = 100./(ONE-EPS1) JR= 0 DO 250 J=1,K IF (YSS(J) .LE. 0.0D0) GOTO 250 IF (T*YSS(J) .LT. YHAT(J)) GOTO 250 T = YHAT(J)/YSS(J) JR= J 250 CONTINUE IF (HALF*T*DELTAW .LT. (EPS2-ONE-DELTA)*VGXAL) GOTO 400 IF (K .EQ. KMAX .AND. JR .NE. 0) GOTO 400 IF (K .EQ. KMAX) GOTO 8250 260 CONTINUE RHOSQ= DMAX1(RHOSQ,DELTA**2+DELTAP) IF(RHOSQ .EQ. 0.0D0) GOTO 8250 270 CONTINUE NAUG = NAUG+1 ISTEP= 3 YANEW=-A(L) YENEW= ONE KDIMR= KDIMR+1 DO 310 J=1,K R(KDIMR)= Y(J) YANEW= YANEW-Y(J)*YA(J) YENEW= YENEW-Y(J)*YE(J) KDIMR= KDIMR+1 310 CONTINUE R(KDIMR)= DSQRT(RHOSQ) K = K+1 JS(K) = L YHAT(K)= 0.0D0 YA(K) = YANEW/R(KDIMR) YE(K) = YENEW/R(KDIMR) JADD = L T = 0.0D0 GOTO 500 400 CONTINUE NEXC = NEXC+1 DO 410 J=1,K YHAT(J)= DMAX1(YHAT(J)-T*YSS(J),0.0D0) 410 CONTINUE ISTEP= 4 415 CONTINUE JDEL= JS(JR) IF (JR .EQ. K) GOTO 446 IDELET= IDELET+1 KMIN1 = K-1 DO 420 J=JR,KMIN1 YHAT(J)= YHAT(J+1) JS(J)=JS(J+1) 420 CONTINUE IF (ISTEP .NE. 4) GOTO 422 DO 421 J=JR,KMIN1 X(J)=X(J+1) 421 CONTINUE 422 CONTINUE JROTAT= JR+1 KR= JR*(JR+1)/2-1 DO 431 JROT=JROTAT,K KR= KR+JROT R1= R(KR) R2= R(KR+1) DMU= DMAX1(DABS(R1),DABS(R2)) DTAU= DMU*DSQRT((R1/DMU)**2+(R2/DMU)**2) R(KR)= DTAU DGAMMA= R1/DTAU DSIGMA= R2/DTAU DNU = DSIGMA/(ONE+DGAMMA) DALPHA= DGAMMA*YA(JROT-1)+DSIGMA*YA(JROT) YA(JROT) = DNU*(YA(JROT-1)+DALPHA)-YA(JROT) YA(JROT-1)= DALPHA DALPHA= DGAMMA*YE(JROT-1)+DSIGMA*YE(JROT) YE(JROT)= DNU*(YE(JROT-1)+DALPHA)-YE(JROT) YE(JROT-1)= DALPHA IF (JROT .EQ. K) GOTO 435 KA= KR DO 430 J=JROT,KMIN1 KA= KA+J DALPHA = DGAMMA*R(KA)+DSIGMA*R(KA+1) R(KA+1)= DNU*(R(KA)+DALPHA)-R(KA+1) R(KA) = DALPHA 430 CONTINUE 431 CONTINUE 435 KR= (JR-1)*JR/2+1 KA= KR+JR DO 445 J=JR,KMIN1 DO 440 J1=1,J R(KR)= R(KA) KR= KR+1 KA= KA+1 440 CONTINUE KA=KA+1 445 CONTINUE 446 CONTINUE 6 KDIMR= KDIMR-K K = K-1 IF (DABS(A(JDEL)) .GT. ONE .AND. JR .LE. K) ISETYA= 1 IF(ISTEP .EQ. 6) GOTO 500 450 CONTINUE YANEW=-A(L) YENEW= ONE IF (K .GT. 0) GOTO 460 KDIMR = 1 IDELET= 0 S1= 0.0D0 GOTO 480 460 CONTINUE IF(ISTEP.NE.4) GOTO 462 DO 461 J=1,K YSS(J)=ONE+X(J) 461 CONTINUE GOTO 472 462 CONTINUE DO 471 J=1,K JA= JS(J) IG= MIN0(L,JA) JG= MAX0(L,JA) KG= (JG-1)*JG/2+IG YSS(J)=ONE+G(KG) 471 CONTINUE 472 CONTINUE CALL QUAD04(R,YSS,KDIMR,K,1) KDIMR= KDIMR+1 S1 = 0.0D0 DO 473 J=1,K S1= S1+YSS(J)**2 R(KDIMR)= YSS(J) YANEW= YANEW-YSS(J)*YA(J) YENEW= YENEW-YSS(J)*YE(J) KDIMR= KDIMR+1 473 CONTINUE 480 R(KDIMR)= DSQRT(DMAX1((ONE+G(LG))-S1,0.0D0)) K= K+1 JS(K)= L IF (ISTEP .EQ. 4) YHAT(K)= T*(ONE+DELTA) IF (R(KDIMR) .EQ. 0.0D0) GOTO 8480 YA(K)= YANEW/R(KDIMR) YE(K)= YENEW/R(KDIMR) IF(ISTEP .EQ. 7) GOTO 750 JADD= L 500 CONTINUE 505 ITER= ITER+1 IF (ISETYA .EQ. 0 .AND. ISETYE .EQ. 0) GOTO 513 DO 510 J=1,K JA=JS(J) IF(ISETYA.EQ.1) YA(J)=-A(JA) IF(ISETYE.EQ.1) YE(J)=ONE 510 CONTINUE IF(ISETYA .EQ. 1) CALL QUAD04(R,YA,KDIMR,K,1) IF(ISETYE .EQ. 1) CALL QUAD04(R,YE,KDIMR,K,1) ISETYA= 0 ISETYE= 0 513 CONTINUE YESQR= 0.0D0 YAYE = 0.0D0 DO 515 J=1,K YESQR= YESQR+YE(J)**2 YAYE=YAYE+YE(J)*YA(J) 515 CONTINUE 516 CONTINUE V= (ONE-U*YAYE)/YESQR DO 517 J=1,K Y(J)=U*YA(J)+V*YE(J) 517 CONTINUE 540 CALL QUAD04(R,Y,KDIMR,K,0) V= ONE-V 550 CONTINUE JR= 0 DO 551 J=1,K IF(Y(J).LE.0.0D0) JR=J 551 CONTINUE JADD= 0 IF (JR .GT. 0) GOTO 600 DO 560 J=1,K YHAT(J)= Y(J) 560 CONTINUE JDEL= 0 T = ONE GOTO 100 600 CONTINUE NDEL= NDEL+1 T = ONE DO 610 J=1,K S1= YHAT(J)-Y(J) IF (S1 .LE. 0.0D0) GOTO 610 S1= YHAT(J)/S1 IF (T .LT. S1) GOTO 610 T = S1 JR= J 610 CONTINUE DO 620 J=1,K YHAT(J)=DMAX1(YHAT(J)+T*(Y(J)-YHAT(J)),0.0D0) 620 CONTINUE ISTEP= 6 GOTO 415 700 CONTINUE NREF = NREF+1 ISTEP = 7 IDELET= 0 KREFAC= K K = 0 750 CONTINUE IF (K .EQ. KREFAC) GOTO 760 L = JS(K+1) LG= L*(L+1)/2 GOTO 450 760 CONTINUE ISETYA= 1 ISETYE= 1 JDEL =-1000 JADD = 1000 T = 0.0D0 GOTO 500 800 INFORM= 0 805 DO 810 J=1,M X(J)=0.0D0 810 CONTINUE ERRX= 0.0D0 DO 820 J=1,K JA = JS(J) ERRX = ERRX+YHAT(J) X(JA)=YHAT(J) 820 CONTINUE ERRX= ONE-ERRX RETURN 8160 IF(IDELET .GT. IDELMX) GOTO 700 INFORM= 1 GOTO 805 8161 INFORM= 4 GOTO 805 8162 INFORM= 2 GOTO 805 8170 INFORM= 2 GOTO 805 8250 IF(IDELET .GT. IDELMX) GOTO 700 INFORM= 3 GOTO 805 8480 CONTINUE IF(IDELET .GT. 0) GOTO 700 NDEL= NDEL+1 IF (INCRW .EQ. NINCRW-1) INCRW=INCRW-1 JDEL =-JS(K) JADD = 0 T = 0.0D0 KDIMR= KDIMR-K K = K-1 S1 = 0.0D0 DO 8484 J=1,K S1=S1+YHAT(J) 8484 CONTINUE IF (S1 .EQ. 0.0D0) GOTO 8486 DO 8485 J=1,K YHAT(J)=YHAT(J)/S1 8485 CONTINUE ISETYA= 1 ISETYE= 1 GOTO 500 8486 YHAT(1)= ONE ISETYA= 1 ISETYE= 1 GOTO 500 END SUBROUTINE QUAD04(R,Y,MDIMR,M,ITRANS) INTEGER MDIMR,M,ITRANS DOUBLE PRECISION R(MDIMR),Y(M) C INTEGER I,J,K,KA,I1,IA,IOLD DOUBLE PRECISION S C C IF (ITRANS.EQ.1) GOTO 3 I= M K= MDIMR Y(I)= Y(I)/R(K) 1 IF (I .EQ. 1) RETURN I1= I K = K-I KA= K I = I-1 J = I S = Y(I) DO 2 IA=I1,M KA= KA+J S = S-R(KA)*Y(IA) J = J+1 2 CONTINUE Y(I)= S/R(K) GOTO 1 3 Y(1)= Y(1)/R(1) I= 1 K= 1 4 IF (I .EQ. M) RETURN IOLD= I I= I+1 K= K+1 S= Y(I) DO 5 IA=1,IOLD S= S-R(K)*Y(IA) K=K+1 5 CONTINUE Y(I)= S/R(K) GOTO 4 END SUBROUTINE MLIS4 (SIMUL,N,XN,FN,FPN,T,TMIN,TMAX,D,D2,G, * AMD,AMF,LOGIC,NAP,NAPMAX,X,TOL,A,TPS,TNC) C INTEGER N,LOGIC,NAP,NAPMAX DOUBLE PRECISION T,TMIN,TMAX,FN,FPN,A,D2,TOL,AMF,AMD, * TPS,TNC DOUBLE PRECISION XN(N),D(N),G(N),X(N) EXTERNAL SIMUL C C Modification of mlis2 (C. Lemarechal) C INTEGER INDIC,INDICA,INDICD,I DOUBLE PRECISION TA,TG,TD,TC,TP,TESF,TESD,FPG,PS,Z1,FFN,FG,FA,FP, * FPA,FPD,SIGN,DEN,F,BARMIN,BARMUL,BARMAX,BARR, * FD,P,Z,TEST,ANUM,MC,MP C C INDIC=4 TESF=AMF*FPN TESD=AMD*FPN TD=0. TG=0. FG=FN FPG=FPN TA=0. FA=FN FPA=FPN INDICA=1 LOGIC=0 BARMIN=0.01 BARMUL=3. BARMAX=.3 BARR=BARMIN IF (T.GT.TMIN) GO TO 20 T=TMIN IF (T.LE.TMAX) GO TO 20 TMIN=TMAX 20 IF (FN+T*FPN.LT.FN+0.9*T*FPN) GO TO 30 T=2.*T GO TO 20 30 IF (T.GE.TMAX) THEN T=TMAX LOGIC=1 ENDIF DO 50 I=1,N 50 X(I)=XN(I)+T*D(I) 100 NAP=NAP+1 IF (NAP.GT.NAPMAX) THEN LOGIC=4 IF (TG.EQ.0.) GO TO 999 DO 120 I=1,N 120 XN(I)=XN(I)+TG*D(I) INDIC=4 CALL SIMUL (N,XN,FN,G) GO TO 999 ENDIF INDIC=4 CALL SIMUL(N,X,F,G) IF (INDIC.EQ.0) THEN LOGIC=5 FN=F DO 170 I=1,N 170 XN(I)=X(I) GO TO 999 ENDIF 200 IF(INDIC.GT.0) GO TO 210 TD=T INDICD=INDIC LOGIC=0 T=TG+0.1*(TD-TG) GO TO 905 210 PS= 0.0 DO 211 I=1,N PS= PS+G(I)*D(I) 211 CONTINUE FP=PS FFN=F-FN IF(FFN.LT.T*TESF) GO TO 300 TD=T FD=F FPD=FP INDICD=INDIC LOGIC=0 IF(TG.NE.0.) GO TO 500 IF(FPD.LT.TESD) GO TO 500 TPS=(FN-F)+TD*FPD TNC=D2*TD*TD P=DMAX1(A*TNC,TPS) IF(P.GT.TOL) GO TO 500 LOGIC=3 FPN=FPD GO TO 999 300 CONTINUE IF(FP.LT.TESD) GO TO 320 LOGIC=0 FN=F FPN=FP DO 310 I=1,N 310 XN(I)=X(I) GO TO 999 320 IF (LOGIC.EQ.0) GO TO 350 FN=F FPN=FP DO 330 I=1,N 330 XN(I)=X(I) GO TO 999 350 TG=T FG=F FPG=FP TA=T T=9.*TG Z=FPN+3.*FP-4.*FFN/TG IF(Z.GT.0.) T=DMIN1(T,TG*DMAX1(1.0D0,-FP/Z)) T=TG+T IF(T.LT.TMAX) GO TO 900 LOGIC=1 T=TMAX GO TO 900 500 IF(INDICA.GT.0 .AND. INDICD.GT.0) GO TO 510 TA=T T=0.9*TG+0.1*TD GO TO 900 510 TEST=BARR*(TD-TG) PS=FP+FPA-3.D0*(FA-F)/(TA-T) Z1=PS*PS-FP*FPA IF (Z1.GE.0.D0) GO TO 520 IF (FP.LT.0.) TC=TD IF (FP.GE.0.) TC=TG GO TO 600 520 Z1=DSQRT(Z1) IF (T-TA.LT.0.) Z1=-Z1 SIGN=(T-TA)/ABS(T-TA) IF ((PS+FP)*SIGN.GT.0.) GO TO 550 DEN=2.D0*PS+FP+FPA ANUM=Z1-FP-PS IF (ABS((T-TA)*ANUM).GE.(TD-TG)*ABS(DEN)) GO TO 530 TC=T+ANUM*(TA-T)/DEN GO TO 600 530 TC=TD GO TO 600 550 TC=T+FP*(TA-T)/(PS+FP+Z1) 600 MC=0 IF (TC.LT.TG) MC=-1 IF (TC.GT.TD) MC=1 TC=DMAX1(TC,TG+TEST) TC=DMIN1(TC,TD-TEST) PS=FPD-FPG IF (PS.NE.0.D0) GO TO 620 TP=0.5*(TD+TG) GO TO 650 620 TP=((FG-FPG*TG)- * (FD-FPD*TD))/PS 650 MP=0 IF (TP.LT.TG) MP=-1 IF (TP.GT.TD) MP=1 TP=DMAX1(TP,TG+TEST) TP=DMIN1(TP,TD-TEST) TA=T IF (MC.EQ.0 .AND. MP.EQ.0) T=DMIN1(TC,TP) IF (MC.EQ.0 .AND. MP.NE.0) T=TC IF (MC.NE.0 .AND. MP.EQ.0) T=TP IF (MC.EQ.1 .AND. MP.EQ.1) T=TD-TEST IF (MC.EQ.-1 .AND. MP.EQ.-1) T=TG+TEST IF (MC*MP.EQ.-1) T=0.5*(TG+TD) IF (MC.EQ.0 .AND. MP.EQ.0) THEN BARR=BARMIN ELSE BARR=DMIN1(BARMUL*BARR,BARMAX) ENDIF 900 FA=F FPA=FP 905 INDICA=INDIC 920 IF (TD.EQ.0.) GO TO 990 IF (TD-TG.LE.TMIN) GO TO 950 DO 930 I=1,N Z=XN(I)+T*D(I) IF (Z.NE.X(I) .AND. Z.NE.XN(I)) GO TO 990 930 CONTINUE 950 LOGIC=6 IF (INDICD.LT.0) LOGIC=INDICD IF (TG.EQ.0.) GO TO 970 DO 960 I=1,N 960 XN(I)=XN(I)+TG*D(I) INDIC=4 CALL SIMUL(N,XN,FN,G) 970 CONTINUE RETURN 990 DO 995 I=1,N 995 X(I)=XN(I)+T*D(I) GO TO 100 999 RETURN END C VERSION 2 DOES NOT USE EISPACK C C ------------------------------------------------------------------ C SUBROUTINE DNLASO(OP, IOVECT, N, NVAL, NFIG, NPERM, * NMVAL, VAL, NMVEC, VEC, NBLOCK, MAXOP, MAXJ, WORK, * IND, IERR) C INTEGER N, NVAL, NFIG, NPERM, NMVAL, NMVEC, NBLOCK, * MAXOP, MAXJ, IND(1), IERR DOUBLE PRECISION 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 DNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF C THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX. THE SUBROUTINE C DNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE DNWLA WHICH IMPLEMENTS C THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND C SUBROUTINE DNPPLA WHICH POST PROCESSES THE OUTPUT OF DNWLA. C HOWEVER DNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING C PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS. C DNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES. IF C THE LARGEST EIGENVALUES ARE DESIRED DNLASO IMPLICITLY USES THE C NEGATIVE OF THE MATRIX. 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 NVAL NVAL SPECIFIES THE EIGENVALUES TO BE FOUND. C DABS(NVAL) IS THE NUMBER OF EIGENVALUES DESIRED. C IF NVAL < 0 THE ALGEBRAICALLY SMALLEST (LEFTMOST) C EIGENVALUES ARE FOUND. IF NVAL > 0 THE ALGEBRAICALLY C LARGEST (RIGHTMOST) EIGENVALUES ARE FOUND. NVAL MUST NOT C BE ZERO. DABS(NVAL) MUST BE LESS THAN MAXJ/2. 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 DABS(NVAL). 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 VEC A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW C DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST DABS(NVAL). IF C NPERM > 0 THEN THE FIRST NPERM COLUMNS OF VEC MUST C CONTAIN 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. DNLASO 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+NV+2) + 2*NBLOCK*NBLOCK + 3*NV C C + THE MAXIMUM OF C N*NBLOCK C AND C MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1) C C WHERE NV = DABS(NVAL) 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 DABS(NVAL). 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 CONTAINS C ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS. SEE C DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES. C C WORK IF WORK IS TERMINATED BEFORE COMPLETION (IERR = -2) C THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS C FOR RESTARTING THE ALGORITHM AND DNLASO 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, NV, IABS, MAX0 LOGICAL RARITZ, SMALL DOUBLE PRECISION DELTA, EPS, DDOT, DNRM2, DABS, TARR(1) C EXTERNAL DNPPLA, DNWLA, DMVPC, DORTQR, DAXPY, EXTERNAL DNPPLA, DMVPC, DORTQR, DAXPY, * DCOPY, DDOT, DNRM2, DSCAL, DSWAP, DLAEIG, DLABCM, * DLABFC, DLAGER, DLARAN, DVSORT C C NOP RETURNED FROM DNWLA AS THE NUMBER OF CALLS TO THE C SUBROUTINE OP. C C NV SET EQUAL TO DABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED, C AND PASSED TO DNWLA. C C SMALL SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED. C C RARITZ RETURNED FROM DNWLA AND PASSED TO DNPPLA. RARITZ IS .TRUE. C IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED. C C DELTA RETURNED FROM DNWLA AS THE EIGENVALUE OF THE MATRIX C WHICH IS CLOSEST TO THE DESIRED EIGENVALUES. C C DNPPLA A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED C BY DNWLA. C C DNWLA 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 GIVEN RITZ VECTORS. 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 DLARAN A SUBROUTINE TO GENERATE RANDOM VECTORS C C DLAEIG, DLAGER, DLABCM, DLABFC SUBROUTINES FOR BAND EIGENVALUE C CALCULATIONS. C C ------------------------------------------------------------------ C C THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS. C NV = IABS(NVAL) IND(1) = 0 IERR = 0 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 (NV.LT.MAX0(1,NPERM)) IERR = IERR + 32 IF (NV.GT.NMVAL) IERR = IERR + 64 IF (NV.GT.MAXOP) IERR = IERR + 128 IF (NV.GE.MAXJ/2) IERR = IERR + 256 IF (NBLOCK.LT.1) IERR = IERR + 512 IF (IERR.NE.0) RETURN C SMALL = NVAL.LT.0 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 DNLASO TERMINATES. C IF (NPERM.EQ.0) GO TO 110 C C THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST C EIGENVALUES ARE DESIRED, SINCE DNWLA WILL IMPLICITLY USE THE C NEGATIVE OF THE MATRIX. C IF (SMALL) GO TO 20 DO 10 I=1,NPERM VAL(I,1) = -VAL(I,1) 10 CONTINUE C C THIS SORTS THE USER SUPPLIED VALUES AND VECTORS. C 20 CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TARR, NMVEC, N, VEC) C C THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON. C IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE. C DO 60 I=1,NPERM VAL(I,2) = DABS(VAL(I,2)) VAL(I,3) = DNRM2(N,VEC(1,I),1) 60 CONTINUE C C THIS PERFORMS THE ORTHONORMALIZATION. C M = N*NBLOCK + 1 CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M)) M = N*NBLOCK - NPERM DO 70 I = 1, NPERM M = M + NPERM + 1 IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 70 IERR = -1 RETURN C 70 CONTINUE C C THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN C THE ARRAY WORK FOR LATER REFERENCE IN DNWLA. 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 RELATIVE MACHINE C PRECISION C C ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT C ***IN A PRODUCTION CODE C 110 EPS = 1.0D0 c DO 120 I = 1,1000 c EPS = 0.5D0*EPS c TEMP = 1.0D0 + EPS c IF(TEMP.EQ.1.0D0) GO TO 130 c 120 CONTINUE eps = 1.0d-14 C C ------------------------------------------------------------------ C C THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM C WITH SELECTIVE ORTHOGONALIZATION. C 130 NBAND = NBLOCK + 1 I1 = 1 + N*NBLOCK I2 = I1 + N*NBLOCK I3 = I2 + NV I4 = I3 + NV I5 = I4 + NV I6 = I5 + MAXJ*NBAND I7 = I6 + NBLOCK*NBLOCK I8 = I7 + NBLOCK*NBLOCK I9 = I8 + MAXJ*(NV+1) I10 = I9 + NBLOCK I11 = I10 + 2*NV + 6 I12 = I11 + MAXJ*(2*NBLOCK+1) I13 = I12 + MAXJ CALL DNWLA(OP, IOVECT, N, NBAND, NV, NFIG, NPERM, VAL, NMVEC, * 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, SMALL, RARITZ, DELTA, EPS, IERR) C C ------------------------------------------------------------------ C C THIS SECTION CALLS DNPPLA (THE POST PROCESSOR). C IF (NPERM.EQ.0) GO TO 140 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 CALL DNPPLA(OP, IOVECT, N, NPERM, NOP, NMVAL, VAL, NMVEC, * VEC, NBLOCK, WORK(I1), WORK(I2), WORK(I3), WORK(I4), * WORK(I5), WORK(I6), DELTA, SMALL, RARITZ, EPS) C 140 IND(1) = NOP RETURN END C C ------------------------------------------------------------------ C SUBROUTINE DNWLA(OP, IOVECT, N, NBAND, NVAL, NFIG, NPERM, VAL, * NMVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP, P1, P0, * RES, TAU, OTAU, T, ALP, BET, S, P2, BOUND, ATEMP, VTEMP, D, * IND, SMALL, RARITZ, DELTA, EPS, IERR) C INTEGER N, NBAND, NVAL, NFIG, NPERM, NMVEC, NBLOCK, MAXOP, MAXJ, * NOP, IND(1), IERR LOGICAL RARITZ, SMALL DOUBLE PRECISION VAL(1), VEC(NMVEC,1), P0(N,1), P1(N,1), * P2(N,1), RES(1), TAU(1), OTAU(1), T(NBAND,1), * ALP(NBLOCK,1), BET(NBLOCK,1), BOUND(1), ATEMP(1), * VTEMP(1), D(1), S(MAXJ,1), DELTA, EPS EXTERNAL OP, IOVECT C C DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE C ORTHOGONALIZATION. C C NBAND NBLOCK + 1 THE BAND WIDTH OF T. C C NVAL THE NUMBER OF DESIRED EIGENVALUES. C C NPERM THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS C INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE C ALGORITHM IS ITERATED). PERMANENT VECTORS ARE ORTHOGONAL C TO THE CURRENT KRYLOV SUBSPACE. C C NOP THE NUMBER OF CALLS TO OP. C C P0, P1, AND P2 THE CURRENT BLOCKS OF LANCZOS VECTORS. C C RES THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS. C C TAU AND OTAU USED TO MONITOR THE NEED FOR ORTHOGONALIZATION. C C T THE BAND MATRIX. C C ALP THE CURRENT DIAGONAL BLOCK. C C BET THE CURRENT OFF DIAGONAL BLOCK. C C BOUND, ATEMP, VTEMP, D TEMPORARY STORAGE USED BY THE BAND C EIGENVALUE SOLVER DLAEIG. C C S EIGENVECTORS OF T. C C SMALL .TRUE. IF THE SMALL EIGENVALUES ARE DESIRED. C C RARITZ RETURNED AS .TRUE. IF A FINAL RAYLEIGH-RITZ PROCEDURE C IS TO BE DONE. C C DELTA RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE C OF THE MATRIX. USED IN ESTIMATING THE ACCURACY OF THE C COMPUTED EIGENVALUES. C C C INTERNAL VARIABLES USED C INTEGER I, I1, II, J, K, L, M, NG, NGOOD, * NLEFT, NSTART, NTHETA, NUMBER, MIN0, MTEMP LOGICAL ENOUGH, TEST DOUBLE PRECISION ALPMAX, ALPMIN, ANORM, BETMAX, BETMIN, * EPSRT, PNORM, RNORM, TEMP, * TMAX, TMIN, TOLA, TOLG, UTOL, DABS, * DMAX1, DMIN1, DSQRT, DDOT, DNRM2, TARR(1), DZERO(1) EXTERNAL DMVPC, DORTQR, DAXPY, DCOPY, DDOT, * DNRM2, DSCAL, DSWAP, DLAEIG, DLAGER, DLARAN, DVSORT C C J THE CURRENT DIMENSION OF T. (THE DIMENSION OF THE CURRENT C KRYLOV SUBSPACE. C C NGOOD THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS C LIE IN THE CURRENT KRYLOV SUBSPACE). C C NLEFT THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED, C I.E. NLEFT = NVAL - NPERM. C C NUMBER = NPERM + NGOOD. C C ANORM AN ESTIMATE OF THE NORM OF THE MATRIX. C C EPS THE RELATIVE MACHINE PRECISION. C C UTOL THE USER TOLERANCE. C C TARR AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN CALLS TO C DLAEIG C C DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE C CONSISTENCY IN CALLS TO DCOPY C DZERO(1) = 0.0D0 RNORM = 0.0D0 IF (NPERM.NE.0) RNORM = DMAX1(-VAL(1),VAL(NPERM)) PNORM = RNORM DELTA = 10.D28 EPSRT = DSQRT(EPS) NLEFT = NVAL - NPERM NOP = 0 NUMBER = NPERM RARITZ = .FALSE. UTOL = DMAX1(DBLE(FLOAT(N))*EPS,10.0D0**DBLE((-FLOAT(NFIG)))) J = MAXJ C C ------------------------------------------------------------------ C C ANY ITERATION OF THE ALGORITHM BEGINS HERE. C 30 DO 50 I=1,NBLOCK TEMP = DNRM2(N,P1(1,I),1) IF (TEMP.EQ.0D0) CALL DLARAN(N, P1(1,I)) 50 CONTINUE IF (NPERM.EQ.0) GO TO 70 DO 60 I=1,NPERM TAU(I) = 1.0D0 OTAU(I) = 0.0D0 60 CONTINUE 70 CALL DCOPY(N*NBLOCK, DZERO, 0, P0, 1) CALL DCOPY(NBLOCK*NBLOCK, DZERO, 0, BET, 1) CALL DCOPY(J*NBAND, DZERO, 0, T, 1) MTEMP = NVAL + 1 DO 75 I = 1, MTEMP CALL DCOPY(J, DZERO, 0, S(1,I), 1) 75 CONTINUE NGOOD = 0 TMIN = 1.0D28 TMAX = -1.0D28 TEST = .TRUE. ENOUGH = .FALSE. BETMAX = 0.0D0 J = 0 C C ------------------------------------------------------------------ C C THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. C 80 J = J + NBLOCK C C THIS IS THE SELECTIVE ORTHOGONALIZATION. C IF (NUMBER.EQ.0) GO TO 110 DO 100 I=1,NUMBER IF (TAU(I).LT.EPSRT) GO TO 100 TEST = .TRUE. TAU(I) = 0.0D0 IF (OTAU(I).NE.0.0D0) OTAU(I) = 1.0D0 DO 90 K=1,NBLOCK TEMP = -DDOT(N,VEC(1,I),1,P1(1,K),1) CALL DAXPY(N, TEMP, VEC(1,I), 1, P1(1,K), 1) C C THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A C NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR. THE ALGORITHM IS C TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. C IF (DABS(TEMP*BET(K,K)).GT.DBLE(FLOAT(N))*EPSRT* * ANORM .AND. I.GT.NPERM) GO TO 380 90 CONTINUE 100 CONTINUE C C IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. C 110 IF(.NOT. TEST) GO TO 160 CALL DORTQR(N, N, NBLOCK, P1, ALP) TEST = .FALSE. IF(J .EQ. NBLOCK) GO TO 160 DO 130 I = 1,NBLOCK IF(ALP(I,I) .GT. 0.0D0) GO TO 130 M = J - 2*NBLOCK + I L = NBLOCK + 1 DO 120 K = I,NBLOCK BET(I,K) = -BET(I,K) T(L,M) = -T(L,M) L = L - 1 M = M + 1 120 CONTINUE 130 CONTINUE C C THIS IS THE LANCZOS STEP. C 160 CALL OP(N, NBLOCK, P1, P2) NOP = NOP + 1 CALL IOVECT(N, NBLOCK, P1, J, 0) C C THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) C DO 180 I=1,NBLOCK DO 170 K=I,NBLOCK CALL DAXPY(N, -BET(I,K), P0(1,K), 1, P2(1,I), 1) 170 CONTINUE 180 CONTINUE C C THIS COMPUTES ALP AND P2=P2-P1*ALP. C DO 200 I=1,NBLOCK DO 190 K=1,I II = I - K + 1 ALP(II,K) = DDOT(N,P1(1,I),1,P2(1,K),1) CALL DAXPY(N, -ALP(II,K), P1(1,I), 1, P2(1,K), 1) IF (K.NE.I) CALL DAXPY(N, -ALP(II,K), P1(1,K), * 1, P2(1,I), 1) 190 CONTINUE 200 CONTINUE C C REORTHOGONALIZATION OF THE SECOND BLOCK C IF(J .NE. NBLOCK) GO TO 220 DO 215 I=1,NBLOCK DO 210 K=1,I TEMP = DDOT(N,P1(1,I),1,P2(1,K),1) CALL DAXPY(N, -TEMP, P1(1,I), 1, P2(1,K), 1) IF (K.NE.I) CALL DAXPY(N, -TEMP, P1(1,K), * 1, P2(1,I), 1) II = I - K + 1 ALP(II,K) = ALP(II,K) + TEMP 210 CONTINUE 215 CONTINUE C C THIS ORTHONORMALIZES THE NEXT BLOCK C 220 CALL DORTQR(N, N, NBLOCK, P2, BET) C C THIS STORES ALP AND BET IN T. C DO 250 I=1,NBLOCK M = J - NBLOCK + I DO 230 K=I,NBLOCK L = K - I + 1 T(L,M) = ALP(L,I) 230 CONTINUE DO 240 K=1,I L = NBLOCK - I + K + 1 T(L,M) = BET(K,I) 240 CONTINUE 250 CONTINUE C C THIS NEGATES T IF SMALL IS FALSE. C IF (SMALL) GO TO 280 M = J - NBLOCK + 1 DO 270 I=M,J DO 260 K=1,L T(K,I) = -T(K,I) 260 CONTINUE 270 CONTINUE C C THIS SHIFTS THE LANCZOS VECTORS C 280 CALL DCOPY(NBLOCK*N, P1, 1, P0, 1) CALL DCOPY(NBLOCK*N, P2, 1, P1, 1) CALL DLAGER(J, NBAND, J-NBLOCK+1, T, TMIN, TMAX) ANORM = DMAX1(RNORM, TMAX, -TMIN) IF (NUMBER.EQ.0) GO TO 305 C C THIS COMPUTES THE EXTREME EIGENVALUES OF ALP. C CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK, 1 P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) ALPMIN = TARR(1) CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR, 1 NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) ALPMAX = TARR(1) C C THIS COMPUTES ALP = BET(TRANSPOSE)*BET. C 305 DO 310 I = 1, NBLOCK DO 300 K = 1, I L = I - K + 1 ALP(L,K) = DDOT(NBLOCK-I+1, BET(I,I), NBLOCK, BET(K,I), 1 NBLOCK) 300 CONTINUE 310 CONTINUE IF(NUMBER .EQ. 0) GO TO 330 C C THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET. C CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK, 1 P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, ANORM*ANORM) BETMIN = DSQRT(TARR(1)) C C THIS UPDATES TAU AND OTAU. C DO 320 I=1,NUMBER TEMP = (TAU(I)*DMAX1(ALPMAX-VAL(I),VAL(I)-ALPMIN) * +OTAU(I)*BETMAX+EPS*ANORM)/BETMIN IF (I.LE.NPERM) TEMP = TEMP + RES(I)/BETMIN OTAU(I) = TAU(I) TAU(I) = TEMP 320 CONTINUE C C THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET. C 330 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR, 1 NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, 2 ANORM*ANORM) BETMAX = DSQRT(TARR(1)) IF (J.LE.2*NBLOCK) GO TO 80 C C ------------------------------------------------------------------ C C THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND C LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK C IS JUSTIFIED. C TOLG = EPSRT*ANORM TOLA = UTOL*RNORM IF(MAXJ-J .LT. NBLOCK .OR. (NOP .GE. MAXOP .AND. 1 NLEFT .NE. 0)) GO TO 390 GO TO 400 C C ------------------------------------------------------------------ C C THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO C SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN C ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY C IS AT 400. C 380 J = J - NBLOCK IERR = -8 390 IF (NLEFT.EQ.0) RETURN TEST = .TRUE. 400 NTHETA = MIN0(J/2,NLEFT+1) CALL DLAEIG(J, NBAND, 1, NTHETA, T, VAL(NUMBER+1), 1 MAXJ, S, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) CALL DMVPC(NBLOCK, BET, MAXJ, J, S, NTHETA, ATEMP, VTEMP, D) C C THIS CHECKS FOR TERMINATION OF A CHECK RUN C IF(NLEFT .NE. 0 .OR. J .LT. 6*NBLOCK) GO TO 410 IF(VAL(NUMBER+1)-ATEMP(1) .GT. VAL(NPERM) - TOLA) GO TO 790 C C THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T C TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED. C 410 IF (NTHETA.LE.NLEFT) GO TO 470 IF (NPERM.EQ.0) GO TO 430 M = NUMBER + NLEFT + 1 IF (VAL(M).GE.VAL(NPERM)) GO TO 430 NPERM = NPERM - 1 NGOOD = 0 NUMBER = NPERM NLEFT = NLEFT + 1 GO TO 400 C C THIS UPDATES DELTA. C 430 M = NUMBER + NLEFT + 1 DELTA = DMIN1(DELTA,VAL(M)) ENOUGH = .TRUE. IF(NLEFT .EQ. 0) GO TO 80 NTHETA = NLEFT VTEMP(NTHETA+1) = 1 C C ------------------------------------------------------------------ C C THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL. C C THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES. C IF (.NOT.(TEST .OR. ENOUGH)) GO TO 470 DELTA = DMIN1(DELTA,ANORM) PNORM = DMAX1(RNORM,DMAX1(-VAL(NUMBER+1),DELTA)) TOLA = UTOL*PNORM NSTART = 0 DO 460 I=1,NTHETA M = NUMBER + I IF (DMIN1(ATEMP(I)*ATEMP(I)/(DELTA-VAL(M)),ATEMP(I)) * .GT.TOLA) GO TO 450 IND(I) = -1 GO TO 460 C 450 ENOUGH = .FALSE. IF (.NOT.TEST) GO TO 470 IND(I) = 1 NSTART = NSTART + 1 460 CONTINUE C C COPY VALUES OF IND INTO VTEMP C DO 465 I = 1,NTHETA VTEMP(I) = DBLE(FLOAT(IND(I))) 465 CONTINUE GO TO 500 C C THIS CHECKS FOR NEW GOOD VECTORS. C 470 NG = 0 DO 490 I=1,NTHETA IF (VTEMP(I).GT.TOLG) GO TO 480 NG = NG + 1 VTEMP(I) = -1 GO TO 490 C 480 VTEMP(I) = 1 490 CONTINUE C IF (NG.LE.NGOOD) GO TO 80 NSTART = NTHETA - NG C C ------------------------------------------------------------------ C C THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS. C IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED. C 500 TEST = TEST .AND. .NOT.ENOUGH NGOOD = NTHETA - NSTART NSTART = NSTART + 1 NTHETA = NTHETA + 1 C C THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND C EIGENVECTORS OF T. THE OTHER EIGENVECTORS ARE SAVED FOR C FORMING STARTING VECTORS, IF NECESSARY. IT ALSO SHIFTS THE C EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS C PAUSE. C CALL DCOPY(NTHETA, VAL(NUMBER+1), 1, VAL(NPERM+1), 1) IF (NSTART.EQ.0) GO TO 580 IF (NSTART.EQ.NTHETA) GO TO 530 CALL DVSORT(NTHETA, VTEMP, ATEMP, 1, VAL(NPERM+1), MAXJ, * J, S) C C THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING C VECTORS. C 530 IF (.NOT.TEST) NSTART = 0 IF (.NOT.TEST) GO TO 580 C C FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW C TEMP = ATEMP(1) DO 535 I = 1, NSTART TEMP = DMIN1(TEMP, ATEMP(I)) 535 CONTINUE M = NGOOD + 1 L = NGOOD + MIN0(NSTART,NBLOCK) DO 540 I=M,L CALL DSCAL(J, TEMP/ATEMP(I), S(1,I), 1) 540 CONTINUE M = (NSTART-1)/NBLOCK IF (M.EQ.0) GO TO 570 L = NGOOD + NBLOCK DO 560 I=1,M DO 550 K=1,NBLOCK L = L + 1 IF (L.GT.NTHETA) GO TO 570 I1 = NGOOD + K CALL DAXPY(J, TEMP/ATEMP(L), S(1,L), 1, S(1,I1), 1) 550 CONTINUE 560 CONTINUE 570 NSTART = MIN0(NSTART,NBLOCK) C C THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS. C 580 IF (NGOOD.EQ.0 .OR. .NOT.(TEST .OR. ENOUGH)) GO TO 600 DO 590 I=1,NGOOD M = NPERM + I RES(M) = ATEMP(I) 590 CONTINUE C C THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE C LANCZOS VECTORS. C 600 NUMBER = NPERM + NGOOD IF (TEST .OR. ENOUGH) CALL DCOPY(N*NBLOCK, DZERO, 0, P1, 1) IF (NGOOD.EQ.0) GO TO 620 M = NPERM + 1 DO 610 I=M,NUMBER CALL DCOPY(N, DZERO, 0, VEC(1,I), 1) 610 CONTINUE 620 DO 670 I=NBLOCK,J,NBLOCK CALL IOVECT(N, NBLOCK, P2, I, 1) DO 660 K=1,NBLOCK M = I - NBLOCK + K IF (NSTART.EQ.0) GO TO 640 DO 630 L=1,NSTART I1 = NGOOD + L CALL DAXPY(N, S(M,I1), P2(1,K), 1, P1(1,L), 1) 630 CONTINUE 640 IF (NGOOD.EQ.0) GO TO 660 DO 650 L=1,NGOOD I1 = L + NPERM CALL DAXPY(N, S(M,L), P2(1,K), 1, VEC(1,I1), 1) 650 CONTINUE 660 CONTINUE 670 CONTINUE IF (TEST .OR. ENOUGH) GO TO 690 C C THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE C TAU RECURRENCE. C M = NPERM + 1 DO 680 I=M,NUMBER TEMP = 1.0D0/DNRM2(N,VEC(1,I),1) CALL DSCAL(N, TEMP, VEC(1,I), 1) TAU(I) = 1.0D0 OTAU(I) = 1.0D0 680 CONTINUE C C SHIFT S VECTORS TO ALIGN FOR LATER CALL TO DLAEIG C CALL DCOPY(NTHETA, VAL(NPERM+1), 1, VTEMP, 1) CALL DVSORT(NTHETA, VTEMP, ATEMP, 0, TARR, MAXJ, J, S) GO TO 80 C C ------------------------------------------------------------------ C C THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE C PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING C THE PERMANENT VECTORS. C 690 IF (NGOOD.EQ.0 .AND. NOP.GE.MAXOP) GO TO 810 IF (NGOOD.EQ.0) GO TO 30 C C THIS ORTHONORMALIZES THE VECTORS C CALL DORTQR(NMVEC, N, NPERM+NGOOD, VEC, S) C C THIS SORTS THE VALUES AND VECTORS. C IF(NPERM .NE. 0) CALL DVSORT(NPERM+NGOOD, VAL, RES, 0, TEMP, * NMVEC, N, VEC) NPERM = NPERM + NGOOD NLEFT = NLEFT - NGOOD RNORM = DMAX1(-VAL(1),VAL(NPERM)) C C THIS DECIDES WHERE TO GO NEXT. C IF (NOP.GE.MAXOP .AND. NLEFT.NE.0) GO TO 810 IF (NLEFT.NE.0) GO TO 30 IF (VAL(NVAL)-VAL(1).LT.TOLA) GO TO 790 C C THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED C TO LOOK FOR UNDISCLOSED MULTIPLICITIES. C M = NPERM - NBLOCK + 1 IF (M.LE.0) RETURN DO 780 I=1,M L = I + NBLOCK - 1 IF (VAL(L)-VAL(I).LT.TOLA) GO TO 30 780 CONTINUE C C THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ C PROCEDURE IS NEEDED. C 790 M = NPERM - NBLOCK IF (M.LE.0) RETURN DO 800 I=1,M L = I + NBLOCK IF (VAL(L)-VAL(I).GE.TOLA) GO TO 800 RARITZ = .TRUE. RETURN 800 CONTINUE C RETURN C C ------------------------------------------------------------------ C C THIS REPORTS THAT MAXOP WAS EXCEEDED. C 810 IERR = -2 GO TO 790 C END C C ------------------------------------------------------------------ C SUBROUTINE DNPPLA(OP, IOVECT, N, NPERM, NOP, NMVAL, VAL, * NMVEC, VEC, NBLOCK, H, HV, P, Q, BOUND, D, DELTA, SMALL, * RARITZ, EPS) C INTEGER N, NPERM, NOP, NMVAL, NMVEC, NBLOCK LOGICAL SMALL, RARITZ DOUBLE PRECISION VAL(NMVAL,1), VEC(NMVEC,1), H(NPERM,1), * HV(NPERM,1), P(N,1), Q(N,1), BOUND(1), D(1), DELTA, EPS EXTERNAL OP, IOVECT C C THIS SUBROUTINE POST PROCESSES THE EIGENVECTORS. BLOCK MATRIX C VECTOR PRODUCTS ARE USED TO MINIMIZED THE NUMBER OF CALLS TO OP. C INTEGER I, J, JJ, K, KK, L, M, MOD DOUBLE PRECISION HMIN, HMAX, TEMP, DDOT, DNRM2, DZERO(1) EXTERNAL DAXPY, DCOPY, DDOT, DLAGER, DLAEIG C C IF RARITZ IS .TRUE. A FINAL RAYLEIGH-RITZ PROCEDURE IS APPLIED C TO THE EIGENVECTORS. C DZERO(1) = 0.0D0 IF (.NOT.RARITZ) GO TO 190 C C ------------------------------------------------------------------ C C THIS CONSTRUCTS H=Q*AQ, WHERE THE COLUMNS OF Q ARE THE C APPROXIMATE EIGENVECTORS. TEMP = -1 IS USED WHEN SMALL IS C FALSE TO AVOID HAVING TO RESORT THE EIGENVALUES AND EIGENVECTORS C COMPUTED BY DLAEIG. C CALL DCOPY(NPERM*NPERM, DZERO, 0, H, 1) TEMP = -1.0D0 IF (SMALL) TEMP = 1.0D0 M = MOD(NPERM,NBLOCK) IF (M.EQ.0) GO TO 40 DO 10 I=1,M CALL DCOPY(N, VEC(1,I), 1, P(1,I), 1) 10 CONTINUE CALL IOVECT(N, M, P, M, 0) CALL OP(N, M, P, Q) NOP = NOP + 1 DO 30 I=1,M DO 20 J=I,NPERM JJ = J - I + 1 H(JJ,I) = TEMP*DDOT(N,VEC(1,J),1,Q(1,I),1) 20 CONTINUE 30 CONTINUE IF (NPERM.LT.NBLOCK) GO TO 90 40 M = M + NBLOCK DO 80 I=M,NPERM,NBLOCK DO 50 J=1,NBLOCK L = I - NBLOCK + J CALL DCOPY(N, VEC(1,L), 1, P(1,J), 1) 50 CONTINUE CALL IOVECT(N, NBLOCK, P, I, 0) CALL OP(N, NBLOCK, P, Q) NOP = NOP + 1 DO 70 J=1,NBLOCK L = I - NBLOCK + J DO 60 K=L,NPERM KK = K - L + 1 H(KK,L) = TEMP*DDOT(N,VEC(1,K),1,Q(1,J),1) 60 CONTINUE 70 CONTINUE 80 CONTINUE C C THIS COMPUTES THE SPECTRAL DECOMPOSITION OF H. C 90 HMIN = H(1,1) HMAX = H(1,1) CALL DLAGER(NPERM, NPERM, 1, H, HMIN, HMAX) CALL DLAEIG(NPERM, NPERM, 1, NPERM, H, VAL, NPERM, 1 HV, BOUND, P, D, Q, EPS, HMIN, HMAX) C C THIS COMPUTES THE RITZ VECTORS--THE COLUMNS OF C Y = QS WHERE S IS THE MATRIX OF EIGENVECTORS OF H. C DO 120 I=1,NPERM CALL DCOPY(N, DZERO, 0, VEC(1,I), 1) 120 CONTINUE M = MOD(NPERM,NBLOCK) IF (M.EQ.0) GO TO 150 CALL IOVECT(N, M, P, M, 1) DO 140 I=1,M DO 130 J=1,NPERM CALL DAXPY(N, HV(I,J), P(1,I), 1, VEC(1,J), 1) 130 CONTINUE 140 CONTINUE IF (NPERM.LT.NBLOCK) GO TO 190 150 M = M + NBLOCK DO 180 I=M,NPERM,NBLOCK CALL IOVECT(N, NBLOCK, P, I, 1) DO 170 J=1,NBLOCK L = I - NBLOCK + J DO 160 K=1,NPERM CALL DAXPY(N, HV(L,K), P(1,J), 1, VEC(1,K), 1) 160 CONTINUE 170 CONTINUE 180 CONTINUE C C ------------------------------------------------------------------ C C THIS SECTION COMPUTES THE RAYLEIGH QUOTIENTS (IN VAL(*,1)) C AND RESIDUAL NORMS (IN VAL(*,2)) OF THE EIGENVECTORS. C 190 IF (.NOT.SMALL) DELTA = -DELTA M = MOD(NPERM,NBLOCK) IF (M.EQ.0) GO TO 220 DO 200 I=1,M CALL DCOPY(N, VEC(1,I), 1, P(1,I), 1) 200 CONTINUE CALL OP(N, M, P, Q) NOP = NOP + 1 DO 210 I=1,M VAL(I,1) = DDOT(N,P(1,I),1,Q(1,I),1) CALL DAXPY(N, -VAL(I,1), P(1,I), 1, Q(1,I), 1) VAL(I,2) = DNRM2(N,Q(1,I),1) 210 CONTINUE IF (NPERM.LT.NBLOCK) GO TO 260 220 M = M + 1 DO 250 I=M,NPERM,NBLOCK DO 230 J=1,NBLOCK L = I - 1 + J CALL DCOPY(N, VEC(1,L), 1, P(1,J), 1) 230 CONTINUE CALL OP(N, NBLOCK, P, Q) NOP = NOP + 1 DO 240 J=1,NBLOCK L = I - 1 + J VAL(L,1) = DDOT(N,P(1,J),1,Q(1,J),1) CALL DAXPY(N, -VAL(L,1), P(1,J), 1, Q(1,J), 1) VAL(L,2) = DNRM2(N,Q(1,J),1) 240 CONTINUE 250 CONTINUE C C THIS COMPUTES THE ACCURACY ESTIMATES. FOR CONSISTENCY WITH DILASO C A DO LOOP IS NOT USED. C 260 I = 0 270 I = I + 1 IF (I.GT.NPERM) RETURN TEMP = DELTA - VAL(I,1) IF (.NOT.SMALL) TEMP = -TEMP VAL(I,4) = 0.0D0 IF (TEMP.GT.0.0D0) VAL(I,4) = VAL(I,2)/TEMP VAL(I,3) = VAL(I,4)*VAL(I,2) GO TO 270 C END DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER I, INCX, J, N, NEXT, NN DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END C C ------------------------------------------------------------------ C SUBROUTINE DORTQR(NZ, N, NBLOCK, Z, B) C INTEGER NZ, N, NBLOCK DOUBLE PRECISION Z(NZ,1), B(NBLOCK,1) C C THIS SUBROUTINE COMPUTES THE QR FACTORIZATION OF THE N X NBLOCK C MATRIX Z. Q IS FORMED IN PLACE AND RETURNED IN Z. R IS C RETURNED IN B. C INTEGER I, J, K, LENGTH, M DOUBLE PRECISION SIGMA, TAU, TEMP, DDOT, DNRM2, DSIGN EXTERNAL DAXPY, DDOT, DNRM2, DSCAL C C THIS SECTION REDUCES Z TO TRIANGULAR FORM. C DO 30 I=1,NBLOCK C C THIS FORMS THE ITH REFLECTION. C LENGTH = N - I + 1 SIGMA = DSIGN(DNRM2(LENGTH,Z(I,I),1),Z(I,I)) B(I,I) = -SIGMA Z(I,I) = Z(I,I) + SIGMA TAU = SIGMA*Z(I,I) IF (I.EQ.NBLOCK) GO TO 30 J = I + 1 C C THIS APPLIES THE ROTATION TO THE REST OF THE COLUMNS. C DO 20 K=J,NBLOCK IF (TAU.EQ.0.0D0) GO TO 10 TEMP = -DDOT(LENGTH,Z(I,I),1,Z(I,K),1)/TAU CALL DAXPY(LENGTH, TEMP, Z(I,I), 1, Z(I,K), 1) 10 B(I,K) = Z(I,K) Z(I,K) = 0.0D0 20 CONTINUE 30 CONTINUE C C THIS ACCUMULATES THE REFLECTIONS IN REVERSE ORDER. C DO 70 M=1,NBLOCK C C THIS RECREATES THE ITH = NBLOCK-M+1)TH REFLECTION. C I = NBLOCK + 1 - M SIGMA = -B(I,I) TAU = Z(I,I)*SIGMA IF (TAU.EQ.0.0D0) GO TO 60 LENGTH = N - NBLOCK + M IF (I.EQ.NBLOCK) GO TO 50 J = I + 1 C C THIS APPLIES IT TO THE LATER COLUMNS. C DO 40 K=J,NBLOCK TEMP = -DDOT(LENGTH,Z(I,I),1,Z(I,K),1)/TAU CALL DAXPY(LENGTH, TEMP, Z(I,I), 1, Z(I,K), 1) 40 CONTINUE 50 CALL DSCAL(LENGTH, -1.0D0/SIGMA, Z(I,I), 1) 60 Z(I,I) = 1.0D0 + Z(I,I) 70 CONTINUE RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DA,DX(1) INTEGER I,INCX,M,MP1,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END C C------------------------------------------------------------------- C SUBROUTINE DVSORT(NUM, VAL, RES, IFLAG, V, NMVEC, N, VEC) INTEGER NUM, IFLAG, NMVEC, N DOUBLE PRECISION VAL(1), RES(1), V(1), VEC(NMVEC,1) C C THIS SUBROUTINE SORTS THE EIGENVALUES (VAL) IN ASCENDING ORDER C WHILE CONCURRENTLY SWAPPING THE RESIDUALS AND VECTORS. INTEGER I, K, M DOUBLE PRECISION TEMP IF(NUM .LE. 1) RETURN DO 20 I = 2, NUM M = NUM - I + 1 DO 10 K = 1, M IF(VAL(K) .LE. VAL(K+1)) GO TO 10 TEMP = VAL(K) VAL(K) = VAL(K+1) VAL(K+1) = TEMP TEMP = RES(K) RES(K) = RES(K+1) RES(K+1) = TEMP CALL DSWAP(N, VEC(1,K), 1, VEC(1,K+1), 1) IF(IFLAG .EQ. 0) GO TO 10 TEMP = V(K) V(K) = V(K+1) V(K+1) = TEMP 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END C C SUBROUTINE DLAEIG(N, NBAND, NL, NR, A, EIGVAL, LDE, 1 EIGVEC, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) C C THIS IS A SPECIALIZED VERSION OF THE SUBROUTINE BNDEIG TAILORED C SPECIFICALLY FOR USE BY THE LASO PACKAGE. C INTEGER N, NBAND, NL, NR, LDE DOUBLE PRECISION A(NBAND,1), EIGVAL(1), 1 EIGVEC(LDE,1), BOUND(2,1), ATEMP(1), D(1), VTEMP(1), 2 EPS, TMIN, TMAX C C LOCAL VARIABLES C INTEGER I, M, NVAL DOUBLE PRECISION ARTOL, ATOL C C FUNCTIONS CALLED C DOUBLE PRECISION DMAX1 C C SUBROUTINES CALLED C C DLABCM, DLABFC, DLAGER, DCOPY C C SET PARAMETERS C ATOL = DBLE(FLOAT(N))*EPS*DMAX1(TMAX,-TMIN) ARTOL = ATOL/DSQRT(EPS) NVAL = NR - NL + 1 C C CHECK FOR SPECIAL CASE OF N = 1 C IF(N .NE. 1) GO TO 30 EIGVAL(1) = A(1,1) EIGVEC(1,1) = 1.0D0 RETURN C C SET UP INITIAL EIGENVALUE BOUNDS C 30 M = NVAL + 1 DO 50 I = 2, M BOUND(1,I) = TMIN BOUND(2,I) = TMAX 50 CONTINUE BOUND(2,1) = TMAX BOUND(1,NVAL + 2) = TMIN IF(NL .EQ. 1) BOUND(2,1) = TMIN IF(NR .EQ. N) BOUND(1,NVAL + 2) = TMAX C 60 CALL DLABCM(N, NBAND, NL, NR, A, EIGVAL, LDE, 1 EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP) RETURN END C C*********************************************************************** C SUBROUTINE DLAGER(N, NBAND, NSTART, A, TMIN, TMAX) C C THIS SUBROUTINE COMPUTES BOUNDS ON THE SPECTRUM OF A BY C EXAMINING THE GERSCHGORIN CIRCLES. ONLY THE NEWLY CREATED C CIRCLES ARE EXAMINED C C FORMAL PARAMETERS C INTEGER N, NBAND, NSTART DOUBLE PRECISION A(NBAND,1), TMIN, TMAX C C LOCAL VARIABLES C INTEGER I, K, L, M DOUBLE PRECISION TEMP C C FUNCTIONS CALLED C INTEGER MIN0 DOUBLE PRECISION DMIN1, DMAX1 C DO 50 K = NSTART, N TEMP = 0.0D0 DO 10 I = 2, NBAND TEMP = TEMP + DABS(A(I,K)) 10 CONTINUE 20 L = MIN0(K,NBAND) IF(L .EQ. 1) GO TO 40 DO 30 I = 2, L M = K - I + 1 TEMP = TEMP + DABS(A(I,M)) 30 CONTINUE 40 TMIN = DMIN1(TMIN, A(1,K)-TEMP) TMAX = DMAX1(TMAX, A(1,K)+TEMP) 50 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE DLARAN(N, X) C C THIS SUBROUTINE SETS THE VECTOR X TO RANDOM NUMBERS C C FORMAL PARAMETERS C INTEGER N DOUBLE PRECISION X(N) C C LOCAL VARIABLES C INTEGER I, IURAND C C FUNCTIONS CALLED C REAL URAND DOUBLE PRECISION DBLE C C SUBROUTINES CALLED C C NONE C C INITIALIZE SEED C c julie 14.3.92- temp or maybe not c common /random/iurand DATA IURAND /0/ C DO 10 I = 1, N X(I) = DBLE(URAND(IURAND)) - 0.5D0 10 CONTINUE RETURN END C C ------------------------------------------------------------------ C SUBROUTINE DMVPC(NBLOCK, BET, MAXJ, J, S, NUMBER, RESNRM, * ORTHCF, RV) C INTEGER NBLOCK, MAXJ, J, NUMBER DOUBLE PRECISION BET(NBLOCK,1), S(MAXJ,1), RESNRM(1), * ORTHCF(1), RV(1) C C THIS SUBROUTINE COMPUTES THE NORM AND THE SMALLEST ELEMENT C (IN ABSOLUTE VALUE) OF THE VECTOR BET*SJI, WHERE SJI C IS AN NBLOCK VECTOR OF THE LAST NBLOCK ELEMENTS OF THE ITH C EIGENVECTOR OF T. THESE QUANTITIES ARE THE RESIDUAL NORM C AND THE ORTHOGONALITY COEFFICIENT RESPECTIVELY FOR THE C CORRESPONDING RITZ PAIR. THE ORTHOGONALITY COEFFICIENT IS C NORMALIZED TO ACCOUNT FOR THE LOCAL REORTHOGONALIZATION. C INTEGER I, K, M DOUBLE PRECISION DDOT, DNRM2, DABS, DMIN1 C M = J - NBLOCK + 1 DO 20 I=1,NUMBER DO 10 K=1,NBLOCK RV(K) = DDOT(NBLOCK,S(M,I),1,BET(K,1),NBLOCK) IF (K.EQ.1) ORTHCF(I) = DABS(RV(K)) ORTHCF(I) = DMIN1(ORTHCF(I),DABS(RV(K))) 10 CONTINUE RESNRM(I) = DNRM2(NBLOCK,RV,1) 20 CONTINUE RETURN END SUBROUTINE DSWAP (N,DX,INCX,DY,INCY) C C INTERCHANGES TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP DTEMP = DX(I + 1) DX(I + 1) = DY(I + 1) DY(I + 1) = DTEMP DTEMP = DX(I + 2) DX(I + 2) = DY(I + 2) DY(I + 2) = DTEMP 50 CONTINUE RETURN END REAL FUNCTION URAND(IY) INTEGER IY C C URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2. THE INTEGER IY C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL C TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED C IN THE INTERVAL (0,1). C INTEGER IA,IC,ITWO,M2,M,MIC DOUBLE PRECISION HALFM REAL S DOUBLE PRECISION DATAN,DSQRT DATA M2/0/,ITWO/2/ IF (M2 .NE. 0) GO TO 20 C C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH C M = 1 10 M2 = M M = ITWO*M2 IF (M .GT. M2) GO TO 10 HALFM = M2 C C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD C IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5 IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1 MIC = (M2 - IC) + M2 C C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT C S = 0.5/HALFM C C COMPUTE NEXT RANDOM NUMBER C 20 IY = IY*IA C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW C INTEGER OVERFLOW ON ADDITION C IF (IY .GT. MIC) IY = (IY - M2) - M2 C IY = IY + IC C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE C WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION C IF (IY/2 .GT. M2) IY = (IY - M2) - M2 C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER C OVERFLOW AFFECTS THE SIGN BIT C IF (IY .LT. 0) IY = (IY + M2) + M2 URAND = FLOAT(IY)*S RETURN END C C*********************************************************************** C SUBROUTINE DLABCM(N, NBAND, NL, NR, A, EIGVAL, 1 LDE, EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP) C C THIS SUBROUTINE ORGANIZES THE CALCULATION OF THE EIGENVALUES C FOR THE BNDEIG PACKAGE. EIGENVALUES ARE COMPUTED BY C A MODIFIED RAYLEIGH QUOTIENT ITERATION. THE EIGENVALUE COUNT C OBTAINED BY EACH FACTORIZATION IS USED TO OCCASIONALLY OVERRIDE C THE COMPUTED RAYLEIGH QUOTIENT WITH A DIFFERENT SHIFT TO C INSURE CONVERGENCE TO THE DESIRED EIGENVALUES. C C FORMAL PARAMETERS. C INTEGER N, NBAND, NL, NR, LDE DOUBLE PRECISION A(NBAND,1), EIGVAL(1), 1 EIGVEC(LDE,1), ATOL, ARTOL, BOUND(2,1), ATEMP(1), 2 D(1), VTEMP(1) C C C LOCAL VARIABLES C LOGICAL FLAG INTEGER I, J, L, M, NUML, NUMVEC, NVAL DOUBLE PRECISION ERRB, GAP, RESID, RQ, SIGMA, VNORM C C C FUNCTIONS CALLED C INTEGER MIN0 DOUBLE PRECISION DMAX1, DMIN1, DDOT, DNRM2 C C SUBROUTINES CALLED C C DLABAX, DLABFC, DLARAN, DAXPY, DCOPY, DSCAL C C REPLACE ZERO VECTORS BY RANDOM C NVAL = NR - NL + 1 FLAG = .FALSE. DO 5 I = 1, NVAL IF(DDOT(N, EIGVEC(1,I), 1, EIGVEC(1,I), 1) .EQ. 0.0D0) 1 CALL DLARAN(N,EIGVEC(1,I)) 5 CONTINUE C C LOOP OVER EIGENVALUES C SIGMA = BOUND(2,NVAL+1) DO 400 J = 1, NVAL NUML = J C C PREPARE TO COMPUTE FIRST RAYLEIGH QUOTIENT C 10 CALL DLABAX(N, NBAND, A, EIGVEC(1,J), VTEMP) VNORM = DNRM2(N, VTEMP, 1) IF(VNORM .EQ. 0.0D0) GO TO 20 CALL DSCAL(N, 1.0D0/VNORM, VTEMP, 1) CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1) CALL DAXPY(N, -SIGMA, EIGVEC(1,J), 1, VTEMP, 1) C C LOOP OVER SHIFTS C C COMPUTE RAYLEIGH QUOTIENT, RESIDUAL NORM, AND CURRENT TOLERANCE C 20 VNORM = DNRM2(N, EIGVEC(1,J), 1) IF(VNORM .NE. 0.0D0) GO TO 30 CALL DLARAN(N, EIGVEC(1,J)) GO TO 10 C 30 RQ = SIGMA + DDOT(N, EIGVEC(1,J), 1, VTEMP, 1) 1 /VNORM/VNORM CALL DAXPY(N, SIGMA-RQ, EIGVEC(1,J), 1, VTEMP, 1) RESID = DMAX1(ATOL, DNRM2(N, VTEMP, 1)/VNORM) CALL DSCAL(N, 1.0/VNORM, EIGVEC(1,J), 1) C C ACCEPT EIGENVALUE IF THE INTERVAL IS SMALL ENOUGH C IF(BOUND(2,J+1) - BOUND(1,J+1) .LT. 3.0D0*ATOL) GO TO 300 C C COMPUTE MINIMAL ERROR BOUND C ERRB = RESID GAP = DMIN1(BOUND(1,J+2) - RQ, RQ - BOUND(2,J)) IF(GAP .GT. RESID) ERRB = DMAX1(ATOL, RESID*RESID/GAP) C C TENTATIVE NEW SHIFT C SIGMA = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1)) C C CHECK FOR TERMINALTION C IF(RESID .GT. 2.0D0*ATOL) GO TO 40 IF(RQ - ERRB .GT. BOUND(2,J) .AND. 1 RQ + ERRB .LT. BOUND(1,J+2)) GO TO 310 C C RQ IS TO THE LEFT OF THE INTERVAL C 40 IF(RQ .GE. BOUND(1,J+1)) GO TO 50 IF(RQ - ERRB .GT. BOUND(2,J)) GO TO 100 IF(RQ + ERRB .LT. BOUND(1,J+1)) CALL DLARAN(N,EIGVEC(1,J)) GO TO 200 C C RQ IS TO THE RIGHT OF THE INTERVAL C 50 IF(RQ .LE. BOUND(2,J+1)) GO TO 100 IF(RQ + ERRB .LT. BOUND(1,J+2)) GO TO 100 C C SAVE THE REJECTED VECTOR IF INDICATED C IF(RQ - ERRB .LE. BOUND(2,J+1)) GO TO 200 DO 60 I = J, NVAL IF(BOUND(2,I+1) .GT. RQ) GO TO 70 60 CONTINUE GO TO 80 C 70 CALL DCOPY(N, EIGVEC(1,J), 1, EIGVEC(1,I), 1) C 80 CALL DLARAN(N, EIGVEC(1,J)) GO TO 200 C C PERTURB RQ TOWARD THE MIDDLE C 100 IF(SIGMA .LT. RQ) SIGMA = DMAX1(SIGMA, RQ-ERRB) IF(SIGMA .GT. RQ) SIGMA = DMIN1(SIGMA, RQ+ERRB) C C FACTOR AND SOLVE C 200 DO 210 I = J, NVAL IF(SIGMA .LT. BOUND(1,I+1)) GO TO 220 210 CONTINUE I = NVAL + 1 220 NUMVEC = I - J NUMVEC = MIN0(NUMVEC, NBAND + 2) IF(RESID .LT. ARTOL) NUMVEC = MIN0(1,NUMVEC) CALL DCOPY(N, EIGVEC(1,J), 1, VTEMP, 1) CALL DLABFC(N, NBAND, A, SIGMA, NUMVEC, LDE, 1 EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL) C C PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW C IF(NUMVEC .EQ. 1) GO TO 227 L = NUMVEC - 1 DO 225 I = 1,L M = J + I CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,M), 1) 225 CONTINUE C C UPDATE INTERVALS C 227 NUML = NUML - NL + 1 IF(NUML .GE. 0) BOUND(2,1) = DMIN1(BOUND(2,1), SIGMA) DO 230 I = J, NVAL IF(SIGMA .LT. BOUND(1,I+1)) GO TO 20 IF(NUML .LT. I) BOUND(1,I+1) = SIGMA IF(NUML .GE. I) BOUND(2,I+1) = SIGMA 230 CONTINUE IF(NUML .LT. NVAL + 1) BOUND(1,NVAL+2) = DMAX1(SIGMA, 1 BOUND(1,NVAL+2)) GO TO 20 C C ACCEPT AN EIGENPAIR C 300 CALL DLARAN(N, EIGVEC(1,J)) FLAG = .TRUE. GO TO 310 C 305 FLAG = .FALSE. RQ = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1)) CALL DLABFC(N, NBAND, A, RQ, NUMVEC, LDE, 1 EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL) VNORM = DNRM2(N, EIGVEC(1,J), 1) IF(VNORM .NE. 0.0) CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1) C C ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES C 310 EIGVAL(J) = RQ IF(J .EQ. 1) GO TO 330 M = J - 1 DO 320 I = 1, M CALL DAXPY(N, -DDOT(N,EIGVEC(1,I),1,EIGVEC(1,J),1), 1 EIGVEC(1,I), 1, EIGVEC(1,J), 1) 320 CONTINUE 330 VNORM = DNRM2(N, EIGVEC(1,J), 1) IF(VNORM .EQ. 0.0D0) GO TO 305 CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1) C C ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE C IF(FLAG) GO TO 305 IF(J .EQ. NVAL) RETURN M = J + 1 DO 340 I = M, NVAL CALL DAXPY(N, -DDOT(N,EIGVEC(1,J),1,EIGVEC(1,I),1), 1 EIGVEC(1,J), 1, EIGVEC(1,I), 1) 340 CONTINUE 400 CONTINUE RETURN C 500 CONTINUE END C C*********************************************************************** C SUBROUTINE DLABFC(N, NBAND, A, SIGMA, NUMBER, LDE, 1 EIGVEC, NUML, LDAD, ATEMP, D, ATOL) C C THIS SUBROUTINE FACTORS (A-SIGMA*I) WHERE A IS A GIVEN BAND C MATRIX AND SIGMA IS AN INPUT PARAMETER. IT ALSO SOLVES ZERO C OR MORE SYSTEMS OF LINEAR EQUATIONS. IT RETURNS THE NUMBER C OF EIGENVALUES OF A LESS THAN SIGMA BY COUNTING THE STURM C SEQUENCE DURING THE FACTORIZATION. TO OBTAIN THE STURM C SEQUENCE COUNT WHILE ALLOWING NON-SYMMETRIC PIVOTING FOR C STABILITY, THE CODE USES A GUPTA'S MULTIPLE PIVOTING C ALGORITHM. C C FORMAL PARAMETERS C INTEGER N, NBAND, NUMBER, LDE, NUML, LDAD DOUBLE PRECISION A(NBAND,1), SIGMA, EIGVEC(LDE,1), 1 ATEMP(LDAD,1), D(LDAD,1), ATOL C C LOCAL VARIABLES C INTEGER I, J, K, KK, L, LA, LD, LPM, M, NB1 DOUBLE PRECISION ZERO(1) C C FUNCTIONS CALLED C INTEGER MIN0 DOUBLE PRECISION DABS C C SUBROUTINES CALLED C C DAXPY, DCOPY, DSWAP C C C INITIALIZE C ZERO(1) = 0.0D0 NB1 = NBAND - 1 NUML = 0 CALL DCOPY(LDAD*NBAND, ZERO, 0, D, 1) C C LOOP OVER COLUMNS OF A C DO 100 K = 1, N C C ADD A COLUMN OF A TO D C D(NBAND, NBAND) = A(1,K) - SIGMA M = MIN0(K, NBAND) - 1 IF(M .EQ. 0) GO TO 20 DO 10 I = 1, M LA = K - I LD = NBAND - I D(LD,NBAND) = A(I+1, LA) 10 CONTINUE C 20 M = MIN0(N-K, NB1) IF(M .EQ. 0) GO TO 40 DO 30 I = 1, M LD = NBAND + I D(LD, NBAND) = A(I+1, K) 30 CONTINUE C C TERMINATE C 40 LPM = 1 IF(NB1 .EQ. 0) GO TO 70 DO 60 I = 1, NB1 L = K - NBAND + I IF(D(I,NBAND) .EQ. 0.0D0) GO TO 60 IF(DABS(D(I,I)) .GE. DABS(D(I,NBAND))) GO TO 50 IF((D(I,NBAND) .LT. 0.0D0 .AND. D(I,I) .LT. 0.0D0) 1 .OR. (D(I,NBAND) .GT. 0.0D0 .AND. D(I,I) .GE. 0.0D0)) 2 LPM = -LPM CALL DSWAP(LDAD-I+1, D(I,I), 1, D(I,NBAND), 1) CALL DSWAP(NUMBER, EIGVEC(L,1), LDE, EIGVEC(K,1), LDE) 50 CALL DAXPY(LDAD-I, -D(I,NBAND)/D(I,I), D(I+1,I), 1, 1 D(I+1,NBAND), 1) CALL DAXPY(NUMBER, -D(I,NBAND)/D(I,I), EIGVEC(L,1), 1 LDE, EIGVEC(K,1), LDE) 60 CONTINUE C C UPDATE STURM SEQUENCE COUNT C 70 IF(D(NBAND,NBAND) .LT. 0.0D0) LPM = -LPM IF(LPM .LT. 0) NUML = NUML + 1 IF(K .EQ. N) GO TO 110 C C COPY FIRST COLUMN OF D INTO ATEMP IF(K .LT. NBAND) GO TO 80 L = K - NB1 CALL DCOPY(LDAD, D, 1, ATEMP(1,L), 1) C C SHIFT THE COLUMNS OF D OVER AND UP C IF(NB1 .EQ. 0) GO TO 100 80 DO 90 I = 1, NB1 CALL DCOPY(LDAD-I, D(I+1,I+1), 1, D(I,I), 1) D(LDAD,I) = 0.0D0 90 CONTINUE 100 CONTINUE C C TRANSFER D TO ATEMP C 110 DO 120 I = 1, NBAND L = N - NBAND + I CALL DCOPY(NBAND-I+1, D(I,I), 1, ATEMP(1,L), 1) 120 CONTINUE C C BACK SUBSTITUTION C IF(NUMBER .EQ. 0) RETURN DO 160 KK = 1, N K = N - KK + 1 IF(DABS(ATEMP(1,K)) .LE. ATOL) 1 ATEMP(1,K) = DSIGN(ATOL,ATEMP(1,K)) C 130 DO 150 I = 1, NUMBER EIGVEC(K,I) = EIGVEC(K,I)/ATEMP(1,K) M = MIN0(LDAD, K) - 1 IF(M .EQ. 0) GO TO 150 DO 140 J = 1, M L = K - J EIGVEC(L,I) = EIGVEC(L,I) - ATEMP(J+1,L)*EIGVEC(K,I) 140 CONTINUE 150 CONTINUE 160 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE DLABAX(N, NBAND, A, X, Y) C C THIS SUBROUTINE SETS Y = A*X C WHERE X AND Y ARE VECTORS OF LENGTH N C AND A IS AN N X NBAND SYMMETRIC BAND MATRIX C C FORMAL PARAMETERS C INTEGER N, NBAND DOUBLE PRECISION A(NBAND,1), X(1), Y(1) C C LOCAL VARIABLES C INTEGER I, K, L, M DOUBLE PRECISION ZERO(1) C C FUNCTIONS CALLED C INTEGER MIN0 C C SUBROUTINES CALLED C C DCOPY C ZERO(1) = 0.0D0 CALL DCOPY(N, ZERO, 0, Y, 1) DO 20 K = 1, N Y(K) = Y(K) + A(1,K)*X(K) M = MIN0(N-K+1, NBAND) IF(M .LT. 2) GO TO 20 DO 10 I = 2, M L = K + I - 1 Y(L) = Y(L) + A(I,K)*X(K) Y(K) = Y(K) + A(I,K)*X(L) 10 CONTINUE 20 CONTINUE RETURN END