+++ /dev/null
- subroutine sumsl(n, d, x, calcf, calcg, iv, liv, lv, v,
- 1 uiparm, urparm, ufparm)
-c
-c *** minimize general unconstrained objective function using ***
-c *** analytic gradient and hessian approx. from secant update ***
-c
- integer n, liv, lv
- integer iv(liv), uiparm(1)
- double precision d(n), x(n), v(lv), urparm(1)
-c dimension v(71 + n*(n+15)/2), uiparm(*), urparm(*)
- external calcf, calcg, ufparm
-c
-c *** purpose ***
-c
-c this routine interacts with subroutine sumit in an attempt
-c to find an n-vector x* that minimizes the (unconstrained)
-c objective function computed by calcf. (often the x* found is
-c a local minimizer rather than a global one.)
-c
-c-------------------------- parameter usage --------------------------
-c
-c n........ (input) the number of variables on which f depends, i.e.,
-c the number of components in x.
-c d........ (input/output) a scale vector such that d(i)*x(i),
-c i = 1,2,...,n, are all in comparable units.
-c d can strongly affect the behavior of sumsl.
-c finding the best choice of d is generally a trial-
-c and-error process. choosing d so that d(i)*x(i)
-c has about the same value for all i often works well.
-c the defaults provided by subroutine deflt (see i
-c below) require the caller to supply d.
-c x........ (input/output) before (initially) calling sumsl, the call-
-c er should set x to an initial guess at x*. when
-c sumsl returns, x contains the best point so far
-c found, i.e., the one that gives the least value so
-c far seen for f(x).
-c calcf.... (input) a subroutine that, given x, computes f(x). calcf
-c must be declared external in the calling program.
-c it is invoked by
-c call calcf(n, x, nf, f, uiparm, urparm, ufparm)
-c when calcf is called, nf is the invocation
-c count for calcf. nf is included for possible use
-c with calcg. if x is out of bounds (e.g., if it
-c would cause overflow in computing f(x)), then calcf
-c should set nf to 0. this will cause a shorter step
-c to be attempted. (if x is in bounds, then calcf
-c should not change nf.) the other parameters are as
-c described above and below. calcf should not change
-c n, p, or x.
-c calcg.... (input) a subroutine that, given x, computes g(x), the gra-
-c dient of f at x. calcg must be declared external in
-c the calling program. it is invoked by
-c call calcg(n, x, nf, g, uiparm, urparm, ufaprm)
-c when calcg is called, nf is the invocation
-c count for calcf at the time f(x) was evaluated. the
-c x passed to calcg is usually the one passed to calcf
-c on either its most recent invocation or the one
-c prior to it. if calcf saves intermediate results
-c for use by calcg, then it is possible to tell from
-c nf whether they are valid for the current x (or
-c which copy is valid if two copies are kept). if g
-c cannot be computed at x, then calcg should set nf to
-c 0. in this case, sumsl will return with iv(1) = 65.
-c (if g can be computed at x, then calcg should not
-c changed nf.) the other parameters to calcg are as
-c described above and below. calcg should not change
-c n or x.
-c iv....... (input/output) an integer value array of length liv (see
-c below) that helps control the sumsl algorithm and
-c that is used to store various intermediate quanti-
-c ties. of particular interest are the initialization/
-c return code iv(1) and the entries in iv that control
-c printing and limit the number of iterations and func-
-c tion evaluations. see the section on iv input
-c values below.
-c liv...... (input) length of iv array. must be at least 60. if li
-c is too small, then sumsl returns with iv(1) = 15.
-c when sumsl returns, the smallest allowed value of
-c liv is stored in iv(lastiv) -- see the section on
-c iv output values below. (this is intended for use
-c with extensions of sumsl that handle constraints.)
-c lv....... (input) length of v array. must be at least 71+n*(n+15)/2.
-c (at least 77+n*(n+17)/2 for smsno, at least
-c 78+n*(n+12) for humsl). if lv is too small, then
-c sumsl returns with iv(1) = 16. when sumsl returns,
-c the smallest allowed value of lv is stored in
-c iv(lastv) -- see the section on iv output values
-c below.
-c v........ (input/output) a floating-point value array of length l
-c (see below) that helps control the sumsl algorithm
-c and that is used to store various intermediate
-c quantities. of particular interest are the entries
-c in v that limit the length of the first step
-c attempted (lmax0) and specify convergence tolerances
-c (afctol, lmaxs, rfctol, sctol, xctol, xftol).
-c uiparm... (input) user integer parameter array passed without change
-c to calcf and calcg.
-c urparm... (input) user floating-point parameter array passed without
-c change to calcf and calcg.
-c ufparm... (input) user external subroutine or function passed without
-c change to calcf and calcg.
-c
-c *** iv input values (from subroutine deflt) ***
-c
-c iv(1)... on input, iv(1) should have a value between 0 and 14......
-c 0 and 12 mean this is a fresh start. 0 means that
-c deflt(2, iv, liv, lv, v)
-c is to be called to provide all default values to iv and
-c v. 12 (the value that deflt assigns to iv(1)) means the
-c caller has already called deflt and has possibly changed
-c some iv and/or v entries to non-default values.
-c 13 means deflt has been called and that sumsl (and
-c sumit) should only do their storage allocation. that is,
-c they should set the output components of iv that tell
-c where various subarrays arrays of v begin, such as iv(g)
-c (and, for humsl and humit only, iv(dtol)), and return.
-c 14 means that a storage has been allocated (by a call
-c with iv(1) = 13) and that the algorithm should be
-c started. when called with iv(1) = 13, sumsl returns
-c iv(1) = 14 unless liv or lv is too small (or n is not
-c positive). default = 12.
-c iv(inith).... iv(25) tells whether the hessian approximation h should
-c be initialized. 1 (the default) means sumit should
-c initialize h to the diagonal matrix whose i-th diagonal
-c element is d(i)**2. 0 means the caller has supplied a
-c cholesky factor l of the initial hessian approximation
-c h = l*(l**t) in v, starting at v(iv(lmat)) = v(iv(42))
-c (and stored compactly by rows). note that iv(lmat) may
-c be initialized by calling sumsl with iv(1) = 13 (see
-c the iv(1) discussion above). default = 1.
-c iv(mxfcal)... iv(17) gives the maximum number of function evaluations
-c (calls on calcf) allowed. if this number does not suf-
-c fice, then sumsl returns with iv(1) = 9. default = 200.
-c iv(mxiter)... iv(18) gives the maximum number of iterations allowed.
-c it also indirectly limits the number of gradient evalua-
-c tions (calls on calcg) to iv(mxiter) + 1. if iv(mxiter)
-c iterations do not suffice, then sumsl returns with
-c iv(1) = 10. default = 150.
-c iv(outlev)... iv(19) controls the number and length of iteration sum-
-c mary lines printed (by itsum). iv(outlev) = 0 means do
-c not print any summary lines. otherwise, print a summary
-c line after each abs(iv(outlev)) iterations. if iv(outlev)
-c is positive, then summary lines of length 78 (plus carri-
-c age control) are printed, including the following... the
-c iteration and function evaluation counts, f = the current
-c function value, relative difference in function values
-c achieved by the latest step (i.e., reldf = (f0-v(f))/f01,
-c where f01 is the maximum of abs(v(f)) and abs(v(f0)) and
-c v(f0) is the function value from the previous itera-
-c tion), the relative function reduction predicted for the
-c step just taken (i.e., preldf = v(preduc) / f01, where
-c v(preduc) is described below), the scaled relative change
-c in x (see v(reldx) below), the step parameter for the
-c step just taken (stppar = 0 means a full newton step,
-c between 0 and 1 means a relaxed newton step, between 1
-c and 2 means a double dogleg step, greater than 2 means
-c a scaled down cauchy step -- see subroutine dbldog), the
-c 2-norm of the scale vector d times the step just taken
-c (see v(dstnrm) below), and npreldf, i.e.,
-c v(nreduc)/f01, where v(nreduc) is described below -- if
-c npreldf is positive, then it is the relative function
-c reduction predicted for a newton step (one with
-c stppar = 0). if npreldf is negative, then it is the
-c negative of the relative function reduction predicted
-c for a step computed with step bound v(lmaxs) for use in
-c testing for singular convergence.
-c if iv(outlev) is negative, then lines of length 50
-c are printed, including only the first 6 items listed
-c above (through reldx).
-c default = 1.
-c iv(parprt)... iv(20) = 1 means print any nondefault v values on a
-c fresh start or any changed v values on a restart.
-c iv(parprt) = 0 means skip this printing. default = 1.
-c iv(prunit)... iv(21) is the output unit number on which all printing
-c is done. iv(prunit) = 0 means suppress all printing.
-c default = standard output unit (unit 6 on most systems).
-c iv(solprt)... iv(22) = 1 means print out the value of x returned (as
-c well as the gradient and the scale vector d).
-c iv(solprt) = 0 means skip this printing. default = 1.
-c iv(statpr)... iv(23) = 1 means print summary statistics upon return-
-c ing. these consist of the function value, the scaled
-c relative change in x caused by the most recent step (see
-c v(reldx) below), the number of function and gradient
-c evaluations (calls on calcf and calcg), and the relative
-c function reductions predicted for the last step taken and
-c for a newton step (or perhaps a step bounded by v(lmaxs)
-c -- see the descriptions of preldf and npreldf under
-c iv(outlev) above).
-c iv(statpr) = 0 means skip this printing.
-c iv(statpr) = -1 means skip this printing as well as that
-c of the one-line termination reason message. default = 1.
-c iv(x0prt).... iv(24) = 1 means print the initial x and scale vector d
-c (on a fresh start only). iv(x0prt) = 0 means skip this
-c printing. default = 1.
-c
-c *** (selected) iv output values ***
-c
-c iv(1)........ on output, iv(1) is a return code....
-c 3 = x-convergence. the scaled relative difference (see
-c v(reldx)) between the current parameter vector x and
-c a locally optimal parameter vector is very likely at
-c most v(xctol).
-c 4 = relative function convergence. the relative differ-
-c ence between the current function value and its lo-
-c cally optimal value is very likely at most v(rfctol).
-c 5 = both x- and relative function convergence (i.e., the
-c conditions for iv(1) = 3 and iv(1) = 4 both hold).
-c 6 = absolute function convergence. the current function
-c value is at most v(afctol) in absolute value.
-c 7 = singular convergence. the hessian near the current
-c iterate appears to be singular or nearly so, and a
-c step of length at most v(lmaxs) is unlikely to yield
-c a relative function decrease of more than v(sctol).
-c 8 = false convergence. the iterates appear to be converg-
-c ing to a noncritical point. this may mean that the
-c convergence tolerances (v(afctol), v(rfctol),
-c v(xctol)) are too small for the accuracy to which
-c the function and gradient are being computed, that
-c there is an error in computing the gradient, or that
-c the function or gradient is discontinuous near x.
-c 9 = function evaluation limit reached without other con-
-c vergence (see iv(mxfcal)).
-c 10 = iteration limit reached without other convergence
-c (see iv(mxiter)).
-c 11 = stopx returned .true. (external interrupt). see the
-c usage notes below.
-c 14 = storage has been allocated (after a call with
-c iv(1) = 13).
-c 17 = restart attempted with n changed.
-c 18 = d has a negative component and iv(dtype) .le. 0.
-c 19...43 = v(iv(1)) is out of range.
-c 63 = f(x) cannot be computed at the initial x.
-c 64 = bad parameters passed to assess (which should not
-c occur).
-c 65 = the gradient could not be computed at x (see calcg
-c above).
-c 67 = bad first parameter to deflt.
-c 80 = iv(1) was out of range.
-c 81 = n is not positive.
-c iv(g)........ iv(28) is the starting subscript in v of the current
-c gradient vector (the one corresponding to x).
-c iv(lastiv)... iv(44) is the least acceptable value of liv. (it is
-c only set if liv is at least 44.)
-c iv(lastv).... iv(45) is the least acceptable value of lv. (it is
-c only set if liv is large enough, at least iv(lastiv).)
-c iv(nfcall)... iv(6) is the number of calls so far made on calcf (i.e.,
-c function evaluations).
-c iv(ngcall)... iv(30) is the number of gradient evaluations (calls on
-c calcg).
-c iv(niter).... iv(31) is the number of iterations performed.
-c
-c *** (selected) v input values (from subroutine deflt) ***
-c
-c v(bias)..... v(43) is the bias parameter used in subroutine dbldog --
-c see that subroutine for details. default = 0.8.
-c v(afctol)... v(31) is the absolute function convergence tolerance.
-c if sumsl finds a point where the function value is less
-c than v(afctol) in absolute value, and if sumsl does not
-c return with iv(1) = 3, 4, or 5, then it returns with
-c iv(1) = 6. this test can be turned off by setting
-c v(afctol) to zero. default = max(10**-20, machep**2),
-c where machep is the unit roundoff.
-c v(dinit).... v(38), if nonnegative, is the value to which the scale
-c vector d is initialized. default = -1.
-c v(lmax0).... v(35) gives the maximum 2-norm allowed for d times the
-c very first step that sumsl attempts. this parameter can
-c markedly affect the performance of sumsl.
-c v(lmaxs).... v(36) is used in testing for singular convergence -- if
-c the function reduction predicted for a step of length
-c bounded by v(lmaxs) is at most v(sctol) * abs(f0), where
-c f0 is the function value at the start of the current
-c iteration, and if sumsl does not return with iv(1) = 3,
-c 4, 5, or 6, then it returns with iv(1) = 7. default = 1.
-c v(rfctol)... v(32) is the relative function convergence tolerance.
-c if the current model predicts a maximum possible function
-c reduction (see v(nreduc)) of at most v(rfctol)*abs(f0)
-c at the start of the current iteration, where f0 is the
-c then current function value, and if the last step attempt-
-c ed achieved no more than twice the predicted function
-c decrease, then sumsl returns with iv(1) = 4 (or 5).
-c default = max(10**-10, machep**(2/3)), where machep is
-c the unit roundoff.
-c v(sctol).... v(37) is the singular convergence tolerance -- see the
-c description of v(lmaxs) above.
-c v(tuner1)... v(26) helps decide when to check for false convergence.
-c this is done if the actual function decrease from the
-c current step is no more than v(tuner1) times its predict-
-c ed value. default = 0.1.
-c v(xctol).... v(33) is the x-convergence tolerance. if a newton step
-c (see v(nreduc)) is tried that has v(reldx) .le. v(xctol)
-c and if this step yields at most twice the predicted func-
-c tion decrease, then sumsl returns with iv(1) = 3 (or 5).
-c (see the description of v(reldx) below.)
-c default = machep**0.5, where machep is the unit roundoff.
-c v(xftol).... v(34) is the false convergence tolerance. if a step is
-c tried that gives no more than v(tuner1) times the predict-
-c ed function decrease and that has v(reldx) .le. v(xftol),
-c and if sumsl does not return with iv(1) = 3, 4, 5, 6, or
-c 7, then it returns with iv(1) = 8. (see the description
-c of v(reldx) below.) default = 100*machep, where
-c machep is the unit roundoff.
-c v(*)........ deflt supplies to v a number of tuning constants, with
-c which it should ordinarily be unnecessary to tinker. see
-c section 17 of version 2.2 of the nl2sol usage summary
-c (i.e., the appendix to ref. 1) for details on v(i),
-c i = decfac, incfac, phmnfc, phmxfc, rdfcmn, rdfcmx,
-c tuner2, tuner3, tuner4, tuner5.
-c
-c *** (selected) v output values ***
-c
-c v(dgnorm)... v(1) is the 2-norm of (diag(d)**-1)*g, where g is the
-c most recently computed gradient.
-c v(dstnrm)... v(2) is the 2-norm of diag(d)*step, where step is the
-c current step.
-c v(f)........ v(10) is the current function value.
-c v(f0)....... v(13) is the function value at the start of the current
-c iteration.
-c v(nreduc)... v(6), if positive, is the maximum function reduction
-c possible according to the current model, i.e., the func-
-c tion reduction predicted for a newton step (i.e.,
-c step = -h**-1 * g, where g is the current gradient and
-c h is the current hessian approximation).
-c if v(nreduc) is negative, then it is the negative of
-c the function reduction predicted for a step computed with
-c a step bound of v(lmaxs) for use in testing for singular
-c convergence.
-c v(preduc)... v(7) is the function reduction predicted (by the current
-c quadratic model) for the current step. this (divided by
-c v(f0)) is used in testing for relative function
-c convergence.
-c v(reldx).... v(17) is the scaled relative change in x caused by the
-c current step, computed as
-c max(abs(d(i)*(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 where x = x0 + step.
-c
-c------------------------------- notes -------------------------------
-c
-c *** algorithm notes ***
-c
-c this routine uses a hessian approximation computed from the
-c bfgs update (see ref 3). only a cholesky factor of the hessian
-c approximation is stored, and this is updated using ideas from
-c ref. 4. steps are computed by the double dogleg scheme described
-c in ref. 2. the steps are assessed as in ref. 1.
-c
-c *** usage notes ***
-c
-c after a return with iv(1) .le. 11, it is possible to restart,
-c i.e., to change some of the iv and v input values described above
-c and continue the algorithm from the point where it was interrupt-
-c ed. iv(1) should not be changed, nor should any entries of i
-c and v other than the input values (those supplied by deflt).
-c those who do not wish to write a calcg which computes the
-c gradient analytically should call smsno rather than sumsl.
-c smsno uses finite differences to compute an approximate gradient.
-c those who would prefer to provide f and g (the function and
-c gradient) by reverse communication rather than by writing subrou-
-c tines calcf and calcg may call on sumit directly. see the com-
-c ments at the beginning of sumit.
-c those who use sumsl interactively may wish to supply their
-c own stopx function, which should return .true. if the break key
-c has been pressed since stopx was last invoked. this makes it
-c possible to externally interrupt sumsl (which will return with
-c iv(1) = 11 if stopx returns .true.).
-c storage for g is allocated at the end of v. thus the caller
-c may make v longer than specified above and may allow calcg to use
-c elements of g beyond the first n as scratch storage.
-c
-c *** portability notes ***
-c
-c the sumsl distribution tape contains both single- and double-
-c precision versions of the sumsl source code, so it should be un-
-c necessary to change precisions.
-c only the functions imdcon and rmdcon contain machine-dependent
-c constants. to change from one machine to another, it should
-c suffice to change the (few) relevant lines in these functions.
-c intrinsic functions are explicitly declared. on certain com-
-c puters (e.g. univac), it may be necessary to comment out these
-c declarations. so that this may be done automatically by a simple
-c program, such declarations are preceded by a comment having c/+
-c in columns 1-3 and blanks in columns 4-72 and are followed by
-c a comment having c/ in columns 1 and 2 and blanks in columns 3-72.
-c the sumsl source code is expressed in 1966 ansi standard
-c fortran. it may be converted to fortran 77 by commenting out all
-c lines that fall between a line having c/6 in columns 1-3 and a
-c line having c/7 in columns 1-3 and by removing (i.e., replacing
-c by a blank) the c in column 1 of the lines that follow the c/7
-c line and precede a line having c/ in columns 1-2 and blanks in
-c columns 3-72. these changes convert some data statements into
-c parameter statements, convert some variables from real to
-c character*4, and make the data statements that initialize these
-c variables use character strings delimited by primes instead
-c of hollerith constants. (such variables and data statements
-c appear only in modules itsum and parck. parameter statements
-c appear nearly everywhere.) these changes also add save state-
-c ments for variables given machine-dependent constants by rmdcon.
-c
-c *** references ***
-c
-c 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), algorithm 573 --
-c an adaptive nonlinear least-squares algorithm, acm trans.
-c math. software 7, pp. 369-383.
-c
-c 2. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
-c mization algorithms which use function and gradient
-c values, j. optim. theory applic. 28, pp. 453-482.
-c
-c 3. dennis, j.e., and more, j.j. (1977), quasi-newton methods, motiva-
-c tion and theory, siam rev. 19, pp. 46-89.
-c
-c 4. goldfarb, d. (1976), factorized variable metric methods for uncon-
-c strained optimization, math. comput. 30, pp. 796-811.
-c
-c *** general ***
-c
-c coded by david m. gay (winter 1980). revised summer 1982.
-c this subroutine was written in connection with research
-c supported in part by the national science foundation under
-c grants mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989,
-c and mcs-7906671.
-c.
-c
-c---------------------------- declarations ---------------------------
-c
- external deflt, sumit
-c
-c deflt... supplies default iv and v input components.
-c sumit... reverse-communication routine that carries out sumsl algo-
-c rithm.
-c
- integer g1, iv1, nf
- double precision f
-c
-c *** subscripts for iv ***
-c
- integer nextv, nfcall, nfgcal, g, toobig, vneed
-c
-c/6
-c data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, toobig/2/, vneed/4/
-c/7
- parameter (nextv=47, nfcall=6, nfgcal=7, g=28, toobig=2, vneed=4)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
- iv1 = iv(1)
- if (iv1 .eq. 12 .or. iv1 .eq. 13) iv(vneed) = iv(vneed) + n
- if (iv1 .eq. 14) go to 10
- if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
- g1 = 1
- if (iv1 .eq. 12) iv(1) = 13
- go to 20
-c
- 10 g1 = iv(g)
-c
- 20 call sumit(d, f, v(g1), iv, 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 calcg(n, x, iv(nfgcal), v(g1), uiparm, urparm, ufparm)
- go to 20
-c
- 50 if (iv(1) .ne. 14) go to 999
-c
-c *** storage allocation
-c
- iv(g) = iv(nextv)
- iv(nextv) = iv(g) + n
- if (iv1 .ne. 13) go to 10
-c
- 999 return
-c *** last card of sumsl follows ***
- end
- subroutine sumit(d, fx, g, iv, liv, lv, n, v, x)
-c
-c *** carry out sumsl (unconstrained minimization) iterations, using
-c *** double-dogleg/bfgs steps.
-c
-c *** parameter declarations ***
-c
- integer liv, lv, n
- integer iv(liv)
- double precision d(n), fx, g(n), v(lv), x(n)
-c
-c-------------------------- parameter usage --------------------------
-c
-c d.... scale vector.
-c fx... function value.
-c g.... gradient vector.
-c iv... integer value array.
-c liv.. length of iv (at least 60).
-c lv... length of v (at least 71 + n*(n+13)/2).
-c n.... number of variables (components in x and g).
-c v.... floating-point value array.
-c x.... vector of parameters to be optimized.
-c
-c *** discussion ***
-c
-c parameters iv, n, v, and x are the same as the corresponding
-c ones to sumsl (which see), except that v can be shorter (since
-c the part of v that sumsl uses for storing g is not needed).
-c moreover, compared with sumsl, 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 sumsl (and smsno), is not referenced by
-c sumit or the subroutines it calls.
-c fx and g need not have been initialized when sumit is called
-c with iv(1) = 12, 13, or 14.
-c
-c iv(1) = 1 means the caller should set fx to f(x), the function value
-c at x, and call sumit again, having changed none of the
-c other parameters. an exception occurs if f(x) cannot be
-c (e.g. if overflow would occur), which may happen because
-c of an oversized step. in this case the caller should set
-c iv(toobig) = iv(2) to 1, which will cause sumit to ig-
-c nore fx and try a smaller step. the parameter nf that
-c sumsl passes to calcf (for possible use by calcg) is a
-c copy of iv(nfcall) = iv(6).
-c iv(1) = 2 means the caller should set g to g(x), the gradient vector
-c of f at x, and call sumit again, having changed none of
-c the other parameters except possibly the scale vector d
-c when iv(dtype) = 0. the parameter nf that sumsl passes
-c to calcg is iv(nfgcal) = iv(7). if g(x) cannot be
-c evaluated, then the caller may set iv(nfgcal) to 0, in
-c which case sumit will return with iv(1) = 65.
-c.
-c *** general ***
-c
-c coded by david m. gay (december 1979). 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 for references.)
-c
-c+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
-c
-c *** local variables ***
-c
- integer dg1, dummy, g01, i, k, l, lstgst, nwtst1, step1,
- 1 temp1, w, x01, z
- double precision t
-c
-c *** constants ***
-c
- double precision half, negone, one, onep2, zero
-c
-c *** no intrinsic functions ***
-c
-c *** external functions and subroutines ***
-c
- external assst, dbdog, deflt, dotprd, itsum, litvmu, livmul,
- 1 ltvmul, lupdat, lvmul, parck, reldst, stopx, vaxpy,
- 2 vcopy, vscopy, vvmulp, v2norm, wzbfgs
- logical stopx
- double precision dotprd, reldst, v2norm
-c
-c assst.... assesses candidate step.
-c dbdog.... computes double-dogleg (candidate) step.
-c deflt.... supplies default iv and v input components.
-c dotprd... returns inner product of two vectors.
-c itsum.... prints iteration summary and info on initial and final x.
-c litvmu... multiplies inverse transpose of lower triangle times vector.
-c livmul... multiplies inverse of lower triangle times vector.
-c ltvmul... multiplies transpose of lower triangle times vector.
-c lupdt.... updates cholesky factor of hessian approximation.
-c lvmul.... multiplies lower triangle times vector.
-c parck.... checks validity of input iv and v values.
-c reldst... computes v(reldx) = relative step size.
-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 vvmulp... multiplies vector by vector raised to power (componentwise).
-c v2norm... returns the 2-norm of a vector.
-c wzbfgs... computes w and z for lupdat corresponding to bfgs update.
-c
-c *** subscripts for iv and v ***
-c
- integer afctol
- integer cnvcod, dg, dgnorm, dinit, dstnrm, dst0, f, f0, fdif,
- 1 gthg, gtstep, g0, incfac, inith, irc, kagqt, lmat, lmax0,
- 2 lmaxs, mode, model, mxfcal, mxiter, nextv, nfcall, nfgcal,
- 3 ngcall, niter, nreduc, nwtstp, preduc, radfac, radinc,
- 4 radius, rad0, reldx, restor, step, stglim, stlstg, toobig,
- 5 tuner4, tuner5, vneed, xirc, x0
-c
-c *** iv subscript values ***
-c
-c/6
-c data cnvcod/55/, dg/37/, g0/48/, inith/25/, irc/29/, kagqt/33/,
-c 1 mode/35/, model/5/, mxfcal/17/, mxiter/18/, nfcall/6/,
-c 2 nfgcal/7/, ngcall/30/, niter/31/, nwtstp/34/, radinc/8/,
-c 3 restor/9/, step/40/, stglim/11/, stlstg/41/, toobig/2/,
-c 4 vneed/4/, xirc/13/, x0/43/
-c/7
- parameter (cnvcod=55, dg=37, g0=48, inith=25, irc=29, kagqt=33,
- 1 mode=35, model=5, mxfcal=17, mxiter=18, nfcall=6,
- 2 nfgcal=7, ngcall=30, niter=31, nwtstp=34, radinc=8,
- 3 restor=9, step=40, stglim=11, stlstg=41, toobig=2,
- 4 vneed=4, xirc=13, x0=43)
-c/
-c
-c *** v subscript values ***
-c
-c/6
-c data afctol/31/
-c data dgnorm/1/, dinit/38/, dstnrm/2/, dst0/3/, f/10/, f0/13/,
-c 1 fdif/11/, gthg/44/, gtstep/4/, incfac/23/, lmat/42/,
-c 2 lmax0/35/, lmaxs/36/, nextv/47/, nreduc/6/, preduc/7/,
-c 3 radfac/16/, radius/8/, rad0/9/, reldx/17/, tuner4/29/,
-c 4 tuner5/30/
-c/7
- parameter (afctol=31)
- parameter (dgnorm=1, dinit=38, dstnrm=2, dst0=3, f=10, f0=13,
- 1 fdif=11, gthg=44, gtstep=4, incfac=23, lmat=42,
- 2 lmax0=35, lmaxs=36, nextv=47, nreduc=6, preduc=7,
- 3 radfac=16, radius=8, rad0=9, reldx=17, tuner4=29,
- 4 tuner5=30)
-c/
-c
-c/6
-c data half/0.5d+0/, negone/-1.d+0/, one/1.d+0/, onep2/1.2d+0/,
-c 1 zero/0.d+0/
-c/7
- parameter (half=0.5d+0, negone=-1.d+0, one=1.d+0, onep2=1.2d+0,
- 1 zero=0.d+0)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
-C Following SAVE statement inserted.
- save l
- i = iv(1)
- if (i .eq. 1) go to 50
- if (i .eq. 2) go to 60
-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+13)/2
- call parck(2, d, iv, liv, lv, n, v)
- i = iv(1) - 2
- if (i .gt. 12) go to 999
- go to (180, 180, 180, 180, 180, 180, 120, 90, 120, 10, 10, 20), i
-c
-c *** storage allocation ***
-c
-10 l = iv(lmat)
- iv(x0) = l + n*(n+1)/2
- iv(step) = iv(x0) + n
- iv(stlstg) = iv(step) + n
- iv(g0) = iv(stlstg) + n
- iv(nwtstp) = iv(g0) + n
- iv(dg) = iv(nwtstp) + n
- iv(nextv) = iv(dg) + n
- 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
- if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
- if (iv(inith) .ne. 1) go to 40
-c
-c *** set the initial hessian approximation to diag(d)**-2 ***
-c
- l = iv(lmat)
- call vscopy(n*(n+1)/2, v(l), zero)
- k = l - 1
- do 30 i = 1, n
- k = k + i
- t = d(i)
- if (t .le. zero) t = one
- v(k) = t
- 30 continue
-c
-c *** compute initial function value ***
-c
- 40 iv(1) = 1
- go to 999
-c
- 50 v(f) = fx
- if (iv(mode) .ge. 0) go to 180
- iv(1) = 2
- if (iv(toobig) .eq. 0) go to 999
- iv(1) = 63
- go to 300
-c
-c *** make sure gradient could be computed ***
-c
- 60 if (iv(nfgcal) .ne. 0) go to 70
- iv(1) = 65
- go to 300
-c
- 70 dg1 = iv(dg)
- call vvmulp(n, v(dg1), g, d, -1)
- v(dgnorm) = v2norm(n, v(dg1))
-c
-c *** test norm of gradient ***
-c
- if (v(dgnorm) .gt. v(afctol)) go to 75
- iv(irc) = 10
- iv(cnvcod) = iv(irc) - 4
-c
- 75 if (iv(cnvcod) .ne. 0) go to 290
- if (iv(mode) .eq. 0) go to 250
-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
- 80 call itsum(d, g, iv, liv, lv, n, v, x)
- 90 k = iv(niter)
- if (k .lt. iv(mxiter)) go to 100
- iv(1) = 10
- go to 300
-c
-c *** update radius ***
-c
- 100 iv(niter) = k + 1
- if(k.gt.0)v(radius) = v(radfac) * v(dstnrm)
-c
-c *** initialize for start of next iteration ***
-c
- g01 = iv(g0)
- x01 = iv(x0)
- v(f0) = v(f)
- iv(irc) = 4
- iv(kagqt) = -1
-c
-c *** copy x to x0, g to g0 ***
-c
- call vcopy(n, v(x01), x)
- call vcopy(n, v(g01), g)
-c
-c *** check stopx and function evaluation limit ***
-c
-C AL 4/30/95
- dummy=iv(nfcall)
- 110 if (.not. stopx(dummy)) go to 130
- iv(1) = 11
- go to 140
-c
-c *** come here when restarting after func. eval. limit or stopx.
-c
- 120 if (v(f) .ge. v(f0)) go to 130
- v(radfac) = one
- k = iv(niter)
- go to 100
-c
- 130 if (iv(nfcall) .lt. iv(mxfcal)) go to 150
- iv(1) = 9
- 140 if (v(f) .ge. v(f0)) go to 300
-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 240
-c
-c. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
-c
- 150 step1 = iv(step)
- dg1 = iv(dg)
- nwtst1 = iv(nwtstp)
- if (iv(kagqt) .ge. 0) go to 160
- l = iv(lmat)
- call livmul(n, v(nwtst1), v(l), g)
- v(nreduc) = half * dotprd(n, v(nwtst1), v(nwtst1))
- call litvmu(n, v(nwtst1), v(l), v(nwtst1))
- call vvmulp(n, v(step1), v(nwtst1), d, 1)
- v(dst0) = v2norm(n, v(step1))
- call vvmulp(n, v(dg1), v(dg1), d, -1)
- call ltvmul(n, v(step1), v(l), v(dg1))
- v(gthg) = v2norm(n, v(step1))
- iv(kagqt) = 0
- 160 call dbdog(v(dg1), lv, n, v(nwtst1), v(step1), v)
- if (iv(irc) .eq. 6) go to 180
-c
-c *** check whether evaluating f(x0 + step) looks worthwhile ***
-c
- if (v(dstnrm) .le. zero) go to 180
- if (iv(irc) .ne. 5) go to 170
- if (v(radfac) .le. one) go to 170
- if (v(preduc) .le. onep2 * v(fdif)) go to 180
-c
-c *** compute f(x0 + step) ***
-c
- 170 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
- 180 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 190
- 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
- 190 k = iv(irc)
- go to (200,230,230,230,200,210,220,220,220,220,220,220,280,250), k
-c
-c *** recompute step with changed radius ***
-c
- 200 v(radius) = v(radfac) * v(dstnrm)
- go to 110
-c
-c *** compute step of length v(lmaxs) for singular convergence test.
-c
- 210 v(radius) = v(lmaxs)
- go to 150
-c
-c *** convergence or false convergence ***
-c
- 220 iv(cnvcod) = k - 4
- if (v(f) .ge. v(f0)) go to 290
- if (iv(xirc) .eq. 14) go to 290
- iv(xirc) = 14
-c
-c. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
-c
- 230 if (iv(irc) .ne. 3) go to 240
- step1 = iv(step)
- temp1 = iv(stlstg)
-c
-c *** set temp1 = hessian * step for use in gradient tests ***
-c
- l = iv(lmat)
- call ltvmul(n, v(temp1), v(l), v(step1))
- call lvmul(n, v(temp1), v(l), v(temp1))
-c
-c *** compute gradient ***
-c
- 240 iv(ngcall) = iv(ngcall) + 1
- iv(1) = 2
- go to 999
-c
-c *** initializations -- g0 = g - g0, etc. ***
-c
- 250 g01 = iv(g0)
- call vaxpy(n, v(g01), negone, v(g01), g)
- step1 = iv(step)
- temp1 = iv(stlstg)
- if (iv(irc) .ne. 3) go to 270
-c
-c *** set v(radfac) by gradient tests ***
-c
-c *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
-c
- call vaxpy(n, v(temp1), negone, v(g01), v(temp1))
- call vvmulp(n, v(temp1), v(temp1), d, -1)
-c
-c *** do gradient tests ***
-c
- if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4))
- 1 go to 260
- if (dotprd(n, g, v(step1))
- 1 .ge. v(gtstep) * v(tuner5)) go to 270
- 260 v(radfac) = v(incfac)
-c
-c *** update h, loop ***
-c
- 270 w = iv(nwtstp)
- z = iv(x0)
- l = iv(lmat)
- call wzbfgs(v(l), n, v(step1), v(w), v(g01), v(z))
-c
-c ** use the n-vectors starting at v(step1) and v(g01) for scratch..
- call lupdat(v(temp1), v(step1), v(l), v(g01), v(l), n, v(w), v(z))
- iv(1) = 2
- go to 80
-c
-c. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
-c
-c *** bad parameters to assess ***
-c
- 280 iv(1) = 64
- go to 300
-c
-c *** print summary of final iteration and other requested items ***
-c
- 290 iv(1) = iv(cnvcod)
- iv(cnvcod) = 0
- 300 call itsum(d, g, iv, liv, lv, n, v, x)
-c
- 999 return
-c
-c *** last line of sumit follows ***
- end
- subroutine dbdog(dig, lv, n, nwtstp, step, v)
-c
-c *** compute double dogleg step ***
-c
-c *** parameter declarations ***
-c
- integer lv, n
- double precision dig(n), nwtstp(n), step(n), v(lv)
-c
-c *** purpose ***
-c
-c this subroutine computes a candidate step (for use in an uncon-
-c strained minimization code) by the double dogleg algorithm of
-c dennis and mei (ref. 1), which is a variation on powell*s dogleg
-c scheme (ref. 2, p. 95).
-c
-c-------------------------- parameter usage --------------------------
-c
-c dig (input) diag(d)**-2 * g -- see algorithm notes.
-c g (input) the current gradient vector.
-c lv (input) length of v.
-c n (input) number of components in dig, g, nwtstp, and step.
-c nwtstp (input) negative newton step -- see algorithm notes.
-c step (output) the computed step.
-c v (i/o) values array, the following components of which are
-c used here...
-c v(bias) (input) bias for relaxed newton step, which is v(bias) of
-c the way from the full newton to the fully relaxed newton
-c step. recommended value = 0.8 .
-c v(dgnorm) (input) 2-norm of diag(d)**-1 * g -- see algorithm notes.
-c v(dstnrm) (output) 2-norm of diag(d) * step, which is v(radius)
-c unless v(stppar) = 0 -- see algorithm notes.
-c v(dst0) (input) 2-norm of diag(d) * nwtstp -- see algorithm notes.
-c v(grdfac) (output) the coefficient of dig in the step returned --
-c step(i) = v(grdfac)*dig(i) + v(nwtfac)*nwtstp(i).
-c v(gthg) (input) square-root of (dig**t) * (hessian) * dig -- see
-c algorithm notes.
-c v(gtstep) (output) inner product between g and step.
-c v(nreduc) (output) function reduction predicted for the full newton
-c step.
-c v(nwtfac) (output) the coefficient of nwtstp in the step returned --
-c see v(grdfac) above.
-c v(preduc) (output) function reduction predicted for the step returned.
-c v(radius) (input) the trust region radius. d times the step returned
-c has 2-norm v(radius) unless v(stppar) = 0.
-c v(stppar) (output) code telling how step was computed... 0 means a
-c full newton step. between 0 and 1 means v(stppar) of the
-c way from the newton to the relaxed newton step. between
-c 1 and 2 means a true double dogleg step, v(stppar) - 1 of
-c the way from the relaxed newton to the cauchy step.
-c greater than 2 means 1 / (v(stppar) - 1) times the cauchy
-c step.
-c
-c------------------------------- notes -------------------------------
-c
-c *** algorithm notes ***
-c
-c let g and h be the current gradient and hessian approxima-
-c tion respectively and let d be the current scale vector. this
-c routine assumes dig = diag(d)**-2 * g and nwtstp = h**-1 * g.
-c the step computed is the same one would get by replacing g and h
-c by diag(d)**-1 * g and diag(d)**-1 * h * diag(d)**-1,
-c computing step, and translating step back to the original
-c variables, i.e., premultiplying it by diag(d)**-1.
-c
-c *** references ***
-c
-c 1. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
-c mization algorithms which use function and gradient
-c values, j. optim. theory applic. 28, pp. 453-482.
-c 2. powell, m.j.d. (1970), a hybrid method for non-linear equations,
-c in numerical methods for non-linear equations, edited by
-c p. rabinowitz, gordon and breach, london.
-c
-c *** general ***
-c
-c coded by david m. gay.
-c this subroutine was written in connection with research supported
-c by the national science foundation under grants mcs-7600324 and
-c mcs-7906671.
-c
-c------------------------ external quantities ------------------------
-c
-c *** functions and subroutines called ***
-c
- external dotprd, v2norm
- double precision dotprd, v2norm
-c
-c dotprd... returns inner product of two vectors.
-c v2norm... returns 2-norm of a vector.
-c
-c *** intrinsic functions ***
-c/+
- double precision dsqrt
-c/
-c-------------------------- local variables --------------------------
-c
- integer i
- double precision cfact, cnorm, ctrnwt, ghinvg, femnsq, gnorm,
- 1 nwtnrm, relax, rlambd, t, t1, t2
- double precision half, one, two, zero
-c
-c *** v subscripts ***
-c
- integer bias, dgnorm, dstnrm, dst0, grdfac, gthg, gtstep,
- 1 nreduc, nwtfac, preduc, radius, stppar
-c
-c *** data initializations ***
-c
-c/6
-c data half/0.5d+0/, one/1.d+0/, two/2.d+0/, zero/0.d+0/
-c/7
- parameter (half=0.5d+0, one=1.d+0, two=2.d+0, zero=0.d+0)
-c/
-c
-c/6
-c data bias/43/, dgnorm/1/, dstnrm/2/, dst0/3/, grdfac/45/,
-c 1 gthg/44/, gtstep/4/, nreduc/6/, nwtfac/46/, preduc/7/,
-c 2 radius/8/, stppar/5/
-c/7
- parameter (bias=43, dgnorm=1, dstnrm=2, dst0=3, grdfac=45,
- 1 gthg=44, gtstep=4, nreduc=6, nwtfac=46, preduc=7,
- 2 radius=8, stppar=5)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- nwtnrm = v(dst0)
- rlambd = one
- if (nwtnrm .gt. zero) rlambd = v(radius) / nwtnrm
- gnorm = v(dgnorm)
- ghinvg = two * v(nreduc)
- v(grdfac) = zero
- v(nwtfac) = zero
- if (rlambd .lt. one) go to 30
-c
-c *** the newton step is inside the trust region ***
-c
- v(stppar) = zero
- v(dstnrm) = nwtnrm
- v(gtstep) = -ghinvg
- v(preduc) = v(nreduc)
- v(nwtfac) = -one
- do 20 i = 1, n
- 20 step(i) = -nwtstp(i)
- go to 999
-c
- 30 v(dstnrm) = v(radius)
- cfact = (gnorm / v(gthg))**2
-c *** cauchy step = -cfact * g.
- cnorm = gnorm * cfact
- relax = one - v(bias) * (one - gnorm*cnorm/ghinvg)
- if (rlambd .lt. relax) go to 50
-c
-c *** step is between relaxed newton and full newton steps ***
-c
- v(stppar) = one - (rlambd - relax) / (one - relax)
- t = -rlambd
- v(gtstep) = t * ghinvg
- v(preduc) = rlambd * (one - half*rlambd) * ghinvg
- v(nwtfac) = t
- do 40 i = 1, n
- 40 step(i) = t * nwtstp(i)
- go to 999
-c
- 50 if (cnorm .lt. v(radius)) go to 70
-c
-c *** the cauchy step lies outside the trust region --
-c *** step = scaled cauchy step ***
-c
- t = -v(radius) / gnorm
- v(grdfac) = t
- v(stppar) = one + cnorm / v(radius)
- v(gtstep) = -v(radius) * gnorm
- v(preduc) = v(radius)*(gnorm - half*v(radius)*(v(gthg)/gnorm)**2)
- do 60 i = 1, n
- 60 step(i) = t * dig(i)
- go to 999
-c
-c *** compute dogleg step between cauchy and relaxed newton ***
-c *** femur = relaxed newton step minus cauchy step ***
-c
- 70 ctrnwt = cfact * relax * ghinvg / gnorm
-c *** ctrnwt = inner prod. of cauchy and relaxed newton steps,
-c *** scaled by gnorm**-1.
- t1 = ctrnwt - gnorm*cfact**2
-c *** t1 = inner prod. of femur and cauchy step, scaled by
-c *** gnorm**-1.
- t2 = v(radius)*(v(radius)/gnorm) - gnorm*cfact**2
- t = relax * nwtnrm
- femnsq = (t/gnorm)*t - ctrnwt - t1
-c *** femnsq = square of 2-norm of femur, scaled by gnorm**-1.
- t = t2 / (t1 + dsqrt(t1**2 + femnsq*t2))
-c *** dogleg step = cauchy step + t * femur.
- t1 = (t - one) * cfact
- v(grdfac) = t1
- t2 = -t * relax
- v(nwtfac) = t2
- v(stppar) = two - t
- v(gtstep) = t1*gnorm**2 + t2*ghinvg
- v(preduc) = -t1*gnorm * ((t2 + one)*gnorm)
- 1 - t2 * (one + half*t2)*ghinvg
- 2 - half * (v(gthg)*t1)**2
- do 80 i = 1, n
- 80 step(i) = t1*dig(i) + t2*nwtstp(i)
-c
- 999 return
-c *** last line of dbdog follows ***
- end
- subroutine ltvmul(n, x, l, y)
-c
-c *** compute x = (l**t)*y, where l is an n x n lower
-c *** triangular matrix stored compactly by rows. x and y may
-c *** occupy the same storage. ***
-c
- integer n
-cal double precision x(n), l(1), y(n)
- double precision x(n), l(n*(n+1)/2), y(n)
-c dimension l(n*(n+1)/2)
- integer i, ij, i0, j
- double precision yi, zero
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c
- i0 = 0
- do 20 i = 1, n
- yi = y(i)
- x(i) = zero
- do 10 j = 1, i
- ij = i0 + j
- x(j) = x(j) + yi*l(ij)
- 10 continue
- i0 = i0 + i
- 20 continue
- 999 return
-c *** last card of ltvmul follows ***
- end
- subroutine lupdat(beta, gamma, l, lambda, lplus, n, w, z)
-c
-c *** compute lplus = secant update of l ***
-c
-c *** parameter declarations ***
-c
- integer n
-cal double precision beta(n), gamma(n), l(1), lambda(n), lplus(1),
- double precision beta(n), gamma(n), l(n*(n+1)/2), lambda(n),
- 1 lplus(n*(n+1)/2),w(n), z(n)
-c dimension l(n*(n+1)/2), lplus(n*(n+1)/2)
-c
-c-------------------------- parameter usage --------------------------
-c
-c beta = scratch vector.
-c gamma = scratch vector.
-c l (input) lower triangular matrix, stored rowwise.
-c lambda = scratch vector.
-c lplus (output) lower triangular matrix, stored rowwise, which may
-c occupy the same storage as l.
-c n (input) length of vector parameters and order of matrices.
-c w (input, destroyed on output) right singular vector of rank 1
-c correction to l.
-c z (input, destroyed on output) left singular vector of rank 1
-c correction to l.
-c
-c------------------------------- notes -------------------------------
-c
-c *** application and usage restrictions ***
-c
-c this routine updates the cholesky factor l of a symmetric
-c positive definite matrix to which a secant update is being
-c applied -- it computes a cholesky factor lplus of
-c l * (i + z*w**t) * (i + w*z**t) * l**t. it is assumed that w
-c and z have been chosen so that the updated matrix is strictly
-c positive definite.
-c
-c *** algorithm notes ***
-c
-c this code uses recurrence 3 of ref. 1 (with d(j) = 1 for all j)
-c to compute lplus of the form l * (i + z*w**t) * q, where q
-c is an orthogonal matrix that makes the result lower triangular.
-c lplus may have some negative diagonal elements.
-c
-c *** references ***
-c
-c 1. goldfarb, d. (1976), factorized variable metric methods for uncon-
-c strained optimization, math. comput. 30, pp. 796-811.
-c
-c *** general ***
-c
-c coded by david m. gay (fall 1979).
-c this subroutine was written in connection with research supported
-c by the national science foundation under grants mcs-7600324 and
-c mcs-7906671.
-c
-c------------------------ external quantities ------------------------
-c
-c *** intrinsic functions ***
-c/+
- double precision dsqrt
-c/
-c-------------------------- local variables --------------------------
-c
- integer i, ij, j, jj, jp1, k, nm1, np1
- double precision a, b, bj, eta, gj, lj, lij, ljj, nu, s, theta,
- 1 wj, zj
- double precision one, zero
-c
-c *** data initializations ***
-c
-c/6
-c data one/1.d+0/, zero/0.d+0/
-c/7
- parameter (one=1.d+0, zero=0.d+0)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- nu = one
- eta = zero
- if (n .le. 1) go to 30
- nm1 = n - 1
-c
-c *** temporarily store s(j) = sum over k = j+1 to n of w(k)**2 in
-c *** lambda(j).
-c
- s = zero
- do 10 i = 1, nm1
- j = n - i
- s = s + w(j+1)**2
- lambda(j) = s
- 10 continue
-c
-c *** compute lambda, gamma, and beta by goldfarb*s recurrence 3.
-c
- do 20 j = 1, nm1
- wj = w(j)
- a = nu*z(j) - eta*wj
- theta = one + a*wj
- s = a*lambda(j)
- lj = dsqrt(theta**2 + a*s)
- if (theta .gt. zero) lj = -lj
- lambda(j) = lj
- b = theta*wj + s
- gamma(j) = b * nu / lj
- beta(j) = (a - b*eta) / lj
- nu = -nu / lj
- eta = -(eta + (a**2)/(theta - lj)) / lj
- 20 continue
- 30 lambda(n) = one + (nu*z(n) - eta*w(n))*w(n)
-c
-c *** update l, gradually overwriting w and z with l*w and l*z.
-c
- np1 = n + 1
- jj = n * (n + 1) / 2
- do 60 k = 1, n
- j = np1 - k
- lj = lambda(j)
- ljj = l(jj)
- lplus(jj) = lj * ljj
- wj = w(j)
- w(j) = ljj * wj
- zj = z(j)
- z(j) = ljj * zj
- if (k .eq. 1) go to 50
- bj = beta(j)
- gj = gamma(j)
- ij = jj + j
- jp1 = j + 1
- do 40 i = jp1, n
- lij = l(ij)
- lplus(ij) = lj*lij + bj*w(i) + gj*z(i)
- w(i) = w(i) + lij*wj
- z(i) = z(i) + lij*zj
- ij = ij + i
- 40 continue
- 50 jj = jj - j
- 60 continue
-c
- 999 return
-c *** last card of lupdat follows ***
- end
- subroutine lvmul(n, x, l, y)
-c
-c *** compute x = l*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)
-c dimension l(n*(n+1)/2)
- integer i, ii, ij, i0, j, np1
- double precision t, zero
-c/6
-c data zero/0.d+0/
-c/7
- parameter (zero=0.d+0)
-c/
-c
- np1 = n + 1
- i0 = n*(n+1)/2
- do 20 ii = 1, n
- i = np1 - ii
- i0 = i0 - i
- t = zero
- do 10 j = 1, i
- ij = i0 + j
- t = t + l(ij)*y(j)
- 10 continue
- x(i) = t
- 20 continue
- 999 return
-c *** last card of lvmul follows ***
- end
- subroutine vvmulp(n, x, y, z, k)
-c
-c *** set x(i) = y(i) * z(i)**k, 1 .le. i .le. n (for k = 1 or -1) ***
-c
- integer n, k
- double precision x(n), y(n), z(n)
- integer i
-c
- if (k .ge. 0) go to 20
- do 10 i = 1, n
- 10 x(i) = y(i) / z(i)
- go to 999
-c
- 20 do 30 i = 1, n
- 30 x(i) = y(i) * z(i)
- 999 return
-c *** last card of vvmulp follows ***
- end
- subroutine wzbfgs (l, n, s, w, y, z)
-c
-c *** compute y and z for lupdat corresponding to bfgs update.
-c
- integer n
-cal double precision l(1), s(n), w(n), y(n), z(n)
- double precision l(n*(n+1)/2), s(n), w(n), y(n), z(n)
-c dimension l(n*(n+1)/2)
-c
-c-------------------------- parameter usage --------------------------
-c
-c l (i/o) cholesky factor of hessian, a lower triang. matrix stored
-c compactly by rows.
-c n (input) order of l and length of s, w, y, z.
-c s (input) the step just taken.
-c w (output) right singular vector of rank 1 correction to l.
-c y (input) change in gradients corresponding to s.
-c z (output) left singular vector of rank 1 correction to l.
-c
-c------------------------------- notes -------------------------------
-c
-c *** algorithm notes ***
-c
-c when s is computed in certain ways, e.g. by gqtstp or
-c dbldog, it is possible to save n**2/2 operations since (l**t)*s
-c or l*(l**t)*s is then known.
-c if the bfgs update to l*(l**t) would reduce its determinant to
-c less than eps times its old value, then this routine in effect
-c replaces y by theta*y + (1 - theta)*l*(l**t)*s, where theta
-c (between 0 and 1) is chosen to make the reduction factor = eps.
-c
-c *** general ***
-c
-c coded by david m. gay (fall 1979).
-c this subroutine was written in connection with research supported
-c by the national science foundation under grants mcs-7600324 and
-c mcs-7906671.
-c
-c------------------------ external quantities ------------------------
-c
-c *** functions and subroutines called ***
-c
- external dotprd, livmul, ltvmul
- double precision dotprd
-c dotprd returns inner product of two vectors.
-c livmul multiplies l**-1 times a vector.
-c ltvmul multiplies l**t times a vector.
-c
-c *** intrinsic functions ***
-c/+
- double precision dsqrt
-c/
-c-------------------------- local variables --------------------------
-c
- integer i
- double precision cs, cy, eps, epsrt, one, shs, ys, theta
-c
-c *** data initializations ***
-c
-c/6
-c data eps/0.1d+0/, one/1.d+0/
-c/7
- parameter (eps=0.1d+0, one=1.d+0)
-c/
-c
-c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
-c
- call ltvmul(n, w, l, s)
- shs = dotprd(n, w, w)
- ys = dotprd(n, y, s)
- if (ys .ge. eps*shs) go to 10
- theta = (one - eps) * shs / (shs - ys)
- epsrt = dsqrt(eps)
- cy = theta / (shs * epsrt)
- cs = (one + (theta-one)/epsrt) / shs
- go to 20
- 10 cy = one / (dsqrt(ys) * dsqrt(shs))
- cs = one / shs
- 20 call livmul(n, z, l, y)
- do 30 i = 1, n
- 30 z(i) = cy * z(i) - cs * w(i)
-c
- 999 return
-c *** last card of wzbfgs follows ***
- end
projects / unres.git / blobdiff