MODULE nlsol ! Non-linear least-squares fitting using the updated versions of ! TOMS 573 (NL2SOL) routines from the Bell Labs. public domain PORT directory. ! nl2sol -- an adaptive nonlinear least-squares algorithm ! ! authors = John E. Dennis, jr., David M. Gay, and Roy E. Welsch ! ! ACM Transactions on Mathematical Software, September, 1981. ! Translation from Fortran 66/77 format to compatibility with ELF90 ! by Alan Miller ! amiller@bigpond.net.au ! http://www.ozemail.com.au/~milleraj ! Latest revision - 31 May 2002 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) ! *** V SUBSCRIPT VALUES *** INTEGER, SAVE :: afctol=31, bias=43, cosmin=47, decfac=22, delta=52, & delta0=44, dfac=41, dgnorm=1, dinit=38, dltfdc=42, & dltfdj=43, dstnrm=2, dst0=3, dstsav=18, dtinit=39, & d0init=40, epslon=19, eta0=42, f=10, fdif=11, & flstgd=12, f0=13, fuzz=45, gtslst=14, gtstep=4, & huberc=48, incfac=23, lmax0=35, lmaxs=36, nreduc=6, & phmnfc=20, phmxfc=21, plstgd=15, preduc=7, radfac=16, & radius=8, rad0=9, rdfcmn=24, rdfcmx=25, reldx=17, & rfctol=32, rlimit=46, rsptol=49, sctol=37, sigmin=50, & stppar=5, tuner1=26, tuner2=27, tuner3=28, tuner4=29, & tuner5=30, xctol=33, xftol=34 ! *** IV SUBSCRIPT VALUES *** INTEGER, SAVE :: algsav=51, cnvcod=55, covmat=26, covprt=14, covreq=15, & dradpr=101, dtype=16, dtype0=54, fdh=74, fx=53, h=56, & hc=71, ierr=75, inith=25, inits=25, ipivot=76, irc=29, & ivneed=3, jcn=66, jtol=59, kagqt=33, lastiv=44, & lastv=45, lmat=42, mlstgd=32, mode=35, model=5, & mxfcal=17, mxiter=18, needhd=36, nextiv=46, nextv=47, & nfcall=6, nfcov=52, nfgcal=7, ngcall=30, ngcov=53, & niter=31, nvdflt=50, nvsave=9, oldn=38, outlev=19, & parprt=20, parsav=49, perm=58, prntit=39, prunit=21, & qrtyp=80, radinc=8, rcond=53, rdreq=57, regd=67, & restor=9, rmat=78, s=62, savei=63, solprt=22, stage=10, & statpr=23, step=40, stglim=11, sused=64, switch=12, & toobig=2, vneed=4, vsave=60, w0=65, xirc=13, xmsave=51, & x0prt=24 REAL (dp), PARAMETER :: one=1.0_dp, zero=0.0_dp REAL (dp), ALLOCATABLE :: scalef(:) CONTAINS FUNCTION i7mdcn(k) RESULT(ival) INTEGER, INTENT(IN) :: k INTEGER :: ival ! *** RETURN INTEGER MACHINE-DEPENDENT CONSTANTS *** ! *** K = 1 MEANS RETURN STANDARD OUTPUT UNIT NUMBER. *** ! *** K = 2 MEANS RETURN ALTERNATE OUTPUT UNIT NUMBER. *** ! *** K = 3 MEANS RETURN INPUT UNIT NUMBER. *** ! (NOTE -- K = 2, 3 ARE USED ONLY BY TEST PROGRAMS.) INTEGER, PARAMETER :: mdcon(3) = (/ 6, 0, 5 /) ival = mdcon(k) RETURN ! *** LAST CARD OF I7MDCN FOLLOWS *** END FUNCTION i7mdcn FUNCTION r7mdc(k) RESULT(fn_val) ! *** RETURN MACHINE DEPENDENT CONSTANTS USED BY NL2SOL *** INTEGER, INTENT(IN) :: k REAL (dp) :: fn_val ! *** the constant returned depends upon k... ! *** k = 1... smallest pos. eta such that -eta exists. ! *** k = 2... square root of eta. ! *** k = 3... unit roundoff = smallest pos. no. machep such ! *** that 1 + machep > 1 .and. 1 - machep < 1. ! *** k = 4... square root of machep. ! *** k = 5... square root of big (see k = 6). ! *** k = 6... largest machine no. big such that -big exists. SELECT CASE (k) CASE(1) fn_val = TINY(1._dp) CASE(2) fn_val = SQRT( TINY(1._dp) ) CASE(3) fn_val = EPSILON(1._dp) CASE(4) fn_val = SQRT( EPSILON(1._dp) ) CASE(5) fn_val = SQRT( HUGE(1._dp) ) CASE(6) fn_val = HUGE(1._dp) CASE DEFAULT WRITE(*, *) '** Illegal argument, k = ', k, ' to FUNCTION rmdcon **' END SELECT RETURN ! *** LAST CARD OF R7MDC FOLLOWS *** END FUNCTION r7mdc SUBROUTINE v7dfl(alg, v) ! *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V *** ! *** ALG = 1 MEANS REGRESSION CONSTANTS. ! *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. INTEGER, INTENT(IN) :: alg REAL (dp), INTENT(OUT) :: v(:) ! REAL :: r7mdc ! EXTERNAL r7mdc ! R7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS REAL (dp) :: machep, mepcrt, sqteps ! *** SUBSCRIPTS FOR V *** REAL (dp), PARAMETER :: three=3.0_dp !------------------------------- BODY -------------------------------- machep = r7mdc(3) v(afctol) = 1.e-20_dp IF (machep > 1.e-10_dp) v(afctol) = machep**2 v(decfac) = 0.5_dp sqteps = r7mdc(4) v(dfac) = 0.6_dp v(dtinit) = 1.e-6_dp mepcrt = machep ** (one/three) v(d0init) = 1.0_dp v(epslon) = 0.1_dp v(incfac) = 2.0_dp v(lmax0) = 1.0_dp v(lmaxs) = 1.0_dp v(phmnfc) = -0.1_dp v(phmxfc) = 0.1_dp v(rdfcmn) = 0.1_dp v(rdfcmx) = 4.0_dp v(rfctol) = MAX(1.e-10_dp, mepcrt**2) v(sctol) = v(rfctol) v(tuner1) = 0.1_dp v(tuner2) = 1.e-4_dp v(tuner3) = 0.75_dp v(tuner4) = 0.5_dp v(tuner5) = 0.75_dp v(xctol) = sqteps v(xftol) = 100.0_dp * machep IF (alg >= 2) GO TO 10 ! *** REGRESSION VALUES v(cosmin) = MAX(1.e-6_dp, 100._dp * machep) v(dinit) = 0.0_dp v(delta0) = sqteps v(dltfdc) = mepcrt v(dltfdj) = sqteps v(fuzz) = 1.5_dp v(huberc) = 0.7_dp v(rlimit) = r7mdc(5) v(rsptol) = 1.e-3_dp v(sigmin) = 1.e-4_dp GO TO 999 ! *** GENERAL OPTIMIZATION VALUES 10 v(bias) = 0.8_dp v(dinit) = -1.0_dp v(eta0) = 1000._dp * machep 999 RETURN ! *** LAST CARD OF V7DFL FOLLOWS *** END SUBROUTINE v7dfl SUBROUTINE ivset(alg, iv, liv, lv, v) ! *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V *** ! *** ALG = 1 MEANS REGRESSION CONSTANTS. ! *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. INTEGER, INTENT(IN) :: alg INTEGER, INTENT(OUT) :: iv(:) INTEGER, INTENT(IN) :: liv INTEGER, INTENT(IN) :: lv REAL (dp), INTENT(OUT) :: v(:) ! INTEGER :: i7mdcn ! EXTERNAL i7mdcn, v7dfl ! I7MDCN... RETURNS MACHINE-DEPENDENT INTEGER CONSTANTS. ! V7DFL.... PROVIDES DEFAULT VALUES TO V. INTEGER :: alg1, miv, mv INTEGER, PARAMETER :: miniv(4) = (/ 82, 59, 103, 103 /), & minv(4) = (/ 98, 71, 101, 85 /) !------------------------------- BODY -------------------------------- IF (prunit <= liv) iv(prunit) = i7mdcn(1) IF (algsav <= liv) iv(algsav) = alg IF (alg < 1 .OR. alg > 4) GO TO 40 miv = miniv(alg) IF (liv < miv) GO TO 20 mv = minv(alg) IF (lv < mv) GO TO 30 alg1 = MOD(alg-1,2) + 1 CALL v7dfl(alg1, v) iv(1) = 12 IF (alg > 2) iv(dradpr) = 1 iv(ivneed) = 0 iv(lastiv) = miv iv(lastv) = mv iv(lmat) = mv + 1 iv(mxfcal) = 200 iv(mxiter) = 150 iv(outlev) = 1 iv(parprt) = 1 iv(perm) = miv + 1 iv(solprt) = 1 iv(statpr) = 1 iv(vneed) = 0 iv(x0prt) = 1 IF (alg1 >= 2) GO TO 10 ! *** REGRESSION VALUES iv(covprt) = 3 iv(covreq) = 1 iv(dtype) = 1 iv(hc) = 0 iv(ierr) = 0 iv(inits) = 0 iv(ipivot) = 0 iv(nvdflt) = 32 iv(vsave) = 58 IF (alg > 2) iv(vsave) = iv(vsave) + 3 iv(parsav) = iv(vsave) + nvsave iv(qrtyp) = 1 iv(rdreq) = 3 iv(rmat) = 0 GO TO 999 ! *** GENERAL OPTIMIZATION VALUES 10 iv(dtype) = 0 iv(inith) = 1 iv(nfcov) = 0 iv(ngcov) = 0 iv(nvdflt) = 25 iv(parsav) = 47 IF (alg > 2) iv(parsav) = 61 GO TO 999 20 iv(1) = 15 GO TO 999 30 iv(1) = 16 GO TO 999 40 iv(1) = 67 999 RETURN ! *** LAST CARD OF IVSET FOLLOWS *** END SUBROUTINE ivset SUBROUTINE a7sst(iv, v) ! *** ASSESS CANDIDATE STEP (***SOL VERSION 2.3) *** INTEGER, INTENT(IN OUT) :: iv(:) REAL (dp), INTENT(IN OUT) :: v(:) ! *** PURPOSE *** ! THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION ROUTINE ! TO ASSESS THE NEXT CANDIDATE STEP. IT MAY RECOMMEND ONE OF SEVERAL ! COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOMPUTING IT USING THE ! SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE TO CONVERGENCE OR FALSE ! CONVERGENCE. SEE THE RETURN CODE LISTING BELOW. !-------------------------- PARAMETER USAGE -------------------------- ! IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION ! BELOW OF IV VALUES REFERENCED. ! LIV (IN) LENGTH OF IV ARRAY. ! LV (IN) LENGTH OF V ARRAY. ! V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION ! BELOW OF V VALUES REFERENCED. ! *** IV VALUES REFERENCED *** ! IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION, ! IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS ! SET WHEN STEP IS DEFINITELY TO BE ACCEPTED). ON INPUT ! AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE ! UNCHANGED SINCE THE PREVIOUS RETURN OF A7SST. ! ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE ! FOLLOWING VALUES... ! 1 = SWITCH MODELS OR TRY SMALLER STEP. ! 2 = SWITCH MODELS OR ACCEPT STEP. ! 3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT TESTS. ! 4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED. ! 5 = RECOMPUTE STEP (USING THE SAME MODEL). ! 6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT ! EVALUATE THE OBJECTIVE FUNCTION. ! 7 = X-CONVERGENCE (SEE V(XCTOL)). ! 8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)). ! 9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE. ! 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)). ! 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)). ! 12 = FALSE CONVERGENCE (SEE V(XFTOL)). ! 13 = IV(IRC) WAS OUT OF RANGE ON INPUT. ! RETURN CODE I HAS PRECEDENCE OVER I+1 FOR I = 9, 10, 11. ! IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL). ! IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING ! THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION. ! IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION, ! THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT. ! IV(NFCALL) (IN) INVOCATION COUNT FOR THE OBJECTIVE FUNCTION. ! IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST ! FUNCTION REDUCTION THIS ITERATION. IV(NFGCAL) REMAINS ! UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED. ! IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER ! OF DECREASES) SO FAR THIS ITERATION. ! IV(RESTOR) (OUT) SET TO 1 IF V(F) HAS BEEN RESTORED AND X SHOULD BE RESTORED ! TO ITS INITIAL VALUE, TO 2 IF X SHOULD BE SAVED, ! TO 3 IF X SHOULD BE RESTORED FROM THE SAVED VALUE, AND TO ! 0 OTHERWISE. ! IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE ! CURRENT ITERATION. ! IV(STGLIM) (IN) MAXIMUM NUMBER OF MODELS TO CONSIDER. ! IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT ! GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL, ! IN WHICH CASE A7SST SETS IV(SWITCH) = 1. ! IV(TOOBIG) (I/O) IS NONZERO ON INPUT IF STEP WAS TOO BIG (E.G., IF ! IT WOULD CAUSE OVERFLOW). IT IS SET TO 0 ON RETURN. ! IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF ! CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS. ! *** V VALUES REFERENCED *** ! V(AFCTOL) (IN) ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. IF THE ! ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS ! THAN V(AFCTOL) AND A7SST DOES NOT RETURN WITH ! IV(IRC) = 11, THEN A7SST RETURNS WITH IV(IRC) = 10. ! V(DECFAC) (IN) FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS ! NONZERO. ! V(DSTNRM) (IN) THE 2-NORM OF D*STEP. ! V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP. ! V(DST0) (IN) THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED, ! I.E., FOR V(NREDUC) >= 0). ! V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC- ! TION VALUE AT X. IF X IS RESTORED TO A PREVIOUS VALUE, ! THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE. ! V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT ! VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION ! DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE). ! V(FLSTGD) (I/O) SAVED VALUE OF V(F). ! V(F0) (IN) OBJECTIVE FUNCTION VALUE AT START OF ITERATION. ! V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP. ! V(GTSTEP) (IN) INNER PRODUCT BETWEEN STEP AND GRADIENT. ! V(INCFAC) (IN) MINIMUM FACTOR BY WHICH TO INCREASE RADIUS. ! V(LMAXS) (IN) MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND). ! IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE ! WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, OR 9 ! DOES NOT OCCUR, IF V(DSTNRM) > V(LMAXS) OR THE CURRENT ! STEP IS A NEWTON STEP, AND IF ! V(PREDUC) <= V(SCTOL) * ABS(V(F0)), THEN A7SST RETURNS ! WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE, THEN ! A7SST REPEATS THIS TEST (DISALLOWING A FULL NEWTON STEP) ! WITH V(PREDUC) COMPUTED FOR A STEP OF LENGTH V(LMAXS) ! (BY A RETURN WITH IV(IRC) = 6). ! V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR ! NEWTON STEP. IF A7SST IS CALLED WITH IV(IRC) = 6, I.E., ! IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR ! USE IN THE SINGULAR CONVERGENCE TEST, THEN V(NREDUC) IS ! SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED. ! V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP. ! V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR ! CURRENT STEP. ! V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS, ! WHICH SHOULD BE V(RADFAC)*DST, WHERE DST IS EITHER THE ! OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF ! DIAG(NEWD)*STEP FOR THE OUTPUT VALUE OF STEP AND THE ! UPDATED VERSION, NEWD, OF THE SCALE VECTOR D. FOR ! IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED. ! V(RDFCMN) (IN) MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT ! VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1. ! V(RDFCMX) (IN) MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0. ! V(RELDX) (IN) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED ! (E.G.) BY FUNCTION RLDST AS ! MAX (D(I)*ABS(X(I)-X0(I)), 1 <= I <= P) / ! MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 <= I <= P). ! V(RFCTOL) (IN) RELATIVE FUNCTION CONVERGENCE TOLERANCE. IF THE ! ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE- ! DICTED AND V(NREDUC) <= V(RFCTOL)*ABS(V(F0)), THEN ! A7SST RETURNS WITH IV(IRC) = 8 OR 9. ! V(SCTOL) (IN) SINGULAR CONVERGENCE TOLERANCE -- SEE V(LMAXS). ! V(STPPAR) (IN) MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP. ! V(TUNER1) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION ! REDUCTION WAS MUCH LESS THAN EXPECTED. SUGGESTED ! VALUE = 0.1. ! V(TUNER2) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION ! REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP. SUGGESTED ! VALUE = 10**-4. ! V(TUNER3) (IN) TUNING CONSTANT USED TO DECIDE IF THE RADIUS ! SHOULD BE INCREASED. SUGGESTED VALUE = 0.75. ! V(XCTOL) (IN) X-CONVERGENCE CRITERION. IF STEP IS A NEWTON STEP ! (V(STPPAR) = 0) HAVING V(RELDX) <= V(XCTOL) AND GIVING ! AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN ! A7SST RETURNS IV(IRC) = 7 OR 9. ! V(XFTOL) (IN) FALSE CONVERGENCE TOLERANCE. IF STEP GAVE NO OR ONLY ! A SMALL FUNCTION DECREASE AND V(RELDX) <= V(XFTOL), ! THEN A7SST RETURNS WITH IV(IRC) = 12. !------------------------------- NOTES ------------------------------- ! *** APPLICATION AND USAGE RESTRICTIONS *** ! THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR LEAST-SQUARES) ! PACKAGE. IT MAY BE USED IN ANY UNCONSTRAINED MINIMIZATION SOLVER THAT ! USES DOGLEG, GOLDFELD-QUANDT-TROTTER, OR LEVENBERG-MARQUARDT STEPS. ! *** ALGORITHM NOTES *** ! SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL ! SWITCHING STRATEGIES. WHILE NL2SOL CONSIDERS ONLY TWO MODELS, ! A7SST IS DESIGNED TO HANDLE ANY NUMBER OF MODELS. ! *** USAGE NOTES *** ! ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES ! STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND ! V(PREDUC) NEED HAVE BEEN INITIALIZED. BETWEEN CALLS, NO I/O ! VALUES EXCEPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER- ! ANCES SHOULD BE CHANGED. ! AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN ! CHANGE THE STOPPING TOLERANCES AND CALL A7SST AGAIN, IN WHICH ! CASE THE STOPPING TESTS WILL BE REPEATED. ! *** REFERENCES *** ! (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981), ! AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, ! ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3. ! (2) POWELL, M.J.D. (1970) A FORTRAN SUBROUTINE FOR SOLVING ! SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL ! METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY ! P. RABINOWITZ, GORDON AND BREACH, LONDON. ! *** HISTORY *** ! JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH ! IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY. ! DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE ! PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS ! PRESENT FORM (FALL 1978), WITH MINOR CHANGES TO THE SINGULAR ! CONVERGENCE TEST IN MAY, 1984 (TO DEAL WITH FULL NEWTON STEPS). ! *** GENERAL *** ! THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED BY ! THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324, DCR75-10143, ! 76-14311DSS, MCS76-11989, AND MCS-7906671. !------------------------ EXTERNAL QUANTITIES ------------------------ ! *** NO EXTERNAL FUNCTIONS AND SUBROUTINES *** !-------------------------- LOCAL VARIABLES -------------------------- LOGICAL :: goodx INTEGER :: i, nfc REAL (dp) :: emax, emaxs, gts, rfac1, xmax ! *** SUBSCRIPTS FOR IV AND V *** ! *** DATA INITIALIZATIONS *** REAL (dp), PARAMETER :: half=0.5_dp, onep2=1.2_dp, two=2.0_dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ nfc = iv(nfcall) iv(switch) = 0 iv(restor) = 0 rfac1 = one goodx = .true. i = iv(irc) SELECT CASE (i) CASE (1) GO TO 20 CASE (2) GO TO 30 CASE (3, 4) GO TO 10 CASE (5) GO TO 40 CASE (6) GO TO 280 CASE (7:11) GO TO 220 CASE (12) GO TO 170 CASE DEFAULT iv(irc) = 13 GO TO 999 END SELECT ! *** INITIALIZE FOR NEW ITERATION *** 10 iv(stage) = 1 iv(radinc) = 0 v(flstgd) = v(f0) IF (iv(toobig) == 0) GO TO 110 iv(stage) = -1 iv(xirc) = i GO TO 60 ! *** STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS *** ! *** FIRST DECIDE WHICH *** 20 IF (iv(model) /= iv(mlstgd)) GO TO 30 ! *** OLD MODEL RETAINED, SMALLER RADIUS TRIED *** ! *** DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION *** iv(stage) = iv(stglim) iv(radinc) = -1 GO TO 110 ! *** A NEW MODEL IS BEING TRIED. DECIDE WHETHER TO KEEP IT. *** 30 iv(stage) = iv(stage) + 1 ! *** NOW WE ADD THE POSSIBILITY THAT STEP WAS RECOMPUTED WITH *** ! *** THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP. *** 40 IF (iv(stage) > 0) GO TO 50 ! *** STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG. *** IF (iv(toobig) /= 0) GO TO 60 ! *** RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF. *** iv(stage) = -iv(stage) i = iv(xirc) SELECT CASE (i) CASE (1) GO TO 20 CASE (2) GO TO 30 CASE (3:4) GO TO 110 CASE (5) GO TO 70 END SELECT 50 IF (iv(toobig) == 0) GO TO 70 ! *** HANDLE OVERSIZE STEP *** iv(toobig) = 0 IF (iv(radinc) > 0) GO TO 80 iv(stage) = -iv(stage) iv(xirc) = iv(irc) 60 iv(toobig) = 0 v(radfac) = v(decfac) iv(radinc) = iv(radinc) - 1 iv(irc) = 5 iv(restor) = 1 v(f) = v(flstgd) GO TO 999 70 IF (v(f) < v(flstgd)) GO TO 110 ! *** THE NEW STEP IS A LOSER. RESTORE OLD MODEL. *** IF (iv(model) == iv(mlstgd)) GO TO 80 iv(model) = iv(mlstgd) iv(switch) = 1 ! *** RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F). 80 IF (v(flstgd) >= v(f0)) GO TO 110 iv(restor) = 1 v(f) = v(flstgd) v(preduc) = v(plstgd) v(gtstep) = v(gtslst) IF (iv(switch) == 0) rfac1 = v(dstnrm) / v(dstsav) v(dstnrm) = v(dstsav) nfc = iv(nfgcal) goodx = .false. 110 v(fdif) = v(f0) - v(f) IF (v(fdif) > v(tuner2) * v(preduc)) GO TO 140 IF (iv(radinc) > 0) GO TO 140 ! *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE ! *** -- SO TRY NEW MODEL OR SMALLER RADIUS IF (v(f) < v(f0)) GO TO 120 iv(mlstgd) = iv(model) v(flstgd) = v(f) v(f) = v(f0) iv(restor) = 1 GO TO 130 120 iv(nfgcal) = nfc 130 iv(irc) = 1 IF (iv(stage) < iv(stglim)) GO TO 160 iv(irc) = 5 iv(radinc) = iv(radinc) - 1 GO TO 160 ! *** NON-TRIVIAL FUNCTION DECREASE ACHIEVED *** 140 iv(nfgcal) = nfc rfac1 = one v(dstsav) = v(dstnrm) IF (v(fdif) > v(preduc)*v(tuner1)) GO TO 190 ! *** DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS ! *** OR ACCEPT STEP WITH DECREASED RADIUS. IF (iv(stage) >= iv(stglim)) GO TO 150 ! *** CONSIDER SWITCHING MODELS *** iv(irc) = 2 GO TO 160 ! *** ACCEPT STEP WITH DECREASED RADIUS *** 150 iv(irc) = 4 ! *** SET V(RADFAC) TO FLETCHER'S DECREASE FACTOR *** 160 iv(xirc) = iv(irc) emax = v(gtstep) + v(fdif) v(radfac) = half * rfac1 IF (emax < v(gtstep)) v(radfac) = rfac1 * MAX(v(rdfcmn), half * v(gtstep)/emax) ! *** DO FALSE CONVERGENCE TEST *** 170 IF (v(reldx) <= v(xftol)) GO TO 180 iv(irc) = iv(xirc) IF (v(f) < v(f0)) GO TO 200 GO TO 230 180 iv(irc) = 12 GO TO 240 ! *** HANDLE GOOD FUNCTION DECREASE *** 190 IF (v(fdif) < (-v(tuner3) * v(gtstep))) GO TO 210 ! *** INCREASING RADIUS LOOKS WORTHWHILE. SEE IF WE JUST ! *** RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP ! *** AFTER RECOMPUTING IT WITH A LARGER RADIUS. IF (iv(radinc) < 0) GO TO 210 IF (iv(restor) == 1) GO TO 210 ! *** WE DID NOT. TRY A LONGER STEP UNLESS THIS WAS A NEWTON STEP. v(radfac) = v(rdfcmx) gts = v(gtstep) IF (v(fdif) < (half/v(radfac) - one) * gts) & v(radfac) = MAX(v(incfac), half*gts/(gts + v(fdif))) iv(irc) = 4 IF (v(stppar) == zero) GO TO 230 IF (v(dst0) >= zero .AND. (v(dst0) < two*v(dstnrm) & .OR. v(nreduc) < onep2*v(fdif))) GO TO 230 ! *** STEP WAS NOT A NEWTON STEP. RECOMPUTE IT WITH ! *** A LARGER RADIUS. iv(irc) = 5 iv(radinc) = iv(radinc) + 1 ! *** SAVE VALUES CORRESPONDING TO GOOD STEP *** 200 v(flstgd) = v(f) iv(mlstgd) = iv(model) IF (iv(restor) /= 1) iv(restor) = 2 v(dstsav) = v(dstnrm) iv(nfgcal) = nfc v(plstgd) = v(preduc) v(gtslst) = v(gtstep) GO TO 230 ! *** ACCEPT STEP WITH RADIUS UNCHANGED *** 210 v(radfac) = one iv(irc) = 3 GO TO 230 ! *** COME HERE FOR A RESTART AFTER CONVERGENCE *** 220 iv(irc) = iv(xirc) IF (v(dstsav) >= zero) GO TO 240 iv(irc) = 12 GO TO 240 ! *** PERFORM CONVERGENCE TESTS *** 230 iv(xirc) = iv(irc) 240 IF (iv(restor) == 1 .AND. v(flstgd) < v(f0)) iv(restor) = 3 IF ( ABS(v(f)) < v(afctol)) iv(irc) = 10 IF (half * v(fdif) > v(preduc)) GO TO 999 emax = v(rfctol) * ABS(v(f0)) emaxs = v(sctol) * ABS(v(f0)) IF (v(preduc) <= emaxs .AND. (v(dstnrm) > v(lmaxs) .OR. & v(stppar) == zero)) iv(irc) = 11 IF (v(dst0) < zero) GO TO 250 i = 0 IF ((v(nreduc) > zero .AND. v(nreduc) <= emax) .OR. & (v(nreduc) == zero .AND. v(preduc) == zero)) i = 2 IF (v(stppar) == zero .AND. v(reldx) <= v(xctol) .AND. goodx) i = i + 1 IF (i > 0) iv(irc) = i + 6 ! *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR ! *** CONVERGENCE TEST. 250 IF (iv(irc) > 5 .AND. iv(irc) /= 12) GO TO 999 IF (v(stppar) == zero) GO TO 999 IF (v(dstnrm) > v(lmaxs)) GO TO 260 IF (v(preduc) >= emaxs) GO TO 999 IF (v(dst0) <= zero) GO TO 270 IF (half * v(dst0) <= v(lmaxs)) GO TO 999 GO TO 270 260 IF (half * v(dstnrm) <= v(lmaxs)) GO TO 999 xmax = v(lmaxs) / v(dstnrm) IF (xmax * (two - xmax) * v(preduc) >= emaxs) GO TO 999 270 IF (v(nreduc) < zero) GO TO 290 ! *** RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST *** v(gtslst) = v(gtstep) v(dstsav) = v(dstnrm) IF (iv(irc) == 12) v(dstsav) = -v(dstsav) v(plstgd) = v(preduc) i = iv(restor) iv(restor) = 2 IF (i == 3) iv(restor) = 0 iv(irc) = 6 GO TO 999 ! *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) *** 280 v(gtstep) = v(gtslst) v(dstnrm) = ABS(v(dstsav)) iv(irc) = iv(xirc) IF (v(dstsav) <= zero) iv(irc) = 12 v(nreduc) = -v(preduc) v(preduc) = v(plstgd) iv(restor) = 3 290 IF (-v(nreduc) <= v(sctol) * ABS(v(f0))) iv(irc) = 11 999 RETURN ! *** LAST LINE OF A7SST FOLLOWS *** END SUBROUTINE a7sst FUNCTION d7tpr(p, x, y) RESULT(fn_val) ! *** RETURN THE INNER PRODUCT OF THE P-VECTORS X AND Y. *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: x(:) REAL (dp), INTENT(IN) :: y(:) REAL (dp) :: fn_val INTEGER :: i REAL (dp) :: t ! REAL :: r7mdc ! EXTERNAL r7mdc ! *** R7MDC(2) RETURNS A MACHINE-DEPENDENT CONSTANT, SQTETA, WHICH ! *** IS SLIGHTLY LARGER THAN THE SMALLEST POSITIVE NUMBER THAT ! *** CAN BE SQUARED WITHOUT UNDERFLOWING. REAL (dp), SAVE :: sqteta = 0.0_dp fn_val = zero IF (p <= 0) GO TO 999 IF (sqteta == zero) sqteta = r7mdc(2) DO i = 1, p t = MAX( ABS(x(i)), ABS(y(i))) IF (t > one) GO TO 10 IF (t < sqteta) CYCLE t = (x(i)/sqteta)*y(i) IF ( ABS(t) < sqteta) CYCLE 10 fn_val = fn_val + x(i)*y(i) END DO 999 RETURN ! *** LAST LINE OF D7TPR FOLLOWS *** END FUNCTION d7tpr SUBROUTINE f7hes(d, g, irt, iv, p, v, x) ! *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING ! *** AT V(IV(FDH)) = V(-IV(H)). ! *** IF IV(COVREQ) >= 0 THEN F7HES USES GRADIENT DIFFERENCES, ! *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN G7LIT. ! IRT VALUES... ! 1 = COMPUTE FUNCTION VALUE, I.E., V(F). ! 2 = COMPUTE G. ! 3 = DONE. ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(IN OUT) :: g(:) INTEGER, INTENT(OUT) :: irt INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: v(:) REAL (dp), INTENT(IN OUT) :: x(:) ! *** LOCAL VARIABLES *** INTEGER :: gsave1, hes, hmi, hpi, hpm, i, k, kind, l, m, mm1, mm1o2, & pp1o2, stpi, stpm, stp0 REAL (dp) :: del REAL (dp), PARAMETER :: half=0.5_dp, negpt5=-0.5_dp, two=2._dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ irt = 4 kind = iv(covreq) m = iv(mode) IF (m > 0) GO TO 10 iv(h) = -ABS(iv(h)) iv(fdh) = 0 iv(kagqt) = -1 v(fx) = v(f) 10 IF (m > p) GO TO 999 IF (kind < 0) GO TO 110 ! *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND ! *** GRADIENT VALUES. gsave1 = iv(w0) + p IF (m > 0) GO TO 20 ! *** FIRST CALL ON F7HES. SET GSAVE = G, TAKE FIRST STEP *** v(gsave1:gsave1+p-1) = g(1:p) iv(switch) = iv(nfgcal) GO TO 90 20 del = v(delta) x(m) = v(xmsave) IF (iv(toobig) == 0) GO TO 40 ! *** HANDLE OVERSIZE V(DELTA) *** IF (del*x(m) > zero) GO TO 30 ! *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** iv(fdh) = -2 GO TO 220 ! *** TRY SHRINKING V(DELTA) *** 30 del = negpt5 * del GO TO 100 40 hes = -iv(h) ! *** SET G = (G - GSAVE)/DEL *** DO i = 1, p g(i) = (g(i) - v(gsave1)) / del gsave1 = gsave1 + 1 END DO ! *** ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** k = hes + m*(m-1)/2 l = k + m - 2 IF (m == 1) GO TO 70 ! *** SET H(I,M) = 0.5 * (H(I,M) + G(I)) FOR I = 1 TO M-1 *** mm1 = m - 1 DO i = 1, mm1 v(k) = half * (v(k) + g(i)) k = k + 1 END DO ! *** ADD H(I,M) = G(I) FOR I = M TO P *** 70 l = l + 1 DO i = m, p v(l) = g(i) l = l + i END DO 90 m = m + 1 iv(mode) = m IF (m > p) GO TO 210 ! *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE *** del = v(delta0) * MAX(one/d(m), ABS(x(m))) IF (x(m) < zero) del = -del v(xmsave) = x(m) 100 x(m) = x(m) + del v(delta) = del irt = 2 GO TO 999 ! *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. 110 stp0 = iv(w0) + p - 1 mm1 = m - 1 mm1o2 = m*mm1/2 IF (m > 0) GO TO 120 ! *** FIRST CALL ON F7HES. *** iv(savei) = 0 GO TO 200 120 i = iv(savei) hes = -iv(h) IF (i > 0) GO TO 180 IF (iv(toobig) == 0) GO TO 140 ! *** HANDLE OVERSIZE STEP *** stpm = stp0 + m del = v(stpm) IF (del*x(xmsave) > zero) GO TO 130 ! *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** iv(fdh) = -2 GO TO 220 ! *** TRY SHRINKING THE STEP *** 130 del = negpt5 * del x(m) = x(xmsave) + del v(stpm) = del irt = 1 GO TO 999 ! *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** 140 pp1o2 = p * (p-1) / 2 hpm = hes + pp1o2 + mm1 v(hpm) = v(f) ! *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** hmi = hes + mm1o2 IF (mm1 == 0) GO TO 160 hpi = hes + pp1o2 DO i = 1, mm1 v(hmi) = v(fx) - (v(f) + v(hpi)) hmi = hmi + 1 hpi = hpi + 1 END DO 160 v(hmi) = v(f) - two*v(fx) ! *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** i = 1 170 iv(savei) = i stpi = stp0 + i v(delta) = x(i) x(i) = x(i) + v(stpi) IF (i == m) x(i) = v(xmsave) - v(stpi) irt = 1 GO TO 999 180 x(i) = v(delta) IF (iv(toobig) == 0) GO TO 190 ! *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** iv(fdh) = -2 GO TO 220 ! *** FINISH COMPUTING H(M,I) *** 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 <= m) GO TO 170 iv(savei) = 0 x(m) = v(xmsave) 200 m = m + 1 iv(mode) = m IF (m > p) GO TO 210 ! *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. ! *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN ! *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STANDARD UNIT VECTOR. del = v(dltfdc) * MAX(one/d(m), ABS(x(m))) IF (x(m) < zero) del = -del v(xmsave) = x(m) x(m) = x(m) + del stpm = stp0 + m v(stpm) = del irt = 1 GO TO 999 ! *** RESTORE V(F), ETC. *** 210 iv(fdh) = hes 220 v(f) = v(fx) irt = 3 IF (kind < 0) GO TO 999 iv(nfgcal) = iv(switch) gsave1 = iv(w0) + p g(1:p) = v(gsave1:gsave1+p-1) 999 RETURN ! *** LAST CARD OF F7HES FOLLOWS *** END SUBROUTINE f7hes FUNCTION v2nrm(p, x) RESULT(fn_val) ! *** RETURN THE 2-NORM OF THE P-VECTOR X, TAKING *** ! *** CARE TO AVOID THE MOST LIKELY UNDERFLOWS. *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: x(:) REAL (dp) :: fn_val INTEGER :: i, j REAL (dp) :: r, scale, t, xi ! REAL :: SQRT ! REAL :: r7mdc ! EXTERNAL r7mdc REAL (dp), SAVE :: sqteta = 0.0_dp IF (p > 0) GO TO 10 fn_val = zero GO TO 999 10 DO i = 1, p IF (x(i) /= zero) GO TO 30 END DO fn_val = zero GO TO 999 30 scale = ABS(x(i)) IF (i < p) GO TO 40 fn_val = scale GO TO 999 40 t = one IF (sqteta == zero) sqteta = r7mdc(2) ! *** SQTETA IS (SLIGHTLY LARGER THAN) THE SQUARE ROOT OF THE ! *** SMALLEST POSITIVE FLOATING POINT NUMBER ON THE MACHINE. ! *** THE TESTS INVOLVING SQTETA ARE DONE TO PREVENT UNDERFLOWS. j = i + 1 DO i = j, p xi = ABS(x(i)) IF (xi > scale) GO TO 50 r = xi / scale IF (r > sqteta) t = t + r*r CYCLE 50 r = scale / xi IF (r <= sqteta) r = zero t = one + t * r*r scale = xi END DO fn_val = scale * SQRT(t) 999 RETURN ! *** LAST LINE OF V2NRM FOLLOWS *** END FUNCTION v2nrm SUBROUTINE v2axy(p, w, a, x, y) ! *** SET W = A*X + Y -- W, X, Y = P-VECTORS, A = SCALAR *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: w(:) REAL (dp), INTENT(IN) :: a REAL (dp), INTENT(IN) :: x(:) REAL (dp), INTENT(IN) :: y(:) w(1:p) = a*x(1:p) + y(1:p) RETURN END SUBROUTINE v2axy SUBROUTINE l7itv(n, x, l, y) ! *** SOLVE (L**T)*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR MATRIX ! *** STORED COMPACTLY BY ROWS. *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: l(:) REAL (dp), INTENT(IN) :: y(:) INTEGER :: i, ii, ij, im1, i0, j, np1 REAL (dp) :: xi x(1:n) = y(1:n) np1 = n + 1 i0 = n*(n+1)/2 DO ii = 1, n i = np1 - ii xi = x(i)/l(i0) x(i) = xi IF (i <= 1) EXIT i0 = i0 - i IF (xi == zero) CYCLE im1 = i - 1 DO j = 1, im1 ij = i0 + j x(j) = x(j) - xi*l(ij) END DO END DO RETURN ! *** LAST CARD OF L7ITV FOLLOWS *** END SUBROUTINE l7itv SUBROUTINE l7ivm(n, x, l, y) ! *** SOLVE L*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR ! *** MATRIX STORED COMPACTLY BY ROWS. *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: l(:) REAL (dp), INTENT(IN) :: y(:) ! REAL :: d7tpr ! EXTERNAL d7tpr INTEGER :: i, j, k REAL (dp) :: t DO k = 1, n IF (y(k) /= zero) GO TO 20 x(k) = zero END DO GO TO 999 20 j = k*(k+1)/2 x(k) = y(k) / l(j) IF (k >= n) GO TO 999 k = k + 1 DO i = k, n t = d7tpr(i-1, l(j+1:), x) j = j + i x(i) = (y(i) - t)/l(j) END DO 999 RETURN ! *** LAST CARD OF L7IVM FOLLOWS *** END SUBROUTINE l7ivm SUBROUTINE l7srt(n1, n, l, a, irc) ! *** COMPUTE ROWS N1 THROUGH N OF THE CHOLESKY FACTOR L OF ! *** A = L*(L**T), WHERE L AND THE LOWER TRIANGLE OF A ARE BOTH ! *** STORED COMPACTLY BY ROWS (AND MAY OCCUPY THE SAME STORAGE). ! *** IRC = 0 MEANS ALL WENT WELL. IRC = J MEANS THE LEADING ! *** PRINCIPAL J X J SUBMATRIX OF A IS NOT POSITIVE DEFINITE -- ! *** AND L(J*(J+1)/2) CONTAINS THE (NONPOS.) REDUCED J-TH DIAGONAL. ! *** PARAMETERS *** INTEGER, INTENT(IN) :: n1 INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: l(:) REAL (dp), INTENT(IN) :: a(:) INTEGER, INTENT(OUT) :: irc ! DIMENSION L(N*(N+1)/2), A(N*(N+1)/2) ! *** LOCAL VARIABLES *** INTEGER :: i, ij, ik, im1, i0, j, jk, jm1, j0, k REAL (dp) :: t, td ! *** INTRINSIC FUNCTIONS *** ! REAL :: SQRT ! *** BODY *** i0 = n1 * (n1 - 1) / 2 DO i = n1, n td = zero IF (i == 1) GO TO 40 j0 = 0 im1 = i - 1 DO j = 1, im1 t = zero IF (j == 1) GO TO 20 jm1 = j - 1 DO k = 1, jm1 ik = i0 + k jk = j0 + k t = t + l(ik)*l(jk) END DO 20 ij = i0 + j j0 = j0 + j t = (a(ij) - t) / l(j0) l(ij) = t td = td + t*t END DO 40 i0 = i0 + i t = a(i0) - td IF (t <= zero) GO TO 60 l(i0) = SQRT(t) END DO irc = 0 GO TO 999 60 l(i0) = t irc = i 999 RETURN ! *** LAST CARD OF L7SRT *** END SUBROUTINE l7srt SUBROUTINE l7svn(p, l, x, y, smallest) ! *** ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L ! *** PARAMETER DECLARATIONS *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: l(:) REAL (dp), INTENT(OUT) :: x(:), y(:), smallest ! DIMENSION L(P*(P+1)/2) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** PURPOSE *** ! THIS FUNCTION RETURNS A GOOD OVER-ESTIMATE OF THE SMALLEST ! SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L. ! *** PARAMETER DESCRIPTION *** ! P (IN) = THE ORDER OF L. L IS A P X P LOWER TRIANGULAR MATRIX. ! L (IN) = ARRAY HOLDING THE ELEMENTS OF L IN ROW ORDER, I.E. ! L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC. ! X (OUT) IF smallest RETURNS A POSITIVE VALUE, THEN X IS A NORMALIZED ! APPROXIMATE LEFT SINGULAR VECTOR CORRESPONDING TO THE ! SMALLEST SINGULAR VALUE. THIS APPROXIMATION MAY BE VERY ! CRUDE. IF smallest RETURNS ZERO, THEN SOME COMPONENTS OF X ! ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES. ! Y (OUT) IF smallest RETURNS A POSITIVE VALUE, THEN Y = (L**-1)*X IS AN ! UNNORMALIZED APPROXIMATE RIGHT SINGULAR VECTOR CORRESPOND- ! ING TO THE SMALLEST SINGULAR VALUE. THIS APPROXIMATION ! MAY BE CRUDE. IF smallest RETURNS ZERO, THEN Y RETAINS ITS ! INPUT VALUE. THE CALLER MAY PASS THE SAME VECTOR FOR X ! AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE Y OVER- ! WRITES X (FOR NONZERO smallest RETURNS). ! *** ALGORITHM NOTES *** ! THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT ! smallest = 0 IS RETURNED IF THE SMALLEST DIAGONAL ELEMENT OF L ! (IN MAGNITUDE) IS NOT MORE THAN THE UNIT ROUNDOFF TIMES THE LARGEST. ! THE ALGORITHM USES A RANDOM NUMBER GENERATOR PROPOSED IN (4), WHICH ! PASSES THE SPECTRAL TEST WITH FLYING COLORS -- SEE (2) AND (3). ! *** SUBROUTINES AND FUNCTIONS CALLED *** ! V2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR. ! *** REFERENCES *** ! (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977), ! AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT ! TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY. ! (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL ! RANDOM-NUMBER GENERATORS -- AN EMPIRICAL VIEW, ! MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. ! (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 ! (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. ! (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER GENERATORS ! WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, PP. 586-593. ! *** HISTORY *** ! DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978). ! *** GENERAL *** ! THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH ! SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS ! MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** LOCAL VARIABLES *** INTEGER :: i, ii, ix, j, ji, jj, jm1, j0, pm1 REAL (dp) :: b, sminus, splus, t, xminus, xplus ! *** EXTERNAL FUNCTIONS AND SUBROUTINES *** ! REAL :: d7tpr, v2nrm ! EXTERNAL d7tpr, v2nrm, v2axy REAL (dp), PARAMETER :: half=0.5_dp, r9973=9973._dp ! *** BODY *** ix = 2 pm1 = p - 1 ! *** FIRST CHECK WHETHER TO RETURN smallest = 0 AND INITIALIZE X *** ii = 0 j0 = p*pm1/2 jj = j0 + p IF (l(jj) == zero) GO TO 110 ix = MOD(3432*ix, 9973) b = half*(one + REAL(ix)/r9973) xplus = b / l(jj) x(p) = xplus IF (p <= 1) GO TO 60 DO i = 1, pm1 ii = ii + i IF (l(ii) == zero) CYCLE ji = j0 + i x(i) = xplus * l(ji) END DO ! *** SOLVE (L**T)*X = B, WHERE THE COMPONENTS OF B HAVE RANDOMLY ! *** CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE. ! DO J = P-1 TO 1 BY -1... DO j = p-1, 1, -1 ! *** DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J ! *** THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I. ix = MOD(3432*ix, 9973) b = half*(one + REAL(ix)/r9973) xplus = (b - x(j)) xminus = (-b - x(j)) splus = ABS(xplus) sminus = ABS(xminus) jm1 = j - 1 j0 = j*jm1/2 jj = j0 + j xplus = xplus/l(jj) xminus = xminus/l(jj) DO i = 1, jm1 ji = j0 + i splus = splus + ABS(x(i) + l(ji)*xplus) sminus = sminus + ABS(x(i) + l(ji)*xminus) END DO IF (sminus > splus) xplus = xminus x(j) = xplus ! *** UPDATE PARTIAL SUMS *** IF (jm1 > 0) CALL v2axy(jm1, x, xplus, l(j0+1:), x) END DO ! *** NORMALIZE X *** 60 t = one/ v2nrm(p, x) DO i = 1, p x(i) = t*x(i) END DO ! *** SOLVE L*Y = X AND RETURN smallest = 1/TWONORM(Y) *** DO j = 1, p jm1 = j - 1 j0 = j*jm1/2 jj = j0 + j t = zero IF (jm1 > 0) t = d7tpr(jm1, l(j0+1:), y) y(j) = (x(j) - t) / l(jj) END DO smallest = one/ v2nrm(p, y) GO TO 999 110 smallest = zero 999 RETURN ! *** LAST CARD OF L7SVN FOLLOWS *** END SUBROUTINE l7svn SUBROUTINE g7qts(d, dig, dihdi, ka, l, p, step, v) ! *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE *** ! *** (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN *** ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(IN) :: dig(:) REAL (dp), INTENT(IN OUT) :: dihdi(:) INTEGER, INTENT(IN OUT) :: ka REAL (dp), INTENT(IN OUT) :: l(:) INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: step(:) REAL (dp), INTENT(IN OUT) :: v(:) ! DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** PURPOSE *** ! GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED ! HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR, ! THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF ! APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE. IN ! OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE ! PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP SUCH THAT THE ! 2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE ! G IS THE GRADIENT, H IS THE HESSIAN, AND D IS A DIAGONAL ! SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D. ! (G7QTS ASSUMES DIG = D**-1 * G AND DIHDI = D**-1 * H * D**-1.) ! *** PARAMETER DESCRIPTION *** ! D (IN) = THE SCALE VECTOR, I.E. THE DIAGONAL OF THE SCALE ! MATRIX D MENTIONED ABOVE UNDER PURPOSE. ! DIG (IN) = THE SCALED GRADIENT VECTOR, D**-1 * G. IF G = 0, THEN ! STEP = 0 AND V(STPPAR) = 0 ARE RETURNED. ! DIHDI (IN) = LOWER TRIANGLE OF THE SCALED HESSIAN (APPROXIMATION), ! I.E., D**-1 * H * D**-1, STORED COMPACTLY BY ROWS., I.E., ! IN THE ORDER (1,1), (2,1), (2,2), (3,1), (3,2), ETC. ! KA (I/O) = THE NUMBER OF HEBDEN ITERATIONS (SO FAR) TAKEN TO DETER- ! MINE STEP. KA < 0 ON INPUT MEANS THIS IS THE FIRST ! ATTEMPT TO DETERMINE STEP (FOR THE PRESENT DIG AND DIHDI) ! -- KA IS INITIALIZED TO 0 IN THIS CASE. OUTPUT WITH ! KA = 0 (OR V(STPPAR) = 0) MEANS STEP = -(H**-1)*G. ! L (I/O) = WORKSPACE OF LENGTH P*(P+1)/2 FOR CHOLESKY FACTORS. ! P (IN) = NUMBER OF PARAMETERS -- THE HESSIAN IS A P X P MATRIX. ! STEP (I/O) = THE STEP COMPUTED. ! V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW. ! W (I/O) = WORKSPACE OF LENGTH 4*P + 6. ! *** ENTRIES IN V *** ! V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G. ! V(DSTNRM) (OUTPUT) = 2-NORM OF D*STEP. ! V(DST0) (I/O) = 2-NORM OF D*(H**-1)*G (FOR POS. DEF. H ONLY), OR ! OVERESTIMATE OF SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1). ! V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED FOR PSI(STEP). FOR THE ! STEP RETURNED, PSI(STEP) WILL EXCEED ITS OPTIMAL VALUE ! BY LESS THAN -V(EPSLON)*PSI(STEP). SUGGESTED VALUE = 0.1. ! V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP. ! V(NREDUC) (OUT) = PSI(-(H**-1)*G) = PSI(NEWTON STEP) (FOR POS. DEF. ! H ONLY -- V(NREDUC) IS SET TO ZERO OTHERWISE). ! V(PHMNFC) (IN) = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP ! (MORE*S SIGMA). THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE ! BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS). ! V(PHMXFC) (IN) (SEE V(PHMNFC).) ! SUGGESTED VALUES -- V(PHMNFC) = -0.25, V(PHMXFC) = 0.5. ! V(PREDUC) (OUT) = PSI(STEP) = PREDICTED OBJ. FUNC. REDUCTION FOR STEP. ! V(RADIUS) (IN) = RADIUS OF CURRENT (SCALED) TRUST REGION. ! V(RAD0) (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL. ! V(STPPAR) (I/O) IS NORMALLY THE MARQUARDT PARAMETER, I.E. THE ALPHA ! DESCRIBED BELOW UNDER ALGORITHM NOTES. IF H + ALPHA*D**2 ! (SEE ALGORITHM NOTES) IS (NEARLY) SINGULAR, HOWEVER, ! THEN V(STPPAR) = -ALPHA. ! *** USAGE NOTES *** ! IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF V(RADIUS), ! THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT WITH ALL PARAMETERS ! UNCHANGED EXCEPT V(RADIUS). (THIS EXPLAINS WHY STEP AND W ARE LISTED AS ! I/O). ON AN INITIAL CALL (ONE WITH KA < 0), STEP AND W NEED NOT BE ! INITIALIZED AND ONLY COMPONENTS V(EPSLON), V(STPPAR), V(PHMNFC), ! V(PHMXFC), V(RADIUS), AND V(RAD0) OF V MUST BE INITIALIZED. ! *** ALGORITHM NOTES *** ! THE DESIRED G-Q-T STEP (REF. 2, 3, 4, 6) SATISFIES ! (H + ALPHA*D**2)*STEP = -G FOR SOME NONNEGATIVE ALPHA SUCH THAT ! H + ALPHA*D**2 IS POSITIVE SEMIDEFINITE. ALPHA AND STEP ARE ! COMPUTED BY A SCHEME ANALOGOUS TO THE ONE DESCRIBED IN REF. 5. ! ESTIMATES OF THE SMALLEST AND LARGEST EIGENVALUES OF THE HESSIAN ! ARE OBTAINED FROM THE GERSCHGORIN CIRCLE THEOREM ENHANCED BY A ! SIMPLE FORM OF THE SCALING DESCRIBED IN REF. 7. CASES IN WHICH ! H + ALPHA*D**2 IS NEARLY (OR EXACTLY) SINGULAR ARE HANDLED BY ! THE TECHNIQUE DISCUSSED IN REF. 2. IN THESE CASES, A STEP OF ! (EXACT) LENGTH V(RADIUS) IS RETURNED FOR WHICH PSI(STEP) EXCEEDS ! ITS OPTIMAL VALUE BY LESS THAN -V(EPSLON)*PSI(STEP). THE TEST ! SUGGESTED IN REF. 6 FOR DETECTING THE SPECIAL CASE IS PERFORMED ! ONCE TWO MATRIX FACTORIZATIONS HAVE BEEN DONE -- DOING SO SOONER ! SEEMS TO DEGRADE THE PERFORMANCE OF OPTIMIZATION ROUTINES THAT ! CALL THIS ROUTINE. ! *** FUNCTIONS AND SUBROUTINES CALLED *** ! D7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. ! L7ITV - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. ! L7IVM - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX. ! L7SRT - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.). ! L7SVN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX. ! R7MDC - RETURNS MACHINE-DEPENDENT CONSTANTS. ! V2NRM - RETURNS 2-NORM OF A VECTOR. ! *** REFERENCES *** ! 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE ! NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. ! SOFTWARE, VOL. 7, NO. 3. ! 2. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS, ! SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP. ! 186-197. ! 3. GOLDFELD, S.M., QUANDT, R.E., AND TROTTER, H.F. (1966), ! MAXIMIZATION BY QUADRATIC HILL-CLIMBING, ECONOMETRICA 34, ! PP. 541-551. ! 4. HEBDEN, M.D. (1973), AN ALGORITHM FOR MINIMIZATION USING EXACT ! SECOND DERIVATIVES, REPORT T.P. 515, THEORETICAL PHYSICS ! DIV., A.E.R.E. HARWELL, OXON., ENGLAND. ! 5. MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN- ! TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES ! IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER- ! VERLAG, BERLIN AND NEW YORK. ! 6. MORE, J.J., AND SORENSEN, D.C. (1981), COMPUTING A TRUST REGION ! STEP, TECHNICAL REPORT ANL-81-83, ARGONNE NATIONAL LAB. ! 7. VARGA, R.S. (1965), MINIMAL GERSCHGORIN SETS, PACIFIC J. MATH. 15, ! PP. 719-729. ! *** GENERAL *** ! CODED BY DAVID M. GAY. ! THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH ! SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS ! MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND MCS-7906671. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** LOCAL VARIABLES *** LOGICAL :: restrt INTEGER :: dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc, & j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x REAL (dp) :: alphak, aki, akk, delta, dst, eps, gtsta, lk, & oldphi, phi, phimax, phimin, psifac, rad, radsq, & root, si, sk, sw, t, twopsi, t1, t2, uk, wi REAL (dp) :: w(4*p+7) ! *** INTRINSIC FUNCTIONS *** ! REAL :: SQRT ! *** EXTERNAL FUNCTIONS AND SUBROUTINES *** ! REAL :: d7tpr, l7svn, r7mdc, v2nrm ! EXTERNAL d7tpr, l7itv, l7ivm, l7srt, l7svn, r7mdc, v2nrm REAL (dp), PARAMETER :: epsfac=50.0_dp, four=4.0_dp, half=0.5_dp, & kappa=2.0_dp, negone=-1.0_dp, & p001=0.001_dp, six=6.0_dp, three=3.0_dp, two=2.0_dp REAL (dp), SAVE :: dgxfac=0.0_dp, big=0.0_dp ! *** BODY *** IF (big <= zero) big = r7mdc(6) ! *** STORE LARGEST ABS. ENTRY IN (D**-1)*H*(D**-1) AT W(DGGDMX). dggdmx = p + 1 ! *** STORE GERSCHGORIN OVER- AND UNDERESTIMATES OF THE LARGEST ! *** AND SMALLEST EIGENVALUES OF (D**-1)*H*(D**-1) AT W(EMAX) ! *** AND W(EMIN) RESPECTIVELY. emax = dggdmx + 1 emin = emax + 1 ! *** FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK, UK, DST, AND ! *** THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR POS. DEF. H) ARE ! *** STORED IN W(LK0), W(UK0), W(DSTSAV), AND W(PHIPIN) RESPECTIVELY. lk0 = emin + 1 phipin = lk0 + 1 uk0 = phipin + 1 dstsav = uk0 + 1 ! *** STORE DIAG OF (D**-1)*H*(D**-1) IN W(DIAG),...,W(DIAG0+P). diag0 = dstsav diag = diag0 + 1 ! *** STORE -D*STEP IN W(Q),...,W(Q0+P). q0 = diag0 + p q = q0 + 1 ! *** ALLOCATE STORAGE FOR SCRATCH VECTOR X *** x = q + p rad = v(radius) radsq = rad**2 ! *** PHITOL = MAX. ERROR ALLOWED IN DST = V(DSTNRM) = 2-NORM OF ! *** D*STEP. phimax = v(phmxfc) * rad phimin = v(phmnfc) * rad psifac = big t1 = two * v(epslon) / (three * (four * (v(phmnfc) + one) * & (kappa + one) + kappa + two) * rad) IF (t1 < big*MIN(rad,one)) psifac = t1 / rad ! *** OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY. IF ! *** WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT. oldphi = zero eps = v(epslon) irc = 0 restrt = .false. kalim = ka + 50 ! *** START OR RESTART, DEPENDING ON KA *** IF (ka >= 0) GO TO 290 ! *** FRESH START *** k = 0 uk = negone ka = 0 kalim = 50 v(dgnorm) = v2nrm(p, dig) v(nreduc) = zero v(dst0) = zero kamin = 3 IF (v(dgnorm) == zero) kamin = 0 ! *** STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P) *** j = 0 DO i = 1, p j = j + i k1 = diag0 + i w(k1) = dihdi(j) END DO ! *** DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI *** t1 = zero j = p * (p + 1) / 2 DO i = 1, j t = ABS(dihdi(i)) IF (t1 < t) t1 = t END DO w(dggdmx) = t1 ! *** TRY ALPHA = 0 *** 30 CALL l7srt(1, p, l, dihdi, irc) IF (irc == 0) GO TO 50 ! *** INDEF. H -- UNDERESTIMATE SMALLEST EIGENVALUE, USE THIS ! *** ESTIMATE TO INITIALIZE LOWER BOUND LK ON ALPHA. j = irc*(irc+1)/2 t = l(j) l(j) = one DO i = 1, irc w(i) = zero END DO w(irc) = one CALL l7itv(irc, w, l, w) t1 = v2nrm(irc, w) lk = -t / t1 / t1 v(dst0) = -lk IF (restrt) GO TO 210 GO TO 70 ! *** POSITIVE DEFINITE H -- COMPUTE UNMODIFIED NEWTON STEP. *** 50 lk = zero CALL l7svn(p, l, w(q:), w(q:), t) IF (t >= one) GO TO 60 IF (v(dgnorm) >= t*t*big) GO TO 70 60 CALL l7ivm(p, w(q:), l, dig) gtsta = d7tpr(p, w(q:), w(q:)) v(nreduc) = half * gtsta CALL l7itv(p, w(q:), l, w(q:)) dst = v2nrm(p, w(q:)) v(dst0) = dst phi = dst - rad IF (phi <= phimax) GO TO 260 IF (restrt) GO TO 210 ! *** PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND ! *** SMALLEST) EIGENVALUES. *** 70 k = 0 DO i = 1, p wi = zero IF (i == 1) GO TO 90 im1 = i - 1 DO j = 1, im1 k = k + 1 t = ABS(dihdi(k)) wi = wi + t w(j) = w(j) + t END DO 90 w(i) = wi k = k + 1 END DO ! *** (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1) *** k = 1 t1 = w(diag) - w(1) IF (p <= 1) GO TO 120 DO i = 2, p j = diag0 + i t = w(j) - w(i) IF (t >= t1) CYCLE t1 = t k = i END DO 120 sk = w(k) j = diag0 + k akk = w(j) k1 = k*(k-1)/2 + 1 inc = 1 t = zero DO i = 1, p IF (i == k) GO TO 130 aki = ABS(dihdi(k1)) si = w(i) j = diag0 + i t1 = half * (akk - w(j) + si - aki) t1 = t1 + SQRT(t1*t1 + sk*aki) IF (t < t1) t = t1 IF (i < k) GO TO 140 130 inc = i 140 k1 = k1 + inc END DO w(emin) = akk - t uk = v(dgnorm)/rad - w(emin) IF (v(dgnorm) == zero) uk = uk + p001 + p001*uk IF (uk <= zero) uk = p001 ! *** COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE *** k = 1 t1 = w(diag) + w(1) IF (p <= 1) GO TO 170 DO i = 2, p j = diag0 + i t = w(j) + w(i) IF (t <= t1) CYCLE t1 = t k = i END DO 170 sk = w(k) j = diag0 + k akk = w(j) k1 = k*(k-1)/2 + 1 inc = 1 t = zero DO i = 1, p IF (i == k) GO TO 180 aki = ABS(dihdi(k1)) si = w(i) j = diag0 + i t1 = half * (w(j) + si - aki - akk) t1 = t1 + SQRT(t1*t1 + sk*aki) IF (t < t1) t = t1 IF (i < k) GO TO 190 180 inc = i 190 k1 = k1 + inc END DO w(emax) = akk + t lk = MAX(lk, v(dgnorm)/rad - w(emax)) ! *** ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE). WE ! *** USE MORE*S SCHEME FOR INITIALIZING IT. alphak = ABS(v(stppar)) * v(rad0)/rad alphak = MIN(uk, MAX(alphak, lk)) IF (irc /= 0) GO TO 210 ! *** COMPUTE L0 FOR POSITIVE DEFINITE H *** CALL l7ivm(p, w, l, w(q:)) t = v2nrm(p, w) w(phipin) = rad / t / t lk = MAX(lk, phi*w(phipin)) ! *** SAFEGUARD ALPHAK AND ADD ALPHAK*I TO (D**-1)*H*(D**-1) *** 210 ka = ka + 1 IF (-v(dst0) >= alphak .OR. alphak < lk .OR. alphak >= uk) & alphak = uk * MAX(p001, SQRT(lk/uk)) IF (alphak <= zero) alphak = half * uk IF (alphak <= zero) alphak = uk k = 0 DO i = 1, p k = k + i j = diag0 + i dihdi(k) = w(j) + alphak END DO ! *** TRY COMPUTING CHOLESKY DECOMPOSITION *** CALL l7srt(1, p, l, dihdi, irc) IF (irc == 0) GO TO 240 ! *** (D**-1)*H*(D**-1) + ALPHAK*I IS INDEFINITE -- OVERESTIMATE ! *** SMALLEST EIGENVALUE FOR USE IN UPDATING LK *** j = (irc*(irc+1))/2 t = l(j) l(j) = one DO i = 1, irc w(i) = zero END DO w(irc) = one CALL l7itv(irc, w, l, w) t1 = v2nrm(irc, w) lk = alphak - t/t1/t1 v(dst0) = -lk IF (uk < lk) uk = lk IF (alphak < lk) GO TO 210 ! *** NASTY CASE -- EXACT GERSCHGORIN BOUNDS. FUDGE LK, UK... t = p001 * alphak IF (t <= zero) t = p001 lk = alphak + t IF (uk <= lk) uk = lk + t GO TO 210 ! *** ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE. ! *** COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE. *** 240 CALL l7ivm(p, w(q:), l, dig) gtsta = d7tpr(p, w(q:), w(q:)) CALL l7itv(p, w(q:), l, w(q:)) dst = v2nrm(p, w(q:)) phi = dst - rad IF (phi <= phimax .AND. phi >= phimin) GO TO 270 IF (phi == oldphi) GO TO 270 oldphi = phi IF (phi < zero) GO TO 330 ! *** UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK *** 250 IF (ka >= kalim) GO TO 270 ! *** THE FOLLOWING MIN IS NECESSARY BECAUSE OF RESTARTS *** IF (phi < zero) uk = MIN(uk, alphak) ! *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES *** IF (kamin == 0) GO TO 210 CALL l7ivm(p, w, l, w(q:)) ! *** THE FOLLOWING, COMMENTED CALCULATION OF ALPHAK IS SOMETIMES ! *** SAFER BUT WORSE IN PERFORMANCE... ! T1 = DST / V2NRM(P, W) ! ALPHAK = ALPHAK + T1 * (PHI/RAD) * T1 t1 = v2nrm(p, w) alphak = alphak + (phi/t1) * (dst/t1) * (dst/rad) lk = MAX(lk, alphak) alphak = lk GO TO 210 ! *** ACCEPTABLE STEP ON FIRST TRY *** 260 alphak = zero ! *** SUCCESSFUL STEP IN GENERAL. COMPUTE STEP = -(D**-1)*Q *** 270 DO i = 1, p j = q0 + i step(i) = -w(j)/d(i) END DO v(gtstep) = -gtsta v(preduc) = half * ( ABS(alphak)*dst*dst + gtsta) GO TO 410 ! *** RESTART WITH NEW RADIUS *** 290 IF (v(dst0) <= zero .OR. v(dst0) - rad > phimax) GO TO 310 ! *** PREPARE TO RETURN NEWTON STEP *** restrt = .true. ka = ka + 1 k = 0 DO i = 1, p k = k + i j = diag0 + i dihdi(k) = w(j) END DO uk = negone GO TO 30 310 kamin = ka + 3 IF (v(dgnorm) == zero) kamin = 0 IF (ka == 0) GO TO 50 dst = w(dstsav) alphak = ABS(v(stppar)) phi = dst - rad t = v(dgnorm)/rad uk = t - w(emin) IF (v(dgnorm) == zero) uk = uk + p001 + p001*uk IF (uk <= zero) uk = p001 IF (rad > v(rad0)) GO TO 320 ! *** SMALLER RADIUS *** lk = zero IF (alphak > zero) lk = w(lk0) lk = MAX(lk, t - w(emax)) IF (v(dst0) > zero) lk = MAX(lk, (v(dst0)-rad)*w(phipin)) GO TO 250 ! *** BIGGER RADIUS *** 320 IF (alphak > zero) uk = MIN(uk, w(uk0)) lk = MAX(zero, -v(dst0), t - w(emax)) IF (v(dst0) > zero) lk = MAX(lk, (v(dst0)-rad)*w(phipin)) GO TO 250 ! *** DECIDE WHETHER TO CHECK FOR SPECIAL CASE... IN PRACTICE (FROM THE ! *** STANDPOINT OF THE CALLING OPTIMIZATION CODE) IT SEEMS BEST NOT TO CHECK ! *** UNTIL A FEW ITERATIONS HAVE FAILED -- HENCE THE TEST ON KAMIN BELOW. 330 delta = alphak + MIN(zero, v(dst0)) twopsi = alphak*dst*dst + gtsta IF (ka >= kamin) GO TO 340 ! *** IF THE TEST IN REF. 2 IS SATISFIED, FALL THROUGH TO HANDLE THE ! *** SPECIAL CASE (AS SOON AS THE MORE-SORENSEN TEST DETECTS IT). IF (psifac >= big) GO TO 340 IF (delta >= psifac*twopsi) GO TO 370 ! *** CHECK FOR THE SPECIAL CASE OF H + ALPHA*D**2 (NEARLY) SINGULAR. ! *** USE ONE STEP OF INVERSE POWER METHOD WITH START FROM L7SVN TO OBTAIN ! *** APPROXIMATE EIGENVECTOR CORRESPONDING TO SMALLEST EIGENVALUE OF ! *** (D**-1)*H*(D**-1). L7SVN RETURNS X AND W WITH L*W = X. 340 CALL l7svn(p, l, w(x:), w, t) ! *** NORMALIZE W *** w(1:p) = t*w(1:p) ! *** COMPLETE CURRENT INV. POWER ITER. -- REPLACE W BY (L**-T)*W. CALL l7itv(p, w, l, w) t2 = one/ v2nrm(p, w) w(1:p) = t2*w(1:p) t = t2 * t ! *** NOW W IS THE DESIRED APPROXIMATE (UNIT) EIGENVECTOR AND ! *** T*X = ((D**-1)*H*(D**-1) + ALPHAK*I)*W. sw = d7tpr(p, w(q:), w) t1 = (rad + dst) * (rad - dst) root = SQRT(sw*sw + t1) IF (sw < zero) root = -root si = t1 / (sw + root) ! *** THE ACTUAL TEST FOR THE SPECIAL CASE... IF ((t2*si)**2 <= eps*(dst**2 + alphak*radsq)) GO TO 380 ! *** UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE) ! *** (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE... IF (v(dst0) <= zero) v(dst0) = MIN(v(dst0), t2**2 - alphak) lk = MAX(lk, -v(dst0)) ! *** CHECK WHETHER WE CAN HOPE TO DETECT THE SPECIAL CASE IN ! *** THE AVAILABLE ARITHMETIC. ACCEPT STEP AS IT IS IF NOT. ! *** IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC. 370 IF (dgxfac == zero) dgxfac = epsfac * r7mdc(3) IF (delta > dgxfac*w(dggdmx)) GO TO 250 GO TO 270 ! *** SPECIAL CASE DETECTED... NEGATE ALPHAK TO INDICATE SPECIAL CASE 380 alphak = -alphak v(preduc) = half * twopsi ! *** ACCEPT CURRENT STEP IF ADDING SI*W WOULD LEAD TO A ! *** FURTHER RELATIVE REDUCTION IN PSI OF LESS THAN V(EPSLON)/3. t1 = zero t = si*(alphak*sw - half*si*(alphak + t*d7tpr(p,w(x:),w))) IF (t < eps*twopsi/six) GO TO 390 v(preduc) = v(preduc) + t dst = rad t1 = -si 390 DO i = 1, p j = q0 + i w(j) = t1*w(i) - w(j) step(i) = w(j) / d(i) END DO v(gtstep) = d7tpr(p, dig, w(q:)) ! *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** 410 v(dstnrm) = dst v(stppar) = alphak w(lk0) = lk w(uk0) = uk v(rad0) = rad w(dstsav) = dst ! *** RESTORE DIAGONAL OF DIHDI *** j = 0 DO i = 1, p j = j + i k = diag0 + i dihdi(j) = w(k) END DO RETURN ! *** LAST CARD OF G7QTS FOLLOWS *** END SUBROUTINE g7qts SUBROUTINE itsum (d, g, iv, p, v, x, b) ! *** PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3) *** ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(IN) :: g(:) INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: v(:) REAL (dp), INTENT(IN) :: x(:) REAL (dp), INTENT(IN), OPTIONAL :: b(:,:) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** LOCAL VARIABLES *** INTEGER :: alg, i, iv1, m, nf, ng, ol, pu REAL (dp) :: nreldf, oldf, preldf, reldf ! *** NO EXTERNAL FUNCTIONS OR SUBROUTINES *** CHARACTER (LEN=4), PARAMETER :: model1(6) = (/ ' ', ' ', ' ', & ' ', ' G ', ' S ' /), & model2(6) = (/ ' G ', ' S ', 'G-S ', & 'S-G ', '-S-G', '-G-S' /) !------------------------------- BODY -------------------------------- pu = iv(prunit) IF (pu == 0) GO TO 520 iv1 = iv(1) IF (iv1 > 62) iv1 = iv1 - 51 ol = iv(outlev) alg = MOD(iv(algsav)-1,2) + 1 IF (iv1 < 2 .OR. iv1 > 15) GO TO 370 IF (iv1 >= 12) GO TO 120 IF (iv1 == 2 .AND. iv(niter) == 0) GO TO 390 IF (ol == 0) GO TO 120 IF (iv1 >= 10 .AND. iv(prntit) == 0) GO TO 120 IF (iv1 > 2) GO TO 10 iv(prntit) = iv(prntit) + 1 IF (iv(prntit) < ABS(ol)) GO TO 520 10 nf = iv(nfcall) - ABS(iv(nfcov)) iv(prntit) = 0 reldf = zero preldf = zero oldf = MAX(ABS(v(f0)),ABS(v(f))) IF (oldf <= zero) GO TO 20 reldf = v(fdif)/oldf preldf = v(preduc)/oldf 20 IF (ol > 0) GO TO 60 ! *** PRINT SHORT SUMMARY LINE *** IF (iv(needhd) == 1 .AND. alg == 1) WRITE (pu,30) 30 FORMAT (/' IT NF F RELDF PRELDF RELDX MODEL STPPAR') IF (iv(needhd) == 1.AND.alg == 2) WRITE (pu,40) 40 FORMAT (/' IT NF F RELDF PRELDF RELDX STPPAR') iv(needhd)=0 IF (alg == 2) GO TO 50 m=iv(sused) WRITE (pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx), model1(m), & model2(m), v(stppar) GO TO 120 50 WRITE (pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx), v(stppar) GO TO 120 ! *** PRINT LONG SUMMARY LINE *** 60 IF (iv(needhd) == 1 .AND. alg == 1) WRITE (pu,70) 70 FORMAT (/' IT NF F RELDF PRELDF RELDX MODEL STPPAR', & ' D*STEP NPRELDF') IF (iv(needhd) == 1 .AND. alg == 2) WRITE (pu,80) 80 FORMAT (/' IT NF F RELDF PRELDF RELDX STPPAR', & ' D*STEP NPRELDF') iv(needhd)=0 nreldf=zero IF (oldf > zero) nreldf=v(nreduc)/oldf IF (alg == 2) GO TO 90 m=iv(sused) WRITE (pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx), model1(m), & model2(m), v(stppar), v(dstnrm), nreldf GO TO 120 90 WRITE (pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx), v(stppar), & v(dstnrm), nreldf 100 FORMAT (i6, i5, e10.3, 2E9.2, e8.1, a3, a4, 2E8.1, e9.2) 110 FORMAT (i6, i5, e11.3, 2E10.2, 3E9.1, e10.2) 120 IF (iv1 <= 2) GO TO 520 i=iv(statpr) IF (i == -1) GO TO 460 IF (i+iv1 < 0) GO TO 460 SELECT CASE (iv1) CASE (1:2) GO TO 520 CASE (3) WRITE (pu,140) 140 FORMAT (/' ***** X-CONVERGENCE *****') GO TO 430 CASE (4) WRITE (pu,160) 160 FORMAT (/' ***** RELATIVE FUNCTION CONVERGENCE *****') GO TO 430 CASE (5) WRITE (pu,180) 180 FORMAT (/' ***** X- AND RELATIVE FUNCTION CONVERGENCE *****') GO TO 430 CASE (6) WRITE (pu,200) 200 FORMAT (/' ***** ABSOLUTE FUNCTION CONVERGENCE *****') GO TO 430 CASE (7) WRITE (pu,220) 220 FORMAT (/' ***** SINGULAR CONVERGENCE *****') GO TO 430 CASE (8) WRITE (pu,240) 240 FORMAT (/' ***** FALSE CONVERGENCE *****') GO TO 430 CASE (9) WRITE (pu,260) 260 FORMAT (/' ***** FUNCTION EVALUATION LIMIT *****') GO TO 430 CASE (10) WRITE (pu,280) 280 FORMAT (/' ***** ITERATION LIMIT *****') GO TO 430 CASE (11) WRITE (pu,300) 300 FORMAT (/' ***** STOPX *****') GO TO 430 CASE (12) WRITE (pu,320) 320 FORMAT (/' ***** INITIAL F(X) CANNOT BE COMPUTED *****') GO TO 390 CASE (13) WRITE (pu,340) 340 FORMAT (/' ***** BAD PARAMETERS TO ASSESS *****') GO TO 520 CASE (14) WRITE (pu,360) 360 FORMAT (/' ***** GRADIENT COULD NOT BE COMPUTED *****') IF (iv(niter) > 0) GO TO 460 GO TO 390 CASE (15) GO TO 500 END SELECT 370 WRITE (pu,380) iv(1) 380 FORMAT (/' ***** IV(1) =', i5, ' *****') GO TO 520 ! *** INITIAL CALL ON ITSUM *** 390 IF (iv(x0prt) /= 0) THEN IF (PRESENT(b)) THEN WRITE (pu,401) (i,x(i),b(1:2,i),d(i),i=1,p) 401 FORMAT (/' I Initial X(I) Lower Bound Upper Bound', & ' D(I)'// (' ', i5, 3g17.6, g14.3)) ELSE WRITE (pu,400) (i,x(i),d(i),i=1,p) 400 FORMAT (/' I Initial X(I) D(I)'// & (' ', i5, g17.6, g14.3)) END IF END IF ! *** THE FOLLOWING ARE TO AVOID UNDEFINED VARIABLES WHEN THE ! *** FUNCTION EVALUATION LIMIT IS 1... v(dstnrm) = zero v(fdif) = zero v(nreduc) = zero v(preduc) = zero v(reldx) = zero IF (iv1 >= 12) GO TO 520 iv(needhd) = 0 iv(prntit) = 0 IF (ol == 0) GO TO 520 IF (ol < 0 .AND. alg == 1) WRITE (pu,30) IF (ol < 0 .AND. alg == 2) WRITE (pu,40) IF (ol > 0 .AND. alg == 1) WRITE (pu,70) IF (ol > 0 .AND. alg == 2) WRITE (pu,80) IF (alg == 1) WRITE (pu,410) iv(nfcall), v(f) IF (alg == 2) WRITE (pu,420) iv(nfcall), v(f) 410 FORMAT (/' 0', i5, e10.3) 420 FORMAT (/' 0', i5, e11.3) GO TO 520 ! *** PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION *** 430 iv(needhd) = 1 IF (iv(statpr) <= 0) GO TO 460 oldf=MAX(ABS(v(f0)),ABS(v(f))) preldf=zero nreldf=zero IF (oldf <= zero) GO TO 440 preldf=v(preduc)/oldf nreldf=v(nreduc)/oldf 440 nf=iv(nfcall)-iv(nfcov) ng=iv(ngcall)-iv(ngcov) WRITE (pu,450) v(f), v(reldx), nf, ng, preldf, nreldf 450 FORMAT (/' FUNCTION', e17.6, ' RELDX', e17.3/ ' FUNC. EVALS', i8, & ' grad. evals', I8/ ' preldf', E16.3, ' npreldf', & E15.3) 460 IF (iv(solprt) == 0) GO TO 520 iv(needhd)=1 IF (iv(algsav) > 2) GO TO 520 WRITE (pu,470) 470 FORMAT (/' I FINAL X(I) D(I) G(I)'/) DO i=1,p WRITE (pu,490) i, x(i), d(i), g(i) END DO 490 FORMAT (' ', i5, e16.6, 2E14.3) GO TO 520 500 WRITE (pu,510) 510 FORMAT (/' INCONSISTENT DIMENSIONS') 520 RETURN ! *** LAST CARD OF ITSUM FOLLOWS *** END SUBROUTINE itsum SUBROUTINE v7cpy(p, y, x) ! *** SET Y = X, WHERE X AND Y ARE P-VECTORS *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: y(:) REAL (dp), INTENT(IN) :: x(:) y(1:p) = x(1:p) RETURN END SUBROUTINE v7cpy SUBROUTINE l7mst(d, g, ierr, ipivot, ka, p, qtr, r, step, v, w) ! *** COMPUTE LEVENBERG-MARQUARDT STEP USING MORE-HEBDEN TECHNIQUE ** ! *** NL2SOL VERSION 2.2. *** ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(IN) :: g(:) INTEGER, INTENT(IN OUT) :: ierr INTEGER, INTENT(IN OUT) :: ipivot(:) INTEGER, INTENT(IN OUT) :: ka INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: qtr(:) REAL (dp), INTENT(IN) :: r(:) REAL (dp), INTENT(OUT) :: step(:) REAL (dp), INTENT(IN OUT) :: v(:), w(:) ! DIMENSION W(P*(P+5)/2 + 4) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** PURPOSE *** ! GIVEN THE R MATRIX FROM THE QR DECOMPOSITION OF A JACOBIAN MATRIX, J, ! AS WELL AS Q-TRANSPOSE TIMES THE CORRESPONDING RESIDUAL VECTOR, RESID, ! THIS SUBROUTINE COMPUTES A LEVENBERG-MARQUARDT STEP OF APPROXIMATE ! LENGTH V(RADIUS) BY THE MORE-TECHNIQUE. ! *** PARAMETER DESCRIPTION *** ! D (IN) = THE SCALE VECTOR. ! G (IN) = THE GRADIENT VECTOR (J**T)*R. ! IERR (I/O) = RETURN CODE FROM QRFACT OR Q7RGS -- 0 MEANS R HAS ! FULL RANK. ! IPIVOT (I/O) = PERMUTATION ARRAY FROM QRFACT OR Q7RGS, WHICH COMPUTE ! QR DECOMPOSITIONS WITH COLUMN PIVOTING. ! KA (I/O). KA < 0 ON INPUT MEANS THIS IS THE FIRST CALL ON ! L7MST FOR THE CURRENT R AND QTR. ON OUTPUT KA CON- ! TAINS THE NUMBER OF HEBDEN ITERATIONS NEEDED TO DETERMINE ! STEP. KA = 0 MEANS A GAUSS-NEWTON STEP. ! P (IN) = NUMBER OF PARAMETERS. ! QTR (IN) = (Q**T)*RESID = Q-TRANSPOSE TIMES THE RESIDUAL VECTOR. ! R (IN) = THE R MATRIX, STORED COMPACTLY BY COLUMNS. ! STEP (OUT) = THE LEVENBERG-MARQUARDT STEP COMPUTED. ! V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW. ! W (I/O) = WORKSPACE OF LENGTH P*(P+5)/2 + 4. ! *** ENTRIES IN V *** ! V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G. ! V(DSTNRM) (I/O) = 2-NORM OF D*STEP. ! V(DST0) (I/O) = 2-NORM OF GAUSS-NEWTON STEP (FOR NONSING. J). ! V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED IN TWONORM(R)**2 MINUS ! TWONORM(R - J*STEP)**2. (SEE ALGORITHM NOTES BELOW.) ! V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP. ! V(NREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED ! FOR A GAUSS-NEWTON STEP. ! V(PHMNFC) (IN) = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP ! (MORE*S SIGMA). THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE ! BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS). ! V(PHMXFC) (IN) (SEE V(PHMNFC).) ! V(PREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED ! BY THE STEP RETURNED. ! V(RADIUS) (IN) = RADIUS OF CURRENT (SCALED) TRUST REGION. ! V(RAD0) (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL. ! V(STPPAR) (I/O) = MARQUARDT PARAMETER (OR ITS NEGATIVE IF THE SPECIAL ! CASE MENTIONED BELOW IN THE ALGORITHM NOTES OCCURS). ! NOTE -- SEE DATA STATEMENT BELOW FOR VALUES OF ABOVE SUBSCRIPTS. ! *** USAGE NOTES *** ! IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF ! V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT ! WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS). (THIS EXPLAINS ! WHY MANY PARAMETERS ARE LISTED AS I/O). ON AN INTIIAL CALL (ONE ! WITH KA = -1), THE CALLER NEED ONLY HAVE INITIALIZED D, G, KA, P, ! QTR, R, V(EPSLON), V(PHMNFC), V(PHMXFC), V(RADIUS), AND V(RAD0). ! *** APPLICATION AND USAGE RESTRICTIONS *** ! THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR LEAST- ! SQUARES) PACKAGE (REF. 1). ! *** ALGORITHM NOTES *** ! THIS CODE IMPLEMENTS THE STEP COMPUTATION SCHEME DESCRIBED IN ! REFS. 2 AND 4. FAST GIVENS TRANSFORMATIONS (SEE REF. 3, PP. 60- ! 62) ARE USED TO COMPUTE STEP WITH A NONZERO MARQUARDT PARAMETER. ! A SPECIAL CASE OCCURS IF J IS (NEARLY) SINGULAR AND V(RADIUS) ! IS SUFFICIENTLY LARGE. IN THIS CASE THE STEP RETURNED IS SUCH ! THAT TWONORM(R)**2 - TWONORM(R - J*STEP)**2 DIFFERS FROM ITS ! OPTIMAL VALUE BY LESS THAN V(EPSLON) TIMES THIS OPTIMAL VALUE, ! WHERE J AND R DENOTE THE ORIGINAL JACOBIAN AND RESIDUAL. (SEE ! REF. 2 FOR MORE DETAILS.) ! *** FUNCTIONS AND SUBROUTINES CALLED *** ! D7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. ! L7ITV - APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. ! L7IVM - APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. ! V7CPY - COPIES ONE VECTOR TO ANOTHER. ! V2NRM - RETURNS 2-NORM OF A VECTOR. ! *** REFERENCES *** ! 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE ! NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. ! SOFTWARE, VOL. 7, NO. 3. ! 2. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS, ! SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP. ! 186-197. ! 3. LAWSON, C.L., AND HANSON, R.J. (1974), SOLVING LEAST SQUARES ! PROBLEMS, PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J. ! 4. MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMENTATION AND ! THEORY, PP.105-116 OF SPRINGER LECTURE NOTES IN MATHEMATICS NO. ! 630, EDITED BY G.A. WATSON, SPRINGER-VERLAG, BERLIN & NEW YORK. ! *** GENERAL *** ! CODED BY DAVID M. GAY. ! THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH ! SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS ! MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND MCS-7906671. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** LOCAL VARIABLES *** INTEGER :: distsav, i, ip1, i1, j1, k, kalim, l, lk0, phipin, & pp1o2, res, res0, rmat, rmat0, uk0 REAL (dp) :: a, adi, alphak, b, dfacsq, dst, dtol, d1, d2, lk, oldphi, phi, & phimax, phimin, psifac, rad, si, sj, sqrtak, t, twopsi, uk, wl ! *** INTRINSIC FUNCTIONS *** ! REAL :: SQRT ! *** EXTERNAL FUNCTIONS AND SUBROUTINES *** ! REAL :: d7tpr, l7svn, r7mdc, v2nrm ! EXTERNAL d7tpr, l7itv, l7ivm, l7svn, r7mdc, v7cpy, v2nrm REAL (dp), PARAMETER :: dfactor=256._dp, eight=8._dp, half=0.5_dp, & negone=-1._dp, p001=0.001_dp, three=3._dp, ttol=2.5_dp REAL (dp), SAVE :: big=0.0_dp ! *** BODY *** ! *** FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK AND UK, ! *** THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR NONSING. J) ! *** AND THE VALUE RETURNED AS V(DSTNRM) ARE STORED AT W(LK0), ! *** W(UK0), W(PHIPIN), AND W(distsav) RESPECTIVELY. lk0 = p + 1 phipin = lk0 + 1 uk0 = phipin + 1 distsav = uk0 + 1 rmat0 = distsav ! *** A COPY OF THE R-MATRIX FROM THE QR DECOMPOSITION OF J IS STORED IN ! *** W STARTING AT W(RMAT), AND A COPY OF THE RESIDUAL VECTOR IS STORED ! *** IN W STARTING AT W(RES). THE LOOPS BELOW THAT UPDATE THE QR DECOMP. ! *** FOR A NONZERO MARQUARDT PARAMETER WORK ON THESE COPIES. ! N.B. The Fortran 77 version did not save the Givens scale factors which were ! stored in W(1:P). In certain circumstances, these values were ! overwritten in routine S7BQN. They have been relocated to scalef(). rmat = rmat0 + 1 pp1o2 = p * (p + 1) / 2 res0 = pp1o2 + rmat0 res = res0 + 1 rad = v(radius) IF (rad > zero) psifac = v(epslon)/((eight*(v(phmnfc) + one) + three) * rad**2) IF (big <= zero) big = r7mdc(6) phimax = v(phmxfc) * rad phimin = v(phmnfc) * rad ! *** DTOL, dfactor, AND DFACSQ ARE USED IN RESCALING THE FAST GIVENS ! *** REPRESENTATION OF THE UPDATED QR DECOMPOSITION. dtol = one/dfactor dfacsq = dfactor*dfactor ! *** OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY. IF ! *** WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT. oldphi = zero lk = zero uk = zero kalim = ka + 12 ! *** START OR RESTART, DEPENDING ON KA *** IF (.NOT. ALLOCATED(scalef)) THEN ALLOCATE( scalef(p) ) ELSE IF (p > SIZE(scalef)) THEN DEALLOCATE( scalef ) ALLOCATE( scalef(p) ) END IF IF (ka == 0) THEN GO TO 20 ELSE IF(ka > 0) THEN GO TO 370 END IF ! *** FRESH START -- COMPUTE V(NREDUC) *** ka = 0 kalim = 12 k = p IF (ierr /= 0) k = ABS(ierr) - 1 v(nreduc) = half* d7tpr(k, qtr, qtr) ! *** SET UP TO TRY INITIAL GAUSS-NEWTON STEP *** 20 v(dst0) = negone IF (ierr /= 0) GO TO 90 CALL l7svn(p, r, step, w(res:), t) IF (t >= one) GO TO 30 IF ( v2nrm(p, qtr) >= big*t) GO TO 90 ! *** COMPUTE GAUSS-NEWTON STEP *** ! *** NOTE -- THE R-MATRIX IS STORED COMPACTLY BY COLUMNS IN ! *** R(1), R(2), R(3), ... IT IS THE TRANSPOSE OF A ! *** LOWER TRIANGULAR MATRIX STORED COMPACTLY BY ROWS, AND WE ! *** TREAT IT AS SUCH WHEN USING L7ITV AND L7IVM. 30 CALL l7itv(p, w, r, qtr) ! *** TEMPORARILY STORE PERMUTED -D*STEP IN STEP. DO i = 1, p j1 = ipivot(i) step(i) = d(j1)*w(i) END DO dst = v2nrm(p, step) v(dst0) = dst phi = dst - rad IF (phi <= phimax) GO TO 410 ! *** IF THIS IS A RESTART, GO TO 110 *** IF (ka > 0) GO TO 110 ! *** GAUSS-NEWTON STEP WAS UNACCEPTABLE. COMPUTE L0 *** DO i = 1, p j1 = ipivot(i) step(i) = d(j1)*(step(i)/dst) END DO CALL l7ivm(p, step, r, step) t = one / v2nrm(p, step) w(phipin) = (t/rad)*t lk = phi*w(phipin) ! *** COMPUTE U0 *** 90 w(1:p) = g(1:p)/d(1:p) v(dgnorm) = v2nrm(p, w) uk = v(dgnorm)/rad IF (uk <= zero) GO TO 390 ! *** ALPHAK WILL BE USED AS THE CURRENT MARQUARDT PARAMETER. WE ! *** USE MORE*S SCHEME FOR INITIALIZING IT. alphak = ABS(v(stppar)) * v(rad0)/rad alphak = MIN(uk, MAX(alphak, lk)) ! *** TOP OF LOOP -- INCREMENT KA, COPY R TO RMAT, QTR TO RES *** 110 ka = ka + 1 CALL v7cpy(pp1o2, w(rmat:), r) CALL v7cpy(p, w(res:), qtr) ! *** SAFEGUARD ALPHAK AND INITIALIZE FAST GIVENS SCALE VECTOR. *** IF (alphak <= zero .OR. alphak < lk .OR. alphak >= uk) & alphak = uk * MAX(p001, SQRT(lk/uk)) IF (alphak <= zero) alphak = half * uk sqrtak = SQRT(alphak) scalef(1:p) = one ! *** ADD ALPHAK*D AND UPDATE QR DECOMP. USING FAST GIVENS TRANS. *** DO i = 1, p ! *** GENERATE, APPLY 1ST GIVENS TRANS. FOR ROW I OF ALPHAK*D. ! *** (USE STEP TO STORE TEMPORARY ROW) *** l = i*(i+1)/2 + rmat0 wl = w(l) d2 = one d1 = scalef(i) j1 = ipivot(i) adi = sqrtak*d(j1) IF (adi >= ABS(wl)) GO TO 150 130 a = adi/wl b = d2*a/d1 t = a*b + one IF (t > ttol) GO TO 150 scalef(i) = d1/t d2 = d2/t w(l) = t*wl a = -a DO j1 = i, p l = l + j1 step(j1) = a*w(l) END DO GO TO 170 150 b = wl/adi a = d1*b/d2 t = a*b + one IF (t > ttol) GO TO 130 scalef(i) = d2/t d2 = d1/t w(l) = t*adi DO j1 = i, p l = l + j1 wl = w(l) step(j1) = -wl w(l) = a*wl END DO 170 IF (i == p) EXIT ! *** NOW USE GIVENS TRANS. TO ZERO ELEMENTS OF TEMP. ROW *** ip1 = i + 1 DO i1 = ip1, p si = step(i1-1) IF (si == zero) CYCLE l = i1*(i1+1)/2 + rmat0 wl = w(l) d1 = scalef(i1) ! *** RESCALE ROW I1 IF NECESSARY *** IF (d1 >= dtol) GO TO 190 d1 = d1*dfacsq wl = wl/dfactor k = l DO j1 = i1, p k = k + j1 w(k) = w(k)/dfactor END DO ! *** USE GIVENS TRANS. TO ZERO NEXT ELEMENT OF TEMP. ROW 190 IF ( ABS(si) > ABS(wl)) GO TO 220 200 a = si/wl b = d2*a/d1 t = a*b + one IF (t > ttol) GO TO 220 w(l) = t*wl scalef(i1) = d1/t d2 = d2/t DO j1 = i1, p l = l + j1 wl = w(l) sj = step(j1) w(l) = wl + b*sj step(j1) = sj - a*wl END DO GO TO 240 220 b = wl/si a = d1*b/d2 t = a*b + one IF (t > ttol) GO TO 200 scalef(i1) = d2/t d2 = d1/t w(l) = t*si DO j1 = i1, p l = l + j1 wl = w(l) sj = step(j1) w(l) = a*wl + sj step(j1) = b*sj - wl END DO ! *** RESCALE TEMP. ROW IF NECESSARY *** 240 IF (d2 >= dtol) CYCLE d2 = d2*dfacsq step(i1:p) = step(i1:p)/dfactor END DO END DO ! *** COMPUTE STEP *** CALL l7itv(p, w(res:), w(rmat:), w(res:)) ! *** RECOVER STEP AND STORE PERMUTED -D*STEP AT W(RES) *** DO i = 1, p j1 = ipivot(i) k = res0 + i t = w(k) step(j1) = -t w(k) = t*d(j1) END DO dst = v2nrm(p, w(res:)) phi = dst - rad IF (phi <= phimax .AND. phi >= phimin) GO TO 430 IF (oldphi == phi) GO TO 430 oldphi = phi ! *** CHECK FOR (AND HANDLE) SPECIAL CASE *** IF (phi > zero) GO TO 310 IF (ka >= kalim) GO TO 430 twopsi = alphak*dst*dst - d7tpr(p, step, g) IF (alphak >= twopsi*psifac) GO TO 310 v(stppar) = -alphak GO TO 440 ! *** UNACCEPTABLE STEP -- UPDATE LK, UK, ALPHAK, AND TRY AGAIN *** 300 IF (phi < zero) uk = MIN(uk, alphak) GO TO 320 310 IF (phi < zero) uk = alphak 320 DO i = 1, p j1 = ipivot(i) k = res0 + i step(i) = d(j1) * (w(k)/dst) END DO CALL l7ivm(p, step, w(rmat:), step) step(1:p) = step(1:p) / SQRT(scalef(1:p)) t = one / v2nrm(p, step) alphak = alphak + t*phi*t/rad lk = MAX(lk, alphak) alphak = lk GO TO 110 ! *** RESTART *** 370 lk = w(lk0) uk = w(uk0) IF (v(dst0) > zero .AND. v(dst0) - rad <= phimax) GO TO 20 alphak = ABS(v(stppar)) dst = w(distsav) phi = dst - rad t = v(dgnorm)/rad IF (rad > v(rad0)) GO TO 380 ! *** SMALLER RADIUS *** uk = t IF (alphak <= zero) lk = zero IF (v(dst0) > zero) lk = MAX(lk, (v(dst0)-rad)*w(phipin)) GO TO 300 ! *** BIGGER RADIUS *** 380 IF (alphak <= zero .OR. uk > t) uk = t lk = zero IF (v(dst0) > zero) lk = MAX(lk, (v(dst0)-rad)*w(phipin)) GO TO 300 ! *** SPECIAL CASE -- RAD <= 0 OR (G = 0 AND J IS SINGULAR) *** 390 v(stppar) = zero dst = zero lk = zero uk = zero v(gtstep) = zero v(preduc) = zero step(1:p) = zero GO TO 450 ! *** ACCEPTABLE GAUSS-NEWTON STEP -- RECOVER STEP FROM W *** 410 alphak = zero DO i = 1, p j1 = ipivot(i) step(j1) = -w(i) END DO ! *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** 430 v(stppar) = alphak 440 v(gtstep) = MIN( d7tpr(p,step,g), zero) v(preduc) = half * (alphak*dst*dst - v(gtstep)) 450 v(dstnrm) = dst w(distsav) = dst w(lk0) = lk w(uk0) = uk v(rad0) = rad RETURN ! *** LAST CARD OF L7MST FOLLOWS *** END SUBROUTINE l7mst SUBROUTINE d7upd(d, dr, iv, n, nn, n2, p, v) ! *** UPDATE SCALE VECTOR D FOR NL2IT *** ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(IN OUT) :: d(:) REAL (dp), INTENT(IN) :: dr(:,:) ! dr(nd,p) INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: nn INTEGER, INTENT(IN) :: n2 INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: v(:) ! DIMENSION V(*) ! *** LOCAL VARIABLES *** INTEGER :: d0, i, jcn0, jcn1, jcni, jtol0, jtoli, k, sii REAL (dp) :: t, vdfac ! *** CONSTANTS *** REAL (dp), PARAMETER :: zero = 0.0_dp ! *** INTRINSIC FUNCTIONS *** ! ! REAL :: SQRT !------------------------------- BODY -------------------------------- IF (iv(dtype) /= 1 .AND. iv(niter) > 0) GO TO 999 jcn1 = iv(jcn) jcn0 = ABS(jcn1) - 1 IF (jcn1 < 0) GO TO 10 iv(jcn) = -jcn1 v(jcn1:jcn1+p-1) = zero 10 DO i = 1, p jcni = jcn0 + i t = v(jcni) DO k = 1, nn t = MAX(t, ABS(dr(k,i))) END DO v(jcni) = t END DO IF (n2 < n) GO TO 999 vdfac = v(dfac) jtol0 = iv(jtol) - 1 d0 = jtol0 + p sii = iv(s) - 1 DO i = 1, p sii = sii + i jcni = jcn0 + i t = v(jcni) IF (v(sii) > zero) t = MAX( SQRT(v(sii)), t) jtoli = jtol0 + i d0 = d0 + 1 IF (t < v(jtoli)) t = MAX(v(d0), v(jtoli)) d(i) = MAX(vdfac*d(i), t) END DO 999 RETURN ! *** LAST CARD OF D7UPD FOLLOWS *** END SUBROUTINE d7upd SUBROUTINE l7nvr(n, lin, l) ! *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED *** ! *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. *** ! *** PARAMETERS *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: lin(:) REAL (dp), INTENT(IN) :: l(:) ! DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2) ! *** LOCAL VARIABLES *** INTEGER :: i, ii, im1, jj, j0, j1, k, k0, np1 REAL (dp) :: t ! *** BODY *** np1 = n + 1 j0 = n*(np1)/2 DO ii = 1, n i = np1 - ii lin(j0) = one/l(j0) IF (i <= 1) GO TO 999 j1 = j0 im1 = i - 1 DO jj = 1, im1 t = zero j0 = j1 k0 = j1 - jj DO k = 1, jj t = t - l(k0)*lin(j0) j0 = j0 - 1 k0 = k0 + k - i END DO lin(j0) = t/l(k0) END DO j0 = j0 - 1 END DO 999 RETURN ! *** LAST CARD OF L7NVR FOLLOWS *** END SUBROUTINE l7nvr SUBROUTINE l7sqr(n, a, l) ! *** COMPUTE A = LOWER TRIANGLE OF L*(L**T), WITH BOTH ! *** L AND A STORED COMPACTLY BY ROWS. (BOTH MAY OCCUPY THE ! *** SAME STORAGE. ! *** PARAMETERS *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: a(:) REAL (dp), INTENT(IN) :: l(:) ! DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) ! *** LOCAL VARIABLES *** INTEGER :: i, ii, ij, ik, ip1, i0, j, jj, jk, j0, k, np1 REAL (dp) :: t np1 = n + 1 i0 = n*(n+1)/2 DO ii = 1, n i = np1 - ii ip1 = i + 1 i0 = i0 - i j0 = i*(i+1)/2 DO jj = 1, i j = ip1 - jj j0 = j0 - j t = 0.0_dp DO k = 1, j ik = i0 + k jk = j0 + k t = t + l(ik)*l(jk) END DO ij = i0 + j a(ij) = t END DO END DO RETURN END SUBROUTINE l7sqr SUBROUTINE l7svx(p, l, x, y, largest) ! *** ESTIMATE LARGEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L ! *** PARAMETER DECLARATIONS *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: l(:) REAL (dp), INTENT(OUT) :: x(:), y(:), largest ! DIMENSION L(P*(P+1)/2) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** PURPOSE *** ! THIS FUNCTION RETURNS A GOOD UNDER-ESTIMATE OF THE LARGEST ! SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L. ! *** PARAMETER DESCRIPTION *** ! P (IN) = THE ORDER OF L. L IS A P X P LOWER TRIANGULAR MATRIX. ! L (IN) = ARRAY HOLDING THE ELEMENTS OF L IN ROW ORDER, I.E. ! L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC. ! X (OUT) IF largest RETURNS A POSITIVE VALUE, THEN X = (L**T)*Y IS AN ! (UNNORMALIZED) APPROXIMATE RIGHT SINGULAR VECTOR ! CORRESPONDING TO THE LARGEST SINGULAR VALUE. THIS ! APPROXIMATION MAY BE CRUDE. ! Y (OUT) IF largest RETURNS A POSITIVE VALUE, THEN Y = L*X IS A NORMALIZED ! APPROXIMATE LEFT SINGULAR VECTOR CORRESPONDING TO THE LARGEST ! SINGULAR VALUE. THIS APPROXIMATION MAY BE VERY CRUDE. ! THE CALLER MAY PASS THE SAME VECTOR FOR X AND Y (NONSTANDARD ! FORTRAN USAGE), IN WHICH CASE X OVER-WRITES Y. ! *** ALGORITHM NOTES *** ! THE ALGORITHM IS BASED ON ANALOGY WITH (1). IT USES A ! RANDOM NUMBER GENERATOR PROPOSED IN (4), WHICH PASSES THE ! SPECTRAL TEST WITH FLYING COLORS -- SEE (2) AND (3). ! *** SUBROUTINES AND FUNCTIONS CALLED *** ! V2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR. ! *** REFERENCES *** ! (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977), ! AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT ! TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY. ! (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL ! RANDOM-NUMBER GENERATORS -- AN EMPIRICAL VIEW, ! MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. ! (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 ! (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. ! (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER GENERATORS ! WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, PP. 586-593. ! *** HISTORY *** ! DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978). ! *** GENERAL *** ! THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH ! SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS ! MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! *** LOCAL VARIABLES *** INTEGER :: i, ix, j, ji, jj, jjj, jm1, j0, pm1, pplus1 REAL (dp) :: b, blji, sminus, splus, t, yi ! *** EXTERNAL FUNCTIONS AND SUBROUTINES *** ! REAL :: d7tpr, v2nrm ! EXTERNAL d7tpr, v2nrm, v2axy REAL (dp), PARAMETER :: half=0.5_dp, r9973=9973._dp ! *** BODY *** ix = 2 pplus1 = p + 1 pm1 = p - 1 ! *** FIRST INITIALIZE X TO PARTIAL SUMS *** j0 = p*pm1/2 jj = j0 + p ix = MOD(3432*ix, 9973) b = half*(one + REAL(ix)/r9973) x(p) = b * l(jj) IF (p <= 1) GO TO 40 DO i = 1, pm1 ji = j0 + i x(i) = b * l(ji) END DO ! *** COMPUTE X = (L**T)*B, WHERE THE COMPONENTS OF B HAVE RANDOMLY ! *** CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE. ! DO J = P-1 TO 1 BY -1... DO j = p-1, 1, -1 ! *** DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J ! *** THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I. ix = MOD(3432*ix, 9973) b = half*(one + REAL(ix)/r9973) jm1 = j - 1 j0 = j*jm1/2 splus = zero sminus = zero DO i = 1, j ji = j0 + i blji = b * l(ji) splus = splus + ABS(blji + x(i)) sminus = sminus + ABS(blji - x(i)) END DO IF (sminus > splus) b = -b x(j) = zero ! *** UPDATE PARTIAL SUMS *** CALL v2axy(j, x, b, l(j0+1:), x) END DO ! *** NORMALIZE X *** 40 t = v2nrm(p, x) IF (t <= zero) GO TO 80 t = one / t x(1:p) = t*x(1:p) ! *** COMPUTE L*X = Y AND RETURN SVMAX = TWONORM(Y) *** DO jjj = 1, p j = pplus1 - jjj ji = j*(j-1)/2 + 1 y(j) = d7tpr(j, l(ji:), x) END DO ! *** NORMALIZE Y AND SET X = (L**T)*Y *** t = one / v2nrm(p, y) ji = 1 DO i = 1, p yi = t * y(i) x(i) = zero CALL v2axy(i, x, yi, l(ji:), x) ji = ji + i END DO largest = v2nrm(p, x) GO TO 999 80 largest = zero 999 RETURN ! *** LAST CARD OF L7SVX FOLLOWS *** END SUBROUTINE l7svx SUBROUTINE l7tsq(n, a, l) ! *** SET A TO LOWER TRIANGLE OF (L**T) * L *** ! *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. *** ! *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: a(:) REAL (dp), INTENT(IN) :: l(:) ! DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) INTEGER :: i, ii, iim1, i1, j, k, m REAL (dp) :: lii, lj ii = 0 DO i = 1, n i1 = ii + 1 ii = ii + i m = 1 IF (i == 1) GO TO 30 iim1 = ii - 1 DO j = i1, iim1 lj = l(j) DO k = i1, j a(m) = a(m) + lj*l(k) m = m + 1 END DO END DO 30 lii = l(ii) DO j = i1, ii a(j) = lii * l(j) END DO END DO RETURN ! *** LAST CARD OF L7TSQ FOLLOWS *** END SUBROUTINE l7tsq SUBROUTINE l7tvm(n, x, l, y) ! *** COMPUTE X = (L**T)*Y, WHERE L IS AN N X N LOWER ! *** TRIANGULAR MATRIX STORED COMPACTLY BY ROWS. INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: l(:) REAL (dp), INTENT(IN) :: y(:) ! DIMENSION L(N*(N+1)/2) INTEGER :: i, ij, i0, j REAL (dp) :: yi i0 = 0 DO i = 1, n yi = y(i) x(i) = zero DO j = 1, i ij = i0 + j x(j) = x(j) + yi*l(ij) END DO i0 = i0 + i END DO RETURN ! *** LAST CARD OF L7TVM FOLLOWS *** END SUBROUTINE l7tvm SUBROUTINE l7vml(n, x, l, y) ! *** COMPUTE X = L*Y, WHERE L IS AN N X N LOWER TRIANGULAR ! *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME ! *** STORAGE. *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: l(:) REAL (dp), INTENT(IN) :: y(:) ! DIMENSION L(N*(N+1)/2) INTEGER :: i, ii, ij, i0, j, np1 REAL (dp) :: t np1 = n + 1 i0 = n*(n+1)/2 DO ii = 1, n i = np1 - ii i0 = i0 - i t = zero DO j = 1, i ij = i0 + j t = t + l(ij)*y(j) END DO x(i) = t END DO RETURN ! *** LAST CARD OF L7VML FOLLOWS *** END SUBROUTINE l7vml SUBROUTINE q7apl(n, p, j, r, ierr) ! *****PARAMETERS. INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: j(:,:) ! j(nn,p) REAL (dp), INTENT(IN OUT) :: r(:) INTEGER, INTENT(IN) :: ierr ! .................................................................. ! *****PURPOSE. ! THIS SUBROUTINE APPLIES TO R THE ORTHOGONAL TRANSFORMATIONS ! STORED IN J BY QRFACT ! *****PARAMETER DESCRIPTION. ! ON INPUT. ! N IS THE NUMBER OF ROWS OF J AND THE SIZE OF THE VECTOR R ! P IS THE NUMBER OF COLUMNS OF J AND THE SIZE OF SIGMA ! J CONTAINS ON AND BELOW ITS DIAGONAL THE COLUMN VECTORS ! U WHICH DETERMINE THE HOUSEHOLDER TRANSFORMATIONS ! IDENT - U*U.TRANSPOSE ! R IS THE RIGHT HAND SIDE VECTOR TO WHICH THE ORTHOGONAL ! TRANSFORMATIONS WILL BE APPLIED ! IERR IF NON-ZERO INDICATES THAT NOT ALL THE TRANSFORMATIONS ! WERE SUCCESSFULLY DETERMINED AND ONLY THE FIRST ! ABS(IERR) - 1 TRANSFORMATIONS WILL BE USED ! ON OUTPUT. ! R HAS BEEN OVERWRITTEN BY ITS TRANSFORMED IMAGE ! *****APPLICATION AND USAGE RESTRICTIONS. ! NONE ! *****ALGORITHM NOTES. ! THE VECTORS U WHICH DETERMINE THE HOUSEHOLDER TRANSFORMATIONS ! ARE NORMALIZED SO THAT THEIR 2-NORM SQUARED IS 2. THE USE OF ! THESE TRANSFORMATIONS HERE IS IN THE SPIRIT OF (1). ! *****SUBROUTINES AND FUNCTIONS CALLED. ! D7TPR - FUNCTION, RETURNS THE INNER PRODUCT OF VECTORS ! *****REFERENCES. ! (1) BUSINGER, P. A., AND GOLUB, G. H. (1965), LINEAR LEAST SQUARES ! SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, NUMER. MATH. 7, PP. 269-276. ! *****HISTORY. ! DESIGNED BY DAVID M. GAY, CODED BY STEPHEN C. PETERS (WINTER 1977) ! CALL ON V2AXY SUBSTITUTED FOR DO LOOP, FALL 1983. ! *****GENERAL. ! THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH ! SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS ! MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. ! .................................................................. ! *****LOCAL VARIABLES. INTEGER :: k, l, nl1 ! *****FUNCTIONS. ! REAL :: d7tpr ! EXTERNAL d7tpr, v2axy ! *** BODY *** k = p IF (ierr /= 0) k = ABS(ierr) - 1 IF ( k == 0) GO TO 999 DO l = 1, k nl1 = n - l + 1 CALL v2axy(nl1, r(l:), - d7tpr(nl1, j(l:,l), r(l:)), j(l:,l), r(l:)) END DO 999 RETURN ! *** LAST LINE OF Q7APL FOLLOWS *** END SUBROUTINE q7apl SUBROUTINE q7rad(n, p, qtr, qtrset, rmat, w, y) ! *** ADD ROWS W TO QR FACTORIZATION WITH R MATRIX RMAT AND ! *** Q**T * RESIDUAL = QTR. Y = NEW COMPONENTS OF RESIDUAL ! *** CORRESPONDING TO W. QTR, Y REFERENCED ONLY IF QTRSET = .TRUE. INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: qtr(:) LOGICAL, INTENT(IN) :: qtrset REAL (dp), INTENT(IN OUT) :: rmat(:) REAL (dp), INTENT(IN OUT) :: w(:,:) ! w(nn,p) REAL (dp), INTENT(IN OUT) :: y(:) ! DIMENSION RMAT(P*(P+1)/2) ! REAL :: SQRT ! REAL :: d7tpr, r7mdc, v2nrm ! EXTERNAL d7tpr, r7mdc, v2axy, v7scl, v2nrm ! *** LOCAL VARIABLES *** INTEGER :: i, ii, ij, ip1, j, k, nk REAL (dp) :: ari, qri, ri, s, t, wi REAL (dp), PARAMETER :: one = 1._dp, zero = 0._dp REAL (dp), SAVE :: big = -1._dp, bigrt = -1._dp, tiny = 0._dp, & tinyrt = 0._dp !------------------------------ BODY ----------------------------------- IF (tiny > zero) GO TO 10 tiny = r7mdc(1) big = r7mdc(6) IF (tiny*big < one) tiny = one / big 10 k = 1 nk = n ii = 0 DO i = 1, p ii = ii + i ip1 = i + 1 ij = ii + i IF (nk <= 1) t = ABS(w(k,i)) IF (nk > 1) t = v2nrm(nk, w(k:,i)) IF (t < tiny) CYCLE ri = rmat(ii) IF (ri /= zero) GO TO 100 IF (nk > 1) GO TO 30 ij = ii DO j = i, p rmat(ij) = w(k,j) ij = ij + j END DO IF (qtrset) qtr(i) = y(k) w(k,i) = zero GO TO 999 30 wi = w(k,i) IF (bigrt > zero) GO TO 40 bigrt = r7mdc(5) tinyrt = r7mdc(2) 40 IF (t <= tinyrt) GO TO 50 IF (t >= bigrt) GO TO 50 IF (wi < zero) t = -t wi = wi + t s = SQRT(t * wi) GO TO 70 50 s = SQRT(t) IF (wi < zero) GO TO 60 wi = wi + t s = s * SQRT(wi) GO TO 70 60 t = -t wi = wi + t s = s * SQRT(-wi) 70 w(k,i) = wi w(k:k+nk-1, i) = ( one/s ) * w(k:k+nk-1, i) rmat(ii) = -t IF (.NOT. qtrset) GO TO 80 CALL v2axy(nk, y(k:), - d7tpr(nk, y(k:), w(k:,i)), w(k:,i), y(k:)) qtr(i) = y(k) 80 IF (ip1 > p) GO TO 999 DO j = ip1, p CALL v2axy(nk, w(k:,j), - d7tpr(nk, w(k:,j), w(k:,i)), w(k:,i), w(k:,j)) rmat(ij) = w(k,j) ij = ij + j END DO IF (nk <= 1) GO TO 999 k = k + 1 nk = nk - 1 CYCLE 100 ari = ABS(ri) IF (ari > t) GO TO 110 t = t * SQRT(one + (ari/t)**2) GO TO 120 110 t = ari * SQRT(one + (t/ari)**2) 120 IF (ri < zero) t = -t ri = ri + t rmat(ii) = -t s = -ri / t IF (nk <= 1) GO TO 150 w(k:k+nk-1, i) = ( one/ri ) * w(k:k+nk-1, i) IF (.NOT. qtrset) GO TO 130 qri = qtr(i) t = s * ( qri + d7tpr(nk, y(k:), w(k:,i)) ) qtr(i) = qri + t 130 IF (ip1 > p) GO TO 999 IF (qtrset) CALL v2axy(nk, y(k:), t, w(k:,i), y(k:)) DO j = ip1, p ri = rmat(ij) t = s * ( ri + d7tpr(nk, w(k:,j), w(k:,i)) ) CALL v2axy(nk, w(k:,j), t, w(k:,i), w(k:,j)) rmat(ij) = ri + t ij = ij + j END DO CYCLE 150 wi = w(k,i) / ri w(k,i) = wi IF (.NOT. qtrset) GO TO 160 qri = qtr(i) t = s * ( qri + y(k)*wi ) qtr(i) = qri + t 160 IF (ip1 > p) GO TO 999 IF (qtrset) y(k) = t*wi + y(k) DO j = ip1, p ri = rmat(ij) t = s * (ri + w(k,j)*wi) w(k,j) = w(k,j) + t*wi rmat(ij) = ri + t ij = ij + j END DO END DO 999 RETURN ! *** LAST LINE OF Q7RAD FOLLOWS *** END SUBROUTINE q7rad SUBROUTINE r7tvm(n, p, y, d, u, x) ! *** SET Y TO R*X, WHERE R IS THE UPPER TRIANGULAR MATRIX WHOSE ! *** DIAGONAL IS IN D AND WHOSE STRICT UPPER TRIANGLE IS IN U. INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: y(:) REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(IN) :: u(:,:) ! u(n,p) REAL (dp), INTENT(IN) :: x(:) ! REAL :: d7tpr ! EXTERNAL d7tpr ! *** LOCAL VARIABLES *** INTEGER :: i, ii, pl, pp1 REAL (dp) :: t ! *** BODY *** pl = MIN(n-1, p) pp1 = pl + 1 DO ii = 1, pl i = pp1 - ii t = x(i) * d(i) IF (i > 1) t = t + d7tpr(i-1, u(:,i), x) y(i) = t END DO RETURN ! *** LAST LINE OF R7TVM FOLLOWS *** END SUBROUTINE r7tvm FUNCTION rldst(p, d, x, x0) RESULT(fn_val) ! *** COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0 *** ! *** NL2SOL VERSION 2.2 *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: d(:), x(:), x0(:) REAL (dp) :: fn_val INTEGER :: i REAL (dp) :: emax, t, xmax ! *** BODY *** emax = zero xmax = zero DO i = 1, p t = ABS(d(i) * (x(i) - x0(i))) IF (emax < t) emax = t t = d(i) * ( ABS(x(i)) + ABS(x0(i))) IF (xmax < t) xmax = t END DO fn_val = zero IF (xmax > zero) fn_val = emax / xmax RETURN ! *** LAST CARD OF RLDST FOLLOWS *** END FUNCTION rldst SUBROUTINE s7lup(a, cosmin, p, size, step, u, w, wchmtd, wscale, y) ! *** UPDATE SYMMETRIC A SO THAT A * STEP = Y *** ! *** (LOWER TRIANGLE OF A STORED ROWWISE *** ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(OUT) :: a(:) REAL (dp), INTENT(IN) :: cosmin INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: size REAL (dp), INTENT(IN) :: step(:) REAL (dp), INTENT(OUT) :: u(:) REAL (dp), INTENT(OUT) :: w(:) REAL (dp), INTENT(IN) :: wchmtd(:) REAL (dp), INTENT(OUT) :: wscale REAL (dp), INTENT(IN) :: y(:) ! DIMENSION A(P*(P+1)/2) ! *** LOCAL VARIABLES *** INTEGER :: i, j, k REAL (dp) :: denmin, sdotwm, t, ui, wi ! *** EXTERNAL FUNCTIONS AND SUBROUTINES *** ! REAL :: d7tpr, v2nrm ! EXTERNAL d7tpr, s7lvm, v2nrm REAL (dp), PARAMETER :: half=0.5_dp !----------------------------------------------------------------------- sdotwm = d7tpr(p, step, wchmtd) denmin = cosmin * v2nrm(p,step) * v2nrm(p, wchmtd) wscale = one IF (denmin /= zero) wscale = MIN(one, ABS(sdotwm/denmin)) t = zero IF (sdotwm /= zero) t = wscale / sdotwm w(1:p) = t * wchmtd(1:p) CALL s7lvm(p, u, a, step) t = half * (size * d7tpr(p, step, u) - d7tpr(p, step, y)) DO i = 1, p u(i) = t*w(i) + y(i) - size*u(i) END DO ! *** SET A = A + U*(W**T) + W*(U**T) *** k = 1 DO i = 1, p ui = u(i) wi = w(i) DO j = 1, i a(k) = size*a(k) + ui*w(j) + wi*u(j) k = k + 1 END DO END DO RETURN ! *** LAST CARD OF S7LUP FOLLOWS *** END SUBROUTINE s7lup SUBROUTINE s7lvm(p, y, s, x) ! *** SET Y = S * X, S = P X P SYMMETRIC MATRIX. *** ! *** LOWER TRIANGLE OF S STORED ROWWISE. *** ! *** PARAMETER DECLARATIONS *** INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: y(:) REAL (dp), INTENT(IN) :: s(:) REAL (dp), INTENT(IN) :: x(:) ! DIMENSION S(P*(P+1)/2) ! *** LOCAL VARIABLES *** INTEGER :: i, im1, j, k REAL (dp) :: xi ! *** EXTERNAL FUNCTION *** ! REAL :: d7tpr ! EXTERNAL d7tpr !----------------------------------------------------------------------- j = 1 DO i = 1, p y(i) = d7tpr(i, s(j:), x) j = j + i END DO IF (p <= 1) GO TO 999 j = 1 DO i = 2, p xi = x(i) im1 = i - 1 j = j + 1 DO k = 1, im1 y(k) = y(k) + s(j)*xi j = j + 1 END DO END DO 999 RETURN ! *** LAST CARD OF S7LVM FOLLOWS *** END SUBROUTINE s7lvm SUBROUTINE parck(alg, d, iv, liv, lv, n, v) ! *** CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES *** ! *** ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT. INTEGER, INTENT(IN) :: alg REAL (dp), INTENT(IN) :: d(:) INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: liv INTEGER, INTENT(IN) :: lv INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: v(:) ! REAL :: r7mdc ! EXTERNAL ivset, r7mdc, v7cpy, v7dfl ! IVSET -- SUPPLIES DEFAULT VALUES TO BOTH IV AND V. ! R7MDC -- RETURNS MACHINE-DEPENDENT CONSTANTS. ! V7CPY -- COPIES ONE VECTOR TO ANOTHER. ! V7DFL -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE. ! *** LOCAL VARIABLES *** INTEGER :: alg1, i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu CHARACTER (LEN=1), PARAMETER :: varnm(2) = (/ 'P', 'P' /), & sh(2) = (/ 'S', 'H' /) CHARACTER (LEN=12) :: which REAL (dp) :: vk REAL (dp), SAVE :: big=0.0_dp, machep=-1.0_dp, tiny=1.0_dp CHARACTER (LEN=6), PARAMETER :: vn(34) = (/ 'EPSLON', 'PHMNFC', 'PHMXFC', & 'DECFAC', 'INCFAC', 'RDFCMN', 'RDFCMX', & 'TUNER1', 'TUNER2', 'TUNER3', 'TUNER4', & 'TUNER5', 'AFCTOL', 'RFCTOL', 'XCTOL.', & 'XFTOL.', 'LMAX0.', 'LMAXS.', 'SCTOL.', & 'DINIT.', 'DTINIT', 'D0INIT', 'DFAC..', & 'DLTFDC', 'DLTFDJ', 'DELTA0', 'FUZZ..', & 'RLIMIT', 'COSMIN', 'HUBERC', 'RSPTOL', & 'SIGMIN', 'ETA0..', 'BIAS..' /) REAL (dp), SAVE :: vm(34) = (/ 0.001_dp, -0.99_dp, 0.001_dp, 0.01_dp, & 1.2_dp, 0.01_dp, 1.2_dp, 0._dp, 0._dp, & 0.001_dp, -1._dp, 0._dp, 0._dp, 0._dp, 0._dp, & 0._dp, 0._dp, 0._dp, 0._dp, -10._dp, 0._dp, & 0._dp, 0._dp, 0._dp, 0._dp, 0._dp, & 1.01_dp, 1.e+10_dp, 0._dp, 0._dp, 0._dp, 0._dp, & 0._dp, 0._dp /) REAL (dp), SAVE :: vx(34) = (/ 0.9_dp, -0.001_dp, 10._dp, 0.8_dp, 100._dp, & 0.8_dp, 100._dp, 0.5_dp, 0.5_dp, 1._dp, & 1._dp, 1._dp, 1._dp, 0.1_dp, 1._dp, 1._dp, & 1._dp, 1._dp, 1._dp, 1._dp, 1._dp, 1._dp, & 1._dp, 1._dp, 1._dp, 1._dp, 1.e+10_dp, & 1._dp, 1._dp, 1._dp, 1._dp, 1._dp, 1._dp, & 1._dp /) CHARACTER (LEN=12), PARAMETER :: cngd = '---CHANGED V', dflt = 'NONDEFAULT V' INTEGER, PARAMETER :: ijmp = 33, jlim(4) = (/ 0, 24, 0, 24 /) INTEGER, PARAMETER :: ndflt(4) = (/ 32, 25, 32, 25 /) INTEGER, PARAMETER :: miniv(4) = (/ 82, 59, 103, 103 /) !............................... BODY ................................ pu = 0 IF (prunit <= liv) pu = iv(prunit) IF (algsav > liv) GO TO 20 IF (alg == iv(algsav)) GO TO 20 IF (pu /= 0) WRITE(pu,10) alg, iv(algsav) 10 FORMAT(/' The first PARAMETER TO ivset should be', i3, ' rather than', i3) iv(1) = 67 GO TO 999 20 IF (alg < 1 .OR. alg > 4) GO TO 340 miv1 = miniv(alg) IF (iv(1) == 15) GO TO 360 alg1 = MOD(alg-1,2) + 1 IF (iv(1) == 0) CALL ivset(alg, iv, liv, lv, v) iv1 = iv(1) IF (iv1 /= 13 .AND. iv1 /= 12) GO TO 30 IF (perm <= liv) miv1 = MAX(miv1, iv(perm) - 1) IF (ivneed <= liv) miv2 = miv1 + MAX(iv(ivneed), 0) IF (lastiv <= liv) iv(lastiv) = miv2 IF (liv < miv1) GO TO 300 iv(ivneed) = 0 iv(lastv) = MAX(iv(vneed), 0) + iv(lmat) - 1 iv(vneed) = 0 IF (liv < miv2) GO TO 300 IF (lv < iv(lastv)) GO TO 320 30 IF (iv1 < 12 .OR. iv1 > 14) GO TO 60 IF (n >= 1) GO TO 50 iv(1) = 81 IF (pu == 0) GO TO 999 WRITE(pu,40) varnm(alg1), n 40 FORMAT(/' /// bad', a1, ' =', i5) GO TO 999 50 IF (iv1 /= 14) iv(nextiv) = iv(perm) IF (iv1 /= 14) iv(nextv) = iv(lmat) IF (iv1 == 13) GO TO 999 k = iv(parsav) - epslon CALL v7dfl(alg1, v(k+1:)) iv(dtype0) = 2 - alg1 iv(oldn) = n which = dflt GO TO 110 60 IF (n == iv(oldn)) GO TO 80 iv(1) = 17 IF (pu == 0) GO TO 999 WRITE(pu,70) varnm(alg1), iv(oldn), n 70 FORMAT(/' /// ', A, ' changed from ', i5, ' TO ', i5) GO TO 999 80 IF (iv1 <= 11 .AND. iv1 >= 1) GO TO 100 iv(1) = 80 IF (pu /= 0) WRITE(pu,90) iv1 90 FORMAT(/' /// iv(1) =', i5, ' should be between 0 AND 14.') GO TO 999 100 which = cngd 110 IF (iv1 == 14) iv1 = 12 IF (big > tiny) GO TO 120 tiny = r7mdc(1) machep = r7mdc(3) big = r7mdc(6) vm(12) = machep vx(12) = big vx(13) = big vm(14) = machep vm(17) = tiny vx(17) = big vm(18) = tiny vx(18) = big vx(20) = big vx(21) = big vx(22) = big vm(24) = machep vm(25) = machep vm(26) = machep vx(28) = r7mdc(5) vm(29) = machep vx(30) = big vm(33) = machep 120 m = 0 i = 1 j = jlim(alg1) k = epslon ndfalt = ndflt(alg1) DO l = 1, ndfalt vk = v(k) WRITE(*, '(a, i4, g13.5)') ' k, v(k) =', k, v(k) IF (vk >= vm(i) .AND. vk <= vx(i)) GO TO 140 m = k IF (pu /= 0) WRITE(pu,130) vn(i), k, vk, vm(i), vx(i) 130 FORMAT(/' /// ', A6, '.. v(', i2, ') =', e11.3, ' should be between', & e11.3, ' AND', e11.3) 140 k = k + 1 i = i + 1 IF (i == j) i = ijmp END DO IF (iv(nvdflt) == ndfalt) GO TO 170 iv(1) = 51 IF (pu == 0) GO TO 999 WRITE(pu,160) iv(nvdflt), ndfalt 160 FORMAT(/' iv(nvdflt) =', i5, ' rather than ', i5) GO TO 999 170 IF ((iv(dtype) > 0 .OR. v(dinit) > zero) .AND. iv1 == 12) GO TO 200 DO i = 1, n IF (d(i) > zero) CYCLE m = 18 IF (pu /= 0) WRITE(pu,180) i, d(i) 180 FORMAT(/' /// d(', i3, ') =', e11.3, ' should be positive') END DO 200 IF (m == 0) GO TO 210 iv(1) = m GO TO 999 210 IF (pu == 0 .OR. iv(parprt) == 0) GO TO 999 IF (iv1 /= 12 .OR. iv(inits) == alg1-1) GO TO 230 m = 1 WRITE(pu,220) sh(alg1), iv(inits) 220 FORMAT(/' Nondefault values....'/ ' init', a1, '..... iv(25) =', i3) 230 IF (iv(dtype) == iv(dtype0)) GO TO 250 IF (m == 0) WRITE(pu,260) which m = 1 WRITE(pu,240) iv(dtype) 240 FORMAT(' dtype..... iv(16) =', i3) 250 i = 1 j = jlim(alg1) k = epslon l = iv(parsav) ndfalt = ndflt(alg1) DO ii = 1, ndfalt IF (v(k) == v(l)) GO TO 280 IF (m == 0) WRITE(pu,260) which 260 FORMAT(/' ' , 3A4, 'VALUES....'/) m = 1 WRITE(pu,270) vn(i), k, v(k) 270 FORMAT(' ', A6, '.. v(', i2, ') =', e15.7) 280 k = k + 1 l = l + 1 i = i + 1 IF (i == j) i = ijmp END DO iv(dtype0) = iv(dtype) parsv1 = iv(parsav) CALL v7cpy(iv(nvdflt), v(parsv1:), v(epslon:)) GO TO 999 300 iv(1) = 15 IF (pu == 0) GO TO 999 WRITE(pu,310) liv, miv2 310 FORMAT(/' /// liv =', i5, ' must be at least', i5) IF (liv < miv1) GO TO 999 IF (lv < iv(lastv)) GO TO 320 GO TO 999 320 iv(1) = 16 IF (pu /= 0) WRITE(pu,330) lv, iv(lastv) 330 FORMAT(/' /// lv =', i5, ' must be at least', i5) GO TO 999 340 iv(1) = 67 IF (pu /= 0) WRITE(pu,350) alg 350 FORMAT(/' /// alg =', i5, ' must be 1 2, 3, OR 4') GO TO 999 360 IF (pu /= 0) WRITE(pu,370) liv, miv1 370 FORMAT(/' /// liv =', i5, ' must be at least', i5, & ' to compute true MIN. liv AND MIN. lv') IF (lastiv <= liv) iv(lastiv) = miv1 IF (lastv <= liv) iv(lastv) = 0 999 RETURN ! *** LAST LINE OF PARCK FOLLOWS *** END SUBROUTINE parck SUBROUTINE n2cvp(iv, p, v) ! *** PRINT COVARIANCE MATRIX FOR RN2G *** INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: v(:) ! *** LOCAL VARIABLES *** INTEGER :: cov1, i, ii, i1, pu REAL (dp) :: t ! *** BODY *** IF (iv(1) > 8) GO TO 999 pu = iv(prunit) IF (pu == 0) GO TO 999 IF (iv(statpr) == 0) GO TO 30 IF (iv(nfcov) > 0) WRITE(pu,10) iv(nfcov) 10 FORMAT(/' ', i4, ' Extra func. evals for covariance AND diagnostics.') IF (iv(ngcov) > 0) WRITE(pu,20) iv(ngcov) 20 FORMAT(' ', i4, ' Extra grad. evals for covariance AND diagnostics.') 30 IF (iv(covprt) <= 0) GO TO 999 cov1 = iv(covmat) IF (iv(regd) <= 0 .AND. cov1 <= 0) GO TO 70 iv(needhd) = 1 t = v(rcond)**2 IF (ABS(iv(covreq)) > 2) GO TO 50 WRITE(pu,40) t 40 FORMAT(/' Reciprocal condition of f.d. hessian = at most', e10.2) GO TO 70 50 WRITE(pu,60) t 60 FORMAT(/' Reciprocal condition of (j**t)*j = at least', e10.2) 70 IF (MOD(iv(covprt),2) == 0) GO TO 999 iv(needhd) = 1 IF (cov1 < 0.0) THEN GO TO 80 ELSE IF (cov1 == 0.0_dp) THEN GO TO 110 ELSE GO TO 130 END IF 80 IF (-1 == cov1) WRITE(pu,90) 90 FORMAT(/' ++++++ indefinite covariance matrix ++++++') IF (-2 == cov1) WRITE(pu,100) 100 FORMAT(/' ++++++ oversize steps in computing covariance +++++') GO TO 999 110 WRITE(pu,120) 120 FORMAT(/' ++++++ covariance matrix NOT computed ++++++') GO TO 999 130 i = ABS(iv(covreq)) IF (i <= 1) WRITE(pu,140) 140 FORMAT(/' covariance = scale * h**-1 * (j**t * j) * h**-1'/ & ' where h = f.d. hessian'/) IF (i == 2) WRITE(pu,150) 150 FORMAT(/' covariance = h**-1, where h = finite-difference hessian'/) IF (i > 2) WRITE(pu,160) 160 FORMAT(/' covariance = scale * j**t * j'/) ii = cov1 - 1 DO i = 1, p i1 = ii + 1 ii = ii + i WRITE(pu,180) i, v(i1:ii) END DO 180 FORMAT(' Row', i3, ' ', 5E12.3/ (t10, 5E12.3)) 999 RETURN ! *** LAST CARD OF N2CVP FOLLOWS *** END SUBROUTINE n2cvp SUBROUTINE o7prd(l, p, s, w, y, z) ! *** FOR I = 1..L, SET S = S + W(I)*Y(.,I)*(Z(.,I)**T), I.E., ! *** ADD W(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I). INTEGER, INTENT(IN) :: l INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: s(:) REAL (dp), INTENT(IN) :: w(:) REAL (dp), INTENT(IN) :: y(:,:) ! y(p,l) REAL (dp), INTENT(IN) :: z(:,:) ! z(p,l) ! DIMENSION S(P*(P+1)/2) INTEGER :: i, j, k, m REAL (dp) :: wk, yi DO k = 1, l wk = w(k) IF (wk == zero) CYCLE m = 1 DO i = 1, p yi = wk * y(i,k) DO j = 1, i s(m) = s(m) + yi*z(j,k) m = m + 1 END DO END DO END DO RETURN ! *** LAST CARD OF O7PRD FOLLOWS *** END SUBROUTINE o7prd SUBROUTINE n2lrd(dr, iv, l, nn, p, r, rd, v) ! *** COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR ! RN2G *** ! *** PARAMETERS *** REAL (dp), INTENT(IN) :: dr(:,:) ! dr(nd,p) INTEGER, INTENT(IN) :: iv(:) REAL (dp), INTENT(IN OUT) :: l(:) INTEGER, INTENT(IN) :: nn INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN) :: r(:) REAL (dp), INTENT(OUT) :: rd(:) REAL (dp), INTENT(OUT) :: v(:) ! *** CODED BY DAVID M. GAY (WINTER 1982, FALL 1983) *** ! *** EXTERNAL FUNCTIONS AND SUBROUTINES *** ! REAL :: d7tpr ! EXTERNAL d7tpr, l7itv, l7ivm, o7prd, v7scp ! *** LOCAL VARIABLES *** INTEGER :: cov, i, j, m, step1 REAL (dp) :: a, ff, s, t, stepsize(p,1) ! *** INTRINSIC FUNCTIONS *** ! REAL :: SQRT REAL (dp), PARAMETER :: negone=-1._dp REAL (dp), SAVE :: onev(1) = 1.0_dp !++++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++ step1 = iv(step) i = iv(rdreq) IF (i <= 0) GO TO 999 IF (MOD(i,4) < 2) GO TO 30 ff = one IF (v(f) /= zero) ff = one / SQRT( ABS(v(f))) rd(1:nn) = negone DO i = 1, nn a = r(i)**2 m = step1 DO j = 1, p v(m) = dr(i,j) m = m + 1 END DO CALL l7ivm(p, v(step1:), l, v(step1:)) s = d7tpr(p, v(step1:), v(step1:)) t = one - s IF (t <= zero) CYCLE a = a * s / t rd(i) = SQRT(a) * ff END DO 30 IF (iv(mode) - p < 2) GO TO 999 ! *** COMPUTE DEFAULT COVARIANCE MATRIX *** cov = ABS(iv(h)) DO i = 1, nn m = step1 DO j = 1, p v(m) = dr(i,j) m = m + 1 END DO CALL l7ivm(p, v(step1:), l, v(step1:)) CALL l7itv(p, v(step1:), l, v(step1:)) stepsize(1:p,1) = v(step1:step1+p-1) CALL o7prd(1, p, v(cov:), onev, stepsize, stepsize) END DO 999 RETURN ! *** LAST LINE OF N2LRD FOLLOWS *** END SUBROUTINE n2lrd SUBROUTINE n2rdp(iv, rd, n, v) ! *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 *** INTEGER, INTENT(IN OUT) :: iv(:) REAL (dp), INTENT(IN) :: rd(:) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: v(:) ! *** NOTE -- V IS PASSED FOR POSSIBLE USE BY REVISED VERSIONS OF ! *** THIS ROUTINE. INTEGER :: pu !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ pu = iv(prunit) IF (pu == 0) GO TO 999 IF (iv(covprt) < 2) GO TO 999 IF (iv(regd) <= 0) GO TO 999 iv(needhd) = 1 IF (v(f) /= 0.0) THEN WRITE(pu,20) rd(1:n) 20 FORMAT(/ & ' regression diagnostic = SQRT( g(i)**t * h(i)**-1 * g(i) / ABS(f) )...'/ & (6E12.3)) ELSE WRITE(pu,40) rd(1:n) 40 FORMAT(/' regression diagnostic = SQRT( g(i)**t * h(i)**-1 * g(i) )...'/ & (6E12.3)) END IF 999 RETURN ! *** LAST LINE OF N2RDP FOLLOWS *** END SUBROUTINE n2rdp SUBROUTINE c7vfn(iv, l, lh, n, p, v) ! *** FINISH COVARIANCE COMPUTATION FOR RN2G, RNSG *** INTEGER, INTENT(IN OUT) :: iv(:) REAL (dp), INTENT(IN OUT) :: l(:) INTEGER, INTENT(IN) :: lh INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: v(:) ! EXTERNAL l7nvr, l7tsq, v7scl ! *** LOCAL VARIABLES *** INTEGER :: cov, i REAL (dp), PARAMETER :: half = 0.5_dp ! *** BODY *** iv(1) = iv(cnvcod) i = iv(mode) - p iv(mode) = 0 iv(cnvcod) = 0 IF (iv(fdh) <= 0) GO TO 999 IF ((i-2)**2 == 1) iv(regd) = 1 IF (MOD(iv(rdreq),2) /= 1) GO TO 999 ! *** FINISH COMPUTING COVARIANCE MATRIX = INVERSE OF F.D. HESSIAN. cov = ABS(iv(h)) iv(fdh) = 0 IF (iv(covmat) /= 0) GO TO 999 IF (i >= 2) GO TO 10 CALL l7nvr(p, v(cov:), l) CALL l7tsq(p, v(cov:), v(cov:)) 10 v(cov:cov+lh-1) = v(f)/(half * REAL(MAX(1,n-p))) * v(cov:cov+lh-1) iv(covmat) = cov 999 RETURN ! *** LAST LINE OF C7VFN FOLLOWS *** END SUBROUTINE c7vfn SUBROUTINE f7dhb(b, d, g, irt, iv, p, v, x) ! *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING ! *** AT V(IV(FDH)) = V(-IV(H)). HONOR SIMPLE BOUNDS IN B. ! *** IF IV(COVREQ) >= 0 THEN F7DHB USES GRADIENT DIFFERENCES, ! *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN G7LIT. ! IRT VALUES... ! 1 = COMPUTE FUNCTION VALUE, I.E., V(F). ! 2 = COMPUTE G. ! 3 = DONE. ! *** PARAMETER DECLARATIONS *** REAL (dp), INTENT(IN) :: b(:,:) ! b(2,p) REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(OUT) :: g(:) INTEGER, INTENT(OUT) :: irt INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: p REAL (dp), INTENT(OUT) :: v(:) REAL (dp), INTENT(OUT) :: x(:) ! *** LOCAL VARIABLES *** LOGICAL :: offsid INTEGER :: gsave1, hes, hmi, hpi, hpm, i, k, kind, l, m, mm1, mm1o2, & newm1, pp1o2, stpi, stpm, stp0 REAL (dp) :: del, del0, t, xm, xm1 ! *** EXTERNAL SUBROUTINES *** ! EXTERNAL v7cpy, v7scp ! V7CPY.... COPY ONE VECTOR TO ANOTHER. ! V7SCP... COPY SCALAR TO ALL COMPONENTS OF A VECTOR. ! *** SUBSCRIPTS FOR IV AND V *** REAL (dp), PARAMETER :: half=0.5_dp, hlim=0.1_dp, two=2._dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ irt = 4 kind = iv(covreq) m = iv(mode) IF (m > 0) GO TO 10 hes = ABS(iv(h)) iv(h) = -hes iv(fdh) = 0 iv(kagqt) = -1 v(fx) = v(f) ! *** SUPPLY ZEROS IN CASE B(1,I) = B(2,I) FOR SOME I *** v(hes:hes+p*(p+1)/2-1) = zero 10 IF (m > p) GO TO 999 IF (kind < 0) GO TO 120 ! *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND ! *** GRADIENT VALUES. gsave1 = iv(w0) + p IF (m > 0) GO TO 20 ! *** FIRST CALL ON F7DHB. SET GSAVE = G, TAKE FIRST STEP *** CALL v7cpy(p, v(gsave1:), g) iv(switch) = iv(nfgcal) GO TO 80 20 del = v(delta) x(m) = v(xmsave) IF (iv(toobig) == 0) GO TO 30 ! *** HANDLE OVERSIZE V(DELTA) *** del0 = v(delta0) * MAX(one/d(m), ABS(x(m))) del = half * del IF ( ABS(del/del0) <= hlim) GO TO 140 30 hes = -iv(h) ! *** SET G = (G - GSAVE)/DEL *** del = one / del DO i = 1, p g(i) = del * (g(i) - v(gsave1)) gsave1 = gsave1 + 1 END DO ! *** ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** k = hes + m*(m-1)/2 l = k + m - 2 IF (m == 1) GO TO 60 ! *** SET H(I,M) = 0.5 * (H(I,M) + G(I)) FOR I = 1 TO M-1 *** mm1 = m - 1 DO i = 1, mm1 IF (b(1,i) < b(2,i)) v(k) = half * (v(k) + g(i)) k = k + 1 END DO ! *** ADD H(I,M) = G(I) FOR I = M TO P *** 60 l = l + 1 DO i = m, p IF (b(1,i) < b(2,i)) v(l) = g(i) l = l + i END DO 80 m = m + 1 iv(mode) = m IF (m > p) GO TO 340 IF (b(1,m) >= b(2,m)) GO TO 80 ! *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE *** del = v(delta0) * MAX(one/d(m), ABS(x(m))) xm = x(m) IF (xm < zero) GO TO 90 xm1 = xm + del IF (xm1 <= b(2,m)) GO TO 110 xm1 = xm - del IF (xm1 >= b(1,m)) GO TO 100 GO TO 280 90 xm1 = xm - del IF (xm1 >= b(1,m)) GO TO 100 xm1 = xm + del IF (xm1 <= b(2,m)) GO TO 110 GO TO 280 100 del = -del 110 v(xmsave) = xm x(m) = xm1 v(delta) = del irt = 2 GO TO 999 ! *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. 120 stp0 = iv(w0) + p - 1 mm1 = m - 1 mm1o2 = m*mm1/2 hes = -iv(h) IF (m > 0) GO TO 130 ! *** FIRST CALL ON F7DHB. *** iv(savei) = 0 GO TO 240 130 IF (iv(toobig) == 0) GO TO 150 ! *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** 140 iv(fdh) = -2 GO TO 350 150 i = iv(savei) IF (i > 0) GO TO 190 ! *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** pp1o2 = p * (p-1) / 2 hpm = hes + pp1o2 + mm1 v(hpm) = v(f) ! *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** newm1 = 1 GO TO 260 160 hmi = hes + mm1o2 IF (mm1 == 0) GO TO 180 hpi = hes + pp1o2 DO i = 1, mm1 t = zero IF (b(1,i) < b(2,i)) t = v(fx) - (v(f) + v(hpi)) v(hmi) = t hmi = hmi + 1 hpi = hpi + 1 END DO 180 v(hmi) = v(f) - two*v(fx) IF (offsid) v(hmi) = v(fx) - two*v(f) ! *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** i = 0 GO TO 200 190 x(i) = v(delta) ! *** FINISH COMPUTING H(M,I) *** 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 > m) GO TO 230 IF (b(1,i) < b(2,i)) GO TO 210 GO TO 200 210 iv(savei) = i stpi = stp0 + i v(delta) = x(i) x(i) = x(i) + v(stpi) irt = 1 IF (i < 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 230 iv(savei) = 0 x(m) = v(xmsave) 240 m = m + 1 iv(mode) = m IF (m > p) GO TO 330 IF (b(1,m) < b(2,m)) GO TO 250 GO TO 240 ! *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. ! *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN ! *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. 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 <= b(2,m)) GO TO 270 offsid = .true. xm1 = xm - del IF (xm - two*del >= b(1,m)) GO TO 300 GO TO 280 270 IF (xm-del >= b(1,m)) GO TO 290 offsid = .true. IF (xm + two*del <= b(2,m)) GO TO 310 280 iv(fdh) = -2 GO TO 350 290 IF (xm >= zero) GO TO 310 xm1 = xm - del 300 del = -del 310 IF (newm1 == 1) GO TO 160 IF (newm1 == 2) GO TO 220 x(m) = xm1 stpm = stp0 + m v(stpm) = del irt = 1 GO TO 999 ! *** HANDLE SPECIAL CASE OF B(1,P) = B(2,P) -- CLEAR SCRATCH VALUES ! *** FROM LAST ROW OF FDH... 330 IF (b(1,p) < b(2,p)) GO TO 340 i = hes + p*(p-1)/2 v(i:i+p-1) = zero ! *** RESTORE V(F), ETC. *** 340 iv(fdh) = hes 350 v(f) = v(fx) irt = 3 IF (kind < 0) GO TO 999 iv(nfgcal) = iv(switch) gsave1 = iv(w0) + p CALL v7cpy(p, g, v(gsave1:)) GO TO 999 999 RETURN ! *** LAST LINE OF F7DHB FOLLOWS *** END SUBROUTINE f7dhb SUBROUTINE v7shf(n, k, x) ! *** SHIFT X(K),...,X(N) LEFT CIRCULARLY ONE POSITION *** INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: k REAL (dp), INTENT(IN OUT) :: x(:) INTEGER :: i, nm1 REAL (dp) :: t IF (k >= n) GO TO 999 nm1 = n - 1 t = x(k) DO i = k, nm1 x(i) = x(i+1) END DO x(n) = t 999 RETURN END SUBROUTINE v7shf SUBROUTINE s7bqn(b, d, dst, ipiv, ipiv1, ipiv2, kb, l, ns, & p, p1, step, td, tg, v, w, x, x0) ! *** COMPUTE BOUNDED MODIFIED NEWTON STEP *** REAL (dp), INTENT(IN) :: b(:,:) ! b(2,p) REAL (dp), INTENT(IN) :: d(:) REAL (dp), INTENT(IN OUT) :: dst(:) INTEGER, INTENT(IN OUT) :: ipiv(:) INTEGER, INTENT(OUT) :: ipiv1(:) INTEGER, INTENT(OUT) :: ipiv2(:) INTEGER, INTENT(IN OUT) :: kb REAL (dp), INTENT(IN OUT) :: l(:) INTEGER, INTENT(OUT) :: ns INTEGER, INTENT(IN) :: p INTEGER, INTENT(IN OUT) :: p1 REAL (dp), INTENT(IN OUT) :: step(:) REAL (dp), INTENT(IN OUT) :: td(:) REAL (dp), INTENT(OUT) :: tg(:) REAL (dp), INTENT(IN OUT) :: v(:) REAL (dp), INTENT(OUT) :: w(:) REAL (dp), INTENT(IN OUT) :: x(:) REAL (dp), INTENT(IN OUT) :: x0(:) ! DIMENSION L(P*(P+1)/2) ! REAL :: d7tpr, r7mdc, v2nrm ! EXTERNAL d7tpr, i7shft, l7itv, l7ivm, q7rsh, r7mdc, v2nrm, & ! v2axy, v7cpy, v7ipr, v7scp, v7shf ! *** LOCAL VARIABLES *** INTEGER :: i, j, k, p0, p1m1 REAL (dp) :: alpha, dst0, dst1, dstmax, dstmin, dx, gts, t, ti, t1, xi REAL (dp), SAVE :: fudge=1.0001_dp, half=0.5_dp, meps2=0._dp, two=2._dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ dstmax = fudge * (one + v(phmxfc)) * v(radius) dstmin = (one + v(phmnfc)) * v(radius) dst1 = zero IF (meps2 <= zero) meps2 = two * r7mdc(3) p0 = p1 ns = 0 DO i = 1, p ipiv1(i) = i ipiv2(i) = i END DO w(1:p1) = -step(1:p1) * td(1:p1) alpha = ABS(v(stppar)) v(preduc) = zero gts = -v(gtstep) IF (kb < 0) dst(1:p) = zero kb = 1 ! *** -W = D TIMES RESTRICTED NEWTON STEP FROM X + DST/D. ! *** FIND T SUCH THAT X - T*W IS STILL FEASIBLE. 30 t = one k = 0 DO i = 1, p1 j = ipiv(i) dx = w(i) / d(j) xi = x(j) - dx IF (xi < b(1,j)) GO TO 40 IF (xi <= b(2,j)) CYCLE ti = ( x(j) - b(2,j) ) / dx k = i GO TO 50 40 ti = ( x(j) - b(1,j) ) / dx k = -i 50 IF (t <= ti) CYCLE t = ti END DO IF (p > p1) step(p1+1:p) = dst(p1+1:p) CALL v2axy(p1, step, -t, w, dst) dst0 = dst1 dst1 = v2nrm(p, step) ! *** CHECK FOR OVERSIZE STEP *** IF (dst1 <= dstmax) GO TO 80 IF (p1 >= p0) GO TO 70 IF (dst0 < dstmin) kb = 0 GO TO 110 70 k = 0 ! *** UPDATE DST, TG, AND V(PREDUC) *** 80 v(dstnrm) = dst1 dst(1:p1) = step(1:p1) t1 = one - t tg(1:p1) = t1 * tg(1:p1) IF (alpha > zero) CALL v2axy(p1, tg, t*alpha, w, tg) v(preduc) = v(preduc) + t*((one - half*t)*gts + half*alpha*t* d7tpr(p1, w, w)) IF (k == 0) GO TO 110 ! *** PERMUTE L, ETC. IF NECESSARY *** p1m1 = p1 - 1 j = ABS(k) IF (j == p1) GO TO 100 ns = ns + 1 ipiv2(p1) = j CALL q7rsh(j, p1, .false., tg, l, w) CALL i7shft(p1, j, ipiv) CALL i7shft(p1, j, ipiv1) CALL v7shf(p1, j, tg) CALL v7shf(p1, j, dst) 100 IF (k < 0) ipiv(p1) = -ipiv(p1) p1 = p1m1 IF (p1 <= 0) GO TO 110 CALL l7ivm(p1, w, l, tg) gts = d7tpr(p1, w, w) CALL l7itv(p1, w, l, w) GO TO 30 ! *** UNSCALE STEP *** 110 DO i = 1, p j = ABS(ipiv(i)) step(j) = dst(i) / d(j) END DO ! *** FUDGE STEP TO ENSURE THAT IT FORCES APPROPRIATE COMPONENTS ! *** TO THEIR BOUNDS *** IF (p1 >= p0) GO TO 150 k = p1 + 1 DO i = k, p0 j = ipiv(i) t = meps2 IF (j > 0) GO TO 130 t = -t j = -j ipiv(i) = j 130 t = t * MAX( ABS(x(j)), ABS(x0(j))) step(j) = step(j) + t END DO 150 CALL v2axy(p, x, one, step, x0) IF (ns > 0) CALL v7ipr(p0, ipiv1, td) RETURN ! *** LAST LINE OF S7BQN FOLLOWS *** END SUBROUTINE s7bqn SUBROUTINE d7mlp(n, x, y, z, k) ! *** SET X = DIAG(Y)**K * Z ! *** FOR X, Z = LOWER TRIANG. MATRICES STORED COMPACTLY BY ROW ! *** K = 1 OR -1. INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: y(:), z(:) INTEGER, INTENT(IN) :: k INTEGER :: i, j, l REAL (dp) :: t REAL (dp), PARAMETER :: one = 1.0_dp l = 1 IF (k < 0) THEN DO i = 1, n t = one / y(i) DO j = 1, i x(l) = t * z(l) l = l + 1 END DO END DO ELSE DO i = 1, n t = y(i) DO j = 1, i x(l) = t * z(l) l = l + 1 END DO END DO END IF RETURN ! *** LAST CARD OF D7MLP FOLLOWS *** END SUBROUTINE d7mlp SUBROUTINE g7qsb(b, d, dihdi, g, ipiv, ipiv1, ipiv2, ka, l, & p, p0, pc, step, td, tg, v, w, x, x0) ! *** COMPUTE HEURISTIC BOUNDED NEWTON STEP *** REAL (dp), INTENT(IN) :: b(:,:) ! b(2,p) REAL (dp), INTENT(IN OUT) :: d(:) REAL (dp), INTENT(IN OUT) :: dihdi(:) REAL (dp), INTENT(IN) :: g(:) INTEGER, INTENT(IN OUT) :: ipiv(:) INTEGER, INTENT(IN OUT) :: ipiv1(:) INTEGER, INTENT(IN OUT) :: ipiv2(:) INTEGER, INTENT(OUT) :: ka REAL (dp), INTENT(IN OUT) :: l(:) INTEGER, INTENT(IN) :: p INTEGER, INTENT(OUT) :: p0 INTEGER, INTENT(IN) :: pc REAL (dp), INTENT(IN OUT) :: step(:,:) ! step(p,2) REAL (dp), INTENT(IN OUT) :: td(:) REAL (dp), INTENT(IN OUT) :: tg(:) REAL (dp), INTENT(IN OUT) :: v(:) REAL (dp), INTENT(IN OUT) :: w(:) REAL (dp), INTENT(IN OUT) :: x(:) REAL (dp), INTENT(IN OUT) :: x0(:) ! DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2) ! REAL :: d7tpr ! EXTERNAL d7tpr, g7qts, s7bqn, s7ipr, v7cpy, v7ipr, v7scp, v7vmp ! *** LOCAL VARIABLES *** INTEGER :: k, kb, kinit, ns, p1, p10 REAL (dp) :: ds0, nred, pred, rad REAL (dp), PARAMETER :: zero = 0.0_dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ p1 = pc IF (ka < 0) GO TO 10 nred = v(nreduc) ds0 = v(dst0) GO TO 20 10 p0 = 0 ka = -1 20 kinit = -1 IF (p0 == p1) kinit = ka CALL v7cpy(p, x, x0) pred = zero rad = v(radius) kb = -1 v(dstnrm) = zero IF (p1 > 0) GO TO 30 nred = zero ds0 = zero step(1:p,1) = zero GO TO 60 30 CALL v7cpy(p, td, d) CALL v7ipr(p, ipiv, td) CALL v7vmp(p, tg, g, d, -1) CALL v7ipr(p, ipiv, tg) 40 k = kinit kinit = -1 v(radius) = rad - v(dstnrm) CALL g7qts(td, tg, dihdi, k, l, p1, step(:,1), v) p0 = p1 IF (ka >= 0) GO TO 50 nred = v(nreduc) ds0 = v(dst0) 50 ka = k v(radius) = rad p10 = p1 CALL s7bqn(b, d, step(:,2), ipiv, ipiv1, ipiv2, kb, l, & ns, p, p1, step(:,1), td, tg, v, w, x, x0) IF (ns > 0) CALL s7ipr(p10, ipiv1, dihdi) pred = pred + v(preduc) IF (ns /= 0) p0 = 0 IF (kb <= 0) GO TO 40 60 v(dst0) = ds0 v(nreduc) = nred v(preduc) = pred v(gtstep) = d7tpr(p, g, step(:,1)) RETURN ! *** LAST LINE OF G7QSB FOLLOWS *** END SUBROUTINE g7qsb SUBROUTINE i7pnvr(n, x, y) ! *** SET PERMUTATION VECTOR X TO INVERSE OF Y *** INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: x(:) INTEGER, INTENT(IN) :: y(:) INTEGER :: i, j DO i = 1, n j = y(i) x(j) = i END DO RETURN ! *** LAST LINE OF I7PNVR FOLLOWS *** END SUBROUTINE i7pnvr SUBROUTINE i7shft(n, k, x) ! *** SHIFT X(K),...,X(N) LEFT CIRCULARLY ONE POSITION *** INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: k INTEGER, INTENT(IN OUT) :: x(:) INTEGER :: i, nm1, t IF (k >= n) GO TO 999 nm1 = n - 1 t = x(k) DO i = k, nm1 x(i) = x(i+1) END DO x(n) = t 999 RETURN END SUBROUTINE i7shft SUBROUTINE l7msb(b, d, g, ierr, ipiv, ipiv1, ipiv2, ka, lmat, & p, p0, pc, qtr, rmat, step, td, tg, v, w, wlm, x, x0) ! *** COMPUTE HEURISTIC BOUNDED NEWTON STEP *** REAL (dp), INTENT(IN) :: b(:,:) ! b(2,p) REAL (dp), INTENT(IN OUT) :: d(:) REAL (dp), INTENT(IN) :: g(:) INTEGER, INTENT(IN OUT) :: ierr INTEGER, INTENT(IN OUT) :: ipiv(:) INTEGER, INTENT(IN OUT) :: ipiv1(:) INTEGER, INTENT(IN OUT) :: ipiv2(:) INTEGER, INTENT(OUT) :: ka REAL (dp), INTENT(IN OUT) :: lmat(:) INTEGER, INTENT(IN) :: p INTEGER, INTENT(OUT) :: p0 INTEGER, INTENT(IN) :: pc REAL (dp), INTENT(IN OUT) :: qtr(:) REAL (dp), INTENT(IN OUT) :: rmat(:) REAL (dp), INTENT(IN OUT) :: step(:,:) ! step(p,3) REAL (dp), INTENT(IN OUT) :: td(:) REAL (dp), INTENT(IN OUT) :: tg(:) REAL (dp), INTENT(IN OUT) :: v(:) REAL (dp), INTENT(IN OUT) :: w(:) REAL (dp), INTENT(IN OUT) :: wlm(:) REAL (dp), INTENT(IN OUT) :: x(:) REAL (dp), INTENT(IN OUT) :: x0(:) ! DIMENSION LMAT(P*(P+1)/2), RMAT(P*(P+1)/2), WLM(P*(P+5)/2 + 4) ! REAL :: d7tpr ! EXTERNAL d7mlp, d7tpr, l7mst, l7tvm, q7rsh, s7bqn, & ! v2axy, v7cpy, v7ipr, v7scp, v7vmp ! *** LOCAL VARIABLES *** INTEGER :: i, j, k, k0, kb, kinit, l, ns, p1, p10, p11 REAL (dp) :: ds0, nred, pred, rad REAL (dp), PARAMETER :: one = 1._dp, zero = 0._dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ p1 = pc IF (ka < 0) GO TO 10 nred = v(nreduc) ds0 = v(dst0) GO TO 20 10 p0 = 0 ka = -1 20 kinit = -1 IF (p0 == p1) kinit = ka CALL v7cpy(p, x, x0) CALL v7cpy(p, td, d) ! *** USE STEP(1,3) AS TEMP. COPY OF QTR *** CALL v7cpy(p, step(:,3), qtr) CALL v7ipr(p, ipiv, td) pred = zero rad = v(radius) kb = -1 v(dstnrm) = zero IF (p1 > 0) GO TO 30 nred = zero ds0 = zero step(1:p,1) = zero GO TO 90 30 CALL v7vmp(p, tg, g, d, -1) CALL v7ipr(p, ipiv, tg) p10 = p1 40 k = kinit kinit = -1 v(radius) = rad - v(dstnrm) CALL v7vmp(p1, tg, tg, td, 1) DO i = 1, p1 ipiv1(i) = i END DO k0 = MAX(0, k) CALL l7mst(td, tg, ierr, ipiv1, k, p1, step(:,3), rmat, step(:,1), v, w) CALL v7vmp(p1, tg, tg, td, -1) p0 = p1 IF (ka >= 0) GO TO 60 nred = v(nreduc) ds0 = v(dst0) 60 ka = k v(radius) = rad l = p1 + 5 IF (k <= k0) CALL d7mlp(p1, lmat, td, rmat, -1) IF (k > k0) CALL d7mlp(p1, lmat, td, wlm(l:), -1) CALL s7bqn(b, d, step(:,2), ipiv, ipiv1, ipiv2, kb, lmat, & ns, p, p1, step(:,1), td, tg, v, w, x, x0) pred = pred + v(preduc) IF (ns == 0) GO TO 80 p0 = 0 ! *** UPDATE RMAT AND QTR *** p11 = p1 + 1 l = p10 + p11 DO k = p11, p10 j = l - k i = ipiv2(j) IF (i < j) CALL q7rsh(i, j, .true., qtr, rmat, w) END DO 80 IF (kb > 0) GO TO 90 ! *** UPDATE LOCAL COPY OF QTR *** CALL v7vmp(p10, w, step(:,2), td, -1) CALL l7tvm(p10, w, lmat, w) CALL v2axy(p10, step(:,3), one, w, qtr) GO TO 40 90 v(dst0) = ds0 v(nreduc) = nred v(preduc) = pred v(gtstep) = d7tpr(p, g, step(:,1)) RETURN ! *** LAST LINE OF L7MSB FOLLOWS *** END SUBROUTINE l7msb SUBROUTINE h2rfa(n, a, b, x, y, z) ! *** APPLY 2X2 HOUSEHOLDER REFLECTION DETERMINED BY X, Y, Z TO ! *** N-VECTORS A, B *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: a(:), b(:) REAL (dp), INTENT(IN) :: x, y, z INTEGER :: i REAL (dp) :: t DO i = 1, n t = a(i)*x + b(i)*y a(i) = a(i) + t b(i) = b(i) + t*z END DO RETURN ! *** LAST LINE OF H2RFA FOLLOWS *** END SUBROUTINE h2rfa SUBROUTINE h2rfg(a, b, x, y, z, fn_val) ! *** DETERMINE X, Y, Z SO I + (1,Z)**T * (X,Y) IS A 2X2 ! *** HOUSEHOLDER REFLECTION SENDING (A,B)**T INTO (C,0)**T, ! *** WHERE C = -SIGN(A)*SQRT(A**2 + B**2) IS THE VALUE H2RFG RETURNS. REAL (dp), INTENT(IN) :: a, b REAL (dp), INTENT(OUT) :: x, y, z, fn_val REAL (dp) :: a1, b1, c, t ! REAL :: SQRT ! *** BODY *** IF (b == zero) THEN x = zero y = zero z = zero fn_val = a ELSE t = ABS(a) + ABS(b) a1 = a / t b1 = b / t c = SQRT(a1**2 + b1**2) IF (a1 > zero) c = -c a1 = a1 - c z = b1 / a1 x = a1 / c y = b1 / c fn_val = t * c END IF RETURN ! *** LAST LINE OF H2RFG FOLLOWS *** END SUBROUTINE h2rfg SUBROUTINE q7rsh(k, p, havqtr, qtr, r, w) ! *** PERMUTE COLUMN K OF R TO COLUMN P, MODIFY QTR ACCORDINGLY *** INTEGER, INTENT(IN) :: k INTEGER, INTENT(IN) :: p LOGICAL, INTENT(IN) :: havqtr REAL (dp), INTENT(IN OUT) :: qtr(:) REAL (dp), INTENT(IN OUT) :: r(:) REAL (dp), INTENT(IN OUT) :: w(:) ! DIMENSION R(P*(P+1)/2) ! REAL :: h2rfg ! EXTERNAL h2rfa, h2rfg, v7cpy ! *** LOCAL VARIABLES *** INTEGER :: i, i1, j, jm1, jp1, j1, km1, k1, pm1 REAL (dp) :: a, b, t, wj, x, y, z REAL (dp), PARAMETER :: zero = 0.0_dp !+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ IF (k >= p) GO TO 999 km1 = k - 1 k1 = k * km1 / 2 w(1:k) = r(k1+1:k1+k) wj = w(k) pm1 = p - 1 j1 = k1 + km1 DO j = k, pm1 jm1 = j - 1 jp1 = j + 1 IF (jm1 > 0) r(k1+1:k1+jm1) = r(j1+2:j1+jm1+1) j1 = j1 + jp1 k1 = k1 + j a = r(j1) b = r(j1+1) IF (b /= zero) GO TO 10 r(k1) = a x = zero z = zero GO TO 40 10 CALL h2rfg(a, b, x, y, z, r(k1)) IF (j == pm1) GO TO 30 i1 = j1 DO i = jp1, pm1 i1 = i1 + i CALL h2rfa(1, r(i1:), r(i1+1:), x, y, z) END DO 30 IF (havqtr) CALL h2rfa(1, qtr(j:), qtr(jp1:), x, y, z) 40 t = x * wj w(j) = wj + t wj = t * z END DO w(p) = wj r(k1+1:k1+p) = w(1:p) 999 RETURN END SUBROUTINE q7rsh SUBROUTINE s7dmp(n, x, y, z, k) ! *** SET X = DIAG(Z)**K * Y * DIAG(Z)**K ! *** FOR X, Y = COMPACTLY STORED LOWER TRIANG. MATRICES ! *** K = 1 OR -1. INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(IN) :: z(:) INTEGER, INTENT(IN) :: k INTEGER :: i, j, l REAL (dp) :: t REAL (dp), PARAMETER :: one = 1.0_dp l = 1 IF (k < 0) THEN DO i = 1, n t = one / z(i) DO j = 1, i x(l) = t * y(l) / z(j) l = l + 1 END DO END DO ELSE DO i = 1, n t = z(i) DO j = 1, i x(l) = t * y(l) * z(j) l = l + 1 END DO END DO END IF RETURN ! *** LAST CARD OF S7DMP FOLLOWS *** END SUBROUTINE s7dmp SUBROUTINE s7ipr(p, ip, h) ! APPLY THE PERMUTATION DEFINED BY IP TO THE ROWS AND COLUMNS OF THE ! P X P SYMMETRIC MATRIX WHOSE LOWER TRIANGLE IS STORED COMPACTLY IN H. ! THUS H.OUTPUT(I,J) = H.INPUT(IP(I), IP(J)). INTEGER, INTENT(IN) :: p INTEGER, INTENT(IN OUT) :: ip(:) REAL (dp), INTENT(IN OUT) :: h(:) INTEGER :: i, j, j1, jm, k, k1, kk, km, kmj, l, m REAL (dp) :: t ! *** BODY *** DO i = 1, p j = ip(i) IF (j == i) CYCLE ip(i) = ABS(j) IF (j < 0) CYCLE k = i 10 j1 = j k1 = k IF (j <= k) GO TO 20 j1 = k k1 = j 20 kmj = k1-j1 l = j1-1 jm = j1*l/2 km = k1*(k1-1)/2 IF (l <= 0) GO TO 40 DO m = 1, l jm = jm+1 t = h(jm) km = km+1 h(jm) = h(km) h(km) = t END DO 40 km = km+1 kk = km+kmj jm = jm+1 t = h(jm) h(jm) = h(kk) h(kk) = t j1 = l l = kmj-1 IF (l <= 0) GO TO 60 DO m = 1, l jm = jm+j1+m t = h(jm) km = km+1 h(jm) = h(km) h(km) = t END DO 60 IF (k1 >= p) GO TO 80 l = p-k1 k1 = k1-1 km = kk DO m = 1, l km = km+k1+m jm = km-kmj t = h(jm) h(jm) = h(km) h(km) = t END DO 80 k = j j = ip(k) ip(k) = -j IF (j > i) GO TO 10 END DO RETURN ! *** LAST LINE OF S7IPR FOLLOWS *** END SUBROUTINE s7ipr SUBROUTINE v7ipr(n, ip, x) ! PERMUTE X SO THAT X.OUTPUT(I) = X.INPUT(IP(I)). ! IP IS UNCHANGED ON OUTPUT. INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN OUT) :: ip(:) REAL (dp), INTENT(IN OUT) :: x(:) INTEGER :: i, j, k REAL (dp) :: t DO i = 1, n j = ip(i) IF (j == i) CYCLE IF (j > 0) GO TO 10 ip(i) = -j CYCLE 10 t = x(i) k = i 20 x(k) = x(j) k = j j = ip(k) ip(k) = -j IF (j > i) GO TO 20 x(k) = t END DO RETURN ! *** LAST LINE OF V7IPR FOLLOWS *** END SUBROUTINE v7ipr SUBROUTINE v7vmp(n, x, y, z, k) ! *** SET X(I) = Y(I) * Z(I)**K, 1 <= I <= N (FOR K = 1 OR -1) *** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: x(:) REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(IN) :: z(:) INTEGER, INTENT(IN) :: k IF (k < 0) THEN x(1:n) = y(1:n) / z(1:n) ELSE x(1:n) = y(1:n) * z(1:n) END IF RETURN ! *** LAST CARD OF V7VMP FOLLOWS *** END SUBROUTINE v7vmp END MODULE nlsol