2 !-----------------------------------------------------------------------------
12 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 !-----------------------------------------------------------------------------
19 !-----------------------------------------------------------------------------
20 subroutine assst(iv, liv, lv, v)
22 ! *** assess candidate step (***sol version 2.3) ***
30 ! this subroutine is called by an unconstrained minimization
31 ! routine to assess the next candidate step. it may recommend one
32 ! of several courses of action, such as accepting the step, recom-
33 ! puting it using the same or a new quadratic model, or halting due
34 ! to convergence or false convergence. see the return code listing
37 !-------------------------- parameter usage --------------------------
39 ! iv (i/o) integer parameter and scratch vector -- see description
40 ! below of iv values referenced.
41 ! liv (in) length of iv array.
42 ! lv (in) length of v array.
43 ! v (i/o) real parameter and scratch vector -- see description
44 ! below of v values referenced.
46 ! *** iv values referenced ***
48 ! iv(irc) (i/o) on input for the first step tried in a new iteration,
49 ! iv(irc) should be set to 3 or 4 (the value to which it is
50 ! set when step is definitely to be accepted). on input
51 ! after step has been recomputed, iv(irc) should be
52 ! unchanged since the previous return of assst.
53 ! on output, iv(irc) is a return code having one of the
55 ! 1 = switch models or try smaller step.
56 ! 2 = switch models or accept step.
57 ! 3 = accept step and determine v(radfac) by gradient
59 ! 4 = accept step, v(radfac) has been determined.
60 ! 5 = recompute step (using the same model).
61 ! 6 = recompute step with radius = v(lmaxs) but do not
62 ! evaulate the objective function.
63 ! 7 = x-convergence (see v(xctol)).
64 ! 8 = relative function convergence (see v(rfctol)).
65 ! 9 = both x- and relative function convergence.
66 ! 10 = absolute function convergence (see v(afctol)).
67 ! 11 = singular convergence (see v(lmaxs)).
68 ! 12 = false convergence (see v(xftol)).
69 ! 13 = iv(irc) was out of range on input.
70 ! return code i has precdence over i+1 for i = 9, 10, 11.
71 ! iv(mlstgd) (i/o) saved value of iv(model).
72 ! iv(model) (i/o) on input, iv(model) should be an integer identifying
73 ! the current quadratic model of the objective function.
74 ! if a previous step yielded a better function reduction,
75 ! then iv(model) will be set to iv(mlstgd) on output.
76 ! iv(nfcall) (in) invocation count for the objective function.
77 ! iv(nfgcal) (i/o) value of iv(nfcall) at step that gave the biggest
78 ! function reduction this iteration. iv(nfgcal) remains
79 ! unchanged until a function reduction is obtained.
80 ! iv(radinc) (i/o) the number of radius increases (or minus the number
81 ! of decreases) so far this iteration.
82 ! iv(restor) (out) set to 1 if v(f) has been restored and x should be
83 ! restored to its initial value, to 2 if x should be saved,
84 ! to 3 if x should be restored from the saved value, and to
86 ! iv(stage) (i/o) count of the number of models tried so far in the
88 ! iv(stglim) (in) maximum number of models to consider.
89 ! iv(switch) (out) set to 0 unless a new model is being tried and it
90 ! gives a smaller function value than the previous model,
91 ! in which case assst sets iv(switch) = 1.
92 ! iv(toobig) (in) is nonzero if step was too big (e.g. if it caused
94 ! iv(xirc) (i/o) value that iv(irc) would have in the absence of
95 ! convergence, false convergence, and oversized steps.
97 ! *** v values referenced ***
99 ! v(afctol) (in) absolute function convergence tolerance. if the
100 ! absolute value of the current function value v(f) is less
101 ! than v(afctol), then assst returns with iv(irc) = 10.
102 ! v(decfac) (in) factor by which to decrease radius when iv(toobig) is
104 ! v(dstnrm) (in) the 2-norm of d*step.
105 ! v(dstsav) (i/o) value of v(dstnrm) on saved step.
106 ! v(dst0) (in) the 2-norm of d times the newton step (when defined,
107 ! i.e., for v(nreduc) .ge. 0).
108 ! v(f) (i/o) on both input and output, v(f) is the objective func-
109 ! tion value at x. if x is restored to a previous value,
110 ! then v(f) is restored to the corresponding value.
111 ! v(fdif) (out) the function reduction v(f0) - v(f) (for the output
112 ! value of v(f) if an earlier step gave a bigger function
113 ! decrease, and for the input value of v(f) otherwise).
114 ! v(flstgd) (i/o) saved value of v(f).
115 ! v(f0) (in) objective function value at start of iteration.
116 ! v(gtslst) (i/o) value of v(gtstep) on saved step.
117 ! v(gtstep) (in) inner product between step and gradient.
118 ! v(incfac) (in) minimum factor by which to increase radius.
119 ! v(lmaxs) (in) maximum reasonable step size (and initial step bound).
120 ! if the actual function decrease is no more than twice
121 ! what was predicted, if a return with iv(irc) = 7, 8, 9,
122 ! or 10 does not occur, if v(dstnrm) .gt. v(lmaxs), and if
123 ! v(preduc) .le. v(sctol) * abs(v(f0)), then assst re-
124 ! turns with iv(irc) = 11. if so doing appears worthwhile,
125 ! then assst repeats this test with v(preduc) computed for
126 ! a step of length v(lmaxs) (by a return with iv(irc) = 6).
127 ! v(nreduc) (i/o) function reduction predicted by quadratic model for
128 ! newton step. if assst is called with iv(irc) = 6, i.e.,
129 ! if v(preduc) has been computed with radius = v(lmaxs) for
130 ! use in the singular convervence test, then v(nreduc) is
131 ! set to -v(preduc) before the latter is restored.
132 ! v(plstgd) (i/o) value of v(preduc) on saved step.
133 ! v(preduc) (i/o) function reduction predicted by quadratic model for
135 ! v(radfac) (out) factor to be used in determining the new radius,
136 ! which should be v(radfac)*dst, where dst is either the
137 ! output value of v(dstnrm) or the 2-norm of
138 ! diag(newd)*step for the output value of step and the
139 ! updated version, newd, of the scale vector d. for
140 ! iv(irc) = 3, v(radfac) = 1.0 is returned.
141 ! v(rdfcmn) (in) minimum value for v(radfac) in terms of the input
142 ! value of v(dstnrm) -- suggested value = 0.1.
143 ! v(rdfcmx) (in) maximum value for v(radfac) -- suggested value = 4.0.
144 ! v(reldx) (in) scaled relative change in x caused by step, computed
145 ! (e.g.) by function reldst as
146 ! max (d(i)*abs(x(i)-x0(i)), 1 .le. i .le. p) /
147 ! max (d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p).
148 ! v(rfctol) (in) relative function convergence tolerance. if the
149 ! actual function reduction is at most twice what was pre-
150 ! dicted and v(nreduc) .le. v(rfctol)*abs(v(f0)), then
151 ! assst returns with iv(irc) = 8 or 9.
152 ! v(stppar) (in) marquardt parameter -- 0 means full newton step.
153 ! v(tuner1) (in) tuning constant used to decide if the function
154 ! reduction was much less than expected. suggested
156 ! v(tuner2) (in) tuning constant used to decide if the function
157 ! reduction was large enough to accept step. suggested
159 ! v(tuner3) (in) tuning constant used to decide if the radius
160 ! should be increased. suggested value = 0.75.
161 ! v(xctol) (in) x-convergence criterion. if step is a newton step
162 ! (v(stppar) = 0) having v(reldx) .le. v(xctol) and giving
163 ! at most twice the predicted function decrease, then
164 ! assst returns iv(irc) = 7 or 9.
165 ! v(xftol) (in) false convergence tolerance. if step gave no or only
166 ! a small function decrease and v(reldx) .le. v(xftol),
167 ! then assst returns with iv(irc) = 12.
169 !------------------------------- notes -------------------------------
171 ! *** application and usage restrictions ***
173 ! this routine is called as part of the nl2sol (nonlinear
174 ! least-squares) package. it may be used in any unconstrained
175 ! minimization solver that uses dogleg, goldfeld-quandt-trotter,
176 ! or levenberg-marquardt steps.
178 ! *** algorithm notes ***
180 ! see (1) for further discussion of the assessing and model
181 ! switching strategies. while nl2sol considers only two models,
182 ! assst is designed to handle any number of models.
184 ! *** usage notes ***
186 ! on the first call of an iteration, only the i/o variables
187 ! step, x, iv(irc), iv(model), v(f), v(dstnrm), v(gtstep), and
188 ! v(preduc) need have been initialized. between calls, no i/o
189 ! values execpt step, x, iv(model), v(f) and the stopping toler-
190 ! ances should be changed.
191 ! after a return for convergence or false convergence, one can
192 ! change the stopping tolerances and call assst again, in which
193 ! case the stopping tests will be repeated.
197 ! (1) dennis, j.e., jr., gay, d.m., and welsch, r.e. (1981),
198 ! an adaptive nonlinear least-squares algorithm,
199 ! acm trans. math. software, vol. 7, no. 3.
201 ! (2) powell, m.j.d. (1970) a fortran subroutine for solving
202 ! systems of nonlinear algebraic equations, in numerical
203 ! methods for nonlinear algebraic equations, edited by
204 ! p. rabinowitz, gordon and breach, london.
208 ! john dennis designed much of this routine, starting with
209 ! ideas in (2). roy welsch suggested the model switching strategy.
210 ! david gay and stephen peters cast this subroutine into a more
211 ! portable form (winter 1977), and david gay cast it into its
212 ! present form (fall 1978).
216 ! this subroutine was written in connection with research
217 ! supported by the national science foundation under grants
218 ! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
221 !------------------------ external quantities ------------------------
223 ! *** no external functions and subroutines ***
225 ! *** intrinsic functions ***
227 !el real(kind=8) :: dabs, dmax1
229 ! *** no common blocks ***
231 !-------------------------- local variables --------------------------
235 real(kind=8) :: emax, emaxs, gts, rfac1, xmax
236 !el real(kind=8) :: half, one, onep2, two, zero
238 ! *** subscripts for iv and v ***
240 !el integer :: afctol, decfac, dstnrm, dstsav, dst0, f, fdif, flstgd, f0,&
241 !el gtslst, gtstep, incfac, irc, lmaxs, mlstgd, model, nfcall,&
242 !el nfgcal, nreduc, plstgd, preduc, radfac, radinc, rdfcmn,&
243 !el rdfcmx, reldx, restor, rfctol, sctol, stage, stglim,&
244 !el stppar, switch, toobig, tuner1, tuner2, tuner3, xctol,&
248 ! *** data initializations ***
251 ! data half/0.5d+0/, one/1.d+0/, onep2/1.2d+0/, two/2.d+0/,
254 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, onep2=1.2d+0, two=2.d+0,&
259 ! data irc/29/, mlstgd/32/, model/5/, nfcall/6/, nfgcal/7/,
260 ! 1 radinc/8/, restor/9/, stage/10/, stglim/11/, switch/12/,
261 ! 2 toobig/2/, xirc/13/
263 integer,parameter :: irc=29, mlstgd=32, model=5, nfcall=6, nfgcal=7,&
264 radinc=8, restor=9, stage=10, stglim=11, switch=12,&
268 ! data afctol/31/, decfac/22/, dstnrm/2/, dst0/3/, dstsav/18/,
269 ! 1 f/10/, fdif/11/, flstgd/12/, f0/13/, gtslst/14/, gtstep/4/,
270 ! 2 incfac/23/, lmaxs/36/, nreduc/6/, plstgd/15/, preduc/7/,
271 ! 3 radfac/16/, rdfcmn/24/, rdfcmx/25/, reldx/17/, rfctol/32/,
272 ! 4 sctol/37/, stppar/5/, tuner1/26/, tuner2/27/, tuner3/28/,
273 ! 5 xctol/33/, xftol/34/
275 integer,parameter :: afctol=31, decfac=22, dstnrm=2, dst0=3, dstsav=18,&
276 f=10, fdif=11, flstgd=12, f0=13, gtslst=14, gtstep=4,&
277 incfac=23, lmaxs=36, nreduc=6, plstgd=15, preduc=7,&
278 radfac=16, rdfcmn=24, rdfcmx=25, reldx=17, rfctol=32,&
279 sctol=37, stppar=5, tuner1=26, tuner2=27, tuner3=28,&
283 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
291 if (i .ge. 1 .and. i .le. 12) &
292 go to (20,30,10,10,40,280,220,220,220,220,220,170), i
296 ! *** initialize for new iteration ***
301 if (iv(toobig) .eq. 0) go to 110
306 ! *** step was recomputed with new model or smaller radius ***
307 ! *** first decide which ***
309 20 if (iv(model) .ne. iv(mlstgd)) go to 30
310 ! *** old model retained, smaller radius tried ***
311 ! *** do not consider any more new models this iteration ***
312 iv(stage) = iv(stglim)
316 ! *** a new model is being tried. decide whether to keep it. ***
318 30 iv(stage) = iv(stage) + 1
320 ! *** now we add the possibiltiy that step was recomputed with ***
321 ! *** the same model, perhaps because of an oversized step. ***
323 40 if (iv(stage) .gt. 0) go to 50
325 ! *** step was recomputed because it was too big. ***
327 if (iv(toobig) .ne. 0) go to 60
329 ! *** restore iv(stage) and pick up where we left off. ***
331 iv(stage) = -iv(stage)
333 go to (20, 30, 110, 110, 70), i
335 50 if (iv(toobig) .eq. 0) go to 70
337 ! *** handle oversize step ***
339 if (iv(radinc) .gt. 0) go to 80
340 iv(stage) = -iv(stage)
343 60 v(radfac) = v(decfac)
344 iv(radinc) = iv(radinc) - 1
349 70 if (v(f) .lt. v(flstgd)) go to 110
351 ! *** the new step is a loser. restore old model. ***
353 if (iv(model) .eq. iv(mlstgd)) go to 80
354 iv(model) = iv(mlstgd)
357 ! *** restore step, etc. only if a previous step decreased v(f).
359 80 if (v(flstgd) .ge. v(f0)) go to 110
362 v(preduc) = v(plstgd)
363 v(gtstep) = v(gtslst)
364 if (iv(switch) .eq. 0) rfac1 = v(dstnrm) / v(dstsav)
365 v(dstnrm) = v(dstsav)
369 110 v(fdif) = v(f0) - v(f)
370 if (v(fdif) .gt. v(tuner2) * v(preduc)) go to 140
371 if(iv(radinc).gt.0) go to 140
373 ! *** no (or only a trivial) function decrease
374 ! *** -- so try new model or smaller radius
376 if (v(f) .lt. v(f0)) go to 120
377 iv(mlstgd) = iv(model)
384 if (iv(stage) .lt. iv(stglim)) go to 160
386 iv(radinc) = iv(radinc) - 1
389 ! *** nontrivial function decrease achieved ***
393 v(dstsav) = v(dstnrm)
394 if (v(fdif) .gt. v(preduc)*v(tuner1)) go to 190
396 ! *** decrease was much less than predicted -- either change models
397 ! *** or accept step with decreased radius.
399 if (iv(stage) .ge. iv(stglim)) go to 150
400 ! *** consider switching models ***
404 ! *** accept step with decreased radius ***
408 ! *** set v(radfac) to fletcher*s decrease factor ***
410 160 iv(xirc) = iv(irc)
411 emax = v(gtstep) + v(fdif)
412 v(radfac) = half * rfac1
413 if (emax .lt. v(gtstep)) v(radfac) = rfac1 * dmax1(v(rdfcmn),&
414 half * v(gtstep)/emax)
416 ! *** do false convergence test ***
418 170 if (v(reldx) .le. v(xftol)) go to 180
420 if (v(f) .lt. v(f0)) go to 200
426 ! *** handle good function decrease ***
428 190 if (v(fdif) .lt. (-v(tuner3) * v(gtstep))) go to 210
430 ! *** increasing radius looks worthwhile. see if we just
431 ! *** recomputed step with a decreased radius or restored step
432 ! *** after recomputing it with a larger radius.
434 if (iv(radinc) .lt. 0) go to 210
435 if (iv(restor) .eq. 1) go to 210
437 ! *** we did not. try a longer step unless this was a newton
440 v(radfac) = v(rdfcmx)
442 if (v(fdif) .lt. (half/v(radfac) - one) * gts) &
443 v(radfac) = dmax1(v(incfac), half*gts/(gts + v(fdif)))
445 if (v(stppar) .eq. zero) go to 230
446 if (v(dst0) .ge. zero .and. (v(dst0) .lt. two*v(dstnrm) &
447 .or. v(nreduc) .lt. onep2*v(fdif))) go to 230
448 ! *** step was not a newton step. recompute it with
449 ! *** a larger radius.
451 iv(radinc) = iv(radinc) + 1
453 ! *** save values corresponding to good step ***
456 iv(mlstgd) = iv(model)
457 if (iv(restor) .ne. 1) iv(restor) = 2
458 v(dstsav) = v(dstnrm)
460 v(plstgd) = v(preduc)
461 v(gtslst) = v(gtstep)
464 ! *** accept step with radius unchanged ***
470 ! *** come here for a restart after convergence ***
472 220 iv(irc) = iv(xirc)
473 if (v(dstsav) .ge. zero) go to 240
477 ! *** perform convergence tests ***
479 230 iv(xirc) = iv(irc)
480 240 if (iv(restor) .eq. 1 .and. v(flstgd) .lt. v(f0)) iv(restor) = 3
481 if (half * v(fdif) .gt. v(preduc)) go to 999
482 emax = v(rfctol) * dabs(v(f0))
483 emaxs = v(sctol) * dabs(v(f0))
484 if (v(dstnrm) .gt. v(lmaxs) .and. v(preduc) .le. emaxs) &
486 if (v(dst0) .lt. zero) go to 250
488 if ((v(nreduc) .gt. zero .and. v(nreduc) .le. emax) .or. &
489 (v(nreduc) .eq. zero .and. v(preduc) .eq. zero)) i = 2
490 if (v(stppar) .eq. zero .and. v(reldx) .le. v(xctol) &
491 .and. goodx) i = i + 1
492 if (i .gt. 0) iv(irc) = i + 6
494 ! *** consider recomputing step of length v(lmaxs) for singular
495 ! *** convergence test.
497 250 if (iv(irc) .gt. 5 .and. iv(irc) .ne. 12) go to 999
498 if (v(dstnrm) .gt. v(lmaxs)) go to 260
499 if (v(preduc) .ge. emaxs) go to 999
500 if (v(dst0) .le. zero) go to 270
501 if (half * v(dst0) .le. v(lmaxs)) go to 999
503 260 if (half * v(dstnrm) .le. v(lmaxs)) go to 999
504 xmax = v(lmaxs) / v(dstnrm)
505 if (xmax * (two - xmax) * v(preduc) .ge. emaxs) go to 999
506 270 if (v(nreduc) .lt. zero) go to 290
508 ! *** recompute v(preduc) for use in singular convergence test ***
510 v(gtslst) = v(gtstep)
511 v(dstsav) = v(dstnrm)
512 if (iv(irc) .eq. 12) v(dstsav) = -v(dstsav)
513 v(plstgd) = v(preduc)
516 if (i .eq. 3) iv(restor) = 0
520 ! *** perform singular convergence test with recomputed v(preduc) ***
522 280 v(gtstep) = v(gtslst)
523 v(dstnrm) = dabs(v(dstsav))
525 if (v(dstsav) .le. zero) iv(irc) = 12
526 v(nreduc) = -v(preduc)
527 v(preduc) = v(plstgd)
529 290 if (-v(nreduc) .le. v(sctol) * dabs(v(f0))) iv(irc) = 11
533 ! *** last card of assst follows ***
535 !-----------------------------------------------------------------------------
536 subroutine deflt(alg, iv, liv, lv, v)
538 ! *** supply ***sol (version 2.3) default values to iv and v ***
540 ! *** alg = 1 means regression constants.
541 ! *** alg = 2 means general unconstrained optimization constants.
544 integer :: alg, iv(liv)
545 real(kind=8) :: v(lv)
547 !el external imdcon, vdflt
549 ! imdcon... returns machine-dependent integer constants.
550 ! vdflt.... provides default values to v.
553 integer :: miniv(2), minv(2)
555 ! *** subscripts for iv ***
557 !el integer algsav, covprt, covreq, dtype, hc, ierr, inith, inits,
558 !el 1 ipivot, ivneed, lastiv, lastv, lmat, mxfcal, mxiter,
559 !el 2 nfcov, ngcov, nvdflt, outlev, parprt, parsav, perm,
560 !el 3 prunit, qrtyp, rdreq, rmat, solprt, statpr, vneed,
563 ! *** iv subscript values ***
566 ! data algsav/51/, covprt/14/, covreq/15/, dtype/16/, hc/71/,
567 ! 1 ierr/75/, inith/25/, inits/25/, ipivot/76/, ivneed/3/,
568 ! 2 lastiv/44/, lastv/45/, lmat/42/, mxfcal/17/, mxiter/18/,
569 ! 3 nfcov/52/, ngcov/53/, nvdflt/50/, outlev/19/, parprt/20/,
570 ! 4 parsav/49/, perm/58/, prunit/21/, qrtyp/80/, rdreq/57/,
571 ! 5 rmat/78/, solprt/22/, statpr/23/, vneed/4/, vsave/60/,
574 integer,parameter :: algsav=51, covprt=14, covreq=15, dtype=16, hc=71,&
575 ierr=75, inith=25, inits=25, ipivot=76, ivneed=3,&
576 lastiv=44, lastv=45, lmat=42, mxfcal=17, mxiter=18,&
577 nfcov=52, ngcov=53, nvdflt=50, outlev=19, parprt=20,&
578 parsav=49, perm=58, prunit=21, qrtyp=80, rdreq=57,&
579 rmat=78, solprt=22, statpr=23, vneed=4, vsave=60,&
582 data miniv(1)/80/, miniv(2)/59/, minv(1)/98/, minv(2)/71/
586 !------------------------------- body --------------------------------
588 if (alg .lt. 1 .or. alg .gt. 2) go to 40
590 if (liv .lt. miv) go to 20
592 if (lv .lt. mv) go to 30
593 call vdflt(alg, lv, v)
605 iv(prunit) = imdcon(1)
611 if (alg .ge. 2) go to 10
613 ! *** regression values
630 ! *** general optimization values
649 ! *** last card of deflt follows ***
651 !-----------------------------------------------------------------------------
652 real(kind=8) function dotprd(p,x,y)
654 ! *** return the inner product of the p-vectors x and y. ***
657 real(kind=8) :: x(p), y(p)
660 !el real(kind=8) :: one, zero
661 real(kind=8) :: sqteta, t
663 !el real(kind=8) :: dmax1, dabs
666 !el real(kind=8) :: rmdcon
668 ! *** rmdcon(2) returns a machine-dependent constant, sqteta, which
669 ! *** is slightly larger than the smallest positive number that
670 ! *** can be squared without underflowing.
673 ! data one/1.d+0/, sqteta/0.d+0/, zero/0.d+0/
675 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
680 if (p .le. 0) go to 999
681 !rc if (sqteta .eq. zero) sqteta = rmdcon(2)
683 !rc t = dmax1(dabs(x(i)), dabs(y(i)))
684 !rc if (t .gt. one) go to 10
685 !rc if (t .lt. sqteta) go to 20
686 !rc t = (x(i)/sqteta)*y(i)
687 !rc if (dabs(t) .lt. sqteta) go to 20
688 10 dotprd = dotprd + x(i)*y(i)
692 ! *** last card of dotprd follows ***
694 !-----------------------------------------------------------------------------
695 subroutine itsum(d, g, iv, liv, lv, p, v, x)
697 ! *** print iteration summary for ***sol (version 2.3) ***
699 ! *** parameter declarations ***
701 integer :: liv, lv, p
703 real(kind=8) :: d(p), g(p), v(lv), x(p)
705 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
707 ! *** local variables ***
709 integer :: alg, i, iv1, m, nf, ng, ol, pu
711 ! real model1(6), model2(6)
713 character(len=4) :: model1(6), model2(6)
715 real(kind=8) :: nreldf, oldf, preldf, reldf !el, zero
717 ! *** intrinsic functions ***
720 !el real(kind=8) :: dabs, dmax1
722 ! *** no external functions or subroutines ***
724 ! *** subscripts for iv and v ***
726 !el integer algsav, dstnrm, f, fdif, f0, needhd, nfcall, nfcov, ngcov,
727 !el 1 ngcall, niter, nreduc, outlev, preduc, prntit, prunit,
728 !el 2 reldx, solprt, statpr, stppar, sused, x0prt
730 ! *** iv subscript values ***
733 ! data algsav/51/, needhd/36/, nfcall/6/, nfcov/52/, ngcall/30/,
734 ! 1 ngcov/53/, niter/31/, outlev/19/, prntit/39/, prunit/21/,
735 ! 2 solprt/22/, statpr/23/, sused/64/, x0prt/24/
737 integer,parameter :: algsav=51, needhd=36, nfcall=6, nfcov=52, ngcall=30,&
738 ngcov=53, niter=31, outlev=19, prntit=39, prunit=21,&
739 solprt=22, statpr=23, sused=64, x0prt=24
742 ! *** v subscript values ***
745 ! data dstnrm/2/, f/10/, f0/13/, fdif/11/, nreduc/6/, preduc/7/,
746 ! 1 reldx/17/, stppar/5/
748 integer,parameter :: dstnrm=2, f=10, f0=13, fdif=11, nreduc=6, preduc=7,&
755 real(kind=8),parameter :: zero=0.d+0
758 ! data model1(1)/4h /, model1(2)/4h /, model1(3)/4h /,
759 ! 1 model1(4)/4h /, model1(5)/4h g /, model1(6)/4h s /,
760 ! 2 model2(1)/4h g /, model2(2)/4h s /, model2(3)/4hg-s /,
761 ! 3 model2(4)/4hs-g /, model2(5)/4h-s-g/, model2(6)/4h-g-s/
763 data model1/' ',' ',' ',' ',' g ',' s '/,&
764 model2/' g ',' s ','g-s ','s-g ','-s-g','-g-s'/
767 !------------------------------- body --------------------------------
770 if (pu .eq. 0) go to 999
772 if (iv1 .gt. 62) iv1 = iv1 - 51
775 if (iv1 .lt. 2 .or. iv1 .gt. 15) go to 370
776 if (iv1 .ge. 12) go to 120
777 if (iv1 .eq. 2 .and. iv(niter) .eq. 0) go to 390
778 if (ol .eq. 0) go to 120
779 if (iv1 .ge. 10 .and. iv(prntit) .eq. 0) go to 120
780 if (iv1 .gt. 2) go to 10
781 iv(prntit) = iv(prntit) + 1
782 if (iv(prntit) .lt. iabs(ol)) go to 999
783 10 nf = iv(nfcall) - iabs(iv(nfcov))
787 oldf = dmax1(dabs(v(f0)), dabs(v(f)))
788 if (oldf .le. zero) go to 20
789 reldf = v(fdif) / oldf
790 preldf = v(preduc) / oldf
791 20 if (ol .gt. 0) go to 60
793 ! *** print short summary line ***
795 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,30)
796 30 format(/10h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,&
798 if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,40)
799 40 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,&
802 if (alg .eq. 2) go to 50
804 write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
805 model1(m), model2(m), v(stppar)
808 50 write(pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
812 ! *** print long summary line ***
814 60 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,70)
815 70 format(/11h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,&
816 2x,13hmodel stppar,2x,6hd*step,2x,7hnpreldf)
817 if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,80)
818 80 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,&
819 3x,6hstppar,3x,6hd*step,3x,7hnpreldf)
822 if (oldf .gt. zero) nreldf = v(nreduc) / oldf
823 if (alg .eq. 2) go to 90
825 write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
826 model1(m), model2(m), v(stppar), v(dstnrm), nreldf
829 90 write(pu,110) iv(niter), nf, v(f), reldf, preldf,&
830 v(reldx), v(stppar), v(dstnrm), nreldf
831 100 format(i6,i5,d10.3,2d9.2,d8.1,a3,a4,2d8.1,d9.2)
832 110 format(i6,i5,d11.3,2d10.2,3d9.1,d10.2)
834 120 if (iv(statpr) .lt. 0) go to 430
835 go to (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,&
839 140 format(/26h ***** x-convergence *****)
843 160 format(/42h ***** relative function convergence *****)
847 180 format(/49h ***** x- and relative function convergence *****)
851 200 format(/42h ***** absolute function convergence *****)
855 220 format(/33h ***** singular convergence *****)
859 240 format(/30h ***** false convergence *****)
863 260 format(/38h ***** function evaluation limit *****)
867 280 format(/28h ***** iteration limit *****)
871 300 format(/18h ***** stopx *****)
875 320 format(/44h ***** initial f(x) cannot be computed *****)
880 340 format(/37h ***** bad parameters to assess *****)
884 360 format(/43h ***** gradient could not be computed *****)
885 if (iv(niter) .gt. 0) go to 480
888 370 write(pu,380) iv(1)
889 380 format(/14h ***** iv(1) =,i5,6h *****)
892 ! *** initial call on itsum ***
894 390 if (iv(x0prt) .ne. 0) write(pu,400) (i, x(i), d(i), i = 1, p)
895 400 format(/23h i initial x(i),8x,4hd(i)//(1x,i5,d17.6,d14.3))
896 ! *** the following are to avoid undefined variables when the
897 ! *** function evaluation limit is 1...
903 if (iv1 .ge. 12) go to 999
906 if (ol .eq. 0) go to 999
907 if (ol .lt. 0 .and. alg .eq. 1) write(pu,30)
908 if (ol .lt. 0 .and. alg .eq. 2) write(pu,40)
909 if (ol .gt. 0 .and. alg .eq. 1) write(pu,70)
910 if (ol .gt. 0 .and. alg .eq. 2) write(pu,80)
911 if (alg .eq. 1) write(pu,410) v(f)
912 if (alg .eq. 2) write(pu,420) v(f)
913 410 format(/11h 0 1,d10.3)
914 !365 format(/11h 0 1,e11.3)
915 420 format(/11h 0 1,d11.3)
918 ! *** print various information requested on solution ***
921 if (iv(statpr) .eq. 0) go to 480
922 oldf = dmax1(dabs(v(f0)), dabs(v(f)))
925 if (oldf .le. zero) go to 440
926 preldf = v(preduc) / oldf
927 nreldf = v(nreduc) / oldf
928 440 nf = iv(nfcall) - iv(nfcov)
929 ng = iv(ngcall) - iv(ngcov)
930 write(pu,450) v(f), v(reldx), nf, ng, preldf, nreldf
931 450 format(/9h function,d17.6,8h reldx,d17.3/12h func. evals,&
932 i8,9x,11hgrad. evals,i8/7h preldf,d16.3,6x,7hnpreldf,d15.3)
934 if (iv(nfcov) .gt. 0) write(pu,460) iv(nfcov)
935 460 format(/1x,i4,50h extra func. evals for covariance and diagnostics.)
936 if (iv(ngcov) .gt. 0) write(pu,470) iv(ngcov)
937 470 format(1x,i4,50h extra grad. evals for covariance and diagnostics.)
939 480 if (iv(solprt) .eq. 0) go to 999
942 490 format(/22h i final x(i),8x,4hd(i),10x,4hg(i)/)
944 write(pu,510) i, x(i), d(i), g(i)
946 510 format(1x,i5,d16.6,2d14.3)
950 530 format(/24h inconsistent dimensions)
952 ! *** last card of itsum follows ***
954 !-----------------------------------------------------------------------------
955 subroutine litvmu(n, x, l, y)
957 ! *** solve (l**t)*x = y, where l is an n x n lower triangular
958 ! *** matrix stored compactly by rows. x and y may occupy the same
962 !al real(kind=8) :: x(n), l(1), y(n)
963 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
964 integer :: i, ii, ij, im1, i0, j, np1
965 real(kind=8) :: xi !el, zero
969 real(kind=8),parameter :: zero=0.d+0
980 if (i .le. 1) go to 999
982 if (xi .eq. zero) go to 30
986 x(j) = x(j) - xi*l(ij)
990 ! *** last card of litvmu follows ***
991 end subroutine litvmu
992 !-----------------------------------------------------------------------------
993 subroutine livmul(n, x, l, y)
995 ! *** solve l*x = y, where l is an n x n lower triangular
996 ! *** matrix stored compactly by rows. x and y may occupy the same
1000 !al real(kind=8) :: x(n), l(1), y(n)
1001 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
1003 !el real(kind=8) :: dotprd
1005 real(kind=8) :: t !el, zero
1009 real(kind=8),parameter :: zero=0.d+0
1013 if (y(k) .ne. zero) go to 20
1019 if (k .ge. n) go to 999
1022 t = dotprd(i-1, l(j+1), x)
1024 x(i) = (y(i) - t)/l(j)
1027 ! *** last card of livmul follows ***
1028 end subroutine livmul
1029 !-----------------------------------------------------------------------------
1030 subroutine parck(alg, d, iv, liv, lv, n, v)
1032 ! *** check ***sol (version 2.3) parameters, print changed values ***
1034 ! *** alg = 1 for regression, alg = 2 for general unconstrained opt.
1036 integer :: alg, liv, lv, n
1038 real(kind=8) :: d(n), v(lv)
1040 !el external rmdcon, vcopy, vdflt
1041 !el real(kind=8) :: rmdcon
1042 ! rmdcon -- returns machine-dependent constants.
1043 ! vcopy -- copies one vector to another.
1044 ! vdflt -- supplies default parameter values to v alone.
1049 ! *** local variables ***
1051 integer :: i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu
1052 integer :: ijmp, jlim(2), miniv(2), ndflt(2)
1054 ! integer varnm(2), sh(2)
1055 ! real cngd(3), dflt(3), vn(2,34), which(3)
1057 character(len=1) :: varnm(2), sh(2)
1058 character(len=4) :: cngd(3), dflt(3), vn(2,34), which(3)
1060 real(kind=8) :: big, machep, tiny, vk, vm(34), vx(34), zero
1062 ! *** iv and v subscripts ***
1064 !el integer algsav, dinit, dtype, dtype0, epslon, inits, ivneed,
1065 !el 1 lastiv, lastv, lmat, nextiv, nextv, nvdflt, oldn,
1066 !el 2 parprt, parsav, perm, prunit, vneed
1070 ! data algsav/51/, dinit/38/, dtype/16/, dtype0/54/, epslon/19/,
1071 ! 1 inits/25/, ivneed/3/, lastiv/44/, lastv/45/, lmat/42/,
1072 ! 2 nextiv/46/, nextv/47/, nvdflt/50/, oldn/38/, parprt/20/,
1073 ! 3 parsav/49/, perm/58/, prunit/21/, vneed/4/
1075 integer,parameter :: algsav=51, dinit=38, dtype=16, dtype0=54, epslon=19,&
1076 inits=25, ivneed=3, lastiv=44, lastv=45, lmat=42,&
1077 nextiv=46, nextv=47, nvdflt=50, oldn=38, parprt=20,&
1078 parsav=49, perm=58, prunit=21, vneed=4
1079 save big, machep, tiny
1082 data big/0.d+0/, machep/-1.d+0/, tiny/1.d+0/, zero/0.d+0/
1084 ! data vn(1,1),vn(2,1)/4hepsl,4hon../
1085 ! data vn(1,2),vn(2,2)/4hphmn,4hfc../
1086 ! data vn(1,3),vn(2,3)/4hphmx,4hfc../
1087 ! data vn(1,4),vn(2,4)/4hdecf,4hac../
1088 ! data vn(1,5),vn(2,5)/4hincf,4hac../
1089 ! data vn(1,6),vn(2,6)/4hrdfc,4hmn../
1090 ! data vn(1,7),vn(2,7)/4hrdfc,4hmx../
1091 ! data vn(1,8),vn(2,8)/4htune,4hr1../
1092 ! data vn(1,9),vn(2,9)/4htune,4hr2../
1093 ! data vn(1,10),vn(2,10)/4htune,4hr3../
1094 ! data vn(1,11),vn(2,11)/4htune,4hr4../
1095 ! data vn(1,12),vn(2,12)/4htune,4hr5../
1096 ! data vn(1,13),vn(2,13)/4hafct,4hol../
1097 ! data vn(1,14),vn(2,14)/4hrfct,4hol../
1098 ! data vn(1,15),vn(2,15)/4hxcto,4hl.../
1099 ! data vn(1,16),vn(2,16)/4hxfto,4hl.../
1100 ! data vn(1,17),vn(2,17)/4hlmax,4h0.../
1101 ! data vn(1,18),vn(2,18)/4hlmax,4hs.../
1102 ! data vn(1,19),vn(2,19)/4hscto,4hl.../
1103 ! data vn(1,20),vn(2,20)/4hdini,4ht.../
1104 ! data vn(1,21),vn(2,21)/4hdtin,4hit../
1105 ! data vn(1,22),vn(2,22)/4hd0in,4hit../
1106 ! data vn(1,23),vn(2,23)/4hdfac,4h..../
1107 ! data vn(1,24),vn(2,24)/4hdltf,4hdc../
1108 ! data vn(1,25),vn(2,25)/4hdltf,4hdj../
1109 ! data vn(1,26),vn(2,26)/4hdelt,4ha0../
1110 ! data vn(1,27),vn(2,27)/4hfuzz,4h..../
1111 ! data vn(1,28),vn(2,28)/4hrlim,4hit../
1112 ! data vn(1,29),vn(2,29)/4hcosm,4hin../
1113 ! data vn(1,30),vn(2,30)/4hhube,4hrc../
1114 ! data vn(1,31),vn(2,31)/4hrspt,4hol../
1115 ! data vn(1,32),vn(2,32)/4hsigm,4hin../
1116 ! data vn(1,33),vn(2,33)/4heta0,4h..../
1117 ! data vn(1,34),vn(2,34)/4hbias,4h..../
1119 data vn(1,1),vn(2,1)/'epsl','on..'/
1120 data vn(1,2),vn(2,2)/'phmn','fc..'/
1121 data vn(1,3),vn(2,3)/'phmx','fc..'/
1122 data vn(1,4),vn(2,4)/'decf','ac..'/
1123 data vn(1,5),vn(2,5)/'incf','ac..'/
1124 data vn(1,6),vn(2,6)/'rdfc','mn..'/
1125 data vn(1,7),vn(2,7)/'rdfc','mx..'/
1126 data vn(1,8),vn(2,8)/'tune','r1..'/
1127 data vn(1,9),vn(2,9)/'tune','r2..'/
1128 data vn(1,10),vn(2,10)/'tune','r3..'/
1129 data vn(1,11),vn(2,11)/'tune','r4..'/
1130 data vn(1,12),vn(2,12)/'tune','r5..'/
1131 data vn(1,13),vn(2,13)/'afct','ol..'/
1132 data vn(1,14),vn(2,14)/'rfct','ol..'/
1133 data vn(1,15),vn(2,15)/'xcto','l...'/
1134 data vn(1,16),vn(2,16)/'xfto','l...'/
1135 data vn(1,17),vn(2,17)/'lmax','0...'/
1136 data vn(1,18),vn(2,18)/'lmax','s...'/
1137 data vn(1,19),vn(2,19)/'scto','l...'/
1138 data vn(1,20),vn(2,20)/'dini','t...'/
1139 data vn(1,21),vn(2,21)/'dtin','it..'/
1140 data vn(1,22),vn(2,22)/'d0in','it..'/
1141 data vn(1,23),vn(2,23)/'dfac','....'/
1142 data vn(1,24),vn(2,24)/'dltf','dc..'/
1143 data vn(1,25),vn(2,25)/'dltf','dj..'/
1144 data vn(1,26),vn(2,26)/'delt','a0..'/
1145 data vn(1,27),vn(2,27)/'fuzz','....'/
1146 data vn(1,28),vn(2,28)/'rlim','it..'/
1147 data vn(1,29),vn(2,29)/'cosm','in..'/
1148 data vn(1,30),vn(2,30)/'hube','rc..'/
1149 data vn(1,31),vn(2,31)/'rspt','ol..'/
1150 data vn(1,32),vn(2,32)/'sigm','in..'/
1151 data vn(1,33),vn(2,33)/'eta0','....'/
1152 data vn(1,34),vn(2,34)/'bias','....'/
1155 data vm(1)/1.0d-3/, vm(2)/-0.99d+0/, vm(3)/1.0d-3/, vm(4)/1.0d-2/,&
1156 vm(5)/1.2d+0/, vm(6)/1.d-2/, vm(7)/1.2d+0/, vm(8)/0.d+0/,&
1157 vm(9)/0.d+0/, vm(10)/1.d-3/, vm(11)/-1.d+0/, vm(13)/0.d+0/,&
1158 vm(15)/0.d+0/, vm(16)/0.d+0/, vm(19)/0.d+0/, vm(20)/-10.d+0/,&
1159 vm(21)/0.d+0/, vm(22)/0.d+0/, vm(23)/0.d+0/, vm(27)/1.01d+0/,&
1160 vm(28)/1.d+10/, vm(30)/0.d+0/, vm(31)/0.d+0/, vm(32)/0.d+0/,&
1162 data vx(1)/0.9d+0/, vx(2)/-1.d-3/, vx(3)/1.d+1/, vx(4)/0.8d+0/,&
1163 vx(5)/1.d+2/, vx(6)/0.8d+0/, vx(7)/1.d+2/, vx(8)/0.5d+0/,&
1164 vx(9)/0.5d+0/, vx(10)/1.d+0/, vx(11)/1.d+0/, vx(14)/0.1d+0/,&
1165 vx(15)/1.d+0/, vx(16)/1.d+0/, vx(19)/1.d+0/, vx(23)/1.d+0/,&
1166 vx(24)/1.d+0/, vx(25)/1.d+0/, vx(26)/1.d+0/, vx(27)/1.d+10/,&
1167 vx(29)/1.d+0/, vx(31)/1.d+0/, vx(32)/1.d+0/, vx(33)/1.d+0/,&
1171 ! data varnm(1)/1hp/, varnm(2)/1hn/, sh(1)/1hs/, sh(2)/1hh/
1172 ! data cngd(1),cngd(2),cngd(3)/4h---c,4hhang,4hed v/,
1173 ! 1 dflt(1),dflt(2),dflt(3)/4hnond,4hefau,4hlt v/
1175 data varnm(1)/'p'/, varnm(2)/'n'/, sh(1)/'s'/, sh(2)/'h'/
1176 data cngd(1),cngd(2),cngd(3)/'---c','hang','ed v'/,&
1177 dflt(1),dflt(2),dflt(3)/'nond','efau','lt v'/
1179 data ijmp/33/, jlim(1)/0/, jlim(2)/24/, ndflt(1)/32/, ndflt(2)/25/
1180 data miniv(1)/80/, miniv(2)/59/
1182 !............................... body ................................
1185 if (prunit .le. liv) pu = iv(prunit)
1186 if (alg .lt. 1 .or. alg .gt. 2) go to 340
1187 if (iv(1) .eq. 0) call deflt(alg, iv, liv, lv, v)
1189 if (iv1 .ne. 13 .and. iv1 .ne. 12) go to 10
1191 if (perm .le. liv) miv1 = max0(miv1, iv(perm) - 1)
1192 if (ivneed .le. liv) miv2 = miv1 + max0(iv(ivneed), 0)
1193 if (lastiv .le. liv) iv(lastiv) = miv2
1194 if (liv .lt. miv1) go to 300
1196 iv(lastv) = max0(iv(vneed), 0) + iv(lmat) - 1
1198 if (liv .lt. miv2) go to 300
1199 if (lv .lt. iv(lastv)) go to 320
1200 10 if (alg .eq. iv(algsav)) go to 30
1201 if (pu .ne. 0) write(pu,20) alg, iv(algsav)
1202 20 format(/39h the first parameter to deflt should be,i3,&
1206 30 if (iv1 .lt. 12 .or. iv1 .gt. 14) go to 60
1207 if (n .ge. 1) go to 50
1209 if (pu .eq. 0) go to 999
1210 write(pu,40) varnm(alg), n
1211 40 format(/8h /// bad,a1,2h =,i5)
1213 50 if (iv1 .ne. 14) iv(nextiv) = iv(perm)
1214 if (iv1 .ne. 14) iv(nextv) = iv(lmat)
1215 if (iv1 .eq. 13) go to 999
1216 k = iv(parsav) - epslon
1217 call vdflt(alg, lv-k, v(k+1))
1218 iv(dtype0) = 2 - alg
1224 60 if (n .eq. iv(oldn)) go to 80
1226 if (pu .eq. 0) go to 999
1227 write(pu,70) varnm(alg), iv(oldn), n
1228 70 format(/5h /// ,1a1,14h changed from ,i5,4h to ,i5)
1231 80 if (iv1 .le. 11 .and. iv1 .ge. 1) go to 100
1233 if (pu .ne. 0) write(pu,90) iv1
1234 90 format(/13h /// iv(1) =,i5,28h should be between 0 and 14.)
1237 100 which(1) = cngd(1)
1241 110 if (iv1 .eq. 14) iv1 = 12
1242 if (big .gt. tiny) go to 120
1269 do 150 l = 1, ndfalt
1271 if (vk .ge. vm(i) .and. vk .le. vx(i)) go to 140
1273 if (pu .ne. 0) write(pu,130) vn(1,i), vn(2,i), k, vk,&
1275 130 format(/6h /// ,2a4,5h.. v(,i2,3h) =,d11.3,7h should,&
1276 11h be between,d11.3,4h and,d11.3)
1279 if (i .eq. j) i = ijmp
1282 if (iv(nvdflt) .eq. ndfalt) go to 170
1284 if (pu .eq. 0) go to 999
1285 write(pu,160) iv(nvdflt), ndfalt
1286 160 format(/13h iv(nvdflt) =,i5,13h rather than ,i5)
1288 170 if ((iv(dtype) .gt. 0 .or. v(dinit) .gt. zero) .and. iv1 .eq. 12) &
1291 if (d(i) .gt. zero) go to 190
1293 if (pu .ne. 0) write(pu,180) i, d(i)
1294 180 format(/8h /// d(,i3,3h) =,d11.3,19h should be positive)
1296 200 if (m .eq. 0) go to 210
1300 210 if (pu .eq. 0 .or. iv(parprt) .eq. 0) go to 999
1301 if (iv1 .ne. 12 .or. iv(inits) .eq. alg-1) go to 230
1303 write(pu,220) sh(alg), iv(inits)
1304 220 format(/22h nondefault values..../5h init,a1,14h..... iv(25) =,&
1306 230 if (iv(dtype) .eq. iv(dtype0)) go to 250
1307 if (m .eq. 0) write(pu,260) which
1309 write(pu,240) iv(dtype)
1310 240 format(20h dtype..... iv(16) =,i3)
1316 do 290 ii = 1, ndfalt
1317 if (v(k) .eq. v(l)) go to 280
1318 if (m .eq. 0) write(pu,260) which
1319 260 format(/1h ,3a4,9halues..../)
1321 write(pu,270) vn(1,i), vn(2,i), k, v(k)
1322 270 format(1x,2a4,5h.. v(,i2,3h) =,d15.7)
1326 if (i .eq. j) i = ijmp
1329 iv(dtype0) = iv(dtype)
1331 call vcopy(iv(nvdflt), v(parsv1), v(epslon))
1335 if (pu .eq. 0) go to 999
1336 write(pu,310) liv, miv2
1337 310 format(/10h /// liv =,i5,17h must be at least,i5)
1338 if (liv .lt. miv1) go to 999
1339 if (lv .lt. iv(lastv)) go to 320
1343 if (pu .eq. 0) go to 999
1344 write(pu,330) lv, iv(lastv)
1345 330 format(/9h /// lv =,i5,17h must be at least,i5)
1349 if (pu .eq. 0) go to 999
1351 350 format(/10h /// alg =,i5,15h must be 1 or 2)
1354 ! *** last card of parck follows ***
1355 end subroutine parck
1356 !-----------------------------------------------------------------------------
1357 real(kind=8) function reldst(p, d, x, x0)
1359 ! *** compute and return relative difference between x and x0 ***
1360 ! *** nl2sol version 2.2 ***
1363 real(kind=8) :: d(p), x(p), x0(p)
1365 !el real(kind=8) :: dabs
1368 real(kind=8) :: emax, t, xmax !el, zero
1372 real(kind=8),parameter :: zero=0.d+0
1378 t = dabs(d(i) * (x(i) - x0(i)))
1379 if (emax .lt. t) emax = t
1380 t = d(i) * (dabs(x(i)) + dabs(x0(i)))
1381 if (xmax .lt. t) xmax = t
1384 if (xmax .gt. zero) reldst = emax / xmax
1386 ! *** last card of reldst follows ***
1388 !-----------------------------------------------------------------------------
1389 subroutine vaxpy(p, w, a, x, y)
1391 ! *** set w = a*x + y -- w, x, y = p-vectors, a = scalar ***
1394 real(kind=8) :: a, w(p), x(p), y(p)
1399 10 w(i) = a*x(i) + y(i)
1401 end subroutine vaxpy
1402 !-----------------------------------------------------------------------------
1403 subroutine vcopy(p, y, x)
1405 ! *** set y = x, where x and y are p-vectors ***
1408 real(kind=8) :: x(p), y(p)
1415 end subroutine vcopy
1416 !-----------------------------------------------------------------------------
1417 subroutine vdflt(alg, lv, v)
1419 ! *** supply ***sol (version 2.3) default values to v ***
1421 ! *** alg = 1 means regression constants.
1422 ! *** alg = 2 means general unconstrained optimization constants.
1424 integer :: alg, l,lv
1425 real(kind=8) :: v(lv)
1427 !el real(kind=8) :: dmax1
1430 !el real(kind=8) :: rmdcon
1431 ! rmdcon... returns machine-dependent constants
1433 real(kind=8) :: machep, mepcrt, sqteps !el one, three
1435 ! *** subscripts for v ***
1437 !el integer afctol, bias, cosmin, decfac, delta0, dfac, dinit, dltfdc,
1438 !el 1 dltfdj, dtinit, d0init, epslon, eta0, fuzz, huberc,
1439 !el 2 incfac, lmax0, lmaxs, phmnfc, phmxfc, rdfcmn, rdfcmx,
1440 !el 3 rfctol, rlimit, rsptol, sctol, sigmin, tuner1, tuner2,
1441 !el 4 tuner3, tuner4, tuner5, xctol, xftol
1444 ! data one/1.d+0/, three/3.d+0/
1446 real(kind=8),parameter :: one=1.d+0, three=3.d+0
1449 ! *** v subscript values ***
1452 ! data afctol/31/, bias/43/, cosmin/47/, decfac/22/, delta0/44/,
1453 ! 1 dfac/41/, dinit/38/, dltfdc/42/, dltfdj/43/, dtinit/39/,
1454 ! 2 d0init/40/, epslon/19/, eta0/42/, fuzz/45/, huberc/48/,
1455 ! 3 incfac/23/, lmax0/35/, lmaxs/36/, phmnfc/20/, phmxfc/21/,
1456 ! 4 rdfcmn/24/, rdfcmx/25/, rfctol/32/, rlimit/46/, rsptol/49/,
1457 ! 5 sctol/37/, sigmin/50/, tuner1/26/, tuner2/27/, tuner3/28/,
1458 ! 6 tuner4/29/, tuner5/30/, xctol/33/, xftol/34/
1460 integer,parameter :: afctol=31, bias=43, cosmin=47, decfac=22, delta0=44,&
1461 dfac=41, dinit=38, dltfdc=42, dltfdj=43, dtinit=39,&
1462 d0init=40, epslon=19, eta0=42, fuzz=45, huberc=48,&
1463 incfac=23, lmax0=35, lmaxs=36, phmnfc=20, phmxfc=21,&
1464 rdfcmn=24, rdfcmx=25, rfctol=32, rlimit=46, rsptol=49,&
1465 sctol=37, sigmin=50, tuner1=26, tuner2=27, tuner3=28,&
1466 tuner4=29, tuner5=30, xctol=33, xftol=34
1469 !------------------------------- body --------------------------------
1473 if (machep .gt. 1.d-10) v(afctol) = machep**2
1479 mepcrt = machep ** (one/three)
1489 v(rfctol) = dmax1(1.d-10, mepcrt**2)
1490 v(sctol) = v(rfctol)
1497 v(xftol) = 1.d+2 * machep
1499 if (alg .ge. 2) go to 10
1501 ! *** regression values
1503 v(cosmin) = dmax1(1.d-6, 1.d+2 * machep)
1509 v(rlimit) = rmdcon(5)
1514 ! *** general optimization values
1518 v(eta0) = 1.0d+3 * machep
1521 ! *** last card of vdflt follows ***
1522 end subroutine vdflt
1523 !-----------------------------------------------------------------------------
1524 subroutine vscopy(p, y, s)
1526 ! *** set p-vector y to scalar s ***
1529 real(kind=8) :: s, y(p)
1536 end subroutine vscopy
1537 !-----------------------------------------------------------------------------
1538 real(kind=8) function v2norm(p, x)
1540 ! *** return the 2-norm of the p-vector x, taking ***
1541 ! *** care to avoid the most likely underflows. ***
1544 real(kind=8) :: x(p)
1547 real(kind=8) :: r, scale, sqteta, t, xi !el, one, zero
1549 !el real(kind=8) :: dabs, dsqrt
1552 !el real(kind=8) :: rmdcon
1555 ! data one/1.d+0/, zero/0.d+0/
1557 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
1562 if (p .gt. 0) go to 10
1566 if (x(i) .ne. zero) go to 30
1571 30 scale = dabs(x(i))
1572 if (i .lt. p) go to 40
1576 if (sqteta .eq. zero) sqteta = rmdcon(2)
1578 ! *** sqteta is (slightly larger than) the square root of the
1579 ! *** smallest positive floating point number on the machine.
1580 ! *** the tests involving sqteta are done to prevent underflows.
1585 if (xi .gt. scale) go to 50
1587 if (r .gt. sqteta) t = t + r*r
1590 if (r .le. sqteta) r = zero
1595 v2norm = scale * dsqrt(t)
1597 ! *** last card of v2norm follows ***
1599 !-----------------------------------------------------------------------------
1600 subroutine humsl(n,d,x,calcf,calcgh,iv,liv,lv,v,uiparm,urparm,ufparm)
1602 ! *** minimize general unconstrained objective function using ***
1603 ! *** (analytic) gradient and hessian provided by the caller. ***
1605 integer :: liv, lv, n
1606 integer :: iv(liv), uiparm(1)
1607 real(kind=8) :: d(n), x(n), v(lv), urparm(1)
1608 real(kind=8),external :: ufparm
1609 ! dimension v(78 + n*(n+12)), uiparm(*), urparm(*)
1610 external :: calcf, calcgh
1612 !------------------------------ discussion ---------------------------
1614 ! this routine is like sumsl, except that the subroutine para-
1615 ! meter calcg of sumsl (which computes the gradient of the objec-
1616 ! tive function) is replaced by the subroutine parameter calcgh,
1617 ! which computes both the gradient and (lower triangle of the)
1618 ! hessian of the objective function. the calling sequence is...
1619 ! call calcgh(n, x, nf, g, h, uiparm, urparm, ufparm)
1620 ! parameters n, x, nf, g, uiparm, urparm, and ufparm are the same
1621 ! as for sumsl, while h is an array of length n*(n+1)/2 in which
1622 ! calcgh must store the lower triangle of the hessian at x. start-
1623 ! ing at h(1), calcgh must store the hessian entries in the order
1624 ! (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ...
1625 ! the value printed (by itsum) in the column labelled stppar
1626 ! is the levenberg-marquardt used in computing the current step.
1627 ! zero means a full newton step. if the special case described in
1628 ! ref. 1 is detected, then stppar is negated. the value printed
1629 ! in the column labelled npreldf is zero if the current hessian
1630 ! is not positive definite.
1631 ! it sometimes proves worthwhile to let d be determined from the
1632 ! diagonal of the hessian matrix by setting iv(dtype) = 1 and
1633 ! v(dinit) = 0. the following iv and v components are relevant...
1635 ! iv(dtol)..... iv(59) gives the starting subscript in v of the dtol
1636 ! array used when d is updated. (iv(dtol) can be
1637 ! initialized by calling humsl with iv(1) = 13.)
1638 ! iv(dtype).... iv(16) tells how the scale vector d should be chosen.
1639 ! iv(dtype) .le. 0 means that d should not be updated, and
1640 ! iv(dtype) .ge. 1 means that d should be updated as
1641 ! described below with v(dfac). default = 0.
1642 ! v(dfac)..... v(41) and the dtol and d0 arrays (see v(dtinit) and
1643 ! v(d0init)) are used in updating the scale vector d when
1644 ! iv(dtype) .gt. 0. (d is initialized according to
1645 ! v(dinit), described in sumsl.) let
1646 ! d1(i) = max(sqrt(abs(h(i,i))), v(dfac)*d(i)),
1647 ! where h(i,i) is the i-th diagonal element of the current
1648 ! hessian. if iv(dtype) = 1, then d(i) is set to d1(i)
1649 ! unless d1(i) .lt. dtol(i), in which case d(i) is set to
1650 ! max(d0(i), dtol(i)).
1651 ! if iv(dtype) .ge. 2, then d is updated during the first
1652 ! iteration as for iv(dtype) = 1 (after any initialization
1653 ! due to v(dinit)) and is left unchanged thereafter.
1655 ! v(dtinit)... v(39), if positive, is the value to which all components
1656 ! of the dtol array (see v(dfac)) are initialized. if
1657 ! v(dtinit) = 0, then it is assumed that the caller has
1658 ! stored dtol in v starting at v(iv(dtol)).
1660 ! v(d0init)... v(40), if positive, is the value to which all components
1661 ! of the d0 vector (see v(dfac)) are initialized. if
1662 ! v(dfac) = 0, then it is assumed that the caller has
1663 ! stored d0 in v starting at v(iv(dtol)+n). default = 1.0.
1667 ! 1. gay, d.m. (1981), computing optimal locally constrained steps,
1668 ! siam j. sci. statist. comput. 2, pp. 186-197.
1672 ! coded by david m. gay (winter 1980). revised sept. 1982.
1673 ! this subroutine was written in connection with research supported
1674 ! in part by the national science foundation under grants
1675 ! mcs-7600324 and mcs-7906671.
1677 !---------------------------- declarations ---------------------------
1679 !el external deflt, humit
1681 ! deflt... provides default input values for iv and v.
1682 ! humit... reverse-communication routine that does humsl algorithm.
1684 integer :: g1, h1, iv1, lh, nf
1687 ! *** subscripts for iv ***
1689 !el integer g, h, nextv, nfcall, nfgcal, toobig, vneed
1692 ! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, h/56/, toobig/2/,
1695 integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28, h=56,&
1699 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
1701 lh = n * (n + 1) / 2
1702 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1703 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
1704 iv(vneed) = iv(vneed) + n*(n+3)/2
1706 if (iv1 .eq. 14) go to 10
1707 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
1710 if (iv1 .eq. 12) iv(1) = 13
1716 20 call humit(d, f, v(g1), v(h1), iv, lh, liv, lv, n, v, x)
1717 if (iv(1) - 2) 30, 40, 50
1720 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
1721 if (nf .le. 0) iv(toobig) = 1
1724 40 call calcgh(n, x, iv(nfgcal), v(g1), v(h1), uiparm, urparm,&
1728 50 if (iv(1) .ne. 14) go to 999
1730 ! *** storage allocation
1734 iv(nextv) = iv(h) + n*(n+1)/2
1735 if (iv1 .ne. 13) go to 10
1738 ! *** last card of humsl follows ***
1739 end subroutine humsl
1740 !-----------------------------------------------------------------------------
1741 subroutine humit(d, fx, g, h, iv, lh, liv, lv, n, v, x)
1743 ! *** carry out humsl (unconstrained minimization) iterations, using
1744 ! *** hessian matrix provided by the caller.
1746 use control, only:stopx
1748 ! *** parameter declarations ***
1750 integer :: lh, liv, lv, n
1752 real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n)
1754 !-------------------------- parameter usage --------------------------
1756 ! d.... scale vector.
1757 ! fx... function value.
1758 ! g.... gradient vector.
1759 ! h.... lower triangle of the hessian, stored rowwise.
1760 ! iv... integer value array.
1761 ! lh... length of h = p*(p+1)/2.
1762 ! liv.. length of iv (at least 60).
1763 ! lv... length of v (at least 78 + n*(n+21)/2).
1764 ! n.... number of variables (components in x and g).
1765 ! v.... floating-point value array.
1766 ! x.... parameter vector.
1768 ! *** discussion ***
1770 ! parameters iv, n, v, and x are the same as the corresponding
1771 ! ones to humsl (which see), except that v can be shorter (since
1772 ! the part of v that humsl uses for storing g and h is not needed).
1773 ! moreover, compared with humsl, iv(1) may have the two additional
1774 ! output values 1 and 2, which are explained below, as is the use
1775 ! of iv(toobig) and iv(nfgcal). the value iv(g), which is an
1776 ! output value from humsl, is not referenced by humit or the
1777 ! subroutines it calls.
1779 ! iv(1) = 1 means the caller should set fx to f(x), the function value
1780 ! at x, and call humit again, having changed none of the
1781 ! other parameters. an exception occurs if f(x) cannot be
1782 ! computed (e.g. if overflow would occur), which may happen
1783 ! because of an oversized step. in this case the caller
1784 ! should set iv(toobig) = iv(2) to 1, which will cause
1785 ! humit to ignore fx and try a smaller step. the para-
1786 ! meter nf that humsl passes to calcf (for possible use by
1787 ! calcgh) is a copy of iv(nfcall) = iv(6).
1788 ! iv(1) = 2 means the caller should set g to g(x), the gradient of f at
1789 ! x, and h to the lower triangle of h(x), the hessian of f
1790 ! at x, and call humit again, having changed none of the
1791 ! other parameters except perhaps the scale vector d.
1792 ! the parameter nf that humsl passes to calcg is
1793 ! iv(nfgcal) = iv(7). if g(x) and h(x) cannot be evaluated,
1794 ! then the caller may set iv(nfgcal) to 0, in which case
1795 ! humit will return with iv(1) = 65.
1796 ! note -- humit overwrites h with the lower triangle
1797 ! of diag(d)**-1 * h(x) * diag(d)**-1.
1801 ! coded by david m. gay (winter 1980). revised sept. 1982.
1802 ! this subroutine was written in connection with research supported
1803 ! in part by the national science foundation under grants
1804 ! mcs-7600324 and mcs-7906671.
1806 ! (see sumsl and humsl for references.)
1808 !+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
1810 ! *** local variables ***
1812 integer :: dg1, dummy, i, j, k, l, lstgst, nn1o2, step1,&
1818 !el real(kind=8) :: one, onep2, zero
1820 ! *** no intrinsic functions ***
1822 ! *** external functions and subroutines ***
1824 !el external assst, deflt, dotprd, dupdu, gqtst, itsum, parck,
1825 !el 1 reldst, slvmul, stopx, vaxpy, vcopy, vscopy, v2norm
1827 !el real(kind=8) :: dotprd, reldst, v2norm
1829 ! assst.... assesses candidate step.
1830 ! deflt.... provides default iv and v input values.
1831 ! dotprd... returns inner product of two vectors.
1832 ! dupdu.... updates scale vector d.
1833 ! gqtst.... computes optimally locally constrained step.
1834 ! itsum.... prints iteration summary and info on initial and final x.
1835 ! parck.... checks validity of input iv and v values.
1836 ! reldst... computes v(reldx) = relative step size.
1837 ! slvmul... multiplies symmetric matrix times vector, given the lower
1838 ! triangle of the matrix.
1839 ! stopx.... returns .true. if the break key has been pressed.
1840 ! vaxpy.... computes scalar times one vector plus another.
1841 ! vcopy.... copies one vector to another.
1842 ! vscopy... sets all elements of a vector to a scalar.
1843 ! v2norm... returns the 2-norm of a vector.
1845 ! *** subscripts for iv and v ***
1847 !el integer cnvcod, dg, dgnorm, dinit, dstnrm, dtinit, dtol,
1848 !el 1 dtype, d0init, f, f0, fdif, gtstep, incfac, irc, kagqt,
1849 !el 2 lmat, lmax0, lmaxs, mode, model, mxfcal, mxiter, nextv,
1850 !el 3 nfcall, nfgcal, ngcall, niter, preduc, radfac, radinc,
1851 !el 4 radius, rad0, reldx, restor, step, stglim, stlstg, stppar,
1852 !el 5 toobig, tuner4, tuner5, vneed, w, xirc, x0
1854 ! *** iv subscript values ***
1857 ! data cnvcod/55/, dg/37/, dtol/59/, dtype/16/, irc/29/, kagqt/33/,
1858 ! 1 lmat/42/, mode/35/, model/5/, mxfcal/17/, mxiter/18/,
1859 ! 2 nextv/47/, nfcall/6/, nfgcal/7/, ngcall/30/, niter/31/,
1860 ! 3 radinc/8/, restor/9/, step/40/, stglim/11/, stlstg/41/,
1861 ! 4 toobig/2/, vneed/4/, w/34/, xirc/13/, x0/43/
1863 integer,parameter :: cnvcod=55, dg=37, dtol=59, dtype=16, irc=29, kagqt=33,&
1864 lmat=42, mode=35, model=5, mxfcal=17, mxiter=18,&
1865 nextv=47, nfcall=6, nfgcal=7, ngcall=30, niter=31,&
1866 radinc=8, restor=9, step=40, stglim=11, stlstg=41,&
1867 toobig=2, vneed=4, w=34, xirc=13, x0=43
1870 ! *** v subscript values ***
1873 ! data dgnorm/1/, dinit/38/, dstnrm/2/, dtinit/39/, d0init/40/,
1874 ! 1 f/10/, f0/13/, fdif/11/, gtstep/4/, incfac/23/, lmax0/35/,
1875 ! 2 lmaxs/36/, preduc/7/, radfac/16/, radius/8/, rad0/9/,
1876 ! 3 reldx/17/, stppar/5/, tuner4/29/, tuner5/30/
1878 integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dtinit=39, d0init=40,&
1879 f=10, f0=13, fdif=11, gtstep=4, incfac=23, lmax0=35,&
1880 lmaxs=36, preduc=7, radfac=16, radius=8, rad0=9,&
1881 reldx=17, stppar=5, tuner4=29, tuner5=30
1885 ! data one/1.d+0/, onep2/1.2d+0/, zero/0.d+0/
1887 real(kind=8),parameter :: one=1.d+0, onep2=1.2d+0, zero=0.d+0
1890 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
1893 if (i .eq. 1) go to 30
1894 if (i .eq. 2) go to 40
1896 ! *** check validity of iv and v input values ***
1898 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1899 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
1900 iv(vneed) = iv(vneed) + n*(n+21)/2 + 7
1901 call parck(2, d, iv, liv, lv, n, v)
1903 if (i .gt. 12) go to 999
1904 nn1o2 = n * (n + 1) / 2
1905 if (lh .ge. nn1o2) go to (210,210,210,210,210,210,160,120,160,&
1910 ! *** storage allocation ***
1912 10 iv(dtol) = iv(lmat) + nn1o2
1913 iv(x0) = iv(dtol) + 2*n
1914 iv(step) = iv(x0) + n
1915 iv(stlstg) = iv(step) + n
1916 iv(dg) = iv(stlstg) + n
1918 iv(nextv) = iv(w) + 4*n + 7
1919 if (iv(1) .ne. 13) go to 20
1923 ! *** initialization ***
1937 if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
1939 if (v(dtinit) .gt. zero) call vscopy(n, v(k), v(dtinit))
1941 if (v(d0init) .gt. zero) call vscopy(n, v(k), v(d0init))
1946 if (iv(mode) .ge. 0) go to 210
1948 if (iv(toobig) .eq. 0) go to 999
1952 ! *** make sure gradient could be computed ***
1954 40 if (iv(nfgcal) .ne. 0) go to 50
1958 ! *** update the scale vector d ***
1961 if (iv(dtype) .le. 0) go to 70
1969 call dupdu(d, v(dg1), iv, liv, lv, n, v)
1971 ! *** compute scaled gradient and its norm ***
1979 v(dgnorm) = v2norm(n, v(dg1))
1981 ! *** compute scaled hessian ***
1987 h(k) = t * h(k) / d(j)
1992 if (iv(cnvcod) .ne. 0) go to 340
1993 if (iv(mode) .eq. 0) go to 300
1995 ! *** allow first step to have scaled 2-norm at most v(lmax0) ***
1997 v(radius) = v(lmax0)
2002 !----------------------------- main loop -----------------------------
2005 ! *** print iteration summary, check iteration limit ***
2007 110 call itsum(d, g, iv, liv, lv, n, v, x)
2009 if (k .lt. iv(mxiter)) go to 130
2013 130 iv(niter) = k + 1
2015 ! *** initialize for start of next iteration ***
2023 ! *** copy x to x0 ***
2025 call vcopy(n, v(x01), x)
2027 ! *** update radius ***
2029 if (k .eq. 0) go to 150
2036 v(radius) = v(radfac) * v2norm(n, v(step1))
2038 ! *** check stopx and function evaluation limit ***
2042 150 if (.not. stopx(dummy)) go to 170
2046 ! *** come here when restarting after func. eval. limit or stopx.
2048 160 if (v(f) .ge. v(f0)) go to 170
2053 170 if (iv(nfcall) .lt. iv(mxfcal)) go to 190
2055 180 if (v(f) .ge. v(f0)) go to 350
2057 ! *** in case of stopx or function evaluation limit with
2058 ! *** improved v(f), evaluate the gradient at x.
2063 !. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
2065 190 step1 = iv(step)
2069 call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1))
2070 if (iv(irc) .eq. 6) go to 210
2072 ! *** check whether evaluating f(x0 + step) looks worthwhile ***
2074 if (v(dstnrm) .le. zero) go to 210
2075 if (iv(irc) .ne. 5) go to 200
2076 if (v(radfac) .le. one) go to 200
2077 if (v(preduc) .le. onep2 * v(fdif)) go to 210
2079 ! *** compute f(x0 + step) ***
2083 call vaxpy(n, x, one, v(step1), v(x01))
2084 iv(nfcall) = iv(nfcall) + 1
2089 !. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
2092 v(reldx) = reldst(n, d, x, v(x01))
2093 call assst(iv, liv, lv, v)
2096 if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
2097 if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
2098 if (iv(restor) .ne. 3) go to 220
2099 call vcopy(n, v(step1), v(lstgst))
2100 call vaxpy(n, x, one, v(step1), v(x01))
2101 v(reldx) = reldst(n, d, x, v(x01))
2104 go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k
2106 ! *** recompute step with new radius ***
2108 230 v(radius) = v(radfac) * v(dstnrm)
2111 ! *** compute step of length v(lmaxs) for singular convergence test.
2113 240 v(radius) = v(lmaxs)
2116 ! *** convergence or false convergence ***
2118 250 iv(cnvcod) = k - 4
2119 if (v(f) .ge. v(f0)) go to 340
2120 if (iv(xirc) .eq. 14) go to 340
2123 !. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
2125 260 if (iv(irc) .ne. 3) go to 290
2128 ! *** prepare for gradient tests ***
2129 ! *** set temp1 = hessian * step + g(x0)
2130 ! *** = diag(d) * (h * step + g(x0))
2132 ! use x0 vector as temporary.
2135 v(k) = d(i) * v(step1)
2139 call slvmul(n, v(temp1), h, v(x01))
2141 v(temp1) = d(i) * v(temp1) + g(i)
2145 ! *** compute gradient and hessian ***
2147 290 iv(ngcall) = iv(ngcall) + 1
2152 if (iv(irc) .ne. 3) go to 110
2154 ! *** set v(radfac) by gradient tests ***
2159 ! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
2163 v(k) = (v(k) - g(i)) / d(i)
2167 ! *** do gradient tests ***
2169 if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320
2170 if (dotprd(n, g, v(step1)) &
2171 .ge. v(gtstep) * v(tuner5)) go to 110
2172 320 v(radfac) = v(incfac)
2175 !. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
2177 ! *** bad parameters to assess ***
2182 ! *** print summary of final iteration and other requested items ***
2184 340 iv(1) = iv(cnvcod)
2186 350 call itsum(d, g, iv, liv, lv, n, v, x)
2190 ! *** last card of humit follows ***
2191 end subroutine humit
2192 !-----------------------------------------------------------------------------
2193 subroutine dupdu(d, hdiag, iv, liv, lv, n, v)
2195 ! *** update scale vector d for humsl ***
2197 ! *** parameter declarations ***
2199 integer :: liv, lv, n
2201 real(kind=8) :: d(n), hdiag(n), v(lv)
2203 ! *** local variables ***
2205 integer :: dtoli, d0i, i
2206 real(kind=8) :: t, vdfac
2208 ! *** intrinsic functions ***
2210 !el real(kind=8) :: dabs, dmax1, dsqrt
2212 ! *** subscripts for iv and v ***
2214 !el integer :: dfac, dtol, dtype, niter
2216 ! data dfac/41/, dtol/59/, dtype/16/, niter/31/
2218 integer,parameter :: dfac=41, dtol=59, dtype=16, niter=31
2221 !------------------------------- body --------------------------------
2224 if (i .eq. 1) go to 10
2225 if (iv(niter) .gt. 0) go to 999
2231 t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i))
2232 if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i))
2239 ! *** last card of dupdu follows ***
2240 end subroutine dupdu
2241 !-----------------------------------------------------------------------------
2242 subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w)
2244 ! *** compute goldfeld-quandt-trotter step by more-hebden technique ***
2245 ! *** (nl2sol version 2.2), modified a la more and sorensen ***
2247 ! *** parameter declarations ***
2250 !al real(kind=8) :: d(p), dig(p), dihdi(1), l(1), v(21), step(p),
2252 real(kind=8) :: d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),&
2253 v(21), step(p),w(4*p+7)
2254 ! dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7)
2256 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2260 ! given the (compactly stored) lower triangle of a scaled
2261 ! hessian (approximation) and a nonzero scaled gradient vector,
2262 ! this subroutine computes a goldfeld-quandt-trotter step of
2263 ! approximate length v(radius) by the more-hebden technique. in
2264 ! other words, step is computed to (approximately) minimize
2265 ! psi(step) = (g**t)*step + 0.5*(step**t)*h*step such that the
2266 ! 2-norm of d*step is at most (approximately) v(radius), where
2267 ! g is the gradient, h is the hessian, and d is a diagonal
2268 ! scale matrix whose diagonal is stored in the parameter d.
2269 ! (gqtst assumes dig = d**-1 * g and dihdi = d**-1 * h * d**-1.)
2271 ! *** parameter description ***
2273 ! d (in) = the scale vector, i.e. the diagonal of the scale
2274 ! matrix d mentioned above under purpose.
2275 ! dig (in) = the scaled gradient vector, d**-1 * g. if g = 0, then
2276 ! step = 0 and v(stppar) = 0 are returned.
2277 ! dihdi (in) = lower triangle of the scaled hessian (approximation),
2278 ! i.e., d**-1 * h * d**-1, stored compactly by rows., i.e.,
2279 ! in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc.
2280 ! ka (i/o) = the number of hebden iterations (so far) taken to deter-
2281 ! mine step. ka .lt. 0 on input means this is the first
2282 ! attempt to determine step (for the present dig and dihdi)
2283 ! -- ka is initialized to 0 in this case. output with
2284 ! ka = 0 (or v(stppar) = 0) means step = -(h**-1)*g.
2285 ! l (i/o) = workspace of length p*(p+1)/2 for cholesky factors.
2286 ! p (in) = number of parameters -- the hessian is a p x p matrix.
2287 ! step (i/o) = the step computed.
2288 ! v (i/o) contains various constants and variables described below.
2289 ! w (i/o) = workspace of length 4*p + 6.
2291 ! *** entries in v ***
2293 ! v(dgnorm) (i/o) = 2-norm of (d**-1)*g.
2294 ! v(dstnrm) (output) = 2-norm of d*step.
2295 ! v(dst0) (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or
2296 ! overestimate of smallest eigenvalue of (d**-1)*h*(d**-1).
2297 ! v(epslon) (in) = max. rel. error allowed for psi(step). for the
2298 ! step returned, psi(step) will exceed its optimal value
2299 ! by less than -v(epslon)*psi(step). suggested value = 0.1.
2300 ! v(gtstep) (out) = inner product between g and step.
2301 ! v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step) (for pos. def.
2302 ! h only -- v(nreduc) is set to zero otherwise).
2303 ! v(phmnfc) (in) = tol. (together with v(phmxfc)) for accepting step
2304 ! (more*s sigma). the error v(dstnrm) - v(radius) must lie
2305 ! between v(phmnfc)*v(radius) and v(phmxfc)*v(radius).
2306 ! v(phmxfc) (in) (see v(phmnfc).)
2307 ! suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5.
2308 ! v(preduc) (out) = psi(step) = predicted obj. func. reduction for step.
2309 ! v(radius) (in) = radius of current (scaled) trust region.
2310 ! v(rad0) (i/o) = value of v(radius) from previous call.
2311 ! v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha
2312 ! described below under algorithm notes. if h + alpha*d**2
2313 ! (see algorithm notes) is (nearly) singular, however,
2314 ! then v(stppar) = -alpha.
2316 ! *** usage notes ***
2318 ! if it is desired to recompute step using a different value of
2319 ! v(radius), then this routine may be restarted by calling it
2320 ! with all parameters unchanged except v(radius). (this explains
2321 ! why step and w are listed as i/o). on an initial call (one with
2322 ! ka .lt. 0), step and w need not be initialized and only compo-
2323 ! nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and
2324 ! v(rad0) of v must be initialized.
2326 ! *** algorithm notes ***
2328 ! the desired g-q-t step (ref. 2, 3, 4, 6) satisfies
2329 ! (h + alpha*d**2)*step = -g for some nonnegative alpha such that
2330 ! h + alpha*d**2 is positive semidefinite. alpha and step are
2331 ! computed by a scheme analogous to the one described in ref. 5.
2332 ! estimates of the smallest and largest eigenvalues of the hessian
2333 ! are obtained from the gerschgorin circle theorem enhanced by a
2334 ! simple form of the scaling described in ref. 7. cases in which
2335 ! h + alpha*d**2 is nearly (or exactly) singular are handled by
2336 ! the technique discussed in ref. 2. in these cases, a step of
2337 ! (exact) length v(radius) is returned for which psi(step) exceeds
2338 ! its optimal value by less than -v(epslon)*psi(step). the test
2339 ! suggested in ref. 6 for detecting the special case is performed
2340 ! once two matrix factorizations have been done -- doing so sooner
2341 ! seems to degrade the performance of optimization routines that
2342 ! call this routine.
2344 ! *** functions and subroutines called ***
2346 ! dotprd - returns inner product of two vectors.
2347 ! litvmu - applies inverse-transpose of compact lower triang. matrix.
2348 ! livmul - applies inverse of compact lower triang. matrix.
2349 ! lsqrt - finds cholesky factor (of compactly stored lower triang.).
2350 ! lsvmin - returns approx. to min. sing. value of lower triang. matrix.
2351 ! rmdcon - returns machine-dependent constants.
2352 ! v2norm - returns 2-norm of a vector.
2354 ! *** references ***
2356 ! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive
2357 ! nonlinear least-squares algorithm, acm trans. math.
2358 ! software, vol. 7, no. 3.
2359 ! 2. gay, d.m. (1981), computing optimal locally constrained steps,
2360 ! siam j. sci. statist. computing, vol. 2, no. 2, pp.
2362 ! 3. goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966),
2363 ! maximization by quadratic hill-climbing, econometrica 34,
2365 ! 4. hebden, m.d. (1973), an algorithm for minimization using exact
2366 ! second derivatives, report t.p. 515, theoretical physics
2367 ! div., a.e.r.e. harwell, oxon., england.
2368 ! 5. more, j.j. (1978), the levenberg-marquardt algorithm, implemen-
2369 ! tation and theory, pp.105-116 of springer lecture notes
2370 ! in mathematics no. 630, edited by g.a. watson, springer-
2371 ! verlag, berlin and new york.
2372 ! 6. more, j.j., and sorensen, d.c. (1981), computing a trust region
2373 ! step, technical report anl-81-83, argonne national lab.
2374 ! 7. varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15,
2379 ! coded by david m. gay.
2380 ! this subroutine was written in connection with research
2381 ! supported by the national science foundation under grants
2382 ! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
2385 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2387 ! *** local variables ***
2390 integer :: dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,&
2391 j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x
2392 real(kind=8) :: alphak, aki, akk, delta, dst, eps, gtsta, lk,&
2393 oldphi, phi, phimax, phimin, psifac, rad, radsq,&
2394 root, si, sk, sw, t, twopsi, t1, t2, uk, wi
2397 real(kind=8) :: big, dgxfac !el, epsfac, four, half, kappa, negone,
2398 !el 1 one, p001, six, three, two, zero
2400 ! *** intrinsic functions ***
2402 !el real(kind=8) :: dabs, dmax1, dmin1, dsqrt
2404 ! *** external functions and subroutines ***
2406 !el external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm
2407 !el real(kind=8) :: dotprd, lsvmin, rmdcon, v2norm
2409 ! *** subscripts for v ***
2411 !el integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc,
2412 !el 1 phmnfc, phmxfc, preduc, radius, rad0
2414 ! data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/,
2415 ! 1 nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/,
2416 ! 2 rad0/9/, stppar/5/
2418 integer,parameter :: dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,&
2419 nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,&
2424 ! data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/,
2425 ! 1 kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/,
2426 ! 2 six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/
2428 real(kind=8), parameter :: epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,&
2429 kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,&
2430 six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0
2433 data big/0.d+0/, dgxfac/0.d+0/
2437 ! *** store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx).
2439 ! *** store gerschgorin over- and underestimates of the largest
2440 ! *** and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax)
2441 ! *** and w(emin) respectively.
2444 ! *** for use in recomputing step, the final values of lk, uk, dst,
2445 ! *** and the inverse derivative of more*s phi at 0 (for pos. def.
2446 ! *** h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin)
2452 ! *** store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p).
2455 ! *** store -d*step in w(q),...,w(q0+p).
2458 ! *** allocate storage for scratch vector x ***
2462 ! *** phitol = max. error allowed in dst = v(dstnrm) = 2-norm of
2464 phimax = v(phmxfc) * rad
2465 phimin = v(phmnfc) * rad
2466 psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) * &
2467 (kappa + one) + kappa + two) * rad**2)
2468 ! *** oldphi is used to detect limits of numerical accuracy. if
2469 ! *** we recompute step and it does not change, then we accept it.
2476 ! *** start or restart, depending on ka ***
2478 if (ka .ge. 0) go to 290
2480 ! *** fresh start ***
2486 v(dgnorm) = v2norm(p, dig)
2490 if (v(dgnorm) .eq. zero) kamin = 0
2492 ! *** store diag(dihdi) in w(diag0+1),...,w(diag0+p) ***
2501 ! *** determine w(dggdmx), the largest element of dihdi ***
2507 if (t1 .lt. t) t1 = t
2511 ! *** try alpha = 0 ***
2513 30 call lsqrt(1, p, l, dihdi, irc)
2514 if (irc .eq. 0) go to 50
2515 ! *** indef. h -- underestimate smallest eigenvalue, use this
2516 ! *** estimate to initialize lower bound lk on alpha.
2523 call litvmu(irc, w, l, w)
2527 if (restrt) go to 210
2530 ! *** positive definite h -- compute unmodified newton step. ***
2532 t = lsvmin(p, l, w(q), w(q))
2533 if (t .ge. one) go to 60
2534 if (big .le. zero) big = rmdcon(6)
2535 if (v(dgnorm) .ge. t*t*big) go to 70
2536 60 call livmul(p, w(q), l, dig)
2537 gtsta = dotprd(p, w(q), w(q))
2538 v(nreduc) = half * gtsta
2539 call litvmu(p, w(q), l, w(q))
2540 dst = v2norm(p, w(q))
2543 if (phi .le. phimax) go to 260
2544 if (restrt) go to 210
2546 ! *** prepare to compute gerschgorin estimates of largest (and
2547 ! *** smallest) eigenvalues. ***
2552 if (i .eq. 1) go to 90
2564 ! *** (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1) ***
2568 if (p .le. 1) go to 120
2572 if (t .ge. t1) go to 110
2584 if (i .eq. k) go to 130
2585 aki = dabs(dihdi(k1))
2588 t1 = half * (akk - w(j) + si - aki)
2589 t1 = t1 + dsqrt(t1*t1 + sk*aki)
2590 if (t .lt. t1) t = t1
2591 if (i .lt. k) go to 140
2597 uk = v(dgnorm)/rad - w(emin)
2598 if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2599 if (uk .le. zero) uk = p001
2601 ! *** compute gerschgorin (over-)estimate of largest eigenvalue ***
2605 if (p .le. 1) go to 170
2609 if (t .le. t1) go to 160
2621 if (i .eq. k) go to 180
2622 aki = dabs(dihdi(k1))
2625 t1 = half * (w(j) + si - aki - akk)
2626 t1 = t1 + dsqrt(t1*t1 + sk*aki)
2627 if (t .lt. t1) t = t1
2628 if (i .lt. k) go to 190
2634 lk = dmax1(lk, v(dgnorm)/rad - w(emax))
2636 ! *** alphak = current value of alpha (see alg. notes above). we
2637 ! *** use more*s scheme for initializing it.
2638 alphak = dabs(v(stppar)) * v(rad0)/rad
2640 if (irc .ne. 0) go to 210
2642 ! *** compute l0 for positive definite h ***
2644 call livmul(p, w, l, w(q))
2646 w(phipin) = dst / t / t
2647 lk = dmax1(lk, phi*w(phipin))
2649 ! *** safeguard alphak and add alphak*i to (d**-1)*h*(d**-1) ***
2652 if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk) &
2653 alphak = uk * dmax1(p001, dsqrt(lk/uk))
2654 if (alphak .le. zero) alphak = half * uk
2655 if (alphak .le. zero) alphak = uk
2660 dihdi(k) = w(j) + alphak
2663 ! *** try computing cholesky decomposition ***
2665 call lsqrt(1, p, l, dihdi, irc)
2666 if (irc .eq. 0) go to 240
2668 ! *** (d**-1)*h*(d**-1) + alphak*i is indefinite -- overestimate
2669 ! *** smallest eigenvalue for use in updating lk ***
2677 call litvmu(irc, w, l, w)
2679 lk = alphak - t/t1/t1
2683 ! *** alphak makes (d**-1)*h*(d**-1) positive definite.
2684 ! *** compute q = -d*step, check for convergence. ***
2686 240 call livmul(p, w(q), l, dig)
2687 gtsta = dotprd(p, w(q), w(q))
2688 call litvmu(p, w(q), l, w(q))
2689 dst = v2norm(p, w(q))
2691 if (phi .le. phimax .and. phi .ge. phimin) go to 270
2692 if (phi .eq. oldphi) go to 270
2694 if (phi .lt. zero) go to 330
2696 ! *** unacceptable alphak -- update lk, uk, alphak ***
2698 250 if (ka .ge. kalim) go to 270
2699 ! *** the following dmin1 is necessary because of restarts ***
2700 if (phi .lt. zero) uk = dmin1(uk, alphak)
2701 ! *** kamin = 0 only iff the gradient vanishes ***
2702 if (kamin .eq. 0) go to 210
2703 call livmul(p, w, l, w(q))
2705 alphak = alphak + (phi/t1) * (dst/t1) * (dst/rad)
2706 lk = dmax1(lk, alphak)
2709 ! *** acceptable step on first try ***
2713 ! *** successful step in general. compute step = -(d**-1)*q ***
2717 step(i) = -w(j)/d(i)
2720 v(preduc) = half * (dabs(alphak)*dst*dst + gtsta)
2724 ! *** restart with new radius ***
2726 290 if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310
2728 ! *** prepare to return newton step ***
2742 if (v(dgnorm) .eq. zero) kamin = 0
2743 if (ka .eq. 0) go to 50
2746 alphak = dabs(v(stppar))
2750 if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2751 if (uk .le. zero) uk = p001
2752 if (rad .gt. v(rad0)) go to 320
2754 ! *** smaller radius ***
2756 if (alphak .gt. zero) lk = w(lk0)
2757 lk = dmax1(lk, t - w(emax))
2758 if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2761 ! *** bigger radius ***
2762 320 if (alphak .gt. zero) uk = dmin1(uk, w(uk0))
2763 lk = dmax1(zero, -v(dst0), t - w(emax))
2764 if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2767 ! *** decide whether to check for special case... in practice (from
2768 ! *** the standpoint of the calling optimization code) it seems best
2769 ! *** not to check until a few iterations have failed -- hence the
2770 ! *** test on kamin below.
2772 330 delta = alphak + dmin1(zero, v(dst0))
2773 twopsi = alphak*dst*dst + gtsta
2774 if (ka .ge. kamin) go to 340
2775 ! *** if the test in ref. 2 is satisfied, fall through to handle
2776 ! *** the special case (as soon as the more-sorensen test detects
2778 if (delta .ge. psifac*twopsi) go to 370
2780 ! *** check for the special case of h + alpha*d**2 (nearly)
2781 ! *** singular. use one step of inverse power method with start
2782 ! *** from lsvmin to obtain approximate eigenvector corresponding
2783 ! *** to smallest eigenvalue of (d**-1)*h*(d**-1). lsvmin returns
2784 ! *** x and w with l*w = x.
2786 340 t = lsvmin(p, l, w(x), w)
2788 ! *** normalize w ***
2791 ! *** complete current inv. power iter. -- replace w by (l**-t)*w.
2792 call litvmu(p, w, l, w)
2793 t2 = one/v2norm(p, w)
2798 ! *** now w is the desired approximate (unit) eigenvector and
2799 ! *** t*x = ((d**-1)*h*(d**-1) + alphak*i)*w.
2801 sw = dotprd(p, w(q), w)
2802 t1 = (rad + dst) * (rad - dst)
2803 root = dsqrt(sw*sw + t1)
2804 if (sw .lt. zero) root = -root
2805 si = t1 / (sw + root)
2807 ! *** the actual test for the special case...
2809 if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380
2811 ! *** update upper bound on smallest eigenvalue (when not positive)
2812 ! *** (as recommended by more and sorensen) and continue...
2814 if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak)
2815 lk = dmax1(lk, -v(dst0))
2817 ! *** check whether we can hope to detect the special case in
2818 ! *** the available arithmetic. accept step as it is if not.
2820 ! *** if not yet available, obtain machine dependent value dgxfac.
2821 370 if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3)
2823 if (delta .gt. dgxfac*w(dggdmx)) go to 250
2826 ! *** special case detected... negate alphak to indicate special case
2828 380 alphak = -alphak
2829 v(preduc) = half * twopsi
2831 ! *** accept current step if adding si*w would lead to a
2832 ! *** further relative reduction in psi of less than v(epslon)/3.
2835 t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w)))
2836 if (t .lt. eps*twopsi/six) go to 390
2837 v(preduc) = v(preduc) + t
2842 w(j) = t1*w(i) - w(j)
2843 step(i) = w(j) / d(i)
2845 v(gtstep) = dotprd(p, dig, w(q))
2847 ! *** save values for use in a possible restart ***
2856 ! *** restore diagonal of dihdi ***
2867 ! *** last card of gqtst follows ***
2868 end subroutine gqtst
2869 !-----------------------------------------------------------------------------
2870 subroutine lsqrt(n1, n, l, a, irc)
2872 ! *** compute rows n1 through n of the cholesky factor l of
2873 ! *** a = l*(l**t), where l and the lower triangle of a are both
2874 ! *** stored compactly by rows (and may occupy the same storage).
2875 ! *** irc = 0 means all went well. irc = j means the leading
2876 ! *** principal j x j submatrix of a is not positive definite --
2877 ! *** and l(j*(j+1)/2) contains the (nonpos.) reduced j-th diagonal.
2879 ! *** parameters ***
2881 integer :: n1, n, irc
2882 !al real(kind=8) :: l(1), a(1)
2883 real(kind=8) :: l(n*(n+1)/2), a(n*(n+1)/2)
2884 ! dimension l(n*(n+1)/2), a(n*(n+1)/2)
2886 ! *** local variables ***
2888 integer :: i, ij, ik, im1, i0, j, jk, jm1, j0, k
2889 real(kind=8) :: t, td !el, zero
2891 ! *** intrinsic functions ***
2893 !el real(kind=8) :: dsqrt
2898 real(kind=8),parameter :: zero=0.d+0
2903 i0 = n1 * (n1 - 1) / 2
2906 if (i .eq. 1) go to 40
2911 if (j .eq. 1) go to 20
2920 t = (a(ij) - t) / l(j0)
2926 if (t .le. zero) go to 60
2938 ! *** last card of lsqrt ***
2939 end subroutine lsqrt
2940 !-----------------------------------------------------------------------------
2941 real(kind=8) function lsvmin(p, l, x, y)
2943 ! *** estimate smallest sing. value of packed lower triang. matrix l
2945 ! *** parameter declarations ***
2948 !al real(kind=8) :: l(1), x(p), y(p)
2949 real(kind=8) :: l(p*(p+1)/2), x(p), y(p)
2950 ! dimension l(p*(p+1)/2)
2952 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2956 ! this function returns a good over-estimate of the smallest
2957 ! singular value of the packed lower triangular matrix l.
2959 ! *** parameter description ***
2961 ! p (in) = the order of l. l is a p x p lower triangular matrix.
2962 ! l (in) = array holding the elements of l in row order, i.e.
2963 ! l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc.
2964 ! x (out) if lsvmin returns a positive value, then x is a normalized
2965 ! approximate left singular vector corresponding to the
2966 ! smallest singular value. this approximation may be very
2967 ! crude. if lsvmin returns zero, then some components of x
2968 ! are zero and the rest retain their input values.
2969 ! y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an
2970 ! unnormalized approximate right singular vector correspond-
2971 ! ing to the smallest singular value. this approximation
2972 ! may be crude. if lsvmin returns zero, then y retains its
2973 ! input value. the caller may pass the same vector for x
2974 ! and y (nonstandard fortran usage), in which case y over-
2975 ! writes x (for nonzero lsvmin returns).
2977 ! *** algorithm notes ***
2979 ! the algorithm is based on (1), with the additional provision that
2980 ! lsvmin = 0 is returned if the smallest diagonal element of l
2981 ! (in magnitude) is not more than the unit roundoff times the
2982 ! largest. the algorithm uses a random number generator proposed
2983 ! in (4), which passes the spectral test with flying colors -- see
2986 ! *** subroutines and functions called ***
2988 ! v2norm - function, returns the 2-norm of a vector.
2990 ! *** references ***
2992 ! (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977),
2993 ! an estimate for the condition number of a matrix, report
2994 ! tm-310, applied math. div., argonne national laboratory.
2996 ! (2) hoaglin, d.c. (1976), theoretical properties of congruential
2997 ! random-number generators -- an empirical view,
2998 ! memorandum ns-340, dept. of statistics, harvard univ.
3000 ! (3) knuth, d.e. (1969), the art of computer programming, vol. 2
3001 ! (seminumerical algorithms), addison-wesley, reading, mass.
3003 ! (4) smith, c.s. (1971), multiplicative pseudo-random number
3004 ! generators with prime modulus, j. assoc. comput. mach. 18,
3009 ! designed and coded by david m. gay (winter 1977/summer 1978).
3013 ! this subroutine was written in connection with research
3014 ! supported by the national science foundation under grants
3015 ! mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989.
3017 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3019 ! *** local variables ***
3021 integer :: i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1
3022 real(kind=8) :: b, sminus, splus, t, xminus, xplus
3026 !el real(kind=8) :: half, one, r9973, zero
3028 ! *** intrinsic functions ***
3032 !el real(kind=8) :: dabs
3034 ! *** external functions and subroutines ***
3036 !el external dotprd, v2norm, vaxpy
3037 !el real(kind=8) :: dotprd, v2norm
3040 ! data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/
3042 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0
3050 ! *** first check whether to return lsvmin = 0 and initialize x ***
3055 if (l(jj) .eq. zero) go to 110
3056 ix = mod(3432*ix, 9973)
3057 b = half*(one + float(ix)/r9973)
3060 if (p .le. 1) go to 60
3063 if (l(ii) .eq. zero) go to 110
3065 x(i) = xplus * l(ji)
3068 ! *** solve (l**t)*x = b, where the components of b have randomly
3069 ! *** chosen magnitudes in (.5,1) with signs chosen to make x large.
3071 ! do j = p-1 to 1 by -1...
3074 ! *** determine x(j) in this iteration. note for i = 1,2,...,j
3075 ! *** that x(i) holds the current partial sum for row i.
3076 ix = mod(3432*ix, 9973)
3077 b = half*(one + float(ix)/r9973)
3079 xminus = (-b - x(j))
3081 sminus = dabs(xminus)
3086 xminus = xminus/l(jj)
3087 if (jm1 .eq. 0) go to 30
3090 splus = splus + dabs(x(i) + l(ji)*xplus)
3091 sminus = sminus + dabs(x(i) + l(ji)*xminus)
3093 30 if (sminus .gt. splus) xplus = xminus
3095 ! *** update partial sums ***
3096 if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x)
3099 ! *** normalize x ***
3101 60 t = one/v2norm(p, x)
3105 ! *** solve l*y = x and return lsvmin = 1/twonorm(y) ***
3112 if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y)
3113 y(j) = (x(j) - t) / l(jj)
3116 lsvmin = one/v2norm(p, y)
3121 ! *** last card of lsvmin follows ***
3123 !-----------------------------------------------------------------------------
3124 subroutine slvmul(p, y, s, x)
3126 ! *** set y = s * x, s = p x p symmetric matrix. ***
3127 ! *** lower triangle of s stored rowwise. ***
3129 ! *** parameter declarations ***
3132 !al real(kind=8) :: s(1), x(p), y(p)
3133 real(kind=8) :: s(p*(p+1)/2), x(p), y(p)
3134 ! dimension s(p*(p+1)/2)
3136 ! *** local variables ***
3138 integer :: i, im1, j, k
3141 ! *** no intrinsic functions ***
3143 ! *** external function ***
3146 !el real(kind=8) :: dotprd
3148 !-----------------------------------------------------------------------
3152 y(i) = dotprd(i, s(j), x)
3156 if (p .le. 1) go to 999
3163 y(k) = y(k) + s(j)*xi
3169 ! *** last card of slvmul follows ***
3170 end subroutine slvmul
3171 !-----------------------------------------------------------------------------
3173 !-----------------------------------------------------------------------------
3174 subroutine minimize(etot,x,iretcode,nfun)
3176 use energy, only: func,gradient,fdum!,etotal,enerprint
3178 ! implicit real*8 (a-h,o-z)
3179 ! include 'DIMENSIONS'
3180 integer,parameter :: liv=60
3181 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3182 !********************************************************************
3183 ! OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
3184 ! the calling subprogram. *
3185 ! when d(i)=1.0, then v(35) is the length of the initial step, *
3186 ! calculated in the usual pythagorean way. *
3187 ! absolute convergence occurs when the function is within v(31) of *
3188 ! zero. unless you know the minimum value in advance, abs convg *
3189 ! is probably not useful. *
3190 ! relative convergence is when the model predicts that the function *
3191 ! will decrease by less than v(32)*abs(fun). *
3192 !********************************************************************
3193 ! include 'COMMON.IOUNITS'
3194 ! include 'COMMON.VAR'
3195 ! include 'COMMON.GEO'
3196 ! include 'COMMON.MINIM'
3198 !el common /srutu/ icall
3199 integer,dimension(liv) :: iv
3200 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
3201 !el real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3202 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3203 real(kind=8) :: energia(0:n_ene)
3204 ! external func,gradient,fdum
3205 ! external func_restr,grad_restr
3206 logical :: not_done,change,reduce
3207 !el common /przechowalnia/ v
3209 integer :: iretcode,nfun,lv,nvar_restr,idum(1),j
3210 real(kind=8) :: etot,rdum(1) !,fdum
3212 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3214 if (.not.allocated(v)) allocate(v(1:lv))
3220 ! DO WHILE (NOT_DONE)
3222 call deflt(2,iv,liv,lv,v)
3223 ! 12 means fresh start, dont call deflt
3225 ! max num of fun calls
3226 if (maxfun.eq.0) maxfun=500
3228 ! max num of iterations
3229 if (maxmin.eq.0) maxmin=1000
3233 ! selects output unit
3235 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
3236 ! 1 means to print out result
3237 iv(22)=print_min_res
3238 ! 1 means to print out summary stats
3239 iv(23)=print_min_stat
3240 ! 1 means to print initial x and d
3241 iv(24)=print_min_ini
3242 ! min val for v(radfac) default is 0.1
3244 ! max val for v(radfac) default is 4.0
3247 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3248 ! the sumsl default is 0.1
3250 ! false conv if (act fnctn decrease) .lt. v(34)
3251 ! the sumsl default is 100*machep
3253 ! absolute convergence
3254 if (tolf.eq.0.0D0) tolf=1.0D-4
3256 ! relative convergence
3257 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3259 ! controls initial step size
3261 ! large vals of d correspond to small components of step
3268 !d print *,'Calling SUMSL'
3269 ! call var_to_geom(nvar,x)
3271 ! call etotal(energia(0))
3275 call x2xx(x,xx,nvar_restr)
3276 call sumsl(nvar_restr,d,xx,func_restr,grad_restr,&
3277 iv,liv,lv,v,idum,rdum,fdum)
3280 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
3284 !d print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
3285 !d write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
3288 call var_to_geom(nvar,x)
3290 ! write (iout,'(a)') 'Reduction worked, minimizing again...'
3296 !el---------------------
3297 ! write (iout,'(/a)') &
3298 ! "Cartesian coordinates of the reference structure after SUMSL"
3299 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3300 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3302 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
3303 ! restyp(itype(i)),i,(c(j,i),j=1,3),&
3304 ! (c(j,i+nres),j=1,3)
3306 !el----------------------------
3307 ! call etotal(energia) !sp
3309 ! call enerprint(energia) !sp
3312 ! write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
3317 end subroutine minimize
3318 !-----------------------------------------------------------------------------
3320 !-----------------------------------------------------------------------------
3321 subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
3323 use energy, only: cartder,zerograd,etotal,sum_gradient
3324 ! implicit real*8 (a-h,o-z)
3325 ! include 'DIMENSIONS'
3326 ! include 'COMMON.CHAIN'
3327 ! include 'COMMON.DERIV'
3328 ! include 'COMMON.VAR'
3329 ! include 'COMMON.INTERACT'
3330 ! include 'COMMON.FFIELD'
3331 ! include 'COMMON.IOUNITS'
3333 integer :: uiparm(1)
3334 real(kind=8) :: urparm(1)
3335 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
3336 integer :: n,nf,ig,ind,i,j,ij,k,igall
3337 real(kind=8) :: f,gphii,gthetai,galphai,gomegai
3338 real(kind=8),external :: ufparm
3341 if (nf-nfl+1) 20,30,40
3342 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
3343 ! write (iout,*) 'grad 20'
3348 ! Intercept NaNs in the coordinates
3349 ! write(iout,*) (var(i),i=1,nvar)
3354 if (x_sum.ne.x_sum) then
3355 write(iout,*)" *** grad_restr : Found NaN in coordinates"
3357 print *," *** grad_restr : Found NaN in coordinates"
3361 call var_to_geom_restr(n,x)
3364 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
3368 ! Convert the Cartesian gradient into internal-coordinate gradient.
3374 IF (mask_phi(i+2).eq.1) THEN
3379 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
3380 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
3393 IF (mask_theta(i+2).eq.1) THEN
3399 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
3400 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
3410 if (itype(i).ne.10) then
3411 IF (mask_side(i).eq.1) THEN
3415 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
3424 if (itype(i).ne.10) then
3425 IF (mask_side(i).eq.1) THEN
3429 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
3437 ! Add the components corresponding to local energy terms.
3444 if (mask_phi(i).eq.1) then
3446 g(ig)=g(ig)+gloc(igall,icg)
3452 if (mask_theta(i).eq.1) then
3454 g(ig)=g(ig)+gloc(igall,icg)
3460 if (itype(i).ne.10) then
3462 if (mask_side(i).eq.1) then
3464 g(ig)=g(ig)+gloc(igall,icg)
3471 !d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
3474 end subroutine grad_restr
3475 !-----------------------------------------------------------------------------
3476 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F
3479 use energy, only: zerograd,etotal,sum_gradient
3480 ! implicit real*8 (a-h,o-z)
3481 ! include 'DIMENSIONS'
3482 ! include 'COMMON.DERIV'
3483 ! include 'COMMON.IOUNITS'
3484 ! include 'COMMON.GEO'
3487 !el common /chuju/ jjj
3488 real(kind=8) :: energia(0:n_ene)
3490 real(kind=8),external :: ufparm
3491 integer :: uiparm(1)
3492 real(kind=8) :: urparm(1)
3493 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3494 ! if (jjj.gt.0) then
3495 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
3499 call var_to_geom_restr(n,x)
3502 !d write (iout,*) 'ETOTAL called from FUNC'
3503 call etotal(energia)
3506 ! if (jjj.gt.0) then
3507 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
3508 ! write (iout,*) 'f=',etot
3512 end subroutine func_restr
3513 !-----------------------------------------------------------------------------
3514 ! subroutine func(n,x,nf,f,uiparm,urparm,ufparm) in module energy
3515 !-----------------------------------------------------------------------------
3516 subroutine x2xx(x,xx,n)
3518 ! implicit real*8 (a-h,o-z)
3519 ! include 'DIMENSIONS'
3520 ! include 'COMMON.VAR'
3521 ! include 'COMMON.CHAIN'
3522 ! include 'COMMON.INTERACT'
3523 integer :: n,i,ij,ig,igall
3524 real(kind=8),dimension(6*nres) :: xx,x !(maxvar) (maxvar=6*maxres)
3534 if (mask_phi(i).eq.1) then
3542 if (mask_theta(i).eq.1) then
3550 if (itype(i).ne.10) then
3552 if (mask_side(i).eq.1) then
3564 !-----------------------------------------------------------------------------
3565 !el subroutine xx2x(x,xx) in module math
3566 !-----------------------------------------------------------------------------
3567 subroutine minim_dc(etot,iretcode,nfun)
3570 use energy, only: fdum,check_ecartint
3571 ! implicit real*8 (a-h,o-z)
3572 ! include 'DIMENSIONS'
3576 integer,parameter :: liv=60
3577 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3578 ! include 'COMMON.SETUP'
3579 ! include 'COMMON.IOUNITS'
3580 ! include 'COMMON.VAR'
3581 ! include 'COMMON.GEO'
3582 ! include 'COMMON.MINIM'
3583 ! include 'COMMON.CHAIN'
3584 integer :: iretcode,nfun,k,i,j,lv,idum(1)
3585 integer,dimension(liv) :: iv
3586 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
3587 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3588 !el common /przechowalnia/ v
3590 real(kind=8) :: energia(0:n_ene)
3591 ! external func_dc,grad_dc ,fdum
3592 logical :: not_done,change,reduce
3593 real(kind=8) :: g(6*nres),f1,etot,rdum(1) !,fdum
3595 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3597 if (.not. allocated(v)) allocate(v(1:lv))
3599 call deflt(2,iv,liv,lv,v)
3600 ! 12 means fresh start, dont call deflt
3602 ! max num of fun calls
3603 if (maxfun.eq.0) maxfun=500
3605 ! max num of iterations
3606 if (maxmin.eq.0) maxmin=1000
3610 ! selects output unit
3612 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
3613 ! 1 means to print out result
3614 iv(22)=print_min_res
3615 ! 1 means to print out summary stats
3616 iv(23)=print_min_stat
3617 ! 1 means to print initial x and d
3618 iv(24)=print_min_ini
3619 ! min val for v(radfac) default is 0.1
3621 ! max val for v(radfac) default is 4.0
3624 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3625 ! the sumsl default is 0.1
3627 ! false conv if (act fnctn decrease) .lt. v(34)
3628 ! the sumsl default is 100*machep
3630 ! absolute convergence
3631 if (tolf.eq.0.0D0) tolf=1.0D-4
3633 ! relative convergence
3634 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3636 ! controls initial step size
3638 ! large vals of d correspond to small components of step
3651 if (ialph(i,1).gt.0) then
3659 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
3669 if (ialph(i,1).gt.0) then
3676 call chainbuild_cart
3680 !d call func_dc(k,x,nf,f,idum,rdum,fdum)
3681 !d call grad_dc(k,x,nf,g,idum,rdum,fdum)
3685 !d call func_dc(k,x,nf,f1,idum,rdum,fdum)
3687 !d print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
3689 !el---------------------
3690 ! write (iout,'(/a)') &
3691 ! "Cartesian coordinates of the reference structure after SUMSL"
3692 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3693 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3695 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
3696 ! restyp(itype(i)),i,(c(j,i),j=1,3),&
3697 ! (c(j,i+nres),j=1,3)
3699 !el----------------------------
3704 end subroutine minim_dc
3705 !-----------------------------------------------------------------------------
3706 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
3709 use energy, only: zerograd,etotal
3710 ! implicit real*8 (a-h,o-z)
3711 ! include 'DIMENSIONS'
3715 ! include 'COMMON.SETUP'
3716 ! include 'COMMON.DERIV'
3717 ! include 'COMMON.IOUNITS'
3718 ! include 'COMMON.GEO'
3719 ! include 'COMMON.CHAIN'
3720 ! include 'COMMON.VAR'
3721 integer :: n,nf,k,i,j
3722 real(kind=8) :: energia(0:n_ene)
3723 real(kind=8),external :: ufparm
3724 integer :: uiparm(1)
3725 real(kind=8) :: urparm(1)
3726 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3729 !bad icg=mod(nf,2)+1
3740 if (ialph(i,1).gt.0) then
3747 call chainbuild_cart
3750 call etotal(energia)
3753 !d print *,'func_dc ',nf,nfl,f
3756 end subroutine func_dc
3757 !-----------------------------------------------------------------------------
3758 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
3761 use energy, only: cartgrad,zerograd,etotal
3762 ! implicit real*8 (a-h,o-z)
3763 ! include 'DIMENSIONS'
3767 ! include 'COMMON.SETUP'
3768 ! include 'COMMON.CHAIN'
3769 ! include 'COMMON.DERIV'
3770 ! include 'COMMON.VAR'
3771 ! include 'COMMON.INTERACT'
3772 ! include 'COMMON.FFIELD'
3773 ! include 'COMMON.MD'
3774 ! include 'COMMON.IOUNITS'
3775 real(kind=8),external :: ufparm
3776 integer :: n,nf,i,j,k
3777 integer :: uiparm(1)
3778 real(kind=8) :: urparm(1)
3779 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
3782 !elwrite(iout,*) "jestesmy w grad dc"
3784 !bad icg=mod(nf,2)+1
3786 !d print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
3787 !elwrite(iout,*) "jestesmy w grad dc nf-nfl+1", nf-nfl+1
3788 if (nf-nfl+1) 20,30,40
3789 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
3803 if (ialph(i,1).gt.0) then
3810 !elwrite(iout,*) "jestesmy w grad dc"
3811 call chainbuild_cart
3814 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
3818 !elwrite(iout,*) "jestesmy w grad dc"
3827 if (ialph(i,1).gt.0) then
3834 !elwrite(iout,*) "jestesmy w grad dc"
3837 end subroutine grad_dc
3838 !-----------------------------------------------------------------------------
3840 !-----------------------------------------------------------------------------
3842 subroutine minim_mcmf
3846 use energy, only: func,gradient,fdum
3847 ! implicit real*8 (a-h,o-z)
3848 ! include 'DIMENSIONS'
3850 integer,parameter :: liv=60
3851 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3852 ! include 'COMMON.VAR'
3853 ! include 'COMMON.IOUNITS'
3854 ! include 'COMMON.MINIM'
3855 ! real(kind=8) :: fdum
3856 ! external func,gradient,fdum
3857 !el real(kind=4) :: ran1,ran2,ran3
3858 ! include 'COMMON.SETUP'
3859 ! include 'COMMON.GEO'
3860 ! include 'COMMON.CHAIN'
3861 ! include 'COMMON.FFIELD'
3862 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3863 real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg
3864 real(kind=8),dimension(6*nres) :: d,garbage !(maxvar) (maxvar=6*maxres)
3865 !el real(kind=8) :: v(1:77+(6*nres)*(6*nres+17)/2+1)
3866 integer,dimension(6) :: indx
3867 integer,dimension(liv) :: iv
3868 integer :: lv,idum(1),nf !
3869 real(kind=8) :: rdum(1)
3870 real(kind=8) :: przes(3),obrot(3,3),eee
3873 integer,dimension(MPI_STATUS_SIZE) :: muster
3875 integer :: ichuj,i,ierr
3876 real(kind=8) :: rad,ene0
3877 data rad /1.745329252d-2/
3878 !el common /przechowalnia/ v
3880 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3881 if (.not. allocated(v)) allocate(v(1:lv))
3886 call mpi_recv(indx,6,mpi_integer,king,idint,CG_COMM,&
3888 if (indx(1).eq.0) return
3889 ! print *, 'worker ',me,' received order ',n,ichuj
3890 call mpi_recv(var,nvar,mpi_double_precision,&
3891 king,idreal,CG_COMM,muster,ierr)
3892 call mpi_recv(ene0,1,mpi_double_precision,&
3893 king,idreal,CG_COMM,muster,ierr)
3894 ! print *, 'worker ',me,' var read '
3897 call deflt(2,iv,liv,lv,v)
3898 ! 12 means fresh start, dont call deflt
3900 ! max num of fun calls
3901 if (maxfun.eq.0) maxfun=500
3903 ! max num of iterations
3904 if (maxmin.eq.0) maxmin=1000
3908 ! selects output unit
3911 ! 1 means to print out result
3913 ! 1 means to print out summary stats
3915 ! 1 means to print initial x and d
3917 ! min val for v(radfac) default is 0.1
3919 ! max val for v(radfac) default is 4.0
3921 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3922 ! the sumsl default is 0.1
3924 ! false conv if (act fnctn decrease) .lt. v(34)
3925 ! the sumsl default is 100*machep
3927 ! absolute convergence
3928 if (tolf.eq.0.0D0) tolf=1.0D-4
3930 ! relative convergence
3931 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3933 ! controls initial step size
3935 ! large vals of d correspond to small components of step
3944 call func(nvar,var,nf,eee,idum,rdum,fdum)
3945 if(eee.gt.1.0d18) then
3946 ! print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
3947 ! print *,' energy before SUMSL =',eee
3948 ! print *,' aborting local minimization'
3955 call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
3956 ! find which conformation was returned from sumsl
3959 ! total # of ftn evaluations (for iwf=0, it includes all minimizations).
3964 call mpi_send(indx,6,mpi_integer,king,idint,CG_COMM,&
3966 ! print '(a5,i3,15f10.5)', 'ENEX0',indx(1),v(10)
3967 call mpi_send(var,nvar,mpi_double_precision,&
3968 king,idreal,CG_COMM,ierr)
3969 call mpi_send(eee,1,mpi_double_precision,king,idreal,&
3971 call mpi_send(ene0,1,mpi_double_precision,king,idreal,&
3975 end subroutine minim_mcmf
3977 !-----------------------------------------------------------------------------
3979 !-----------------------------------------------------------------------------
3980 ! algorithm 611, collected algorithms from acm.
3981 ! algorithm appeared in acm-trans. math. software, vol.9, no. 4,
3982 ! dec., 1983, p. 503-524.
3983 integer function imdcon(k)
3987 ! *** return integer machine-dependent constants ***
3989 ! *** k = 1 means return standard output unit number. ***
3990 ! *** k = 2 means return alternate output unit number. ***
3991 ! *** k = 3 means return input unit number. ***
3992 ! (note -- k = 2, 3 are used only by test programs.)
3994 ! +++ port version follows...
3998 ! data mdperm(1)/2/, mdperm(2)/4/, mdperm(3)/1/
3999 ! imdcon = i1mach(mdperm(k))
4000 ! +++ end of port version +++
4002 ! +++ non-port version follows...
4004 data mdcon(1)/6/, mdcon(2)/8/, mdcon(3)/5/
4006 ! +++ end of non-port version +++
4009 ! *** last card of imdcon follows ***
4011 !-----------------------------------------------------------------------------
4012 real(kind=8) function rmdcon(k)
4014 ! *** return machine dependent constants used by nl2sol ***
4016 ! +++ comments below contain data statements for various machines. +++
4017 ! +++ to convert to another machine, place a c in column 1 of the +++
4018 ! +++ data statement line(s) that correspond to the current machine +++
4019 ! +++ and remove the c from column 1 of the data statement line(s) +++
4020 ! +++ that correspond to the new machine. +++
4024 ! *** the constant returned depends on k...
4026 ! *** k = 1... smallest pos. eta such that -eta exists.
4027 ! *** k = 2... square root of eta.
4028 ! *** k = 3... unit roundoff = smallest pos. no. machep such
4029 ! *** that 1 + machep .gt. 1 .and. 1 - machep .lt. 1.
4030 ! *** k = 4... square root of machep.
4031 ! *** k = 5... square root of big (see k = 6).
4032 ! *** k = 6... largest machine no. big such that -big exists.
4034 real(kind=8) :: big, eta, machep
4035 integer :: bigi(4), etai(4), machei(4)
4037 !el real(kind=8) :: dsqrt
4039 equivalence (big,bigi(1)), (eta,etai(1)), (machep,machei(1))
4041 ! +++ ibm 360, ibm 370, or xerox +++
4043 ! data big/z7fffffffffffffff/, eta/z0010000000000000/,
4044 ! 1 machep/z3410000000000000/
4046 ! +++ data general +++
4048 ! data big/0.7237005577d+76/, eta/0.5397605347d-78/,
4049 ! 1 machep/2.22044605d-16/
4053 ! data big/1.7d+38/, eta/2.938735878d-39/, machep/2.775557562d-17/
4057 ! data big/1.157920892d+77/, eta/8.636168556d-78/,
4058 ! 1 machep/5.551115124d-17/
4062 ! data big/1.69d+38/, eta/5.9d-39/, machep/2.1680435d-19/
4066 ! data big/"377777100000000000000000/,
4067 ! 1 eta/"002400400000000000000000/,
4068 ! 2 machep/"104400000000000000000000/
4072 ! data big/o0777777777777777,o7777777777777777/,
4073 ! 1 eta/o1771000000000000,o7770000000000000/,
4074 ! 2 machep/o1451000000000000,o0000000000000000/
4076 ! +++ control data +++
4078 ! data big/37767777777777777777b,37167777777777777777b/,
4079 ! 1 eta/00014000000000000000b,00000000000000000000b/,
4080 ! 2 machep/15614000000000000000b,15010000000000000000b/
4084 ! data big/1.0d+9786/, eta/1.0d-9860/, machep/1.4210855d-14/
4088 ! data big/8.988d+307/, eta/1.2d-308/, machep/1.734723476d-18/
4092 data big/1.7d+38/, eta/2.939d-39/, machep/1.3877788d-17/
4096 ! data bigi(1)/577767777777777777777b/,
4097 ! 1 bigi(2)/000007777777777777776b/,
4098 ! 2 etai(1)/200004000000000000000b/,
4099 ! 3 etai(2)/000000000000000000000b/,
4100 ! 4 machei(1)/377224000000000000000b/,
4101 ! 5 machei(2)/000000000000000000000b/
4103 ! +++ port library -- requires more than just a data statement... +++
4106 ! double precision d1mach, zero
4107 ! data big/0.d+0/, eta/0.d+0/, machep/0.d+0/, zero/0.d+0/
4108 ! if (big .gt. zero) go to 1
4111 ! machep = d1mach(4)
4114 ! +++ end of port +++
4116 !------------------------------- body --------------------------------
4118 go to (10, 20, 30, 40, 50, 60), k
4123 20 rmdcon = dsqrt(256.d+0*eta)/16.d+0
4129 40 rmdcon = dsqrt(machep)
4132 50 rmdcon = dsqrt(big/256.d+0)*16.d+0
4138 ! *** last card of rmdcon follows ***
4140 !-----------------------------------------------------------------------------
4142 !-----------------------------------------------------------------------------
4143 subroutine sc_move(n_start,n_end,n_maxtry,e_drop,n_fun,etot)
4146 use random, only: iran_num
4147 use energy, only: esc
4148 ! Perform a quick search over side-chain arrangments (over
4149 ! residues n_start to n_end) for a given (frozen) CA trace
4150 ! Only side-chains are minimized (at most n_maxtry times each),
4152 ! Stops if energy drops by e_drop, otherwise tries all residues
4153 ! in the given range
4154 ! If there is an energy drop, full minimization may be useful
4155 ! n_start, n_end CAN be modified by this routine, but only if
4156 ! out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
4157 ! NOTE: this move should never increase the energy
4161 ! implicit real*8 (a-h,o-z)
4162 ! include 'DIMENSIONS'
4164 ! include 'COMMON.GEO'
4165 ! include 'COMMON.VAR'
4166 ! include 'COMMON.HEADER'
4167 ! include 'COMMON.IOUNITS'
4168 ! include 'COMMON.CHAIN'
4169 ! include 'COMMON.FFIELD'
4171 ! External functions
4172 !el integer iran_num
4173 !el external iran_num
4176 integer :: n_start,n_end,n_maxtry
4177 real(kind=8) :: e_drop
4181 real(kind=8) :: etot
4184 ! real(kind=8) :: energy(0:n_ene)
4185 real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1)
4186 real(kind=8) :: orig_e,cur_e
4187 integer :: n,n_steps,n_first,n_cur,n_tot !,i
4188 real(kind=8) :: orig_w(0:n_ene)
4189 real(kind=8) :: wtime
4191 !elwrite(iout,*) "in sc_move etot= ", etot
4192 ! Set non side-chain weights to zero (minimization is faster)
4193 ! NOTE: e(2) does not actually depend on the side-chain, only CA
4222 ! Make sure n_start, n_end are within proper range
4223 if (n_start.lt.2) n_start=2
4224 if (n_end.gt.nres-1) n_end=nres-1
4225 !rc if (n_start.lt.n_end) then
4226 if (n_start.gt.n_end) then
4231 ! Save the initial values of energy and coordinates
4233 !d call etotal(energy)
4234 !d write (iout,*) 'start sc ene',energy(0)
4235 !d call enerprint(energy(0))
4241 !rc cur_alph(i)=alph(i)
4242 !rc cur_omeg(i)=omeg(i)
4245 !t wtime=MPI_WTIME()
4246 ! Try (one by one) all specified residues, starting from a
4247 ! random position in sequence
4248 ! Stop early if the energy has decreased by at least e_drop
4249 n_tot=n_end-n_start+1
4250 n_first=iran_num(0,n_tot-1)
4253 !rc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
4254 do while (n.lt.n_tot)
4255 n_cur=n_start+mod(n_first+n,n_tot)
4256 call single_sc_move(n_cur,n_maxtry,e_drop,&
4258 !elwrite(iout,*) "after msingle sc_move etot= ", etot
4259 ! If a lower energy was found, update the current structure...
4260 !rc if (etot.lt.cur_e) then
4263 !rc cur_alph(i)=alph(i)
4264 !rc cur_omeg(i)=omeg(i)
4267 ! ...else revert to the previous one
4270 !rc alph(i)=cur_alph(i)
4271 !rc omeg(i)=cur_omeg(i)
4277 !d call etotal(energy)
4278 !d print *,'running',n,energy(0)
4282 !d call etotal(energy)
4283 !d write (iout,*) 'end sc ene',energy(0)
4285 ! Put the original weights back to calculate the full energy
4301 !t write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
4303 end subroutine sc_move
4304 !-----------------------------------------------------------------------------
4305 subroutine single_sc_move(res_pick,n_maxtry,e_drop,n_steps,n_fun,e_sc)
4307 ! Perturb one side-chain (res_pick) and minimize the
4308 ! neighbouring region, keeping all CA's and non-neighbouring
4310 ! Try until e_drop energy improvement is achieved, or n_maxtry
4311 ! attempts have been made
4312 ! At the start, e_sc should contain the side-chain-only energy(0)
4313 ! nsteps and nfun for this move are ADDED to n_steps and n_fun
4315 use energy, only: esc
4316 use geometry, only:dist
4318 ! implicit real*8 (a-h,o-z)
4319 ! include 'DIMENSIONS'
4320 ! include 'COMMON.VAR'
4321 ! include 'COMMON.INTERACT'
4322 ! include 'COMMON.CHAIN'
4323 ! include 'COMMON.MINIM'
4324 ! include 'COMMON.FFIELD'
4325 ! include 'COMMON.IOUNITS'
4327 ! External functions
4328 !el double precision dist
4332 integer :: res_pick,n_maxtry
4333 real(kind=8) :: e_drop
4335 ! Input/Output arguments
4336 integer :: n_steps,n_fun
4337 real(kind=8) :: e_sc
4342 integer :: nres_moved
4343 integer :: iretcode,loc_nfun,orig_maxfun,n_try
4344 real(kind=8) :: sc_dist,sc_dist_cutoff
4345 ! real(kind=8) :: energy_(0:n_ene)
4346 real(kind=8) :: evdw,escloc,orig_e,cur_e
4347 real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1)
4348 real(kind=8) :: var(6*nres) !(maxvar) (maxvar=6*maxres)
4350 real(kind=8) :: orig_theta(1:nres),orig_phi(1:nres),&
4351 orig_alph(1:nres),orig_omeg(1:nres)
4353 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4354 ! Define what is meant by "neighbouring side-chain"
4355 sc_dist_cutoff=5.0D0
4357 ! Don't do glycine or ends
4359 if (i.eq.10 .or. i.eq.ntyp1) return
4361 ! Freeze everything (later will relax only selected side-chains)
4369 ! Find the neighbours of the side-chain to move
4370 ! and save initial variables
4375 ! Don't do glycine (itype(j)==10)
4376 if (itype(i).ne.10) then
4377 sc_dist=dist(nres+i,nres+res_pick)
4379 sc_dist=sc_dist_cutoff
4381 if (sc_dist.lt.sc_dist_cutoff) then
4382 nres_moved=nres_moved+1
4392 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4393 !elwrite(iout,*) "in sinle wsc=",wsc,"evdw",evdw,"wscloc",wscloc,"escloc",escloc
4394 e_sc=wsc*evdw+wscloc*escloc
4395 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4396 !d call etotal(energy)
4397 !d print *,'new ',(energy(k),k=0,n_ene)
4402 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
4403 ! Move the selected residue (don't worry if it fails)
4404 call gen_side(iabs(itype(res_pick)),theta(res_pick+1),&
4405 alph(res_pick),omeg(res_pick),fail)
4407 ! Minimize the side-chains starting from the new arrangement
4408 call geom_to_var(nvar,var)
4413 !rc orig_theta(i)=theta(i)
4414 !rc orig_phi(i)=phi(i)
4415 !rc orig_alph(i)=alph(i)
4416 !rc orig_omeg(i)=omeg(i)
4419 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4420 call minimize_sc1(e_sc,var,iretcode,loc_nfun)
4422 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4423 !v write(*,'(2i3,2f12.5,2i3)')
4424 !v & res_pick,nres_moved,orig_e,e_sc-cur_e,
4425 !v & iretcode,loc_nfun
4427 !$$$ if (iretcode.eq.8) then
4428 !$$$ write(iout,*)'Coordinates just after code 8'
4429 !$$$ call chainbuild
4430 !$$$ call all_varout
4431 !$$$ call flush(iout)
4433 !$$$ theta(i)=orig_theta(i)
4434 !$$$ phi(i)=orig_phi(i)
4435 !$$$ alph(i)=orig_alph(i)
4436 !$$$ omeg(i)=orig_omeg(i)
4438 !$$$ write(iout,*)'Coordinates just before code 8'
4439 !$$$ call chainbuild
4440 !$$$ call all_varout
4441 !$$$ call flush(iout)
4444 n_fun=n_fun+loc_nfun
4446 call var_to_geom(nvar,var)
4448 ! If a lower energy was found, update the current structure...
4449 if (e_sc.lt.cur_e) then
4451 !v call etotal(energy)
4454 !d e_sc1=wsc*evdw+wscloc*escloc
4455 !d print *,' new',e_sc1,energy(0)
4456 !v print *,'new ',energy(0)
4457 !d call enerprint(energy(0))
4460 if (mask_side(i).eq.1) then
4466 ! ...else revert to the previous one
4469 if (mask_side(i).eq.1) then
4478 n_steps=n_steps+n_try
4480 ! Reset the minimization mask_r to false
4484 end subroutine single_sc_move
4485 !-----------------------------------------------------------------------------
4486 subroutine sc_minimize(etot,iretcode,nfun)
4488 ! Minimizes side-chains only, leaving backbone frozen
4490 use energy, only: etotal
4492 ! implicit real*8 (a-h,o-z)
4493 ! include 'DIMENSIONS'
4494 ! include 'COMMON.VAR'
4495 ! include 'COMMON.CHAIN'
4496 ! include 'COMMON.FFIELD'
4499 real(kind=8) :: etot
4500 integer :: iretcode,nfun
4504 real(kind=8) :: orig_w(0:n_ene),energy_(0:n_ene)
4505 real(kind=8) :: var(6*nres) !(maxvar)(maxvar=6*maxres)
4508 ! Set non side-chain weights to zero (minimization is faster)
4509 ! NOTE: e(2) does not actually depend on the side-chain, only CA
4536 ! Prepare to freeze backbone
4543 ! Minimize the side-chains
4545 call geom_to_var(nvar,var)
4546 call minimize(etot,var,iretcode,nfun)
4547 call var_to_geom(nvar,var)
4550 ! Put the original weights back and calculate the full energy
4565 call etotal(energy_)
4569 end subroutine sc_minimize
4570 !-----------------------------------------------------------------------------
4571 subroutine minimize_sc1(etot,x,iretcode,nfun)
4573 use energy, only: func,gradient,fdum,etotal,enerprint
4575 ! implicit real*8 (a-h,o-z)
4576 ! include 'DIMENSIONS'
4577 integer,parameter :: liv=60
4578 integer :: iretcode,nfun
4579 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4580 ! include 'COMMON.IOUNITS'
4581 ! include 'COMMON.VAR'
4582 ! include 'COMMON.GEO'
4583 ! include 'COMMON.MINIM'
4584 !el integer :: icall
4585 !el common /srutu/ icall
4586 integer,dimension(liv) :: iv
4587 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
4588 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
4589 real(kind=8) :: energia(0:n_ene)
4590 !el real(kind=8) :: fdum
4591 ! external gradient,fdum !func,
4592 ! external func_restr1,grad_restr1
4593 logical :: not_done,change,reduce
4594 !el common /przechowalnia/ v
4596 integer :: nvar_restr,lv,i,j
4598 real(kind=8) :: rdum(1),etot !,fdum
4600 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4601 if (.not. allocated(v)) allocate(v(1:lv))
4603 call deflt(2,iv,liv,lv,v)
4604 ! 12 means fresh start, dont call deflt
4606 ! max num of fun calls
4607 if (maxfun.eq.0) maxfun=500
4609 ! max num of iterations
4610 if (maxmin.eq.0) maxmin=1000
4614 ! selects output unit
4617 ! 1 means to print out result
4619 ! 1 means to print out summary stats
4621 ! 1 means to print initial x and d
4623 ! min val for v(radfac) default is 0.1
4625 ! max val for v(radfac) default is 4.0
4628 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
4629 ! the sumsl default is 0.1
4631 ! false conv if (act fnctn decrease) .lt. v(34)
4632 ! the sumsl default is 100*machep
4634 ! absolute convergence
4635 if (tolf.eq.0.0D0) tolf=1.0D-4
4637 ! relative convergence
4638 if (rtolf.eq.0.0D0) rtolf=1.0D-4
4640 ! controls initial step size
4642 ! large vals of d correspond to small components of step
4651 write(iout,*) "mask_r",mask_r,"petla if minimize_sc1"
4652 call x2xx(x,xx,nvar_restr)
4653 call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,&
4654 iv,liv,lv,v,idum,rdum,fdum)
4657 write(iout,*) "mask_r",mask_r,"petla else minimize_sc1"
4658 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
4660 !el---------------------
4661 ! write (iout,'(/a)') &
4662 ! "Cartesian coordinates of the reference structure after SUMSL"
4663 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
4664 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4666 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
4667 ! restyp(itype(i)),i,(c(j,i),j=1,3),&
4668 ! (c(j,i+nres),j=1,3)
4670 ! call etotal(energia)
4671 ! call enerprint(energia)
4672 !el----------------------------
4678 end subroutine minimize_sc1
4679 !-----------------------------------------------------------------------------
4680 subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
4683 use energy, only: zerograd,esc,sc_grad
4684 ! implicit real*8 (a-h,o-z)
4685 ! include 'DIMENSIONS'
4686 ! include 'COMMON.DERIV'
4687 ! include 'COMMON.IOUNITS'
4688 ! include 'COMMON.GEO'
4689 ! include 'COMMON.FFIELD'
4690 ! include 'COMMON.INTERACT'
4691 ! include 'COMMON.TIME1'
4693 !el common /chuju/ jjj
4694 real(kind=8) :: energia(0:n_ene),evdw,escloc
4695 real(kind=8) :: e1,e2,f
4696 real(kind=8),external :: ufparm
4697 integer :: uiparm(1)
4698 real(kind=8) :: urparm(1)
4699 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
4704 ! Intercept NaNs in the coordinates, before calling etotal
4710 if (x_sum.ne.x_sum) then
4711 write(iout,*)" *** func_restr1 : Found NaN in coordinates"
4718 call var_to_geom_restr(n,x)
4721 !d write (iout,*) 'ETOTAL called from FUNC'
4724 f=wsc*evdw+wscloc*escloc
4725 !d call etotal(energia(0))
4726 !d f=wsc*energia(1)+wscloc*energia(12)
4727 !d print *,f,evdw,escloc,energia(0)
4729 ! Sum up the components of the Cartesian gradient.
4733 gradx(j,i,icg)=wsc*gvdwx(j,i)
4738 end subroutine func_restr1
4739 !-----------------------------------------------------------------------------
4740 subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
4742 use energy, only: cartder,zerograd,esc,sc_grad
4743 ! implicit real*8 (a-h,o-z)
4744 ! include 'DIMENSIONS'
4745 ! include 'COMMON.CHAIN'
4746 ! include 'COMMON.DERIV'
4747 ! include 'COMMON.VAR'
4748 ! include 'COMMON.INTERACT'
4749 ! include 'COMMON.FFIELD'
4750 ! include 'COMMON.IOUNITS'
4752 integer :: i,j,k,ind,n,nf,uiparm(1)
4753 real(kind=8) :: f,urparm(1)
4754 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
4755 integer :: ig,igall,ij
4756 real(kind=8) :: gphii,gthetai,galphai,gomegai
4757 real(kind=8),external :: ufparm
4760 if (nf-nfl+1) 20,30,40
4761 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
4762 ! write (iout,*) 'grad 20'
4765 30 call var_to_geom_restr(n,x)
4768 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
4772 ! Convert the Cartesian gradient into internal-coordinate gradient.
4778 IF (mask_phi(i+2).eq.1) THEN
4783 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
4784 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
4797 IF (mask_theta(i+2).eq.1) THEN
4803 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
4804 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
4814 if (itype(i).ne.10) then
4815 IF (mask_side(i).eq.1) THEN
4819 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
4828 if (itype(i).ne.10) then
4829 IF (mask_side(i).eq.1) THEN
4833 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
4841 ! Add the components corresponding to local energy terms.
4848 if (mask_phi(i).eq.1) then
4850 g(ig)=g(ig)+gloc(igall,icg)
4856 if (mask_theta(i).eq.1) then
4858 g(ig)=g(ig)+gloc(igall,icg)
4864 if (itype(i).ne.10) then
4866 if (mask_side(i).eq.1) then
4868 g(ig)=g(ig)+gloc(igall,icg)
4875 !d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
4878 end subroutine grad_restr1
4879 !-----------------------------------------------------------------------------
4880 subroutine egb1(evdw)
4882 ! This subroutine calculates the interaction energy of nonbonded side chains
4883 ! assuming the Gay-Berne potential of interaction.
4886 use energy, only: sc_grad
4887 ! implicit real*8 (a-h,o-z)
4888 ! include 'DIMENSIONS'
4889 ! include 'COMMON.GEO'
4890 ! include 'COMMON.VAR'
4891 ! include 'COMMON.LOCAL'
4892 ! include 'COMMON.CHAIN'
4893 ! include 'COMMON.DERIV'
4894 ! include 'COMMON.NAMES'
4895 ! include 'COMMON.INTERACT'
4896 ! include 'COMMON.IOUNITS'
4897 ! include 'COMMON.CALC'
4898 ! include 'COMMON.CONTROL'
4900 real(kind=8) :: evdw
4902 integer :: iint,ind,itypi,itypi1,itypj
4903 real(kind=8) :: xi,yi,zi,rrij,sig,sig0ij,rij_shift,fac,e1,e2,&
4905 !elwrite(iout,*) "check evdw"
4906 ! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
4909 ! if (icall.eq.0) lprn=.true.
4911 do i=iatsc_s,iatsc_e
4913 itypi=iabs(itype(i))
4914 itypi1=iabs(itype(i+1))
4918 dxi=dc_norm(1,nres+i)
4919 dyi=dc_norm(2,nres+i)
4920 dzi=dc_norm(3,nres+i)
4921 dsci_inv=dsc_inv(itypi)
4922 !elwrite(iout,*) itypi,itypi1,xi,yi,zi,dxi,dyi,dzi,dsci_inv
4924 ! Calculate SC interaction energy.
4926 do iint=1,nint_gr(i)
4927 do j=istart(i,iint),iend(i,iint)
4928 IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
4930 itypj=iabs(itype(j))
4931 dscj_inv=dsc_inv(itypj)
4932 sig0ij=sigma(itypi,itypj)
4933 chi1=chi(itypi,itypj)
4934 chi2=chi(itypj,itypi)
4941 alf12=0.5D0*(alf1+alf2)
4942 ! For diagnostics only!!!
4955 dxj=dc_norm(1,nres+j)
4956 dyj=dc_norm(2,nres+j)
4957 dzj=dc_norm(3,nres+j)
4958 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
4960 ! Calculate angle-dependent terms of energy and contributions to their
4964 sig=sig0ij*dsqrt(sigsq)
4965 rij_shift=1.0D0/rij-sig+sig0ij
4966 ! I hate to put IF's in the loops, but here don't have another choice!!!!
4967 if (rij_shift.le.0.0D0) then
4969 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
4970 !d restyp(itypi),i,restyp(itypj),j, &
4971 !d rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
4975 !---------------------------------------------------------------
4976 rij_shift=1.0D0/rij_shift
4977 fac=rij_shift**expon
4978 e1=fac*fac*aa(itypi,itypj)
4979 e2=fac*bb(itypi,itypj)
4980 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
4981 eps2der=evdwij*eps3rt
4982 eps3der=evdwij*eps2rt
4983 evdwij=evdwij*eps2rt*eps3rt
4986 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
4987 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
4988 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
4989 !d restyp(itypi),i,restyp(itypj),j, &
4990 !d epsi,sigm,chi1,chi2,chip1,chip2, &
4991 !d eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
4992 !d om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
4996 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
4999 ! Calculate gradient components.
5000 e1=e1*eps1*eps2rt**2*eps3rt**2
5001 fac=-expon*(e1+evdwij)*rij_shift
5004 ! Calculate the radial part of the gradient
5008 ! Calculate angular part of the gradient.
5010 !elwrite(iout,*) evdw
5012 !elwrite(iout,*) "evdw=",evdw,j,iint,i
5014 !elwrite(iout,*) evdw
5016 !elwrite(iout,*) evdw
5018 !elwrite(iout,*) evdw
5020 !elwrite(iout,*) evdw,i
5022 !-----------------------------------------------------------------------------
5024 !-----------------------------------------------------------------------------
5025 subroutine sumsl(n,d,x,calcf,calcg,iv,liv,lv,v,uiparm,urparm,ufparm)
5027 ! *** minimize general unconstrained objective function using ***
5028 ! *** analytic gradient and hessian approx. from secant update ***
5030 integer :: n, liv, lv
5031 integer :: iv(liv), uiparm(1)
5032 real(kind=8) :: d(n), x(n), v(lv), urparm(1)
5033 real(kind=8),external :: ufparm !funtion name as an argument
5035 ! dimension v(71 + n*(n+15)/2), uiparm(*), urparm(*)
5036 external :: calcf, calcg !subroutine name as an argument
5040 ! this routine interacts with subroutine sumit in an attempt
5041 ! to find an n-vector x* that minimizes the (unconstrained)
5042 ! objective function computed by calcf. (often the x* found is
5043 ! a local minimizer rather than a global one.)
5045 !-------------------------- parameter usage --------------------------
5047 ! n........ (input) the number of variables on which f depends, i.e.,
5048 ! the number of components in x.
5049 ! d........ (input/output) a scale vector such that d(i)*x(i),
5050 ! i = 1,2,...,n, are all in comparable units.
5051 ! d can strongly affect the behavior of sumsl.
5052 ! finding the best choice of d is generally a trial-
5053 ! and-error process. choosing d so that d(i)*x(i)
5054 ! has about the same value for all i often works well.
5055 ! the defaults provided by subroutine deflt (see i
5056 ! below) require the caller to supply d.
5057 ! x........ (input/output) before (initially) calling sumsl, the call-
5058 ! er should set x to an initial guess at x*. when
5059 ! sumsl returns, x contains the best point so far
5060 ! found, i.e., the one that gives the least value so
5061 ! far seen for f(x).
5062 ! calcf.... (input) a subroutine that, given x, computes f(x). calcf
5063 ! must be declared external in the calling program.
5065 ! call calcf(n, x, nf, f, uiparm, urparm, ufparm)
5066 ! when calcf is called, nf is the invocation
5067 ! count for calcf. nf is included for possible use
5068 ! with calcg. if x is out of bounds (e.g., if it
5069 ! would cause overflow in computing f(x)), then calcf
5070 ! should set nf to 0. this will cause a shorter step
5071 ! to be attempted. (if x is in bounds, then calcf
5072 ! should not change nf.) the other parameters are as
5073 ! described above and below. calcf should not change
5075 ! calcg.... (input) a subroutine that, given x, computes g(x), the gra-
5076 ! dient of f at x. calcg must be declared external in
5077 ! the calling program. it is invoked by
5078 ! call calcg(n, x, nf, g, uiparm, urparm, ufaprm)
5079 ! when calcg is called, nf is the invocation
5080 ! count for calcf at the time f(x) was evaluated. the
5081 ! x passed to calcg is usually the one passed to calcf
5082 ! on either its most recent invocation or the one
5083 ! prior to it. if calcf saves intermediate results
5084 ! for use by calcg, then it is possible to tell from
5085 ! nf whether they are valid for the current x (or
5086 ! which copy is valid if two copies are kept). if g
5087 ! cannot be computed at x, then calcg should set nf to
5088 ! 0. in this case, sumsl will return with iv(1) = 65.
5089 ! (if g can be computed at x, then calcg should not
5090 ! changed nf.) the other parameters to calcg are as
5091 ! described above and below. calcg should not change
5093 ! iv....... (input/output) an integer value array of length liv (see
5094 ! below) that helps control the sumsl algorithm and
5095 ! that is used to store various intermediate quanti-
5096 ! ties. of particular interest are the initialization/
5097 ! return code iv(1) and the entries in iv that control
5098 ! printing and limit the number of iterations and func-
5099 ! tion evaluations. see the section on iv input
5101 ! liv...... (input) length of iv array. must be at least 60. if li
5102 ! is too small, then sumsl returns with iv(1) = 15.
5103 ! when sumsl returns, the smallest allowed value of
5104 ! liv is stored in iv(lastiv) -- see the section on
5105 ! iv output values below. (this is intended for use
5106 ! with extensions of sumsl that handle constraints.)
5107 ! lv....... (input) length of v array. must be at least 71+n*(n+15)/2.
5108 ! (at least 77+n*(n+17)/2 for smsno, at least
5109 ! 78+n*(n+12) for humsl). if lv is too small, then
5110 ! sumsl returns with iv(1) = 16. when sumsl returns,
5111 ! the smallest allowed value of lv is stored in
5112 ! iv(lastv) -- see the section on iv output values
5114 ! v........ (input/output) a floating-point value array of length l
5115 ! (see below) that helps control the sumsl algorithm
5116 ! and that is used to store various intermediate
5117 ! quantities. of particular interest are the entries
5118 ! in v that limit the length of the first step
5119 ! attempted (lmax0) and specify convergence tolerances
5120 ! (afctol, lmaxs, rfctol, sctol, xctol, xftol).
5121 ! uiparm... (input) user integer parameter array passed without change
5122 ! to calcf and calcg.
5123 ! urparm... (input) user floating-point parameter array passed without
5124 ! change to calcf and calcg.
5125 ! ufparm... (input) user external subroutine or function passed without
5126 ! change to calcf and calcg.
5128 ! *** iv input values (from subroutine deflt) ***
5130 ! iv(1)... on input, iv(1) should have a value between 0 and 14......
5131 ! 0 and 12 mean this is a fresh start. 0 means that
5132 ! deflt(2, iv, liv, lv, v)
5133 ! is to be called to provide all default values to iv and
5134 ! v. 12 (the value that deflt assigns to iv(1)) means the
5135 ! caller has already called deflt and has possibly changed
5136 ! some iv and/or v entries to non-default values.
5137 ! 13 means deflt has been called and that sumsl (and
5138 ! sumit) should only do their storage allocation. that is,
5139 ! they should set the output components of iv that tell
5140 ! where various subarrays arrays of v begin, such as iv(g)
5141 ! (and, for humsl and humit only, iv(dtol)), and return.
5142 ! 14 means that a storage has been allocated (by a call
5143 ! with iv(1) = 13) and that the algorithm should be
5144 ! started. when called with iv(1) = 13, sumsl returns
5145 ! iv(1) = 14 unless liv or lv is too small (or n is not
5146 ! positive). default = 12.
5147 ! iv(inith).... iv(25) tells whether the hessian approximation h should
5148 ! be initialized. 1 (the default) means sumit should
5149 ! initialize h to the diagonal matrix whose i-th diagonal
5150 ! element is d(i)**2. 0 means the caller has supplied a
5151 ! cholesky factor l of the initial hessian approximation
5152 ! h = l*(l**t) in v, starting at v(iv(lmat)) = v(iv(42))
5153 ! (and stored compactly by rows). note that iv(lmat) may
5154 ! be initialized by calling sumsl with iv(1) = 13 (see
5155 ! the iv(1) discussion above). default = 1.
5156 ! iv(mxfcal)... iv(17) gives the maximum number of function evaluations
5157 ! (calls on calcf) allowed. if this number does not suf-
5158 ! fice, then sumsl returns with iv(1) = 9. default = 200.
5159 ! iv(mxiter)... iv(18) gives the maximum number of iterations allowed.
5160 ! it also indirectly limits the number of gradient evalua-
5161 ! tions (calls on calcg) to iv(mxiter) + 1. if iv(mxiter)
5162 ! iterations do not suffice, then sumsl returns with
5163 ! iv(1) = 10. default = 150.
5164 ! iv(outlev)... iv(19) controls the number and length of iteration sum-
5165 ! mary lines printed (by itsum). iv(outlev) = 0 means do
5166 ! not print any summary lines. otherwise, print a summary
5167 ! line after each abs(iv(outlev)) iterations. if iv(outlev)
5168 ! is positive, then summary lines of length 78 (plus carri-
5169 ! age control) are printed, including the following... the
5170 ! iteration and function evaluation counts, f = the current
5171 ! function value, relative difference in function values
5172 ! achieved by the latest step (i.e., reldf = (f0-v(f))/f01,
5173 ! where f01 is the maximum of abs(v(f)) and abs(v(f0)) and
5174 ! v(f0) is the function value from the previous itera-
5175 ! tion), the relative function reduction predicted for the
5176 ! step just taken (i.e., preldf = v(preduc) / f01, where
5177 ! v(preduc) is described below), the scaled relative change
5178 ! in x (see v(reldx) below), the step parameter for the
5179 ! step just taken (stppar = 0 means a full newton step,
5180 ! between 0 and 1 means a relaxed newton step, between 1
5181 ! and 2 means a double dogleg step, greater than 2 means
5182 ! a scaled down cauchy step -- see subroutine dbldog), the
5183 ! 2-norm of the scale vector d times the step just taken
5184 ! (see v(dstnrm) below), and npreldf, i.e.,
5185 ! v(nreduc)/f01, where v(nreduc) is described below -- if
5186 ! npreldf is positive, then it is the relative function
5187 ! reduction predicted for a newton step (one with
5188 ! stppar = 0). if npreldf is negative, then it is the
5189 ! negative of the relative function reduction predicted
5190 ! for a step computed with step bound v(lmaxs) for use in
5191 ! testing for singular convergence.
5192 ! if iv(outlev) is negative, then lines of length 50
5193 ! are printed, including only the first 6 items listed
5194 ! above (through reldx).
5196 ! iv(parprt)... iv(20) = 1 means print any nondefault v values on a
5197 ! fresh start or any changed v values on a restart.
5198 ! iv(parprt) = 0 means skip this printing. default = 1.
5199 ! iv(prunit)... iv(21) is the output unit number on which all printing
5200 ! is done. iv(prunit) = 0 means suppress all printing.
5201 ! default = standard output unit (unit 6 on most systems).
5202 ! iv(solprt)... iv(22) = 1 means print out the value of x returned (as
5203 ! well as the gradient and the scale vector d).
5204 ! iv(solprt) = 0 means skip this printing. default = 1.
5205 ! iv(statpr)... iv(23) = 1 means print summary statistics upon return-
5206 ! ing. these consist of the function value, the scaled
5207 ! relative change in x caused by the most recent step (see
5208 ! v(reldx) below), the number of function and gradient
5209 ! evaluations (calls on calcf and calcg), and the relative
5210 ! function reductions predicted for the last step taken and
5211 ! for a newton step (or perhaps a step bounded by v(lmaxs)
5212 ! -- see the descriptions of preldf and npreldf under
5213 ! iv(outlev) above).
5214 ! iv(statpr) = 0 means skip this printing.
5215 ! iv(statpr) = -1 means skip this printing as well as that
5216 ! of the one-line termination reason message. default = 1.
5217 ! iv(x0prt).... iv(24) = 1 means print the initial x and scale vector d
5218 ! (on a fresh start only). iv(x0prt) = 0 means skip this
5219 ! printing. default = 1.
5221 ! *** (selected) iv output values ***
5223 ! iv(1)........ on output, iv(1) is a return code....
5224 ! 3 = x-convergence. the scaled relative difference (see
5225 ! v(reldx)) between the current parameter vector x and
5226 ! a locally optimal parameter vector is very likely at
5228 ! 4 = relative function convergence. the relative differ-
5229 ! ence between the current function value and its lo-
5230 ! cally optimal value is very likely at most v(rfctol).
5231 ! 5 = both x- and relative function convergence (i.e., the
5232 ! conditions for iv(1) = 3 and iv(1) = 4 both hold).
5233 ! 6 = absolute function convergence. the current function
5234 ! value is at most v(afctol) in absolute value.
5235 ! 7 = singular convergence. the hessian near the current
5236 ! iterate appears to be singular or nearly so, and a
5237 ! step of length at most v(lmaxs) is unlikely to yield
5238 ! a relative function decrease of more than v(sctol).
5239 ! 8 = false convergence. the iterates appear to be converg-
5240 ! ing to a noncritical point. this may mean that the
5241 ! convergence tolerances (v(afctol), v(rfctol),
5242 ! v(xctol)) are too small for the accuracy to which
5243 ! the function and gradient are being computed, that
5244 ! there is an error in computing the gradient, or that
5245 ! the function or gradient is discontinuous near x.
5246 ! 9 = function evaluation limit reached without other con-
5247 ! vergence (see iv(mxfcal)).
5248 ! 10 = iteration limit reached without other convergence
5250 ! 11 = stopx returned .true. (external interrupt). see the
5251 ! usage notes below.
5252 ! 14 = storage has been allocated (after a call with
5254 ! 17 = restart attempted with n changed.
5255 ! 18 = d has a negative component and iv(dtype) .le. 0.
5256 ! 19...43 = v(iv(1)) is out of range.
5257 ! 63 = f(x) cannot be computed at the initial x.
5258 ! 64 = bad parameters passed to assess (which should not
5260 ! 65 = the gradient could not be computed at x (see calcg
5262 ! 67 = bad first parameter to deflt.
5263 ! 80 = iv(1) was out of range.
5264 ! 81 = n is not positive.
5265 ! iv(g)........ iv(28) is the starting subscript in v of the current
5266 ! gradient vector (the one corresponding to x).
5267 ! iv(lastiv)... iv(44) is the least acceptable value of liv. (it is
5268 ! only set if liv is at least 44.)
5269 ! iv(lastv).... iv(45) is the least acceptable value of lv. (it is
5270 ! only set if liv is large enough, at least iv(lastiv).)
5271 ! iv(nfcall)... iv(6) is the number of calls so far made on calcf (i.e.,
5272 ! function evaluations).
5273 ! iv(ngcall)... iv(30) is the number of gradient evaluations (calls on
5275 ! iv(niter).... iv(31) is the number of iterations performed.
5277 ! *** (selected) v input values (from subroutine deflt) ***
5279 ! v(bias)..... v(43) is the bias parameter used in subroutine dbldog --
5280 ! see that subroutine for details. default = 0.8.
5281 ! v(afctol)... v(31) is the absolute function convergence tolerance.
5282 ! if sumsl finds a point where the function value is less
5283 ! than v(afctol) in absolute value, and if sumsl does not
5284 ! return with iv(1) = 3, 4, or 5, then it returns with
5285 ! iv(1) = 6. this test can be turned off by setting
5286 ! v(afctol) to zero. default = max(10**-20, machep**2),
5287 ! where machep is the unit roundoff.
5288 ! v(dinit).... v(38), if nonnegative, is the value to which the scale
5289 ! vector d is initialized. default = -1.
5290 ! v(lmax0).... v(35) gives the maximum 2-norm allowed for d times the
5291 ! very first step that sumsl attempts. this parameter can
5292 ! markedly affect the performance of sumsl.
5293 ! v(lmaxs).... v(36) is used in testing for singular convergence -- if
5294 ! the function reduction predicted for a step of length
5295 ! bounded by v(lmaxs) is at most v(sctol) * abs(f0), where
5296 ! f0 is the function value at the start of the current
5297 ! iteration, and if sumsl does not return with iv(1) = 3,
5298 ! 4, 5, or 6, then it returns with iv(1) = 7. default = 1.
5299 ! v(rfctol)... v(32) is the relative function convergence tolerance.
5300 ! if the current model predicts a maximum possible function
5301 ! reduction (see v(nreduc)) of at most v(rfctol)*abs(f0)
5302 ! at the start of the current iteration, where f0 is the
5303 ! then current function value, and if the last step attempt-
5304 ! ed achieved no more than twice the predicted function
5305 ! decrease, then sumsl returns with iv(1) = 4 (or 5).
5306 ! default = max(10**-10, machep**(2/3)), where machep is
5307 ! the unit roundoff.
5308 ! v(sctol).... v(37) is the singular convergence tolerance -- see the
5309 ! description of v(lmaxs) above.
5310 ! v(tuner1)... v(26) helps decide when to check for false convergence.
5311 ! this is done if the actual function decrease from the
5312 ! current step is no more than v(tuner1) times its predict-
5313 ! ed value. default = 0.1.
5314 ! v(xctol).... v(33) is the x-convergence tolerance. if a newton step
5315 ! (see v(nreduc)) is tried that has v(reldx) .le. v(xctol)
5316 ! and if this step yields at most twice the predicted func-
5317 ! tion decrease, then sumsl returns with iv(1) = 3 (or 5).
5318 ! (see the description of v(reldx) below.)
5319 ! default = machep**0.5, where machep is the unit roundoff.
5320 ! v(xftol).... v(34) is the false convergence tolerance. if a step is
5321 ! tried that gives no more than v(tuner1) times the predict-
5322 ! ed function decrease and that has v(reldx) .le. v(xftol),
5323 ! and if sumsl does not return with iv(1) = 3, 4, 5, 6, or
5324 ! 7, then it returns with iv(1) = 8. (see the description
5325 ! of v(reldx) below.) default = 100*machep, where
5326 ! machep is the unit roundoff.
5327 ! v(*)........ deflt supplies to v a number of tuning constants, with
5328 ! which it should ordinarily be unnecessary to tinker. see
5329 ! section 17 of version 2.2 of the nl2sol usage summary
5330 ! (i.e., the appendix to ref. 1) for details on v(i),
5331 ! i = decfac, incfac, phmnfc, phmxfc, rdfcmn, rdfcmx,
5332 ! tuner2, tuner3, tuner4, tuner5.
5334 ! *** (selected) v output values ***
5336 ! v(dgnorm)... v(1) is the 2-norm of (diag(d)**-1)*g, where g is the
5337 ! most recently computed gradient.
5338 ! v(dstnrm)... v(2) is the 2-norm of diag(d)*step, where step is the
5340 ! v(f)........ v(10) is the current function value.
5341 ! v(f0)....... v(13) is the function value at the start of the current
5343 ! v(nreduc)... v(6), if positive, is the maximum function reduction
5344 ! possible according to the current model, i.e., the func-
5345 ! tion reduction predicted for a newton step (i.e.,
5346 ! step = -h**-1 * g, where g is the current gradient and
5347 ! h is the current hessian approximation).
5348 ! if v(nreduc) is negative, then it is the negative of
5349 ! the function reduction predicted for a step computed with
5350 ! a step bound of v(lmaxs) for use in testing for singular
5352 ! v(preduc)... v(7) is the function reduction predicted (by the current
5353 ! quadratic model) for the current step. this (divided by
5354 ! v(f0)) is used in testing for relative function
5356 ! v(reldx).... v(17) is the scaled relative change in x caused by the
5357 ! current step, computed as
5358 ! max(abs(d(i)*(x(i)-x0(i)), 1 .le. i .le. p) /
5359 ! max(d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p),
5360 ! where x = x0 + step.
5362 !------------------------------- notes -------------------------------
5364 ! *** algorithm notes ***
5366 ! this routine uses a hessian approximation computed from the
5367 ! bfgs update (see ref 3). only a cholesky factor of the hessian
5368 ! approximation is stored, and this is updated using ideas from
5369 ! ref. 4. steps are computed by the double dogleg scheme described
5370 ! in ref. 2. the steps are assessed as in ref. 1.
5372 ! *** usage notes ***
5374 ! after a return with iv(1) .le. 11, it is possible to restart,
5375 ! i.e., to change some of the iv and v input values described above
5376 ! and continue the algorithm from the point where it was interrupt-
5377 ! ed. iv(1) should not be changed, nor should any entries of i
5378 ! and v other than the input values (those supplied by deflt).
5379 ! those who do not wish to write a calcg which computes the
5380 ! gradient analytically should call smsno rather than sumsl.
5381 ! smsno uses finite differences to compute an approximate gradient.
5382 ! those who would prefer to provide f and g (the function and
5383 ! gradient) by reverse communication rather than by writing subrou-
5384 ! tines calcf and calcg may call on sumit directly. see the com-
5385 ! ments at the beginning of sumit.
5386 ! those who use sumsl interactively may wish to supply their
5387 ! own stopx function, which should return .true. if the break key
5388 ! has been pressed since stopx was last invoked. this makes it
5389 ! possible to externally interrupt sumsl (which will return with
5390 ! iv(1) = 11 if stopx returns .true.).
5391 ! storage for g is allocated at the end of v. thus the caller
5392 ! may make v longer than specified above and may allow calcg to use
5393 ! elements of g beyond the first n as scratch storage.
5395 ! *** portability notes ***
5397 ! the sumsl distribution tape contains both single- and double-
5398 ! precision versions of the sumsl source code, so it should be un-
5399 ! necessary to change precisions.
5400 ! only the functions imdcon and rmdcon contain machine-dependent
5401 ! constants. to change from one machine to another, it should
5402 ! suffice to change the (few) relevant lines in these functions.
5403 ! intrinsic functions are explicitly declared. on certain com-
5404 ! puters (e.g. univac), it may be necessary to comment out these
5405 ! declarations. so that this may be done automatically by a simple
5406 ! program, such declarations are preceded by a comment having c/+
5407 ! in columns 1-3 and blanks in columns 4-72 and are followed by
5408 ! a comment having c/ in columns 1 and 2 and blanks in columns 3-72.
5409 ! the sumsl source code is expressed in 1966 ansi standard
5410 ! fortran. it may be converted to fortran 77 by commenting out all
5411 ! lines that fall between a line having c/6 in columns 1-3 and a
5412 ! line having c/7 in columns 1-3 and by removing (i.e., replacing
5413 ! by a blank) the c in column 1 of the lines that follow the c/7
5414 ! line and precede a line having c/ in columns 1-2 and blanks in
5415 ! columns 3-72. these changes convert some data statements into
5416 ! parameter statements, convert some variables from real to
5417 ! character*4, and make the data statements that initialize these
5418 ! variables use character strings delimited by primes instead
5419 ! of hollerith constants. (such variables and data statements
5420 ! appear only in modules itsum and parck. parameter statements
5421 ! appear nearly everywhere.) these changes also add save state-
5422 ! ments for variables given machine-dependent constants by rmdcon.
5424 ! *** references ***
5426 ! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), algorithm 573 --
5427 ! an adaptive nonlinear least-squares algorithm, acm trans.
5428 ! math. software 7, pp. 369-383.
5430 ! 2. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
5431 ! mization algorithms which use function and gradient
5432 ! values, j. optim. theory applic. 28, pp. 453-482.
5434 ! 3. dennis, j.e., and more, j.j. (1977), quasi-newton methods, motiva-
5435 ! tion and theory, siam rev. 19, pp. 46-89.
5437 ! 4. goldfarb, d. (1976), factorized variable metric methods for uncon-
5438 ! strained optimization, math. comput. 30, pp. 796-811.
5442 ! coded by david m. gay (winter 1980). revised summer 1982.
5443 ! this subroutine was written in connection with research
5444 ! supported in part by the national science foundation under
5445 ! grants mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989,
5449 !---------------------------- declarations ---------------------------
5451 !el external deflt, sumit
5453 ! deflt... supplies default iv and v input components.
5454 ! sumit... reverse-communication routine that carries out sumsl algo-
5457 integer :: g1, iv1, nf
5460 ! *** subscripts for iv ***
5462 !el integer nextv, nfcall, nfgcal, g, toobig, vneed
5465 ! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, toobig/2/, vneed/4/
5467 integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28,&
5471 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
5473 !elwrite(iout,*) "in sumsl"
5474 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
5476 if (iv1 .eq. 12 .or. iv1 .eq. 13) iv(vneed) = iv(vneed) + n
5477 if (iv1 .eq. 14) go to 10
5478 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
5480 if (iv1 .eq. 12) iv(1) = 13
5484 !elwrite(iout,*) "in sumsl go to 10"
5487 !elwrite(iout,*) "in sumsl"
5488 20 call sumit(d, f, v(g1), iv, liv, lv, n, v, x)
5489 !elwrite(iout,*) "in sumsl, go to 20"
5491 !elwrite(iout,*) "in sumsl, go to 20, po sumit"
5492 !elwrite(iout,*) "in sumsl iv()", iv(1)-2
5493 if (iv(1) - 2) 30, 40, 50
5496 !elwrite(iout,*) "in sumsl iv",iv(nfcall)
5497 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
5498 !elwrite(iout,*) "in sumsl"
5499 if (nf .le. 0) iv(toobig) = 1
5502 !elwrite(iout,*) "in sumsl"
5503 40 call calcg(n, x, iv(nfgcal), v(g1), uiparm, urparm, ufparm)
5504 !elwrite(iout,*) "in sumsl"
5507 50 if (iv(1) .ne. 14) go to 999
5509 ! *** storage allocation
5512 iv(nextv) = iv(g) + n
5513 if (iv1 .ne. 13) go to 10
5514 !elwrite(iout,*) "in sumsl"
5517 ! *** last card of sumsl follows ***
5518 end subroutine sumsl
5519 !-----------------------------------------------------------------------------
5520 subroutine sumit(d,fx,g,iv,liv,lv,n,v,x)
5522 use control, only:stopx
5524 ! *** carry out sumsl (unconstrained minimization) iterations, using
5525 ! *** double-dogleg/bfgs steps.
5527 ! *** parameter declarations ***
5529 integer :: liv, lv, n
5531 real(kind=8) :: d(n), fx, g(n), v(lv), x(n)
5533 !-------------------------- parameter usage --------------------------
5535 ! d.... scale vector.
5536 ! fx... function value.
5537 ! g.... gradient vector.
5538 ! iv... integer value array.
5539 ! liv.. length of iv (at least 60).
5540 ! lv... length of v (at least 71 + n*(n+13)/2).
5541 ! n.... number of variables (components in x and g).
5542 ! v.... floating-point value array.
5543 ! x.... vector of parameters to be optimized.
5545 ! *** discussion ***
5547 ! parameters iv, n, v, and x are the same as the corresponding
5548 ! ones to sumsl (which see), except that v can be shorter (since
5549 ! the part of v that sumsl uses for storing g is not needed).
5550 ! moreover, compared with sumsl, iv(1) may have the two additional
5551 ! output values 1 and 2, which are explained below, as is the use
5552 ! of iv(toobig) and iv(nfgcal). the value iv(g), which is an
5553 ! output value from sumsl (and smsno), is not referenced by
5554 ! sumit or the subroutines it calls.
5555 ! fx and g need not have been initialized when sumit is called
5556 ! with iv(1) = 12, 13, or 14.
5558 ! iv(1) = 1 means the caller should set fx to f(x), the function value
5559 ! at x, and call sumit again, having changed none of the
5560 ! other parameters. an exception occurs if f(x) cannot be
5561 ! (e.g. if overflow would occur), which may happen because
5562 ! of an oversized step. in this case the caller should set
5563 ! iv(toobig) = iv(2) to 1, which will cause sumit to ig-
5564 ! nore fx and try a smaller step. the parameter nf that
5565 ! sumsl passes to calcf (for possible use by calcg) is a
5566 ! copy of iv(nfcall) = iv(6).
5567 ! iv(1) = 2 means the caller should set g to g(x), the gradient vector
5568 ! of f at x, and call sumit again, having changed none of
5569 ! the other parameters except possibly the scale vector d
5570 ! when iv(dtype) = 0. the parameter nf that sumsl passes
5571 ! to calcg is iv(nfgcal) = iv(7). if g(x) cannot be
5572 ! evaluated, then the caller may set iv(nfgcal) to 0, in
5573 ! which case sumit will return with iv(1) = 65.
5577 ! coded by david m. gay (december 1979). revised sept. 1982.
5578 ! this subroutine was written in connection with research supported
5579 ! in part by the national science foundation under grants
5580 ! mcs-7600324 and mcs-7906671.
5582 ! (see sumsl for references.)
5584 !+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
5586 ! *** local variables ***
5588 integer :: dg1, dummy, g01, i, k, l, lstgst, nwtst1, step1,&
5591 !el logical :: lstopx
5595 !el real(kind=8) :: half, negone, one, onep2, zero
5597 ! *** no intrinsic functions ***
5599 ! *** external functions and subroutines ***
5601 !el external assst, dbdog, deflt, dotprd, itsum, litvmu, livmul,
5602 !el 1 ltvmul, lupdat, lvmul, parck, reldst, stopx, vaxpy,
5603 !el 2 vcopy, vscopy, vvmulp, v2norm, wzbfgs
5605 !el real(kind=8) :: dotprd, reldst, v2norm
5607 ! assst.... assesses candidate step.
5608 ! dbdog.... computes double-dogleg (candidate) step.
5609 ! deflt.... supplies default iv and v input components.
5610 ! dotprd... returns inner product of two vectors.
5611 ! itsum.... prints iteration summary and info on initial and final x.
5612 ! litvmu... multiplies inverse transpose of lower triangle times vector.
5613 ! livmul... multiplies inverse of lower triangle times vector.
5614 ! ltvmul... multiplies transpose of lower triangle times vector.
5615 ! lupdt.... updates cholesky factor of hessian approximation.
5616 ! lvmul.... multiplies lower triangle times vector.
5617 ! parck.... checks validity of input iv and v values.
5618 ! reldst... computes v(reldx) = relative step size.
5619 ! stopx.... returns .true. if the break key has been pressed.
5620 ! vaxpy.... computes scalar times one vector plus another.
5621 ! vcopy.... copies one vector to another.
5622 ! vscopy... sets all elements of a vector to a scalar.
5623 ! vvmulp... multiplies vector by vector raised to power (componentwise).
5624 ! v2norm... returns the 2-norm of a vector.
5625 ! wzbfgs... computes w and z for lupdat corresponding to bfgs update.
5627 ! *** subscripts for iv and v ***
5630 !el integer cnvcod, dg, dgnorm, dinit, dstnrm, dst0, f, f0, fdif,
5631 !el 1 gthg, gtstep, g0, incfac, inith, irc, kagqt, lmat, lmax0,
5632 !el 2 lmaxs, mode, model, mxfcal, mxiter, nextv, nfcall, nfgcal,
5633 !el 3 ngcall, niter, nreduc, nwtstp, preduc, radfac, radinc,
5634 !el 4 radius, rad0, reldx, restor, step, stglim, stlstg, toobig,
5635 !el 5 tuner4, tuner5, vneed, xirc, x0
5637 ! *** iv subscript values ***
5640 ! data cnvcod/55/, dg/37/, g0/48/, inith/25/, irc/29/, kagqt/33/,
5641 ! 1 mode/35/, model/5/, mxfcal/17/, mxiter/18/, nfcall/6/,
5642 ! 2 nfgcal/7/, ngcall/30/, niter/31/, nwtstp/34/, radinc/8/,
5643 ! 3 restor/9/, step/40/, stglim/11/, stlstg/41/, toobig/2/,
5644 ! 4 vneed/4/, xirc/13/, x0/43/
5646 integer,parameter :: cnvcod=55, dg=37, g0=48, inith=25, irc=29, kagqt=33,&
5647 mode=35, model=5, mxfcal=17, mxiter=18, nfcall=6,&
5648 nfgcal=7, ngcall=30, niter=31, nwtstp=34, radinc=8,&
5649 restor=9, step=40, stglim=11, stlstg=41, toobig=2,&
5650 vneed=4, xirc=13, x0=43
5653 ! *** v subscript values ***
5657 ! data dgnorm/1/, dinit/38/, dstnrm/2/, dst0/3/, f/10/, f0/13/,
5658 ! 1 fdif/11/, gthg/44/, gtstep/4/, incfac/23/, lmat/42/,
5659 ! 2 lmax0/35/, lmaxs/36/, nextv/47/, nreduc/6/, preduc/7/,
5660 ! 3 radfac/16/, radius/8/, rad0/9/, reldx/17/, tuner4/29/,
5663 integer,parameter :: afctol=31
5664 integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dst0=3, f=10, f0=13,&
5665 fdif=11, gthg=44, gtstep=4, incfac=23, lmat=42,&
5666 lmax0=35, lmaxs=36, nextv=47, nreduc=6, preduc=7,&
5667 radfac=16, radius=8, rad0=9, reldx=17, tuner4=29,&
5672 ! data half/0.5d+0/, negone/-1.d+0/, one/1.d+0/, onep2/1.2d+0/,
5675 real(kind=8),parameter :: half=0.5d+0, negone=-1.d+0, one=1.d+0,&
5676 onep2=1.2d+0,zero=0.d+0
5679 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
5681 ! Following SAVE statement inserted.
5684 if (i .eq. 1) go to 50
5685 if (i .eq. 2) go to 60
5687 ! *** check validity of iv and v input values ***
5689 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
5690 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
5691 iv(vneed) = iv(vneed) + n*(n+13)/2
5692 call parck(2, d, iv, liv, lv, n, v)
5694 if (i .gt. 12) go to 999
5695 go to (180, 180, 180, 180, 180, 180, 120, 90, 120, 10, 10, 20), i
5697 ! *** storage allocation ***
5700 iv(x0) = l + n*(n+1)/2
5701 iv(step) = iv(x0) + n
5702 iv(stlstg) = iv(step) + n
5703 iv(g0) = iv(stlstg) + n
5704 iv(nwtstp) = iv(g0) + n
5705 iv(dg) = iv(nwtstp) + n
5706 iv(nextv) = iv(dg) + n
5707 if (iv(1) .ne. 13) go to 20
5711 ! *** initialization ***
5724 if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
5725 if (iv(inith) .ne. 1) go to 40
5727 ! *** set the initial hessian approximation to diag(d)**-2 ***
5730 call vscopy(n*(n+1)/2, v(l), zero)
5735 if (t .le. zero) t = one
5739 ! *** compute initial function value ***
5745 if (iv(mode) .ge. 0) go to 180
5747 if (iv(toobig) .eq. 0) go to 999
5751 ! *** make sure gradient could be computed ***
5753 60 if (iv(nfgcal) .ne. 0) go to 70
5758 call vvmulp(n, v(dg1), g, d, -1)
5759 v(dgnorm) = v2norm(n, v(dg1))
5761 ! *** test norm of gradient ***
5763 if (v(dgnorm) .gt. v(afctol)) go to 75
5765 iv(cnvcod) = iv(irc) - 4
5767 75 if (iv(cnvcod) .ne. 0) go to 290
5768 if (iv(mode) .eq. 0) go to 250
5770 ! *** allow first step to have scaled 2-norm at most v(lmax0) ***
5772 v(radius) = v(lmax0)
5777 !----------------------------- main loop -----------------------------
5780 ! *** print iteration summary, check iteration limit ***
5782 80 call itsum(d, g, iv, liv, lv, n, v, x)
5784 if (k .lt. iv(mxiter)) go to 100
5788 ! *** update radius ***
5790 100 iv(niter) = k + 1
5791 if(k.gt.0)v(radius) = v(radfac) * v(dstnrm)
5793 ! *** initialize for start of next iteration ***
5801 ! *** copy x to x0, g to g0 ***
5803 call vcopy(n, v(x01), x)
5804 call vcopy(n, v(g01), g)
5806 ! *** check stopx and function evaluation limit ***
5810 !el lstopx = stopx(dummy)
5811 !elwrite(iout,*) "lstopx",lstopx,dummy
5812 110 if (.not. stopx(dummy)) go to 130
5814 ! write (iout,*) "iv(1)=11 !!!!"
5817 ! *** come here when restarting after func. eval. limit or stopx.
5819 120 if (v(f) .ge. v(f0)) go to 130
5824 130 if (iv(nfcall) .lt. iv(mxfcal)) go to 150
5826 140 if (v(f) .ge. v(f0)) go to 300
5828 ! *** in case of stopx or function evaluation limit with
5829 ! *** improved v(f), evaluate the gradient at x.
5834 !. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
5836 150 step1 = iv(step)
5839 if (iv(kagqt) .ge. 0) go to 160
5841 call livmul(n, v(nwtst1), v(l), g)
5842 v(nreduc) = half * dotprd(n, v(nwtst1), v(nwtst1))
5843 call litvmu(n, v(nwtst1), v(l), v(nwtst1))
5844 call vvmulp(n, v(step1), v(nwtst1), d, 1)
5845 v(dst0) = v2norm(n, v(step1))
5846 call vvmulp(n, v(dg1), v(dg1), d, -1)
5847 call ltvmul(n, v(step1), v(l), v(dg1))
5848 v(gthg) = v2norm(n, v(step1))
5850 160 call dbdog(v(dg1), lv, n, v(nwtst1), v(step1), v)
5851 if (iv(irc) .eq. 6) go to 180
5853 ! *** check whether evaluating f(x0 + step) looks worthwhile ***
5855 if (v(dstnrm) .le. zero) go to 180
5856 if (iv(irc) .ne. 5) go to 170
5857 if (v(radfac) .le. one) go to 170
5858 if (v(preduc) .le. onep2 * v(fdif)) go to 180
5860 ! *** compute f(x0 + step) ***
5864 call vaxpy(n, x, one, v(step1), v(x01))
5865 iv(nfcall) = iv(nfcall) + 1
5870 !. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
5873 v(reldx) = reldst(n, d, x, v(x01))
5874 call assst(iv, liv, lv, v)
5877 if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
5878 if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
5879 if (iv(restor) .ne. 3) go to 190
5880 call vcopy(n, v(step1), v(lstgst))
5881 call vaxpy(n, x, one, v(step1), v(x01))
5882 v(reldx) = reldst(n, d, x, v(x01))
5885 go to (200,230,230,230,200,210,220,220,220,220,220,220,280,250), k
5887 ! *** recompute step with changed radius ***
5889 200 v(radius) = v(radfac) * v(dstnrm)
5892 ! *** compute step of length v(lmaxs) for singular convergence test.
5894 210 v(radius) = v(lmaxs)
5897 ! *** convergence or false convergence ***
5899 220 iv(cnvcod) = k - 4
5900 if (v(f) .ge. v(f0)) go to 290
5901 if (iv(xirc) .eq. 14) go to 290
5904 !. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
5906 230 if (iv(irc) .ne. 3) go to 240
5910 ! *** set temp1 = hessian * step for use in gradient tests ***
5913 call ltvmul(n, v(temp1), v(l), v(step1))
5914 call lvmul(n, v(temp1), v(l), v(temp1))
5916 ! *** compute gradient ***
5918 240 iv(ngcall) = iv(ngcall) + 1
5922 ! *** initializations -- g0 = g - g0, etc. ***
5925 call vaxpy(n, v(g01), negone, v(g01), g)
5928 if (iv(irc) .ne. 3) go to 270
5930 ! *** set v(radfac) by gradient tests ***
5932 ! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
5934 call vaxpy(n, v(temp1), negone, v(g01), v(temp1))
5935 call vvmulp(n, v(temp1), v(temp1), d, -1)
5937 ! *** do gradient tests ***
5939 if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) &
5941 if (dotprd(n, g, v(step1)) &
5942 .ge. v(gtstep) * v(tuner5)) go to 270
5943 260 v(radfac) = v(incfac)
5945 ! *** update h, loop ***
5950 call wzbfgs(v(l), n, v(step1), v(w), v(g01), v(z))
5952 ! ** use the n-vectors starting at v(step1) and v(g01) for scratch..
5953 call lupdat(v(temp1), v(step1), v(l), v(g01), v(l), n, v(w), v(z))
5957 !. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
5959 ! *** bad parameters to assess ***
5964 ! *** print summary of final iteration and other requested items ***
5966 290 iv(1) = iv(cnvcod)
5968 300 call itsum(d, g, iv, liv, lv, n, v, x)
5972 ! *** last line of sumit follows ***
5973 end subroutine sumit
5974 !-----------------------------------------------------------------------------
5975 subroutine dbdog(dig,lv,n,nwtstp,step,v)
5977 ! *** compute double dogleg step ***
5979 ! *** parameter declarations ***
5982 real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv)
5986 ! this subroutine computes a candidate step (for use in an uncon-
5987 ! strained minimization code) by the double dogleg algorithm of
5988 ! dennis and mei (ref. 1), which is a variation on powell*s dogleg
5989 ! scheme (ref. 2, p. 95).
5991 !-------------------------- parameter usage --------------------------
5993 ! dig (input) diag(d)**-2 * g -- see algorithm notes.
5994 ! g (input) the current gradient vector.
5995 ! lv (input) length of v.
5996 ! n (input) number of components in dig, g, nwtstp, and step.
5997 ! nwtstp (input) negative newton step -- see algorithm notes.
5998 ! step (output) the computed step.
5999 ! v (i/o) values array, the following components of which are
6001 ! v(bias) (input) bias for relaxed newton step, which is v(bias) of
6002 ! the way from the full newton to the fully relaxed newton
6003 ! step. recommended value = 0.8 .
6004 ! v(dgnorm) (input) 2-norm of diag(d)**-1 * g -- see algorithm notes.
6005 ! v(dstnrm) (output) 2-norm of diag(d) * step, which is v(radius)
6006 ! unless v(stppar) = 0 -- see algorithm notes.
6007 ! v(dst0) (input) 2-norm of diag(d) * nwtstp -- see algorithm notes.
6008 ! v(grdfac) (output) the coefficient of dig in the step returned --
6009 ! step(i) = v(grdfac)*dig(i) + v(nwtfac)*nwtstp(i).
6010 ! v(gthg) (input) square-root of (dig**t) * (hessian) * dig -- see
6012 ! v(gtstep) (output) inner product between g and step.
6013 ! v(nreduc) (output) function reduction predicted for the full newton
6015 ! v(nwtfac) (output) the coefficient of nwtstp in the step returned --
6016 ! see v(grdfac) above.
6017 ! v(preduc) (output) function reduction predicted for the step returned.
6018 ! v(radius) (input) the trust region radius. d times the step returned
6019 ! has 2-norm v(radius) unless v(stppar) = 0.
6020 ! v(stppar) (output) code telling how step was computed... 0 means a
6021 ! full newton step. between 0 and 1 means v(stppar) of the
6022 ! way from the newton to the relaxed newton step. between
6023 ! 1 and 2 means a true double dogleg step, v(stppar) - 1 of
6024 ! the way from the relaxed newton to the cauchy step.
6025 ! greater than 2 means 1 / (v(stppar) - 1) times the cauchy
6028 !------------------------------- notes -------------------------------
6030 ! *** algorithm notes ***
6032 ! let g and h be the current gradient and hessian approxima-
6033 ! tion respectively and let d be the current scale vector. this
6034 ! routine assumes dig = diag(d)**-2 * g and nwtstp = h**-1 * g.
6035 ! the step computed is the same one would get by replacing g and h
6036 ! by diag(d)**-1 * g and diag(d)**-1 * h * diag(d)**-1,
6037 ! computing step, and translating step back to the original
6038 ! variables, i.e., premultiplying it by diag(d)**-1.
6040 ! *** references ***
6042 ! 1. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
6043 ! mization algorithms which use function and gradient
6044 ! values, j. optim. theory applic. 28, pp. 453-482.
6045 ! 2. powell, m.j.d. (1970), a hybrid method for non-linear equations,
6046 ! in numerical methods for non-linear equations, edited by
6047 ! p. rabinowitz, gordon and breach, london.
6051 ! coded by david m. gay.
6052 ! this subroutine was written in connection with research supported
6053 ! by the national science foundation under grants mcs-7600324 and
6056 !------------------------ external quantities ------------------------
6058 ! *** functions and subroutines called ***
6060 !el external dotprd, v2norm
6061 !el real(kind=8) :: dotprd, v2norm
6063 ! dotprd... returns inner product of two vectors.
6064 ! v2norm... returns 2-norm of a vector.
6066 ! *** intrinsic functions ***
6068 !el real(kind=8) :: dsqrt
6070 !-------------------------- local variables --------------------------
6073 real(kind=8) :: cfact, cnorm, ctrnwt, ghinvg, femnsq, gnorm,&
6074 nwtnrm, relax, rlambd, t, t1, t2
6075 !el real(kind=8) :: half, one, two, zero
6077 ! *** v subscripts ***
6079 !el integer bias, dgnorm, dstnrm, dst0, grdfac, gthg, gtstep,
6080 !el 1 nreduc, nwtfac, preduc, radius, stppar
6082 ! *** data initializations ***
6085 ! data half/0.5d+0/, one/1.d+0/, two/2.d+0/, zero/0.d+0/
6087 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, two=2.d+0, zero=0.d+0
6091 ! data bias/43/, dgnorm/1/, dstnrm/2/, dst0/3/, grdfac/45/,
6092 ! 1 gthg/44/, gtstep/4/, nreduc/6/, nwtfac/46/, preduc/7/,
6093 ! 2 radius/8/, stppar/5/
6095 integer,parameter :: bias=43, dgnorm=1, dstnrm=2, dst0=3, grdfac=45,&
6096 gthg=44, gtstep=4, nreduc=6, nwtfac=46, preduc=7,&
6100 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6104 if (nwtnrm .gt. zero) rlambd = v(radius) / nwtnrm
6106 ghinvg = two * v(nreduc)
6109 if (rlambd .lt. one) go to 30
6111 ! *** the newton step is inside the trust region ***
6116 v(preduc) = v(nreduc)
6119 20 step(i) = -nwtstp(i)
6122 30 v(dstnrm) = v(radius)
6123 cfact = (gnorm / v(gthg))**2
6124 ! *** cauchy step = -cfact * g.
6125 cnorm = gnorm * cfact
6126 relax = one - v(bias) * (one - gnorm*cnorm/ghinvg)
6127 if (rlambd .lt. relax) go to 50
6129 ! *** step is between relaxed newton and full newton steps ***
6131 v(stppar) = one - (rlambd - relax) / (one - relax)
6133 v(gtstep) = t * ghinvg
6134 v(preduc) = rlambd * (one - half*rlambd) * ghinvg
6137 40 step(i) = t * nwtstp(i)
6140 50 if (cnorm .lt. v(radius)) go to 70
6142 ! *** the cauchy step lies outside the trust region --
6143 ! *** step = scaled cauchy step ***
6145 t = -v(radius) / gnorm
6147 v(stppar) = one + cnorm / v(radius)
6148 v(gtstep) = -v(radius) * gnorm
6149 v(preduc) = v(radius)*(gnorm - half*v(radius)*(v(gthg)/gnorm)**2)
6151 60 step(i) = t * dig(i)
6154 ! *** compute dogleg step between cauchy and relaxed newton ***
6155 ! *** femur = relaxed newton step minus cauchy step ***
6157 70 ctrnwt = cfact * relax * ghinvg / gnorm
6158 ! *** ctrnwt = inner prod. of cauchy and relaxed newton steps,
6159 ! *** scaled by gnorm**-1.
6160 t1 = ctrnwt - gnorm*cfact**2
6161 ! *** t1 = inner prod. of femur and cauchy step, scaled by
6163 t2 = v(radius)*(v(radius)/gnorm) - gnorm*cfact**2
6165 femnsq = (t/gnorm)*t - ctrnwt - t1
6166 ! *** femnsq = square of 2-norm of femur, scaled by gnorm**-1.
6167 t = t2 / (t1 + dsqrt(t1**2 + femnsq*t2))
6168 ! *** dogleg step = cauchy step + t * femur.
6169 t1 = (t - one) * cfact
6174 v(gtstep) = t1*gnorm**2 + t2*ghinvg
6175 v(preduc) = -t1*gnorm * ((t2 + one)*gnorm) &
6176 - t2 * (one + half*t2)*ghinvg &
6177 - half * (v(gthg)*t1)**2
6179 80 step(i) = t1*dig(i) + t2*nwtstp(i)
6182 ! *** last line of dbdog follows ***
6183 end subroutine dbdog
6184 !-----------------------------------------------------------------------------
6185 subroutine ltvmul(n,x,l,y)
6187 ! *** compute x = (l**t)*y, where l is an n x n lower
6188 ! *** triangular matrix stored compactly by rows. x and y may
6189 ! *** occupy the same storage. ***
6192 !al real(kind=8) :: x(n), l(1), y(n)
6193 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
6194 ! dimension l(n*(n+1)/2)
6195 integer :: i, ij, i0, j
6196 real(kind=8) :: yi !el, zero
6200 real(kind=8),parameter :: zero=0.d+0
6209 x(j) = x(j) + yi*l(ij)
6214 ! *** last card of ltvmul follows ***
6215 end subroutine ltvmul
6216 !-----------------------------------------------------------------------------
6217 subroutine lupdat(beta,gamma,l,lambda,lplus,n,w,z)
6219 ! *** compute lplus = secant update of l ***
6221 ! *** parameter declarations ***
6224 !al double precision beta(n), gamma(n), l(1), lambda(n), lplus(1),
6225 real(kind=8) :: beta(n), gamma(n), l(n*(n+1)/2), lambda(n), &
6226 lplus(n*(n+1)/2),w(n), z(n)
6227 ! dimension l(n*(n+1)/2), lplus(n*(n+1)/2)
6229 !-------------------------- parameter usage --------------------------
6231 ! beta = scratch vector.
6232 ! gamma = scratch vector.
6233 ! l (input) lower triangular matrix, stored rowwise.
6234 ! lambda = scratch vector.
6235 ! lplus (output) lower triangular matrix, stored rowwise, which may
6236 ! occupy the same storage as l.
6237 ! n (input) length of vector parameters and order of matrices.
6238 ! w (input, destroyed on output) right singular vector of rank 1
6240 ! z (input, destroyed on output) left singular vector of rank 1
6243 !------------------------------- notes -------------------------------
6245 ! *** application and usage restrictions ***
6247 ! this routine updates the cholesky factor l of a symmetric
6248 ! positive definite matrix to which a secant update is being
6249 ! applied -- it computes a cholesky factor lplus of
6250 ! l * (i + z*w**t) * (i + w*z**t) * l**t. it is assumed that w
6251 ! and z have been chosen so that the updated matrix is strictly
6252 ! positive definite.
6254 ! *** algorithm notes ***
6256 ! this code uses recurrence 3 of ref. 1 (with d(j) = 1 for all j)
6257 ! to compute lplus of the form l * (i + z*w**t) * q, where q
6258 ! is an orthogonal matrix that makes the result lower triangular.
6259 ! lplus may have some negative diagonal elements.
6261 ! *** references ***
6263 ! 1. goldfarb, d. (1976), factorized variable metric methods for uncon-
6264 ! strained optimization, math. comput. 30, pp. 796-811.
6268 ! coded by david m. gay (fall 1979).
6269 ! this subroutine was written in connection with research supported
6270 ! by the national science foundation under grants mcs-7600324 and
6273 !------------------------ external quantities ------------------------
6275 ! *** intrinsic functions ***
6277 !el real(kind=8) :: dsqrt
6279 !-------------------------- local variables --------------------------
6281 integer :: i, ij, j, jj, jp1, k, nm1, np1
6282 real(kind=8) :: a, b, bj, eta, gj, lj, lij, ljj, nu, s, theta,&
6284 !el real(kind=8) :: one, zero
6286 ! *** data initializations ***
6289 ! data one/1.d+0/, zero/0.d+0/
6291 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
6294 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6298 if (n .le. 1) go to 30
6301 ! *** temporarily store s(j) = sum over k = j+1 to n of w(k)**2 in
6311 ! *** compute lambda, gamma, and beta by goldfarb*s recurrence 3.
6315 a = nu*z(j) - eta*wj
6318 lj = dsqrt(theta**2 + a*s)
6319 if (theta .gt. zero) lj = -lj
6322 gamma(j) = b * nu / lj
6323 beta(j) = (a - b*eta) / lj
6325 eta = -(eta + (a**2)/(theta - lj)) / lj
6327 30 lambda(n) = one + (nu*z(n) - eta*w(n))*w(n)
6329 ! *** update l, gradually overwriting w and z with l*w and l*z.
6332 jj = n * (n + 1) / 2
6337 lplus(jj) = lj * ljj
6342 if (k .eq. 1) go to 50
6349 lplus(ij) = lj*lij + bj*w(i) + gj*z(i)
6350 w(i) = w(i) + lij*wj
6351 z(i) = z(i) + lij*zj
6358 ! *** last card of lupdat follows ***
6359 end subroutine lupdat
6360 !-----------------------------------------------------------------------------
6361 subroutine lvmul(n,x,l,y)
6363 ! *** compute x = l*y, where l is an n x n lower triangular
6364 ! *** matrix stored compactly by rows. x and y may occupy the same
6368 !al double precision x(n), l(1), y(n)
6369 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
6370 ! dimension l(n*(n+1)/2)
6371 integer :: i, ii, ij, i0, j, np1
6372 real(kind=8) :: t !el, zero
6376 real(kind=8),parameter :: zero=0.d+0
6392 ! *** last card of lvmul follows ***
6393 end subroutine lvmul
6394 !-----------------------------------------------------------------------------
6395 subroutine vvmulp(n,x,y,z,k)
6397 ! *** set x(i) = y(i) * z(i)**k, 1 .le. i .le. n (for k = 1 or -1) ***
6400 real(kind=8) :: x(n), y(n), z(n)
6403 if (k .ge. 0) go to 20
6405 10 x(i) = y(i) / z(i)
6409 30 x(i) = y(i) * z(i)
6411 ! *** last card of vvmulp follows ***
6412 end subroutine vvmulp
6413 !-----------------------------------------------------------------------------
6414 subroutine wzbfgs(l,n,s,w,y,z)
6416 ! *** compute y and z for lupdat corresponding to bfgs update.
6419 !al double precision l(1), s(n), w(n), y(n), z(n)
6420 real(kind=8) :: l(n*(n+1)/2), s(n), w(n), y(n), z(n)
6421 ! dimension l(n*(n+1)/2)
6423 !-------------------------- parameter usage --------------------------
6425 ! l (i/o) cholesky factor of hessian, a lower triang. matrix stored
6426 ! compactly by rows.
6427 ! n (input) order of l and length of s, w, y, z.
6428 ! s (input) the step just taken.
6429 ! w (output) right singular vector of rank 1 correction to l.
6430 ! y (input) change in gradients corresponding to s.
6431 ! z (output) left singular vector of rank 1 correction to l.
6433 !------------------------------- notes -------------------------------
6435 ! *** algorithm notes ***
6437 ! when s is computed in certain ways, e.g. by gqtstp or
6438 ! dbldog, it is possible to save n**2/2 operations since (l**t)*s
6439 ! or l*(l**t)*s is then known.
6440 ! if the bfgs update to l*(l**t) would reduce its determinant to
6441 ! less than eps times its old value, then this routine in effect
6442 ! replaces y by theta*y + (1 - theta)*l*(l**t)*s, where theta
6443 ! (between 0 and 1) is chosen to make the reduction factor = eps.
6447 ! coded by david m. gay (fall 1979).
6448 ! this subroutine was written in connection with research supported
6449 ! by the national science foundation under grants mcs-7600324 and
6452 !------------------------ external quantities ------------------------
6454 ! *** functions and subroutines called ***
6456 !el external dotprd, livmul, ltvmul
6457 !el real(kind=8) :: dotprd
6458 ! dotprd returns inner product of two vectors.
6459 ! livmul multiplies l**-1 times a vector.
6460 ! ltvmul multiplies l**t times a vector.
6462 ! *** intrinsic functions ***
6464 !el real(kind=8) :: dsqrt
6466 !-------------------------- local variables --------------------------
6469 real(kind=8) :: cs, cy, epsrt, shs, ys, theta !el, eps, one
6471 ! *** data initializations ***
6474 ! data eps/0.1d+0/, one/1.d+0/
6476 real(kind=8),parameter :: eps=0.1d+0, one=1.d+0
6479 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6481 call ltvmul(n, w, l, s)
6482 shs = dotprd(n, w, w)
6483 ys = dotprd(n, y, s)
6484 if (ys .ge. eps*shs) go to 10
6485 theta = (one - eps) * shs / (shs - ys)
6487 cy = theta / (shs * epsrt)
6488 cs = (one + (theta-one)/epsrt) / shs
6490 10 cy = one / (dsqrt(ys) * dsqrt(shs))
6492 20 call livmul(n, z, l, y)
6494 30 z(i) = cy * z(i) - cs * w(i)
6497 ! *** last card of wzbfgs follows ***
6498 end subroutine wzbfgs
6499 !-----------------------------------------------------------------------------
6500 !-----------------------------------------------------------------------------