2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
18 !-----------------------------------------------------------------------------
20 !-----------------------------------------------------------------------------
22 !-----------------------------------------------------------------------------
23 subroutine assst(iv, liv, lv, v)
25 ! *** assess candidate step (***sol version 2.3) ***
33 ! this subroutine is called by an unconstrained minimization
34 ! routine to assess the next candidate step. it may recommend one
35 ! of several courses of action, such as accepting the step, recom-
36 ! puting it using the same or a new quadratic model, or halting due
37 ! to convergence or false convergence. see the return code listing
40 !-------------------------- parameter usage --------------------------
42 ! iv (i/o) integer parameter and scratch vector -- see description
43 ! below of iv values referenced.
44 ! liv (in) length of iv array.
45 ! lv (in) length of v array.
46 ! v (i/o) real parameter and scratch vector -- see description
47 ! below of v values referenced.
49 ! *** iv values referenced ***
51 ! iv(irc) (i/o) on input for the first step tried in a new iteration,
52 ! iv(irc) should be set to 3 or 4 (the value to which it is
53 ! set when step is definitely to be accepted). on input
54 ! after step has been recomputed, iv(irc) should be
55 ! unchanged since the previous return of assst.
56 ! on output, iv(irc) is a return code having one of the
58 ! 1 = switch models or try smaller step.
59 ! 2 = switch models or accept step.
60 ! 3 = accept step and determine v(radfac) by gradient
62 ! 4 = accept step, v(radfac) has been determined.
63 ! 5 = recompute step (using the same model).
64 ! 6 = recompute step with radius = v(lmaxs) but do not
65 ! evaulate the objective function.
66 ! 7 = x-convergence (see v(xctol)).
67 ! 8 = relative function convergence (see v(rfctol)).
68 ! 9 = both x- and relative function convergence.
69 ! 10 = absolute function convergence (see v(afctol)).
70 ! 11 = singular convergence (see v(lmaxs)).
71 ! 12 = false convergence (see v(xftol)).
72 ! 13 = iv(irc) was out of range on input.
73 ! return code i has precdence over i+1 for i = 9, 10, 11.
74 ! iv(mlstgd) (i/o) saved value of iv(model).
75 ! iv(model) (i/o) on input, iv(model) should be an integer identifying
76 ! the current quadratic model of the objective function.
77 ! if a previous step yielded a better function reduction,
78 ! then iv(model) will be set to iv(mlstgd) on output.
79 ! iv(nfcall) (in) invocation count for the objective function.
80 ! iv(nfgcal) (i/o) value of iv(nfcall) at step that gave the biggest
81 ! function reduction this iteration. iv(nfgcal) remains
82 ! unchanged until a function reduction is obtained.
83 ! iv(radinc) (i/o) the number of radius increases (or minus the number
84 ! of decreases) so far this iteration.
85 ! iv(restor) (out) set to 1 if v(f) has been restored and x should be
86 ! restored to its initial value, to 2 if x should be saved,
87 ! to 3 if x should be restored from the saved value, and to
89 ! iv(stage) (i/o) count of the number of models tried so far in the
91 ! iv(stglim) (in) maximum number of models to consider.
92 ! iv(switch) (out) set to 0 unless a new model is being tried and it
93 ! gives a smaller function value than the previous model,
94 ! in which case assst sets iv(switch) = 1.
95 ! iv(toobig) (in) is nonzero if step was too big (e.g. if it caused
97 ! iv(xirc) (i/o) value that iv(irc) would have in the absence of
98 ! convergence, false convergence, and oversized steps.
100 ! *** v values referenced ***
102 ! v(afctol) (in) absolute function convergence tolerance. if the
103 ! absolute value of the current function value v(f) is less
104 ! than v(afctol), then assst returns with iv(irc) = 10.
105 ! v(decfac) (in) factor by which to decrease radius when iv(toobig) is
107 ! v(dstnrm) (in) the 2-norm of d*step.
108 ! v(dstsav) (i/o) value of v(dstnrm) on saved step.
109 ! v(dst0) (in) the 2-norm of d times the newton step (when defined,
110 ! i.e., for v(nreduc) .ge. 0).
111 ! v(f) (i/o) on both input and output, v(f) is the objective func-
112 ! tion value at x. if x is restored to a previous value,
113 ! then v(f) is restored to the corresponding value.
114 ! v(fdif) (out) the function reduction v(f0) - v(f) (for the output
115 ! value of v(f) if an earlier step gave a bigger function
116 ! decrease, and for the input value of v(f) otherwise).
117 ! v(flstgd) (i/o) saved value of v(f).
118 ! v(f0) (in) objective function value at start of iteration.
119 ! v(gtslst) (i/o) value of v(gtstep) on saved step.
120 ! v(gtstep) (in) inner product between step and gradient.
121 ! v(incfac) (in) minimum factor by which to increase radius.
122 ! v(lmaxs) (in) maximum reasonable step size (and initial step bound).
123 ! if the actual function decrease is no more than twice
124 ! what was predicted, if a return with iv(irc) = 7, 8, 9,
125 ! or 10 does not occur, if v(dstnrm) .gt. v(lmaxs), and if
126 ! v(preduc) .le. v(sctol) * abs(v(f0)), then assst re-
127 ! turns with iv(irc) = 11. if so doing appears worthwhile,
128 ! then assst repeats this test with v(preduc) computed for
129 ! a step of length v(lmaxs) (by a return with iv(irc) = 6).
130 ! v(nreduc) (i/o) function reduction predicted by quadratic model for
131 ! newton step. if assst is called with iv(irc) = 6, i.e.,
132 ! if v(preduc) has been computed with radius = v(lmaxs) for
133 ! use in the singular convervence test, then v(nreduc) is
134 ! set to -v(preduc) before the latter is restored.
135 ! v(plstgd) (i/o) value of v(preduc) on saved step.
136 ! v(preduc) (i/o) function reduction predicted by quadratic model for
138 ! v(radfac) (out) factor to be used in determining the new radius,
139 ! which should be v(radfac)*dst, where dst is either the
140 ! output value of v(dstnrm) or the 2-norm of
141 ! diag(newd)*step for the output value of step and the
142 ! updated version, newd, of the scale vector d. for
143 ! iv(irc) = 3, v(radfac) = 1.0 is returned.
144 ! v(rdfcmn) (in) minimum value for v(radfac) in terms of the input
145 ! value of v(dstnrm) -- suggested value = 0.1.
146 ! v(rdfcmx) (in) maximum value for v(radfac) -- suggested value = 4.0.
147 ! v(reldx) (in) scaled relative change in x caused by step, computed
148 ! (e.g.) by function reldst as
149 ! max (d(i)*abs(x(i)-x0(i)), 1 .le. i .le. p) /
150 ! max (d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p).
151 ! v(rfctol) (in) relative function convergence tolerance. if the
152 ! actual function reduction is at most twice what was pre-
153 ! dicted and v(nreduc) .le. v(rfctol)*abs(v(f0)), then
154 ! assst returns with iv(irc) = 8 or 9.
155 ! v(stppar) (in) marquardt parameter -- 0 means full newton step.
156 ! v(tuner1) (in) tuning constant used to decide if the function
157 ! reduction was much less than expected. suggested
159 ! v(tuner2) (in) tuning constant used to decide if the function
160 ! reduction was large enough to accept step. suggested
162 ! v(tuner3) (in) tuning constant used to decide if the radius
163 ! should be increased. suggested value = 0.75.
164 ! v(xctol) (in) x-convergence criterion. if step is a newton step
165 ! (v(stppar) = 0) having v(reldx) .le. v(xctol) and giving
166 ! at most twice the predicted function decrease, then
167 ! assst returns iv(irc) = 7 or 9.
168 ! v(xftol) (in) false convergence tolerance. if step gave no or only
169 ! a small function decrease and v(reldx) .le. v(xftol),
170 ! then assst returns with iv(irc) = 12.
172 !------------------------------- notes -------------------------------
174 ! *** application and usage restrictions ***
176 ! this routine is called as part of the nl2sol (nonlinear
177 ! least-squares) package. it may be used in any unconstrained
178 ! minimization solver that uses dogleg, goldfeld-quandt-trotter,
179 ! or levenberg-marquardt steps.
181 ! *** algorithm notes ***
183 ! see (1) for further discussion of the assessing and model
184 ! switching strategies. while nl2sol considers only two models,
185 ! assst is designed to handle any number of models.
187 ! *** usage notes ***
189 ! on the first call of an iteration, only the i/o variables
190 ! step, x, iv(irc), iv(model), v(f), v(dstnrm), v(gtstep), and
191 ! v(preduc) need have been initialized. between calls, no i/o
192 ! values execpt step, x, iv(model), v(f) and the stopping toler-
193 ! ances should be changed.
194 ! after a return for convergence or false convergence, one can
195 ! change the stopping tolerances and call assst again, in which
196 ! case the stopping tests will be repeated.
200 ! (1) dennis, j.e., jr., gay, d.m., and welsch, r.e. (1981),
201 ! an adaptive nonlinear least-squares algorithm,
202 ! acm trans. math. software, vol. 7, no. 3.
204 ! (2) powell, m.j.d. (1970) a fortran subroutine for solving
205 ! systems of nonlinear algebraic equations, in numerical
206 ! methods for nonlinear algebraic equations, edited by
207 ! p. rabinowitz, gordon and breach, london.
211 ! john dennis designed much of this routine, starting with
212 ! ideas in (2). roy welsch suggested the model switching strategy.
213 ! david gay and stephen peters cast this subroutine into a more
214 ! portable form (winter 1977), and david gay cast it into its
215 ! present form (fall 1978).
219 ! this subroutine was written in connection with research
220 ! supported by the national science foundation under grants
221 ! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
224 !------------------------ external quantities ------------------------
226 ! *** no external functions and subroutines ***
228 ! *** intrinsic functions ***
230 !el real(kind=8) :: dabs, dmax1
232 ! *** no common blocks ***
234 !-------------------------- local variables --------------------------
238 real(kind=8) :: emax, emaxs, gts, rfac1, xmax
239 !el real(kind=8) :: half, one, onep2, two, zero
241 ! *** subscripts for iv and v ***
243 !el integer :: afctol, decfac, dstnrm, dstsav, dst0, f, fdif, flstgd, f0,&
244 !el gtslst, gtstep, incfac, irc, lmaxs, mlstgd, model, nfcall,&
245 !el nfgcal, nreduc, plstgd, preduc, radfac, radinc, rdfcmn,&
246 !el rdfcmx, reldx, restor, rfctol, sctol, stage, stglim,&
247 !el stppar, switch, toobig, tuner1, tuner2, tuner3, xctol,&
251 ! *** data initializations ***
254 ! data half/0.5d+0/, one/1.d+0/, onep2/1.2d+0/, two/2.d+0/,
257 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, onep2=1.2d+0, two=2.d+0,&
262 ! data irc/29/, mlstgd/32/, model/5/, nfcall/6/, nfgcal/7/,
263 ! 1 radinc/8/, restor/9/, stage/10/, stglim/11/, switch/12/,
264 ! 2 toobig/2/, xirc/13/
266 integer,parameter :: irc=29, mlstgd=32, model=5, nfcall=6, nfgcal=7,&
267 radinc=8, restor=9, stage=10, stglim=11, switch=12,&
271 ! data afctol/31/, decfac/22/, dstnrm/2/, dst0/3/, dstsav/18/,
272 ! 1 f/10/, fdif/11/, flstgd/12/, f0/13/, gtslst/14/, gtstep/4/,
273 ! 2 incfac/23/, lmaxs/36/, nreduc/6/, plstgd/15/, preduc/7/,
274 ! 3 radfac/16/, rdfcmn/24/, rdfcmx/25/, reldx/17/, rfctol/32/,
275 ! 4 sctol/37/, stppar/5/, tuner1/26/, tuner2/27/, tuner3/28/,
276 ! 5 xctol/33/, xftol/34/
278 integer,parameter :: afctol=31, decfac=22, dstnrm=2, dst0=3, dstsav=18,&
279 f=10, fdif=11, flstgd=12, f0=13, gtslst=14, gtstep=4,&
280 incfac=23, lmaxs=36, nreduc=6, plstgd=15, preduc=7,&
281 radfac=16, rdfcmn=24, rdfcmx=25, reldx=17, rfctol=32,&
282 sctol=37, stppar=5, tuner1=26, tuner2=27, tuner3=28,&
286 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
294 if (i .ge. 1 .and. i .le. 12) &
295 go to (20,30,10,10,40,280,220,220,220,220,220,170), i
299 ! *** initialize for new iteration ***
304 if (iv(toobig) .eq. 0) go to 110
309 ! *** step was recomputed with new model or smaller radius ***
310 ! *** first decide which ***
312 20 if (iv(model) .ne. iv(mlstgd)) go to 30
313 ! *** old model retained, smaller radius tried ***
314 ! *** do not consider any more new models this iteration ***
315 iv(stage) = iv(stglim)
319 ! *** a new model is being tried. decide whether to keep it. ***
321 30 iv(stage) = iv(stage) + 1
323 ! *** now we add the possibiltiy that step was recomputed with ***
324 ! *** the same model, perhaps because of an oversized step. ***
326 40 if (iv(stage) .gt. 0) go to 50
328 ! *** step was recomputed because it was too big. ***
330 if (iv(toobig) .ne. 0) go to 60
332 ! *** restore iv(stage) and pick up where we left off. ***
334 iv(stage) = -iv(stage)
336 go to (20, 30, 110, 110, 70), i
338 50 if (iv(toobig) .eq. 0) go to 70
340 ! *** handle oversize step ***
342 if (iv(radinc) .gt. 0) go to 80
343 iv(stage) = -iv(stage)
346 60 v(radfac) = v(decfac)
347 iv(radinc) = iv(radinc) - 1
352 70 if (v(f) .lt. v(flstgd)) go to 110
354 ! *** the new step is a loser. restore old model. ***
356 if (iv(model) .eq. iv(mlstgd)) go to 80
357 iv(model) = iv(mlstgd)
360 ! *** restore step, etc. only if a previous step decreased v(f).
362 80 if (v(flstgd) .ge. v(f0)) go to 110
365 v(preduc) = v(plstgd)
366 v(gtstep) = v(gtslst)
367 if (iv(switch) .eq. 0) rfac1 = v(dstnrm) / v(dstsav)
368 v(dstnrm) = v(dstsav)
372 110 v(fdif) = v(f0) - v(f)
373 if (v(fdif) .gt. v(tuner2) * v(preduc)) go to 140
374 if(iv(radinc).gt.0) go to 140
376 ! *** no (or only a trivial) function decrease
377 ! *** -- so try new model or smaller radius
379 if (v(f) .lt. v(f0)) go to 120
380 iv(mlstgd) = iv(model)
387 if (iv(stage) .lt. iv(stglim)) go to 160
389 iv(radinc) = iv(radinc) - 1
392 ! *** nontrivial function decrease achieved ***
396 v(dstsav) = v(dstnrm)
397 if (v(fdif) .gt. v(preduc)*v(tuner1)) go to 190
399 ! *** decrease was much less than predicted -- either change models
400 ! *** or accept step with decreased radius.
402 if (iv(stage) .ge. iv(stglim)) go to 150
403 ! *** consider switching models ***
407 ! *** accept step with decreased radius ***
411 ! *** set v(radfac) to fletcher*s decrease factor ***
413 160 iv(xirc) = iv(irc)
414 emax = v(gtstep) + v(fdif)
415 v(radfac) = half * rfac1
416 if (emax .lt. v(gtstep)) v(radfac) = rfac1 * dmax1(v(rdfcmn),&
417 half * v(gtstep)/emax)
419 ! *** do false convergence test ***
421 170 if (v(reldx) .le. v(xftol)) go to 180
423 if (v(f) .lt. v(f0)) go to 200
429 ! *** handle good function decrease ***
431 190 if (v(fdif) .lt. (-v(tuner3) * v(gtstep))) go to 210
433 ! *** increasing radius looks worthwhile. see if we just
434 ! *** recomputed step with a decreased radius or restored step
435 ! *** after recomputing it with a larger radius.
437 if (iv(radinc) .lt. 0) go to 210
438 if (iv(restor) .eq. 1) go to 210
440 ! *** we did not. try a longer step unless this was a newton
443 v(radfac) = v(rdfcmx)
445 if (v(fdif) .lt. (half/v(radfac) - one) * gts) &
446 v(radfac) = dmax1(v(incfac), half*gts/(gts + v(fdif)))
448 if (v(stppar) .eq. zero) go to 230
449 if (v(dst0) .ge. zero .and. (v(dst0) .lt. two*v(dstnrm) &
450 .or. v(nreduc) .lt. onep2*v(fdif))) go to 230
451 ! *** step was not a newton step. recompute it with
452 ! *** a larger radius.
454 iv(radinc) = iv(radinc) + 1
456 ! *** save values corresponding to good step ***
459 iv(mlstgd) = iv(model)
460 if (iv(restor) .ne. 1) iv(restor) = 2
461 v(dstsav) = v(dstnrm)
463 v(plstgd) = v(preduc)
464 v(gtslst) = v(gtstep)
467 ! *** accept step with radius unchanged ***
473 ! *** come here for a restart after convergence ***
475 220 iv(irc) = iv(xirc)
476 if (v(dstsav) .ge. zero) go to 240
480 ! *** perform convergence tests ***
482 230 iv(xirc) = iv(irc)
483 240 if (iv(restor) .eq. 1 .and. v(flstgd) .lt. v(f0)) iv(restor) = 3
484 if (half * v(fdif) .gt. v(preduc)) go to 999
485 emax = v(rfctol) * dabs(v(f0))
486 emaxs = v(sctol) * dabs(v(f0))
487 if (v(dstnrm) .gt. v(lmaxs) .and. v(preduc) .le. emaxs) &
489 if (v(dst0) .lt. zero) go to 250
491 if ((v(nreduc) .gt. zero .and. v(nreduc) .le. emax) .or. &
492 (v(nreduc) .eq. zero .and. v(preduc) .eq. zero)) i = 2
493 if (v(stppar) .eq. zero .and. v(reldx) .le. v(xctol) &
494 .and. goodx) i = i + 1
495 if (i .gt. 0) iv(irc) = i + 6
497 ! *** consider recomputing step of length v(lmaxs) for singular
498 ! *** convergence test.
500 250 if (iv(irc) .gt. 5 .and. iv(irc) .ne. 12) go to 999
501 if (v(dstnrm) .gt. v(lmaxs)) go to 260
502 if (v(preduc) .ge. emaxs) go to 999
503 if (v(dst0) .le. zero) go to 270
504 if (half * v(dst0) .le. v(lmaxs)) go to 999
506 260 if (half * v(dstnrm) .le. v(lmaxs)) go to 999
507 xmax = v(lmaxs) / v(dstnrm)
508 if (xmax * (two - xmax) * v(preduc) .ge. emaxs) go to 999
509 270 if (v(nreduc) .lt. zero) go to 290
511 ! *** recompute v(preduc) for use in singular convergence test ***
513 v(gtslst) = v(gtstep)
514 v(dstsav) = v(dstnrm)
515 if (iv(irc) .eq. 12) v(dstsav) = -v(dstsav)
516 v(plstgd) = v(preduc)
519 if (i .eq. 3) iv(restor) = 0
523 ! *** perform singular convergence test with recomputed v(preduc) ***
525 280 v(gtstep) = v(gtslst)
526 v(dstnrm) = dabs(v(dstsav))
528 if (v(dstsav) .le. zero) iv(irc) = 12
529 v(nreduc) = -v(preduc)
530 v(preduc) = v(plstgd)
532 290 if (-v(nreduc) .le. v(sctol) * dabs(v(f0))) iv(irc) = 11
536 ! *** last card of assst follows ***
538 !-----------------------------------------------------------------------------
539 subroutine deflt(alg, iv, liv, lv, v)
541 ! *** supply ***sol (version 2.3) default values to iv and v ***
543 ! *** alg = 1 means regression constants.
544 ! *** alg = 2 means general unconstrained optimization constants.
547 integer :: alg, iv(liv)
548 real(kind=8) :: v(lv)
550 !el external imdcon, vdflt
552 ! imdcon... returns machine-dependent integer constants.
553 ! vdflt.... provides default values to v.
556 integer :: miniv(2), minv(2)
558 ! *** subscripts for iv ***
560 !el integer algsav, covprt, covreq, dtype, hc, ierr, inith, inits,
561 !el 1 ipivot, ivneed, lastiv, lastv, lmat, mxfcal, mxiter,
562 !el 2 nfcov, ngcov, nvdflt, outlev, parprt, parsav, perm,
563 !el 3 prunit, qrtyp, rdreq, rmat, solprt, statpr, vneed,
566 ! *** iv subscript values ***
569 ! data algsav/51/, covprt/14/, covreq/15/, dtype/16/, hc/71/,
570 ! 1 ierr/75/, inith/25/, inits/25/, ipivot/76/, ivneed/3/,
571 ! 2 lastiv/44/, lastv/45/, lmat/42/, mxfcal/17/, mxiter/18/,
572 ! 3 nfcov/52/, ngcov/53/, nvdflt/50/, outlev/19/, parprt/20/,
573 ! 4 parsav/49/, perm/58/, prunit/21/, qrtyp/80/, rdreq/57/,
574 ! 5 rmat/78/, solprt/22/, statpr/23/, vneed/4/, vsave/60/,
577 integer,parameter :: algsav=51, covprt=14, covreq=15, dtype=16, hc=71,&
578 ierr=75, inith=25, inits=25, ipivot=76, ivneed=3,&
579 lastiv=44, lastv=45, lmat=42, mxfcal=17, mxiter=18,&
580 nfcov=52, ngcov=53, nvdflt=50, outlev=19, parprt=20,&
581 parsav=49, perm=58, prunit=21, qrtyp=80, rdreq=57,&
582 rmat=78, solprt=22, statpr=23, vneed=4, vsave=60,&
585 data miniv(1)/80/, miniv(2)/59/, minv(1)/98/, minv(2)/71/
589 !------------------------------- body --------------------------------
591 if (alg .lt. 1 .or. alg .gt. 2) go to 40
593 if (liv .lt. miv) go to 20
595 if (lv .lt. mv) go to 30
596 call vdflt(alg, lv, v)
608 iv(prunit) = imdcon(1)
614 if (alg .ge. 2) go to 10
616 ! *** regression values
633 ! *** general optimization values
652 ! *** last card of deflt follows ***
654 !-----------------------------------------------------------------------------
655 real(kind=8) function dotprd(p,x,y)
657 ! *** return the inner product of the p-vectors x and y. ***
660 real(kind=8) :: x(p), y(p)
663 !el real(kind=8) :: one, zero
664 real(kind=8) :: sqteta, t
666 !el real(kind=8) :: dmax1, dabs
669 !el real(kind=8) :: rmdcon
671 ! *** rmdcon(2) returns a machine-dependent constant, sqteta, which
672 ! *** is slightly larger than the smallest positive number that
673 ! *** can be squared without underflowing.
676 ! data one/1.d+0/, sqteta/0.d+0/, zero/0.d+0/
678 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
683 if (p .le. 0) go to 999
684 !rc if (sqteta .eq. zero) sqteta = rmdcon(2)
686 !rc t = dmax1(dabs(x(i)), dabs(y(i)))
687 !rc if (t .gt. one) go to 10
688 !rc if (t .lt. sqteta) go to 20
689 !rc t = (x(i)/sqteta)*y(i)
690 !rc if (dabs(t) .lt. sqteta) go to 20
691 10 dotprd = dotprd + x(i)*y(i)
695 ! *** last card of dotprd follows ***
697 !-----------------------------------------------------------------------------
698 subroutine itsum(d, g, iv, liv, lv, p, v, x)
700 ! *** print iteration summary for ***sol (version 2.3) ***
702 ! *** parameter declarations ***
704 integer :: liv, lv, p
706 real(kind=8) :: d(p), g(p), v(lv), x(p)
708 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
710 ! *** local variables ***
712 integer :: alg, i, iv1, m, nf, ng, ol, pu
714 ! real model1(6), model2(6)
716 character(len=4) :: model1(6), model2(6)
718 real(kind=8) :: nreldf, oldf, preldf, reldf !el, zero
720 ! *** intrinsic functions ***
723 !el real(kind=8) :: dabs, dmax1
725 ! *** no external functions or subroutines ***
727 ! *** subscripts for iv and v ***
729 !el integer algsav, dstnrm, f, fdif, f0, needhd, nfcall, nfcov, ngcov,
730 !el 1 ngcall, niter, nreduc, outlev, preduc, prntit, prunit,
731 !el 2 reldx, solprt, statpr, stppar, sused, x0prt
733 ! *** iv subscript values ***
736 ! data algsav/51/, needhd/36/, nfcall/6/, nfcov/52/, ngcall/30/,
737 ! 1 ngcov/53/, niter/31/, outlev/19/, prntit/39/, prunit/21/,
738 ! 2 solprt/22/, statpr/23/, sused/64/, x0prt/24/
740 integer,parameter :: algsav=51, needhd=36, nfcall=6, nfcov=52, ngcall=30,&
741 ngcov=53, niter=31, outlev=19, prntit=39, prunit=21,&
742 solprt=22, statpr=23, sused=64, x0prt=24
745 ! *** v subscript values ***
748 ! data dstnrm/2/, f/10/, f0/13/, fdif/11/, nreduc/6/, preduc/7/,
749 ! 1 reldx/17/, stppar/5/
751 integer,parameter :: dstnrm=2, f=10, f0=13, fdif=11, nreduc=6, preduc=7,&
758 real(kind=8),parameter :: zero=0.d+0
761 ! data model1(1)/4h /, model1(2)/4h /, model1(3)/4h /,
762 ! 1 model1(4)/4h /, model1(5)/4h g /, model1(6)/4h s /,
763 ! 2 model2(1)/4h g /, model2(2)/4h s /, model2(3)/4hg-s /,
764 ! 3 model2(4)/4hs-g /, model2(5)/4h-s-g/, model2(6)/4h-g-s/
766 data model1/' ',' ',' ',' ',' g ',' s '/,&
767 model2/' g ',' s ','g-s ','s-g ','-s-g','-g-s'/
770 !------------------------------- body --------------------------------
773 if (pu .eq. 0) go to 999
775 if (iv1 .gt. 62) iv1 = iv1 - 51
778 if (iv1 .lt. 2 .or. iv1 .gt. 15) go to 370
779 if (iv1 .ge. 12) go to 120
780 if (iv1 .eq. 2 .and. iv(niter) .eq. 0) go to 390
781 if (ol .eq. 0) go to 120
782 if (iv1 .ge. 10 .and. iv(prntit) .eq. 0) go to 120
783 if (iv1 .gt. 2) go to 10
784 iv(prntit) = iv(prntit) + 1
785 if (iv(prntit) .lt. iabs(ol)) go to 999
786 10 nf = iv(nfcall) - iabs(iv(nfcov))
790 oldf = dmax1(dabs(v(f0)), dabs(v(f)))
791 if (oldf .le. zero) go to 20
792 reldf = v(fdif) / oldf
793 preldf = v(preduc) / oldf
794 20 if (ol .gt. 0) go to 60
796 ! *** print short summary line ***
798 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,30)
799 30 format(/10h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,&
801 if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,40)
802 40 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,&
805 if (alg .eq. 2) go to 50
807 write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
808 model1(m), model2(m), v(stppar)
811 50 write(pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
815 ! *** print long summary line ***
817 60 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,70)
818 70 format(/11h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,&
819 2x,13hmodel stppar,2x,6hd*step,2x,7hnpreldf)
820 if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,80)
821 80 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,&
822 3x,6hstppar,3x,6hd*step,3x,7hnpreldf)
825 if (oldf .gt. zero) nreldf = v(nreduc) / oldf
826 if (alg .eq. 2) go to 90
828 write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),&
829 model1(m), model2(m), v(stppar), v(dstnrm), nreldf
832 90 write(pu,110) iv(niter), nf, v(f), reldf, preldf,&
833 v(reldx), v(stppar), v(dstnrm), nreldf
834 100 format(i6,i5,d10.3,2d9.2,d8.1,a3,a4,2d8.1,d9.2)
835 110 format(i6,i5,d11.3,2d10.2,3d9.1,d10.2)
837 120 if (iv(statpr) .lt. 0) go to 430
838 go to (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,&
842 140 format(/26h ***** x-convergence *****)
846 160 format(/42h ***** relative function convergence *****)
850 180 format(/49h ***** x- and relative function convergence *****)
854 200 format(/42h ***** absolute function convergence *****)
858 220 format(/33h ***** singular convergence *****)
862 240 format(/30h ***** false convergence *****)
866 260 format(/38h ***** function evaluation limit *****)
870 280 format(/28h ***** iteration limit *****)
874 300 format(/18h ***** stopx *****)
878 320 format(/44h ***** initial f(x) cannot be computed *****)
883 340 format(/37h ***** bad parameters to assess *****)
887 360 format(/43h ***** gradient could not be computed *****)
888 if (iv(niter) .gt. 0) go to 480
891 370 write(pu,380) iv(1)
892 380 format(/14h ***** iv(1) =,i5,6h *****)
895 ! *** initial call on itsum ***
897 390 if (iv(x0prt) .ne. 0) write(pu,400) (i, x(i), d(i), i = 1, p)
898 400 format(/23h i initial x(i),8x,4hd(i)//(1x,i5,d17.6,d14.3))
899 ! *** the following are to avoid undefined variables when the
900 ! *** function evaluation limit is 1...
906 if (iv1 .ge. 12) go to 999
909 if (ol .eq. 0) go to 999
910 if (ol .lt. 0 .and. alg .eq. 1) write(pu,30)
911 if (ol .lt. 0 .and. alg .eq. 2) write(pu,40)
912 if (ol .gt. 0 .and. alg .eq. 1) write(pu,70)
913 if (ol .gt. 0 .and. alg .eq. 2) write(pu,80)
914 if (alg .eq. 1) write(pu,410) v(f)
915 if (alg .eq. 2) write(pu,420) v(f)
916 410 format(/11h 0 1,d10.3)
917 !365 format(/11h 0 1,e11.3)
918 420 format(/11h 0 1,d11.3)
921 ! *** print various information requested on solution ***
924 if (iv(statpr) .eq. 0) go to 480
925 oldf = dmax1(dabs(v(f0)), dabs(v(f)))
928 if (oldf .le. zero) go to 440
929 preldf = v(preduc) / oldf
930 nreldf = v(nreduc) / oldf
931 440 nf = iv(nfcall) - iv(nfcov)
932 ng = iv(ngcall) - iv(ngcov)
933 write(pu,450) v(f), v(reldx), nf, ng, preldf, nreldf
934 450 format(/9h function,d17.6,8h reldx,d17.3/12h func. evals,&
935 i8,9x,11hgrad. evals,i8/7h preldf,d16.3,6x,7hnpreldf,d15.3)
937 if (iv(nfcov) .gt. 0) write(pu,460) iv(nfcov)
938 460 format(/1x,i4,50h extra func. evals for covariance and diagnostics.)
939 if (iv(ngcov) .gt. 0) write(pu,470) iv(ngcov)
940 470 format(1x,i4,50h extra grad. evals for covariance and diagnostics.)
942 480 if (iv(solprt) .eq. 0) go to 999
945 490 format(/22h i final x(i),8x,4hd(i),10x,4hg(i)/)
947 write(pu,510) i, x(i), d(i), g(i)
949 510 format(1x,i5,d16.6,2d14.3)
953 530 format(/24h inconsistent dimensions)
955 ! *** last card of itsum follows ***
957 !-----------------------------------------------------------------------------
958 subroutine litvmu(n, x, l, y)
960 ! *** solve (l**t)*x = y, where l is an n x n lower triangular
961 ! *** matrix stored compactly by rows. x and y may occupy the same
965 !al real(kind=8) :: x(n), l(1), y(n)
966 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
967 integer :: i, ii, ij, im1, i0, j, np1
968 real(kind=8) :: xi !el, zero
972 real(kind=8),parameter :: zero=0.d+0
983 if (i .le. 1) go to 999
985 if (xi .eq. zero) go to 30
989 x(j) = x(j) - xi*l(ij)
993 ! *** last card of litvmu follows ***
994 end subroutine litvmu
995 !-----------------------------------------------------------------------------
996 subroutine livmul(n, x, l, y)
998 ! *** solve l*x = y, where l is an n x n lower triangular
999 ! *** matrix stored compactly by rows. x and y may occupy the same
1003 !al real(kind=8) :: x(n), l(1), y(n)
1004 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
1006 !el real(kind=8) :: dotprd
1008 real(kind=8) :: t !el, zero
1012 real(kind=8),parameter :: zero=0.d+0
1016 if (y(k) .ne. zero) go to 20
1022 if (k .ge. n) go to 999
1025 t = dotprd(i-1, l(j+1), x)
1027 x(i) = (y(i) - t)/l(j)
1030 ! *** last card of livmul follows ***
1031 end subroutine livmul
1032 !-----------------------------------------------------------------------------
1033 subroutine parck(alg, d, iv, liv, lv, n, v)
1035 ! *** check ***sol (version 2.3) parameters, print changed values ***
1037 ! *** alg = 1 for regression, alg = 2 for general unconstrained opt.
1039 integer :: alg, liv, lv, n
1041 real(kind=8) :: d(n), v(lv)
1043 !el external rmdcon, vcopy, vdflt
1044 !el real(kind=8) :: rmdcon
1045 ! rmdcon -- returns machine-dependent constants.
1046 ! vcopy -- copies one vector to another.
1047 ! vdflt -- supplies default parameter values to v alone.
1052 ! *** local variables ***
1054 integer :: i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu
1055 integer :: ijmp, jlim(2), miniv(2), ndflt(2)
1057 ! integer varnm(2), sh(2)
1058 ! real cngd(3), dflt(3), vn(2,34), which(3)
1060 character(len=1) :: varnm(2), sh(2)
1061 character(len=4) :: cngd(3), dflt(3), vn(2,34), which(3)
1063 real(kind=8) :: big, machep, tiny, vk, vm(34), vx(34), zero
1065 ! *** iv and v subscripts ***
1067 !el integer algsav, dinit, dtype, dtype0, epslon, inits, ivneed,
1068 !el 1 lastiv, lastv, lmat, nextiv, nextv, nvdflt, oldn,
1069 !el 2 parprt, parsav, perm, prunit, vneed
1073 ! data algsav/51/, dinit/38/, dtype/16/, dtype0/54/, epslon/19/,
1074 ! 1 inits/25/, ivneed/3/, lastiv/44/, lastv/45/, lmat/42/,
1075 ! 2 nextiv/46/, nextv/47/, nvdflt/50/, oldn/38/, parprt/20/,
1076 ! 3 parsav/49/, perm/58/, prunit/21/, vneed/4/
1078 integer,parameter :: algsav=51, dinit=38, dtype=16, dtype0=54, epslon=19,&
1079 inits=25, ivneed=3, lastiv=44, lastv=45, lmat=42,&
1080 nextiv=46, nextv=47, nvdflt=50, oldn=38, parprt=20,&
1081 parsav=49, perm=58, prunit=21, vneed=4
1082 save big, machep, tiny
1085 data big/0.d+0/, machep/-1.d+0/, tiny/1.d+0/, zero/0.d+0/
1087 ! data vn(1,1),vn(2,1)/4hepsl,4hon../
1088 ! data vn(1,2),vn(2,2)/4hphmn,4hfc../
1089 ! data vn(1,3),vn(2,3)/4hphmx,4hfc../
1090 ! data vn(1,4),vn(2,4)/4hdecf,4hac../
1091 ! data vn(1,5),vn(2,5)/4hincf,4hac../
1092 ! data vn(1,6),vn(2,6)/4hrdfc,4hmn../
1093 ! data vn(1,7),vn(2,7)/4hrdfc,4hmx../
1094 ! data vn(1,8),vn(2,8)/4htune,4hr1../
1095 ! data vn(1,9),vn(2,9)/4htune,4hr2../
1096 ! data vn(1,10),vn(2,10)/4htune,4hr3../
1097 ! data vn(1,11),vn(2,11)/4htune,4hr4../
1098 ! data vn(1,12),vn(2,12)/4htune,4hr5../
1099 ! data vn(1,13),vn(2,13)/4hafct,4hol../
1100 ! data vn(1,14),vn(2,14)/4hrfct,4hol../
1101 ! data vn(1,15),vn(2,15)/4hxcto,4hl.../
1102 ! data vn(1,16),vn(2,16)/4hxfto,4hl.../
1103 ! data vn(1,17),vn(2,17)/4hlmax,4h0.../
1104 ! data vn(1,18),vn(2,18)/4hlmax,4hs.../
1105 ! data vn(1,19),vn(2,19)/4hscto,4hl.../
1106 ! data vn(1,20),vn(2,20)/4hdini,4ht.../
1107 ! data vn(1,21),vn(2,21)/4hdtin,4hit../
1108 ! data vn(1,22),vn(2,22)/4hd0in,4hit../
1109 ! data vn(1,23),vn(2,23)/4hdfac,4h..../
1110 ! data vn(1,24),vn(2,24)/4hdltf,4hdc../
1111 ! data vn(1,25),vn(2,25)/4hdltf,4hdj../
1112 ! data vn(1,26),vn(2,26)/4hdelt,4ha0../
1113 ! data vn(1,27),vn(2,27)/4hfuzz,4h..../
1114 ! data vn(1,28),vn(2,28)/4hrlim,4hit../
1115 ! data vn(1,29),vn(2,29)/4hcosm,4hin../
1116 ! data vn(1,30),vn(2,30)/4hhube,4hrc../
1117 ! data vn(1,31),vn(2,31)/4hrspt,4hol../
1118 ! data vn(1,32),vn(2,32)/4hsigm,4hin../
1119 ! data vn(1,33),vn(2,33)/4heta0,4h..../
1120 ! data vn(1,34),vn(2,34)/4hbias,4h..../
1122 data vn(1,1),vn(2,1)/'epsl','on..'/
1123 data vn(1,2),vn(2,2)/'phmn','fc..'/
1124 data vn(1,3),vn(2,3)/'phmx','fc..'/
1125 data vn(1,4),vn(2,4)/'decf','ac..'/
1126 data vn(1,5),vn(2,5)/'incf','ac..'/
1127 data vn(1,6),vn(2,6)/'rdfc','mn..'/
1128 data vn(1,7),vn(2,7)/'rdfc','mx..'/
1129 data vn(1,8),vn(2,8)/'tune','r1..'/
1130 data vn(1,9),vn(2,9)/'tune','r2..'/
1131 data vn(1,10),vn(2,10)/'tune','r3..'/
1132 data vn(1,11),vn(2,11)/'tune','r4..'/
1133 data vn(1,12),vn(2,12)/'tune','r5..'/
1134 data vn(1,13),vn(2,13)/'afct','ol..'/
1135 data vn(1,14),vn(2,14)/'rfct','ol..'/
1136 data vn(1,15),vn(2,15)/'xcto','l...'/
1137 data vn(1,16),vn(2,16)/'xfto','l...'/
1138 data vn(1,17),vn(2,17)/'lmax','0...'/
1139 data vn(1,18),vn(2,18)/'lmax','s...'/
1140 data vn(1,19),vn(2,19)/'scto','l...'/
1141 data vn(1,20),vn(2,20)/'dini','t...'/
1142 data vn(1,21),vn(2,21)/'dtin','it..'/
1143 data vn(1,22),vn(2,22)/'d0in','it..'/
1144 data vn(1,23),vn(2,23)/'dfac','....'/
1145 data vn(1,24),vn(2,24)/'dltf','dc..'/
1146 data vn(1,25),vn(2,25)/'dltf','dj..'/
1147 data vn(1,26),vn(2,26)/'delt','a0..'/
1148 data vn(1,27),vn(2,27)/'fuzz','....'/
1149 data vn(1,28),vn(2,28)/'rlim','it..'/
1150 data vn(1,29),vn(2,29)/'cosm','in..'/
1151 data vn(1,30),vn(2,30)/'hube','rc..'/
1152 data vn(1,31),vn(2,31)/'rspt','ol..'/
1153 data vn(1,32),vn(2,32)/'sigm','in..'/
1154 data vn(1,33),vn(2,33)/'eta0','....'/
1155 data vn(1,34),vn(2,34)/'bias','....'/
1158 data vm(1)/1.0d-3/, vm(2)/-0.99d+0/, vm(3)/1.0d-3/, vm(4)/1.0d-2/,&
1159 vm(5)/1.2d+0/, vm(6)/1.d-2/, vm(7)/1.2d+0/, vm(8)/0.d+0/,&
1160 vm(9)/0.d+0/, vm(10)/1.d-3/, vm(11)/-1.d+0/, vm(13)/0.d+0/,&
1161 vm(15)/0.d+0/, vm(16)/0.d+0/, vm(19)/0.d+0/, vm(20)/-10.d+0/,&
1162 vm(21)/0.d+0/, vm(22)/0.d+0/, vm(23)/0.d+0/, vm(27)/1.01d+0/,&
1163 vm(28)/1.d+10/, vm(30)/0.d+0/, vm(31)/0.d+0/, vm(32)/0.d+0/,&
1165 data vx(1)/0.9d+0/, vx(2)/-1.d-3/, vx(3)/1.d+1/, vx(4)/0.8d+0/,&
1166 vx(5)/1.d+2/, vx(6)/0.8d+0/, vx(7)/1.d+2/, vx(8)/0.5d+0/,&
1167 vx(9)/0.5d+0/, vx(10)/1.d+0/, vx(11)/1.d+0/, vx(14)/0.1d+0/,&
1168 vx(15)/1.d+0/, vx(16)/1.d+0/, vx(19)/1.d+0/, vx(23)/1.d+0/,&
1169 vx(24)/1.d+0/, vx(25)/1.d+0/, vx(26)/1.d+0/, vx(27)/1.d+10/,&
1170 vx(29)/1.d+0/, vx(31)/1.d+0/, vx(32)/1.d+0/, vx(33)/1.d+0/,&
1174 ! data varnm(1)/1hp/, varnm(2)/1hn/, sh(1)/1hs/, sh(2)/1hh/
1175 ! data cngd(1),cngd(2),cngd(3)/4h---c,4hhang,4hed v/,
1176 ! 1 dflt(1),dflt(2),dflt(3)/4hnond,4hefau,4hlt v/
1178 data varnm(1)/'p'/, varnm(2)/'n'/, sh(1)/'s'/, sh(2)/'h'/
1179 data cngd(1),cngd(2),cngd(3)/'---c','hang','ed v'/,&
1180 dflt(1),dflt(2),dflt(3)/'nond','efau','lt v'/
1182 data ijmp/33/, jlim(1)/0/, jlim(2)/24/, ndflt(1)/32/, ndflt(2)/25/
1183 data miniv(1)/80/, miniv(2)/59/
1185 !............................... body ................................
1188 if (prunit .le. liv) pu = iv(prunit)
1189 if (alg .lt. 1 .or. alg .gt. 2) go to 340
1190 if (iv(1) .eq. 0) call deflt(alg, iv, liv, lv, v)
1192 if (iv1 .ne. 13 .and. iv1 .ne. 12) go to 10
1194 if (perm .le. liv) miv1 = max0(miv1, iv(perm) - 1)
1195 if (ivneed .le. liv) miv2 = miv1 + max0(iv(ivneed), 0)
1196 if (lastiv .le. liv) iv(lastiv) = miv2
1197 if (liv .lt. miv1) go to 300
1199 iv(lastv) = max0(iv(vneed), 0) + iv(lmat) - 1
1201 if (liv .lt. miv2) go to 300
1202 if (lv .lt. iv(lastv)) go to 320
1203 10 if (alg .eq. iv(algsav)) go to 30
1204 if (pu .ne. 0) write(pu,20) alg, iv(algsav)
1205 20 format(/39h the first parameter to deflt should be,i3,&
1209 30 if (iv1 .lt. 12 .or. iv1 .gt. 14) go to 60
1210 if (n .ge. 1) go to 50
1212 if (pu .eq. 0) go to 999
1213 write(pu,40) varnm(alg), n
1214 40 format(/8h /// bad,a1,2h =,i5)
1216 50 if (iv1 .ne. 14) iv(nextiv) = iv(perm)
1217 if (iv1 .ne. 14) iv(nextv) = iv(lmat)
1218 if (iv1 .eq. 13) go to 999
1219 k = iv(parsav) - epslon
1220 call vdflt(alg, lv-k, v(k+1))
1221 iv(dtype0) = 2 - alg
1227 60 if (n .eq. iv(oldn)) go to 80
1229 if (pu .eq. 0) go to 999
1230 write(pu,70) varnm(alg), iv(oldn), n
1231 70 format(/5h /// ,1a1,14h changed from ,i5,4h to ,i5)
1234 80 if (iv1 .le. 11 .and. iv1 .ge. 1) go to 100
1236 if (pu .ne. 0) write(pu,90) iv1
1237 90 format(/13h /// iv(1) =,i5,28h should be between 0 and 14.)
1240 100 which(1) = cngd(1)
1244 110 if (iv1 .eq. 14) iv1 = 12
1245 if (big .gt. tiny) go to 120
1272 do 150 l = 1, ndfalt
1274 if (vk .ge. vm(i) .and. vk .le. vx(i)) go to 140
1276 if (pu .ne. 0) write(pu,130) vn(1,i), vn(2,i), k, vk,&
1278 130 format(/6h /// ,2a4,5h.. v(,i2,3h) =,d11.3,7h should,&
1279 11h be between,d11.3,4h and,d11.3)
1282 if (i .eq. j) i = ijmp
1285 if (iv(nvdflt) .eq. ndfalt) go to 170
1287 if (pu .eq. 0) go to 999
1288 write(pu,160) iv(nvdflt), ndfalt
1289 160 format(/13h iv(nvdflt) =,i5,13h rather than ,i5)
1291 170 if ((iv(dtype) .gt. 0 .or. v(dinit) .gt. zero) .and. iv1 .eq. 12) &
1294 if (d(i) .gt. zero) go to 190
1296 if (pu .ne. 0) write(pu,180) i, d(i)
1297 180 format(/8h /// d(,i3,3h) =,d11.3,19h should be positive)
1299 200 if (m .eq. 0) go to 210
1303 210 if (pu .eq. 0 .or. iv(parprt) .eq. 0) go to 999
1304 if (iv1 .ne. 12 .or. iv(inits) .eq. alg-1) go to 230
1306 write(pu,220) sh(alg), iv(inits)
1307 220 format(/22h nondefault values..../5h init,a1,14h..... iv(25) =,&
1309 230 if (iv(dtype) .eq. iv(dtype0)) go to 250
1310 if (m .eq. 0) write(pu,260) which
1312 write(pu,240) iv(dtype)
1313 240 format(20h dtype..... iv(16) =,i3)
1319 do 290 ii = 1, ndfalt
1320 if (v(k) .eq. v(l)) go to 280
1321 if (m .eq. 0) write(pu,260) which
1322 260 format(/1h ,3a4,9halues..../)
1324 write(pu,270) vn(1,i), vn(2,i), k, v(k)
1325 270 format(1x,2a4,5h.. v(,i2,3h) =,d15.7)
1329 if (i .eq. j) i = ijmp
1332 iv(dtype0) = iv(dtype)
1334 call vcopy(iv(nvdflt), v(parsv1), v(epslon))
1338 if (pu .eq. 0) go to 999
1339 write(pu,310) liv, miv2
1340 310 format(/10h /// liv =,i5,17h must be at least,i5)
1341 if (liv .lt. miv1) go to 999
1342 if (lv .lt. iv(lastv)) go to 320
1346 if (pu .eq. 0) go to 999
1347 write(pu,330) lv, iv(lastv)
1348 330 format(/9h /// lv =,i5,17h must be at least,i5)
1352 if (pu .eq. 0) go to 999
1354 350 format(/10h /// alg =,i5,15h must be 1 or 2)
1357 ! *** last card of parck follows ***
1358 end subroutine parck
1359 !-----------------------------------------------------------------------------
1360 real(kind=8) function reldst(p, d, x, x0)
1362 ! *** compute and return relative difference between x and x0 ***
1363 ! *** nl2sol version 2.2 ***
1366 real(kind=8) :: d(p), x(p), x0(p)
1368 !el real(kind=8) :: dabs
1371 real(kind=8) :: emax, t, xmax !el, zero
1375 real(kind=8),parameter :: zero=0.d+0
1381 t = dabs(d(i) * (x(i) - x0(i)))
1382 if (emax .lt. t) emax = t
1383 t = d(i) * (dabs(x(i)) + dabs(x0(i)))
1384 if (xmax .lt. t) xmax = t
1387 if (xmax .gt. zero) reldst = emax / xmax
1389 ! *** last card of reldst follows ***
1391 !-----------------------------------------------------------------------------
1392 subroutine vaxpy(p, w, a, x, y)
1394 ! *** set w = a*x + y -- w, x, y = p-vectors, a = scalar ***
1397 real(kind=8) :: a, w(p), x(p), y(p)
1402 10 w(i) = a*x(i) + y(i)
1404 end subroutine vaxpy
1405 !-----------------------------------------------------------------------------
1406 subroutine vcopy(p, y, x)
1408 ! *** set y = x, where x and y are p-vectors ***
1411 real(kind=8) :: x(p), y(p)
1418 end subroutine vcopy
1419 !-----------------------------------------------------------------------------
1420 subroutine vdflt(alg, lv, v)
1422 ! *** supply ***sol (version 2.3) default values to v ***
1424 ! *** alg = 1 means regression constants.
1425 ! *** alg = 2 means general unconstrained optimization constants.
1427 integer :: alg, l,lv
1428 real(kind=8) :: v(lv)
1430 !el real(kind=8) :: dmax1
1433 !el real(kind=8) :: rmdcon
1434 ! rmdcon... returns machine-dependent constants
1436 real(kind=8) :: machep, mepcrt, sqteps !el one, three
1438 ! *** subscripts for v ***
1440 !el integer afctol, bias, cosmin, decfac, delta0, dfac, dinit, dltfdc,
1441 !el 1 dltfdj, dtinit, d0init, epslon, eta0, fuzz, huberc,
1442 !el 2 incfac, lmax0, lmaxs, phmnfc, phmxfc, rdfcmn, rdfcmx,
1443 !el 3 rfctol, rlimit, rsptol, sctol, sigmin, tuner1, tuner2,
1444 !el 4 tuner3, tuner4, tuner5, xctol, xftol
1447 ! data one/1.d+0/, three/3.d+0/
1449 real(kind=8),parameter :: one=1.d+0, three=3.d+0
1452 ! *** v subscript values ***
1455 ! data afctol/31/, bias/43/, cosmin/47/, decfac/22/, delta0/44/,
1456 ! 1 dfac/41/, dinit/38/, dltfdc/42/, dltfdj/43/, dtinit/39/,
1457 ! 2 d0init/40/, epslon/19/, eta0/42/, fuzz/45/, huberc/48/,
1458 ! 3 incfac/23/, lmax0/35/, lmaxs/36/, phmnfc/20/, phmxfc/21/,
1459 ! 4 rdfcmn/24/, rdfcmx/25/, rfctol/32/, rlimit/46/, rsptol/49/,
1460 ! 5 sctol/37/, sigmin/50/, tuner1/26/, tuner2/27/, tuner3/28/,
1461 ! 6 tuner4/29/, tuner5/30/, xctol/33/, xftol/34/
1463 integer,parameter :: afctol=31, bias=43, cosmin=47, decfac=22, delta0=44,&
1464 dfac=41, dinit=38, dltfdc=42, dltfdj=43, dtinit=39,&
1465 d0init=40, epslon=19, eta0=42, fuzz=45, huberc=48,&
1466 incfac=23, lmax0=35, lmaxs=36, phmnfc=20, phmxfc=21,&
1467 rdfcmn=24, rdfcmx=25, rfctol=32, rlimit=46, rsptol=49,&
1468 sctol=37, sigmin=50, tuner1=26, tuner2=27, tuner3=28,&
1469 tuner4=29, tuner5=30, xctol=33, xftol=34
1472 !------------------------------- body --------------------------------
1476 if (machep .gt. 1.d-10) v(afctol) = machep**2
1482 mepcrt = machep ** (one/three)
1492 v(rfctol) = dmax1(1.d-10, mepcrt**2)
1493 v(sctol) = v(rfctol)
1500 v(xftol) = 1.d+2 * machep
1502 if (alg .ge. 2) go to 10
1504 ! *** regression values
1506 v(cosmin) = dmax1(1.d-6, 1.d+2 * machep)
1512 v(rlimit) = rmdcon(5)
1517 ! *** general optimization values
1521 v(eta0) = 1.0d+3 * machep
1524 ! *** last card of vdflt follows ***
1525 end subroutine vdflt
1526 !-----------------------------------------------------------------------------
1527 subroutine vscopy(p, y, s)
1529 ! *** set p-vector y to scalar s ***
1532 real(kind=8) :: s, y(p)
1539 end subroutine vscopy
1540 !-----------------------------------------------------------------------------
1541 real(kind=8) function v2norm(p, x)
1543 ! *** return the 2-norm of the p-vector x, taking ***
1544 ! *** care to avoid the most likely underflows. ***
1547 real(kind=8) :: x(p)
1550 real(kind=8) :: r, scale, sqteta, t, xi !el, one, zero
1552 !el real(kind=8) :: dabs, dsqrt
1555 !el real(kind=8) :: rmdcon
1558 ! data one/1.d+0/, zero/0.d+0/
1560 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
1565 if (p .gt. 0) go to 10
1569 if (x(i) .ne. zero) go to 30
1574 30 scale = dabs(x(i))
1575 if (i .lt. p) go to 40
1579 if (sqteta .eq. zero) sqteta = rmdcon(2)
1581 ! *** sqteta is (slightly larger than) the square root of the
1582 ! *** smallest positive floating point number on the machine.
1583 ! *** the tests involving sqteta are done to prevent underflows.
1588 if (xi .gt. scale) go to 50
1590 if (r .gt. sqteta) t = t + r*r
1593 if (r .le. sqteta) r = zero
1598 v2norm = scale * dsqrt(t)
1600 ! *** last card of v2norm follows ***
1602 !-----------------------------------------------------------------------------
1603 subroutine humsl(n,d,x,calcf,calcgh,iv,liv,lv,v,uiparm,urparm,ufparm)
1605 ! *** minimize general unconstrained objective function using ***
1606 ! *** (analytic) gradient and hessian provided by the caller. ***
1608 integer :: liv, lv, n
1609 integer :: iv(liv), uiparm(1)
1610 real(kind=8) :: d(n), x(n), v(lv), urparm(1)
1611 real(kind=8),external :: ufparm
1612 ! dimension v(78 + n*(n+12)), uiparm(*), urparm(*)
1613 external :: calcf, calcgh
1615 !------------------------------ discussion ---------------------------
1617 ! this routine is like sumsl, except that the subroutine para-
1618 ! meter calcg of sumsl (which computes the gradient of the objec-
1619 ! tive function) is replaced by the subroutine parameter calcgh,
1620 ! which computes both the gradient and (lower triangle of the)
1621 ! hessian of the objective function. the calling sequence is...
1622 ! call calcgh(n, x, nf, g, h, uiparm, urparm, ufparm)
1623 ! parameters n, x, nf, g, uiparm, urparm, and ufparm are the same
1624 ! as for sumsl, while h is an array of length n*(n+1)/2 in which
1625 ! calcgh must store the lower triangle of the hessian at x. start-
1626 ! ing at h(1), calcgh must store the hessian entries in the order
1627 ! (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ...
1628 ! the value printed (by itsum) in the column labelled stppar
1629 ! is the levenberg-marquardt used in computing the current step.
1630 ! zero means a full newton step. if the special case described in
1631 ! ref. 1 is detected, then stppar is negated. the value printed
1632 ! in the column labelled npreldf is zero if the current hessian
1633 ! is not positive definite.
1634 ! it sometimes proves worthwhile to let d be determined from the
1635 ! diagonal of the hessian matrix by setting iv(dtype) = 1 and
1636 ! v(dinit) = 0. the following iv and v components are relevant...
1638 ! iv(dtol)..... iv(59) gives the starting subscript in v of the dtol
1639 ! array used when d is updated. (iv(dtol) can be
1640 ! initialized by calling humsl with iv(1) = 13.)
1641 ! iv(dtype).... iv(16) tells how the scale vector d should be chosen.
1642 ! iv(dtype) .le. 0 means that d should not be updated, and
1643 ! iv(dtype) .ge. 1 means that d should be updated as
1644 ! described below with v(dfac). default = 0.
1645 ! v(dfac)..... v(41) and the dtol and d0 arrays (see v(dtinit) and
1646 ! v(d0init)) are used in updating the scale vector d when
1647 ! iv(dtype) .gt. 0. (d is initialized according to
1648 ! v(dinit), described in sumsl.) let
1649 ! d1(i) = max(sqrt(abs(h(i,i))), v(dfac)*d(i)),
1650 ! where h(i,i) is the i-th diagonal element of the current
1651 ! hessian. if iv(dtype) = 1, then d(i) is set to d1(i)
1652 ! unless d1(i) .lt. dtol(i), in which case d(i) is set to
1653 ! max(d0(i), dtol(i)).
1654 ! if iv(dtype) .ge. 2, then d is updated during the first
1655 ! iteration as for iv(dtype) = 1 (after any initialization
1656 ! due to v(dinit)) and is left unchanged thereafter.
1658 ! v(dtinit)... v(39), if positive, is the value to which all components
1659 ! of the dtol array (see v(dfac)) are initialized. if
1660 ! v(dtinit) = 0, then it is assumed that the caller has
1661 ! stored dtol in v starting at v(iv(dtol)).
1663 ! v(d0init)... v(40), if positive, is the value to which all components
1664 ! of the d0 vector (see v(dfac)) are initialized. if
1665 ! v(dfac) = 0, then it is assumed that the caller has
1666 ! stored d0 in v starting at v(iv(dtol)+n). default = 1.0.
1670 ! 1. gay, d.m. (1981), computing optimal locally constrained steps,
1671 ! siam j. sci. statist. comput. 2, pp. 186-197.
1675 ! coded by david m. gay (winter 1980). revised sept. 1982.
1676 ! this subroutine was written in connection with research supported
1677 ! in part by the national science foundation under grants
1678 ! mcs-7600324 and mcs-7906671.
1680 !---------------------------- declarations ---------------------------
1682 !el external deflt, humit
1684 ! deflt... provides default input values for iv and v.
1685 ! humit... reverse-communication routine that does humsl algorithm.
1687 integer :: g1, h1, iv1, lh, nf
1690 ! *** subscripts for iv ***
1692 !el integer g, h, nextv, nfcall, nfgcal, toobig, vneed
1695 ! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, h/56/, toobig/2/,
1698 integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28, h=56,&
1702 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
1704 lh = n * (n + 1) / 2
1705 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1706 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
1707 iv(vneed) = iv(vneed) + n*(n+3)/2
1709 if (iv1 .eq. 14) go to 10
1710 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
1713 if (iv1 .eq. 12) iv(1) = 13
1719 20 call humit(d, f, v(g1), v(h1), iv, lh, liv, lv, n, v, x)
1720 if (iv(1) - 2) 30, 40, 50
1723 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
1724 if (nf .le. 0) iv(toobig) = 1
1727 40 call calcgh(n, x, iv(nfgcal), v(g1), v(h1), uiparm, urparm,&
1731 50 if (iv(1) .ne. 14) go to 999
1733 ! *** storage allocation
1737 iv(nextv) = iv(h) + n*(n+1)/2
1738 if (iv1 .ne. 13) go to 10
1741 ! *** last card of humsl follows ***
1742 end subroutine humsl
1743 !-----------------------------------------------------------------------------
1744 subroutine humit(d, fx, g, h, iv, lh, liv, lv, n, v, x)
1746 ! *** carry out humsl (unconstrained minimization) iterations, using
1747 ! *** hessian matrix provided by the caller.
1750 use control, only:stopx
1752 ! *** parameter declarations ***
1754 integer :: lh, liv, lv, n
1756 real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n)
1758 !-------------------------- parameter usage --------------------------
1760 ! d.... scale vector.
1761 ! fx... function value.
1762 ! g.... gradient vector.
1763 ! h.... lower triangle of the hessian, stored rowwise.
1764 ! iv... integer value array.
1765 ! lh... length of h = p*(p+1)/2.
1766 ! liv.. length of iv (at least 60).
1767 ! lv... length of v (at least 78 + n*(n+21)/2).
1768 ! n.... number of variables (components in x and g).
1769 ! v.... floating-point value array.
1770 ! x.... parameter vector.
1772 ! *** discussion ***
1774 ! parameters iv, n, v, and x are the same as the corresponding
1775 ! ones to humsl (which see), except that v can be shorter (since
1776 ! the part of v that humsl uses for storing g and h is not needed).
1777 ! moreover, compared with humsl, iv(1) may have the two additional
1778 ! output values 1 and 2, which are explained below, as is the use
1779 ! of iv(toobig) and iv(nfgcal). the value iv(g), which is an
1780 ! output value from humsl, is not referenced by humit or the
1781 ! subroutines it calls.
1783 ! iv(1) = 1 means the caller should set fx to f(x), the function value
1784 ! at x, and call humit again, having changed none of the
1785 ! other parameters. an exception occurs if f(x) cannot be
1786 ! computed (e.g. if overflow would occur), which may happen
1787 ! because of an oversized step. in this case the caller
1788 ! should set iv(toobig) = iv(2) to 1, which will cause
1789 ! humit to ignore fx and try a smaller step. the para-
1790 ! meter nf that humsl passes to calcf (for possible use by
1791 ! calcgh) is a copy of iv(nfcall) = iv(6).
1792 ! iv(1) = 2 means the caller should set g to g(x), the gradient of f at
1793 ! x, and h to the lower triangle of h(x), the hessian of f
1794 ! at x, and call humit again, having changed none of the
1795 ! other parameters except perhaps the scale vector d.
1796 ! the parameter nf that humsl passes to calcg is
1797 ! iv(nfgcal) = iv(7). if g(x) and h(x) cannot be evaluated,
1798 ! then the caller may set iv(nfgcal) to 0, in which case
1799 ! humit will return with iv(1) = 65.
1800 ! note -- humit overwrites h with the lower triangle
1801 ! of diag(d)**-1 * h(x) * diag(d)**-1.
1805 ! coded by david m. gay (winter 1980). revised sept. 1982.
1806 ! this subroutine was written in connection with research supported
1807 ! in part by the national science foundation under grants
1808 ! mcs-7600324 and mcs-7906671.
1810 ! (see sumsl and humsl for references.)
1812 !+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
1814 ! *** local variables ***
1816 integer :: dg1, dummy, i, j, k, l, lstgst, nn1o2, step1,&
1822 !el real(kind=8) :: one, onep2, zero
1824 ! *** no intrinsic functions ***
1826 ! *** external functions and subroutines ***
1828 !el external assst, deflt, dotprd, dupdu, gqtst, itsum, parck,
1829 !el 1 reldst, slvmul, stopx, vaxpy, vcopy, vscopy, v2norm
1831 !el real(kind=8) :: dotprd, reldst, v2norm
1833 ! assst.... assesses candidate step.
1834 ! deflt.... provides default iv and v input values.
1835 ! dotprd... returns inner product of two vectors.
1836 ! dupdu.... updates scale vector d.
1837 ! gqtst.... computes optimally locally constrained step.
1838 ! itsum.... prints iteration summary and info on initial and final x.
1839 ! parck.... checks validity of input iv and v values.
1840 ! reldst... computes v(reldx) = relative step size.
1841 ! slvmul... multiplies symmetric matrix times vector, given the lower
1842 ! triangle of the matrix.
1843 ! stopx.... returns .true. if the break key has been pressed.
1844 ! vaxpy.... computes scalar times one vector plus another.
1845 ! vcopy.... copies one vector to another.
1846 ! vscopy... sets all elements of a vector to a scalar.
1847 ! v2norm... returns the 2-norm of a vector.
1849 ! *** subscripts for iv and v ***
1851 !el integer cnvcod, dg, dgnorm, dinit, dstnrm, dtinit, dtol,
1852 !el 1 dtype, d0init, f, f0, fdif, gtstep, incfac, irc, kagqt,
1853 !el 2 lmat, lmax0, lmaxs, mode, model, mxfcal, mxiter, nextv,
1854 !el 3 nfcall, nfgcal, ngcall, niter, preduc, radfac, radinc,
1855 !el 4 radius, rad0, reldx, restor, step, stglim, stlstg, stppar,
1856 !el 5 toobig, tuner4, tuner5, vneed, w, xirc, x0
1858 ! *** iv subscript values ***
1861 ! data cnvcod/55/, dg/37/, dtol/59/, dtype/16/, irc/29/, kagqt/33/,
1862 ! 1 lmat/42/, mode/35/, model/5/, mxfcal/17/, mxiter/18/,
1863 ! 2 nextv/47/, nfcall/6/, nfgcal/7/, ngcall/30/, niter/31/,
1864 ! 3 radinc/8/, restor/9/, step/40/, stglim/11/, stlstg/41/,
1865 ! 4 toobig/2/, vneed/4/, w/34/, xirc/13/, x0/43/
1867 integer,parameter :: cnvcod=55, dg=37, dtol=59, dtype=16, irc=29, kagqt=33,&
1868 lmat=42, mode=35, model=5, mxfcal=17, mxiter=18,&
1869 nextv=47, nfcall=6, nfgcal=7, ngcall=30, niter=31,&
1870 radinc=8, restor=9, step=40, stglim=11, stlstg=41,&
1871 toobig=2, vneed=4, w=34, xirc=13, x0=43
1874 ! *** v subscript values ***
1877 ! data dgnorm/1/, dinit/38/, dstnrm/2/, dtinit/39/, d0init/40/,
1878 ! 1 f/10/, f0/13/, fdif/11/, gtstep/4/, incfac/23/, lmax0/35/,
1879 ! 2 lmaxs/36/, preduc/7/, radfac/16/, radius/8/, rad0/9/,
1880 ! 3 reldx/17/, stppar/5/, tuner4/29/, tuner5/30/
1882 integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dtinit=39, d0init=40,&
1883 f=10, f0=13, fdif=11, gtstep=4, incfac=23, lmax0=35,&
1884 lmaxs=36, preduc=7, radfac=16, radius=8, rad0=9,&
1885 reldx=17, stppar=5, tuner4=29, tuner5=30
1889 ! data one/1.d+0/, onep2/1.2d+0/, zero/0.d+0/
1891 real(kind=8),parameter :: one=1.d+0, onep2=1.2d+0, zero=0.d+0
1894 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
1897 if (i .eq. 1) go to 30
1898 if (i .eq. 2) go to 40
1900 ! *** check validity of iv and v input values ***
1902 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1903 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
1904 iv(vneed) = iv(vneed) + n*(n+21)/2 + 7
1905 call parck(2, d, iv, liv, lv, n, v)
1907 if (i .gt. 12) go to 999
1908 nn1o2 = n * (n + 1) / 2
1909 if (lh .ge. nn1o2) go to (210,210,210,210,210,210,160,120,160,&
1914 ! *** storage allocation ***
1916 10 iv(dtol) = iv(lmat) + nn1o2
1917 iv(x0) = iv(dtol) + 2*n
1918 iv(step) = iv(x0) + n
1919 iv(stlstg) = iv(step) + n
1920 iv(dg) = iv(stlstg) + n
1922 iv(nextv) = iv(w) + 4*n + 7
1923 if (iv(1) .ne. 13) go to 20
1927 ! *** initialization ***
1941 if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
1943 if (v(dtinit) .gt. zero) call vscopy(n, v(k), v(dtinit))
1945 if (v(d0init) .gt. zero) call vscopy(n, v(k), v(d0init))
1950 if (iv(mode) .ge. 0) go to 210
1952 if (iv(toobig) .eq. 0) go to 999
1956 ! *** make sure gradient could be computed ***
1958 40 if (iv(nfgcal) .ne. 0) go to 50
1962 ! *** update the scale vector d ***
1965 if (iv(dtype) .le. 0) go to 70
1973 call dupdu(d, v(dg1), iv, liv, lv, n, v)
1975 ! *** compute scaled gradient and its norm ***
1983 v(dgnorm) = v2norm(n, v(dg1))
1985 ! *** compute scaled hessian ***
1991 h(k) = t * h(k) / d(j)
1996 if (iv(cnvcod) .ne. 0) go to 340
1997 if (iv(mode) .eq. 0) go to 300
1999 ! *** allow first step to have scaled 2-norm at most v(lmax0) ***
2001 v(radius) = v(lmax0)
2006 !----------------------------- main loop -----------------------------
2009 ! *** print iteration summary, check iteration limit ***
2011 110 call itsum(d, g, iv, liv, lv, n, v, x)
2013 if (k .lt. iv(mxiter)) go to 130
2017 130 iv(niter) = k + 1
2019 ! *** initialize for start of next iteration ***
2027 ! *** copy x to x0 ***
2029 call vcopy(n, v(x01), x)
2031 ! *** update radius ***
2033 if (k .eq. 0) go to 150
2040 v(radius) = v(radfac) * v2norm(n, v(step1))
2042 ! *** check stopx and function evaluation limit ***
2046 150 if (.not. stopx(dummy)) go to 170
2050 ! *** come here when restarting after func. eval. limit or stopx.
2052 160 if (v(f) .ge. v(f0)) go to 170
2057 170 if (iv(nfcall) .lt. iv(mxfcal)) go to 190
2059 180 if (v(f) .ge. v(f0)) go to 350
2061 ! *** in case of stopx or function evaluation limit with
2062 ! *** improved v(f), evaluate the gradient at x.
2067 !. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
2069 190 step1 = iv(step)
2073 call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1))
2074 if (iv(irc) .eq. 6) go to 210
2076 ! *** check whether evaluating f(x0 + step) looks worthwhile ***
2078 if (v(dstnrm) .le. zero) go to 210
2079 if (iv(irc) .ne. 5) go to 200
2080 if (v(radfac) .le. one) go to 200
2081 if (v(preduc) .le. onep2 * v(fdif)) go to 210
2083 ! *** compute f(x0 + step) ***
2087 call vaxpy(n, x, one, v(step1), v(x01))
2088 iv(nfcall) = iv(nfcall) + 1
2093 !. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
2096 v(reldx) = reldst(n, d, x, v(x01))
2097 call assst(iv, liv, lv, v)
2100 if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
2101 if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
2102 if (iv(restor) .ne. 3) go to 220
2103 call vcopy(n, v(step1), v(lstgst))
2104 call vaxpy(n, x, one, v(step1), v(x01))
2105 v(reldx) = reldst(n, d, x, v(x01))
2108 go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k
2110 ! *** recompute step with new radius ***
2112 230 v(radius) = v(radfac) * v(dstnrm)
2115 ! *** compute step of length v(lmaxs) for singular convergence test.
2117 240 v(radius) = v(lmaxs)
2120 ! *** convergence or false convergence ***
2122 250 iv(cnvcod) = k - 4
2123 if (v(f) .ge. v(f0)) go to 340
2124 if (iv(xirc) .eq. 14) go to 340
2127 !. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
2129 260 if (iv(irc) .ne. 3) go to 290
2132 ! *** prepare for gradient tests ***
2133 ! *** set temp1 = hessian * step + g(x0)
2134 ! *** = diag(d) * (h * step + g(x0))
2136 ! use x0 vector as temporary.
2139 v(k) = d(i) * v(step1)
2143 call slvmul(n, v(temp1), h, v(x01))
2145 v(temp1) = d(i) * v(temp1) + g(i)
2149 ! *** compute gradient and hessian ***
2151 290 iv(ngcall) = iv(ngcall) + 1
2156 if (iv(irc) .ne. 3) go to 110
2158 ! *** set v(radfac) by gradient tests ***
2163 ! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
2167 v(k) = (v(k) - g(i)) / d(i)
2171 ! *** do gradient tests ***
2173 if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320
2174 if (dotprd(n, g, v(step1)) &
2175 .ge. v(gtstep) * v(tuner5)) go to 110
2176 320 v(radfac) = v(incfac)
2179 !. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
2181 ! *** bad parameters to assess ***
2186 ! *** print summary of final iteration and other requested items ***
2188 340 iv(1) = iv(cnvcod)
2190 350 call itsum(d, g, iv, liv, lv, n, v, x)
2194 ! *** last card of humit follows ***
2195 end subroutine humit
2196 !-----------------------------------------------------------------------------
2197 subroutine dupdu(d, hdiag, iv, liv, lv, n, v)
2199 ! *** update scale vector d for humsl ***
2201 ! *** parameter declarations ***
2203 integer :: liv, lv, n
2205 real(kind=8) :: d(n), hdiag(n), v(lv)
2207 ! *** local variables ***
2209 integer :: dtoli, d0i, i
2210 real(kind=8) :: t, vdfac
2212 ! *** intrinsic functions ***
2214 !el real(kind=8) :: dabs, dmax1, dsqrt
2216 ! *** subscripts for iv and v ***
2218 !el integer :: dfac, dtol, dtype, niter
2220 ! data dfac/41/, dtol/59/, dtype/16/, niter/31/
2222 integer,parameter :: dfac=41, dtol=59, dtype=16, niter=31
2225 !------------------------------- body --------------------------------
2228 if (i .eq. 1) go to 10
2229 if (iv(niter) .gt. 0) go to 999
2235 t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i))
2236 if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i))
2243 ! *** last card of dupdu follows ***
2244 end subroutine dupdu
2245 !-----------------------------------------------------------------------------
2246 subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w)
2248 ! *** compute goldfeld-quandt-trotter step by more-hebden technique ***
2249 ! *** (nl2sol version 2.2), modified a la more and sorensen ***
2251 ! *** parameter declarations ***
2254 !al real(kind=8) :: d(p), dig(p), dihdi(1), l(1), v(21), step(p),
2256 real(kind=8) :: d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),&
2257 v(21), step(p),w(4*p+7)
2258 ! dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7)
2260 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2264 ! given the (compactly stored) lower triangle of a scaled
2265 ! hessian (approximation) and a nonzero scaled gradient vector,
2266 ! this subroutine computes a goldfeld-quandt-trotter step of
2267 ! approximate length v(radius) by the more-hebden technique. in
2268 ! other words, step is computed to (approximately) minimize
2269 ! psi(step) = (g**t)*step + 0.5*(step**t)*h*step such that the
2270 ! 2-norm of d*step is at most (approximately) v(radius), where
2271 ! g is the gradient, h is the hessian, and d is a diagonal
2272 ! scale matrix whose diagonal is stored in the parameter d.
2273 ! (gqtst assumes dig = d**-1 * g and dihdi = d**-1 * h * d**-1.)
2275 ! *** parameter description ***
2277 ! d (in) = the scale vector, i.e. the diagonal of the scale
2278 ! matrix d mentioned above under purpose.
2279 ! dig (in) = the scaled gradient vector, d**-1 * g. if g = 0, then
2280 ! step = 0 and v(stppar) = 0 are returned.
2281 ! dihdi (in) = lower triangle of the scaled hessian (approximation),
2282 ! i.e., d**-1 * h * d**-1, stored compactly by rows., i.e.,
2283 ! in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc.
2284 ! ka (i/o) = the number of hebden iterations (so far) taken to deter-
2285 ! mine step. ka .lt. 0 on input means this is the first
2286 ! attempt to determine step (for the present dig and dihdi)
2287 ! -- ka is initialized to 0 in this case. output with
2288 ! ka = 0 (or v(stppar) = 0) means step = -(h**-1)*g.
2289 ! l (i/o) = workspace of length p*(p+1)/2 for cholesky factors.
2290 ! p (in) = number of parameters -- the hessian is a p x p matrix.
2291 ! step (i/o) = the step computed.
2292 ! v (i/o) contains various constants and variables described below.
2293 ! w (i/o) = workspace of length 4*p + 6.
2295 ! *** entries in v ***
2297 ! v(dgnorm) (i/o) = 2-norm of (d**-1)*g.
2298 ! v(dstnrm) (output) = 2-norm of d*step.
2299 ! v(dst0) (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or
2300 ! overestimate of smallest eigenvalue of (d**-1)*h*(d**-1).
2301 ! v(epslon) (in) = max. rel. error allowed for psi(step). for the
2302 ! step returned, psi(step) will exceed its optimal value
2303 ! by less than -v(epslon)*psi(step). suggested value = 0.1.
2304 ! v(gtstep) (out) = inner product between g and step.
2305 ! v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step) (for pos. def.
2306 ! h only -- v(nreduc) is set to zero otherwise).
2307 ! v(phmnfc) (in) = tol. (together with v(phmxfc)) for accepting step
2308 ! (more*s sigma). the error v(dstnrm) - v(radius) must lie
2309 ! between v(phmnfc)*v(radius) and v(phmxfc)*v(radius).
2310 ! v(phmxfc) (in) (see v(phmnfc).)
2311 ! suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5.
2312 ! v(preduc) (out) = psi(step) = predicted obj. func. reduction for step.
2313 ! v(radius) (in) = radius of current (scaled) trust region.
2314 ! v(rad0) (i/o) = value of v(radius) from previous call.
2315 ! v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha
2316 ! described below under algorithm notes. if h + alpha*d**2
2317 ! (see algorithm notes) is (nearly) singular, however,
2318 ! then v(stppar) = -alpha.
2320 ! *** usage notes ***
2322 ! if it is desired to recompute step using a different value of
2323 ! v(radius), then this routine may be restarted by calling it
2324 ! with all parameters unchanged except v(radius). (this explains
2325 ! why step and w are listed as i/o). on an initial call (one with
2326 ! ka .lt. 0), step and w need not be initialized and only compo-
2327 ! nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and
2328 ! v(rad0) of v must be initialized.
2330 ! *** algorithm notes ***
2332 ! the desired g-q-t step (ref. 2, 3, 4, 6) satisfies
2333 ! (h + alpha*d**2)*step = -g for some nonnegative alpha such that
2334 ! h + alpha*d**2 is positive semidefinite. alpha and step are
2335 ! computed by a scheme analogous to the one described in ref. 5.
2336 ! estimates of the smallest and largest eigenvalues of the hessian
2337 ! are obtained from the gerschgorin circle theorem enhanced by a
2338 ! simple form of the scaling described in ref. 7. cases in which
2339 ! h + alpha*d**2 is nearly (or exactly) singular are handled by
2340 ! the technique discussed in ref. 2. in these cases, a step of
2341 ! (exact) length v(radius) is returned for which psi(step) exceeds
2342 ! its optimal value by less than -v(epslon)*psi(step). the test
2343 ! suggested in ref. 6 for detecting the special case is performed
2344 ! once two matrix factorizations have been done -- doing so sooner
2345 ! seems to degrade the performance of optimization routines that
2346 ! call this routine.
2348 ! *** functions and subroutines called ***
2350 ! dotprd - returns inner product of two vectors.
2351 ! litvmu - applies inverse-transpose of compact lower triang. matrix.
2352 ! livmul - applies inverse of compact lower triang. matrix.
2353 ! lsqrt - finds cholesky factor (of compactly stored lower triang.).
2354 ! lsvmin - returns approx. to min. sing. value of lower triang. matrix.
2355 ! rmdcon - returns machine-dependent constants.
2356 ! v2norm - returns 2-norm of a vector.
2358 ! *** references ***
2360 ! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive
2361 ! nonlinear least-squares algorithm, acm trans. math.
2362 ! software, vol. 7, no. 3.
2363 ! 2. gay, d.m. (1981), computing optimal locally constrained steps,
2364 ! siam j. sci. statist. computing, vol. 2, no. 2, pp.
2366 ! 3. goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966),
2367 ! maximization by quadratic hill-climbing, econometrica 34,
2369 ! 4. hebden, m.d. (1973), an algorithm for minimization using exact
2370 ! second derivatives, report t.p. 515, theoretical physics
2371 ! div., a.e.r.e. harwell, oxon., england.
2372 ! 5. more, j.j. (1978), the levenberg-marquardt algorithm, implemen-
2373 ! tation and theory, pp.105-116 of springer lecture notes
2374 ! in mathematics no. 630, edited by g.a. watson, springer-
2375 ! verlag, berlin and new york.
2376 ! 6. more, j.j., and sorensen, d.c. (1981), computing a trust region
2377 ! step, technical report anl-81-83, argonne national lab.
2378 ! 7. varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15,
2383 ! coded by david m. gay.
2384 ! this subroutine was written in connection with research
2385 ! supported by the national science foundation under grants
2386 ! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
2389 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2391 ! *** local variables ***
2394 integer :: dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,&
2395 j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x
2396 real(kind=8) :: alphak, aki, akk, delta, dst, eps, gtsta, lk,&
2397 oldphi, phi, phimax, phimin, psifac, rad, radsq,&
2398 root, si, sk, sw, t, twopsi, t1, t2, uk, wi
2401 real(kind=8) :: big, dgxfac !el, epsfac, four, half, kappa, negone,
2402 !el 1 one, p001, six, three, two, zero
2404 ! *** intrinsic functions ***
2406 !el real(kind=8) :: dabs, dmax1, dmin1, dsqrt
2408 ! *** external functions and subroutines ***
2410 !el external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm
2411 !el real(kind=8) :: dotprd, lsvmin, rmdcon, v2norm
2413 ! *** subscripts for v ***
2415 !el integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc,
2416 !el 1 phmnfc, phmxfc, preduc, radius, rad0
2418 ! data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/,
2419 ! 1 nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/,
2420 ! 2 rad0/9/, stppar/5/
2422 integer,parameter :: dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,&
2423 nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,&
2428 ! data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/,
2429 ! 1 kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/,
2430 ! 2 six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/
2432 real(kind=8), parameter :: epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,&
2433 kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,&
2434 six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0
2437 data big/0.d+0/, dgxfac/0.d+0/
2441 ! *** store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx).
2443 ! *** store gerschgorin over- and underestimates of the largest
2444 ! *** and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax)
2445 ! *** and w(emin) respectively.
2448 ! *** for use in recomputing step, the final values of lk, uk, dst,
2449 ! *** and the inverse derivative of more*s phi at 0 (for pos. def.
2450 ! *** h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin)
2456 ! *** store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p).
2459 ! *** store -d*step in w(q),...,w(q0+p).
2462 ! *** allocate storage for scratch vector x ***
2466 ! *** phitol = max. error allowed in dst = v(dstnrm) = 2-norm of
2468 phimax = v(phmxfc) * rad
2469 phimin = v(phmnfc) * rad
2470 psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) * &
2471 (kappa + one) + kappa + two) * rad**2)
2472 ! *** oldphi is used to detect limits of numerical accuracy. if
2473 ! *** we recompute step and it does not change, then we accept it.
2480 ! *** start or restart, depending on ka ***
2482 if (ka .ge. 0) go to 290
2484 ! *** fresh start ***
2490 v(dgnorm) = v2norm(p, dig)
2494 if (v(dgnorm) .eq. zero) kamin = 0
2496 ! *** store diag(dihdi) in w(diag0+1),...,w(diag0+p) ***
2505 ! *** determine w(dggdmx), the largest element of dihdi ***
2511 if (t1 .lt. t) t1 = t
2515 ! *** try alpha = 0 ***
2517 30 call lsqrt(1, p, l, dihdi, irc)
2518 if (irc .eq. 0) go to 50
2519 ! *** indef. h -- underestimate smallest eigenvalue, use this
2520 ! *** estimate to initialize lower bound lk on alpha.
2527 call litvmu(irc, w, l, w)
2531 if (restrt) go to 210
2534 ! *** positive definite h -- compute unmodified newton step. ***
2536 t = lsvmin(p, l, w(q), w(q))
2537 if (t .ge. one) go to 60
2538 if (big .le. zero) big = rmdcon(6)
2539 if (v(dgnorm) .ge. t*t*big) go to 70
2540 60 call livmul(p, w(q), l, dig)
2541 gtsta = dotprd(p, w(q), w(q))
2542 v(nreduc) = half * gtsta
2543 call litvmu(p, w(q), l, w(q))
2544 dst = v2norm(p, w(q))
2547 if (phi .le. phimax) go to 260
2548 if (restrt) go to 210
2550 ! *** prepare to compute gerschgorin estimates of largest (and
2551 ! *** smallest) eigenvalues. ***
2556 if (i .eq. 1) go to 90
2568 ! *** (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1) ***
2572 if (p .le. 1) go to 120
2576 if (t .ge. t1) go to 110
2588 if (i .eq. k) go to 130
2589 aki = dabs(dihdi(k1))
2592 t1 = half * (akk - w(j) + si - aki)
2593 t1 = t1 + dsqrt(t1*t1 + sk*aki)
2594 if (t .lt. t1) t = t1
2595 if (i .lt. k) go to 140
2601 uk = v(dgnorm)/rad - w(emin)
2602 if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2603 if (uk .le. zero) uk = p001
2605 ! *** compute gerschgorin (over-)estimate of largest eigenvalue ***
2609 if (p .le. 1) go to 170
2613 if (t .le. t1) go to 160
2625 if (i .eq. k) go to 180
2626 aki = dabs(dihdi(k1))
2629 t1 = half * (w(j) + si - aki - akk)
2630 t1 = t1 + dsqrt(t1*t1 + sk*aki)
2631 if (t .lt. t1) t = t1
2632 if (i .lt. k) go to 190
2638 lk = dmax1(lk, v(dgnorm)/rad - w(emax))
2640 ! *** alphak = current value of alpha (see alg. notes above). we
2641 ! *** use more*s scheme for initializing it.
2642 alphak = dabs(v(stppar)) * v(rad0)/rad
2644 if (irc .ne. 0) go to 210
2646 ! *** compute l0 for positive definite h ***
2648 call livmul(p, w, l, w(q))
2650 w(phipin) = dst / t / t
2651 lk = dmax1(lk, phi*w(phipin))
2653 ! *** safeguard alphak and add alphak*i to (d**-1)*h*(d**-1) ***
2656 if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk) &
2657 alphak = uk * dmax1(p001, dsqrt(lk/uk))
2658 if (alphak .le. zero) alphak = half * uk
2659 if (alphak .le. zero) alphak = uk
2664 dihdi(k) = w(j) + alphak
2667 ! *** try computing cholesky decomposition ***
2669 call lsqrt(1, p, l, dihdi, irc)
2670 if (irc .eq. 0) go to 240
2672 ! *** (d**-1)*h*(d**-1) + alphak*i is indefinite -- overestimate
2673 ! *** smallest eigenvalue for use in updating lk ***
2681 call litvmu(irc, w, l, w)
2683 lk = alphak - t/t1/t1
2687 ! *** alphak makes (d**-1)*h*(d**-1) positive definite.
2688 ! *** compute q = -d*step, check for convergence. ***
2690 240 call livmul(p, w(q), l, dig)
2691 gtsta = dotprd(p, w(q), w(q))
2692 call litvmu(p, w(q), l, w(q))
2693 dst = v2norm(p, w(q))
2695 if (phi .le. phimax .and. phi .ge. phimin) go to 270
2696 if (phi .eq. oldphi) go to 270
2698 if (phi .lt. zero) go to 330
2700 ! *** unacceptable alphak -- update lk, uk, alphak ***
2702 250 if (ka .ge. kalim) go to 270
2703 ! *** the following dmin1 is necessary because of restarts ***
2704 if (phi .lt. zero) uk = dmin1(uk, alphak)
2705 ! *** kamin = 0 only iff the gradient vanishes ***
2706 if (kamin .eq. 0) go to 210
2707 call livmul(p, w, l, w(q))
2709 alphak = alphak + (phi/t1) * (dst/t1) * (dst/rad)
2710 lk = dmax1(lk, alphak)
2713 ! *** acceptable step on first try ***
2717 ! *** successful step in general. compute step = -(d**-1)*q ***
2721 step(i) = -w(j)/d(i)
2724 v(preduc) = half * (dabs(alphak)*dst*dst + gtsta)
2728 ! *** restart with new radius ***
2730 290 if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310
2732 ! *** prepare to return newton step ***
2746 if (v(dgnorm) .eq. zero) kamin = 0
2747 if (ka .eq. 0) go to 50
2750 alphak = dabs(v(stppar))
2754 if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2755 if (uk .le. zero) uk = p001
2756 if (rad .gt. v(rad0)) go to 320
2758 ! *** smaller radius ***
2760 if (alphak .gt. zero) lk = w(lk0)
2761 lk = dmax1(lk, t - w(emax))
2762 if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2765 ! *** bigger radius ***
2766 320 if (alphak .gt. zero) uk = dmin1(uk, w(uk0))
2767 lk = dmax1(zero, -v(dst0), t - w(emax))
2768 if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2771 ! *** decide whether to check for special case... in practice (from
2772 ! *** the standpoint of the calling optimization code) it seems best
2773 ! *** not to check until a few iterations have failed -- hence the
2774 ! *** test on kamin below.
2776 330 delta = alphak + dmin1(zero, v(dst0))
2777 twopsi = alphak*dst*dst + gtsta
2778 if (ka .ge. kamin) go to 340
2779 ! *** if the test in ref. 2 is satisfied, fall through to handle
2780 ! *** the special case (as soon as the more-sorensen test detects
2782 if (delta .ge. psifac*twopsi) go to 370
2784 ! *** check for the special case of h + alpha*d**2 (nearly)
2785 ! *** singular. use one step of inverse power method with start
2786 ! *** from lsvmin to obtain approximate eigenvector corresponding
2787 ! *** to smallest eigenvalue of (d**-1)*h*(d**-1). lsvmin returns
2788 ! *** x and w with l*w = x.
2790 340 t = lsvmin(p, l, w(x), w)
2792 ! *** normalize w ***
2795 ! *** complete current inv. power iter. -- replace w by (l**-t)*w.
2796 call litvmu(p, w, l, w)
2797 t2 = one/v2norm(p, w)
2802 ! *** now w is the desired approximate (unit) eigenvector and
2803 ! *** t*x = ((d**-1)*h*(d**-1) + alphak*i)*w.
2805 sw = dotprd(p, w(q), w)
2806 t1 = (rad + dst) * (rad - dst)
2807 root = dsqrt(sw*sw + t1)
2808 if (sw .lt. zero) root = -root
2809 si = t1 / (sw + root)
2811 ! *** the actual test for the special case...
2813 if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380
2815 ! *** update upper bound on smallest eigenvalue (when not positive)
2816 ! *** (as recommended by more and sorensen) and continue...
2818 if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak)
2819 lk = dmax1(lk, -v(dst0))
2821 ! *** check whether we can hope to detect the special case in
2822 ! *** the available arithmetic. accept step as it is if not.
2824 ! *** if not yet available, obtain machine dependent value dgxfac.
2825 370 if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3)
2827 if (delta .gt. dgxfac*w(dggdmx)) go to 250
2830 ! *** special case detected... negate alphak to indicate special case
2832 380 alphak = -alphak
2833 v(preduc) = half * twopsi
2835 ! *** accept current step if adding si*w would lead to a
2836 ! *** further relative reduction in psi of less than v(epslon)/3.
2839 t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w)))
2840 if (t .lt. eps*twopsi/six) go to 390
2841 v(preduc) = v(preduc) + t
2846 w(j) = t1*w(i) - w(j)
2847 step(i) = w(j) / d(i)
2849 v(gtstep) = dotprd(p, dig, w(q))
2851 ! *** save values for use in a possible restart ***
2860 ! *** restore diagonal of dihdi ***
2871 ! *** last card of gqtst follows ***
2872 end subroutine gqtst
2873 !-----------------------------------------------------------------------------
2874 subroutine lsqrt(n1, n, l, a, irc)
2876 ! *** compute rows n1 through n of the cholesky factor l of
2877 ! *** a = l*(l**t), where l and the lower triangle of a are both
2878 ! *** stored compactly by rows (and may occupy the same storage).
2879 ! *** irc = 0 means all went well. irc = j means the leading
2880 ! *** principal j x j submatrix of a is not positive definite --
2881 ! *** and l(j*(j+1)/2) contains the (nonpos.) reduced j-th diagonal.
2883 ! *** parameters ***
2885 integer :: n1, n, irc
2886 !al real(kind=8) :: l(1), a(1)
2887 real(kind=8) :: l(n*(n+1)/2), a(n*(n+1)/2)
2888 ! dimension l(n*(n+1)/2), a(n*(n+1)/2)
2890 ! *** local variables ***
2892 integer :: i, ij, ik, im1, i0, j, jk, jm1, j0, k
2893 real(kind=8) :: t, td !el, zero
2895 ! *** intrinsic functions ***
2897 !el real(kind=8) :: dsqrt
2902 real(kind=8),parameter :: zero=0.d+0
2907 i0 = n1 * (n1 - 1) / 2
2910 if (i .eq. 1) go to 40
2915 if (j .eq. 1) go to 20
2924 t = (a(ij) - t) / l(j0)
2930 if (t .le. zero) go to 60
2942 ! *** last card of lsqrt ***
2943 end subroutine lsqrt
2944 !-----------------------------------------------------------------------------
2945 real(kind=8) function lsvmin(p, l, x, y)
2947 ! *** estimate smallest sing. value of packed lower triang. matrix l
2949 ! *** parameter declarations ***
2952 !al real(kind=8) :: l(1), x(p), y(p)
2953 real(kind=8) :: l(p*(p+1)/2), x(p), y(p)
2954 ! dimension l(p*(p+1)/2)
2956 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2960 ! this function returns a good over-estimate of the smallest
2961 ! singular value of the packed lower triangular matrix l.
2963 ! *** parameter description ***
2965 ! p (in) = the order of l. l is a p x p lower triangular matrix.
2966 ! l (in) = array holding the elements of l in row order, i.e.
2967 ! l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc.
2968 ! x (out) if lsvmin returns a positive value, then x is a normalized
2969 ! approximate left singular vector corresponding to the
2970 ! smallest singular value. this approximation may be very
2971 ! crude. if lsvmin returns zero, then some components of x
2972 ! are zero and the rest retain their input values.
2973 ! y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an
2974 ! unnormalized approximate right singular vector correspond-
2975 ! ing to the smallest singular value. this approximation
2976 ! may be crude. if lsvmin returns zero, then y retains its
2977 ! input value. the caller may pass the same vector for x
2978 ! and y (nonstandard fortran usage), in which case y over-
2979 ! writes x (for nonzero lsvmin returns).
2981 ! *** algorithm notes ***
2983 ! the algorithm is based on (1), with the additional provision that
2984 ! lsvmin = 0 is returned if the smallest diagonal element of l
2985 ! (in magnitude) is not more than the unit roundoff times the
2986 ! largest. the algorithm uses a random number generator proposed
2987 ! in (4), which passes the spectral test with flying colors -- see
2990 ! *** subroutines and functions called ***
2992 ! v2norm - function, returns the 2-norm of a vector.
2994 ! *** references ***
2996 ! (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977),
2997 ! an estimate for the condition number of a matrix, report
2998 ! tm-310, applied math. div., argonne national laboratory.
3000 ! (2) hoaglin, d.c. (1976), theoretical properties of congruential
3001 ! random-number generators -- an empirical view,
3002 ! memorandum ns-340, dept. of statistics, harvard univ.
3004 ! (3) knuth, d.e. (1969), the art of computer programming, vol. 2
3005 ! (seminumerical algorithms), addison-wesley, reading, mass.
3007 ! (4) smith, c.s. (1971), multiplicative pseudo-random number
3008 ! generators with prime modulus, j. assoc. comput. mach. 18,
3013 ! designed and coded by david m. gay (winter 1977/summer 1978).
3017 ! this subroutine was written in connection with research
3018 ! supported by the national science foundation under grants
3019 ! mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989.
3021 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3023 ! *** local variables ***
3025 integer :: i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1
3026 real(kind=8) :: b, sminus, splus, t, xminus, xplus
3030 !el real(kind=8) :: half, one, r9973, zero
3032 ! *** intrinsic functions ***
3036 !el real(kind=8) :: dabs
3038 ! *** external functions and subroutines ***
3040 !el external dotprd, v2norm, vaxpy
3041 !el real(kind=8) :: dotprd, v2norm
3044 ! data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/
3046 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0
3054 ! *** first check whether to return lsvmin = 0 and initialize x ***
3059 if (l(jj) .eq. zero) go to 110
3060 ix = mod(3432*ix, 9973)
3061 b = half*(one + float(ix)/r9973)
3064 if (p .le. 1) go to 60
3067 if (l(ii) .eq. zero) go to 110
3069 x(i) = xplus * l(ji)
3072 ! *** solve (l**t)*x = b, where the components of b have randomly
3073 ! *** chosen magnitudes in (.5,1) with signs chosen to make x large.
3075 ! do j = p-1 to 1 by -1...
3078 ! *** determine x(j) in this iteration. note for i = 1,2,...,j
3079 ! *** that x(i) holds the current partial sum for row i.
3080 ix = mod(3432*ix, 9973)
3081 b = half*(one + float(ix)/r9973)
3083 xminus = (-b - x(j))
3085 sminus = dabs(xminus)
3090 xminus = xminus/l(jj)
3091 if (jm1 .eq. 0) go to 30
3094 splus = splus + dabs(x(i) + l(ji)*xplus)
3095 sminus = sminus + dabs(x(i) + l(ji)*xminus)
3097 30 if (sminus .gt. splus) xplus = xminus
3099 ! *** update partial sums ***
3100 if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x)
3103 ! *** normalize x ***
3105 60 t = one/v2norm(p, x)
3109 ! *** solve l*y = x and return lsvmin = 1/twonorm(y) ***
3116 if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y)
3117 y(j) = (x(j) - t) / l(jj)
3120 lsvmin = one/v2norm(p, y)
3125 ! *** last card of lsvmin follows ***
3127 !-----------------------------------------------------------------------------
3128 subroutine slvmul(p, y, s, x)
3130 ! *** set y = s * x, s = p x p symmetric matrix. ***
3131 ! *** lower triangle of s stored rowwise. ***
3133 ! *** parameter declarations ***
3136 !al real(kind=8) :: s(1), x(p), y(p)
3137 real(kind=8) :: s(p*(p+1)/2), x(p), y(p)
3138 ! dimension s(p*(p+1)/2)
3140 ! *** local variables ***
3142 integer :: i, im1, j, k
3145 ! *** no intrinsic functions ***
3147 ! *** external function ***
3150 !el real(kind=8) :: dotprd
3152 !-----------------------------------------------------------------------
3156 y(i) = dotprd(i, s(j), x)
3160 if (p .le. 1) go to 999
3167 y(k) = y(k) + s(j)*xi
3173 ! *** last card of slvmul follows ***
3174 end subroutine slvmul
3175 !-----------------------------------------------------------------------------
3177 !-----------------------------------------------------------------------------
3178 subroutine minimize(etot,x,iretcode,nfun)
3180 use energy, only: func,gradient,fdum!,etotal,enerprint
3182 ! implicit real*8 (a-h,o-z)
3183 ! include 'DIMENSIONS'
3184 integer,parameter :: liv=60
3185 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3186 !********************************************************************
3187 ! OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
3188 ! the calling subprogram. *
3189 ! when d(i)=1.0, then v(35) is the length of the initial step, *
3190 ! calculated in the usual pythagorean way. *
3191 ! absolute convergence occurs when the function is within v(31) of *
3192 ! zero. unless you know the minimum value in advance, abs convg *
3193 ! is probably not useful. *
3194 ! relative convergence is when the model predicts that the function *
3195 ! will decrease by less than v(32)*abs(fun). *
3196 !********************************************************************
3197 ! include 'COMMON.IOUNITS'
3198 ! include 'COMMON.VAR'
3199 ! include 'COMMON.GEO'
3200 ! include 'COMMON.MINIM'
3202 !el common /srutu/ icall
3203 integer,dimension(liv) :: iv
3204 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
3205 !el real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3206 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3207 real(kind=8) :: energia(0:n_ene)
3208 ! external func,gradient,fdum
3209 ! external func_restr,grad_restr
3210 logical :: not_done,change,reduce
3211 !el common /przechowalnia/ v
3213 integer :: iretcode,nfun,lv,nvar_restr,idum(1),j
3214 real(kind=8) :: etot,rdum(1) !,fdum
3216 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3218 if (.not.allocated(v)) allocate(v(1:lv))
3224 ! DO WHILE (NOT_DONE)
3226 call deflt(2,iv,liv,lv,v)
3227 ! 12 means fresh start, dont call deflt
3229 ! max num of fun calls
3230 if (maxfun.eq.0) maxfun=500
3232 ! max num of iterations
3233 if (maxmin.eq.0) maxmin=1000
3237 ! selects output unit
3239 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
3240 ! 1 means to print out result
3241 iv(22)=print_min_res
3242 ! 1 means to print out summary stats
3243 iv(23)=print_min_stat
3244 ! 1 means to print initial x and d
3245 iv(24)=print_min_ini
3246 ! min val for v(radfac) default is 0.1
3248 ! max val for v(radfac) default is 4.0
3251 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3252 ! the sumsl default is 0.1
3254 ! false conv if (act fnctn decrease) .lt. v(34)
3255 ! the sumsl default is 100*machep
3257 ! absolute convergence
3258 if (tolf.eq.0.0D0) tolf=1.0D-4
3260 ! relative convergence
3261 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3263 ! controls initial step size
3265 ! large vals of d correspond to small components of step
3272 !d print *,'Calling SUMSL'
3273 ! call var_to_geom(nvar,x)
3275 ! call etotal(energia(0))
3279 call x2xx(x,xx,nvar_restr)
3280 call sumsl(nvar_restr,d,xx,func_restr,grad_restr,&
3281 iv,liv,lv,v,idum,rdum,fdum)
3284 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
3288 !d print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
3289 !d write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
3292 call var_to_geom(nvar,x)
3294 ! write (iout,'(a)') 'Reduction worked, minimizing again...'
3298 write(iout,*) 'Warning calling chainbuild'
3301 !el---------------------
3302 ! write (iout,'(/a)') &
3303 ! "Cartesian coordinates of the reference structure after SUMSL"
3304 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3305 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3307 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
3308 ! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
3309 ! (c(j,i+nres),j=1,3)
3311 !el----------------------------
3312 ! call etotal(energia) !sp
3314 ! call enerprint(energia) !sp
3317 ! write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
3322 end subroutine minimize
3323 !-----------------------------------------------------------------------------
3325 !-----------------------------------------------------------------------------
3326 subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
3328 use energy, only: cartder,zerograd,etotal,sum_gradient
3329 ! implicit real*8 (a-h,o-z)
3330 ! include 'DIMENSIONS'
3331 ! include 'COMMON.CHAIN'
3332 ! include 'COMMON.DERIV'
3333 ! include 'COMMON.VAR'
3334 ! include 'COMMON.INTERACT'
3335 ! include 'COMMON.FFIELD'
3336 ! include 'COMMON.IOUNITS'
3338 integer :: uiparm(1)
3339 real(kind=8) :: urparm(1)
3340 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
3341 integer :: n,nf,ig,ind,i,j,ij,k,igall
3342 real(kind=8) :: f,gphii,gthetai,galphai,gomegai
3343 real(kind=8),external :: ufparm
3346 if (nf-nfl+1) 20,30,40
3347 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
3348 ! write (iout,*) 'grad 20'
3353 ! Intercept NaNs in the coordinates
3354 ! write(iout,*) (var(i),i=1,nvar)
3359 if (x_sum.ne.x_sum) then
3360 write(iout,*)" *** grad_restr : Found NaN in coordinates"
3362 print *," *** grad_restr : Found NaN in coordinates"
3366 call var_to_geom_restr(n,x)
3367 write(iout,*) 'Warning calling chainbuild'
3370 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
3374 ! Convert the Cartesian gradient into internal-coordinate gradient.
3380 IF (mask_phi(i+2).eq.1) THEN
3385 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
3386 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
3399 IF (mask_theta(i+2).eq.1) THEN
3405 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
3406 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
3416 if (itype(i,1).ne.10) then
3417 IF (mask_side(i).eq.1) THEN
3421 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
3430 if (itype(i,1).ne.10) then
3431 IF (mask_side(i).eq.1) THEN
3435 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
3443 ! Add the components corresponding to local energy terms.
3450 if (mask_phi(i).eq.1) then
3452 g(ig)=g(ig)+gloc(igall,icg)
3458 if (mask_theta(i).eq.1) then
3460 g(ig)=g(ig)+gloc(igall,icg)
3466 if (itype(i,1).ne.10) then
3468 if (mask_side(i).eq.1) then
3470 g(ig)=g(ig)+gloc(igall,icg)
3477 !d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
3480 end subroutine grad_restr
3481 !-----------------------------------------------------------------------------
3482 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F
3485 use energy, only: zerograd,etotal,sum_gradient
3486 ! implicit real*8 (a-h,o-z)
3487 ! include 'DIMENSIONS'
3488 ! include 'COMMON.DERIV'
3489 ! include 'COMMON.IOUNITS'
3490 ! include 'COMMON.GEO'
3493 !el common /chuju/ jjj
3494 real(kind=8) :: energia(0:n_ene)
3496 real(kind=8),external :: ufparm
3497 integer :: uiparm(1)
3498 real(kind=8) :: urparm(1)
3499 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3500 ! if (jjj.gt.0) then
3501 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
3505 call var_to_geom_restr(n,x)
3507 write(iout,*) 'Warning calling chainbuild'
3509 !d write (iout,*) 'ETOTAL called from FUNC'
3510 call etotal(energia)
3513 ! if (jjj.gt.0) then
3514 ! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
3515 ! write (iout,*) 'f=',etot
3519 end subroutine func_restr
3520 !-----------------------------------------------------------------------------
3521 ! subroutine func(n,x,nf,f,uiparm,urparm,ufparm) in module energy
3522 !-----------------------------------------------------------------------------
3523 subroutine x2xx(x,xx,n)
3525 ! implicit real*8 (a-h,o-z)
3526 ! include 'DIMENSIONS'
3527 ! include 'COMMON.VAR'
3528 ! include 'COMMON.CHAIN'
3529 ! include 'COMMON.INTERACT'
3530 integer :: n,i,ij,ig,igall
3531 real(kind=8),dimension(6*nres) :: xx,x !(maxvar) (maxvar=6*maxres)
3533 !el allocate(varall(nvar)) allocated in alioc_ener_arrays
3543 if (mask_phi(i).eq.1) then
3551 if (mask_theta(i).eq.1) then
3559 if (itype(i,1).ne.10) then
3561 if (mask_side(i).eq.1) then
3573 !-----------------------------------------------------------------------------
3574 !el subroutine xx2x(x,xx) in module math
3575 !-----------------------------------------------------------------------------
3576 subroutine minim_dc(etot,iretcode,nfun)
3579 use energy, only: fdum,check_ecartint
3580 ! implicit real*8 (a-h,o-z)
3581 ! include 'DIMENSIONS'
3585 integer,parameter :: liv=60
3586 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3587 ! include 'COMMON.SETUP'
3588 ! include 'COMMON.IOUNITS'
3589 ! include 'COMMON.VAR'
3590 ! include 'COMMON.GEO'
3591 ! include 'COMMON.MINIM'
3592 ! include 'COMMON.CHAIN'
3593 integer :: iretcode,nfun,k,i,j,lv,idum(1)
3594 integer,dimension(liv) :: iv
3595 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
3596 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
3597 !el common /przechowalnia/ v
3599 real(kind=8) :: energia(0:n_ene)
3600 ! external func_dc,grad_dc ,fdum
3601 logical :: not_done,change,reduce
3602 real(kind=8) :: g(6*nres),f1,etot,rdum(1) !,fdum
3604 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3606 if (.not. allocated(v)) allocate(v(1:lv))
3608 call deflt(2,iv,liv,lv,v)
3609 ! 12 means fresh start, dont call deflt
3611 ! max num of fun calls
3612 if (maxfun.eq.0) maxfun=500
3614 ! max num of iterations
3615 if (maxmin.eq.0) maxmin=1000
3619 ! selects output unit
3621 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
3622 ! 1 means to print out result
3623 iv(22)=print_min_res
3624 ! 1 means to print out summary stats
3625 iv(23)=print_min_stat
3626 ! 1 means to print initial x and d
3627 iv(24)=print_min_ini
3628 ! min val for v(radfac) default is 0.1
3630 ! max val for v(radfac) default is 4.0
3633 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3634 ! the sumsl default is 0.1
3636 ! false conv if (act fnctn decrease) .lt. v(34)
3637 ! the sumsl default is 100*machep
3639 ! absolute convergence
3640 if (tolf.eq.0.0D0) tolf=1.0D-4
3642 ! relative convergence
3643 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3645 ! controls initial step size
3647 ! large vals of d correspond to small components of step
3660 if (ialph(i,1).gt.0) then
3668 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
3678 if (ialph(i,1).gt.0) then
3685 call chainbuild_cart
3689 !d call func_dc(k,x,nf,f,idum,rdum,fdum)
3690 !d call grad_dc(k,x,nf,g,idum,rdum,fdum)
3694 !d call func_dc(k,x,nf,f1,idum,rdum,fdum)
3696 !d print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
3698 !el---------------------
3699 ! write (iout,'(/a)') &
3700 ! "Cartesian coordinates of the reference structure after SUMSL"
3701 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3702 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3704 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
3705 ! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
3706 ! (c(j,i+nres),j=1,3)
3708 !el----------------------------
3713 end subroutine minim_dc
3714 !-----------------------------------------------------------------------------
3715 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
3718 use energy, only: zerograd,etotal
3719 ! implicit real*8 (a-h,o-z)
3720 ! include 'DIMENSIONS'
3724 ! include 'COMMON.SETUP'
3725 ! include 'COMMON.DERIV'
3726 ! include 'COMMON.IOUNITS'
3727 ! include 'COMMON.GEO'
3728 ! include 'COMMON.CHAIN'
3729 ! include 'COMMON.VAR'
3730 integer :: n,nf,k,i,j
3731 real(kind=8) :: energia(0:n_ene)
3732 real(kind=8),external :: ufparm
3733 integer :: uiparm(1)
3734 real(kind=8) :: urparm(1)
3735 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3738 !bad icg=mod(nf,2)+1
3749 if (ialph(i,1).gt.0) then
3756 call chainbuild_cart
3759 call etotal(energia)
3762 !d print *,'func_dc ',nf,nfl,f
3765 end subroutine func_dc
3766 !-----------------------------------------------------------------------------
3767 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
3770 use energy, only: cartgrad,zerograd,etotal
3772 ! implicit real*8 (a-h,o-z)
3773 ! include 'DIMENSIONS'
3777 ! include 'COMMON.SETUP'
3778 ! include 'COMMON.CHAIN'
3779 ! include 'COMMON.DERIV'
3780 ! include 'COMMON.VAR'
3781 ! include 'COMMON.INTERACT'
3782 ! include 'COMMON.FFIELD'
3783 ! include 'COMMON.MD'
3784 ! include 'COMMON.IOUNITS'
3785 real(kind=8),external :: ufparm
3786 integer :: n,nf,i,j,k
3787 integer :: uiparm(1)
3788 real(kind=8) :: urparm(1)
3789 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
3792 !elwrite(iout,*) "jestesmy w grad dc"
3794 !bad icg=mod(nf,2)+1
3796 !d print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
3797 !elwrite(iout,*) "jestesmy w grad dc nf-nfl+1", nf-nfl+1
3798 if (nf-nfl+1) 20,30,40
3799 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
3813 if (ialph(i,1).gt.0) then
3820 !elwrite(iout,*) "jestesmy w grad dc"
3821 call chainbuild_cart
3824 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
3828 !elwrite(iout,*) "jestesmy w grad dc"
3837 if (ialph(i,1).gt.0) then
3844 !elwrite(iout,*) "jestesmy w grad dc"
3847 end subroutine grad_dc
3848 !-----------------------------------------------------------------------------
3850 !-----------------------------------------------------------------------------
3852 subroutine minim_mcmf
3856 use energy, only: func,gradient,fdum
3857 ! implicit real*8 (a-h,o-z)
3858 ! include 'DIMENSIONS'
3860 integer,parameter :: liv=60
3861 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3862 ! include 'COMMON.VAR'
3863 ! include 'COMMON.IOUNITS'
3864 ! include 'COMMON.MINIM'
3865 ! real(kind=8) :: fdum
3866 ! external func,gradient,fdum
3867 !el real(kind=4) :: ran1,ran2,ran3
3868 ! include 'COMMON.SETUP'
3869 ! include 'COMMON.GEO'
3870 ! include 'COMMON.CHAIN'
3871 ! include 'COMMON.FFIELD'
3872 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3873 real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg
3874 real(kind=8),dimension(6*nres) :: d,garbage !(maxvar) (maxvar=6*maxres)
3875 !el real(kind=8) :: v(1:77+(6*nres)*(6*nres+17)/2+1)
3876 integer,dimension(6) :: indx
3877 integer,dimension(liv) :: iv
3878 integer :: lv,idum(1),nf !
3879 real(kind=8) :: rdum(1)
3880 real(kind=8) :: przes(3),obrot(3,3),eee
3883 integer,dimension(MPI_STATUS_SIZE) :: muster
3885 integer :: ichuj,i,ierr
3886 real(kind=8) :: rad,ene0
3887 data rad /1.745329252d-2/
3888 !el common /przechowalnia/ v
3890 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
3891 if (.not. allocated(v)) allocate(v(1:lv))
3896 call mpi_recv(indx,6,mpi_integer,king,idint,CG_COMM,&
3898 if (indx(1).eq.0) return
3899 ! print *, 'worker ',me,' received order ',n,ichuj
3900 call mpi_recv(var,nvar,mpi_double_precision,&
3901 king,idreal,CG_COMM,muster,ierr)
3902 call mpi_recv(ene0,1,mpi_double_precision,&
3903 king,idreal,CG_COMM,muster,ierr)
3904 ! print *, 'worker ',me,' var read '
3907 call deflt(2,iv,liv,lv,v)
3908 ! 12 means fresh start, dont call deflt
3910 ! max num of fun calls
3911 if (maxfun.eq.0) maxfun=500
3913 ! max num of iterations
3914 if (maxmin.eq.0) maxmin=1000
3918 ! selects output unit
3921 ! 1 means to print out result
3923 ! 1 means to print out summary stats
3925 ! 1 means to print initial x and d
3927 ! min val for v(radfac) default is 0.1
3929 ! max val for v(radfac) default is 4.0
3931 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
3932 ! the sumsl default is 0.1
3934 ! false conv if (act fnctn decrease) .lt. v(34)
3935 ! the sumsl default is 100*machep
3937 ! absolute convergence
3938 if (tolf.eq.0.0D0) tolf=1.0D-4
3940 ! relative convergence
3941 if (rtolf.eq.0.0D0) rtolf=1.0D-4
3943 ! controls initial step size
3945 ! large vals of d correspond to small components of step
3954 call func(nvar,var,nf,eee,idum,rdum,fdum)
3955 if(eee.gt.1.0d18) then
3956 ! print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
3957 ! print *,' energy before SUMSL =',eee
3958 ! print *,' aborting local minimization'
3965 call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
3966 ! find which conformation was returned from sumsl
3969 ! total # of ftn evaluations (for iwf=0, it includes all minimizations).
3974 call mpi_send(indx,6,mpi_integer,king,idint,CG_COMM,&
3976 ! print '(a5,i3,15f10.5)', 'ENEX0',indx(1),v(10)
3977 call mpi_send(var,nvar,mpi_double_precision,&
3978 king,idreal,CG_COMM,ierr)
3979 call mpi_send(eee,1,mpi_double_precision,king,idreal,&
3981 call mpi_send(ene0,1,mpi_double_precision,king,idreal,&
3985 end subroutine minim_mcmf
3987 !-----------------------------------------------------------------------------
3989 !-----------------------------------------------------------------------------
3990 ! algorithm 611, collected algorithms from acm.
3991 ! algorithm appeared in acm-trans. math. software, vol.9, no. 4,
3992 ! dec., 1983, p. 503-524.
3993 integer function imdcon(k)
3997 ! *** return integer machine-dependent constants ***
3999 ! *** k = 1 means return standard output unit number. ***
4000 ! *** k = 2 means return alternate output unit number. ***
4001 ! *** k = 3 means return input unit number. ***
4002 ! (note -- k = 2, 3 are used only by test programs.)
4004 ! +++ port version follows...
4008 ! data mdperm(1)/2/, mdperm(2)/4/, mdperm(3)/1/
4009 ! imdcon = i1mach(mdperm(k))
4010 ! +++ end of port version +++
4012 ! +++ non-port version follows...
4014 data mdcon(1)/6/, mdcon(2)/8/, mdcon(3)/5/
4016 ! +++ end of non-port version +++
4019 ! *** last card of imdcon follows ***
4021 !-----------------------------------------------------------------------------
4022 real(kind=8) function rmdcon(k)
4024 ! *** return machine dependent constants used by nl2sol ***
4026 ! +++ comments below contain data statements for various machines. +++
4027 ! +++ to convert to another machine, place a c in column 1 of the +++
4028 ! +++ data statement line(s) that correspond to the current machine +++
4029 ! +++ and remove the c from column 1 of the data statement line(s) +++
4030 ! +++ that correspond to the new machine. +++
4034 ! *** the constant returned depends on k...
4036 ! *** k = 1... smallest pos. eta such that -eta exists.
4037 ! *** k = 2... square root of eta.
4038 ! *** k = 3... unit roundoff = smallest pos. no. machep such
4039 ! *** that 1 + machep .gt. 1 .and. 1 - machep .lt. 1.
4040 ! *** k = 4... square root of machep.
4041 ! *** k = 5... square root of big (see k = 6).
4042 ! *** k = 6... largest machine no. big such that -big exists.
4044 real(kind=8) :: big, eta, machep
4045 integer :: bigi(4), etai(4), machei(4)
4047 !el real(kind=8) :: dsqrt
4049 equivalence (big,bigi(1)), (eta,etai(1)), (machep,machei(1))
4051 ! +++ ibm 360, ibm 370, or xerox +++
4053 ! data big/z7fffffffffffffff/, eta/z0010000000000000/,
4054 ! 1 machep/z3410000000000000/
4056 ! +++ data general +++
4058 ! data big/0.7237005577d+76/, eta/0.5397605347d-78/,
4059 ! 1 machep/2.22044605d-16/
4063 ! data big/1.7d+38/, eta/2.938735878d-39/, machep/2.775557562d-17/
4067 ! data big/1.157920892d+77/, eta/8.636168556d-78/,
4068 ! 1 machep/5.551115124d-17/
4072 ! data big/1.69d+38/, eta/5.9d-39/, machep/2.1680435d-19/
4076 ! data big/"377777100000000000000000/,
4077 ! 1 eta/"002400400000000000000000/,
4078 ! 2 machep/"104400000000000000000000/
4082 ! data big/o0777777777777777,o7777777777777777/,
4083 ! 1 eta/o1771000000000000,o7770000000000000/,
4084 ! 2 machep/o1451000000000000,o0000000000000000/
4086 ! +++ control data +++
4088 ! data big/37767777777777777777b,37167777777777777777b/,
4089 ! 1 eta/00014000000000000000b,00000000000000000000b/,
4090 ! 2 machep/15614000000000000000b,15010000000000000000b/
4094 ! data big/1.0d+9786/, eta/1.0d-9860/, machep/1.4210855d-14/
4098 ! data big/8.988d+307/, eta/1.2d-308/, machep/1.734723476d-18/
4102 data big/1.7d+38/, eta/2.939d-39/, machep/1.3877788d-17/
4106 ! data bigi(1)/577767777777777777777b/,
4107 ! 1 bigi(2)/000007777777777777776b/,
4108 ! 2 etai(1)/200004000000000000000b/,
4109 ! 3 etai(2)/000000000000000000000b/,
4110 ! 4 machei(1)/377224000000000000000b/,
4111 ! 5 machei(2)/000000000000000000000b/
4113 ! +++ port library -- requires more than just a data statement... +++
4116 ! double precision d1mach, zero
4117 ! data big/0.d+0/, eta/0.d+0/, machep/0.d+0/, zero/0.d+0/
4118 ! if (big .gt. zero) go to 1
4121 ! machep = d1mach(4)
4124 ! +++ end of port +++
4126 !------------------------------- body --------------------------------
4128 go to (10, 20, 30, 40, 50, 60), k
4133 20 rmdcon = dsqrt(256.d+0*eta)/16.d+0
4139 40 rmdcon = dsqrt(machep)
4142 50 rmdcon = dsqrt(big/256.d+0)*16.d+0
4148 ! *** last card of rmdcon follows ***
4150 !-----------------------------------------------------------------------------
4152 !-----------------------------------------------------------------------------
4153 subroutine sc_move(n_start,n_end,n_maxtry,e_drop,n_fun,etot)
4156 use random, only: iran_num
4157 use energy, only: esc
4158 ! Perform a quick search over side-chain arrangments (over
4159 ! residues n_start to n_end) for a given (frozen) CA trace
4160 ! Only side-chains are minimized (at most n_maxtry times each),
4162 ! Stops if energy drops by e_drop, otherwise tries all residues
4163 ! in the given range
4164 ! If there is an energy drop, full minimization may be useful
4165 ! n_start, n_end CAN be modified by this routine, but only if
4166 ! out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
4167 ! NOTE: this move should never increase the energy
4171 ! implicit real*8 (a-h,o-z)
4172 ! include 'DIMENSIONS'
4174 ! include 'COMMON.GEO'
4175 ! include 'COMMON.VAR'
4176 ! include 'COMMON.HEADER'
4177 ! include 'COMMON.IOUNITS'
4178 ! include 'COMMON.CHAIN'
4179 ! include 'COMMON.FFIELD'
4181 ! External functions
4182 !el integer iran_num
4183 !el external iran_num
4186 integer :: n_start,n_end,n_maxtry
4187 real(kind=8) :: e_drop
4191 real(kind=8) :: etot
4194 ! real(kind=8) :: energy(0:n_ene)
4195 real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1)
4196 real(kind=8) :: orig_e,cur_e
4197 integer :: n,n_steps,n_first,n_cur,n_tot !,i
4198 real(kind=8) :: orig_w(0:n_ene)
4199 real(kind=8) :: wtime
4201 !elwrite(iout,*) "in sc_move etot= ", etot
4202 ! Set non side-chain weights to zero (minimization is faster)
4203 ! NOTE: e(2) does not actually depend on the side-chain, only CA
4232 ! Make sure n_start, n_end are within proper range
4233 if (n_start.lt.2) n_start=2
4234 if (n_end.gt.nres-1) n_end=nres-1
4235 !rc if (n_start.lt.n_end) then
4236 if (n_start.gt.n_end) then
4241 ! Save the initial values of energy and coordinates
4243 !d call etotal(energy)
4244 !d write (iout,*) 'start sc ene',energy(0)
4245 !d call enerprint(energy(0))
4251 !rc cur_alph(i)=alph(i)
4252 !rc cur_omeg(i)=omeg(i)
4255 !t wtime=MPI_WTIME()
4256 ! Try (one by one) all specified residues, starting from a
4257 ! random position in sequence
4258 ! Stop early if the energy has decreased by at least e_drop
4259 n_tot=n_end-n_start+1
4260 n_first=iran_num(0,n_tot-1)
4263 !rc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
4264 do while (n.lt.n_tot)
4265 n_cur=n_start+mod(n_first+n,n_tot)
4266 call single_sc_move(n_cur,n_maxtry,e_drop,&
4268 !elwrite(iout,*) "after msingle sc_move etot= ", etot
4269 ! If a lower energy was found, update the current structure...
4270 !rc if (etot.lt.cur_e) then
4273 !rc cur_alph(i)=alph(i)
4274 !rc cur_omeg(i)=omeg(i)
4277 ! ...else revert to the previous one
4280 !rc alph(i)=cur_alph(i)
4281 !rc omeg(i)=cur_omeg(i)
4287 !d call etotal(energy)
4288 !d print *,'running',n,energy(0)
4292 !d call etotal(energy)
4293 !d write (iout,*) 'end sc ene',energy(0)
4295 ! Put the original weights back to calculate the full energy
4311 !t write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
4313 end subroutine sc_move
4314 !-----------------------------------------------------------------------------
4315 subroutine single_sc_move(res_pick,n_maxtry,e_drop,n_steps,n_fun,e_sc)
4317 ! Perturb one side-chain (res_pick) and minimize the
4318 ! neighbouring region, keeping all CA's and non-neighbouring
4320 ! Try until e_drop energy improvement is achieved, or n_maxtry
4321 ! attempts have been made
4322 ! At the start, e_sc should contain the side-chain-only energy(0)
4323 ! nsteps and nfun for this move are ADDED to n_steps and n_fun
4325 use energy, only: esc
4326 use geometry, only:dist
4328 ! implicit real*8 (a-h,o-z)
4329 ! include 'DIMENSIONS'
4330 ! include 'COMMON.VAR'
4331 ! include 'COMMON.INTERACT'
4332 ! include 'COMMON.CHAIN'
4333 ! include 'COMMON.MINIM'
4334 ! include 'COMMON.FFIELD'
4335 ! include 'COMMON.IOUNITS'
4337 ! External functions
4338 !el double precision dist
4342 integer :: res_pick,n_maxtry
4343 real(kind=8) :: e_drop
4345 ! Input/Output arguments
4346 integer :: n_steps,n_fun
4347 real(kind=8) :: e_sc
4352 integer :: nres_moved
4353 integer :: iretcode,loc_nfun,orig_maxfun,n_try
4354 real(kind=8) :: sc_dist,sc_dist_cutoff
4355 ! real(kind=8) :: energy_(0:n_ene)
4356 real(kind=8) :: evdw,escloc,orig_e,cur_e
4357 real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1)
4358 real(kind=8) :: var(6*nres) !(maxvar) (maxvar=6*maxres)
4360 real(kind=8) :: orig_theta(1:nres),orig_phi(1:nres),&
4361 orig_alph(1:nres),orig_omeg(1:nres)
4363 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4364 ! Define what is meant by "neighbouring side-chain"
4365 sc_dist_cutoff=5.0D0
4367 ! Don't do glycine or ends
4369 if (i.eq.10 .or. i.eq.ntyp1) return
4371 ! Freeze everything (later will relax only selected side-chains)
4379 ! Find the neighbours of the side-chain to move
4380 ! and save initial variables
4385 ! Don't do glycine (itype(j,1)==10)
4386 if (itype(i,1).ne.10) then
4387 sc_dist=dist(nres+i,nres+res_pick)
4389 sc_dist=sc_dist_cutoff
4391 if (sc_dist.lt.sc_dist_cutoff) then
4392 nres_moved=nres_moved+1
4398 write(iout,*) 'Warning calling chainbuild'
4402 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4403 !elwrite(iout,*) "in sinle wsc=",wsc,"evdw",evdw,"wscloc",wscloc,"escloc",escloc
4404 e_sc=wsc*evdw+wscloc*escloc
4405 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4406 !d call etotal(energy)
4407 !d print *,'new ',(energy(k),k=0,n_ene)
4412 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
4413 ! Move the selected residue (don't worry if it fails)
4414 call gen_side(iabs(itype(res_pick,1)),theta(res_pick+1),&
4415 alph(res_pick),omeg(res_pick),fail)
4417 ! Minimize the side-chains starting from the new arrangement
4418 call geom_to_var(nvar,var)
4423 !rc orig_theta(i)=theta(i)
4424 !rc orig_phi(i)=phi(i)
4425 !rc orig_alph(i)=alph(i)
4426 !rc orig_omeg(i)=omeg(i)
4429 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4430 call minimize_sc1(e_sc,var,iretcode,loc_nfun)
4432 !elwrite(iout,*) "in sinle etot/ e_sc",e_sc
4433 !v write(*,'(2i3,2f12.5,2i3)')
4434 !v & res_pick,nres_moved,orig_e,e_sc-cur_e,
4435 !v & iretcode,loc_nfun
4437 !$$$ if (iretcode.eq.8) then
4438 !$$$ write(iout,*)'Coordinates just after code 8'
4439 !$$$ call chainbuild
4440 !$$$ call all_varout
4441 !$$$ call flush(iout)
4443 !$$$ theta(i)=orig_theta(i)
4444 !$$$ phi(i)=orig_phi(i)
4445 !$$$ alph(i)=orig_alph(i)
4446 !$$$ omeg(i)=orig_omeg(i)
4448 !$$$ write(iout,*)'Coordinates just before code 8'
4449 !$$$ call chainbuild
4450 !$$$ call all_varout
4451 !$$$ call flush(iout)
4454 n_fun=n_fun+loc_nfun
4456 call var_to_geom(nvar,var)
4458 ! If a lower energy was found, update the current structure...
4459 if (e_sc.lt.cur_e) then
4461 !v call etotal(energy)
4464 !d e_sc1=wsc*evdw+wscloc*escloc
4465 !d print *,' new',e_sc1,energy(0)
4466 !v print *,'new ',energy(0)
4467 !d call enerprint(energy(0))
4470 if (mask_side(i).eq.1) then
4476 ! ...else revert to the previous one
4479 if (mask_side(i).eq.1) then
4488 n_steps=n_steps+n_try
4490 ! Reset the minimization mask_r to false
4494 end subroutine single_sc_move
4495 !-----------------------------------------------------------------------------
4496 subroutine sc_minimize(etot,iretcode,nfun)
4498 ! Minimizes side-chains only, leaving backbone frozen
4500 use energy, only: etotal
4502 ! implicit real*8 (a-h,o-z)
4503 ! include 'DIMENSIONS'
4504 ! include 'COMMON.VAR'
4505 ! include 'COMMON.CHAIN'
4506 ! include 'COMMON.FFIELD'
4509 real(kind=8) :: etot
4510 integer :: iretcode,nfun
4514 real(kind=8) :: orig_w(0:n_ene),energy_(0:n_ene)
4515 real(kind=8) :: var(6*nres) !(maxvar)(maxvar=6*maxres)
4518 ! Set non side-chain weights to zero (minimization is faster)
4519 ! NOTE: e(2) does not actually depend on the side-chain, only CA
4546 ! Prepare to freeze backbone
4553 ! Minimize the side-chains
4555 call geom_to_var(nvar,var)
4556 call minimize(etot,var,iretcode,nfun)
4557 call var_to_geom(nvar,var)
4560 ! Put the original weights back and calculate the full energy
4573 write(iout,*) 'Warning calling chainbuild'
4575 call etotal(energy_)
4579 end subroutine sc_minimize
4580 !-----------------------------------------------------------------------------
4581 subroutine minimize_sc1(etot,x,iretcode,nfun)
4583 use energy, only: func,gradient,fdum,etotal,enerprint
4585 ! implicit real*8 (a-h,o-z)
4586 ! include 'DIMENSIONS'
4587 integer,parameter :: liv=60
4588 integer :: iretcode,nfun
4589 ! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4590 ! include 'COMMON.IOUNITS'
4591 ! include 'COMMON.VAR'
4592 ! include 'COMMON.GEO'
4593 ! include 'COMMON.MINIM'
4594 !el integer :: icall
4595 !el common /srutu/ icall
4596 integer,dimension(liv) :: iv
4597 real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
4598 real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
4599 real(kind=8) :: energia(0:n_ene)
4600 !el real(kind=8) :: fdum
4601 ! external gradient,fdum !func,
4602 ! external func_restr1,grad_restr1
4603 logical :: not_done,change,reduce
4604 !el common /przechowalnia/ v
4606 integer :: nvar_restr,lv,i,j
4608 real(kind=8) :: rdum(1),etot !,fdum
4610 lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
4611 if (.not. allocated(v)) allocate(v(1:lv))
4613 call deflt(2,iv,liv,lv,v)
4614 ! 12 means fresh start, dont call deflt
4616 ! max num of fun calls
4617 if (maxfun.eq.0) maxfun=500
4619 ! max num of iterations
4620 if (maxmin.eq.0) maxmin=1000
4624 ! selects output unit
4627 ! 1 means to print out result
4629 ! 1 means to print out summary stats
4631 ! 1 means to print initial x and d
4633 ! min val for v(radfac) default is 0.1
4635 ! max val for v(radfac) default is 4.0
4638 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
4639 ! the sumsl default is 0.1
4641 ! false conv if (act fnctn decrease) .lt. v(34)
4642 ! the sumsl default is 100*machep
4644 ! absolute convergence
4645 if (tolf.eq.0.0D0) tolf=1.0D-4
4647 ! relative convergence
4648 if (rtolf.eq.0.0D0) rtolf=1.0D-4
4650 ! controls initial step size
4652 ! large vals of d correspond to small components of step
4661 call x2xx(x,xx,nvar_restr)
4662 call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,&
4663 iv,liv,lv,v,idum,rdum,fdum)
4666 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
4668 !el---------------------
4669 ! write (iout,'(/a)') &
4670 ! "Cartesian coordinates of the reference structure after SUMSL"
4671 ! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
4672 ! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4674 ! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
4675 ! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
4676 ! (c(j,i+nres),j=1,3)
4678 ! call etotal(energia)
4679 ! call enerprint(energia)
4680 !el----------------------------
4686 end subroutine minimize_sc1
4687 !-----------------------------------------------------------------------------
4688 subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
4691 use energy, only: zerograd,esc,sc_grad
4692 ! implicit real*8 (a-h,o-z)
4693 ! include 'DIMENSIONS'
4694 ! include 'COMMON.DERIV'
4695 ! include 'COMMON.IOUNITS'
4696 ! include 'COMMON.GEO'
4697 ! include 'COMMON.FFIELD'
4698 ! include 'COMMON.INTERACT'
4699 ! include 'COMMON.TIME1'
4701 !el common /chuju/ jjj
4702 real(kind=8) :: energia(0:n_ene),evdw,escloc
4703 real(kind=8) :: e1,e2,f
4704 real(kind=8),external :: ufparm
4705 integer :: uiparm(1)
4706 real(kind=8) :: urparm(1)
4707 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
4712 ! Intercept NaNs in the coordinates, before calling etotal
4718 if (x_sum.ne.x_sum) then
4719 write(iout,*)" *** func_restr1 : Found NaN in coordinates"
4726 call var_to_geom_restr(n,x)
4728 write(iout,*) 'Warning calling chainbuild'
4730 !d write (iout,*) 'ETOTAL called from FUNC'
4733 f=wsc*evdw+wscloc*escloc
4734 !d call etotal(energia(0))
4735 !d f=wsc*energia(1)+wscloc*energia(12)
4736 !d print *,f,evdw,escloc,energia(0)
4738 ! Sum up the components of the Cartesian gradient.
4742 gradx(j,i,icg)=wsc*gvdwx(j,i)
4747 end subroutine func_restr1
4748 !-----------------------------------------------------------------------------
4749 subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
4751 use energy, only: cartder,zerograd,esc,sc_grad
4752 ! implicit real*8 (a-h,o-z)
4753 ! include 'DIMENSIONS'
4754 ! include 'COMMON.CHAIN'
4755 ! include 'COMMON.DERIV'
4756 ! include 'COMMON.VAR'
4757 ! include 'COMMON.INTERACT'
4758 ! include 'COMMON.FFIELD'
4759 ! include 'COMMON.IOUNITS'
4761 integer :: i,j,k,ind,n,nf,uiparm(1)
4762 real(kind=8) :: f,urparm(1)
4763 real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
4764 integer :: ig,igall,ij
4765 real(kind=8) :: gphii,gthetai,galphai,gomegai
4766 real(kind=8),external :: ufparm
4769 if (nf-nfl+1) 20,30,40
4770 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
4771 ! write (iout,*) 'grad 20'
4774 30 call var_to_geom_restr(n,x)
4775 write(iout,*) 'Warning calling chainbuild'
4778 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
4782 ! Convert the Cartesian gradient into internal-coordinate gradient.
4788 IF (mask_phi(i+2).eq.1) THEN
4793 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
4794 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
4807 IF (mask_theta(i+2).eq.1) THEN
4813 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
4814 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
4824 if (itype(i,1).ne.10) then
4825 IF (mask_side(i).eq.1) THEN
4829 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
4838 if (itype(i,1).ne.10) then
4839 IF (mask_side(i).eq.1) THEN
4843 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
4851 ! Add the components corresponding to local energy terms.
4858 if (mask_phi(i).eq.1) then
4860 g(ig)=g(ig)+gloc(igall,icg)
4866 if (mask_theta(i).eq.1) then
4868 g(ig)=g(ig)+gloc(igall,icg)
4874 if (itype(i,1).ne.10) then
4876 if (mask_side(i).eq.1) then
4878 g(ig)=g(ig)+gloc(igall,icg)
4885 !d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
4888 end subroutine grad_restr1
4889 !-----------------------------------------------------------------------------
4890 subroutine egb1(evdw)
4892 ! This subroutine calculates the interaction energy of nonbonded side chains
4893 ! assuming the Gay-Berne potential of interaction.
4896 use energy, only: sc_grad
4897 ! use control, only:stopx
4898 ! implicit real*8 (a-h,o-z)
4899 ! include 'DIMENSIONS'
4900 ! include 'COMMON.GEO'
4901 ! include 'COMMON.VAR'
4902 ! include 'COMMON.LOCAL'
4903 ! include 'COMMON.CHAIN'
4904 ! include 'COMMON.DERIV'
4905 ! include 'COMMON.NAMES'
4906 ! include 'COMMON.INTERACT'
4907 ! include 'COMMON.IOUNITS'
4908 ! include 'COMMON.CALC'
4909 ! include 'COMMON.CONTROL'
4911 real(kind=8) :: evdw
4913 integer :: iint,ind,itypi,itypi1,itypj
4914 real(kind=8) :: xi,yi,zi,rrij,sig,sig0ij,rij_shift,fac,e1,e2,&
4916 !elwrite(iout,*) "check evdw"
4917 ! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
4920 ! if (icall.eq.0) lprn=.true.
4922 do i=iatsc_s,iatsc_e
4924 itypi=iabs(itype(i,1))
4925 itypi1=iabs(itype(i+1,1))
4929 dxi=dc_norm(1,nres+i)
4930 dyi=dc_norm(2,nres+i)
4931 dzi=dc_norm(3,nres+i)
4932 dsci_inv=dsc_inv(itypi)
4933 !elwrite(iout,*) itypi,itypi1,xi,yi,zi,dxi,dyi,dzi,dsci_inv
4935 ! Calculate SC interaction energy.
4937 do iint=1,nint_gr(i)
4938 do j=istart(i,iint),iend(i,iint)
4939 IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
4941 itypj=iabs(itype(j,1))
4942 dscj_inv=dsc_inv(itypj)
4943 sig0ij=sigma(itypi,itypj)
4944 chi1=chi(itypi,itypj)
4945 chi2=chi(itypj,itypi)
4952 alf12=0.5D0*(alf1+alf2)
4953 ! For diagnostics only!!!
4966 dxj=dc_norm(1,nres+j)
4967 dyj=dc_norm(2,nres+j)
4968 dzj=dc_norm(3,nres+j)
4969 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
4971 ! Calculate angle-dependent terms of energy and contributions to their
4975 sig=sig0ij*dsqrt(sigsq)
4976 rij_shift=1.0D0/rij-sig+sig0ij
4977 ! I hate to put IF's in the loops, but here don't have another choice!!!!
4978 if (rij_shift.le.0.0D0) then
4980 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
4981 !d restyp(itypi),i,restyp(itypj),j, &
4982 !d rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
4986 !---------------------------------------------------------------
4987 rij_shift=1.0D0/rij_shift
4988 fac=rij_shift**expon
4989 e1=fac*fac*aa_aq(itypi,itypj)
4990 e2=fac*bb_aq(itypi,itypj)
4991 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
4992 eps2der=evdwij*eps3rt
4993 eps3der=evdwij*eps2rt
4994 evdwij=evdwij*eps2rt*eps3rt
4996 if (wliptran.gt.0.0) print *,"WARNING eps_aq used!"
4998 sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
4999 epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
5000 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
5001 !d restyp(itypi),i,restyp(itypj),j, &
5002 !d epsi,sigm,chi1,chi2,chip1,chip2, &
5003 !d eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
5004 !d om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
5008 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
5011 ! Calculate gradient components.
5012 e1=e1*eps1*eps2rt**2*eps3rt**2
5013 fac=-expon*(e1+evdwij)*rij_shift
5016 ! Calculate the radial part of the gradient
5020 ! Calculate angular part of the gradient.
5022 !elwrite(iout,*) evdw
5024 !elwrite(iout,*) "evdw=",evdw,j,iint,i
5026 !elwrite(iout,*) evdw
5028 !elwrite(iout,*) evdw
5030 !elwrite(iout,*) evdw
5032 !elwrite(iout,*) evdw,i
5034 !-----------------------------------------------------------------------------
5036 !-----------------------------------------------------------------------------
5037 subroutine sumsl(n,d,x,calcf,calcg,iv,liv,lv,v,uiparm,urparm,ufparm)
5039 ! *** minimize general unconstrained objective function using ***
5040 ! *** analytic gradient and hessian approx. from secant update ***
5043 integer :: n, liv, lv
5044 integer :: iv(liv), uiparm(1)
5045 real(kind=8) :: d(n), x(n), v(lv), urparm(1)
5046 real(kind=8),external :: ufparm !funtion name as an argument
5048 ! dimension v(71 + n*(n+15)/2), uiparm(*), urparm(*)
5049 external :: calcf, calcg !subroutine name as an argument
5053 ! this routine interacts with subroutine sumit in an attempt
5054 ! to find an n-vector x* that minimizes the (unconstrained)
5055 ! objective function computed by calcf. (often the x* found is
5056 ! a local minimizer rather than a global one.)
5058 !-------------------------- parameter usage --------------------------
5060 ! n........ (input) the number of variables on which f depends, i.e.,
5061 ! the number of components in x.
5062 ! d........ (input/output) a scale vector such that d(i)*x(i),
5063 ! i = 1,2,...,n, are all in comparable units.
5064 ! d can strongly affect the behavior of sumsl.
5065 ! finding the best choice of d is generally a trial-
5066 ! and-error process. choosing d so that d(i)*x(i)
5067 ! has about the same value for all i often works well.
5068 ! the defaults provided by subroutine deflt (see i
5069 ! below) require the caller to supply d.
5070 ! x........ (input/output) before (initially) calling sumsl, the call-
5071 ! er should set x to an initial guess at x*. when
5072 ! sumsl returns, x contains the best point so far
5073 ! found, i.e., the one that gives the least value so
5074 ! far seen for f(x).
5075 ! calcf.... (input) a subroutine that, given x, computes f(x). calcf
5076 ! must be declared external in the calling program.
5078 ! call calcf(n, x, nf, f, uiparm, urparm, ufparm)
5079 ! when calcf is called, nf is the invocation
5080 ! count for calcf. nf is included for possible use
5081 ! with calcg. if x is out of bounds (e.g., if it
5082 ! would cause overflow in computing f(x)), then calcf
5083 ! should set nf to 0. this will cause a shorter step
5084 ! to be attempted. (if x is in bounds, then calcf
5085 ! should not change nf.) the other parameters are as
5086 ! described above and below. calcf should not change
5088 ! calcg.... (input) a subroutine that, given x, computes g(x), the gra-
5089 ! dient of f at x. calcg must be declared external in
5090 ! the calling program. it is invoked by
5091 ! call calcg(n, x, nf, g, uiparm, urparm, ufaprm)
5092 ! when calcg is called, nf is the invocation
5093 ! count for calcf at the time f(x) was evaluated. the
5094 ! x passed to calcg is usually the one passed to calcf
5095 ! on either its most recent invocation or the one
5096 ! prior to it. if calcf saves intermediate results
5097 ! for use by calcg, then it is possible to tell from
5098 ! nf whether they are valid for the current x (or
5099 ! which copy is valid if two copies are kept). if g
5100 ! cannot be computed at x, then calcg should set nf to
5101 ! 0. in this case, sumsl will return with iv(1) = 65.
5102 ! (if g can be computed at x, then calcg should not
5103 ! changed nf.) the other parameters to calcg are as
5104 ! described above and below. calcg should not change
5106 ! iv....... (input/output) an integer value array of length liv (see
5107 ! below) that helps control the sumsl algorithm and
5108 ! that is used to store various intermediate quanti-
5109 ! ties. of particular interest are the initialization/
5110 ! return code iv(1) and the entries in iv that control
5111 ! printing and limit the number of iterations and func-
5112 ! tion evaluations. see the section on iv input
5114 ! liv...... (input) length of iv array. must be at least 60. if li
5115 ! is too small, then sumsl returns with iv(1) = 15.
5116 ! when sumsl returns, the smallest allowed value of
5117 ! liv is stored in iv(lastiv) -- see the section on
5118 ! iv output values below. (this is intended for use
5119 ! with extensions of sumsl that handle constraints.)
5120 ! lv....... (input) length of v array. must be at least 71+n*(n+15)/2.
5121 ! (at least 77+n*(n+17)/2 for smsno, at least
5122 ! 78+n*(n+12) for humsl). if lv is too small, then
5123 ! sumsl returns with iv(1) = 16. when sumsl returns,
5124 ! the smallest allowed value of lv is stored in
5125 ! iv(lastv) -- see the section on iv output values
5127 ! v........ (input/output) a floating-point value array of length l
5128 ! (see below) that helps control the sumsl algorithm
5129 ! and that is used to store various intermediate
5130 ! quantities. of particular interest are the entries
5131 ! in v that limit the length of the first step
5132 ! attempted (lmax0) and specify convergence tolerances
5133 ! (afctol, lmaxs, rfctol, sctol, xctol, xftol).
5134 ! uiparm... (input) user integer parameter array passed without change
5135 ! to calcf and calcg.
5136 ! urparm... (input) user floating-point parameter array passed without
5137 ! change to calcf and calcg.
5138 ! ufparm... (input) user external subroutine or function passed without
5139 ! change to calcf and calcg.
5141 ! *** iv input values (from subroutine deflt) ***
5143 ! iv(1)... on input, iv(1) should have a value between 0 and 14......
5144 ! 0 and 12 mean this is a fresh start. 0 means that
5145 ! deflt(2, iv, liv, lv, v)
5146 ! is to be called to provide all default values to iv and
5147 ! v. 12 (the value that deflt assigns to iv(1)) means the
5148 ! caller has already called deflt and has possibly changed
5149 ! some iv and/or v entries to non-default values.
5150 ! 13 means deflt has been called and that sumsl (and
5151 ! sumit) should only do their storage allocation. that is,
5152 ! they should set the output components of iv that tell
5153 ! where various subarrays arrays of v begin, such as iv(g)
5154 ! (and, for humsl and humit only, iv(dtol)), and return.
5155 ! 14 means that a storage has been allocated (by a call
5156 ! with iv(1) = 13) and that the algorithm should be
5157 ! started. when called with iv(1) = 13, sumsl returns
5158 ! iv(1) = 14 unless liv or lv is too small (or n is not
5159 ! positive). default = 12.
5160 ! iv(inith).... iv(25) tells whether the hessian approximation h should
5161 ! be initialized. 1 (the default) means sumit should
5162 ! initialize h to the diagonal matrix whose i-th diagonal
5163 ! element is d(i)**2. 0 means the caller has supplied a
5164 ! cholesky factor l of the initial hessian approximation
5165 ! h = l*(l**t) in v, starting at v(iv(lmat)) = v(iv(42))
5166 ! (and stored compactly by rows). note that iv(lmat) may
5167 ! be initialized by calling sumsl with iv(1) = 13 (see
5168 ! the iv(1) discussion above). default = 1.
5169 ! iv(mxfcal)... iv(17) gives the maximum number of function evaluations
5170 ! (calls on calcf) allowed. if this number does not suf-
5171 ! fice, then sumsl returns with iv(1) = 9. default = 200.
5172 ! iv(mxiter)... iv(18) gives the maximum number of iterations allowed.
5173 ! it also indirectly limits the number of gradient evalua-
5174 ! tions (calls on calcg) to iv(mxiter) + 1. if iv(mxiter)
5175 ! iterations do not suffice, then sumsl returns with
5176 ! iv(1) = 10. default = 150.
5177 ! iv(outlev)... iv(19) controls the number and length of iteration sum-
5178 ! mary lines printed (by itsum). iv(outlev) = 0 means do
5179 ! not print any summary lines. otherwise, print a summary
5180 ! line after each abs(iv(outlev)) iterations. if iv(outlev)
5181 ! is positive, then summary lines of length 78 (plus carri-
5182 ! age control) are printed, including the following... the
5183 ! iteration and function evaluation counts, f = the current
5184 ! function value, relative difference in function values
5185 ! achieved by the latest step (i.e., reldf = (f0-v(f))/f01,
5186 ! where f01 is the maximum of abs(v(f)) and abs(v(f0)) and
5187 ! v(f0) is the function value from the previous itera-
5188 ! tion), the relative function reduction predicted for the
5189 ! step just taken (i.e., preldf = v(preduc) / f01, where
5190 ! v(preduc) is described below), the scaled relative change
5191 ! in x (see v(reldx) below), the step parameter for the
5192 ! step just taken (stppar = 0 means a full newton step,
5193 ! between 0 and 1 means a relaxed newton step, between 1
5194 ! and 2 means a double dogleg step, greater than 2 means
5195 ! a scaled down cauchy step -- see subroutine dbldog), the
5196 ! 2-norm of the scale vector d times the step just taken
5197 ! (see v(dstnrm) below), and npreldf, i.e.,
5198 ! v(nreduc)/f01, where v(nreduc) is described below -- if
5199 ! npreldf is positive, then it is the relative function
5200 ! reduction predicted for a newton step (one with
5201 ! stppar = 0). if npreldf is negative, then it is the
5202 ! negative of the relative function reduction predicted
5203 ! for a step computed with step bound v(lmaxs) for use in
5204 ! testing for singular convergence.
5205 ! if iv(outlev) is negative, then lines of length 50
5206 ! are printed, including only the first 6 items listed
5207 ! above (through reldx).
5209 ! iv(parprt)... iv(20) = 1 means print any nondefault v values on a
5210 ! fresh start or any changed v values on a restart.
5211 ! iv(parprt) = 0 means skip this printing. default = 1.
5212 ! iv(prunit)... iv(21) is the output unit number on which all printing
5213 ! is done. iv(prunit) = 0 means suppress all printing.
5214 ! default = standard output unit (unit 6 on most systems).
5215 ! iv(solprt)... iv(22) = 1 means print out the value of x returned (as
5216 ! well as the gradient and the scale vector d).
5217 ! iv(solprt) = 0 means skip this printing. default = 1.
5218 ! iv(statpr)... iv(23) = 1 means print summary statistics upon return-
5219 ! ing. these consist of the function value, the scaled
5220 ! relative change in x caused by the most recent step (see
5221 ! v(reldx) below), the number of function and gradient
5222 ! evaluations (calls on calcf and calcg), and the relative
5223 ! function reductions predicted for the last step taken and
5224 ! for a newton step (or perhaps a step bounded by v(lmaxs)
5225 ! -- see the descriptions of preldf and npreldf under
5226 ! iv(outlev) above).
5227 ! iv(statpr) = 0 means skip this printing.
5228 ! iv(statpr) = -1 means skip this printing as well as that
5229 ! of the one-line termination reason message. default = 1.
5230 ! iv(x0prt).... iv(24) = 1 means print the initial x and scale vector d
5231 ! (on a fresh start only). iv(x0prt) = 0 means skip this
5232 ! printing. default = 1.
5234 ! *** (selected) iv output values ***
5236 ! iv(1)........ on output, iv(1) is a return code....
5237 ! 3 = x-convergence. the scaled relative difference (see
5238 ! v(reldx)) between the current parameter vector x and
5239 ! a locally optimal parameter vector is very likely at
5241 ! 4 = relative function convergence. the relative differ-
5242 ! ence between the current function value and its lo-
5243 ! cally optimal value is very likely at most v(rfctol).
5244 ! 5 = both x- and relative function convergence (i.e., the
5245 ! conditions for iv(1) = 3 and iv(1) = 4 both hold).
5246 ! 6 = absolute function convergence. the current function
5247 ! value is at most v(afctol) in absolute value.
5248 ! 7 = singular convergence. the hessian near the current
5249 ! iterate appears to be singular or nearly so, and a
5250 ! step of length at most v(lmaxs) is unlikely to yield
5251 ! a relative function decrease of more than v(sctol).
5252 ! 8 = false convergence. the iterates appear to be converg-
5253 ! ing to a noncritical point. this may mean that the
5254 ! convergence tolerances (v(afctol), v(rfctol),
5255 ! v(xctol)) are too small for the accuracy to which
5256 ! the function and gradient are being computed, that
5257 ! there is an error in computing the gradient, or that
5258 ! the function or gradient is discontinuous near x.
5259 ! 9 = function evaluation limit reached without other con-
5260 ! vergence (see iv(mxfcal)).
5261 ! 10 = iteration limit reached without other convergence
5263 ! 11 = stopx returned .true. (external interrupt). see the
5264 ! usage notes below.
5265 ! 14 = storage has been allocated (after a call with
5267 ! 17 = restart attempted with n changed.
5268 ! 18 = d has a negative component and iv(dtype) .le. 0.
5269 ! 19...43 = v(iv(1)) is out of range.
5270 ! 63 = f(x) cannot be computed at the initial x.
5271 ! 64 = bad parameters passed to assess (which should not
5273 ! 65 = the gradient could not be computed at x (see calcg
5275 ! 67 = bad first parameter to deflt.
5276 ! 80 = iv(1) was out of range.
5277 ! 81 = n is not positive.
5278 ! iv(g)........ iv(28) is the starting subscript in v of the current
5279 ! gradient vector (the one corresponding to x).
5280 ! iv(lastiv)... iv(44) is the least acceptable value of liv. (it is
5281 ! only set if liv is at least 44.)
5282 ! iv(lastv).... iv(45) is the least acceptable value of lv. (it is
5283 ! only set if liv is large enough, at least iv(lastiv).)
5284 ! iv(nfcall)... iv(6) is the number of calls so far made on calcf (i.e.,
5285 ! function evaluations).
5286 ! iv(ngcall)... iv(30) is the number of gradient evaluations (calls on
5288 ! iv(niter).... iv(31) is the number of iterations performed.
5290 ! *** (selected) v input values (from subroutine deflt) ***
5292 ! v(bias)..... v(43) is the bias parameter used in subroutine dbldog --
5293 ! see that subroutine for details. default = 0.8.
5294 ! v(afctol)... v(31) is the absolute function convergence tolerance.
5295 ! if sumsl finds a point where the function value is less
5296 ! than v(afctol) in absolute value, and if sumsl does not
5297 ! return with iv(1) = 3, 4, or 5, then it returns with
5298 ! iv(1) = 6. this test can be turned off by setting
5299 ! v(afctol) to zero. default = max(10**-20, machep**2),
5300 ! where machep is the unit roundoff.
5301 ! v(dinit).... v(38), if nonnegative, is the value to which the scale
5302 ! vector d is initialized. default = -1.
5303 ! v(lmax0).... v(35) gives the maximum 2-norm allowed for d times the
5304 ! very first step that sumsl attempts. this parameter can
5305 ! markedly affect the performance of sumsl.
5306 ! v(lmaxs).... v(36) is used in testing for singular convergence -- if
5307 ! the function reduction predicted for a step of length
5308 ! bounded by v(lmaxs) is at most v(sctol) * abs(f0), where
5309 ! f0 is the function value at the start of the current
5310 ! iteration, and if sumsl does not return with iv(1) = 3,
5311 ! 4, 5, or 6, then it returns with iv(1) = 7. default = 1.
5312 ! v(rfctol)... v(32) is the relative function convergence tolerance.
5313 ! if the current model predicts a maximum possible function
5314 ! reduction (see v(nreduc)) of at most v(rfctol)*abs(f0)
5315 ! at the start of the current iteration, where f0 is the
5316 ! then current function value, and if the last step attempt-
5317 ! ed achieved no more than twice the predicted function
5318 ! decrease, then sumsl returns with iv(1) = 4 (or 5).
5319 ! default = max(10**-10, machep**(2/3)), where machep is
5320 ! the unit roundoff.
5321 ! v(sctol).... v(37) is the singular convergence tolerance -- see the
5322 ! description of v(lmaxs) above.
5323 ! v(tuner1)... v(26) helps decide when to check for false convergence.
5324 ! this is done if the actual function decrease from the
5325 ! current step is no more than v(tuner1) times its predict-
5326 ! ed value. default = 0.1.
5327 ! v(xctol).... v(33) is the x-convergence tolerance. if a newton step
5328 ! (see v(nreduc)) is tried that has v(reldx) .le. v(xctol)
5329 ! and if this step yields at most twice the predicted func-
5330 ! tion decrease, then sumsl returns with iv(1) = 3 (or 5).
5331 ! (see the description of v(reldx) below.)
5332 ! default = machep**0.5, where machep is the unit roundoff.
5333 ! v(xftol).... v(34) is the false convergence tolerance. if a step is
5334 ! tried that gives no more than v(tuner1) times the predict-
5335 ! ed function decrease and that has v(reldx) .le. v(xftol),
5336 ! and if sumsl does not return with iv(1) = 3, 4, 5, 6, or
5337 ! 7, then it returns with iv(1) = 8. (see the description
5338 ! of v(reldx) below.) default = 100*machep, where
5339 ! machep is the unit roundoff.
5340 ! v(*)........ deflt supplies to v a number of tuning constants, with
5341 ! which it should ordinarily be unnecessary to tinker. see
5342 ! section 17 of version 2.2 of the nl2sol usage summary
5343 ! (i.e., the appendix to ref. 1) for details on v(i),
5344 ! i = decfac, incfac, phmnfc, phmxfc, rdfcmn, rdfcmx,
5345 ! tuner2, tuner3, tuner4, tuner5.
5347 ! *** (selected) v output values ***
5349 ! v(dgnorm)... v(1) is the 2-norm of (diag(d)**-1)*g, where g is the
5350 ! most recently computed gradient.
5351 ! v(dstnrm)... v(2) is the 2-norm of diag(d)*step, where step is the
5353 ! v(f)........ v(10) is the current function value.
5354 ! v(f0)....... v(13) is the function value at the start of the current
5356 ! v(nreduc)... v(6), if positive, is the maximum function reduction
5357 ! possible according to the current model, i.e., the func-
5358 ! tion reduction predicted for a newton step (i.e.,
5359 ! step = -h**-1 * g, where g is the current gradient and
5360 ! h is the current hessian approximation).
5361 ! if v(nreduc) is negative, then it is the negative of
5362 ! the function reduction predicted for a step computed with
5363 ! a step bound of v(lmaxs) for use in testing for singular
5365 ! v(preduc)... v(7) is the function reduction predicted (by the current
5366 ! quadratic model) for the current step. this (divided by
5367 ! v(f0)) is used in testing for relative function
5369 ! v(reldx).... v(17) is the scaled relative change in x caused by the
5370 ! current step, computed as
5371 ! max(abs(d(i)*(x(i)-x0(i)), 1 .le. i .le. p) /
5372 ! max(d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p),
5373 ! where x = x0 + step.
5375 !------------------------------- notes -------------------------------
5377 ! *** algorithm notes ***
5379 ! this routine uses a hessian approximation computed from the
5380 ! bfgs update (see ref 3). only a cholesky factor of the hessian
5381 ! approximation is stored, and this is updated using ideas from
5382 ! ref. 4. steps are computed by the double dogleg scheme described
5383 ! in ref. 2. the steps are assessed as in ref. 1.
5385 ! *** usage notes ***
5387 ! after a return with iv(1) .le. 11, it is possible to restart,
5388 ! i.e., to change some of the iv and v input values described above
5389 ! and continue the algorithm from the point where it was interrupt-
5390 ! ed. iv(1) should not be changed, nor should any entries of i
5391 ! and v other than the input values (those supplied by deflt).
5392 ! those who do not wish to write a calcg which computes the
5393 ! gradient analytically should call smsno rather than sumsl.
5394 ! smsno uses finite differences to compute an approximate gradient.
5395 ! those who would prefer to provide f and g (the function and
5396 ! gradient) by reverse communication rather than by writing subrou-
5397 ! tines calcf and calcg may call on sumit directly. see the com-
5398 ! ments at the beginning of sumit.
5399 ! those who use sumsl interactively may wish to supply their
5400 ! own stopx function, which should return .true. if the break key
5401 ! has been pressed since stopx was last invoked. this makes it
5402 ! possible to externally interrupt sumsl (which will return with
5403 ! iv(1) = 11 if stopx returns .true.).
5404 ! storage for g is allocated at the end of v. thus the caller
5405 ! may make v longer than specified above and may allow calcg to use
5406 ! elements of g beyond the first n as scratch storage.
5408 ! *** portability notes ***
5410 ! the sumsl distribution tape contains both single- and double-
5411 ! precision versions of the sumsl source code, so it should be un-
5412 ! necessary to change precisions.
5413 ! only the functions imdcon and rmdcon contain machine-dependent
5414 ! constants. to change from one machine to another, it should
5415 ! suffice to change the (few) relevant lines in these functions.
5416 ! intrinsic functions are explicitly declared. on certain com-
5417 ! puters (e.g. univac), it may be necessary to comment out these
5418 ! declarations. so that this may be done automatically by a simple
5419 ! program, such declarations are preceded by a comment having c/+
5420 ! in columns 1-3 and blanks in columns 4-72 and are followed by
5421 ! a comment having c/ in columns 1 and 2 and blanks in columns 3-72.
5422 ! the sumsl source code is expressed in 1966 ansi standard
5423 ! fortran. it may be converted to fortran 77 by commenting out all
5424 ! lines that fall between a line having c/6 in columns 1-3 and a
5425 ! line having c/7 in columns 1-3 and by removing (i.e., replacing
5426 ! by a blank) the c in column 1 of the lines that follow the c/7
5427 ! line and precede a line having c/ in columns 1-2 and blanks in
5428 ! columns 3-72. these changes convert some data statements into
5429 ! parameter statements, convert some variables from real to
5430 ! character*4, and make the data statements that initialize these
5431 ! variables use character strings delimited by primes instead
5432 ! of hollerith constants. (such variables and data statements
5433 ! appear only in modules itsum and parck. parameter statements
5434 ! appear nearly everywhere.) these changes also add save state-
5435 ! ments for variables given machine-dependent constants by rmdcon.
5437 ! *** references ***
5439 ! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), algorithm 573 --
5440 ! an adaptive nonlinear least-squares algorithm, acm trans.
5441 ! math. software 7, pp. 369-383.
5443 ! 2. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
5444 ! mization algorithms which use function and gradient
5445 ! values, j. optim. theory applic. 28, pp. 453-482.
5447 ! 3. dennis, j.e., and more, j.j. (1977), quasi-newton methods, motiva-
5448 ! tion and theory, siam rev. 19, pp. 46-89.
5450 ! 4. goldfarb, d. (1976), factorized variable metric methods for uncon-
5451 ! strained optimization, math. comput. 30, pp. 796-811.
5455 ! coded by david m. gay (winter 1980). revised summer 1982.
5456 ! this subroutine was written in connection with research
5457 ! supported in part by the national science foundation under
5458 ! grants mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989,
5462 !---------------------------- declarations ---------------------------
5464 !el external deflt, sumit
5466 ! deflt... supplies default iv and v input components.
5467 ! sumit... reverse-communication routine that carries out sumsl algo-
5470 integer :: g1, iv1, nf
5473 ! *** subscripts for iv ***
5475 !el integer nextv, nfcall, nfgcal, g, toobig, vneed
5478 ! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, toobig/2/, vneed/4/
5480 integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28,&
5484 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
5486 !elwrite(iout,*) "in sumsl"
5487 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
5489 if (iv1 .eq. 12 .or. iv1 .eq. 13) iv(vneed) = iv(vneed) + n
5490 if (iv1 .eq. 14) go to 10
5491 if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
5493 if (iv1 .eq. 12) iv(1) = 13
5497 !elwrite(iout,*) "in sumsl go to 10"
5500 !elwrite(iout,*) "in sumsl"
5501 20 call sumit(d, f, v(g1), iv, liv, lv, n, v, x)
5502 !elwrite(iout,*) "in sumsl, go to 20"
5504 !elwrite(iout,*) "in sumsl, go to 20, po sumit"
5505 !elwrite(iout,*) "in sumsl iv()", iv(1)-2
5506 if (iv(1) - 2) 30, 40, 50
5509 !elwrite(iout,*) "in sumsl iv",iv(nfcall)
5510 call calcf(n, x, nf, f, uiparm, urparm, ufparm)
5511 !elwrite(iout,*) "in sumsl"
5512 if (nf .le. 0) iv(toobig) = 1
5515 !elwrite(iout,*) "in sumsl"
5516 40 call calcg(n, x, iv(nfgcal), v(g1), uiparm, urparm, ufparm)
5517 !elwrite(iout,*) "in sumsl"
5520 50 if (iv(1) .ne. 14) go to 999
5522 ! *** storage allocation
5525 iv(nextv) = iv(g) + n
5526 if (iv1 .ne. 13) go to 10
5527 !elwrite(iout,*) "in sumsl"
5530 ! *** last card of sumsl follows ***
5531 end subroutine sumsl
5532 !-----------------------------------------------------------------------------
5533 subroutine sumit(d,fx,g,iv,liv,lv,n,v,x)
5535 use control, only:stopx
5537 ! *** carry out sumsl (unconstrained minimization) iterations, using
5538 ! *** double-dogleg/bfgs steps.
5540 ! *** parameter declarations ***
5542 integer :: liv, lv, n
5544 real(kind=8) :: d(n), fx, g(n), v(lv), x(n)
5546 !-------------------------- parameter usage --------------------------
5548 ! d.... scale vector.
5549 ! fx... function value.
5550 ! g.... gradient vector.
5551 ! iv... integer value array.
5552 ! liv.. length of iv (at least 60).
5553 ! lv... length of v (at least 71 + n*(n+13)/2).
5554 ! n.... number of variables (components in x and g).
5555 ! v.... floating-point value array.
5556 ! x.... vector of parameters to be optimized.
5558 ! *** discussion ***
5560 ! parameters iv, n, v, and x are the same as the corresponding
5561 ! ones to sumsl (which see), except that v can be shorter (since
5562 ! the part of v that sumsl uses for storing g is not needed).
5563 ! moreover, compared with sumsl, iv(1) may have the two additional
5564 ! output values 1 and 2, which are explained below, as is the use
5565 ! of iv(toobig) and iv(nfgcal). the value iv(g), which is an
5566 ! output value from sumsl (and smsno), is not referenced by
5567 ! sumit or the subroutines it calls.
5568 ! fx and g need not have been initialized when sumit is called
5569 ! with iv(1) = 12, 13, or 14.
5571 ! iv(1) = 1 means the caller should set fx to f(x), the function value
5572 ! at x, and call sumit again, having changed none of the
5573 ! other parameters. an exception occurs if f(x) cannot be
5574 ! (e.g. if overflow would occur), which may happen because
5575 ! of an oversized step. in this case the caller should set
5576 ! iv(toobig) = iv(2) to 1, which will cause sumit to ig-
5577 ! nore fx and try a smaller step. the parameter nf that
5578 ! sumsl passes to calcf (for possible use by calcg) is a
5579 ! copy of iv(nfcall) = iv(6).
5580 ! iv(1) = 2 means the caller should set g to g(x), the gradient vector
5581 ! of f at x, and call sumit again, having changed none of
5582 ! the other parameters except possibly the scale vector d
5583 ! when iv(dtype) = 0. the parameter nf that sumsl passes
5584 ! to calcg is iv(nfgcal) = iv(7). if g(x) cannot be
5585 ! evaluated, then the caller may set iv(nfgcal) to 0, in
5586 ! which case sumit will return with iv(1) = 65.
5590 ! coded by david m. gay (december 1979). revised sept. 1982.
5591 ! this subroutine was written in connection with research supported
5592 ! in part by the national science foundation under grants
5593 ! mcs-7600324 and mcs-7906671.
5595 ! (see sumsl for references.)
5597 !+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++
5599 ! *** local variables ***
5601 integer :: dg1, dummy, g01, i, k, l, lstgst, nwtst1, step1,&
5604 !el logical :: lstopx
5608 !el real(kind=8) :: half, negone, one, onep2, zero
5610 ! *** no intrinsic functions ***
5612 ! *** external functions and subroutines ***
5614 !el external assst, dbdog, deflt, dotprd, itsum, litvmu, livmul,
5615 !el 1 ltvmul, lupdat, lvmul, parck, reldst, stopx, vaxpy,
5616 !el 2 vcopy, vscopy, vvmulp, v2norm, wzbfgs
5618 !el real(kind=8) :: dotprd, reldst, v2norm
5620 ! assst.... assesses candidate step.
5621 ! dbdog.... computes double-dogleg (candidate) step.
5622 ! deflt.... supplies default iv and v input components.
5623 ! dotprd... returns inner product of two vectors.
5624 ! itsum.... prints iteration summary and info on initial and final x.
5625 ! litvmu... multiplies inverse transpose of lower triangle times vector.
5626 ! livmul... multiplies inverse of lower triangle times vector.
5627 ! ltvmul... multiplies transpose of lower triangle times vector.
5628 ! lupdt.... updates cholesky factor of hessian approximation.
5629 ! lvmul.... multiplies lower triangle times vector.
5630 ! parck.... checks validity of input iv and v values.
5631 ! reldst... computes v(reldx) = relative step size.
5632 ! stopx.... returns .true. if the break key has been pressed.
5633 ! vaxpy.... computes scalar times one vector plus another.
5634 ! vcopy.... copies one vector to another.
5635 ! vscopy... sets all elements of a vector to a scalar.
5636 ! vvmulp... multiplies vector by vector raised to power (componentwise).
5637 ! v2norm... returns the 2-norm of a vector.
5638 ! wzbfgs... computes w and z for lupdat corresponding to bfgs update.
5640 ! *** subscripts for iv and v ***
5643 !el integer cnvcod, dg, dgnorm, dinit, dstnrm, dst0, f, f0, fdif,
5644 !el 1 gthg, gtstep, g0, incfac, inith, irc, kagqt, lmat, lmax0,
5645 !el 2 lmaxs, mode, model, mxfcal, mxiter, nextv, nfcall, nfgcal,
5646 !el 3 ngcall, niter, nreduc, nwtstp, preduc, radfac, radinc,
5647 !el 4 radius, rad0, reldx, restor, step, stglim, stlstg, toobig,
5648 !el 5 tuner4, tuner5, vneed, xirc, x0
5650 ! *** iv subscript values ***
5653 ! data cnvcod/55/, dg/37/, g0/48/, inith/25/, irc/29/, kagqt/33/,
5654 ! 1 mode/35/, model/5/, mxfcal/17/, mxiter/18/, nfcall/6/,
5655 ! 2 nfgcal/7/, ngcall/30/, niter/31/, nwtstp/34/, radinc/8/,
5656 ! 3 restor/9/, step/40/, stglim/11/, stlstg/41/, toobig/2/,
5657 ! 4 vneed/4/, xirc/13/, x0/43/
5659 integer,parameter :: cnvcod=55, dg=37, g0=48, inith=25, irc=29, kagqt=33,&
5660 mode=35, model=5, mxfcal=17, mxiter=18, nfcall=6,&
5661 nfgcal=7, ngcall=30, niter=31, nwtstp=34, radinc=8,&
5662 restor=9, step=40, stglim=11, stlstg=41, toobig=2,&
5663 vneed=4, xirc=13, x0=43
5666 ! *** v subscript values ***
5670 ! data dgnorm/1/, dinit/38/, dstnrm/2/, dst0/3/, f/10/, f0/13/,
5671 ! 1 fdif/11/, gthg/44/, gtstep/4/, incfac/23/, lmat/42/,
5672 ! 2 lmax0/35/, lmaxs/36/, nextv/47/, nreduc/6/, preduc/7/,
5673 ! 3 radfac/16/, radius/8/, rad0/9/, reldx/17/, tuner4/29/,
5676 integer,parameter :: afctol=31
5677 integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dst0=3, f=10, f0=13,&
5678 fdif=11, gthg=44, gtstep=4, incfac=23, lmat=42,&
5679 lmax0=35, lmaxs=36, nextv=47, nreduc=6, preduc=7,&
5680 radfac=16, radius=8, rad0=9, reldx=17, tuner4=29,&
5685 ! data half/0.5d+0/, negone/-1.d+0/, one/1.d+0/, onep2/1.2d+0/,
5688 real(kind=8),parameter :: half=0.5d+0, negone=-1.d+0, one=1.d+0,&
5689 onep2=1.2d+0,zero=0.d+0
5692 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
5694 ! Following SAVE statement inserted.
5697 if (i .eq. 1) go to 50
5698 if (i .eq. 2) go to 60
5700 ! *** check validity of iv and v input values ***
5702 if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
5703 if (iv(1) .eq. 12 .or. iv(1) .eq. 13) &
5704 iv(vneed) = iv(vneed) + n*(n+13)/2
5705 call parck(2, d, iv, liv, lv, n, v)
5707 if (i .gt. 12) go to 999
5708 go to (180, 180, 180, 180, 180, 180, 120, 90, 120, 10, 10, 20), i
5710 ! *** storage allocation ***
5713 iv(x0) = l + n*(n+1)/2
5714 iv(step) = iv(x0) + n
5715 iv(stlstg) = iv(step) + n
5716 iv(g0) = iv(stlstg) + n
5717 iv(nwtstp) = iv(g0) + n
5718 iv(dg) = iv(nwtstp) + n
5719 iv(nextv) = iv(dg) + n
5720 if (iv(1) .ne. 13) go to 20
5724 ! *** initialization ***
5737 if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
5738 if (iv(inith) .ne. 1) go to 40
5740 ! *** set the initial hessian approximation to diag(d)**-2 ***
5743 call vscopy(n*(n+1)/2, v(l), zero)
5748 if (t .le. zero) t = one
5752 ! *** compute initial function value ***
5758 if (iv(mode) .ge. 0) go to 180
5760 if (iv(toobig) .eq. 0) go to 999
5764 ! *** make sure gradient could be computed ***
5766 60 if (iv(nfgcal) .ne. 0) go to 70
5771 call vvmulp(n, v(dg1), g, d, -1)
5772 v(dgnorm) = v2norm(n, v(dg1))
5774 ! *** test norm of gradient ***
5776 if (v(dgnorm) .gt. v(afctol)) go to 75
5778 iv(cnvcod) = iv(irc) - 4
5780 75 if (iv(cnvcod) .ne. 0) go to 290
5781 if (iv(mode) .eq. 0) go to 250
5783 ! *** allow first step to have scaled 2-norm at most v(lmax0) ***
5785 v(radius) = v(lmax0)
5790 !----------------------------- main loop -----------------------------
5793 ! *** print iteration summary, check iteration limit ***
5795 80 call itsum(d, g, iv, liv, lv, n, v, x)
5797 if (k .lt. iv(mxiter)) go to 100
5801 ! *** update radius ***
5803 100 iv(niter) = k + 1
5804 if(k.gt.0)v(radius) = v(radfac) * v(dstnrm)
5806 ! *** initialize for start of next iteration ***
5814 ! *** copy x to x0, g to g0 ***
5816 call vcopy(n, v(x01), x)
5817 call vcopy(n, v(g01), g)
5819 ! *** check stopx and function evaluation limit ***
5823 !el lstopx = stopx(dummy)
5824 !elwrite(iout,*) "lstopx",lstopx,dummy
5825 110 if (.not. stopx(dummy)) go to 130
5827 ! write (iout,*) "iv(1)=11 !!!!"
5830 ! *** come here when restarting after func. eval. limit or stopx.
5832 120 if (v(f) .ge. v(f0)) go to 130
5837 130 if (iv(nfcall) .lt. iv(mxfcal)) go to 150
5839 140 if (v(f) .ge. v(f0)) go to 300
5841 ! *** in case of stopx or function evaluation limit with
5842 ! *** improved v(f), evaluate the gradient at x.
5847 !. . . . . . . . . . . . . compute candidate step . . . . . . . . . .
5849 150 step1 = iv(step)
5852 if (iv(kagqt) .ge. 0) go to 160
5854 call livmul(n, v(nwtst1), v(l), g)
5855 v(nreduc) = half * dotprd(n, v(nwtst1), v(nwtst1))
5856 call litvmu(n, v(nwtst1), v(l), v(nwtst1))
5857 call vvmulp(n, v(step1), v(nwtst1), d, 1)
5858 v(dst0) = v2norm(n, v(step1))
5859 call vvmulp(n, v(dg1), v(dg1), d, -1)
5860 call ltvmul(n, v(step1), v(l), v(dg1))
5861 v(gthg) = v2norm(n, v(step1))
5863 160 call dbdog(v(dg1), lv, n, v(nwtst1), v(step1), v)
5864 if (iv(irc) .eq. 6) go to 180
5866 ! *** check whether evaluating f(x0 + step) looks worthwhile ***
5868 if (v(dstnrm) .le. zero) go to 180
5869 if (iv(irc) .ne. 5) go to 170
5870 if (v(radfac) .le. one) go to 170
5871 if (v(preduc) .le. onep2 * v(fdif)) go to 180
5873 ! *** compute f(x0 + step) ***
5877 call vaxpy(n, x, one, v(step1), v(x01))
5878 iv(nfcall) = iv(nfcall) + 1
5883 !. . . . . . . . . . . . . assess candidate step . . . . . . . . . . .
5886 v(reldx) = reldst(n, d, x, v(x01))
5887 call assst(iv, liv, lv, v)
5890 if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
5891 if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
5892 if (iv(restor) .ne. 3) go to 190
5893 call vcopy(n, v(step1), v(lstgst))
5894 call vaxpy(n, x, one, v(step1), v(x01))
5895 v(reldx) = reldst(n, d, x, v(x01))
5898 go to (200,230,230,230,200,210,220,220,220,220,220,220,280,250), k
5900 ! *** recompute step with changed radius ***
5902 200 v(radius) = v(radfac) * v(dstnrm)
5905 ! *** compute step of length v(lmaxs) for singular convergence test.
5907 210 v(radius) = v(lmaxs)
5910 ! *** convergence or false convergence ***
5912 220 iv(cnvcod) = k - 4
5913 if (v(f) .ge. v(f0)) go to 290
5914 if (iv(xirc) .eq. 14) go to 290
5917 !. . . . . . . . . . . . process acceptable step . . . . . . . . . . .
5919 230 if (iv(irc) .ne. 3) go to 240
5923 ! *** set temp1 = hessian * step for use in gradient tests ***
5926 call ltvmul(n, v(temp1), v(l), v(step1))
5927 call lvmul(n, v(temp1), v(l), v(temp1))
5929 ! *** compute gradient ***
5931 240 iv(ngcall) = iv(ngcall) + 1
5935 ! *** initializations -- g0 = g - g0, etc. ***
5938 call vaxpy(n, v(g01), negone, v(g01), g)
5941 if (iv(irc) .ne. 3) go to 270
5943 ! *** set v(radfac) by gradient tests ***
5945 ! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) ***
5947 call vaxpy(n, v(temp1), negone, v(g01), v(temp1))
5948 call vvmulp(n, v(temp1), v(temp1), d, -1)
5950 ! *** do gradient tests ***
5952 if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) &
5954 if (dotprd(n, g, v(step1)) &
5955 .ge. v(gtstep) * v(tuner5)) go to 270
5956 260 v(radfac) = v(incfac)
5958 ! *** update h, loop ***
5963 call wzbfgs(v(l), n, v(step1), v(w), v(g01), v(z))
5965 ! ** use the n-vectors starting at v(step1) and v(g01) for scratch..
5966 call lupdat(v(temp1), v(step1), v(l), v(g01), v(l), n, v(w), v(z))
5970 !. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . .
5972 ! *** bad parameters to assess ***
5977 ! *** print summary of final iteration and other requested items ***
5979 290 iv(1) = iv(cnvcod)
5981 300 call itsum(d, g, iv, liv, lv, n, v, x)
5985 ! *** last line of sumit follows ***
5986 end subroutine sumit
5987 !-----------------------------------------------------------------------------
5988 subroutine dbdog(dig,lv,n,nwtstp,step,v)
5990 ! *** compute double dogleg step ***
5992 ! *** parameter declarations ***
5995 real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv)
5999 ! this subroutine computes a candidate step (for use in an uncon-
6000 ! strained minimization code) by the double dogleg algorithm of
6001 ! dennis and mei (ref. 1), which is a variation on powell*s dogleg
6002 ! scheme (ref. 2, p. 95).
6004 !-------------------------- parameter usage --------------------------
6006 ! dig (input) diag(d)**-2 * g -- see algorithm notes.
6007 ! g (input) the current gradient vector.
6008 ! lv (input) length of v.
6009 ! n (input) number of components in dig, g, nwtstp, and step.
6010 ! nwtstp (input) negative newton step -- see algorithm notes.
6011 ! step (output) the computed step.
6012 ! v (i/o) values array, the following components of which are
6014 ! v(bias) (input) bias for relaxed newton step, which is v(bias) of
6015 ! the way from the full newton to the fully relaxed newton
6016 ! step. recommended value = 0.8 .
6017 ! v(dgnorm) (input) 2-norm of diag(d)**-1 * g -- see algorithm notes.
6018 ! v(dstnrm) (output) 2-norm of diag(d) * step, which is v(radius)
6019 ! unless v(stppar) = 0 -- see algorithm notes.
6020 ! v(dst0) (input) 2-norm of diag(d) * nwtstp -- see algorithm notes.
6021 ! v(grdfac) (output) the coefficient of dig in the step returned --
6022 ! step(i) = v(grdfac)*dig(i) + v(nwtfac)*nwtstp(i).
6023 ! v(gthg) (input) square-root of (dig**t) * (hessian) * dig -- see
6025 ! v(gtstep) (output) inner product between g and step.
6026 ! v(nreduc) (output) function reduction predicted for the full newton
6028 ! v(nwtfac) (output) the coefficient of nwtstp in the step returned --
6029 ! see v(grdfac) above.
6030 ! v(preduc) (output) function reduction predicted for the step returned.
6031 ! v(radius) (input) the trust region radius. d times the step returned
6032 ! has 2-norm v(radius) unless v(stppar) = 0.
6033 ! v(stppar) (output) code telling how step was computed... 0 means a
6034 ! full newton step. between 0 and 1 means v(stppar) of the
6035 ! way from the newton to the relaxed newton step. between
6036 ! 1 and 2 means a true double dogleg step, v(stppar) - 1 of
6037 ! the way from the relaxed newton to the cauchy step.
6038 ! greater than 2 means 1 / (v(stppar) - 1) times the cauchy
6041 !------------------------------- notes -------------------------------
6043 ! *** algorithm notes ***
6045 ! let g and h be the current gradient and hessian approxima-
6046 ! tion respectively and let d be the current scale vector. this
6047 ! routine assumes dig = diag(d)**-2 * g and nwtstp = h**-1 * g.
6048 ! the step computed is the same one would get by replacing g and h
6049 ! by diag(d)**-1 * g and diag(d)**-1 * h * diag(d)**-1,
6050 ! computing step, and translating step back to the original
6051 ! variables, i.e., premultiplying it by diag(d)**-1.
6053 ! *** references ***
6055 ! 1. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti-
6056 ! mization algorithms which use function and gradient
6057 ! values, j. optim. theory applic. 28, pp. 453-482.
6058 ! 2. powell, m.j.d. (1970), a hybrid method for non-linear equations,
6059 ! in numerical methods for non-linear equations, edited by
6060 ! p. rabinowitz, gordon and breach, london.
6064 ! coded by david m. gay.
6065 ! this subroutine was written in connection with research supported
6066 ! by the national science foundation under grants mcs-7600324 and
6069 !------------------------ external quantities ------------------------
6071 ! *** functions and subroutines called ***
6073 !el external dotprd, v2norm
6074 !el real(kind=8) :: dotprd, v2norm
6076 ! dotprd... returns inner product of two vectors.
6077 ! v2norm... returns 2-norm of a vector.
6079 ! *** intrinsic functions ***
6081 !el real(kind=8) :: dsqrt
6083 !-------------------------- local variables --------------------------
6086 real(kind=8) :: cfact, cnorm, ctrnwt, ghinvg, femnsq, gnorm,&
6087 nwtnrm, relax, rlambd, t, t1, t2
6088 !el real(kind=8) :: half, one, two, zero
6090 ! *** v subscripts ***
6092 !el integer bias, dgnorm, dstnrm, dst0, grdfac, gthg, gtstep,
6093 !el 1 nreduc, nwtfac, preduc, radius, stppar
6095 ! *** data initializations ***
6098 ! data half/0.5d+0/, one/1.d+0/, two/2.d+0/, zero/0.d+0/
6100 real(kind=8),parameter :: half=0.5d+0, one=1.d+0, two=2.d+0, zero=0.d+0
6104 ! data bias/43/, dgnorm/1/, dstnrm/2/, dst0/3/, grdfac/45/,
6105 ! 1 gthg/44/, gtstep/4/, nreduc/6/, nwtfac/46/, preduc/7/,
6106 ! 2 radius/8/, stppar/5/
6108 integer,parameter :: bias=43, dgnorm=1, dstnrm=2, dst0=3, grdfac=45,&
6109 gthg=44, gtstep=4, nreduc=6, nwtfac=46, preduc=7,&
6113 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6117 if (nwtnrm .gt. zero) rlambd = v(radius) / nwtnrm
6119 ghinvg = two * v(nreduc)
6122 if (rlambd .lt. one) go to 30
6124 ! *** the newton step is inside the trust region ***
6129 v(preduc) = v(nreduc)
6132 20 step(i) = -nwtstp(i)
6135 30 v(dstnrm) = v(radius)
6136 cfact = (gnorm / v(gthg))**2
6137 ! *** cauchy step = -cfact * g.
6138 cnorm = gnorm * cfact
6139 relax = one - v(bias) * (one - gnorm*cnorm/ghinvg)
6140 if (rlambd .lt. relax) go to 50
6142 ! *** step is between relaxed newton and full newton steps ***
6144 v(stppar) = one - (rlambd - relax) / (one - relax)
6146 v(gtstep) = t * ghinvg
6147 v(preduc) = rlambd * (one - half*rlambd) * ghinvg
6150 40 step(i) = t * nwtstp(i)
6153 50 if (cnorm .lt. v(radius)) go to 70
6155 ! *** the cauchy step lies outside the trust region --
6156 ! *** step = scaled cauchy step ***
6158 t = -v(radius) / gnorm
6160 v(stppar) = one + cnorm / v(radius)
6161 v(gtstep) = -v(radius) * gnorm
6162 v(preduc) = v(radius)*(gnorm - half*v(radius)*(v(gthg)/gnorm)**2)
6164 60 step(i) = t * dig(i)
6167 ! *** compute dogleg step between cauchy and relaxed newton ***
6168 ! *** femur = relaxed newton step minus cauchy step ***
6170 70 ctrnwt = cfact * relax * ghinvg / gnorm
6171 ! *** ctrnwt = inner prod. of cauchy and relaxed newton steps,
6172 ! *** scaled by gnorm**-1.
6173 t1 = ctrnwt - gnorm*cfact**2
6174 ! *** t1 = inner prod. of femur and cauchy step, scaled by
6176 t2 = v(radius)*(v(radius)/gnorm) - gnorm*cfact**2
6178 femnsq = (t/gnorm)*t - ctrnwt - t1
6179 ! *** femnsq = square of 2-norm of femur, scaled by gnorm**-1.
6180 t = t2 / (t1 + dsqrt(t1**2 + femnsq*t2))
6181 ! *** dogleg step = cauchy step + t * femur.
6182 t1 = (t - one) * cfact
6187 v(gtstep) = t1*gnorm**2 + t2*ghinvg
6188 v(preduc) = -t1*gnorm * ((t2 + one)*gnorm) &
6189 - t2 * (one + half*t2)*ghinvg &
6190 - half * (v(gthg)*t1)**2
6192 80 step(i) = t1*dig(i) + t2*nwtstp(i)
6195 ! *** last line of dbdog follows ***
6196 end subroutine dbdog
6197 !-----------------------------------------------------------------------------
6198 subroutine ltvmul(n,x,l,y)
6200 ! *** compute x = (l**t)*y, where l is an n x n lower
6201 ! *** triangular matrix stored compactly by rows. x and y may
6202 ! *** occupy the same storage. ***
6205 !al real(kind=8) :: x(n), l(1), y(n)
6206 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
6207 ! dimension l(n*(n+1)/2)
6208 integer :: i, ij, i0, j
6209 real(kind=8) :: yi !el, zero
6213 real(kind=8),parameter :: zero=0.d+0
6222 x(j) = x(j) + yi*l(ij)
6227 ! *** last card of ltvmul follows ***
6228 end subroutine ltvmul
6229 !-----------------------------------------------------------------------------
6230 subroutine lupdat(beta,gamma,l,lambda,lplus,n,w,z)
6232 ! *** compute lplus = secant update of l ***
6234 ! *** parameter declarations ***
6237 !al double precision beta(n), gamma(n), l(1), lambda(n), lplus(1),
6238 real(kind=8) :: beta(n), gamma(n), l(n*(n+1)/2), lambda(n), &
6239 lplus(n*(n+1)/2),w(n), z(n)
6240 ! dimension l(n*(n+1)/2), lplus(n*(n+1)/2)
6242 !-------------------------- parameter usage --------------------------
6244 ! beta = scratch vector.
6245 ! gamma = scratch vector.
6246 ! l (input) lower triangular matrix, stored rowwise.
6247 ! lambda = scratch vector.
6248 ! lplus (output) lower triangular matrix, stored rowwise, which may
6249 ! occupy the same storage as l.
6250 ! n (input) length of vector parameters and order of matrices.
6251 ! w (input, destroyed on output) right singular vector of rank 1
6253 ! z (input, destroyed on output) left singular vector of rank 1
6256 !------------------------------- notes -------------------------------
6258 ! *** application and usage restrictions ***
6260 ! this routine updates the cholesky factor l of a symmetric
6261 ! positive definite matrix to which a secant update is being
6262 ! applied -- it computes a cholesky factor lplus of
6263 ! l * (i + z*w**t) * (i + w*z**t) * l**t. it is assumed that w
6264 ! and z have been chosen so that the updated matrix is strictly
6265 ! positive definite.
6267 ! *** algorithm notes ***
6269 ! this code uses recurrence 3 of ref. 1 (with d(j) = 1 for all j)
6270 ! to compute lplus of the form l * (i + z*w**t) * q, where q
6271 ! is an orthogonal matrix that makes the result lower triangular.
6272 ! lplus may have some negative diagonal elements.
6274 ! *** references ***
6276 ! 1. goldfarb, d. (1976), factorized variable metric methods for uncon-
6277 ! strained optimization, math. comput. 30, pp. 796-811.
6281 ! coded by david m. gay (fall 1979).
6282 ! this subroutine was written in connection with research supported
6283 ! by the national science foundation under grants mcs-7600324 and
6286 !------------------------ external quantities ------------------------
6288 ! *** intrinsic functions ***
6290 !el real(kind=8) :: dsqrt
6292 !-------------------------- local variables --------------------------
6294 integer :: i, ij, j, jj, jp1, k, nm1, np1
6295 real(kind=8) :: a, b, bj, eta, gj, lj, lij, ljj, nu, s, theta,&
6297 !el real(kind=8) :: one, zero
6299 ! *** data initializations ***
6302 ! data one/1.d+0/, zero/0.d+0/
6304 real(kind=8),parameter :: one=1.d+0, zero=0.d+0
6307 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6311 if (n .le. 1) go to 30
6314 ! *** temporarily store s(j) = sum over k = j+1 to n of w(k)**2 in
6324 ! *** compute lambda, gamma, and beta by goldfarb*s recurrence 3.
6328 a = nu*z(j) - eta*wj
6331 lj = dsqrt(theta**2 + a*s)
6332 if (theta .gt. zero) lj = -lj
6335 gamma(j) = b * nu / lj
6336 beta(j) = (a - b*eta) / lj
6338 eta = -(eta + (a**2)/(theta - lj)) / lj
6340 30 lambda(n) = one + (nu*z(n) - eta*w(n))*w(n)
6342 ! *** update l, gradually overwriting w and z with l*w and l*z.
6345 jj = n * (n + 1) / 2
6350 lplus(jj) = lj * ljj
6355 if (k .eq. 1) go to 50
6362 lplus(ij) = lj*lij + bj*w(i) + gj*z(i)
6363 w(i) = w(i) + lij*wj
6364 z(i) = z(i) + lij*zj
6371 ! *** last card of lupdat follows ***
6372 end subroutine lupdat
6373 !-----------------------------------------------------------------------------
6374 subroutine lvmul(n,x,l,y)
6376 ! *** compute x = l*y, where l is an n x n lower triangular
6377 ! *** matrix stored compactly by rows. x and y may occupy the same
6381 !al double precision x(n), l(1), y(n)
6382 real(kind=8) :: x(n), l(n*(n+1)/2), y(n)
6383 ! dimension l(n*(n+1)/2)
6384 integer :: i, ii, ij, i0, j, np1
6385 real(kind=8) :: t !el, zero
6389 real(kind=8),parameter :: zero=0.d+0
6405 ! *** last card of lvmul follows ***
6406 end subroutine lvmul
6407 !-----------------------------------------------------------------------------
6408 subroutine vvmulp(n,x,y,z,k)
6410 ! *** set x(i) = y(i) * z(i)**k, 1 .le. i .le. n (for k = 1 or -1) ***
6413 real(kind=8) :: x(n), y(n), z(n)
6416 if (k .ge. 0) go to 20
6418 10 x(i) = y(i) / z(i)
6422 30 x(i) = y(i) * z(i)
6424 ! *** last card of vvmulp follows ***
6425 end subroutine vvmulp
6426 !-----------------------------------------------------------------------------
6427 subroutine wzbfgs(l,n,s,w,y,z)
6429 ! *** compute y and z for lupdat corresponding to bfgs update.
6432 !al double precision l(1), s(n), w(n), y(n), z(n)
6433 real(kind=8) :: l(n*(n+1)/2), s(n), w(n), y(n), z(n)
6434 ! dimension l(n*(n+1)/2)
6436 !-------------------------- parameter usage --------------------------
6438 ! l (i/o) cholesky factor of hessian, a lower triang. matrix stored
6439 ! compactly by rows.
6440 ! n (input) order of l and length of s, w, y, z.
6441 ! s (input) the step just taken.
6442 ! w (output) right singular vector of rank 1 correction to l.
6443 ! y (input) change in gradients corresponding to s.
6444 ! z (output) left singular vector of rank 1 correction to l.
6446 !------------------------------- notes -------------------------------
6448 ! *** algorithm notes ***
6450 ! when s is computed in certain ways, e.g. by gqtstp or
6451 ! dbldog, it is possible to save n**2/2 operations since (l**t)*s
6452 ! or l*(l**t)*s is then known.
6453 ! if the bfgs update to l*(l**t) would reduce its determinant to
6454 ! less than eps times its old value, then this routine in effect
6455 ! replaces y by theta*y + (1 - theta)*l*(l**t)*s, where theta
6456 ! (between 0 and 1) is chosen to make the reduction factor = eps.
6460 ! coded by david m. gay (fall 1979).
6461 ! this subroutine was written in connection with research supported
6462 ! by the national science foundation under grants mcs-7600324 and
6465 !------------------------ external quantities ------------------------
6467 ! *** functions and subroutines called ***
6469 !el external dotprd, livmul, ltvmul
6470 !el real(kind=8) :: dotprd
6471 ! dotprd returns inner product of two vectors.
6472 ! livmul multiplies l**-1 times a vector.
6473 ! ltvmul multiplies l**t times a vector.
6475 ! *** intrinsic functions ***
6477 !el real(kind=8) :: dsqrt
6479 !-------------------------- local variables --------------------------
6482 real(kind=8) :: cs, cy, epsrt, shs, ys, theta !el, eps, one
6484 ! *** data initializations ***
6487 ! data eps/0.1d+0/, one/1.d+0/
6489 real(kind=8),parameter :: eps=0.1d+0, one=1.d+0
6492 !+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++
6494 call ltvmul(n, w, l, s)
6495 shs = dotprd(n, w, w)
6496 ys = dotprd(n, y, s)
6497 if (ys .ge. eps*shs) go to 10
6498 theta = (one - eps) * shs / (shs - ys)
6500 cy = theta / (shs * epsrt)
6501 cs = (one + (theta-one)/epsrt) / shs
6503 10 cy = one / (dsqrt(ys) * dsqrt(shs))
6505 20 call livmul(n, z, l, y)
6507 30 z(i) = cy * z(i) - cs * w(i)
6510 ! *** last card of wzbfgs follows ***
6511 end subroutine wzbfgs
6512 !-----------------------------------------------------------------------------
6513 !-----------------------------------------------------------------------------