C ALGORITHM 717, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 19, NO. 1, MARCH, 1993, PP. 109-130. cat >README <<'//GO.SYSIN DD README' Subroutines for Maximum Likelihood and Quasi-Likelihood Estimation of Parameters in Nonlinear Regression Models by David S. Bunch (UC Davis), David M. Gay (AT&T Bell Laboratories), and Roy E. Welsch (MIT); submission to ACM Transactions on Mathematical Software. _____________ Preliminaries ============= To use the Fortran subroutines and main programs in this large file, you need to split it into its 53 constituent files. The first of these logical files is called README, which includes the text you are reading now, as well as additional information for linking and running the test programs. The instructions for decomposing the files is given next, followed by the remainder of README, and then the other files. ____________________ Splitting into Files ==================== If you are using a UNIX system, just feed this file to /bin/sh (in an empty directory), as in sh thisfilename Alternatively, you can feed this file to the Fortran program shown below. You could also split this large file up by hand: each constituent file is preceded by a line of the form cat >filename <<'//GO.SYSIN DD filename' and followed by a corresponding line of the form //GO.SYSIN DD filename We've indented the above lines here for display purposes, but the real lines start in column 1. Here is the promised Fortran program for decomposing this file. On some systems you will need to make a change in format statement 90, as described in the comment below. PROGRAM UNPACK CHARACTER*100 LINE, TERMIN CHARACTER*16 FNAME INTEGER FNLEN, I, LINENO, OUTFIL PARAMETER (OUTFIL=1) LINENO = 0 10 READ(*,'(A)',END=999) LINE LINENO = LINENO + 1 IF (LINE(1:5) .NE. 'cat >') CALL BADCAT(LINE, LINENO) DO 20 I = 6, 100 IF (LINE(I:I) .EQ. ' ') GO TO 30 20 CONTINUE CALL BADCAT(LINE, LINENO) 30 FNAME = LINE(6:I) FNLEN = I - 5 IF (LINE(I+1:I+17) .NE. '<<''//GO.SYSIN DD ') 1 CALL BADCAT(LINE, LINENO) IF (LINE(I+18:I+16+FNLEN) .NE. FNAME) 1 CALL BADCAT(LINE, LINENO) TERMIN = LINE(I+4:I+FNLEN+16) OPEN(OUTFIL,FILE=FNAME,STATUS='NEW',ERR=40) GO TO 60 40 OPEN(OUTFIL,FILE=FNAME,STATUS='OLD',ERR=50) WRITE(*,*) 'overwriting ',FNAME(1:FNLEN) GO TO 70 50 WRITE(*,*) 'cannot open ',FNAME(1:FNLEN) GO TO 999 60 WRITE(*,*) FNAME(1:FNLEN) 70 LINENO = LINENO + 1 READ(*,'(A)',END=100) LINE IF (LINE .EQ. TERMIN) THEN CLOSE(OUTFIL) GO TO 10 END IF DO 80 I = 100, 1, -1 IF (LINE(I:I) .NE. ' ') GO TO 90 80 CONTINUE ******** On systems where carriage controls end up in written files ******** (to be honored by a printer or subsequent program), such as ******** VAX VMS and most UNIX systems, you need to omit "1X," ******** in the following WRITE statement, changing it to *90 WRITE(OUTFIL,'(A)') LINE(I:I) 90 WRITE(OUTFIL,'(1X,A)') LINE(1:I) GO TO 70 100 WRITE(*,*) 'Premature end of file' 999 END SUBROUTINE BADCAT(LINE, LINENO) CHARACTER*100 LINE INTEGER LINENO WRITE(*,*) 'Line ', LINENO, ': Bad "cat" line:' WRITE(*, '(A)') LINE STOP END ________ Overview ======== Information is given below regarding subroutines and test programs, including how to link and run the test programs. There are five basic test programs: MADSEN (simple test problems, no bounds on parameters) MADSENB (simple test problems, with bounds on parameters) PMAIN (problems from Gay and Welsch 1988, mix of bounds and no bounds) MLMNP (multinomial probit estimation problems from Bunch 1991, no bounds) MLMNPB (" ", with bounds) These exist in both single- and double-precision versions. The following two documents are available from the authors: "Driver PMAIN for DGL[FG][B ]", by David M. Gay. "MLMNP and MLMNPB: Fortran Programs for Maximum Likelihood Estimation of Linear-in-Parameter Multinomial Probit Models", by David S. Bunch. Postscript for these documents is or will be available from netlib. For details, send netlib@research.att.com the one-line electronic mail message send index from opt/nlr MADSEN and MADSENB do not require input. PMAIN has a single sample input file for producing test results, but can also be run interactively or with other input files; see the "Driver PMAIN" document. MLMNP and MLMNPB require input files for Fortran units 1 and 2; four examples are included. ____________________ Machine dependencies ==================== There are two machine-dependent files, dmdc.f and smdc.f (double- and single-precision versions, respectively), which provide machine constants. You must activate (i.e., remove the 'C' from column 1) the lines that pertain to your machine, or that pertain to the PORT routines D1MACH (for dmdc.f) and R1MACH (for smdc.f), if you choose to use those routines (which have constants for a much wider selection of machines than do dmdc.f and smdc.f). __________ Precisions ========== As previously noted, we've provided both single- and double-precision versions of our optimization subroutines and test problems. If you are a referee, you may want to try both; otherwise, unless you are using a Cray or CDC machine (or something similar whose single precision has substantially more accuracy than does binary IEEE aritihmetic), you are probably better off using the double-precision routines. ________ makefile ======== If you are using a UNIX system, you can probably just type make to cause everything to be compiled and all test programs to be run. If you run the single-precision tests on a 32-bit computer, you may get a few instances of FALSE or SINGULAR CONVERGENCE. For comparison purposes, we include *.sgi files, which are *.out files we obtained on an SGI computer (IEEE arithmetic; compilation was with f2c and cc). We note, however, that many of the test problems are very nonlinear, and differences in compilers and machines will often produce slightly different output. ____________________________________ Opening files from the test programs ==================================== To run the MLMNP and MLMNPB programs, you may need to adjust the OPEN and REWIND statements near the beginning of mlmnp.f and mlmnpb.f (for double-precision, or smlmnp.f and smlmnpb.f for single). ________________ Summary of files ================ 1. README This file. 2. makefile For UNIX systems; encodes the information below about what files are needed for what programs. DOUBLE PRECISION SOURCE FILES 3. dmdc.f0 Edit this into dmdc.f . 4. dglfg.f Top-level routines DGLG, DGLF, DRGLG (no bounds), followed by routines they call that are not in dgletc.f . 5. dglfgb.f Top-level routines DGLGB, DGLFB, DRGLGB (simple bounds), followed by routines they call that are not in dgletc.f . 6. dgletc.f Routines needed whether or not there are simple bounds. 7. madsen.f Simple example program, no bounds. Needs dmdc.f, dglfg.f, dgletc.f . 8. madsenb.f Simple test program, variant of madsen.f with bounds. Needs dmdc.f, dglfgb.f, dgletc.f . 9. dpmain.f General test program PMAIN. Needs dmdc.f, dglfg.f, dglfgb.f, dgletc.f, and mecdf.f . 10. mecdf.f Computes approximation to multivariate normal cumulative distribution function (uses Mendell- Elston approximation.) 11. mlmnp.f Program MLMNP for linear-in-parameter multinomial probit models, no bounds. Needs dmdc.f, dglfg.f, dgletc.f, mecdf.f, mnpsubs.f . 12. mlmnpb.f Program MLMNPB for linear-in-parameter multinomial probit models with simple bounds. Needs dmdc.f, dglfgb.f, dgletc.f, mecdf.f, mnpsubs.f . 13. mnpsubs.f Needed by mlmnp.f and mlmnpb.f . TEST DATA FILES 14. pmain.in Input for PMAIN (from Gay & Welsch, 1988). The following *.fu? files are input for MLMNP and MLMNPB. The files named *.fu1 are to be read by Fortran unit 1, and those named *.fu2 are to be read by Fortran unit 2. 15. daganzo.fu2 Choice data set from Daganzo (1979). 16. mnpex1.fu1 Example 1: a model to use with daganzo.fu2 17. mnpex2.fu1 Example 2: another model to use with daganzo.fu2 18. rent.fu2 Choice data set from MBA survey on rental housing 19. rent1.fu1 Example 3: a model to use with rent.fu2 20. rent2.fu1 Example 4: another model to use with rent.fu2 SINGLE PRECISION SOURCE FILES corresponding to files 3-13. 21. smdc.f0 Edit this into smdc.f . 22. sglfg.f 23. sglfgb.f 24. sgletc.f 25. smadsen.f 26. smadsenb.f 27. spmain.f 28. smecdf.f 29. smlmnp.f 30. smlmnpb.f 31. smnpsubs.f SAMPLE OUTPUTS, DOUBLE PRECISION 32. madsen.sgi 33. madsenb.sgi 34. mnpex1.sgi 35. mnpex1b.sgi 36. mnpex2.sgi 37. mnpex2b.sgi 38. pmain.sgi 39. rent1.sgi 40. rent1b.sgi 41. rent2.sgi 42. rent2b.sgi SAMPLE OUTPUTS, SINGLE PRECISION 43. smadsen.sgi 44. smadsenb.sgi 45. smnpex1.sgi 46. smnpex1b.sgi 47. smnpex2.sgi 48. smnpex2b.sgi 49. spmain.sgi 50. srent1.sgi 51. srent1b.sgi 52. srent2.sgi 53. srent2b.sgi //GO.SYSIN DD README cat >makefile <<'//GO.SYSIN DD makefile' .SUFFIXES: .f .o FFLAGS = -u F77 = f77 L = .f.o: $(F77) -c $(FFLAGS) $*.f both: out sout out: madsen.out madsenb.out pmain.out mnpex1.out mnpex1b.out \ mnpex2.out mnpex2b.out rent1.out rent1b.out rent2.out rent2b.out sout: smadsen.out smadsenb.out spmain.out smnpex1.out smnpex1b.out \ smnpex2.out mnpex2b.out srent1.out srent1b.out srent2.out srent2b.out dmdc.f: dmdc.f0 true # Obtain dmcd.f from dmdc.f0 by activating the statements false # appropriate to your machine smdc.f: smdc.f0 true # Obtain dmcd.f from smdc.f0 by activating the statements false # appropriate to your machine madsen.out: madsen.f dglfg.o dgletc.o dmdc.o $(F77) madsen.f dglfg.o dgletc.o dmdc.o $L a.out >$@ madsenb.out: madsenb.f dglfgb.o dgletc.o dmdc.o $(F77) madsenb.f dglfgb.o dgletc.o dmdc.o $L a.out >$@ pmain: dpmain.o mecdf.o dglfg.o dglfgb.o dgletc.o dmdc.o $(F77) dpmain.o mecdf.o dglfg.o dglfgb.o dgletc.o dmdc.o $L -o $@ pmain.out: pmain pmain.in pmain $@ mlmnp: mlmnp.o mecdf.o mnpsubs.o dglfg.o dgletc.o dmdc.o $(F77) mlmnp.o mecdf.o mnpsubs.o dglfg.o dgletc.o dmdc.o $L -o $@ mlmnpb: mlmnpb.o mecdf.o mnpsubs.o dglfgb.o dgletc.o dmdc.o $(F77) mlmnpb.o mecdf.o mnpsubs.o dglfgb.o dgletc.o dmdc.o $L -o $@ mnpex1.out mnpex1b.out: mlmnp mlmnpb mnpex1.fu1 daganzo.fu2 rm -f fort.? ln mnpex1.fu1 fort.1 ln daganzo.fu2 fort.2 mlmnp >mnpex1.out mlmnpb >mnpex1b.out mnpex2.out mnpex2b.out: mlmnp mlmnpb mnpex2.fu1 daganzo.fu2 rm -f fort.? ln mnpex2.fu1 fort.1 ln daganzo.fu2 fort.2 mlmnp >mnpex2.out mlmnpb >mnpex2b.out rent1.out rent1b.out: mlmnp mlmnpb rent1.fu1 rent.fu2 rm -f fort.? ln rent1.fu1 fort.1 ln rent.fu2 fort.2 mlmnp >rent1.out mlmnpb >rent1b.out rent2.out rent2b.out: mlmnp mlmnpb rent2.fu1 rent.fu2 rm -f fort.? ln rent2.fu1 fort.1 ln rent.fu2 fort.2 mlmnp >rent2.out mlmnpb >rent2b.out # single-precision runs... smadsen.out: smadsen.f sglfg.o sgletc.o smdc.o $(F77) smadsen.f sglfg.o sgletc.o smdc.o $L a.out >$@ smadsenb.out: smadsenb.f sglfgb.o sgletc.o smdc.o $(F77) smadsenb.f sglfgb.o sgletc.o smdc.o $L a.out >$@ spmain: spmain.o smecdf.o sglfg.o sglfgb.o sgletc.o smdc.o $(F77) spmain.o smecdf.o sglfg.o sglfgb.o sgletc.o smdc.o $L -o $@ spmain.out: spmain pmain.in spmain $@ smlmnp: smlmnp.o smecdf.o smnpsubs.o sglfg.o sgletc.o smdc.o $(F77) smlmnp.o smecdf.o smnpsubs.o sglfg.o sgletc.o smdc.o $L -o $@ smlmnpb: smlmnpb.o smecdf.o smnpsubs.o sglfgb.o sgletc.o smdc.o $(F77) smlmnpb.o smecdf.o smnpsubs.o sglfgb.o sgletc.o smdc.o $L -o $@ smnpex1.out smnpex1b.out: smlmnp smlmnpb mnpex1.fu1 daganzo.fu2 rm -f fort.? ln mnpex1.fu1 fort.1 ln daganzo.fu2 fort.2 smlmnp >smnpex1.out smlmnpb >smnpex1b.out smnpex2.out smnpex2b.out: smlmnp smlmnpb mnpex2.fu1 daganzo.fu2 rm -f fort.? ln mnpex2.fu1 fort.1 ln daganzo.fu2 fort.2 smlmnp >smnpex2.out smlmnpb >smnpex2b.out srent1.out srent1b.out: smlmnp smlmnpb rent1.fu1 rent.fu2 rm -f fort.? ln rent1.fu1 fort.1 ln rent.fu2 fort.2 smlmnp >srent1.out smlmnpb >srent1b.out srent2.out srent2b.out: smlmnp smlmnpb rent2.fu1 rent.fu2 rm -f fort.? ln rent2.fu1 fort.1 ln rent.fu2 fort.2 smlmnp >srent2.out smlmnpb >srent2b.out clean: rm -f *.out *.o pmain mlmnp mlmnpb spmain smlmnp smlmnpb //GO.SYSIN DD makefile cat >dmdc.f0 <<'//GO.SYSIN DD dmdc.f0' DOUBLE PRECISION FUNCTION DR7MDC(K) C C *** RETURN MACHINE DEPENDENT CONSTANTS USED BY NL2SOL *** C C +++ COMMENTS BELOW CONTAIN DATA STATEMENTS FOR VARIOUS MACHINES. +++ C +++ TO CONVERT TO ANOTHER MACHINE, PLACE A C IN COLUMN 1 OF THE +++ C +++ DATA STATEMENT LINE(S) THAT CORRESPOND TO THE CURRENT MACHINE +++ C +++ AND REMOVE THE C FROM COLUMN 1 OF THE DATA STATEMENT LINE(S) +++ C +++ THAT CORRESPOND TO THE NEW MACHINE. +++ C INTEGER K C C *** THE CONSTANT RETURNED DEPENDS ON K... C C *** K = 1... SMALLEST POS. ETA SUCH THAT -ETA EXISTS. C *** K = 2... SQUARE ROOT OF ETA. C *** K = 3... UNIT ROUNDOFF = SMALLEST POS. NO. MACHEP SUCH C *** THAT 1 + MACHEP .GT. 1 .AND. 1 - MACHEP .LT. 1. C *** K = 4... SQUARE ROOT OF MACHEP. C *** K = 5... SQUARE ROOT OF BIG (SEE K = 6). C *** K = 6... LARGEST MACHINE NO. BIG SUCH THAT -BIG EXISTS. C DOUBLE PRECISION BIG, ETA, MACHEP, ZERO INTEGER BIGI(2), ETAI(2), MACHEI(2) EQUIVALENCE (BIG,BIGI(1)), (ETA,ETAI(1)), (MACHEP,MACHEI(1)) PARAMETER (ZERO=0.D+0) C C +++ IEEE ARITHMETIC MACHINES IN WHICH THE MOST SIGNIFICANT BYTE C +++ IS STORED FIRST, SUCH AS THE AT&T 3B SERIES AND MACHINES C +++ BASED ON SPARC, MIPS, AND MOTOROLA 68XXX PROCESSORS. C C DATA BIGI(1),BIGI(2) / 2146435071, -1 / C DATA ETAI(1),ETAI(2) / 1048576, 0 / C DATA MACHEI(1),MACHEI(2) / 1017118720, 0 / C C +++ IEEE ARITHMETIC MACHINES IN WHICH THE LEAST SIGNIFICANT BYTE C +++ IS STORED FIRST, SUCH AS MACHINES BASED ON INTEL PROCESSORS, C +++ E.G. PERSONAL COMPUTERS WITH AN INTEL 80X87. C C DATA BIGI(1),BIGI(2) / -1, 2146435071 / C DATA ETAI(1),ETAI(2) / 0, 1048576 / C DATA MACHEI(1),MACHEI(2) / 0, 1017118720 / C C +++ IBM, AMDAHL, OR XEROX MAINFRAME +++ C C DATA BIGI(1),BIGI(2)/2147483647, -1/ C DATA ETAI(1),ETAI(2)/1048576, 0/ C DATA MACHEI(1),MACHEI(2)/873463808,0/ C C +++ VAX +++ C C DATA BIGI(1),BIGI(2) / -32769, -1 / C DATA ETAI(1),ETAI(2) / 128, 0 / C DATA MACHEI(1),MACHEI(2) / 9344, 0 / C C +++ CRAY +++ C C DATA BIGI(1)/6917247552664371199/ C DATA BIGI(2)/128891879815246481/ C DATA ETAI(1)/2332160919536140288/ C DATA ETAI(2)/0/ C DATA MACHEI(1)/4585931058058362880/ C DATA MACHEI(2)/0/ C C +++ PORT LIBRARY -- REQUIRES MORE THAN JUST A DATA STATEMENT, +++ C +++ BUT HAS CONSTANTS FOR MANY MORE MACHINES. +++ C C To get the current D1MACH, which has constants for many more C machines, ask netlib@research.att.com to C send d1mach from cor C For machines with rounded arithmetic (e.g., IEEE or VAX arithmetic), C use MACHEP = 0.5D0 * D1MACH(4) below. C C DOUBLE PRECISION D1MACH C EXTERNAL D1MACH C DATA BIG/0.D+0/, ETA/0.D+0/, MACHEP/0.D+0/, ZERO/0.D+0/ C IF (BIG .GT. ZERO) GO TO 1 C BIG = D1MACH(2) C ETA = D1MACH(1) C MACHEP = D1MACH(4) C1 CONTINUE C C +++ END OF PORT +++ C C------------------------------- BODY -------------------------------- C IF (MACHEP .LE. ZERO) THEN WRITE(*,*) 'Edit DR7MDC to activate the appropriate statements' STOP 987 ENDIF GO TO (10, 20, 30, 40, 50, 60), K C 10 DR7MDC = ETA GO TO 999 C 20 DR7MDC = SQRT(256.D+0*ETA)/16.D+0 GO TO 999 C 30 DR7MDC = MACHEP GO TO 999 C 40 DR7MDC = SQRT(MACHEP) GO TO 999 C 50 DR7MDC = SQRT(BIG/256.D+0)*16.D+0 GO TO 999 C 60 DR7MDC = BIG C 999 RETURN C *** LAST LINE OF DR7MDC FOLLOWS *** END INTEGER FUNCTION I7MDCN(K) C INTEGER K C C *** RETURN INTEGER MACHINE-DEPENDENT CONSTANTS *** C C *** K = 1 MEANS RETURN STANDARD OUTPUT UNIT NUMBER. *** C *** K = 2 MEANS RETURN ALTERNATE OUTPUT UNIT NUMBER. *** C *** K = 3 MEANS RETURN INPUT UNIT NUMBER. *** C (NOTE -- K = 2, 3 ARE USED ONLY BY TEST PROGRAMS.) C C +++ PORT VERSION FOLLOWS... C INTEGER I1MACH C EXTERNAL I1MACH C INTEGER MDPERM(3) C DATA MDPERM(1)/2/, MDPERM(2)/4/, MDPERM(3)/1/ C I7MDCN = I1MACH(MDPERM(K)) C +++ END OF PORT VERSION +++ C C +++ NON-PORT VERSION FOLLOWS... INTEGER MDCON(3) DATA MDCON(1)/6/, MDCON(2)/8/, MDCON(3)/5/ I7MDCN = MDCON(K) C +++ END OF NON-PORT VERSION +++ C 999 RETURN C *** LAST LINE OF I7MDCN FOLLOWS *** END //GO.SYSIN DD dmdc.f0 cat >dglfg.f <<'//GO.SYSIN DD dglfg.f' SUBROUTINE DGLG(N, P, PS, X, RHO, RHOI, RHOR, IV, LIV, LV, V, 1 CALCRJ, UI, UR, UF) C C *** GENERALIZED LINEAR REGRESSION A LA NL2SOL *** C C *** PARAMETERS *** C INTEGER N, P, PS, LIV, LV INTEGER IV(LIV), RHOI(*), UI(*) DOUBLE PRECISION X(*), RHOR(*), V(LV), UR(*) EXTERNAL CALCRJ, RHO, UF C C *** PARAMETER USAGE *** C C N....... TOTAL NUMBER OF RESIDUALS. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S). C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS. C SEE DRGLG FOR DETAILS ABOUT RHO. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV, AT LEAST 90 + P. C LV...... LENGTH OF V, AT LEAST C 105 + P*(3*P + 16) + 2*N + 4*PS C + N*(P + 1 + (P-PS+1)*(P-PS+2)/2). C V....... FLOATING-POINT VALUES ARRAY. C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR AND JACOBIAN MATRIX. C UI...... PASSED UNCHANGED TO CALCRJ. C UR...... PASSED UNCHANGED TO CALCRJ. C UF...... PASSED UNCHANGED TO CALCRJ. C C *** CALCRJ CALLING SEQUENCE... C C CALL CALCRJ(N, P, X, NF, NEED, R, RP, UI, UR, UF) C C PARAMETERS N, P, X, UI, UR, AND UF ARE AS ABOVE. C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(P,N). C NEED IS AN INTEGER ARRAY OF LENGTH 2... C NEED(1) = 1 MEANS CALCRJ SHOULD COMPUTE THE RESIDUAL VECTOR R, C AND NEED(2) IS THE VALUE NF HAD AT THE LAST X WHERE C CALCRJ MIGHT BE CALLED WITH NEED(1) = 2. C NEED(1) = 2 MEANS CALCRJ SHOULD COMPUTE THE JACOBIAN MATRIX RP, C WHERE RP(J,I) = DERIVATIVE OF R(I) WITH RESPECT TO X(J). C (CALCRJ SHOULD NOT CHANGE NEED AND SHOULD CHANGE AT MOST ONE OF R C AND RP. IF R OR RP, AS APPROPRIATE, CANNOT BE COMPUTED, THEN CALCRJ C SHOULD SET NF TO 0. OTHERWISE IT SHOULD NOT CHANGE NF.) C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DIVSET, DRGLG C C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DRGLG ... CARRIES OUT OPTIMIZATION ITERATIONS. C C C *** LOCAL VARIABLES *** C INTEGER D1, DR1, I, IV1, NEED1(2), NEED2(2), NF, R1, RD1 C C *** IV COMPONENTS *** C INTEGER D, J, NEXTV, NFCALL, NFGCAL, R, REGD, REGD0, TOOBIG, VNEED PARAMETER (D=27, J=70, NEXTV=47, NFCALL=6, NFGCAL=7, R=61, 1 REGD=67, REGD0=82, TOOBIG=2, VNEED=4) SAVE NEED1, NEED2 DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 I = (P-PS+2)*(P-PS+1)/2 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+1+I) CALL DRGLG(X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, RHOI, 1 RHOR, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + (P - PS + 1)*N IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N IV(NEXTV) = IV(J) + N*PS IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RD1 = IV(REGD0) C 20 CALL DRGLG(V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, V(R1), 1 V(RD1), RHO, RHOI, RHOR, V, X) IF (IV(1)-2) 30, 50, 60 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) NEED1(2) = IV(NFGCAL) CALL CALCRJ(N, PS, X, NF, NEED1, V(R1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE DR = GRADIENT OF R COMPONENTS *** C 50 CALL CALCRJ(N, PS, X, IV(NFGCAL), NEED2, V(R1), V(DR1), UI, UR,UF) IF (IV(NFGCAL) .EQ. 0) IV(TOOBIG) = 1 GO TO 20 C C *** INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED C *** AND PRINT IT IF SO REQUESTED... C 60 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 C 999 RETURN C C *** LAST LINE OF DGLG FOLLOWS *** END SUBROUTINE DGLF(N, P, PS, X, RHO, RHOI, RHOR, IV, LIV, LV, V, 1 CALCRJ, UI, UR, UF) C C *** GENERALIZED LINEAR REGRESSION, FINITE-DIFFERENCE JACOBIAN *** C C *** PARAMETERS *** C INTEGER N, P, PS, LIV, LV INTEGER IV(LIV), RHOI(*), UI(*) DOUBLE PRECISION X(*), V(LV), RHOR(*), UR(*) EXTERNAL CALCRJ, RHO, UF C C *** PARAMETER USAGE *** C C N....... TOTAL NUMBER OF RESIDUALS. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S). C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS. C SEE DRGLG FOR DETAILS ABOUT RHO. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV, AT LEAST 90 + P. C LV...... LENGTH OF V, AT LEAST C 105 + P*(3*P + 16) + 2*N + 4*PS C + N*(P + 3 + (P-PS+1)*(P-PS+2)/2). C V....... FLOATING-POINT VALUES ARRAY. C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR. C UI...... PASSED UNCHANGED TO CALCRJ. C UR...... PASSED UNCHANGED TO CALCRJ. C UF...... PASSED UNCHANGED TO CALCRJ. C C *** CALCRJ CALLING SEQUENCE... C C CALL CALCRJ(N, P, X, NF, NEED, R, RP, UI, UR, UF) C C PARAMETERS N, P, X, UI, UR, AND UF ARE AS ABOVE. C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(P,N). C NEED MAY BE REGARDED AS AN INTEGER THAT ALWAYS HAS THE VALUE 1 C WHEN DGLF CALLS CALCRJ. THIS MEANS CALCRJ SHOULD COMPUTE THE C RESIDUAL VECTOR R. (CALCRJ SHOULD NOT CHANGE NEED OR RP. IF R C CANNOT BE COMPUTED, THEN CALCRJ SHOULD SET NF TO 0. OTHERWISE IT C SHOULD NOT CHANGE NF. FOR COMPATIBILITY WITH DGLG, NEED IS A C VECTOR OF LENGTH 2.) C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DIVSET, DRGLG,DV7CPY C C DIVSET... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DRGLG... CARRIES OUT OPTIMIZATION ITERATIONS. C DV7CPY... COPIES ONE VECTOR TO ANOTHER. C C *** LOCAL VARIABLES *** C INTEGER D1, DK, DR1, I, I1, IV1, J1K, J1K0, K, NEED(2), NF, 1 NG, RD1, R1, R21, RN, RS1 DOUBLE PRECISION H, H0, HLIM, NEGPT5, ONE, XK, ZERO C C *** IV AND V COMPONENTS *** C INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL, 1 NGCALL, NGCOV, R, RDREQ, REGD, REGD0, TOOBIG, VNEED PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35, 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53, 2 R=61, RDREQ=57, REGD=67, REGD0=82, TOOBIG=2, VNEED=4) SAVE NEED DATA HLIM/0.1D+0/, NEGPT5/-0.5D+0/, ONE/1.D+0/, ZERO/0.D+0/ DATA NEED(1)/1/, NEED(2)/0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) IV(COVREQ) = -IABS(IV(COVREQ)) IF (IV(COVREQ) .EQ. 0 .AND. IV(RDREQ) .GT. 0) IV(COVREQ) = -1 IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 I = (P-PS+2)*(P-PS+1)/2 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+3+I) CALL DRGLG(X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, RHOI, 1 RHOR, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + (P - PS + 3)*N IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N IV(NEXTV) = IV(J) + N*PS IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RD1 = IV(REGD0) R21 = RD1 - N RS1 = R21 - N RN = RS1 + N - 1 C 20 CALL DRGLG(V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, V(R1), 1 V(RD1), RHO, RHOI, RHOR, V, X) IF (IV(1)-2) 30, 50, 120 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) CALL CALCRJ(N, PS, X, NF, NEED, V(R1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 CALL DV7CPY(N, V(RS1), V(R1)) IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R *** C C *** INITIALIZE D IF NECESSARY *** C 50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO) 1 CALL DV7SCP(P, V(D1), ONE) C DK = D1 NG = IV(NGCALL) - 1 IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1 J1K0 = DR1 NF = IV(NFCALL) IF (NF .EQ. IV(NFGCAL)) GO TO 70 NG = NG + 1 CALL CALCRJ(N, PS, X, NF, NEED, V(RS1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 70 60 IV(TOOBIG) = 1 IV(NGCALL) = NG GO TO 20 70 DO 110 K = 1, PS XK = X(K) H = V(DLTFDJ) * MAX( ABS(XK), ONE/V(DK)) H0 = H DK = DK + 1 80 X(K) = XK + H NG = NG + 1 NF = -NG CALL CALCRJ(N, PS, X, NF, NEED, V(R21), V(DR1), UI, UR, UF) IF (NF .LT. 0) GO TO 90 H = NEGPT5 * H IF ( ABS(H/H0) .GE. HLIM) GO TO 80 GO TO 60 90 X(K) = XK IV(NGCALL) = NG I1 = R21 J1K = J1K0 J1K0 = J1K0 + 1 DO 100 I = RS1, RN V(J1K) = (V(I1) - V(I)) / H I1 = I1 + 1 J1K = J1K + PS 100 CONTINUE 110 CONTINUE GO TO 20 C 120 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 C 999 RETURN C C *** LAST LINE OF DGLF FOLLOWS *** END SUBROUTINE DRGLG(D, DR, IV, LIV, LV, N, ND, NN, P, PS, R, 1 RD, RHO, RHOI, RHOR, V, X) C C *** ITERATION DRIVER FOR GENERALIZED (NON)LINEAR MODELS (ETC.) C INTEGER LIV, LV, N, ND, NN, P, PS INTEGER IV(LIV), RHOI(*) DOUBLE PRECISION D(P), DR(ND,N), R(*), RD(*), RHOR(*), 1 V(LV), X(*) C DIMENSION RD(N, (P-PS)*(P-PS+1)/2 + 1) EXTERNAL RHO C C-------------------------- PARAMETER USAGE -------------------------- C C D....... SCALE VECTOR. C DR...... DERIVATIVES OF R AT X. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV... LIV MUST BE AT LEAST P + 90. C LV...... LENGTH OF V... LV MUST BE AT LEAST C 105 + P*(2*P+16) + 2*N + 4*PS. C N....... TOTAL NUMBER OF RESIDUALS. C ND...... LEADING DIMENSION OF DR -- MUST BE AT LEAST PS. C NN...... LEAD DIMENSION OF R, RD. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS. C R....... RESIDUALS (OR MEANS -- FUNCTIONS OF X) WHEN DRGLG IS CALLED C WITH IV(1) = 1. C RD...... RD(I) = HALF * (G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN C IV(RDREQ) IS 2, 3, 5, OR 6. DRGLG SETS IV(REGD) = 1 IF RD C IS SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE C TO COMPUTE IT, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN) C WAS INDEFINITE. BEFORE CONVERGENCE, RD IS ALSO USED AS C TEMPORARY STORAGE. C RHO..... COMPUTES INFO ABOUT OBJECTIVE FUNCTION. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C V....... FLOATING-POINT VALUES ARRAY. C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C C *** CALLING SEQUENCE FOR RHO... C C CALL RHO(NEED, F, N, NF, XN, R, RD, RHOI, RHOR, W) C C PARAMETER DECLARATIONS FOR RHO... C C INTEGER NEED(2), N, NF, RHOI(*) C FLOATING-POINT F, XN(*), R(*), RD(N,*), RHOR(*), W(N) C C RHOI AND RHOR ARE FOR RHO TO USE AS IT SEES FIT. THEY ARE PASSED C TO RHO WITHOUT CHANGE. IF IV(RDREQ) IS AT LEAST 4, I.E., IF MORE C THAN THE SIMPLEST REGRESSION DIAGNOSTIC INFORMATION IS TO BE COMPUTED, C THEN SOME COMPONENTS OF RHOI AND RHOR MUST CONVEY SOME EXTRA C DETAILS, AS DESCRIBED BELOW. C F, R, RD, AND W ARE EXPLAINED BELOW WITH NEED. C XN IS THE VECTOR OF NUISANCE PARAMETERS (OF LENGTH P - PS). IF C RHO NEEDS TO KNOW THE LENGTH OF XN, THEN THIS LENGTH SHOULD BE C COMMUNICATED THROUGH RHOI (OR THROUGH COMMON). RHO SHOULD NOT CHANGE C XN. C NEED(1) = 1 MEANS RHO SHOULD SET F TO THE SUM OF THE LOSS FUNCTION C VALUES AT THE RESIDUALS R(I). NF IS THE CURRENT FUNCTION INVOCATION C COUNT (A VALUE THAT IS INCREMENTED EACH TIME A NEW PARAMETER EXTIMATE C X IS CONSIDERED). NEED(2) IS THE VALUE NF HAD AT THE LAST R WHERE C RHO MIGHT BE CALLED WITH NEED(1) = 2. IF RHO SAVES INTERMEDIATE C RESULTS FOR USE IN CALLS WITH NEED(1) = 2, THEN IT CAN USE NF TO TELL C WHICH INTERMEDIATE RESULTS ARE APPROPRIATE, AND IT CAN SAVE SOME OF C THESE RESULTS IN R. C NEED(1) = 2 MEANS RHO SHOULD SET R(I) TO THE LOSS FUNCTION C DERIVATIVE WITH RESPECT TO THE RESIDUALS THAT WERE PASSED TO RHO WHEN C NF HAD THE SAME VALUE IT DOES NOW (AND NEED(1) WAS 1). RHO SHOULD C ALSO SET W(I) TO THE APPROXIMATION OF THE SECOND DERIVATIVE OF THE C LOSS FUNCTION (WITH RESPECT TO THE I-TH RESIDUAL) THAT SHOULD BE USED C IN THE GAUSS-NEWTON MODEL. WHEN THERE ARE NUISANCE PARAMETERS (I.E., C WHEN PS .LT. P) RHO SHOULD ALSO SET R(I+K*N) TO THE DERIVATIVE OF THE C LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL AND XN(K), AND IT C SHOULD SET RD(I,J + K*(K+1)/2 + 1) TO THE SECOND PARTIAL DERIVATIVE C OF THE I-TH RESIDUAL WITH RESPECT TO XN(J) AND XN(K), 0 .LE. J .LE. K C AND 1 .LE. K .LE. P - PS, WHERE XN(0) MEANS THE I-TH RESIDUAL ITSELF. C IN ANY EVENT, RHO SHOULD ALSO SET RD(I,1) TO THE (TRUE) SECOND C DERIVATIVE OF THE LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL. C NF (THE FUNCTION INVOCATION COUNT WHOSE NORMAL USE IS EXPLAINED C ABOVE) SHOULD NOT BE CHANGED UNLESS RHO CANNOT CARRY OUT THE REQUESTED C TASK, IN WHICH CASE RHO SHOULD SET NF TO 0. C C C *** REGRESSION DIAGNOSTICS *** C C IV(RDREQ) INDICATES WHETHER A COVARIANCE MATRIX AND REGRESSION C DIAGNOSTIC VECTOR ARE TO BE COMPUTED. IV(RDREQ) HAS THE FORM C IV(RDREQ) = CVR +2*RDR, WHERE CVR = 0 OR 1 AND RDR = 0, 1, OR 2, C SO THAT C C CVR = MOD(IV(RDREQ), 2) C RDR = MOD(IV(RDREQ)/2, 3). C C CVR = 0 FOR NO COVARIANCE MATRIX C = 1 IF A COVARIANCE MATRIX ESTIMATE IS DESIRED C C RDR = 0 FOR NO LEAVE-ONE-OUT DIAGNOSTIC INFORMATION. C = 1 TO HAVE ONE-STEP ESTIMATES OF F(X(I)) - F(X*) STORED IN RD, C WHERE X(I) MINIMIZES F (THE OBJECTIVE FUNCTION) WITH C COMPONENT I OF R REMOVED AND X* MINIMIZES THE FULL F. C = 2 FOR MORE DETAILED ONE-STEP LEAVE-ONE-OUT INFORMATION, AS C DICTATED BY THE IV COMPONENTS DESCRIBED BELOW. C C FOR RDR = 2, THE FOLLOWING COMPONENTS OF IV ARE RELEVANT... C C NFIX = IV(83) = NUMBER OF TRAILING NUISANCE PARAMETERS TO TREAT C AS FIXED WHEN COMPUTING DIAGNOSTIC VECTORS (0 .LE. NFIX .LE. C P - PS, SO X(I) IS KEPT FIXED FOR P - NFIX .LT. I .LE. P). C C LOO = IV(84) TELLS WHAT TO LEAVE OUT... C = 1 MEANS LEAVE OUT EACH COMPONENT OF R SEPARATELY, AND C = 2 MEANS LEAVE OUT CONTIGUOUS BLOCKS OF R COMPONENTS. C FOR LOO = 2, IV(85) IS THE STARTING SUBSCRIPT IN RHOI C OF AN ARRAY BS OF BLOCK SIZES, IV(86) IS THE STRIDE FOR BS, C AND IV(87) = NB IS THE NUMBER OF BLOCKS, SO THAT C BS(I) = RHOI(IV(85) + (I-1)*IV(86)), 1 .LE. I .LE. NB. C NOTE THAT IF ALL BLOCKS ARE THE SAME SIZE, THEN IT SUFFICES C TO SET RHOI(IV(85)) = BLOCKSIZE AND IV(86) = 0. C NOTE THAT LOO = 1 IS EQUIVALENT TO LOO = 2 WITH C RHOI(IV(85)) = 1, IV(86) = 0, IV(87) = N. C = 3,4 ARE SIMILAR TO LOO = 1,2, RESPECTIVELY, BUT LEAVING A C FRACTION OUT. IN THIS CASE, IV(88) IS THE STARTING C SUBSCRIPT IN RHOR OF AN ARRAY FLO OF FRACTIONS TO LEAVE OUT, C AND IV(89) IS THE STRIDE FOR FLO... C FLO(I) = RHOR(IV(88) + (I-1)*IV(89)), 1 .LE. I .LE. NB. C C XNOTI = IV(90) TELLS WHAT DIAGNOSTIC INFORMATION TO STORE... C = 0 MEANS JUST STORE ONE-STEP ESTIMATES OF F(X(I)) - F(X*) IN C RD(I), 1 .LE. I .LE. NB. C .GT. 0 MEANS ALSO STORE ONE-STEP ESTIMATES OF X(I) ESTIMATES C IN RHOR, STARTING AT RHOR(XNOTI)... C X(I)(J) = RHOR((I-1)*(P-NFIX) + J + XNOTI-1), C 1 .LE. I .LE. NB, 1 .LE. J .LE. P - NFIX. C C SOMETIMES ONE-STEP ESTIMATES OF X(I) DO NOT EXIST, BECAUSE THE C APPROXIMATE UPDATED HESSIAN MATRIX IS INDEFINITE. IN SUCH CASES, C THE CORRESPONDING RD COMPONENT IS SET TO -1, AND, IF XNOTI IS C POSITIVE, THE SOLUTION X IS RETURNED AS X(I). WHEN ONE-STEP ESTIMATES C OF X(I) DO EXIST, THE CORRESPONDING COMPONENT OF RD IS POSITIVE. C C SUMMARY OF RHOI COMPONENTS (FOR RDR = MOD(IV(RDREQ)/2, 3) = 2)... C C IV(83) = NFIX C IV(84) = LOO C IV(85) = START IN RHOI OF BS C IV(86) = STRIDE FOR BS C IV(87) = NB C IV(88) = START IN RHOR OF FLO C IV(89) = STRIDE FOR FLO C IV(90) = XNOTI (START IN RHOR OF X(I)). C C C *** COVARIANCE MATRIX ESTIMATE *** C C IF IV(RDREQ) INDICATES THAT A COVARIANCE MATRIX IS TO BE COMPUTED, C THEN IV(COVREQ) = IV(15) DETERMINES THE FORM OF THE COMPUTED C COVARIANCE MATRIX ESTIMATE AND, SIMULTANEOUSLY, THE FORM OF C APPROXIMATE HESSIAN MATRIX USED IN COMPUTING REGRESSION DIAGNOSTIC C INFORMATION. IN ALL CASES, SOME APPROXIMATE FINAL HESSIAN MATRIX C IS OBTAINED, AND ITS INVERSE IS THE COVARIANCE MATRIX ESTIMATE C (WHICH MAY HAVE TO BE SCALED APPROPRIATELY -- THAT IS UP TO YOU). C IF IV(COVREQ) IS AT MOST 2 IN ABSOLUTE VALUE, THEN THE FINAL C HESSIAN APPROXIMATION IS COMPUTED BY FINITE DIFFERENCES -- GRADIENT C DIFFERENCES IF IV(COVREQ) IS NONNEGATIVE, FUNCTION DIFFERENCES C OTHERWISE. IF (IV(COVREQ)) IS AT LEAST 3 IN ABSOLUTE VALUE, THEN THE C CURRENT GAUSS-NEWTON HESSIAN APPROXIMATION IS TAKEN AS THE FINAL C HESSIAN APPROXIMATION. FOR SOME PROBLEMS THIS SAVES TIME AND YIELDS C THE SAME OR NEARLY THE SAME HESSIAN APPROXIMATION AS DO FINITE C DIFFERENCES. FOR OTHER PROBLEMS, THE TWO KINDS OF HESSIAN C APPROXIMATIONS MAY GIVE DECIDEDLY DIFFERENT REGRESSION DIAGNOSTICS AND C COVARIANCE MATRIX ESTIMATES. C C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL DD7UP5,DIVSET, DG2LRD, DN3RDP, DD7TPR, DQ7ADR, DVSUM, 1 DG7LIT,DITSUM, DL7NVR, DL7ITV, DL7IVM,DL7SRT, DL7SQR, 2 DL7SVX, DL7SVN, DL7TSQ,DL7VML,DO7PRD,DV2AXY,DV7CPY, 3 DV7SCL, DV7SCP DOUBLE PRECISION DD7TPR, DL7SVX, DL7SVN,DVSUM C C DD7UP5... UPDATES SCALE VECTOR D. C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DG2LRD.... COMPUTES REGRESSION DIAGNOSTIC. C DN3RDP... PRINTS REGRESSION DIAGNOSTIC. C DD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C DQ7ADR.... ADDS ROWS TO QR FACTORIZATION. C DVSUM..... RETURNS SUM OF ELEMENTS OF A VECTOR. C DG7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM. C DITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. C DL7NVR... INVERTS COMPACTLY STORED TRIANGULAR MATRIX. C DL7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C DL7IVM... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C DL7SQR... COMPUTES L*(L**T) FOR LOWER TRIANG. MATRIX L. C DL7SVX... UNDERESTIMATES LARGEST SINGULAR VALUE OF TRIANG. MATRIX. C DL7SVN... OVERESTIMATES SMALLEST SINGULAR VALUE OF TRIANG. MATRIX. C DL7TSQ... COMPUTES (L**T)*L FOR LOWER TRIANG. MATRIX L. C DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DO7PRD.... ADDS OUTER PRODUCT OF VECTORS TO A MATRIX. C DV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCL... MULTIPLIES A VECTOR BY A SCALAR. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C C *** LOCAL VARIABLES *** C LOGICAL JUSTG, UPDATD, ZEROG INTEGER G1, HN1, I, II, IV1, J, J1, JTOL1, K, L, LH, 1 NEED1(2), NEED2(2), PMPS, PS1, PSLEN, QTR1, 2 RMAT1, STEP1, TEMP1, TEMP2, TEMP3, TEMP4, W, WI, Y1 DOUBLE PRECISION RHMAX, RHTOL, RHO1, RHO2, T C DOUBLE PRECISION ONE, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COVMAT, DINIT, DTYPE, DTINIT, D0INIT, F, 1 F0, FDH, G, H, HC, IPIVOT, IVNEED, JCN, JTOL, LMAT, 2 MODE, NEXTIV, NEXTV, NF0, NF1, NFCALL, NFCOV, NFGCAL, 3 NGCALL, NGCOV, PERM, QTR, RDREQ, REGD, RESTOR, 4 RMAT, RSPTOL, STEP, TOOBIG, VNEED, XNOTI, Y C C *** IV SUBSCRIPT VALUES *** C PARAMETER (CNVCOD=55, COVMAT=26, DTYPE=16, F0=13, FDH=74, G=28, 1 H=56, HC=71, IPIVOT=76, IVNEED=3, JCN=66, JTOL=59, 2 LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, NFCALL=6, 3 NFCOV=52, NF0=68, NF1=69, NFGCAL=7, NGCALL=30, 4 NGCOV=53, PERM=58, QTR=77, RESTOR=9, RMAT=78, RDREQ=57, 5 REGD=67, STEP=40, TOOBIG=2, VNEED=4, XNOTI=90, Y=48) C C *** V SUBSCRIPT VALUES *** C PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RSPTOL=49) PARAMETER (ONE=1.D+0, ZERO=0.D+0) SAVE NEED1, NEED2 DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C LH = P * (P+1) / 2 IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) PS1 = PS + 1 IV1 = IV(1) IF (IV1 .GT. 2) GO TO 10 W = IV(Y) + P IV(RESTOR) = 0 I = IV1 + 2 IF (IV(TOOBIG) .EQ. 0) GO TO (120, 110, 110, 130), I V(F) = V(F0) IF (I .NE. 3) IV(1) = 2 GO TO 40 C C *** FRESH START OR RESTART -- CHECK INPUT INTEGERS *** C 10 IF (ND .LT. PS) GO TO 360 IF (PS .GT. P) GO TO 360 IF (PS .LE. 0) GO TO 360 IF (N .LE. 0) GO TO 360 IF (IV1 .EQ. 14) GO TO 30 IF (IV1 .GT. 16) GO TO 420 IF (IV1 .LT. 12) GO TO 40 IF (IV1 .EQ. 12) IV(1) = 13 IF (IV(1) .NE. 13) GO TO 20 IV(IVNEED) = IV(IVNEED) + P IV(VNEED) = IV(VNEED) + P*(P+13)/2 + 2*N + 4*PS C *** ADJUST IV(PERM) TO MAKE ROOM FOR IV INPUT COMPONENTS C *** NEEDED WHEN IV(RDREQ) IS 4 OR 5... I = XNOTI + 1 IF (IV(PERM) .LT. I) IV(PERM) = I C 20 CALL DG7LIT(D, X, IV, LIV, LV, P, PS, V, X, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(IPIVOT) = IV(NEXTIV) IV(NEXTIV) = IV(IPIVOT) + P IV(Y) = IV(NEXTV) IV(G) = IV(Y) + P + N IV(RMAT) = IV(G) + P + 4*PS IV(QTR) = IV(RMAT) + LH IV(JTOL) = IV(QTR) + P + N IV(JCN) = IV(JTOL) + 2*P IV(NEXTV) = IV(JCN) + P IF (IV1 .EQ. 13) GO TO 999 C 30 JTOL1 = IV(JTOL) IF (V(DINIT) .GE. ZERO) CALL DV7SCP(P, D, V(DINIT)) IF (V(DTINIT) .GT. ZERO) CALL DV7SCP(P, V(JTOL1), V(DTINIT)) I = JTOL1 + P IF (V(D0INIT) .GT. ZERO) CALL DV7SCP(P, V(I), V(D0INIT)) IV(NF0) = 0 IV(NF1) = 0 C 40 G1 = IV(G) Y1 = IV(Y) CALL DG7LIT(D, V(G1), IV, LIV, LV, P, PS, V, X, V(Y1)) IF (IV(1) - 2) 50, 60, 380 C 50 V(F) = ZERO IF (IV(NF1) .EQ. 0) GO TO 999 IF (IV(RESTOR) .NE. 2) GO TO 999 IV(NF0) = IV(NF1) CALL DV7CPY(N, RD, R) IV(REGD) = 0 GO TO 999 C 60 IF (IV(MODE) .GT. 0) GO TO 370 CALL DV7SCP(P, V(G1), ZERO) RMAT1 = IABS(IV(RMAT)) QTR1 = IABS(IV(QTR)) CALL DV7SCP(PS, V(QTR1), ZERO) IV(REGD) = 0 CALL DV7SCP(PS, V(Y1), ZERO) CALL DV7SCP(LH, V(RMAT1), ZERO) IF (IV(RESTOR) .NE. 3) GO TO 70 CALL DV7CPY(N, R, RD) IV(NF1) = IV(NF0) 70 CALL RHO(NEED2, T, N, IV(NFGCAL), X(PS1), R, RD, RHOI, RHOR, V(W)) IF (IV(NFGCAL) .GT. 0) GO TO 90 80 IV(TOOBIG) = 1 GO TO 40 90 IF (IV(MODE) .LT. 0) GO TO 999 DO 100 I = 1, N 100 CALL DV2AXY(PS, V(Y1), R(I), DR(1,I), V(Y1)) GO TO 999 C C *** COMPUTE F(X) *** C 110 I = IV(NFCALL) NEED1(2) = IV(NFGCAL) CALL RHO(NEED1, V(F), N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IV(NF1) = I IF (I .LE. 0) GO TO 80 GO TO 40 C C *** COMPUTE GRADIENT INFORMATION FOR FINITE-DIFFERENCE HESSIAN *** C 120 IV(1) = 2 JUSTG = .TRUE. I = IV(NFCALL) CALL RHO(NEED1, T, N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IF (I .LE. 0) GO TO 80 CALL RHO(NEED2, T, N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IF (I .LE. 0) GO TO 80 GO TO 250 C C *** PREPARE TO COMPUTE GRADIENT INFORMATION WHILE ITERATING *** C 130 JUSTG = .FALSE. G1 = IV(G) C C *** DECIDE WHETHER TO UPDATE D BELOW *** C I = IV(DTYPE) UPDATD = .FALSE. IF (I .LE. 0) GO TO 140 IF (I .EQ. 1 .OR. IV(MODE) .LT. 0) UPDATD = .TRUE. C C *** COMPUTE RMAT AND QTR *** C 140 QTR1 = IABS(IV(QTR)) RMAT1 = IABS(IV(RMAT)) IV(RMAT) = RMAT1 IV(HC) = 0 IV(NF0) = 0 IV(NF1) = 0 IF (IV(MODE) .LT. 0) GO TO 160 C C *** ADJUST Y *** C Y1 = IV(Y) WI = W STEP1 = IV(STEP) DO 150 I = 1, N T = V(WI) - RD(I) WI = WI + 1 IF (T .NE. ZERO) CALL DV2AXY(PS, V(Y1), 1 T*DD7TPR(PS,V(STEP1),DR(1,I)), DR(1,I), V(Y1)) 150 CONTINUE C C *** CHECK FOR NEGATIVE W COMPONENTS *** C 160 J1 = W + N - 1 DO 170 WI = W, J1 IF (V(WI) .LT. ZERO) GO TO 240 170 CONTINUE C C *** W IS NONNEGATIVE. COMPUTE QR FACTORIZATION *** C *** AND, IF NECESSARY, USE SEMINORMAL EQUATIONS *** C RHMAX = ZERO RHTOL = V(RSPTOL) TEMP1 = G1 + P ZEROG = .TRUE. WI = W DO 200 I = 1, N RHO1 = R(I) RHO2 = V(WI) WI = WI + 1 T = SQRT(RHO2) IF (RHMAX .LT. RHO2) RHMAX = RHO2 IF (RHO2 .GT. RHTOL*RHMAX) GO TO 180 C *** SEMINORMAL EQUATIONS *** CALL DV2AXY(PS, V(G1), RHO1, DR(1,I), V(G1)) RHO1 = ZERO ZEROG = .FALSE. GO TO 190 180 RHO1 = RHO1 / T C *** QR ACCUMULATION *** 190 CALL DV7SCL(PS, V(TEMP1), T, DR(1,I)) CALL DQ7ADR(PS, V(QTR1), V(RMAT1), V(TEMP1), RHO1) 200 CONTINUE C C *** COMPUTE G FROM RMAT AND QTR *** C TEMP2 = TEMP1 + PS CALL DL7VML(PS, V(TEMP1), V(RMAT1), V(QTR1)) IF (ZEROG) GO TO 220 IV(QTR) = -QTR1 IF (DL7SVX(PS, V(RMAT1), V(TEMP2), V(TEMP2)) * RHTOL .GE. 1 DL7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2))) GO TO 230 CALL DL7IVM(PS, V(TEMP2), V(RMAT1), V(G1)) C C *** SEMINORMAL EQUATIONS CORRECTION OF BJOERCK -- C *** ONE CYCLE OF ITERATIVE REFINEMENT... C TEMP3 = TEMP2 + PS TEMP4 = TEMP3 + PS CALL DL7ITV(PS, V(TEMP3), V(RMAT1), V(TEMP2)) CALL DV7SCP(PS, V(TEMP4), ZERO) RHMAX = ZERO WI = W DO 210 I = 1, N RHO2 = V(WI) WI = WI + 1 IF (RHMAX .LT. RHO2) RHMAX = RHO2 RHO1 = ZERO IF (RHO2 .LE. RHTOL*RHMAX) RHO1 = R(I) T = RHO1 - RHO2*DD7TPR(PS, V(TEMP3), DR(1,I)) CALL DV2AXY(PS, V(TEMP4), T, DR(1,I), V(TEMP4)) 210 CONTINUE CALL DL7IVM(PS, V(TEMP3), V(RMAT1), V(TEMP4)) CALL DV2AXY(PS, V(TEMP2), ONE, V(TEMP3), V(TEMP2)) CALL DV2AXY(PS, V(QTR1), ONE, V(TEMP2), V(QTR1)) 220 IV(QTR) = QTR1 230 CALL DV2AXY(PS, V(G1), ONE, V(TEMP1), V(G1)) IF (PS .GE. P) GO TO 350 GO TO 270 C C *** INDEFINITE GN HESSIAN... *** C 240 IV(RMAT) = -RMAT1 IV(HC) = RMAT1 CALL DO7PRD(N, LH, PS, V(RMAT1), V(W), DR, DR) C C *** COMPUTE GRADIENT *** C 250 G1 = IV(G) CALL DV7SCP(P, V(G1), ZERO) DO 260 I = 1, N 260 CALL DV2AXY(PS, V(G1), R(I), DR(1,I), V(G1)) IF (PS .GE. P) GO TO 350 C C *** COMPUTE GRADIENT COMPONENTS OF NUISANCE PARAMETERS *** C 270 K = P - PS J1 = 1 G1 = G1 + PS DO 280 J = 1, K J1 = J1 + NN V(G1) =DVSUM(N, R(J1)) G1 = G1 + 1 280 CONTINUE IF (JUSTG) GO TO 390 C C *** COMPUTE HESSIAN COMPONENTS OF NUISANCE PARAMETERS *** C I = PS*PS1/2 PSLEN = P*(P+1)/2 - I HN1 = RMAT1 + I CALL DV7SCP(PSLEN, V(HN1), ZERO) PMPS = P - PS K = HN1 J1 = 1 DO 310 II = 1, PMPS J1 = J1 + NN J = J1 DO 290 I = 1, N CALL DV2AXY(PS, V(K), RD(J), DR(1,I), V(K)) J = J + 1 290 CONTINUE K = K + PS DO 300 I = 1, II J1 = J1 + NN V(K) =DVSUM(N, RD(J1)) K = K + 1 300 CONTINUE 310 CONTINUE IF (IV(RMAT) .LE. 0) GO TO 350 J = IV(LMAT) CALL DV7CPY(PSLEN, V(J), V(HN1)) IF (DL7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2)) .LE. ZERO) GO TO 320 CALL DL7SRT(PS1, P, V(RMAT1), V(RMAT1), I) IF (I .LE. 0) GO TO 330 C C *** HESSIAN IS NOT POSITIVE DEFINITE *** C 320 CALL DL7SQR(PS, V(RMAT1), V(RMAT1)) CALL DV7CPY(PSLEN, V(HN1), V(J)) IV(HC) = RMAT1 IV(RMAT) = -RMAT1 GO TO 350 C C *** NUISANCE PARS LEAVE HESSIAN POS. DEF. GET REST OF QTR *** C 330 J = QTR1 + PS G1 = IV(G) + PS DO 340 I = PS1, P T = DD7TPR(I-1, V(HN1), V(QTR1)) HN1 = HN1 + I V(J) = (V(G1) - T) / V(HN1-1) J = J + 1 G1 = G1 + 1 340 CONTINUE 350 IF (JUSTG) GO TO 390 IF (UPDATD) CALL DD7UP5(D, IV, LIV, LV, P, PS, V) GO TO 40 C C *** MISC. DETAILS *** C C *** BAD N, ND, OR P *** C 360 IV(1) = 66 GO TO 420 C C *** COVARIANCE OR INITIAL S COMPUTATION *** C 370 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(NFGCAL) = IV(NFCALL) IV(1) = -1 GO TO 999 C C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** C 380 IF (IV(COVMAT) .NE. 0) GO TO 410 IF (IV(REGD) .NE. 0) GO TO 410 C C *** SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE *** C K = IV(FDH) IF (K .LE. 0) GO TO 400 IF (IV(RDREQ) .LE. 0) GO TO 410 C C *** COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF C DESIRED *** C IV(MODE) = P + 1 IV(NGCALL) = IV(NGCALL) + 1 IV(NGCOV) = IV(NGCOV) + 1 IV(CNVCOD) = IV(1) IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(NFGCAL) = IV(NFCALL) IV(1) = -1 GO TO 999 C 390 IF (IV(MODE) .LE. P) GO TO 40 C *** SAVE RD IN W FOR POSSIBLE USE IN OTHER DIAGNOSTICS *** CALL DV7CPY(N, V(W), RD) C *** OVERWRITE RD WITH REGRESSION DIAGNOSTICS *** L = IV(LMAT) I = IV(JCN) STEP1 = IV(STEP) CALL DG2LRD(DR, IV, V(L), LH, LIV, LV, ND, N, P, PS, R, RD, 1 RHOI, RHOR, V, V(STEP1), X, V(I)) IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 IF (MOD(IV(RDREQ),2) .EQ. 0) GO TO 410 C C *** FINISH COVARIANCE COMPUTATION *** C I = IABS(IV(H)) IV(FDH) = 0 CALL DL7NVR(P, V(I), V(L)) CALL DL7TSQ(P, V(I), V(I)) IV(COVMAT) = I GO TO 410 C C *** COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN *** C 400 IV(COVMAT) = K IV(REGD) = K C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 410 G1 = IV(G) 420 CALL DITSUM(D, V(G1), IV, LIV, LV, P, V, X) IF (IV(1) .LE. 6 .AND. IV(RDREQ) .GT. 0) 1 CALL DN3RDP(IV, LIV, LV, N, P, RD, RHOI, RHOR, V) C 999 RETURN C *** LAST LINE OF DRGLG FOLLOWS *** END SUBROUTINE DF7HES(D, G, IRT, IV, LIV, LV, P, V, X) C C *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING C *** AT V(IV(FDH)) = V(-IV(H)). C C *** IF IV(COVREQ) .GE. 0 THEN DF7HES USES GRADIENT DIFFERENCES, C *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN DG7LIT. C C IRT VALUES... C 1 = COMPUTE FUNCTION VALUE, I.E., V(F). C 2 = COMPUTE G. C 3 = DONE. C C C *** PARAMETER DECLARATIONS *** C INTEGER IRT, LIV, LV, P INTEGER IV(LIV) DOUBLE PRECISION D(P), G(P), V(LV), X(P) C C *** LOCAL VARIABLES *** C INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2, 1 PP1O2, STPI, STPM, STP0 DOUBLE PRECISION DEL, HALF, NEGPT5, ONE, TWO, ZERO C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DV7CPY C C DV7CPY.... COPY ONE VECTOR TO ANOTHER. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H, KAGQT, MODE, 1 NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE C PARAMETER (HALF=0.5D+0, NEGPT5=-0.5D+0, ONE=1.D+0, TWO=2.D+0, 1 ZERO=0.D+0) C PARAMETER (COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42, F=10, 1 FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7, 2 SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C IRT = 4 KIND = IV(COVREQ) M = IV(MODE) IF (M .GT. 0) GO TO 10 IV(H) = -IABS(IV(H)) IV(FDH) = 0 IV(KAGQT) = -1 V(FX) = V(F) 10 IF (M .GT. P) GO TO 999 IF (KIND .LT. 0) GO TO 110 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND C *** GRADIENT VALUES. C GSAVE1 = IV(W) + P IF (M .GT. 0) GO TO 20 C *** FIRST CALL ON DF7HES. SET GSAVE = G, TAKE FIRST STEP *** CALL DV7CPY(P, V(GSAVE1), G) IV(SWITCH) = IV(NFGCAL) GO TO 90 C 20 DEL = V(DELTA) X(M) = V(XMSAVE) IF (IV(TOOBIG) .EQ. 0) GO TO 40 C C *** HANDLE OVERSIZE V(DELTA) *** C IF (DEL*X(M) .GT. ZERO) GO TO 30 C *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** IV(FDH) = -2 GO TO 220 C C *** TRY SHRINKING V(DELTA) *** 30 DEL = NEGPT5 * DEL GO TO 100 C 40 HES = -IV(H) C C *** SET G = (G - GSAVE)/DEL *** C DO 50 I = 1, P G(I) = (G(I) - V(GSAVE1)) / DEL GSAVE1 = GSAVE1 + 1 50 CONTINUE C C *** ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** C K = HES + M*(M-1)/2 L = K + M - 2 IF (M .EQ. 1) GO TO 70 C C *** SET H(I,M) = 0.5 * (H(I,M) + G(I)) FOR I = 1 TO M-1 *** C MM1 = M - 1 DO 60 I = 1, MM1 V(K) = HALF * (V(K) + G(I)) K = K + 1 60 CONTINUE C C *** ADD H(I,M) = G(I) FOR I = M TO P *** C 70 L = L + 1 DO 80 I = M, P V(L) = G(I) L = L + I 80 CONTINUE C 90 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 210 C C *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE *** C DEL = V(DELTA0) * MAX(ONE/D(M), ABS(X(M))) IF (X(M) .LT. ZERO) DEL = -DEL V(XMSAVE) = X(M) 100 X(M) = X(M) + DEL V(DELTA) = DEL IRT = 2 GO TO 999 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. C 110 STP0 = IV(W) + P - 1 MM1 = M - 1 MM1O2 = M*MM1/2 IF (M .GT. 0) GO TO 120 C *** FIRST CALL ON DF7HES. *** IV(SAVEI) = 0 GO TO 200 C 120 I = IV(SAVEI) HES = -IV(H) IF (I .GT. 0) GO TO 180 IF (IV(TOOBIG) .EQ. 0) GO TO 140 C C *** HANDLE OVERSIZE STEP *** C STPM = STP0 + M DEL = V(STPM) IF (DEL*X(XMSAVE) .GT. ZERO) GO TO 130 C *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** IV(FDH) = -2 GO TO 220 C C *** TRY SHRINKING THE STEP *** 130 DEL = NEGPT5 * DEL X(M) = X(XMSAVE) + DEL V(STPM) = DEL IRT = 1 GO TO 999 C C *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** C 140 PP1O2 = P * (P-1) / 2 HPM = HES + PP1O2 + MM1 V(HPM) = V(F) C C *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** C HMI = HES + MM1O2 IF (MM1 .EQ. 0) GO TO 160 HPI = HES + PP1O2 DO 150 I = 1, MM1 V(HMI) = V(FX) - (V(F) + V(HPI)) HMI = HMI + 1 HPI = HPI + 1 150 CONTINUE 160 V(HMI) = V(F) - TWO*V(FX) C C *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** C I = 1 C 170 IV(SAVEI) = I STPI = STP0 + I V(DELTA) = X(I) X(I) = X(I) + V(STPI) IF (I .EQ. M) X(I) = V(XMSAVE) - V(STPI) IRT = 1 GO TO 999 C 180 X(I) = V(DELTA) IF (IV(TOOBIG) .EQ. 0) GO TO 190 C *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** IV(FDH) = -2 GO TO 220 C C *** FINISH COMPUTING H(M,I) *** C 190 STPI = STP0 + I HMI = HES + MM1O2 + I - 1 STPM = STP0 + M V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM)) I = I + 1 IF (I .LE. M) GO TO 170 IV(SAVEI) = 0 X(M) = V(XMSAVE) C 200 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 210 C C *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. C *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN C *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. C DEL = V(DLTFDC) * MAX(ONE/D(M), ABS(X(M))) IF (X(M) .LT. ZERO) DEL = -DEL V(XMSAVE) = X(M) X(M) = X(M) + DEL STPM = STP0 + M V(STPM) = DEL IRT = 1 GO TO 999 C C *** RESTORE V(F), ETC. *** C 210 IV(FDH) = HES 220 V(F) = V(FX) IRT = 3 IF (KIND .LT. 0) GO TO 999 IV(NFGCAL) = IV(SWITCH) GSAVE1 = IV(W) + P CALL DV7CPY(P, G, V(GSAVE1)) GO TO 999 C 999 RETURN C *** LAST LINE OF DF7HES FOLLOWS *** END SUBROUTINE DG2LRD(DR, IV, L, LH, LIV, LV, ND, N, P, PS, R, RD, 1 RHOI, RHOR, V, W, X, Z) C C *** COMPUTE REGRESSION DIAGNOSTIC FOR DRGLG *** C C *** PARAMETERS *** C INTEGER LH, LIV, LV, ND, N, P, PS INTEGER IV(LIV), RHOI(*) DOUBLE PRECISION DR(ND,P), L(LH), R(N), RD(N), RHOR(*), V(LV), 1 W(P), X(P), Z(P) C C *** CODED BY DAVID M. GAY (SPRING 1986, SUMMER 1991) *** C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL DD7TPR, DL7ITV, DL7IVM,DL7SRT, DL7SQR, DS7LVM, 1 DV2AXY,DV7CPY, DV7SCP DOUBLE PRECISION DD7TPR C C DD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C DL7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C DL7IVM... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C DL7SQR... COMPUTES L*(L**T) FOR LOWER TRIANG. MATRIX L. C DS7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR. C DV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C C *** LOCAL VARIABLES *** C LOGICAL USEFLO INTEGER BS1, BSINC, FLO1, FLOINC, H1, HPS1, I, 1 J, J1, K, KI, KI1, KID, L1, LE, LL, LOO1, N1, 2 PMPS, PP1O2, PS1, PX, RDR, XNI, ZAP1, ZAPLEN DOUBLE PRECISION FRAC, HI, RI, S, T, T1 C C *** CONSTANTS *** C DOUBLE PRECISION HALF, NEGONE, ONE, ZERO C C C *** IV SUBSCRIPTS *** C INTEGER BS, BSSTR, COVREQ, FDH, FLO, FLOSTR, LOO, NB, NFIX, 1 RDREQ, REGD, XNOTI PARAMETER (BS=85, BSSTR=86, COVREQ=15, FDH=74, FLO=88, FLOSTR=89, 1 LOO=84, NB=87, NFIX=83, RDREQ=57, REGD=67, XNOTI=90) PARAMETER (HALF=0.5D+0, NEGONE=-1.D+0, ONE=1.D+0, ZERO=0.D+0) C C++++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++ C I = IV(RDREQ) RDR = MOD(I/2, 3) IF (RDR .EQ. 0) GO TO 999 H1 = IV(FDH) USEFLO = .FALSE. PX = P N1 = N FRAC = ONE XNI = 0 IF (RDR .EQ. 1) GO TO 120 LOO1 = IV(LOO) IF (LOO1 .LE. 0 .OR. LOO1 .GT. 6) THEN IV(REGD) = -1 GO TO 999 ENDIF IF (LOO1 .GT. 3) THEN USEFLO = .TRUE. FLO1 = IV(FLO) FLOINC = IV(FLOSTR) LOO1 = LOO1 - 3 ENDIF XNI = IV(XNOTI) PX = P - IV(NFIX) IF (PX .LT. PS .OR. PX .GT. P) THEN IV(REGD) = -2 GO TO 999 ENDIF IF (LOO1 .EQ. 1) GO TO 120 N1 = IV(NB) IF (N1 .LE. 0 .OR. N1 .GT. N) THEN IV(REGD) = -3 GO TO 999 ENDIF BS1 = IV(BS) BSINC = IV(BSSTR) IF (H1 .LE. 0) GO TO 190 IF (IABS(IV(COVREQ)) .GE. 3) CALL DL7SQR(P, V(H1), L) PP1O2 = PX*(PX+1)/2 PS1 = PS + 1 ZAP1 = PS*(PS1)/2 + 1 LE = 0 DO 100 I = 1, N1 IF (USEFLO) THEN FRAC = RHOR(FLO1) FLO1 = FLO1 + FLOINC ENDIF L1 = LE + 1 IF (L1 .GT. N) GO TO 110 LE = LE + RHOI(BS1) IF (LE .GT. N) LE = N BS1 = BS1 + BSINC CALL DV7CPY(PP1O2, L, V(H1)) IF (PS .GE. PX) GO TO 50 K = ZAP1 KI = L1 DO 40 J = PS1, P KI = KI + N KI1 = KI DO 10 LL = L1, LE CALL DV2AXY(PS, L(K), -FRAC*RD(KI1), DR(1,LL), L(K)) KI1 = KI1 + 1 10 CONTINUE K = K + PS DO 30 J1 = PS1, J KI = KI + N KI1 = KI T = ZERO DO 20 LL = L1, LE T = T + RD(KI1) KI1 = KI1 + 1 20 CONTINUE L(K) = L(K) - FRAC*T K = K + 1 30 CONTINUE 40 CONTINUE 50 DO 70 LL = L1, LE T = -FRAC*RD(LL) K = 1 DO 60 J = 1, PS CALL DV2AXY(J, L(K), T*DR(J,LL), DR(1,LL), L(K)) K = K + J 60 CONTINUE 70 CONTINUE CALL DL7SRT(1, PX, L, L, J) IF (J .EQ. 0) THEN CALL DV7SCP(PX, W, ZERO) DO 90 LL = L1, LE CALL DV2AXY(PS, W, R(LL), DR(1,LL), W) IF (PS1 .GT. PX) GO TO 90 K = L1 DO 80 J = PS1, P K = K + N W(J) = W(J) + R(K) 80 CONTINUE 90 CONTINUE CALL DL7IVM(PX, W, L, W) CALL DL7ITV(PX, W, L, W) CALL DS7LVM(PX, Z, V(H1), W) RD(I) = HALF * FRAC * DD7TPR(PX, W, Z) IF (XNI .GT. 0) THEN CALL DV2AXY(PX, RHOR(XNI), FRAC, W, X) XNI = XNI + PX ENDIF ELSE RD(I) = NEGONE IF (XNI .GT. 0) THEN CALL DV7CPY(PX, RHOR(XNI), X) XNI = XNI + PX ENDIF ENDIF 100 CONTINUE 110 IV(REGD) = 1 C *** RESTORE L *** CALL DL7SRT(1, P, L, V(H1), J) GO TO 999 C 120 IF (H1 .LE. 0) GO TO 190 IF (IABS(IV(COVREQ)) .GE. 3) CALL DL7SQR(P, V(H1), L) IF (PS .GE. PX) GO TO 170 PS1 = PS + 1 PMPS = PX - PS ZAP1 = PS*(PS1)/2 ZAPLEN = PX*(PX+1)/2 - ZAP1 HPS1 = H1 + ZAP1 ZAP1 = ZAP1 + 1 DO 160 I = 1, N IF (USEFLO) THEN FRAC = RHOR(FLO1) FLO1 = FLO1 + FLOINC ENDIF CALL DV7CPY(ZAPLEN, L(ZAP1), V(HPS1)) CALL DV7SCP(PS, W, ZERO) K = ZAP1 KI = I KID = KI DO 140 J = PS1, PX KI = KI + N CALL DV2AXY(PS, L(K), -FRAC*RD(KI), DR(1,I), L(K)) K = K + PS KID = KID + N W(J) = FRAC*R(KID) DO 130 J1 = PS1, J KI = KI + N L(K) = L(K) - FRAC*RD(KI) K = K + 1 130 CONTINUE 140 CONTINUE CALL DL7SRT(PS1, PX, L, L, J) IF (J .NE. 0) GO TO 150 CALL DV7CPY(PS, Z, DR(1,I)) CALL DV7SCP(PMPS, Z(PS1), ZERO) CALL DL7IVM(PX, Z, L, Z) HI = DD7TPR(PX, Z, Z) CALL DL7IVM(PX, W, L, W) RI = FRAC*R(I) C *** FIRST PS ELEMENTS OF W VANISH *** T = DD7TPR(PMPS, W(PS1), Z(PS1)) S = FRAC*RD(I) T1 = ONE - S*HI IF (T1 .LE. ZERO) GO TO 150 CALL DV2AXY(PX, W, (RI + S*T)/T1, Z, W) CALL DL7ITV(PX, W, L, W) CALL DS7LVM(PX, Z, V(H1), W) RD(I) = HALF * DD7TPR(PX, W, Z) IF (XNI .GT. 0) THEN CALL DV2AXY(PX, RHOR(XNI), ONE, W, X) XNI = XNI + PX ENDIF GO TO 160 150 RD(I) = NEGONE IF (XNI .GT. 0) THEN CALL DV7CPY(PX, RHOR(XNI), X) XNI = XNI + PX ENDIF 160 CONTINUE C C *** RESTORE L *** C CALL DV7CPY(ZAPLEN, L(ZAP1), V(HPS1)) CALL DL7SRT(PS1, PX, L, L, J) GO TO 200 C 170 DO 180 I = 1, N IF (USEFLO) THEN FRAC = RHOR(FLO1) FLO1 = FLO1 + FLOINC ENDIF CALL DL7IVM(PX, Z, L, DR(1,I)) S = DD7TPR(PX, Z, Z) T = ONE - FRAC*RD(I) * S IF (T .LE. ZERO) THEN RD(I) = NEGONE IF (XNI .GT. 0) THEN CALL DV7CPY(PX, RHOR(XNI), X) XNI = XNI + PX ENDIF ELSE RD(I) = HALF * FRAC * (R(I)/T)**2 * S IF (XNI .GT. 0) THEN CALL DL7ITV(PX, Z, L, Z) CALL DV2AXY(PX, RHOR(XNI), FRAC*R(I)/T, Z, X) XNI = XNI + PX ENDIF ENDIF 180 CONTINUE GO TO 200 C 190 CALL DV7SCP(N1, RD, NEGONE) 200 IV(REGD) = 1 C 999 RETURN C *** LAST LINE OF DG2LRD FOLLOWS *** END SUBROUTINE DG7LIT(D, G, IV, LIV, LV, P, PS, V, X, Y) C C *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR *** C *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) *** C C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, P, PS INTEGER IV(LIV) DOUBLE PRECISION D(P), G(P), V(LV), X(P), Y(P) C C-------------------------- PARAMETER USAGE -------------------------- C C D.... SCALE VECTOR. C IV... INTEGER VALUE ARRAY. C LIV.. LENGTH OF IV. MUST BE AT LEAST 82. C LH... LENGTH OF H = P*(P+1)/2. C LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7. C G.... GRADIENT AT X (WHEN IV(1) = 2). C P.... NUMBER OF PARAMETERS (COMPONENTS IN X). C PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. C V.... FLOATING-POINT VALUE ARRAY. C X.... PARAMETER VECTOR. C Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). C C *** DISCUSSION *** C C DG7LIT PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF C REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES C IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED C FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES C THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED C COMPACTLY BY ROWS IN V, STARTING AT IV(HC)), AND DG7LIT BUILDS AN C APPROXIMATION, S, TO THE SECOND-ORDER TERM. THE CALLER ALSO C PROVIDES THE FUNCTION VALUE, GRADIENT, AND PART OF THE YIELD C VECTOR USED IN UPDATING S. DG7LIT DECIDES DYNAMICALLY WHETHER OR C NOT TO USE S WHEN CHOOSING THE NEXT STEP TO TRY... THE HESSIAN C APPROXIMATION USED IS EITHER HC ALONE (GAUSS-NEWTON MODEL) OR C HC + S (AUGMENTED MODEL). C C IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT C CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO C 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS C IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN C COMPUTED HAS NONZERO VALUES IN THESE ROWS. C C IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY C FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS C USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME C WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2, C 1, OR 2). C C FOR UPDATING S,DG7LIT ASSUMES THAT THE GRADIENT HAS THE FORM C OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE C GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS C THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP, C THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF C RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY C PART OF THIS IN Y, NAMELY THE SUM OVER I OF C RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING DG7LIT WITH IV(1) = 2 AND C IV(MODE) = 0 (WHERE MODE = 38). G THEN CONTANS THE OTHER PART, C SO THAT THE DESIRED YIELD VECTOR IS G - Y. IF PS .LT. P, THEN C THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF C GRAD(R(I,X)), STEP, AND Y. C C PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING C ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER C (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS C NOT NEEDED). MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE C TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, C AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). THE VALUES IV(D), C IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND C NL2SNO), ARE NOT REFERENCED BY DG7LIT OR THE SUBROUTINES IT CALLS. C C WHEN DG7LIT IS FIRST CALLED, I.E., WHEN DG7LIT IS CALLED WITH C IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED. TO C OBTAIN THESE STARTING VALUES,DG7LIT RETURNS FIRST WITH IV(1) = 1, C THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES. ON C SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT C Y MUST ALSO BE SUPPLIED. (NOTE THAT Y IS USED FOR SCRATCH -- ITS C INPUT CONTENTS ARE LOST. BY CONTRAST, HC IS NEVER CHANGED.) C ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY C IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE C IN COMPUTING A COVARIANCE MATRIX. IN THIS CASE DG7LIT WILL MAKE A C NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. C WHEN IV(MODE) IS POSITIVE, Y SHOULD NOT BE CHANGED. C C IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE C FUNCTION VALUE AT X, AND CALL DG7LIT AGAIN, HAVING CHANGED C NONE OF THE OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X) C CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH C MAY HAPPEN BECAUSE OF AN OVERSIZED STEP. IN THIS CASE C THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL C CAUSE DG7LIT TO IGNORE V(F) AND TRY A SMALLER STEP. NOTE C THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE C IN IV(NFCALL) = IV(6). THIS MAY BE USED TO IDENTIFY C WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM- C PUTING G, HC, AND Y THE NEXT TIME DG7LIT RETURNS WITH C IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT OF F AT C X. THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON C HESSIAN AT X. IF IV(MODE) = 0, THEN THE CALLER SHOULD C ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE. C THE CALLER SHOULD THEN CALL DG7LIT AGAIN (WITH IV(1) = 2). C THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT C CHANGE X. NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE C VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH C IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS. C IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1. MLPIT C IS AN EXAMPLE WHERE THIS INFORMATION IS USED. IF G OR HC C CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET C IV(TOOBIG) TO 1, IN WHICH CASE DG7LIT WILL RETURN WITH C IV(1) = 15. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS. C C (SEE NL2SOL FOR REFERENCES.) C C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER DUMMY, DIG1, G01, H1, HC1, I, IPIV1, J, K, L, LMAT1, 1 LSTGST, PP1O2, QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1, 2 TEMP1, TEMP2, W1, X01 DOUBLE PRECISION E, STTSST, T, T1 C C *** CONSTANTS *** C DOUBLE PRECISION HALF, NEGONE, ONE, ONEP2, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C LOGICAL STOPX DOUBLE PRECISION DD7TPR, DL7SVX, DL7SVN, DRLDST, DR7MDC, DV2NRM EXTERNAL DA7SST, DD7TPR,DF7HES,DG7QTS,DITSUM, DL7MST,DL7SRT, 1 DL7SQR, DL7SVX, DL7SVN, DL7TVM,DL7VML,DPARCK, DRLDST, 2 DR7MDC, DS7LUP, DS7LVM, STOPX,DV2AXY,DV7CPY, DV7SCP, 3 DV2NRM C C DA7SST.... ASSESSES CANDIDATE STEP. C DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. C DF7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE). C DG7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. C DL7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). C DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C DL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. C DL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. C DL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. C DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. C DRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. C DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. C DS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- C ANGLE OF A SYMMETRIC MATRIX. C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. C DV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C DV2NRM... RETURNS THE 2-NORM OF A VECTOR. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COSMIN, COVMAT, COVREQ, DGNORM, DIG, DSTNRM, F, 1 FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR, INCFAC, INITS, 2 IPIVOT, IRC, KAGQT, KALM, LMAT, LMAX0, LMAXS, MODE, MODEL, 3 MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL, NFCOV, NGCOV, 4 NGCALL, NITER, NVSAVE, PHMXFC, PREDUC, QTR, RADFAC, 5 RADINC, RADIUS, RAD0, RCOND, RDREQ, REGD, RELDX, RESTOR, 6 RMAT, S, SIZE, STEP, STGLIM, STLSTG, STPPAR, SUSED, 7 SWITCH, TOOBIG, TUNER4, TUNER5, VNEED, VSAVE, W, WSCALE, 8 XIRC, X0 C C *** IV SUBSCRIPT VALUES *** C PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DIG=37, FDH=74, H=56, 1 HC=71, IERR=75, INITS=25, IPIVOT=76, IRC=29, KAGQT=33, 2 KALM=34, LMAT=42, MODE=35, MODEL=5, MXFCAL=17, 3 MXITER=18, NEXTV=47, NFCALL=6, NFGCAL=7, NFCOV=52, 4 NGCOV=53, NGCALL=30, NITER=31, QTR=77, RADINC=8, 5 RDREQ=57, REGD=67, RESTOR=9, RMAT=78, S=62, STEP=40, 6 STGLIM=11, STLSTG=41, SUSED=64, SWITCH=12, TOOBIG=2, 7 VNEED=4, VSAVE=60, W=65, XIRC=13, X0=43) C C *** V SUBSCRIPT VALUES *** C PARAMETER (COSMIN=47, DGNORM=1, DSTNRM=2, F=10, FDIF=11, FUZZ=45, 1 F0=13, GTSTEP=4, INCFAC=23, LMAX0=35, LMAXS=36, 2 NVSAVE=9, PHMXFC=21, PREDUC=7, RADFAC=16, RADIUS=8, 3 RAD0=9, RCOND=53, RELDX=17, SIZE=55, STPPAR=5, 4 TUNER4=29, TUNER5=30, WSCALE=56) C C PARAMETER (HALF=0.5D+0, NEGONE=-1.D+0, ONE=1.D+0, ONEP2=1.2D+0, 1 ZERO=0.D+0) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C I = IV(1) IF (I .EQ. 1) GO TO 40 IF (I .EQ. 2) GO TO 50 C IF (I .EQ. 12 .OR. I .EQ. 13) 1 IV(VNEED) = IV(VNEED) + P*(3*P + 19)/2 + 7 CALL DPARCK(1, D, IV, LIV, LV, P, V) I = IV(1) - 2 IF (I .GT. 12) GO TO 999 GO TO (290, 290, 290, 290, 290, 290, 170, 120, 170, 10, 10, 20), I C C *** STORAGE ALLOCATION *** C 10 PP1O2 = P * (P + 1) / 2 IV(S) = IV(LMAT) + PP1O2 IV(X0) = IV(S) + PP1O2 IV(STEP) = IV(X0) + P IV(STLSTG) = IV(STEP) + P IV(DIG) = IV(STLSTG) + P IV(W) = IV(DIG) + P IV(H) = IV(W) + 4*P + 7 IV(NEXTV) = IV(H) + PP1O2 IF (IV(1) .NE. 13) GO TO 20 IV(1) = 14 GO TO 999 C C *** INITIALIZATION *** C 20 IV(NITER) = 0 IV(NFCALL) = 1 IV(NGCALL) = 1 IV(NFGCAL) = 1 IV(MODE) = -1 IV(STGLIM) = 2 IV(TOOBIG) = 0 IV(CNVCOD) = 0 IV(COVMAT) = 0 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(RADINC) = 0 IV(RESTOR) = 0 IV(FDH) = 0 V(RAD0) = ZERO V(STPPAR) = ZERO V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC)) C C *** SET INITIAL MODEL AND S MATRIX *** C IV(MODEL) = 1 IF (IV(S) .LT. 0) GO TO 999 IF (IV(INITS) .GT. 1) IV(MODEL) = 2 S1 = IV(S) IF (IV(INITS) .EQ. 0 .OR. IV(INITS) .GT. 2) 1 CALL DV7SCP(P*(P+1)/2, V(S1), ZERO) IV(1) = 1 J = IV(IPIVOT) IF (J .LE. 0) GO TO 999 DO 30 I = 1, P IV(J) = I J = J + 1 30 CONTINUE GO TO 999 C C *** NEW FUNCTION VALUE *** C 40 IF (IV(MODE) .EQ. 0) GO TO 290 IF (IV(MODE) .GT. 0) GO TO 520 C IV(1) = 2 IF (IV(TOOBIG) .EQ. 0) GO TO 999 IV(1) = 63 GO TO 999 C C *** NEW GRADIENT *** C 50 IV(KALM) = -1 IV(KAGQT) = -1 IV(FDH) = 0 IF (IV(MODE) .GT. 0) GO TO 520 C C *** MAKE SURE GRADIENT COULD BE COMPUTED *** C IF (IV(TOOBIG) .EQ. 0) GO TO 60 IV(1) = 65 GO TO 999 60 IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 610 C C *** COMPUTE D**-1 * GRADIENT *** C DIG1 = IV(DIG) K = DIG1 DO 70 I = 1, P V(K) = G(I) / D(I) K = K + 1 70 CONTINUE V(DGNORM) = DV2NRM(P, V(DIG1)) C IF (IV(CNVCOD) .NE. 0) GO TO 510 IF (IV(MODE) .EQ. 0) GO TO 440 IV(MODE) = 0 V(F0) = V(F) IF (IV(INITS) .LE. 2) GO TO 100 C C *** ARRANGE FOR FINITE-DIFFERENCE INITIAL S *** C IV(XIRC) = IV(COVREQ) IV(COVREQ) = -1 IF (IV(INITS) .GT. 3) IV(COVREQ) = 1 IV(CNVCOD) = 70 GO TO 530 C C *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S *** C 80 IV(CNVCOD) = 0 IV(MODE) = 0 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(COVREQ) = IV(XIRC) S1 = IV(S) PP1O2 = PS * (PS + 1) / 2 HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 90 CALL DV2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1)) GO TO 100 90 RMAT1 = IV(RMAT) CALL DL7SQR(PS, V(S1), V(RMAT1)) CALL DV2AXY(PP1O2, V(S1), NEGONE, V(S1), V(H1)) 100 IV(1) = 2 C C C----------------------------- MAIN LOOP ----------------------------- C C C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** C 110 CALL DITSUM(D, G, IV, LIV, LV, P, V, X) 120 K = IV(NITER) IF (K .LT. IV(MXITER)) GO TO 130 IV(1) = 10 GO TO 999 130 IV(NITER) = K + 1 C C *** UPDATE RADIUS *** C IF (K .EQ. 0) GO TO 150 STEP1 = IV(STEP) DO 140 I = 1, P V(STEP1) = D(I) * V(STEP1) STEP1 = STEP1 + 1 140 CONTINUE STEP1 = IV(STEP) T = V(RADFAC) * DV2NRM(P, V(STEP1)) IF (V(RADFAC) .LT. ONE .OR. T .GT. V(RADIUS)) V(RADIUS) = T C C *** INITIALIZE FOR START OF NEXT ITERATION *** C 150 X01 = IV(X0) V(F0) = V(F) IV(IRC) = 4 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(MODEL) C C *** COPY X TO X0 *** C CALL DV7CPY(P, V(X01), X) C C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** C 160 IF (.NOT. STOPX(DUMMY)) GO TO 180 IV(1) = 11 GO TO 190 C C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. C 170 IF (V(F) .GE. V(F0)) GO TO 180 V(RADFAC) = ONE K = IV(NITER) GO TO 130 C 180 IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 200 IV(1) = 9 190 IF (V(F) .GE. V(F0)) GO TO 999 C C *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH C *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. C IV(CNVCOD) = IV(1) GO TO 430 C C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . C 200 STEP1 = IV(STEP) W1 = IV(W) H1 = IV(H) T1 = ONE IF (IV(MODEL) .EQ. 2) GO TO 210 T1 = ZERO C C *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... C RMAT1 = IV(RMAT) IF (RMAT1 .LE. 0) GO TO 210 QTR1 = IV(QTR) IF (QTR1 .LE. 0) GO TO 210 IPIV1 = IV(IPIVOT) CALL DL7MST(D, G, IV(IERR), IV(IPIV1), IV(KALM), P, V(QTR1), 1 V(RMAT1), V(STEP1), V, V(W1)) C *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN, C *** SO WE MARK IT INVALID... IV(H) = -IABS(H1) C *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO C *** MARK INVALID THE INFORMATION DG7QTS MAY HAVE STORED IN V... IV(KAGQT) = -1 GO TO 260 C 210 IF (H1 .GT. 0) GO TO 250 C C *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** C H1 = -H1 IV(H) = H1 IV(FDH) = 0 J = IV(HC) IF (J .GT. 0) GO TO 220 J = H1 RMAT1 = IV(RMAT) CALL DL7SQR(P, V(H1), V(RMAT1)) 220 S1 = IV(S) DO 240 I = 1, P T = ONE / D(I) DO 230 K = 1, I V(H1) = T * (V(J) + T1*V(S1)) / D(K) J = J + 1 H1 = H1 + 1 S1 = S1 + 1 230 CONTINUE 240 CONTINUE H1 = IV(H) IV(KAGQT) = -1 C C *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** C 250 DIG1 = IV(DIG) LMAT1 = IV(LMAT) CALL DG7QTS(D, V(DIG1), V(H1), IV(KAGQT), V(LMAT1), P, V(STEP1), 1 V, V(W1)) IF (IV(KALM) .GT. 0) IV(KALM) = 0 C 260 IF (IV(IRC) .NE. 6) GO TO 270 IF (IV(RESTOR) .NE. 2) GO TO 290 RSTRST = 2 GO TO 300 C C *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** C 270 IV(TOOBIG) = 0 IF (V(DSTNRM) .LE. ZERO) GO TO 290 IF (IV(IRC) .NE. 5) GO TO 280 IF (V(RADFAC) .LE. ONE) GO TO 280 IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 280 IF (IV(RESTOR) .NE. 2) GO TO 290 RSTRST = 0 GO TO 300 C C *** COMPUTE F(X0 + STEP) *** C 280 X01 = IV(X0) STEP1 = IV(STEP) CALL DV2AXY(P, X, ONE, V(STEP1), V(X01)) IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 999 C C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . C 290 RSTRST = 3 300 X01 = IV(X0) V(RELDX) = DRLDST(P, D, X, V(X01)) CALL DA7SST(IV, LIV, LV, V) STEP1 = IV(STEP) LSTGST = IV(STLSTG) I = IV(RESTOR) + 1 GO TO (340, 310, 320, 330), I 310 CALL DV7CPY(P, X, V(X01)) GO TO 340 320 CALL DV7CPY(P, V(LSTGST), V(STEP1)) GO TO 340 330 CALL DV7CPY(P, V(STEP1), V(LSTGST)) CALL DV2AXY(P, X, ONE, V(STEP1), V(X01)) V(RELDX) = DRLDST(P, D, X, V(X01)) IV(RESTOR) = RSTRST C C *** IF NECESSARY, SWITCH MODELS *** C 340 IF (IV(SWITCH) .EQ. 0) GO TO 350 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL DV7CPY(NVSAVE, V, V(L)) 350 L = IV(IRC) - 4 STPMOD = IV(MODEL) IF (L .GT. 0) GO TO (370,380,390,390,390,390,390,390,500,440), L C C *** DECIDE WHETHER TO CHANGE MODELS *** C E = V(PREDUC) - V(FDIF) S1 = IV(S) CALL DS7LVM(PS, Y, V(S1), V(STEP1)) STTSST = HALF * DD7TPR(PS, V(STEP1), Y) IF (IV(MODEL) .EQ. 1) STTSST = -STTSST IF ( ABS(E + STTSST) * V(FUZZ) .GE. ABS(E)) GO TO 360 C C *** SWITCH MODELS *** C IV(MODEL) = 3 - IV(MODEL) IF (-2 .LT. L) GO TO 400 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL DV7CPY(NVSAVE, V(L), V) GO TO 160 C 360 IF (-3 .LT. L) GO TO 400 C C *** RECOMPUTE STEP WITH NEW RADIUS *** C 370 V(RADIUS) = V(RADFAC) * V(DSTNRM) GO TO 160 C C *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST C 380 V(RADIUS) = V(LMAXS) GO TO 200 C C *** CONVERGENCE OR FALSE CONVERGENCE *** C 390 IV(CNVCOD) = L IF (V(F) .GE. V(F0)) GO TO 510 IF (IV(XIRC) .EQ. 14) GO TO 510 IV(XIRC) = 14 C C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . C 400 IV(COVMAT) = 0 IV(REGD) = 0 C C *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** C IF (IV(IRC) .NE. 3) GO TO 430 STEP1 = IV(STEP) TEMP1 = IV(STLSTG) TEMP2 = IV(W) C C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** C HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 410 CALL DS7LVM(P, V(TEMP1), V(HC1), V(STEP1)) GO TO 420 410 RMAT1 = IV(RMAT) CALL DL7TVM(P, V(TEMP1), V(RMAT1), V(STEP1)) CALL DL7VML(P, V(TEMP1), V(RMAT1), V(TEMP1)) C 420 IF (STPMOD .EQ. 1) GO TO 430 S1 = IV(S) CALL DS7LVM(PS, V(TEMP2), V(S1), V(STEP1)) CALL DV2AXY(PS, V(TEMP1), ONE, V(TEMP2), V(TEMP1)) C C *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** C 430 IV(NGCALL) = IV(NGCALL) + 1 G01 = IV(W) CALL DV7CPY(P, V(G01), G) IV(1) = 2 IV(TOOBIG) = 0 GO TO 999 C C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** C 440 G01 = IV(W) CALL DV2AXY(P, V(G01), NEGONE, V(G01), G) STEP1 = IV(STEP) TEMP1 = IV(STLSTG) TEMP2 = IV(W) IF (IV(IRC) .NE. 3) GO TO 470 C C *** SET V(RADFAC) BY GRADIENT TESTS *** C C *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(X0) - G(X))) *** C K = TEMP1 L = G01 DO 450 I = 1, P V(K) = (V(K) - V(L)) / D(I) K = K + 1 L = L + 1 450 CONTINUE C C *** DO GRADIENT TESTS *** C IF (DV2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) GO TO 460 IF (DD7TPR(P, G, V(STEP1)) 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 470 460 V(RADFAC) = V(INCFAC) C C *** COMPUTE Y VECTOR NEEDED FOR UPDATING S *** C 470 CALL DV2AXY(PS, Y, NEGONE, Y, G) C C *** DETERMINE SIZING FACTOR V(SIZE) *** C C *** SET TEMP1 = S * STEP *** S1 = IV(S) CALL DS7LVM(PS, V(TEMP1), V(S1), V(STEP1)) C T1 = ABS(DD7TPR(PS, V(STEP1), V(TEMP1))) T = ABS(DD7TPR(PS, V(STEP1), Y)) V(SIZE) = ONE IF (T .LT. T1) V(SIZE) = T / T1 C C *** SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI *** C HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 480 CALL DS7LVM(PS, V(G01), V(HC1), V(STEP1)) GO TO 490 C 480 RMAT1 = IV(RMAT) CALL DL7TVM(PS, V(G01), V(RMAT1), V(STEP1)) CALL DL7VML(PS, V(G01), V(RMAT1), V(G01)) C 490 CALL DV2AXY(PS, V(G01), ONE, Y, V(G01)) C C *** UPDATE S *** C CALL DS7LUP(V(S1), V(COSMIN), PS, V(SIZE), V(STEP1), V(TEMP1), 1 V(TEMP2), V(G01), V(WSCALE), Y) IV(1) = 2 GO TO 110 C C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . C C *** BAD PARAMETERS TO ASSESS *** C 500 IV(1) = 64 GO TO 999 C C C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** C 510 IF (IV(RDREQ) .EQ. 0) GO TO 600 IF (IV(FDH) .NE. 0) GO TO 600 IF (IV(CNVCOD) .GE. 7) GO TO 600 IF (IV(REGD) .GT. 0) GO TO 600 IF (IV(COVMAT) .GT. 0) GO TO 600 IF (IABS(IV(COVREQ)) .GE. 3) GO TO 560 IF (IV(RESTOR) .EQ. 0) IV(RESTOR) = 2 GO TO 530 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE *** C 520 IV(RESTOR) = 0 530 CALL DF7HES(D, G, I, IV, LIV, LV, P, V, X) GO TO (540, 550, 580), I 540 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 999 C 550 IV(NGCOV) = IV(NGCOV) + 1 IV(NGCALL) = IV(NGCALL) + 1 IV(NFGCAL) = IV(NFCALL) + IV(NGCOV) IV(1) = 2 GO TO 999 C 560 H1 = IABS(IV(H)) IV(H) = -H1 PP1O2 = P * (P + 1) / 2 RMAT1 = IV(RMAT) IF (RMAT1 .LE. 0) GO TO 570 LMAT1 = IV(LMAT) CALL DV7CPY(PP1O2, V(LMAT1), V(RMAT1)) V(RCOND) = ZERO GO TO 590 570 HC1 = IV(HC) IV(FDH) = H1 CALL DV7CPY(P*(P+1)/2, V(H1), V(HC1)) C C *** COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN C *** FOR USE IN CALLER*S COVARIANCE CALCULATION... C 580 LMAT1 = IV(LMAT) H1 = IV(FDH) IF (H1 .LE. 0) GO TO 600 IF (IV(CNVCOD) .EQ. 70) GO TO 80 CALL DL7SRT(1, P, V(LMAT1), V(H1), I) IV(FDH) = -1 V(RCOND) = ZERO IF (I .NE. 0) GO TO 600 C 590 IV(FDH) = -1 STEP1 = IV(STEP) T = DL7SVN(P, V(LMAT1), V(STEP1), V(STEP1)) IF (T .LE. ZERO) GO TO 600 T = T / DL7SVX(P, V(LMAT1), V(STEP1), V(STEP1)) IF (T .GT. DR7MDC(4)) IV(FDH) = H1 V(RCOND) = T C 600 IV(MODE) = 0 IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 GO TO 999 C C *** SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH C *** IV(HC) .LE. 0 AND IV(RMAT) .LE. 0 C 610 IV(1) = 1400 C 999 RETURN C C *** LAST LINE OF DG7LIT FOLLOWS *** END SUBROUTINE DL7NVR(N, LIN, L) C C *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED *** C *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. *** C C *** PARAMETERS *** C INTEGER N DOUBLE PRECISION L(1), LIN(1) C DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2) C C *** LOCAL VARIABLES *** C INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1 DOUBLE PRECISION ONE, T, ZERO PARAMETER (ONE=1.D+0, ZERO=0.D+0) C C *** BODY *** C NP1 = N + 1 J0 = N*(NP1)/2 DO 30 II = 1, N I = NP1 - II LIN(J0) = ONE/L(J0) IF (I .LE. 1) GO TO 999 J1 = J0 IM1 = I - 1 DO 20 JJ = 1, IM1 T = ZERO J0 = J1 K0 = J1 - JJ DO 10 K = 1, JJ T = T - L(K0)*LIN(J0) J0 = J0 - 1 K0 = K0 + K - I 10 CONTINUE LIN(J0) = T/L(K0) 20 CONTINUE J0 = J0 - 1 30 CONTINUE 999 RETURN C *** LAST LINE OF DL7NVR FOLLOWS *** END SUBROUTINE DL7TSQ(N, A, L) C C *** SET A TO LOWER TRIANGLE OF (L**T) * L *** C C *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. *** C *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. *** C INTEGER N DOUBLE PRECISION A(1), L(1) C DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) C INTEGER I, II, IIM1, I1, J, K, M DOUBLE PRECISION LII, LJ C II = 0 DO 50 I = 1, N I1 = II + 1 II = II + I M = 1 IF (I .EQ. 1) GO TO 30 IIM1 = II - 1 DO 20 J = I1, IIM1 LJ = L(J) DO 10 K = I1, J A(M) = A(M) + LJ*L(K) M = M + 1 10 CONTINUE 20 CONTINUE 30 LII = L(II) DO 40 J = I1, II 40 A(J) = LII * L(J) 50 CONTINUE C 999 RETURN C *** LAST LINE OF DL7TSQ FOLLOWS *** END SUBROUTINE DN3RDP(IV, LIV, LV, N, P, RD, RHOI, RHOR, V) C C *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 *** C INTEGER LIV, LV, N, P INTEGER IV(LIV), RHOI(*) DOUBLE PRECISION RD(N), RHOR(*), V(LV) C C *** NOTE -- V IS PASSED FOR POSSIBLE USE BY REVISED VERSIONS OF C *** THIS ROUTINE. C INTEGER COV1, I, I1, I2, IEND, II, J, K, K1, NI, PU, PX, PX1, XNI C C *** IV AND V SUBSCRIPTS *** C INTEGER BS, BSSTR, COVMAT, COVPRT, COVREQ, LOO, NB, NEEDHD, NFCOV, 1 NFIX, NGCOV, PRUNIT, RDREQ, REGD, RCOND, STATPR, XNOTI C PARAMETER (BS=85, BSSTR=86, COVMAT=26, COVPRT=14, COVREQ=15, 1 LOO=84, NB=87, NEEDHD=36, NFCOV=52, NFIX=83, NGCOV=53, 2 PRUNIT=21, RDREQ=57, REGD=67, RCOND=53, STATPR=23, 3 XNOTI=90) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C PU = IV(PRUNIT) IF (PU .EQ. 0) GO TO 999 IF (IV(STATPR) .EQ. 0) GO TO 30 IF (IV(NFCOV) .GT. 0) WRITE(PU,10) IV(NFCOV) 10 FORMAT(/1X,I4,50H EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOST 1ICS.) IF (IV(NGCOV) .GT. 0) WRITE(PU,20) IV(NGCOV) 20 FORMAT(1X,I4,50H EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTI 1CS.) IF (IV(NFCOV) .GT. 0 .OR. IV(NGCOV) .GT. 0) IV(NEEDHD) = 1 C 30 IF (IV(COVPRT) .LE. 0) GO TO 999 COV1 = IV(COVMAT) IF (IV(REGD) .LE. 0 .AND. COV1 .LE. 0) GO TO 70 IV(NEEDHD) = 1 IF (IABS(IV(COVREQ)) .GT. 2) GO TO 50 C WRITE(PU,40) V(RCOND) 40 FORMAT(/53H SQRT(RECIPROCAL CONDITION OF F.D. HESSIAN) = AT MOST, 1 G10.2) GO TO 70 C 50 WRITE(PU,60) V(RCOND) 60 FORMAT(/54H SQRT(RECIPROCAL CONDITION OF (J**T)*RHO"*J) = AT MOST, 1 G10.2) C 70 IF (MOD(IV(COVPRT),2) .EQ. 0) GO TO 210 IV(NEEDHD) = 1 IF (COV1) 80,110,130 80 IF (-1 .EQ. COV1) WRITE(PU,90) 90 FORMAT(/43H ++++++ INDEFINITE COVARIANCE MATRIX ++++++) IF (-2 .EQ. COV1) WRITE(PU,100) 100 FORMAT(/52H ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++) GO TO 999 C 110 WRITE(PU,120 ) 120 FORMAT(/45H ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++) GO TO 210 C 130 IF (IABS(IV(COVREQ)) .LT. 3) GO TO 150 WRITE(PU,140) 140 FORMAT(/35H COVARIANCE = (J**T * RHO" * J)**-1/) GO TO 170 150 WRITE(PU,160) 160 FORMAT(/56H COVARIANCE = H**-1, WHERE H = FINITE-DIFFERENCE HESSIA 1N/) 170 II = COV1 - 1 DO 180 I = 1, P I1 = II + 1 I2 = II + MIN(I, 5) II = II + I WRITE(PU,190) I, (V(J), J = I1, I2) IF (I .LE. 5) GO TO 180 I2 = I2 + 1 WRITE(PU,200) (V(J), J = I2, II) 180 CONTINUE 190 FORMAT(4H ROW,I3,2X,5G12.3) 200 FORMAT(9X,5G12.3) 210 IF (IV(COVPRT) .LT. 2) GO TO 999 I = IV(REGD) + 4 GO TO (230, 250, 270, 290, 310), I WRITE(PU,220) IV(REGD) 220 FORMAT(/18H BUG... IV(REGD) =,I10) GO TO 999 230 WRITE(PU,240) NB, IV(NB) 240 FORMAT(/17H BAD IV(NB) = IV(,I2,3H) =,I10) GO TO 999 250 WRITE(PU,260) NFIX, IV(NFIX) 260 FORMAT(/19H BAD IV(NFIX) = IV(,I2,3H) =,I10) GO TO 999 270 WRITE(PU,280) LOO, IV(LOO) 280 FORMAT(/18H BAD IV(LOO) = IV(,I2,3H) =,I10) GO TO 999 290 WRITE(PU,300) 300 FORMAT(/42H REGRESSION DIAGNOSTIC VECTOR NOT COMPUTED) GO TO 999 310 IV(NEEDHD) = 1 XNI = 0 I = MOD(IV(RDREQ)/2, 3) + 1 GO TO (999, 330, 320), I 320 XNI = IV(XNOTI) PX = P - IV(NFIX) PX1 = PX - 1 IF (IV(LOO) .GT. 1) GO TO 400 330 WRITE(PU,340) 340 FORMAT (74H REGRESSION DIAGNOSTIC = 0.5 * G(I)**T * H(I)**-1 * H * 1 H(I)**-1 * G(I)...) IF (XNI .LE. 0) GO TO 380 WRITE(PU, 350) 350 FORMAT(29H I RD(I) X(I)) DO 360 I = 1, N WRITE(PU, 370) I, RD(I), (RHOR(J), J = XNI, XNI+PX1) XNI = XNI + PX 360 CONTINUE 370 FORMAT(1X,I5,G13.3,4G15.6/(19X,4G15.6)) GO TO 999 C 380 WRITE(PU,390) RD 390 FORMAT(6G12.3) GO TO 999 C 400 WRITE(PU,410) 410 FORMAT(/77H BLOCK REGRESSION DIAGNOSTIC = 0.5 * G(I)**T * H(I)**-1 1 * H * H(I)**-1 * G(I)) NI = IV(NB) K = IV(BS) K1 = IV(BSSTR) IEND = 0 IF (XNI .GT. 0) GO TO 450 WRITE(PU,420) 420 FORMAT(28H BLOCK FIRST LAST RD(I)) DO 440 I = 1, NI I1 = IEND + 1 IF (I1 .GT. N) GO TO 999 IEND = IEND + RHOI(K) K = K + K1 IF (IEND .GT. N) IEND = N WRITE(PU,430) I, I1, IEND, RD(I) 430 FORMAT(I6,I7,I6,G12.3) 440 CONTINUE GO TO 999 C 450 WRITE(PU, 460) 460 FORMAT(41H BLOCK FIRST LAST RD(I) X(I)) DO 480 I = 1, NI I1 = IEND + 1 IF (I1 .GT. N) GO TO 999 IEND = IEND + RHOI(K) K = K + K1 IF (IEND .GT. N) IEND = N WRITE(PU,470) I, I1, IEND, RD(I), (RHOR(J), J = XNI, XNI+PX1) 470 FORMAT(I6,I7,I6,G12.3,3G15.6/(31X,3G15.6)) XNI = XNI + PX 480 CONTINUE C 999 RETURN C *** LAST LINE OF DN3RDP FOLLOWS *** END //GO.SYSIN DD dglfg.f cat >dglfgb.f <<'//GO.SYSIN DD dglfgb.f' SUBROUTINE DGLGB(N, P, PS, X, B, RHO, RHOI, RHOR, IV, LIV, LV, 1 V, CALCRJ, UI, UR, UF) C C *** GENERALIZED LINEAR REGRESSION A LA NL2SOL, PLUS SIMPLE BOUNDS *** C C *** PARAMETERS *** C INTEGER N, P, PS, LIV, LV INTEGER IV(LIV), RHOI(*), UI(*) DOUBLE PRECISION B(2,P), X(P), RHOR(*), V(LV), UR(*) EXTERNAL CALCRJ, RHO, UF C C *** PARAMETER USAGE *** C C N....... TOTAL NUMBER OF RESIDUALS. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S). C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C B....... BOUNDS TO ENFORCE... B(1,I) .LE. X(I) .LE. B(2,I). C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS. C SEE DRGLG FOR DETAILS ABOUT RHO. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV (SEE DISCUSSION BELOW). C LV...... LENGTH OF V (SEE DISCUSSION BELOW). C V....... FLOATING-POINT VALUES ARRAY. C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR AND JACOBIAN MATRIX. C UI...... PASSED UNCHANGED TO CALCRJ. C UR...... PASSED UNCHANGED TO CALCRJ. C UF...... PASSED UNCHANGED TO CALCRJ. C C *** CALCRJ CALLING SEQUENCE... C C CALL CALCRJ(N, P, X, NF, NEED, R, RP, UI, UR, UF) C C PARAMETERS N, P, X, UI, UR, AND UF ARE AS ABOVE. C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(P,N). C NEED IS AN INTEGER ARRAY OF LENGTH 2... C NEED(1) = 1 MEANS CALCRJ SHOULD COMPUTE THE RESIDUAL VECTOR R, C AND NEED(2) IS THE VALUE NF HAD AT THE LAST X WHERE C CALCRJ MIGHT BE CALLED WITH NEED(1) = 2. C NEED(1) = 2 MEANS CALCRJ SHOULD COMPUTE THE JACOBIAN MATRIX RP, C WHERE RP(J,I) = DERIVATIVE OF R(I) WITH RESPECT TO X(J). C (CALCRJ SHOULD NOT CHANGE NEED AND SHOULD CHANGE AT MOST ONE OF R C AND RP. IF R OR RP, AS APPROPRIATE, CANNOT BE COMPUTED, THEN CALCRJ C SHOULD SET NF TO 0. OTHERWISE IT SHOULD NOT CHANGE NF.) C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DIVSET, DRGLGB C C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DRGLGB... CARRIES OUT OPTIMIZATION ITERATIONS. C C C *** LOCAL VARIABLES *** C INTEGER D1, DR1, I, IV1, NEED1(2), NEED2(2), NF, R1, RD1 C C *** IV COMPONENTS *** C INTEGER D, J, NEXTV, NFCALL, NFGCAL, R, REGD, REGD0, TOOBIG, VNEED PARAMETER (D=27, J=70, NEXTV=47, NFCALL=6, NFGCAL=7, R=61, 1 REGD=67, REGD0=82, TOOBIG=2, VNEED=4) SAVE NEED1, NEED2 DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 I = (P-PS+2)*(P-PS+1)/2 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+1+I) CALL DRGLGB(B, X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, 1 RHO, RHOI,RHOR, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + (P - PS + 1)*N IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N IV(NEXTV) = IV(J) + N*PS IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RD1 = IV(REGD0) C 20 CALL DRGLGB(B, V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, 1 V(R1), V(RD1), RHO, RHOI, RHOR, V, X) IF (IV(1)-2) 30, 50, 60 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) NEED1(2) = IV(NFGCAL) CALL CALCRJ(N, PS, X, NF, NEED1, V(R1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE DR = GRADIENT OF R COMPONENTS *** C 50 CALL CALCRJ(N, PS, X, IV(NFGCAL), NEED2, V(R1), V(DR1), UI, UR,UF) IF (IV(NFGCAL) .EQ. 0) IV(TOOBIG) = 1 GO TO 20 C C *** INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED C *** AND PRINT IT IF SO REQUESTED... C 60 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 C 999 RETURN C C *** LAST LINE OF DGLGB FOLLOWS *** END SUBROUTINE DGLFB(N, P, PS, X, B, RHO, RHOI, RHOR, IV, LIV, LV, V, 1 CALCRJ, UI, UR, UF) C C *** GENERALIZED LINEAR REGRESSION, FINITE-DIFFERENCE JACOBIAN *** C *** WITH SIMPLE BOUNDS ON X *** C C *** PARAMETERS *** C INTEGER N, P, PS, LIV, LV INTEGER IV(LIV), RHOI(*), UI(*) DOUBLE PRECISION B(2,P), X(P), V(LV), RHOR(*), UR(*) EXTERNAL CALCRJ, RHO, UF C C *** PARAMETER USAGE *** C C N....... TOTAL NUMBER OF RESIDUALS. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S). C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C B....... BOUNDS TO ENFORCE... B(1,I) .LE. X(I) .LE. B(2,I). C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS. C SEE DRGLG FOR DETAILS ABOUT RHO. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV (SEE DISCUSSION BELOW). C LV...... LENGTH OF V (SEE DISCUSSION BELOW). C V....... FLOATING-POINT VALUES ARRAY. C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR. C UI...... PASSED UNCHANGED TO CALCRJ. C UR...... PASSED UNCHANGED TO CALCRJ. C UF...... PASSED UNCHANGED TO CALCRJ. C C *** CALCRJ CALLING SEQUENCE... C C CALL CALCRJ(N, P, X, NF, NEED, R, RP, UI, UR, UF) C C PARAMETERS N, P, X, UI, UR, AND UF ARE AS ABOVE. C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(P,N). C NEED MAY BE REGARDED AS AN INTEGER THAT ALWAYS HAS THE VALUE 1 C WHEN DGLFB CALLS CALCRJ. THIS MEANS CALCRJ SHOULD COMPUTE THE C RESIDUAL VECTOR R. (CALCRJ SHOULD NOT CHANGE NEED OR RP. IF R C CANNOT BE COMPUTED, THEN CALCRJ SHOULD SET NF TO 0. OTHERWISE IT C SHOULD NOT CHANGE NF. FOR COMPATIBILITY WITH DGLG, NEED IS A C VECTOR OF LENGTH 2.) C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DIVSET, DRGLGB,DV7CPY C C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DRGLGB... CARRIES OUT OPTIMIZATION ITERATIONS. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C C *** LOCAL VARIABLES *** C INTEGER D1, DK, DR1, I, I1, IV1, J1K, J1K0, K, NEED(2), NF, 1 NG, RD1, R1, R21, RS1, RSN DOUBLE PRECISION H, H0, HLIM, NEGPT5, T, ONE, XK, XK1, ZERO C C *** IV AND V COMPONENTS *** C INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL, 1 NGCALL, NGCOV, R, REGD0, TOOBIG, VNEED PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35, 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53, 2 R=61, REGD0=82, TOOBIG=2, VNEED=4) SAVE NEED DATA HLIM/0.1D+0/, NEGPT5/-0.5D+0/, ONE/1.D+0/, ZERO/0.D+0/ DATA NEED(1)/1/, NEED(2)/0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) IV(COVREQ) = -IABS(IV(COVREQ)) IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 I = (P-PS+2)*(P-PS+1)/2 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+3+I) CALL DRGLGB(B, X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, 1 RHOI, RHOR, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + (P - PS + 3)*N IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N IV(NEXTV) = IV(J) + N*PS IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RD1 = IV(REGD0) R21 = RD1 - N RS1 = R21 - N RSN = RS1 + N - 1 C 20 CALL DRGLGB(B, V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, 1 V(R1), V(RD1), RHO, RHOI, RHOR, V, X) IF (IV(1)-2) 30, 50, 999 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) CALL CALCRJ(N, PS, X, NF, NEED, V(R1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 CALL DV7CPY(N, V(RS1), V(R1)) IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R *** C C *** INITIALIZE D IF NECESSARY *** C 50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO) 1 CALL DV7SCP(P, V(D1), ONE) C DK = D1 NG = IV(NGCALL) - 1 IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1 J1K0 = DR1 NF = IV(NFCALL) IF (NF .EQ. IV(NFGCAL)) GO TO 70 NG = NG + 1 CALL CALCRJ(N, PS, X, NF, NEED, V(RS1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 70 60 IV(TOOBIG) = 1 IV(NGCALL) = NG GO TO 20 70 DO 130 K = 1, PS J1K = J1K0 J1K0 = J1K0 + 1 IF (B(1,K) .GE. B(2,K)) GO TO 120 XK = X(K) H = V(DLTFDJ) * MAX( ABS(XK), ONE/V(DK)) H0 = H DK = DK + 1 T = NEGPT5 XK1 = XK + H IF (XK - H .GE. B(1,K)) GO TO 80 T = -T IF (XK1 .GT. B(2,K)) GO TO 60 80 IF (XK1 .LE. B(2,K)) GO TO 90 T = -T H = -H XK1 = XK + H IF (XK1 .LT. B(1,K)) GO TO 60 90 X(K) = XK1 NF = IV(NFGCAL) CALL CALCRJ(N, PS, X, NF, NEED, V(R21), V(DR1), UI, UR, UF) NG = NG + 1 IF (NF .GT. 0) GO TO 100 H = T * H XK1 = XK + H IF ( ABS(H/H0) .GE. HLIM) GO TO 90 GO TO 60 100 X(K) = XK IV(NGCALL) = NG I1 = R21 DO 110 I = RS1, RSN V(J1K) = (V(I1) - V(I)) / H I1 = I1 + 1 J1K = J1K + PS 110 CONTINUE GO TO 130 C *** SUPPLY A ZERO DERIVATIVE FOR CONSTANT COMPONENTS... 120 DO 125 I = 1, N V(J1K) = ZERO J1K = J1K + PS 125 CONTINUE 130 CONTINUE GO TO 20 C 999 RETURN C C *** LAST LINE OF DGLFB FOLLOWS *** END SUBROUTINE DRGLGB(B, D, DR, IV, LIV, LV, N, ND, NN, P, PS, R, 1 RD, RHO, RHOI, RHOR, V, X) C C *** ITERATION DRIVER FOR GENERALIZED (NON)LINEAR MODELS (ETC.) C INTEGER LIV, LV, N, ND, NN, P, PS INTEGER IV(LIV), RHOI(*) DOUBLE PRECISION B(2,P), D(P), DR(ND,N), R(*), RD(*), RHOR(*), 1 V(LV), X(*) C DIMENSION RD(N, (P-PS)*(P-PS+1)/2 + 1) EXTERNAL RHO C C-------------------------- PARAMETER USAGE -------------------------- C C B........ BOUNDS ON X. C D........ SCALE VECTOR. C DR....... DERIVATIVES OF R AT X. C IV....... INTEGER VALUES ARRAY. C LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 82. C LV....... LENGTH OF V... LV MUST BE AT LEAST 105 + P*(2*P+16). C N........ TOTAL NUMBER OF RESIDUALS. C ND....... LEADING DIMENSION OF DR -- MUST BE AT LEAST PS. C NN....... LEAD DIMENSION OF R, RD. C P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS....... NUMBER OF NON-NUISANCE PARAMETERS. C R........ RESIDUALS (OR MEANS -- FUNCTIONS OF X) WHEN DRGLGB IS CALLED C WITH IV(1) = 1. C RD....... TEMPORARY STORAGE. C RHO...... COMPUTES INFO ABOUT OBJECTIVE FUNCTION. C RHOI..... PASSED WITHOUT CHANGE TO RHO. C RHOR..... PASSED WITHOUT CHANGE TO RHO. C V........ FLOATING-POINT VALUES ARRAY. C X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C C *** CALLING SEQUENCE FOR RHO... C C CALL RHO(NEED, F, N, NF, XN, R, RD, RHOI, RHOR, W) C C PARAMETER DECLARATIONS FOR RHO... C C INTEGER NEED(2), N, NF, RHOI(*) C FLOATING-POINT F, XN(*), R(*), RD(N,*), RHOR(*), W(N) C C RHOI AND RHOR ARE FOR RHO TO USE AS IT SEES FIT. THEY ARE PASSED C TO RHO WITHOUT CHANGE. C F, R, RD, AND W ARE EXPLAINED BELOW WITH NEED. C XN IS THE VECTOR OF NUISANCE PARAMETERS (OF LENGTH P - PS). IF C RHO NEEDS TO KNOW THE LENGTH OF XN, THEN THIS LENGTH SHOULD BE C COMMUNICATED THROUGH RHOI (OR THROUGH COMMON). RHO SHOULD NOT CHANGE C XN. C NEED(1) = 1 MEANS RHO SHOULD SET F TO THE SUM OF THE LOSS FUNCTION C VALUES AT THE RESIDUALS R(I). NF IS THE CURRENT FUNCTION INVOCATION C COUNT (A VALUE THAT IS INCREMENTED EACH TIME A NEW PARAMETER EXTIMATE C X IS CONSIDERED). NEED(2) IS THE VALUE NF HAD AT THE LAST R WHERE C RHO MIGHT BE CALLED WITH NEED(1) = 2. IF RHO SAVES INTERMEDIATE C RESULTS FOR USE IN CALLS WITH NEED(1) = 2, THEN IT CAN USE NF TO TELL C WHICH INTERMEDIATE RESULTS ARE APPROPRIATE, AND IT CAN SAVE SOME OF C THESE RESULTS IN R. C NEED(1) = 2 MEANS RHO SHOULD SET R(I) TO THE LOSS FUNCTION C DERIVATIVE WITH RESPECT TO THE RESIDUALS THAT WERE PASSED TO RHO WHEN C NF HAD THE SAME VALUE IT DOES NOW (AND NEED(1) WAS 1). RHO SHOULD C ALSO SET W(I) TO THE APPROXIMATION OF THE SECOND DERIVATIVE OF THE C LOSS FUNCTION (WITH RESPECT TO THE I-TH RESIDUAL) THAT SHOULD BE USED C IN THE GAUSS-NEWTON MODEL. WHEN THERE ARE NUISANCE PARAMETERS (I.E., C WHEN PS .LT. P) RHO SHOULD ALSO SET R(I+K*N) TO THE DERIVATIVE OF THE C LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL AND XN(K), AND IT C SHOULD SET RD(I,J + K*(K+1)/2 + 1) TO THE SECOND PARTIAL DERIVATIVE C OF THE I-TH RESIDUAL WITH RESPECT TO XN(J) AND XN(K), 0 .LE. J .LE. K C AND 1 .LE. K .LE. P - PS, WHERE XN(0) MEANS THE I-TH RESIDUAL ITSELF. C IN ANY EVENT, RHO SHOULD ALSO SET RD(I,1) TO THE (TRUE) SECOND C DERIVATIVE OF THE LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL. C NF (THE FUNCTION INVOCATION COUNT WHOSE NORMAL USE IS EXPLAINED C ABOVE) SHOULD NOT BE CHANGED UNLESS RHO CANNOT CARRY OUT THE REQUESTED C TASK, IN WHICH CASE RHO SHOULD SET NF TO 0. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL DIVSET, DD7TPR, DD7UP5, DG7ITB,DITSUM, DL7ITV, DL7IVM, 1 DL7SRT, DL7SQR, DL7SVX, DL7SVN,DL7VML,DO7PRD, 2 DQ7ADR,DV2AXY,DV7CPY, DV7SCL, DV7SCP, DVSUM DOUBLE PRECISION DD7TPR, DL7SVX, DL7SVN, DVSUM C C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C DD7UP5... UPDATES SCALE VECTOR D. C DG7ITB... PERFORMS BASIC MINIMIZATION ALGORITHM. C DITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. C DL7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C DL7IVM... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C DL7SQR... COMPUTES L*(L**T) FOR LOWER TRIANG. MATRIX L. C DL7SVX... UNDERESTIMATES LARGEST SINGULAR VALUE OF TRIANG. MATRIX. C DL7SVN... OVERESTIMATES SMALLEST SINGULAR VALUE OF TRIANG. MATRIX. C DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DO7PRD.... ADDS OUTER PRODUCT OF VECTORS TO A MATRIX. C DQ7ADR... ADDS ROWS TO QR FACTORIZATION. C DV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C DV7SCL... MULTIPLIES A VECTOR BY A SCALAR. C DVSUM.... RETURNS SUM OF ELEMENTS OF A VECTOR. C C *** LOCAL VARIABLES *** C LOGICAL UPDATD, ZEROG INTEGER G1, HN1, I, II, IV1, J, J1, JTOL1, K, LH, 1 NEED1(2), NEED2(2), PMPS, PS1, PSLEN, QTR1, 2 RMAT1, STEP1, TEMP1, TEMP2, TEMP3, TEMP4, W, WI, Y1 DOUBLE PRECISION RHMAX, RHTOL, RHO1, RHO2, T C DOUBLE PRECISION ONE, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER DINIT, DTYPE, DTINIT, D0INIT, F, 1 F0, G, HC, IPIVOT, IVNEED, JCN, JTOL, LMAT, 2 MODE, NEXTIV, NEXTV, NF0, NF1, NFCALL, NFGCAL, 3 QTR, RDREQ, REGD, RESTOR, RMAT, 4 RSPTOL, STEP, TOOBIG, VNEED C C *** IV SUBSCRIPT VALUES *** C PARAMETER (DTYPE=16, F0=13, G=28, HC=71, IPIVOT=76, IVNEED=3, 1 JCN=66, JTOL=59, LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, 2 NFCALL=6, NF0=68, NF1=69, NFGCAL=7, QTR=77, RESTOR=9, 3 RMAT=78, RDREQ=57, REGD=67, STEP=40, TOOBIG=2, VNEED=4) C C *** V SUBSCRIPT VALUES *** C PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RSPTOL=49) PARAMETER (ONE=1.D+0, ZERO=0.D+0) SAVE NEED1, NEED2 DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C LH = P * (P+1) / 2 IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) PS1 = PS + 1 IV1 = IV(1) IF (IV1 .GT. 2) GO TO 10 W = IV(G) - N IV(RESTOR) = 0 IF (IV(TOOBIG) .EQ. 0) GO TO (110, 120), IV1 V(F) = V(F0) IF (IV1 .NE. 1) IV(1) = 2 GO TO 40 C C *** FRESH START OR RESTART -- CHECK INPUT INTEGERS *** C 10 IF (ND .LT. PS) GO TO 340 IF (PS .GT. P) GO TO 340 IF (PS .LE. 0) GO TO 340 IF (N .LE. 0) GO TO 340 IF (IV1 .EQ. 14) GO TO 30 IF (IV1 .GT. 16) GO TO 360 IF (IV1 .LT. 12) GO TO 40 IF (IV1 .EQ. 12) IV(1) = 13 IF (IV(1) .NE. 13) GO TO 20 IV(IVNEED) = IV(IVNEED) + P IV(VNEED) = IV(VNEED) + P*(P+13)/2 + 2*N + 4*PS 20 CALL DG7ITB(B, D, X, IV, LIV, LV, P, PS, V, X, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(IPIVOT) = IV(NEXTIV) IV(NEXTIV) = IV(IPIVOT) + P IV(G) = IV(NEXTV) + P + N IV(RMAT) = IV(G) + P + 4*PS IV(QTR) = IV(RMAT) + LH IV(JTOL) = IV(QTR) + P + N IV(JCN) = IV(JTOL) + 2*P IV(NEXTV) = IV(JCN) + P C *** TURN OFF COVARIANCE COMPUTATION *** IV(RDREQ) = 0 IF (IV1 .EQ. 13) GO TO 999 C 30 JTOL1 = IV(JTOL) IF (V(DINIT) .GE. ZERO) CALL DV7SCP(P, D, V(DINIT)) IF (V(DTINIT) .GT. ZERO) CALL DV7SCP(P, V(JTOL1), V(DTINIT)) I = JTOL1 + P IF (V(D0INIT) .GT. ZERO) CALL DV7SCP(P, V(I), V(D0INIT)) IV(NF0) = 0 IV(NF1) = 0 C 40 G1 = IV(G) Y1 = G1 - (P + N) CALL DG7ITB(B, D, V(G1), IV, LIV, LV, P, PS, V, X, V(Y1)) IF (IV(1) - 2) 50, 60, 350 C 50 V(F) = ZERO IF (IV(NF1) .EQ. 0) GO TO 999 IF (IV(RESTOR) .NE. 2) GO TO 999 IV(NF0) = IV(NF1) CALL DV7CPY(N, RD, R) IV(REGD) = 0 GO TO 999 C 60 CALL DV7SCP(P, V(G1), ZERO) RMAT1 = IABS(IV(RMAT)) QTR1 = IABS(IV(QTR)) CALL DV7SCP(PS, V(QTR1), ZERO) IV(REGD) = 0 CALL DV7SCP(PS, V(Y1), ZERO) CALL DV7SCP(LH, V(RMAT1), ZERO) IF (IV(RESTOR) .NE. 3) GO TO 70 CALL DV7CPY(N, R, RD) IV(NF1) = IV(NF0) 70 CALL RHO(NEED2, T, N, IV(NFGCAL), X(PS1), R, RD, RHOI, RHOR, V(W)) IF (IV(NFGCAL) .GT. 0) GO TO 90 80 IV(TOOBIG) = 1 GO TO 40 90 IF (IV(MODE) .LT. 0) GO TO 999 DO 100 I = 1, N 100 CALL DV2AXY(PS, V(Y1), R(I), DR(1,I), V(Y1)) GO TO 999 C C *** COMPUTE F(X) *** C 110 I = IV(NFCALL) NEED1(2) = IV(NFGCAL) CALL RHO(NEED1, V(F), N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IV(NF1) = I IF (I .LE. 0) GO TO 80 GO TO 40 C 120 G1 = IV(G) C C *** DECIDE WHETHER TO UPDATE D BELOW *** C I = IV(DTYPE) UPDATD = .FALSE. IF (I .LE. 0) GO TO 130 IF (I .EQ. 1 .OR. IV(MODE) .LT. 0) UPDATD = .TRUE. C C *** COMPUTE RMAT AND QTR *** C 130 QTR1 = IABS(IV(QTR)) RMAT1 = IABS(IV(RMAT)) IV(RMAT) = RMAT1 IV(HC) = 0 IV(NF0) = 0 IV(NF1) = 0 IF (IV(MODE) .LT. 0) GO TO 150 C C *** ADJUST Y *** C Y1 = IV(G) - (P + N) WI = W STEP1 = IV(STEP) DO 140 I = 1, N T = V(WI) - RD(I) WI = WI + 1 IF (T .NE. ZERO) CALL DV2AXY(PS, V(Y1), 1 T*DD7TPR(PS,V(STEP1),DR(1,I)), DR(1,I), V(Y1)) 140 CONTINUE C C *** CHECK FOR NEGATIVE W COMPONENTS *** C 150 J1 = W + N - 1 DO 160 WI = W, J1 IF (V(WI) .LT. ZERO) GO TO 230 160 CONTINUE C C *** W IS NONNEGATIVE. COMPUTE QR FACTORIZATION *** C *** AND, IF NECESSARY, USE SEMINORMAL EQUATIONS *** C RHMAX = ZERO RHTOL = V(RSPTOL) TEMP1 = G1 + P ZEROG = .TRUE. WI = W DO 190 I = 1, N RHO1 = R(I) RHO2 = V(WI) WI = WI + 1 T = SQRT(RHO2) IF (RHMAX .LT. RHO2) RHMAX = RHO2 IF (RHO2 .GT. RHTOL*RHMAX) GO TO 170 C *** SEMINORMAL EQUATIONS *** CALL DV2AXY(PS, V(G1), RHO1, DR(1,I), V(G1)) RHO1 = ZERO ZEROG = .FALSE. GO TO 180 170 RHO1 = RHO1 / T C *** QR ACCUMULATION *** 180 CALL DV7SCL(PS, V(TEMP1), T, DR(1,I)) CALL DQ7ADR(PS, V(QTR1), V(RMAT1), V(TEMP1), RHO1) 190 CONTINUE C C *** COMPUTE G FROM RMAT AND QTR *** C TEMP2 = TEMP1 + P CALL DL7VML(PS, V(TEMP1), V(RMAT1), V(QTR1)) IF (ZEROG) GO TO 210 IV(QTR) = -QTR1 IF (DL7SVX(PS, V(RMAT1), V(TEMP2), V(TEMP2)) * RHTOL .GE. 1 DL7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2))) GO TO 220 CALL DL7IVM(PS, V(TEMP2), V(RMAT1), V(G1)) C C *** SEMINORMAL EQUATIONS CORRECTION OF BJOERCK -- C *** ONE CYCLE OF ITERATIVE REFINEMENT... C TEMP3 = TEMP2 + PS TEMP4 = TEMP3 + PS CALL DL7ITV(PS, V(TEMP3), V(RMAT1), V(TEMP2)) CALL DV7SCP(PS, V(TEMP4), ZERO) RHMAX = ZERO WI = W DO 200 I = 1, N RHO2 = V(WI) WI = WI + 1 IF (RHMAX .LT. RHO2) RHMAX = RHO2 RHO1 = ZERO IF (RHO2 .LE. RHTOL*RHMAX) RHO1 = R(I) T = RHO1 - RHO2*DD7TPR(PS, V(TEMP3), DR(1,I)) CALL DV2AXY(PS, V(TEMP4), T, DR(1,I), V(TEMP4)) 200 CONTINUE CALL DL7IVM(PS, V(TEMP3), V(RMAT1), V(TEMP4)) CALL DV2AXY(PS, V(TEMP2), ONE, V(TEMP3), V(TEMP2)) CALL DV2AXY(PS, V(QTR1), ONE, V(TEMP2), V(QTR1)) 210 IV(QTR) = QTR1 220 CALL DV2AXY(PS, V(G1), ONE, V(TEMP1), V(G1)) IF (PS .GE. P) GO TO 330 GO TO 250 C C *** INDEFINITE GN HESSIAN... *** C 230 IV(RMAT) = -RMAT1 IV(HC) = RMAT1 CALL DO7PRD(N, LH, PS, V(RMAT1), V(W), DR, DR) C C *** COMPUTE GRADIENT *** C G1 = IV(G) DO 240 I = 1, N 240 CALL DV2AXY(PS, V(G1), R(I), DR(1,I), V(G1)) IF (PS .GE. P) GO TO 330 C C *** COMPUTE GRADIENT COMPONENTS OF NUISANCE PARAMETERS *** C 250 K = P - PS J1 = 1 G1 = G1 + PS DO 260 J = 1, K J1 = J1 + NN V(G1) = DVSUM(N, R(J1)) G1 = G1 + 1 260 CONTINUE C C *** COMPUTE HESSIAN COMPONENTS OF NUISANCE PARAMETERS *** C I = PS*PS1/2 PSLEN = P*(P+1)/2 - I HN1 = RMAT1 + I CALL DV7SCP(PSLEN, V(HN1), ZERO) PMPS = P - PS K = HN1 J1 = 1 DO 290 II = 1, PMPS J1 = J1 + NN J = J1 DO 270 I = 1, N CALL DV2AXY(PS, V(K), RD(J), DR(1,I), V(K)) J = J + 1 270 CONTINUE K = K + PS DO 280 I = 1, II J1 = J1 + NN V(K) = DVSUM(N, RD(J1)) K = K + 1 280 CONTINUE 290 CONTINUE IF (IV(RMAT) .LE. 0) GO TO 330 J = IV(LMAT) CALL DV7CPY(PSLEN, V(J), V(HN1)) IF (DL7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2)) .LE. ZERO) GO TO 300 CALL DL7SRT(PS1, P, V(RMAT1), V(RMAT1), I) IF (I .LE. 0) GO TO 310 C C *** HESSIAN IS NOT POSITIVE DEFINITE *** C 300 CALL DL7SQR(PS, V(RMAT1), V(RMAT1)) CALL DV7CPY(PSLEN, V(HN1), V(J)) IV(HC) = RMAT1 IV(RMAT) = -RMAT1 GO TO 330 C C *** NUISANCE PARS LEAVE HESSIAN POS. DEF. GET REST OF QTR *** C 310 J = QTR1 + PS G1 = IV(G) + PS DO 320 I = PS1, P T = DD7TPR(I-1, V(HN1), V(QTR1)) HN1 = HN1 + I V(J) = (V(G1) - T) / V(HN1-1) J = J + 1 G1 = G1 + 1 320 CONTINUE 330 IF (UPDATD) CALL DD7UP5(D, IV, LIV, LV, P, PS, V) GO TO 40 C C *** MISC. DETAILS *** C C *** BAD N, ND, OR P *** C 340 IV(1) = 66 GO TO 360 C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 350 G1 = IV(G) 360 CALL DITSUM(D, V(G1), IV, LIV, LV, P, V, X) C 999 RETURN C *** LAST LINE OF DRGLGB FOLLOWS *** END SUBROUTINE DD7MLP(N, X, Y, Z, K) C C *** SET X = DIAG(Y)**K * Z C *** FOR X, Z = LOWER TRIANG. MATRICES STORED COMPACTLY BY ROW C *** K = 1 OR -1. C INTEGER N, K DOUBLE PRECISION X(*), Y(N), Z(*) INTEGER I, J, L DOUBLE PRECISION ONE, T DATA ONE/1.D+0/ C L = 1 IF (K .GE. 0) GO TO 30 DO 20 I = 1, N T = ONE / Y(I) DO 10 J = 1, I X(L) = T * Z(L) L = L + 1 10 CONTINUE 20 CONTINUE GO TO 999 C 30 DO 50 I = 1, N T = Y(I) DO 40 J = 1, I X(L) = T * Z(L) L = L + 1 40 CONTINUE 50 CONTINUE 999 RETURN C *** LAST LINE OF DD7MLP FOLLOWS *** END SUBROUTINE DF7DHB(B, D, G, IRT, IV, LIV, LV, P, V, X) C C *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING C *** AT V(IV(FDH)) = V(-IV(H)). HONOR SIMPLE BOUNDS IN B. C C *** IF IV(COVREQ) .GE. 0 THEN DF7DHB USES GRADIENT DIFFERENCES, C *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN DG7LIT. C C IRT VALUES... C 1 = COMPUTE FUNCTION VALUE, I.E., V(F). C 2 = COMPUTE G. C 3 = DONE. C C C *** PARAMETER DECLARATIONS *** C INTEGER IRT, LIV, LV, P INTEGER IV(LIV) DOUBLE PRECISION B(2,P), D(P), G(P), V(LV), X(P) C C *** LOCAL VARIABLES *** C LOGICAL OFFSID INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2, 1 NEWM1, PP1O2, STPI, STPM, STP0 DOUBLE PRECISION DEL, DEL0, T, XM, XM1 DOUBLE PRECISION HALF, HLIM, ONE, TWO, ZERO C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DV7CPY, DV7SCP C C DV7CPY.... COPY ONE VECTOR TO ANOTHER. C DV7SCP... COPY SCALAR TO ALL COMPONENTS OF A VECTOR. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H, KAGQT, MODE, 1 NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE C PARAMETER (HALF=0.5D+0, HLIM=0.1D+0, ONE=1.D+0, TWO=2.D+0, 1 ZERO=0.D+0) C PARAMETER (COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42, F=10, 1 FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7, 2 SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C IRT = 4 KIND = IV(COVREQ) M = IV(MODE) IF (M .GT. 0) GO TO 10 HES = IABS(IV(H)) IV(H) = -HES IV(FDH) = 0 IV(KAGQT) = -1 V(FX) = V(F) C *** SUPPLY ZEROS IN CASE B(1,I) = B(2,I) FOR SOME I *** CALL DV7SCP(P*(P+1)/2, V(HES), ZERO) 10 IF (M .GT. P) GO TO 999 IF (KIND .LT. 0) GO TO 120 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND C *** GRADIENT VALUES. C GSAVE1 = IV(W) + P IF (M .GT. 0) GO TO 20 C *** FIRST CALL ON DF7DHB. SET GSAVE = G, TAKE FIRST STEP *** CALL DV7CPY(P, V(GSAVE1), G) IV(SWITCH) = IV(NFGCAL) GO TO 80 C 20 DEL = V(DELTA) X(M) = V(XMSAVE) IF (IV(TOOBIG) .EQ. 0) GO TO 30 C C *** HANDLE OVERSIZE V(DELTA) *** C DEL0 = V(DELTA0) * MAX(ONE/D(M), ABS(X(M))) DEL = HALF * DEL IF ( ABS(DEL/DEL0) .LE. HLIM) GO TO 140 C 30 HES = -IV(H) C C *** SET G = (G - GSAVE)/DEL *** C DEL = ONE / DEL DO 40 I = 1, P G(I) = DEL * (G(I) - V(GSAVE1)) GSAVE1 = GSAVE1 + 1 40 CONTINUE C C *** ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** C K = HES + M*(M-1)/2 L = K + M - 2 IF (M .EQ. 1) GO TO 60 C C *** SET H(I,M) = 0.5 * (H(I,M) + G(I)) FOR I = 1 TO M-1 *** C MM1 = M - 1 DO 50 I = 1, MM1 IF (B(1,I) .LT. B(2,I)) V(K) = HALF * (V(K) + G(I)) K = K + 1 50 CONTINUE C C *** ADD H(I,M) = G(I) FOR I = M TO P *** C 60 L = L + 1 DO 70 I = M, P IF (B(1,I) .LT. B(2,I)) V(L) = G(I) L = L + I 70 CONTINUE C 80 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 340 IF (B(1,M) .GE. B(2,M)) GO TO 80 C C *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE *** C DEL = V(DELTA0) * MAX(ONE/D(M), ABS(X(M))) XM = X(M) IF (XM .LT. ZERO) GO TO 90 XM1 = XM + DEL IF (XM1 .LE. B(2,M)) GO TO 110 XM1 = XM - DEL IF (XM1 .GE. B(1,M)) GO TO 100 GO TO 280 90 XM1 = XM - DEL IF (XM1 .GE. B(1,M)) GO TO 100 XM1 = XM + DEL IF (XM1 .LE. B(2,M)) GO TO 110 GO TO 280 C 100 DEL = -DEL 110 V(XMSAVE) = XM X(M) = XM1 V(DELTA) = DEL IRT = 2 GO TO 999 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. C 120 STP0 = IV(W) + P - 1 MM1 = M - 1 MM1O2 = M*MM1/2 HES = -IV(H) IF (M .GT. 0) GO TO 130 C *** FIRST CALL ON DF7DHB. *** IV(SAVEI) = 0 GO TO 240 C 130 IF (IV(TOOBIG) .EQ. 0) GO TO 150 C *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** 140 IV(FDH) = -2 GO TO 350 150 I = IV(SAVEI) IF (I .GT. 0) GO TO 190 C C *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** C PP1O2 = P * (P-1) / 2 HPM = HES + PP1O2 + MM1 V(HPM) = V(F) C C *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** C NEWM1 = 1 GO TO 260 160 HMI = HES + MM1O2 IF (MM1 .EQ. 0) GO TO 180 HPI = HES + PP1O2 DO 170 I = 1, MM1 T = ZERO IF (B(1,I) .LT. B(2,I)) T = V(FX) - (V(F) + V(HPI)) V(HMI) = T HMI = HMI + 1 HPI = HPI + 1 170 CONTINUE 180 V(HMI) = V(F) - TWO*V(FX) IF (OFFSID) V(HMI) = V(FX) - TWO*V(F) C C *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** C I = 0 GO TO 200 C 190 X(I) = V(DELTA) C C *** FINISH COMPUTING H(M,I) *** C STPI = STP0 + I HMI = HES + MM1O2 + I - 1 STPM = STP0 + M V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM)) 200 I = I + 1 IF (I .GT. M) GO TO 230 IF (B(1,I) .LT. B(2,I)) GO TO 210 GO TO 200 C 210 IV(SAVEI) = I STPI = STP0 + I V(DELTA) = X(I) X(I) = X(I) + V(STPI) IRT = 1 IF (I .LT. M) GO TO 999 NEWM1 = 2 GO TO 260 220 X(M) = V(XMSAVE) - DEL IF (OFFSID) X(M) = V(XMSAVE) + TWO*DEL GO TO 999 C 230 IV(SAVEI) = 0 X(M) = V(XMSAVE) C 240 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 330 IF (B(1,M) .LT. B(2,M)) GO TO 250 GO TO 240 C C *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. C *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN C *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. C 250 V(XMSAVE) = X(M) NEWM1 = 3 260 XM = V(XMSAVE) DEL = V(DLTFDC) * MAX(ONE/D(M), ABS(XM)) XM1 = XM + DEL OFFSID = .FALSE. IF (XM1 .LE. B(2,M)) GO TO 270 OFFSID = .TRUE. XM1 = XM - DEL IF (XM - TWO*DEL .GE. B(1,M)) GO TO 300 GO TO 280 270 IF (XM-DEL .GE. B(1,M)) GO TO 290 OFFSID = .TRUE. IF (XM + TWO*DEL .LE. B(2,M)) GO TO 310 C 280 IV(FDH) = -2 GO TO 350 C 290 IF (XM .GE. ZERO) GO TO 310 XM1 = XM - DEL 300 DEL = -DEL 310 GO TO (160, 220, 320), NEWM1 320 X(M) = XM1 STPM = STP0 + M V(STPM) = DEL IRT = 1 GO TO 999 C C *** HANDLE SPECIAL CASE OF B(1,P) = B(2,P) -- CLEAR SCRATCH VALUES C *** FROM LAST ROW OF FDH... C 330 IF (B(1,P) .LT. B(2,P)) GO TO 340 I = HES + P*(P-1)/2 CALL DV7SCP(P, V(I), ZERO) C C *** RESTORE V(F), ETC. *** C 340 IV(FDH) = HES 350 V(F) = V(FX) IRT = 3 IF (KIND .LT. 0) GO TO 999 IV(NFGCAL) = IV(SWITCH) GSAVE1 = IV(W) + P CALL DV7CPY(P, G, V(GSAVE1)) GO TO 999 C 999 RETURN C *** LAST LINE OF DF7DHB FOLLOWS *** END SUBROUTINE DG7ITB(B, D, G, IV, LIV, LV, P, PS, V, X, Y) C C *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR *** C *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) *** C *** HAVING SIMPLE BOUNDS ON THE PARAMETERS BEING ESTIMATED. *** C C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, P, PS INTEGER IV(LIV) DOUBLE PRECISION B(2,P), D(P), G(P), V(LV), X(P), Y(P) C C-------------------------- PARAMETER USAGE -------------------------- C C B.... VECTOR OF LOWER AND UPPER BOUNDS ON X. C D.... SCALE VECTOR. C IV... INTEGER VALUE ARRAY. C LIV.. LENGTH OF IV. MUST BE AT LEAST 80. C LH... LENGTH OF H = P*(P+1)/2. C LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7. C G.... GRADIENT AT X (WHEN IV(1) = 2). C HC... GAUSS-NEWTON HESSIAN AT X (WHEN IV(1) = 2). C P.... NUMBER OF PARAMETERS (COMPONENTS IN X). C PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. C V.... FLOATING-POINT VALUE ARRAY. C X.... PARAMETER VECTOR. C Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). C C *** DISCUSSION *** C C DG7ITB IS SIMILAR TO DG7LIT, EXCEPT FOR THE EXTRA PARAMETER B C -- DG7ITB ENFORCES THE BOUNDS B(1,I) .LE. X(I) .LE. B(2,I), C I = 1(1)P. C DG7ITB PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF C REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES C IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED C FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES C THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED C COMPACTLY BY ROWS), AND DG7ITB BUILDS AN APPROXIMATION, S, TO THE C SECOND-ORDER TERM. THE CALLER ALSO PROVIDES THE FUNCTION VALUE, C GRADIENT, AND PART OF THE YIELD VECTOR USED IN UPDATING S. C DG7ITB DECIDES DYNAMICALLY WHETHER OR NOT TO USE S WHEN CHOOSING C THE NEXT STEP TO TRY... THE HESSIAN APPROXIMATION USED IS EITHER C HC ALONE (GAUSS-NEWTON MODEL) OR HC + S (AUGMENTED MODEL). C IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT C CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO C 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS C IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN C COMPUTED HAS NONZERO VALUES IN THESE ROWS. C C IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY C FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS C USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME C WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2, C 1, OR 2). C C FOR UPDATING S, DG7ITB ASSUMES THAT THE GRADIENT HAS THE FORM C OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE C GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS C THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP, C THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF C RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY C PART OF THIS IN Y, NAMELY THE SUM OVER I OF C RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING DG7ITB WITH IV(1) = 2 AND C IV(MODE) = 0 (WHERE MODE = 38). G THEN CONTANS THE OTHER PART, C SO THAT THE DESIRED YIELD VECTOR IS G - Y. IF PS .LT. P, THEN C THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF C GRAD(R(I,X)), STEP, AND Y. C C PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING C ONES TO DN2GB (AND NL2SOL), EXCEPT THAT V CAN BE SHORTER C (SINCE THE PART OF V THAT DN2GB USES FOR STORING D, J, AND R IS C NOT NEEDED). MOREOVER, COMPARED WITH DN2GB (AND NL2SOL), IV(1) C MAY HAVE THE TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE C EXPLAINED BELOW, AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). C THE VALUES IV(D), IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM C DN2GB (AND DN2FB), ARE NOT REFERENCED BY DG7ITB O