1 subroutine sumsl_a(n, d, x, calcf, calcg, iv, liv, lv, v,
2 1 uiparm, urparm, ufparm)
4 c *** minimize general unconstrained objective function using ***
5 c *** analytic gradient and hessian approx. from secant update ***
8 integer iv(liv), uiparm(1)
9 double precision d(n), x(n), v(lv), urparm(1)
10 c dimension v(71 + n*(n+15)/2), uiparm(*), urparm(*)
11 external calcf, calcg, ufparm
15 c this routine interacts with subroutine sumit in an attempt
16 c to find an n-vector x* that minimizes the (unconstrained)
17 c objective function computed by calcf. (often the x* found is
18 c a local minimizer rather than a global one.)
20 c-------------------------- parameter usage --------------------------
22 c n........ (input) the number of variables on which f depends, i.e.,
23 c the number of components in x.
24 c d........ (input/output) a scale vector such that d(i)*x(i),
25 c i = 1,2,...,n, are all in comparable units.
26 c d can strongly affect the behavior of sumsl.
27 c finding the best choice of d is generally a trial-
28 c and-error process. choosing d so that d(i)*x(i)
29 c has about the same value for all i often works well.
30 c the defaults provided by subroutine deflt (see i
31 c below) require the caller to supply d.
32 c x........ (input/output) before (initially) calling sumsl, the call-
33 c er should set x to an initial guess at x*. when
34 c sumsl returns, x contains the best point so far
35 c found, i.e., the one that gives the least value so
37 c calcf.... (input) a subroutine that, given x, computes f(x). calcf
38 c must be declared external in the calling program.
40 c call calcf(n, x, nf, f, uiparm, urparm, ufparm)
41 c when calcf is called, nf is the invocation
42 c count for calcf. nf is included for possible use
43 c with calcg. if x is out of bounds (e.g., if it
44 c would cause overflow in computing f(x)), then calcf
45 c should set nf to 0. this will cause a shorter step
46 c to be attempted. (if x is in bounds, then calcf
47 c should not change nf.) the other parameters are as
48 c described above and below. calcf should not change
50 c calcg.... (input) a subroutine that, given x, computes g(x), the gra-
51 c dient of f at x. calcg must be declared external in
52 c the calling program. it is invoked by
53 c call calcg(n, x, nf, g, uiparm, urparm, ufaprm)
54 c when calcg is called, nf is the invocation
55 c count for calcf at the time f(x) was evaluated. the
56 c x passed to calcg is usually the one passed to calcf
57 c on either its most recent invocation or the one
58 c prior to it. if calcf saves intermediate results
59 c for use by calcg, then it is possible to tell from
60 c nf whether they are valid for the current x (or
61 c which copy is valid if two copies are kept). if g
62 c cannot be computed at x, then calcg should set nf to
63 c 0. in this case, sumsl will return with iv(1) = 65.
64 c (if g can be computed at x, then calcg should not
65 c changed nf.) the other parameters to calcg are as
66 c described above and below. calcg should not change
68 c iv....... (input/output) an integer value array of length liv (see
69 c below) that helps control the sumsl algorithm and
70 c that is used to store various intermediate quanti-
71 c ties. of particular interest are the initialization/
72 c return code iv(1) and the entries in iv that control
73 c printing and limit the number of iterations and func-
74 c tion evaluations. see the section on iv input
76 c liv...... (input) length of iv array. must be at least 60. if li
77 c is too small, then sumsl returns with iv(1) = 15.
78 c when sumsl returns, the smallest allowed value of
79 c liv is stored in iv(lastiv) -- see the section on
80 c iv output values below. (this is intended for use
81 c with extensions of sumsl that handle constraints.)
82 c lv....... (input) length of v array. must be at least 71+n*(n+15)/2.
83 c (at least 77+n*(n+17)/2 for smsno, at least
84 c 78+n*(n+12) for humsl). if lv is too small, then
85 c sumsl returns with iv(1) = 16. when sumsl returns,
86 c the smallest allowed value of lv is stored in
87 c iv(lastv) -- see the section on iv output values
89 c v........ (input/output) a floating-point value array of length l
90 c (see below) that helps control the sumsl algorithm
91 c and that is used to store various intermediate
92 c quantities. of particular interest are the entries
93 c in v that limit the length of the first step
94 c attempted (lmax0) and specify convergence tolerances
95 c (afctol, lmaxs, rfctol, sctol, xctol, xftol).
96 c uiparm... (input) user integer parameter array passed without change
98 c urparm... (input) user floating-point parameter array passed without
99 c change to calcf and calcg.
100 c ufparm... (input) user external subroutine or function passed without
101 c change to calcf and calcg.
103 c *** iv input values (from subroutine deflt) ***
105 c iv(1)... on input, iv(1) should have a value between 0 and 14......
106 c 0 and 12 mean this is a fresh start. 0 means that
107 c deflt(2, iv, liv, lv, v)
108 c is to be called to provide all default values to iv and
109 c v. 12 (the value that deflt assigns to iv(1)) means the
110 c caller has already called deflt and has possibly changed
111 c some iv and/or v entries to non-default values.
112 c 13 means deflt has been called and that sumsl (and
113 c sumit) should only do their storage allocation. that is,
114 c they should set the output components of iv that tell
115 c where various subarrays arrays of v begin, such as iv(g)
116 c (and, for humsl and humit only, iv(dtol)), and return.
117 c 14 means that a storage has been allocated (by a call
118 c with iv(1) = 13) and that the algorithm should be
119 c started. when called with iv(1) = 13, sumsl returns
120 c iv(1) = 14 unless liv or lv is too small (or n is not
121 c positive). default = 12.
122 c iv(inith).... iv(25) tells whether the hessian approximation h should
123 c be initialized. 1 (the default) means sumit should
124 c initialize h to the diagonal matrix whose i-th diagonal
125 c element is d(i)**2. 0 means the caller has supplied a
126 c cholesky factor l of the initial hessian approximation
127 c h = l*(l**t) in v, starting at v(iv(lmat)) = v(iv(42))
128 c (and stored compactly by rows). note that iv(lmat) may
129 c be initialized by calling sumsl with iv(1) = 13 (see
130 c the iv(1) discussion above). default = 1.
131 c iv(mxfcal)... iv(17) gives the maximum number of function evaluations
132 c (calls on calcf) allowed. if this number does not suf-
133 c fice, then sumsl returns with iv(1) = 9. default = 200.
134 c iv(mxiter)... iv(18) gives the maximum number of iterations allowed.
135 c it also indirectly limits the number of gradient evalua-
136 c tions (calls on calcg) to iv(mxiter) + 1. if iv(mxiter)
137 c iterations do not suffice, then sumsl returns with
138 c iv(1) = 10. default = 150.
139 c iv(outlev)... iv(19) controls the number and length of iteration sum-
140 c mary lines printed (by itsum). iv(outlev) = 0 means do
141 c not print any summary lines. otherwise, print a summary
142 c line after each abs(iv(outlev)) iterations. if iv(outlev)
143 c is positive, then summary lines of length 78 (plus carri-
144 c age control) are printed, including the following... the
145 c iteration and function evaluation counts, f = the current
146 c function value, relative difference in function values
147 c achieved by the latest step (i.e., reldf = (f0-v(f))/f01,
148 c where f01 is the maximum of abs(v(f)) and abs(v(f0)) and
149 c v(f0) is the function value from the previous itera-
150 c tion), the relative function reduction predicted for the
151 c step just taken (i.e., preldf = v(preduc) / f01, where
152 c v(preduc) is described below), the scaled relative change
153 c in x (see v(reldx) below), the step parameter for the
154 c step just taken (stppar = 0 means a full newton step,
155 c between 0 and 1 means a relaxed newton step, between 1
156 c and 2 means a double dogleg step, greater than 2 means
157 c a scaled down cauchy step -- see subroutine dbldog), the
158 c 2-norm of the scale vector d times the step just taken
159 c (see v(dstnrm) below), and npreldf, i.e.,
160 c v(nreduc)/f01, where v(nreduc) is described below -- if
161 c npreldf is positive, then it is the relative function
162 c reduction predicted for a newton step (one with
163 c stppar = 0). if npreldf is negative, then it is the
164 c negative of the relative function reduction predicted
165 c for a step computed with step bound v(lmaxs) for use in
166 c testing for singular convergence.
167 c if iv(outlev) is negative, then lines of length 50
168 c are printed, including only the first 6 items listed
169 c above (through reldx).
171 c iv(parprt)... iv(20) = 1 means print any nondefault v values on a
172 c fresh start or any changed v values on a restart.
173 c iv(parprt) = 0 means skip this printing. default = 1.
174 c iv(prunit)... iv(21) is the output unit number on which all printing
175 c is done. iv(prunit) = 0 means suppress all printing.
176 c default = standard output unit (unit 6 on most systems).
177 c iv(solprt)... iv(22) = 1 means print out the value of x returned (as
178 c well as the gradient and the scale vector d).
179 c iv(solprt) = 0 means skip this printing. default = 1.
180 c iv(statpr)... iv(23) = 1 means print summary statistics upon return-
181 c ing. these consist of the function value, the scaled
182 c relative change in x caused by the most recent step (see
183 c v(reldx) below), the number of function and gradient
184 c evaluations (calls on calcf and calcg), and the relative
185 c function reductions predicted for the last step taken and
186 c for a newton step (or perhaps a step bounded by v(lmaxs)
187 c -- see the descriptions of preldf and npreldf under
189 c iv(statpr) = 0 means skip this printing.
190 c iv(statpr) = -1 means skip this printing as well as that
191 c of the one-line termination reason message. default = 1.
192 c iv(x0prt).... iv(24) = 1 means print the initial x and scale vector d
193 c (on a fresh start only). iv(x0prt) = 0 means skip this
194 c printing. default = 1.
196 c *** (selected) iv output values ***
198 c iv(1)........ on output, iv(1) is a return code....
199 c 3 = x-convergence. the scaled relative difference (see
200 c v(reldx)) between the current parameter vector x and
201 c a locally optimal parameter vector is very likely at
203 c 4 = relative function convergence. the relative differ-
204 c ence between the current function value and its lo-
205 c cally optimal value is very likely at most v(rfctol).
206 c 5 = both x- and relative function convergence (i.e., the
207 c conditions for iv(1) = 3 and iv(1) = 4 both hold).
208 c 6 = absolute function convergence. the current function
209 c value is at most v(afctol) in absolute value.
210 c 7 = singular convergence. the hessian near the current
211 c iterate appears to be singular or nearly so, and a
212 c step of length at most v(lmaxs) is unlikely to yield
213 c a relative function decrease of more than v(sctol).
214 c 8 = false convergence. the iterates appear to be converg-
215 c ing to a noncritical point. this may mean that the
216 c convergence tolerances (v(afctol), v(rfctol),
217 c v(xctol)) are too small for the accuracy to which
218 c the function and gradient are being computed, that
219 c there is an error in computing the gradient, or that
220 c the function or gradient is discontinuous near x.
221 c 9 = function evaluation limit reached without other con-
222 c vergence (see iv(mxfcal)).
223 c 10 = iteration limit reached without other convergence
225 c 11 = stopx returned .true. (external interrupt). see the
227 c 14 = storage has been allocated (after a call with
229 c 17 = restart attempted with n changed.
230 c 18 = d has a negative component and iv(dtype) .le. 0.
231 c 19...43 = v(iv(1)) is out of range.
232 c 63 = f(x) cannot be computed at the initial x.
233 c 64 = bad parameters passed to assess (which should not
235 c 65 = the gradient could not be computed at x (see calcg
237 c 67 = bad first parameter to deflt.
238 c 80 = iv(1) was out of range.
239 c 81 = n is not positive.
240 c iv(g)........ iv(28) is the starting subscript in v of the current
241 c gradient vector (the one corresponding to x).
242 c iv(lastiv)... iv(44) is the least acceptable value of liv. (it is
243 c only set if liv is at least 44.)
244 c iv(lastv).... iv(45) is the least acceptable value of lv. (it is
245 c only set if liv is large enough, at least iv(lastiv).)
246 c iv(nfcall)... iv(6) is the number of calls so far made on calcf (i.e.,
247 c function evaluations).
248 c iv(ngcall)... iv(30) is the number of gradient evaluations (calls on
250 c iv(niter).... iv(31) is the number of iterations performed.
252 c *** (selected) v input values (from subroutine deflt) ***
254 c v(bias)..... v(43) is the bias parameter used in subroutine dbldog --
255 c see that subroutine for details. default = 0.8.
256 c v(afctol)... v(31) is the absolute function convergence tolerance.
257 c if sumsl finds a point where the function value is less
258 c than v(afctol) in absolute value, and if sumsl does not
259 c return with iv(1) = 3, 4, or 5, then it returns with
260 c iv(1) = 6. this test can be turned off by setting
261 c v(afctol) to zero. default = max(10**-20, machep**2),
262 c where machep is the unit roundoff.
263 c v(dinit).... v(38), if nonnegative, is the value to which the scale
264 c vector d is initialized. default = -1.
265 c v(lmax0).... v(35) gives the maximum 2-norm allowed for d times the
266 c very first step that sumsl attempts. this parameter can
267 c markedly affect the performance of sumsl.
268 c v(lmaxs).... v(36) is used in testing for singular convergence -- if
269 c the function reduction predicted for a step of length
270 c bounded by v(lmaxs) is at most v(sctol) * abs(f0), where
271 c f0 is the function value at the start of the current
272 c iteration, and if sumsl does not return with iv(1) = 3,
273 c 4, 5, or 6, then it returns with iv(1) = 7. default = 1.
274 c v(rfctol)... v(32) is the relative function convergence tolerance.
275 c if the current model predicts a maximum possible function
276 c reduction (see v(nreduc)) of at most v(rfctol)*abs(f0)
277 c at the start of the current iteration, where f0 is the
278 c then current function value, and if the last step attempt-
279 c ed achieved no more than twice the predicted function
280 c decrease, then sumsl returns with iv(1) = 4 (or 5).
281 c default = max(10**-10, machep**(2/3)), where machep is
283 c v(sctol).... v(37) is the singular convergence tolerance -- see the
284 c description of v(lmaxs) above.
285 c v(tuner1)... v(26) helps decide when to check for false convergence.
286 c this is done if the actual function decrease from the
287 c current step is no more than v(tuner1) times its predict-
288 c ed value. default = 0.1.
289 c v(xctol).... v(33) is the x-convergence tolerance. if a newton step
290 c (see v(nreduc)) is tried that has v(reldx) .le. v(xctol)
291 c and if this step yields at most twice the predicted func-
292 c tion decrease, then sumsl returns with iv(1) = 3 (or 5).
293 c (see the description of v(reldx) below.)
294 c default = machep**0.5, where machep is the unit roundoff.
295 c v(xftol).... v(34) is the false convergence tolerance. if a step is
296 c tried that gives no more than v(tuner1) times the predict-
297 c ed function decrease and that has v(reldx) .le. v(xftol),
298 c and if sumsl does not return with iv(1) = 3, 4, 5, 6, or
299 c 7, then it returns with iv(1) = 8. (see the description
300 c of v(reldx) below.) default = 100*machep, where
301 c machep is the unit roundoff.
302 c v(*)........ deflt supplies to v a number of tuning constants, with
303 c which it should ordinarily be unnecessary to tinker. see
304 c section 17 of version 2.2 of the nl2sol usage summary
305 c (i.e., the appendix to ref. 1) for details on v(i),
306 c i = decfac, incfac, phmnfc, phmxfc, rdfcmn, rdfcmx,
307 c tuner2, tuner3, tuner4, tuner5.
309 c *** (selected) v output values ***
311 c v(dgnorm)... v(1) is the 2-norm of (diag(d)**-1)*g, where g is the
312 c most recently computed gradient.
313 c v(dstnrm)... v(2) is the 2-norm of diag(d)*step, where step is the
315 c v(f)........ v(10) is the current function value.
316 c v(f0)....... v(13) is the function value at the start of the current
318 c v(nreduc)... v(6), if positive, is the maximum function reduction
319 c possible according to the current model, i.e., the func-
320 c tion reduction predicted for a newton step (i.e.,
321 c step = -h**-1 * g, where g is the current gradient and
322 c h is the current hessian approximation).
323 c if v(nreduc) is negative, then it is the negative of
324 c the function reduction predicted for a step computed with
325 c a step bound of v(lmaxs) for use in testing for singular
327 c v(preduc)... v(7) is the function reduction predicted (by the current
328 c quadratic model) for the current step. this (divided by
329 c v(f0)) is used in testing for relative function
331 c v(reldx).... v(17) is the scaled relative change in x caused by the
332 c current step, computed as
333 c max(abs(d(i)*(x(i)-x0(i)), 1 .le. i .le. p) /
334 c max(d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p),
335 c where x = x0 + step.
337 c------------------------------- notes -------------------------------
339 c *** algorithm notes ***
341 c this routine uses a hessian approximation computed from the
342 c bfgs update (see ref 3). only a cholesky factor of the hessian
343 c approximation is stored, and this is updated using ideas from
344 c ref. 4. steps are computed by the double dogleg scheme described
345 c in ref. 2. the steps are assessed as in ref. 1.
347 c *** usage notes ***
349 c after a return with iv(1) .le. 11, it is possible to restart,
350 c i.e., to change some of the iv and v input values described above
351 c and continue the algorithm from the point where it was interrupt-
352 c ed. iv(1) should not be changed, nor should any entries of i
353 c and v other than the input values (those supplied by deflt).
354 c those who do not wish to write a calcg which computes the
355 c gradient analytically should call smsno rather than sumsl.
356 c smsno uses finite differences to compute an approximate gradient.
357 c those who would prefer to provide f and g (the function and
358 c gradient) by reverse communication rather than by writing subrou-
359 c tines calcf and calcg may call on sumit directly. see the com-
360 c ments at the beginning of sumit.
361 c those who use sumsl interactively may wish to supply their
362 c own stopx function, which should return .true. if the break key
363 c has been pressed since stopx was last invoked. this makes it
364 c possible to externally interrupt sumsl (which will return with
365 c iv(1) = 11 if stopx returns .true.).
366 c storage for g is allocated at the end of v. thus the caller
367 c may make v longer than specified above and may allow calcg to use
368 c elements of g beyond the first n as scratch storage.
370 c *** portability notes ***
372 c the sumsl distribution tape contains both single- and double-
373 c precision versions of the sumsl source code, so it should be un-
374 c necessary to change precisions.
375 c only the functions imdcon and rmdcon contain machine-dependent
376 c constants. to change from one machine to another, it should
377 c suffice to change the (few) relevant lines in these functions.
378 c intrinsic functions are explicitly declared. on certain com-
379 c puters (e.g. univac), it may be necessary to comment out these
380 c declarations. so that this may be done automatically by a simple
381 c program, such declarations are preceded by a comment having c/+
382 c in columns 1-3 and blanks in columns 4-72 and are followed by
383 c a comment having c/ in columns 1 and 2 and blanks in columns 3-72.
384 c the sumsl source code is expressed in 1966 ansi standard
385 c fortran. it may be converted to fortran 77 by commenting out all
386 c lines that fall between a line having c/6 in columns 1-3 and a
387 c line having c/7 in columns 1-3 and by removing (i.e., replacing
388 c by a blank) the c in column 1 of the lines that follow the c/7
389 c line and precede a line having c/ in columns 1-2 and blanks in
390 c columns 3-72. these changes convert some data statements into
391 c parameter statements, convert some variables from real to
392 c character*4, and make the data statements that initialize these
393 c variables use character strings delimited by primes instead
394 c of hollerith constants. (such variables and data statements
395 c appear only in modules itsum and parck. parameter statements
396 c appear nearly everywhere.) these changes also add save state-
397 c ments for variables given machine-dependent constants by rmdcon.
401 c 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), algorithm 573 --
402 c an adaptive nonlinear least-squares algorithm, acm trans.
403 c math. software 7, pp. 369-383.
405 c 2. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
406 c mization algorithms which use function and gradient
407 c values, j. optim. theory applic. 28, pp. 453-482.
409 c 3. dennis, j.e., and more, j.j. (1977), quasi-newton methods, motiva-
410 c tion and theory, siam rev. 19, pp. 46-89.
412 c 4. goldfarb, d. (1976), factorized variable metric methods for uncon-
413 c strained optimization, math. comput. 30, pp. 796-811.
417 c coded by david m. gay (winter 1980). revised summer 1982.
418 c this subroutine was written in connection with research
419 c supported in part by the national science foundation under
420 c grants mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989,
424 c---------------------------- declarations ---------------------------
426 external deflt, sumit
428 c deflt... supplies default iv and v input components.
429 c sumit... reverse-communication routine that carries out sumsl algo-
435 c *** subscripts for iv ***
437 integer nextv, nfcall, nfgcal, g, toobig, vneed
440 c data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, toobig/2/, vneed/4/
442 parameter (nextv=47, nfcall=6, nfgcal=7, g=28, toobig=2, vneed=4)
445 c+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
447 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
449 if (iv1 .eq. 12 .or. iv1 .eq. 13) iv(vneed) = iv(vneed) + n
450 if (iv1 .eq. 14) go to 10
451 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
453 if (iv1 .eq. 12) iv(1) = 13
458 20 call sumit(d, f, v(g1), iv, liv, lv, n, v, x)
459 if (iv(1) - 2) 30, 40, 50
462 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
463 if (nf .le. 0) iv(toobig) = 1
466 40 call calcg(n, x, iv(nfgcal), v(g1), uiparm, urparm, ufparm)
469 50 if (iv(1) .ne. 14) go to 999
471 c *** storage allocation
474 iv(nextv) = iv(g) + n
475 if (iv1 .ne. 13) go to 10
478 c *** last card of sumsl follows ***
projects / unres.git / blob