C============================================================================== C The subroutines for double precision laso calculation C propre to the algorithm which find a few eigenvalues C============================================================================== C VERSION 2 DOES NOT USE EISPACK 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,*), VAL(NMVAL,*), WORK(*) 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, TEMP, DDOT, DNRM2, DABS, TARR(1) EXTERNAL DNPPLA, DNWLA, 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 DO 120 I = 1,1000 EPS = 0.5D0*EPS TEMP = 1.0D0 + EPS IF(TEMP.EQ.1.0D0) GO TO 130 120 CONTINUE 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 C END 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(*), VEC(NMVEC,*), P0(N,*), P1(N,*), * P2(N,*), RES(*), TAU(*), OTAU(*), T(NBAND,*), * ALP(NBLOCK,*), BET(NBLOCK,*), BOUND(*), ATEMP(*), * VTEMP(*), D(*), S(MAXJ,*), 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.D30 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.0D30 TMAX = -1.0D30 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------------------------------------------------------------------------------ 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,*), VEC(NMVEC,*), H(NPERM,*), * HV(NPERM,*), P(N,*), Q(N,*), BOUND(*), D(*), 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 C============================================================================== C Common subroutines for double precision laso calculation C============================================================================== 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(*),DY(*),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF (N.LE.0) RETURN IF (DA.EQ.0D0) RETURN IF (INCX.EQ.1.AND.INCY.EQ.1) GOTO 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) GOTO 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 END C------------------------------------------------------------------------------ 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(*),DY(*) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF (N.LE.0) RETURN IF (INCX.EQ.1.AND.INCY.EQ.1)GOTO 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) GOTO 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 C END C------------------------------------------------------------------------------ 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(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0D0 DTEMP = 0D0 IF (N.LE.0) RETURN IF (INCX.EQ.1.AND.INCY.EQ.1)GOTO 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) GOTO 40 DO 30 I = 1, M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF (N.LT.5) GOTO 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 C END 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 (Lower-band-packed Storage Mode) C C FORMAL PARAMETERS C INTEGER N, NBAND DOUBLE PRECISION A(NBAND,*), X(*), Y(*) C C LOCAL VARIABLES C INTEGER I, K, L, M DOUBLE PRECISION ZERO C C FUNCTIONS CALLED C INTEGER MIN0 C C SUBROUTINES CALLED C C DCOPY C ZERO = 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 C END C------------------------------------------------------------------------------ SUBROUTINE DLABCM(N, NBAND, NL, NR, A, EIGVAL, & 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,*), EIGVAL(*), & EIGVEC(LDE,*), ATOL, ARTOL, BOUND(2,*), ATEMP(*), & D(*), VTEMP(*) 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 C END C------------------------------------------------------------------------------ SUBROUTINE DLABFC(N, NBAND, A, SIGMA, NUMBER, LDE, & 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,*), SIGMA, EIGVEC(LDE,*), & ATEMP(LDAD,*), D(LDAD,*), 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 END 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,*), EIGVAL(*), & EIGVEC(LDE,*), BOUND(2,*), ATEMP(*), D(*), VTEMP(*), & 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) GOTO 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, & EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP) END 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,*), 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 C END 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 DATA IURAND /0/ C DO 10 I = 1, N X(I) = DBLE(URAND(IURAND)) - 0.5D0 10 CONTINUE C END C------------------------------------------------------------------------------ SUBROUTINE DMVPC(NBLOCK, BET, MAXJ, J, S, NUMBER, RESNRM, & ORTHCF, RV) C INTEGER NBLOCK, MAXJ, J, NUMBER DOUBLE PRECISION BET(NBLOCK,*), S(MAXJ,*), RESNRM(*), & ORTHCF(*), RV(*) 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 C END C------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER I, INCX, J, N, NEXT, NN DOUBLE PRECISION DX(*), 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 END C------------------------------------------------------------------------------ SUBROUTINE DORTQR(NZ, N, NBLOCK, Z, B) C INTEGER NZ, N, NBLOCK DOUBLE PRECISION Z(NZ,*), B(NBLOCK,*) 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 END C------------------------------------------------------------------------------ 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 C END C------------------------------------------------------------------------------ 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 C END C------------------------------------------------------------------------------ SUBROUTINE DVSORT(NUM, VAL, RES, IFLAG, V, NMVEC, N, VEC) INTEGER NUM, IFLAG, NMVEC, N DOUBLE PRECISION VAL(*), RES(*), V(*), VEC(NMVEC,*) 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 C------------------------------------------------------------------------------ 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 C END C------------------------------------------------------------------------------