+++ /dev/null
- subroutine assst(iv, liv, lv, v)
-c
-c *** assess candidate step (***sol version 2.3) ***
-c
- integer liv, l
- integer iv(liv)
- double precision v(lv)
-c
-c *** purpose ***
-c
-c this subroutine is called by an unconstrained minimization
-c routine to assess the next candidate step. it may recommend one
-c of several courses of action, such as accepting the step, recom-
-c puting it using the same or a new quadratic model, or halting due
-c to convergence or false convergence. see the return code listing
-c below.
-c
-c-------------------------- parameter usage --------------------------
-c
-c iv (i/o) integer parameter and scratch vector -- see description
-c below of iv values referenced.
-c liv (in) length of iv array.
-c lv (in) length of v array.
-c v (i/o) real parameter and scratch vector -- see description
-c below of v values referenced.
-c
-c *** iv values referenced ***
-c
-c iv(irc) (i/o) on input for the first step tried in a new iteration,
-c iv(irc) should be set to 3 or 4 (the value to which it is
-c set when step is definitely to be accepted). on input
-c after step has been recomputed, iv(irc) should be
-c unchanged since the previous return of assst.
-c on output, iv(irc) is a return code having one of the
-c following values...
-c 1 = switch models or try smaller step.
-c 2 = switch models or accept step.
-c 3 = accept step and determine v(radfac) by gradient
-c tests.
-c 4 = accept step, v(radfac) has been determined.
-c 5 = recompute step (using the same model).
-c 6 = recompute step with radius = v(lmaxs) but do not
-c evaulate the objective function.
-c 7 = x-convergence (see v(xctol)).
-c 8 = relative function convergence (see v(rfctol)).
-c 9 = both x- and relative function convergence.
-c 10 = absolute function convergence (see v(afctol)).
-c 11 = singular convergence (see v(lmaxs)).
-c 12 = false convergence (see v(xftol)).
-c 13 = iv(irc) was out of range on input.
-c return code i has precdence over i+1 for i = 9, 10, 11.
-c iv(mlstgd) (i/o) saved value of iv(model).
-c iv(model) (i/o) on input, iv(model) should be an integer identifying
-c the current quadratic model of the objective function.
-c if a previous step yielded a better function reduction,
-c then iv(model) will be set to iv(mlstgd) on output.
-c iv(nfcall) (in) invocation count for the objective function.
-c iv(nfgcal) (i/o) value of iv(nfcall) at step that gave the biggest
-c function reduction this iteration. iv(nfgcal) remains
-c unchanged until a function reduction is obtained.
-c iv(radinc) (i/o) the number of radius increases (or minus the number
-c of decreases) so far this iteration.
-c iv(restor) (out) set to 1 if v(f) has been restored and x should be
-c restored to its initial value, to 2 if x should be saved,
-c to 3 if x should be restored from the saved value, and to
-c 0 otherwise.
-c iv(stage) (i/o) count of the number of models tried so far in the
-c current iteration.
-c iv(stglim) (in) maximum number of models to consider.
-c iv(switch) (out) set to 0 unless a new model is being tried and it
-c gives a smaller function value than the previous model,
-c in which case assst sets iv(switch) = 1.
-c iv(toobig) (in) is nonzero if step was too big (e.g. if it caused
-c overflow).
-c iv(xirc) (i/o) value that iv(irc) would have in the absence of
-c convergence, false convergence, and oversized steps.
-c
-c *** v values referenced ***
-c
-c v(afctol) (in) absolute function convergence tolerance. if the
-c absolute value of the current function value v(f) is less
-c than v(afctol), then assst returns with iv(irc) = 10.
-c v(decfac) (in) factor by which to decrease radius when iv(toobig) is
-c nonzero.
-c v(dstnrm) (in) the 2-norm of d*step.
-c v(dstsav) (i/o) value of v(dstnrm) on saved step.
-c v(dst0) (in) the 2-norm of d times the newton step (when defined,
-c i.e., for v(nreduc) .ge. 0).
-c v(f) (i/o) on both input and output, v(f) is the objective func-
-c tion value at x. if x is restored to a previous value,
-c then v(f) is restored to the corresponding value.
-c v(fdif) (out) the function reduction v(f0) - v(f) (for the output
-c value of v(f) if an earlier step gave a bigger function
-c decrease, and for the input value of v(f) otherwise).
-c v(flstgd) (i/o) saved value of v(f).
-c v(f0) (in) objective function value at start of iteration.
-c v(gtslst) (i/o) value of v(gtstep) on saved step.
-c v(gtstep) (in) inner product between step and gradient.
-c v(incfac) (in) minimum factor by which to increase radius.
-c v(lmaxs) (in) maximum reasonable step size (and initial step bound).
-c if the actual function decrease is no more than twice
-c what was predicted, if a return with iv(irc) = 7, 8, 9,
-c or 10 does not occur, if v(dstnrm) .gt. v(lmaxs), and if
-c v(preduc) .le. v(sctol) * abs(v(f0)), then assst re-
-c turns with iv(irc) = 11. if so doing appears worthwhile,
-c then assst repeats this test with v(preduc) computed for
-c a step of length v(lmaxs) (by a return with iv(irc) = 6).
-c v(nreduc) (i/o) function reduction predicted by quadratic model for
-c newton step. if assst is called with iv(irc) = 6, i.e.,
-c if v(preduc) has been computed with radius = v(lmaxs) for
-c use in the singular convervence test, then v(nreduc) is
-c set to -v(preduc) before the latter is restored.
-c v(plstgd) (i/o) value of v(preduc) on saved step.
-c v(preduc) (i/o) function reduction predicted by quadratic model for
-c current step.
-c v(radfac) (out) factor to be used in determining the new radius,
-c which should be v(radfac)*dst, where dst is either the
-c output value of v(dstnrm) or the 2-norm of
-c diag(newd)*step for the output value of step and the
-c updated version, newd, of the scale vector d. for
-c iv(irc) = 3, v(radfac) = 1.0 is returned.
-c v(rdfcmn) (in) minimum value for v(radfac) in terms of the input
-c value of v(dstnrm) -- suggested value = 0.1.
-c v(rdfcmx) (in) maximum value for v(radfac) -- suggested value = 4.0.
-c v(reldx) (in) scaled relative change in x caused by step, computed
-c (e.g.) by function reldst as
-c max (d(i)*abs(x(i)-x0(i)), 1 .le. i .le. p) /
-c max (d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p).
-c v(rfctol) (in) relative function convergence tolerance. if the
-c actual function reduction is at most twice what was pre-
-c dicted and v(nreduc) .le. v(rfctol)*abs(v(f0)), then
-c assst returns with iv(irc) = 8 or 9.
-c v(stppar) (in) marquardt parameter -- 0 means full newton step.
-c v(tuner1) (in) tuning constant used to decide if the function
-c reduction was much less than expected. suggested
-c value = 0.1.
-c v(tuner2) (in) tuning constant used to decide if the function
-c reduction was large enough to accept step. suggested
-c value = 10**-4.
-c v(tuner3) (in) tuning constant used to decide if the radius
-c should be increased. suggested value = 0.75.
-c v(xctol) (in) x-convergence criterion. if step is a newton step
-c (v(stppar) = 0) having v(reldx) .le. v(xctol) and giving
-c at most twice the predicted function decrease, then
-c assst returns iv(irc) = 7 or 9.
-c v(xftol) (in) false convergence tolerance. if step gave no or only
-c a small function decrease and v(reldx) .le. v(xftol),
-c then assst returns with iv(irc) = 12.
-c
-c------------------------------- notes -------------------------------
-c
-c *** application and usage restrictions ***
-c
-c this routine is called as part of the nl2sol (nonlinear
-c least-squares) package. it may be used in any unconstrained
-c minimization solver that uses dogleg, goldfeld-quandt-trotter,
-c or levenberg-marquardt steps.
-c
-c *** algorithm notes ***
-c
-c see (1) for further discussion of the assessing and model
-c switching strategies. while nl2sol considers only two models,
-c assst is designed to handle any number of models.
-c
-c *** usage notes ***
-c
-c on the first call of an iteration, only the i/o variables
-c step, x, iv(irc), iv(model), v(f), v(dstnrm), v(gtstep), and
-c v(preduc) need have been initialized. between calls, no i/o
-c values execpt step, x, iv(model), v(f) and the stopping toler-
-c ances should be changed.
-c after a return for convergence or false convergence, one can
-c change the stopping tolerances and call assst again, in which
-c case the stopping tests will be repeated.
-c
-c *** references ***
-c
-c (1) dennis, j.e., jr., gay, d.m., and welsch, r.e. (1981),
-c an adaptive nonlinear least-squares algorithm,
-c acm trans. math. software, vol. 7, no. 3.
-c
-c (2) powell, m.j.d. (1970) a fortran subroutine for solving
-c systems of nonlinear algebraic equations, in numerical
-c methods for nonlinear algebraic equations, edited by
-c p. rabinowitz, gordon and breach, london.
-c
-c *** history ***
-c
-c john dennis designed much of this routine, starting with
-c ideas in (2). roy welsch suggested the model switching strategy.
-c david gay and stephen peters cast this subroutine into a more
-c portable form (winter 1977), and david gay cast it into its
-c present form (fall 1978).
-c
-c *** general ***
-c
-c this subroutine was written in connection with research
-c supported by the national science foundation under grants
-c mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
-c mcs-7906671.
-c
-c------------------------ external quantities ------------------------
-c
-c *** no external functions and subroutines ***
-c
-c *** intrinsic functions ***
-c/+
- double precision dabs, dmax1
-c/
-c *** no common blocks ***
-c
-c-------------------------- local variables --------------------------
-c
- logical goodx
- integer i, nfc
- double precision emax, emaxs, gts, rfac1, xmax
- double precision half, one, onep2, two, zero
-c
-c *** subscripts for iv and v ***
-c
- integer afctol, decfac, dstnrm, dstsav, dst0, f, fdif, flstgd, f0,
- 1 gtslst, gtstep, incfac, irc, lmaxs, mlstgd, model, nfcall,
- 2 nfgcal, nreduc, plstgd, preduc, radfac, radinc, rdfcmn,
- 3 rdfcmx, reldx, restor, rfctol, sctol, stage, stglim,
- 4 stppar, switch, toobig, tuner1, tuner2, tuner3, xctol,
- 5 xftol, xirc
-c
-c *** data initializations ***
-c
-c/6
-c data half/0.5d+0/, one/1.d+0/, onep2/1.2d+0/, two/2.d+0/,
-c 1 zero/0.d+0/
-c/7
- parameter (half=0.5d+0, one=1.d+0, onep2=1.2d+0, two=2.d+0,
- 1 zero=0.d+0)
-c/
-c
-c/6
-c data irc/29/, mlstgd/32/, model/5/, nfcall/6/, nfgcal/7/,
-c 1 radinc/8/, restor/9/, stage/10/, stglim/11/, switch/12/,
-c 2 toobig/2/, xirc/13/
-c/7
- parameter (irc=29, mlstgd=32, model=5, nfcall=6, nfgcal=7,
- 1 radinc=8, restor=9, stage=10, stglim=11, switch=12,
- 2 toobig=2, xirc=13)
-c/
-c/6
-c data afctol/31/, decfac/22/, dstnrm/2/, dst0/3/, dstsav/18/,
-c 1 f/10/, fdif/11/, flstgd/12/, f0/13/, gtslst/14/, gtstep/4/,
-c 2 incfac/23/, lmaxs/36/, nreduc/6/, plstgd/15/, preduc/7/,
-c 3 radfac/16/, rdfcmn/24/, rdfcmx/25/, reldx/17/, rfctol/32/,
-c 4 sctol/37/, stppar/5/, tuner1/26/, tuner2/27/, tuner3/28/,
-c 5 xctol/33/, xftol/34/
-c/7
- parameter (afctol=31, decfac=22, dstnrm=2, dst0=3, dstsav=18,
- 1 f=10, fdif=11, flstgd=12, f0=13, gtslst=14, gtstep=4,
- 2 incfac=23, lmaxs=36, nreduc=6, plstgd=15, preduc=7,
- 3 radfac=16, rdfcmn=24, rdfcmx=25, reldx=17, rfctol=32,
- 4 sctol=37, stppar=5, tuner1=26, tuner2=27, tuner3=28,
- 5 xctol=33, xftol=34)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- nfc = iv(nfcall)
- iv(switch) = 0
- iv(restor) = 0
- rfac1 = one
- goodx = .true.
- i = iv(irc)
- if (i .ge. 1 .and. i .le. 12)
- 1 go to (20,30,10,10,40,280,220,220,220,220,220,170), i
- iv(irc) = 13
- go to 999
-c
-c *** initialize for new iteration ***
-c
- 10 iv(stage) = 1
- iv(radinc) = 0
- v(flstgd) = v(f0)
- if (iv(toobig) .eq. 0) go to 110
- iv(stage) = -1
- iv(xirc) = i
- go to 60
-c
-c *** step was recomputed with new model or smaller radius ***
-c *** first decide which ***
-c
- 20 if (iv(model) .ne. iv(mlstgd)) go to 30
-c *** old model retained, smaller radius tried ***
-c *** do not consider any more new models this iteration ***
- iv(stage) = iv(stglim)
- iv(radinc) = -1
- go to 110
-c
-c *** a new model is being tried. decide whether to keep it. ***
-c
- 30 iv(stage) = iv(stage) + 1
-c
-c *** now we add the possibiltiy that step was recomputed with ***
-c *** the same model, perhaps because of an oversized step. ***
-c
- 40 if (iv(stage) .gt. 0) go to 50
-c
-c *** step was recomputed because it was too big. ***
-c
- if (iv(toobig) .ne. 0) go to 60
-c
-c *** restore iv(stage) and pick up where we left off. ***
-c
- iv(stage) = -iv(stage)
- i = iv(xirc)
- go to (20, 30, 110, 110, 70), i
-c
- 50 if (iv(toobig) .eq. 0) go to 70
-c
-c *** handle oversize step ***
-c
- if (iv(radinc) .gt. 0) go to 80
- iv(stage) = -iv(stage)
- iv(xirc) = iv(irc)
-c
- 60 v(radfac) = v(decfac)
- iv(radinc) = iv(radinc) - 1
- iv(irc) = 5
- iv(restor) = 1
- go to 999
-c
- 70 if (v(f) .lt. v(flstgd)) go to 110
-c
-c *** the new step is a loser. restore old model. ***
-c
- if (iv(model) .eq. iv(mlstgd)) go to 80
- iv(model) = iv(mlstgd)
- iv(switch) = 1
-c
-c *** restore step, etc. only if a previous step decreased v(f).
-c
- 80 if (v(flstgd) .ge. v(f0)) go to 110
- iv(restor) = 1
- v(f) = v(flstgd)
- v(preduc) = v(plstgd)
- v(gtstep) = v(gtslst)
- if (iv(switch) .eq. 0) rfac1 = v(dstnrm) / v(dstsav)
- v(dstnrm) = v(dstsav)
- nfc = iv(nfgcal)
- goodx = .false.
-c
- 110 v(fdif) = v(f0) - v(f)
- if (v(fdif) .gt. v(tuner2) * v(preduc)) go to 140
- if(iv(radinc).gt.0) go to 140
-c
-c *** no (or only a trivial) function decrease
-c *** -- so try new model or smaller radius
-c
- if (v(f) .lt. 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) .lt. iv(stglim)) go to 160
- iv(irc) = 5
- iv(radinc) = iv(radinc) - 1
- go to 160
-c
-c *** nontrivial function decrease achieved ***
-c
- 140 iv(nfgcal) = nfc
- rfac1 = one
- v(dstsav) = v(dstnrm)
- if (v(fdif) .gt. v(preduc)*v(tuner1)) go to 190
-c
-c *** decrease was much less than predicted -- either change models
-c *** or accept step with decreased radius.
-c
- if (iv(stage) .ge. iv(stglim)) go to 150
-c *** consider switching models ***
- iv(irc) = 2
- go to 160
-c
-c *** accept step with decreased radius ***
-c
- 150 iv(irc) = 4
-c
-c *** set v(radfac) to fletcher*s decrease factor ***
-c
- 160 iv(xirc) = iv(irc)
- emax = v(gtstep) + v(fdif)
- v(radfac) = half * rfac1
- if (emax .lt. v(gtstep)) v(radfac) = rfac1 * dmax1(v(rdfcmn),
- 1 half * v(gtstep)/emax)
-c
-c *** do false convergence test ***
-c
- 170 if (v(reldx) .le. v(xftol)) go to 180
- iv(irc) = iv(xirc)
- if (v(f) .lt. v(f0)) go to 200
- go to 230
-c
- 180 iv(irc) = 12
- go to 240
-c
-c *** handle good function decrease ***
-c
- 190 if (v(fdif) .lt. (-v(tuner3) * v(gtstep))) go to 210
-c
-c *** increasing radius looks worthwhile. see if we just
-c *** recomputed step with a decreased radius or restored step
-c *** after recomputing it with a larger radius.
-c
- if (iv(radinc) .lt. 0) go to 210
- if (iv(restor) .eq. 1) go to 210
-c
-c *** we did not. try a longer step unless this was a newton
-c *** step.
-c
- v(radfac) = v(rdfcmx)
- gts = v(gtstep)
- if (v(fdif) .lt. (half/v(radfac) - one) * gts)
- 1 v(radfac) = dmax1(v(incfac), half*gts/(gts + v(fdif)))
- iv(irc) = 4
- if (v(stppar) .eq. zero) go to 230
- if (v(dst0) .ge. zero .and. (v(dst0) .lt. two*v(dstnrm)
- 1 .or. v(nreduc) .lt. onep2*v(fdif))) go to 230
-c *** step was not a newton step. recompute it with
-c *** a larger radius.
- iv(irc) = 5
- iv(radinc) = iv(radinc) + 1
-c
-c *** save values corresponding to good step ***
-c
- 200 v(flstgd) = v(f)
- iv(mlstgd) = iv(model)
- if (iv(restor) .ne. 1) iv(restor) = 2
- v(dstsav) = v(dstnrm)
- iv(nfgcal) = nfc
- v(plstgd) = v(preduc)
- v(gtslst) = v(gtstep)
- go to 230
-c
-c *** accept step with radius unchanged ***
-c
- 210 v(radfac) = one
- iv(irc) = 3
- go to 230
-c
-c *** come here for a restart after convergence ***
-c
- 220 iv(irc) = iv(xirc)
- if (v(dstsav) .ge. zero) go to 240
- iv(irc) = 12
- go to 240
-c
-c *** perform convergence tests ***
-c
- 230 iv(xirc) = iv(irc)
- 240 if (iv(restor) .eq. 1 .and. v(flstgd) .lt. v(f0)) iv(restor) = 3
- if (half * v(fdif) .gt. v(preduc)) go to 999
- emax = v(rfctol) * dabs(v(f0))
- emaxs = v(sctol) * dabs(v(f0))
- if (v(dstnrm) .gt. v(lmaxs) .and. v(preduc) .le. emaxs)
- 1 iv(irc) = 11
- if (v(dst0) .lt. zero) go to 250
- i = 0
- if ((v(nreduc) .gt. zero .and. v(nreduc) .le. emax) .or.
- 1 (v(nreduc) .eq. zero. and. v(preduc) .eq. zero)) i = 2
- if (v(stppar) .eq. zero .and. v(reldx) .le. v(xctol)
- 1 .and. goodx) i = i + 1
- if (i .gt. 0) iv(irc) = i + 6
-c
-c *** consider recomputing step of length v(lmaxs) for singular
-c *** convergence test.
-c
- 250 if (iv(irc) .gt. 5 .and. iv(irc) .ne. 12) go to 999
- if (v(dstnrm) .gt. v(lmaxs)) go to 260
- if (v(preduc) .ge. emaxs) go to 999
- if (v(dst0) .le. zero) go to 270
- if (half * v(dst0) .le. v(lmaxs)) go to 999
- go to 270
- 260 if (half * v(dstnrm) .le. v(lmaxs)) go to 999
- xmax = v(lmaxs) / v(dstnrm)
- if (xmax * (two - xmax) * v(preduc) .ge. emaxs) go to 999
- 270 if (v(nreduc) .lt. zero) go to 290
-c
-c *** recompute v(preduc) for use in singular convergence test ***
-c
- v(gtslst) = v(gtstep)
- v(dstsav) = v(dstnrm)
- if (iv(irc) .eq. 12) v(dstsav) = -v(dstsav)
- v(plstgd) = v(preduc)
- i = iv(restor)
- iv(restor) = 2
- if (i .eq. 3) iv(restor) = 0
- iv(irc) = 6
- go to 999
-c
-c *** perform singular convergence test with recomputed v(preduc) ***
-c
- 280 v(gtstep) = v(gtslst)
- v(dstnrm) = dabs(v(dstsav))
- iv(irc) = iv(xirc)
- if (v(dstsav) .le. zero) iv(irc) = 12
- v(nreduc) = -v(preduc)
- v(preduc) = v(plstgd)
- iv(restor) = 3
- 290 if (-v(nreduc) .le. v(sctol) * dabs(v(f0))) iv(irc) = 11
-c
- 999 return
-c
-c *** last card of assst follows ***
- end
- subroutine deflt(alg, iv, liv, lv, v)
-c
-c *** supply ***sol (version 2.3) default values to iv and v ***
-c
-c *** alg = 1 means regression constants.
-c *** alg = 2 means general unconstrained optimization constants.
-c
- integer liv, l
- integer alg, iv(liv)
- double precision v(lv)
-c
- external imdcon, vdflt
- integer imdcon
-c imdcon... returns machine-dependent integer constants.
-c vdflt.... provides default values to v.
-c
- integer miv, m
- integer miniv(2), minv(2)
-c
-c *** subscripts for iv ***
-c
- integer algsav, covprt, covreq, dtype, hc, ierr, inith, inits,
- 1 ipivot, ivneed, lastiv, lastv, lmat, mxfcal, mxiter,
- 2 nfcov, ngcov, nvdflt, outlev, parprt, parsav, perm,
- 3 prunit, qrtyp, rdreq, rmat, solprt, statpr, vneed,
- 4 vsave, x0prt
-c
-c *** iv subscript values ***
-c
-c/6
-c data algsav/51/, covprt/14/, covreq/15/, dtype/16/, hc/71/,
-c 1 ierr/75/, inith/25/, inits/25/, ipivot/76/, ivneed/3/,
-c 2 lastiv/44/, lastv/45/, lmat/42/, mxfcal/17/, mxiter/18/,
-c 3 nfcov/52/, ngcov/53/, nvdflt/50/, outlev/19/, parprt/20/,
-c 4 parsav/49/, perm/58/, prunit/21/, qrtyp/80/, rdreq/57/,
-c 5 rmat/78/, solprt/22/, statpr/23/, vneed/4/, vsave/60/,
-c 6 x0prt/24/
-c/7
- parameter (algsav=51, covprt=14, covreq=15, dtype=16, hc=71,
- 1 ierr=75, inith=25, inits=25, ipivot=76, ivneed=3,
- 2 lastiv=44, lastv=45, lmat=42, mxfcal=17, mxiter=18,
- 3 nfcov=52, ngcov=53, nvdflt=50, outlev=19, parprt=20,
- 4 parsav=49, perm=58, prunit=21, qrtyp=80, rdreq=57,
- 5 rmat=78, solprt=22, statpr=23, vneed=4, vsave=60,
- 6 x0prt=24)
-c/
- data miniv(1)/80/, miniv(2)/59/, minv(1)/98/, minv(2)/71/
-c
-c------------------------------- body --------------------------------
-c
- if (alg .lt. 1 .or. alg .gt. 2) go to 40
- miv = miniv(alg)
- if (liv .lt. miv) go to 20
- mv = minv(alg)
- if (lv .lt. mv) go to 30
- call vdflt(alg, lv, v)
- iv(1) = 12
- iv(algsav) = alg
- 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(prunit) = imdcon(1)
- iv(solprt) = 1
- iv(statpr) = 1
- iv(vneed) = 0
- iv(x0prt) = 1
-c
- if (alg .ge. 2) go to 10
-c
-c *** regression values
-c
- 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(parsav) = 67
- iv(qrtyp) = 1
- iv(rdreq) = 3
- iv(rmat) = 0
- iv(vsave) = 58
- go to 999
-c
-c *** general optimization values
-c
- 10 iv(dtype) = 0
- iv(inith) = 1
- iv(nfcov) = 0
- iv(ngcov) = 0
- iv(nvdflt) = 25
- iv(parsav) = 47
- go to 999
-c
- 20 iv(1) = 15
- go to 999
-c
- 30 iv(1) = 16
- go to 999
-c
- 40 iv(1) = 67
-c
- 999 return
-c *** last card of deflt follows ***
- end
- double precision function dotprd(p, x, y)
-c
-c *** return the inner product of the p-vectors x and y. ***
-c
- integer p
- double precision x(p), y(p)
-c
- integer i
- double precision one, sqteta, t, zero
-c/+
- double precision dmax1, dabs
-c/
- external rmdcon
- double precision rmdcon
-c
-c *** rmdcon(2) returns a machine-dependent constant, sqteta, which
-c *** is slightly larger than the smallest positive number that
-c *** can be squared without underflowing.
-c
-c/6
-c data one/1.d+0/, sqteta/0.d+0/, zero/0.d+0/
-c/7
- parameter (one=1.d+0, zero=0.d+0)
- data sqteta/0.d+0/
-c/
-c
- dotprd = zero
- if (p .le. 0) go to 999
-crc if (sqteta .eq. zero) sqteta = rmdcon(2)
- do 20 i = 1, p
-crc t = dmax1(dabs(x(i)), dabs(y(i)))
-crc if (t .gt. one) go to 10
-crc if (t .lt. sqteta) go to 20
-crc t = (x(i)/sqteta)*y(i)
-crc if (dabs(t) .lt. sqteta) go to 20
- 10 dotprd = dotprd + x(i)*y(i)
- 20 continue
-c
- 999 return
-c *** last card of dotprd follows ***
- end
- subroutine itsum(d, g, iv, liv, lv, p, v, x)
-c
-c *** print iteration summary for ***sol (version 2.3) ***
-c
-c *** parameter declarations ***
-c
- integer liv, lv, p
- integer iv(liv)
- double precision d(p), g(p), v(lv), x(p)
-c
-c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-c
-c *** local variables ***
-c
- integer alg, i, iv1, m, nf, ng, ol, pu
-c/6
-c real model1(6), model2(6)
-c/7
- character*4 model1(6), model2(6)
-c/
- double precision nreldf, oldf, preldf, reldf, zero
-c
-c *** intrinsic functions ***
-c/+
- integer iabs
- double precision dabs, dmax1
-c/
-c *** no external functions or subroutines ***
-c
-c *** subscripts for iv and v ***
-c
- integer algsav, dstnrm, f, fdif, f0, needhd, nfcall, nfcov, ngcov,
- 1 ngcall, niter, nreduc, outlev, preduc, prntit, prunit,
- 2 reldx, solprt, statpr, stppar, sused, x0prt
-c
-c *** iv subscript values ***
-c
-c/6
-c data algsav/51/, needhd/36/, nfcall/6/, nfcov/52/, ngcall/30/,
-c 1 ngcov/53/, niter/31/, outlev/19/, prntit/39/, prunit/21/,
-c 2 solprt/22/, statpr/23/, sused/64/, x0prt/24/
-c/7
- parameter (algsav=51, needhd=36, nfcall=6, nfcov=52, ngcall=30,
- 1 ngcov=53, niter=31, outlev=19, prntit=39, prunit=21,
- 2 solprt=22, statpr=23, sused=64, x0prt=24)
-c/
-c
-c *** v subscript values ***
-c
-c/6
-c data dstnrm/2/, f/10/, f0/13/, fdif/11/, nreduc/6/, preduc/7/,
-c 1 reldx/17/, stppar/5/
-c/7
- parameter (dstnrm=2, f=10, f0=13, fdif=11, nreduc=6, preduc=7,
- 1 reldx=17, stppar=5)
-c/
-c
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c/6
-c data model1(1)/4h /, model1(2)/4h /, model1(3)/4h /,
-c 1 model1(4)/4h /, model1(5)/4h g /, model1(6)/4h s /,
-c 2 model2(1)/4h g /, model2(2)/4h s /, model2(3)/4hg-s /,
-c 3 model2(4)/4hs-g /, model2(5)/4h-s-g/, model2(6)/4h-g-s/
-c/7
- data model1/' ',' ',' ',' ',' g ',' s '/,
- 1 model2/' g ',' s ','g-s ','s-g ','-s-g','-g-s'/
-c/
-c
-c------------------------------- body --------------------------------
-c
- pu = iv(prunit)
- if (pu .eq. 0) go to 999
- iv1 = iv(1)
- if (iv1 .gt. 62) iv1 = iv1 - 51
- ol = iv(outlev)
- alg = iv(algsav)
- if (iv1 .lt. 2 .or. iv1 .gt. 15) go to 370
- if (iv1 .ge. 12) go to 120
- if (iv1 .eq. 2 .and. iv(niter) .eq. 0) go to 390
- if (ol .eq. 0) go to 120
- if (iv1 .ge. 10 .and. iv(prntit) .eq. 0) go to 120
- if (iv1 .gt. 2) go to 10
- iv(prntit) = iv(prntit) + 1
- if (iv(prntit) .lt. iabs(ol)) go to 999
- 10 nf = iv(nfcall) - iabs(iv(nfcov))
- iv(prntit) = 0
- reldf = zero
- preldf = zero
- oldf = dmax1(dabs(v(f0)), dabs(v(f)))
- if (oldf .le. zero) go to 20
- reldf = v(fdif) / oldf
- preldf = v(preduc) / oldf
- 20 if (ol .gt. 0) go to 60
-c
-c *** print short summary line ***
-c
- if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,30)
- 30 format(/10h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,
- 1 2x,13hmodel stppar)
- if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,40)
- 40 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,
- 1 3x,6hstppar)
- iv(needhd) = 0
- if (alg .eq. 2) go to 50
- m = iv(sused)
- write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),
- 1 model1(m), model2(m), v(stppar)
- go to 120
-c
- 50 write(pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx),
- 1 v(stppar)
- go to 120
-c
-c *** print long summary line ***
-c
- 60 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,70)
- 70 format(/11h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,
- 1 2x,13hmodel stppar,2x,6hd*step,2x,7hnpreldf)
- if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,80)
- 80 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,
- 1 3x,6hstppar,3x,6hd*step,3x,7hnpreldf)
- iv(needhd) = 0
- nreldf = zero
- if (oldf .gt. zero) nreldf = v(nreduc) / oldf
- if (alg .eq. 2) go to 90
- m = iv(sused)
- write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),
- 1 model1(m), model2(m), v(stppar), v(dstnrm), nreldf
- go to 120
-c
- 90 write(pu,110) iv(niter), nf, v(f), reldf, preldf,
- 1 v(reldx), v(stppar), v(dstnrm), nreldf
- 100 format(i6,i5,d10.3,2d9.2,d8.1,a3,a4,2d8.1,d9.2)
- 110 format(i6,i5,d11.3,2d10.2,3d9.1,d10.2)
-c
- 120 if (iv(statpr) .lt. 0) go to 430
- go to (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,
- 1 330, 350, 520), iv1
-c
- 130 write(pu,140)
- 140 format(/26h ***** x-convergence *****)
- go to 430
-c
- 150 write(pu,160)
- 160 format(/42h ***** relative function convergence *****)
- go to 430
-c
- 170 write(pu,180)
- 180 format(/49h ***** x- and relative function convergence *****)
- go to 430
-c
- 190 write(pu,200)
- 200 format(/42h ***** absolute function convergence *****)
- go to 430
-c
- 210 write(pu,220)
- 220 format(/33h ***** singular convergence *****)
- go to 430
-c
- 230 write(pu,240)
- 240 format(/30h ***** false convergence *****)
- go to 430
-c
- 250 write(pu,260)
- 260 format(/38h ***** function evaluation limit *****)
- go to 430
-c
- 270 write(pu,280)
- 280 format(/28h ***** iteration limit *****)
- go to 430
-c
- 290 write(pu,300)
- 300 format(/18h ***** stopx *****)
- go to 430
-c
- 310 write(pu,320)
- 320 format(/44h ***** initial f(x) cannot be computed *****)
-c
- go to 390
-c
- 330 write(pu,340)
- 340 format(/37h ***** bad parameters to assess *****)
- go to 999
-c
- 350 write(pu,360)
- 360 format(/43h ***** gradient could not be computed *****)
- if (iv(niter) .gt. 0) go to 480
- go to 390
-c
- 370 write(pu,380) iv(1)
- 380 format(/14h ***** iv(1) =,i5,6h *****)
- go to 999
-c
-c *** initial call on itsum ***
-c
- 390 if (iv(x0prt) .ne. 0) write(pu,400) (i, x(i), d(i), i = 1, p)
- 400 format(/23h i initial x(i),8x,4hd(i)//(1x,i5,d17.6,d14.3))
-c *** the following are to avoid undefined variables when the
-c *** function evaluation limit is 1...
- v(dstnrm) = zero
- v(fdif) = zero
- v(nreduc) = zero
- v(preduc) = zero
- v(reldx) = zero
- if (iv1 .ge. 12) go to 999
- iv(needhd) = 0
- iv(prntit) = 0
- if (ol .eq. 0) go to 999
- if (ol .lt. 0 .and. alg .eq. 1) write(pu,30)
- if (ol .lt. 0 .and. alg .eq. 2) write(pu,40)
- if (ol .gt. 0 .and. alg .eq. 1) write(pu,70)
- if (ol .gt. 0 .and. alg .eq. 2) write(pu,80)
- if (alg .eq. 1) write(pu,410) v(f)
- if (alg .eq. 2) write(pu,420) v(f)
- 410 format(/11h 0 1,d10.3)
-c365 format(/11h 0 1,e11.3)
- 420 format(/11h 0 1,d11.3)
- go to 999
-c
-c *** print various information requested on solution ***
-c
- 430 iv(needhd) = 1
- if (iv(statpr) .eq. 0) go to 480
- oldf = dmax1(dabs(v(f0)), dabs(v(f)))
- preldf = zero
- nreldf = zero
- if (oldf .le. 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(/9h function,d17.6,8h reldx,d17.3/12h func. evals,
- 1 i8,9x,11hgrad. evals,i8/7h preldf,d16.3,6x,7hnpreldf,d15.3)
-c
- if (iv(nfcov) .gt. 0) write(pu,460) iv(nfcov)
- 460 format(/1x,i4,50h extra func. evals for covariance and diagnost
- 1ics.)
- if (iv(ngcov) .gt. 0) write(pu,470) iv(ngcov)
- 470 format(1x,i4,50h extra grad. evals for covariance and diagnosti
- 1cs.)
-c
- 480 if (iv(solprt) .eq. 0) go to 999
- iv(needhd) = 1
- write(pu,490)
- 490 format(/22h i final x(i),8x,4hd(i),10x,4hg(i)/)
- do 500 i = 1, p
- write(pu,510) i, x(i), d(i), g(i)
- 500 continue
- 510 format(1x,i5,d16.6,2d14.3)
- go to 999
-c
- 520 write(pu,530)
- 530 format(/24h inconsistent dimensions)
- 999 return
-c *** last card of itsum follows ***
- end
- subroutine litvmu(n, x, l, y)
-c
-c *** solve (l**t)*x = y, where l is an n x n lower triangular
-c *** matrix stored compactly by rows. x and y may occupy the same
-c *** storage. ***
-c
- integer n
-cal double precision x(n), l(1), y(n)
- double precision x(n), l(n*(n+1)/2), y(n)
- integer i, ii, ij, im1, i0, j, np1
- double precision xi, zero
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c
- do 10 i = 1, n
- 10 x(i) = y(i)
- np1 = n + 1
- i0 = n*(n+1)/2
- do 30 ii = 1, n
- i = np1 - ii
- xi = x(i)/l(i0)
- x(i) = xi
- if (i .le. 1) go to 999
- i0 = i0 - i
- if (xi .eq. zero) go to 30
- im1 = i - 1
- do 20 j = 1, im1
- ij = i0 + j
- x(j) = x(j) - xi*l(ij)
- 20 continue
- 30 continue
- 999 return
-c *** last card of litvmu follows ***
- end
- subroutine livmul(n, x, l, y)
-c
-c *** solve l*x = y, where l is an n x n lower triangular
-c *** matrix stored compactly by rows. x and y may occupy the same
-c *** storage. ***
-c
- integer n
-cal double precision x(n), l(1), y(n)
- double precision x(n), l(n*(n+1)/2), y(n)
- external dotprd
- double precision dotprd
- integer i, j, k
- double precision t, zero
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c
- do 10 k = 1, n
- if (y(k) .ne. zero) go to 20
- x(k) = zero
- 10 continue
- go to 999
- 20 j = k*(k+1)/2
- x(k) = y(k) / l(j)
- if (k .ge. n) go to 999
- k = k + 1
- do 30 i = k, n
- t = dotprd(i-1, l(j+1), x)
- j = j + i
- x(i) = (y(i) - t)/l(j)
- 30 continue
- 999 return
-c *** last card of livmul follows ***
- end
- subroutine parck(alg, d, iv, liv, lv, n, v)
-c
-c *** check ***sol (version 2.3) parameters, print changed values ***
-c
-c *** alg = 1 for regression, alg = 2 for general unconstrained opt.
-c
- integer alg, liv, lv, n
- integer iv(liv)
- double precision d(n), v(lv)
-c
- external rmdcon, vcopy, vdflt
- double precision rmdcon
-c rmdcon -- returns machine-dependent constants.
-c vcopy -- copies one vector to another.
-c vdflt -- supplies default parameter values to v alone.
-c/+
- integer max0
-c/
-c
-c *** local variables ***
-c
- integer i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu
- integer ijmp, jlim(2), miniv(2), ndflt(2)
-c/6
-c integer varnm(2), sh(2)
-c real cngd(3), dflt(3), vn(2,34), which(3)
-c/7
- character*1 varnm(2), sh(2)
- character*4 cngd(3), dflt(3), vn(2,34), which(3)
-c/
- double precision big, machep, tiny, vk, vm(34), vx(34), zero
-c
-c *** iv and v subscripts ***
-c
- integer algsav, dinit, dtype, dtype0, epslon, inits, ivneed,
- 1 lastiv, lastv, lmat, nextiv, nextv, nvdflt, oldn,
- 2 parprt, parsav, perm, prunit, vneed
-c
-c
-c/6
-c data algsav/51/, dinit/38/, dtype/16/, dtype0/54/, epslon/19/,
-c 1 inits/25/, ivneed/3/, lastiv/44/, lastv/45/, lmat/42/,
-c 2 nextiv/46/, nextv/47/, nvdflt/50/, oldn/38/, parprt/20/,
-c 3 parsav/49/, perm/58/, prunit/21/, vneed/4/
-c/7
- parameter (algsav=51, dinit=38, dtype=16, dtype0=54, epslon=19,
- 1 inits=25, ivneed=3, lastiv=44, lastv=45, lmat=42,
- 2 nextiv=46, nextv=47, nvdflt=50, oldn=38, parprt=20,
- 3 parsav=49, perm=58, prunit=21, vneed=4)
- save big, machep, tiny
-c/
-c
- data big/0.d+0/, machep/-1.d+0/, tiny/1.d+0/, zero/0.d+0/
-c/6
-c data vn(1,1),vn(2,1)/4hepsl,4hon../
-c data vn(1,2),vn(2,2)/4hphmn,4hfc../
-c data vn(1,3),vn(2,3)/4hphmx,4hfc../
-c data vn(1,4),vn(2,4)/4hdecf,4hac../
-c data vn(1,5),vn(2,5)/4hincf,4hac../
-c data vn(1,6),vn(2,6)/4hrdfc,4hmn../
-c data vn(1,7),vn(2,7)/4hrdfc,4hmx../
-c data vn(1,8),vn(2,8)/4htune,4hr1../
-c data vn(1,9),vn(2,9)/4htune,4hr2../
-c data vn(1,10),vn(2,10)/4htune,4hr3../
-c data vn(1,11),vn(2,11)/4htune,4hr4../
-c data vn(1,12),vn(2,12)/4htune,4hr5../
-c data vn(1,13),vn(2,13)/4hafct,4hol../
-c data vn(1,14),vn(2,14)/4hrfct,4hol../
-c data vn(1,15),vn(2,15)/4hxcto,4hl.../
-c data vn(1,16),vn(2,16)/4hxfto,4hl.../
-c data vn(1,17),vn(2,17)/4hlmax,4h0.../
-c data vn(1,18),vn(2,18)/4hlmax,4hs.../
-c data vn(1,19),vn(2,19)/4hscto,4hl.../
-c data vn(1,20),vn(2,20)/4hdini,4ht.../
-c data vn(1,21),vn(2,21)/4hdtin,4hit../
-c data vn(1,22),vn(2,22)/4hd0in,4hit../
-c data vn(1,23),vn(2,23)/4hdfac,4h..../
-c data vn(1,24),vn(2,24)/4hdltf,4hdc../
-c data vn(1,25),vn(2,25)/4hdltf,4hdj../
-c data vn(1,26),vn(2,26)/4hdelt,4ha0../
-c data vn(1,27),vn(2,27)/4hfuzz,4h..../
-c data vn(1,28),vn(2,28)/4hrlim,4hit../
-c data vn(1,29),vn(2,29)/4hcosm,4hin../
-c data vn(1,30),vn(2,30)/4hhube,4hrc../
-c data vn(1,31),vn(2,31)/4hrspt,4hol../
-c data vn(1,32),vn(2,32)/4hsigm,4hin../
-c data vn(1,33),vn(2,33)/4heta0,4h..../
-c data vn(1,34),vn(2,34)/4hbias,4h..../
-c/7
- data vn(1,1),vn(2,1)/'epsl','on..'/
- data vn(1,2),vn(2,2)/'phmn','fc..'/
- data vn(1,3),vn(2,3)/'phmx','fc..'/
- data vn(1,4),vn(2,4)/'decf','ac..'/
- data vn(1,5),vn(2,5)/'incf','ac..'/
- data vn(1,6),vn(2,6)/'rdfc','mn..'/
- data vn(1,7),vn(2,7)/'rdfc','mx..'/
- data vn(1,8),vn(2,8)/'tune','r1..'/
- data vn(1,9),vn(2,9)/'tune','r2..'/
- data vn(1,10),vn(2,10)/'tune','r3..'/
- data vn(1,11),vn(2,11)/'tune','r4..'/
- data vn(1,12),vn(2,12)/'tune','r5..'/
- data vn(1,13),vn(2,13)/'afct','ol..'/
- data vn(1,14),vn(2,14)/'rfct','ol..'/
- data vn(1,15),vn(2,15)/'xcto','l...'/
- data vn(1,16),vn(2,16)/'xfto','l...'/
- data vn(1,17),vn(2,17)/'lmax','0...'/
- data vn(1,18),vn(2,18)/'lmax','s...'/
- data vn(1,19),vn(2,19)/'scto','l...'/
- data vn(1,20),vn(2,20)/'dini','t...'/
- data vn(1,21),vn(2,21)/'dtin','it..'/
- data vn(1,22),vn(2,22)/'d0in','it..'/
- data vn(1,23),vn(2,23)/'dfac','....'/
- data vn(1,24),vn(2,24)/'dltf','dc..'/
- data vn(1,25),vn(2,25)/'dltf','dj..'/
- data vn(1,26),vn(2,26)/'delt','a0..'/
- data vn(1,27),vn(2,27)/'fuzz','....'/
- data vn(1,28),vn(2,28)/'rlim','it..'/
- data vn(1,29),vn(2,29)/'cosm','in..'/
- data vn(1,30),vn(2,30)/'hube','rc..'/
- data vn(1,31),vn(2,31)/'rspt','ol..'/
- data vn(1,32),vn(2,32)/'sigm','in..'/
- data vn(1,33),vn(2,33)/'eta0','....'/
- data vn(1,34),vn(2,34)/'bias','....'/
-c/
-c
- data vm(1)/1.0d-3/, vm(2)/-0.99d+0/, vm(3)/1.0d-3/, vm(4)/1.0d-2/,
- 1 vm(5)/1.2d+0/, vm(6)/1.d-2/, vm(7)/1.2d+0/, vm(8)/0.d+0/,
- 2 vm(9)/0.d+0/, vm(10)/1.d-3/, vm(11)/-1.d+0/, vm(13)/0.d+0/,
- 3 vm(15)/0.d+0/, vm(16)/0.d+0/, vm(19)/0.d+0/, vm(20)/-10.d+0/,
- 4 vm(21)/0.d+0/, vm(22)/0.d+0/, vm(23)/0.d+0/, vm(27)/1.01d+0/,
- 5 vm(28)/1.d+10/, vm(30)/0.d+0/, vm(31)/0.d+0/, vm(32)/0.d+0/,
- 6 vm(34)/0.d+0/
- data vx(1)/0.9d+0/, vx(2)/-1.d-3/, vx(3)/1.d+1/, vx(4)/0.8d+0/,
- 1 vx(5)/1.d+2/, vx(6)/0.8d+0/, vx(7)/1.d+2/, vx(8)/0.5d+0/,
- 2 vx(9)/0.5d+0/, vx(10)/1.d+0/, vx(11)/1.d+0/, vx(14)/0.1d+0/,
- 3 vx(15)/1.d+0/, vx(16)/1.d+0/, vx(19)/1.d+0/, vx(23)/1.d+0/,
- 4 vx(24)/1.d+0/, vx(25)/1.d+0/, vx(26)/1.d+0/, vx(27)/1.d+10/,
- 5 vx(29)/1.d+0/, vx(31)/1.d+0/, vx(32)/1.d+0/, vx(33)/1.d+0/,
- 6 vx(34)/1.d+0/
-c
-c/6
-c data varnm(1)/1hp/, varnm(2)/1hn/, sh(1)/1hs/, sh(2)/1hh/
-c data cngd(1),cngd(2),cngd(3)/4h---c,4hhang,4hed v/,
-c 1 dflt(1),dflt(2),dflt(3)/4hnond,4hefau,4hlt v/
-c/7
- data varnm(1)/'p'/, varnm(2)/'n'/, sh(1)/'s'/, sh(2)/'h'/
- data cngd(1),cngd(2),cngd(3)/'---c','hang','ed v'/,
- 1 dflt(1),dflt(2),dflt(3)/'nond','efau','lt v'/
-c/
- data ijmp/33/, jlim(1)/0/, jlim(2)/24/, ndflt(1)/32/, ndflt(2)/25/
- data miniv(1)/80/, miniv(2)/59/
-c
-c............................... body ................................
-c
- pu = 0
- if (prunit .le. liv) pu = iv(prunit)
- if (alg .lt. 1 .or. alg .gt. 2) go to 340
- if (iv(1) .eq. 0) call deflt(alg, iv, liv, lv, v)
- iv1 = iv(1)
- if (iv1 .ne. 13 .and. iv1 .ne. 12) go to 10
- miv1 = miniv(alg)
- if (perm .le. liv) miv1 = max0(miv1, iv(perm) - 1)
- if (ivneed .le. liv) miv2 = miv1 + max0(iv(ivneed), 0)
- if (lastiv .le. liv) iv(lastiv) = miv2
- if (liv .lt. miv1) go to 300
- iv(ivneed) = 0
- iv(lastv) = max0(iv(vneed), 0) + iv(lmat) - 1
- iv(vneed) = 0
- if (liv .lt. miv2) go to 300
- if (lv .lt. iv(lastv)) go to 320
- 10 if (alg .eq. iv(algsav)) go to 30
- if (pu .ne. 0) write(pu,20) alg, iv(algsav)
- 20 format(/39h the first parameter to deflt should be,i3,
- 1 12h rather than,i3)
- iv(1) = 82
- go to 999
- 30 if (iv1 .lt. 12 .or. iv1 .gt. 14) go to 60
- if (n .ge. 1) go to 50
- iv(1) = 81
- if (pu .eq. 0) go to 999
- write(pu,40) varnm(alg), n
- 40 format(/8h /// bad,a1,2h =,i5)
- go to 999
- 50 if (iv1 .ne. 14) iv(nextiv) = iv(perm)
- if (iv1 .ne. 14) iv(nextv) = iv(lmat)
- if (iv1 .eq. 13) go to 999
- k = iv(parsav) - epslon
- call vdflt(alg, lv-k, v(k+1))
- iv(dtype0) = 2 - alg
- iv(oldn) = n
- which(1) = dflt(1)
- which(2) = dflt(2)
- which(3) = dflt(3)
- go to 110
- 60 if (n .eq. iv(oldn)) go to 80
- iv(1) = 17
- if (pu .eq. 0) go to 999
- write(pu,70) varnm(alg), iv(oldn), n
- 70 format(/5h /// ,1a1,14h changed from ,i5,4h to ,i5)
- go to 999
-c
- 80 if (iv1 .le. 11 .and. iv1 .ge. 1) go to 100
- iv(1) = 80
- if (pu .ne. 0) write(pu,90) iv1
- 90 format(/13h /// iv(1) =,i5,28h should be between 0 and 14.)
- go to 999
-c
- 100 which(1) = cngd(1)
- which(2) = cngd(2)
- which(3) = cngd(3)
-c
- 110 if (iv1 .eq. 14) iv1 = 12
- if (big .gt. tiny) go to 120
- tiny = rmdcon(1)
- machep = rmdcon(3)
- big = rmdcon(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) = rmdcon(5)
- vm(29) = machep
- vx(30) = big
- vm(33) = machep
- 120 m = 0
- i = 1
- j = jlim(alg)
- k = epslon
- ndfalt = ndflt(alg)
- do 150 l = 1, ndfalt
- vk = v(k)
- if (vk .ge. vm(i) .and. vk .le. vx(i)) go to 140
- m = k
- if (pu .ne. 0) write(pu,130) vn(1,i), vn(2,i), k, vk,
- 1 vm(i), vx(i)
- 130 format(/6h /// ,2a4,5h.. v(,i2,3h) =,d11.3,7h should,
- 1 11h be between,d11.3,4h and,d11.3)
- 140 k = k + 1
- i = i + 1
- if (i .eq. j) i = ijmp
- 150 continue
-c
- if (iv(nvdflt) .eq. ndfalt) go to 170
- iv(1) = 51
- if (pu .eq. 0) go to 999
- write(pu,160) iv(nvdflt), ndfalt
- 160 format(/13h iv(nvdflt) =,i5,13h rather than ,i5)
- go to 999
- 170 if ((iv(dtype) .gt. 0 .or. v(dinit) .gt. zero) .and. iv1 .eq. 12)
- 1 go to 200
- do 190 i = 1, n
- if (d(i) .gt. zero) go to 190
- m = 18
- if (pu .ne. 0) write(pu,180) i, d(i)
- 180 format(/8h /// d(,i3,3h) =,d11.3,19h should be positive)
- 190 continue
- 200 if (m .eq. 0) go to 210
- iv(1) = m
- go to 999
-c
- 210 if (pu .eq. 0 .or. iv(parprt) .eq. 0) go to 999
- if (iv1 .ne. 12 .or. iv(inits) .eq. alg-1) go to 230
- m = 1
- write(pu,220) sh(alg), iv(inits)
- 220 format(/22h nondefault values..../5h init,a1,14h..... iv(25) =,
- 1 i3)
- 230 if (iv(dtype) .eq. iv(dtype0)) go to 250
- if (m .eq. 0) write(pu,260) which
- m = 1
- write(pu,240) iv(dtype)
- 240 format(20h dtype..... iv(16) =,i3)
- 250 i = 1
- j = jlim(alg)
- k = epslon
- l = iv(parsav)
- ndfalt = ndflt(alg)
- do 290 ii = 1, ndfalt
- if (v(k) .eq. v(l)) go to 280
- if (m .eq. 0) write(pu,260) which
- 260 format(/1h ,3a4,9halues..../)
- m = 1
- write(pu,270) vn(1,i), vn(2,i), k, v(k)
- 270 format(1x,2a4,5h.. v(,i2,3h) =,d15.7)
- 280 k = k + 1
- l = l + 1
- i = i + 1
- if (i .eq. j) i = ijmp
- 290 continue
-c
- iv(dtype0) = iv(dtype)
- parsv1 = iv(parsav)
- call vcopy(iv(nvdflt), v(parsv1), v(epslon))
- go to 999
-c
- 300 iv(1) = 15
- if (pu .eq. 0) go to 999
- write(pu,310) liv, miv2
- 310 format(/10h /// liv =,i5,17h must be at least,i5)
- if (liv .lt. miv1) go to 999
- if (lv .lt. iv(lastv)) go to 320
- go to 999
-c
- 320 iv(1) = 16
- if (pu .eq. 0) go to 999
- write(pu,330) lv, iv(lastv)
- 330 format(/9h /// lv =,i5,17h must be at least,i5)
- go to 999
-c
- 340 iv(1) = 67
- if (pu .eq. 0) go to 999
- write(pu,350) alg
- 350 format(/10h /// alg =,i5,15h must be 1 or 2)
-c
- 999 return
-c *** last card of parck follows ***
- end
- double precision function reldst(p, d, x, x0)
-c
-c *** compute and return relative difference between x and x0 ***
-c *** nl2sol version 2.2 ***
-c
- integer p
- double precision d(p), x(p), x0(p)
-c/+
- double precision dabs
-c/
- integer i
- double precision emax, t, xmax, zero
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c
- emax = zero
- xmax = zero
- do 10 i = 1, p
- t = dabs(d(i) * (x(i) - x0(i)))
- if (emax .lt. t) emax = t
- t = d(i) * (dabs(x(i)) + dabs(x0(i)))
- if (xmax .lt. t) xmax = t
- 10 continue
- reldst = zero
- if (xmax .gt. zero) reldst = emax / xmax
- 999 return
-c *** last card of reldst follows ***
- end
-c logical function stopx(idummy)
-c *****parameters...
-c integer idummy
-c
-c ..................................................................
-c
-c *****purpose...
-c this function may serve as the stopx (asynchronous interruption)
-c function for the nl2sol (nonlinear least-squares) package at
-c those installations which do not wish to implement a
-c dynamic stopx.
-c
-c *****algorithm notes...
-c at installations where the nl2sol system is used
-c interactively, this dummy stopx should be replaced by a
-c function that returns .true. if and only if the interrupt
-c (break) key has been pressed since the last call on stopx.
-c
-c ..................................................................
-c
-c stopx = .false.
-c return
-c end
- subroutine vaxpy(p, w, a, x, y)
-c
-c *** set w = a*x + y -- w, x, y = p-vectors, a = scalar ***
-c
- integer p
- double precision a, w(p), x(p), y(p)
-c
- integer i
-c
- do 10 i = 1, p
- 10 w(i) = a*x(i) + y(i)
- return
- end
- subroutine vcopy(p, y, x)
-c
-c *** set y = x, where x and y are p-vectors ***
-c
- integer p
- double precision x(p), y(p)
-c
- integer i
-c
- do 10 i = 1, p
- 10 y(i) = x(i)
- return
- end
- subroutine vdflt(alg, lv, v)
-c
-c *** supply ***sol (version 2.3) default values to v ***
-c
-c *** alg = 1 means regression constants.
-c *** alg = 2 means general unconstrained optimization constants.
-c
- integer alg, l
- double precision v(lv)
-c/+
- double precision dmax1
-c/
- external rmdcon
- double precision rmdcon
-c rmdcon... returns machine-dependent constants
-c
- double precision machep, mepcrt, one, sqteps, three
-c
-c *** subscripts for v ***
-c
- integer afctol, bias, cosmin, decfac, delta0, dfac, dinit, dltfdc,
- 1 dltfdj, dtinit, d0init, epslon, eta0, fuzz, huberc,
- 2 incfac, lmax0, lmaxs, phmnfc, phmxfc, rdfcmn, rdfcmx,
- 3 rfctol, rlimit, rsptol, sctol, sigmin, tuner1, tuner2,
- 4 tuner3, tuner4, tuner5, xctol, xftol
-c
-c/6
-c data one/1.d+0/, three/3.d+0/
-c/7
- parameter (one=1.d+0, three=3.d+0)
-c/
-c
-c *** v subscript values ***
-c
-c/6
-c data afctol/31/, bias/43/, cosmin/47/, decfac/22/, delta0/44/,
-c 1 dfac/41/, dinit/38/, dltfdc/42/, dltfdj/43/, dtinit/39/,
-c 2 d0init/40/, epslon/19/, eta0/42/, fuzz/45/, huberc/48/,
-c 3 incfac/23/, lmax0/35/, lmaxs/36/, phmnfc/20/, phmxfc/21/,
-c 4 rdfcmn/24/, rdfcmx/25/, rfctol/32/, rlimit/46/, rsptol/49/,
-c 5 sctol/37/, sigmin/50/, tuner1/26/, tuner2/27/, tuner3/28/,
-c 6 tuner4/29/, tuner5/30/, xctol/33/, xftol/34/
-c/7
- parameter (afctol=31, bias=43, cosmin=47, decfac=22, delta0=44,
- 1 dfac=41, dinit=38, dltfdc=42, dltfdj=43, dtinit=39,
- 2 d0init=40, epslon=19, eta0=42, fuzz=45, huberc=48,
- 3 incfac=23, lmax0=35, lmaxs=36, phmnfc=20, phmxfc=21,
- 4 rdfcmn=24, rdfcmx=25, rfctol=32, rlimit=46, rsptol=49,
- 5 sctol=37, sigmin=50, tuner1=26, tuner2=27, tuner3=28,
- 6 tuner4=29, tuner5=30, xctol=33, xftol=34)
-c/
-c
-c------------------------------- body --------------------------------
-c
- machep = rmdcon(3)
- v(afctol) = 1.d-20
- if (machep .gt. 1.d-10) v(afctol) = machep**2
- v(decfac) = 0.5d+0
- sqteps = rmdcon(4)
- v(dfac) = 0.6d+0
- v(delta0) = sqteps
- v(dtinit) = 1.d-6
- mepcrt = machep ** (one/three)
- v(d0init) = 1.d+0
- v(epslon) = 0.1d+0
- v(incfac) = 2.d+0
- v(lmax0) = 1.d+0
- v(lmaxs) = 1.d+0
- v(phmnfc) = -0.1d+0
- v(phmxfc) = 0.1d+0
- v(rdfcmn) = 0.1d+0
- v(rdfcmx) = 4.d+0
- v(rfctol) = dmax1(1.d-10, mepcrt**2)
- v(sctol) = v(rfctol)
- v(tuner1) = 0.1d+0
- v(tuner2) = 1.d-4
- v(tuner3) = 0.75d+0
- v(tuner4) = 0.5d+0
- v(tuner5) = 0.75d+0
- v(xctol) = sqteps
- v(xftol) = 1.d+2 * machep
-c
- if (alg .ge. 2) go to 10
-c
-c *** regression values
-c
- v(cosmin) = dmax1(1.d-6, 1.d+2 * machep)
- v(dinit) = 0.d+0
- v(dltfdc) = mepcrt
- v(dltfdj) = sqteps
- v(fuzz) = 1.5d+0
- v(huberc) = 0.7d+0
- v(rlimit) = rmdcon(5)
- v(rsptol) = 1.d-3
- v(sigmin) = 1.d-4
- go to 999
-c
-c *** general optimization values
-c
- 10 v(bias) = 0.8d+0
- v(dinit) = -1.0d+0
- v(eta0) = 1.0d+3 * machep
-c
- 999 return
-c *** last card of vdflt follows ***
- end
- subroutine vscopy(p, y, s)
-c
-c *** set p-vector y to scalar s ***
-c
- integer p
- double precision s, y(p)
-c
- integer i
-c
- do 10 i = 1, p
- 10 y(i) = s
- return
- end
- double precision function v2norm(p, x)
-c
-c *** return the 2-norm of the p-vector x, taking ***
-c *** care to avoid the most likely underflows. ***
-c
- integer p
- double precision x(p)
-c
- integer i, j
- double precision one, r, scale, sqteta, t, xi, zero
-c/+
- double precision dabs, dsqrt
-c/
- external rmdcon
- double precision rmdcon
-c
-c/6
-c data one/1.d+0/, zero/0.d+0/
-c/7
- parameter (one=1.d+0, zero=0.d+0)
- save sqteta
-c/
- data sqteta/0.d+0/
-c
- if (p .gt. 0) go to 10
- v2norm = zero
- go to 999
- 10 do 20 i = 1, p
- if (x(i) .ne. zero) go to 30
- 20 continue
- v2norm = zero
- go to 999
-c
- 30 scale = dabs(x(i))
- if (i .lt. p) go to 40
- v2norm = scale
- go to 999
- 40 t = one
- if (sqteta .eq. zero) sqteta = rmdcon(2)
-c
-c *** sqteta is (slightly larger than) the square root of the
-c *** smallest positive floating point number on the machine.
-c *** the tests involving sqteta are done to prevent underflows.
-c
- j = i + 1
- do 60 i = j, p
- xi = dabs(x(i))
- if (xi .gt. scale) go to 50
- r = xi / scale
- if (r .gt. sqteta) t = t + r*r
- go to 60
- 50 r = scale / xi
- if (r .le. sqteta) r = zero
- t = one + t * r*r
- scale = xi
- 60 continue
-c
- v2norm = scale * dsqrt(t)
- 999 return
-c *** last card of v2norm follows ***
- end
- subroutine humsl(n, d, x, calcf, calcgh, iv, liv, lv, v,
- 1 uiparm, urparm, ufparm)
-c
-c *** minimize general unconstrained objective function using ***
-c *** (analytic) gradient and hessian provided by the caller. ***
-c
- integer liv, lv, n
- integer iv(liv), uiparm(1)
- double precision d(n), x(n), v(lv), urparm(1)
-c dimension v(78 + n*(n+12)), uiparm(*), urparm(*)
- external calcf, calcgh, ufparm
-c
-c------------------------------ discussion ---------------------------
-c
-c this routine is like sumsl, except that the subroutine para-
-c meter calcg of sumsl (which computes the gradient of the objec-
-c tive function) is replaced by the subroutine parameter calcgh,
-c which computes both the gradient and (lower triangle of the)
-c hessian of the objective function. the calling sequence is...
-c call calcgh(n, x, nf, g, h, uiparm, urparm, ufparm)
-c parameters n, x, nf, g, uiparm, urparm, and ufparm are the same
-c as for sumsl, while h is an array of length n*(n+1)/2 in which
-c calcgh must store the lower triangle of the hessian at x. start-
-c ing at h(1), calcgh must store the hessian entries in the order
-c (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ...
-c the value printed (by itsum) in the column labelled stppar
-c is the levenberg-marquardt used in computing the current step.
-c zero means a full newton step. if the special case described in
-c ref. 1 is detected, then stppar is negated. the value printed
-c in the column labelled npreldf is zero if the current hessian
-c is not positive definite.
-c it sometimes proves worthwhile to let d be determined from the
-c diagonal of the hessian matrix by setting iv(dtype) = 1 and
-c v(dinit) = 0. the following iv and v components are relevant...
-c
-c iv(dtol)..... iv(59) gives the starting subscript in v of the dtol
-c array used when d is updated. (iv(dtol) can be
-c initialized by calling humsl with iv(1) = 13.)
-c iv(dtype).... iv(16) tells how the scale vector d should be chosen.
-c iv(dtype) .le. 0 means that d should not be updated, and
-c iv(dtype) .ge. 1 means that d should be updated as
-c described below with v(dfac). default = 0.
-c v(dfac)..... v(41) and the dtol and d0 arrays (see v(dtinit) and
-c v(d0init)) are used in updating the scale vector d when
-c iv(dtype) .gt. 0. (d is initialized according to
-c v(dinit), described in sumsl.) let
-c d1(i) = max(sqrt(abs(h(i,i))), v(dfac)*d(i)),
-c where h(i,i) is the i-th diagonal element of the current
-c hessian. if iv(dtype) = 1, then d(i) is set to d1(i)
-c unless d1(i) .lt. dtol(i), in which case d(i) is set to
-c max(d0(i), dtol(i)).
-c if iv(dtype) .ge. 2, then d is updated during the first
-c iteration as for iv(dtype) = 1 (after any initialization
-c due to v(dinit)) and is left unchanged thereafter.
-c default = 0.6.
-c v(dtinit)... v(39), if positive, is the value to which all components
-c of the dtol array (see v(dfac)) are initialized. if
-c v(dtinit) = 0, then it is assumed that the caller has
-c stored dtol in v starting at v(iv(dtol)).
-c default = 10**-6.
-c v(d0init)... v(40), if positive, is the value to which all components
-c of the d0 vector (see v(dfac)) are initialized. if
-c v(dfac) = 0, then it is assumed that the caller has
-c stored d0 in v starting at v(iv(dtol)+n). default = 1.0.
-c
-c *** reference ***
-c
-c 1. gay, d.m. (1981), computing optimal locally constrained steps,
-c siam j. sci. statist. comput. 2, pp. 186-197.
-c.
-c *** general ***
-c
-c coded by david m. gay (winter 1980). revised sept. 1982.
-c this subroutine was written in connection with research supported
-c in part by the national science foundation under grants
-c mcs-7600324 and mcs-7906671.
-c
-c---------------------------- declarations ---------------------------
-c
- external deflt, humit
-c
-c deflt... provides default input values for iv and v.
-c humit... reverse-communication routine that does humsl algorithm.
-c
- integer g1, h1, iv1, lh, nf
- double precision f
-c
-c *** subscripts for iv ***
-c
- integer g, h, nextv, nfcall, nfgcal, toobig, vneed
-c
-c/6
-c data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, h/56/, toobig/2/,
-c 1 vneed/4/
-c/7
- parameter (nextv=47, nfcall=6, nfgcal=7, g=28, h=56, toobig=2,
- 1 vneed=4)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- lh = n * (n + 1) / 2
- if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
- if (iv(1) .eq. 12 .or. iv(1) .eq. 13)
- 1 iv(vneed) = iv(vneed) + n*(n+3)/2
- iv1 = iv(1)
- if (iv1 .eq. 14) go to 10
- if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
- g1 = 1
- h1 = 1
- if (iv1 .eq. 12) iv(1) = 13
- go to 20
-c
- 10 g1 = iv(g)
- h1 = iv(h)
-c
- 20 call humit(d, f, v(g1), v(h1), iv, lh, liv, lv, n, v, x)
- if (iv(1) - 2) 30, 40, 50
-c
- 30 nf = iv(nfcall)
- call calcf(n, x, nf, f, uiparm, urparm, ufparm)
- if (nf .le. 0) iv(toobig) = 1
- go to 20
-c
- 40 call calcgh(n, x, iv(nfgcal), v(g1), v(h1), uiparm, urparm,
- 1 ufparm)
- go to 20
-c
- 50 if (iv(1) .ne. 14) go to 999
-c
-c *** storage allocation
-c
- iv(g) = iv(nextv)
- iv(h) = iv(g) + n
- iv(nextv) = iv(h) + n*(n+1)/2
- if (iv1 .ne. 13) go to 10
-c
- 999 return
-c *** last card of humsl follows ***
- end
- subroutine humit(d, fx, g, h, iv, lh, liv, lv, n, v, x)
-c
-c *** carry out humsl (unconstrained minimization) iterations, using
-c *** hessian matrix provided by the caller.
-c
-c *** parameter declarations ***
-c
- integer lh, liv, lv, n
- integer iv(liv)
- double precision d(n), fx, g(n), h(lh), v(lv), x(n)
-c
-c-------------------------- parameter usage --------------------------
-c
-c d.... scale vector.
-c fx... function value.
-c g.... gradient vector.
-c h.... lower triangle of the hessian, stored rowwise.
-c iv... integer value array.
-c lh... length of h = p*(p+1)/2.
-c liv.. length of iv (at least 60).
-c lv... length of v (at least 78 + n*(n+21)/2).
-c n.... number of variables (components in x and g).
-c v.... floating-point value array.
-c x.... parameter vector.
-c
-c *** discussion ***
-c
-c parameters iv, n, v, and x are the same as the corresponding
-c ones to humsl (which see), except that v can be shorter (since
-c the part of v that humsl uses for storing g and h is not needed).
-c moreover, compared with humsl, iv(1) may have the two additional
-c output values 1 and 2, which are explained below, as is the use
-c of iv(toobig) and iv(nfgcal). the value iv(g), which is an
-c output value from humsl, is not referenced by humit or the
-c subroutines it calls.
-c
-c iv(1) = 1 means the caller should set fx to f(x), the function value
-c at x, and call humit again, having changed none of the
-c other parameters. an exception occurs if f(x) cannot be
-c computed (e.g. if overflow would occur), which may happen
-c because of an oversized step. in this case the caller
-c should set iv(toobig) = iv(2) to 1, which will cause
-c humit to ignore fx and try a smaller step. the para-
-c meter nf that humsl passes to calcf (for possible use by
-c calcgh) is a copy of iv(nfcall) = iv(6).
-c iv(1) = 2 means the caller should set g to g(x), the gradient of f at
-c x, and h to the lower triangle of h(x), the hessian of f
-c at x, and call humit again, having changed none of the
-c other parameters except perhaps the scale vector d.
-c the parameter nf that humsl passes to calcg is
-c iv(nfgcal) = iv(7). if g(x) and h(x) cannot be evaluated,
-c then the caller may set iv(nfgcal) to 0, in which case
-c humit will return with iv(1) = 65.
-c note -- humit overwrites h with the lower triangle
-c of diag(d)**-1 * h(x) * diag(d)**-1.
-c.
-c *** general ***
-c
-c coded by david m. gay (winter 1980). revised sept. 1982.
-c this subroutine was written in connection with research supported
-c in part by the national science foundation under grants
-c mcs-7600324 and mcs-7906671.
-c
-c (see sumsl and humsl for references.)
-c
-c+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
-c
-c *** local variables ***
-c
- integer dg1, dummy, i, j, k, l, lstgst, nn1o2, step1,
- 1 temp1, w1, x01
- double precision t
-c
-c *** constants ***
-c
- double precision one, onep2, zero
-c
-c *** no intrinsic functions ***
-c
-c *** external functions and subroutines ***
-c
- external assst, deflt, dotprd, dupdu, gqtst, itsum, parck,
- 1 reldst, slvmul, stopx, vaxpy, vcopy, vscopy, v2norm
- logical stopx
- double precision dotprd, reldst, v2norm
-c
-c assst.... assesses candidate step.
-c deflt.... provides default iv and v input values.
-c dotprd... returns inner product of two vectors.
-c dupdu.... updates scale vector d.
-c gqtst.... computes optimally locally constrained step.
-c itsum.... prints iteration summary and info on initial and final x.
-c parck.... checks validity of input iv and v values.
-c reldst... computes v(reldx) = relative step size.
-c slvmul... multiplies symmetric matrix times vector, given the lower
-c triangle of the matrix.
-c stopx.... returns .true. if the break key has been pressed.
-c vaxpy.... computes scalar times one vector plus another.
-c vcopy.... copies one vector to another.
-c vscopy... sets all elements of a vector to a scalar.
-c v2norm... returns the 2-norm of a vector.
-c
-c *** subscripts for iv and v ***
-c
- integer cnvcod, dg, dgnorm, dinit, dstnrm, dtinit, dtol,
- 1 dtype, d0init, f, f0, fdif, gtstep, incfac, irc, kagqt,
- 2 lmat, lmax0, lmaxs, mode, model, mxfcal, mxiter, nextv,
- 3 nfcall, nfgcal, ngcall, niter, preduc, radfac, radinc,
- 4 radius, rad0, reldx, restor, step, stglim, stlstg, stppar,
- 5 toobig, tuner4, tuner5, vneed, w, xirc, x0
-c
-c *** iv subscript values ***
-c
-c/6
-c data cnvcod/55/, dg/37/, dtol/59/, dtype/16/, irc/29/, kagqt/33/,
-c 1 lmat/42/, mode/35/, model/5/, mxfcal/17/, mxiter/18/,
-c 2 nextv/47/, nfcall/6/, nfgcal/7/, ngcall/30/, niter/31/,
-c 3 radinc/8/, restor/9/, step/40/, stglim/11/, stlstg/41/,
-c 4 toobig/2/, vneed/4/, w/34/, xirc/13/, x0/43/
-c/7
- parameter (cnvcod=55, dg=37, dtol=59, dtype=16, irc=29, kagqt=33,
- 1 lmat=42, mode=35, model=5, mxfcal=17, mxiter=18,
- 2 nextv=47, nfcall=6, nfgcal=7, ngcall=30, niter=31,
- 3 radinc=8, restor=9, step=40, stglim=11, stlstg=41,
- 4 toobig=2, vneed=4, w=34, xirc=13, x0=43)
-c/
-c
-c *** v subscript values ***
-c
-c/6
-c data dgnorm/1/, dinit/38/, dstnrm/2/, dtinit/39/, d0init/40/,
-c 1 f/10/, f0/13/, fdif/11/, gtstep/4/, incfac/23/, lmax0/35/,
-c 2 lmaxs/36/, preduc/7/, radfac/16/, radius/8/, rad0/9/,
-c 3 reldx/17/, stppar/5/, tuner4/29/, tuner5/30/
-c/7
- parameter (dgnorm=1, dinit=38, dstnrm=2, dtinit=39, d0init=40,
- 1 f=10, f0=13, fdif=11, gtstep=4, incfac=23, lmax0=35,
- 2 lmaxs=36, preduc=7, radfac=16, radius=8, rad0=9,
- 3 reldx=17, stppar=5, tuner4=29, tuner5=30)
-c/
-c
-c/6
-c data one/1.d+0/, onep2/1.2d+0/, zero/0.d+0/
-c/7
- parameter (one=1.d+0, onep2=1.2d+0, zero=0.d+0)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- i = iv(1)
- if (i .eq. 1) go to 30
- if (i .eq. 2) go to 40
-c
-c *** check validity of iv and v input values ***
-c
- if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
- if (iv(1) .eq. 12 .or. iv(1) .eq. 13)
- 1 iv(vneed) = iv(vneed) + n*(n+21)/2 + 7
- call parck(2, d, iv, liv, lv, n, v)
- i = iv(1) - 2
- if (i .gt. 12) go to 999
- nn1o2 = n * (n + 1) / 2
- if (lh .ge. nn1o2) go to (210,210,210,210,210,210,160,120,160,
- 1 10,10,20), i
- iv(1) = 66
- go to 350
-c
-c *** storage allocation ***
-c
- 10 iv(dtol) = iv(lmat) + nn1o2
- iv(x0) = iv(dtol) + 2*n
- iv(step) = iv(x0) + n
- iv(stlstg) = iv(step) + n
- iv(dg) = iv(stlstg) + n
- iv(w) = iv(dg) + n
- iv(nextv) = iv(w) + 4*n + 7
- if (iv(1) .ne. 13) go to 20
- iv(1) = 14
- go to 999
-c
-c *** initialization ***
-c
- 20 iv(niter) = 0
- iv(nfcall) = 1
- iv(ngcall) = 1
- iv(nfgcal) = 1
- iv(mode) = -1
- iv(model) = 1
- iv(stglim) = 1
- iv(toobig) = 0
- iv(cnvcod) = 0
- iv(radinc) = 0
- v(rad0) = zero
- v(stppar) = zero
- if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
- k = iv(dtol)
- if (v(dtinit) .gt. zero) call vscopy(n, v(k), v(dtinit))
- k = k + n
- if (v(d0init) .gt. zero) call vscopy(n, v(k), v(d0init))
- iv(1) = 1
- go to 999
-c
- 30 v(f) = fx
- if (iv(mode) .ge. 0) go to 210
- iv(1) = 2
- if (iv(toobig) .eq. 0) go to 999
- iv(1) = 63
- go to 350
-c
-c *** make sure gradient could be computed ***
-c
- 40 if (iv(nfgcal) .ne. 0) go to 50
- iv(1) = 65
- go to 350
-c
-c *** update the scale vector d ***
-c
- 50 dg1 = iv(dg)
- if (iv(dtype) .le. 0) go to 70
- k = dg1
- j = 0
- do 60 i = 1, n
- j = j + i
- v(k) = h(j)
- k = k + 1
- 60 continue
- call dupdu(d, v(dg1), iv, liv, lv, n, v)
-c
-c *** compute scaled gradient and its norm ***
-c
- 70 dg1 = iv(dg)
- k = dg1
- do 80 i = 1, n
- v(k) = g(i) / d(i)
- k = k + 1
- 80 continue
- v(dgnorm) = v2norm(n, v(dg1))
-c
-c *** compute scaled hessian ***
-c
- k = 1
- do 100 i = 1, n
- t = one / d(i)
- do 90 j = 1, i
- h(k) = t * h(k) / d(j)
- k = k + 1
- 90 continue
- 100 continue
-c
- if (iv(cnvcod) .ne. 0) go to 340
- if (iv(mode) .eq. 0) go to 300
-c
-c *** allow first step to have scaled 2-norm at most v(lmax0) ***
-c
- v(radius) = v(lmax0)
-c
- iv(mode) = 0
-c
-c
-c----------------------------- main loop -----------------------------
-c
-c
-c *** print iteration summary, check iteration limit ***
-c
- 110 call itsum(d, g, iv, liv, lv, n, v, x)
- 120 k = iv(niter)
- if (k .lt. iv(mxiter)) go to 130
- iv(1) = 10
- go to 350
-c
- 130 iv(niter) = k + 1
-c
-c *** initialize for start of next iteration ***
-c
- dg1 = iv(dg)
- x01 = iv(x0)
- v(f0) = v(f)
- iv(irc) = 4
- iv(kagqt) = -1
-c
-c *** copy x to x0 ***
-c
- call vcopy(n, v(x01), x)
-c
-c *** update radius ***
-c
- if (k .eq. 0) go to 150
- step1 = iv(step)
- k = step1
- do 140 i = 1, n
- v(k) = d(i) * v(k)
- k = k + 1
- 140 continue
- v(radius) = v(radfac) * v2norm(n, v(step1))
-c
-c *** check stopx and function evaluation limit ***
-c
-C AL 4/30/95
- dummy=iv(nfcall)
- 150 if (.not. stopx(dummy)) go to 170
- iv(1) = 11
- go to 180
-c
-c *** come here when restarting after func. eval. limit or stopx.
-c
- 160 if (v(f) .ge. v(f0)) go to 170
- v(radfac) = one
- k = iv(niter)
- go to 130
-c
- 170 if (iv(nfcall) .lt. iv(mxfcal)) go to 190
- iv(1) = 9
- 180 if (v(f) .ge. v(f0)) go to 350
-c
-c *** in case of stopx or function evaluation limit with
-c *** improved v(f), evaluate the gradient at x.
-c
- iv(cnvcod) = iv(1)
- go to 290
-c
-c. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
-c
- 190 step1 = iv(step)
- dg1 = iv(dg)
- l = iv(lmat)
- w1 = iv(w)
- call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1))
- if (iv(irc) .eq. 6) go to 210
-c
-c *** check whether evaluating f(x0 + step) looks worthwhile ***
-c
- if (v(dstnrm) .le. zero) go to 210
- if (iv(irc) .ne. 5) go to 200
- if (v(radfac) .le. one) go to 200
- if (v(preduc) .le. onep2 * v(fdif)) go to 210
-c
-c *** compute f(x0 + step) ***
-c
- 200 x01 = iv(x0)
- step1 = iv(step)
- call vaxpy(n, x, one, v(step1), v(x01))
- iv(nfcall) = iv(nfcall) + 1
- iv(1) = 1
- iv(toobig) = 0
- go to 999
-c
-c. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
-c
- 210 x01 = iv(x0)
- v(reldx) = reldst(n, d, x, v(x01))
- call assst(iv, liv, lv, v)
- step1 = iv(step)
- lstgst = iv(stlstg)
- if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
- if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
- if (iv(restor) .ne. 3) go to 220
- call vcopy(n, v(step1), v(lstgst))
- call vaxpy(n, x, one, v(step1), v(x01))
- v(reldx) = reldst(n, d, x, v(x01))
-c
- 220 k = iv(irc)
- go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k
-c
-c *** recompute step with new radius ***
-c
- 230 v(radius) = v(radfac) * v(dstnrm)
- go to 150
-c
-c *** compute step of length v(lmaxs) for singular convergence test.
-c
- 240 v(radius) = v(lmaxs)
- go to 190
-c
-c *** convergence or false convergence ***
-c
- 250 iv(cnvcod) = k - 4
- if (v(f) .ge. v(f0)) go to 340
- if (iv(xirc) .eq. 14) go to 340
- iv(xirc) = 14
-c
-c. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
-c
- 260 if (iv(irc) .ne. 3) go to 290
- temp1 = lstgst
-c
-c *** prepare for gradient tests ***
-c *** set temp1 = hessian * step + g(x0)
-c *** = diag(d) * (h * step + g(x0))
-c
-c use x0 vector as temporary.
- k = x01
- do 270 i = 1, n
- v(k) = d(i) * v(step1)
- k = k + 1
- step1 = step1 + 1
- 270 continue
- call slvmul(n, v(temp1), h, v(x01))
- do 280 i = 1, n
- v(temp1) = d(i) * v(temp1) + g(i)
- temp1 = temp1 + 1
- 280 continue
-c
-c *** compute gradient and hessian ***
-c
- 290 iv(ngcall) = iv(ngcall) + 1
- iv(1) = 2
- go to 999
-c
- 300 iv(1) = 2
- if (iv(irc) .ne. 3) go to 110
-c
-c *** set v(radfac) by gradient tests ***
-c
- temp1 = iv(stlstg)
- step1 = iv(step)
-c
-c *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
-c
- k = temp1
- do 310 i = 1, n
- v(k) = (v(k) - g(i)) / d(i)
- k = k + 1
- 310 continue
-c
-c *** do gradient tests ***
-c
- if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320
- if (dotprd(n, g, v(step1))
- 1 .ge. v(gtstep) * v(tuner5)) go to 110
- 320 v(radfac) = v(incfac)
- go to 110
-c
-c. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
-c
-c *** bad parameters to assess ***
-c
- 330 iv(1) = 64
- go to 350
-c
-c *** print summary of final iteration and other requested items ***
-c
- 340 iv(1) = iv(cnvcod)
- iv(cnvcod) = 0
- 350 call itsum(d, g, iv, liv, lv, n, v, x)
-c
- 999 return
-c
-c *** last card of humit follows ***
- end
- subroutine dupdu(d, hdiag, iv, liv, lv, n, v)
-c
-c *** update scale vector d for humsl ***
-c
-c *** parameter declarations ***
-c
- integer liv, lv, n
- integer iv(liv)
- double precision d(n), hdiag(n), v(lv)
-c
-c *** local variables ***
-c
- integer dtoli, d0i, i
- double precision t, vdfac
-c
-c *** intrinsic functions ***
-c/+
- double precision dabs, dmax1, dsqrt
-c/
-c *** subscripts for iv and v ***
-c
- integer dfac, dtol, dtype, niter
-c/6
-c data dfac/41/, dtol/59/, dtype/16/, niter/31/
-c/7
- parameter (dfac=41, dtol=59, dtype=16, niter=31)
-c/
-c
-c------------------------------- body --------------------------------
-c
- i = iv(dtype)
- if (i .eq. 1) go to 10
- if (iv(niter) .gt. 0) go to 999
-c
- 10 dtoli = iv(dtol)
- d0i = dtoli + n
- vdfac = v(dfac)
- do 20 i = 1, n
- t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i))
- if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i))
- d(i) = t
- dtoli = dtoli + 1
- d0i = d0i + 1
- 20 continue
-c
- 999 return
-c *** last card of dupdu follows ***
- end
- subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w)
-c
-c *** compute goldfeld-quandt-trotter step by more-hebden technique ***
-c *** (nl2sol version 2.2), modified a la more and sorensen ***
-c
-c *** parameter declarations ***
-c
- integer ka, p
-cal double precision d(p), dig(p), dihdi(1), l(1), v(21), step(p),
-cal 1 w(1)
- double precision d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),
- 1 v(21), step(p),w(4*p+7)
-c dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7)
-c
-c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-c
-c *** purpose ***
-c
-c given the (compactly stored) lower triangle of a scaled
-c hessian (approximation) and a nonzero scaled gradient vector,
-c this subroutine computes a goldfeld-quandt-trotter step of
-c approximate length v(radius) by the more-hebden technique. in
-c other words, step is computed to (approximately) minimize
-c psi(step) = (g**t)*step + 0.5*(step**t)*h*step such that the
-c 2-norm of d*step is at most (approximately) v(radius), where
-c g is the gradient, h is the hessian, and d is a diagonal
-c scale matrix whose diagonal is stored in the parameter d.
-c (gqtst assumes dig = d**-1 * g and dihdi = d**-1 * h * d**-1.)
-c
-c *** parameter description ***
-c
-c d (in) = the scale vector, i.e. the diagonal of the scale
-c matrix d mentioned above under purpose.
-c dig (in) = the scaled gradient vector, d**-1 * g. if g = 0, then
-c step = 0 and v(stppar) = 0 are returned.
-c dihdi (in) = lower triangle of the scaled hessian (approximation),
-c i.e., d**-1 * h * d**-1, stored compactly by rows., i.e.,
-c in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc.
-c ka (i/o) = the number of hebden iterations (so far) taken to deter-
-c mine step. ka .lt. 0 on input means this is the first
-c attempt to determine step (for the present dig and dihdi)
-c -- ka is initialized to 0 in this case. output with
-c ka = 0 (or v(stppar) = 0) means step = -(h**-1)*g.
-c l (i/o) = workspace of length p*(p+1)/2 for cholesky factors.
-c p (in) = number of parameters -- the hessian is a p x p matrix.
-c step (i/o) = the step computed.
-c v (i/o) contains various constants and variables described below.
-c w (i/o) = workspace of length 4*p + 6.
-c
-c *** entries in v ***
-c
-c v(dgnorm) (i/o) = 2-norm of (d**-1)*g.
-c v(dstnrm) (output) = 2-norm of d*step.
-c v(dst0) (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or
-c overestimate of smallest eigenvalue of (d**-1)*h*(d**-1).
-c v(epslon) (in) = max. rel. error allowed for psi(step). for the
-c step returned, psi(step) will exceed its optimal value
-c by less than -v(epslon)*psi(step). suggested value = 0.1.
-c v(gtstep) (out) = inner product between g and step.
-c v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step) (for pos. def.
-c h only -- v(nreduc) is set to zero otherwise).
-c v(phmnfc) (in) = tol. (together with v(phmxfc)) for accepting step
-c (more*s sigma). the error v(dstnrm) - v(radius) must lie
-c between v(phmnfc)*v(radius) and v(phmxfc)*v(radius).
-c v(phmxfc) (in) (see v(phmnfc).)
-c suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5.
-c v(preduc) (out) = psi(step) = predicted obj. func. reduction for step.
-c v(radius) (in) = radius of current (scaled) trust region.
-c v(rad0) (i/o) = value of v(radius) from previous call.
-c v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha
-c described below under algorithm notes. if h + alpha*d**2
-c (see algorithm notes) is (nearly) singular, however,
-c then v(stppar) = -alpha.
-c
-c *** usage notes ***
-c
-c if it is desired to recompute step using a different value of
-c v(radius), then this routine may be restarted by calling it
-c with all parameters unchanged except v(radius). (this explains
-c why step and w are listed as i/o). on an initial call (one with
-c ka .lt. 0), step and w need not be initialized and only compo-
-c nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and
-c v(rad0) of v must be initialized.
-c
-c *** algorithm notes ***
-c
-c the desired g-q-t step (ref. 2, 3, 4, 6) satisfies
-c (h + alpha*d**2)*step = -g for some nonnegative alpha such that
-c h + alpha*d**2 is positive semidefinite. alpha and step are
-c computed by a scheme analogous to the one described in ref. 5.
-c estimates of the smallest and largest eigenvalues of the hessian
-c are obtained from the gerschgorin circle theorem enhanced by a
-c simple form of the scaling described in ref. 7. cases in which
-c h + alpha*d**2 is nearly (or exactly) singular are handled by
-c the technique discussed in ref. 2. in these cases, a step of
-c (exact) length v(radius) is returned for which psi(step) exceeds
-c its optimal value by less than -v(epslon)*psi(step). the test
-c suggested in ref. 6 for detecting the special case is performed
-c once two matrix factorizations have been done -- doing so sooner
-c seems to degrade the performance of optimization routines that
-c call this routine.
-c
-c *** functions and subroutines called ***
-c
-c dotprd - returns inner product of two vectors.
-c litvmu - applies inverse-transpose of compact lower triang. matrix.
-c livmul - applies inverse of compact lower triang. matrix.
-c lsqrt - finds cholesky factor (of compactly stored lower triang.).
-c lsvmin - returns approx. to min. sing. value of lower triang. matrix.
-c rmdcon - returns machine-dependent constants.
-c v2norm - returns 2-norm of a vector.
-c
-c *** references ***
-c
-c 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive
-c nonlinear least-squares algorithm, acm trans. math.
-c software, vol. 7, no. 3.
-c 2. gay, d.m. (1981), computing optimal locally constrained steps,
-c siam j. sci. statist. computing, vol. 2, no. 2, pp.
-c 186-197.
-c 3. goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966),
-c maximization by quadratic hill-climbing, econometrica 34,
-c pp. 541-551.
-c 4. hebden, m.d. (1973), an algorithm for minimization using exact
-c second derivatives, report t.p. 515, theoretical physics
-c div., a.e.r.e. harwell, oxon., england.
-c 5. more, j.j. (1978), the levenberg-marquardt algorithm, implemen-
-c tation and theory, pp.105-116 of springer lecture notes
-c in mathematics no. 630, edited by g.a. watson, springer-
-c verlag, berlin and new york.
-c 6. more, j.j., and sorensen, d.c. (1981), computing a trust region
-c step, technical report anl-81-83, argonne national lab.
-c 7. varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15,
-c pp. 719-729.
-c
-c *** general ***
-c
-c coded by david m. gay.
-c this subroutine was written in connection with research
-c supported by the national science foundation under grants
-c mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
-c mcs-7906671.
-c
-c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-c
-c *** local variables ***
-c
- logical restrt
- integer dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,
- 1 j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x
- double precision alphak, aki, akk, delta, dst, eps, gtsta, lk,
- 1 oldphi, phi, phimax, phimin, psifac, rad, radsq,
- 2 root, si, sk, sw, t, twopsi, t1, t2, uk, wi
-c
-c *** constants ***
- double precision big, dgxfac, epsfac, four, half, kappa, negone,
- 1 one, p001, six, three, two, zero
-c
-c *** intrinsic functions ***
-c/+
- double precision dabs, dmax1, dmin1, dsqrt
-c/
-c *** external functions and subroutines ***
-c
- external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm
- double precision dotprd, lsvmin, rmdcon, v2norm
-c
-c *** subscripts for v ***
-c
- integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc,
- 1 phmnfc, phmxfc, preduc, radius, rad0
-c/6
-c data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/,
-c 1 nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/,
-c 2 rad0/9/, stppar/5/
-c/7
- parameter (dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,
- 1 nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,
- 2 rad0=9, stppar=5)
-c/
-c
-c/6
-c data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/,
-c 1 kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/,
-c 2 six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/
-c/7
- parameter (epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,
- 1 kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,
- 2 six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0)
- save dgxfac
-c/
- data big/0.d+0/, dgxfac/0.d+0/
-c
-c *** body ***
-c
-c *** store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx).
- dggdmx = p + 1
-c *** store gerschgorin over- and underestimates of the largest
-c *** and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax)
-c *** and w(emin) respectively.
- emax = dggdmx + 1
- emin = emax + 1
-c *** for use in recomputing step, the final values of lk, uk, dst,
-c *** and the inverse derivative of more*s phi at 0 (for pos. def.
-c *** h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin)
-c *** respectively.
- lk0 = emin + 1
- phipin = lk0 + 1
- uk0 = phipin + 1
- dstsav = uk0 + 1
-c *** store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p).
- diag0 = dstsav
- diag = diag0 + 1
-c *** store -d*step in w(q),...,w(q0+p).
- q0 = diag0 + p
- q = q0 + 1
-c *** allocate storage for scratch vector x ***
- x = q + p
- rad = v(radius)
- radsq = rad**2
-c *** phitol = max. error allowed in dst = v(dstnrm) = 2-norm of
-c *** d*step.
- phimax = v(phmxfc) * rad
- phimin = v(phmnfc) * rad
- psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) *
- 1 (kappa + one) + kappa + two) * rad**2)
-c *** oldphi is used to detect limits of numerical accuracy. if
-c *** we recompute step and it does not change, then we accept it.
- oldphi = zero
- eps = v(epslon)
- irc = 0
- restrt = .false.
- kalim = ka + 50
-c
-c *** start or restart, depending on ka ***
-c
- if (ka .ge. 0) go to 290
-c
-c *** fresh start ***
-c
- k = 0
- uk = negone
- ka = 0
- kalim = 50
- v(dgnorm) = v2norm(p, dig)
- v(nreduc) = zero
- v(dst0) = zero
- kamin = 3
- if (v(dgnorm) .eq. zero) kamin = 0
-c
-c *** store diag(dihdi) in w(diag0+1),...,w(diag0+p) ***
-c
- j = 0
- do 10 i = 1, p
- j = j + i
- k1 = diag0 + i
- w(k1) = dihdi(j)
- 10 continue
-c
-c *** determine w(dggdmx), the largest element of dihdi ***
-c
- t1 = zero
- j = p * (p + 1) / 2
- do 20 i = 1, j
- t = dabs(dihdi(i))
- if (t1 .lt. t) t1 = t
- 20 continue
- w(dggdmx) = t1
-c
-c *** try alpha = 0 ***
-c
- 30 call lsqrt(1, p, l, dihdi, irc)
- if (irc .eq. 0) go to 50
-c *** indef. h -- underestimate smallest eigenvalue, use this
-c *** estimate to initialize lower bound lk on alpha.
- j = irc*(irc+1)/2
- t = l(j)
- l(j) = one
- do 40 i = 1, irc
- 40 w(i) = zero
- w(irc) = one
- call litvmu(irc, w, l, w)
- t1 = v2norm(irc, w)
- lk = -t / t1 / t1
- v(dst0) = -lk
- if (restrt) go to 210
- go to 70
-c
-c *** positive definite h -- compute unmodified newton step. ***
- 50 lk = zero
- t = lsvmin(p, l, w(q), w(q))
- if (t .ge. one) go to 60
- if (big .le. zero) big = rmdcon(6)
- if (v(dgnorm) .ge. t*t*big) go to 70
- 60 call livmul(p, w(q), l, dig)
- gtsta = dotprd(p, w(q), w(q))
- v(nreduc) = half * gtsta
- call litvmu(p, w(q), l, w(q))
- dst = v2norm(p, w(q))
- v(dst0) = dst
- phi = dst - rad
- if (phi .le. phimax) go to 260
- if (restrt) go to 210
-c
-c *** prepare to compute gerschgorin estimates of largest (and
-c *** smallest) eigenvalues. ***
-c
- 70 k = 0
- do 100 i = 1, p
- wi = zero
- if (i .eq. 1) go to 90
- im1 = i - 1
- do 80 j = 1, im1
- k = k + 1
- t = dabs(dihdi(k))
- wi = wi + t
- w(j) = w(j) + t
- 80 continue
- 90 w(i) = wi
- k = k + 1
- 100 continue
-c
-c *** (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1) ***
-c
- k = 1
- t1 = w(diag) - w(1)
- if (p .le. 1) go to 120
- do 110 i = 2, p
- j = diag0 + i
- t = w(j) - w(i)
- if (t .ge. t1) go to 110
- t1 = t
- k = i
- 110 continue
-c
- 120 sk = w(k)
- j = diag0 + k
- akk = w(j)
- k1 = k*(k-1)/2 + 1
- inc = 1
- t = zero
- do 150 i = 1, p
- if (i .eq. k) go to 130
- aki = dabs(dihdi(k1))
- si = w(i)
- j = diag0 + i
- t1 = half * (akk - w(j) + si - aki)
- t1 = t1 + dsqrt(t1*t1 + sk*aki)
- if (t .lt. t1) t = t1
- if (i .lt. k) go to 140
- 130 inc = i
- 140 k1 = k1 + inc
- 150 continue
-c
- w(emin) = akk - t
- uk = v(dgnorm)/rad - w(emin)
- if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
- if (uk .le. zero) uk = p001
-c
-c *** compute gerschgorin (over-)estimate of largest eigenvalue ***
-c
- k = 1
- t1 = w(diag) + w(1)
- if (p .le. 1) go to 170
- do 160 i = 2, p
- j = diag0 + i
- t = w(j) + w(i)
- if (t .le. t1) go to 160
- t1 = t
- k = i
- 160 continue
-c
- 170 sk = w(k)
- j = diag0 + k
- akk = w(j)
- k1 = k*(k-1)/2 + 1
- inc = 1
- t = zero
- do 200 i = 1, p
- if (i .eq. k) go to 180
- aki = dabs(dihdi(k1))
- si = w(i)
- j = diag0 + i
- t1 = half * (w(j) + si - aki - akk)
- t1 = t1 + dsqrt(t1*t1 + sk*aki)
- if (t .lt. t1) t = t1
- if (i .lt. k) go to 190
- 180 inc = i
- 190 k1 = k1 + inc
- 200 continue
-c
- w(emax) = akk + t
- lk = dmax1(lk, v(dgnorm)/rad - w(emax))
-c
-c *** alphak = current value of alpha (see alg. notes above). we
-c *** use more*s scheme for initializing it.
- alphak = dabs(v(stppar)) * v(rad0)/rad
-c
- if (irc .ne. 0) go to 210
-c
-c *** compute l0 for positive definite h ***
-c
- call livmul(p, w, l, w(q))
- t = v2norm(p, w)
- w(phipin) = dst / t / t
- lk = dmax1(lk, phi*w(phipin))
-c
-c *** safeguard alphak and add alphak*i to (d**-1)*h*(d**-1) ***
-c
- 210 ka = ka + 1
- if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk)
- 1 alphak = uk * dmax1(p001, dsqrt(lk/uk))
- if (alphak .le. zero) alphak = half * uk
- if (alphak .le. zero) alphak = uk
- k = 0
- do 220 i = 1, p
- k = k + i
- j = diag0 + i
- dihdi(k) = w(j) + alphak
- 220 continue
-c
-c *** try computing cholesky decomposition ***
-c
- call lsqrt(1, p, l, dihdi, irc)
- if (irc .eq. 0) go to 240
-c
-c *** (d**-1)*h*(d**-1) + alphak*i is indefinite -- overestimate
-c *** smallest eigenvalue for use in updating lk ***
-c
- j = (irc*(irc+1))/2
- t = l(j)
- l(j) = one
- do 230 i = 1, irc
- 230 w(i) = zero
- w(irc) = one
- call litvmu(irc, w, l, w)
- t1 = v2norm(irc, w)
- lk = alphak - t/t1/t1
- v(dst0) = -lk
- go to 210
-c
-c *** alphak makes (d**-1)*h*(d**-1) positive definite.
-c *** compute q = -d*step, check for convergence. ***
-c
- 240 call livmul(p, w(q), l, dig)
- gtsta = dotprd(p, w(q), w(q))
- call litvmu(p, w(q), l, w(q))
- dst = v2norm(p, w(q))
- phi = dst - rad
- if (phi .le. phimax .and. phi .ge. phimin) go to 270
- if (phi .eq. oldphi) go to 270
- oldphi = phi
- if (phi .lt. zero) go to 330
-c
-c *** unacceptable alphak -- update lk, uk, alphak ***
-c
- 250 if (ka .ge. kalim) go to 270
-c *** the following dmin1 is necessary because of restarts ***
- if (phi .lt. zero) uk = dmin1(uk, alphak)
-c *** kamin = 0 only iff the gradient vanishes ***
- if (kamin .eq. 0) go to 210
- call livmul(p, w, l, w(q))
- t1 = v2norm(p, w)
- alphak = alphak + (phi/t1) * (dst/t1) * (dst/rad)
- lk = dmax1(lk, alphak)
- go to 210
-c
-c *** acceptable step on first try ***
-c
- 260 alphak = zero
-c
-c *** successful step in general. compute step = -(d**-1)*q ***
-c
- 270 do 280 i = 1, p
- j = q0 + i
- step(i) = -w(j)/d(i)
- 280 continue
- v(gtstep) = -gtsta
- v(preduc) = half * (dabs(alphak)*dst*dst + gtsta)
- go to 410
-c
-c
-c *** restart with new radius ***
-c
- 290 if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310
-c
-c *** prepare to return newton step ***
-c
- restrt = .true.
- ka = ka + 1
- k = 0
- do 300 i = 1, p
- k = k + i
- j = diag0 + i
- dihdi(k) = w(j)
- 300 continue
- uk = negone
- go to 30
-c
- 310 kamin = ka + 3
- if (v(dgnorm) .eq. zero) kamin = 0
- if (ka .eq. 0) go to 50
-c
- dst = w(dstsav)
- alphak = dabs(v(stppar))
- phi = dst - rad
- t = v(dgnorm)/rad
- uk = t - w(emin)
- if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
- if (uk .le. zero) uk = p001
- if (rad .gt. v(rad0)) go to 320
-c
-c *** smaller radius ***
- lk = zero
- if (alphak .gt. zero) lk = w(lk0)
- lk = dmax1(lk, t - w(emax))
- if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
- go to 250
-c
-c *** bigger radius ***
- 320 if (alphak .gt. zero) uk = dmin1(uk, w(uk0))
- lk = dmax1(zero, -v(dst0), t - w(emax))
- if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
- go to 250
-c
-c *** decide whether to check for special case... in practice (from
-c *** the standpoint of the calling optimization code) it seems best
-c *** not to check until a few iterations have failed -- hence the
-c *** test on kamin below.
-c
- 330 delta = alphak + dmin1(zero, v(dst0))
- twopsi = alphak*dst*dst + gtsta
- if (ka .ge. kamin) go to 340
-c *** if the test in ref. 2 is satisfied, fall through to handle
-c *** the special case (as soon as the more-sorensen test detects
-c *** it).
- if (delta .ge. psifac*twopsi) go to 370
-c
-c *** check for the special case of h + alpha*d**2 (nearly)
-c *** singular. use one step of inverse power method with start
-c *** from lsvmin to obtain approximate eigenvector corresponding
-c *** to smallest eigenvalue of (d**-1)*h*(d**-1). lsvmin returns
-c *** x and w with l*w = x.
-c
- 340 t = lsvmin(p, l, w(x), w)
-c
-c *** normalize w ***
- do 350 i = 1, p
- 350 w(i) = t*w(i)
-c *** complete current inv. power iter. -- replace w by (l**-t)*w.
- call litvmu(p, w, l, w)
- t2 = one/v2norm(p, w)
- do 360 i = 1, p
- 360 w(i) = t2*w(i)
- t = t2 * t
-c
-c *** now w is the desired approximate (unit) eigenvector and
-c *** t*x = ((d**-1)*h*(d**-1) + alphak*i)*w.
-c
- sw = dotprd(p, w(q), w)
- t1 = (rad + dst) * (rad - dst)
- root = dsqrt(sw*sw + t1)
- if (sw .lt. zero) root = -root
- si = t1 / (sw + root)
-c
-c *** the actual test for the special case...
-c
- if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380
-c
-c *** update upper bound on smallest eigenvalue (when not positive)
-c *** (as recommended by more and sorensen) and continue...
-c
- if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak)
- lk = dmax1(lk, -v(dst0))
-c
-c *** check whether we can hope to detect the special case in
-c *** the available arithmetic. accept step as it is if not.
-c
-c *** if not yet available, obtain machine dependent value dgxfac.
- 370 if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3)
-c
- if (delta .gt. dgxfac*w(dggdmx)) go to 250
- go to 270
-c
-c *** special case detected... negate alphak to indicate special case
-c
- 380 alphak = -alphak
- v(preduc) = half * twopsi
-c
-c *** accept current step if adding si*w would lead to a
-c *** further relative reduction in psi of less than v(epslon)/3.
-c
- t1 = zero
- t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w)))
- if (t .lt. eps*twopsi/six) go to 390
- v(preduc) = v(preduc) + t
- dst = rad
- t1 = -si
- 390 do 400 i = 1, p
- j = q0 + i
- w(j) = t1*w(i) - w(j)
- step(i) = w(j) / d(i)
- 400 continue
- v(gtstep) = dotprd(p, dig, w(q))
-c
-c *** save values for use in a possible restart ***
-c
- 410 v(dstnrm) = dst
- v(stppar) = alphak
- w(lk0) = lk
- w(uk0) = uk
- v(rad0) = rad
- w(dstsav) = dst
-c
-c *** restore diagonal of dihdi ***
-c
- j = 0
- do 420 i = 1, p
- j = j + i
- k = diag0 + i
- dihdi(j) = w(k)
- 420 continue
-c
- 999 return
-c
-c *** last card of gqtst follows ***
- end
- subroutine lsqrt(n1, n, l, a, irc)
-c
-c *** compute rows n1 through n of the cholesky factor l of
-c *** a = l*(l**t), where l and the lower triangle of a are both
-c *** stored compactly by rows (and may occupy the same storage).
-c *** irc = 0 means all went well. irc = j means the leading
-c *** principal j x j submatrix of a is not positive definite --
-c *** and l(j*(j+1)/2) contains the (nonpos.) reduced j-th diagonal.
-c
-c *** parameters ***
-c
- integer n1, n, irc
-cal double precision l(1), a(1)
- double precision l(n*(n+1)/2), a(n*(n+1)/2)
-c dimension l(n*(n+1)/2), a(n*(n+1)/2)
-c
-c *** local variables ***
-c
- integer i, ij, ik, im1, i0, j, jk, jm1, j0, k
- double precision t, td, zero
-c
-c *** intrinsic functions ***
-c/+
- double precision dsqrt
-c/
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c
-c *** body ***
-c
- i0 = n1 * (n1 - 1) / 2
- do 50 i = n1, n
- td = zero
- if (i .eq. 1) go to 40
- j0 = 0
- im1 = i - 1
- do 30 j = 1, im1
- t = zero
- if (j .eq. 1) go to 20
- jm1 = j - 1
- do 10 k = 1, jm1
- ik = i0 + k
- jk = j0 + k
- t = t + l(ik)*l(jk)
- 10 continue
- 20 ij = i0 + j
- j0 = j0 + j
- t = (a(ij) - t) / l(j0)
- l(ij) = t
- td = td + t*t
- 30 continue
- 40 i0 = i0 + i
- t = a(i0) - td
- if (t .le. zero) go to 60
- l(i0) = dsqrt(t)
- 50 continue
-c
- irc = 0
- go to 999
-c
- 60 l(i0) = t
- irc = i
-c
- 999 return
-c
-c *** last card of lsqrt ***
- end
- double precision function lsvmin(p, l, x, y)
-c
-c *** estimate smallest sing. value of packed lower triang. matrix l
-c
-c *** parameter declarations ***
-c
- integer p
-cal double precision l(1), x(p), y(p)
- double precision l(p*(p+1)/2), x(p), y(p)
-c dimension l(p*(p+1)/2)
-c
-c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-c
-c *** purpose ***
-c
-c this function returns a good over-estimate of the smallest
-c singular value of the packed lower triangular matrix l.
-c
-c *** parameter description ***
-c
-c p (in) = the order of l. l is a p x p lower triangular matrix.
-c l (in) = array holding the elements of l in row order, i.e.
-c l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc.
-c x (out) if lsvmin returns a positive value, then x is a normalized
-c approximate left singular vector corresponding to the
-c smallest singular value. this approximation may be very
-c crude. if lsvmin returns zero, then some components of x
-c are zero and the rest retain their input values.
-c y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an
-c unnormalized approximate right singular vector correspond-
-c ing to the smallest singular value. this approximation
-c may be crude. if lsvmin returns zero, then y retains its
-c input value. the caller may pass the same vector for x
-c and y (nonstandard fortran usage), in which case y over-
-c writes x (for nonzero lsvmin returns).
-c
-c *** algorithm notes ***
-c
-c the algorithm is based on (1), with the additional provision that
-c lsvmin = 0 is returned if the smallest diagonal element of l
-c (in magnitude) is not more than the unit roundoff times the
-c largest. the algorithm uses a random number generator proposed
-c in (4), which passes the spectral test with flying colors -- see
-c (2) and (3).
-c
-c *** subroutines and functions called ***
-c
-c v2norm - function, returns the 2-norm of a vector.
-c
-c *** references ***
-c
-c (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977),
-c an estimate for the condition number of a matrix, report
-c tm-310, applied math. div., argonne national laboratory.
-c
-c (2) hoaglin, d.c. (1976), theoretical properties of congruential
-c random-number generators -- an empirical view,
-c memorandum ns-340, dept. of statistics, harvard univ.
-c
-c (3) knuth, d.e. (1969), the art of computer programming, vol. 2
-c (seminumerical algorithms), addison-wesley, reading, mass.
-c
-c (4) smith, c.s. (1971), multiplicative pseudo-random number
-c generators with prime modulus, j. assoc. comput. mach. 18,
-c pp. 586-593.
-c
-c *** history ***
-c
-c designed and coded by david m. gay (winter 1977/summer 1978).
-c
-c *** general ***
-c
-c this subroutine was written in connection with research
-c supported by the national science foundation under grants
-c mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989.
-c
-c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-c
-c *** local variables ***
-c
- integer i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1
- double precision b, sminus, splus, t, xminus, xplus
-c
-c *** constants ***
-c
- double precision half, one, r9973, zero
-c
-c *** intrinsic functions ***
-c/+
- integer mod
- real float
- double precision dabs
-c/
-c *** external functions and subroutines ***
-c
- external dotprd, v2norm, vaxpy
- double precision dotprd, v2norm
-c
-c/6
-c data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/
-c/7
- parameter (half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0)
-c/
-c
-c *** body ***
-c
- ix = 2
- pm1 = p - 1
-c
-c *** first check whether to return lsvmin = 0 and initialize x ***
-c
- ii = 0
- j0 = p*pm1/2
- jj = j0 + p
- if (l(jj) .eq. zero) go to 110
- ix = mod(3432*ix, 9973)
- b = half*(one + float(ix)/r9973)
- xplus = b / l(jj)
- x(p) = xplus
- if (p .le. 1) go to 60
- do 10 i = 1, pm1
- ii = ii + i
- if (l(ii) .eq. zero) go to 110
- ji = j0 + i
- x(i) = xplus * l(ji)
- 10 continue
-c
-c *** solve (l**t)*x = b, where the components of b have randomly
-c *** chosen magnitudes in (.5,1) with signs chosen to make x large.
-c
-c do j = p-1 to 1 by -1...
- do 50 jjj = 1, pm1
- j = p - jjj
-c *** determine x(j) in this iteration. note for i = 1,2,...,j
-c *** that x(i) holds the current partial sum for row i.
- ix = mod(3432*ix, 9973)
- b = half*(one + float(ix)/r9973)
- xplus = (b - x(j))
- xminus = (-b - x(j))
- splus = dabs(xplus)
- sminus = dabs(xminus)
- jm1 = j - 1
- j0 = j*jm1/2
- jj = j0 + j
- xplus = xplus/l(jj)
- xminus = xminus/l(jj)
- if (jm1 .eq. 0) go to 30
- do 20 i = 1, jm1
- ji = j0 + i
- splus = splus + dabs(x(i) + l(ji)*xplus)
- sminus = sminus + dabs(x(i) + l(ji)*xminus)
- 20 continue
- 30 if (sminus .gt. splus) xplus = xminus
- x(j) = xplus
-c *** update partial sums ***
- if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x)
- 50 continue
-c
-c *** normalize x ***
-c
- 60 t = one/v2norm(p, x)
- do 70 i = 1, p
- 70 x(i) = t*x(i)
-c
-c *** solve l*y = x and return lsvmin = 1/twonorm(y) ***
-c
- do 100 j = 1, p
- jm1 = j - 1
- j0 = j*jm1/2
- jj = j0 + j
- t = zero
- if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y)
- y(j) = (x(j) - t) / l(jj)
- 100 continue
-c
- lsvmin = one/v2norm(p, y)
- go to 999
-c
- 110 lsvmin = zero
- 999 return
-c *** last card of lsvmin follows ***
- end
- subroutine slvmul(p, y, s, x)
-c
-c *** set y = s * x, s = p x p symmetric matrix. ***
-c *** lower triangle of s stored rowwise. ***
-c
-c *** parameter declarations ***
-c
- integer p
-cal double precision s(1), x(p), y(p)
- double precision s(p*(p+1)/2), x(p), y(p)
-c dimension s(p*(p+1)/2)
-c
-c *** local variables ***
-c
- integer i, im1, j, k
- double precision xi
-c
-c *** no intrinsic functions ***
-c
-c *** external function ***
-c
- external dotprd
- double precision dotprd
-c
-c-----------------------------------------------------------------------
-c
- j = 1
- do 10 i = 1, p
- y(i) = dotprd(i, s(j), x)
- j = j + i
- 10 continue
-c
- if (p .le. 1) go to 999
- j = 1
- do 40 i = 2, p
- xi = x(i)
- im1 = i - 1
- j = j + 1
- do 30 k = 1, im1
- y(k) = y(k) + s(j)*xi
- j = j + 1
- 30 continue
- 40 continue
-c
- 999 return
-c *** last card of slvmul follows ***
- end
projects / unres.git / blobdiff