#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'README' <<'END_OF_FILE' XSince some compiler optimizations cause MACHAR to misbehave, XI've replaced the call on MACHAR in setup.f with corresponding Xcalls on PORT routines I1MACH and D1MACH, which now adapt Xthemselves automatically to most currently used machines. XThe original setup.f is now setup.f0, and the original machar.f Xis now machar.f0. X XAs of this writing, we do not know how to contact Rod Bain, the Xauthor of NNES. X X-- David M. Gay (dmg@bell-labs.com) X Bell Labs, Murray Hill X 8 February 1999 END_OF_FILE if test 485 -ne `wc -c <'README'`; then echo shar: \"'README'\" unpacked with wrong size! fi # end of 'README' fi if test -f 'abmul.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'abmul.f'\" else echo shar: Extracting \"'abmul.f'\" \(7666 characters\) sed "s/^X//" >'abmul.f' <<'END_OF_FILE' X SUBROUTINE ABMUL(NRADEC,NRAACT,NCBDEC,NCBACT,NCDEC ,NCACT , X $ AMAT ,BMAT ,CMAT ,AROW) XC XC FEB. 8, 1991 XC XC MATRIX MULTIPLICATION AB=C XC XC VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 XC EACH ROW OF MATRIX A IS SAVED AS A COLUMN, AROW, BEFORE USE. XC XC NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX XC NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT FOR 2ND INDEX XC NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT XC XC I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED XC NCBDEC IS THE NUMBER OF COLUMNS OF B DECLARED XC NCDEC IS THE COMMON DECLARED DIMENSION XC XC MODIFICATION OF MATRIX MULTIPLIER DONATED BY PROF. JAMES XC MACKINNON, QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION AMAT(NRADEC,NCDEC) ,BMAT(NCDEC,NCBDEC), X $ CMAT(NRADEC,NCBDEC) ,AROW(NCDEC) X DATA ZERO /0.0D0/ XC XC FIND NUMBER OF GROUPS OF SIZE 32, 16 ... XC X NCC32=NCACT/32 X NCC32R=NCACT-32*NCC32 X NCC16=NCC32R/16 X NCC16R=NCC32R-16*NCC16 X NCC8=NCC16R/8 X NCC8R=NCC16R-8*NCC8 X NCC4=NCC8R/4 X NCC4R=NCC8R-4*NCC4 XC XC REASSIGN ROWS TO VECTOR AROW. XC X DO 100 I=1,NRAACT X K=0 X IF(NCC32.GT.0) THEN X DO 200 KK=1,NCC32 X K=K+32 X AROW(K-31)=AMAT(I,K-31) X AROW(K-30)=AMAT(I,K-30) X AROW(K-29)=AMAT(I,K-29) X AROW(K-28)=AMAT(I,K-28) X AROW(K-27)=AMAT(I,K-27) X AROW(K-26)=AMAT(I,K-26) X AROW(K-25)=AMAT(I,K-25) X AROW(K-24)=AMAT(I,K-24) X AROW(K-23)=AMAT(I,K-23) X AROW(K-22)=AMAT(I,K-22) X AROW(K-21)=AMAT(I,K-21) X AROW(K-20)=AMAT(I,K-20) X AROW(K-19)=AMAT(I,K-19) X AROW(K-18)=AMAT(I,K-18) X AROW(K-17)=AMAT(I,K-17) X AROW(K-16)=AMAT(I,K-16) X AROW(K-15)=AMAT(I,K-15) X AROW(K-14)=AMAT(I,K-14) X AROW(K-13)=AMAT(I,K-13) X AROW(K-12)=AMAT(I,K-12) X AROW(K-11)=AMAT(I,K-11) X AROW(K-10)=AMAT(I,K-10) X AROW(K-9)=AMAT(I,K-9) X AROW(K-8)=AMAT(I,K-8) X AROW(K-7)=AMAT(I,K-7) X AROW(K-6)=AMAT(I,K-6) X AROW(K-5)=AMAT(I,K-5) X AROW(K-4)=AMAT(I,K-4) X AROW(K-3)=AMAT(I,K-3) X AROW(K-2)=AMAT(I,K-2) X AROW(K-1)=AMAT(I,K-1) X AROW(K)=AMAT(I,K) X200 CONTINUE X END IF X IF(NCC16.GT.0) THEN X DO 300 KK=1,NCC16 X K=K+16 X AROW(K-15)=AMAT(I,K-15) X AROW(K-14)=AMAT(I,K-14) X AROW(K-13)=AMAT(I,K-13) X AROW(K-12)=AMAT(I,K-12) X AROW(K-11)=AMAT(I,K-11) X AROW(K-10)=AMAT(I,K-10) X AROW(K-9)=AMAT(I,K-9) X AROW(K-8)=AMAT(I,K-8) X AROW(K-7)=AMAT(I,K-7) X AROW(K-6)=AMAT(I,K-6) X AROW(K-5)=AMAT(I,K-5) X AROW(K-4)=AMAT(I,K-4) X AROW(K-3)=AMAT(I,K-3) X AROW(K-2)=AMAT(I,K-2) X AROW(K-1)=AMAT(I,K-1) X AROW(K)=AMAT(I,K) X300 CONTINUE X END IF X IF(NCC8.GT.0) THEN X DO 400 KK=1,NCC8 X K=K+8 X AROW(K-7)=AMAT(I,K-7) X AROW(K-6)=AMAT(I,K-6) X AROW(K-5)=AMAT(I,K-5) X AROW(K-4)=AMAT(I,K-4) X AROW(K-3)=AMAT(I,K-3) X AROW(K-2)=AMAT(I,K-2) X AROW(K-1)=AMAT(I,K-1) X AROW(K)=AMAT(I,K) X400 CONTINUE X END IF X IF(NCC4.GT.0) THEN X DO 500 KK=1,NCC4 X K=K+4 X AROW(K-3)=AMAT(I,K-3) X AROW(K-2)=AMAT(I,K-2) X AROW(K-1)=AMAT(I,K-1) X AROW(K)=AMAT(I,K) X500 CONTINUE X END IF X IF(NCC4R.GT.0) THEN X DO 600 KK=1,NCC4R X K=K+1 X AROW(K)=AMAT(I,K) X600 CONTINUE X END IF XC XC FIND ENTRY FOR MATRIX C USING COLUMN VECTOR AROW. XC X DO 700 J=1,NCBACT X SUM=ZERO X K=0 X IF(NCC32.GT.0) THEN X DO 800 KK=1,NCC32 X K=K+32 X SUM=SUM X $ +AROW(K-31)*BMAT(K-31,J)+AROW(K-30)*BMAT(K-30,J) X $ +AROW(K-29)*BMAT(K-29,J)+AROW(K-28)*BMAT(K-28,J) X $ +AROW(K-27)*BMAT(K-27,J)+AROW(K-26)*BMAT(K-26,J) X $ +AROW(K-25)*BMAT(K-25,J)+AROW(K-24)*BMAT(K-24,J) X SUM=SUM X $ +AROW(K-23)*BMAT(K-23,J)+AROW(K-22)*BMAT(K-22,J) X $ +AROW(K-21)*BMAT(K-21,J)+AROW(K-20)*BMAT(K-20,J) X $ +AROW(K-19)*BMAT(K-19,J)+AROW(K-18)*BMAT(K-18,J) X $ +AROW(K-17)*BMAT(K-17,J)+AROW(K-16)*BMAT(K-16,J) X SUM=SUM X $ +AROW(K-15)*BMAT(K-15,J)+AROW(K-14)*BMAT(K-14,J) X $ +AROW(K-13)*BMAT(K-13,J)+AROW(K-12)*BMAT(K-12,J) X $ +AROW(K-11)*BMAT(K-11,J)+AROW(K-10)*BMAT(K-10,J) X $ +AROW(K-9) *BMAT(K-9,J) +AROW(K-8) *BMAT(K-8,J) X SUM=SUM X $ +AROW(K-7)*BMAT(K-7,J)+AROW(K-6)*BMAT(K-6,J) X $ +AROW(K-5)*BMAT(K-5,J)+AROW(K-4)*BMAT(K-4,J) X $ +AROW(K-3)*BMAT(K-3,J)+AROW(K-2)*BMAT(K-2,J) X $ +AROW(K-1)*BMAT(K-1,J)+AROW(K) *BMAT(K,J) X800 CONTINUE X END IF X IF(NCC16.GT.0) THEN X DO 900 KK=1,NCC16 X K=K+16 X SUM=SUM X $ +AROW(K-15)*BMAT(K-15,J)+AROW(K-14)*BMAT(K-14,J) X $ +AROW(K-13)*BMAT(K-13,J)+AROW(K-12)*BMAT(K-12,J) X $ +AROW(K-11)*BMAT(K-11,J)+AROW(K-10)*BMAT(K-10,J) X $ +AROW(K-9) *BMAT(K-9,J) +AROW(K-8) *BMAT(K-8,J) X SUM=SUM X $ +AROW(K-7)*BMAT(K-7,J)+AROW(K-6)*BMAT(K-6,J) X $ +AROW(K-5)*BMAT(K-5,J)+AROW(K-4)*BMAT(K-4,J) X $ +AROW(K-3)*BMAT(K-3,J)+AROW(K-2)*BMAT(K-2,J) X $ +AROW(K-1)*BMAT(K-1,J)+AROW(K) *BMAT(K,J) X900 CONTINUE X END IF X IF(NCC8.GT.0) THEN X DO 1000 KK=1,NCC8 X K=K+8 X SUM=SUM X $ +AROW(K-7)*BMAT(K-7,J)+AROW(K-6)*BMAT(K-6,J) X $ +AROW(K-5)*BMAT(K-5,J)+AROW(K-4)*BMAT(K-4,J) X $ +AROW(K-3)*BMAT(K-3,J)+AROW(K-2)*BMAT(K-2,J) X $ +AROW(K-1)*BMAT(K-1,J)+AROW(K) *BMAT(K,J) X1000 CONTINUE X END IF X IF(NCC4.GT.0) THEN X DO 1100 KK=1,NCC4 X K=K+4 X SUM=SUM X $ +AROW(K-3)*BMAT(K-3,J)+AROW(K-2)*BMAT(K-2,J) X $ +AROW(K-1)*BMAT(K-1,J)+AROW(K) *BMAT(K,J) X1100 CONTINUE X END IF X IF(NCC4R.GT.0) THEN X DO 1200 KK=1,NCC4R X K=K+1 X SUM=SUM+AROW(K)*BMAT(K,J) X1200 CONTINUE X END IF X CMAT(I,J)=SUM X700 CONTINUE X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ABMUL. XC X END END_OF_FILE if test 7666 -ne `wc -c <'abmul.f'`; then echo shar: \"'abmul.f'\" unpacked with wrong size! fi # end of 'abmul.f' fi if test -f 'ascalf.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ascalf.f'\" else echo shar: Extracting \"'ascalf.f'\" \(1363 characters\) sed "s/^X//" >'ascalf.f' <<'END_OF_FILE' X SUBROUTINE ASCALF(N,EPSMCH,FVECC,JAC,SCALEF) XC XC FEB. 13, 1991 XC XC THIS SUBROUTINE ESTABLISHES SCALING FACTORS FOR THE XC RESIDUAL VECTOR IF FUNCTION ADAPTIVE SCALING IS CHOSEN XC USING INTEGER VARIABLE ITSCLF. XC XC NOTE: IN QUASI-NEWTON METHODS THE SCALING FACTORS ARE XC UPDATED ONLY WHEN THE JACOBIAN IS EVALUATED EXPLI- XC CITLY. XC XC SCALING FACTORS ARE DETERMINED FROM THE INFINITY NORMS XC OF THE ROWS OF THE JACOBIAN AND THE VALUES OF THE CURRENT XC FUNCTION VECTOR, FVECC. XC XC A MINIMUM TOLERANCE ON THE SCALING FACTOR IS THE SQUARE XC ROOT OF THE MACHINE PRECISION, SQRTEP. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X DIMENSION FVECC(N) ,SCALEF(N) X DATA ZERO,ONE /0.0D0,1.0D0/ XC X SQRTEP=SQRT(EPSMCH) XC XC I COUNTS THE ROWS. XC X DO 100 I=1,N X AMAX=ZERO XC XC FIND MAXIMUM ENTRY IN ROW I. XC X DO 200 J=1,N X AMAX=MAX(AMAX,ABS(JAC(I,J))) X200 CONTINUE XC X AMAX=MAX(AMAX,FVECC(I)) XC XC SET SCALING FACTOR TO A DEFAULT OF ONE IF ITH ROW IS ZEROS. XC X IF(AMAX.EQ.ZERO) AMAX=ONE X AMAX=MAX(AMAX,SQRTEP) X SCALEF(I)=ONE/AMAX X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ASCALF. XC X END END_OF_FILE if test 1363 -ne `wc -c <'ascalf.f'`; then echo shar: \"'ascalf.f'\" unpacked with wrong size! fi # end of 'ascalf.f' fi if test -f 'ascalx.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ascalx.f'\" else echo shar: Extracting \"'ascalx.f'\" \(1196 characters\) sed "s/^X//" >'ascalx.f' <<'END_OF_FILE' X SUBROUTINE ASCALX(N,EPSMCH,JAC,SCALEX) XC XC FEB. 13, 1991 XC XC THIS SUBROUTINE ESTABLISHES SCALING FACTORS FOR THE XC COMPONENET VECTOR IF ADAPTIVE SCALING IS CHOSEN USING XC INTEGER ITSCLX. XC XC NOTE: IN QUASI-NEWTON METHODS THE SCALING FACTORS ARE XC UPDATED ONLY WHEN THE JACOBIAN IS EVALUATED EXPLI- XC CITLY. XC XC SCALING FACTORS ARE DETERMINED FROM THE INFINITY NORMS XC OF THE COLUMNS OF THE JACOBIAN. XC XC A MINIMUM TOLERANCE ON THE SCALING FACTOR IS THE SQUARE XC ROOT OF THE MACHINE PRECISION, SQRTEP. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X DIMENSION SCALEX(N) X DATA ZERO,ONE /0.0D0,1.0D0/ XC X SQRTEP=SQRT(EPSMCH) XC XC J COUNTS COLUMNS. XC X DO 100 J=1,N X AMAX=ZERO XC XC FIND MAXIMUM ENTRY IN JTH COLUMN. XC X DO 200 I=1,N X AMAX=MAX(AMAX,ABS(JAC(I,J))) X200 CONTINUE XC XC IF A COLUMN IS ALL ZEROS SET AMAX TO ONE. XC X IF(AMAX.EQ.ZERO) AMAX=ONE X SCALEX(J)=MAX(AMAX,SQRTEP) X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ASCALX. XC X END END_OF_FILE if test 1196 -ne `wc -c <'ascalx.f'`; then echo shar: \"'ascalx.f'\" unpacked with wrong size! fi # end of 'ascalx.f' fi if test -f 'atamul.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'atamul.f'\" else echo shar: Extracting \"'atamul.f'\" \(4482 characters\) sed "s/^X//" >'atamul.f' <<'END_OF_FILE' X SUBROUTINE ATAMUL(NRADEC,NCADEC,NRAACT,NCAACT,NRBDEC,NCBDEC, X $ AMAT ,BMAT) XC XC FEB. 8, 1991 XC XC MATRIX MULTIPLICATION: A^A=B XC XC VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. XC XC NRADEC IS NUMBER OF ROWS IN A DECLARED XC NCADEC IS NUMBER OF COLUMNS IN A DECLARED XC NRAACT IS THE LIMIT FOR THE 1ST INDEX IN A XC NCAACT IS THE LIMIT FOR THE 2ND INDEX IN A XC NRBDEC IS NUMBER OF ROWS IN B DECLARED XC NCBDEC IS NUMBER OF COLUMNS IN B DECLARED XC XC MODIFIED VERSION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES XC MACKINNON, QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION AMAT(NRADEC,NCADEC), BMAT(NRBDEC,NCBDEC) X DATA ZERO /0.0D0/ XC XC FIND NUMBER OF GROUPS OF SIZE 32, 16 ... XC X NCC32=NRAACT/32 X NCC32R=NRAACT-32*NCC32 X NCC16=NCC32R/16 X NCC16R=NCC32R-16*NCC16 X NCC8=NCC16R/8 X NCC8R=NCC16R-8*NCC8 X NCC4=NCC8R/4 X NCC4R=NCC8R-4*NCC4 XC XC FIND ENTRY IN MATRIX B. XC X DO 100 I=1,NCAACT X DO 200 J=I,NCAACT X SUM=ZERO X K=0 X IF(NCC32.GT.0) THEN X DO 300 KK=1,NCC32 X K=K+32 X SUM=SUM X $ +AMAT(K-31,I)*AMAT(K-31,J)+AMAT(K-30,I)*AMAT(K-30,J) X $ +AMAT(K-29,I)*AMAT(K-29,J)+AMAT(K-28,I)*AMAT(K-28,J) X $ +AMAT(K-27,I)*AMAT(K-27,J)+AMAT(K-26,I)*AMAT(K-26,J) X $ +AMAT(K-25,I)*AMAT(K-25,J)+AMAT(K-24,I)*AMAT(K-24,J) X SUM=SUM X $ +AMAT(K-23,I)*AMAT(K-23,J)+AMAT(K-22,I)*AMAT(K-22,J) X $ +AMAT(K-21,I)*AMAT(K-21,J)+AMAT(K-20,I)*AMAT(K-20,J) X $ +AMAT(K-19,I)*AMAT(K-19,J)+AMAT(K-18,I)*AMAT(K-18,J) X $ +AMAT(K-17,I)*AMAT(K-17,J)+AMAT(K-16,I)*AMAT(K-16,J) X SUM=SUM X $ +AMAT(K-15,I)*AMAT(K-15,J)+AMAT(K-14,I)*AMAT(K-14,J) X $ +AMAT(K-13,I)*AMAT(K-13,J)+AMAT(K-12,I)*AMAT(K-12,J) X $ +AMAT(K-11,I)*AMAT(K-11,J)+AMAT(K-10,I)*AMAT(K-10,J) X $ +AMAT(K-9,I)*AMAT(K-9,J) +AMAT(K-8,I)*AMAT(K-8,J) X SUM=SUM X $ +AMAT(K-7,I)*AMAT(K-7,J)+AMAT(K-6,I)*AMAT(K-6,J) X $ +AMAT(K-5,I)*AMAT(K-5,J)+AMAT(K-4,I)*AMAT(K-4,J) X $ +AMAT(K-3,I)*AMAT(K-3,J)+AMAT(K-2,I)*AMAT(K-2,J) X $ +AMAT(K-1,I)*AMAT(K-1,J)+AMAT(K,I) *AMAT(K,J) X300 CONTINUE X END IF X IF(NCC16.GT.0) THEN X DO 400 KK=1,NCC16 X K=K+16 X SUM=SUM X $ +AMAT(K-15,I)*AMAT(K-15,J)+AMAT(K-14,I)*AMAT(K-14,J) X $ +AMAT(K-13,I)*AMAT(K-13,J)+AMAT(K-12,I)*AMAT(K-12,J) X $ +AMAT(K-11,I)*AMAT(K-11,J)+AMAT(K-10,I)*AMAT(K-10,J) X $ +AMAT(K-9,I)*AMAT(K-9,J) +AMAT(K-8,I) *AMAT(K-8,J) X SUM=SUM X $ +AMAT(K-7,I)*AMAT(K-7,J)+AMAT(K-6,I)*AMAT(K-6,J) X $ +AMAT(K-5,I)*AMAT(K-5,J)+AMAT(K-4,I)*AMAT(K-4,J) X $ +AMAT(K-3,I)*AMAT(K-3,J)+AMAT(K-2,I)*AMAT(K-2,J) X $ +AMAT(K-1,I)*AMAT(K-1,J)+AMAT(K,I)*AMAT(K,J) X400 CONTINUE X END IF X IF(NCC8.GT.0) THEN X DO 500 KK=1,NCC8 X K=K+8 X SUM=SUM X $ +AMAT(K-7,I)*AMAT(K-7,J)+AMAT(K-6,I)*AMAT(K-6,J) X $ +AMAT(K-5,I)*AMAT(K-5,J)+AMAT(K-4,I)*AMAT(K-4,J) X $ +AMAT(K-3,I)*AMAT(K-3,J)+AMAT(K-2,I)*AMAT(K-2,J) X $ +AMAT(K-1,I)*AMAT(K-1,J)+AMAT(K,I)*AMAT(K,J) X500 CONTINUE X END IF X IF(NCC4.GT.0) THEN X DO 600 KK=1,NCC4 X K=K+4 X SUM=SUM X $ +AMAT(K-3,I)*AMAT(K-3,J)+AMAT(K-2,I)*AMAT(K-2,J) X $ +AMAT(K-1,I)*AMAT(K-1,J)+AMAT(K,I)*AMAT(K,J) X600 CONTINUE X END IF X IF(NCC4R.GT.0) THEN X DO 700 KK=1,NCC4R X K=K+1 X SUM=SUM+AMAT(K,I)*AMAT(K,J) X700 CONTINUE X END IF X BMAT(I,J)=SUM X IF(I.NE.J) BMAT(J,I)=BMAT(I,J) X200 CONTINUE X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ATAMUL. XC X END END_OF_FILE if test 4482 -ne `wc -c <'atamul.f'`; then echo shar: \"'atamul.f'\" unpacked with wrong size! fi # end of 'atamul.f' fi if test -f 'ataov.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ataov.f'\" else echo shar: Extracting \"'ataov.f'\" \(1786 characters\) sed "s/^X//" >'ataov.f' <<'END_OF_FILE' X SUBROUTINE ATAOV(OVERFL,MAXEXP,N,NUNIT,OUTPUT,A,B,SCALEF) XC XC SEPT. 8, 1991 XC XC THIS SUBROUTINE FINDS THE PRODUCT OF THE TRANSPOSE OF THE XC MATRIX A AND MATRIX A. EACH ENTRY IS CHECKED BEFORE BEING XC ACCEPTED. IF IT WOULD CAUSE AN OVERFLOW 10**MAXEXP IS XC INSERTED IN ITS PLACE. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X INTEGER OUTPUT X DIMENSION A(N,N) ,B(N,N) ,SCALEF(N) X LOGICAL OVERFL ,WRNSUP X COMMON/NNES_2/WRNSUP X DATA ZERO,ONE,TWO,TEN /0.0D0,1.0D0,2.0D0,10.0D0/ XC X EPS=TEN**(-MAXEXP) X OVERFL=.FALSE. XC X DO 100 I=1,N X DO 200 J=I+1,N X SUM=ZERO X DO 300 K=1,N X IF(LOG10(ABS(A(K,I))+EPS)+LOG10(ABS(A(K,J))+EPS) X $ +TWO*LOG10(SCALEF(K)).GT.MAXEXP) THEN X OVERFL=.TRUE. X B(I,J)=SIGN(TEN**MAXEXP,A(K,I))*SIGN(ONE,A(K,J)) X IF(OUTPUT.GT.2.AND.(.NOT.WRNSUP)) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,2) B(I,J) X2 FORMAT(T3,'*',4X,'WARNING: COMPONENT IN', X $ ' MATRIX-MATRIX PRODUCT SET TO ',1PD12.3, X $ T74,'*') X END IF X GO TO 201 X END IF X SUM=SUM+A(K,I)*A(K,J)*SCALEF(K)*SCALEF(K) X300 CONTINUE X B(I,J)=SUM X B(J,I)=SUM X201 CONTINUE X200 CONTINUE X SUM=ZERO X DO 400 K=1,N X IF(TWO*(LOG10(ABS(A(K,I))+EPS)+LOG10(SCALEF(K))). X $ GT.MAXEXP) THEN X OVERFL=.TRUE. X B(I,I)=TEN**MAXEXP X IF(OUTPUT.GT.2.AND.(.NOT.WRNSUP)) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) B(I,I) X END IF X GO TO 401 X END IF X SUM=SUM+A(K,I)*A(K,I)*SCALEF(K)*SCALEF(K) X400 CONTINUE X B(I,I)=SUM X401 CONTINUE X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ATAOV. XC X END END_OF_FILE if test 1786 -ne `wc -c <'ataov.f'`; then echo shar: \"'ataov.f'\" unpacked with wrong size! fi # end of 'ataov.f' fi if test -f 'atbmul.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'atbmul.f'\" else echo shar: Extracting \"'atbmul.f'\" \(4592 characters\) sed "s/^X//" >'atbmul.f' <<'END_OF_FILE' X SUBROUTINE ATBMUL(NCADEC,NCAACT,NCBDEC,NCBACT,NCDEC,NCACT, X $ AMAT ,BMAT ,CMAT) XC XC FEB. 8, 1991 XC XC MATRIX MULTIPLICATION: A^B=C XC XC VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. XC XC NCADEC IS 2ND DIM. OF AMAT; NCAACT IS ACTUAL LIMIT FOR 2ND INDEX XC NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT FOR 2ND INDEX XC NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT XC XC I.E. NCADEC IS NUMBER OF COLUMNS OF A DECLARED XC NCBDEC IS NUMBER OF COLUMNS OF B DECLARED XC NCDEC IS THE NUMBER OF ROWS IN BOTH A AND B DECLARED XC XC MODIFIED VERSION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES XC MACKINNON, QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION AMAT(NCDEC,NCADEC), BMAT(NCDEC,NCBDEC), X $ CMAT(NCADEC,NCBDEC) X DATA ZERO /0.0D0/ XC XC FIND NUMBER OF GROUPS OF SIZE 32, 16 ... XC X NCC32=NCACT/32 X NCC32R=NCACT-32*NCC32 X NCC16=NCC32R/16 X NCC16R=NCC32R-16*NCC16 X NCC8=NCC16R/8 X NCC8R=NCC16R-8*NCC8 X NCC4=NCC8R/4 X NCC4R=NCC8R-4*NCC4 XC XC FIND ENTRY IN MATRIX C. XC X DO 100 I=1,NCAACT X DO 200 J=1,NCBACT X SUM=ZERO X K=0 X IF(NCC32.GT.0) THEN X DO 300 KK=1,NCC32 X K=K+32 X SUM=SUM X $ +AMAT(K-31,I)*BMAT(K-31,J)+AMAT(K-30,I)*BMAT(K-30,J) X $ +AMAT(K-29,I)*BMAT(K-29,J)+AMAT(K-28,I)*BMAT(K-28,J) X $ +AMAT(K-27,I)*BMAT(K-27,J)+AMAT(K-26,I)*BMAT(K-26,J) X $ +AMAT(K-25,I)*BMAT(K-25,J)+AMAT(K-24,I)*BMAT(K-24,J) X SUM=SUM X $ +AMAT(K-23,I)*BMAT(K-23,J)+AMAT(K-22,I)*BMAT(K-22,J) X $ +AMAT(K-21,I)*BMAT(K-21,J)+AMAT(K-20,I)*BMAT(K-20,J) X $ +AMAT(K-19,I)*BMAT(K-19,J)+AMAT(K-18,I)*BMAT(K-18,J) X $ +AMAT(K-17,I)*BMAT(K-17,J)+AMAT(K-16,I)*BMAT(K-16,J) X SUM=SUM X $ +AMAT(K-15,I)*BMAT(K-15,J)+AMAT(K-14,I)*BMAT(K-14,J) X $ +AMAT(K-13,I)*BMAT(K-13,J)+AMAT(K-12,I)*BMAT(K-12,J) X $ +AMAT(K-11,I)*BMAT(K-11,J)+AMAT(K-10,I)*BMAT(K-10,J) X $ +AMAT(K-9,I)*BMAT(K-9,J) +AMAT(K-8,I) *BMAT(K-8,J) X SUM=SUM X $ +AMAT(K-7,I)*BMAT(K-7,J)+AMAT(K-6,I)*BMAT(K-6,J) X $ +AMAT(K-5,I)*BMAT(K-5,J)+AMAT(K-4,I)*BMAT(K-4,J) X $ +AMAT(K-3,I)*BMAT(K-3,J)+AMAT(K-2,I)*BMAT(K-2,J) X $ +AMAT(K-1,I)*BMAT(K-1,J)+AMAT(K,I) *BMAT(K,J) X300 CONTINUE X END IF X IF(NCC16.GT.0) THEN X DO 400 KK=1,NCC16 X K=K+16 X SUM=SUM X $ +AMAT(K-15,I)*BMAT(K-15,J)+AMAT(K-14,I)*BMAT(K-14,J) X $ +AMAT(K-13,I)*BMAT(K-13,J)+AMAT(K-12,I)*BMAT(K-12,J) X $ +AMAT(K-11,I)*BMAT(K-11,J)+AMAT(K-10,I)*BMAT(K-10,J) X $ +AMAT(K-9,I)*BMAT(K-9,J) +AMAT(K-8,I) *BMAT(K-8,J) X SUM=SUM X $ +AMAT(K-7,I)*BMAT(K-7,J)+AMAT(K-6,I)*BMAT(K-6,J) X $ +AMAT(K-5,I)*BMAT(K-5,J)+AMAT(K-4,I)*BMAT(K-4,J) X $ +AMAT(K-3,I)*BMAT(K-3,J)+AMAT(K-2,I)*BMAT(K-2,J) X $ +AMAT(K-1,I)*BMAT(K-1,J)+AMAT(K,I) *BMAT(K,J) X400 CONTINUE X END IF X IF(NCC8.GT.0) THEN X DO 500 KK=1,NCC8 X K=K+8 X SUM=SUM X $ +AMAT(K-7,I)*BMAT(K-7,J)+AMAT(K-6,I)*BMAT(K-6,J) X $ +AMAT(K-5,I)*BMAT(K-5,J)+AMAT(K-4,I)*BMAT(K-4,J) X $ +AMAT(K-3,I)*BMAT(K-3,J)+AMAT(K-2,I)*BMAT(K-2,J) X $ +AMAT(K-1,I)*BMAT(K-1,J)+AMAT(K,I) *BMAT(K,J) X500 CONTINUE X END IF X IF(NCC4.GT.0) THEN X DO 600 KK=1,NCC4 X K=K+4 X SUM=SUM X $ +AMAT(K-3,I)*BMAT(K-3,J)+AMAT(K-2,I)*BMAT(K-2,J) X $ +AMAT(K-1,I)*BMAT(K-1,J)+AMAT(K,I) *BMAT(K,J) X600 CONTINUE X END IF X IF(NCC4R.GT.0) THEN X DO 700 KK=1,NCC4R X K=K+1 X SUM=SUM+AMAT(K,I)*BMAT(K,J) X700 CONTINUE X END IF X CMAT(I,J)=SUM X200 CONTINUE X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ATBMUL. XC X END X END_OF_FILE if test 4592 -ne `wc -c <'atbmul.f'`; then echo shar: \"'atbmul.f'\" unpacked with wrong size! fi # end of 'atbmul.f' fi if test -f 'atvov.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'atvov.f'\" else echo shar: Extracting \"'atvov.f'\" \(1198 characters\) sed "s/^X//" >'atvov.f' <<'END_OF_FILE' X SUBROUTINE ATVOV(OVERFL,MAXEXP,N,NUNIT,OUTPUT,AMAT,BVEC,CVEC) XC XC FEB. 8 ,1991 XC XC THIS SUBROUTINE FINDS THE PRODUCT OF THE TRANSPOSE OF THE XC MATRIX A AND THE VECTOR B WHERE EACH ENTRY IS CHECKED TO XC PREVENT OVERFLOWS. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X INTEGER OUTPUT X DIMENSION AMAT(N,N) ,BVEC(N) ,CVEC(N) X LOGICAL OVERFL ,WRNSUP X COMMON/NNES_2/WRNSUP X DATA ZERO,ONE,TEN /0.0D0,1.0D0,10.0D0/ XC X EPS=TEN**(-MAXEXP) X OVERFL=.FALSE. XC X DO 100 I=1,N X SUM=ZERO X DO 200 J=1,N X IF(LOG10(ABS(AMAT(J,I))+EPS)+LOG10(ABS(BVEC(J))+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X CVEC(I)=SIGN(TEN**MAXEXP,AMAT(J,I)) X $ *SIGN(ONE,BVEC(J)) X IF(OUTPUT.GT.2.AND.(.NOT.WRNSUP)) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,2) CVEC(I) X2 FORMAT(T3,'*',4X,'WARNING: COMPONENT IN', X $ ' MATRIX-VECTOR PRODUCT SET TO ',1PD12.3,T74,'*') X END IF X GO TO 101 X END IF X SUM=SUM+AMAT(J,I)*BVEC(J) X200 CONTINUE X CVEC(I)=SUM X101 CONTINUE X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE ATVOV. XC X END END_OF_FILE if test 1198 -ne `wc -c <'atvov.f'`; then echo shar: \"'atvov.f'\" unpacked with wrong size! fi # end of 'atvov.f' fi if test -f 'avmul.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'avmul.f'\" else echo shar: Extracting \"'avmul.f'\" \(4104 characters\) sed "s/^X//" >'avmul.f' <<'END_OF_FILE' X SUBROUTINE AVMUL(NRADEC,NRAACT,NCDEC ,NCACT ,AMAT ,BVEC ,CVEC) XC XC FEB. 8, 1991 XC XC MATRIX-VECTOR MULTIPLICATION AB=C XC XC VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 XC EACH ROW OF MATRIX A IS SAVED AS A COLUMN BEFORE USE. XC XC NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX XC NCDEC IS COMMON DIMENSION OF AMAT & BVEC; NCACT IS ACTUAL LIMIT XC XC I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED XC NCDEC IS THE COMMON DECLARED DIMENSION (COLUMNS OF A AND XC ROWS OF B) XC XC MODIFICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES XC MACKINNON, QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION AMAT(NRADEC,NCDEC), BVEC(NCDEC), CVEC(NRADEC) X DATA ZERO /0.0D0/ XC XC FIND NUMBER OF GROUPS OF SIZE 32, 16 ... XC X NCC32=NCACT/32 X NCC32R=NCACT-32*NCC32 X NCC16=NCC32R/16 X NCC16R=NCC32R-16*NCC16 X NCC8=NCC16R/8 X NCC8R=NCC16R-8*NCC8 X NCC4=NCC8R/4 X NCC4R=NCC8R-4*NCC4 X DO 100 I=1,NRAACT XC XC FIND ENTRY FOR VECTOR C. XC X SUM=ZERO X K=0 X IF(NCC32.GT.0) THEN X DO 200 KK=1,NCC32 X K=K+32 X SUM=SUM X $ +AMAT(I,K-31)*BVEC(K-31)+AMAT(I,K-30)*BVEC(K-30) X $ +AMAT(I,K-29)*BVEC(K-29)+AMAT(I,K-28)*BVEC(K-28) X $ +AMAT(I,K-27)*BVEC(K-27)+AMAT(I,K-26)*BVEC(K-26) X $ +AMAT(I,K-25)*BVEC(K-25)+AMAT(I,K-24)*BVEC(K-24) X SUM=SUM X $ +AMAT(I,K-23)*BVEC(K-23)+AMAT(I,K-22)*BVEC(K-22) X $ +AMAT(I,K-21)*BVEC(K-21)+AMAT(I,K-20)*BVEC(K-20) X $ +AMAT(I,K-19)*BVEC(K-19)+AMAT(I,K-18)*BVEC(K-18) X $ +AMAT(I,K-17)*BVEC(K-17)+AMAT(I,K-16)*BVEC(K-16) X SUM=SUM X $ +AMAT(I,K-15)*BVEC(K-15)+AMAT(I,K-14)*BVEC(K-14) X $ +AMAT(I,K-13)*BVEC(K-13)+AMAT(I,K-12)*BVEC(K-12) X $ +AMAT(I,K-11)*BVEC(K-11)+AMAT(I,K-10)*BVEC(K-10) X $ +AMAT(I,K-9)*BVEC(K-9) +AMAT(I,K-8) *BVEC(K-8) X SUM=SUM X $ +AMAT(I,K-7)*BVEC(K-7)+AMAT(I,K-6)*BVEC(K-6) X $ +AMAT(I,K-5)*BVEC(K-5)+AMAT(I,K-4)*BVEC(K-4) X $ +AMAT(I,K-3)*BVEC(K-3)+AMAT(I,K-2)*BVEC(K-2) X $ +AMAT(I,K-1)*BVEC(K-1)+AMAT(I,K) *BVEC(K) X200 CONTINUE X END IF X IF(NCC16.GT.0) THEN X DO 300 KK=1,NCC16 X K=K+16 X SUM=SUM X $ +AMAT(I,K-15)*BVEC(K-15)+AMAT(I,K-14)*BVEC(K-14) X $ +AMAT(I,K-13)*BVEC(K-13)+AMAT(I,K-12)*BVEC(K-12) X $ +AMAT(I,K-11)*BVEC(K-11)+AMAT(I,K-10)*BVEC(K-10) X $ +AMAT(I,K-9)*BVEC(K-9) +AMAT(I,K-8) *BVEC(K-8) X SUM=SUM X $ +AMAT(I,K-7)*BVEC(K-7)+AMAT(I,K-6)*BVEC(K-6) X $ +AMAT(I,K-5)*BVEC(K-5)+AMAT(I,K-4)*BVEC(K-4) X $ +AMAT(I,K-3)*BVEC(K-3)+AMAT(I,K-2)*BVEC(K-2) X $ +AMAT(I,K-1)*BVEC(K-1)+AMAT(I,K) *BVEC(K) X300 CONTINUE X END IF X IF(NCC8.GT.0) THEN X DO 400 KK=1,NCC8 X K=K+8 X SUM=SUM X $ +AMAT(I,K-7)*BVEC(K-7)+AMAT(I,K-6)*BVEC(K-6) X $ +AMAT(I,K-5)*BVEC(K-5)+AMAT(I,K-4)*BVEC(K-4) X $ +AMAT(I,K-3)*BVEC(K-3)+AMAT(I,K-2)*BVEC(K-2) X $ +AMAT(I,K-1)*BVEC(K-1)+AMAT(I,K) *BVEC(K) X400 CONTINUE X END IF X IF(NCC4.GT.0) THEN X DO 500 KK=1,NCC4 X K=K+4 X SUM=SUM X $ +AMAT(I,K-3)*BVEC(K-3)+AMAT(I,K-2)*BVEC(K-2) X $ +AMAT(I,K-1)*BVEC(K-1)+AMAT(I,K) *BVEC(K) X500 CONTINUE X END IF X IF(NCC4R.GT.0) THEN X DO 600 KK=1,NCC4R X K=K+1 X SUM=SUM+AMAT(I,K)*BVEC(K) X600 CONTINUE X END IF X CVEC(I)=SUM X100 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE AVMUL. XC X END END_OF_FILE if test 4104 -ne `wc -c <'avmul.f'`; then echo shar: \"'avmul.f'\" unpacked with wrong size! fi # end of 'avmul.f' fi if test -f 'bakdif.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'bakdif.f'\" else echo shar: Extracting \"'bakdif.f'\" \(510 characters\) sed "s/^X//" >'bakdif.f' <<'END_OF_FILE' X SUBROUTINE BAKDIF(OVERFL,J ,N ,DELTAJ,TEMPJ ,FVEC , X $ FVECJ1,JACFDM,XC ,FVECEV) XC XC FEB. 6, 1991 XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JACFDM(N,N) X DIMENSION FVEC(N),FVECJ1(N),XC(N) X LOGICAL OVERFL X EXTERNAL FVECEV X DELTAJ=TEMPJ-XC(J) X CALL FVECEV(OVERFL,N,FVECJ1,XC) X IF(.NOT.OVERFL) THEN X DO 100 I=1,N X JACFDM(I,J)=(FVEC(I)-FVECJ1(I))/DELTAJ X100 CONTINUE X END IF X RETURN XC XC LAST CARD OF SUBROUTINE BAKDIF. XC X END END_OF_FILE if test 510 -ne `wc -c <'bakdif.f'`; then echo shar: \"'bakdif.f'\" unpacked with wrong size! fi # end of 'bakdif.f' fi if test -f 'bnddif.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'bnddif.f'\" else echo shar: Extracting \"'bnddif.f'\" \(1208 characters\) sed "s/^X//" >'bnddif.f' <<'END_OF_FILE' X SUBROUTINE BNDDIF(OVERFL,J ,N ,EPSMCH,BOUNDL,BOUNDU, X $ FVECC ,FVECJ1,JACFDM,WV3 ,XC ,FVECEV) XC XC FEB. 15, 1991 XC XC FINITE DIFFERENCE CALCULATION WHEN THE BOUNDS FOR COMPONENT J XC ARE SO CLOSE THAT NEITHER A FORWARD NOR BACKWARD DIFFERENCE XC CAN BE PERFORMED. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JACFDM(N,N) X DIMENSION BOUNDL(N) ,BOUNDU(N), FVECC(N) ,FVECJ1(N) , X $ WV3(N) ,XC(N) X LOGICAL OVERFL X EXTERNAL FVECEV X DATA ZERO /0.0D0/ XC X EPS3Q=EPSMCH**0.75 XC XC STORE CURRENT X DO 100 I=1,N X WV3(I)=FVECC(I) X100 CONTINUE X XC(J)=BOUNDU(J) X CALL FVECEV(OVERFL,N,FVECJ1,XC) X IF(.NOT.OVERFL) THEN X XC(J)=BOUNDL(J) X CALL FVECEV(OVERFL,N,FVECC,XC) X IF(.NOT.OVERFL) THEN X DO 200 I=1,N XC XC ENSURE THAT THE JACOBIAN CALCULATION ISN'T JUST NOISE. XC X IF(FVECJ1(I)-FVECC(I).GT.EPS3Q) THEN X JACFDM(I,J)=(FVECJ1(I)-FVECC(I))/ X $ (BOUNDU(J)-BOUNDL(J)) X ELSE X JACFDM(I,J)=ZERO X END IF X200 CONTINUE X END IF X END IF X DO 300 I=1,N X FVECC(I)=WV3(I) X300 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE BNDDIF. XC X END END_OF_FILE if test 1208 -ne `wc -c <'bnddif.f'`; then echo shar: \"'bnddif.f'\" unpacked with wrong size! fi # end of 'bnddif.f' fi if test -f 'broyfa.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'broyfa.f'\" else echo shar: Extracting \"'broyfa.f'\" \(4367 characters\) sed "s/^X//" >'broyfa.f' <<'END_OF_FILE' X SUBROUTINE BROYFA(OVERCH,OVERFL,SCLFCH,SCLXCH,MAXEXP, X $ N ,NUNIT ,OUTPUT,EPSMCH,A , X $ DELF ,FVEC ,FVECC ,JAC ,RDIAG , X $ S ,SCALEF,SCALEX,T ,W , X $ XC ,XPLUS ) XC XC FEB. 23, 1992 XC XC THE BROYDEN QUASI-NEWTON METHOD IS APPLIED TO THE FACTORED XC FORM OF THE JACOBIAN. XC XC NOTE: T AND W ARE TEMPORARY WORKING VECTORS ONLY. XC XC THE UPDATE OCCURS ONLY IF A SIGNIFICANT CHANGE IN THE XC JACOBIAN WOULD RESULT, I.E., NOT ALL THE VALUES IN VECTOR W XC ARE LESS THAN THE THRESHOLD IN MAGNITUDE. IF THE VECTOR XC W IS ESSENTIALLY ZERO THEN THE LOGICAL VARIABLE SKIPUP XC REMAINS SET AT TRUE. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X INTEGER OUTPUT X DIMENSION A(N,N) ,DELF(N) ,FVEC(N) ,FVECC(N) , X $ RDIAG(N) ,S(N) ,SCALEF(N) ,SCALEX(N), X $ T(N) ,W(N) ,XC(N) ,XPLUS(N) X LOGICAL OVERCH ,OVERFL ,SCLFCH ,SCLXCH , X $ SKIPUP X DATA ZERO,TEN /0.0D0,10.0D0/ XC X OVERFL=.FALSE. X EPS=TEN**(-MAXEXP) X SQRTEP=SQRT(EPSMCH) XC X DO 100 I=1,N X A(I,I)=RDIAG(I) X S(I)=XPLUS(I)-XC(I) X100 CONTINUE XC XC R IS NOW IN THE UPPER TRIANGLE OF A. XC X SKIPUP=.TRUE. XC XC THE BROYDEN UPDATE IS CONDENSED INTO THE FORM XC XC A(NEW) = A(OLD) + T S^ XC XC THE PRODUCT A*S IS FORMED IN TWO STAGES AS R IS IN THE UPPER XC TRIANGLE OF MATRIX A AND Q^ IS IN JAC. XC XC FIRST MULTIPLY R*S (A IS CONSIDERED UPPER TRIANGULAR) XC X CALL UVMUL(N,N,N,N,A,S,T) XC XC NOTE: THIS T IS TEMPORARY - IT IS THE T FROM BELOW WHICH XC IS SENT TO SUBROUTINE QRUPDA. XC X DO 200 I=1,N X CALL INNERP(OVERCH,OVERFL,MAXEXP,N ,N ,N , X $ NUNIT ,OUTPUT,SUM ,JAC(N*(I-1)+1,1),T) X W(I)=SCALEF(I)*(FVEC(I)-FVECC(I))-SUM XC XC TEST TO ENSURE VECTOR W IS NONZERO. ANY VALUE GREATER XC THAN THE THRESHOLD WILL SET SKIPUP TO FALSE. XC X IF(ABS(W(I)).GT.SQRTEP*SCALEF(I)*(ABS(FVEC(I))+ X $ ABS(FVECC(I)))) THEN X SKIPUP=.FALSE. X ELSE X W(I)=ZERO X END IF X200 CONTINUE XC XC IF W(I)=0 FOR ALL I THEN THE UPDATE IS SKIPPED. XC X IF(.NOT.SKIPUP) THEN XC XC T=Q^W; Q^ IS IN JAC. XC X CALL AVMUL(N,N,N,N,JAC,W,T) X IF(SCLXCH) THEN X DO 300 I=1,N X W(I)=S(I)*SCALEX(I) X300 CONTINUE X ELSE X CALL MATCOP(N,N,1,1,N,1,S,W) X END IF X CALL TWONRM(OVERFL,MAXEXP,N,EPSMCH,DENOM,W) XC XC IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN. XC X IF(OVERFL.OR.LOG10(DENOM+EPS).GT.MAXEXP/2) THEN X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,2) X2 FORMAT(T3,'*',4X,'WARNING: JACOBIAN NOT UPDATED', X $ ' TO AVOID OVERFLOW IN DENOMINATOR OF',T74,'*') X WRITE(NUNIT,3) X3 FORMAT(T3,'*',4X,'BROYDEN UPDATE',T74,'*') X END IF X RETURN X ELSE X DENOM=DENOM*DENOM X END IF XC XC IF DENOM IS ZERO AVOID DIVIDE BY ZERO AND CONTINUE WITH XC SAME JACOBIAN. XC X IF(DENOM.EQ.ZERO) RETURN XC XC THE SCALED VERSION OF S REPLACES THE ORIGINAL BEFORE XC BEING SENT TO QRUPDA. XC X DO 400 I=1,N X S(I)=S(I)*SCALEX(I)*SCALEX(I)/DENOM X400 CONTINUE XC XC UPDATE THE QR DECOMPOSITION USING A SERIES OF GIVENS XC ROTATIONS. XC X CALL QRUPDA(OVERFL,MAXEXP,N,EPSMCH,A,JAC,T,S) XC XC RESET RDIAG AS DIAGONAL OF CURRENT R WHICH IS IN XC THE UPPER TRIANGE OF A. XC X DO 500 I=1,N X RDIAG(I)=A(I,I) X500 CONTINUE X END IF XC XC UPDATE THE GRADIENT VECTOR, DELF. THE NEW Q^ IS IN JAC. XC XC DELF = (QR)^F = R^Q^F = R^JAC F XC X IF(SCLFCH) THEN X DO 600 I=1,N X W(I)=FVEC(I)*SCALEF(I) X600 CONTINUE X ELSE X CALL MATCOP(N,N,1,1,N,1,FVEC,W) X END IF X CALL AVMUL(N,N,N,N,JAC,W,T) X CALL UTBMUL(N,N,1,1,N,N,A,T,DELF) X RETURN XC XC LAST CARD OF SUBROUTINE BROYFA. XC X END END_OF_FILE if test 4367 -ne `wc -c <'broyfa.f'`; then echo shar: \"'broyfa.f'\" unpacked with wrong size! fi # end of 'broyfa.f' fi if test -f 'broyun.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'broyun.f'\" else echo shar: Extracting \"'broyun.f'\" \(2115 characters\) sed "s/^X//" >'broyun.f' <<'END_OF_FILE' X SUBROUTINE BROYUN(OVERFL,MAXEXP,N ,NUNIT ,OUTPUT, X $ EPSMCH,FVEC ,FVECC ,JAC ,SCALEX, X $ WV1 ,XC ,XPLUS) XC XC FEB. 23, 1992 XC XC UPDATE THE JACOBIAN USING BROYDEN'S METHOD. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X INTEGER OUTPUT X DIMENSION FVEC(N) ,FVECC(N) ,SCALEX(N),WV1(N) , X $ XC(N) ,XPLUS(N) X LOGICAL OVERFL X DATA ZERO,TEN /0.0D0,10.0D0/ XC X EPS=TEN**(-MAXEXP) X SQRTEP=SQRT(EPSMCH) XC X DO 100 I=1,N X WV1(I)=(XPLUS(I)-XC(I))*SCALEX(I) X100 CONTINUE X CALL TWONRM(OVERFL,MAXEXP,N,EPSMCH,DENOM,WV1) XC XC IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN. XC X IF(OVERFL.OR.LOG10(DENOM+EPS).GT.MAXEXP/2) THEN X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,2) X2 FORMAT(T3,'*',4X,'WARNING: JACOBIAN NOT UPDATED', X $ ' TO AVOID OVERFLOW IN DENOMINATOR OF',T74,'*') X WRITE(NUNIT,3) X3 FORMAT(T3,'*',4X,'BROYDEN UPDATE',T74,'*') X END IF X RETURN X ELSE X DENOM=DENOM*DENOM X END IF XC XC IF DENOM IS ZERO, AVOID OVERFLOW, CONTINUE WITH SAME JACOBIAN. XC X IF(DENOM.EQ.ZERO) RETURN XC XC UPDATE JACOBIAN BY ROWS. XC X DO 200 I=1,N X SUM=ZERO X DO 300 J=1,N X SUM=SUM+JAC(I,J)*(XPLUS(J)-XC(J)) X300 CONTINUE X TEMPI=FVEC(I)-FVECC(I)-SUM XC XC CHECK TO ENSURE THAT SOME MEANINGFUL CHANGE IS BEING MADE XC TO THE APPROXIMATE JACOBIAN; IF NOT, SKIP UPDATING ROW I. XC X IF(ABS(TEMPI).GE.SQRTEP*(ABS(FVEC(I))+ABS(FVECC(I)))) X $ THEN X TEMPI=TEMPI/DENOM X DO 400 J=1,N X JAC(I,J)=JAC(I,J)+TEMPI* X $ (XPLUS(J)-XC(J))*SCALEX(J)*SCALEX(J) X400 CONTINUE X END IF X200 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE BROYUN. XC X END END_OF_FILE if test 2115 -ne `wc -c <'broyun.f'`; then echo shar: \"'broyun.f'\" unpacked with wrong size! fi # end of 'broyun.f' fi if test -f 'cholde.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'cholde.f'\" else echo shar: Extracting \"'cholde.f'\" \(1982 characters\) sed "s/^X//" >'cholde.f' <<'END_OF_FILE' X SUBROUTINE CHOLDE(N,MAXADD,MAXFFL,SQRTEP,H,L) XC XC FEB. 23, 1992 XC XC THIS SUBROUTINE FINDS THE CHOLESKY DECOMPOSITION OF THE XC MATRIX, H, AND RETURNS IT IN THE LOWER TRIANGLE OF XC MATRIX, L. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION L(N,N) ,MAXADD ,MAXFFL ,MINL ,MINL2 , X $ MINLJJ X DIMENSION H(N,N) X DATA ZERO /0.0D0/ XC X MINL=SQRT(SQRTEP)*MAXFFL XC XC MAXFFL EQUALS 0 WHEN THE MATRIX IS KNOWN TO BE POSITIVE XC DEFINITE. XC X IF(MAXFFL.EQ.ZERO) THEN XC XC FIND SQUARE ROOT OF LARGEST MAGNITUDE DIAGONAL ELEMENT XC AND SET MINL2. XC X DO 100 I=1,N X MAXFFL=MAX(MAXFFL,ABS(H(I,I))) X100 CONTINUE X MAXFFL=SQRT(MAXFFL) X MINL2=SQRTEP*MAXFFL X END IF XC XC MAXADD CONTAINS THE MAXIMUM AMOUNT WHICH IS IMPLICITLY ADDED XC TO ANY DIAGONAL ELEMENT OF MATRIX H. XC X MAXADD=ZERO X DO 200 J=1,N X SUM=ZERO X DO 300 I=1,J-1 X SUM=SUM+L(J,I)*L(J,I) X300 CONTINUE X L(J,J)=H(J,J)-SUM X MINLJJ=ZERO X DO 400 I=J+1,N X SUM=ZERO X DO 500 K=1,J-1 X SUM=SUM+L(I,K)*L(J,K) X500 CONTINUE X L(I,J)=H(J,I)-SUM X MINLJJ=MAX(MINLJJ,ABS(L(I,J))) X400 CONTINUE X MINLJJ=MAX(MINLJJ/MAXFFL,MINL) X IF(L(J,J).GT.MINLJJ*MINLJJ) THEN XC XC NORMAL CHOLESKY DECOMPOSITION. XC X L(J,J)=SQRT(L(J,J)) X ELSE XC XC IMPLICITLY PERTURB DIAGONAL OF H. XC X IF(MINLJJ.LT.MINL2) THEN X MINLJJ=MINL2 X END IF X MAXADD=MAX(MAXADD,MINLJJ*MINLJJ-L(J,J)) X L(J,J)=MINLJJ X END IF X DO 600 I=J+1,N X L(I,J)=L(I,J)/L(J,J) X600 CONTINUE X200 CONTINUE X RETURN XC XC LAST CARD OF SUBROUTINE CHOLDE. XC X END END_OF_FILE if test 1982 -ne `wc -c <'cholde.f'`; then echo shar: \"'cholde.f'\" unpacked with wrong size! fi # end of 'cholde.f' fi if test -f 'chsolv.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chsolv.f'\" else echo shar: Extracting \"'chsolv.f'\" \(692 characters\) sed "s/^X//" >'chsolv.f' <<'END_OF_FILE' X SUBROUTINE CHSOLV(OVERCH,OVERFL,MAXEXP,N ,NUNIT , X $ OUTPUT,L ,RHS ,S ,WV2 ) XC XC FEB. 14, 1991 XC XC THIS SUBROUTINE USES FORWARD/BACKWARD SUBSTITUTION TO SOLVE THE XC SYSTEM OF LINEAR EQUATIONS: XC XC (LL^)S=RHS XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION L(N,N) X INTEGER OUTPUT X DIMENSION RHS(N) ,S(N) ,WV2(N) X LOGICAL OVERCH ,OVERFL XC X CALL LSOLV(OVERCH,OVERFL,MAXEXP,N,NUNIT,OUTPUT,L,WV2,RHS) X CALL LTSOLV(OVERCH,OVERFL,MAXEXP,N,NUNIT,OUTPUT,L,S,WV2) XC X RETURN XC XC LAST CARD OF SUBROUTINE CHSOLV. XC X END END_OF_FILE if test 692 -ne `wc -c <'chsolv.f'`; then echo shar: \"'chsolv.f'\" unpacked with wrong size! fi # end of 'chsolv.f' fi if test -f 'condno.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'condno.f'\" else echo shar: Extracting \"'condno.f'\" \(4495 characters\) sed "s/^X//" >'condno.f' <<'END_OF_FILE' X SUBROUTINE CONDNO(OVERCH,OVERFL,MAXEXP,N ,NUNIT , X $ OUTPUT,CONNUM,A ,P ,PM , X $ Q ,RDIAG ) XC XC FEB. 14, 1991 XC XC THIS SUBROUTINE ESTIMATES THE CONDITION NUMBER OF A XC QR-DECOMPOSED MATRIX USING THE METHOD OF CLINE, MOLER, XC STEWART AND WILKINSON (SIAM J. N.A. 16 P368 (1979) ). XC XC IF A POTENTIAL OVERFLOW IS DETECTED AT ANY POINT THEN A XC CONDITION NUMBER EQUIVALENT TO THAT OF A SINGULAR MATRIX XC IS ASSIGNED BY THE CALLING SUBROUTINE. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION A(N,N) ,P(N) ,PM(N) ,RDIAG(N) ,Q(N) X LOGICAL OVERCH ,OVERFL X DATA ZERO,ONE,TEN /0.0D0,1.0D0,10.0D0/ XC X OVERFL=.FALSE. X EPS=TEN**(-MAXEXP) XC X CONNUM=ABS(RDIAG(1)) X DO 100 J=2,N X TEMP=ZERO X DO 200 I=1,J-1 X IF(OVERCH) THEN X IF(ABS(A(I,J)).GT.TEN**(MAXEXP-1)) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X TEMP=TEMP+ABS(A(I,J)) X200 CONTINUE X TEMP=TEMP+ABS(RDIAG(J)) X CONNUM=MAX(CONNUM,TEMP) X100 CONTINUE X Q(1)=ONE/RDIAG(1) X DO 300 I=2,N X IF(OVERCH) THEN X IF(LOG10(ABS(Q(1))+EPS)+LOG10(ABS(A(1,I))+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X P(I)=A(1,I)*Q(1) X300 CONTINUE X DO 400 J=2,N X IF(OVERCH) THEN X IF(LOG10(ABS(P(J))+EPS)-LOG10(ABS(RDIAG(J))+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X QP=(ONE-P(J))/RDIAG(J) X QM=(-ONE-P(J))/RDIAG(J) X TEMP=ABS(QP) X TEMPM=ABS(QM) X DO 500 I=J+1,N X IF(OVERCH) THEN X IF(LOG10(ABS(A(J,I))+EPS)+LOG10(ABS(QM)+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X PM(I)=P(I)+A(J,I)*QM X IF(OVERCH) THEN X IF(LOG10(ABS(PM(I))+EPS)-LOG10(ABS(RDIAG(I))+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X TEMPM=TEMPM+(ABS(PM(I))/ABS(RDIAG(I))) X IF(OVERCH) THEN X IF(TEMPM.GT.TEN**(MAXEXP-1)) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X IF(OVERCH) THEN X IF(LOG10(ABS(A(J,I))+EPS)+LOG10(ABS(QP)+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X P(I)=P(I)+A(J,I)*QP X IF(OVERCH) THEN X IF(LOG10(ABS(P(I))+EPS)-LOG10(ABS(RDIAG(I))+EPS) X $ .GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X TEMP=TEMP+(ABS(P(I))/ABS(RDIAG(I))) X IF(OVERCH) THEN X IF(TEMP.GT.TEN**(MAXEXP-1)) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X500 CONTINUE X IF(TEMP.GE.TEMPM) THEN X Q(J)=QP X ELSE X Q(J)=QM X DO 600 I=J+1,N X P(I)=PM(I) X600 CONTINUE X END IF X400 CONTINUE X QNORM=ZERO X DO 700 J=1,N X QNORM=QNORM+ABS(Q(J)) X IF(OVERCH) THEN X IF(QNORM.GT.TEN**(MAXEXP-1)) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X700 CONTINUE X IF(LOG10(CONNUM)-LOG10(QNORM).GT.MAXEXP) THEN X OVERFL=.TRUE. X RETURN X END IF X CONNUM=CONNUM/QNORM X CALL RSOLV(OVERCH,OVERFL,MAXEXP,N,NUNIT,OUTPUT,A,RDIAG,Q) X IF(OVERFL) RETURN X QNORM=ZERO X DO 800 J=1,N X QNORM=QNORM+ABS(Q(J)) X IF(OVERCH) THEN X IF(QNORM.GT.TEN**(MAXEXP-1)) THEN X OVERFL=.TRUE. X RETURN X END IF X END IF X800 CONTINUE X CONNUM=CONNUM*QNORM X RETURN XC XC LAST CARD OF SUBROUTINE CONDNO. XC X END END_OF_FILE if test 4495 -ne `wc -c <'condno.f'`; then echo shar: \"'condno.f'\" unpacked with wrong size! fi # end of 'condno.f' fi if test -f 'd1mach.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'d1mach.f'\" else echo shar: Extracting \"'d1mach.f'\" \(7515 characters\) sed "s/^X//" >'d1mach.f' <<'END_OF_FILE' X DOUBLE PRECISION FUNCTION D1MACH(I) X INTEGER I XC XC DOUBLE-PRECISION MACHINE CONSTANTS XC D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. XC D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. XC D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. XC D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. XC D1MACH( 5) = LOG10(B) XC X INTEGER SMALL(2) X INTEGER LARGE(2) X INTEGER RIGHT(2) X INTEGER DIVER(2) X INTEGER LOG10(2) X INTEGER SC, CRAY1(38), J X COMMON /D9MACH/ CRAY1 X SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC X DOUBLE PRECISION DMACH(5) X EQUIVALENCE (DMACH(1),SMALL(1)) X EQUIVALENCE (DMACH(2),LARGE(1)) X EQUIVALENCE (DMACH(3),RIGHT(1)) X EQUIVALENCE (DMACH(4),DIVER(1)) X EQUIVALENCE (DMACH(5),LOG10(1)) XC THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES. XC R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF XC D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR XC MANY MACHINES YET. XC TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 XC ON THE NEXT LINE X DATA SC/0/ XC AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. XC CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY XC mail netlib@research.bell-labs.com XC send old1mach from blas XC PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. XC XC MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. XC DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / XC DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / XC DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / XC DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / XC DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/ XC XC MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING XC 32-BIT INTEGERS. XC DATA SMALL(1),SMALL(2) / 8388608, 0 / XC DATA LARGE(1),LARGE(2) / 2147483647, -1 / XC DATA RIGHT(1),RIGHT(2) / 612368384, 0 / XC DATA DIVER(1),DIVER(2) / 620756992, 0 / XC DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ XC XC MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. XC DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / XC DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / XC DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / XC DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / XC DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/ XC XC ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES. X IF (SC .NE. 987) THEN X DMACH(1) = 1.D13 X IF ( SMALL(1) .EQ. 1117925532 X * .AND. SMALL(2) .EQ. -448790528) THEN X* *** IEEE BIG ENDIAN *** X SMALL(1) = 1048576 X SMALL(2) = 0 X LARGE(1) = 2146435071 X LARGE(2) = -1 X RIGHT(1) = 1017118720 X RIGHT(2) = 0 X DIVER(1) = 1018167296 X DIVER(2) = 0 X LOG10(1) = 1070810131 X LOG10(2) = 1352628735 X ELSE IF ( SMALL(2) .EQ. 1117925532 X * .AND. SMALL(1) .EQ. -448790528) THEN X* *** IEEE LITTLE ENDIAN *** X SMALL(2) = 1048576 X SMALL(1) = 0 X LARGE(2) = 2146435071 X LARGE(1) = -1 X RIGHT(2) = 1017118720 X RIGHT(1) = 0 X DIVER(2) = 1018167296 X DIVER(1) = 0 X LOG10(2) = 1070810131 X LOG10(1) = 1352628735 X ELSE IF ( SMALL(1) .EQ. -2065213935 X * .AND. SMALL(2) .EQ. 10752) THEN X* *** VAX WITH D_FLOATING *** X SMALL(1) = 128 X SMALL(2) = 0 X LARGE(1) = -32769 X LARGE(2) = -1 X RIGHT(1) = 9344 X RIGHT(2) = 0 X DIVER(1) = 9472 X DIVER(2) = 0 X LOG10(1) = 546979738 X LOG10(2) = -805796613 X ELSE IF ( SMALL(1) .EQ. 1267827943 X * .AND. SMALL(2) .EQ. 704643072) THEN X* *** IBM MAINFRAME *** X SMALL(1) = 1048576 X SMALL(2) = 0 X LARGE(1) = 2147483647 X LARGE(2) = -1 X RIGHT(1) = 856686592 X RIGHT(2) = 0 X DIVER(1) = 873463808 X DIVER(2) = 0 X LOG10(1) = 1091781651 X LOG10(2) = 1352628735 X ELSE IF ( SMALL(1) .EQ. 1120022684 X * .AND. SMALL(2) .EQ. -448790528) THEN X* *** CONVEX C-1 *** X SMALL(1) = 1048576 X SMALL(2) = 0 X LARGE(1) = 2147483647 X LARGE(2) = -1 X RIGHT(1) = 1019215872 X RIGHT(2) = 0 X DIVER(1) = 1020264448 X DIVER(2) = 0 X LOG10(1) = 1072907283 X LOG10(2) = 1352628735 X ELSE IF ( SMALL(1) .EQ. 815547074 X * .AND. SMALL(2) .EQ. 58688) THEN X* *** VAX G-FLOATING *** X SMALL(1) = 16 X SMALL(2) = 0 X LARGE(1) = -32769 X LARGE(2) = -1 X RIGHT(1) = 15552 X RIGHT(2) = 0 X DIVER(1) = 15568 X DIVER(2) = 0 X LOG10(1) = 1142112243 X LOG10(2) = 2046775455 X ELSE X DMACH(2) = 1.D27 + 1 X DMACH(3) = 1.D27 X LARGE(2) = LARGE(2) - RIGHT(2) X IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN X CRAY1(1) = 67291416 X DO 10 J = 1, 20 X CRAY1(J+1) = CRAY1(J) + CRAY1(J) X 10 CONTINUE X CRAY1(22) = CRAY1(21) + 321322 X DO 20 J = 22, 37 X CRAY1(J+1) = CRAY1(J) + CRAY1(J) X 20 CONTINUE X IF (CRAY1(38) .EQ. SMALL(1)) THEN X* *** CRAY *** X CALL I1MCRY(SMALL(1), J, 8285, 8388608, 0) X SMALL(2) = 0 X CALL I1MCRY(LARGE(1), J, 24574, 16777215, 16777215) X CALL I1MCRY(LARGE(2), J, 0, 16777215, 16777214) X CALL I1MCRY(RIGHT(1), J, 16291, 8388608, 0) X RIGHT(2) = 0 X CALL I1MCRY(DIVER(1), J, 16292, 8388608, 0) X DIVER(2) = 0 X CALL I1MCRY(LOG10(1), J, 16383, 10100890, 8715215) X CALL I1MCRY(LOG10(2), J, 0, 16226447, 9001388) X ELSE X WRITE(*,9000) X STOP 779 X END IF X ELSE X WRITE(*,9000) X STOP 779 X END IF X END IF X SC = 987 X END IF X* SANITY CHECK X IF (DMACH(4) .GE. 1.0D0) STOP 778 X IF (I .LT. 1 .OR. I .GT. 5) THEN X WRITE(*,*) 'D1MACH(I): I =',I,' is out of bounds.' X STOP X END IF X D1MACH = DMACH(I) X RETURN X 9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/ X *' appropriate for your machine.') X* /* Standard C source for D1MACH -- remove the * in column 1 */ X*#include X*#include X*#include X*double d1mach_(long *i) X*{ X* switch(*i){ X* case 1: return DBL_MIN; X* case 2: return DBL_MAX; X* case 3: return DBL_EPSILON/FLT_RADIX; X* case 4: return DBL_EPSILON; X* case 5: return log10(FLT_RADIX); X* } X* fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i); X* exit(1); return 0; /* some compilers demand return values */ X*} X END X SUBROUTINE I1MCRY(A, A1, B, C, D) X**** SPECIAL COMPUTATION FOR OLD CRAY MACHINES **** X INTEGER A, A1, B, C, D X A1 = 16777216*B + C X A = 16777216*A1 + D X END END_OF_FILE if test 7515 -ne `wc -c <'d1mach.f'`; then echo shar: \"'d1mach.f'\" unpacked with wrong size! fi # end of 'd1mach.f' fi if test -f 'delcau.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'delcau.f'\" else echo shar: Extracting \"'delcau.f'\" \(4979 characters\) sed "s/^X//" >'delcau.f' <<'END_OF_FILE' X SUBROUTINE DELCAU(CAUCHY,OVERCH,OVERFL,ITNUM ,MAXEXP, X $ N ,NUNIT ,OUTPUT,BETA ,CAULEN, X $ DELTA ,EPSMCH,MAXSTP,NEWLEN,SQRTZ , X $ A ,DELF ,SCALEX,WV1 ) XC XC FEB. 23, 1992 XC XC THIS SUBROUTINE ESTABLISHES AN INITIAL TRUST REGION, DELTA, XC IF ONE IS NOT SPECIFIED BY THE USER AND FINDS THE LENGTH OF XC THE SCALED CAUCHY STEP, CAULEN, AT EACH STEP IF THE DOUBLE XC DOGLEG OPTION IS BEING USED. XC XC THE USER HAS TWO CHOICES FOR THE INITIAL TRUST REGION: XC 1) BASED ON THE SCALED CAUCHY STEP XC 2) BASED ON THE SCALED NEWTON STEP XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION MAXSTP ,NEWLEN X INTEGER OUTPUT X DIMENSION A(N,N) ,DELF(N) ,SCALEX(N) ,WV1(N) X LOGICAL CAUCHY ,OVERCH ,OVERFL X DATA ZERO,THREE,TEN/0.0D0,3.0D0,10.0D0/ XC X OVERFL=.FALSE. X EPS=TEN**(-MAXEXP) XC XC IF DELTA IS NOT GIVEN EVALUATE IT USING EITHER THE CAUCHY XC STEP OR THE NEWTON STEP AS SPECIFIED BY THE USER. XC XC THE SCALED CAUCHY LENGTH, CAULEN, IS REQUIRED IN TWO CASES. XC 1) WHEN SELECTED AS THE CRITERION FOR THE INITIAL DELTA XC 2) IN THE DOUBLE DOGLEG STEP REGARDLESS OF (1) XC X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,2) X2 FORMAT(T3,'*',4X,'DETERMINATION OF SCALED CAUCHY', X $ ' STEP LENGTH, CAULEN',T74,'*') X END IF XC XC FIND FACTOR WHICH GIVES CAUCHY POINT WHEN MULTPLYING XC STEEPEST DESCENT DIRECTION, DELF. XC XC CAULEN= ZETA**1.5/BETA XC = SQRTZ**3/BETA XC X DO 100 I=1,N X WV1(I)=DELF(I)/SCALEX(I) X100 CONTINUE X CALL TWONRM(OVERFL,MAXEXP,N,EPSMCH,SQRTZ,WV1) X IF(OVERFL) THEN X CAULEN=TEN**MAXEXP X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,3) SQRTZ X3 FORMAT(T3,'*',7X,'ZETA SET TO ',1PD11.3,' TO' X $ ' AVOID OVERFLOW',T74,'*') X WRITE(NUNIT,4) CAULEN X4 FORMAT(T3,'*',7X,'SCALED CAUCHY LENGTH, CAULEN ', X $ 'SET TO',1PD9.2,' TO AVOID OVERFLOW',T74,'*') X IF(ITNUM.EQ.1) THEN X WRITE(NUNIT,5) X5 FORMAT(T3,'*',7X,'THE PROBLEM SHOULD BE RESCALED', X $ ' OR A NEW STARTING POINT CHOSEN',T74,'*') X WRITE(NUNIT,6) X6 FORMAT(T3,'*',7X,'EXECUTION CONTINUES WITH', X $ ' SUBSTITUTIONS AS LISTED ABOVE',T74,'*') X END IF X END IF X ELSE X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,7) SQRTZ X7 FORMAT(T3,'*',7X,'SQUARE ROOT OF ZETA, SQRTZ: ', X $ 1PD12.3,T74,'*') X END IF X END IF XC XC NOTE: THE LOWER TRIANGLE OF MATRIX A NOW CONTAINS THE XC TRANSPOSE OF R WHERE A=QR. XC X BETA=ZERO X DO 200 I=1,N X TEMP=ZERO X DO 300 J=I,N X IF(OVERCH) THEN X IF(LOG10(ABS(A(J,I))+EPS)+LOG10(ABS(DELF(J))+EPS). X $ GT.MAXEXP) THEN X CAULEN=SQRT(EPSMCH) X GO TO 301 X END IF X END IF X TEMP=TEMP+A(J,I)*DELF(J)/(SCALEX(J)*SCALEX(J)) X300 CONTINUE X BETA=BETA+TEMP*TEMP X200 CONTINUE X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,8) BETA X8 FORMAT(T3,'*',7X,'BETA: ',1PD11.3,6X,'NOTE: ', X $ 'CAULEN=ZETA**1.5/BETA',T74,'*') X WRITE(NUNIT,1) X END IF XC XC AVOID OVERFLOWS IN FINDING CAULEN. XC X IF(THREE*LOG10(SQRTZ+EPS)-LOG10(BETA+EPS).LT.MAXEXP X $ .AND.(.NOT.OVERFL).AND.BETA.NE.ZERO) THEN XC XC NORMAL DETERMINATION. XC X CAULEN=SQRTZ*SQRTZ*SQRTZ/BETA XC XC THIS STEP AVOIDS DIVIDE BY ZERO IN DOGLEG IN THE XC (RARE) CASE WHERE DELF(I)=0 FOR ALL I BUT THE XC POINT IS NOT YET A SOLUTION - MOST LIKELY A BAD XC STARTING ESTIMATE. XC X CAULEN=MAX(CAULEN,TEN**(-MAXEXP)) XC X ELSE XC XC SUBSTITUTION TO AVOID OVERFLOW. XC X CAULEN=TEN**MAXEXP X END IF X301 CONTINUE X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,9) CAULEN X9 FORMAT(T3,'*',7X,'SCALED CAUCHY LENGTH, CAULEN: ', X $ 1PD12.3,T74,'*') X END IF XC XC ESTABLISH INITIAL TRUST REGION IF NEEDED. XC X IF(DELTA.LT.ZERO) THEN XC XC USE DISTANCE TO CAUCHY POINT OR LENGTH OF NEWTON STEP. XC X IF(CAUCHY) THEN X DELTA=MIN(CAULEN,MAXSTP) X ELSE X DELTA=MIN(NEWLEN,MAXSTP) X END IF X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,10) DELTA X10 FORMAT(T3,'*',7X,'INITIAL TRUST REGION SIZE, DELTA: ', X $ 1PD12.3,T74,'*') X END IF X END IF X RETURN XC XC LAST CARD OF SUBROUTINE DELCAU. XC X END END_OF_FILE if test 4979 -ne `wc -c <'delcau.f'`; then echo shar: \"'delcau.f'\" unpacked with wrong size! fi # end of 'delcau.f' fi if test -f 'deufls.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'deufls.f'\" else echo shar: Extracting \"'deufls.f'\" \(23007 characters\) sed "s/^X//" >'deufls.f' <<'END_OF_FILE' X SUBROUTINE DEUFLS(ABORT ,DEUFLH,GEOMS ,OVERCH,OVERFL, X $ QNFAIL,QRSING,RESTRT,SCLFCH,SCLXCH, X $ ACPCOD,ACPTCR,CONTYP,ITNUM ,JUPDM , X $ MAXEXP,MAXLIN,N ,NFUNC ,NUNIT , X $ OUTPUT,QNUPDM,STOPCR,ALPHA ,CONFAC, X $ DELFTS,EPSMCH,FCNMAX,FCNNEW,FCNOLD, X $ LAMBDA,NEWMAX,SBRNRM,SIGMA ,A , X $ ASTORE,BOUNDL,BOUNDU,DELF ,FVEC , X $ HHPI ,JAC ,RDIAG ,RHS ,S , X $ SBAR ,SCALEF,SCALEX,SN ,WV2 , X $ XC ,XPLUS ,FVECEV) XC XC FEB. 23, 1992 XC XC THIS SUBROUTINE CONDUCTS A LINE SEARCH IN THE NEWTON XC DIRECTION IF NO CONSTRAINTS ARE VIOLATED. IF THE FIRST XC TRIAL IS A FAILURE THERE ARE TWO TYPES OF LINE SEARCH XC AVAILABLE. XC 1) REDUCE THE RELAXATION FACTOR, LAMBDA, TO XC SIGMA*LAMBDA WHERE SIGMA IS USER-SPECIFIED XC (GEOMETRIC LINE SEARCH) XC 2) AT THE FIRST STEP MINIMIZE A QUADRATIC THROUGH XC THE OBJECTIVE FUNCTION AT THE CURRENT POINT AND XC TRIAL ESTIMATE (WHICH MUST BE A FAILURE) WITH XC INITIAL SLOPE DELFTS. AT SUBSEQUENT STEPS MINI- XC MIZE A CUBIC THROUGH THE OBJECTIVE FUNCTION AT XC THE TWO MOST RECENT FAILURES AND THE CURRENT XC POINT, AGAIN USING THE INITIAL SLOPE, DELFTS. XC XC CONVIO INDICATES A CONSTRAINT VIOLATION BY ONE OR MORE XC COMPONENTS XC FRSTST INDICATES THE FIRST STEP IN THE LINE SEARCH. XC XC RATIO RATIO OF PROPOSED STEP LENGTH IN (I)TH DIRECTION XC TO DISTANCE FROM (I)TH COMPONENT TO BOUNDARY XC VIOLATED XC RATIOM MINIMUM OF RATIOS FOR ALL CONSTAINTS VIOLATED XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) ,LAMBDA ,LAMPRE ,LAMTMP , X $ NEWMAX X INTEGER ACPCOD ,ACPTCR ,CONTYP ,OUTPUT , X $ STOPCR ,QNUPDM X DIMENSION A(N,N) ,ASTORE(N,N),BOUNDL(N),BOUNDU(N), X $ FVEC(N) ,HHPI(N) ,RDIAG(N) ,RHS(N) , X $ S(N) ,DELF(N) ,SBAR(N) ,SCALEF(N), X $ SCALEX(N),SN(N) ,WV2(N) ,XC(N) , X $ XPLUS(N) X LOGICAL ABORT ,CONVIO ,DEUFLH ,FRSTST , X $ GEOMS ,OVERCH ,OVERFL ,QNFAIL , X $ QRSING ,RESTRT ,SCLFCH ,SCLXCH , X $ WRNSUP X COMMON/NNES_2/WRNSUP X EXTERNAL FVECEV X DATA ZERO,POINT1,POINT5,ONE,TWO,THREE,TEN /0.0D0,0.1D0,0.5D0, X $ 1.0D0,2.0D0,3.0D0,10.0D0/ XC X FRSTST=.TRUE. X OVERFL=.FALSE. X EPS=TEN**(-MAXEXP) X DO 100 K=1,MAXLIN X RATIOM=ONE XC X CONVIO=.FALSE. XC XC FIND TRIAL POINT AND CHECK IF CONSTRAINTS VIOLATED (IF XC CONTYP IS NOT EQUAL TO ZERO). XC X DO 200 I=1,N XC XC NOTE: WV2 MARKS VIOLATIONS. WV2(I) CHANGES TO XC 1 FOR LOWER BOUND VIOLATIONS AND TO 2 FOR XC FOR UPPER BOUND VIOLATIONS. CONSTRAINT VIOL- XC ATIONS CAN ONLY OCCUR AT THE FIRST STEP. XC X WV2(I)=-ONE X XPLUS(I)=XC(I)+LAMBDA*SN(I) X IF(CONTYP.GT.0.AND.FRSTST) THEN X IF(XPLUS(I).LT.BOUNDL(I)) THEN X CONVIO=.TRUE. X WV2(I)=ONE X ELSEIF(XPLUS(I).GT.BOUNDU(I)) THEN X CONVIO=.TRUE. X WV2(I)=TWO X ELSE X WV2(I)=-ONE X END IF X END IF X200 CONTINUE XC XC IF CONSTRAINTS ARE VIOLATED FIRST REDUCE THE STEP XC SIZES FOR THE VIOLATING COMPONENTS TO OBTAIN A XC FEASIBLE POINT. IF THE DIRECTION TO THIS MODIFIED XC POINT IS NOT A DESCENT DIRECTION OR IF THE MODIFIED XC STEP DOES NOT LEAD TO AN ACCEPTABLE POINT THEN RETURN XC TO THE NEWTON DIRECTION AND START A LINE SEARCH AT A XC FEASIBLE POINT WHERE THE COMPONENT WHICH HAS THE XC SMALLEST VALUE OF RATIO (DEFINED BELOW) IS TAKEN TO XC "CONFAC" OF THE DISTANCE TO THE BOUNDARY. DEFAULT XC VALUE OF CONFAC IS 0.95. XC X IF(CONVIO) THEN X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,2) K X2 FORMAT(T3,'*',10X,'LINE SEARCH STEP:',I3,T74,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,3) LAMBDA X3 FORMAT(T3,'*',10X,'LAMBDA FOR ATTEMPTED STEP: ', X $ 1PD12.3,T74,'*') X WRITE(NUNIT,4) X4 FORMAT(T3,'*',10X,'CONSTRAINT VIOLATED',T74,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,5) X5 FORMAT(T3,'*',10X,'TRIAL ESTIMATES (VIOLATIONS', X $ ' MARKED)',T74,'*') X WRITE(NUNIT,1) X DO 300 I=1,N X IF(WV2(I).GT.ZERO) THEN X WRITE(NUNIT,6) I,XPLUS(I) X6 FORMAT(T3,'*',13X,'XPLUS(',I3,') = ',1PD12.3, X $ 2X,'*',T74,'*') X ELSE X WRITE(NUNIT,7) I,XPLUS(I) X7 FORMAT(T3,'*',13X,'XPLUS(',I3,') = ',1PD12.3, X $ T74,'*') X END IF X300 CONTINUE X END IF X DO 400 I=1,N X IF(WV2(I).GT.ZERO) THEN XC XC FIND RATIO FOR THIS VIOLATING COMPONENT. XC X IF(WV2(I).EQ.ONE) THEN X RATIO=-(XC(I)-BOUNDL(I))/ X $ (XPLUS(I)-XC(I)) X ELSEIF(WV2(I).EQ.TWO) THEN X RATIO=(BOUNDU(I)-XC(I))/ X $ (XPLUS(I)-XC(I)) X END IF XC XC NOTE: THIS LINE IS FOR OUTPUT PURPOSES ONLY. XC X WV2(I)=RATIO XC X RATIOM=MIN(RATIOM,RATIO) X IF(RATIO.GT.POINT5) THEN X S(I)=CONFAC*RATIO*LAMBDA*SN(I) X ELSE XC XC WITHIN BUFFER ZONE - ONLY TAKE 1/2 XC OF THE STEP YOU WOULD TAKE OTHERWISE. XC X S(I)=POINT5*CONFAC*RATIO*LAMBDA*SN(I) X END IF XC XC ESTABLISH MODIFIED TRIAL POINT. XC X XPLUS(I)=XC(I)+S(I) X ELSE XC XC FOR NONVIOLATORS XPLUS REMAINS UNCHANGED BUT XC THE COMPONENT OF S IS LOADED TO CHECK THE XC DIRECTIONAL DERIVATIVE. XC X S(I)=LAMBDA*SN(I) X END IF X400 CONTINUE X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,8) X8 FORMAT(T3,'*',7X,'NEW S AND XPLUS VECTORS', X $ ' (WITH RATIOS FOR VIOLATIONS)',T74,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,9) X9 FORMAT(T3,'*',7X,'NOTE: RATIOS ARE RATIO OF', X $ ' LENGTH TO BOUNDARY FROM CURRENT',T74,'*') X WRITE(NUNIT,10) X10 FORMAT(T3,'*',13X,'X VECTOR TO MAGNITUDE OF', X $ ' CORRESPONDING PROPOSED STEP',T74,'*') X WRITE(NUNIT,1) X DO 500 I=1,N X IF(WV2(I).LT.ZERO) THEN X WRITE(NUNIT,11) I,S(I),I,XPLUS(I) X11 FORMAT(T3,'*',7X,'S(',I3,') = ',1PD12.3,4X, X $ 'XPLUS(',I3,') = ',1PD12.3,T74,'*') X ELSE X WRITE(NUNIT,12) I,S(I),I,XPLUS(I),WV2(I) X12 FORMAT(T3,'*',7X,'S(',I3,') = ',1PD12.3,4X, X $ 'XPLUS(',I3,') = ',1PD12.3,1X,1PD11.3,T74,'*') X END IF X500 CONTINUE X WRITE(NUNIT,1) X WRITE(NUNIT,13) RATIOM X13 FORMAT(T3,'*',7X,'MINIMUM OF RATIOS, RATIOM: ', X $ 1PD12.3,T74,'*') X END IF XC XC CHECK DIRECTIONAL DERIVATIVE FOR MODIFIED POINT, DLFTSM. XC X CALL INNERP(OVERCH,OVERFL,MAXEXP,N ,N ,N , X $ NUNIT ,OUTPUT,DLFTSM,DELF ,S ) X OVERFL=.FALSE. X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,14) DLFTSM X14 FORMAT(T3,'*',7X,'INNER PRODUCT OF DELF AND S FOR', X $ ' MODIFIED S: ',1PD12.3,T74,'*') X END IF XC XC IF INNER PRODUCT IS POSITIVE RETURN TO NEWTON DIRECTION XC AND CONDUCT A LINE SEARCH WITHIN THE FEASIBLE REGION. XC X IF(DLFTSM.GT.ZERO) THEN X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,15) X15 FORMAT(T3,'*',7X,'DELFTS > 0',' START LINE', X $ ' SEARCH AT LAMBDA=CONFAC*LAMBDA*RATIOM',T74,'*') X WRITE(NUNIT,16) X16 FORMAT(T3,'*',7X,'NOTE: NO TRIAL POINT WAS', X $ ' EVALUATED AT THIS STEP OF LINE SEARCH',T74,'*') X END IF XC XC THE STARTING POINT IS SET AT JUST INSIDE THE MOST XC VIOLATED BOUNDARY. XC X LAMBDA=CONFAC*RATIOM*LAMBDA XC XC LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE. XC X GO TO 101 X END IF XC X END IF XC XC NO CONSTRAINTS VIOLATED - EVALUATE RESIDUAL VECTOR XC AT NEW POINT. XC X CALL FVECEV(OVERFL,N,FVEC,XPLUS) X NFUNC=NFUNC+1 XC XC CHECK FOR OVERFLOW IN FUNCTION VECTOR EVALUATION. XC IF SO, REDUCE STEP LENGTH AND CONTINUE LINE SEARCH. XC X IF(OVERFL) THEN XC X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,17) X17 FORMAT(T3,'*',7X,'OVERFLOW IN FUNCTION VECTOR', X $ ' - STEP LENGTH REDUCED',T74,'*') XC XC FORCE STEP TO BE WITHIN CONSTRAINTS - DON'T CALL XC THIS THE FIRST STEP, I.E. FRSTST STAYS AT TRUE. XC X END IF X IF(CONVIO) THEN X LAMBDA=RATIOM*CONFAC*LAMBDA X ELSE X LAMBDA=SIGMA*LAMBDA X END IF XC XC LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE. XC X GO TO 101 X END IF XC XC EVALUATE (POSSIBLY SCALED) OBJECTIVE FUNCTION AT NEW POINT. XC X CALL FCNEVL(OVERFL,MAXEXP,N ,NUNIT ,OUTPUT, X $ EPSMCH,FCNNEW,FVEC ,SCALEF,WV2 ) X IF(OVERFL) THEN XC X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,18) X18 FORMAT(T3,'*',7X,'OVERFLOW IN OBJECTIVE FUNCTION', X $ ' - STEP LENGTH REDUCED',T74,'*') XC XC FORCE STEP TO BE WITHIN CONSTRAINTS - DON'T CALL XC THIS THE FIRST STEP, I.E. FRSTST STAYS AT TRUE. XC X END IF X IF(CONVIO) THEN X LAMBDA=RATIOM*CONFAC*LAMBDA X ELSE X LAMBDA=SIGMA*LAMBDA X END IF X GO TO 101 X END IF XC XC IF DEUFLHARD'S METHOD IS BEING USED FOR EITHER XC RELAXATION FACTOR INITIALIZATION OR THE SECOND XC ACCEPTANCE CRITERION THEN EVALUATE SBAR. EVALU- XC ATION METHOD DEPENDS UPON WHETHER THE JACOBIAN XC WAS PERTURBED IN THE SOLUTION OF THE LINEAR SYSTEM. XC LOGICAL VARIABLE QRSING IS TRUE IF PERTURBATION XC TOOK PLACE. XC X IF(DEUFLH.OR.ACPTCR.EQ.12) THEN X IF(QRSING) THEN XC XC FORM -J^F AS RIGHT HAND SIDE - METHOD DEPENDS ON XC WHETHER QNUPDM EQUALS 0 OR 1 IF A QUASI-NEWTON XC UPDATE IS BEING USED. IF JUPDM IS 0 THEN THE NEWTON XC STEP HAS BEEN FOUND IN SUBROUTINE NSTPUN. XC X IF(JUPDM.EQ.0.OR.QNUPDM.EQ.0) THEN XC XC UNSCALED JACOBIAN IN MATRIX JAC. XC X DO 600 I=1,N X IF(SCLFCH) THEN X WV2(I)=-FVEC(I)*SCALEF(I)*SCALEF(I) X ELSE X WV2(I)=-FVEC(I) X END IF X600 CONTINUE X CALL ATBMUL(N,N,1,1,N,N,JAC,WV2,RHS) X ELSE XC XC R IN UPPER TRIANGLE OF A PLUS RDIAG AND Q^ IN JAC XC - FROM QR DECOMPOSITION OF SCALED JACOBIAN. XC X DO 700 I=1,N X WV2(I)=ZERO X DO 800 J=1,N X WV2(I)=WV2(I)-JAC(I,J)*FVEC(J)*SCALEF(J) X800 CONTINUE X700 CONTINUE X RHS(1)=RDIAG(1)*WV2(1) X DO 900 J=2,N X RHS(J)=ZERO X DO 1000 I=1,J-1 X RHS(J)=RHS(J)+A(I,J)*WV2(I) X1000 CONTINUE X RHS(J)=RHS(J)+RDIAG(J)*WV2(J) X900 CONTINUE X END IF X CALL CHSOLV(OVERCH,OVERFL,MAXEXP,N ,NUNIT , X $ OUTPUT,A ,RHS ,SBAR ,WV2 ) X ELSE XC XC RIGHT HAND SIDE IS -FVEC. XC X IF(QNUPDM.EQ.0.OR.JUPDM.EQ.0) THEN XC XC QR DECOMPOSITION OF SCALED JACOBIAN STORED IN XC ASTORE. XC X DO 1100 I=1,N X SBAR(I)=-FVEC(I)*SCALEF(I) X1100 CONTINUE X CALL QRSOLV(OVERCH,OVERFL,MAXEXP,N ,NUNIT , X $ OUTPUT,ASTORE,HHPI ,RDIAG ,SBAR ) X ELSE XC XC SET UP RIGHT HAND SIDE - MULTIPLY -FVEC BY Q^ XC (STORED IN JAC). XC X DO 1200 I=1,N X WV2(I)=-FVEC(I)*SCALEF(I) X1200 CONTINUE X CALL AVMUL(N,N,N,N,JAC,WV2,SBAR) X CALL RSOLV(OVERCH,OVERFL,MAXEXP,N ,NUNIT , X $ OUTPUT,A ,RDIAG ,SBAR ) X END IF X END IF XC XC NORM OF SCALED SBAR IS NEEDED FOR SECOND ACCEPTANCE TEST. XC X IF(ACPTCR.EQ.12) THEN X DO 1300 I=1,N X WV2(I)=SCALEX(I)*SBAR(I) X1300 CONTINUE X CALL TWONRM(OVERFL,MAXEXP,N,EPSMCH,SBRNRM,WV2) X END IF X END IF XC X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) K X WRITE(NUNIT,1) X WRITE(NUNIT,3) LAMBDA X WRITE(NUNIT,1) X IF(.NOT.CONVIO) THEN X WRITE(NUNIT,19) X19 FORMAT(T3,'*',10X,'NEW COMPONENT/FCN VECTORS', X $ ' (XPLUS(I)=XC(I)+LAMBDA*SN(I))',T74,'*') X ELSE X WRITE(NUNIT,20) X20 FORMAT(T3,'*',10X,'NEW FUNCTION VECTORS', X $ ' AT MODIFIED POINT',T74,'*') X END IF X WRITE(NUNIT,1) X DO 1400 I=1,N X WRITE(NUNIT,21) I,XPLUS(I),I,FVEC(I) X21 FORMAT(T3,'*',10X,'XPLUS(',I3,') = ',1PD12.3, X $ 5X,'FVEC(',I3,') = ',1PD12.3,T74,'*') X1400 CONTINUE X WRITE(NUNIT,1) X IF(.NOT.SCLFCH) THEN X WRITE(NUNIT,22) FCNNEW X22 FORMAT(T3,'*',10X,'OBJECTIVE FUNCTION VALUE', X $ ' AT XPLUS: .........'1PD12.3,T74,'*') X ELSE X WRITE(10,23) FCNNEW X23 FORMAT(T3,'*',10X,'SCALED OBJECTIVE FUNCTION VALUE', X $ ' AT XPLUS: ..'1PD12.3,T74,'*') X END IF X WRITE(NUNIT,24) FCNMAX+ALPHA*LAMBDA*DELFTS X24 FORMAT(T3,'*',10X,'FCNMAX+ALPHA*LAMBDA*DELFTS:', X $ ' ................',1PD12.3,T74,'*') X IF(DEUFLH.OR.ACPTCR.EQ.12) THEN X IF(ITNUM .GT.0) THEN X IF(.NOT.SCLXCH) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,25) X25 FORMAT(T3,'*',10X,'DEUFLHARD SBAR VECTOR', X $ T74,'*') X WRITE(NUNIT,1) X DO 1500 I=1,N X WRITE(NUNIT,26) I,SBAR(I) X26 FORMAT(T3,'*',10X,'SBAR(',I3,') = ', X $ 1PD12.3,T74,'*') X1500 CONTINUE X ELSE X WRITE(NUNIT,1) X WRITE(NUNIT,27) X27 FORMAT(T3,'*',10X,'DEUFLHARD SBAR VECTOR', X $ 14X,'IN SCALED X UNITS',T74,'*') X WRITE(NUNIT,1) X DO 1600 I=1,N X WRITE(NUNIT,28) I,SBAR(I),I,SCALEX(I)*SBAR(I) X28 FORMAT(T3,'*',10X,'SBAR(',I3,') = ',1PD12.3, X $ 8X,'SBAR(',I3,') = ',1PD12.3,T74,'*') X1600 CONTINUE X END IF X END IF X END IF X IF(ACPTCR.EQ.12) THEN X WRITE(NUNIT,1) X IF(.NOT.SCLXCH) THEN X WRITE(NUNIT,29) SBRNRM X29 FORMAT(T3,'*',10X,'VALUE OF SBRNRM', X $ ' AT XPLUS: ..................'1PD12.3,T74,'*') X ELSE X WRITE(NUNIT,30) SBRNRM X30 FORMAT(T3,'*',10X,'VALUE OF SCALED SBRNRM', X $ ' AT XPLUS: ...........'1PD12.3,T74,'*') X END IF X WRITE(NUNIT,31) NEWMAX X31 FORMAT(T3,'*',10X,'NEWMAX:',' ..............', X $ '......................',1PD12.3,T74,'*') X END IF X END IF XC XC CHECK FOR ACCEPTABLE STEP. XC X IF(FCNNEW.LT.FCNMAX+ALPHA*LAMBDA*DELFTS) THEN X ACPCOD=1 XC XC NOTE: STEP WILL BE ACCEPTED REGARDLESS OF NEXT TEST. XC THIS SECTION IS FOR BOOKKEEPING ONLY. XC X IF(ACPTCR.EQ.12) THEN X IF(SBRNRM.LT.NEWMAX) THEN X ACPCOD=12 X END IF X END IF XC X RETURN X END IF X IF(ACPTCR.EQ.12.AND.SBRNRM.LT.NEWMAX) THEN X ACPCOD=2 X RETURN X END IF XC XC FAILURE OF STEP ACCEPTANCE TEST. XC X IF(CONVIO) THEN X LAMBDA=CONFAC*RATIOM*LAMBDA XC XC LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE. XC X GO TO 101 X END IF X IF(LAMBDA.EQ.ZERO) THEN X IF(OUTPUT.GT.0) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,32) X32 FORMAT(T3,'*',7X,'LAMBDA IS 0.0: NO PROGRESS', X $ ' POSSIBLE - CHECK BOUNDS OR START',T74,'*') X END IF X ABORT=.TRUE. X RETURN X END IF X IF(GEOMS) THEN XC XC GEOMETRIC LINE SEARCH XC X LAMBDA=SIGMA*LAMBDA XC X ELSE XC X IF(FRSTST) THEN X FRSTST=.FALSE. XC XC FIND MINIMUM OF QUADRATIC AT FIRST STEP. XC X LAMTMP=-(LAMBDA*LAMBDA)*DELFTS/ X $ (TWO*(FCNNEW-FCNOLD-LAMBDA*DELFTS)) X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,33) LAMTMP X33 FORMAT(T3,'*',13X,'TEMPORARY LAMBDA FROM', X $ ' QUADRATIC MODEL: ',1PD11.3,T74,'*') X END IF XC X ELSE XC XC FIND MINIMUM OF CUBIC AT SUBSEQUENT STEPS. XC X FACTOR=ONE/(LAMBDA-LAMPRE) X IF(LAMBDA*LAMBDA.EQ.ZERO) THEN X LAMBDA=SIGMA*LAMBDA XC XC NOTE: IF THIS LAMBDA**2 WAS ZERO ANY SUBSEQUENT XC LAMBDA**2 WILL ALSO BE ZERO. XC X GO TO 101 X END IF X ACUBIC=FACTOR*((ONE/LAMBDA*LAMBDA)*(FCNNEW-FCNOLD- X $ LAMBDA*DELFTS)-((ONE/LAMPRE*LAMBDA)*(FPLPRE- X $ FCNOLD-LAMPRE*DELFTS))) X BCUBIC=FACTOR*((-LAMPRE/LAMBDA*LAMBDA)*(FCNNEW-FCNOLD- X $ LAMBDA*DELFTS)+((LAMBDA/LAMPRE*LAMBDA)*(FPLPRE- X $ FCNOLD-LAMPRE*DELFTS))) X IF(TWO*LOG10(ABS(BCUBIC)+EPS).GT.DBLE(MAXEXP)) X $ THEN X LAMTMP=SIGMA*LAMBDA X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,34) X34 FORMAT(T3,'*',13X,'POTENTIAL OVERFLOW IN' X $ ' CALCULATING TRIAL LAMBDA FROM',T74,'*') X WRITE(NUNIT,35) X35 FORMAT(T3,'*',13X,'CUBIC MODEL - LAMBDA', X $ ' SET TO SIGMA*LAMBDA',T74,'*') X END IF X ELSE X DISC=BCUBIC*BCUBIC-THREE*ACUBIC*DELFTS X IF(ABS(ACUBIC).LE.EPSMCH) THEN X LAMTMP=-DELFTS/(TWO*BCUBIC) X ELSE X IF(DISC.LT.ZERO) THEN X LAMTMP=SIGMA*LAMBDA X ELSE X LAMTMP=(-BCUBIC+SQRT(DISC))/(THREE*ACUBIC) X END IF X END IF X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,36) LAMTMP X36 FORMAT(T3,'*',13X,'TEMPORARY LAMBDA FROM', X $ ' CUBIC MODEL: .....',1PD11.3,T74,'*') X END IF X END IF X IF(LAMTMP.GT.SIGMA*LAMBDA) THEN X LAMTMP=SIGMA*LAMBDA X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,37) X37 FORMAT(T3,'*',13X,'LAMTMP TOO LARGE - REDUCED', X $ ' TO SIGMA*LAMBDA',T74,'*') X END IF X END IF X END IF X LAMPRE=LAMBDA X FPLPRE=FCNNEW X IF(LAMTMP.LT.POINT1*LAMBDA) THEN X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,38) X38 FORMAT(T3,'*',13X,'LAMTMP TOO SMALL - INCREASED', X $ ' TO 0.1*LAMBDA',T74,'*') X END IF X LAMBDA=POINT1*LAMBDA X ELSE X IF(OUTPUT.GT.4.AND.LAMTMP.NE.SIGMA*LAMBDA) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,39) X39 FORMAT(T3,'*',13X,'LAMTMP WITHIN LIMITS - ', X $ 'LAMBDA SET TO LAMTMP',T74,'*') X END IF X LAMBDA=LAMTMP X END IF X END IF X101 CONTINUE X100 CONTINUE XC XC FAILURE IN LINE SEARCH XC X ACPCOD=0 XC XC IF A QUASI-NEWTON STEP HAS FAILED IN THE LINE SEARCH THEN XC SET QNFAIL TO TRUE ANS RETURN TO SUBROUTINE NNES. THIS WILL XC CAUSE THE JACOBIAN TO BE RE-EVALUATED EXPLICITLY AND A LINE XC SEARCH IN THE NEW DIRECTION CONDUCTED. XC X IF(.NOT.RESTRT) THEN X QNFAIL=.TRUE. X RETURN X END IF XC XC FALL THROUGH MAIN LOOP - WARNING GIVEN. XC X IF(OUTPUT.GT.2.AND.(.NOT.WRNSUP)) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,40) MAXLIN X40 FORMAT(T3,'*',7X,'WARNING: ',I3,' CYCLES COMPLETED', X $ ' IN LINE SEARCH WITHOUT SUCCESS',T74,'*') X END IF X IF(STOPCR.EQ.2.OR.STOPCR.EQ.3) THEN X STOPCR=12 X WRITE(NUNIT,1) X WRITE(NUNIT,41) X41 FORMAT(T3,'*',7X,'STOPPING CRITERION RESET FROM ', X $ '2 TO 12 TO AVOID HANGING',T74,'*') X END IF X RETURN XC XC LAST CARD OF SUBROUTINE DEUFLS. XC X END END_OF_FILE if test 23013 -ne `wc -c <'deufls.f'`; then echo shar: \"'deufls.f'\" unpacked with wrong size! fi # end of 'deufls.f' fi if test -f 'dogleg.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dogleg.f'\" else echo shar: Extracting \"'dogleg.f'\" \(7772 characters\) sed "s/^X//" >'dogleg.f' <<'END_OF_FILE' X SUBROUTINE DOGLEG(FRSTDG,NEWTKN,OVERCH,OVERFL,MAXEXP, X $ N ,NOTRST,NUNIT ,OUTPUT,BETA , X $ CAULEN,DELTA ,ETAFAC,NEWLEN,SQRTZ , X $ DELF ,S ,SCALEX,SN ,SSDHAT, X $ VHAT ) XC XC FEB. 23, 1992 XC XC THIS SUBROUTINE FINDS A TRUST REGION STEP USING THE XC (DOUBLE) DOGLEG METHOD. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION LAMBDA ,NEWLEN X INTEGER OUTPUT X DIMENSION DELF(N) ,S(N) ,SCALEX(N),SN(N) , X $ SSDHAT(N),VHAT(N) X LOGICAL FRSTDG ,NEWTKN ,OVERCH ,OVERFL , X $ WRNSUP X COMMON/NNES_2/WRNSUP X DATA ZERO,ONE /0.0D0,1.0D0/ XC X OVERFL=.FALSE. XC X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',70X,'*') X WRITE(NUNIT,2) NOTRST,DELTA X2 FORMAT(T3,'*',4X,'TRUST REGION STEP:',I6,2X, X $ 'TRUST REGION LENGTH, DELTA:',1PD11.3,2X,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,3) NEWLEN X3 FORMAT(T3,'*',7X,'LENGTH OF NEWTON STEP, NEWLEN: ', X $ 1PD11.3,21X,'*') X END IF XC XC CHECK FOR NEWTON STEP WITHIN TRUST REGION - IF SO USE XC NEWTON STEP. XC X IF(NEWLEN.LE.DELTA) THEN X DO 100 I=1,N X S(I)=SN(I) X100 CONTINUE X NEWTKN=.TRUE. X TEMP=DELTA X DELTA=NEWLEN X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,4) X4 FORMAT(T3,'*',7X,'NEWTON STEP WITHIN ACCEPTABLE RANGE', X $ ' ( <= THAN DELTA)',11X,'*') X IF(TEMP.EQ.DELTA) THEN X WRITE(NUNIT,5) DELTA X5 FORMAT(T3,'*',7X,'DELTA STAYS AT LENGTH OF NEWTON', X $ ' STEP: ',1PD11.3,14X,'*') X ELSE X WRITE(NUNIT,6) X6 FORMAT(T3,'*',7X,'DELTA SET TO LENGTH OF NEWTON', X $ ' STEP',29X,'*') X END IF X WRITE(NUNIT,1) X WRITE(NUNIT,7) X7 FORMAT(T3,'*',7X,'FULL NEWTON STEP ATTEMPTED',37X,'*') X END IF X RETURN X ELSE XC XC NEWTON STEP NOT WITHIN TRUST REGION - APPLY (DOUBLE) XC DOGLEG PROCEDURE. (IF ETAFAC EQUALS 1.0 THEN THE SINGLE XC DOGLEG PROCEDURE IS BEING APPLIED). XC X NEWTKN=.FALSE. X IF(FRSTDG) THEN XC XC SPECIAL SECTION FOR FIRST DOGLEG STEP - CALCULATES XC CAUCHY POINT (MINIMIZER OF MODEL FUNCTION IN STEEPEST XC DESCENT DIRECTION OF THE OBJECTIVE FUNCTION). XC X FRSTDG=.FALSE. X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X IF(ETAFAC.EQ.ONE) THEN X WRITE(NUNIT,8) X8 FORMAT(T3,'*',7X,'FIRST SINGLE DOGLEG STEP', X $ 39X,'*') X ELSE X WRITE(NUNIT,9) X9 FORMAT(T3,'*',7X,'FIRST DOUBLE DOGLEG STEP', X $ 39X,'*') X END IF X WRITE(NUNIT,1) X WRITE(NUNIT,10) X10 FORMAT(T3,'*',10X,'SCALED CAUCHY STEP',42X,'*') X WRITE(NUNIT,1) X END IF XC XC NOTE: BETA AND SQRTZ WERE CALCULATED IN SUBROUTINE DELCAU. XC X ZETA=SQRTZ*SQRTZ XC XC FIND STEP TO CAUCHY POINT. XC X FACTOR=-(ZETA/BETA) X DO 200 I=1,N X SSDHAT(I)=FACTOR*(DELF(I)/SCALEX(I)) X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,11) I,SSDHAT(I) X11 FORMAT(T3,'*',13X,'SSDHAT(',I3,') = ',1PD12.3, X $ 31X,'*') X END IF X200 CONTINUE X CALL INNERP(OVERCH,OVERFL,MAXEXP,N,N,N ,NUNIT , X $ OUTPUT,DELFTS,DELF ,SN ) X OVERFL=.FALSE. XC XC PROTECT AGAINST (RARE) CASE WHEN CALCULATED DIRECTIONAL XC DERIVATIVE EQUALS ZERO. XC X IF(DELFTS.NE.ZERO) THEN XC XC STANDARD EXECUTION. XC X GAMMA=(ZETA/ABS(DELFTS))*(ZETA/BETA) X ETA=ETAFAC+(ONE-ETAFAC)*GAMMA X ELSE X IF(OUTPUT.GT.1.AND.(.NOT.WRNSUP)) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,12) X12 FORMAT(T3,'*',4X,'WARNING: DELFTS=0; ETA SET', X $ ' TO 1.0 TO AVOID DIVISION BY ZERO',8X,'*') X END IF X ETA=ONE X END IF X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,13) ETA X13 FORMAT(T3,'*',10X,'ETA = ',1PD11.3,43X,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,14) X14 FORMAT(T3,'*',10X,'VHAT VECTOR VHAT(I)=ETA*', X $ 'SN(I)*SCALEX(I)-SSDHAT(I)',7X,'*') X WRITE(NUNIT,1) X END IF X DO 300 I=1,N X VHAT(I)=ETA*SCALEX(I)*SN(I)-SSDHAT(I) X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,15) I,VHAT(I) X15 FORMAT(T3,'*',13X,'VHAT(',I3,') = ',1PD12.3, X $ 33X,'*') X END IF X300 CONTINUE X END IF XC XC ETA*NEWLEN <= DELTA MEANS TAKE STEP IN NEWTON DIRECTION XC TO TRUST REGION BOUNDARY. XC X IF(ETA*NEWLEN.LE.DELTA) THEN X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,16) X16 FORMAT(T3,'*',10X,'ETA*NEWLEN <= DELTA S(I)', X $ '= (DELTA/NEWLEN)*SN(I)',10X,'*') X END IF X DO 400 I=1,N X S(I)=(DELTA/NEWLEN)*SN(I) X400 CONTINUE X ELSE XC XC DISTANCE TO CAUCHY POINT >= DELTA MEANS TAKE STEP IN XC STEEPEST DESCENT DIRECTION TO TRUST REGION BOUNDARY. XC X IF(CAULEN.GE.DELTA) THEN X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,17) X17 FORMAT(T3,'*',10X,'CAULEN >= DELTA S(I)', X $ '=(DELTA/CAULEN)*(SSDHAT(I)/SCALEX(I))',1X,'*') X END IF X DO 500 I=1,N X S(I)=(DELTA/CAULEN)*(SSDHAT(I)/SCALEX(I)) X500 CONTINUE X ELSE XC XC TAKE (DOUBLE) DOGLEG STEP. XC X CALL INNERP(OVERCH,OVERFL,MAXEXP,N,N,N ,NUNIT , X $ OUTPUT,TEMP ,SSDHAT,VHAT ) X CALL INNERP(OVERCH,OVERFL,MAXEXP,N,N,N ,NUNIT , X $ OUTPUT,TEMPV ,VHAT ,VHAT ) X OVERFL=.FALSE. X LAMBDA=(-TEMP+SQRT(TEMP*TEMP-TEMPV*(CAULEN*CAULEN X $ -DELTA*DELTA)))/TEMPV X IF(OUTPUT.GT.4) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,18) X18 FORMAT(T3,'*',10X,'S(I)=(SSDHAT(I)+LAMBDA*VHAT(I))', X $ '/SCALEX(I)',19X,'*') X WRITE(NUNIT,1) X WRITE(NUNIT,19) LAMBDA X19 FORMAT(T3,'*',10X,'WHERE LAMBDA = ',1PD12.3,33X,'*') X END IF X DO 600 I=1,N X S(I)=(SSDHAT(I)+LAMBDA*VHAT(I))/SCALEX(I) X600 CONTINUE X END IF X END IF X IF(OUTPUT.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,1) X WRITE(NUNIT,20) X20 FORMAT(T3,'*',10X,'REVISED STEP FROM SUBROUTINE', X $ ' DOGLEG',25X,'*') X WRITE(NUNIT,1) X DO 700 I=1,N X WRITE(NUNIT,21) I,S(I) X21 FORMAT(T3,'*',13X,'S(',I3,') = ',1PD12.3,36X,'*') X700 CONTINUE X END IF X END IF X RETURN X END END_OF_FILE if test 7772 -ne `wc -c <'dogleg.f'`; then echo shar: \"'dogleg.f'\" unpacked with wrong size! fi # end of 'dogleg.f' fi if test -f 'ex1.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ex1.f'\" else echo shar: Extracting \"'ex1.f'\" \(3556 characters\) sed "s/^X//" >'ex1.f' <<'END_OF_FILE' X PROGRAM EX1 XC XC ROSENBROCK'S BANANA XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X PARAMETER ( MGLL= 10, X $ N= 2, X $ NUNIT= 10) X DOUBLE PRECISION X $ JAC(N,N) ,LAM0 ,MSTPF ,NSTTOL X INTEGER ACPTCR ,OUTPUT ,QNUPDM ,STOPCR , X $ SUPPRS ,TRMCOD ,TRUPDM X DIMENSION A(N,N) ,BOUNDL(N) ,BOUNDU(N) ,DELF(N) , X $ FSAVE(N) ,FTRACK(0:MGLL-1) ,FVEC(N) , X $ FVECC(N) ,H(N,N) ,HHPI(N) ,PLEE(N,N), X $ RDIAG(N) ,S(N) ,SBAR(N) ,SCALEF(N), X $ SCALEX(N) ,SN(N) ,SSDHAT(N) , X $ STRACK(0:MGLL-1) ,VHAT(N) ,WV1(N) , X $ WV2(N) ,WV3(N) ,WV4(N) ,XC(N) , X $ XPLUS(N) ,XSAVE(N) X LOGICAL ABSNEW ,BYPASS ,CAUCHY ,DEUFLH , X $ GEOMS ,LINESR ,MATSUP ,NEWTON , X $ OVERCH ,WRNSUP X CHARACTER*6 HELP X EXTERNAL FCN,JACOB X COMMON/NNES_1/MATSUP X COMMON/NNES_2/WRNSUP X COMMON/NNES_3/BYPASS X COMMON/NNES_4/NFETOT X OPEN(UNIT=NUNIT,FILE='EX1.OUT',STATUS='UNKNOWN') X XC(1)=-1.2D0 X XC(2)=1.0D0 X CALL SETUP(ABSNEW ,CAUCHY ,DEUFLH ,GEOMS ,LINESR , X $ NEWTON ,OVERCH ,ACPTCR ,ITSCLF ,ITSCLX , X $ JACTYP ,JUPDM ,MAXEXP ,MAXIT ,MAXNS , X $ MAXQNS ,MINQNS ,N ,NARMIJ ,NIEJEV , X $ NJACCH ,OUTPUT ,QNUPDM ,STOPCR ,SUPPRS , X $ TRUPDM ,ALPHA ,CONFAC ,DELTA ,DELFAC , X $ EPSMCH ,ETAFAC ,FDTOLJ ,FTOL ,LAM0 , X $ MSTPF ,NSTTOL ,OMEGA ,RATIOF ,SIGMA , X $ STPTOL ,BOUNDL ,BOUNDU ,SCALEF ,SCALEX , X $ HELP) X CALL NNES(ABSNEW ,CAUCHY ,DEUFLH ,GEOMS ,LINESR , X $ NEWTON ,OVERCH ,ACPTCR ,ITSCLF ,ITSCLX , X $ JACTYP ,JUPDM ,MAXEXP ,MAXIT ,MAXNS , X $ MAXQNS ,MGLL ,MINQNS ,N ,NARMIJ , X $ NIEJEV ,NJACCH ,NJETOT ,NUNIT ,OUTPUT , X $ QNUPDM ,STOPCR ,SUPPRS ,TRMCOD ,TRUPDM , X $ ALPHA ,CONFAC ,DELTA ,DELFAC ,EPSMCH , X $ ETAFAC ,FCNNEW ,FDTOLJ ,FTOL ,LAM0 , X $ MSTPF ,NSTTOL ,OMEGA ,RATIOF ,SIGMA , X $ STPTOL ,A ,BOUNDL ,BOUNDU ,DELF , X $ FSAVE ,FTRACK ,FVEC ,FVECC ,H , X $ HHPI ,JAC ,PLEE ,RDIAG ,S , X $ SBAR ,SCALEF ,SCALEX ,SN ,SSDHAT , X $ STRACK ,VHAT ,WV1 ,WV2 ,WV3 , X $ WV4 ,XC ,XPLUS ,XSAVE ,HELP , X $ FCN ,JACOB ) X STOP X END X SUBROUTINE FCN(OVERFL,N,FVEC,XC) X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION FVEC(N),XC(N) X COMMON/NNES_4/NFETOT X LOGICAL OVERFL X OVERFL=.FALSE. X NFETOT=NFETOT+1 Xc Xc Rosenbrock Banana Function Xc Xc f1 = 10(x2 - x1**2) Xc f2 = 1 - x1 Xc Xc start: (-1.2,1) Xc (6.39,-0.221) Xc Xc sol'n: (1,1) Xc X FVEC(1)=10.0D0*(XC(2)-XC(1)**2) X FVEC(2)=1.0D0-XC(1) X RETURN X END X SUBROUTINE JACOB(OVERFL,N,JAC,XC) X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X DIMENSION XC(N) X LOGICAL OVERFL X OVERFL=.FALSE. X RETURN X END X END_OF_FILE if test 3556 -ne `wc -c <'ex1.f'`; then echo shar: \"'ex1.f'\" unpacked with wrong size! fi # end of 'ex1.f' fi if test -f 'ex2.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ex2.f'\" else echo shar: Extracting \"'ex2.f'\" \(4486 characters\) sed "s/^X//" >'ex2.f' <<'END_OF_FILE' X PROGRAM EX2 XC XC XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X PARAMETER ( MGLL= 10, X $ N= 6, X $ NUNIT= 10) X DOUBLE PRECISION X $ JAC(N,N) ,LAM0 ,MSTPF ,NSTTOL X INTEGER ACPTCR ,OUTPUT ,QNUPDM ,STOPCR , X $ SUPPRS ,TRMCOD ,TRUPDM X DIMENSION A(N,N) ,BOUNDL(N) ,BOUNDU(N) ,DELF(N) , X $ FSAVE(N) ,FTRACK(0:MGLL-1) ,FVEC(N) , X $ FVECC(N) ,H(N,N) ,HHPI(N) ,PLEE(N,N), X $ RDIAG(N) ,S(N) ,SBAR(N) ,SCALEF(N), X $ SCALEX(N) ,SN(N) ,SSDHAT(N) , X $ STRACK(0:MGLL-1) ,VHAT(N) ,WV1(N) , X $ WV2(N) ,WV3(N) ,WV4(N) ,XC(N) , X $ XPLUS(N) ,XSAVE(N) X LOGICAL ABSNEW ,BYPASS ,CAUCHY ,DEUFLH , X $ GEOMS ,LINESR ,MATSUP ,NEWTON , X $ OVERCH ,WRNSUP X CHARACTER*6 HELP X EXTERNAL FCN,JACOB X COMMON/NNES_1/MATSUP X COMMON/NNES_2/WRNSUP X COMMON/NNES_3/BYPASS X COMMON/NNES_4/NFETOT X OPEN(UNIT=NUNIT,FILE='EX2.OUT',STATUS='UNKNOWN') X DO 100 I=1,N X XC(I)=1.0D0 X100 CONTINUE X CALL SETUP(ABSNEW ,CAUCHY ,DEUFLH ,GEOMS ,LINESR , X $ NEWTON ,OVERCH ,ACPTCR ,ITSCLF ,ITSCLX , X $ JACTYP ,JUPDM ,MAXEXP ,MAXIT ,MAXNS , X $ MAXQNS ,MINQNS ,N ,NARMIJ ,NIEJEV , X $ NJACCH ,OUTPUT ,QNUPDM ,STOPCR ,SUPPRS , X $ TRUPDM ,ALPHA ,CONFAC ,DELTA ,DELFAC , X $ EPSMCH ,ETAFAC ,FDTOLJ ,FTOL ,LAM0 , X $ MSTPF ,NSTTOL ,OMEGA ,RATIOF ,SIGMA , X $ STPTOL ,BOUNDL ,BOUNDU ,SCALEF ,SCALEX , X $ HELP) X OUTPUT=5 X BOUNDL(3)=0.0D0 X BOUNDL(4)=0.0D0 X BOUNDL(6)=0.0D0 XC XC STOPS PREMATURELY UNLESS THIS ARE SET VERY STIFF XC X NSTTOL=1.0D-16 X STPTOL=1.0D-16 XC XC THIS PARAMETER AFFECTS ILL-CONDITIONED JACOBIAN MECHANISM XC XC THIS WILL WORK, =.FALSE. WILL FAIL XC X BYPASS=.TRUE. X CALL NNES(ABSNEW ,CAUCHY ,DEUFLH ,GEOMS ,LINESR , X $ NEWTON ,OVERCH ,ACPTCR ,ITSCLF ,ITSCLX , X $ JACTYP ,JUPDM ,MAXEXP ,MAXIT ,MAXNS , X $ MAXQNS ,MGLL ,MINQNS ,N ,NARMIJ , X $ NIEJEV ,NJACCH ,NJETOT ,NUNIT ,OUTPUT , X $ QNUPDM ,STOPCR ,SUPPRS ,TRMCOD ,TRUPDM , X $ ALPHA ,CONFAC ,DELTA ,DELFAC ,EPSMCH , X $ ETAFAC ,FCNNEW ,FDTOLJ ,FTOL ,LAM0 , X $ MSTPF ,NSTTOL ,OMEGA ,RATIOF ,SIGMA , X $ STPTOL ,A ,BOUNDL ,BOUNDU ,DELF , X $ FSAVE ,FTRACK ,FVEC ,FVECC ,H , X $ HHPI ,JAC ,PLEE ,RDIAG ,S , X $ SBAR ,SCALEF ,SCALEX ,SN ,SSDHAT , X $ STRACK ,VHAT ,WV1 ,WV2 ,WV3 , X $ WV4 ,XC ,XPLUS ,XSAVE ,HELP , X $ FCN ,JACOB ) X STOP X END X SUBROUTINE FCN(OVERFL,N,FVEC,XC) X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DIMENSION FVEC(N),XC(N) X COMMON/NNES_4/NFETOT X LOGICAL OVERFL X OVERFL=.FALSE. X NFETOT=NFETOT+1 Xc Xc Hiebert's 2nd Chemical Engineering Problem Xc Xc source: Hiebert; Sandia Technical Report #SAND80-0181 Xc Sandia National Laboratories, Albuquerque, NM (1980) Xc Xc f1 = X1 + X2 + X4 - .001 Xc f2 = X5 + X6 -55 Xc f3 = X1 + X2 + X3 +2X5 + X6 - 110.001 Xc f4 = X1 - 0.1X2 Xc f5 = X1 - 10000X3X4 Xc f6 = X5 - 5.5e15X3X6 Xc Xc start: (0,0,0,0,0,0) Xc (1,1,1,1,1,1) Xc (1e-4,1e-3,0,1e-4,55,1e-4) Xc (10,10,10,10,10,10,10) Xc Xc sol'n (8.264e-5,8.264e-4,9.091e-5,9.091e-5,55,1.1e-10) Xc X FVEC(1)=XC(1)+XC(2)+XC(4)-0.001D0 X FVEC(2)=XC(5)+XC(6)-55.0D0 X FVEC(3)=XC(1)+XC(2)+XC(3)+2.0D0*XC(5)+XC(6)-110.001D0 X FVEC(4)=XC(1)-0.1D0*XC(2) X FVEC(5)=XC(1)-1.0D04*XC(3)*XC(4) X FVEC(6)=XC(5)-5.5D15*XC(3)*XC(6) X RETURN X END X SUBROUTINE JACOB(OVERFL,N,JAC,XC) X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X DIMENSION XC(N) X LOGICAL OVERFL X OVERFL=.FALSE. X RETURN X END X END_OF_FILE if test 4486 -ne `wc -c <'ex2.f'`; then echo shar: \"'ex2.f'\" unpacked with wrong size! fi # end of 'ex2.f' fi if test -f 'fcnevl.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'fcnevl.f'\" else echo shar: Extracting \"'fcnevl.f'\" \(1240 characters\) sed "s/^X//" >'fcnevl.f' <<'END_OF_FILE' X SUBROUTINE FCNEVL(OVERFL,MAXEXP,N ,NUNIT ,OUTPUT, X $ EPSMCH,FCNNEW,FVEC ,SCALEF,WV1 ) XC XC FEB. 23, 1992 XC XC THE OBJECTIVE FUNCTION, FCNNEW, DEFINED BY: XC XC FCNNEW:=1/2(SCALEF*FVEC^SCALEF*FVEC) XC XC IS EVALUATED BY THIS SUBROUTINE. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X INTEGER OUTPUT X DIMENSION FVEC(N) ,SCALEF(N) ,WV1(N) X LOGICAL OVERFL ,WRNSUP X COMMON/NNES_2/WRNSUP X DATA TWO,TEN /2.0D0,10.0D0/ XC X OVERFL=.FALSE. X EPS=TEN**(-MAXEXP) XC X DO 100 I=1,N X WV1(I)=FVEC(I)*SCALEF(I) X100 CONTINUE X CALL TWONRM(OVERFL,MAXEXP,N,EPSMCH,FCNNEW,WV1) XC XC IF AN OVERFLOW WOULD OCCUR SUBSTITUTE A LARGE VALUE XC FOR FCNNEW. XC X IF(OVERFL.OR.TWO*LOG10(FCNNEW+EPS).GT.MAXEXP) THEN X OVERFL=.TRUE. X FCNNEW=TEN**MAXEXP X IF(OUTPUT.GT.2.AND.(.NOT.WRNSUP)) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,2) FCNNEW X2 FORMAT(T3,'*',4X,'WARNING: TO AVOID OVERFLOW', X $ ' OBJECTIVE FUNCTION SET TO: ',1PD11.3,T74,'*') X END IF X RETURN X END IF X FCNNEW=FCNNEW*FCNNEW/TWO X RETURN X END END_OF_FILE if test 1240 -ne `wc -c <'fcnevl.f'`; then echo shar: \"'fcnevl.f'\" unpacked with wrong size! fi # end of 'fcnevl.f' fi if test -f 'fordif.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'fordif.f'\" else echo shar: Extracting \"'fordif.f'\" \(484 characters\) sed "s/^X//" >'fordif.f' <<'END_OF_FILE' X SUBROUTINE FORDIF(OVERFL,J ,N ,DELTAJ,FVEC , X $ FVECJ1,JACFDM,XC ,FVECEV) XC XC FEB. 6, 1991 XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JACFDM(N,N) X DIMENSION FVEC(N),FVECJ1(N),XC(N) X LOGICAL OVERFL X EXTERNAL FVECEV X CALL FVECEV(OVERFL,N,FVECJ1,XC) X IF(.NOT.OVERFL) THEN X DO 100 I=1,N X JACFDM(I,J)=(FVECJ1(I)-FVEC(I))/DELTAJ X100 CONTINUE X END IF X X RETURN XC XC LAST CARD OF SUBROUTINE FORDIF. XC X END END_OF_FILE if test 484 -ne `wc -c <'fordif.f'`; then echo shar: \"'fordif.f'\" unpacked with wrong size! fi # end of 'fordif.f' fi if test -f 'gradf.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'gradf.f'\" else echo shar: Extracting \"'gradf.f'\" \(2490 characters\) sed "s/^X//" >'gradf.f' <<'END_OF_FILE' X SUBROUTINE GRADF(OVERCH,OVERFL,RESTRT,SCLFCH,SCLXCH,JUPDM , X $ MAXEXP,N ,NUNIT ,OUTPUT,QNUPDM,DELF , X $ FVECC ,JAC ,SCALEF,SCALEX,WV1 ) XC XC FEB. 23, 1992 XC XC THIS SUBROUTINE COMPUTES THE GRADIENT OF THE FUNCTION XC XC F=1/2{SCALEF*FVECC)^(SCALEF*FVECC} XC XC WHICH IS USED AS THE OBJECTIVE FUNCTION FOR MINIMIZATION. XC XC NOTE: WHEN THE FACTORED FORM OF THE JACOBIAN IS UPDATED IN XC QUASI-NEWTON METHODS THE GRADIENT IS UPDATED AS WELL XC IN THE SAME SUBROUTINE - IT IS PRINTED HERE THOUGH. XC IN THESE CASES QNUPDM > 0. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X DOUBLE PRECISION JAC(N,N) X INTEGER OUTPUT ,QNUPDM X DIMENSION DELF(N) ,FVECC(N),SCALEF(N),SCALEX(N),WV1(N) X LOGICAL OVERCH ,OVERFL ,RESTRT ,SCLFCH ,SCLXCH XC X IF(RESTRT.OR.JUPDM.EQ.0.OR.QNUPDM.EQ.0) THEN XC XC GRADIENT NOT ALREADY UPDATED: FIND DELF = J^F. XC X DO 100 I=1,N X IF(SCLFCH) THEN X WV1(I)=FVECC(I)*SCALEF(I)*SCALEF(I) X ELSE X WV1(I)=FVECC(I) X END IF X100 CONTINUE X IF(OVERCH) THEN ! CHECK EACH ENTRY INDIVIDUALLY X CALL ATVOV(OVERFL,MAXEXP,N ,NUNIT , X $ OUTPUT,JAC ,WV1 ,DELF ) X ELSE X CALL ATBMUL(N,N,1,1,N,N,JAC,WV1,DELF) X END IF X END IF XC XC PRINT GRADIENT VECTOR, DELF. XC X IF(OUTPUT.GT.3) THEN X IF(.NOT.SCLXCH) THEN X WRITE(NUNIT,1) X1 FORMAT(T3,'*',T74,'*') X IF(.NOT.SCLFCH) THEN X WRITE(NUNIT,2) X2 FORMAT(T3,'*',4X,'GRADIENT OF OBJECTIVE FUNCTION', X $ T74,'*') X ELSE X WRITE(NUNIT,3) X3 FORMAT(T3,'*',4X,'GRADIENT OF SCALED OBJECTIVE', X $ ' FUNCTION',T74,'*') X END IF X WRITE(NUNIT,1) X DO 300 I=1,N X WRITE(NUNIT,4) I,DELF(I) X4 FORMAT(T3,'*',6X,'DELF(',I3,') = ',1PD12.3,T74,'*') X300 CONTINUE X ELSE X WRITE(NUNIT,1) X WRITE(NUNIT,5) X5 FORMAT(T3,'*',4X,'GRADIENT OF OBJECTIVE FUNCTION',9X, X $ 'IN SCALED X UNITS',T74,'*') X WRITE(NUNIT,1) X DO 400 I=1,N X WRITE(NUNIT,6) I,DELF(I),I,SCALEF(I)*SCALEF(I)* X $ DELF(I)/SCALEX(I) X6 FORMAT(T3,'*',6X,'DELF(',I3,') = ',1PD12.3, X $ 9X,'DELF(',I3,') = ',1PD12.3,T74,'*') X400 CONTINUE X END IF X END IF X RETURN XC XC LAST CARD OF SUBROUTINE GRADF. XC X END END_OF_FILE if test 2490 -ne `wc -c <'gradf.f'`; then echo shar: \"'gradf.f'\" unpacked with wrong size! fi # end of 'gradf.f' fi if test -f 'hlpchk.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'hlpchk.f'\" else echo shar: Extracting \"'hlpchk.f'\" \(177 characters\) sed "s/^X//" >'hlpchk.f' <<'END_OF_FILE' X PROGRAM HLPCHK X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X CHARACTER*6 HELP X NUNIT=10 X OPEN(UNIT=NUNIT,FILE='H.OUT',STATUS='UNKNOWN') X HELP='ALL' X CALL OLHELP(NUNIT,HELP) X STOP X END END_OF_FILE if test 177 -ne `wc -c <'hlpchk.f'`; then echo shar: \"'hlpchk.f'\" unpacked with wrong size! fi # end of 'hlpchk.f' fi if test -f 'i1mach.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'i1mach.f'\" else echo shar: Extracting \"'i1mach.f'\" \(9699 characters\) sed "s/^X//" >'i1mach.f' <<'END_OF_FILE' X INTEGER FUNCTION I1MACH(I) X INTEGER I XC XC I1MACH( 1) = THE STANDARD INPUT UNIT. XC I1MACH( 2) = THE STANDARD OUTPUT UNIT. XC I1MACH( 3) = THE STANDARD PUNCH UNIT. XC I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT. XC I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT. XC I1MACH( 6) = THE NUMBER OF CHARACTERS PER CHARACTER STORAGE UNIT. XC INTEGERS HAVE FORM SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) XC I1MACH( 7) = A, THE BASE. XC I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS. XC I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE. XC FLOATS HAVE FORM SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) XC WHERE EMIN .LE. E .LE. EMAX. XC I1MACH(10) = B, THE BASE. XC SINGLE-PRECISION XC I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS. XC I1MACH(12) = EMIN, THE SMALLEST EXPONENT E. XC I1MACH(13) = EMAX, THE LARGEST EXPONENT E. XC DOUBLE-PRECISION XC I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS. XC I1MACH(15) = EMIN, THE SMALLEST EXPONENT E. XC I1MACH(16) = EMAX, THE LARGEST EXPONENT E. XC X INTEGER IMACH(16), OUTPUT, SC, SMALL(2) X SAVE IMACH, SC X REAL RMACH X EQUIVALENCE (IMACH(4),OUTPUT), (RMACH,SMALL(1)) X INTEGER I3, J, K, T3E(3) X DATA T3E(1) / 9777664 / X DATA T3E(2) / 5323660 / X DATA T3E(3) / 46980 / XC THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES, XC INCLUDING AUTO-DOUBLE COMPILERS. XC TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 XC ON THE NEXT LINE X DATA SC/0/ XC AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. XC CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY XC mail netlib@research.bell-labs.com XC send old1mach from blas XC PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. XC XC MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. XC XC DATA IMACH( 1) / 5 / XC DATA IMACH( 2) / 6 / XC DATA IMACH( 3) / 43 / XC DATA IMACH( 4) / 6 / XC DATA IMACH( 5) / 36 / XC DATA IMACH( 6) / 4 / XC DATA IMACH( 7) / 2 / XC DATA IMACH( 8) / 35 / XC DATA IMACH( 9) / O377777777777 / XC DATA IMACH(10) / 2 / XC DATA IMACH(11) / 27 / XC DATA IMACH(12) / -127 / XC DATA IMACH(13) / 127 / XC DATA IMACH(14) / 63 / XC DATA IMACH(15) / -127 / XC DATA IMACH(16) / 127 /, SC/987/ XC XC MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING XC 32-BIT INTEGER ARITHMETIC. XC XC DATA IMACH( 1) / 5 / XC DATA IMACH( 2) / 6 / XC DATA IMACH( 3) / 7 / XC DATA IMACH( 4) / 6 / XC DATA IMACH( 5) / 32 / XC DATA IMACH( 6) / 4 / XC DATA IMACH( 7) / 2 / XC DATA IMACH( 8) / 31 / XC DATA IMACH( 9) / 2147483647 / XC DATA IMACH(10) / 2 / XC DATA IMACH(11) / 24 / XC DATA IMACH(12) / -127 / XC DATA IMACH(13) / 127 / XC DATA IMACH(14) / 56 / XC DATA IMACH(15) / -127 / XC DATA IMACH(16) / 127 /, SC/987/ XC XC MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. XC XC NOTE THAT THE PUNCH UNIT, I1MACH(3), HAS BEEN SET TO 7 XC WHICH IS APPROPRIATE FOR THE UNIVAC-FOR SYSTEM. XC IF YOU HAVE THE UNIVAC-FTN SYSTEM, SET IT TO 1. XC XC DATA IMACH( 1) / 5 / XC DATA IMACH( 2) / 6 / XC DATA IMACH( 3) / 7 / XC DATA IMACH( 4) / 6 / XC DATA IMACH( 5) / 36 / XC DATA IMACH( 6) / 6 / XC DATA IMACH( 7) / 2 / XC DATA IMACH( 8) / 35 / XC DATA IMACH( 9) / O377777777777 / XC DATA IMACH(10) / 2 / XC DATA IMACH(11) / 27 / XC DATA IMACH(12) / -128 / XC DATA IMACH(13) / 127 / XC DATA IMACH(14) / 60 / XC DATA IMACH(15) /-1024 / XC DATA IMACH(16) / 1023 /, SC/987/ XC X IF (SC .NE. 987) THEN X* *** CHECK FOR AUTODOUBLE *** X SMALL(2) = 0 X RMACH = 1E13 X IF (SMALL(2) .NE. 0) THEN X* *** AUTODOUBLED *** X IF ( (SMALL(1) .EQ. 1117925532 X * .AND. SMALL(2) .EQ. -448790528) X * .OR. (SMALL(2) .EQ. 1117925532 X * .AND. SMALL(1) .EQ. -448790528)) THEN X* *** IEEE *** X IMACH(10) = 2 X IMACH(14) = 53 X IMACH(15) = -1021 X IMACH(16) = 1024 X ELSE IF ( SMALL(1) .EQ. -2065213935 X * .AND. SMALL(2) .EQ. 10752) THEN X* *** VAX WITH D_FLOATING *** X IMACH(10) = 2 X IMACH(14) = 56 X IMACH(15) = -127 X IMACH(16) = 127 X ELSE IF ( SMALL(1) .EQ. 1267827943 X * .AND. SMALL(2) .EQ. 704643072) THEN X* *** IBM MAINFRAME *** X IMACH(10) = 16 X IMACH(14) = 14 X IMACH(15) = -64 X IMACH(16) = 63 X ELSE X WRITE(*,9010) X STOP 777 X END IF X IMACH(11) = IMACH(14) X IMACH(12) = IMACH(15) X IMACH(13) = IMACH(16) X ELSE X RMACH = 1234567. X IF (SMALL(1) .EQ. 1234613304) THEN X* *** IEEE *** X IMACH(10) = 2 X IMACH(11) = 24 X IMACH(12) = -125 X IMACH(13) = 128 X IMACH(14) = 53 X IMACH(15) = -1021 X IMACH(16) = 1024 X SC = 987 X ELSE IF (SMALL(1) .EQ. -1271379306) THEN X* *** VAX *** X IMACH(10) = 2 X IMACH(11) = 24 X IMACH(12) = -127 X IMACH(13) = 127 X IMACH(14) = 56 X IMACH(15) = -127 X IMACH(16) = 127 X SC = 987 X ELSE IF (SMALL(1) .EQ. 1175639687) THEN X* *** IBM MAINFRAME *** X IMACH(10) = 16 X IMACH(11) = 6 X IMACH(12) = -64 X IMACH(13) = 63 X IMACH(14) = 14 X IMACH(15) = -64 X IMACH(16) = 63 X SC = 987 X ELSE IF (SMALL(1) .EQ. 1251390520) THEN X* *** CONVEX C-1 *** X IMACH(10) = 2 X IMACH(11) = 24 X IMACH(12) = -128 X IMACH(13) = 127 X IMACH(14) = 53 X IMACH(15) = -1024 X IMACH(16) = 1023 X ELSE X DO 10 I3 = 1, 3 X J = SMALL(1) / 10000000 X K = SMALL(1) - 10000000*J X IF (K .NE. T3E(I3)) GO TO 20 X SMALL(1) = J X 10 CONTINUE X* *** CRAY T3E *** X IMACH( 1) = 5 X IMACH( 2) = 6 X IMACH( 3) = 0 X IMACH( 4) = 0 X IMACH( 5) = 64 X IMACH( 6) = 8 X IMACH( 7) = 2 X IMACH( 8) = 63 X CALL I1MCR1(IMACH(9), K, 32767, 16777215, 16777215) X IMACH(10) = 2 X IMACH(11) = 53 X IMACH(12) = -1021 X IMACH(13) = 1024 X IMACH(14) = 53 X IMACH(15) = -1021 X IMACH(16) = 1024 X GO TO 35 X 20 CALL I1MCR1(J, K, 16405, 9876536, 0) X IF (SMALL(1) .NE. J) THEN X WRITE(*,9020) X STOP 777 X END IF X* *** CRAY 1, XMP, 2, AND 3 *** X IMACH(1) = 5 X IMACH(2) = 6 X IMACH(3) = 102 X IMACH(4) = 6 X IMACH(5) = 46 X IMACH(6) = 8 X IMACH(7) = 2 X IMACH(8) = 45 X CALL I1MCR1(IMACH(9), K, 0, 4194303, 16777215) X IMACH(10) = 2 X IMACH(11) = 47 X IMACH(12) = -8188 X IMACH(13) = 8189 X IMACH(14) = 94 X IMACH(15) = -8141 X IMACH(16) = 8189 X GO TO 35 X END IF X END IF X IMACH( 1) = 5 X IMACH( 2) = 6 X IMACH( 3) = 7 X IMACH( 4) = 6 X IMACH( 5) = 32 X IMACH( 6) = 4 X IMACH( 7) = 2 X IMACH( 8) = 31 X IMACH( 9) = 2147483647 X 35 SC = 987 X END IF X 9010 FORMAT(/' Adjust autodoubled I1MACH by uncommenting data'/ X * ' statements appropriate for your machine and setting'/ X * ' IMACH(I) = IMACH(I+3) for I = 11, 12, and 13.') X 9020 FORMAT(/' Adjust I1MACH by uncommenting data statements'/ X * ' appropriate for your machine.') X IF (I .LT. 1 .OR. I .GT. 16) GO TO 40 X I1MACH = IMACH(I) X RETURN X 40 WRITE(*,*) 'I1MACH(I): I =',I,' is out of bounds.' X STOP X* /* C source for I1MACH -- remove the * in column 1 */ X* /* Note that some values may need changing. */ X*#include X*#include X*#include X*#include X* X*long i1mach_(long *i) X*{ X* switch(*i){ X* case 1: return 5; /* standard input */ X* case 2: return 6; /* standard output */ X* case 3: return 7; /* standard punch */ X* case 4: return 0; /* standard error */ X* case 5: return 32; /* bits per integer */ X* case 6: return sizeof(int); X* case 7: return 2; /* base for integers */ X* case 8: return 31; /* digits of integer base */ X* case 9: return LONG_MAX; X* case 10: return FLT_RADIX; X* case 11: return FLT_MANT_DIG; X* case 12: return FLT_MIN_EXP; X* case 13: return FLT_MAX_EXP; X* case 14: return DBL_MANT_DIG; X* case 15: return DBL_MIN_EXP; X* case 16: return DBL_MAX_EXP; X* } X* fprintf(stderr, "invalid argument: i1mach(%ld)\n", *i); X* exit(1);return 0; /* some compilers demand return values */ X*} X END X SUBROUTINE I1MCR1(A, A1, B, C, D) X**** SPECIAL COMPUTATION FOR OLD CRAY MACHINES **** X INTEGER A, A1, B, C, D X A1 = 16777216*B + C X A = 16777216*A1 + D X END END_OF_FILE if test 9699 -ne `wc -c <'i1mach.f'`; then echo shar: \"'i1mach.f'\" unpacked with wrong size! fi # end of 'i1mach.f' fi if test -f 'initch.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'initch.f'\" else echo shar: Extracting \"'initch.f'\" \(9587 characters\) sed "s/^X//" >'initch.f' <<'END_OF_FILE' X SUBROUTINE INITCH(INSTOP,LINESR,NEWTON,OVERFL,SCLFCH, X $ SCLXCH,ACPTCR,CONTYP,JACTYP,JUPDM , X $ MAXEXP,N ,NUNIT ,OUTPUT,QNUPDM, X $ STOPCR,TRUPDM,EPSMCH,FCNOLD,FTOL , X $ BOUNDL,BOUNDU,FVECC ,SCALEF,SCALEX, X $ WV1 ,XC ,FVECEV) XC XC AUG. 27, 1991 XC XC THIS SUBROUTINE FIRST CHECKS TO SEE IF N IS WITHIN THE XC ACCEPTABLE RANGE. XC XC THE SECOND CHECK IS TO SEE IF THE INITIAL ESTIMATE IS XC ALREADY A SOLUTION BY THE FUNCTION VALUE CRITERION, FTOL. XC XC THE THIRD CHECK IS MADE TO SEE IF THE NEWTON OPTION IS BEING XC USED WITH THE LINE SEARCH. IF NOT A WARNING IS GIVEN AND XC THE LINE SEARCH OPTION IS INVOKED. XC XC THE FOURTH CHECK IS TO ENSURE APPLICABILITY OF SELECTED XC VALUES FOR INTEGER CONSTANTS. XC XC THE FIFTH CHECK IS TO WARN THE USER IF INITIAL ESTIMATES XC ARE NOT WITHIN THE RANGES SET BY THE BOUNDL AND BOUNDU XC VECTORS. CONTYP IS CHANGED FROM 0 TO 1 IF ANY BOUND HAS XC BEEN SET BY THE USER XC XC THE SIXTH CHECK ENSURES BOUNDL(I) < BOUNDU(I) FOR ALL I. XC X IMPLICIT DOUBLE PRECISION (A-H,O-Z) X INTEGER ACPTCR ,CONTYP ,OUTPUT ,QNUPDM , X $ STOPCR ,TRUPDM X DIMENSION BOUNDL(N),BOUNDU(N),FVECC(N) ,SCALEF(N), X $ SCALEX(N),WV1(N) ,XC(N) X LOGICAL FRSTER ,INSTOP ,LINESR ,NEWTON , X $ OVERFL ,SCLFCH ,SCLXCH X EXTERNAL FVECEV X DATA ZERO,ONE,TEN /0.0D0,1.0D0,10.0D0/ XC X INSTOP=.FALSE. X TEMP1=-TEN**MAXEXP X TEMP2=TEN**MAXEXP XC XC CHECK FOR N IN RANGE. XC X IF(N.LE.0) THEN X INSTOP=.TRUE. X WRITE(NUNIT,1) X1 FORMAT(T3,72('*')) X WRITE(NUNIT,2) X2 FORMAT(T3,'*',T74,'*') X WRITE(NUNIT,3) X3 FORMAT(T3,'*',2X,'N IS OUT OF RANGE - RESET TO POSITIVE', X $ ' INTEGER',T74,'*') X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF XC XC CHECK FOR SCALING FACTORS POSITIVE. XC X FRSTER=.TRUE. X SCLFCH=.FALSE. X SCLXCH=.FALSE. X DO 100 I=1,N X IF(SCALEF(I).LE.ZERO) THEN X IF(FRSTER) THEN X INSTOP=.TRUE. X FRSTER=.FALSE. X WRITE(NUNIT,1) X END IF X WRITE(NUNIT,2) X WRITE(NUNIT,4) I,SCALEF(I) X4 FORMAT(T3,'*',7X,'SCALEF(',I3,') = ',1PD12.3, X $ 4X,'SHOULD BE POSITIVE',T74,'*') X END IF X IF(SCALEF(I).NE.ONE) SCLFCH=.TRUE. X IF(SCALEX(I).LE.ZERO) THEN X IF(FRSTER) THEN X INSTOP=.TRUE. X FRSTER=.FALSE. X WRITE(NUNIT,1) X END IF X WRITE(NUNIT,2) X WRITE(NUNIT,5) I,SCALEX(I) X5 FORMAT(T3,'*',7X,'SCALEX(',I3,') = ',1PD12.3, X $ 4X,'SHOULD BE POSITIVE',T74,'*') X END IF X IF(SCALEX(I).NE.ONE) SCLXCH=.TRUE. X100 CONTINUE X IF(.NOT.FRSTER) THEN X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF XC XC EVALUATE INITIAL RESIDUAL VECTOR AND OBJECTIVE FUNCTION AND XC CHECK TO SEE IF THE INITIAL GUESS IS ALREADY A SOLUTION. XC X CALL FVECEV(OVERFL,N,FVECC,XC) XC XC NOTE: NUMBER OF LINE SEARCH FUNCTION EVALUATIONS, NFUNC, XC INITIALIZED AT 1 WHICH REPRESENTS THIS EVALUATION. XC X IF(OVERFL) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,6) X6 FORMAT(T3,'*',7X,'OVERFLOW IN INITIAL FUNCTION', X $ ' VECTOR EVALUATION',T74,'*') X WRITE(NUNIT,2) X WRITE(NUNIT,1) X INSTOP=.TRUE. X RETURN X END IF X CALL FCNEVL(OVERFL,MAXEXP,N ,NUNIT ,OUTPUT,EPSMCH, X $ FCNOLD,FVECC ,SCALEF,WV1 ) X IF(OVERFL) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,7) X7 FORMAT(T3,'*',7X,'OVERFLOW IN INITIAL OBJECTIVE ', X $ ' FUNCTION EVALUATION',T74,'*') X WRITE(NUNIT,2) X WRITE(NUNIT,1) X INSTOP=.TRUE. X RETURN X END IF XC XC CHECK FOR SOLUTION USING SECOND STOPPING CRITERION. XC X DO 200 I=1,N X IF(ABS(FVECC(I)).GT.FTOL) GO TO 201 X200 CONTINUE X INSTOP=.TRUE. X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,8) X8 FORMAT(T3,'*',2X,'WARNING: THIS IS ALREADY A SOLUTION', X $ ' BY THE CRITERIA OF THE SOLVER',T74,'*') X WRITE(NUNIT,2) XC XC IF THE PROBLEM IS BADLY SCALED THE OBJECTIVE FUNCTION XC MAY MEET THE TOLERANCE ALTHOUGH THE INITIAL ESTIMATE XC IS NOT THE SOLUTION. XC X WRITE(NUNIT,9) X9 FORMAT(T3,'*',2X,'THIS MAY POSSIBLY BE ALLEVIATED BY ', X $ 'RESCALING THE PROBLEM IF THE',T74,'*') X WRITE(NUNIT,10) X10 FORMAT(T3,'*',2X,'INITIAL ESTIMATE IS KNOWN NOT TO BE', X $ ' A SOLUTION',T74,'*') X WRITE(NUNIT,2) X WRITE(NUNIT,1) X201 CONTINUE ! FVEC(I) > FTOL FOR SOME I FROM 200 LOOP XC XC CHECK FOR NEWTON'S METHOD REQUESTED BUT LINE SEARCH NOT XC BEING USED. XC X IF(NEWTON.AND.(.NOT.LINESR)) THEN X LINESR=.TRUE. X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,11) X11 FORMAT(T3,'*',2X,'WARNING: INCOMPATIBLE OPTIONS', X $ ': NEWTON=.TRUE. AND LINESR=.FALSE.',T74,'*') X WRITE(NUNIT,12) X12 FORMAT(T3,'*',2X,'LINESR SET TO .TRUE.; EXECUTION' X $ ' OF NEWTON METHOD CONTINUING',T74,'*') X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF XC XC CHECK INTEGER CONSTANTS. XC X IF(ACPTCR.NE.1.AND.ACPTCR.NE.12) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,13) ACPTCR X13 FORMAT(T3,'*',2X,'ACPTCR NOT AN ACCEPTABLE VALUE: ', X $ I5,T74,'*') X INSTOP=.TRUE. X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF X IF(JACTYP.LT.0.OR.JACTYP.GT.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,14) JACTYP X14 FORMAT(T3,'*',2X,'JACTYP:',I5,' - NOT IN PROPER RANGE', X $ T74,'*') X INSTOP=.TRUE. X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF X IF(STOPCR.NE.1.AND.STOPCR.NE.12.AND.STOPCR.NE.2. X $ AND.STOPCR.NE.3) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,15) STOPCR X15 FORMAT(T3,'*',2X,'STOPCR NOT AN ACCEPTABLE VALUE: ', X $ I5,T74,'*') X INSTOP=.TRUE. X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF X IF(QNUPDM.LT.0.OR.QNUPDM.GT.1) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,16) QNUPDM X16 FORMAT(T3,'*',2X,'QNUPDM:',I5,' - NOT IN PROPER RANGE', X $ T74,'*') X INSTOP=.TRUE. X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF X IF(TRUPDM.LT.0.OR.TRUPDM.GT.1) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,17) TRUPDM X17 FORMAT(T3,'*',2X,'TRUPDM:',I5,' - NOT IN PROPER RANGE', X $ T74,'*') X INSTOP=.TRUE. X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF X IF(JUPDM.LT.0.OR.JUPDM.GT.2) THEN X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,18) JUPDM X18 FORMAT(T3,'*',2X,'JUPDM:',I5,' - NOT IN PROPER RANGE', X $ T74,'*') X INSTOP=.TRUE. X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF XC XC CHECK FOR INITIAL ESTIMATES NOT WITHIN SPECIFIED BOUNDS AND XC SET CONTYP TO 1 => AT LEAST ONE BOUND IS IN EFFECT. XC X CONTYP=0 X DO 300 I=1,N X IF((BOUNDL(I).NE.TEMP1.OR.BOUNDU(I).NE.TEMP2)) THEN X CONTYP=1 X GO TO 301 X END IF X300 CONTINUE X301 CONTINUE X FRSTER=.TRUE. X IF(CONTYP.NE.0) THEN X DO 400 I=1,N XC XC CHECK FOR INITIAL ESTIMATES OUT OF RANGE AND LOWER XC BOUND GREATER THAN OR EQUAL TO THE UPPER BOUND. XC X IF(XC(I).LT.BOUNDL(I).OR.XC(I).GT.BOUNDU(I)) THEN X IF(FRSTER) THEN X INSTOP=.TRUE. X FRSTER=.FALSE. X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,19) X19 FORMAT(T3,'*',7X,'COMPONENTS MUST BE WITHIN', X $ ' BOUNDS',T74,'*') X WRITE(NUNIT,2) X WRITE(NUNIT,20) X20 FORMAT(T3,'*',8X,'NO.',9X,'XC',16X,'BOUNDL',10X, X $ 'BOUNDU',T74,'*') X WRITE(NUNIT,2) X END IF X WRITE(NUNIT,21) I,XC(I),BOUNDL(I),BOUNDU(I) X21 FORMAT(T3,'*',7X,I3,3X,1PD12.3,9X,1PD12.3,4X, X $ 1PD12.3,T74,'*') X END IF X400 CONTINUE X IF(.NOT.FRSTER) THEN X WRITE(NUNIT,2) X WRITE(NUNIT,1) X END IF XC X FRSTER=.TRUE. X DO 500 I=1,N X IF(BOUNDL(I).GE.BOUNDU(I)) THEN X IF(FRSTER) THEN X FRSTER=.FALSE. X WRITE(NUNIT,1) X WRITE(NUNIT,2) X WRITE(NUNIT,22) X22 FORMAT(T3,'*',7X,'LOWER BOUND MUST BE LESS THAN', X $