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