C************************************************************************* C * C SPARSE SVD(A) VIA EIGENSYSTEM OF A'A SYMMETRIC MATRIX * C * C * C************************************************************************* C * C (C) COPYRIGHT 1991 * C MICHAEL W. BERRY * C ALL RIGHTS RESERVED * C * C PERMISSION TO COPY ALL OR PART OF ANY OF THIS * C SOFTWARE IS ONLY GRANTED UPON APPROVAL FROM THE * C AUTHOR: MICHAEL W. BERRY, DEPT. OF COMPUTER * C SCIENCE, UNIVERSITY OF TENNESSEE, 107 AYRES HALL, * C KNOXVILLE, TN, 37996-1301 (BERRY@CS.UTK.EDU) * C * C * C * C************************************************************************* C PROGRAM LAS2 C C.... THIS SAMPLE PROGRAM USES LANDR TO COMPUTE SINGULAR TRIPLETS OF A VIA C.... THE EQUIVALENT SYMMETRIC EIGENVALUE PROBLEM C.... C.... B X = LAMBDA X, WHERE X' = (U',V'), LAMBDA = SIGMA**2, C.... C.... C.... B = A'A , AND A IS M (NROW) BY N (NCOL) (NROW>>NCOL), C.... C.... C.... SO THAT {U,SQRT(LAMBDA),V} IS A SINGULAR TRIPLET OF A. C.... (A' = TRANSPOSE OF A) C.... C.... USER SUPPLIED ROUTINES: OPA,OPB C.... C.... OPA( X,Y) TAKES AN N-VECTOR X AND SHOULD RETURN A*X IN Y. C.... OPB(NCOL,X,Y) TAKES AN N-VECTOR X AND SHOULD RETURN B*X IN Y. C C.... THE COMMON BLOCK "GETPUT" IS USED SOLELY BY THE USER-SUPPLIED C.... SUBROUTINE STORE AND HOLDS THE LANCZOS VECTORS. C C.... USER SHOULD EDIT THE FUNCTION TIMER WITH AN APPROPRIATE CALL TO C.... AN INSTRINSIC TIMING ROUTINE THAT RETURNS ELAPSED USER TIME C.... REAL*8 A REAL*4 T0,TIME COMMON/GETPUT/ A(500000) C C.... SEE WRITE-UP FOR DETAILS ON MEMORY ALLOCATION FOR W AND RITZ. C.... N <= LMTN, MAXPRS <= LMTPRS, LANMAX <= LMTLAN. C,T0,TIME INTEGER LMTN,LMTPRS,LMTLAN,LMTNW PARAMETER (LMTN = 100,LMTPRS = 50,LMTLAN = LMTPRS) PARAMETER (LMTNW = 600000) C C.....NOTE: PARAMETER (LMTNW = 6*LMTN+4*LMTLAN+1+LMTLAN*LMTLAN) C INTEGER I,EV,NW,N,LANMAX,MAXPRS,J,NEIG,IERR,MSGLVL,NSIG REAL*8 ENDL,ENDR,KAPPA,W(LMTNW),RITZ(LMTN),BND(LMTN) REAL*8 D(LMTN),TMP,TMP1,XNORM,DDOT LOGICAL VECTORS C COMMON /COUNT/ MCOUNT C C.... COMMON BLOCK FOR S-VECTORS (SHOULD BE CONTAINED IN RITVEC ALSO) C C.... XV1, XV2 ARE TEMP ARRAYS NEEDED COMPUTING FOR SINGULAR VECTORS C.... (IF VECTORS.EQ.TRUE) C.... LEADING DIMENSIONS OF XV1 MUST BE .GE. (NROW+NCOL) C.... LEADING DIMENSIONS OF XV2 MUST BE THE ORDER OF A'A (NCOL) C.... NSIG IS THE NUMBER OF ACCEPTED RITZ VALUES BASED ON KAPPA C.... (RELATIVE ACCURACY) C REAL*8 XV1,XV2 COMMON /VECTORS/ XV1(500,500),XV2(500),NSIG C C.... THIS PROGRAM USES THE HARWELL-BOEING SPARSE MATRIX C FORMAT FOR THE SPARSE MATRIX A C C NMAX= (BOUND ON (NUMBER OF COLUMNS IN A) + 1 C NZMAX= (BOUND ON NUMBER OF NONZEROS IN A) C C.... COMMON FOR HARWELL-BOEING FORMAT (SHOULD BE CONTAINED IN OPA, OPB, AND C.... RITVEC ALSO) C PARAMETER(NMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NMAX), ROWIND(NZMAX) COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND C C.... INTEGER TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD INTEGER NROW, NCOL, NNZERO, NRHS C CHARACTER TITLE*72, KEY*8, TYPE*3, PTRFMT*16, INDFMT*16, 1 VALFMT*20, RHSFMT*20 C CHARACTER*40 NAME C INTEGER NCELLS,MCOUNT,INPUT,MATRIX,RESULTS C COMMON /OUTPUT/ RESULTS C NCELLS(N,LANMAX,I) = 6*N+4*LANMAX+1+I*LANMAX*LANMAX INPUT=1 MATRIX=7 RESULTS=2 MCOUNT=0 C C OPEN FILES FOR INPUT/OUTPUT C OPEN (INPUT,FILE='LAI1') OPEN (MATRIX,FILE='LAI7') OPEN (RESULTS,FILE='LAO2') C C.... READ (MATRIX,10) TITLE, KEY, TOTCRD, PTRCRD, C....1 INDCRD, VALCRD, RHSCRD, C C.....FOR DATSETS MISSING TOTCRD THROUGH RHSCRD FIELDS C READ (MATRIX,11) TITLE, KEY, 1 TYPE, NROW, NCOL, NNZERO, NRHS, PTRFMT, INDFMT, VALFMT, RHSFMT C 10 FORMAT (A72, A8 / 5I14 / A3, 11X, 4I14 / 2A16, 2A20) 11 FORMAT (A72, A8 // A3, 11X, 4I14 / 2A16, 2A20) C C LEAVE IF MATRIX TOO BIG ---- C IF (NCOL .GE. NMAX .OR. NNZERO .GT. NZMAX) THEN WRITE(RESULTS,*) ' SORRY, YOUR MATRIX IS TOO BIG ' STOP ENDIF C C READ DATA... C READ (MATRIX,PTRFMT) (POINTR (I), I = 1, NCOL+1) READ (MATRIX,INDFMT) (ROWIND (I), I = 1, NNZERO) READ (MATRIX,VALFMT) (VALUE (I), I = 1, NNZERO) C C DEFINE LAST ELEMENT OF POINTR IN CASE IT IS NOT. C POINTR(NCOL+1) = NNZERO + 1 C N=NCOL C C------------------------- C C.... NW IS THE SIZE OF THE WORKING ARRAY W C NW = LMTNW C READ(INPUT,*) NAME,LANMAX,MAXPRS,ENDL,ENDR,VECTORS,KAPPA,MSGLVL C WRITE(RESULTS,2000) N,LANMAX,MAXPRS,ENDL,ENDR,VECTORS,KAPPA, * TITLE,NAME,NROW,NCOL,N IF(VECTORS) THEN EV=3 ELSE EV=0 ENDIF IF(EV.GT.0) OPEN(EV,FILE='LAO3',FORM='UNFORMATTED') C C.... DATA VALIDATION C IF (ENDL.GE.ENDR) THEN WRITE(RESULTS,*)'ENDL MUST BE STRICTLY LESS THAN ENDR' STOP ELSE IF (MAXPRS.GT.LANMAX) THEN WRITE(RESULTS,*)'MAXPRS CANNOT EXCEED LANMAX ', * '(MAXPRS =',MAXPRS,' LANMAX =',LANMAX,')' STOP ELSE IF (EV.LE.0.AND.NCELLS(N,LANMAX,0).GT.NW) THEN WRITE(RESULTS,*)'6*N+4*LANMAX+1 CANNOT EXCEED ', * NW,'(N =',N,' LANMAX =',LANMAX,')' STOP ELSE IF (EV.GT.0.AND.NCELLS(N,LANMAX,1).GT.NW) THEN WRITE(RESULTS,*)'6*N+4*LANMAX+1+LANMAX*LANMAX CANNOT EXCEED ', * NW,'(N =',N,' LANMAX =',LANMAX,')' STOP ENDIF C C.... TO GET A RANDOM STARTING VECTOR, THE FIRST N CELLS C.... MUST BE INITIALIZED TO ZERO. C CALL DSCAL(N,0.0D0,W,1) C C.... MAKE A LANCZOS RUN; SEE LANDR FOR MEANINGS OF PARAMETERS. C T0=TIMER() CALL LANDR(N,LANMAX,MAXPRS,ENDL,ENDR,EV,KAPPA, C CALL LANDR(...,J,NEIG,RITZ,BND,W,NW,IERR,MSGLVL) * J,NEIG,RITZ,BND,W,NW,IERR,MSGLVL) TIME=TIMER()-T0 IF (IERR.NE.0) WRITE(RESULTS,2999)IERR C C.... IF YOU PERMUTE EIGENVALUES, DON'T FORGET TO PERMUTE THE C.... EIGENVECTORS WITH THEM. C WRITE(RESULTS,9998) TIME WRITE(RESULTS,9999)J,NEIG,(I,RITZ(I),BND(I),I = 1,J) C C C 9998 FORMAT(/1X,'...... LANSO EXECUTION TIME=',1PE10.2) 9994 FORMAT(/1X,'...... LASVD EXECUTION TIME=',1PE10.2) 9995 FORMAT(1X,'...... ' * /1X,'...... ',6X,' MXV =',I12,3X,' NSIG =',I3 * /1X,'...... ' * /1X,'...... ',3X,3X,' COMPUTED S-VALUES ',2X, * '(','ERROR BNDS',')' * /1X,'...... ' * /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')')) 9997 FORMAT(1X,'...... ' * /1X,'...... ',6X,' MXV =',I12,3X,' NSIG =',I3 * /1X,'...... ' * /1X,'...... ',3X,3X,' COMPUTED S-VALUES ',2X, * '(','RES. NORMS',')' * /1X,'...... ' * /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')')) 9999 FORMAT(1X,'...... ' * /1X,'...... ',6X,' J =',I3,3X,' NEIG =',I3 * /1X,'...... ' * /1X,'...... ',3X,3X,' COMPUTED RITZ VALUES',2X, * '(','ERROR BNDS',')' * /1X,'...... ' * /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')')) 2000 FORMAT( * 1X,'... ' * /1X,'... SOLVE THE [A^TA] EIGENPROBLEM' * /1X,'... NO. OF EQUATIONS =',I5 * /1X,'... MAX. NO. OF LANCZOS STEPS =',I4 * /1X,'... MAX. NO. OF EIGENPAIRS =',I5 * /1X,'... LEFT END OF THE INTERVAL =',1PE22.14 * /1X,'... RIGHT END OF THE INTERVAL =',1PE22.14 * /1X,'... WANT S-VECTORS? [T/F] =',L4 * /1X,'... KAPPA =',1PE22.14 * //1X,A50 * /1X,A50 * /1X,'... NO. OF DOCUMENTS (ROWS) = ',I8 * /1X,'... NO. OF TERMS (COLS) = ',I8 * /1X,'... ORDER OF MATRIX A = ',I8 * /1X,'... ') 2999 FORMAT(1X,'... RETURN FLAG =',I9,4(' '),'...') C C.... COMPUTE RESIDUAL ERROR WHEN SINGULAR VALUES AND VECTORS ARE COMPUTED C IF(VECTORS) THEN C T0=TIMER() DO 35 JJ = 1,NSIG C C...........MULTIPLY BY MATRIX B FIRST C CALL OPB(N,XV1(1,JJ),XV2) TMP =DDOT(N,XV1(1,JJ),1,XV2,1) CALL DAXPY(N,-TMP,XV1(1,JJ),1,XV2,1) TMP =DSQRT(TMP) XNORM=DSQRT(DDOT(N,XV2,1,XV2,1)) C C...........MULTIPLY BY MATRIX A TO GET (SCALED) LEFT S-VECTOR C CALL OPA(XV1(1,JJ),XV1(NCOL+1,JJ)) TMP1=1.0D0/TMP CALL DSCAL(NROW,TMP1,XV1(NCOL+1,JJ),1) XNORM=XNORM*TMP1 BND(JJ)=XNORM D(JJ)=TMP C C...........WRITE LEFT S-VECTOR TO UNFORMATTED OUTPUT FILE C WRITE(EV) (XV1(NCOL+I,JJ),I=1,NROW) 35 CONTINUE C TIME=TIME+TIMER()-T0 WRITE(RESULTS,9994) TIME WRITE(RESULTS,9997) MCOUNT,NSIG,(I,D(I),BND(I),I = 1,NSIG) ELSE NSIG=0 DO 45 I=1,J K=J-I+1 IF (BND(K).LE.KAPPA*DABS(RITZ(K))) THEN NSIG=NSIG+1 ELSE GO TO 46 ENDIF 45 CONTINUE 46 WRITE(RESULTS,9994) TIME WRITE(RESULTS,9995) MCOUNT,NSIG, * (I,DSQRT(RITZ(J-NSIG+I)),BND(J-NSIG+I),I = 1,NSIG) ENDIF C STOP END C SUBROUTINE OPM(N,A,B) INTEGER N REAL*8 A(1),B(1) INTEGER I C DO 100 I = 1,N B(I) = A(I) 100 CONTINUE RETURN END C SUBROUTINE OPB(N, X, Y) C C.....MULTIPLICATION OF MATRIX B BY A VECTOR X , WHERE C C B = A'A, WHERE A IS NROW BY NCOL (NROW >> NCOL) C C HENCE, B IS OF ORDER N:=NCOL (Y STORES PRODUCT VECTOR) C PARAMETER(NMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NMAX), ROWIND(NZMAX) COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND C COMMON /COUNT/ MCOUNT INTEGER N,MCOUNT REAL*8 X(1), Y(1), ZTEMP(NMAX) C C--------------------- C MCOUNT=MCOUNT+2 C DO 55 I=1,NCOL Y(I)=0.0 55 CONTINUE C DO 85 I=1,NROW ZTEMP(I)=0.0 85 CONTINUE C C.....MULTIPLY BY SPARSE C C DO 15 I=1,NCOL C DO 10 J=POINTR(I),POINTR(I+1)-1 ZTEMP(ROWIND(J))=ZTEMP(ROWIND(J))+VALUE(J)*X(I) 10 CONTINUE C 15 CONTINUE C C.....MULTIPLY BY SPARSE C' C DO 25 I =1,NCOL C DO 20 J=POINTR(I),POINTR(I+1)-1 Y(I)=Y(I)+VALUE(J)*ZTEMP(ROWIND(J)) 20 CONTINUE 25 CONTINUE C RETURN END C SUBROUTINE OPA (X, Y) C C.....MULTIPLICATION OF MATRIX A BY VECTOR X , WHERE C C A IS NROW BY NCOL (NROW >> NCOL) C C (Y STORES PRODUCT VECTOR) C PARAMETER(NMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NMAX), ROWIND(NZMAX) COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND C COMMON /COUNT/ MCOUNT INTEGER MCOUNT REAL*8 X(1), Y(1) C C--------------------- C MCOUNT=MCOUNT+1 C DO 55 I=1,NROW Y(I)=0.0 55 CONTINUE C DO 15 I=1,NCOL C DO 10 J=POINTR(I),POINTR(I+1)-1 Y(ROWIND(J))=Y(ROWIND(J))+VALUE(J)*X(I) 10 CONTINUE C 15 CONTINUE C RETURN END C C.... USER SUPPLIED ELAPSED CPU TIME ROUTINE C REAL*4 FUNCTION TIMER() REAL*4 TMP(2), ETIME TIMER=ETIME(TMP(1)) RETURN END C C @(#)LANDR.F 3.11 (BNP) 4/29/89; FROM LANDR.F 2.16 6/7/88 C SUBROUTINE LANDR(N,LANMAX,MAXPRS,ENDL,ENDR,EV,KAPPA, C SUBROUTINE LANDR(...,J,NEIG,RITZ,BND,W,NW,IERR,MSGLVL) * J,NEIG,RITZ,BND,W,NW,IERR,MSGLVL) INTEGER N,LANMAX,MAXPRS,EV,J,NEIG,NW,IERR,MSGLVL REAL*8 ENDL,ENDR,KAPPA,RITZ(LANMAX),BND(LANMAX),W(NW) C C.... THE PROGRAM MAKE A LANCZOS RUN USING A LINEAR OPERATOR THAT ACTS C.... THROUGH A USER SUPPLIED SUBROUTINE CALLED OPB IN ORDER TO COMPUTE C.... EITHER ALL EIGENVALUES OUTSIDE AN INTERVAL [ENDL,ENDR] OR THE C.... FIRST MAXPRS EIGENVALUES, WHICHEVER OCCURS FIRST. THE INNER C.... PRODUCT IS DEFINED IMPLICITLY BY A USER SUPPLIED SUBROUTINE OPM. C C ***************************************** C * * C * LANCZOS ALORITHM WITH * C * SELECTIVE ORTHOGONALIZATION * C * L A N S O * C * (USING SIMON'S RECURRENCE) * C * * C ***************************************** C C.... INPUTS C.... N DIMENSION OF THE EIGENPROBLEM C.... LANMAX UPPER LIMIT TO THE NUMBER OF LANCZOS STEPS C.... MAXPRS UPPER LIMIT TO THE NUMBER OF WANTED EIGENPAIRS C.... ENDL LEFT END OF THE INTERVAL CONTAINING THE UNWANTED EIGENVALUES C.... ENDR RIGHT END OF THE INTERVAL CONTAINING THE UNWANTED EIGENVALUES C.... EV .LE.0 MEANS EIGENVALUES ONLY, C.... .GT.0 MEANS BOTH EIGENVALUES AND EIGENVECTORS ARE WANTED C.... AND EV BECOMES THE OUTPUT CHANNEL FOR EIGENVECTORS C.... KAPPA RELATIVE ACCURACY OF RITZ VALUES ACCEPTABLE AS EIGENVALUES C.... OF THE MATRIX A'A C.... NW LENGTH OF THE WORK ARRAY W C.... W WORK ARRAY OF LENGTH NW C C.... OUTPUTS C.... J NUMBER OF LANCZOS STEPS ACTUALLY TAKEN C.... NEIG NUMBER OF RITZ VALUES STABILIZED C.... RITZ ARRAY TO HOLD THE RITZ VALUES C.... BND ARRAY TO HOLD THE ERROR BOUNDS C.... IERR ERROR FLAG C C.... SUBROUTINES: MACHAR,LANSO,RITVEC C REAL*8 RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 COMMON/RDATA/RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 C INTEGER MT,I,NQ(4),N1,L2,L3,L4,L5,LSQ6 C C.... MACHAR SPECIFIC (EPS DECLARED IN COMMON/RDATA/) C INTEGER IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP REAL*8 EPSNEG,XMIN,XMAX C C.... CHECK INPUT DATA C MT = 6*N+1+4*LANMAX IF (EV.GT.0) MT = MT+LANMAX*LANMAX IERR = 0 IF (N.LE.0) IERR = IERR+1 IF (LANMAX.LE.0) IERR = IERR+2 IF (ENDR.LE.ENDL) IERR = IERR+4 IF (MAXPRS.LE.0) IERR = IERR+8 IF (MAXPRS.GT.LANMAX) IERR = IERR+16 IF (LANMAX.GT.N) IERR = IERR+32 IF (MT.GT.NW) IERR = IERR+64 IF (IERR.GT.0) RETURN C C.... COMPUTE MACHINE PRECISION C CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, * EPS,EPSNEG,XMIN,XMAX) EPS1 = EPS*DSQRT(DBLE(N)) REPS = DSQRT(EPS) EPS34 = REPS*DSQRT(REPS) C C.... SET POINTERS AND INITIALIZE C NQ(1) = N+1 DO 20 I = 2,4 NQ(I) = NQ(I-1)+N 20 CONTINUE N1 = 1+5*N L2 = N+N1 L3 = LANMAX+L2 L4 = 1+LANMAX+L3 L5 = LANMAX+L4 CALL LANSO(N,LANMAX,MAXPRS,ENDL,ENDR,J,NEIG,RITZ,BND, C CALL LANSO(...,W,W(N1),W(L2),W(L3),W(L4),W(L5),NQ,IERR,MSGLVL) * W,W(N1),W(L2),W(L3),W(L4),W(L5),NQ,IERR,MSGLVL) C C.... COMPUTE EIGENVECTORS C IF (EV.GT.0) THEN LSQ6 = LANMAX+L5 C C.... RATIONALIZE KAPPA C KAPPA = DMAX1(DABS(KAPPA),EPS34) C C.... W(NQ(4)) BECOMES WORKING STORAGE FROM THIS POINT ON C CALL RITVEC(N,J,EV,KAPPA,RNM,RITZ,BND,W(L2),W(L3), C CALL RITVEC(...,W(LSQ6),W(NQ(4)),W(N1),IERR,MSGLVL) * W(LSQ6),W(NQ(4)),W(N1),IERR,MSGLVL) ENDIF RETURN END C C @(#)LANSO.F 3.14 (BNP) 4/29/89; FROM LANSO.F 2.18 7/7/88 C SUBROUTINE LANSO(N,LANMAX,MAXPRS,ENDL,ENDR,J,NEIG,RITZ,BND, C SUBROUTINE LANSO(...,R,WRK,ALF,BET,ETA,OLDETA,NQ,IERR,MSGLVL) * R,WRK,ALF,BET,ETA,OLDETA,NQ,IERR,MSGLVL) INTEGER N,LANMAX,MAXPRS,J,NEIG,NQ(4),IERR,MSGLVL INTEGER MCOUNT REAL*8 ENDL,ENDR,R(5*N),WRK(N), * ALF(LANMAX),BET(LANMAX+1),RITZ(LANMAX),BND(LANMAX), * ETA(LANMAX),OLDETA(LANMAX) C COMMON /MXV/ MCOUNT C C.... INPUTS C.... N DIMENSION OF THE EIGENPROBLEM C.... LANMAX UPPER LIMIT TO THE NUMBER OF LANCZOS STEPS C.... MAXPRS UPPER LIMIT TO THE NUMBER OF WANTED EIGENPAIRS C.... ENDL LEFT END OF THE INTERVAL CONTAINING THE WANTED EIGENVALUES C.... ENDR RIGHT END OF THE INTERVAL CONTAINING THE WANTED EIGENVALUES C C.... WORK SPACE C.... R HOLDS 5 VECTORS OF LENGTH N. SEE THE TEXT FOR DETAILS. C.... NQ(4) CONTAINS THE POINTERS TO THE BEGINING OF EACH VECTOR IN R. C.... ALF ARRAY TO HOLD DIAGONAL OF THE TRIDIAGONAL T C.... BET ARRAY TO HOLD OFF-DIAGONAL OF T C.... ETA ORTHOGONALITY ESTIMATE OF LANCZOS VECTORS AT STEP J C.... OLDETA ORTHOGONALITY ESTIMATE OF LANCZOS VECTORS AT STEP J-1 C C.... OUTPUTS C.... J ACTUAL NUMBER OF LANCZOS STEPS TAKEN C.... NEIG NUMBER OF COMPUTED EIGENPAIRS C.... RITZ ARRAY TO HOLD THE CONVERGED RITZ VALUES C.... BND ARRAY TO HOLD THE ERROR BOUNDS C.... IERR ERROR FLAG C C.... BLAS ROUTINES: DATX,DAXPY,DCOPY,DDOT,DSCAL,IDAMAX C.... SUBROUTINES: DSORT2,IMTQLB,ORTBND,PURGE,STARTV,STPONE C.... USER-SUPPLIED: OPB,OPM,STORE C INTEGER STORQ,RETRQ,STORP,RETRP PARAMETER (STORQ = 1,RETRQ = 2,STORP = 3,RETRP = 4) C INTEGER MAXLL REAL*8 FOUR PARAMETER (MAXLL = 2,FOUR = 4.0D0) C REAL*8 RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 COMMON/RDATA/RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 C LOGICAL ENOUGH INTEGER LL,I,L,FIRST,LAST,MID,ID1,ID2,ID3,IDAMAX REAL*8 T,GAPL,GAP,DDOT,STARTV C REAL*8 ONE,ZERO DATA ONE,ZERO/1.0D0,0.0D0/ C CALL STPONE(N,R,WRK,ALF,NQ,IERR,MSGLVL) IF (RNM.EQ.ZERO.OR.IERR.NE.0) RETURN C ETA(1) = EPS1 OLDETA(1) = EPS1 C LL = 0 FIRST = 2 LAST = MIN(MAXPRS+MAX(8,MAXPRS),LANMAX) ENOUGH = .FALSE. DO 100 ID1 = 1,MAXPRS IF (ENOUGH) GOTO 200 IF (RNM.LE.TOL) RNM = ZERO C C.... LANCZOS LOOP C DO 10 J = FIRST,LAST MID = NQ(2) NQ(2) = NQ(1) NQ(1) = MID MID = NQ(3) NQ(3) = NQ(4) NQ(4) = MID CALL STORE(N,STORQ,J-1,R(NQ(2))) IF (J-1.LE.MAXLL) CALL STORE(N,STORP,J-1,R(NQ(4))) BET(J) = RNM C C... RESTART IF INVARIANT SUBSPACE IS FOUND C IF (BET(J).EQ.ZERO) THEN RNM = STARTV(N,J,R,WRK,NQ,EPS,IERR,MSGLVL) IF (IERR.NE.0) RETURN ENOUGH = RNM.EQ.ZERO IF (ENOUGH) GOTO 15 ENDIF C C.... TAKE A LANCZOS STEP C T = ONE/RNM CALL DATX(N,T,R,1,R(NQ(1)),1) CALL DSCAL(N,T,R(NQ(3)),1) C CALL OPB(N,R(NQ(3)),R) CALL DAXPY(N,-RNM,R(NQ(2)),1,R,1) C ALF(J) = DDOT(N,R,1,R(NQ(3)),1) CALL DAXPY(N,-ALF(J),R(NQ(1)),1,R,1) C C.... ORTHOGONALIZE AGAINST INITIAL LANCZOS VECTORS C IF (J.LE.MAXLL+1.AND.DABS(ALF(J-1)).GT.FOUR*DABS(ALF(J))) THEN LL = J-1 ENDIF DO 5 I = 1,MIN(LL,J-2) CALL STORE(N,RETRP,I,WRK) T = DDOT(N,WRK,1,R,1) CALL STORE(N,RETRQ,I,WRK) CALL DAXPY(N,-T,WRK,1,R,1) ETA(I) = EPS1 OLDETA(I) = EPS1 5 CONTINUE C C.... EXTENDED LOCAL REORTHOGONALIZATION C T = DDOT(N,R,1,R(NQ(4)),1) CALL DAXPY(N,-T,R(NQ(2)),1,R,1) IF (BET(J).GT.ZERO) BET(J) = BET(J)+T T = DDOT(N,R,1,R(NQ(3)),1) CALL DAXPY(N,-T,R(NQ(1)),1,R,1) ALF(J) = ALF(J)+T CALL OPM(N,R,R(NQ(4))) RNM = DSQRT(DDOT(N,R,1,R(NQ(4)),1)) ANORM = BET(J)+DABS(ALF(J))+RNM TOL = REPS*ANORM C C.... UPDATE THE ORTHOGONALITY BOUNDS C CALL ORTBND(N,J,ALF,BET,ETA,OLDETA,RNM,EPS1) C C.... RESTORE THE ORTHOGONALITY STATE WHEN NEEDED C CALL PURGE(N,LL,J,R,R(NQ(1)),R(NQ(4)),R(NQ(3)),WRK, C CALL PURGE(...,ETA,OLDETA,MSGLVL) * ETA,OLDETA,MSGLVL) IF (RNM.LE.TOL) RNM = ZERO 10 CONTINUE J = LAST 15 IF (ENOUGH) J = J-1 FIRST = J+1 BET(J+1) = RNM C C.... NOW ANALYZE T C L = 1 DO 40 ID2 = 1,J IF (L.GT.J) GOTO 50 DO 20 I = L,J IF (BET(I+1).EQ.ZERO) GOTO 30 20 CONTINUE I = J C.... NOW I IS AT THE END OF AN UNREDUCED SUBMATRIX C 30 CALL DCOPY(I-L+1,ALF(L),1,RITZ(L),-1) CALL DCOPY(I-L,BET(L+1),1,WRK(L+1),-1) CALL IMTQLB(I-L+1,RITZ(L),WRK(L),BND(L),IERR) IF (IERR.NE.0) THEN PRINT *,' IMTQLB FAILED TO CONVERGE (IERR =',IERR,')' PRINT *,' L =',L,' I =',I PRINT *,(ID3,RITZ(ID3),WRK(ID3),BND(ID3),ID3 = L,I) ENDIF DO 35 ID3 = L,I BND(ID3) = RNM*DABS(BND(ID3)) 35 CONTINUE L = I+1 40 CONTINUE C C.... SORT EIGENVALUES INTO INCREASING ORDER C 50 CALL DSORT2(J,RITZ,BND) C C.... MASSAGE ERROR BOUNDS FOR VERY CLOSE RITZ VALUES C MID = IDAMAX(J,BND,1) DO 70 L = -1,1,2 DO 60 I = ((J+1)-L*(J-1))/2,MID-L,L IF (DABS(RITZ(I+L)-RITZ(I)).LT.EPS34*DABS(RITZ(I))) THEN IF (BND(I).GT.TOL.AND.BND(I+L).GT.TOL) THEN BND(I+L) = DSQRT(BND(I)**2+BND(I+L)**2) BND(I) = ZERO ENDIF ENDIF 60 CONTINUE 70 CONTINUE C C.... REFINE THE ERROR BOUNDS C NEIG = 0 GAPL = RITZ(J)-RITZ(1) DO 80 I = 1,J GAP = GAPL IF (I.LT.J) GAPL = RITZ(I+1)-RITZ(I) GAP = MIN(GAP,GAPL) IF (GAP.GT.BND(I)) THEN BND(I) = BND(I)*(BND(I)/GAP) ENDIF IF (BND(I).LE.16.0D0*EPS*DABS(RITZ(I))) THEN NEIG = NEIG+1 C ENOUGH = ENOUGH.OR.ENDL.LT.RITZ(I).AND.RITZ(I).LT.ENDR IF(.NOT.(ENOUGH) ) THEN ENOUGH = ENDL.LT.RITZ(I).AND.RITZ(I).LT.ENDR ENDIF ENDIF 80 CONTINUE C C.... SHOULD WE STOP? C IF (NEIG.LT.MAXPRS) THEN IF (NEIG.EQ.0) THEN LAST = FIRST+8 ELSE LAST = FIRST+MAX(2,((J-6)*(MAXPRS-NEIG))/NEIG) ENDIF LAST = MIN(LAST,LANMAX) ELSE ENOUGH = .TRUE. ENDIF ENOUGH = ENOUGH.OR.FIRST.GT.LANMAX 100 CONTINUE 200 CALL STORE(N,STORQ,J,R(NQ(1))) RETURN END SUBROUTINE DATX(N,DA,DX,INCX,DY,INCY) C C DY := DA*DX C REAL*8 DX(*),DY(*),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N IF (N.LE.0) RETURN IF (DA.EQ.0.0D0) RETURN IF (INCX.EQ.1.AND.INCY.EQ.1) GO TO 20 C C UNEQUAL INCREMENTS OR EQUAL INCREMENTS .NE. ONE 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) = DA*DX(IX) IX = IX+INCX IY = IY+INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 M = MOD(N,4) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M DY(I) = DA*DX(I) 30 CONTINUE IF (N.LT.4) RETURN 40 MP1 = M+1 DO 50 I = MP1,N,4 DY(I) = DA*DX(I) DY(I+1) = DA*DX(I+1) DY(I+2) = DA*DX(I+2) DY(I+3) = DA*DX(I+3) 50 CONTINUE RETURN END C C @(#)DSORT2.F 3.2 (BNP) 12/9/88 C SUBROUTINE DSORT2(N,ARRAY1,ARRAY2) INTEGER N REAL*8 ARRAY1(0:N-1),ARRAY2(0:N-1) C C.... SORT ARRAY1 AND ARRAY2 INTO INCREASING ORDER FOR ARRAY1 C INTEGER IGAP,I,J REAL*8 TEMP C IGAP = N/2 10 IF (IGAP.GT.0) THEN DO 200 I = IGAP,N-1 J = I-IGAP 50 IF (J.GE.0.AND.ARRAY1(J).GT.ARRAY1(J+IGAP)) THEN TEMP = ARRAY1(J) ARRAY1(J) = ARRAY1(J+IGAP) ARRAY1(J+IGAP) = TEMP TEMP = ARRAY2(J) ARRAY2(J) = ARRAY2(J+IGAP) ARRAY2(J+IGAP) = TEMP ELSE GO TO 200 ENDIF J = J-IGAP GO TO 50 200 CONTINUE ELSE RETURN ENDIF IGAP = IGAP/2 GO TO 10 END SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,NM,MML,IERR REAL*8 D(N),E(N),Z(NM,N) REAL*8 B,C,F,G,P,R,S,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1. C C E HAS BEEN DESTROYED. C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0D0 C DO 240 L = 1, N J = 0 C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 TST1 = DABS(D(M)) + DABS(D(M+1)) TST2 = TST1 + DABS(E(M)) IF (TST2 .EQ. TST1) GO TO 120 110 CONTINUE C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... G = (D(L+1) - P) / (2.0D0 * E(L)) R = PYTHAG(G,1.0D0) G = D(M) - P + E(L) / (G + DSIGN(R,G)) S = 1.0D0 C = 1.0D0 P = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML I = M - II F = S * E(I) B = C * E(I) R = PYTHAG(F,G) E(I+1) = R IF (R .EQ. 0.0D0) GO TO 210 S = F / R C = G / R G = D(I+1) - P R = (D(I) - G) * S + 2.0D0 * C * B P = S * R D(I+1) = G + P G = C * R - B C .......... FORM VECTOR .......... DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0D0 GO TO 105 C .......... RECOVER FROM UNDERFLOW .......... 210 D(I+1) = D(I+1) - P E(M) = 0.0D0 GO TO 105 240 CONTINUE C .......... ORDER EIGENVALUES AND EIGENVECTORS .......... DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END C C @(#)IMTQLB.F 3.5 (BNP) 4/28/89 C SUBROUTINE IMTQLB(N,D,E,BND,IERR) C INTEGER I,J,L,M,N,II,MML,IERR REAL*8 D(N),E(N),BND(N) REAL*8 B,C,F,G,P,R,S,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A MODIFICATION OF THE ALGOL PROCEDURE IMTQL1, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E HAS BEEN DESTROYED. C C BND WILL HOLD THE TOP ELEMENTS OF THE NORMALIZED EIGENVECTORS. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C THIS VERSION DATED AUGUST 1983 (NOVEMBER 1988). C IERR = 0 BND(1) = 1.0D0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N BND(I) = 0.0D0 100 E(I-1) = E(I) C E(N) = 0.0D0 C DO 290 L = 1, N J = 0 C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 TST1 = DABS(D(M)) + DABS(D(M+1)) TST2 = TST1 + DABS(E(M)) IF (TST2 .EQ. TST1) GO TO 120 110 CONTINUE C 120 P = D(L) F = BND(L) IF (M .EQ. L) GO TO 215 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... G = (D(L+1) - P) / (2.0D0 * E(L)) R = PYTHAG(G,1.0D0) G = D(M) - P + E(L) / (G + DSIGN(R,G)) S = 1.0D0 C = 1.0D0 P = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML I = M - II F = S * E(I) B = C * E(I) R = PYTHAG(F,G) E(I+1) = R IF (R .EQ. 0.0D0) GO TO 210 S = F / R C = G / R G = D(I+1) - P R = (D(I) - G) * S + 2.0D0 * C * B P = S * R D(I+1) = G + P G = C * R - B F = BND(I+1) BND(I+1) = S*BND(I)+C*F BND(I) = C*BND(I)-S*F 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0D0 GO TO 105 C .......... RECOVER FROM UNDERFLOW .......... 210 D(I+1) = D(I+1) - P E(M) = 0.0D0 GO TO 105 C .......... ORDER EIGENVALUES .......... 215 IF (L .EQ. 1) GO TO 250 C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... DO 230 II = 2, L I = L + 2 - II IF (P .GE. D(I-1)) GO TO 270 D(I) = D(I-1) BND(I) = BND(I-1) 230 CONTINUE C 250 I = 1 270 D(I) = P BND(I) = F 290 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) C----------------------------------------------------------------------- C THIS FORTRAN 77 SUBROUTINE IS INTENDED TO DETERMINE THE PARAMETERS C OF THE FLOATING-POINT ARITHMETIC SYSTEM SPECIFIED BELOW. THE C DETERMINATION OF THE FIRST THREE USES AN EXTENSION OF AN ALGORITHM C DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, INCORPORATING SOME, C BUT NOT ALL, OF THE IMPROVEMENTS SUGGESTED BY M. GENTLEMAN AND S. C MAROVICH, CACM 17 (1974), PP. 276-277. AN EARLIER VERSION OF THIS C PROGRAM WAS PUBLISHED IN THE BOOK SOFTWARE MANUAL FOR THE C ELEMENTARY FUNCTIONS BY W. J. CODY AND W. WAITE, PRENTICE-HALL, C ENGLEWOOD CLIFFS, NJ, 1980. THE PRESENT VERSION IS DOCUMENTED IN C W. J. CODY, "MACHAR: A SUBROUTINE TO DYNAMICALLY DETERMINE MACHINE C PARAMETERS," TOMS 14, DECEMBER, 1988. C C PARAMETER VALUES REPORTED ARE AS FOLLOWS: C C IBETA - THE RADIX FOR THE FLOATING-POINT REPRESENTATION C IT - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING-POINT C SIGNIFICAND C IRND - 0 IF FLOATING-POINT ADDITION CHOPS C 1 IF FLOATING-POINT ADDITION ROUNDS, BUT NOT IN THE C IEEE STYLE C 2 IF FLOATING-POINT ADDITION ROUNDS IN THE IEEE STYLE C 3 IF FLOATING-POINT ADDITION CHOPS, AND THERE IS C PARTIAL UNDERFLOW C 4 IF FLOATING-POINT ADDITION ROUNDS, BUT NOT IN THE C IEEE STYLE, AND THERE IS PARTIAL UNDERFLOW C 5 IF FLOATING-POINT ADDITION ROUNDS IN THE IEEE STYLE, C AND THERE IS PARTIAL UNDERFLOW C NGRD - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION WITH C TRUNCATING ARITHMETIC. IT IS C 0 IF FLOATING-POINT ARITHMETIC ROUNDS, OR IF IT C TRUNCATES AND ONLY IT BASE IBETA DIGITS C PARTICIPATE IN THE POST-NORMALIZATION SHIFT OF THE C FLOATING-POINT SIGNIFICAND IN MULTIPLICATION; C 1 IF FLOATING-POINT ARITHMETIC TRUNCATES AND MORE C THAN IT BASE IBETA DIGITS PARTICIPATE IN THE C POST-NORMALIZATION SHIFT OF THE FLOATING-POINT C SIGNIFICAND IN MULTIPLICATION. C MACHEP - THE LARGEST NEGATIVE INTEGER SUCH THAT C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, EXCEPT THAT C MACHEP IS BOUNDED BELOW BY -(IT+3) C NEGEPS - THE LARGEST NEGATIVE INTEGER SUCH THAT C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, EXCEPT THAT C NEGEPS IS BOUNDED BELOW BY -(IT+3) C IEXP - THE NUMBER OF BITS (DECIMAL PLACES IF IBETA = 10) C RESERVED FOR THE REPRESENTATION OF THE EXPONENT C (INCLUDING THE BIAS OR SIGN) OF A FLOATING-POINT C NUMBER C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT C FLOAT(IBETA)**MINEXP IS POSITIVE AND NORMALIZED C MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS C EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH C THAT 1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER C IBETA = 2 OR IRND = 0, EPS = FLOAT(IBETA)**MACHEP. C OTHERWISE, EPS = (FLOAT(IBETA)**MACHEP)/2 C EPSNEG - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT C 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2 C OR IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS. C OTHERWISE, EPSNEG = (IBETA**NEGEPS)/2. BECAUSE C NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT C BE THE SMALLEST NUMBER THAT CAN ALTER 1.0 BY C SUBTRACTION. C XMIN - THE SMALLEST NON-VANISHING NORMALIZED FLOATING-POINT C POWER OF THE RADIX, I.E., XMIN = FLOAT(IBETA)**MINEXP C XMAX - THE LARGEST FINITE FLOATING-POINT NUMBER. IN C PARTICULAR XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP C NOTE - ON SOME MACHINES XMAX WILL BE ONLY THE C SECOND, OR PERHAPS THIRD, LARGEST NUMBER, BEING C TOO SMALL BY 1 OR 2 UNITS IN THE LAST DIGIT OF C THE SIGNIFICAND. C C LATEST REVISION - DECEMBER 4, 1987 C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C----------------------------------------------------------------------- INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP, 1 MINEXP,MX,NEGEP,NGRD,NXRES REAL*8 1 A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA, 2 TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO C----------------------------------------------------------------------- CONV(I) = DBLE(I) ONE = CONV(1) TWO = ONE + ONE ZERO = ONE - ONE C----------------------------------------------------------------------- C DETERMINE IBETA, BETA ALA MALCOLM. C----------------------------------------------------------------------- A = ONE 10 A = A + A TEMP = A+ONE TEMP1 = TEMP-A IF (TEMP1-ONE .EQ. ZERO) GO TO 10 B = ONE 20 B = B + B TEMP = A+B ITEMP = INT(TEMP-A) IF (ITEMP .EQ. 0) GO TO 20 IBETA = ITEMP BETA = CONV(IBETA) C----------------------------------------------------------------------- C DETERMINE IT, IRND. C----------------------------------------------------------------------- IT = 0 B = ONE 100 IT = IT + 1 B = B * BETA TEMP = B+ONE TEMP1 = TEMP-B IF (TEMP1-ONE .EQ. ZERO) GO TO 100 IRND = 0 BETAH = BETA / TWO TEMP = A+BETAH IF (TEMP-A .NE. ZERO) IRND = 1 TEMPA = A + BETA TEMP = TEMPA+BETAH IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 C----------------------------------------------------------------------- C DETERMINE NEGEP, EPSNEG. C----------------------------------------------------------------------- NEGEP = IT + 3 BETAIN = ONE / BETA A = ONE DO 200 I = 1, NEGEP A = A * BETAIN 200 CONTINUE B = A 210 TEMP = ONE-A IF (TEMP-ONE .NE. ZERO) GO TO 220 A = A * BETA NEGEP = NEGEP - 1 GO TO 210 220 NEGEP = -NEGEP EPSNEG = A C----------------------------------------------------------------------- C DETERMINE MACHEP, EPS. C----------------------------------------------------------------------- MACHEP = -IT - 3 A = B 300 TEMP = ONE+A IF (TEMP-ONE .NE. ZERO) GO TO 320 A = A * BETA MACHEP = MACHEP + 1 GO TO 300 320 EPS = A C----------------------------------------------------------------------- C DETERMINE NGRD. C----------------------------------------------------------------------- NGRD = 0 TEMP = ONE+EPS IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1 C----------------------------------------------------------------------- C DETERMINE IEXP, MINEXP, XMIN. C C LOOP TO DETERMINE LARGEST I AND K = 2**I SUCH THAT C (1/BETA) ** (2**(I)) C DOES NOT UNDERFLOW. C EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW. C----------------------------------------------------------------------- I = 0 K = 1 Z = BETAIN T = ONE + EPS NXRES = 0 400 Y = Z Z = Y * Y C----------------------------------------------------------------------- C CHECK FOR UNDERFLOW HERE. C----------------------------------------------------------------------- A = Z * ONE TEMP = Z * T IF ((A+A .EQ. ZERO) .OR. (DABS(Z) .GE. Y)) GO TO 410 TEMP1 = TEMP * BETAIN IF (TEMP1*BETA .EQ. Z) GO TO 410 I = I + 1 K = K + K GO TO 400 410 IF (IBETA .EQ. 10) GO TO 420 IEXP = I + 1 MX = K + K GO TO 450 C----------------------------------------------------------------------- C THIS SEGMENT IS FOR DECIMAL MACHINES ONLY. C----------------------------------------------------------------------- 420 IEXP = 2 IZ = IBETA 430 IF (K .LT. IZ) GO TO 440 IZ = IZ * IBETA IEXP = IEXP + 1 GO TO 430 440 MX = IZ + IZ - 1 C----------------------------------------------------------------------- C LOOP TO DETERMINE MINEXP, XMIN. C EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW. C----------------------------------------------------------------------- 450 XMIN = Y Y = Y * BETAIN C----------------------------------------------------------------------- C CHECK FOR UNDERFLOW HERE. C----------------------------------------------------------------------- A = Y * ONE TEMP = Y * T IF (((A+A) .EQ. ZERO) .OR. (DABS(Y) .GE. XMIN)) GO TO 460 K = K + 1 TEMP1 = TEMP * BETAIN IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN GO TO 450 ELSE NXRES = 3 XMIN = Y END IF 460 MINEXP = -K C----------------------------------------------------------------------- C DETERMINE MAXEXP, XMAX. C----------------------------------------------------------------------- IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500 MX = MX + MX IEXP = IEXP + 1 500 MAXEXP = MX + MINEXP C----------------------------------------------------------------- C ADJUST IRND TO REFLECT PARTIAL UNDERFLOW. C----------------------------------------------------------------- IRND = IRND + NXRES C----------------------------------------------------------------- C ADJUST FOR IEEE-STYLE MACHINES. C----------------------------------------------------------------- IF (IRND .GE. 2) MAXEXP = MAXEXP - 2 C----------------------------------------------------------------- C ADJUST FOR MACHINES WITH IMPLICIT LEADING BIT IN BINARY C SIGNIFICAND, AND MACHINES WITH RADIX POINT AT EXTREME C RIGHT OF SIGNIFICAND. C----------------------------------------------------------------- I = MAXEXP + MINEXP IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1 IF (I .GT. 20) MAXEXP = MAXEXP - 1 IF (A .NE. Y) MAXEXP = MAXEXP - 2 XMAX = ONE - EPSNEG IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG XMAX = XMAX / (BETA * BETA * BETA * XMIN) I = MAXEXP + MINEXP + 3 IF (I .LE. 0) GO TO 520 DO 510 J = 1, I IF (IBETA .EQ. 2) XMAX = XMAX + XMAX IF (IBETA .NE. 2) XMAX = XMAX * BETA 510 CONTINUE 520 RETURN C---------- LAST CARD OF MACHAR ---------- END C C @(#)ORTBND.F 3.4 (BNP) 4/28/89; FROM ORTBND.F 2.7 10/17/87 C SUBROUTINE ORTBND(N,J,ALF,BET,ETA,OLDETA,RNM,EPS1) INTEGER N,J REAL*8 ALF(J),BET(J),ETA(J),OLDETA(J),RNM,EPS1 C C.... UPDATE THE ETA RECURRENCE. C C.... INPUTS C.... N DIMENSION OF THE EIGENPROBLEM C.... J DIMENSION OF T C.... ALF DIAGONAL ELEMENTS OF THE TRIDIAGONAL T C.... BET OFF-DIAGONAL ELEMENTS OF T C.... ETA ORTHOGONALITY ESTIMATE OF LANCZOS VECTORS AT STEP J C.... OLDETA ORTHOGONALITY ESTIMATE OF LANCZOS VECTORS AT STEP J-1 C.... RNM NORM OF THE NEXT RESIDUAL VECTOR C.... EPS1 ROUNDOFF ESTIMATE FOR DOT PRODUCT OF TWO UNIT VECTORS C C.... OUTPUTS C.... ETA ORTHOGONALITY ESTIMATE OF LANCZOS VECTORS AT STEP J+1 C.... OLDETA ORTHOGONALITY ESTIMATE OF LANCZOS VECTORS AT STEP J C C.... BLAS ROUTINES: DSWAP C INTEGER I C IF (J.LE.1) RETURN IF (RNM.EQ.0.0D0) GOTO 200 IF (J.GT.2) THEN OLDETA(1) = (BET(2)*ETA(2)+(ALF(1)-ALF(J))* * ETA(1)-BET(J)*OLDETA(1))/RNM+EPS1 ENDIF DO 100 I = 2,J-2 OLDETA(I) = (BET(I+1)*ETA(I+1)+(ALF(I)-ALF(J))*ETA(I)+ * BET(I)*ETA(I-1)-BET(J)*OLDETA(I))/RNM+EPS1 100 CONTINUE 200 OLDETA(J-1) = EPS1 CALL DSWAP(J-1,OLDETA,1,ETA,1) ETA(J) = EPS1 RETURN END C C @(#)PURGE.F 3.12 (BNP) 4/29/89; FROM PURGE.F 2.13 6/25/88 C SUBROUTINE PURGE(N,LL,J,R,Q,RA,QA,WRK,ETA,OLDETA,MSGLVL) INTEGER N,LL,J,MSGLVL REAL*8 R(N),Q(N),RA(N),QA(N),WRK(N),ETA(J),OLDETA(J) C C.... THIS ROUTINE EXAMINES ETA TO DECIDE WHETHER C.... RE-ORTHOGONALIZATION SHOULD BE PERFORMED. C C.... N DIMENSION OF THE EIGENPROBLEM C.... LL NO. OF INITIAL LANCZOS VECTORS IN LOCAL ORTHOG. C.... J CURRENT LANCZOS STEP C.... R THE RESIDUAL VECTOR TO BECOME THE NEXT LANCZOS VECTOR C.... Q THE CURRENT LANCZOS VECTOR C.... RA THE PRODUCT OF THE MASS MATRIX AND R C.... QA THE PRODUCT OF THE MASS MATRIX AND Q C.... WRK A TEMPORARY VECTOR TO HOLD THE PREVIOUS LANCZOS VECTORS C.... ETA STATE OF ORTHOGONALITY BETWEEN R AND PREVIOUS LANCZOS VECTORS C.... OLDETA STATE OF ORTHOGONALITY BETWEEN Q AND PREVIOUS LANCZOS VECTORS C C.... BLAS ROUTINES: DAXPY,DDOT,IDAMAX C.... SUBROUTINES: NONE C.... USER-SUPPLIED: OPM,STORE C INTEGER RETRQ PARAMETER (RETRQ = 2) C REAL*8 RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 COMMON/RDATA/RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 INTEGER K,I,LOOP,IDAMAX REAL*8 T,TQ,TR,REPS1,DDOT C IF (J.LE.LL+2) RETURN K = IDAMAX(J-(LL+2),ETA(LL+1),1)+LL C IF (DABS(ETA(K)).GT.REPS) THEN REPS1 = EPS1/REPS DO 55 LOOP = 1,2 IF (RNM.GT.TOL) THEN C C.... BRING IN A LANCZOS VECTOR T AND ORTHOGONALIZE BOTH R AND Q C.... AGAINST IT C TQ = 0.0D0 TR = 0.0D0 DO 50 I = LL+1,J-1 CALL STORE(N,RETRQ,I,WRK) T = -DDOT(N,QA,1,WRK,1) TQ = TQ+DABS(T) CALL DAXPY(N,T,WRK,1,Q,1) T = -DDOT(N,RA,1,WRK,1) TR = TR+DABS(T) CALL DAXPY(N,T,WRK,1,R,1) 50 CONTINUE CALL OPM(N,Q,QA) C C.... RESTORE LOCAL ORTHOGONALITY C T = -DDOT(N,R,1,QA,1) TR = TR+DABS(T) CALL DAXPY(N,T,Q,1,R,1) C CALL OPM(N,R,RA) RNM = DSQRT(DDOT(N,RA,1,R,1)) IF (TQ.LE.REPS1.AND.TR.LE.REPS1*RNM) GOTO 58 ENDIF 55 CONTINUE 58 DO 60 I = LL+1,J ETA(I) = EPS1 OLDETA(I) = EPS1 60 CONTINUE ENDIF RETURN END REAL*8 FUNCTION PYTHAG(A,B) REAL*8 A,B C C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C REAL*8 P,R,S,T,U P = DMAX1(DABS(A),DABS(B)) IF (P .EQ. 0.0D0) GO TO 20 R = (DMIN1(DABS(A),DABS(B))/P)**2 10 CONTINUE T = 4.0D0 + R IF (T .EQ. 4.0D0) GO TO 20 S = R/T U = 1.0D0 + 2.0D0*S P = U*P R = (S/U)**2 * R GO TO 10 20 PYTHAG = P RETURN END C C @(#)RANDOM.F 3.2 (BNP) 12/9/88; FROM RANDOM.F 2.2 10/13/87 C REAL*8 FUNCTION RANDOM(IY) INTEGER IY C C RANDOM 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 RANDOM.THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY C BETWEEN SUBSEQUENT CALLS TO RANDOM.VALUES OF RANDOM WILL BE RETURNED C IN THE INTERVAL (0,1). C INTEGER M,IA,IC,MIC REAL*8 HALFM,S C INTEGER M2,ITWO DATA M2,ITWO/0,2/ IF (M2.EQ.0) THEN 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.5D0/HALFM C C COMPUTE NEXT RANDOM NUMBER C ENDIF 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 RANDOM = DBLE(IY)*S RETURN END C C @(#)RITVEC.F 1.3 (BNP) 4/5/89 C SUBROUTINE RITVEC(N,J,EV,KAPPA,RNM,RITZ,BND,ALF,BET,S,WRK1,WRK2, * IERR,MSGLVL) C C.... COMMON BLOCK FOR S-VECTORS (SHOULD BE CONTAINED IN RITVEC ALSO) C C.... XV1, XV2 ARE TEMP ARRAYS NEEDED COMPUTING FOR SINGULAR VECTORS C.... (IF VECTORS.EQ.TRUE) C.... LEADING DIMENSIONS OF XV1 MUST BE .GE. (NROW+NCOL) C.... LEADING DIMENSIONS OF XV2 MUST BE THE ORDER OF A'A (NCOL) C REAL*8 XV1,XV2 COMMON /VECTORS/ XV1(500,500),XV2(500),NSIG C INTEGER N,J,EV,IERR,MSGLVL,NSIG REAL*8 KAPPA,RNM,RITZ(J),BND(J),ALF(J),BET(J), * S(J,J),WRK1(N),WRK2(N) C C.... SUBROUTINES: IMTQL2,STORE C.... BLAS ROUTINES: DAXPY,DCOPY,DSCAL C INTEGER RETRQ PARAMETER (RETRQ = 2) C INTEGER I,K,KOUNT C CALL DSCAL(J*J,0.0D0,S,1) DO 10 I = 1,J S(I,I) = 1.0D0 10 CONTINUE CALL DCOPY(J,ALF,1,WRK1,-1) CALL DCOPY(J-1,BET(2),1,WRK2(2),-1) CALL IMTQL2(J,J,WRK1,WRK2,S,IERR) IF (IERR.NE.0) RETURN C C.... ON RETURN WRK1 CONTAINS EIGENVALUES IN ASCENDING ORDER C.... AND S CONTAINS THE CORRESPONDING EIGENVECTORS C C.... REWIND(EV) C.... WRITE(EV)N,J,KAPPA KOUNT=0 DO 50 K = 1,J C IF (BND(K).LE.KAPPA*DABS(RITZ(K))) THEN C CALL DSCAL(N,0.0D0,WRK1,1) C DO 20 I = 1,J CALL STORE(N,RETRQ,I,WRK2) CALL DAXPY(N,S(J-I+1,K),WRK2,1,WRK1,1) 20 CONTINUE C C.... WRITE(EV)RITZ(K),BND(K),(WRK1(I),I=1,N) IF(EV.GT.0) WRITE(EV) (WRK1(I),I=1,N) KOUNT=KOUNT+1 DO 60 II=1,N XV1(II,KOUNT)=WRK1(II) 60 CONTINUE C ENDIF C 50 CONTINUE NSIG=KOUNT C.... CLOSE(EV) RETURN END C C @(#)STARTV.F 1.4 (BNP) 12/9/88 C REAL*8 FUNCTION STARTV(N,J,R,WRK,NQ,EPS,IERR,MSGLVL) INTEGER N,J,NQ(4),IERR,MSGLVL REAL*8 R(5*N),WRK(N),EPS C C.... THIS ROUTINE DELIVERS A STARTING VECTOR IN R AND RETURNS |R|; C.... IT RETURNS ZERO IF RANGE IS SPANNED. IERR GETS SET TO NON-ZERO C.... IF NO STARTING VECTOR WITHIN RANGE OF OPERATOR CAN BE FOUND. C C.... N DIMENSION OF THE EIGENPROBLEM C.... J STARTING INDEX FOR A LANCZOS RUN C.... R AN ARRAY CONTAINING [R(J),Q(J),Q(J-1),P(J),P(J-1)/MR(J)] C.... NQ(4) LOCATION POINTERS FOR THE ARRAY R C C.... BLAS ROUTINES: DAXPY,DDOT C.... SUBROUTINES: RANDOM C.... USER-SUPPLIED: OPB,OPM,STORE C INTEGER RETRQ PARAMETER (RETRQ = 2) C INTEGER IRAND,I,ID REAL*8 RNM2,T,RANDOM,DDOT C REAL*8 ONE,ZERO DATA ONE,ZERO/1.0D0,0.0D0/ C C.... GET INITIAL VECTOR, DEFAULT IS RANDOM C RNM2 = DDOT(N,R,1,R,1) IRAND = 918272+J DO 20 ID = 1,3 IF (ID.GT.1.OR.J.GT.1.OR.RNM2.EQ.ZERO) THEN DO 10 I = 1,N R(I) = RANDOM(IRAND) 10 CONTINUE ENDIF CALL OPM(N,R,R(NQ(3))) C C.... APPLY OPERATOR TO PUT R IN RANGE (ESSENTIAL IF M SINGULAR) C CALL OPB(N,R(NQ(3)),R) CALL OPM(N,R,R(NQ(3))) RNM2 = DDOT(N,R,1,R(NQ(3)),1) IF (RNM2.GT.ZERO) GOTO 30 20 CONTINUE C C.... FATAL ERROR C IERR = 8192 STARTV = -ONE RETURN 30 IF (J.GT.1) THEN DO 40 I = 1,J-1 CALL STORE(N,RETRQ,I,WRK) T = -DDOT(N,R(NQ(3)),1,WRK,1) CALL DAXPY(N,T,WRK,1,R,1) 40 CONTINUE C C.... MAKE SURE Q(J) IS ORTHOGONAL TO Q(J-1) C T = DDOT(N,R(NQ(4)),1,R,1) CALL DAXPY(N,-T,R(NQ(2)),1,R,1) C CALL OPM(N,R,R(NQ(3))) T = DDOT(N,R(NQ(3)),1,R,1) IF (T.LE.EPS*RNM2) T = ZERO RNM2 = T ENDIF STARTV = DSQRT(RNM2) RETURN END C C @(#)STORE.F 3.3 (BNP) 3/16/89; FROM STORE.F 2.2 10/13/87 C SUBROUTINE STORE(N,ISW,J,S) INTEGER N,ISW,J REAL*8 S(N) C INTEGER MAXLL PARAMETER (MAXLL = 2) C INTEGER STORQ,RETRQ,STORP,RETRP PARAMETER (STORQ = 1,RETRQ = 2,STORP = 3,RETRP = 4) C REAL*8 A COMMON/GETPUT/A(500000) C IF (ISW.EQ.STORQ) THEN CALL DCOPY(N,S,1,A((J+MAXLL-1)*N+1),1) ELSE IF (ISW.EQ.RETRQ) THEN CALL DCOPY(N,A((J+MAXLL-1)*N+1),1,S,1) ELSE IF (ISW.EQ.STORP) THEN IF (J.GT.MAXLL) STOP 'STORE: (STORP) J.GT.MAXLL' CALL DCOPY(N,S,1,A((J-1)*N+1),1) ELSE IF (ISW.EQ.RETRP) THEN IF (J.GT.MAXLL) STOP 'STORE: (RETRP) J.GT.MAXLL' CALL DCOPY(N,A((J-1)*N+1),1,S,1) ENDIF RETURN END C C @(#)STPONE.F 3.6 (BNP) 12/9/88; FROM STPONE.F 2.8 6/17/88 C SUBROUTINE STPONE(N,R,WRK,ALF,NQ,IERR,MSGLVL) INTEGER N,NQ(4),IERR,MSGLVL REAL*8 R(5*N),WRK(N),ALF(1) C C.... THIS ROUTINE PERFORMS THE FIRST STEP OF THE LANCZOS ALGORITHM. C.... IT PERFORMS A STEP OF EXTENDED LOCAL RE-ORTHOGONALIZATION. C C.... N DIMENSION OF THE EIGENPROBLEM C.... R AN ARRAY CONTAINING [R(J),Q(J),Q(J-1),P(J),P(J-1)/MR(J)] C.... ALF DIAGONAL ELEMENTS OF T C.... NQ(4) LOCATION POINTERS FOR THE ARRAY R C C.... BLAS ROUTINES: DATX,DAXPY,DDOT,DSCAL C.... SUBROUTINES: STARTV C.... USER-SUPPLIED: OPB,OPM,STORE C REAL*8 RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 COMMON/RDATA/RNM,ANORM,TOL,EPS,EPS1,REPS,EPS34 C REAL*8 T,STARTV,DDOT C REAL*8 ONE,ZERO DATA ONE,ZERO/1.0D0,0.0D0/ C C.... GET INITIAL VECTOR, DEFAULT IS RANDOM C RNM = STARTV(N,1,R,WRK,NQ,EPS,IERR,MSGLVL) IF (RNM.EQ.ZERO.OR.IERR.NE.0) RETURN C C.... NORMALIZE STARTING VECTOR C T = ONE/RNM CALL DATX(N,T,R,1,R(NQ(1)),1) CALL DSCAL(N,T,R(NQ(3)),1) C C.... TAKE THE FIRST STEP C CALL OPB(N,R(NQ(3)),R) ALF(1) = DDOT(N,R,1,R(NQ(3)),1) CALL DAXPY(N,-ALF(1),R(NQ(1)),1,R,1) T = DDOT(N,R,1,R(NQ(3)),1) CALL DAXPY(N,-T,R(NQ(1)),1,R,1) ALF(1) = ALF(1)+T CALL OPM(N,R,R(NQ(4))) RNM = DSQRT(DDOT(N,R,1,R(NQ(4)),1)) ANORM = RNM+DABS(ALF(1)) TOL = REPS*ANORM C RETURN END C SUBROUTINE DAXPY( N, DA, DX, INCX, DY, INCY ) * * CONSTANT TIMES A VECTOR PLUS A VECTOR. * USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. * JACK DONGARRA, LINPACK, 3/11/78. * * .. SCALAR ARGUMENTS .. INTEGER INCX, INCY, N REAL*8 DA * .. * .. ARRAY ARGUMENTS .. REAL*8 DX( 1 ), DY( 1 ) * .. * .. LOCAL SCALARS .. INTEGER I, IX, IY, M, MP1 * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MOD * .. * .. EXECUTABLE STATEMENTS .. * IF( N.LE.0 ) $ RETURN IF( DA.EQ.0.0D0 ) $ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) $ GO TO 20 * * CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS * NOT EQUAL TO 1 * 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 * * CODE FOR BOTH INCREMENTS EQUAL TO 1 * * * CLEAN-UP LOOP * 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 ) * * COPIES A VECTOR, X, TO A VECTOR, Y. * USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. * JACK DONGARRA, LINPACK, 3/11/78. * * .. SCALAR ARGUMENTS .. INTEGER INCX, INCY, N * .. * .. ARRAY ARGUMENTS .. REAL*8 DX( 1 ), DY( 1 ) * .. * .. LOCAL SCALARS .. INTEGER I, IX, IY, M, MP1 * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MOD * .. * .. EXECUTABLE STATEMENTS .. * IF( N.LE.0 ) $ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) $ GO TO 20 * * CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS * NOT EQUAL TO 1 * 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 * * CODE FOR BOTH INCREMENTS EQUAL TO 1 * * * CLEAN-UP LOOP * 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 REAL*8 FUNCTION DDOT( N, DX, INCX, DY, INCY ) * * FORMS THE DOT PRODUCT OF TWO VECTORS. * USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. * JACK DONGARRA, LINPACK, 3/11/78. * * .. SCALAR ARGUMENTS .. INTEGER INCX, INCY, N * .. * .. ARRAY ARGUMENTS .. REAL*8 DX( 1 ), DY( 1 ) * .. * .. LOCAL SCALARS .. INTEGER I, IX, IY, M, MP1 REAL*8 DTEMP * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MOD * .. * .. EXECUTABLE STATEMENTS .. * DDOT = 0.0D0 DTEMP = 0.0D0 IF( N.LE.0 ) $ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) $ GO TO 20 * * CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS * NOT EQUAL TO 1 * 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 * * CODE FOR BOTH INCREMENTS EQUAL TO 1 * * * CLEAN-UP LOOP * 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 SUBROUTINE DSCAL( N, DA, DX, INCX ) * * SCALES A VECTOR BY A CONSTANT. * USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. * JACK DONGARRA, LINPACK, 3/11/78. * * .. SCALAR ARGUMENTS .. INTEGER INCX, N REAL*8 DA * .. * .. ARRAY ARGUMENTS .. REAL*8 DX( 1 ) * .. * .. LOCAL SCALARS .. INTEGER I, IX, M, MP1, NINCX * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MOD * .. * .. EXECUTABLE STATEMENTS .. * IF( N.LE.0 ) $ RETURN IF( INCX.EQ.1 ) $ GO TO 20 * * CODE FOR INCREMENT NOT EQUAL TO 1 * IX = 1 IF( INCX.LT.0 ) $ IX = 1 - ( N-1 )*INCX NINCX = IX + ( N-1 )*INCX DO 10 I = IX, NINCX, INCX DX( I ) = DA*DX( I ) 10 CONTINUE RETURN * * CODE FOR INCREMENT EQUAL TO 1 * * * CLEAN-UP LOOP * 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 SUBROUTINE DSWAP( N, DX, INCX, DY, INCY ) * * INTERCHANGES TWO VECTORS. * USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. * JACK DONGARRA, LINPACK, 3/11/78. * * .. SCALAR ARGUMENTS .. INTEGER INCX, INCY, N * .. * .. ARRAY ARGUMENTS .. REAL*8 DX( 1 ), DY( 1 ) * .. * .. LOCAL SCALARS .. INTEGER I, IX, IY, M, MP1 REAL*8 DTEMP * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC MOD * .. * .. EXECUTABLE STATEMENTS .. * IF( N.LE.0 ) $ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) $ GO TO 20 * * CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS * NOT EQUAL TO 1 * 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 * * CODE FOR BOTH INCREMENTS EQUAL TO 1 * * * CLEAN-UP LOOP * 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 INTEGER FUNCTION IDAMAX( N, DX, INCX ) * * FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. * JACK DONGARRA, LINPACK, 3/11/78. * * .. SCALAR ARGUMENTS .. INTEGER INCX, N * .. * .. ARRAY ARGUMENTS .. REAL*8 DX( 1 ) * .. * .. LOCAL SCALARS .. INTEGER I, IX REAL*8 DMAX * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC DABS * .. * .. EXECUTABLE STATEMENTS .. * IDAMAX = 0 IF( N.LT.1 ) $ RETURN IDAMAX = 1 IF( N.EQ.1 ) $ RETURN IF( INCX.EQ.1 ) $ GO TO 30 * * CODE FOR INCREMENT NOT EQUAL TO 1 * IX = 1 IF( INCX.LT.0 ) $ IX = 1 - ( N-1 )*INCX DMAX = DABS( DX( IX ) ) IX = IX + INCX DO 20 I = 2, N IF( DABS( DX( IX ) ).LE.DMAX ) $ GO TO 10 IDAMAX = I DMAX = DABS( DX( IX ) ) 10 IX = IX + INCX 20 CONTINUE IF( INCX.LT.0 ) $ IDAMAX = N - IDAMAX + 1 RETURN * * CODE FOR INCREMENT EQUAL TO 1 * 30 DMAX = DABS( DX( 1 ) ) DO 40 I = 2, N IF( DABS( DX( I ) ).LE.DMAX ) $ GO TO 40 IDAMAX = I DMAX = DABS( DX( I ) ) 40 CONTINUE RETURN END