2 !-----------------------------------------------------------------------------
16 ! double precision funcgrad_dc
18 !-----------------------------------------------------------------------------
21 !-----------------------------------------------------------------------------
23 !-----------------------------------------------------------------------------
25 !-----------------------------------------------------------------------------
26 subroutine assst(iv, liv, lv, v)
28 ! *** assess candidate step (***sol version 2.3) ***
37 ! this subroutine is called by an unconstrained minimization
38 ! routine to assess the next candidate step. it may recommend one
39 ! of several courses of action, such as accepting the step, recom-
40 ! puting it using the same or a new quadratic model, or halting due
41 ! to convergence or false convergence. see the return code listing
44 !-------------------------- parameter usage --------------------------
46 ! iv (i/o) integer parameter and scratch vector -- see description
47 ! below of iv values referenced.
48 ! liv (in) length of iv array.
49 ! lv (in) length of v array.
50 ! v (i/o) real parameter and scratch vector -- see description
51 ! below of v values referenced.
53 ! *** iv values referenced ***
55 ! iv(irc) (i/o) on input for the first step tried in a new iteration,
56 ! iv(irc) should be set to 3 or 4 (the value to which it is
57 ! set when step is definitely to be accepted). on input
58 ! after step has been recomputed, iv(irc) should be
59 ! unchanged since the previous return of assst.
60 ! on output, iv(irc) is a return code having one of the
62 ! 1 = switch models or try smaller step.
63 ! 2 = switch models or accept step.
64 ! 3 = accept step and determine v(radfac) by gradient
66 ! 4 = accept step, v(radfac) has been determined.
67 ! 5 = recompute step (using the same model).
68 ! 6 = recompute step with radius = v(lmaxs) but do not
69 ! evaulate the objective function.
70 ! 7 = x-convergence (see v(xctol)).
71 ! 8 = relative function convergence (see v(rfctol)).
72 ! 9 = both x- and relative function convergence.
73 ! 10 = absolute function convergence (see v(afctol)).
74 ! 11 = singular convergence (see v(lmaxs)).
75 ! 12 = false convergence (see v(xftol)).
76 ! 13 = iv(irc) was out of range on input.
77 ! return code i has precdence over i+1 for i = 9, 10, 11.
78 ! iv(mlstgd) (i/o) saved value of iv(model).
79 ! iv(model) (i/o) on input, iv(model) should be an integer identifying
80 ! the current quadratic model of the objective function.
81 ! if a previous step yielded a better function reduction,
82 ! then iv(model) will be set to iv(mlstgd) on output.
83 ! iv(nfcall) (in) invocation count for the objective function.
84 ! iv(nfgcal) (i/o) value of iv(nfcall) at step that gave the biggest
85 ! function reduction this iteration. iv(nfgcal) remains
86 ! unchanged until a function reduction is obtained.
87 ! iv(radinc) (i/o) the number of radius increases (or minus the number
88 ! of decreases) so far this iteration.
89 ! iv(restor) (out) set to 1 if v(f) has been restored and x should be
90 ! restored to its initial value, to 2 if x should be saved,
91 ! to 3 if x should be restored from the saved value, and to
93 ! iv(stage) (i/o) count of the number of models tried so far in the
95 ! iv(stglim) (in) maximum number of models to consider.
96 ! iv(switch) (out) set to 0 unless a new model is being tried and it
97 ! gives a smaller function value than the previous model,
98 ! in which case assst sets iv(switch) = 1.
99 ! iv(toobig) (in) is nonzero if step was too big (e.g. if it caused
101 ! iv(xirc) (i/o) value that iv(irc) would have in the absence of
102 ! convergence, false convergence, and oversized steps.
104 ! *** v values referenced ***
106 ! v(afctol) (in) absolute function convergence tolerance. if the
107 ! absolute value of the current function value v(f) is less
108 ! than v(afctol), then assst returns with iv(irc) = 10.
109 ! v(decfac) (in) factor by which to decrease radius when iv(toobig) is
111 ! v(dstnrm) (in) the 2-norm of d*step.
112 ! v(dstsav) (i/o) value of v(dstnrm) on saved step.
113 ! v(dst0) (in) the 2-norm of d times the newton step (when defined,
114 ! i.e., for v(nreduc) .ge. 0).
115 ! v(f) (i/o) on both input and output, v(f) is the objective func-
116 ! tion value at x. if x is restored to a previous value,
117 ! then v(f) is restored to the corresponding value.
118 ! v(fdif) (out) the function reduction v(f0) - v(f) (for the output
119 ! value of v(f) if an earlier step gave a bigger function
120 ! decrease, and for the input value of v(f) otherwise).
121 ! v(flstgd) (i/o) saved value of v(f).
122 ! v(f0) (in) objective function value at start of iteration.
123 ! v(gtslst) (i/o) value of v(gtstep) on saved step.
124 ! v(gtstep) (in) inner product between step and gradient.
125 ! v(incfac) (in) minimum factor by which to increase radius.
126 ! v(lmaxs) (in) maximum reasonable step size (and initial step bound).
127 ! if the actual function decrease is no more than twice
128 ! what was predicted, if a return with iv(irc) = 7, 8, 9,
129 ! or 10 does not occur, if v(dstnrm) .gt. v(lmaxs), and if
130 ! v(preduc) .le. v(sctol) * abs(v(f0)), then assst re-
131 ! turns with iv(irc) = 11. if so doing appears worthwhile,
132 ! then assst repeats this test with v(preduc) computed for
133 ! a step of length v(lmaxs) (by a return with iv(irc) = 6).
134 ! v(nreduc) (i/o) function reduction predicted by quadratic model for
135 ! newton step. if assst is called with iv(irc) = 6, i.e.,
136 ! if v(preduc) has been computed with radius = v(lmaxs) for
137 ! use in the singular convervence test, then v(nreduc) is
138 ! set to -v(preduc) before the latter is restored.
139 ! v(plstgd) (i/o) value of v(preduc) on saved step.
140 ! v(preduc) (i/o) function reduction predicted by quadratic model for
142 ! v(radfac) (out) factor to be used in determining the new radius,
143 ! which should be v(radfac)*dst, where dst is either the
144 ! output value of v(dstnrm) or the 2-norm of
145 ! diag(newd)*step for the output value of step and the
146 ! updated version, newd, of the scale vector d. for
147 ! iv(irc) = 3, v(radfac) = 1.0 is returned.
148 ! v(rdfcmn) (in) minimum value for v(radfac) in terms of the input
149 ! value of v(dstnrm) -- suggested value = 0.1.
150 ! v(rdfcmx) (in) maximum value for v(radfac) -- suggested value = 4.0.
151 ! v(reldx) (in) scaled relative change in x caused by step, computed
152 ! (e.g.) by function reldst as
153 ! max (d(i)*abs(x(i)-x0(i)), 1 .le. i .le. p) /
154 ! max (d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p).
155 ! v(rfctol) (in) relative function convergence tolerance. if the
156 ! actual function reduction is at most twice what was pre-
157 ! dicted and v(nreduc) .le. v(rfctol)*abs(v(f0)), then
158 ! assst returns with iv(irc) = 8 or 9.
159 ! v(stppar) (in) marquardt parameter -- 0 means full newton step.
160 ! v(tuner1) (in) tuning constant used to decide if the function
161 ! reduction was much less than expected. suggested
163 ! v(tuner2) (in) tuning constant used to decide if the function
164 ! reduction was large enough to accept step. suggested
166 ! v(tuner3) (in) tuning constant used to decide if the radius
167 ! should be increased. suggested value = 0.75.
168 ! v(xctol) (in) x-convergence criterion. if step is a newton step
169 ! (v(stppar) = 0) having v(reldx) .le. v(xctol) and giving
170 ! at most twice the predicted function decrease, then
171 ! assst returns iv(irc) = 7 or 9.
172 ! v(xftol) (in) false convergence tolerance. if step gave no or only
173 ! a small function decrease and v(reldx) .le. v(xftol),
174 ! then assst returns with iv(irc) = 12.
176 !------------------------------- notes -------------------------------
178 ! *** application and usage restrictions ***
180 ! this routine is called as part of the nl2sol (nonlinear
181 ! least-squares) package. it may be used in any unconstrained
182 ! minimization solver that uses dogleg, goldfeld-quandt-trotter,
183 ! or levenberg-marquardt steps.
185 ! *** algorithm notes ***
187 ! see (1) for further discussion of the assessing and model
188 ! switching strategies. while nl2sol considers only two models,
189 ! assst is designed to handle any number of models.
191 ! *** usage notes ***
193 ! on the first call of an iteration, only the i/o variables
194 ! step, x, iv(irc), iv(model), v(f), v(dstnrm), v(gtstep), and
195 ! v(preduc) need have been initialized. between calls, no i/o
196 ! values execpt step, x, iv(model), v(f) and the stopping toler-
197 ! ances should be changed.
198 ! after a return for convergence or false convergence, one can
199 ! change the stopping tolerances and call assst again, in which
200 ! case the stopping tests will be repeated.
204 ! (1) dennis, j.e., jr., gay, d.m., and welsch, r.e. (1981),
205 ! an adaptive nonlinear least-squares algorithm,
206 ! acm trans. math. software, vol. 7, no. 3.
208 ! (2) powell, m.j.d. (1970) a fortran subroutine for solving
209 ! systems of nonlinear algebraic equations, in numerical
210 ! methods for nonlinear algebraic equations, edited by
211 ! p. rabinowitz, gordon and breach, london.
215 ! john dennis designed much of this routine, starting with
216 ! ideas in (2). roy welsch suggested the model switching strategy.
217 ! david gay and stephen peters cast this subroutine into a more
218 ! portable form (winter 1977), and david gay cast it into its
219 ! present form (fall 1978).
223 ! this subroutine was written in connection with research
224 ! supported by the national science foundation under grants
225 ! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
228 !------------------------ external quantities ------------------------
230 ! *** no external functions and subroutines ***
232 ! *** intrinsic functions ***
234 !el real(kind=8) :: dabs, dmax1
236 ! *** no common blocks ***
238 !-------------------------- local variables --------------------------
242 real(kind=8) :: emax, emaxs, gts, rfac1, xmax
243 !el real(kind=8) :: half, one, onep2, two, zero
245 ! *** subscripts for iv and v ***
247 !el integer :: afctol, decfac, dstnrm, dstsav, dst0, f, fdif, flstgd, f0,&
248 !el gtslst, gtstep, incfac, irc, lmaxs, mlstgd, model, nfcall,&
249 !el nfgcal, nreduc, plstgd, preduc, radfac, radinc, rdfcmn,&
250 !el rdfcmx, reldx, restor, rfctol, sctol, stage, stglim,&
251 !el stppar, switch, toobig, tuner1, tuner2, tuner3, xctol,&
255 ! *** data initializations ***
258 ! data half/0.5d+0/, one/1.d+0/, onep2/1.2d+0/, two/2.d+0/,
261 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, onep2=1.2d+0, two=2.d+0,&
266 ! data irc/29/, mlstgd/32/, model/5/, nfcall/6/, nfgcal/7/,
267 ! 1 radinc/8/, restor/9/, stage/10/, stglim/11/, switch/12/,
268 ! 2 toobig/2/, xirc/13/
270 integer,parameter :: irc=29, mlstgd=32, model=5, nfcall=6, nfgcal=7,&
271 radinc=8, restor=9, stage=10, stglim=11, switch=12,&
275 ! data afctol/31/, decfac/22/, dstnrm/2/, dst0/3/, dstsav/18/,
276 ! 1 f/10/, fdif/11/, flstgd/12/, f0/13/, gtslst/14/, gtstep/4/,
277 ! 2 incfac/23/, lmaxs/36/, nreduc/6/, plstgd/15/, preduc/7/,
278 ! 3 radfac/16/, rdfcmn/24/, rdfcmx/25/, reldx/17/, rfctol/32/,
279 ! 4 sctol/37/, stppar/5/, tuner1/26/, tuner2/27/, tuner3/28/,
280 ! 5 xctol/33/, xftol/34/
282 integer,parameter :: afctol=31, decfac=22, dstnrm=2, dst0=3, dstsav=18,&
283 f=10, fdif=11, flstgd=12, f0=13, gtslst=14, gtstep=4,&
284 incfac=23, lmaxs=36, nreduc=6, plstgd=15, preduc=7,&
285 radfac=16, rdfcmn=24, rdfcmx=25, reldx=17, rfctol=32,&
286 sctol=37, stppar=5, tuner1=26, tuner2=27, tuner3=28,&
290 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
298 if (i .ge. 1 .and. i .le. 12) &
299 go to (20,30,10,10,40,280,220,220,220,220,220,170), i
303 ! *** initialize for new iteration ***
308 if (iv(toobig) .eq. 0) go to 110
313 ! *** step was recomputed with new model or smaller radius ***
314 ! *** first decide which ***
316 20 if (iv(model) .ne. iv(mlstgd)) go to 30
317 ! *** old model retained, smaller radius tried ***
318 ! *** do not consider any more new models this iteration ***
319 iv(stage) = iv(stglim)
323 ! *** a new model is being tried. decide whether to keep it. ***
325 30 iv(stage) = iv(stage) + 1
327 ! *** now we add the possibiltiy that step was recomputed with ***
328 ! *** the same model, perhaps because of an oversized step. ***
330 40 if (iv(stage) .gt. 0) go to 50
332 ! *** step was recomputed because it was too big. ***
334 if (iv(toobig) .ne. 0) go to 60
336 ! *** restore iv(stage) and pick up where we left off. ***
338 iv(stage) = -iv(stage)
340 go to (20, 30, 110, 110, 70), i
342 50 if (iv(toobig) .eq. 0) go to 70
344 ! *** handle oversize step ***
346 if (iv(radinc) .gt. 0) go to 80
347 iv(stage) = -iv(stage)
350 60 v(radfac) = v(decfac)
351 iv(radinc) = iv(radinc) - 1
356 70 if (v(f) .lt. v(flstgd)) go to 110
358 ! *** the new step is a loser. restore old model. ***
360 if (iv(model) .eq. iv(mlstgd)) go to 80
361 iv(model) = iv(mlstgd)
364 ! *** restore step, etc. only if a previous step decreased v(f).
366 80 if (v(flstgd) .ge. v(f0)) go to 110
369 v(preduc) = v(plstgd)
370 v(gtstep) = v(gtslst)
371 if (iv(switch) .eq. 0) rfac1 = v(dstnrm) / v(dstsav)
372 v(dstnrm) = v(dstsav)
376 110 v(fdif) = v(f0) - v(f)
377 if (v(fdif) .gt. v(tuner2) * v(preduc)) go to 140
378 if(iv(radinc).gt.0) go to 140
380 ! *** no (or only a trivial) function decrease
381 ! *** -- so try new model or smaller radius
383 if (v(f) .lt. v(f0)) go to 120
384 iv(mlstgd) = iv(model)
391 if (iv(stage) .lt. iv(stglim)) go to 160
393 iv(radinc) = iv(radinc) - 1
396 ! *** nontrivial function decrease achieved ***
400 v(dstsav) = v(dstnrm)
401 if (v(fdif) .gt. v(preduc)*v(tuner1)) go to 190
403 ! *** decrease was much less than predicted -- either change models
404 ! *** or accept step with decreased radius.
406 if (iv(stage) .ge. iv(stglim)) go to 150
407 ! *** consider switching models ***
411 ! *** accept step with decreased radius ***
415 ! *** set v(radfac) to fletcher*s decrease factor ***
417 160 iv(xirc) = iv(irc)
418 emax = v(gtstep) + v(fdif)
419 v(radfac) = half * rfac1
420 if (emax .lt. v(gtstep)) v(radfac) = rfac1 * dmax1(v(rdfcmn),&
421 half * v(gtstep)/emax)
423 ! *** do false convergence test ***
425 170 if (v(reldx) .le. v(xftol)) go to 180
427 if (v(f) .lt. v(f0)) go to 200
433 ! *** handle good function decrease ***
435 190 if (v(fdif) .lt. (-v(tuner3) * v(gtstep))) go to 210
437 ! *** increasing radius looks worthwhile. see if we just
438 ! *** recomputed step with a decreased radius or restored step
439 ! *** after recomputing it with a larger radius.
441 if (iv(radinc) .lt. 0) go to 210
442 if (iv(restor) .eq. 1) go to 210
444 ! *** we did not. try a longer step unless this was a newton
447 v(radfac) = v(rdfcmx)
449 if (v(fdif) .lt. (half/v(radfac) - one) * gts) &
450 v(radfac) = dmax1(v(incfac), half*gts/(gts + v(fdif)))
452 if (v(stppar) .eq. zero) go to 230
453 if (v(dst0) .ge. zero .and. (v(dst0) .lt. two*v(dstnrm) &
454 .or. v(nreduc) .lt. onep2*v(fdif))) go to 230
455 ! *** step was not a newton step. recompute it with
456 ! *** a larger radius.
458 iv(radinc) = iv(radinc) + 1
460 ! *** save values corresponding to good step ***
463 iv(mlstgd) = iv(model)
464 if (iv(restor) .ne. 1) iv(restor) = 2
465 v(dstsav) = v(dstnrm)
467 v(plstgd) = v(preduc)
468 v(gtslst) = v(gtstep)
471 ! *** accept step with radius unchanged ***
477 ! *** come here for a restart after convergence ***
479 220 iv(irc) = iv(xirc)
480 if (v(dstsav) .ge. zero) go to 240
484 ! *** perform convergence tests ***
486 230 iv(xirc) = iv(irc)
487 240 if (iv(restor) .eq. 1 .and. v(flstgd) .lt. v(f0)) iv(restor) = 3
488 if (half * v(fdif) .gt. v(preduc)) go to 999
489 emax = v(rfctol) * dabs(v(f0))
490 emaxs = v(sctol) * dabs(v(f0))
491 if (v(dstnrm) .gt. v(lmaxs) .and. v(preduc) .le. emaxs) &
493 if (v(dst0) .lt. zero) go to 250
495 if ((v(nreduc) .gt. zero .and. v(nreduc) .le. emax) .or. &
496 (v(nreduc) .eq. zero .and. v(preduc) .eq. zero)) i = 2
497 if (v(stppar) .eq. zero .and. v(reldx) .le. v(xctol) &
498 .and. goodx) i = i + 1
499 if (i .gt. 0) iv(irc) = i + 6
501 ! *** consider recomputing step of length v(lmaxs) for singular
502 ! *** convergence test.
504 250 if (iv(irc) .gt. 5 .and. iv(irc) .ne. 12) go to 999
505 if (v(dstnrm) .gt. v(lmaxs)) go to 260
506 if (v(preduc) .ge. emaxs) go to 999
507 if (v(dst0) .le. zero) go to 270
508 if (half * v(dst0) .le. v(lmaxs)) go to 999
510 260 if (half * v(dstnrm) .le. v(lmaxs)) go to 999
511 xmax = v(lmaxs) / v(dstnrm)
512 if (xmax * (two - xmax) * v(preduc) .ge. emaxs) go to 999
513 270 if (v(nreduc) .lt. zero) go to 290
515 ! *** recompute v(preduc) for use in singular convergence test ***
517 v(gtslst) = v(gtstep)
518 v(dstsav) = v(dstnrm)
519 if (iv(irc) .eq. 12) v(dstsav) = -v(dstsav)
520 v(plstgd) = v(preduc)
523 if (i .eq. 3) iv(restor) = 0
527 ! *** perform singular convergence test with recomputed v(preduc) ***
529 280 v(gtstep) = v(gtslst)
530 v(dstnrm) = dabs(v(dstsav))
532 if (v(dstsav) .le. zero) iv(irc) = 12
533 v(nreduc) = -v(preduc)
534 v(preduc) = v(plstgd)
536 290 if (-v(nreduc) .le. v(sctol) * dabs(v(f0))) iv(irc) = 11
540 ! *** last card of assst follows ***
542 !-----------------------------------------------------------------------------
543 subroutine deflt(alg, iv, liv, lv, v)
545 ! *** supply ***sol (version 2.3) default values to iv and v ***
547 ! *** alg = 1 means regression constants.
548 ! *** alg = 2 means general unconstrained optimization constants.
554 real(kind=8) :: v(lv)
556 !el external imdcon, vdflt
558 ! imdcon... returns machine-dependent integer constants.
559 ! vdflt.... provides default values to v.
562 integer :: miniv(2), minv(2)
564 ! *** subscripts for iv ***
566 !el integer algsav, covprt, covreq, dtype, hc, ierr, inith, inits,
567 !el 1 ipivot, ivneed, lastiv, lastv, lmat, mxfcal, mxiter,
568 !el 2 nfcov, ngcov, nvdflt, outlev, parprt, parsav, perm,
569 !el 3 prunit, qrtyp, rdreq, rmat, solprt, statpr, vneed,
572 ! *** iv subscript values ***
575 ! data algsav/51/, covprt/14/, covreq/15/, dtype/16/, hc/71/,
576 ! 1 ierr/75/, inith/25/, inits/25/, ipivot/76/, ivneed/3/,
577 ! 2 lastiv/44/, lastv/45/, lmat/42/, mxfcal/17/, mxiter/18/,
578 ! 3 nfcov/52/, ngcov/53/, nvdflt/50/, outlev/19/, parprt/20/,
579 ! 4 parsav/49/, perm/58/, prunit/21/, qrtyp/80/, rdreq/57/,
580 ! 5 rmat/78/, solprt/22/, statpr/23/, vneed/4/, vsave/60/,
583 integer,parameter :: algsav=51, covprt=14, covreq=15, dtype=16, hc=71,&
584 ierr=75, inith=25, inits=25, ipivot=76, ivneed=3,&
585 lastiv=44, lastv=45, lmat=42, mxfcal=17, mxiter=18,&
586 nfcov=52, ngcov=53, nvdflt=50, outlev=19, parprt=20,&
587 parsav=49, perm=58, prunit=21, qrtyp=80, rdreq=57,&
588 rmat=78, solprt=22, statpr=23, vneed=4, vsave=60,&
591 data miniv(1)/80/, miniv(2)/59/, minv(1)/98/, minv(2)/71/
595 !------------------------------- body --------------------------------
597 if (alg .lt. 1 .or. alg .gt. 2) go to 40
599 if (liv .lt. miv) go to 20
601 if (lv .lt. mv) go to 30
602 call vdflt(alg, lv, v)
614 iv(prunit) = imdcon(1)
620 if (alg .ge. 2) go to 10
622 ! *** regression values
639 ! *** general optimization values
658 ! *** last card of deflt follows ***
660 !-----------------------------------------------------------------------------
661 real(kind=8) function dotprd(p,x,y)
663 ! *** return the inner product of the p-vectors x and y. ***
666 real(kind=8) :: x(p), y(p)
669 !el real(kind=8) :: one, zero
670 real(kind=8) :: sqteta, t
672 !el real(kind=8) :: dmax1, dabs
675 !el real(kind=8) :: rmdcon
677 ! *** rmdcon(2) returns a machine-dependent constant, sqteta, which
678 ! *** is slightly larger than the smallest positive number that
679 ! *** can be squared without underflowing.
682 ! data one/1.d+0/, sqteta/0.d+0/, zero/0.d+0/
684 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
689 if (p .le. 0) go to 999
690 !rc if (sqteta .eq. zero) sqteta = rmdcon(2)
692 !rc t = dmax1(dabs(x(i)), dabs(y(i)))
693 !rc if (t .gt. one) go to 10
694 !rc if (t .lt. sqteta) go to 20
695 !rc t = (x(i)/sqteta)*y(i)
696 !rc if (dabs(t) .lt. sqteta) go to 20
697 10 dotprd = dotprd + x(i)*y(i)
701 ! *** last card of dotprd follows ***
703 !-----------------------------------------------------------------------------
704 subroutine itsum(d, g, iv, liv, lv, p, v, x)
706 ! *** print iteration summary for ***sol (version 2.3) ***
708 ! *** parameter declarations ***
713 real(kind=8) :: d(p), g(p), v(lv), x(p)
715 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
717 ! *** local variables ***
719 integer :: alg, i, iv1, m, nf, ng, ol, pu
721 ! real model1(6), model2(6)
723 character(len=4) :: model1(6), model2(6)
725 real(kind=8) :: nreldf, oldf, preldf, reldf !el, zero
727 ! *** intrinsic functions ***
730 !el real(kind=8) :: dabs, dmax1
732 ! *** no external functions or subroutines ***
734 ! *** subscripts for iv and v ***
736 !el integer algsav, dstnrm, f, fdif, f0, needhd, nfcall, nfcov, ngcov,
737 !el 1 ngcall, niter, nreduc, outlev, preduc, prntit, prunit,
738 !el 2 reldx, solprt, statpr, stppar, sused, x0prt
740 ! *** iv subscript values ***
743 ! data algsav/51/, needhd/36/, nfcall/6/, nfcov/52/, ngcall/30/,
744 ! 1 ngcov/53/, niter/31/, outlev/19/, prntit/39/, prunit/21/,
745 ! 2 solprt/22/, statpr/23/, sused/64/, x0prt/24/
747 integer,parameter :: algsav=51, needhd=36, nfcall=6, nfcov=52, ngcall=30,&
748 ngcov=53, niter=31, outlev=19, prntit=39, prunit=21,&
749 solprt=22, statpr=23, sused=64, x0prt=24
752 ! *** v subscript values ***
755 ! data dstnrm/2/, f/10/, f0/13/, fdif/11/, nreduc/6/, preduc/7/,
756 ! 1 reldx/17/, stppar/5/
758 integer,parameter :: dstnrm=2, f=10, f0=13, fdif=11, nreduc=6, preduc=7,&
765 real(kind=8),parameter :: zero=0.d+0
768 ! data model1(1)/4h /, model1(2)/4h /, model1(3)/4h /,
769 ! 1 model1(4)/4h /, model1(5)/4h g /, model1(6)/4h s /,
770 ! 2 model2(1)/4h g /, model2(2)/4h s /, model2(3)/4hg-s /,
771 ! 3 model2(4)/4hs-g /, model2(5)/4h-s-g/, model2(6)/4h-g-s/
773 data model1/' ',' ',' ',' ',' g ',' s '/,&
774 model2/' g ',' s ','g-s ','s-g ','-s-g','-g-s'/
777 !------------------------------- body --------------------------------
780 if (pu .eq. 0) go to 999
782 if (iv1 .gt. 62) iv1 = iv1 - 51
785 if (iv1 .lt. 2 .or. iv1 .gt. 15) go to 370
786 if (iv1 .ge. 12) go to 120
787 if (iv1 .eq. 2 .and. iv(niter) .eq. 0) go to 390
788 if (ol .eq. 0) go to 120
789 if (iv1 .ge. 10 .and. iv(prntit) .eq. 0) go to 120
790 if (iv1 .gt. 2) go to 10
791 iv(prntit) = iv(prntit) + 1
792 if (iv(prntit) .lt. iabs(ol)) go to 999
793 10 nf = iv(nfcall) - iabs(iv(nfcov))
797 oldf = dmax1(dabs(v(f0)), dabs(v(f)))
798 if (oldf .le. zero) go to 20
799 reldf = v(fdif) / oldf
800 preldf = v(preduc) / oldf
801 20 if (ol .gt. 0) go to 60
803 ! *** print short summary line ***
805 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,30)
806 30 format(/10h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,&
808 if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,40)
809 40 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,&
812 if (alg .eq. 2) go to 50
814 write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
815 model1(m), model2(m), v(stppar)
818 50 write(pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
822 ! *** print long summary line ***
824 60 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,70)
825 70 format(/11h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,&
826 2x,13hmodel stppar,2x,6hd*step,2x,7hnpreldf)
827 if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,80)
828 80 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,&
829 3x,6hstppar,3x,6hd*step,3x,7hnpreldf)
832 if (oldf .gt. zero) nreldf = v(nreduc) / oldf
833 if (alg .eq. 2) go to 90
835 write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
836 model1(m), model2(m), v(stppar), v(dstnrm), nreldf
839 90 write(pu,110) iv(niter), nf, v(f), reldf, preldf,&
840 v(reldx), v(stppar), v(dstnrm), nreldf
841 100 format(i6,i5,d10.3,2d9.2,d8.1,a3,a4,2d8.1,d9.2)
842 110 format(i6,i5,d11.3,2d10.2,3d9.1,d10.2)
844 120 if (iv(statpr) .lt. 0) go to 430
845 go to (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,&
849 140 format(/26h ***** x-convergence *****)
853 160 format(/42h ***** relative function convergence *****)
857 180 format(/49h ***** x- and relative function convergence *****)
861 200 format(/42h ***** absolute function convergence *****)
865 220 format(/33h ***** singular convergence *****)
869 240 format(/30h ***** false convergence *****)
873 260 format(/38h ***** function evaluation limit *****)
877 280 format(/28h ***** iteration limit *****)
881 300 format(/18h ***** stopx *****)
885 320 format(/44h ***** initial f(x) cannot be computed *****)
890 340 format(/37h ***** bad parameters to assess *****)
894 360 format(/43h ***** gradient could not be computed *****)
895 if (iv(niter) .gt. 0) go to 480
898 370 write(pu,380) iv(1)
899 380 format(/14h ***** iv(1) =,i5,6h *****)
902 ! *** initial call on itsum ***
904 390 if (iv(x0prt) .ne. 0) write(pu,400) (i, x(i), d(i), i = 1, p)
905 400 format(/23h i initial x(i),8x,4hd(i)//(1x,i5,d17.6,d14.3))
906 ! *** the following are to avoid undefined variables when the
907 ! *** function evaluation limit is 1...
913 if (iv1 .ge. 12) go to 999
916 if (ol .eq. 0) go to 999
917 if (ol .lt. 0 .and. alg .eq. 1) write(pu,30)
918 if (ol .lt. 0 .and. alg .eq. 2) write(pu,40)
919 if (ol .gt. 0 .and. alg .eq. 1) write(pu,70)
920 if (ol .gt. 0 .and. alg .eq. 2) write(pu,80)
921 if (alg .eq. 1) write(pu,410) v(f)
922 if (alg .eq. 2) write(pu,420) v(f)
923 410 format(/11h 0 1,d10.3)
924 !365 format(/11h 0 1,e11.3)
925 420 format(/11h 0 1,d11.3)
928 ! *** print various information requested on solution ***
931 if (iv(statpr) .eq. 0) go to 480
932 oldf = dmax1(dabs(v(f0)), dabs(v(f)))
935 if (oldf .le. zero) go to 440
936 preldf = v(preduc) / oldf
937 nreldf = v(nreduc) / oldf
938 440 nf = iv(nfcall) - iv(nfcov)
939 ng = iv(ngcall) - iv(ngcov)
940 write(pu,450) v(f), v(reldx), nf, ng, preldf, nreldf
941 450 format(/9h function,d17.6,8h reldx,d17.3/12h func. evals,&
942 i8,9x,11hgrad. evals,i8/7h preldf,d16.3,6x,7hnpreldf,d15.3)
944 if (iv(nfcov) .gt. 0) write(pu,460) iv(nfcov)
945 460 format(/1x,i4,50h extra func. evals for covariance and diagnostics.)
946 if (iv(ngcov) .gt. 0) write(pu,470) iv(ngcov)
947 470 format(1x,i4,50h extra grad. evals for covariance and diagnostics.)
949 480 if (iv(solprt) .eq. 0) go to 999
952 490 format(/22h i final x(i),8x,4hd(i),10x,4hg(i)/)
954 write(pu,510) i, x(i), d(i), g(i)
956 510 format(1x,i5,d16.6,2d14.3)
960 530 format(/24h inconsistent dimensions)
962 ! *** last card of itsum follows ***
964 !-----------------------------------------------------------------------------
965 subroutine litvmu(n, x, l, y)
967 ! *** solve (l**t)*x = y, where l is an n x n lower triangular
968 ! *** matrix stored compactly by rows. x and y may occupy the same
972 !al real(kind=8) :: x(n), l(1), y(n)
973 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
974 integer :: i, ii, ij, im1, i0, j, np1
975 real(kind=8) :: xi !el, zero
979 real(kind=8),parameter :: zero=0.d+0
990 if (i .le. 1) go to 999
992 if (xi .eq. zero) go to 30
996 x(j) = x(j) - xi*l(ij)
1000 ! *** last card of litvmu follows ***
1001 end subroutine litvmu
1002 !-----------------------------------------------------------------------------
1003 subroutine livmul(n, x, l, y)
1005 ! *** solve l*x = y, where l is an n x n lower triangular
1006 ! *** matrix stored compactly by rows. x and y may occupy the same
1010 !al real(kind=8) :: x(n), l(1), y(n)
1011 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
1013 !el real(kind=8) :: dotprd
1015 real(kind=8) :: t !el, zero
1019 real(kind=8),parameter :: zero=0.d+0
1023 if (y(k) .ne. zero) go to 20
1029 if (k .ge. n) go to 999
1032 t = dotprd(i-1, l(j+1), x)
1034 x(i) = (y(i) - t)/l(j)
1037 ! *** last card of livmul follows ***
1038 end subroutine livmul
1039 !-----------------------------------------------------------------------------
1040 subroutine parck(alg, d, iv, liv, lv, n, v)
1042 ! *** check ***sol (version 2.3) parameters, print changed values ***
1044 ! *** alg = 1 for regression, alg = 2 for general unconstrained opt.
1046 integer :: alg, liv, n
1049 real(kind=8) :: d(n), v(lv)
1051 !el external rmdcon, vcopy, vdflt
1052 !el real(kind=8) :: rmdcon
1053 ! rmdcon -- returns machine-dependent constants.
1054 ! vcopy -- copies one vector to another.
1055 ! vdflt -- supplies default parameter values to v alone.
1060 ! *** local variables ***
1062 integer :: i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu
1063 integer :: ijmp, jlim(2), miniv(2), ndflt(2)
1065 ! integer varnm(2), sh(2)
1066 ! real cngd(3), dflt(3), vn(2,34), which(3)
1068 character(len=1) :: varnm(2), sh(2)
1069 character(len=4) :: cngd(3), dflt(3), vn(2,34), which(3)
1071 real(kind=8) :: big, machep, tiny, vk, vm(34), vx(34), zero
1073 ! *** iv and v subscripts ***
1075 !el integer algsav, dinit, dtype, dtype0, epslon, inits, ivneed,
1076 !el 1 lastiv, lastv, lmat, nextiv, nextv, nvdflt, oldn,
1077 !el 2 parprt, parsav, perm, prunit, vneed
1081 ! data algsav/51/, dinit/38/, dtype/16/, dtype0/54/, epslon/19/,
1082 ! 1 inits/25/, ivneed/3/, lastiv/44/, lastv/45/, lmat/42/,
1083 ! 2 nextiv/46/, nextv/47/, nvdflt/50/, oldn/38/, parprt/20/,
1084 ! 3 parsav/49/, perm/58/, prunit/21/, vneed/4/
1086 integer,parameter :: algsav=51, dinit=38, dtype=16, dtype0=54, epslon=19,&
1087 inits=25, ivneed=3, lastiv=44, lastv=45, lmat=42,&
1088 nextiv=46, nextv=47, nvdflt=50, oldn=38, parprt=20,&
1089 parsav=49, perm=58, prunit=21, vneed=4
1090 save big, machep, tiny
1093 data big/0.d+0/, machep/-1.d+0/, tiny/1.d+0/, zero/0.d+0/
1095 ! data vn(1,1),vn(2,1)/4hepsl,4hon../
1096 ! data vn(1,2),vn(2,2)/4hphmn,4hfc../
1097 ! data vn(1,3),vn(2,3)/4hphmx,4hfc../
1098 ! data vn(1,4),vn(2,4)/4hdecf,4hac../
1099 ! data vn(1,5),vn(2,5)/4hincf,4hac../
1100 ! data vn(1,6),vn(2,6)/4hrdfc,4hmn../
1101 ! data vn(1,7),vn(2,7)/4hrdfc,4hmx../
1102 ! data vn(1,8),vn(2,8)/4htune,4hr1../
1103 ! data vn(1,9),vn(2,9)/4htune,4hr2../
1104 ! data vn(1,10),vn(2,10)/4htune,4hr3../
1105 ! data vn(1,11),vn(2,11)/4htune,4hr4../
1106 ! data vn(1,12),vn(2,12)/4htune,4hr5../
1107 ! data vn(1,13),vn(2,13)/4hafct,4hol../
1108 ! data vn(1,14),vn(2,14)/4hrfct,4hol../
1109 ! data vn(1,15),vn(2,15)/4hxcto,4hl.../
1110 ! data vn(1,16),vn(2,16)/4hxfto,4hl.../
1111 ! data vn(1,17),vn(2,17)/4hlmax,4h0.../
1112 ! data vn(1,18),vn(2,18)/4hlmax,4hs.../
1113 ! data vn(1,19),vn(2,19)/4hscto,4hl.../
1114 ! data vn(1,20),vn(2,20)/4hdini,4ht.../
1115 ! data vn(1,21),vn(2,21)/4hdtin,4hit../
1116 ! data vn(1,22),vn(2,22)/4hd0in,4hit../
1117 ! data vn(1,23),vn(2,23)/4hdfac,4h..../
1118 ! data vn(1,24),vn(2,24)/4hdltf,4hdc../
1119 ! data vn(1,25),vn(2,25)/4hdltf,4hdj../
1120 ! data vn(1,26),vn(2,26)/4hdelt,4ha0../
1121 ! data vn(1,27),vn(2,27)/4hfuzz,4h..../
1122 ! data vn(1,28),vn(2,28)/4hrlim,4hit../
1123 ! data vn(1,29),vn(2,29)/4hcosm,4hin../
1124 ! data vn(1,30),vn(2,30)/4hhube,4hrc../
1125 ! data vn(1,31),vn(2,31)/4hrspt,4hol../
1126 ! data vn(1,32),vn(2,32)/4hsigm,4hin../
1127 ! data vn(1,33),vn(2,33)/4heta0,4h..../
1128 ! data vn(1,34),vn(2,34)/4hbias,4h..../
1130 data vn(1,1),vn(2,1)/'epsl','on..'/
1131 data vn(1,2),vn(2,2)/'phmn','fc..'/
1132 data vn(1,3),vn(2,3)/'phmx','fc..'/
1133 data vn(1,4),vn(2,4)/'decf','ac..'/
1134 data vn(1,5),vn(2,5)/'incf','ac..'/
1135 data vn(1,6),vn(2,6)/'rdfc','mn..'/
1136 data vn(1,7),vn(2,7)/'rdfc','mx..'/
1137 data vn(1,8),vn(2,8)/'tune','r1..'/
1138 data vn(1,9),vn(2,9)/'tune','r2..'/
1139 data vn(1,10),vn(2,10)/'tune','r3..'/
1140 data vn(1,11),vn(2,11)/'tune','r4..'/
1141 data vn(1,12),vn(2,12)/'tune','r5..'/
1142 data vn(1,13),vn(2,13)/'afct','ol..'/
1143 data vn(1,14),vn(2,14)/'rfct','ol..'/
1144 data vn(1,15),vn(2,15)/'xcto','l...'/
1145 data vn(1,16),vn(2,16)/'xfto','l...'/
1146 data vn(1,17),vn(2,17)/'lmax','0...'/
1147 data vn(1,18),vn(2,18)/'lmax','s...'/
1148 data vn(1,19),vn(2,19)/'scto','l...'/
1149 data vn(1,20),vn(2,20)/'dini','t...'/
1150 data vn(1,21),vn(2,21)/'dtin','it..'/
1151 data vn(1,22),vn(2,22)/'d0in','it..'/
1152 data vn(1,23),vn(2,23)/'dfac','....'/
1153 data vn(1,24),vn(2,24)/'dltf','dc..'/
1154 data vn(1,25),vn(2,25)/'dltf','dj..'/
1155 data vn(1,26),vn(2,26)/'delt','a0..'/
1156 data vn(1,27),vn(2,27)/'fuzz','....'/
1157 data vn(1,28),vn(2,28)/'rlim','it..'/
1158 data vn(1,29),vn(2,29)/'cosm','in..'/
1159 data vn(1,30),vn(2,30)/'hube','rc..'/
1160 data vn(1,31),vn(2,31)/'rspt','ol..'/
1161 data vn(1,32),vn(2,32)/'sigm','in..'/
1162 data vn(1,33),vn(2,33)/'eta0','....'/
1163 data vn(1,34),vn(2,34)/'bias','....'/
1166 data vm(1)/1.0d-3/, vm(2)/-0.99d+0/, vm(3)/1.0d-3/, vm(4)/1.0d-2/,&
1167 vm(5)/1.2d+0/, vm(6)/1.d-2/, vm(7)/1.2d+0/, vm(8)/0.d+0/,&
1168 vm(9)/0.d+0/, vm(10)/1.d-3/, vm(11)/-1.d+0/, vm(13)/0.d+0/,&
1169 vm(15)/0.d+0/, vm(16)/0.d+0/, vm(19)/0.d+0/, vm(20)/-10.d+0/,&
1170 vm(21)/0.d+0/, vm(22)/0.d+0/, vm(23)/0.d+0/, vm(27)/1.01d+0/,&
1171 vm(28)/1.d+10/, vm(30)/0.d+0/, vm(31)/0.d+0/, vm(32)/0.d+0/,&
1173 data vx(1)/0.9d+0/, vx(2)/-1.d-3/, vx(3)/1.d+1/, vx(4)/0.8d+0/,&
1174 vx(5)/1.d+2/, vx(6)/0.8d+0/, vx(7)/1.d+2/, vx(8)/0.5d+0/,&
1175 vx(9)/0.5d+0/, vx(10)/1.d+0/, vx(11)/1.d+0/, vx(14)/0.1d+0/,&
1176 vx(15)/1.d+0/, vx(16)/1.d+0/, vx(19)/1.d+0/, vx(23)/1.d+0/,&
1177 vx(24)/1.d+0/, vx(25)/1.d+0/, vx(26)/1.d+0/, vx(27)/1.d+10/,&
1178 vx(29)/1.d+0/, vx(31)/1.d+0/, vx(32)/1.d+0/, vx(33)/1.d+0/,&
1182 ! data varnm(1)/1hp/, varnm(2)/1hn/, sh(1)/1hs/, sh(2)/1hh/
1183 ! data cngd(1),cngd(2),cngd(3)/4h---c,4hhang,4hed v/,
1184 ! 1 dflt(1),dflt(2),dflt(3)/4hnond,4hefau,4hlt v/
1186 data varnm(1)/'p'/, varnm(2)/'n'/, sh(1)/'s'/, sh(2)/'h'/
1187 data cngd(1),cngd(2),cngd(3)/'---c','hang','ed v'/,&
1188 dflt(1),dflt(2),dflt(3)/'nond','efau','lt v'/
1190 data ijmp/33/, jlim(1)/0/, jlim(2)/24/, ndflt(1)/32/, ndflt(2)/25/
1191 data miniv(1)/80/, miniv(2)/59/
1193 !............................... body ................................
1196 if (prunit .le. liv) pu = iv(prunit)
1197 if (alg .lt. 1 .or. alg .gt. 2) go to 340
1198 if (iv(1) .eq. 0) call deflt(alg, iv, liv, lv, v)
1200 if (iv1 .ne. 13 .and. iv1 .ne. 12) go to 10
1202 if (perm .le. liv) miv1 = max0(miv1, iv(perm) - 1)
1203 if (ivneed .le. liv) miv2 = miv1 + max0(iv(ivneed), 0)
1204 if (lastiv .le. liv) iv(lastiv) = miv2
1205 if (liv .lt. miv1) go to 300
1207 iv(lastv) = max0(iv(vneed), 0) + iv(lmat) - 1
1209 if (liv .lt. miv2) go to 300
1210 if (lv .lt. iv(lastv)) go to 320
1211 10 if (alg .eq. iv(algsav)) go to 30
1212 if (pu .ne. 0) write(pu,20) alg, iv(algsav)
1213 20 format(/39h the first parameter to deflt should be,i3,&
1217 30 if (iv1 .lt. 12 .or. iv1 .gt. 14) go to 60
1218 if (n .ge. 1) go to 50
1220 if (pu .eq. 0) go to 999
1221 write(pu,40) varnm(alg), n
1222 40 format(/8h /// bad,a1,2h =,i5)
1224 50 if (iv1 .ne. 14) iv(nextiv) = iv(perm)
1225 if (iv1 .ne. 14) iv(nextv) = iv(lmat)
1226 if (iv1 .eq. 13) go to 999
1227 k = iv(parsav) - epslon
1228 call vdflt(alg, lv-k, v(k+1))
1229 iv(dtype0) = 2 - alg
1235 60 if (n .eq. iv(oldn)) go to 80
1237 if (pu .eq. 0) go to 999
1238 write(pu,70) varnm(alg), iv(oldn), n
1239 70 format(/5h /// ,1a1,14h changed from ,i5,4h to ,i5)
1242 80 if (iv1 .le. 11 .and. iv1 .ge. 1) go to 100
1244 if (pu .ne. 0) write(pu,90) iv1
1245 90 format(/13h /// iv(1) =,i5,28h should be between 0 and 14.)
1248 100 which(1) = cngd(1)
1252 110 if (iv1 .eq. 14) iv1 = 12
1253 if (big .gt. tiny) go to 120
1280 do 150 l = 1, ndfalt
1282 if (vk .ge. vm(i) .and. vk .le. vx(i)) go to 140
1284 if (pu .ne. 0) write(pu,130) vn(1,i), vn(2,i), k, vk,&
1286 130 format(/6h /// ,2a4,5h.. v(,i2,3h) =,d11.3,7h should,&
1287 11h be between,d11.3,4h and,d11.3)
1290 if (i .eq. j) i = ijmp
1293 if (iv(nvdflt) .eq. ndfalt) go to 170
1295 if (pu .eq. 0) go to 999
1296 write(pu,160) iv(nvdflt), ndfalt
1297 160 format(/13h iv(nvdflt) =,i5,13h rather than ,i5)
1299 170 if ((iv(dtype) .gt. 0 .or. v(dinit) .gt. zero) .and. iv1 .eq. 12) &
1302 if (d(i) .gt. zero) go to 190
1304 if (pu .ne. 0) write(pu,180) i, d(i)
1305 180 format(/8h /// d(,i3,3h) =,d11.3,19h should be positive)
1307 200 if (m .eq. 0) go to 210
1311 210 if (pu .eq. 0 .or. iv(parprt) .eq. 0) go to 999
1312 if (iv1 .ne. 12 .or. iv(inits) .eq. alg-1) go to 230
1314 write(pu,220) sh(alg), iv(inits)
1315 220 format(/22h nondefault values..../5h init,a1,14h..... iv(25) =,&
1317 230 if (iv(dtype) .eq. iv(dtype0)) go to 250
1318 if (m .eq. 0) write(pu,260) which
1320 write(pu,240) iv(dtype)
1321 240 format(20h dtype..... iv(16) =,i3)
1327 do 290 ii = 1, ndfalt
1328 if (v(k) .eq. v(l)) go to 280
1329 if (m .eq. 0) write(pu,260) which
1330 260 format(/1h ,3a4,9halues..../)
1332 write(pu,270) vn(1,i), vn(2,i), k, v(k)
1333 270 format(1x,2a4,5h.. v(,i2,3h) =,d15.7)
1337 if (i .eq. j) i = ijmp
1340 iv(dtype0) = iv(dtype)
1342 call vcopy(iv(nvdflt), v(parsv1), v(epslon))
1346 if (pu .eq. 0) go to 999
1347 write(pu,310) liv, miv2
1348 310 format(/10h /// liv =,i5,17h must be at least,i5)
1349 if (liv .lt. miv1) go to 999
1350 if (lv .lt. iv(lastv)) go to 320
1354 if (pu .eq. 0) go to 999
1355 write(pu,330) lv, iv(lastv)
1356 330 format(/9h /// lv =,i5,17h must be at least,i5)
1360 if (pu .eq. 0) go to 999
1362 350 format(/10h /// alg =,i5,15h must be 1 or 2)
1365 ! *** last card of parck follows ***
1366 end subroutine parck
1367 !-----------------------------------------------------------------------------
1368 real(kind=8) function reldst(p, d, x, x0)
1370 ! *** compute and return relative difference between x and x0 ***
1371 ! *** nl2sol version 2.2 ***
1374 real(kind=8) :: d(p), x(p), x0(p)
1376 !el real(kind=8) :: dabs
1379 real(kind=8) :: emax, t, xmax !el, zero
1383 real(kind=8),parameter :: zero=0.d+0
1389 t = dabs(d(i) * (x(i) - x0(i)))
1390 if (emax .lt. t) emax = t
1391 t = d(i) * (dabs(x(i)) + dabs(x0(i)))
1392 if (xmax .lt. t) xmax = t
1395 if (xmax .gt. zero) reldst = emax / xmax
1397 ! *** last card of reldst follows ***
1399 !-----------------------------------------------------------------------------
1400 subroutine vaxpy(p, w, a, x, y)
1402 ! *** set w = a*x + y -- w, x, y = p-vectors, a = scalar ***
1405 real(kind=8) :: a, w(p), x(p), y(p)
1410 10 w(i) = a*x(i) + y(i)
1412 end subroutine vaxpy
1413 !-----------------------------------------------------------------------------
1414 subroutine vcopy(p, y, x)
1416 ! *** set y = x, where x and y are p-vectors ***
1419 real(kind=8) :: x(p), y(p)
1426 end subroutine vcopy
1427 !-----------------------------------------------------------------------------
1428 subroutine vdflt(alg, lv, v)
1430 ! *** supply ***sol (version 2.3) default values to v ***
1432 ! *** alg = 1 means regression constants.
1433 ! *** alg = 2 means general unconstrained optimization constants.
1437 real(kind=8) :: v(lv)
1439 !el real(kind=8) :: dmax1
1442 !el real(kind=8) :: rmdcon
1443 ! rmdcon... returns machine-dependent constants
1445 real(kind=8) :: machep, mepcrt, sqteps !el one, three
1447 ! *** subscripts for v ***
1449 !el integer afctol, bias, cosmin, decfac, delta0, dfac, dinit, dltfdc,
1450 !el 1 dltfdj, dtinit, d0init, epslon, eta0, fuzz, huberc,
1451 !el 2 incfac, lmax0, lmaxs, phmnfc, phmxfc, rdfcmn, rdfcmx,
1452 !el 3 rfctol, rlimit, rsptol, sctol, sigmin, tuner1, tuner2,
1453 !el 4 tuner3, tuner4, tuner5, xctol, xftol
1456 ! data one/1.d+0/, three/3.d+0/
1458 real(kind=8),parameter :: one=1.d+0, three=3.d+0
1461 ! *** v subscript values ***
1464 ! data afctol/31/, bias/43/, cosmin/47/, decfac/22/, delta0/44/,
1465 ! 1 dfac/41/, dinit/38/, dltfdc/42/, dltfdj/43/, dtinit/39/,
1466 ! 2 d0init/40/, epslon/19/, eta0/42/, fuzz/45/, huberc/48/,
1467 ! 3 incfac/23/, lmax0/35/, lmaxs/36/, phmnfc/20/, phmxfc/21/,
1468 ! 4 rdfcmn/24/, rdfcmx/25/, rfctol/32/, rlimit/46/, rsptol/49/,
1469 ! 5 sctol/37/, sigmin/50/, tuner1/26/, tuner2/27/, tuner3/28/,
1470 ! 6 tuner4/29/, tuner5/30/, xctol/33/, xftol/34/
1472 integer,parameter :: afctol=31, bias=43, cosmin=47, decfac=22, delta0=44,&
1473 dfac=41, dinit=38, dltfdc=42, dltfdj=43, dtinit=39,&
1474 d0init=40, epslon=19, eta0=42, fuzz=45, huberc=48,&
1475 incfac=23, lmax0=35, lmaxs=36, phmnfc=20, phmxfc=21,&
1476 rdfcmn=24, rdfcmx=25, rfctol=32, rlimit=46, rsptol=49,&
1477 sctol=37, sigmin=50, tuner1=26, tuner2=27, tuner3=28,&
1478 tuner4=29, tuner5=30, xctol=33, xftol=34
1481 !------------------------------- body --------------------------------
1485 if (machep .gt. 1.d-10) v(afctol) = machep**2
1491 mepcrt = machep ** (one/three)
1501 v(rfctol) = dmax1(1.d-10, mepcrt**2)
1502 v(sctol) = v(rfctol)
1509 v(xftol) = 1.d+2 * machep
1511 if (alg .ge. 2) go to 10
1513 ! *** regression values
1515 v(cosmin) = dmax1(1.d-6, 1.d+2 * machep)
1521 v(rlimit) = rmdcon(5)
1526 ! *** general optimization values
1530 v(eta0) = 1.0d+3 * machep
1533 ! *** last card of vdflt follows ***
1534 end subroutine vdflt
1535 !-----------------------------------------------------------------------------
1536 subroutine vscopy(p, y, s)
1538 ! *** set p-vector y to scalar s ***
1541 real(kind=8) :: s, y(p)
1548 end subroutine vscopy
1549 !-----------------------------------------------------------------------------
1550 real(kind=8) function v2norm(p, x)
1552 ! *** return the 2-norm of the p-vector x, taking ***
1553 ! *** care to avoid the most likely underflows. ***
1556 real(kind=8) :: x(p)
1559 real(kind=8) :: r, scale, sqteta, t, xi !el, one, zero
1561 !el real(kind=8) :: dabs, dsqrt
1564 !el real(kind=8) :: rmdcon
1567 ! data one/1.d+0/, zero/0.d+0/
1569 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
1574 if (p .gt. 0) go to 10
1578 if (x(i) .ne. zero) go to 30
1583 30 scale = dabs(x(i))
1584 if (i .lt. p) go to 40
1588 if (sqteta .eq. zero) sqteta = rmdcon(2)
1590 ! *** sqteta is (slightly larger than) the square root of the
1591 ! *** smallest positive floating point number on the machine.
1592 ! *** the tests involving sqteta are done to prevent underflows.
1597 if (xi .gt. scale) go to 50
1599 if (r .gt. sqteta) t = t + r*r
1602 if (r .le. sqteta) r = zero
1607 v2norm = scale * dsqrt(t)
1609 ! *** last card of v2norm follows ***
1611 !-----------------------------------------------------------------------------
1612 subroutine humsl(n,d,x,calcf,calcgh,iv,liv,lv,v,uiparm,urparm,ufparm)
1614 ! *** minimize general unconstrained objective function using ***
1615 ! *** (analytic) gradient and hessian provided by the caller. ***
1620 integer :: uiparm(1)
1621 real(kind=8) :: d(n), x(n), v(lv), urparm(1)
1622 real(kind=8),external :: ufparm
1623 ! dimension v(78 + n*(n+12)), uiparm(*), urparm(*)
1624 external :: calcf, calcgh
1626 !------------------------------ discussion ---------------------------
1628 ! this routine is like sumsl, except that the subroutine para-
1629 ! meter calcg of sumsl (which computes the gradient of the objec-
1630 ! tive function) is replaced by the subroutine parameter calcgh,
1631 ! which computes both the gradient and (lower triangle of the)
1632 ! hessian of the objective function. the calling sequence is...
1633 ! call calcgh(n, x, nf, g, h, uiparm, urparm, ufparm)
1634 ! parameters n, x, nf, g, uiparm, urparm, and ufparm are the same
1635 ! as for sumsl, while h is an array of length n*(n+1)/2 in which
1636 ! calcgh must store the lower triangle of the hessian at x. start-
1637 ! ing at h(1), calcgh must store the hessian entries in the order
1638 ! (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ...
1639 ! the value printed (by itsum) in the column labelled stppar
1640 ! is the levenberg-marquardt used in computing the current step.
1641 ! zero means a full newton step. if the special case described in
1642 ! ref. 1 is detected, then stppar is negated. the value printed
1643 ! in the column labelled npreldf is zero if the current hessian
1644 ! is not positive definite.
1645 ! it sometimes proves worthwhile to let d be determined from the
1646 ! diagonal of the hessian matrix by setting iv(dtype) = 1 and
1647 ! v(dinit) = 0. the following iv and v components are relevant...
1649 ! iv(dtol)..... iv(59) gives the starting subscript in v of the dtol
1650 ! array used when d is updated. (iv(dtol) can be
1651 ! initialized by calling humsl with iv(1) = 13.)
1652 ! iv(dtype).... iv(16) tells how the scale vector d should be chosen.
1653 ! iv(dtype) .le. 0 means that d should not be updated, and
1654 ! iv(dtype) .ge. 1 means that d should be updated as
1655 ! described below with v(dfac). default = 0.
1656 ! v(dfac)..... v(41) and the dtol and d0 arrays (see v(dtinit) and
1657 ! v(d0init)) are used in updating the scale vector d when
1658 ! iv(dtype) .gt. 0. (d is initialized according to
1659 ! v(dinit), described in sumsl.) let
1660 ! d1(i) = max(sqrt(abs(h(i,i))), v(dfac)*d(i)),
1661 ! where h(i,i) is the i-th diagonal element of the current
1662 ! hessian. if iv(dtype) = 1, then d(i) is set to d1(i)
1663 ! unless d1(i) .lt. dtol(i), in which case d(i) is set to
1664 ! max(d0(i), dtol(i)).
1665 ! if iv(dtype) .ge. 2, then d is updated during the first
1666 ! iteration as for iv(dtype) = 1 (after any initialization
1667 ! due to v(dinit)) and is left unchanged thereafter.
1669 ! v(dtinit)... v(39), if positive, is the value to which all components
1670 ! of the dtol array (see v(dfac)) are initialized. if
1671 ! v(dtinit) = 0, then it is assumed that the caller has
1672 ! stored dtol in v starting at v(iv(dtol)).
1674 ! v(d0init)... v(40), if positive, is the value to which all components
1675 ! of the d0 vector (see v(dfac)) are initialized. if
1676 ! v(dfac) = 0, then it is assumed that the caller has
1677 ! stored d0 in v starting at v(iv(dtol)+n). default = 1.0.
1681 ! 1. gay, d.m. (1981), computing optimal locally constrained steps,
1682 ! siam j. sci. statist. comput. 2, pp. 186-197.
1686 ! coded by david m. gay (winter 1980). revised sept. 1982.
1687 ! this subroutine was written in connection with research supported
1688 ! in part by the national science foundation under grants
1689 ! mcs-7600324 and mcs-7906671.
1691 !---------------------------- declarations ---------------------------
1693 !el external deflt, humit
1695 ! deflt... provides default input values for iv and v.
1696 ! humit... reverse-communication routine that does humsl algorithm.
1698 integer :: g1, h1, iv1, lh, nf
1701 ! *** subscripts for iv ***
1703 !el integer g, h, nextv, nfcall, nfgcal, toobig, vneed
1706 ! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, h/56/, toobig/2/,
1709 integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28, h=56,&
1713 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
1715 lh = n * (n + 1) / 2
1716 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1717 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
1718 iv(vneed) = iv(vneed) + n*(n+3)/2
1720 if (iv1 .eq. 14) go to 10
1721 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
1724 if (iv1 .eq. 12) iv(1) = 13
1730 20 call humit(d, f, v(g1), v(h1), iv, lh, liv, lv, n, v, x)
1731 if (iv(1) - 2) 30, 40, 50
1734 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
1735 if (nf .le. 0) iv(toobig) = 1
1738 40 call calcgh(n, x, iv(nfgcal), v(g1), v(h1), uiparm, urparm,&
1742 50 if (iv(1) .ne. 14) go to 999
1744 ! *** storage allocation
1748 iv(nextv) = iv(h) + n*(n+1)/2
1749 if (iv1 .ne. 13) go to 10
1752 ! *** last card of humsl follows ***
1753 end subroutine humsl
1754 !-----------------------------------------------------------------------------
1755 subroutine humit(d, fx, g, h, iv, lh, liv, lv, n, v, x)
1757 ! *** carry out humsl (unconstrained minimization) iterations, using
1758 ! *** hessian matrix provided by the caller.
1761 use control, only:stopx
1763 ! *** parameter declarations ***
1765 integer :: lh, liv, n
1768 real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n)
1770 !-------------------------- parameter usage --------------------------
1772 ! d.... scale vector.
1773 ! fx... function value.
1774 ! g.... gradient vector.
1775 ! h.... lower triangle of the hessian, stored rowwise.
1776 ! iv... integer value array.
1777 ! lh... length of h = p*(p+1)/2.
1778 ! liv.. length of iv (at least 60).
1779 ! lv... length of v (at least 78 + n*(n+21)/2).
1780 ! n.... number of variables (components in x and g).
1781 ! v.... floating-point value array.
1782 ! x.... parameter vector.
1784 ! *** discussion ***
1786 ! parameters iv, n, v, and x are the same as the corresponding
1787 ! ones to humsl (which see), except that v can be shorter (since
1788 ! the part of v that humsl uses for storing g and h is not needed).
1789 ! moreover, compared with humsl, iv(1) may have the two additional
1790 ! output values 1 and 2, which are explained below, as is the use
1791 ! of iv(toobig) and iv(nfgcal). the value iv(g), which is an
1792 ! output value from humsl, is not referenced by humit or the
1793 ! subroutines it calls.
1795 ! iv(1) = 1 means the caller should set fx to f(x), the function value
1796 ! at x, and call humit again, having changed none of the
1797 ! other parameters. an exception occurs if f(x) cannot be
1798 ! computed (e.g. if overflow would occur), which may happen
1799 ! because of an oversized step. in this case the caller
1800 ! should set iv(toobig) = iv(2) to 1, which will cause
1801 ! humit to ignore fx and try a smaller step. the para-
1802 ! meter nf that humsl passes to calcf (for possible use by
1803 ! calcgh) is a copy of iv(nfcall) = iv(6).
1804 ! iv(1) = 2 means the caller should set g to g(x), the gradient of f at
1805 ! x, and h to the lower triangle of h(x), the hessian of f
1806 ! at x, and call humit again, having changed none of the
1807 ! other parameters except perhaps the scale vector d.
1808 ! the parameter nf that humsl passes to calcg is
1809 ! iv(nfgcal) = iv(7). if g(x) and h(x) cannot be evaluated,
1810 ! then the caller may set iv(nfgcal) to 0, in which case
1811 ! humit will return with iv(1) = 65.
1812 ! note -- humit overwrites h with the lower triangle
1813 ! of diag(d)**-1 * h(x) * diag(d)**-1.
1817 ! coded by david m. gay (winter 1980). revised sept. 1982.
1818 ! this subroutine was written in connection with research supported
1819 ! in part by the national science foundation under grants
1820 ! mcs-7600324 and mcs-7906671.
1822 ! (see sumsl and humsl for references.)
1824 !+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
1826 ! *** local variables ***
1828 integer :: dg1, dummy, i, j, k, l, lstgst, nn1o2, step1,&
1834 !el real(kind=8) :: one, onep2, zero
1836 ! *** no intrinsic functions ***
1838 ! *** external functions and subroutines ***
1840 !el external assst, deflt, dotprd, dupdu, gqtst, itsum, parck,
1841 !el 1 reldst, slvmul, stopx, vaxpy, vcopy, vscopy, v2norm
1843 !el real(kind=8) :: dotprd, reldst, v2norm
1845 ! assst.... assesses candidate step.
1846 ! deflt.... provides default iv and v input values.
1847 ! dotprd... returns inner product of two vectors.
1848 ! dupdu.... updates scale vector d.
1849 ! gqtst.... computes optimally locally constrained step.
1850 ! itsum.... prints iteration summary and info on initial and final x.
1851 ! parck.... checks validity of input iv and v values.
1852 ! reldst... computes v(reldx) = relative step size.
1853 ! slvmul... multiplies symmetric matrix times vector, given the lower
1854 ! triangle of the matrix.
1855 ! stopx.... returns .true. if the break key has been pressed.
1856 ! vaxpy.... computes scalar times one vector plus another.
1857 ! vcopy.... copies one vector to another.
1858 ! vscopy... sets all elements of a vector to a scalar.
1859 ! v2norm... returns the 2-norm of a vector.
1861 ! *** subscripts for iv and v ***
1863 !el integer cnvcod, dg, dgnorm, dinit, dstnrm, dtinit, dtol,
1864 !el 1 dtype, d0init, f, f0, fdif, gtstep, incfac, irc, kagqt,
1865 !el 2 lmat, lmax0, lmaxs, mode, model, mxfcal, mxiter, nextv,
1866 !el 3 nfcall, nfgcal, ngcall, niter, preduc, radfac, radinc,
1867 !el 4 radius, rad0, reldx, restor, step, stglim, stlstg, stppar,
1868 !el 5 toobig, tuner4, tuner5, vneed, w, xirc, x0
1870 ! *** iv subscript values ***
1873 ! data cnvcod/55/, dg/37/, dtol/59/, dtype/16/, irc/29/, kagqt/33/,
1874 ! 1 lmat/42/, mode/35/, model/5/, mxfcal/17/, mxiter/18/,
1875 ! 2 nextv/47/, nfcall/6/, nfgcal/7/, ngcall/30/, niter/31/,
1876 ! 3 radinc/8/, restor/9/, step/40/, stglim/11/, stlstg/41/,
1877 ! 4 toobig/2/, vneed/4/, w/34/, xirc/13/, x0/43/
1879 integer,parameter :: cnvcod=55, dg=37, dtol=59, dtype=16, irc=29, kagqt=33,&
1880 lmat=42, mode=35, model=5, mxfcal=17, mxiter=18,&
1881 nextv=47, nfcall=6, nfgcal=7, ngcall=30, niter=31,&
1882 radinc=8, restor=9, step=40, stglim=11, stlstg=41,&
1883 toobig=2, vneed=4, w=34, xirc=13, x0=43
1886 ! *** v subscript values ***
1889 ! data dgnorm/1/, dinit/38/, dstnrm/2/, dtinit/39/, d0init/40/,
1890 ! 1 f/10/, f0/13/, fdif/11/, gtstep/4/, incfac/23/, lmax0/35/,
1891 ! 2 lmaxs/36/, preduc/7/, radfac/16/, radius/8/, rad0/9/,
1892 ! 3 reldx/17/, stppar/5/, tuner4/29/, tuner5/30/
1894 integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dtinit=39, d0init=40,&
1895 f=10, f0=13, fdif=11, gtstep=4, incfac=23, lmax0=35,&
1896 lmaxs=36, preduc=7, radfac=16, radius=8, rad0=9,&
1897 reldx=17, stppar=5, tuner4=29, tuner5=30
1901 ! data one/1.d+0/, onep2/1.2d+0/, zero/0.d+0/
1903 real(kind=8),parameter :: one=1.d+0, onep2=1.2d+0, zero=0.d+0
1906 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
1909 if (i .eq. 1) go to 30
1910 if (i .eq. 2) go to 40
1912 ! *** check validity of iv and v input values ***
1914 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1915 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
1916 iv(vneed) = iv(vneed) + n*(n+21)/2 + 7
1917 call parck(2, d, iv, liv, lv, n, v)
1919 if (i .gt. 12) go to 999
1920 nn1o2 = n * (n + 1) / 2
1921 if (lh .ge. nn1o2) go to (210,210,210,210,210,210,160,120,160,&
1926 ! *** storage allocation ***
1928 10 iv(dtol) = iv(lmat) + nn1o2
1929 iv(x0) = iv(dtol) + 2*n
1930 iv(step) = iv(x0) + n
1931 iv(stlstg) = iv(step) + n
1932 iv(dg) = iv(stlstg) + n
1934 iv(nextv) = iv(w) + 4*n + 7
1935 if (iv(1) .ne. 13) go to 20
1939 ! *** initialization ***
1953 if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
1955 if (v(dtinit) .gt. zero) call vscopy(n, v(k), v(dtinit))
1957 if (v(d0init) .gt. zero) call vscopy(n, v(k), v(d0init))
1962 if (iv(mode) .ge. 0) go to 210
1964 if (iv(toobig) .eq. 0) go to 999
1968 ! *** make sure gradient could be computed ***
1970 40 if (iv(nfgcal) .ne. 0) go to 50
1974 ! *** update the scale vector d ***
1977 if (iv(dtype) .le. 0) go to 70
1985 call dupdu(d, v(dg1), iv, liv, lv, n, v)
1987 ! *** compute scaled gradient and its norm ***
1995 v(dgnorm) = v2norm(n, v(dg1))
1997 ! *** compute scaled hessian ***
2003 h(k) = t * h(k) / d(j)
2008 if (iv(cnvcod) .ne. 0) go to 340
2009 if (iv(mode) .eq. 0) go to 300
2011 ! *** allow first step to have scaled 2-norm at most v(lmax0) ***
2013 v(radius) = v(lmax0)
2018 !----------------------------- main loop -----------------------------
2021 ! *** print iteration summary, check iteration limit ***
2023 110 call itsum(d, g, iv, liv, lv, n, v, x)
2025 if (k .lt. iv(mxiter)) go to 130
2029 130 iv(niter) = k + 1
2031 ! *** initialize for start of next iteration ***
2039 ! *** copy x to x0 ***
2041 call vcopy(n, v(x01), x)
2043 ! *** update radius ***
2045 if (k .eq. 0) go to 150
2052 v(radius) = v(radfac) * v2norm(n, v(step1))
2054 ! *** check stopx and function evaluation limit ***
2058 150 if (.not. stopx(dummy)) go to 170
2062 ! *** come here when restarting after func. eval. limit or stopx.
2064 160 if (v(f) .ge. v(f0)) go to 170
2069 170 if (iv(nfcall) .lt. iv(mxfcal)) go to 190
2071 180 if (v(f) .ge. v(f0)) go to 350
2073 ! *** in case of stopx or function evaluation limit with
2074 ! *** improved v(f), evaluate the gradient at x.
2079 !. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
2081 190 step1 = iv(step)
2085 call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1))
2086 if (iv(irc) .eq. 6) go to 210
2088 ! *** check whether evaluating f(x0 + step) looks worthwhile ***
2090 if (v(dstnrm) .le. zero) go to 210
2091 if (iv(irc) .ne. 5) go to 200
2092 if (v(radfac) .le. one) go to 200
2093 if (v(preduc) .le. onep2 * v(fdif)) go to 210
2095 ! *** compute f(x0 + step) ***
2099 call vaxpy(n, x, one, v(step1), v(x01))
2100 iv(nfcall) = iv(nfcall) + 1
2105 !. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
2108 v(reldx) = reldst(n, d, x, v(x01))
2109 call assst(iv, liv, lv, v)
2112 if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
2113 if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
2114 if (iv(restor) .ne. 3) go to 220
2115 call vcopy(n, v(step1), v(lstgst))
2116 call vaxpy(n, x, one, v(step1), v(x01))
2117 v(reldx) = reldst(n, d, x, v(x01))
2120 go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k
2122 ! *** recompute step with new radius ***
2124 230 v(radius) = v(radfac) * v(dstnrm)
2127 ! *** compute step of length v(lmaxs) for singular convergence test.
2129 240 v(radius) = v(lmaxs)
2132 ! *** convergence or false convergence ***
2134 250 iv(cnvcod) = k - 4
2135 if (v(f) .ge. v(f0)) go to 340
2136 if (iv(xirc) .eq. 14) go to 340
2139 !. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
2141 260 if (iv(irc) .ne. 3) go to 290
2144 ! *** prepare for gradient tests ***
2145 ! *** set temp1 = hessian * step + g(x0)
2146 ! *** = diag(d) * (h * step + g(x0))
2148 ! use x0 vector as temporary.
2151 v(k) = d(i) * v(step1)
2155 call slvmul(n, v(temp1), h, v(x01))
2157 v(temp1) = d(i) * v(temp1) + g(i)
2161 ! *** compute gradient and hessian ***
2163 290 iv(ngcall) = iv(ngcall) + 1
2168 if (iv(irc) .ne. 3) go to 110
2170 ! *** set v(radfac) by gradient tests ***
2175 ! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
2179 v(k) = (v(k) - g(i)) / d(i)
2183 ! *** do gradient tests ***
2185 if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320
2186 if (dotprd(n, g, v(step1)) &
2187 .ge. v(gtstep) * v(tuner5)) go to 110
2188 320 v(radfac) = v(incfac)
2191 !. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
2193 ! *** bad parameters to assess ***
2198 ! *** print summary of final iteration and other requested items ***
2200 340 iv(1) = iv(cnvcod)
2202 350 call itsum(d, g, iv, liv, lv, n, v, x)
2206 ! *** last card of humit follows ***
2207 end subroutine humit
2208 !-----------------------------------------------------------------------------
2209 subroutine dupdu(d, hdiag, iv, liv, lv, n, v)
2211 ! *** update scale vector d for humsl ***
2213 ! *** parameter declarations ***
2218 real(kind=8) :: d(n), hdiag(n), v(lv)
2220 ! *** local variables ***
2222 integer :: dtoli, d0i, i
2223 real(kind=8) :: t, vdfac
2225 ! *** intrinsic functions ***
2227 !el real(kind=8) :: dabs, dmax1, dsqrt
2229 ! *** subscripts for iv and v ***
2231 !el integer :: dfac, dtol, dtype, niter
2233 ! data dfac/41/, dtol/59/, dtype/16/, niter/31/
2235 integer,parameter :: dfac=41, dtol=59, dtype=16, niter=31
2238 !------------------------------- body --------------------------------
2241 if (i .eq. 1) go to 10
2242 if (iv(niter) .gt. 0) go to 999
2248 t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i))
2249 if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i))
2256 ! *** last card of dupdu follows ***
2257 end subroutine dupdu
2258 !-----------------------------------------------------------------------------
2259 subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w)
2261 ! *** compute goldfeld-quandt-trotter step by more-hebden technique ***
2262 ! *** (nl2sol version 2.2), modified a la more and sorensen ***
2264 ! *** parameter declarations ***
2267 !al real(kind=8) :: d(p), dig(p), dihdi(1), l(1), v(21), step(p),
2269 real(kind=8) :: d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),&
2270 v(21), step(p),w(4*p+7)
2271 ! dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7)
2273 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2277 ! given the (compactly stored) lower triangle of a scaled
2278 ! hessian (approximation) and a nonzero scaled gradient vector,
2279 ! this subroutine computes a goldfeld-quandt-trotter step of
2280 ! approximate length v(radius) by the more-hebden technique. in
2281 ! other words, step is computed to (approximately) minimize
2282 ! psi(step) = (g**t)*step + 0.5*(step**t)*h*step such that the
2283 ! 2-norm of d*step is at most (approximately) v(radius), where
2284 ! g is the gradient, h is the hessian, and d is a diagonal
2285 ! scale matrix whose diagonal is stored in the parameter d.
2286 ! (gqtst assumes dig = d**-1 * g and dihdi = d**-1 * h * d**-1.)
2288 ! *** parameter description ***
2290 ! d (in) = the scale vector, i.e. the diagonal of the scale
2291 ! matrix d mentioned above under purpose.
2292 ! dig (in) = the scaled gradient vector, d**-1 * g. if g = 0, then
2293 ! step = 0 and v(stppar) = 0 are returned.
2294 ! dihdi (in) = lower triangle of the scaled hessian (approximation),
2295 ! i.e., d**-1 * h * d**-1, stored compactly by rows., i.e.,
2296 ! in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc.
2297 ! ka (i/o) = the number of hebden iterations (so far) taken to deter-
2298 ! mine step. ka .lt. 0 on input means this is the first
2299 ! attempt to determine step (for the present dig and dihdi)
2300 ! -- ka is initialized to 0 in this case. output with
2301 ! ka = 0 (or v(stppar) = 0) means step = -(h**-1)*g.
2302 ! l (i/o) = workspace of length p*(p+1)/2 for cholesky factors.
2303 ! p (in) = number of parameters -- the hessian is a p x p matrix.
2304 ! step (i/o) = the step computed.
2305 ! v (i/o) contains various constants and variables described below.
2306 ! w (i/o) = workspace of length 4*p + 6.
2308 ! *** entries in v ***
2310 ! v(dgnorm) (i/o) = 2-norm of (d**-1)*g.
2311 ! v(dstnrm) (output) = 2-norm of d*step.
2312 ! v(dst0) (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or
2313 ! overestimate of smallest eigenvalue of (d**-1)*h*(d**-1).
2314 ! v(epslon) (in) = max. rel. error allowed for psi(step). for the
2315 ! step returned, psi(step) will exceed its optimal value
2316 ! by less than -v(epslon)*psi(step). suggested value = 0.1.
2317 ! v(gtstep) (out) = inner product between g and step.
2318 ! v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step) (for pos. def.
2319 ! h only -- v(nreduc) is set to zero otherwise).
2320 ! v(phmnfc) (in) = tol. (together with v(phmxfc)) for accepting step
2321 ! (more*s sigma). the error v(dstnrm) - v(radius) must lie
2322 ! between v(phmnfc)*v(radius) and v(phmxfc)*v(radius).
2323 ! v(phmxfc) (in) (see v(phmnfc).)
2324 ! suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5.
2325 ! v(preduc) (out) = psi(step) = predicted obj. func. reduction for step.
2326 ! v(radius) (in) = radius of current (scaled) trust region.
2327 ! v(rad0) (i/o) = value of v(radius) from previous call.
2328 ! v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha
2329 ! described below under algorithm notes. if h + alpha*d**2
2330 ! (see algorithm notes) is (nearly) singular, however,
2331 ! then v(stppar) = -alpha.
2333 ! *** usage notes ***
2335 ! if it is desired to recompute step using a different value of
2336 ! v(radius), then this routine may be restarted by calling it
2337 ! with all parameters unchanged except v(radius). (this explains
2338 ! why step and w are listed as i/o). on an initial call (one with
2339 ! ka .lt. 0), step and w need not be initialized and only compo-
2340 ! nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and
2341 ! v(rad0) of v must be initialized.
2343 ! *** algorithm notes ***
2345 ! the desired g-q-t step (ref. 2, 3, 4, 6) satisfies
2346 ! (h + alpha*d**2)*step = -g for some nonnegative alpha such that
2347 ! h + alpha*d**2 is positive semidefinite. alpha and step are
2348 ! computed by a scheme analogous to the one described in ref. 5.
2349 ! estimates of the smallest and largest eigenvalues of the hessian
2350 ! are obtained from the gerschgorin circle theorem enhanced by a
2351 ! simple form of the scaling described in ref. 7. cases in which
2352 ! h + alpha*d**2 is nearly (or exactly) singular are handled by
2353 ! the technique discussed in ref. 2. in these cases, a step of
2354 ! (exact) length v(radius) is returned for which psi(step) exceeds
2355 ! its optimal value by less than -v(epslon)*psi(step). the test
2356 ! suggested in ref. 6 for detecting the special case is performed
2357 ! once two matrix factorizations have been done -- doing so sooner
2358 ! seems to degrade the performance of optimization routines that
2359 ! call this routine.
2361 ! *** functions and subroutines called ***
2363 ! dotprd - returns inner product of two vectors.
2364 ! litvmu - applies inverse-transpose of compact lower triang. matrix.
2365 ! livmul - applies inverse of compact lower triang. matrix.
2366 ! lsqrt - finds cholesky factor (of compactly stored lower triang.).
2367 ! lsvmin - returns approx. to min. sing. value of lower triang. matrix.
2368 ! rmdcon - returns machine-dependent constants.
2369 ! v2norm - returns 2-norm of a vector.
2371 ! *** references ***
2373 ! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive
2374 ! nonlinear least-squares algorithm, acm trans. math.
2375 ! software, vol. 7, no. 3.
2376 ! 2. gay, d.m. (1981), computing optimal locally constrained steps,
2377 ! siam j. sci. statist. computing, vol. 2, no. 2, pp.
2379 ! 3. goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966),
2380 ! maximization by quadratic hill-climbing, econometrica 34,
2382 ! 4. hebden, m.d. (1973), an algorithm for minimization using exact
2383 ! second derivatives, report t.p. 515, theoretical physics
2384 ! div., a.e.r.e. harwell, oxon., england.
2385 ! 5. more, j.j. (1978), the levenberg-marquardt algorithm, implemen-
2386 ! tation and theory, pp.105-116 of springer lecture notes
2387 ! in mathematics no. 630, edited by g.a. watson, springer-
2388 ! verlag, berlin and new york.
2389 ! 6. more, j.j., and sorensen, d.c. (1981), computing a trust region
2390 ! step, technical report anl-81-83, argonne national lab.
2391 ! 7. varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15,
2396 ! coded by david m. gay.
2397 ! this subroutine was written in connection with research
2398 ! supported by the national science foundation under grants
2399 ! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
2402 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2404 ! *** local variables ***
2407 integer :: dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,&
2408 j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x
2409 real(kind=8) :: alphak, aki, akk, delta, dst, eps, gtsta, lk,&
2410 oldphi, phi, phimax, phimin, psifac, rad, radsq,&
2411 root, si, sk, sw, t, twopsi, t1, t2, uk, wi
2414 real(kind=8) :: big, dgxfac !el, epsfac, four, half, kappa, negone,
2415 !el 1 one, p001, six, three, two, zero
2417 ! *** intrinsic functions ***
2419 !el real(kind=8) :: dabs, dmax1, dmin1, dsqrt
2421 ! *** external functions and subroutines ***
2423 !el external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm
2424 !el real(kind=8) :: dotprd, lsvmin, rmdcon, v2norm
2426 ! *** subscripts for v ***
2428 !el integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc,
2429 !el 1 phmnfc, phmxfc, preduc, radius, rad0
2431 ! data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/,
2432 ! 1 nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/,
2433 ! 2 rad0/9/, stppar/5/
2435 integer,parameter :: dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,&
2436 nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,&
2441 ! data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/,
2442 ! 1 kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/,
2443 ! 2 six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/
2445 real(kind=8), parameter :: epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,&
2446 kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,&
2447 six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0
2450 data big/0.d+0/, dgxfac/0.d+0/
2454 ! *** store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx).
2456 ! *** store gerschgorin over- and underestimates of the largest
2457 ! *** and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax)
2458 ! *** and w(emin) respectively.
2461 ! *** for use in recomputing step, the final values of lk, uk, dst,
2462 ! *** and the inverse derivative of more*s phi at 0 (for pos. def.
2463 ! *** h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin)
2469 ! *** store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p).
2472 ! *** store -d*step in w(q),...,w(q0+p).
2475 ! *** allocate storage for scratch vector x ***
2479 ! *** phitol = max. error allowed in dst = v(dstnrm) = 2-norm of
2481 phimax = v(phmxfc) * rad
2482 phimin = v(phmnfc) * rad
2483 psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) * &
2484 (kappa + one) + kappa + two) * rad**2)
2485 ! *** oldphi is used to detect limits of numerical accuracy. if
2486 ! *** we recompute step and it does not change, then we accept it.
2493 ! *** start or restart, depending on ka ***
2495 if (ka .ge. 0) go to 290
2497 ! *** fresh start ***
2503 v(dgnorm) = v2norm(p, dig)
2507 if (v(dgnorm) .eq. zero) kamin = 0
2509 ! *** store diag(dihdi) in w(diag0+1),...,w(diag0+p) ***
2518 ! *** determine w(dggdmx), the largest element of dihdi ***
2524 if (t1 .lt. t) t1 = t
2528 ! *** try alpha = 0 ***
2530 30 call lsqrt(1, p, l, dihdi, irc)
2531 if (irc .eq. 0) go to 50
2532 ! *** indef. h -- underestimate smallest eigenvalue, use this
2533 ! *** estimate to initialize lower bound lk on alpha.
2540 call litvmu(irc, w, l, w)
2544 if (restrt) go to 210
2547 ! *** positive definite h -- compute unmodified newton step. ***
2549 t = lsvmin(p, l, w(q), w(q))
2550 if (t .ge. one) go to 60
2551 if (big .le. zero) big = rmdcon(6)
2552 if (v(dgnorm) .ge. t*t*big) go to 70
2553 60 call livmul(p, w(q), l, dig)
2554 gtsta = dotprd(p, w(q), w(q))
2555 v(nreduc) = half * gtsta
2556 call litvmu(p, w(q), l, w(q))
2557 dst = v2norm(p, w(q))
2560 if (phi .le. phimax) go to 260
2561 if (restrt) go to 210
2563 ! *** prepare to compute gerschgorin estimates of largest (and
2564 ! *** smallest) eigenvalues. ***
2569 if (i .eq. 1) go to 90
2581 ! *** (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1) ***
2585 if (p .le. 1) go to 120
2589 if (t .ge. t1) go to 110
2601 if (i .eq. k) go to 130
2602 aki = dabs(dihdi(k1))
2605 t1 = half * (akk - w(j) + si - aki)
2606 t1 = t1 + dsqrt(t1*t1 + sk*aki)
2607 if (t .lt. t1) t = t1
2608 if (i .lt. k) go to 140
2614 uk = v(dgnorm)/rad - w(emin)
2615 if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2616 if (uk .le. zero) uk = p001
2618 ! *** compute gerschgorin (over-)estimate of largest eigenvalue ***
2622 if (p .le. 1) go to 170
2626 if (t .le. t1) go to 160
2638 if (i .eq. k) go to 180
2639 aki = dabs(dihdi(k1))
2642 t1 = half * (w(j) + si - aki - akk)
2643 t1 = t1 + dsqrt(t1*t1 + sk*aki)
2644 if (t .lt. t1) t = t1
2645 if (i .lt. k) go to 190
2651 lk = dmax1(lk, v(dgnorm)/rad - w(emax))
2653 ! *** alphak = current value of alpha (see alg. notes above). we
2654 ! *** use more*s scheme for initializing it.
2655 alphak = dabs(v(stppar)) * v(rad0)/rad
2657 if (irc .ne. 0) go to 210
2659 ! *** compute l0 for positive definite h ***
2661 call livmul(p, w, l, w(q))
2663 w(phipin) = dst / t / t
2664 lk = dmax1(lk, phi*w(phipin))
2666 ! *** safeguard alphak and add alphak*i to (d**-1)*h*(d**-1) ***
2669 if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk) &
2670 alphak = uk * dmax1(p001, dsqrt(lk/uk))
2671 if (alphak .le. zero) alphak = half * uk
2672 if (alphak .le. zero) alphak = uk
2677 dihdi(k) = w(j) + alphak
2680 ! *** try computing cholesky decomposition ***
2682 call lsqrt(1, p, l, dihdi, irc)
2683 if (irc .eq. 0) go to 240
2685 ! *** (d**-1)*h*(d**-1) + alphak*i is indefinite -- overestimate
2686 ! *** smallest eigenvalue for use in updating lk ***
2694 call litvmu(irc, w, l, w)
2696 lk = alphak - t/t1/t1
2700 ! *** alphak makes (d**-1)*h*(d**-1) positive definite.
2701 ! *** compute q = -d*step, check for convergence. ***
2703 240 call livmul(p, w(q), l, dig)
2704 gtsta = dotprd(p, w(q), w(q))
2705 call litvmu(p, w(q), l, w(q))
2706 dst = v2norm(p, w(q))
2708 if (phi .le. phimax .and. phi .ge. phimin) go to 270
2709 if (phi .eq. oldphi) go to 270
2711 if (phi .lt. zero) go to 330
2713 ! *** unacceptable alphak -- update lk, uk, alphak ***
2715 250 if (ka .ge. kalim) go to 270
2716 ! *** the following dmin1 is necessary because of restarts ***
2717 if (phi .lt. zero) uk = dmin1(uk, alphak)
2718 ! *** kamin = 0 only iff the gradient vanishes ***
2719 if (kamin .eq. 0) go to 210
2720 call livmul(p, w, l, w(q))
2722 alphak = alphak + (phi/t1) * (dst/t1) * (dst/rad)
2723 lk = dmax1(lk, alphak)
2726 ! *** acceptable step on first try ***
2730 ! *** successful step in general. compute step = -(d**-1)*q ***
2734 step(i) = -w(j)/d(i)
2737 v(preduc) = half * (dabs(alphak)*dst*dst + gtsta)
2741 ! *** restart with new radius ***
2743 290 if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310
2745 ! *** prepare to return newton step ***
2759 if (v(dgnorm) .eq. zero) kamin = 0
2760 if (ka .eq. 0) go to 50
2763 alphak = dabs(v(stppar))
2767 if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2768 if (uk .le. zero) uk = p001
2769 if (rad .gt. v(rad0)) go to 320
2771 ! *** smaller radius ***
2773 if (alphak .gt. zero) lk = w(lk0)
2774 lk = dmax1(lk, t - w(emax))
2775 if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2778 ! *** bigger radius ***
2779 320 if (alphak .gt. zero) uk = dmin1(uk, w(uk0))
2780 lk = dmax1(zero, -v(dst0), t - w(emax))
2781 if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2784 ! *** decide whether to check for special case... in practice (from
2785 ! *** the standpoint of the calling optimization code) it seems best
2786 ! *** not to check until a few iterations have failed -- hence the
2787 ! *** test on kamin below.
2789 330 delta = alphak + dmin1(zero, v(dst0))
2790 twopsi = alphak*dst*dst + gtsta
2791 if (ka .ge. kamin) go to 340
2792 ! *** if the test in ref. 2 is satisfied, fall through to handle
2793 ! *** the special case (as soon as the more-sorensen test detects
2795 if (delta .ge. psifac*twopsi) go to 370
2797 ! *** check for the special case of h + alpha*d**2 (nearly)
2798 ! *** singular. use one step of inverse power method with start
2799 ! *** from lsvmin to obtain approximate eigenvector corresponding
2800 ! *** to smallest eigenvalue of (d**-1)*h*(d**-1). lsvmin returns
2801 ! *** x and w with l*w = x.
2803 340 t = lsvmin(p, l, w(x), w)
2805 ! *** normalize w ***
2808 ! *** complete current inv. power iter. -- replace w by (l**-t)*w.
2809 call litvmu(p, w, l, w)
2810 t2 = one/v2norm(p, w)
2815 ! *** now w is the desired approximate (unit) eigenvector and
2816 ! *** t*x = ((d**-1)*h*(d**-1) + alphak*i)*w.
2818 sw = dotprd(p, w(q), w)
2819 t1 = (rad + dst) * (rad - dst)
2820 root = dsqrt(sw*sw + t1)
2821 if (sw .lt. zero) root = -root
2822 si = t1 / (sw + root)
2824 ! *** the actual test for the special case...
2826 if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380
2828 ! *** update upper bound on smallest eigenvalue (when not positive)
2829 ! *** (as recommended by more and sorensen) and continue...
2831 if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak)
2832 lk = dmax1(lk, -v(dst0))
2834 ! *** check whether we can hope to detect the special case in
2835 ! *** the available arithmetic. accept step as it is if not.
2837 ! *** if not yet available, obtain machine dependent value dgxfac.
2838 370 if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3)
2840 if (delta .gt. dgxfac*w(dggdmx)) go to 250
2843 ! *** special case detected... negate alphak to indicate special case
2845 380 alphak = -alphak
2846 v(preduc) = half * twopsi
2848 ! *** accept current step if adding si*w would lead to a
2849 ! *** further relative reduction in psi of less than v(epslon)/3.
2852 t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w)))
2853 if (t .lt. eps*twopsi/six) go to 390
2854 v(preduc) = v(preduc) + t
2859 w(j) = t1*w(i) - w(j)
2860 step(i) = w(j) / d(i)
2862 v(gtstep) = dotprd(p, dig, w(q))
2864 ! *** save values for use in a possible restart ***
2873 ! *** restore diagonal of dihdi ***
2884 ! *** last card of gqtst follows ***
2885 end subroutine gqtst
2886 !-----------------------------------------------------------------------------
2887 subroutine lsqrt(n1, n, l, a, irc)
2889 ! *** compute rows n1 through n of the cholesky factor l of
2890 ! *** a = l*(l**t), where l and the lower triangle of a are both
2891 ! *** stored compactly by rows (and may occupy the same storage).
2892 ! *** irc = 0 means all went well. irc = j means the leading
2893 ! *** principal j x j submatrix of a is not positive definite --
2894 ! *** and l(j*(j+1)/2) contains the (nonpos.) reduced j-th diagonal.
2896 ! *** parameters ***
2898 integer :: n1, n, irc
2899 !al real(kind=8) :: l(1), a(1)
2900 real(kind=8) :: l(n*(n+1)/2), a(n*(n+1)/2)
2901 ! dimension l(n*(n+1)/2), a(n*(n+1)/2)
2903 ! *** local variables ***
2905 integer :: i, ij, ik, im1, i0, j, jk, jm1, j0, k
2906 real(kind=8) :: t, td !el, zero
2908 ! *** intrinsic functions ***
2910 !el real(kind=8) :: dsqrt
2915 real(kind=8),parameter :: zero=0.d+0
2920 i0 = n1 * (n1 - 1) / 2
2923 if (i .eq. 1) go to 40
2928 if (j .eq. 1) go to 20
2937 t = (a(ij) - t) / l(j0)
2943 if (t .le. zero) go to 60
2955 ! *** last card of lsqrt ***
2956 end subroutine lsqrt
2957 !-----------------------------------------------------------------------------
2958 real(kind=8) function lsvmin(p, l, x, y)
2960 ! *** estimate smallest sing. value of packed lower triang. matrix l
2962 ! *** parameter declarations ***
2965 !al real(kind=8) :: l(1), x(p), y(p)
2966 real(kind=8) :: l(p*(p+1)/2), x(p), y(p)
2967 ! dimension l(p*(p+1)/2)
2969 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2973 ! this function returns a good over-estimate of the smallest
2974 ! singular value of the packed lower triangular matrix l.
2976 ! *** parameter description ***
2978 ! p (in) = the order of l. l is a p x p lower triangular matrix.
2979 ! l (in) = array holding the elements of l in row order, i.e.
2980 ! l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc.
2981 ! x (out) if lsvmin returns a positive value, then x is a normalized
2982 ! approximate left singular vector corresponding to the
2983 ! smallest singular value. this approximation may be very
2984 ! crude. if lsvmin returns zero, then some components of x
2985 ! are zero and the rest retain their input values.
2986 ! y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an
2987 ! unnormalized approximate right singular vector correspond-
2988 ! ing to the smallest singular value. this approximation
2989 ! may be crude. if lsvmin returns zero, then y retains its
2990 ! input value. the caller may pass the same vector for x
2991 ! and y (nonstandard fortran usage), in which case y over-
2992 ! writes x (for nonzero lsvmin returns).
2994 ! *** algorithm notes ***
2996 ! the algorithm is based on (1), with the additional provision that
2997 ! lsvmin = 0 is returned if the smallest diagonal element of l
2998 ! (in magnitude) is not more than the unit roundoff times the
2999 ! largest. the algorithm uses a random number generator proposed
3000 ! in (4), which passes the spectral test with flying colors -- see
3003 ! *** subroutines and functions called ***
3005 ! v2norm - function, returns the 2-norm of a vector.
3007 ! *** references ***
3009 ! (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977),
3010 ! an estimate for the condition number of a matrix, report
3011 ! tm-310, applied math. div., argonne national laboratory.
3013 ! (2) hoaglin, d.c. (1976), theoretical properties of congruential
3014 ! random-number generators -- an empirical view,
3015 ! memorandum ns-340, dept. of statistics, harvard univ.
3017 ! (3) knuth, d.e. (1969), the art of computer programming, vol. 2
3018 ! (seminumerical algorithms), addison-wesley, reading, mass.
3020 ! (4) smith, c.s. (1971), multiplicative pseudo-random number
3021 ! generators with prime modulus, j. assoc. comput. mach. 18,
3026 ! designed and coded by david m. gay (winter 1977/summer 1978).
3030 ! this subroutine was written in connection with research
3031 ! supported by the national science foundation under grants
3032 ! mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989.
3034 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3036 ! *** local variables ***
3038 integer :: i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1
3039 real(kind=8) :: b, sminus, splus, t, xminus, xplus
3043 !el real(kind=8) :: half, one, r9973, zero
3045 ! *** intrinsic functions ***
3049 !el real(kind=8) :: dabs
3051 ! *** external functions and subroutines ***
3053 !el external dotprd, v2norm, vaxpy
3054 !el real(kind=8) :: dotprd, v2norm
3057 ! data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/
3059 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0
3067 ! *** first check whether to return lsvmin = 0 and initialize x ***
3072 if (l(jj) .eq. zero) go to 110
3073 ix = mod(3432*ix, 9973)
3074 b = half*(one + float(ix)/r9973)
3077 if (p .le. 1) go to 60
3080 if (l(ii) .eq. zero) go to 110
3082 x(i) = xplus * l(ji)
3085 ! *** solve (l**t)*x = b, where the components of b have randomly
3086 ! *** chosen magnitudes in (.5,1) with signs chosen to make x large.
3088 ! do j = p-1 to 1 by -1...
3091 ! *** determine x(j) in this iteration. note for i = 1,2,...,j
3092 ! *** that x(i) holds the current partial sum for row i.
3093 ix = mod(3432*ix, 9973)
3094 b = half*(one + float(ix)/r9973)
3096 xminus = (-b - x(j))
3098 sminus = dabs(xminus)
3103 xminus = xminus/l(jj)
3104 if (jm1 .eq. 0) go to 30
3107 splus = splus + dabs(x(i) + l(ji)*xplus)
3108 sminus = sminus + dabs(x(i) + l(ji)*xminus)
3110 30 if (sminus .gt. splus) xplus = xminus
3112 ! *** update partial sums ***
3113 if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x)
3116 ! *** normalize x ***
3118 60 t = one/v2norm(p, x)
3122 ! *** solve l*y = x and return lsvmin = 1/twonorm(y) ***
3129 if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y)
3130 y(j) = (x(j) - t) / l(jj)
3133 lsvmin = one/v2norm(p, y)
3138 ! *** last card of lsvmin follows ***
3140 !-----------------------------------------------------------------------------
3141 subroutine slvmul(p, y, s, x)
3143 ! *** set y = s * x, s = p x p symmetric matrix. ***
3144 ! *** lower triangle of s stored rowwise. ***
3146 ! *** parameter declarations ***
3149 !al real(kind=8) :: s(1), x(p), y(p)
3150 real(kind=8) :: s(p*(p+1)/2), x(p), y(p)
3151 ! dimension s(p*(p+1)/2)
3153 ! *** local variables ***
3155 integer :: i, im1, j, k
3158 ! *** no intrinsic functions ***
3160 ! *** external function ***
3163 !el real(kind=8) :: dotprd
3165 !-----------------------------------------------------------------------
3169 y(i) = dotprd(i, s(j), x)
3173 if (p .le. 1) go to 999
3180 y(k) = y(k) + s(j)*xi
3186 ! *** last card of slvmul follows ***
3187 end subroutine slvmul
3188 !-----------------------------------------------------------------------------
3190 !-----------------------------------------------------------------------------
3191 subroutine minimize(etot,x,iretcode,nfun)
3198 use energy, only: funcgrad,etotal,enerprint
3200 use energy, only: func,gradient,fdum!,etotal,enerprint
3203 ! implicit real*8 (a-h,o-z)
3204 ! include 'DIMENSIONS'
3205 integer,parameter :: liv=60
3206 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3207 !********************************************************************
3208 ! OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
3209 ! the calling subprogram. *
3210 ! when d(i)=1.0, then v(35) is the length of the initial step, *
3211 ! calculated in the usual pythagorean way. *
3212 ! absolute convergence occurs when the function is within v(31) of *
3213 ! zero. unless you know the minimum value in advance, abs convg *
3214 ! is probably not useful. *
3215 ! relative convergence is when the model predicts that the function *
3216 ! will decrease by less than v(32)*abs(fun). *
3217 !********************************************************************
3218 ! include 'COMMON.IOUNITS'
3219 ! include 'COMMON.VAR'
3220 ! include 'COMMON.GEO'
3221 ! include 'COMMON.MINIM'
3223 !el common /srutu/ icall
3224 integer,dimension(liv) :: iv
3225 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
3226 !el real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3227 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3228 real(kind=8) :: energia(0:n_ene),grdmin
3229 ! external func,gradient,fdum
3230 ! external func_restr,grad_restr
3231 logical :: not_done,change,reduce
3232 !el common /przechowalnia/ v
3234 integer :: iretcode,nfun,nvar_restr,idum(1),j
3236 real(kind=8) :: etot,rdum(1) !,fdum
3240 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3242 if (.not.allocated(v)) allocate(v(1:lv))
3246 coordtype='RIGIDBODY'
3249 jprint=print_min_stat
3251 if (.not. allocated(scale)) allocate (scale(nvar))
3253 !c set scaling parameter for function and derivative values;
3254 !c use square root of median eigenvalue of typical Hessian
3266 !c write (iout,*) "Calling lbfgs"
3267 write (iout,*) 'Calling LBFGS, minimization in angles'
3268 call var_to_geom(nvar,x)
3270 call etotal(energia(0))
3271 call enerprint(energia(0))
3272 call lbfgs (nvar,x,etot,grdmin,funcgrad,optsave)
3274 write (iout,*) "Minimized energy",etot
3281 ! DO WHILE (NOT_DONE)
3283 call deflt(2,iv,liv,lv,v)
3284 ! 12 means fresh start, dont call deflt
3286 ! max num of fun calls
3287 if (maxfun.eq.0) maxfun=500
3289 ! max num of iterations
3290 if (maxmin.eq.0) maxmin=1000
3294 ! selects output unit
3296 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
3297 ! 1 means to print out result
3298 iv(22)=print_min_res
3299 ! 1 means to print out summary stats
3300 iv(23)=print_min_stat
3301 ! 1 means to print initial x and d
3302 iv(24)=print_min_ini
3303 ! min val for v(radfac) default is 0.1
3305 ! max val for v(radfac) default is 4.0
3308 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3309 ! the sumsl default is 0.1
3311 ! false conv if (act fnctn decrease) .lt. v(34)
3312 ! the sumsl default is 100*machep
3314 ! absolute convergence
3315 if (tolf.eq.0.0D0) tolf=1.0D-4
3317 ! relative convergence
3318 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3320 ! controls initial step size
3322 ! large vals of d correspond to small components of step
3329 !d print *,'Calling SUMSL'
3330 ! call var_to_geom(nvar,x)
3332 ! call etotal(energia(0))
3336 call x2xx(x,xx,nvar_restr)
3337 call sumsl(nvar_restr,d,xx,func_restr,grad_restr,&
3338 iv,liv,lv,v,idum,rdum,fdum)
3341 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
3345 !d print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
3346 !d write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
3349 call var_to_geom(nvar,x)
3351 ! write (iout,'(a)') 'Reduction worked, minimizing again...'
3355 write(iout,*) 'Warning calling chainbuild'
3358 !el---------------------
3359 ! write (iout,'(/a)') &
3360 ! "Cartesian coordinates of the reference structure after SUMSL"
3361 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3362 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3364 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
3365 ! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
3366 ! (c(j,i+nres),j=1,3)
3368 !el----------------------------
3369 ! call etotal(energia) !sp
3371 ! call enerprint(energia) !sp
3374 ! write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
3379 end subroutine minimize
3380 !-----------------------------------------------------------------------------
3382 !-----------------------------------------------------------------------------
3384 subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
3386 use energy, only: cartder,zerograd,etotal,sum_gradient
3387 ! implicit real*8 (a-h,o-z)
3388 ! include 'DIMENSIONS'
3389 ! include 'COMMON.CHAIN'
3390 ! include 'COMMON.DERIV'
3391 ! include 'COMMON.VAR'
3392 ! include 'COMMON.INTERACT'
3393 ! include 'COMMON.FFIELD'
3394 ! include 'COMMON.IOUNITS'
3396 integer :: uiparm(1)
3397 real(kind=8) :: urparm(1)
3398 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
3399 integer :: n,nf,ig,ind,i,j,ij,k,igall
3400 real(kind=8) :: f,gphii,gthetai,galphai,gomegai
3401 real(kind=8),external :: ufparm
3404 if (nf-nfl+1) 20,30,40
3405 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
3406 ! write (iout,*) 'grad 20'
3411 ! Intercept NaNs in the coordinates
3412 ! write(iout,*) (var(i),i=1,nvar)
3417 if (x_sum.ne.x_sum) then
3418 write(iout,*)" *** grad_restr : Found NaN in coordinates"
3420 print *," *** grad_restr : Found NaN in coordinates"
3424 call var_to_geom_restr(n,x)
3425 write(iout,*) 'Warning calling chainbuild'
3428 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
3432 ! Convert the Cartesian gradient into internal-coordinate gradient.
3438 IF (mask_phi(i+2).eq.1) THEN
3443 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
3444 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
3457 IF (mask_theta(i+2).eq.1) THEN
3463 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
3464 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
3474 if (itype(i,1).ne.10) then
3475 IF (mask_side(i).eq.1) THEN
3479 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
3488 if (itype(i,1).ne.10) then
3489 IF (mask_side(i).eq.1) THEN
3493 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
3501 ! Add the components corresponding to local energy terms.
3508 if (mask_phi(i).eq.1) then
3510 g(ig)=g(ig)+gloc(igall,icg)
3516 if (mask_theta(i).eq.1) then
3518 g(ig)=g(ig)+gloc(igall,icg)
3524 if (itype(i,1).ne.10) then
3526 if (mask_side(i).eq.1) then
3528 g(ig)=g(ig)+gloc(igall,icg)
3535 !d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
3538 end subroutine grad_restr
3541 !-----------------------------------------------------------------------------
3542 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F
3545 use energy, only: zerograd,etotal,sum_gradient
3546 ! implicit real*8 (a-h,o-z)
3547 ! include 'DIMENSIONS'
3548 ! include 'COMMON.DERIV'
3549 ! include 'COMMON.IOUNITS'
3550 ! include 'COMMON.GEO'
3553 !el common /chuju/ jjj
3554 real(kind=8) :: energia(0:n_ene)
3556 real(kind=8),external :: ufparm
3557 integer :: uiparm(1)
3558 real(kind=8) :: urparm(1)
3559 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3560 ! if (jjj.gt.0) then
3561 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
3565 call var_to_geom_restr(n,x)
3567 write(iout,*) 'Warning calling chainbuild'
3569 !d write (iout,*) 'ETOTAL called from FUNC'
3570 call etotal(energia)
3573 ! if (jjj.gt.0) then
3574 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
3575 ! write (iout,*) 'f=',etot
3579 end subroutine func_restr
3581 !-----------------------------------------------------------------------------
3582 ! subroutine func(n,x,nf,f,uiparm,urparm,ufparm) in module energy
3583 !-----------------------------------------------------------------------------
3584 subroutine x2xx(x,xx,n)
3586 ! implicit real*8 (a-h,o-z)
3587 ! include 'DIMENSIONS'
3588 ! include 'COMMON.VAR'
3589 ! include 'COMMON.CHAIN'
3590 ! include 'COMMON.INTERACT'
3591 integer :: n,i,ij,ig,igall
3592 real(kind=8),dimension(6*nres) :: xx,x !(maxvar) (maxvar=6*maxres)
3594 !el allocate(varall(nvar)) allocated in alioc_ener_arrays
3604 if (mask_phi(i).eq.1) then
3612 if (mask_theta(i).eq.1) then
3620 if (itype(i,1).ne.10) then
3622 if (mask_side(i).eq.1) then
3634 !-----------------------------------------------------------------------------
3635 !el subroutine xx2x(x,xx) in module math
3636 !-----------------------------------------------------------------------------
3637 subroutine minim_dc(etot,iretcode,nfun)
3647 use energy, only: fdum,check_ecartint,etotal,enerprint
3648 ! implicit real*8 (a-h,o-z)
3649 ! include 'DIMENSIONS'
3653 integer,parameter :: liv=60
3654 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3655 ! include 'COMMON.SETUP'
3656 ! include 'COMMON.IOUNITS'
3657 ! include 'COMMON.VAR'
3658 ! include 'COMMON.GEO'
3659 ! include 'COMMON.MINIM'
3660 ! include 'COMMON.CHAIN'
3661 integer :: iretcode,nfun,k,i,j,idum(1)
3662 integer :: lv,nvarx,nf
3663 integer,dimension(liv) :: iv
3664 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
3665 real(kind=8),dimension(:),allocatable :: x,d,xx !(maxvar) (maxvar=6*maxres)
3666 !el common /przechowalnia/ v
3668 real(kind=8) :: energia(0:n_ene),grdmin
3669 ! external func_dc,grad_dc ,fdum
3670 logical :: not_done,change,reduce
3671 real(kind=8) :: g(6*nres),f1,etot,rdum(1) !,fdum
3675 coordtype='CARTESIAN'
3678 jprint=print_min_stat
3680 write(iout,*) "in minim_dc LBFGS"
3683 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3684 print *,"lv value",lv
3685 if (.not. allocated(v)) allocate(v(1:lv))
3687 call deflt(2,iv,liv,lv,v)
3688 print *,"after delft"
3689 ! 12 means fresh start, dont call deflt
3691 ! max num of fun calls
3692 if (maxfun.eq.0) maxfun=500
3694 ! max num of iterations
3695 if (maxmin.eq.0) maxmin=1000
3699 ! selects output unit
3701 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
3702 ! 1 means to print out result
3703 iv(22)=print_min_res
3704 ! 1 means to print out summary stats
3705 iv(23)=print_min_stat
3706 ! 1 means to print initial x and d
3707 iv(24)=print_min_ini
3708 ! min val for v(radfac) default is 0.1
3710 ! max val for v(radfac) default is 4.0
3713 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3714 ! the sumsl default is 0.1
3716 ! false conv if (act fnctn decrease) .lt. v(34)
3717 ! the sumsl default is 100*machep
3719 ! absolute convergence
3720 if (tolf.eq.0.0D0) tolf=1.0D-4
3722 ! relative convergence
3723 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3725 ! controls initial step size
3727 ! large vals of d correspond to small components of step
3732 if (.not.allocated(x)) then
3734 allocate(xx(6*nres))
3745 if (ialph(i,1).gt.0) then
3753 call chainbuild_cart
3754 call etotal(energia(0))
3755 call enerprint(energia(0))
3760 !c perform dynamic allocation of some global arrays
3762 if (.not. allocated(scale)) allocate (scale(nvarx))
3764 !c set scaling parameter for function and derivative values;
3765 !c use square root of median eigenvalue of typical Hessian
3777 write (iout,*) "minim_dc Calling lbfgs"
3778 call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
3781 ! call check_ecartint
3782 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
3783 ! call check_ecartint
3793 if (ialph(i,1).gt.0) then
3800 call chainbuild_cart
3804 !d call func_dc(k,x,nf,f,idum,rdum,fdum)
3805 !d call grad_dc(k,x,nf,g,idum,rdum,fdum)
3809 !d call func_dc(k,x,nf,f1,idum,rdum,fdum)
3811 !d print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
3813 !el---------------------
3814 ! write (iout,'(/a)') &
3815 ! "Cartesian coordinates of the reference structure after SUMSL"
3816 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3817 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3819 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
3820 ! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
3821 ! (c(j,i+nres),j=1,3)
3823 !el----------------------------
3830 end subroutine minim_dc
3834 real(kind=8) function funcgrad_dc(x,g)
3835 use energy, only: etotal, zerograd,cartgrad
3836 ! implicit real*8 (a-h,o-z)
3841 real(kind=8),dimension(nres*6):: x,g
3842 double precision energia(0:n_ene)
3854 if (ialph(i,1).gt.0) then
3861 call chainbuild_cart
3863 call etotal(energia(0))
3864 !c write (iout,*) "energia",energia(0)
3865 funcgrad_dc=energia(0)
3875 if (ialph(i,1).gt.0) then
3886 !-----------------------------------------------------------------------------
3887 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
3890 use energy, only: zerograd,etotal
3891 ! implicit real*8 (a-h,o-z)
3892 ! include 'DIMENSIONS'
3896 ! include 'COMMON.SETUP'
3897 ! include 'COMMON.DERIV'
3898 ! include 'COMMON.IOUNITS'
3899 ! include 'COMMON.GEO'
3900 ! include 'COMMON.CHAIN'
3901 ! include 'COMMON.VAR'
3902 integer :: n,nf,k,i,j
3903 real(kind=8) :: energia(0:n_ene)
3904 real(kind=8),external :: ufparm
3905 integer :: uiparm(1)
3906 real(kind=8) :: urparm(1)
3907 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3910 !bad icg=mod(nf,2)+1
3921 if (ialph(i,1).gt.0) then
3928 call chainbuild_cart
3931 call etotal(energia)
3934 !d print *,'func_dc ',nf,nfl,f
3937 end subroutine func_dc
3938 !-----------------------------------------------------------------------------
3939 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
3942 use energy, only: cartgrad,zerograd,etotal
3944 ! implicit real*8 (a-h,o-z)
3945 ! include 'DIMENSIONS'
3949 ! include 'COMMON.SETUP'
3950 ! include 'COMMON.CHAIN'
3951 ! include 'COMMON.DERIV'
3952 ! include 'COMMON.VAR'
3953 ! include 'COMMON.INTERACT'
3954 ! include 'COMMON.FFIELD'
3955 ! include 'COMMON.MD'
3956 ! include 'COMMON.IOUNITS'
3957 real(kind=8),external :: ufparm
3958 integer :: n,nf,i,j,k
3959 integer :: uiparm(1)
3960 real(kind=8) :: urparm(1)
3961 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
3964 !elwrite(iout,*) "jestesmy w grad dc"
3966 !bad icg=mod(nf,2)+1
3968 !d print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
3969 !elwrite(iout,*) "jestesmy w grad dc nf-nfl+1", nf-nfl+1
3970 if (nf-nfl+1) 20,30,40
3971 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
3985 if (ialph(i,1).gt.0) then
3992 !elwrite(iout,*) "jestesmy w grad dc"
3993 call chainbuild_cart
3996 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
4000 !elwrite(iout,*) "jestesmy w grad dc"
4009 if (ialph(i,1).gt.0) then
4016 !elwrite(iout,*) "jestesmy w grad dc"
4019 end subroutine grad_dc
4021 !-----------------------------------------------------------------------------
4023 !-----------------------------------------------------------------------------
4025 subroutine minim_mcmf
4030 use energy, only: func,gradient,fdum
4031 ! implicit real*8 (a-h,o-z)
4032 ! include 'DIMENSIONS'
4034 integer,parameter :: liv=60
4035 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4036 ! include 'COMMON.VAR'
4037 ! include 'COMMON.IOUNITS'
4038 ! include 'COMMON.MINIM'
4039 ! real(kind=8) :: fdum
4040 ! external func,gradient,fdum
4041 !el real(kind=4) :: ran1,ran2,ran3
4042 ! include 'COMMON.SETUP'
4043 ! include 'COMMON.GEO'
4044 ! include 'COMMON.CHAIN'
4045 ! include 'COMMON.FFIELD'
4046 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
4047 real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg
4048 real(kind=8),dimension(6*nres) :: d,garbage !(maxvar) (maxvar=6*maxres)
4049 !el real(kind=8) :: v(1:77+(6*nres)*(6*nres+17)/2+1)
4050 integer,dimension(6) :: indx
4051 integer,dimension(liv) :: iv
4052 integer :: idum(1),nf !
4054 real(kind=8) :: rdum(1)
4055 real(kind=8) :: przes(3),obrot(3,3),eee
4058 integer,dimension(MPI_STATUS_SIZE) :: muster
4060 integer :: ichuj,i,ierr
4061 real(kind=8) :: rad,ene0
4062 data rad /1.745329252d-2/
4063 !el common /przechowalnia/ v
4065 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4066 if (.not. allocated(v)) allocate(v(1:lv))
4071 call mpi_recv(indx,6,mpi_integer,king,idint,CG_COMM,&
4073 if (indx(1).eq.0) return
4074 ! print *, 'worker ',me,' received order ',n,ichuj
4075 call mpi_recv(var,nvar,mpi_double_precision,&
4076 king,idreal,CG_COMM,muster,ierr)
4077 call mpi_recv(ene0,1,mpi_double_precision,&
4078 king,idreal,CG_COMM,muster,ierr)
4079 ! print *, 'worker ',me,' var read '
4082 call deflt(2,iv,liv,lv,v)
4083 ! 12 means fresh start, dont call deflt
4085 ! max num of fun calls
4086 if (maxfun.eq.0) maxfun=500
4088 ! max num of iterations
4089 if (maxmin.eq.0) maxmin=1000
4093 ! selects output unit
4096 ! 1 means to print out result
4098 ! 1 means to print out summary stats
4100 ! 1 means to print initial x and d
4102 ! min val for v(radfac) default is 0.1
4104 ! max val for v(radfac) default is 4.0
4106 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
4107 ! the sumsl default is 0.1
4109 ! false conv if (act fnctn decrease) .lt. v(34)
4110 ! the sumsl default is 100*machep
4112 ! absolute convergence
4113 if (tolf.eq.0.0D0) tolf=1.0D-4
4115 ! relative convergence
4116 if (rtolf.eq.0.0D0) rtolf=1.0D-4
4118 ! controls initial step size
4120 ! large vals of d correspond to small components of step
4129 call func(nvar,var,nf,eee,idum,rdum,fdum)
4130 if(eee.gt.1.0d18) then
4131 ! print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
4132 ! print *,' energy before SUMSL =',eee
4133 ! print *,' aborting local minimization'
4140 call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
4141 ! find which conformation was returned from sumsl
4144 ! total # of ftn evaluations (for iwf=0, it includes all minimizations).
4149 call mpi_send(indx,6,mpi_integer,king,idint,CG_COMM,&
4151 ! print '(a5,i3,15f10.5)', 'ENEX0',indx(1),v(10)
4152 call mpi_send(var,nvar,mpi_double_precision,&
4153 king,idreal,CG_COMM,ierr)
4154 call mpi_send(eee,1,mpi_double_precision,king,idreal,&
4156 call mpi_send(ene0,1,mpi_double_precision,king,idreal,&
4161 end subroutine minim_mcmf
4163 !-----------------------------------------------------------------------------
4165 !-----------------------------------------------------------------------------
4166 ! algorithm 611, collected algorithms from acm.
4167 ! algorithm appeared in acm-trans. math. software, vol.9, no. 4,
4168 ! dec., 1983, p. 503-524.
4169 integer function imdcon(k)
4173 ! *** return integer machine-dependent constants ***
4175 ! *** k = 1 means return standard output unit number. ***
4176 ! *** k = 2 means return alternate output unit number. ***
4177 ! *** k = 3 means return input unit number. ***
4178 ! (note -- k = 2, 3 are used only by test programs.)
4180 ! +++ port version follows...
4184 ! data mdperm(1)/2/, mdperm(2)/4/, mdperm(3)/1/
4185 ! imdcon = i1mach(mdperm(k))
4186 ! +++ end of port version +++
4188 ! +++ non-port version follows...
4190 data mdcon(1)/6/, mdcon(2)/8/, mdcon(3)/5/
4192 ! +++ end of non-port version +++
4195 ! *** last card of imdcon follows ***
4197 !-----------------------------------------------------------------------------
4198 real(kind=8) function rmdcon(k)
4200 ! *** return machine dependent constants used by nl2sol ***
4202 ! +++ comments below contain data statements for various machines. +++
4203 ! +++ to convert to another machine, place a c in column 1 of the +++
4204 ! +++ data statement line(s) that correspond to the current machine +++
4205 ! +++ and remove the c from column 1 of the data statement line(s) +++
4206 ! +++ that correspond to the new machine. +++
4210 ! *** the constant returned depends on k...
4212 ! *** k = 1... smallest pos. eta such that -eta exists.
4213 ! *** k = 2... square root of eta.
4214 ! *** k = 3... unit roundoff = smallest pos. no. machep such
4215 ! *** that 1 + machep .gt. 1 .and. 1 - machep .lt. 1.
4216 ! *** k = 4... square root of machep.
4217 ! *** k = 5... square root of big (see k = 6).
4218 ! *** k = 6... largest machine no. big such that -big exists.
4220 real(kind=8) :: big, eta, machep
4221 integer :: bigi(4), etai(4), machei(4)
4223 !el real(kind=8) :: dsqrt
4225 equivalence (big,bigi(1)), (eta,etai(1)), (machep,machei(1))
4227 ! +++ ibm 360, ibm 370, or xerox +++
4229 ! data big/z7fffffffffffffff/, eta/z0010000000000000/,
4230 ! 1 machep/z3410000000000000/
4232 ! +++ data general +++
4234 ! data big/0.7237005577d+76/, eta/0.5397605347d-78/,
4235 ! 1 machep/2.22044605d-16/
4239 ! data big/1.7d+38/, eta/2.938735878d-39/, machep/2.775557562d-17/
4243 ! data big/1.157920892d+77/, eta/8.636168556d-78/,
4244 ! 1 machep/5.551115124d-17/
4248 ! data big/1.69d+38/, eta/5.9d-39/, machep/2.1680435d-19/
4252 ! data big/"377777100000000000000000/,
4253 ! 1 eta/"002400400000000000000000/,
4254 ! 2 machep/"104400000000000000000000/
4258 ! data big/o0777777777777777,o7777777777777777/,
4259 ! 1 eta/o1771000000000000,o7770000000000000/,
4260 ! 2 machep/o1451000000000000,o0000000000000000/
4262 ! +++ control data +++
4264 ! data big/37767777777777777777b,37167777777777777777b/,
4265 ! 1 eta/00014000000000000000b,00000000000000000000b/,
4266 ! 2 machep/15614000000000000000b,15010000000000000000b/
4270 ! data big/1.0d+9786/, eta/1.0d-9860/, machep/1.4210855d-14/
4274 ! data big/8.988d+307/, eta/1.2d-308/, machep/1.734723476d-18/
4278 data big/1.7d+38/, eta/2.939d-39/, machep/1.3877788d-17/
4282 ! data bigi(1)/577767777777777777777b/,
4283 ! 1 bigi(2)/000007777777777777776b/,
4284 ! 2 etai(1)/200004000000000000000b/,
4285 ! 3 etai(2)/000000000000000000000b/,
4286 ! 4 machei(1)/377224000000000000000b/,
4287 ! 5 machei(2)/000000000000000000000b/
4289 ! +++ port library -- requires more than just a data statement... +++
4292 ! double precision d1mach, zero
4293 ! data big/0.d+0/, eta/0.d+0/, machep/0.d+0/, zero/0.d+0/
4294 ! if (big .gt. zero) go to 1
4297 ! machep = d1mach(4)
4300 ! +++ end of port +++
4302 !------------------------------- body --------------------------------
4304 go to (10, 20, 30, 40, 50, 60), k
4309 20 rmdcon = dsqrt(256.d+0*eta)/16.d+0
4315 40 rmdcon = dsqrt(machep)
4318 50 rmdcon = dsqrt(big/256.d+0)*16.d+0
4324 ! *** last card of rmdcon follows ***
4326 !-----------------------------------------------------------------------------
4328 !-----------------------------------------------------------------------------
4329 subroutine sc_move(n_start,n_end,n_maxtry,e_drop,n_fun,etot)
4332 use random, only: iran_num
4333 use energy, only: esc
4334 ! Perform a quick search over side-chain arrangments (over
4335 ! residues n_start to n_end) for a given (frozen) CA trace
4336 ! Only side-chains are minimized (at most n_maxtry times each),
4338 ! Stops if energy drops by e_drop, otherwise tries all residues
4339 ! in the given range
4340 ! If there is an energy drop, full minimization may be useful
4341 ! n_start, n_end CAN be modified by this routine, but only if
4342 ! out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
4343 ! NOTE: this move should never increase the energy
4347 ! implicit real*8 (a-h,o-z)
4348 ! include 'DIMENSIONS'
4350 ! include 'COMMON.GEO'
4351 ! include 'COMMON.VAR'
4352 ! include 'COMMON.HEADER'
4353 ! include 'COMMON.IOUNITS'
4354 ! include 'COMMON.CHAIN'
4355 ! include 'COMMON.FFIELD'
4357 ! External functions
4358 !el integer iran_num
4359 !el external iran_num
4362 integer :: n_start,n_end,n_maxtry
4363 real(kind=8) :: e_drop
4367 real(kind=8) :: etot
4370 ! real(kind=8) :: energy(0:n_ene)
4371 real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1)
4372 real(kind=8) :: orig_e,cur_e
4373 integer :: n,n_steps,n_first,n_cur,n_tot !,i
4374 real(kind=8) :: orig_w(0:n_ene)
4375 real(kind=8) :: wtime
4377 !elwrite(iout,*) "in sc_move etot= ", etot
4378 ! Set non side-chain weights to zero (minimization is faster)
4379 ! NOTE: e(2) does not actually depend on the side-chain, only CA
4408 ! Make sure n_start, n_end are within proper range
4409 if (n_start.lt.2) n_start=2
4410 if (n_end.gt.nres-1) n_end=nres-1
4411 !rc if (n_start.lt.n_end) then
4412 if (n_start.gt.n_end) then
4417 ! Save the initial values of energy and coordinates
4419 !d call etotal(energy)
4420 !d write (iout,*) 'start sc ene',energy(0)
4421 !d call enerprint(energy(0))
4427 !rc cur_alph(i)=alph(i)
4428 !rc cur_omeg(i)=omeg(i)
4431 !t wtime=MPI_WTIME()
4432 ! Try (one by one) all specified residues, starting from a
4433 ! random position in sequence
4434 ! Stop early if the energy has decreased by at least e_drop
4435 n_tot=n_end-n_start+1
4436 n_first=iran_num(0,n_tot-1)
4439 !rc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
4440 do while (n.lt.n_tot)
4441 n_cur=n_start+mod(n_first+n,n_tot)
4442 call single_sc_move(n_cur,n_maxtry,e_drop,&
4444 !elwrite(iout,*) "after msingle sc_move etot= ", etot
4445 ! If a lower energy was found, update the current structure...
4446 !rc if (etot.lt.cur_e) then
4449 !rc cur_alph(i)=alph(i)
4450 !rc cur_omeg(i)=omeg(i)
4453 ! ...else revert to the previous one
4456 !rc alph(i)=cur_alph(i)
4457 !rc omeg(i)=cur_omeg(i)
4463 !d call etotal(energy)
4464 !d print *,'running',n,energy(0)
4468 !d call etotal(energy)
4469 !d write (iout,*) 'end sc ene',energy(0)
4471 ! Put the original weights back to calculate the full energy
4487 !t write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
4489 end subroutine sc_move
4490 !-----------------------------------------------------------------------------
4491 subroutine single_sc_move(res_pick,n_maxtry,e_drop,n_steps,n_fun,e_sc)
4493 ! Perturb one side-chain (res_pick) and minimize the
4494 ! neighbouring region, keeping all CA's and non-neighbouring
4496 ! Try until e_drop energy improvement is achieved, or n_maxtry
4497 ! attempts have been made
4498 ! At the start, e_sc should contain the side-chain-only energy(0)
4499 ! nsteps and nfun for this move are ADDED to n_steps and n_fun
4501 use energy, only: esc
4502 use geometry, only:dist
4504 ! implicit real*8 (a-h,o-z)
4505 ! include 'DIMENSIONS'
4506 ! include 'COMMON.VAR'
4507 ! include 'COMMON.INTERACT'
4508 ! include 'COMMON.CHAIN'
4509 ! include 'COMMON.MINIM'
4510 ! include 'COMMON.FFIELD'
4511 ! include 'COMMON.IOUNITS'
4513 ! External functions
4514 !el double precision dist
4518 integer :: res_pick,n_maxtry
4519 real(kind=8) :: e_drop
4521 ! Input/Output arguments
4522 integer :: n_steps,n_fun
4523 real(kind=8) :: e_sc
4528 integer :: nres_moved
4529 integer :: iretcode,loc_nfun,orig_maxfun,n_try
4530 real(kind=8) :: sc_dist,sc_dist_cutoff
4531 ! real(kind=8) :: energy_(0:n_ene)
4532 real(kind=8) :: evdw,escloc,orig_e,cur_e
4533 real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1)
4534 real(kind=8) :: var(6*nres) !(maxvar) (maxvar=6*maxres)
4536 real(kind=8) :: orig_theta(1:nres),orig_phi(1:nres),&
4537 orig_alph(1:nres),orig_omeg(1:nres)
4539 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4540 ! Define what is meant by "neighbouring side-chain"
4541 sc_dist_cutoff=5.0D0
4543 ! Don't do glycine or ends
4545 if (i.eq.10 .or. i.eq.ntyp1 .or. molnum(res_pick).eq.5) return
4547 ! Freeze everything (later will relax only selected side-chains)
4555 ! Find the neighbours of the side-chain to move
4556 ! and save initial variables
4561 ! Don't do glycine (itype(j,1)==10)
4562 if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1) &
4563 .and.(molnum(i).ne.5)) then
4564 sc_dist=dist(nres+i,nres+res_pick)
4566 sc_dist=sc_dist_cutoff
4568 if (sc_dist.lt.sc_dist_cutoff) then
4569 nres_moved=nres_moved+1
4575 ! write(iout,*) 'Warning calling chainbuild'
4577 ! write(iout,*) "before egb1"
4580 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4581 !elwrite(iout,*) "in sinle wsc=",wsc,"evdw",evdw,"wscloc",wscloc,"escloc",escloc
4582 e_sc=wsc*evdw+wscloc*escloc
4583 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4584 !d call etotal(energy)
4585 !d print *,'new ',(energy(k),k=0,n_ene)
4590 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
4591 ! Move the selected residue (don't worry if it fails)
4592 call gen_side(iabs(itype(res_pick,molnum(res_pick))),theta(res_pick+1),&
4593 alph(res_pick),omeg(res_pick),fail,molnum(res_pick))
4595 ! Minimize the side-chains starting from the new arrangement
4596 call geom_to_var(nvar,var)
4601 !rc orig_theta(i)=theta(i)
4602 !rc orig_phi(i)=phi(i)
4603 !rc orig_alph(i)=alph(i)
4604 !rc orig_omeg(i)=omeg(i)
4607 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4608 call minimize_sc1(e_sc,var,iretcode,loc_nfun)
4610 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4611 !v write(*,'(2i3,2f12.5,2i3)')
4612 !v & res_pick,nres_moved,orig_e,e_sc-cur_e,
4613 !v & iretcode,loc_nfun
4615 !$$$ if (iretcode.eq.8) then
4616 !$$$ write(iout,*)'Coordinates just after code 8'
4617 !$$$ call chainbuild
4618 !$$$ call all_varout
4619 !$$$ call flush(iout)
4621 !$$$ theta(i)=orig_theta(i)
4622 !$$$ phi(i)=orig_phi(i)
4623 !$$$ alph(i)=orig_alph(i)
4624 !$$$ omeg(i)=orig_omeg(i)
4626 !$$$ write(iout,*)'Coordinates just before code 8'
4627 !$$$ call chainbuild
4628 !$$$ call all_varout
4629 !$$$ call flush(iout)
4632 n_fun=n_fun+loc_nfun
4634 call var_to_geom(nvar,var)
4636 ! If a lower energy was found, update the current structure...
4637 if (e_sc.lt.cur_e) then
4639 !v call etotal(energy)
4642 !d e_sc1=wsc*evdw+wscloc*escloc
4643 !d print *,' new',e_sc1,energy(0)
4644 !v print *,'new ',energy(0)
4645 !d call enerprint(energy(0))
4648 if (mask_side(i).eq.1) then
4654 ! ...else revert to the previous one
4657 if (mask_side(i).eq.1) then
4666 n_steps=n_steps+n_try
4668 ! Reset the minimization mask_r to false
4672 end subroutine single_sc_move
4673 !-----------------------------------------------------------------------------
4674 subroutine sc_minimize(etot,iretcode,nfun)
4676 ! Minimizes side-chains only, leaving backbone frozen
4678 use energy, only: etotal
4680 ! implicit real*8 (a-h,o-z)
4681 ! include 'DIMENSIONS'
4682 ! include 'COMMON.VAR'
4683 ! include 'COMMON.CHAIN'
4684 ! include 'COMMON.FFIELD'
4687 real(kind=8) :: etot
4688 integer :: iretcode,nfun
4692 real(kind=8) :: orig_w(0:n_ene),energy_(0:n_ene)
4693 real(kind=8) :: var(6*nres) !(maxvar)(maxvar=6*maxres)
4696 ! Set non side-chain weights to zero (minimization is faster)
4697 ! NOTE: e(2) does not actually depend on the side-chain, only CA
4724 ! Prepare to freeze backbone
4731 ! Minimize the side-chains
4733 call geom_to_var(nvar,var)
4734 call minimize(etot,var,iretcode,nfun)
4735 call var_to_geom(nvar,var)
4738 ! Put the original weights back and calculate the full energy
4751 write(iout,*) 'Warning calling chainbuild'
4753 call etotal(energy_)
4757 end subroutine sc_minimize
4758 !-----------------------------------------------------------------------------
4759 subroutine minimize_sc1(etot,x,iretcode,nfun)
4761 use energy, only: func,gradient,fdum,etotal,enerprint
4764 ! implicit real*8 (a-h,o-z)
4765 ! include 'DIMENSIONS'
4766 integer,parameter :: liv=60
4767 integer :: iretcode,nfun
4768 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4769 ! include 'COMMON.IOUNITS'
4770 ! include 'COMMON.VAR'
4771 ! include 'COMMON.GEO'
4772 ! include 'COMMON.MINIM'
4773 !el integer :: icall
4774 !el common /srutu/ icall
4775 integer,dimension(liv) :: iv
4776 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
4777 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
4778 real(kind=8) :: energia(0:n_ene)
4779 !el real(kind=8) :: fdum
4780 ! external gradient,fdum !func,
4781 ! external func_restr1,grad_restr1
4782 logical :: not_done,change,reduce
4783 !el common /przechowalnia/ v
4785 integer :: nvar_restr,i,j
4787 real(kind=8) :: rdum(1),etot !,fdum
4789 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4790 if (.not. allocated(v)) allocate(v(1:lv))
4792 call deflt(2,iv,liv,lv,v)
4793 ! 12 means fresh start, dont call deflt
4795 ! max num of fun calls
4796 if (maxfun.eq.0) maxfun=500
4798 ! max num of iterations
4799 if (maxmin.eq.0) maxmin=1000
4803 ! selects output unit
4806 ! 1 means to print out result
4808 ! 1 means to print out summary stats
4810 ! 1 means to print initial x and d
4812 ! min val for v(radfac) default is 0.1
4814 ! max val for v(radfac) default is 4.0
4817 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
4818 ! the sumsl default is 0.1
4820 ! false conv if (act fnctn decrease) .lt. v(34)
4821 ! the sumsl default is 100*machep
4823 ! absolute convergence
4824 if (tolf.eq.0.0D0) tolf=1.0D-4
4826 ! relative convergence
4827 if (rtolf.eq.0.0D0) rtolf=1.0D-4
4829 ! controls initial step size
4831 ! large vals of d correspond to small components of step
4840 call x2xx(x,xx,nvar_restr)
4841 call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,&
4842 iv,liv,lv,v,idum,rdum,fdum)
4845 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
4847 !el---------------------
4848 ! write (iout,'(/a)') &
4849 ! "Cartesian coordinates of the reference structure after SUMSL"
4850 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
4851 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4853 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
4854 ! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
4855 ! (c(j,i+nres),j=1,3)
4857 ! call etotal(energia)
4858 ! call enerprint(energia)
4859 !el----------------------------
4865 end subroutine minimize_sc1
4866 !-----------------------------------------------------------------------------
4867 subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
4870 use energy, only: zerograd,esc,sc_grad
4871 ! implicit real*8 (a-h,o-z)
4872 ! include 'DIMENSIONS'
4873 ! include 'COMMON.DERIV'
4874 ! include 'COMMON.IOUNITS'
4875 ! include 'COMMON.GEO'
4876 ! include 'COMMON.FFIELD'
4877 ! include 'COMMON.INTERACT'
4878 ! include 'COMMON.TIME1'
4880 !el common /chuju/ jjj
4881 real(kind=8) :: energia(0:n_ene),evdw,escloc
4882 real(kind=8) :: e1,e2,f
4883 real(kind=8),external :: ufparm
4884 integer :: uiparm(1)
4885 real(kind=8) :: urparm(1)
4886 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
4891 ! Intercept NaNs in the coordinates, before calling etotal
4897 if (x_sum.ne.x_sum) then
4898 write(iout,*)" *** func_restr1 : Found NaN in coordinates"
4905 call var_to_geom_restr(n,x)
4907 write(iout,*) 'Warning calling chainbuild'
4909 !d write (iout,*) 'ETOTAL called from FUNC'
4912 f=wsc*evdw+wscloc*escloc
4913 !d call etotal(energia(0))
4914 !d f=wsc*energia(1)+wscloc*energia(12)
4915 !d print *,f,evdw,escloc,energia(0)
4917 ! Sum up the components of the Cartesian gradient.
4921 gradx(j,i,icg)=wsc*gvdwx(j,i)
4926 end subroutine func_restr1
4927 !-----------------------------------------------------------------------------
4929 subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
4931 use energy, only: cartder,zerograd,esc,sc_grad
4932 ! implicit real*8 (a-h,o-z)
4933 ! include 'DIMENSIONS'
4934 ! include 'COMMON.CHAIN'
4935 ! include 'COMMON.DERIV'
4936 ! include 'COMMON.VAR'
4937 ! include 'COMMON.INTERACT'
4938 ! include 'COMMON.FFIELD'
4939 ! include 'COMMON.IOUNITS'
4941 integer :: i,j,k,ind,n,nf,uiparm(1)
4942 real(kind=8) :: f,urparm(1)
4943 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
4944 integer :: ig,igall,ij
4945 real(kind=8) :: gphii,gthetai,galphai,gomegai
4946 real(kind=8),external :: ufparm
4949 if (nf-nfl+1) 20,30,40
4950 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
4951 ! write (iout,*) 'grad 20'
4954 30 call var_to_geom_restr(n,x)
4955 write(iout,*) 'Warning calling chainbuild'
4958 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
4962 ! Convert the Cartesian gradient into internal-coordinate gradient.
4968 IF (mask_phi(i+2).eq.1) THEN
4973 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
4974 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
4987 IF (mask_theta(i+2).eq.1) THEN
4993 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
4994 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
5004 if (itype(i,1).ne.10) then
5005 IF (mask_side(i).eq.1) THEN
5009 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
5018 if (itype(i,1).ne.10) then
5019 IF (mask_side(i).eq.1) THEN
5023 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
5031 ! Add the components corresponding to local energy terms.
5038 if (mask_phi(i).eq.1) then
5040 g(ig)=g(ig)+gloc(igall,icg)
5046 if (mask_theta(i).eq.1) then
5048 g(ig)=g(ig)+gloc(igall,icg)
5054 if (itype(i,1).ne.10) then
5056 if (mask_side(i).eq.1) then
5058 g(ig)=g(ig)+gloc(igall,icg)
5065 !d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
5068 end subroutine grad_restr1
5069 !-----------------------------------------------------------------------------
5070 subroutine egb1(evdw)
5072 ! This subroutine calculates the interaction energy of nonbonded side chains
5073 ! assuming the Gay-Berne potential of interaction.
5076 use energy, only: sc_grad
5077 ! use control, only:stopx
5078 ! implicit real*8 (a-h,o-z)
5079 ! include 'DIMENSIONS'
5080 ! include 'COMMON.GEO'
5081 ! include 'COMMON.VAR'
5082 ! include 'COMMON.LOCAL'
5083 ! include 'COMMON.CHAIN'
5084 ! include 'COMMON.DERIV'
5085 ! include 'COMMON.NAMES'
5086 ! include 'COMMON.INTERACT'
5087 ! include 'COMMON.IOUNITS'
5088 ! include 'COMMON.CALC'
5089 ! include 'COMMON.CONTROL'
5091 real(kind=8) :: evdw
5093 integer :: iint,ind,itypi,itypi1,itypj
5094 real(kind=8) :: xi,yi,zi,rrij,sig,sig0ij,rij_shift,fac,e1,e2,&
5096 !elwrite(iout,*) "check evdw"
5097 ! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
5100 ! if (icall.eq.0) lprn=.true.
5102 do i=iatsc_s,iatsc_e
5103 if ((itype(i,1).eq.ntyp1).or.(molnum(i).gt.1)) cycle
5104 itypi=iabs(itype(i,1))
5105 itypi1=iabs(itype(i+1,1))
5106 ! print *,"ebg1",i,itypi,itypi1
5110 dxi=dc_norm(1,nres+i)
5111 dyi=dc_norm(2,nres+i)
5112 dzi=dc_norm(3,nres+i)
5113 dsci_inv=dsc_inv(itypi)
5114 !elwrite(iout,*) itypi,itypi1,xi,yi,zi,dxi,dyi,dzi,dsci_inv
5116 ! Calculate SC interaction energy.
5118 do iint=1,nint_gr(i)
5119 do j=istart(i,iint),iend(i,iint)
5120 IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
5122 itypj=iabs(itype(j,1))
5123 if ((itype(j,1).eq.ntyp1).or.(molnum(j).gt.1)) cycle
5124 ! print *,"ebg1",j,itypj
5126 dscj_inv=dsc_inv(itypj)
5127 sig0ij=sigma(itypi,itypj)
5128 chi1=chi(itypi,itypj)
5129 chi2=chi(itypj,itypi)
5136 alf12=0.5D0*(alf1+alf2)
5137 ! For diagnostics only!!!
5150 dxj=dc_norm(1,nres+j)
5151 dyj=dc_norm(2,nres+j)
5152 dzj=dc_norm(3,nres+j)
5153 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
5155 ! Calculate angle-dependent terms of energy and contributions to their
5159 sig=sig0ij*dsqrt(sigsq)
5160 rij_shift=1.0D0/rij-sig+sig0ij
5161 ! I hate to put IF's in the loops, but here don't have another choice!!!!
5162 if (rij_shift.le.0.0D0) then
5164 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
5165 !d restyp(itypi),i,restyp(itypj),j, &
5166 !d rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
5170 !---------------------------------------------------------------
5171 rij_shift=1.0D0/rij_shift
5172 fac=rij_shift**expon
5173 e1=fac*fac*aa_aq(itypi,itypj)
5174 e2=fac*bb_aq(itypi,itypj)
5175 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
5176 eps2der=evdwij*eps3rt
5177 eps3der=evdwij*eps2rt
5178 evdwij=evdwij*eps2rt*eps3rt
5180 ! if (wliptran.gt.0.0) print *,"WARNING eps_aq used!"
5182 sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
5183 epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
5184 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
5185 !d restyp(itypi),i,restyp(itypj),j, &
5186 !d epsi,sigm,chi1,chi2,chip1,chip2, &
5187 !d eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
5188 !d om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
5192 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
5195 ! Calculate gradient components.
5196 e1=e1*eps1*eps2rt**2*eps3rt**2
5197 fac=-expon*(e1+evdwij)*rij_shift
5200 ! Calculate the radial part of the gradient
5204 ! Calculate angular part of the gradient.
5206 !elwrite(iout,*) evdw
5208 !elwrite(iout,*) "evdw=",evdw,j,iint,i
5210 !elwrite(iout,*) evdw
5212 !elwrite(iout,*) evdw
5214 !elwrite(iout,*) evdw
5216 !elwrite(iout,*) evdw,i
5218 !-----------------------------------------------------------------------------
5220 !-----------------------------------------------------------------------------
5221 subroutine sumsl(n,d,x,calcf,calcg,iv,liv,lv,v,uiparm,urparm,ufparm)
5223 ! *** minimize general unconstrained objective function using ***
5224 ! *** analytic gradient and hessian approx. from secant update ***
5230 integer :: uiparm(1)
5231 real(kind=8) :: d(n), x(n), v(lv), urparm(1)
5232 real(kind=8),external :: ufparm !funtion name as an argument
5234 ! dimension v(71 + n*(n+15)/2), uiparm(*), urparm(*)
5235 external :: calcf, calcg !subroutine name as an argument
5239 ! this routine interacts with subroutine sumit in an attempt
5240 ! to find an n-vector x* that minimizes the (unconstrained)
5241 ! objective function computed by calcf. (often the x* found is
5242 ! a local minimizer rather than a global one.)
5244 !-------------------------- parameter usage --------------------------
5246 ! n........ (input) the number of variables on which f depends, i.e.,
5247 ! the number of components in x.
5248 ! d........ (input/output) a scale vector such that d(i)*x(i),
5249 ! i = 1,2,...,n, are all in comparable units.
5250 ! d can strongly affect the behavior of sumsl.
5251 ! finding the best choice of d is generally a trial-
5252 ! and-error process. choosing d so that d(i)*x(i)
5253 ! has about the same value for all i often works well.
5254 ! the defaults provided by subroutine deflt (see i
5255 ! below) require the caller to supply d.
5256 ! x........ (input/output) before (initially) calling sumsl, the call-
5257 ! er should set x to an initial guess at x*. when
5258 ! sumsl returns, x contains the best point so far
5259 ! found, i.e., the one that gives the least value so
5260 ! far seen for f(x).
5261 ! calcf.... (input) a subroutine that, given x, computes f(x). calcf
5262 ! must be declared external in the calling program.
5264 ! call calcf(n, x, nf, f, uiparm, urparm, ufparm)
5265 ! when calcf is called, nf is the invocation
5266 ! count for calcf. nf is included for possible use
5267 ! with calcg. if x is out of bounds (e.g., if it
5268 ! would cause overflow in computing f(x)), then calcf
5269 ! should set nf to 0. this will cause a shorter step
5270 ! to be attempted. (if x is in bounds, then calcf
5271 ! should not change nf.) the other parameters are as
5272 ! described above and below. calcf should not change
5274 ! calcg.... (input) a subroutine that, given x, computes g(x), the gra-
5275 ! dient of f at x. calcg must be declared external in
5276 ! the calling program. it is invoked by
5277 ! call calcg(n, x, nf, g, uiparm, urparm, ufaprm)
5278 ! when calcg is called, nf is the invocation
5279 ! count for calcf at the time f(x) was evaluated. the
5280 ! x passed to calcg is usually the one passed to calcf
5281 ! on either its most recent invocation or the one
5282 ! prior to it. if calcf saves intermediate results
5283 ! for use by calcg, then it is possible to tell from
5284 ! nf whether they are valid for the current x (or
5285 ! which copy is valid if two copies are kept). if g
5286 ! cannot be computed at x, then calcg should set nf to
5287 ! 0. in this case, sumsl will return with iv(1) = 65.
5288 ! (if g can be computed at x, then calcg should not
5289 ! changed nf.) the other parameters to calcg are as
5290 ! described above and below. calcg should not change
5292 ! iv....... (input/output) an integer value array of length liv (see
5293 ! below) that helps control the sumsl algorithm and
5294 ! that is used to store various intermediate quanti-
5295 ! ties. of particular interest are the initialization/
5296 ! return code iv(1) and the entries in iv that control
5297 ! printing and limit the number of iterations and func-
5298 ! tion evaluations. see the section on iv input
5300 ! liv...... (input) length of iv array. must be at least 60. if li
5301 ! is too small, then sumsl returns with iv(1) = 15.
5302 ! when sumsl returns, the smallest allowed value of
5303 ! liv is stored in iv(lastiv) -- see the section on
5304 ! iv output values below. (this is intended for use
5305 ! with extensions of sumsl that handle constraints.)
5306 ! lv....... (input) length of v array. must be at least 71+n*(n+15)/2.
5307 ! (at least 77+n*(n+17)/2 for smsno, at least
5308 ! 78+n*(n+12) for humsl). if lv is too small, then
5309 ! sumsl returns with iv(1) = 16. when sumsl returns,
5310 ! the smallest allowed value of lv is stored in
5311 ! iv(lastv) -- see the section on iv output values
5313 ! v........ (input/output) a floating-point value array of length l
5314 ! (see below) that helps control the sumsl algorithm
5315 ! and that is used to store various intermediate
5316 ! quantities. of particular interest are the entries
5317 ! in v that limit the length of the first step
5318 ! attempted (lmax0) and specify convergence tolerances
5319 ! (afctol, lmaxs, rfctol, sctol, xctol, xftol).
5320 ! uiparm... (input) user integer parameter array passed without change
5321 ! to calcf and calcg.
5322 ! urparm... (input) user floating-point parameter array passed without
5323 ! change to calcf and calcg.
5324 ! ufparm... (input) user external subroutine or function passed without
5325 ! change to calcf and calcg.
5327 ! *** iv input values (from subroutine deflt) ***
5329 ! iv(1)... on input, iv(1) should have a value between 0 and 14......
5330 ! 0 and 12 mean this is a fresh start. 0 means that
5331 ! deflt(2, iv, liv, lv, v)
5332 ! is to be called to provide all default values to iv and
5333 ! v. 12 (the value that deflt assigns to iv(1)) means the
5334 ! caller has already called deflt and has possibly changed
5335 ! some iv and/or v entries to non-default values.
5336 ! 13 means deflt has been called and that sumsl (and
5337 ! sumit) should only do their storage allocation. that is,
5338 ! they should set the output components of iv that tell
5339 ! where various subarrays arrays of v begin, such as iv(g)
5340 ! (and, for humsl and humit only, iv(dtol)), and return.
5341 ! 14 means that a storage has been allocated (by a call
5342 ! with iv(1) = 13) and that the algorithm should be
5343 ! started. when called with iv(1) = 13, sumsl returns
5344 ! iv(1) = 14 unless liv or lv is too small (or n is not
5345 ! positive). default = 12.
5346 ! iv(inith).... iv(25) tells whether the hessian approximation h should
5347 ! be initialized. 1 (the default) means sumit should
5348 ! initialize h to the diagonal matrix whose i-th diagonal
5349 ! element is d(i)**2. 0 means the caller has supplied a
5350 ! cholesky factor l of the initial hessian approximation
5351 ! h = l*(l**t) in v, starting at v(iv(lmat)) = v(iv(42))
5352 ! (and stored compactly by rows). note that iv(lmat) may
5353 ! be initialized by calling sumsl with iv(1) = 13 (see
5354 ! the iv(1) discussion above). default = 1.
5355 ! iv(mxfcal)... iv(17) gives the maximum number of function evaluations
5356 ! (calls on calcf) allowed. if this number does not suf-
5357 ! fice, then sumsl returns with iv(1) = 9. default = 200.
5358 ! iv(mxiter)... iv(18) gives the maximum number of iterations allowed.
5359 ! it also indirectly limits the number of gradient evalua-
5360 ! tions (calls on calcg) to iv(mxiter) + 1. if iv(mxiter)
5361 ! iterations do not suffice, then sumsl returns with
5362 ! iv(1) = 10. default = 150.
5363 ! iv(outlev)... iv(19) controls the number and length of iteration sum-
5364 ! mary lines printed (by itsum). iv(outlev) = 0 means do
5365 ! not print any summary lines. otherwise, print a summary
5366 ! line after each abs(iv(outlev)) iterations. if iv(outlev)
5367 ! is positive, then summary lines of length 78 (plus carri-
5368 ! age control) are printed, including the following... the
5369 ! iteration and function evaluation counts, f = the current
5370 ! function value, relative difference in function values
5371 ! achieved by the latest step (i.e., reldf = (f0-v(f))/f01,
5372 ! where f01 is the maximum of abs(v(f)) and abs(v(f0)) and
5373 ! v(f0) is the function value from the previous itera-
5374 ! tion), the relative function reduction predicted for the
5375 ! step just taken (i.e., preldf = v(preduc) / f01, where
5376 ! v(preduc) is described below), the scaled relative change
5377 ! in x (see v(reldx) below), the step parameter for the
5378 ! step just taken (stppar = 0 means a full newton step,
5379 ! between 0 and 1 means a relaxed newton step, between 1
5380 ! and 2 means a double dogleg step, greater than 2 means
5381 ! a scaled down cauchy step -- see subroutine dbldog), the
5382 ! 2-norm of the scale vector d times the step just taken
5383 ! (see v(dstnrm) below), and npreldf, i.e.,
5384 ! v(nreduc)/f01, where v(nreduc) is described below -- if
5385 ! npreldf is positive, then it is the relative function
5386 ! reduction predicted for a newton step (one with
5387 ! stppar = 0). if npreldf is negative, then it is the
5388 ! negative of the relative function reduction predicted
5389 ! for a step computed with step bound v(lmaxs) for use in
5390 ! testing for singular convergence.
5391 ! if iv(outlev) is negative, then lines of length 50
5392 ! are printed, including only the first 6 items listed
5393 ! above (through reldx).
5395 ! iv(parprt)... iv(20) = 1 means print any nondefault v values on a
5396 ! fresh start or any changed v values on a restart.
5397 ! iv(parprt) = 0 means skip this printing. default = 1.
5398 ! iv(prunit)... iv(21) is the output unit number on which all printing
5399 ! is done. iv(prunit) = 0 means suppress all printing.
5400 ! default = standard output unit (unit 6 on most systems).
5401 ! iv(solprt)... iv(22) = 1 means print out the value of x returned (as
5402 ! well as the gradient and the scale vector d).
5403 ! iv(solprt) = 0 means skip this printing. default = 1.
5404 ! iv(statpr)... iv(23) = 1 means print summary statistics upon return-
5405 ! ing. these consist of the function value, the scaled
5406 ! relative change in x caused by the most recent step (see
5407 ! v(reldx) below), the number of function and gradient
5408 ! evaluations (calls on calcf and calcg), and the relative
5409 ! function reductions predicted for the last step taken and
5410 ! for a newton step (or perhaps a step bounded by v(lmaxs)
5411 ! -- see the descriptions of preldf and npreldf under
5412 ! iv(outlev) above).
5413 ! iv(statpr) = 0 means skip this printing.
5414 ! iv(statpr) = -1 means skip this printing as well as that
5415 ! of the one-line termination reason message. default = 1.
5416 ! iv(x0prt).... iv(24) = 1 means print the initial x and scale vector d
5417 ! (on a fresh start only). iv(x0prt) = 0 means skip this
5418 ! printing. default = 1.
5420 ! *** (selected) iv output values ***
5422 ! iv(1)........ on output, iv(1) is a return code....
5423 ! 3 = x-convergence. the scaled relative difference (see
5424 ! v(reldx)) between the current parameter vector x and
5425 ! a locally optimal parameter vector is very likely at
5427 ! 4 = relative function convergence. the relative differ-
5428 ! ence between the current function value and its lo-
5429 ! cally optimal value is very likely at most v(rfctol).
5430 ! 5 = both x- and relative function convergence (i.e., the
5431 ! conditions for iv(1) = 3 and iv(1) = 4 both hold).
5432 ! 6 = absolute function convergence. the current function
5433 ! value is at most v(afctol) in absolute value.
5434 ! 7 = singular convergence. the hessian near the current
5435 ! iterate appears to be singular or nearly so, and a
5436 ! step of length at most v(lmaxs) is unlikely to yield
5437 ! a relative function decrease of more than v(sctol).
5438 ! 8 = false convergence. the iterates appear to be converg-
5439 ! ing to a noncritical point. this may mean that the
5440 ! convergence tolerances (v(afctol), v(rfctol),
5441 ! v(xctol)) are too small for the accuracy to which
5442 ! the function and gradient are being computed, that
5443 ! there is an error in computing the gradient, or that
5444 ! the function or gradient is discontinuous near x.
5445 ! 9 = function evaluation limit reached without other con-
5446 ! vergence (see iv(mxfcal)).
5447 ! 10 = iteration limit reached without other convergence
5449 ! 11 = stopx returned .true. (external interrupt). see the
5450 ! usage notes below.
5451 ! 14 = storage has been allocated (after a call with
5453 ! 17 = restart attempted with n changed.
5454 ! 18 = d has a negative component and iv(dtype) .le. 0.
5455 ! 19...43 = v(iv(1)) is out of range.
5456 ! 63 = f(x) cannot be computed at the initial x.
5457 ! 64 = bad parameters passed to assess (which should not
5459 ! 65 = the gradient could not be computed at x (see calcg
5461 ! 67 = bad first parameter to deflt.
5462 ! 80 = iv(1) was out of range.
5463 ! 81 = n is not positive.
5464 ! iv(g)........ iv(28) is the starting subscript in v of the current
5465 ! gradient vector (the one corresponding to x).
5466 ! iv(lastiv)... iv(44) is the least acceptable value of liv. (it is
5467 ! only set if liv is at least 44.)
5468 ! iv(lastv).... iv(45) is the least acceptable value of lv. (it is
5469 ! only set if liv is large enough, at least iv(lastiv).)
5470 ! iv(nfcall)... iv(6) is the number of calls so far made on calcf (i.e.,
5471 ! function evaluations).
5472 ! iv(ngcall)... iv(30) is the number of gradient evaluations (calls on
5474 ! iv(niter).... iv(31) is the number of iterations performed.
5476 ! *** (selected) v input values (from subroutine deflt) ***
5478 ! v(bias)..... v(43) is the bias parameter used in subroutine dbldog --
5479 ! see that subroutine for details. default = 0.8.
5480 ! v(afctol)... v(31) is the absolute function convergence tolerance.
5481 ! if sumsl finds a point where the function value is less
5482 ! than v(afctol) in absolute value, and if sumsl does not
5483 ! return with iv(1) = 3, 4, or 5, then it returns with
5484 ! iv(1) = 6. this test can be turned off by setting
5485 ! v(afctol) to zero. default = max(10**-20, machep**2),
5486 ! where machep is the unit roundoff.
5487 ! v(dinit).... v(38), if nonnegative, is the value to which the scale
5488 ! vector d is initialized. default = -1.
5489 ! v(lmax0).... v(35) gives the maximum 2-norm allowed for d times the
5490 ! very first step that sumsl attempts. this parameter can
5491 ! markedly affect the performance of sumsl.
5492 ! v(lmaxs).... v(36) is used in testing for singular convergence -- if
5493 ! the function reduction predicted for a step of length
5494 ! bounded by v(lmaxs) is at most v(sctol) * abs(f0), where
5495 ! f0 is the function value at the start of the current
5496 ! iteration, and if sumsl does not return with iv(1) = 3,
5497 ! 4, 5, or 6, then it returns with iv(1) = 7. default = 1.
5498 ! v(rfctol)... v(32) is the relative function convergence tolerance.
5499 ! if the current model predicts a maximum possible function
5500 ! reduction (see v(nreduc)) of at most v(rfctol)*abs(f0)
5501 ! at the start of the current iteration, where f0 is the
5502 ! then current function value, and if the last step attempt-
5503 ! ed achieved no more than twice the predicted function
5504 ! decrease, then sumsl returns with iv(1) = 4 (or 5).
5505 ! default = max(10**-10, machep**(2/3)), where machep is
5506 ! the unit roundoff.
5507 ! v(sctol).... v(37) is the singular convergence tolerance -- see the
5508 ! description of v(lmaxs) above.
5509 ! v(tuner1)... v(26) helps decide when to check for false convergence.
5510 ! this is done if the actual function decrease from the
5511 ! current step is no more than v(tuner1) times its predict-
5512 ! ed value. default = 0.1.
5513 ! v(xctol).... v(33) is the x-convergence tolerance. if a newton step
5514 ! (see v(nreduc)) is tried that has v(reldx) .le. v(xctol)
5515 ! and if this step yields at most twice the predicted func-
5516 ! tion decrease, then sumsl returns with iv(1) = 3 (or 5).
5517 ! (see the description of v(reldx) below.)
5518 ! default = machep**0.5, where machep is the unit roundoff.
5519 ! v(xftol).... v(34) is the false convergence tolerance. if a step is
5520 ! tried that gives no more than v(tuner1) times the predict-
5521 ! ed function decrease and that has v(reldx) .le. v(xftol),
5522 ! and if sumsl does not return with iv(1) = 3, 4, 5, 6, or
5523 ! 7, then it returns with iv(1) = 8. (see the description
5524 ! of v(reldx) below.) default = 100*machep, where
5525 ! machep is the unit roundoff.
5526 ! v(*)........ deflt supplies to v a number of tuning constants, with
5527 ! which it should ordinarily be unnecessary to tinker. see
5528 ! section 17 of version 2.2 of the nl2sol usage summary
5529 ! (i.e., the appendix to ref. 1) for details on v(i),
5530 ! i = decfac, incfac, phmnfc, phmxfc, rdfcmn, rdfcmx,
5531 ! tuner2, tuner3, tuner4, tuner5.
5533 ! *** (selected) v output values ***
5535 ! v(dgnorm)... v(1) is the 2-norm of (diag(d)**-1)*g, where g is the
5536 ! most recently computed gradient.
5537 ! v(dstnrm)... v(2) is the 2-norm of diag(d)*step, where step is the
5539 ! v(f)........ v(10) is the current function value.
5540 ! v(f0)....... v(13) is the function value at the start of the current
5542 ! v(nreduc)... v(6), if positive, is the maximum function reduction
5543 ! possible according to the current model, i.e., the func-
5544 ! tion reduction predicted for a newton step (i.e.,
5545 ! step = -h**-1 * g, where g is the current gradient and
5546 ! h is the current hessian approximation).
5547 ! if v(nreduc) is negative, then it is the negative of
5548 ! the function reduction predicted for a step computed with
5549 ! a step bound of v(lmaxs) for use in testing for singular
5551 ! v(preduc)... v(7) is the function reduction predicted (by the current
5552 ! quadratic model) for the current step. this (divided by
5553 ! v(f0)) is used in testing for relative function
5555 ! v(reldx).... v(17) is the scaled relative change in x caused by the
5556 ! current step, computed as
5557 ! max(abs(d(i)*(x(i)-x0(i)), 1 .le. i .le. p) /
5558 ! max(d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p),
5559 ! where x = x0 + step.
5561 !------------------------------- notes -------------------------------
5563 ! *** algorithm notes ***
5565 ! this routine uses a hessian approximation computed from the
5566 ! bfgs update (see ref 3). only a cholesky factor of the hessian
5567 ! approximation is stored, and this is updated using ideas from
5568 ! ref. 4. steps are computed by the double dogleg scheme described
5569 ! in ref. 2. the steps are assessed as in ref. 1.
5571 ! *** usage notes ***
5573 ! after a return with iv(1) .le. 11, it is possible to restart,
5574 ! i.e., to change some of the iv and v input values described above
5575 ! and continue the algorithm from the point where it was interrupt-
5576 ! ed. iv(1) should not be changed, nor should any entries of i
5577 ! and v other than the input values (those supplied by deflt).
5578 ! those who do not wish to write a calcg which computes the
5579 ! gradient analytically should call smsno rather than sumsl.
5580 ! smsno uses finite differences to compute an approximate gradient.
5581 ! those who would prefer to provide f and g (the function and
5582 ! gradient) by reverse communication rather than by writing subrou-
5583 ! tines calcf and calcg may call on sumit directly. see the com-
5584 ! ments at the beginning of sumit.
5585 ! those who use sumsl interactively may wish to supply their
5586 ! own stopx function, which should return .true. if the break key
5587 ! has been pressed since stopx was last invoked. this makes it
5588 ! possible to externally interrupt sumsl (which will return with
5589 ! iv(1) = 11 if stopx returns .true.).
5590 ! storage for g is allocated at the end of v. thus the caller
5591 ! may make v longer than specified above and may allow calcg to use
5592 ! elements of g beyond the first n as scratch storage.
5594 ! *** portability notes ***
5596 ! the sumsl distribution tape contains both single- and double-
5597 ! precision versions of the sumsl source code, so it should be un-
5598 ! necessary to change precisions.
5599 ! only the functions imdcon and rmdcon contain machine-dependent
5600 ! constants. to change from one machine to another, it should
5601 ! suffice to change the (few) relevant lines in these functions.
5602 ! intrinsic functions are explicitly declared. on certain com-
5603 ! puters (e.g. univac), it may be necessary to comment out these
5604 ! declarations. so that this may be done automatically by a simple
5605 ! program, such declarations are preceded by a comment having c/+
5606 ! in columns 1-3 and blanks in columns 4-72 and are followed by
5607 ! a comment having c/ in columns 1 and 2 and blanks in columns 3-72.
5608 ! the sumsl source code is expressed in 1966 ansi standard
5609 ! fortran. it may be converted to fortran 77 by commenting out all
5610 ! lines that fall between a line having c/6 in columns 1-3 and a
5611 ! line having c/7 in columns 1-3 and by removing (i.e., replacing
5612 ! by a blank) the c in column 1 of the lines that follow the c/7
5613 ! line and precede a line having c/ in columns 1-2 and blanks in
5614 ! columns 3-72. these changes convert some data statements into
5615 ! parameter statements, convert some variables from real to
5616 ! character*4, and make the data statements that initialize these
5617 ! variables use character strings delimited by primes instead
5618 ! of hollerith constants. (such variables and data statements
5619 ! appear only in modules itsum and parck. parameter statements
5620 ! appear nearly everywhere.) these changes also add save state-
5621 ! ments for variables given machine-dependent constants by rmdcon.
5623 ! *** references ***
5625 ! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), algorithm 573 --
5626 ! an adaptive nonlinear least-squares algorithm, acm trans.
5627 ! math. software 7, pp. 369-383.
5629 ! 2. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
5630 ! mization algorithms which use function and gradient
5631 ! values, j. optim. theory applic. 28, pp. 453-482.
5633 ! 3. dennis, j.e., and more, j.j. (1977), quasi-newton methods, motiva-
5634 ! tion and theory, siam rev. 19, pp. 46-89.
5636 ! 4. goldfarb, d. (1976), factorized variable metric methods for uncon-
5637 ! strained optimization, math. comput. 30, pp. 796-811.
5641 ! coded by david m. gay (winter 1980). revised summer 1982.
5642 ! this subroutine was written in connection with research
5643 ! supported in part by the national science foundation under
5644 ! grants mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989,
5648 !---------------------------- declarations ---------------------------
5650 !el external deflt, sumit
5652 ! deflt... supplies default iv and v input components.
5653 ! sumit... reverse-communication routine that carries out sumsl algo-
5660 ! *** subscripts for iv ***
5662 !el integer nextv, nfcall, nfgcal, g, toobig, vneed
5665 ! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, toobig/2/, vneed/4/
5667 integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28,&
5671 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
5673 !elwrite(iout,*) "in sumsl"
5674 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
5676 if (iv1 .eq. 12 .or. iv1 .eq. 13) iv(vneed) = iv(vneed) + n
5677 if (iv1 .eq. 14) go to 10
5678 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
5680 if (iv1 .eq. 12) iv(1) = 13
5684 print *,"my new g1",g1
5685 !elwrite(iout,*) "in sumsl go to 10"
5688 !elwrite(iout,*) "in sumsl"
5689 20 call sumit(d, f, v(g1), iv, liv, lv, n, v, x)
5690 !elwrite(iout,*) "in sumsl, go to 20"
5692 !elwrite(iout,*) "in sumsl, go to 20, po sumit"
5693 !elwrite(iout,*) "in sumsl iv()", iv(1)-2
5694 if (iv(1) - 2) 30, 40, 50
5697 !elwrite(iout,*) "in sumsl iv",iv(nfcall)
5698 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
5699 !elwrite(iout,*) "in sumsl"
5700 if (nf .le. 0) iv(toobig) = 1
5703 !elwrite(iout,*) "in sumsl"
5704 40 call calcg(n, x, iv(nfgcal), v(g1), uiparm, urparm, ufparm)
5705 !elwrite(iout,*) "in sumsl"
5708 50 if (iv(1) .ne. 14) go to 999
5710 ! *** storage allocation
5713 iv(nextv) = iv(g) + n
5714 if (iv1 .ne. 13) go to 10
5715 !elwrite(iout,*) "in sumsl"
5718 ! *** last card of sumsl follows ***
5719 end subroutine sumsl
5720 !-----------------------------------------------------------------------------
5721 subroutine sumit(d,fx,g,iv,liv,lv,n,v,x)
5723 use control, only:stopx
5724 use MD_data, only: itime_mat
5726 ! *** carry out sumsl (unconstrained minimization) iterations, using
5727 ! *** double-dogleg/bfgs steps.
5729 ! *** parameter declarations ***
5734 real(kind=8) :: d(n), fx, g(n), v(lv), x(n)
5736 !-------------------------- parameter usage --------------------------
5738 ! d.... scale vector.
5739 ! fx... function value.
5740 ! g.... gradient vector.
5741 ! iv... integer value array.
5742 ! liv.. length of iv (at least 60).
5743 ! lv... length of v (at least 71 + n*(n+13)/2).
5744 ! n.... number of variables (components in x and g).
5745 ! v.... floating-point value array.
5746 ! x.... vector of parameters to be optimized.
5748 ! *** discussion ***
5750 ! parameters iv, n, v, and x are the same as the corresponding
5751 ! ones to sumsl (which see), except that v can be shorter (since
5752 ! the part of v that sumsl uses for storing g is not needed).
5753 ! moreover, compared with sumsl, iv(1) may have the two additional
5754 ! output values 1 and 2, which are explained below, as is the use
5755 ! of iv(toobig) and iv(nfgcal). the value iv(g), which is an
5756 ! output value from sumsl (and smsno), is not referenced by
5757 ! sumit or the subroutines it calls.
5758 ! fx and g need not have been initialized when sumit is called
5759 ! with iv(1) = 12, 13, or 14.
5761 ! iv(1) = 1 means the caller should set fx to f(x), the function value
5762 ! at x, and call sumit again, having changed none of the
5763 ! other parameters. an exception occurs if f(x) cannot be
5764 ! (e.g. if overflow would occur), which may happen because
5765 ! of an oversized step. in this case the caller should set
5766 ! iv(toobig) = iv(2) to 1, which will cause sumit to ig-
5767 ! nore fx and try a smaller step. the parameter nf that
5768 ! sumsl passes to calcf (for possible use by calcg) is a
5769 ! copy of iv(nfcall) = iv(6).
5770 ! iv(1) = 2 means the caller should set g to g(x), the gradient vector
5771 ! of f at x, and call sumit again, having changed none of
5772 ! the other parameters except possibly the scale vector d
5773 ! when iv(dtype) = 0. the parameter nf that sumsl passes
5774 ! to calcg is iv(nfgcal) = iv(7). if g(x) cannot be
5775 ! evaluated, then the caller may set iv(nfgcal) to 0, in
5776 ! which case sumit will return with iv(1) = 65.
5780 ! coded by david m. gay (december 1979). revised sept. 1982.
5781 ! this subroutine was written in connection with research supported
5782 ! in part by the national science foundation under grants
5783 ! mcs-7600324 and mcs-7906671.
5785 ! (see sumsl for references.)
5787 !+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
5789 ! *** local variables ***
5791 integer :: dg1, dummy, g01, i, k, l, lstgst, nwtst1, step1,&
5794 !el logical :: lstopx
5798 !el real(kind=8) :: half, negone, one, onep2, zero
5800 ! *** no intrinsic functions ***
5802 ! *** external functions and subroutines ***
5804 !el external assst, dbdog, deflt, dotprd, itsum, litvmu, livmul,
5805 !el 1 ltvmul, lupdat, lvmul, parck, reldst, stopx, vaxpy,
5806 !el 2 vcopy, vscopy, vvmulp, v2norm, wzbfgs
5808 !el real(kind=8) :: dotprd, reldst, v2norm
5810 ! assst.... assesses candidate step.
5811 ! dbdog.... computes double-dogleg (candidate) step.
5812 ! deflt.... supplies default iv and v input components.
5813 ! dotprd... returns inner product of two vectors.
5814 ! itsum.... prints iteration summary and info on initial and final x.
5815 ! litvmu... multiplies inverse transpose of lower triangle times vector.
5816 ! livmul... multiplies inverse of lower triangle times vector.
5817 ! ltvmul... multiplies transpose of lower triangle times vector.
5818 ! lupdt.... updates cholesky factor of hessian approximation.
5819 ! lvmul.... multiplies lower triangle times vector.
5820 ! parck.... checks validity of input iv and v values.
5821 ! reldst... computes v(reldx) = relative step size.
5822 ! stopx.... returns .true. if the break key has been pressed.
5823 ! vaxpy.... computes scalar times one vector plus another.
5824 ! vcopy.... copies one vector to another.
5825 ! vscopy... sets all elements of a vector to a scalar.
5826 ! vvmulp... multiplies vector by vector raised to power (componentwise).
5827 ! v2norm... returns the 2-norm of a vector.
5828 ! wzbfgs... computes w and z for lupdat corresponding to bfgs update.
5830 ! *** subscripts for iv and v ***
5833 !el integer cnvcod, dg, dgnorm, dinit, dstnrm, dst0, f, f0, fdif,
5834 !el 1 gthg, gtstep, g0, incfac, inith, irc, kagqt, lmat, lmax0,
5835 !el 2 lmaxs, mode, model, mxfcal, mxiter, nextv, nfcall, nfgcal,
5836 !el 3 ngcall, niter, nreduc, nwtstp, preduc, radfac, radinc,
5837 !el 4 radius, rad0, reldx, restor, step, stglim, stlstg, toobig,
5838 !el 5 tuner4, tuner5, vneed, xirc, x0
5840 ! *** iv subscript values ***
5843 ! data cnvcod/55/, dg/37/, g0/48/, inith/25/, irc/29/, kagqt/33/,
5844 ! 1 mode/35/, model/5/, mxfcal/17/, mxiter/18/, nfcall/6/,
5845 ! 2 nfgcal/7/, ngcall/30/, niter/31/, nwtstp/34/, radinc/8/,
5846 ! 3 restor/9/, step/40/, stglim/11/, stlstg/41/, toobig/2/,
5847 ! 4 vneed/4/, xirc/13/, x0/43/
5849 integer,parameter :: cnvcod=55, dg=37, g0=48, inith=25, irc=29, kagqt=33,&
5850 mode=35, model=5, mxfcal=17, mxiter=18, nfcall=6,&
5851 nfgcal=7, ngcall=30, niter=31, nwtstp=34, radinc=8,&
5852 restor=9, step=40, stglim=11, stlstg=41, toobig=2,&
5853 vneed=4, xirc=13, x0=43
5856 ! *** v subscript values ***
5860 ! data dgnorm/1/, dinit/38/, dstnrm/2/, dst0/3/, f/10/, f0/13/,
5861 ! 1 fdif/11/, gthg/44/, gtstep/4/, incfac/23/, lmat/42/,
5862 ! 2 lmax0/35/, lmaxs/36/, nextv/47/, nreduc/6/, preduc/7/,
5863 ! 3 radfac/16/, radius/8/, rad0/9/, reldx/17/, tuner4/29/,
5866 integer,parameter :: afctol=31
5867 integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dst0=3, f=10, f0=13,&
5868 fdif=11, gthg=44, gtstep=4, incfac=23, lmat=42,&
5869 lmax0=35, lmaxs=36, nextv=47, nreduc=6, preduc=7,&
5870 radfac=16, radius=8, rad0=9, reldx=17, tuner4=29,&
5875 ! data half/0.5d+0/, negone/-1.d+0/, one/1.d+0/, onep2/1.2d+0/,
5878 real(kind=8),parameter :: half=0.5d+0, negone=-1.d+0, one=1.d+0,&
5879 onep2=1.2d+0,zero=0.d+0
5882 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
5884 ! Following SAVE statement inserted.
5887 if (i .eq. 1) go to 50
5888 if (i .eq. 2) go to 60
5890 ! *** check validity of iv and v input values ***
5892 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
5893 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
5894 iv(vneed) = iv(vneed) + n*(n+13)/2
5895 call parck(2, d, iv, liv, lv, n, v)
5897 if (i .gt. 12) go to 999
5898 go to (180, 180, 180, 180, 180, 180, 120, 90, 120, 10, 10, 20), i
5900 ! *** storage allocation ***
5903 iv(x0) = l + n*(n+1)/2
5904 iv(step) = iv(x0) + n
5905 iv(stlstg) = iv(step) + n
5906 iv(g0) = iv(stlstg) + n
5907 iv(nwtstp) = iv(g0) + n
5908 iv(dg) = iv(nwtstp) + n
5909 iv(nextv) = iv(dg) + n
5910 if (iv(1) .ne. 13) go to 20
5914 ! *** initialization ***
5927 if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
5928 if (iv(inith) .ne. 1) go to 40
5930 ! *** set the initial hessian approximation to diag(d)**-2 ***
5933 call vscopy(n*(n+1)/2, v(l), zero)
5938 if (t .le. zero) t = one
5942 ! *** compute initial function value ***
5948 if (iv(mode) .ge. 0) go to 180
5950 if (iv(toobig) .eq. 0) go to 999
5954 ! *** make sure gradient could be computed ***
5956 60 if (iv(nfgcal) .ne. 0) go to 70
5961 call vvmulp(n, v(dg1), g, d, -1)
5962 v(dgnorm) = v2norm(n, v(dg1))
5964 ! *** test norm of gradient ***
5966 if (v(dgnorm) .gt. v(afctol)) go to 75
5968 iv(cnvcod) = iv(irc) - 4
5970 75 if (iv(cnvcod) .ne. 0) go to 290
5971 if (iv(mode) .eq. 0) go to 250
5973 ! *** allow first step to have scaled 2-norm at most v(lmax0) ***
5975 v(radius) = v(lmax0)
5980 !----------------------------- main loop -----------------------------
5983 ! *** print iteration summary, check iteration limit ***
5985 80 call itsum(d, g, iv, liv, lv, n, v, x)
5989 if (k .lt. iv(mxiter)) go to 100
5993 ! *** update radius ***
5995 100 iv(niter) = k + 1
5996 if(k.gt.0)v(radius) = v(radfac) * v(dstnrm)
5998 ! *** initialize for start of next iteration ***
6006 ! *** copy x to x0, g to g0 ***
6008 call vcopy(n, v(x01), x)
6009 call vcopy(n, v(g01), g)
6011 ! *** check stopx and function evaluation limit ***
6015 !el lstopx = stopx(dummy)
6016 !elwrite(iout,*) "lstopx",lstopx,dummy
6017 110 if (.not. stopx(dummy)) go to 130
6019 ! write (iout,*) "iv(1)=11 !!!!"
6022 ! *** come here when restarting after func. eval. limit or stopx.
6024 120 if (v(f) .ge. v(f0)) go to 130
6029 130 if (iv(nfcall) .lt. iv(mxfcal)) go to 150
6031 140 if (v(f) .ge. v(f0)) go to 300
6033 ! *** in case of stopx or function evaluation limit with
6034 ! *** improved v(f), evaluate the gradient at x.
6039 !. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
6041 150 step1 = iv(step)
6044 if (iv(kagqt) .ge. 0) go to 160
6046 call livmul(n, v(nwtst1), v(l), g)
6047 v(nreduc) = half * dotprd(n, v(nwtst1), v(nwtst1))
6048 call litvmu(n, v(nwtst1), v(l), v(nwtst1))
6049 call vvmulp(n, v(step1), v(nwtst1), d, 1)
6050 v(dst0) = v2norm(n, v(step1))
6051 call vvmulp(n, v(dg1), v(dg1), d, -1)
6052 call ltvmul(n, v(step1), v(l), v(dg1))
6053 v(gthg) = v2norm(n, v(step1))
6055 160 call dbdog(v(dg1), lv, n, v(nwtst1), v(step1), v)
6056 if (iv(irc) .eq. 6) go to 180
6058 ! *** check whether evaluating f(x0 + step) looks worthwhile ***
6060 if (v(dstnrm) .le. zero) go to 180
6061 if (iv(irc) .ne. 5) go to 170
6062 if (v(radfac) .le. one) go to 170
6063 if (v(preduc) .le. onep2 * v(fdif)) go to 180
6065 ! *** compute f(x0 + step) ***
6069 call vaxpy(n, x, one, v(step1), v(x01))
6070 iv(nfcall) = iv(nfcall) + 1
6075 !. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
6078 v(reldx) = reldst(n, d, x, v(x01))
6079 call assst(iv, liv, lv, v)
6082 if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
6083 if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
6084 if (iv(restor) .ne. 3) go to 190
6085 call vcopy(n, v(step1), v(lstgst))
6086 call vaxpy(n, x, one, v(step1), v(x01))
6087 v(reldx) = reldst(n, d, x, v(x01))
6090 go to (200,230,230,230,200,210,220,220,220,220,220,220,280,250), k
6092 ! *** recompute step with changed radius ***
6094 200 v(radius) = v(radfac) * v(dstnrm)
6097 ! *** compute step of length v(lmaxs) for singular convergence test.
6099 210 v(radius) = v(lmaxs)
6102 ! *** convergence or false convergence ***
6104 220 iv(cnvcod) = k - 4
6105 if (v(f) .ge. v(f0)) go to 290
6106 if (iv(xirc) .eq. 14) go to 290
6109 !. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
6111 230 if (iv(irc) .ne. 3) go to 240
6115 ! *** set temp1 = hessian * step for use in gradient tests ***
6118 call ltvmul(n, v(temp1), v(l), v(step1))
6119 call lvmul(n, v(temp1), v(l), v(temp1))
6121 ! *** compute gradient ***
6123 240 iv(ngcall) = iv(ngcall) + 1
6127 ! *** initializations -- g0 = g - g0, etc. ***
6130 call vaxpy(n, v(g01), negone, v(g01), g)
6133 if (iv(irc) .ne. 3) go to 270
6135 ! *** set v(radfac) by gradient tests ***
6137 ! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
6139 call vaxpy(n, v(temp1), negone, v(g01), v(temp1))
6140 call vvmulp(n, v(temp1), v(temp1), d, -1)
6142 ! *** do gradient tests ***
6144 if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) &
6146 if (dotprd(n, g, v(step1)) &
6147 .ge. v(gtstep) * v(tuner5)) go to 270
6148 260 v(radfac) = v(incfac)
6150 ! *** update h, loop ***
6155 call wzbfgs(v(l), n, v(step1), v(w), v(g01), v(z))
6157 ! ** use the n-vectors starting at v(step1) and v(g01) for scratch..
6158 call lupdat(v(temp1), v(step1), v(l), v(g01), v(l), n, v(w), v(z))
6162 !. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
6164 ! *** bad parameters to assess ***
6169 ! *** print summary of final iteration and other requested items ***
6171 290 iv(1) = iv(cnvcod)
6173 300 call itsum(d, g, iv, liv, lv, n, v, x)
6177 ! *** last line of sumit follows ***
6178 end subroutine sumit
6179 !-----------------------------------------------------------------------------
6180 subroutine dbdog(dig,lv,n,nwtstp,step,v)
6182 ! *** compute double dogleg step ***
6184 ! *** parameter declarations ***
6188 real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv)
6192 ! this subroutine computes a candidate step (for use in an uncon-
6193 ! strained minimization code) by the double dogleg algorithm of
6194 ! dennis and mei (ref. 1), which is a variation on powell*s dogleg
6195 ! scheme (ref. 2, p. 95).
6197 !-------------------------- parameter usage --------------------------
6199 ! dig (input) diag(d)**-2 * g -- see algorithm notes.
6200 ! g (input) the current gradient vector.
6201 ! lv (input) length of v.
6202 ! n (input) number of components in dig, g, nwtstp, and step.
6203 ! nwtstp (input) negative newton step -- see algorithm notes.
6204 ! step (output) the computed step.
6205 ! v (i/o) values array, the following components of which are
6207 ! v(bias) (input) bias for relaxed newton step, which is v(bias) of
6208 ! the way from the full newton to the fully relaxed newton
6209 ! step. recommended value = 0.8 .
6210 ! v(dgnorm) (input) 2-norm of diag(d)**-1 * g -- see algorithm notes.
6211 ! v(dstnrm) (output) 2-norm of diag(d) * step, which is v(radius)
6212 ! unless v(stppar) = 0 -- see algorithm notes.
6213 ! v(dst0) (input) 2-norm of diag(d) * nwtstp -- see algorithm notes.
6214 ! v(grdfac) (output) the coefficient of dig in the step returned --
6215 ! step(i) = v(grdfac)*dig(i) + v(nwtfac)*nwtstp(i).
6216 ! v(gthg) (input) square-root of (dig**t) * (hessian) * dig -- see
6218 ! v(gtstep) (output) inner product between g and step.
6219 ! v(nreduc) (output) function reduction predicted for the full newton
6221 ! v(nwtfac) (output) the coefficient of nwtstp in the step returned --
6222 ! see v(grdfac) above.
6223 ! v(preduc) (output) function reduction predicted for the step returned.
6224 ! v(radius) (input) the trust region radius. d times the step returned
6225 ! has 2-norm v(radius) unless v(stppar) = 0.
6226 ! v(stppar) (output) code telling how step was computed... 0 means a
6227 ! full newton step. between 0 and 1 means v(stppar) of the
6228 ! way from the newton to the relaxed newton step. between
6229 ! 1 and 2 means a true double dogleg step, v(stppar) - 1 of
6230 ! the way from the relaxed newton to the cauchy step.
6231 ! greater than 2 means 1 / (v(stppar) - 1) times the cauchy
6234 !------------------------------- notes -------------------------------
6236 ! *** algorithm notes ***
6238 ! let g and h be the current gradient and hessian approxima-
6239 ! tion respectively and let d be the current scale vector. this
6240 ! routine assumes dig = diag(d)**-2 * g and nwtstp = h**-1 * g.
6241 ! the step computed is the same one would get by replacing g and h
6242 ! by diag(d)**-1 * g and diag(d)**-1 * h * diag(d)**-1,
6243 ! computing step, and translating step back to the original
6244 ! variables, i.e., premultiplying it by diag(d)**-1.
6246 ! *** references ***
6248 ! 1. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
6249 ! mization algorithms which use function and gradient
6250 ! values, j. optim. theory applic. 28, pp. 453-482.
6251 ! 2. powell, m.j.d. (1970), a hybrid method for non-linear equations,
6252 ! in numerical methods for non-linear equations, edited by
6253 ! p. rabinowitz, gordon and breach, london.
6257 ! coded by david m. gay.
6258 ! this subroutine was written in connection with research supported
6259 ! by the national science foundation under grants mcs-7600324 and
6262 !------------------------ external quantities ------------------------
6264 ! *** functions and subroutines called ***
6266 !el external dotprd, v2norm
6267 !el real(kind=8) :: dotprd, v2norm
6269 ! dotprd... returns inner product of two vectors.
6270 ! v2norm... returns 2-norm of a vector.
6272 ! *** intrinsic functions ***
6274 !el real(kind=8) :: dsqrt
6276 !-------------------------- local variables --------------------------
6279 real(kind=8) :: cfact, cnorm, ctrnwt, ghinvg, femnsq, gnorm,&
6280 nwtnrm, relax, rlambd, t, t1, t2
6281 !el real(kind=8) :: half, one, two, zero
6283 ! *** v subscripts ***
6285 !el integer bias, dgnorm, dstnrm, dst0, grdfac, gthg, gtstep,
6286 !el 1 nreduc, nwtfac, preduc, radius, stppar
6288 ! *** data initializations ***
6291 ! data half/0.5d+0/, one/1.d+0/, two/2.d+0/, zero/0.d+0/
6293 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, two=2.d+0, zero=0.d+0
6297 ! data bias/43/, dgnorm/1/, dstnrm/2/, dst0/3/, grdfac/45/,
6298 ! 1 gthg/44/, gtstep/4/, nreduc/6/, nwtfac/46/, preduc/7/,
6299 ! 2 radius/8/, stppar/5/
6301 integer,parameter :: bias=43, dgnorm=1, dstnrm=2, dst0=3, grdfac=45,&
6302 gthg=44, gtstep=4, nreduc=6, nwtfac=46, preduc=7,&
6306 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6310 if (nwtnrm .gt. zero) rlambd = v(radius) / nwtnrm
6312 ghinvg = two * v(nreduc)
6315 if (rlambd .lt. one) go to 30
6317 ! *** the newton step is inside the trust region ***
6322 v(preduc) = v(nreduc)
6325 20 step(i) = -nwtstp(i)
6328 30 v(dstnrm) = v(radius)
6329 cfact = (gnorm / v(gthg))**2
6330 ! *** cauchy step = -cfact * g.
6331 cnorm = gnorm * cfact
6332 relax = one - v(bias) * (one - gnorm*cnorm/ghinvg)
6333 if (rlambd .lt. relax) go to 50
6335 ! *** step is between relaxed newton and full newton steps ***
6337 v(stppar) = one - (rlambd - relax) / (one - relax)
6339 v(gtstep) = t * ghinvg
6340 v(preduc) = rlambd * (one - half*rlambd) * ghinvg
6343 40 step(i) = t * nwtstp(i)
6346 50 if (cnorm .lt. v(radius)) go to 70
6348 ! *** the cauchy step lies outside the trust region --
6349 ! *** step = scaled cauchy step ***
6351 t = -v(radius) / gnorm
6353 v(stppar) = one + cnorm / v(radius)
6354 v(gtstep) = -v(radius) * gnorm
6355 v(preduc) = v(radius)*(gnorm - half*v(radius)*(v(gthg)/gnorm)**2)
6357 60 step(i) = t * dig(i)
6360 ! *** compute dogleg step between cauchy and relaxed newton ***
6361 ! *** femur = relaxed newton step minus cauchy step ***
6363 70 ctrnwt = cfact * relax * ghinvg / gnorm
6364 ! *** ctrnwt = inner prod. of cauchy and relaxed newton steps,
6365 ! *** scaled by gnorm**-1.
6366 t1 = ctrnwt - gnorm*cfact**2
6367 ! *** t1 = inner prod. of femur and cauchy step, scaled by
6369 t2 = v(radius)*(v(radius)/gnorm) - gnorm*cfact**2
6371 femnsq = (t/gnorm)*t - ctrnwt - t1
6372 ! *** femnsq = square of 2-norm of femur, scaled by gnorm**-1.
6373 t = t2 / (t1 + dsqrt(t1**2 + femnsq*t2))
6374 ! *** dogleg step = cauchy step + t * femur.
6375 t1 = (t - one) * cfact
6380 v(gtstep) = t1*gnorm**2 + t2*ghinvg
6381 v(preduc) = -t1*gnorm * ((t2 + one)*gnorm) &
6382 - t2 * (one + half*t2)*ghinvg &
6383 - half * (v(gthg)*t1)**2
6385 80 step(i) = t1*dig(i) + t2*nwtstp(i)
6388 ! *** last line of dbdog follows ***
6389 end subroutine dbdog
6390 !-----------------------------------------------------------------------------
6391 subroutine ltvmul(n,x,l,y)
6393 ! *** compute x = (l**t)*y, where l is an n x n lower
6394 ! *** triangular matrix stored compactly by rows. x and y may
6395 ! *** occupy the same storage. ***
6398 !al real(kind=8) :: x(n), l(1), y(n)
6399 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
6400 ! dimension l(n*(n+1)/2)
6401 integer :: i, ij, i0, j
6402 real(kind=8) :: yi !el, zero
6406 real(kind=8),parameter :: zero=0.d+0
6415 x(j) = x(j) + yi*l(ij)
6420 ! *** last card of ltvmul follows ***
6421 end subroutine ltvmul
6422 !-----------------------------------------------------------------------------
6423 subroutine lupdat(beta,gamma,l,lambda,lplus,n,w,z)
6425 ! *** compute lplus = secant update of l ***
6427 ! *** parameter declarations ***
6430 !al double precision beta(n), gamma(n), l(1), lambda(n), lplus(1),
6431 real(kind=8) :: beta(n), gamma(n), l(n*(n+1)/2), lambda(n), &
6432 lplus(n*(n+1)/2),w(n), z(n)
6433 ! dimension l(n*(n+1)/2), lplus(n*(n+1)/2)
6435 !-------------------------- parameter usage --------------------------
6437 ! beta = scratch vector.
6438 ! gamma = scratch vector.
6439 ! l (input) lower triangular matrix, stored rowwise.
6440 ! lambda = scratch vector.
6441 ! lplus (output) lower triangular matrix, stored rowwise, which may
6442 ! occupy the same storage as l.
6443 ! n (input) length of vector parameters and order of matrices.
6444 ! w (input, destroyed on output) right singular vector of rank 1
6446 ! z (input, destroyed on output) left singular vector of rank 1
6449 !------------------------------- notes -------------------------------
6451 ! *** application and usage restrictions ***
6453 ! this routine updates the cholesky factor l of a symmetric
6454 ! positive definite matrix to which a secant update is being
6455 ! applied -- it computes a cholesky factor lplus of
6456 ! l * (i + z*w**t) * (i + w*z**t) * l**t. it is assumed that w
6457 ! and z have been chosen so that the updated matrix is strictly
6458 ! positive definite.
6460 ! *** algorithm notes ***
6462 ! this code uses recurrence 3 of ref. 1 (with d(j) = 1 for all j)
6463 ! to compute lplus of the form l * (i + z*w**t) * q, where q
6464 ! is an orthogonal matrix that makes the result lower triangular.
6465 ! lplus may have some negative diagonal elements.
6467 ! *** references ***
6469 ! 1. goldfarb, d. (1976), factorized variable metric methods for uncon-
6470 ! strained optimization, math. comput. 30, pp. 796-811.
6474 ! coded by david m. gay (fall 1979).
6475 ! this subroutine was written in connection with research supported
6476 ! by the national science foundation under grants mcs-7600324 and
6479 !------------------------ external quantities ------------------------
6481 ! *** intrinsic functions ***
6483 !el real(kind=8) :: dsqrt
6485 !-------------------------- local variables --------------------------
6487 integer :: i, ij, j, jj, jp1, k, nm1, np1
6488 real(kind=8) :: a, b, bj, eta, gj, lj, lij, ljj, nu, s, theta,&
6490 !el real(kind=8) :: one, zero
6492 ! *** data initializations ***
6495 ! data one/1.d+0/, zero/0.d+0/
6497 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
6500 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6504 if (n .le. 1) go to 30
6507 ! *** temporarily store s(j) = sum over k = j+1 to n of w(k)**2 in
6517 ! *** compute lambda, gamma, and beta by goldfarb*s recurrence 3.
6521 a = nu*z(j) - eta*wj
6524 lj = dsqrt(theta**2 + a*s)
6525 if (theta .gt. zero) lj = -lj
6528 gamma(j) = b * nu / lj
6529 beta(j) = (a - b*eta) / lj
6531 eta = -(eta + (a**2)/(theta - lj)) / lj
6533 30 lambda(n) = one + (nu*z(n) - eta*w(n))*w(n)
6535 ! *** update l, gradually overwriting w and z with l*w and l*z.
6538 jj = n * (n + 1) / 2
6543 lplus(jj) = lj * ljj
6548 if (k .eq. 1) go to 50
6555 lplus(ij) = lj*lij + bj*w(i) + gj*z(i)
6556 w(i) = w(i) + lij*wj
6557 z(i) = z(i) + lij*zj
6564 ! *** last card of lupdat follows ***
6565 end subroutine lupdat
6566 !-----------------------------------------------------------------------------
6567 subroutine lvmul(n,x,l,y)
6569 ! *** compute x = l*y, where l is an n x n lower triangular
6570 ! *** matrix stored compactly by rows. x and y may occupy the same
6574 !al double precision x(n), l(1), y(n)
6575 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
6576 ! dimension l(n*(n+1)/2)
6577 integer :: i, ii, ij, i0, j, np1
6578 real(kind=8) :: t !el, zero
6582 real(kind=8),parameter :: zero=0.d+0
6598 ! *** last card of lvmul follows ***
6599 end subroutine lvmul
6600 !-----------------------------------------------------------------------------
6601 subroutine vvmulp(n,x,y,z,k)
6603 ! *** set x(i) = y(i) * z(i)**k, 1 .le. i .le. n (for k = 1 or -1) ***
6606 real(kind=8) :: x(n), y(n), z(n)
6609 if (k .ge. 0) go to 20
6611 10 x(i) = y(i) / z(i)
6615 30 x(i) = y(i) * z(i)
6617 ! *** last card of vvmulp follows ***
6618 end subroutine vvmulp
6619 !-----------------------------------------------------------------------------
6620 subroutine wzbfgs(l,n,s,w,y,z)
6622 ! *** compute y and z for lupdat corresponding to bfgs update.
6625 !al double precision l(1), s(n), w(n), y(n), z(n)
6626 real(kind=8) :: l(n*(n+1)/2), s(n), w(n), y(n), z(n)
6627 ! dimension l(n*(n+1)/2)
6629 !-------------------------- parameter usage --------------------------
6631 ! l (i/o) cholesky factor of hessian, a lower triang. matrix stored
6632 ! compactly by rows.
6633 ! n (input) order of l and length of s, w, y, z.
6634 ! s (input) the step just taken.
6635 ! w (output) right singular vector of rank 1 correction to l.
6636 ! y (input) change in gradients corresponding to s.
6637 ! z (output) left singular vector of rank 1 correction to l.
6639 !------------------------------- notes -------------------------------
6641 ! *** algorithm notes ***
6643 ! when s is computed in certain ways, e.g. by gqtstp or
6644 ! dbldog, it is possible to save n**2/2 operations since (l**t)*s
6645 ! or l*(l**t)*s is then known.
6646 ! if the bfgs update to l*(l**t) would reduce its determinant to
6647 ! less than eps times its old value, then this routine in effect
6648 ! replaces y by theta*y + (1 - theta)*l*(l**t)*s, where theta
6649 ! (between 0 and 1) is chosen to make the reduction factor = eps.
6653 ! coded by david m. gay (fall 1979).
6654 ! this subroutine was written in connection with research supported
6655 ! by the national science foundation under grants mcs-7600324 and
6658 !------------------------ external quantities ------------------------
6660 ! *** functions and subroutines called ***
6662 !el external dotprd, livmul, ltvmul
6663 !el real(kind=8) :: dotprd
6664 ! dotprd returns inner product of two vectors.
6665 ! livmul multiplies l**-1 times a vector.
6666 ! ltvmul multiplies l**t times a vector.
6668 ! *** intrinsic functions ***
6670 !el real(kind=8) :: dsqrt
6672 !-------------------------- local variables --------------------------
6675 real(kind=8) :: cs, cy, epsrt, shs, ys, theta !el, eps, one
6677 ! *** data initializations ***
6680 ! data eps/0.1d+0/, one/1.d+0/
6682 real(kind=8),parameter :: eps=0.1d+0, one=1.d+0
6685 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6687 call ltvmul(n, w, l, s)
6688 shs = dotprd(n, w, w)
6689 ys = dotprd(n, y, s)
6690 if (ys .ge. eps*shs) go to 10
6691 theta = (one - eps) * shs / (shs - ys)
6693 cy = theta / (shs * epsrt)
6694 cs = (one + (theta-one)/epsrt) / shs
6696 10 cy = one / (dsqrt(ys) * dsqrt(shs))
6698 20 call livmul(n, z, l, y)
6700 30 z(i) = cy * z(i) - cs * w(i)
6703 ! *** last card of wzbfgs follows ***
6704 end subroutine wzbfgs
6705 !-----------------------------------------------------------------------------
6709 ! ###################################################
6710 ! ## COPYRIGHT (C) 1999 by Jay William Ponder ##
6711 ! ## All Rights Reserved ##
6712 ! ###################################################
6714 ! ##############################################################
6716 ! ## subroutine lbfgs -- limited memory BFGS optimization ##
6718 ! ##############################################################
6721 ! "lbfgs" is a limited memory BFGS quasi-newton nonlinear
6722 ! optimization routine
6724 ! literature references:
6726 ! J. Nocedal, "Updating Quasi-Newton Matrices with Limited
6727 ! Storage", Mathematics of Computation, 35, 773-782 (1980)
6729 ! D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method
6730 ! for Large Scale Optimization", Mathematical Programming,
6731 ! 45, 503-528 (1989)
6733 ! J. Nocedal and S. J. Wright, "Numerical Optimization",
6734 ! Springer-Verlag, New York, 1999, Section 9.1
6736 ! variables and parameters:
6738 ! nvar number of parameters in the objective function
6739 ! x0 contains starting point upon input, upon return
6740 ! contains the best point found
6741 ! minimum during optimization contains best current function
6742 ! value; returns final best function value
6743 ! grdmin normal exit if rms gradient gets below this value
6744 ! ncalls total number of function/gradient evaluations
6746 ! required external routines:
6748 ! fgvalue function to evaluate function and gradient values
6749 ! optsave subroutine to write out info about current status
6752 subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave)
6761 use MD_data, only: itime_mat
6763 use energy_data,only:statusbf,niter,ncalls
6770 integer niter,ncalls
6771 character*9 statusbf
6774 real*8 f,f_old,fgvalue
6775 real*8 f_move,x_move
6777 real*8 minimum,grdmin
6778 real*8 angle,rms,beta
6781 real*8, allocatable :: rho(:)
6782 real*8, allocatable :: alpha(:)
6783 real*8, allocatable :: x_old(:)
6784 real*8, allocatable :: g(:)
6785 real*8, allocatable :: g_old(:)
6786 real*8, allocatable :: p(:)
6787 real*8, allocatable :: q(:)
6788 real*8, allocatable :: r(:)
6789 real*8, allocatable :: h0(:)
6790 real*8, allocatable :: s(:,:)
6791 real*8, allocatable :: y(:,:)
6794 character*20 keyword
6795 character*240 record
6796 character*240 string
6797 external fgvalue,optsave
6798 ! common /lbfgstat/ status,niter,ncalls
6801 !c initialize some values to be used below
6804 rms = sqrt(dble(nvar))
6805 if (coordtype .eq. 'CARTESIAN') then
6806 rms = rms / sqrt(3.0d0)
6807 else if (coordtype .eq. 'RIGIDBODY') then
6808 rms = rms / sqrt(6.0d0)
6815 !c perform dynamic allocation of some global arrays
6817 if (.not. allocated(scale)) allocate (scale(nvar))
6819 !c set default values for variable scale factors
6821 if (.not. set_scale) then
6823 if (scale(i) .eq. 0.0d0) scale(i) = 1.0d0
6827 !c set default parameters for the optimization
6830 if (fctmin .eq. 0.0d0) fctmin = -100000000.0d0
6831 if (maxiter .eq. 0) maxiter = 1000000
6832 if (nextiter .eq. 0) nextiter = 1
6833 if (jprint .lt. 0) jprint = 1
6834 if (iwrite .lt. 0) iwrite = 1
6835 write(iout,*) "maxiter",maxiter
6837 !c set default parameters for the line search
6839 if (stpmax .eq. 0.0d0) stpmax = 5.0d0
6846 !c search the keywords for optimization parameters
6852 call gettext (record,keyword,next)
6853 call upcase (keyword)
6854 string = record(next:240)
6855 if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then
6856 read (string,*,err=10,end=10) msav
6857 msav = max(0,min(msav,nvar))
6858 else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then
6860 else if (keyword(1:7) .eq. 'FCTMIN ') then
6861 read (string,*,err=10,end=10) fctmin
6862 else if (keyword(1:8) .eq. 'MAXITER ') then
6863 read (string,*,err=10,end=10) maxiter
6864 else if (keyword(1:8) .eq. 'STEPMAX ') then
6865 read (string,*,err=10,end=10) stpmax
6866 else if (keyword(1:8) .eq. 'STEPMIN ') then
6867 read (string,*,err=10,end=10) stpmin
6868 else if (keyword(1:6) .eq. 'CAPPA ') then
6869 read (string,*,err=10,end=10) cappa
6870 else if (keyword(1:9) .eq. 'SLOPEMAX ') then
6871 read (string,*,err=10,end=10) slpmax
6872 else if (keyword(1:7) .eq. 'ANGMAX ') then
6873 read (string,*,err=10,end=10) angmax
6874 else if (keyword(1:7) .eq. 'INTMAX ') then
6875 read (string,*,err=10,end=10) intmax
6881 !c print header information about the optimization method
6883 if (jprint .gt. 0) then
6884 if (msav .eq. 0) then
6886 20 format (/,' Steepest Descent Gradient Optimization :')
6888 30 format (/,' SD Iter F Value G RMS F Move',&
6889 ' X Move Angle FG Call Comment',/)
6892 40 format (/,' Limited Memory BFGS Quasi-Newton',&
6895 50 format (/,' QN Iter F Value G RMS F Move',&
6896 ' X Move Angle FG Call Comment',/)
6901 !c perform dynamic allocation of some local arrays
6903 allocate (x_old(nvar))
6905 allocate (g_old(nvar))
6910 if (msav .ne. 0) then
6911 allocate (rho(msav))
6912 allocate (alpha(msav))
6913 allocate (s(nvar,msav))
6914 allocate (y(nvar,msav))
6917 !c evaluate the function and get the initial gradient
6919 niter = nextiter - 1
6920 maxiter = niter + maxiter
6930 g_norm = g_norm + g(i)*g(i)
6931 g_rms = g_rms + (g(i)*scale(i))**2
6933 g_norm = sqrt(g_norm)
6934 g_rms = sqrt(g_rms) / rms
6935 f_move = 0.5d0 * stpmax * g_norm
6937 !c print initial information prior to first iteration
6939 write(jout,*) "before first"
6940 if (jprint .gt. 0) then
6941 if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then
6942 write (jout,60) niter,f,g_rms,ncalls
6943 60 format (i6,f14.4,f11.4,2x,i7)
6945 write (jout,70) niter,f,g_rms,ncalls
6946 70 format (i6,d14.4,d11.4,2x,i7)
6951 !c write initial intermediate prior to first iteration
6953 if (iwrite .gt. 0) call optsave (niter,f,x0)
6955 !c tests of the various termination criteria
6957 if (niter .ge. maxiter) then
6958 statusbf = 'IterLimit'
6961 if (f .le. fctmin) then
6962 statusbf = 'SmallFct '
6965 if (g_rms .le. grdmin) then
6966 statusbf = 'SmallGrad'
6970 !c start of a new limited memory BFGS iteration
6972 do while (.not. done)
6974 !c write (jout,*) "LBFGS niter",niter
6975 muse = min(niter-1,msav)
6977 if (m .gt. msav) m = 1
6979 !c estimate Hessian diagonal and compute the Hg product
6988 if (k .eq. 0) k = msav
6991 alpha(k) = alpha(k) + s(i,k)*q(i)
6993 alpha(k) = alpha(k) * rho(k)
6995 q(i) = q(i) - alpha(k)*y(i,k)
7004 beta = beta + y(i,k)*r(i)
7006 beta = beta * rho(k)
7008 r(i) = r(i) + s(i,k)*(alpha(k)-beta)
7011 if (k .gt. msav) k = 1
7014 !c set search direction and store current point and gradient
7022 !c perform line search along the new conjugate direction
7025 !c write (jout,*) "Before search"
7026 call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,statusbf)
7027 !c write (jout,*) "After search"
7029 !c update variables based on results of this iteration
7031 if (msav .ne. 0) then
7035 s(i,m) = x0(i) - x_old(i)
7036 y(i,m) = g(i) - g_old(i)
7037 ys = ys + y(i,m)*s(i,m)
7038 yy = yy + y(i,m)*y(i,m)
7044 !c get the sizes of the moves made during this iteration
7047 ! if (f_move.eq.0) f_move=1.0d0
7051 x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
7053 x_move = sqrt(x_move) / rms
7054 if (coordtype .eq. 'INTERNAL') then
7055 x_move = radian * x_move
7058 !c compute the rms gradient per optimization parameter
7062 g_rms = g_rms + (g(i)*scale(i))**2
7064 g_rms = sqrt(g_rms) / rms
7066 !c test for error due to line search problems
7068 if (statusbf.eq.'BadIntpln' .or. statusbf.eq.'IntplnErr') then
7070 if (nerr .ge. maxerr) done = .true.
7075 !c test for too many total iterations
7077 if (niter .ge. maxiter) then
7078 statusbf = 'IterLimit'
7082 !c test the normal termination criteria
7084 if (f .le. fctmin) then
7085 statusbf = 'SmallFct '
7088 if (g_rms .le. grdmin) then
7089 statusbf = 'SmallGrad'
7093 !c print intermediate results for the current iteration
7095 ! write(iout,*) "niter1",niter
7097 if (jprint .gt. 0) then
7098 if (done .or. mod(niter,jprint).eq.0) then
7099 if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.&
7100 g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.&
7101 f_move.gt.-1.0d5) then
7102 write (jout,80) niter,f,g_rms,f_move,x_move,&
7103 angle,ncalls,statusbf
7104 80 format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9)
7106 write (jout,90) niter,f,g_rms,f_move,x_move,&
7107 angle,ncalls,statusbf
7108 90 format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9)
7114 !c write intermediate results for the current iteration
7116 if (iwrite .gt. 0) then
7117 if (done .or. mod(niter,iwrite).eq.0) then
7118 call optsave (niter,f,x0)
7123 !c perform deallocation of some local arrays
7132 if (msav .ne. 0) then
7139 !c set final value of the objective function
7142 if (jprint .gt. 0) then
7143 if (statusbf.eq.'SmallGrad' .or. statusbf.eq.'SmallFct ') then
7144 write (jout,100) statusbf
7145 100 format (/,' LBFGS -- Normal Termination due to ',a9)
7147 write (jout,110) statusbf
7148 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9)
7155 ! double precision function funcgrad(x,g)
7157 ! include 'DIMENSIONS'
7158 ! include 'COMMON.CONTROL'
7159 ! include 'COMMON.CHAIN'
7160 ! include 'COMMON.DERIV'
7161 ! include 'COMMON.VAR'
7162 ! include 'COMMON.INTERACT'
7163 ! include 'COMMON.FFIELD'
7164 ! include 'COMMON.MD'
7165 ! include 'COMMON.QRESTR'
7166 ! include 'COMMON.IOUNITS'
7167 ! include 'COMMON.GEO'
7168 ! double precision energia(0:n_ene)
7169 ! double precision x(nvar),g(nvar)
7171 ! if (jjj.gt.0) then
7172 ! write (iout,*) "in func x"
7173 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
7175 ! call var_to_geom(nvar,x)
7177 ! call chainbuild_extconf
7178 ! call etotal(energia(0))
7180 ! funcgrad=energia(0)
7181 ! call cart2intgrad(nvar,g)
7183 ! Add the components corresponding to local energy terms.
7185 ! Add the usampl contributions
7188 ! gloc(i,icg)=gloc(i,icg)+dugamma(i)
7191 ! gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
7195 !d write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
7196 ! g(i)=g(i)+gloc(i,icg)
7201 !-----------------------------------------------------------------------------