X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fminim.F90;fp=source%2Funres%2Fminim.F90;h=43056404767f7122e86f0098a611a1d5492b8c6f;hb=2611e37f82b576a1366c2d78fce87c1a55852037;hp=0000000000000000000000000000000000000000;hpb=a1e9d5a6b704c90ebd10f86a655780212a880dce;p=unres4.git diff --git a/source/unres/minim.F90 b/source/unres/minim.F90 new file mode 100644 index 0000000..4305640 --- /dev/null +++ b/source/unres/minim.F90 @@ -0,0 +1,6508 @@ + module minimm +!----------------------------------------------------------------------------- + use io_units + use names + use math +! use MPI_data + use geometry_data + use energy_data + use control_data + use minim_data + use geometry +! use csa_data +! use energy + implicit none +!----------------------------------------------------------------------------- +! +! +!----------------------------------------------------------------------------- + contains +!----------------------------------------------------------------------------- +! cored.f +!----------------------------------------------------------------------------- + subroutine assst(iv, liv, lv, v) +! +! *** assess candidate step (***sol version 2.3) *** +! + integer :: liv, l,lv + integer :: iv(liv) + real(kind=8) :: v(lv) +! +! *** purpose *** +! +! this subroutine is called by an unconstrained minimization +! routine to assess the next candidate step. it may recommend one +! of several courses of action, such as accepting the step, recom- +! puting it using the same or a new quadratic model, or halting due +! to convergence or false convergence. see the return code listing +! below. +! +!-------------------------- parameter usage -------------------------- +! +! iv (i/o) integer parameter and scratch vector -- see description +! below of iv values referenced. +! liv (in) length of iv array. +! lv (in) length of v array. +! v (i/o) real parameter and scratch vector -- see description +! below of v values referenced. +! +! *** iv values referenced *** +! +! iv(irc) (i/o) on input for the first step tried in a new iteration, +! iv(irc) should be set to 3 or 4 (the value to which it is +! set when step is definitely to be accepted). on input +! after step has been recomputed, iv(irc) should be +! unchanged since the previous return of assst. +! on output, iv(irc) is a return code having one of the +! following values... +! 1 = switch models or try smaller step. +! 2 = switch models or accept step. +! 3 = accept step and determine v(radfac) by gradient +! tests. +! 4 = accept step, v(radfac) has been determined. +! 5 = recompute step (using the same model). +! 6 = recompute step with radius = v(lmaxs) but do not +! evaulate the objective function. +! 7 = x-convergence (see v(xctol)). +! 8 = relative function convergence (see v(rfctol)). +! 9 = both x- and relative function convergence. +! 10 = absolute function convergence (see v(afctol)). +! 11 = singular convergence (see v(lmaxs)). +! 12 = false convergence (see v(xftol)). +! 13 = iv(irc) was out of range on input. +! return code i has precdence over i+1 for i = 9, 10, 11. +! iv(mlstgd) (i/o) saved value of iv(model). +! iv(model) (i/o) on input, iv(model) should be an integer identifying +! the current quadratic model of the objective function. +! if a previous step yielded a better function reduction, +! then iv(model) will be set to iv(mlstgd) on output. +! iv(nfcall) (in) invocation count for the objective function. +! iv(nfgcal) (i/o) value of iv(nfcall) at step that gave the biggest +! function reduction this iteration. iv(nfgcal) remains +! unchanged until a function reduction is obtained. +! iv(radinc) (i/o) the number of radius increases (or minus the number +! of decreases) so far this iteration. +! iv(restor) (out) set to 1 if v(f) has been restored and x should be +! restored to its initial value, to 2 if x should be saved, +! to 3 if x should be restored from the saved value, and to +! 0 otherwise. +! iv(stage) (i/o) count of the number of models tried so far in the +! current iteration. +! iv(stglim) (in) maximum number of models to consider. +! iv(switch) (out) set to 0 unless a new model is being tried and it +! gives a smaller function value than the previous model, +! in which case assst sets iv(switch) = 1. +! iv(toobig) (in) is nonzero if step was too big (e.g. if it caused +! overflow). +! iv(xirc) (i/o) value that iv(irc) would have in the absence of +! convergence, false convergence, and oversized steps. +! +! *** v values referenced *** +! +! v(afctol) (in) absolute function convergence tolerance. if the +! absolute value of the current function value v(f) is less +! than v(afctol), then assst returns with iv(irc) = 10. +! v(decfac) (in) factor by which to decrease radius when iv(toobig) is +! nonzero. +! v(dstnrm) (in) the 2-norm of d*step. +! v(dstsav) (i/o) value of v(dstnrm) on saved step. +! v(dst0) (in) the 2-norm of d times the newton step (when defined, +! i.e., for v(nreduc) .ge. 0). +! v(f) (i/o) on both input and output, v(f) is the objective func- +! tion value at x. if x is restored to a previous value, +! then v(f) is restored to the corresponding value. +! v(fdif) (out) the function reduction v(f0) - v(f) (for the output +! value of v(f) if an earlier step gave a bigger function +! decrease, and for the input value of v(f) otherwise). +! v(flstgd) (i/o) saved value of v(f). +! v(f0) (in) objective function value at start of iteration. +! v(gtslst) (i/o) value of v(gtstep) on saved step. +! v(gtstep) (in) inner product between step and gradient. +! v(incfac) (in) minimum factor by which to increase radius. +! v(lmaxs) (in) maximum reasonable step size (and initial step bound). +! if the actual function decrease is no more than twice +! what was predicted, if a return with iv(irc) = 7, 8, 9, +! or 10 does not occur, if v(dstnrm) .gt. v(lmaxs), and if +! v(preduc) .le. v(sctol) * abs(v(f0)), then assst re- +! turns with iv(irc) = 11. if so doing appears worthwhile, +! then assst repeats this test with v(preduc) computed for +! a step of length v(lmaxs) (by a return with iv(irc) = 6). +! v(nreduc) (i/o) function reduction predicted by quadratic model for +! newton step. if assst is called with iv(irc) = 6, i.e., +! if v(preduc) has been computed with radius = v(lmaxs) for +! use in the singular convervence test, then v(nreduc) is +! set to -v(preduc) before the latter is restored. +! v(plstgd) (i/o) value of v(preduc) on saved step. +! v(preduc) (i/o) function reduction predicted by quadratic model for +! current step. +! v(radfac) (out) factor to be used in determining the new radius, +! which should be v(radfac)*dst, where dst is either the +! output value of v(dstnrm) or the 2-norm of +! diag(newd)*step for the output value of step and the +! updated version, newd, of the scale vector d. for +! iv(irc) = 3, v(radfac) = 1.0 is returned. +! v(rdfcmn) (in) minimum value for v(radfac) in terms of the input +! value of v(dstnrm) -- suggested value = 0.1. +! v(rdfcmx) (in) maximum value for v(radfac) -- suggested value = 4.0. +! v(reldx) (in) scaled relative change in x caused by step, computed +! (e.g.) by function reldst as +! max (d(i)*abs(x(i)-x0(i)), 1 .le. i .le. p) / +! max (d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p). +! v(rfctol) (in) relative function convergence tolerance. if the +! actual function reduction is at most twice what was pre- +! dicted and v(nreduc) .le. v(rfctol)*abs(v(f0)), then +! assst returns with iv(irc) = 8 or 9. +! v(stppar) (in) marquardt parameter -- 0 means full newton step. +! v(tuner1) (in) tuning constant used to decide if the function +! reduction was much less than expected. suggested +! value = 0.1. +! v(tuner2) (in) tuning constant used to decide if the function +! reduction was large enough to accept step. suggested +! value = 10**-4. +! v(tuner3) (in) tuning constant used to decide if the radius +! should be increased. suggested value = 0.75. +! v(xctol) (in) x-convergence criterion. if step is a newton step +! (v(stppar) = 0) having v(reldx) .le. v(xctol) and giving +! at most twice the predicted function decrease, then +! assst returns iv(irc) = 7 or 9. +! v(xftol) (in) false convergence tolerance. if step gave no or only +! a small function decrease and v(reldx) .le. v(xftol), +! then assst returns with iv(irc) = 12. +! +!------------------------------- notes ------------------------------- +! +! *** application and usage restrictions *** +! +! this routine is called as part of the nl2sol (nonlinear +! least-squares) package. it may be used in any unconstrained +! minimization solver that uses dogleg, goldfeld-quandt-trotter, +! or levenberg-marquardt steps. +! +! *** algorithm notes *** +! +! see (1) for further discussion of the assessing and model +! switching strategies. while nl2sol considers only two models, +! assst is designed to handle any number of models. +! +! *** usage notes *** +! +! on the first call of an iteration, only the i/o variables +! step, x, iv(irc), iv(model), v(f), v(dstnrm), v(gtstep), and +! v(preduc) need have been initialized. between calls, no i/o +! values execpt step, x, iv(model), v(f) and the stopping toler- +! ances should be changed. +! after a return for convergence or false convergence, one can +! change the stopping tolerances and call assst again, in which +! case the stopping tests will be repeated. +! +! *** references *** +! +! (1) dennis, j.e., jr., gay, d.m., and welsch, r.e. (1981), +! an adaptive nonlinear least-squares algorithm, +! acm trans. math. software, vol. 7, no. 3. +! +! (2) powell, m.j.d. (1970) a fortran subroutine for solving +! systems of nonlinear algebraic equations, in numerical +! methods for nonlinear algebraic equations, edited by +! p. rabinowitz, gordon and breach, london. +! +! *** history *** +! +! john dennis designed much of this routine, starting with +! ideas in (2). roy welsch suggested the model switching strategy. +! david gay and stephen peters cast this subroutine into a more +! portable form (winter 1977), and david gay cast it into its +! present form (fall 1978). +! +! *** general *** +! +! this subroutine was written in connection with research +! supported by the national science foundation under grants +! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and +! mcs-7906671. +! +!------------------------ external quantities ------------------------ +! +! *** no external functions and subroutines *** +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dabs, dmax1 +!/ +! *** no common blocks *** +! +!-------------------------- local variables -------------------------- +! + logical :: goodx + integer :: i, nfc + real(kind=8) :: emax, emaxs, gts, rfac1, xmax +!el real(kind=8) :: half, one, onep2, two, zero +! +! *** subscripts for iv and v *** +! +!el integer :: afctol, decfac, dstnrm, dstsav, dst0, f, fdif, flstgd, f0,& +!el gtslst, gtstep, incfac, irc, lmaxs, mlstgd, model, nfcall,& +!el nfgcal, nreduc, plstgd, preduc, radfac, radinc, rdfcmn,& +!el rdfcmx, reldx, restor, rfctol, sctol, stage, stglim,& +!el stppar, switch, toobig, tuner1, tuner2, tuner3, xctol,& +!el xftol, xirc +! +! +! *** data initializations *** +! +!/6 +! data half/0.5d+0/, one/1.d+0/, onep2/1.2d+0/, two/2.d+0/, +! 1 zero/0.d+0/ +!/7 + real(kind=8),parameter :: half=0.5d+0, one=1.d+0, onep2=1.2d+0, two=2.d+0,& + zero=0.d+0 +!/ +! +!/6 +! data irc/29/, mlstgd/32/, model/5/, nfcall/6/, nfgcal/7/, +! 1 radinc/8/, restor/9/, stage/10/, stglim/11/, switch/12/, +! 2 toobig/2/, xirc/13/ +!/7 + integer,parameter :: irc=29, mlstgd=32, model=5, nfcall=6, nfgcal=7,& + radinc=8, restor=9, stage=10, stglim=11, switch=12,& + toobig=2, xirc=13 +!/ +!/6 +! data afctol/31/, decfac/22/, dstnrm/2/, dst0/3/, dstsav/18/, +! 1 f/10/, fdif/11/, flstgd/12/, f0/13/, gtslst/14/, gtstep/4/, +! 2 incfac/23/, lmaxs/36/, nreduc/6/, plstgd/15/, preduc/7/, +! 3 radfac/16/, rdfcmn/24/, rdfcmx/25/, reldx/17/, rfctol/32/, +! 4 sctol/37/, stppar/5/, tuner1/26/, tuner2/27/, tuner3/28/, +! 5 xctol/33/, xftol/34/ +!/7 + integer,parameter :: afctol=31, decfac=22, dstnrm=2, dst0=3, dstsav=18,& + f=10, fdif=11, flstgd=12, f0=13, gtslst=14, gtstep=4,& + incfac=23, lmaxs=36, nreduc=6, plstgd=15, preduc=7,& + radfac=16, rdfcmn=24, rdfcmx=25, reldx=17, rfctol=32,& + sctol=37, stppar=5, tuner1=26, tuner2=27, tuner3=28,& + xctol=33, xftol=34 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! + nfc = iv(nfcall) + iv(switch) = 0 + iv(restor) = 0 + rfac1 = one + goodx = .true. + i = iv(irc) + if (i .ge. 1 .and. i .le. 12) & + go to (20,30,10,10,40,280,220,220,220,220,220,170), i + iv(irc) = 13 + go to 999 +! +! *** initialize for new iteration *** +! + 10 iv(stage) = 1 + iv(radinc) = 0 + v(flstgd) = v(f0) + if (iv(toobig) .eq. 0) go to 110 + iv(stage) = -1 + iv(xirc) = i + go to 60 +! +! *** step was recomputed with new model or smaller radius *** +! *** first decide which *** +! + 20 if (iv(model) .ne. iv(mlstgd)) go to 30 +! *** old model retained, smaller radius tried *** +! *** do not consider any more new models this iteration *** + iv(stage) = iv(stglim) + iv(radinc) = -1 + go to 110 +! +! *** a new model is being tried. decide whether to keep it. *** +! + 30 iv(stage) = iv(stage) + 1 +! +! *** now we add the possibiltiy that step was recomputed with *** +! *** the same model, perhaps because of an oversized step. *** +! + 40 if (iv(stage) .gt. 0) go to 50 +! +! *** step was recomputed because it was too big. *** +! + if (iv(toobig) .ne. 0) go to 60 +! +! *** restore iv(stage) and pick up where we left off. *** +! + iv(stage) = -iv(stage) + i = iv(xirc) + go to (20, 30, 110, 110, 70), i +! + 50 if (iv(toobig) .eq. 0) go to 70 +! +! *** handle oversize step *** +! + if (iv(radinc) .gt. 0) go to 80 + iv(stage) = -iv(stage) + iv(xirc) = iv(irc) +! + 60 v(radfac) = v(decfac) + iv(radinc) = iv(radinc) - 1 + iv(irc) = 5 + iv(restor) = 1 + go to 999 +! + 70 if (v(f) .lt. v(flstgd)) go to 110 +! +! *** the new step is a loser. restore old model. *** +! + if (iv(model) .eq. iv(mlstgd)) go to 80 + iv(model) = iv(mlstgd) + iv(switch) = 1 +! +! *** restore step, etc. only if a previous step decreased v(f). +! + 80 if (v(flstgd) .ge. v(f0)) go to 110 + iv(restor) = 1 + v(f) = v(flstgd) + v(preduc) = v(plstgd) + v(gtstep) = v(gtslst) + if (iv(switch) .eq. 0) rfac1 = v(dstnrm) / v(dstsav) + v(dstnrm) = v(dstsav) + nfc = iv(nfgcal) + goodx = .false. +! + 110 v(fdif) = v(f0) - v(f) + if (v(fdif) .gt. v(tuner2) * v(preduc)) go to 140 + if(iv(radinc).gt.0) go to 140 +! +! *** no (or only a trivial) function decrease +! *** -- so try new model or smaller radius +! + if (v(f) .lt. v(f0)) go to 120 + iv(mlstgd) = iv(model) + v(flstgd) = v(f) + v(f) = v(f0) + iv(restor) = 1 + go to 130 + 120 iv(nfgcal) = nfc + 130 iv(irc) = 1 + if (iv(stage) .lt. iv(stglim)) go to 160 + iv(irc) = 5 + iv(radinc) = iv(radinc) - 1 + go to 160 +! +! *** nontrivial function decrease achieved *** +! + 140 iv(nfgcal) = nfc + rfac1 = one + v(dstsav) = v(dstnrm) + if (v(fdif) .gt. v(preduc)*v(tuner1)) go to 190 +! +! *** decrease was much less than predicted -- either change models +! *** or accept step with decreased radius. +! + if (iv(stage) .ge. iv(stglim)) go to 150 +! *** consider switching models *** + iv(irc) = 2 + go to 160 +! +! *** accept step with decreased radius *** +! + 150 iv(irc) = 4 +! +! *** set v(radfac) to fletcher*s decrease factor *** +! + 160 iv(xirc) = iv(irc) + emax = v(gtstep) + v(fdif) + v(radfac) = half * rfac1 + if (emax .lt. v(gtstep)) v(radfac) = rfac1 * dmax1(v(rdfcmn),& + half * v(gtstep)/emax) +! +! *** do false convergence test *** +! + 170 if (v(reldx) .le. v(xftol)) go to 180 + iv(irc) = iv(xirc) + if (v(f) .lt. v(f0)) go to 200 + go to 230 +! + 180 iv(irc) = 12 + go to 240 +! +! *** handle good function decrease *** +! + 190 if (v(fdif) .lt. (-v(tuner3) * v(gtstep))) go to 210 +! +! *** increasing radius looks worthwhile. see if we just +! *** recomputed step with a decreased radius or restored step +! *** after recomputing it with a larger radius. +! + if (iv(radinc) .lt. 0) go to 210 + if (iv(restor) .eq. 1) go to 210 +! +! *** we did not. try a longer step unless this was a newton +! *** step. + + v(radfac) = v(rdfcmx) + gts = v(gtstep) + if (v(fdif) .lt. (half/v(radfac) - one) * gts) & + v(radfac) = dmax1(v(incfac), half*gts/(gts + v(fdif))) + iv(irc) = 4 + if (v(stppar) .eq. zero) go to 230 + if (v(dst0) .ge. zero .and. (v(dst0) .lt. two*v(dstnrm) & + .or. v(nreduc) .lt. onep2*v(fdif))) go to 230 +! *** step was not a newton step. recompute it with +! *** a larger radius. + iv(irc) = 5 + iv(radinc) = iv(radinc) + 1 +! +! *** save values corresponding to good step *** +! + 200 v(flstgd) = v(f) + iv(mlstgd) = iv(model) + if (iv(restor) .ne. 1) iv(restor) = 2 + v(dstsav) = v(dstnrm) + iv(nfgcal) = nfc + v(plstgd) = v(preduc) + v(gtslst) = v(gtstep) + go to 230 +! +! *** accept step with radius unchanged *** +! + 210 v(radfac) = one + iv(irc) = 3 + go to 230 +! +! *** come here for a restart after convergence *** +! + 220 iv(irc) = iv(xirc) + if (v(dstsav) .ge. zero) go to 240 + iv(irc) = 12 + go to 240 +! +! *** perform convergence tests *** +! + 230 iv(xirc) = iv(irc) + 240 if (iv(restor) .eq. 1 .and. v(flstgd) .lt. v(f0)) iv(restor) = 3 + if (half * v(fdif) .gt. v(preduc)) go to 999 + emax = v(rfctol) * dabs(v(f0)) + emaxs = v(sctol) * dabs(v(f0)) + if (v(dstnrm) .gt. v(lmaxs) .and. v(preduc) .le. emaxs) & + iv(irc) = 11 + if (v(dst0) .lt. zero) go to 250 + i = 0 + if ((v(nreduc) .gt. zero .and. v(nreduc) .le. emax) .or. & + (v(nreduc) .eq. zero .and. v(preduc) .eq. zero)) i = 2 + if (v(stppar) .eq. zero .and. v(reldx) .le. v(xctol) & + .and. goodx) i = i + 1 + if (i .gt. 0) iv(irc) = i + 6 +! +! *** consider recomputing step of length v(lmaxs) for singular +! *** convergence test. +! + 250 if (iv(irc) .gt. 5 .and. iv(irc) .ne. 12) go to 999 + if (v(dstnrm) .gt. v(lmaxs)) go to 260 + if (v(preduc) .ge. emaxs) go to 999 + if (v(dst0) .le. zero) go to 270 + if (half * v(dst0) .le. v(lmaxs)) go to 999 + go to 270 + 260 if (half * v(dstnrm) .le. v(lmaxs)) go to 999 + xmax = v(lmaxs) / v(dstnrm) + if (xmax * (two - xmax) * v(preduc) .ge. emaxs) go to 999 + 270 if (v(nreduc) .lt. zero) go to 290 +! +! *** recompute v(preduc) for use in singular convergence test *** +! + v(gtslst) = v(gtstep) + v(dstsav) = v(dstnrm) + if (iv(irc) .eq. 12) v(dstsav) = -v(dstsav) + v(plstgd) = v(preduc) + i = iv(restor) + iv(restor) = 2 + if (i .eq. 3) iv(restor) = 0 + iv(irc) = 6 + go to 999 +! +! *** perform singular convergence test with recomputed v(preduc) *** +! + 280 v(gtstep) = v(gtslst) + v(dstnrm) = dabs(v(dstsav)) + iv(irc) = iv(xirc) + if (v(dstsav) .le. zero) iv(irc) = 12 + v(nreduc) = -v(preduc) + v(preduc) = v(plstgd) + iv(restor) = 3 + 290 if (-v(nreduc) .le. v(sctol) * dabs(v(f0))) iv(irc) = 11 +! + 999 return +! +! *** last card of assst follows *** + end subroutine assst +!----------------------------------------------------------------------------- + subroutine deflt(alg, iv, liv, lv, v) +! +! *** supply ***sol (version 2.3) default values to iv and v *** +! +! *** alg = 1 means regression constants. +! *** alg = 2 means general unconstrained optimization constants. +! + integer :: liv, l,lv + integer :: alg, iv(liv) + real(kind=8) :: v(lv) +! +!el external imdcon, vdflt +!el integer imdcon +! imdcon... returns machine-dependent integer constants. +! vdflt.... provides default values to v. +! + integer :: miv, m + integer :: miniv(2), minv(2) +! +! *** subscripts for iv *** +! +!el integer algsav, covprt, covreq, dtype, hc, ierr, inith, inits, +!el 1 ipivot, ivneed, lastiv, lastv, lmat, mxfcal, mxiter, +!el 2 nfcov, ngcov, nvdflt, outlev, parprt, parsav, perm, +!el 3 prunit, qrtyp, rdreq, rmat, solprt, statpr, vneed, +!el 4 vsave, x0prt +! +! *** iv subscript values *** +! +!/6 +! data algsav/51/, covprt/14/, covreq/15/, dtype/16/, hc/71/, +! 1 ierr/75/, inith/25/, inits/25/, ipivot/76/, ivneed/3/, +! 2 lastiv/44/, lastv/45/, lmat/42/, mxfcal/17/, mxiter/18/, +! 3 nfcov/52/, ngcov/53/, nvdflt/50/, outlev/19/, parprt/20/, +! 4 parsav/49/, perm/58/, prunit/21/, qrtyp/80/, rdreq/57/, +! 5 rmat/78/, solprt/22/, statpr/23/, vneed/4/, vsave/60/, +! 6 x0prt/24/ +!/7 + integer,parameter :: algsav=51, covprt=14, covreq=15, dtype=16, hc=71,& + ierr=75, inith=25, inits=25, ipivot=76, ivneed=3,& + lastiv=44, lastv=45, lmat=42, mxfcal=17, mxiter=18,& + nfcov=52, ngcov=53, nvdflt=50, outlev=19, parprt=20,& + parsav=49, perm=58, prunit=21, qrtyp=80, rdreq=57,& + rmat=78, solprt=22, statpr=23, vneed=4, vsave=60,& + x0prt=24 +!/ + data miniv(1)/80/, miniv(2)/59/, minv(1)/98/, minv(2)/71/ +!el local variables + integer :: mv +! +!------------------------------- body -------------------------------- +! + if (alg .lt. 1 .or. alg .gt. 2) go to 40 + miv = miniv(alg) + if (liv .lt. miv) go to 20 + mv = minv(alg) + if (lv .lt. mv) go to 30 + call vdflt(alg, lv, v) + iv(1) = 12 + iv(algsav) = alg + iv(ivneed) = 0 + iv(lastiv) = miv + iv(lastv) = mv + iv(lmat) = mv + 1 + iv(mxfcal) = 200 + iv(mxiter) = 150 + iv(outlev) = 1 + iv(parprt) = 1 + iv(perm) = miv + 1 + iv(prunit) = imdcon(1) + iv(solprt) = 1 + iv(statpr) = 1 + iv(vneed) = 0 + iv(x0prt) = 1 +! + if (alg .ge. 2) go to 10 +! +! *** regression values +! + iv(covprt) = 3 + iv(covreq) = 1 + iv(dtype) = 1 + iv(hc) = 0 + iv(ierr) = 0 + iv(inits) = 0 + iv(ipivot) = 0 + iv(nvdflt) = 32 + iv(parsav) = 67 + iv(qrtyp) = 1 + iv(rdreq) = 3 + iv(rmat) = 0 + iv(vsave) = 58 + go to 999 +! +! *** general optimization values +! + 10 iv(dtype) = 0 + iv(inith) = 1 + iv(nfcov) = 0 + iv(ngcov) = 0 + iv(nvdflt) = 25 + iv(parsav) = 47 + go to 999 +! + 20 iv(1) = 15 + go to 999 +! + 30 iv(1) = 16 + go to 999 +! + 40 iv(1) = 67 +! + 999 return +! *** last card of deflt follows *** + end subroutine deflt +!----------------------------------------------------------------------------- + real(kind=8) function dotprd(p,x,y) +! +! *** return the inner product of the p-vectors x and y. *** +! + integer :: p + real(kind=8) :: x(p), y(p) +! + integer :: i +!el real(kind=8) :: one, zero + real(kind=8) :: sqteta, t +!/+ +!el real(kind=8) :: dmax1, dabs +!/ +!el external rmdcon +!el real(kind=8) :: rmdcon +! +! *** rmdcon(2) returns a machine-dependent constant, sqteta, which +! *** is slightly larger than the smallest positive number that +! *** can be squared without underflowing. +! +!/6 +! data one/1.d+0/, sqteta/0.d+0/, zero/0.d+0/ +!/7 + real(kind=8),parameter :: one=1.d+0, zero=0.d+0 + data sqteta/0.d+0/ +!/ +! + dotprd = zero + if (p .le. 0) go to 999 +!rc if (sqteta .eq. zero) sqteta = rmdcon(2) + do 20 i = 1, p +!rc t = dmax1(dabs(x(i)), dabs(y(i))) +!rc if (t .gt. one) go to 10 +!rc if (t .lt. sqteta) go to 20 +!rc t = (x(i)/sqteta)*y(i) +!rc if (dabs(t) .lt. sqteta) go to 20 + 10 dotprd = dotprd + x(i)*y(i) + 20 continue +! + 999 return +! *** last card of dotprd follows *** + end function dotprd +!----------------------------------------------------------------------------- + subroutine itsum(d, g, iv, liv, lv, p, v, x) +! +! *** print iteration summary for ***sol (version 2.3) *** +! +! *** parameter declarations *** +! + integer :: liv, lv, p + integer :: iv(liv) + real(kind=8) :: d(p), g(p), v(lv), x(p) +! +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! +! *** local variables *** +! + integer :: alg, i, iv1, m, nf, ng, ol, pu +!/6 +! real model1(6), model2(6) +!/7 + character(len=4) :: model1(6), model2(6) +!/ + real(kind=8) :: nreldf, oldf, preldf, reldf !el, zero +! +! *** intrinsic functions *** +!/+ +!el integer :: iabs +!el real(kind=8) :: dabs, dmax1 +!/ +! *** no external functions or subroutines *** +! +! *** subscripts for iv and v *** +! +!el integer algsav, dstnrm, f, fdif, f0, needhd, nfcall, nfcov, ngcov, +!el 1 ngcall, niter, nreduc, outlev, preduc, prntit, prunit, +!el 2 reldx, solprt, statpr, stppar, sused, x0prt +! +! *** iv subscript values *** +! +!/6 +! data algsav/51/, needhd/36/, nfcall/6/, nfcov/52/, ngcall/30/, +! 1 ngcov/53/, niter/31/, outlev/19/, prntit/39/, prunit/21/, +! 2 solprt/22/, statpr/23/, sused/64/, x0prt/24/ +!/7 + integer,parameter :: algsav=51, needhd=36, nfcall=6, nfcov=52, ngcall=30,& + ngcov=53, niter=31, outlev=19, prntit=39, prunit=21,& + solprt=22, statpr=23, sused=64, x0prt=24 +!/ +! +! *** v subscript values *** +! +!/6 +! data dstnrm/2/, f/10/, f0/13/, fdif/11/, nreduc/6/, preduc/7/, +! 1 reldx/17/, stppar/5/ +!/7 + integer,parameter :: dstnrm=2, f=10, f0=13, fdif=11, nreduc=6, preduc=7,& + reldx=17, stppar=5 +!/ +! +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +!/6 +! data model1(1)/4h /, model1(2)/4h /, model1(3)/4h /, +! 1 model1(4)/4h /, model1(5)/4h g /, model1(6)/4h s /, +! 2 model2(1)/4h g /, model2(2)/4h s /, model2(3)/4hg-s /, +! 3 model2(4)/4hs-g /, model2(5)/4h-s-g/, model2(6)/4h-g-s/ +!/7 + data model1/' ',' ',' ',' ',' g ',' s '/,& + model2/' g ',' s ','g-s ','s-g ','-s-g','-g-s'/ +!/ +! +!------------------------------- body -------------------------------- +! + pu = iv(prunit) + if (pu .eq. 0) go to 999 + iv1 = iv(1) + if (iv1 .gt. 62) iv1 = iv1 - 51 + ol = iv(outlev) + alg = iv(algsav) + if (iv1 .lt. 2 .or. iv1 .gt. 15) go to 370 + if (iv1 .ge. 12) go to 120 + if (iv1 .eq. 2 .and. iv(niter) .eq. 0) go to 390 + if (ol .eq. 0) go to 120 + if (iv1 .ge. 10 .and. iv(prntit) .eq. 0) go to 120 + if (iv1 .gt. 2) go to 10 + iv(prntit) = iv(prntit) + 1 + if (iv(prntit) .lt. iabs(ol)) go to 999 + 10 nf = iv(nfcall) - iabs(iv(nfcov)) + iv(prntit) = 0 + reldf = zero + preldf = zero + oldf = dmax1(dabs(v(f0)), dabs(v(f))) + if (oldf .le. zero) go to 20 + reldf = v(fdif) / oldf + preldf = v(preduc) / oldf + 20 if (ol .gt. 0) go to 60 +! +! *** print short summary line *** +! + if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,30) + 30 format(/10h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,& + 2x,13hmodel stppar) + if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,40) + 40 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,& + 3x,6hstppar) + iv(needhd) = 0 + if (alg .eq. 2) go to 50 + m = iv(sused) + write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),& + model1(m), model2(m), v(stppar) + go to 120 +! + 50 write(pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx),& + v(stppar) + go to 120 +! +! *** print long summary line *** +! + 60 if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,70) + 70 format(/11h it nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,& + 2x,13hmodel stppar,2x,6hd*step,2x,7hnpreldf) + if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,80) + 80 format(/11h it nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,& + 3x,6hstppar,3x,6hd*step,3x,7hnpreldf) + iv(needhd) = 0 + nreldf = zero + if (oldf .gt. zero) nreldf = v(nreduc) / oldf + if (alg .eq. 2) go to 90 + m = iv(sused) + write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),& + model1(m), model2(m), v(stppar), v(dstnrm), nreldf + go to 120 +! + 90 write(pu,110) iv(niter), nf, v(f), reldf, preldf,& + v(reldx), v(stppar), v(dstnrm), nreldf + 100 format(i6,i5,d10.3,2d9.2,d8.1,a3,a4,2d8.1,d9.2) + 110 format(i6,i5,d11.3,2d10.2,3d9.1,d10.2) +! + 120 if (iv(statpr) .lt. 0) go to 430 + go to (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,& + 330, 350, 520), iv1 +! + 130 write(pu,140) + 140 format(/26h ***** x-convergence *****) + go to 430 +! + 150 write(pu,160) + 160 format(/42h ***** relative function convergence *****) + go to 430 +! + 170 write(pu,180) + 180 format(/49h ***** x- and relative function convergence *****) + go to 430 +! + 190 write(pu,200) + 200 format(/42h ***** absolute function convergence *****) + go to 430 +! + 210 write(pu,220) + 220 format(/33h ***** singular convergence *****) + go to 430 +! + 230 write(pu,240) + 240 format(/30h ***** false convergence *****) + go to 430 +! + 250 write(pu,260) + 260 format(/38h ***** function evaluation limit *****) + go to 430 +! + 270 write(pu,280) + 280 format(/28h ***** iteration limit *****) + go to 430 +! + 290 write(pu,300) + 300 format(/18h ***** stopx *****) + go to 430 +! + 310 write(pu,320) + 320 format(/44h ***** initial f(x) cannot be computed *****) +! + go to 390 +! + 330 write(pu,340) + 340 format(/37h ***** bad parameters to assess *****) + go to 999 +! + 350 write(pu,360) + 360 format(/43h ***** gradient could not be computed *****) + if (iv(niter) .gt. 0) go to 480 + go to 390 +! + 370 write(pu,380) iv(1) + 380 format(/14h ***** iv(1) =,i5,6h *****) + go to 999 +! +! *** initial call on itsum *** +! + 390 if (iv(x0prt) .ne. 0) write(pu,400) (i, x(i), d(i), i = 1, p) + 400 format(/23h i initial x(i),8x,4hd(i)//(1x,i5,d17.6,d14.3)) +! *** the following are to avoid undefined variables when the +! *** function evaluation limit is 1... + v(dstnrm) = zero + v(fdif) = zero + v(nreduc) = zero + v(preduc) = zero + v(reldx) = zero + if (iv1 .ge. 12) go to 999 + iv(needhd) = 0 + iv(prntit) = 0 + if (ol .eq. 0) go to 999 + if (ol .lt. 0 .and. alg .eq. 1) write(pu,30) + if (ol .lt. 0 .and. alg .eq. 2) write(pu,40) + if (ol .gt. 0 .and. alg .eq. 1) write(pu,70) + if (ol .gt. 0 .and. alg .eq. 2) write(pu,80) + if (alg .eq. 1) write(pu,410) v(f) + if (alg .eq. 2) write(pu,420) v(f) + 410 format(/11h 0 1,d10.3) +!365 format(/11h 0 1,e11.3) + 420 format(/11h 0 1,d11.3) + go to 999 +! +! *** print various information requested on solution *** +! + 430 iv(needhd) = 1 + if (iv(statpr) .eq. 0) go to 480 + oldf = dmax1(dabs(v(f0)), dabs(v(f))) + preldf = zero + nreldf = zero + if (oldf .le. zero) go to 440 + preldf = v(preduc) / oldf + nreldf = v(nreduc) / oldf + 440 nf = iv(nfcall) - iv(nfcov) + ng = iv(ngcall) - iv(ngcov) + write(pu,450) v(f), v(reldx), nf, ng, preldf, nreldf + 450 format(/9h function,d17.6,8h reldx,d17.3/12h func. evals,& + i8,9x,11hgrad. evals,i8/7h preldf,d16.3,6x,7hnpreldf,d15.3) +! + if (iv(nfcov) .gt. 0) write(pu,460) iv(nfcov) + 460 format(/1x,i4,50h extra func. evals for covariance and diagnostics.) + if (iv(ngcov) .gt. 0) write(pu,470) iv(ngcov) + 470 format(1x,i4,50h extra grad. evals for covariance and diagnostics.) +! + 480 if (iv(solprt) .eq. 0) go to 999 + iv(needhd) = 1 + write(pu,490) + 490 format(/22h i final x(i),8x,4hd(i),10x,4hg(i)/) + do 500 i = 1, p + write(pu,510) i, x(i), d(i), g(i) + 500 continue + 510 format(1x,i5,d16.6,2d14.3) + go to 999 +! + 520 write(pu,530) + 530 format(/24h inconsistent dimensions) + 999 return +! *** last card of itsum follows *** + end subroutine itsum +!----------------------------------------------------------------------------- + subroutine litvmu(n, x, l, y) +! +! *** solve (l**t)*x = y, where l is an n x n lower triangular +! *** matrix stored compactly by rows. x and y may occupy the same +! *** storage. *** +! + integer :: n +!al real(kind=8) :: x(n), l(1), y(n) + real(kind=8) :: x(n), l(n*(n+1)/2), y(n) + integer :: i, ii, ij, im1, i0, j, np1 + real(kind=8) :: xi !el, zero +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +! + do 10 i = 1, n + 10 x(i) = y(i) + np1 = n + 1 + i0 = n*(n+1)/2 + do 30 ii = 1, n + i = np1 - ii + xi = x(i)/l(i0) + x(i) = xi + if (i .le. 1) go to 999 + i0 = i0 - i + if (xi .eq. zero) go to 30 + im1 = i - 1 + do 20 j = 1, im1 + ij = i0 + j + x(j) = x(j) - xi*l(ij) + 20 continue + 30 continue + 999 return +! *** last card of litvmu follows *** + end subroutine litvmu +!----------------------------------------------------------------------------- + subroutine livmul(n, x, l, y) +! +! *** solve l*x = y, where l is an n x n lower triangular +! *** matrix stored compactly by rows. x and y may occupy the same +! *** storage. *** +! + integer :: n +!al real(kind=8) :: x(n), l(1), y(n) + real(kind=8) :: x(n), l(n*(n+1)/2), y(n) +!el external dotprd +!el real(kind=8) :: dotprd + integer :: i, j, k + real(kind=8) :: t !el, zero +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +! + do 10 k = 1, n + if (y(k) .ne. zero) go to 20 + x(k) = zero + 10 continue + go to 999 + 20 j = k*(k+1)/2 + x(k) = y(k) / l(j) + if (k .ge. n) go to 999 + k = k + 1 + do 30 i = k, n + t = dotprd(i-1, l(j+1), x) + j = j + i + x(i) = (y(i) - t)/l(j) + 30 continue + 999 return +! *** last card of livmul follows *** + end subroutine livmul +!----------------------------------------------------------------------------- + subroutine parck(alg, d, iv, liv, lv, n, v) +! +! *** check ***sol (version 2.3) parameters, print changed values *** +! +! *** alg = 1 for regression, alg = 2 for general unconstrained opt. +! + integer :: alg, liv, lv, n + integer :: iv(liv) + real(kind=8) :: d(n), v(lv) +! +!el external rmdcon, vcopy, vdflt +!el real(kind=8) :: rmdcon +! rmdcon -- returns machine-dependent constants. +! vcopy -- copies one vector to another. +! vdflt -- supplies default parameter values to v alone. +!/+ +!el integer :: max0 +!/ +! +! *** local variables *** +! + integer :: i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu + integer :: ijmp, jlim(2), miniv(2), ndflt(2) +!/6 +! integer varnm(2), sh(2) +! real cngd(3), dflt(3), vn(2,34), which(3) +!/7 + character(len=1) :: varnm(2), sh(2) + character(len=4) :: cngd(3), dflt(3), vn(2,34), which(3) +!/ + real(kind=8) :: big, machep, tiny, vk, vm(34), vx(34), zero +! +! *** iv and v subscripts *** +! +!el integer algsav, dinit, dtype, dtype0, epslon, inits, ivneed, +!el 1 lastiv, lastv, lmat, nextiv, nextv, nvdflt, oldn, +!el 2 parprt, parsav, perm, prunit, vneed +! +! +!/6 +! data algsav/51/, dinit/38/, dtype/16/, dtype0/54/, epslon/19/, +! 1 inits/25/, ivneed/3/, lastiv/44/, lastv/45/, lmat/42/, +! 2 nextiv/46/, nextv/47/, nvdflt/50/, oldn/38/, parprt/20/, +! 3 parsav/49/, perm/58/, prunit/21/, vneed/4/ +!/7 + integer,parameter :: algsav=51, dinit=38, dtype=16, dtype0=54, epslon=19,& + inits=25, ivneed=3, lastiv=44, lastv=45, lmat=42,& + nextiv=46, nextv=47, nvdflt=50, oldn=38, parprt=20,& + parsav=49, perm=58, prunit=21, vneed=4 + save big, machep, tiny +!/ +! + data big/0.d+0/, machep/-1.d+0/, tiny/1.d+0/, zero/0.d+0/ +!/6 +! data vn(1,1),vn(2,1)/4hepsl,4hon../ +! data vn(1,2),vn(2,2)/4hphmn,4hfc../ +! data vn(1,3),vn(2,3)/4hphmx,4hfc../ +! data vn(1,4),vn(2,4)/4hdecf,4hac../ +! data vn(1,5),vn(2,5)/4hincf,4hac../ +! data vn(1,6),vn(2,6)/4hrdfc,4hmn../ +! data vn(1,7),vn(2,7)/4hrdfc,4hmx../ +! data vn(1,8),vn(2,8)/4htune,4hr1../ +! data vn(1,9),vn(2,9)/4htune,4hr2../ +! data vn(1,10),vn(2,10)/4htune,4hr3../ +! data vn(1,11),vn(2,11)/4htune,4hr4../ +! data vn(1,12),vn(2,12)/4htune,4hr5../ +! data vn(1,13),vn(2,13)/4hafct,4hol../ +! data vn(1,14),vn(2,14)/4hrfct,4hol../ +! data vn(1,15),vn(2,15)/4hxcto,4hl.../ +! data vn(1,16),vn(2,16)/4hxfto,4hl.../ +! data vn(1,17),vn(2,17)/4hlmax,4h0.../ +! data vn(1,18),vn(2,18)/4hlmax,4hs.../ +! data vn(1,19),vn(2,19)/4hscto,4hl.../ +! data vn(1,20),vn(2,20)/4hdini,4ht.../ +! data vn(1,21),vn(2,21)/4hdtin,4hit../ +! data vn(1,22),vn(2,22)/4hd0in,4hit../ +! data vn(1,23),vn(2,23)/4hdfac,4h..../ +! data vn(1,24),vn(2,24)/4hdltf,4hdc../ +! data vn(1,25),vn(2,25)/4hdltf,4hdj../ +! data vn(1,26),vn(2,26)/4hdelt,4ha0../ +! data vn(1,27),vn(2,27)/4hfuzz,4h..../ +! data vn(1,28),vn(2,28)/4hrlim,4hit../ +! data vn(1,29),vn(2,29)/4hcosm,4hin../ +! data vn(1,30),vn(2,30)/4hhube,4hrc../ +! data vn(1,31),vn(2,31)/4hrspt,4hol../ +! data vn(1,32),vn(2,32)/4hsigm,4hin../ +! data vn(1,33),vn(2,33)/4heta0,4h..../ +! data vn(1,34),vn(2,34)/4hbias,4h..../ +!/7 + data vn(1,1),vn(2,1)/'epsl','on..'/ + data vn(1,2),vn(2,2)/'phmn','fc..'/ + data vn(1,3),vn(2,3)/'phmx','fc..'/ + data vn(1,4),vn(2,4)/'decf','ac..'/ + data vn(1,5),vn(2,5)/'incf','ac..'/ + data vn(1,6),vn(2,6)/'rdfc','mn..'/ + data vn(1,7),vn(2,7)/'rdfc','mx..'/ + data vn(1,8),vn(2,8)/'tune','r1..'/ + data vn(1,9),vn(2,9)/'tune','r2..'/ + data vn(1,10),vn(2,10)/'tune','r3..'/ + data vn(1,11),vn(2,11)/'tune','r4..'/ + data vn(1,12),vn(2,12)/'tune','r5..'/ + data vn(1,13),vn(2,13)/'afct','ol..'/ + data vn(1,14),vn(2,14)/'rfct','ol..'/ + data vn(1,15),vn(2,15)/'xcto','l...'/ + data vn(1,16),vn(2,16)/'xfto','l...'/ + data vn(1,17),vn(2,17)/'lmax','0...'/ + data vn(1,18),vn(2,18)/'lmax','s...'/ + data vn(1,19),vn(2,19)/'scto','l...'/ + data vn(1,20),vn(2,20)/'dini','t...'/ + data vn(1,21),vn(2,21)/'dtin','it..'/ + data vn(1,22),vn(2,22)/'d0in','it..'/ + data vn(1,23),vn(2,23)/'dfac','....'/ + data vn(1,24),vn(2,24)/'dltf','dc..'/ + data vn(1,25),vn(2,25)/'dltf','dj..'/ + data vn(1,26),vn(2,26)/'delt','a0..'/ + data vn(1,27),vn(2,27)/'fuzz','....'/ + data vn(1,28),vn(2,28)/'rlim','it..'/ + data vn(1,29),vn(2,29)/'cosm','in..'/ + data vn(1,30),vn(2,30)/'hube','rc..'/ + data vn(1,31),vn(2,31)/'rspt','ol..'/ + data vn(1,32),vn(2,32)/'sigm','in..'/ + data vn(1,33),vn(2,33)/'eta0','....'/ + data vn(1,34),vn(2,34)/'bias','....'/ +!/ +! + data vm(1)/1.0d-3/, vm(2)/-0.99d+0/, vm(3)/1.0d-3/, vm(4)/1.0d-2/,& + vm(5)/1.2d+0/, vm(6)/1.d-2/, vm(7)/1.2d+0/, vm(8)/0.d+0/,& + vm(9)/0.d+0/, vm(10)/1.d-3/, vm(11)/-1.d+0/, vm(13)/0.d+0/,& + vm(15)/0.d+0/, vm(16)/0.d+0/, vm(19)/0.d+0/, vm(20)/-10.d+0/,& + vm(21)/0.d+0/, vm(22)/0.d+0/, vm(23)/0.d+0/, vm(27)/1.01d+0/,& + vm(28)/1.d+10/, vm(30)/0.d+0/, vm(31)/0.d+0/, vm(32)/0.d+0/,& + vm(34)/0.d+0/ + data vx(1)/0.9d+0/, vx(2)/-1.d-3/, vx(3)/1.d+1/, vx(4)/0.8d+0/,& + vx(5)/1.d+2/, vx(6)/0.8d+0/, vx(7)/1.d+2/, vx(8)/0.5d+0/,& + vx(9)/0.5d+0/, vx(10)/1.d+0/, vx(11)/1.d+0/, vx(14)/0.1d+0/,& + vx(15)/1.d+0/, vx(16)/1.d+0/, vx(19)/1.d+0/, vx(23)/1.d+0/,& + vx(24)/1.d+0/, vx(25)/1.d+0/, vx(26)/1.d+0/, vx(27)/1.d+10/,& + vx(29)/1.d+0/, vx(31)/1.d+0/, vx(32)/1.d+0/, vx(33)/1.d+0/,& + vx(34)/1.d+0/ +! +!/6 +! data varnm(1)/1hp/, varnm(2)/1hn/, sh(1)/1hs/, sh(2)/1hh/ +! data cngd(1),cngd(2),cngd(3)/4h---c,4hhang,4hed v/, +! 1 dflt(1),dflt(2),dflt(3)/4hnond,4hefau,4hlt v/ +!/7 + data varnm(1)/'p'/, varnm(2)/'n'/, sh(1)/'s'/, sh(2)/'h'/ + data cngd(1),cngd(2),cngd(3)/'---c','hang','ed v'/,& + dflt(1),dflt(2),dflt(3)/'nond','efau','lt v'/ +!/ + data ijmp/33/, jlim(1)/0/, jlim(2)/24/, ndflt(1)/32/, ndflt(2)/25/ + data miniv(1)/80/, miniv(2)/59/ +! +!............................... body ................................ +! + pu = 0 + if (prunit .le. liv) pu = iv(prunit) + if (alg .lt. 1 .or. alg .gt. 2) go to 340 + if (iv(1) .eq. 0) call deflt(alg, iv, liv, lv, v) + iv1 = iv(1) + if (iv1 .ne. 13 .and. iv1 .ne. 12) go to 10 + miv1 = miniv(alg) + if (perm .le. liv) miv1 = max0(miv1, iv(perm) - 1) + if (ivneed .le. liv) miv2 = miv1 + max0(iv(ivneed), 0) + if (lastiv .le. liv) iv(lastiv) = miv2 + if (liv .lt. miv1) go to 300 + iv(ivneed) = 0 + iv(lastv) = max0(iv(vneed), 0) + iv(lmat) - 1 + iv(vneed) = 0 + if (liv .lt. miv2) go to 300 + if (lv .lt. iv(lastv)) go to 320 + 10 if (alg .eq. iv(algsav)) go to 30 + if (pu .ne. 0) write(pu,20) alg, iv(algsav) + 20 format(/39h the first parameter to deflt should be,i3,& + 12h rather than,i3) + iv(1) = 82 + go to 999 + 30 if (iv1 .lt. 12 .or. iv1 .gt. 14) go to 60 + if (n .ge. 1) go to 50 + iv(1) = 81 + if (pu .eq. 0) go to 999 + write(pu,40) varnm(alg), n + 40 format(/8h /// bad,a1,2h =,i5) + go to 999 + 50 if (iv1 .ne. 14) iv(nextiv) = iv(perm) + if (iv1 .ne. 14) iv(nextv) = iv(lmat) + if (iv1 .eq. 13) go to 999 + k = iv(parsav) - epslon + call vdflt(alg, lv-k, v(k+1)) + iv(dtype0) = 2 - alg + iv(oldn) = n + which(1) = dflt(1) + which(2) = dflt(2) + which(3) = dflt(3) + go to 110 + 60 if (n .eq. iv(oldn)) go to 80 + iv(1) = 17 + if (pu .eq. 0) go to 999 + write(pu,70) varnm(alg), iv(oldn), n + 70 format(/5h /// ,1a1,14h changed from ,i5,4h to ,i5) + go to 999 +! + 80 if (iv1 .le. 11 .and. iv1 .ge. 1) go to 100 + iv(1) = 80 + if (pu .ne. 0) write(pu,90) iv1 + 90 format(/13h /// iv(1) =,i5,28h should be between 0 and 14.) + go to 999 +! + 100 which(1) = cngd(1) + which(2) = cngd(2) + which(3) = cngd(3) +! + 110 if (iv1 .eq. 14) iv1 = 12 + if (big .gt. tiny) go to 120 + tiny = rmdcon(1) + machep = rmdcon(3) + big = rmdcon(6) + vm(12) = machep + vx(12) = big + vx(13) = big + vm(14) = machep + vm(17) = tiny + vx(17) = big + vm(18) = tiny + vx(18) = big + vx(20) = big + vx(21) = big + vx(22) = big + vm(24) = machep + vm(25) = machep + vm(26) = machep + vx(28) = rmdcon(5) + vm(29) = machep + vx(30) = big + vm(33) = machep + 120 m = 0 + i = 1 + j = jlim(alg) + k = epslon + ndfalt = ndflt(alg) + do 150 l = 1, ndfalt + vk = v(k) + if (vk .ge. vm(i) .and. vk .le. vx(i)) go to 140 + m = k + if (pu .ne. 0) write(pu,130) vn(1,i), vn(2,i), k, vk,& + vm(i), vx(i) + 130 format(/6h /// ,2a4,5h.. v(,i2,3h) =,d11.3,7h should,& + 11h be between,d11.3,4h and,d11.3) + 140 k = k + 1 + i = i + 1 + if (i .eq. j) i = ijmp + 150 continue +! + if (iv(nvdflt) .eq. ndfalt) go to 170 + iv(1) = 51 + if (pu .eq. 0) go to 999 + write(pu,160) iv(nvdflt), ndfalt + 160 format(/13h iv(nvdflt) =,i5,13h rather than ,i5) + go to 999 + 170 if ((iv(dtype) .gt. 0 .or. v(dinit) .gt. zero) .and. iv1 .eq. 12) & + go to 200 + do 190 i = 1, n + if (d(i) .gt. zero) go to 190 + m = 18 + if (pu .ne. 0) write(pu,180) i, d(i) + 180 format(/8h /// d(,i3,3h) =,d11.3,19h should be positive) + 190 continue + 200 if (m .eq. 0) go to 210 + iv(1) = m + go to 999 +! + 210 if (pu .eq. 0 .or. iv(parprt) .eq. 0) go to 999 + if (iv1 .ne. 12 .or. iv(inits) .eq. alg-1) go to 230 + m = 1 + write(pu,220) sh(alg), iv(inits) + 220 format(/22h nondefault values..../5h init,a1,14h..... iv(25) =,& + i3) + 230 if (iv(dtype) .eq. iv(dtype0)) go to 250 + if (m .eq. 0) write(pu,260) which + m = 1 + write(pu,240) iv(dtype) + 240 format(20h dtype..... iv(16) =,i3) + 250 i = 1 + j = jlim(alg) + k = epslon + l = iv(parsav) + ndfalt = ndflt(alg) + do 290 ii = 1, ndfalt + if (v(k) .eq. v(l)) go to 280 + if (m .eq. 0) write(pu,260) which + 260 format(/1h ,3a4,9halues..../) + m = 1 + write(pu,270) vn(1,i), vn(2,i), k, v(k) + 270 format(1x,2a4,5h.. v(,i2,3h) =,d15.7) + 280 k = k + 1 + l = l + 1 + i = i + 1 + if (i .eq. j) i = ijmp + 290 continue +! + iv(dtype0) = iv(dtype) + parsv1 = iv(parsav) + call vcopy(iv(nvdflt), v(parsv1), v(epslon)) + go to 999 +! + 300 iv(1) = 15 + if (pu .eq. 0) go to 999 + write(pu,310) liv, miv2 + 310 format(/10h /// liv =,i5,17h must be at least,i5) + if (liv .lt. miv1) go to 999 + if (lv .lt. iv(lastv)) go to 320 + go to 999 +! + 320 iv(1) = 16 + if (pu .eq. 0) go to 999 + write(pu,330) lv, iv(lastv) + 330 format(/9h /// lv =,i5,17h must be at least,i5) + go to 999 +! + 340 iv(1) = 67 + if (pu .eq. 0) go to 999 + write(pu,350) alg + 350 format(/10h /// alg =,i5,15h must be 1 or 2) +! + 999 return +! *** last card of parck follows *** + end subroutine parck +!----------------------------------------------------------------------------- + real(kind=8) function reldst(p, d, x, x0) +! +! *** compute and return relative difference between x and x0 *** +! *** nl2sol version 2.2 *** +! + integer :: p + real(kind=8) :: d(p), x(p), x0(p) +!/+ +!el real(kind=8) :: dabs +!/ + integer :: i + real(kind=8) :: emax, t, xmax !el, zero +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +! + emax = zero + xmax = zero + do 10 i = 1, p + t = dabs(d(i) * (x(i) - x0(i))) + if (emax .lt. t) emax = t + t = d(i) * (dabs(x(i)) + dabs(x0(i))) + if (xmax .lt. t) xmax = t + 10 continue + reldst = zero + if (xmax .gt. zero) reldst = emax / xmax + 999 return +! *** last card of reldst follows *** + end function reldst +!----------------------------------------------------------------------------- + subroutine vaxpy(p, w, a, x, y) +! +! *** set w = a*x + y -- w, x, y = p-vectors, a = scalar *** +! + integer :: p + real(kind=8) :: a, w(p), x(p), y(p) +! + integer :: i +! + do 10 i = 1, p + 10 w(i) = a*x(i) + y(i) + return + end subroutine vaxpy +!----------------------------------------------------------------------------- + subroutine vcopy(p, y, x) +! +! *** set y = x, where x and y are p-vectors *** +! + integer :: p + real(kind=8) :: x(p), y(p) +! + integer :: i +! + do 10 i = 1, p + 10 y(i) = x(i) + return + end subroutine vcopy +!----------------------------------------------------------------------------- + subroutine vdflt(alg, lv, v) +! +! *** supply ***sol (version 2.3) default values to v *** +! +! *** alg = 1 means regression constants. +! *** alg = 2 means general unconstrained optimization constants. +! + integer :: alg, l,lv + real(kind=8) :: v(lv) +!/+ +!el real(kind=8) :: dmax1 +!/ +!el external rmdcon +!el real(kind=8) :: rmdcon +! rmdcon... returns machine-dependent constants +! + real(kind=8) :: machep, mepcrt, sqteps !el one, three +! +! *** subscripts for v *** +! +!el integer afctol, bias, cosmin, decfac, delta0, dfac, dinit, dltfdc, +!el 1 dltfdj, dtinit, d0init, epslon, eta0, fuzz, huberc, +!el 2 incfac, lmax0, lmaxs, phmnfc, phmxfc, rdfcmn, rdfcmx, +!el 3 rfctol, rlimit, rsptol, sctol, sigmin, tuner1, tuner2, +!el 4 tuner3, tuner4, tuner5, xctol, xftol +! +!/6 +! data one/1.d+0/, three/3.d+0/ +!/7 + real(kind=8),parameter :: one=1.d+0, three=3.d+0 +!/ +! +! *** v subscript values *** +! +!/6 +! data afctol/31/, bias/43/, cosmin/47/, decfac/22/, delta0/44/, +! 1 dfac/41/, dinit/38/, dltfdc/42/, dltfdj/43/, dtinit/39/, +! 2 d0init/40/, epslon/19/, eta0/42/, fuzz/45/, huberc/48/, +! 3 incfac/23/, lmax0/35/, lmaxs/36/, phmnfc/20/, phmxfc/21/, +! 4 rdfcmn/24/, rdfcmx/25/, rfctol/32/, rlimit/46/, rsptol/49/, +! 5 sctol/37/, sigmin/50/, tuner1/26/, tuner2/27/, tuner3/28/, +! 6 tuner4/29/, tuner5/30/, xctol/33/, xftol/34/ +!/7 + integer,parameter :: afctol=31, bias=43, cosmin=47, decfac=22, delta0=44,& + dfac=41, dinit=38, dltfdc=42, dltfdj=43, dtinit=39,& + d0init=40, epslon=19, eta0=42, fuzz=45, huberc=48,& + incfac=23, lmax0=35, lmaxs=36, phmnfc=20, phmxfc=21,& + rdfcmn=24, rdfcmx=25, rfctol=32, rlimit=46, rsptol=49,& + sctol=37, sigmin=50, tuner1=26, tuner2=27, tuner3=28,& + tuner4=29, tuner5=30, xctol=33, xftol=34 +!/ +! +!------------------------------- body -------------------------------- +! + machep = rmdcon(3) + v(afctol) = 1.d-20 + if (machep .gt. 1.d-10) v(afctol) = machep**2 + v(decfac) = 0.5d+0 + sqteps = rmdcon(4) + v(dfac) = 0.6d+0 + v(delta0) = sqteps + v(dtinit) = 1.d-6 + mepcrt = machep ** (one/three) + v(d0init) = 1.d+0 + v(epslon) = 0.1d+0 + v(incfac) = 2.d+0 + v(lmax0) = 1.d+0 + v(lmaxs) = 1.d+0 + v(phmnfc) = -0.1d+0 + v(phmxfc) = 0.1d+0 + v(rdfcmn) = 0.1d+0 + v(rdfcmx) = 4.d+0 + v(rfctol) = dmax1(1.d-10, mepcrt**2) + v(sctol) = v(rfctol) + v(tuner1) = 0.1d+0 + v(tuner2) = 1.d-4 + v(tuner3) = 0.75d+0 + v(tuner4) = 0.5d+0 + v(tuner5) = 0.75d+0 + v(xctol) = sqteps + v(xftol) = 1.d+2 * machep +! + if (alg .ge. 2) go to 10 +! +! *** regression values +! + v(cosmin) = dmax1(1.d-6, 1.d+2 * machep) + v(dinit) = 0.d+0 + v(dltfdc) = mepcrt + v(dltfdj) = sqteps + v(fuzz) = 1.5d+0 + v(huberc) = 0.7d+0 + v(rlimit) = rmdcon(5) + v(rsptol) = 1.d-3 + v(sigmin) = 1.d-4 + go to 999 +! +! *** general optimization values +! + 10 v(bias) = 0.8d+0 + v(dinit) = -1.0d+0 + v(eta0) = 1.0d+3 * machep +! + 999 return +! *** last card of vdflt follows *** + end subroutine vdflt +!----------------------------------------------------------------------------- + subroutine vscopy(p, y, s) +! +! *** set p-vector y to scalar s *** +! + integer :: p + real(kind=8) :: s, y(p) +! + integer :: i +! + do 10 i = 1, p + 10 y(i) = s + return + end subroutine vscopy +!----------------------------------------------------------------------------- + real(kind=8) function v2norm(p, x) +! +! *** return the 2-norm of the p-vector x, taking *** +! *** care to avoid the most likely underflows. *** +! + integer :: p + real(kind=8) :: x(p) +! + integer :: i, j + real(kind=8) :: r, scale, sqteta, t, xi !el, one, zero +!/+ +!el real(kind=8) :: dabs, dsqrt +!/ +!el external rmdcon +!el real(kind=8) :: rmdcon +! +!/6 +! data one/1.d+0/, zero/0.d+0/ +!/7 + real(kind=8),parameter :: one=1.d+0, zero=0.d+0 + save sqteta +!/ + data sqteta/0.d+0/ +! + if (p .gt. 0) go to 10 + v2norm = zero + go to 999 + 10 do 20 i = 1, p + if (x(i) .ne. zero) go to 30 + 20 continue + v2norm = zero + go to 999 +! + 30 scale = dabs(x(i)) + if (i .lt. p) go to 40 + v2norm = scale + go to 999 + 40 t = one + if (sqteta .eq. zero) sqteta = rmdcon(2) +! +! *** sqteta is (slightly larger than) the square root of the +! *** smallest positive floating point number on the machine. +! *** the tests involving sqteta are done to prevent underflows. +! + j = i + 1 + do 60 i = j, p + xi = dabs(x(i)) + if (xi .gt. scale) go to 50 + r = xi / scale + if (r .gt. sqteta) t = t + r*r + go to 60 + 50 r = scale / xi + if (r .le. sqteta) r = zero + t = one + t * r*r + scale = xi + 60 continue +! + v2norm = scale * dsqrt(t) + 999 return +! *** last card of v2norm follows *** + end function v2norm +!----------------------------------------------------------------------------- + subroutine humsl(n,d,x,calcf,calcgh,iv,liv,lv,v,uiparm,urparm,ufparm) +! +! *** minimize general unconstrained objective function using *** +! *** (analytic) gradient and hessian provided by the caller. *** +! + integer :: liv, lv, n + integer :: iv(liv), uiparm(1) + real(kind=8) :: d(n), x(n), v(lv), urparm(1) + real(kind=8),external :: ufparm +! dimension v(78 + n*(n+12)), uiparm(*), urparm(*) + external :: calcf, calcgh +! +!------------------------------ discussion --------------------------- +! +! this routine is like sumsl, except that the subroutine para- +! meter calcg of sumsl (which computes the gradient of the objec- +! tive function) is replaced by the subroutine parameter calcgh, +! which computes both the gradient and (lower triangle of the) +! hessian of the objective function. the calling sequence is... +! call calcgh(n, x, nf, g, h, uiparm, urparm, ufparm) +! parameters n, x, nf, g, uiparm, urparm, and ufparm are the same +! as for sumsl, while h is an array of length n*(n+1)/2 in which +! calcgh must store the lower triangle of the hessian at x. start- +! ing at h(1), calcgh must store the hessian entries in the order +! (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ... +! the value printed (by itsum) in the column labelled stppar +! is the levenberg-marquardt used in computing the current step. +! zero means a full newton step. if the special case described in +! ref. 1 is detected, then stppar is negated. the value printed +! in the column labelled npreldf is zero if the current hessian +! is not positive definite. +! it sometimes proves worthwhile to let d be determined from the +! diagonal of the hessian matrix by setting iv(dtype) = 1 and +! v(dinit) = 0. the following iv and v components are relevant... +! +! iv(dtol)..... iv(59) gives the starting subscript in v of the dtol +! array used when d is updated. (iv(dtol) can be +! initialized by calling humsl with iv(1) = 13.) +! iv(dtype).... iv(16) tells how the scale vector d should be chosen. +! iv(dtype) .le. 0 means that d should not be updated, and +! iv(dtype) .ge. 1 means that d should be updated as +! described below with v(dfac). default = 0. +! v(dfac)..... v(41) and the dtol and d0 arrays (see v(dtinit) and +! v(d0init)) are used in updating the scale vector d when +! iv(dtype) .gt. 0. (d is initialized according to +! v(dinit), described in sumsl.) let +! d1(i) = max(sqrt(abs(h(i,i))), v(dfac)*d(i)), +! where h(i,i) is the i-th diagonal element of the current +! hessian. if iv(dtype) = 1, then d(i) is set to d1(i) +! unless d1(i) .lt. dtol(i), in which case d(i) is set to +! max(d0(i), dtol(i)). +! if iv(dtype) .ge. 2, then d is updated during the first +! iteration as for iv(dtype) = 1 (after any initialization +! due to v(dinit)) and is left unchanged thereafter. +! default = 0.6. +! v(dtinit)... v(39), if positive, is the value to which all components +! of the dtol array (see v(dfac)) are initialized. if +! v(dtinit) = 0, then it is assumed that the caller has +! stored dtol in v starting at v(iv(dtol)). +! default = 10**-6. +! v(d0init)... v(40), if positive, is the value to which all components +! of the d0 vector (see v(dfac)) are initialized. if +! v(dfac) = 0, then it is assumed that the caller has +! stored d0 in v starting at v(iv(dtol)+n). default = 1.0. +! +! *** reference *** +! +! 1. gay, d.m. (1981), computing optimal locally constrained steps, +! siam j. sci. statist. comput. 2, pp. 186-197. +!. +! *** general *** +! +! coded by david m. gay (winter 1980). revised sept. 1982. +! this subroutine was written in connection with research supported +! in part by the national science foundation under grants +! mcs-7600324 and mcs-7906671. +! +!---------------------------- declarations --------------------------- +! +!el external deflt, humit +! +! deflt... provides default input values for iv and v. +! humit... reverse-communication routine that does humsl algorithm. +! + integer :: g1, h1, iv1, lh, nf + real(kind=8) :: f +! +! *** subscripts for iv *** +! +!el integer g, h, nextv, nfcall, nfgcal, toobig, vneed +! +!/6 +! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, h/56/, toobig/2/, +! 1 vneed/4/ +!/7 + integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28, h=56,& + toobig=2,vneed=4 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! + lh = n * (n + 1) / 2 + if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v) + if (iv(1) .eq. 12 .or. iv(1) .eq. 13) & + iv(vneed) = iv(vneed) + n*(n+3)/2 + iv1 = iv(1) + if (iv1 .eq. 14) go to 10 + if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10 + g1 = 1 + h1 = 1 + if (iv1 .eq. 12) iv(1) = 13 + go to 20 +! + 10 g1 = iv(g) + h1 = iv(h) +! + 20 call humit(d, f, v(g1), v(h1), iv, lh, liv, lv, n, v, x) + if (iv(1) - 2) 30, 40, 50 +! + 30 nf = iv(nfcall) + call calcf(n, x, nf, f, uiparm, urparm, ufparm) + if (nf .le. 0) iv(toobig) = 1 + go to 20 +! + 40 call calcgh(n, x, iv(nfgcal), v(g1), v(h1), uiparm, urparm,& + ufparm) + go to 20 +! + 50 if (iv(1) .ne. 14) go to 999 +! +! *** storage allocation +! + iv(g) = iv(nextv) + iv(h) = iv(g) + n + iv(nextv) = iv(h) + n*(n+1)/2 + if (iv1 .ne. 13) go to 10 +! + 999 return +! *** last card of humsl follows *** + end subroutine humsl +!----------------------------------------------------------------------------- + subroutine humit(d, fx, g, h, iv, lh, liv, lv, n, v, x) +! +! *** carry out humsl (unconstrained minimization) iterations, using +! *** hessian matrix provided by the caller. +! +!el use control + use control, only:stopx + +! *** parameter declarations *** +! + integer :: lh, liv, lv, n + integer :: iv(liv) + real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n) +! +!-------------------------- parameter usage -------------------------- +! +! d.... scale vector. +! fx... function value. +! g.... gradient vector. +! h.... lower triangle of the hessian, stored rowwise. +! iv... integer value array. +! lh... length of h = p*(p+1)/2. +! liv.. length of iv (at least 60). +! lv... length of v (at least 78 + n*(n+21)/2). +! n.... number of variables (components in x and g). +! v.... floating-point value array. +! x.... parameter vector. +! +! *** discussion *** +! +! parameters iv, n, v, and x are the same as the corresponding +! ones to humsl (which see), except that v can be shorter (since +! the part of v that humsl uses for storing g and h is not needed). +! moreover, compared with humsl, iv(1) may have the two additional +! output values 1 and 2, which are explained below, as is the use +! of iv(toobig) and iv(nfgcal). the value iv(g), which is an +! output value from humsl, is not referenced by humit or the +! subroutines it calls. +! +! iv(1) = 1 means the caller should set fx to f(x), the function value +! at x, and call humit again, having changed none of the +! other parameters. an exception occurs if f(x) cannot be +! computed (e.g. if overflow would occur), which may happen +! because of an oversized step. in this case the caller +! should set iv(toobig) = iv(2) to 1, which will cause +! humit to ignore fx and try a smaller step. the para- +! meter nf that humsl passes to calcf (for possible use by +! calcgh) is a copy of iv(nfcall) = iv(6). +! iv(1) = 2 means the caller should set g to g(x), the gradient of f at +! x, and h to the lower triangle of h(x), the hessian of f +! at x, and call humit again, having changed none of the +! other parameters except perhaps the scale vector d. +! the parameter nf that humsl passes to calcg is +! iv(nfgcal) = iv(7). if g(x) and h(x) cannot be evaluated, +! then the caller may set iv(nfgcal) to 0, in which case +! humit will return with iv(1) = 65. +! note -- humit overwrites h with the lower triangle +! of diag(d)**-1 * h(x) * diag(d)**-1. +!. +! *** general *** +! +! coded by david m. gay (winter 1980). revised sept. 1982. +! this subroutine was written in connection with research supported +! in part by the national science foundation under grants +! mcs-7600324 and mcs-7906671. +! +! (see sumsl and humsl for references.) +! +!+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++ +! +! *** local variables *** +! + integer :: dg1, dummy, i, j, k, l, lstgst, nn1o2, step1,& + temp1, w1, x01 + real(kind=8) :: t +! +! *** constants *** +! +!el real(kind=8) :: one, onep2, zero +! +! *** no intrinsic functions *** +! +! *** external functions and subroutines *** +! +!el external assst, deflt, dotprd, dupdu, gqtst, itsum, parck, +!el 1 reldst, slvmul, stopx, vaxpy, vcopy, vscopy, v2norm +!el logical stopx +!el real(kind=8) :: dotprd, reldst, v2norm +! +! assst.... assesses candidate step. +! deflt.... provides default iv and v input values. +! dotprd... returns inner product of two vectors. +! dupdu.... updates scale vector d. +! gqtst.... computes optimally locally constrained step. +! itsum.... prints iteration summary and info on initial and final x. +! parck.... checks validity of input iv and v values. +! reldst... computes v(reldx) = relative step size. +! slvmul... multiplies symmetric matrix times vector, given the lower +! triangle of the matrix. +! stopx.... returns .true. if the break key has been pressed. +! vaxpy.... computes scalar times one vector plus another. +! vcopy.... copies one vector to another. +! vscopy... sets all elements of a vector to a scalar. +! v2norm... returns the 2-norm of a vector. +! +! *** subscripts for iv and v *** +! +!el integer cnvcod, dg, dgnorm, dinit, dstnrm, dtinit, dtol, +!el 1 dtype, d0init, f, f0, fdif, gtstep, incfac, irc, kagqt, +!el 2 lmat, lmax0, lmaxs, mode, model, mxfcal, mxiter, nextv, +!el 3 nfcall, nfgcal, ngcall, niter, preduc, radfac, radinc, +!el 4 radius, rad0, reldx, restor, step, stglim, stlstg, stppar, +!el 5 toobig, tuner4, tuner5, vneed, w, xirc, x0 +! +! *** iv subscript values *** +! +!/6 +! data cnvcod/55/, dg/37/, dtol/59/, dtype/16/, irc/29/, kagqt/33/, +! 1 lmat/42/, mode/35/, model/5/, mxfcal/17/, mxiter/18/, +! 2 nextv/47/, nfcall/6/, nfgcal/7/, ngcall/30/, niter/31/, +! 3 radinc/8/, restor/9/, step/40/, stglim/11/, stlstg/41/, +! 4 toobig/2/, vneed/4/, w/34/, xirc/13/, x0/43/ +!/7 + integer,parameter :: cnvcod=55, dg=37, dtol=59, dtype=16, irc=29, kagqt=33,& + lmat=42, mode=35, model=5, mxfcal=17, mxiter=18,& + nextv=47, nfcall=6, nfgcal=7, ngcall=30, niter=31,& + radinc=8, restor=9, step=40, stglim=11, stlstg=41,& + toobig=2, vneed=4, w=34, xirc=13, x0=43 +!/ +! +! *** v subscript values *** +! +!/6 +! data dgnorm/1/, dinit/38/, dstnrm/2/, dtinit/39/, d0init/40/, +! 1 f/10/, f0/13/, fdif/11/, gtstep/4/, incfac/23/, lmax0/35/, +! 2 lmaxs/36/, preduc/7/, radfac/16/, radius/8/, rad0/9/, +! 3 reldx/17/, stppar/5/, tuner4/29/, tuner5/30/ +!/7 + integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dtinit=39, d0init=40,& + f=10, f0=13, fdif=11, gtstep=4, incfac=23, lmax0=35,& + lmaxs=36, preduc=7, radfac=16, radius=8, rad0=9,& + reldx=17, stppar=5, tuner4=29, tuner5=30 +!/ +! +!/6 +! data one/1.d+0/, onep2/1.2d+0/, zero/0.d+0/ +!/7 + real(kind=8),parameter :: one=1.d+0, onep2=1.2d+0, zero=0.d+0 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! + i = iv(1) + if (i .eq. 1) go to 30 + if (i .eq. 2) go to 40 +! +! *** check validity of iv and v input values *** +! + if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v) + if (iv(1) .eq. 12 .or. iv(1) .eq. 13) & + iv(vneed) = iv(vneed) + n*(n+21)/2 + 7 + call parck(2, d, iv, liv, lv, n, v) + i = iv(1) - 2 + if (i .gt. 12) go to 999 + nn1o2 = n * (n + 1) / 2 + if (lh .ge. nn1o2) go to (210,210,210,210,210,210,160,120,160,& + 10,10,20), i + iv(1) = 66 + go to 350 +! +! *** storage allocation *** +! + 10 iv(dtol) = iv(lmat) + nn1o2 + iv(x0) = iv(dtol) + 2*n + iv(step) = iv(x0) + n + iv(stlstg) = iv(step) + n + iv(dg) = iv(stlstg) + n + iv(w) = iv(dg) + n + iv(nextv) = iv(w) + 4*n + 7 + if (iv(1) .ne. 13) go to 20 + iv(1) = 14 + go to 999 +! +! *** initialization *** +! + 20 iv(niter) = 0 + iv(nfcall) = 1 + iv(ngcall) = 1 + iv(nfgcal) = 1 + iv(mode) = -1 + iv(model) = 1 + iv(stglim) = 1 + iv(toobig) = 0 + iv(cnvcod) = 0 + iv(radinc) = 0 + v(rad0) = zero + v(stppar) = zero + if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit)) + k = iv(dtol) + if (v(dtinit) .gt. zero) call vscopy(n, v(k), v(dtinit)) + k = k + n + if (v(d0init) .gt. zero) call vscopy(n, v(k), v(d0init)) + iv(1) = 1 + go to 999 +! + 30 v(f) = fx + if (iv(mode) .ge. 0) go to 210 + iv(1) = 2 + if (iv(toobig) .eq. 0) go to 999 + iv(1) = 63 + go to 350 +! +! *** make sure gradient could be computed *** +! + 40 if (iv(nfgcal) .ne. 0) go to 50 + iv(1) = 65 + go to 350 +! +! *** update the scale vector d *** +! + 50 dg1 = iv(dg) + if (iv(dtype) .le. 0) go to 70 + k = dg1 + j = 0 + do 60 i = 1, n + j = j + i + v(k) = h(j) + k = k + 1 + 60 continue + call dupdu(d, v(dg1), iv, liv, lv, n, v) +! +! *** compute scaled gradient and its norm *** +! + 70 dg1 = iv(dg) + k = dg1 + do 80 i = 1, n + v(k) = g(i) / d(i) + k = k + 1 + 80 continue + v(dgnorm) = v2norm(n, v(dg1)) +! +! *** compute scaled hessian *** +! + k = 1 + do 100 i = 1, n + t = one / d(i) + do 90 j = 1, i + h(k) = t * h(k) / d(j) + k = k + 1 + 90 continue + 100 continue +! + if (iv(cnvcod) .ne. 0) go to 340 + if (iv(mode) .eq. 0) go to 300 +! +! *** allow first step to have scaled 2-norm at most v(lmax0) *** +! + v(radius) = v(lmax0) +! + iv(mode) = 0 +! +! +!----------------------------- main loop ----------------------------- +! +! +! *** print iteration summary, check iteration limit *** +! + 110 call itsum(d, g, iv, liv, lv, n, v, x) + 120 k = iv(niter) + if (k .lt. iv(mxiter)) go to 130 + iv(1) = 10 + go to 350 +! + 130 iv(niter) = k + 1 +! +! *** initialize for start of next iteration *** +! + dg1 = iv(dg) + x01 = iv(x0) + v(f0) = v(f) + iv(irc) = 4 + iv(kagqt) = -1 +! +! *** copy x to x0 *** +! + call vcopy(n, v(x01), x) +! +! *** update radius *** +! + if (k .eq. 0) go to 150 + step1 = iv(step) + k = step1 + do 140 i = 1, n + v(k) = d(i) * v(k) + k = k + 1 + 140 continue + v(radius) = v(radfac) * v2norm(n, v(step1)) +! +! *** check stopx and function evaluation limit *** +! +! AL 4/30/95 + dummy=iv(nfcall) + 150 if (.not. stopx(dummy)) go to 170 + iv(1) = 11 + go to 180 +! +! *** come here when restarting after func. eval. limit or stopx. +! + 160 if (v(f) .ge. v(f0)) go to 170 + v(radfac) = one + k = iv(niter) + go to 130 +! + 170 if (iv(nfcall) .lt. iv(mxfcal)) go to 190 + iv(1) = 9 + 180 if (v(f) .ge. v(f0)) go to 350 +! +! *** in case of stopx or function evaluation limit with +! *** improved v(f), evaluate the gradient at x. +! + iv(cnvcod) = iv(1) + go to 290 +! +!. . . . . . . . . . . . . compute candidate step . . . . . . . . . . +! + 190 step1 = iv(step) + dg1 = iv(dg) + l = iv(lmat) + w1 = iv(w) + call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1)) + if (iv(irc) .eq. 6) go to 210 +! +! *** check whether evaluating f(x0 + step) looks worthwhile *** +! + if (v(dstnrm) .le. zero) go to 210 + if (iv(irc) .ne. 5) go to 200 + if (v(radfac) .le. one) go to 200 + if (v(preduc) .le. onep2 * v(fdif)) go to 210 +! +! *** compute f(x0 + step) *** +! + 200 x01 = iv(x0) + step1 = iv(step) + call vaxpy(n, x, one, v(step1), v(x01)) + iv(nfcall) = iv(nfcall) + 1 + iv(1) = 1 + iv(toobig) = 0 + go to 999 +! +!. . . . . . . . . . . . . assess candidate step . . . . . . . . . . . +! + 210 x01 = iv(x0) + v(reldx) = reldst(n, d, x, v(x01)) + call assst(iv, liv, lv, v) + step1 = iv(step) + lstgst = iv(stlstg) + if (iv(restor) .eq. 1) call vcopy(n, x, v(x01)) + if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1)) + if (iv(restor) .ne. 3) go to 220 + call vcopy(n, v(step1), v(lstgst)) + call vaxpy(n, x, one, v(step1), v(x01)) + v(reldx) = reldst(n, d, x, v(x01)) +! + 220 k = iv(irc) + go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k +! +! *** recompute step with new radius *** +! + 230 v(radius) = v(radfac) * v(dstnrm) + go to 150 +! +! *** compute step of length v(lmaxs) for singular convergence test. +! + 240 v(radius) = v(lmaxs) + go to 190 +! +! *** convergence or false convergence *** +! + 250 iv(cnvcod) = k - 4 + if (v(f) .ge. v(f0)) go to 340 + if (iv(xirc) .eq. 14) go to 340 + iv(xirc) = 14 +! +!. . . . . . . . . . . . process acceptable step . . . . . . . . . . . +! + 260 if (iv(irc) .ne. 3) go to 290 + temp1 = lstgst +! +! *** prepare for gradient tests *** +! *** set temp1 = hessian * step + g(x0) +! *** = diag(d) * (h * step + g(x0)) +! +! use x0 vector as temporary. + k = x01 + do 270 i = 1, n + v(k) = d(i) * v(step1) + k = k + 1 + step1 = step1 + 1 + 270 continue + call slvmul(n, v(temp1), h, v(x01)) + do 280 i = 1, n + v(temp1) = d(i) * v(temp1) + g(i) + temp1 = temp1 + 1 + 280 continue +! +! *** compute gradient and hessian *** +! + 290 iv(ngcall) = iv(ngcall) + 1 + iv(1) = 2 + go to 999 +! + 300 iv(1) = 2 + if (iv(irc) .ne. 3) go to 110 +! +! *** set v(radfac) by gradient tests *** +! + temp1 = iv(stlstg) + step1 = iv(step) +! +! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) *** +! + k = temp1 + do 310 i = 1, n + v(k) = (v(k) - g(i)) / d(i) + k = k + 1 + 310 continue +! +! *** do gradient tests *** +! + if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320 + if (dotprd(n, g, v(step1)) & + .ge. v(gtstep) * v(tuner5)) go to 110 + 320 v(radfac) = v(incfac) + go to 110 +! +!. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . . +! +! *** bad parameters to assess *** +! + 330 iv(1) = 64 + go to 350 +! +! *** print summary of final iteration and other requested items *** +! + 340 iv(1) = iv(cnvcod) + iv(cnvcod) = 0 + 350 call itsum(d, g, iv, liv, lv, n, v, x) +! + 999 return +! +! *** last card of humit follows *** + end subroutine humit +!----------------------------------------------------------------------------- + subroutine dupdu(d, hdiag, iv, liv, lv, n, v) +! +! *** update scale vector d for humsl *** +! +! *** parameter declarations *** +! + integer :: liv, lv, n + integer :: iv(liv) + real(kind=8) :: d(n), hdiag(n), v(lv) +! +! *** local variables *** +! + integer :: dtoli, d0i, i + real(kind=8) :: t, vdfac +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dabs, dmax1, dsqrt +!/ +! *** subscripts for iv and v *** +! +!el integer :: dfac, dtol, dtype, niter +!/6 +! data dfac/41/, dtol/59/, dtype/16/, niter/31/ +!/7 + integer,parameter :: dfac=41, dtol=59, dtype=16, niter=31 +!/ +! +!------------------------------- body -------------------------------- +! + i = iv(dtype) + if (i .eq. 1) go to 10 + if (iv(niter) .gt. 0) go to 999 +! + 10 dtoli = iv(dtol) + d0i = dtoli + n + vdfac = v(dfac) + do 20 i = 1, n + t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i)) + if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i)) + d(i) = t + dtoli = dtoli + 1 + d0i = d0i + 1 + 20 continue +! + 999 return +! *** last card of dupdu follows *** + end subroutine dupdu +!----------------------------------------------------------------------------- + subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w) +! +! *** compute goldfeld-quandt-trotter step by more-hebden technique *** +! *** (nl2sol version 2.2), modified a la more and sorensen *** +! +! *** parameter declarations *** +! + integer :: ka, p +!al real(kind=8) :: d(p), dig(p), dihdi(1), l(1), v(21), step(p), +!al 1 w(1) + real(kind=8) :: d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),& + v(21), step(p),w(4*p+7) +! dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7) +! +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! +! *** purpose *** +! +! given the (compactly stored) lower triangle of a scaled +! hessian (approximation) and a nonzero scaled gradient vector, +! this subroutine computes a goldfeld-quandt-trotter step of +! approximate length v(radius) by the more-hebden technique. in +! other words, step is computed to (approximately) minimize +! psi(step) = (g**t)*step + 0.5*(step**t)*h*step such that the +! 2-norm of d*step is at most (approximately) v(radius), where +! g is the gradient, h is the hessian, and d is a diagonal +! scale matrix whose diagonal is stored in the parameter d. +! (gqtst assumes dig = d**-1 * g and dihdi = d**-1 * h * d**-1.) +! +! *** parameter description *** +! +! d (in) = the scale vector, i.e. the diagonal of the scale +! matrix d mentioned above under purpose. +! dig (in) = the scaled gradient vector, d**-1 * g. if g = 0, then +! step = 0 and v(stppar) = 0 are returned. +! dihdi (in) = lower triangle of the scaled hessian (approximation), +! i.e., d**-1 * h * d**-1, stored compactly by rows., i.e., +! in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc. +! ka (i/o) = the number of hebden iterations (so far) taken to deter- +! mine step. ka .lt. 0 on input means this is the first +! attempt to determine step (for the present dig and dihdi) +! -- ka is initialized to 0 in this case. output with +! ka = 0 (or v(stppar) = 0) means step = -(h**-1)*g. +! l (i/o) = workspace of length p*(p+1)/2 for cholesky factors. +! p (in) = number of parameters -- the hessian is a p x p matrix. +! step (i/o) = the step computed. +! v (i/o) contains various constants and variables described below. +! w (i/o) = workspace of length 4*p + 6. +! +! *** entries in v *** +! +! v(dgnorm) (i/o) = 2-norm of (d**-1)*g. +! v(dstnrm) (output) = 2-norm of d*step. +! v(dst0) (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or +! overestimate of smallest eigenvalue of (d**-1)*h*(d**-1). +! v(epslon) (in) = max. rel. error allowed for psi(step). for the +! step returned, psi(step) will exceed its optimal value +! by less than -v(epslon)*psi(step). suggested value = 0.1. +! v(gtstep) (out) = inner product between g and step. +! v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step) (for pos. def. +! h only -- v(nreduc) is set to zero otherwise). +! v(phmnfc) (in) = tol. (together with v(phmxfc)) for accepting step +! (more*s sigma). the error v(dstnrm) - v(radius) must lie +! between v(phmnfc)*v(radius) and v(phmxfc)*v(radius). +! v(phmxfc) (in) (see v(phmnfc).) +! suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5. +! v(preduc) (out) = psi(step) = predicted obj. func. reduction for step. +! v(radius) (in) = radius of current (scaled) trust region. +! v(rad0) (i/o) = value of v(radius) from previous call. +! v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha +! described below under algorithm notes. if h + alpha*d**2 +! (see algorithm notes) is (nearly) singular, however, +! then v(stppar) = -alpha. +! +! *** usage notes *** +! +! if it is desired to recompute step using a different value of +! v(radius), then this routine may be restarted by calling it +! with all parameters unchanged except v(radius). (this explains +! why step and w are listed as i/o). on an initial call (one with +! ka .lt. 0), step and w need not be initialized and only compo- +! nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and +! v(rad0) of v must be initialized. +! +! *** algorithm notes *** +! +! the desired g-q-t step (ref. 2, 3, 4, 6) satisfies +! (h + alpha*d**2)*step = -g for some nonnegative alpha such that +! h + alpha*d**2 is positive semidefinite. alpha and step are +! computed by a scheme analogous to the one described in ref. 5. +! estimates of the smallest and largest eigenvalues of the hessian +! are obtained from the gerschgorin circle theorem enhanced by a +! simple form of the scaling described in ref. 7. cases in which +! h + alpha*d**2 is nearly (or exactly) singular are handled by +! the technique discussed in ref. 2. in these cases, a step of +! (exact) length v(radius) is returned for which psi(step) exceeds +! its optimal value by less than -v(epslon)*psi(step). the test +! suggested in ref. 6 for detecting the special case is performed +! once two matrix factorizations have been done -- doing so sooner +! seems to degrade the performance of optimization routines that +! call this routine. +! +! *** functions and subroutines called *** +! +! dotprd - returns inner product of two vectors. +! litvmu - applies inverse-transpose of compact lower triang. matrix. +! livmul - applies inverse of compact lower triang. matrix. +! lsqrt - finds cholesky factor (of compactly stored lower triang.). +! lsvmin - returns approx. to min. sing. value of lower triang. matrix. +! rmdcon - returns machine-dependent constants. +! v2norm - returns 2-norm of a vector. +! +! *** references *** +! +! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive +! nonlinear least-squares algorithm, acm trans. math. +! software, vol. 7, no. 3. +! 2. gay, d.m. (1981), computing optimal locally constrained steps, +! siam j. sci. statist. computing, vol. 2, no. 2, pp. +! 186-197. +! 3. goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966), +! maximization by quadratic hill-climbing, econometrica 34, +! pp. 541-551. +! 4. hebden, m.d. (1973), an algorithm for minimization using exact +! second derivatives, report t.p. 515, theoretical physics +! div., a.e.r.e. harwell, oxon., england. +! 5. more, j.j. (1978), the levenberg-marquardt algorithm, implemen- +! tation and theory, pp.105-116 of springer lecture notes +! in mathematics no. 630, edited by g.a. watson, springer- +! verlag, berlin and new york. +! 6. more, j.j., and sorensen, d.c. (1981), computing a trust region +! step, technical report anl-81-83, argonne national lab. +! 7. varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15, +! pp. 719-729. +! +! *** general *** +! +! coded by david m. gay. +! this subroutine was written in connection with research +! supported by the national science foundation under grants +! mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and +! mcs-7906671. +! +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! +! *** local variables *** +! + logical :: restrt + integer :: dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,& + j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x + real(kind=8) :: alphak, aki, akk, delta, dst, eps, gtsta, lk,& + oldphi, phi, phimax, phimin, psifac, rad, radsq,& + root, si, sk, sw, t, twopsi, t1, t2, uk, wi +! +! *** constants *** + real(kind=8) :: big, dgxfac !el, epsfac, four, half, kappa, negone, +!el 1 one, p001, six, three, two, zero +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dabs, dmax1, dmin1, dsqrt +!/ +! *** external functions and subroutines *** +! +!el external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm +!el real(kind=8) :: dotprd, lsvmin, rmdcon, v2norm +! +! *** subscripts for v *** +! +!el integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc, +!el 1 phmnfc, phmxfc, preduc, radius, rad0 +!/6 +! data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/, +! 1 nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/, +! 2 rad0/9/, stppar/5/ +!/7 + integer,parameter :: dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,& + nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,& + rad0=9, stppar=5 +!/ +! +!/6 +! data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/, +! 1 kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/, +! 2 six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/ +!/7 + real(kind=8), parameter :: epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,& + kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,& + six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0 + save dgxfac +!/ + data big/0.d+0/, dgxfac/0.d+0/ +! +! *** body *** +! +! *** store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx). + dggdmx = p + 1 +! *** store gerschgorin over- and underestimates of the largest +! *** and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax) +! *** and w(emin) respectively. + emax = dggdmx + 1 + emin = emax + 1 +! *** for use in recomputing step, the final values of lk, uk, dst, +! *** and the inverse derivative of more*s phi at 0 (for pos. def. +! *** h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin) +! *** respectively. + lk0 = emin + 1 + phipin = lk0 + 1 + uk0 = phipin + 1 + dstsav = uk0 + 1 +! *** store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p). + diag0 = dstsav + diag = diag0 + 1 +! *** store -d*step in w(q),...,w(q0+p). + q0 = diag0 + p + q = q0 + 1 +! *** allocate storage for scratch vector x *** + x = q + p + rad = v(radius) + radsq = rad**2 +! *** phitol = max. error allowed in dst = v(dstnrm) = 2-norm of +! *** d*step. + phimax = v(phmxfc) * rad + phimin = v(phmnfc) * rad + psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) * & + (kappa + one) + kappa + two) * rad**2) +! *** oldphi is used to detect limits of numerical accuracy. if +! *** we recompute step and it does not change, then we accept it. + oldphi = zero + eps = v(epslon) + irc = 0 + restrt = .false. + kalim = ka + 50 +! +! *** start or restart, depending on ka *** +! + if (ka .ge. 0) go to 290 +! +! *** fresh start *** +! + k = 0 + uk = negone + ka = 0 + kalim = 50 + v(dgnorm) = v2norm(p, dig) + v(nreduc) = zero + v(dst0) = zero + kamin = 3 + if (v(dgnorm) .eq. zero) kamin = 0 +! +! *** store diag(dihdi) in w(diag0+1),...,w(diag0+p) *** +! + j = 0 + do 10 i = 1, p + j = j + i + k1 = diag0 + i + w(k1) = dihdi(j) + 10 continue +! +! *** determine w(dggdmx), the largest element of dihdi *** +! + t1 = zero + j = p * (p + 1) / 2 + do 20 i = 1, j + t = dabs(dihdi(i)) + if (t1 .lt. t) t1 = t + 20 continue + w(dggdmx) = t1 +! +! *** try alpha = 0 *** +! + 30 call lsqrt(1, p, l, dihdi, irc) + if (irc .eq. 0) go to 50 +! *** indef. h -- underestimate smallest eigenvalue, use this +! *** estimate to initialize lower bound lk on alpha. + j = irc*(irc+1)/2 + t = l(j) + l(j) = one + do 40 i = 1, irc + 40 w(i) = zero + w(irc) = one + call litvmu(irc, w, l, w) + t1 = v2norm(irc, w) + lk = -t / t1 / t1 + v(dst0) = -lk + if (restrt) go to 210 + go to 70 +! +! *** positive definite h -- compute unmodified newton step. *** + 50 lk = zero + t = lsvmin(p, l, w(q), w(q)) + if (t .ge. one) go to 60 + if (big .le. zero) big = rmdcon(6) + if (v(dgnorm) .ge. t*t*big) go to 70 + 60 call livmul(p, w(q), l, dig) + gtsta = dotprd(p, w(q), w(q)) + v(nreduc) = half * gtsta + call litvmu(p, w(q), l, w(q)) + dst = v2norm(p, w(q)) + v(dst0) = dst + phi = dst - rad + if (phi .le. phimax) go to 260 + if (restrt) go to 210 +! +! *** prepare to compute gerschgorin estimates of largest (and +! *** smallest) eigenvalues. *** +! + 70 k = 0 + do 100 i = 1, p + wi = zero + if (i .eq. 1) go to 90 + im1 = i - 1 + do 80 j = 1, im1 + k = k + 1 + t = dabs(dihdi(k)) + wi = wi + t + w(j) = w(j) + t + 80 continue + 90 w(i) = wi + k = k + 1 + 100 continue +! +! *** (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1) *** +! + k = 1 + t1 = w(diag) - w(1) + if (p .le. 1) go to 120 + do 110 i = 2, p + j = diag0 + i + t = w(j) - w(i) + if (t .ge. t1) go to 110 + t1 = t + k = i + 110 continue +! + 120 sk = w(k) + j = diag0 + k + akk = w(j) + k1 = k*(k-1)/2 + 1 + inc = 1 + t = zero + do 150 i = 1, p + if (i .eq. k) go to 130 + aki = dabs(dihdi(k1)) + si = w(i) + j = diag0 + i + t1 = half * (akk - w(j) + si - aki) + t1 = t1 + dsqrt(t1*t1 + sk*aki) + if (t .lt. t1) t = t1 + if (i .lt. k) go to 140 + 130 inc = i + 140 k1 = k1 + inc + 150 continue +! + w(emin) = akk - t + uk = v(dgnorm)/rad - w(emin) + if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk + if (uk .le. zero) uk = p001 +! +! *** compute gerschgorin (over-)estimate of largest eigenvalue *** +! + k = 1 + t1 = w(diag) + w(1) + if (p .le. 1) go to 170 + do 160 i = 2, p + j = diag0 + i + t = w(j) + w(i) + if (t .le. t1) go to 160 + t1 = t + k = i + 160 continue +! + 170 sk = w(k) + j = diag0 + k + akk = w(j) + k1 = k*(k-1)/2 + 1 + inc = 1 + t = zero + do 200 i = 1, p + if (i .eq. k) go to 180 + aki = dabs(dihdi(k1)) + si = w(i) + j = diag0 + i + t1 = half * (w(j) + si - aki - akk) + t1 = t1 + dsqrt(t1*t1 + sk*aki) + if (t .lt. t1) t = t1 + if (i .lt. k) go to 190 + 180 inc = i + 190 k1 = k1 + inc + 200 continue +! + w(emax) = akk + t + lk = dmax1(lk, v(dgnorm)/rad - w(emax)) +! +! *** alphak = current value of alpha (see alg. notes above). we +! *** use more*s scheme for initializing it. + alphak = dabs(v(stppar)) * v(rad0)/rad +! + if (irc .ne. 0) go to 210 +! +! *** compute l0 for positive definite h *** +! + call livmul(p, w, l, w(q)) + t = v2norm(p, w) + w(phipin) = dst / t / t + lk = dmax1(lk, phi*w(phipin)) +! +! *** safeguard alphak and add alphak*i to (d**-1)*h*(d**-1) *** +! + 210 ka = ka + 1 + if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk) & + alphak = uk * dmax1(p001, dsqrt(lk/uk)) + if (alphak .le. zero) alphak = half * uk + if (alphak .le. zero) alphak = uk + k = 0 + do 220 i = 1, p + k = k + i + j = diag0 + i + dihdi(k) = w(j) + alphak + 220 continue +! +! *** try computing cholesky decomposition *** +! + call lsqrt(1, p, l, dihdi, irc) + if (irc .eq. 0) go to 240 +! +! *** (d**-1)*h*(d**-1) + alphak*i is indefinite -- overestimate +! *** smallest eigenvalue for use in updating lk *** +! + j = (irc*(irc+1))/2 + t = l(j) + l(j) = one + do 230 i = 1, irc + 230 w(i) = zero + w(irc) = one + call litvmu(irc, w, l, w) + t1 = v2norm(irc, w) + lk = alphak - t/t1/t1 + v(dst0) = -lk + go to 210 +! +! *** alphak makes (d**-1)*h*(d**-1) positive definite. +! *** compute q = -d*step, check for convergence. *** +! + 240 call livmul(p, w(q), l, dig) + gtsta = dotprd(p, w(q), w(q)) + call litvmu(p, w(q), l, w(q)) + dst = v2norm(p, w(q)) + phi = dst - rad + if (phi .le. phimax .and. phi .ge. phimin) go to 270 + if (phi .eq. oldphi) go to 270 + oldphi = phi + if (phi .lt. zero) go to 330 +! +! *** unacceptable alphak -- update lk, uk, alphak *** +! + 250 if (ka .ge. kalim) go to 270 +! *** the following dmin1 is necessary because of restarts *** + if (phi .lt. zero) uk = dmin1(uk, alphak) +! *** kamin = 0 only iff the gradient vanishes *** + if (kamin .eq. 0) go to 210 + call livmul(p, w, l, w(q)) + t1 = v2norm(p, w) + alphak = alphak + (phi/t1) * (dst/t1) * (dst/rad) + lk = dmax1(lk, alphak) + go to 210 +! +! *** acceptable step on first try *** +! + 260 alphak = zero +! +! *** successful step in general. compute step = -(d**-1)*q *** +! + 270 do 280 i = 1, p + j = q0 + i + step(i) = -w(j)/d(i) + 280 continue + v(gtstep) = -gtsta + v(preduc) = half * (dabs(alphak)*dst*dst + gtsta) + go to 410 +! +! +! *** restart with new radius *** +! + 290 if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310 +! +! *** prepare to return newton step *** +! + restrt = .true. + ka = ka + 1 + k = 0 + do 300 i = 1, p + k = k + i + j = diag0 + i + dihdi(k) = w(j) + 300 continue + uk = negone + go to 30 +! + 310 kamin = ka + 3 + if (v(dgnorm) .eq. zero) kamin = 0 + if (ka .eq. 0) go to 50 +! + dst = w(dstsav) + alphak = dabs(v(stppar)) + phi = dst - rad + t = v(dgnorm)/rad + uk = t - w(emin) + if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk + if (uk .le. zero) uk = p001 + if (rad .gt. v(rad0)) go to 320 +! +! *** smaller radius *** + lk = zero + if (alphak .gt. zero) lk = w(lk0) + lk = dmax1(lk, t - w(emax)) + if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin)) + go to 250 +! +! *** bigger radius *** + 320 if (alphak .gt. zero) uk = dmin1(uk, w(uk0)) + lk = dmax1(zero, -v(dst0), t - w(emax)) + if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin)) + go to 250 +! +! *** decide whether to check for special case... in practice (from +! *** the standpoint of the calling optimization code) it seems best +! *** not to check until a few iterations have failed -- hence the +! *** test on kamin below. +! + 330 delta = alphak + dmin1(zero, v(dst0)) + twopsi = alphak*dst*dst + gtsta + if (ka .ge. kamin) go to 340 +! *** if the test in ref. 2 is satisfied, fall through to handle +! *** the special case (as soon as the more-sorensen test detects +! *** it). + if (delta .ge. psifac*twopsi) go to 370 +! +! *** check for the special case of h + alpha*d**2 (nearly) +! *** singular. use one step of inverse power method with start +! *** from lsvmin to obtain approximate eigenvector corresponding +! *** to smallest eigenvalue of (d**-1)*h*(d**-1). lsvmin returns +! *** x and w with l*w = x. +! + 340 t = lsvmin(p, l, w(x), w) +! +! *** normalize w *** + do 350 i = 1, p + 350 w(i) = t*w(i) +! *** complete current inv. power iter. -- replace w by (l**-t)*w. + call litvmu(p, w, l, w) + t2 = one/v2norm(p, w) + do 360 i = 1, p + 360 w(i) = t2*w(i) + t = t2 * t +! +! *** now w is the desired approximate (unit) eigenvector and +! *** t*x = ((d**-1)*h*(d**-1) + alphak*i)*w. +! + sw = dotprd(p, w(q), w) + t1 = (rad + dst) * (rad - dst) + root = dsqrt(sw*sw + t1) + if (sw .lt. zero) root = -root + si = t1 / (sw + root) +! +! *** the actual test for the special case... +! + if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380 +! +! *** update upper bound on smallest eigenvalue (when not positive) +! *** (as recommended by more and sorensen) and continue... +! + if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak) + lk = dmax1(lk, -v(dst0)) +! +! *** check whether we can hope to detect the special case in +! *** the available arithmetic. accept step as it is if not. +! +! *** if not yet available, obtain machine dependent value dgxfac. + 370 if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3) +! + if (delta .gt. dgxfac*w(dggdmx)) go to 250 + go to 270 +! +! *** special case detected... negate alphak to indicate special case +! + 380 alphak = -alphak + v(preduc) = half * twopsi +! +! *** accept current step if adding si*w would lead to a +! *** further relative reduction in psi of less than v(epslon)/3. +! + t1 = zero + t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w))) + if (t .lt. eps*twopsi/six) go to 390 + v(preduc) = v(preduc) + t + dst = rad + t1 = -si + 390 do 400 i = 1, p + j = q0 + i + w(j) = t1*w(i) - w(j) + step(i) = w(j) / d(i) + 400 continue + v(gtstep) = dotprd(p, dig, w(q)) +! +! *** save values for use in a possible restart *** +! + 410 v(dstnrm) = dst + v(stppar) = alphak + w(lk0) = lk + w(uk0) = uk + v(rad0) = rad + w(dstsav) = dst +! +! *** restore diagonal of dihdi *** +! + j = 0 + do 420 i = 1, p + j = j + i + k = diag0 + i + dihdi(j) = w(k) + 420 continue +! + 999 return +! +! *** last card of gqtst follows *** + end subroutine gqtst +!----------------------------------------------------------------------------- + subroutine lsqrt(n1, n, l, a, irc) +! +! *** compute rows n1 through n of the cholesky factor l of +! *** a = l*(l**t), where l and the lower triangle of a are both +! *** stored compactly by rows (and may occupy the same storage). +! *** irc = 0 means all went well. irc = j means the leading +! *** principal j x j submatrix of a is not positive definite -- +! *** and l(j*(j+1)/2) contains the (nonpos.) reduced j-th diagonal. +! +! *** parameters *** +! + integer :: n1, n, irc +!al real(kind=8) :: l(1), a(1) + real(kind=8) :: l(n*(n+1)/2), a(n*(n+1)/2) +! dimension l(n*(n+1)/2), a(n*(n+1)/2) +! +! *** local variables *** +! + integer :: i, ij, ik, im1, i0, j, jk, jm1, j0, k + real(kind=8) :: t, td !el, zero +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dsqrt +!/ +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +! +! *** body *** +! + i0 = n1 * (n1 - 1) / 2 + do 50 i = n1, n + td = zero + if (i .eq. 1) go to 40 + j0 = 0 + im1 = i - 1 + do 30 j = 1, im1 + t = zero + if (j .eq. 1) go to 20 + jm1 = j - 1 + do 10 k = 1, jm1 + ik = i0 + k + jk = j0 + k + t = t + l(ik)*l(jk) + 10 continue + 20 ij = i0 + j + j0 = j0 + j + t = (a(ij) - t) / l(j0) + l(ij) = t + td = td + t*t + 30 continue + 40 i0 = i0 + i + t = a(i0) - td + if (t .le. zero) go to 60 + l(i0) = dsqrt(t) + 50 continue +! + irc = 0 + go to 999 +! + 60 l(i0) = t + irc = i +! + 999 return +! +! *** last card of lsqrt *** + end subroutine lsqrt +!----------------------------------------------------------------------------- + real(kind=8) function lsvmin(p, l, x, y) +! +! *** estimate smallest sing. value of packed lower triang. matrix l +! +! *** parameter declarations *** +! + integer :: p +!al real(kind=8) :: l(1), x(p), y(p) + real(kind=8) :: l(p*(p+1)/2), x(p), y(p) +! dimension l(p*(p+1)/2) +! +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! +! *** purpose *** +! +! this function returns a good over-estimate of the smallest +! singular value of the packed lower triangular matrix l. +! +! *** parameter description *** +! +! p (in) = the order of l. l is a p x p lower triangular matrix. +! l (in) = array holding the elements of l in row order, i.e. +! l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc. +! x (out) if lsvmin returns a positive value, then x is a normalized +! approximate left singular vector corresponding to the +! smallest singular value. this approximation may be very +! crude. if lsvmin returns zero, then some components of x +! are zero and the rest retain their input values. +! y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an +! unnormalized approximate right singular vector correspond- +! ing to the smallest singular value. this approximation +! may be crude. if lsvmin returns zero, then y retains its +! input value. the caller may pass the same vector for x +! and y (nonstandard fortran usage), in which case y over- +! writes x (for nonzero lsvmin returns). +! +! *** algorithm notes *** +! +! the algorithm is based on (1), with the additional provision that +! lsvmin = 0 is returned if the smallest diagonal element of l +! (in magnitude) is not more than the unit roundoff times the +! largest. the algorithm uses a random number generator proposed +! in (4), which passes the spectral test with flying colors -- see +! (2) and (3). +! +! *** subroutines and functions called *** +! +! v2norm - function, returns the 2-norm of a vector. +! +! *** references *** +! +! (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977), +! an estimate for the condition number of a matrix, report +! tm-310, applied math. div., argonne national laboratory. +! +! (2) hoaglin, d.c. (1976), theoretical properties of congruential +! random-number generators -- an empirical view, +! memorandum ns-340, dept. of statistics, harvard univ. +! +! (3) knuth, d.e. (1969), the art of computer programming, vol. 2 +! (seminumerical algorithms), addison-wesley, reading, mass. +! +! (4) smith, c.s. (1971), multiplicative pseudo-random number +! generators with prime modulus, j. assoc. comput. mach. 18, +! pp. 586-593. +! +! *** history *** +! +! designed and coded by david m. gay (winter 1977/summer 1978). +! +! *** general *** +! +! this subroutine was written in connection with research +! supported by the national science foundation under grants +! mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989. +! +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! +! *** local variables *** +! + integer :: i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1 + real(kind=8) :: b, sminus, splus, t, xminus, xplus +! +! *** constants *** +! +!el real(kind=8) :: half, one, r9973, zero +! +! *** intrinsic functions *** +!/+ +!el integer mod +!el real float +!el real(kind=8) :: dabs +!/ +! *** external functions and subroutines *** +! +!el external dotprd, v2norm, vaxpy +!el real(kind=8) :: dotprd, v2norm +! +!/6 +! data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/ +!/7 + real(kind=8),parameter :: half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0 +!/ +! +! *** body *** +! + ix = 2 + pm1 = p - 1 +! +! *** first check whether to return lsvmin = 0 and initialize x *** +! + ii = 0 + j0 = p*pm1/2 + jj = j0 + p + if (l(jj) .eq. zero) go to 110 + ix = mod(3432*ix, 9973) + b = half*(one + float(ix)/r9973) + xplus = b / l(jj) + x(p) = xplus + if (p .le. 1) go to 60 + do 10 i = 1, pm1 + ii = ii + i + if (l(ii) .eq. zero) go to 110 + ji = j0 + i + x(i) = xplus * l(ji) + 10 continue +! +! *** solve (l**t)*x = b, where the components of b have randomly +! *** chosen magnitudes in (.5,1) with signs chosen to make x large. +! +! do j = p-1 to 1 by -1... + do 50 jjj = 1, pm1 + j = p - jjj +! *** determine x(j) in this iteration. note for i = 1,2,...,j +! *** that x(i) holds the current partial sum for row i. + ix = mod(3432*ix, 9973) + b = half*(one + float(ix)/r9973) + xplus = (b - x(j)) + xminus = (-b - x(j)) + splus = dabs(xplus) + sminus = dabs(xminus) + jm1 = j - 1 + j0 = j*jm1/2 + jj = j0 + j + xplus = xplus/l(jj) + xminus = xminus/l(jj) + if (jm1 .eq. 0) go to 30 + do 20 i = 1, jm1 + ji = j0 + i + splus = splus + dabs(x(i) + l(ji)*xplus) + sminus = sminus + dabs(x(i) + l(ji)*xminus) + 20 continue + 30 if (sminus .gt. splus) xplus = xminus + x(j) = xplus +! *** update partial sums *** + if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x) + 50 continue +! +! *** normalize x *** +! + 60 t = one/v2norm(p, x) + do 70 i = 1, p + 70 x(i) = t*x(i) +! +! *** solve l*y = x and return lsvmin = 1/twonorm(y) *** +! + do 100 j = 1, p + jm1 = j - 1 + j0 = j*jm1/2 + jj = j0 + j + t = zero + if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y) + y(j) = (x(j) - t) / l(jj) + 100 continue +! + lsvmin = one/v2norm(p, y) + go to 999 +! + 110 lsvmin = zero + 999 return +! *** last card of lsvmin follows *** + end function lsvmin +!----------------------------------------------------------------------------- + subroutine slvmul(p, y, s, x) +! +! *** set y = s * x, s = p x p symmetric matrix. *** +! *** lower triangle of s stored rowwise. *** +! +! *** parameter declarations *** +! + integer :: p +!al real(kind=8) :: s(1), x(p), y(p) + real(kind=8) :: s(p*(p+1)/2), x(p), y(p) +! dimension s(p*(p+1)/2) +! +! *** local variables *** +! + integer :: i, im1, j, k + real(kind=8) :: xi +! +! *** no intrinsic functions *** +! +! *** external function *** +! +!el external dotprd +!el real(kind=8) :: dotprd +! +!----------------------------------------------------------------------- +! + j = 1 + do 10 i = 1, p + y(i) = dotprd(i, s(j), x) + j = j + i + 10 continue +! + if (p .le. 1) go to 999 + j = 1 + do 40 i = 2, p + xi = x(i) + im1 = i - 1 + j = j + 1 + do 30 k = 1, im1 + y(k) = y(k) + s(j)*xi + j = j + 1 + 30 continue + 40 continue +! + 999 return +! *** last card of slvmul follows *** + end subroutine slvmul +!----------------------------------------------------------------------------- +! minimize_p.F +!----------------------------------------------------------------------------- + subroutine minimize(etot,x,iretcode,nfun) + + use energy, only: func,gradient,fdum!,etotal,enerprint + use comm_srutu +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' + integer,parameter :: liv=60 +! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) +!******************************************************************** +! OPTIMIZE sets up SUMSL or DFP and provides a simple interface for * +! the calling subprogram. * +! when d(i)=1.0, then v(35) is the length of the initial step, * +! calculated in the usual pythagorean way. * +! absolute convergence occurs when the function is within v(31) of * +! zero. unless you know the minimum value in advance, abs convg * +! is probably not useful. * +! relative convergence is when the model predicts that the function * +! will decrease by less than v(32)*abs(fun). * +!******************************************************************** +! include 'COMMON.IOUNITS' +! include 'COMMON.VAR' +! include 'COMMON.GEO' +! include 'COMMON.MINIM' + integer :: i +!el common /srutu/ icall + integer,dimension(liv) :: iv + real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv) +!el real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) + real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) + real(kind=8) :: energia(0:n_ene) +! external func,gradient,fdum +! external func_restr,grad_restr + logical :: not_done,change,reduce +!el common /przechowalnia/ v +!el local variables + integer :: iretcode,nfun,lv,nvar_restr,idum(1),j + real(kind=8) :: etot,rdum(1) !,fdum + + lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) + + if (.not.allocated(v)) allocate(v(1:lv)) + + icall = 1 + + NOT_DONE=.TRUE. + +! DO WHILE (NOT_DONE) + + call deflt(2,iv,liv,lv,v) +! 12 means fresh start, dont call deflt + iv(1)=12 +! max num of fun calls + if (maxfun.eq.0) maxfun=500 + iv(17)=maxfun +! max num of iterations + if (maxmin.eq.0) maxmin=1000 + iv(18)=maxmin +! controls output + iv(19)=2 +! selects output unit + iv(21)=0 + if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout +! 1 means to print out result + iv(22)=print_min_res +! 1 means to print out summary stats + iv(23)=print_min_stat +! 1 means to print initial x and d + iv(24)=print_min_ini +! min val for v(radfac) default is 0.1 + v(24)=0.1D0 +! max val for v(radfac) default is 4.0 + v(25)=2.0D0 +! v(25)=4.0D0 +! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +! the sumsl default is 0.1 + v(26)=0.1D0 +! false conv if (act fnctn decrease) .lt. v(34) +! the sumsl default is 100*machep + v(34)=v(34)/100.0D0 +! absolute convergence + if (tolf.eq.0.0D0) tolf=1.0D-4 + v(31)=tolf +! relative convergence + if (rtolf.eq.0.0D0) rtolf=1.0D-4 + v(32)=rtolf +! controls initial step size + v(35)=1.0D-1 +! large vals of d correspond to small components of step + do i=1,nphi + d(i)=1.0D-1 + enddo + do i=nphi+1,nvar + d(i)=1.0D-1 + enddo +!d print *,'Calling SUMSL' +! call var_to_geom(nvar,x) +! call chainbuild +! call etotal(energia(0)) +! etot = energia(0) +!elmask_r=.true. + IF (mask_r) THEN + call x2xx(x,xx,nvar_restr) + call sumsl(nvar_restr,d,xx,func_restr,grad_restr,& + iv,liv,lv,v,idum,rdum,fdum) + call xx2x(x,xx) + ELSE + call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum) + ENDIF + etot=v(10) + iretcode=iv(1) +!d print *,'Exit SUMSL; return code:',iretcode,' energy:',etot +!d write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1) +! call intout +! change=reduce(x) + call var_to_geom(nvar,x) +! if (change) then +! write (iout,'(a)') 'Reduction worked, minimizing again...' +! else +! not_done=.false. +! endif + call chainbuild + +!el--------------------- +! write (iout,'(/a)') & +! "Cartesian coordinates of the reference structure after SUMSL" +! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & +! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" +! do i=1,nres +! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') & +! restyp(itype(i)),i,(c(j,i),j=1,3),& +! (c(j,i+nres),j=1,3) +! enddo +!el---------------------------- +! call etotal(energia) !sp +! etot=energia(0) +! call enerprint(energia) !sp + nfun=iv(6) + +! write (*,*) 'Processor',MyID,' leaves MINIMIZE.' + +! ENDDO ! NOT_DONE + + return + end subroutine minimize +!----------------------------------------------------------------------------- +! gradient_p.F +!----------------------------------------------------------------------------- + subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm) + + use energy, only: cartder,zerograd,etotal,sum_gradient +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.VAR' +! include 'COMMON.INTERACT' +! include 'COMMON.FFIELD' +! include 'COMMON.IOUNITS' +!EL external ufparm + integer :: uiparm(1) + real(kind=8) :: urparm(1) + real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres) + integer :: n,nf,ig,ind,i,j,ij,k,igall + real(kind=8) :: f,gphii,gthetai,galphai,gomegai + real(kind=8),external :: ufparm + + icg=mod(nf,2)+1 + if (nf-nfl+1) 20,30,40 + 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm) +! write (iout,*) 'grad 20' + if (nf.eq.0) return + goto 40 + 30 continue +#ifdef OSF +! Intercept NaNs in the coordinates +! write(iout,*) (var(i),i=1,nvar) + x_sum=0.D0 + do i=1,n + x_sum=x_sum+x(i) + enddo + if (x_sum.ne.x_sum) then + write(iout,*)" *** grad_restr : Found NaN in coordinates" + call flush(iout) + print *," *** grad_restr : Found NaN in coordinates" + return + endif +#endif + call var_to_geom_restr(n,x) + call chainbuild +! +! Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +! + 40 call cartder +! +! Convert the Cartesian gradient into internal-coordinate gradient. +! + + ig=0 + ind=nres-2 + do i=2,nres-2 + IF (mask_phi(i+2).eq.1) THEN + gphii=0.0D0 + do j=i+1,nres-1 + ind=ind+1 + do k=1,3 + gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) + gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) + enddo + enddo + ig=ig+1 + g(ig)=gphii + ELSE + ind=ind+nres-1-i + ENDIF + enddo + + + ind=0 + do i=1,nres-2 + IF (mask_theta(i+2).eq.1) THEN + ig=ig+1 + gthetai=0.0D0 + do j=i+1,nres-1 + ind=ind+1 + do k=1,3 + gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) + gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) + enddo + enddo + g(ig)=gthetai + ELSE + ind=ind+nres-1-i + ENDIF + enddo + + do i=2,nres-1 + if (itype(i).ne.10) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + galphai=0.0D0 + do k=1,3 + galphai=galphai+dxds(k,i)*gradx(k,i,icg) + enddo + g(ig)=galphai + ENDIF + endif + enddo + + + do i=2,nres-1 + if (itype(i).ne.10) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + gomegai=0.0D0 + do k=1,3 + gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) + enddo + g(ig)=gomegai + ENDIF + endif + enddo + +! +! Add the components corresponding to local energy terms. +! + + ig=0 + igall=0 + do i=4,nres + igall=igall+1 + if (mask_phi(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do i=3,nres + igall=igall+1 + if (mask_theta(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do ij=1,2 + do i=2,nres-1 + if (itype(i).ne.10) then + igall=igall+1 + if (mask_side(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + endif + enddo + enddo + +!d do i=1,ig +!d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) +!d enddo + return + end subroutine grad_restr +!----------------------------------------------------------------------------- + subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F + + use comm_chu + use energy, only: zerograd,etotal,sum_gradient +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.DERIV' +! include 'COMMON.IOUNITS' +! include 'COMMON.GEO' + integer :: n,nf +!el integer :: jjj +!el common /chuju/ jjj + real(kind=8) :: energia(0:n_ene) + real(kind=8) :: f + real(kind=8),external :: ufparm + integer :: uiparm(1) + real(kind=8) :: urparm(1) + real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres) +! if (jjj.gt.0) then +! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n) +! endif + nfl=nf + icg=mod(nf,2)+1 + call var_to_geom_restr(n,x) + call zerograd + call chainbuild +!d write (iout,*) 'ETOTAL called from FUNC' + call etotal(energia) + call sum_gradient + f=energia(0) +! if (jjj.gt.0) then +! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n) +! write (iout,*) 'f=',etot +! jjj=0 +! endif + return + end subroutine func_restr +!----------------------------------------------------------------------------- +! subroutine func(n,x,nf,f,uiparm,urparm,ufparm) in module energy +!----------------------------------------------------------------------------- + subroutine x2xx(x,xx,n) + +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.VAR' +! include 'COMMON.CHAIN' +! include 'COMMON.INTERACT' + integer :: n,i,ij,ig,igall + real(kind=8),dimension(6*nres) :: xx,x !(maxvar) (maxvar=6*maxres) + +!el allocate(varall(nvar)) allocated in alioc_ener_arrays + + do i=1,nvar + varall(i)=x(i) + enddo + + ig=0 + igall=0 + do i=4,nres + igall=igall+1 + if (mask_phi(i).eq.1) then + ig=ig+1 + xx(ig)=x(igall) + endif + enddo + + do i=3,nres + igall=igall+1 + if (mask_theta(i).eq.1) then + ig=ig+1 + xx(ig)=x(igall) + endif + enddo + + do ij=1,2 + do i=2,nres-1 + if (itype(i).ne.10) then + igall=igall+1 + if (mask_side(i).eq.1) then + ig=ig+1 + xx(ig)=x(igall) + endif + endif + enddo + enddo + + n=ig + + return + end subroutine x2xx +!----------------------------------------------------------------------------- +!el subroutine xx2x(x,xx) in module math +!----------------------------------------------------------------------------- + subroutine minim_dc(etot,iretcode,nfun) + + use MPI_data + use energy, only: fdum,check_ecartint +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + integer,parameter :: liv=60 +! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) +! include 'COMMON.SETUP' +! include 'COMMON.IOUNITS' +! include 'COMMON.VAR' +! include 'COMMON.GEO' +! include 'COMMON.MINIM' +! include 'COMMON.CHAIN' + integer :: iretcode,nfun,k,i,j,lv,idum(1) + integer,dimension(liv) :: iv + real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv) + real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) +!el common /przechowalnia/ v + + real(kind=8) :: energia(0:n_ene) +! external func_dc,grad_dc ,fdum + logical :: not_done,change,reduce + real(kind=8) :: g(6*nres),f1,etot,rdum(1) !,fdum + + lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) + + if (.not. allocated(v)) allocate(v(1:lv)) + + call deflt(2,iv,liv,lv,v) +! 12 means fresh start, dont call deflt + iv(1)=12 +! max num of fun calls + if (maxfun.eq.0) maxfun=500 + iv(17)=maxfun +! max num of iterations + if (maxmin.eq.0) maxmin=1000 + iv(18)=maxmin +! controls output + iv(19)=2 +! selects output unit + iv(21)=0 + if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout +! 1 means to print out result + iv(22)=print_min_res +! 1 means to print out summary stats + iv(23)=print_min_stat +! 1 means to print initial x and d + iv(24)=print_min_ini +! min val for v(radfac) default is 0.1 + v(24)=0.1D0 +! max val for v(radfac) default is 4.0 + v(25)=2.0D0 +! v(25)=4.0D0 +! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +! the sumsl default is 0.1 + v(26)=0.1D0 +! false conv if (act fnctn decrease) .lt. v(34) +! the sumsl default is 100*machep + v(34)=v(34)/100.0D0 +! absolute convergence + if (tolf.eq.0.0D0) tolf=1.0D-4 + v(31)=tolf +! relative convergence + if (rtolf.eq.0.0D0) rtolf=1.0D-4 + v(32)=rtolf +! controls initial step size + v(35)=1.0D-1 +! large vals of d correspond to small components of step + do i=1,6*nres + d(i)=1.0D-1 + enddo + + k=0 + do i=1,nres-1 + do j=1,3 + k=k+1 + x(k)=dc(j,i) + enddo + enddo + do i=2,nres-1 + if (ialph(i,1).gt.0) then + do j=1,3 + k=k+1 + x(k)=dc(j,i+nres) + enddo + endif + enddo + call check_ecartint + call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum) + call check_ecartint + k=0 + do i=1,nres-1 + do j=1,3 + k=k+1 + dc(j,i)=x(k) + enddo + enddo + do i=2,nres-1 + if (ialph(i,1).gt.0) then + do j=1,3 + k=k+1 + dc(j,i+nres)=x(k) + enddo + endif + enddo + call chainbuild_cart + +!d call zerograd +!d nf=0 +!d call func_dc(k,x,nf,f,idum,rdum,fdum) +!d call grad_dc(k,x,nf,g,idum,rdum,fdum) +!d +!d do i=1,k +!d x(i)=x(i)+1.0D-5 +!d call func_dc(k,x,nf,f1,idum,rdum,fdum) +!d x(i)=x(i)-1.0D-5 +!d print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5 +!d enddo +!el--------------------- +! write (iout,'(/a)') & +! "Cartesian coordinates of the reference structure after SUMSL" +! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & +! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" +! do i=1,nres +! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') & +! restyp(itype(i)),i,(c(j,i),j=1,3),& +! (c(j,i+nres),j=1,3) +! enddo +!el---------------------------- + etot=v(10) + iretcode=iv(1) + nfun=iv(6) + return + end subroutine minim_dc +!----------------------------------------------------------------------------- + subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm) + + use MPI_data + use energy, only: zerograd,etotal +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif +! include 'COMMON.SETUP' +! include 'COMMON.DERIV' +! include 'COMMON.IOUNITS' +! include 'COMMON.GEO' +! include 'COMMON.CHAIN' +! include 'COMMON.VAR' + integer :: n,nf,k,i,j + real(kind=8) :: energia(0:n_ene) + real(kind=8),external :: ufparm + integer :: uiparm(1) + real(kind=8) :: urparm(1) + real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres) + real(kind=8) :: f + nfl=nf +!bad icg=mod(nf,2)+1 + icg=1 + + k=0 + do i=1,nres-1 + do j=1,3 + k=k+1 + dc(j,i)=x(k) + enddo + enddo + do i=2,nres-1 + if (ialph(i,1).gt.0) then + do j=1,3 + k=k+1 + dc(j,i+nres)=x(k) + enddo + endif + enddo + call chainbuild_cart + + call zerograd + call etotal(energia) + f=energia(0) + +!d print *,'func_dc ',nf,nfl,f + + return + end subroutine func_dc +!----------------------------------------------------------------------------- + subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm) + + use MPI_data + use energy, only: cartgrad,zerograd,etotal +! use MD_data +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif +! include 'COMMON.SETUP' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.VAR' +! include 'COMMON.INTERACT' +! include 'COMMON.FFIELD' +! include 'COMMON.MD' +! include 'COMMON.IOUNITS' + real(kind=8),external :: ufparm + integer :: n,nf,i,j,k + integer :: uiparm(1) + real(kind=8) :: urparm(1) + real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres) + real(kind=8) :: f +! +!elwrite(iout,*) "jestesmy w grad dc" +! +!bad icg=mod(nf,2)+1 + icg=1 +!d print *,'grad_dc ',nf,nfl,nf-nfl+1,icg +!elwrite(iout,*) "jestesmy w grad dc nf-nfl+1", nf-nfl+1 + if (nf-nfl+1) 20,30,40 + 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm) +!d print *,20 + if (nf.eq.0) return + goto 40 + 30 continue +!d print *,30 + k=0 + do i=1,nres-1 + do j=1,3 + k=k+1 + dc(j,i)=x(k) + enddo + enddo + do i=2,nres-1 + if (ialph(i,1).gt.0) then + do j=1,3 + k=k+1 + dc(j,i+nres)=x(k) + enddo + endif + enddo +!elwrite(iout,*) "jestesmy w grad dc" + call chainbuild_cart + +! +! Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +! + 40 call cartgrad +!d print *,40 +!elwrite(iout,*) "jestesmy w grad dc" + k=0 + do i=1,nres-1 + do j=1,3 + k=k+1 + g(k)=gcart(j,i) + enddo + enddo + do i=2,nres-1 + if (ialph(i,1).gt.0) then + do j=1,3 + k=k+1 + g(k)=gxcart(j,i) + enddo + endif + enddo +!elwrite(iout,*) "jestesmy w grad dc" + + return + end subroutine grad_dc +!----------------------------------------------------------------------------- +! minim_mcmf.F +!----------------------------------------------------------------------------- +#ifdef MPI + subroutine minim_mcmf + + use MPI_data + use csa_data + use energy, only: func,gradient,fdum +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' + include 'mpif.h' + integer,parameter :: liv=60 +! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) +! include 'COMMON.VAR' +! include 'COMMON.IOUNITS' +! include 'COMMON.MINIM' +! real(kind=8) :: fdum +! external func,gradient,fdum +!el real(kind=4) :: ran1,ran2,ran3 +! include 'COMMON.SETUP' +! include 'COMMON.GEO' +! include 'COMMON.CHAIN' +! include 'COMMON.FFIELD' + real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) + real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg + real(kind=8),dimension(6*nres) :: d,garbage !(maxvar) (maxvar=6*maxres) +!el real(kind=8) :: v(1:77+(6*nres)*(6*nres+17)/2+1) + integer,dimension(6) :: indx + integer,dimension(liv) :: iv + integer :: lv,idum(1),nf ! + real(kind=8) :: rdum(1) + real(kind=8) :: przes(3),obrot(3,3),eee + logical :: non_conv + + integer,dimension(MPI_STATUS_SIZE) :: muster + + integer :: ichuj,i,ierr + real(kind=8) :: rad,ene0 + data rad /1.745329252d-2/ +!el common /przechowalnia/ v + + lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) + if (.not. allocated(v)) allocate(v(1:lv)) + + ichuj=0 + 10 continue + ichuj = ichuj + 1 + call mpi_recv(indx,6,mpi_integer,king,idint,CG_COMM,& + muster,ierr) + if (indx(1).eq.0) return +! print *, 'worker ',me,' received order ',n,ichuj + call mpi_recv(var,nvar,mpi_double_precision,& + king,idreal,CG_COMM,muster,ierr) + call mpi_recv(ene0,1,mpi_double_precision,& + king,idreal,CG_COMM,muster,ierr) +! print *, 'worker ',me,' var read ' + + + call deflt(2,iv,liv,lv,v) +! 12 means fresh start, dont call deflt + iv(1)=12 +! max num of fun calls + if (maxfun.eq.0) maxfun=500 + iv(17)=maxfun +! max num of iterations + if (maxmin.eq.0) maxmin=1000 + iv(18)=maxmin +! controls output + iv(19)=2 +! selects output unit +! iv(21)=iout + iv(21)=0 +! 1 means to print out result + iv(22)=0 +! 1 means to print out summary stats + iv(23)=0 +! 1 means to print initial x and d + iv(24)=0 +! min val for v(radfac) default is 0.1 + v(24)=0.1D0 +! max val for v(radfac) default is 4.0 + v(25)=2.0D0 +! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +! the sumsl default is 0.1 + v(26)=0.1D0 +! false conv if (act fnctn decrease) .lt. v(34) +! the sumsl default is 100*machep + v(34)=v(34)/100.0D0 +! absolute convergence + if (tolf.eq.0.0D0) tolf=1.0D-4 + v(31)=tolf +! relative convergence + if (rtolf.eq.0.0D0) rtolf=1.0D-4 + v(32)=rtolf +! controls initial step size + v(35)=1.0D-1 +! large vals of d correspond to small components of step + do i=1,nphi + d(i)=1.0D-1 + enddo + do i=nphi+1,nvar + d(i)=1.0D-1 + enddo +! minimize energy + + call func(nvar,var,nf,eee,idum,rdum,fdum) + if(eee.gt.1.0d18) then +! print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL' +! print *,' energy before SUMSL =',eee +! print *,' aborting local minimization' + iv(1)=-1 + v(10)=eee + nf=1 + go to 201 + endif + + call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum) +! find which conformation was returned from sumsl + nf=iv(7)+1 + 201 continue +! total # of ftn evaluations (for iwf=0, it includes all minimizations). + indx(4)=nf + indx(5)=iv(1) + eee=v(10) + + call mpi_send(indx,6,mpi_integer,king,idint,CG_COMM,& + ierr) +! print '(a5,i3,15f10.5)', 'ENEX0',indx(1),v(10) + call mpi_send(var,nvar,mpi_double_precision,& + king,idreal,CG_COMM,ierr) + call mpi_send(eee,1,mpi_double_precision,king,idreal,& + CG_COMM,ierr) + call mpi_send(ene0,1,mpi_double_precision,king,idreal,& + CG_COMM,ierr) + go to 10 + return + end subroutine minim_mcmf +#endif +!----------------------------------------------------------------------------- +! rmdd.f +!----------------------------------------------------------------------------- +! algorithm 611, collected algorithms from acm. +! algorithm appeared in acm-trans. math. software, vol.9, no. 4, +! dec., 1983, p. 503-524. + integer function imdcon(k) +! + integer :: k +! +! *** return integer machine-dependent constants *** +! +! *** k = 1 means return standard output unit number. *** +! *** k = 2 means return alternate output unit number. *** +! *** k = 3 means return input unit number. *** +! (note -- k = 2, 3 are used only by test programs.) +! +! +++ port version follows... +! external i1mach +! integer i1mach +! integer mdperm(3) +! data mdperm(1)/2/, mdperm(2)/4/, mdperm(3)/1/ +! imdcon = i1mach(mdperm(k)) +! +++ end of port version +++ +! +! +++ non-port version follows... + integer :: mdcon(3) + data mdcon(1)/6/, mdcon(2)/8/, mdcon(3)/5/ + imdcon = mdcon(k) +! +++ end of non-port version +++ +! + 999 return +! *** last card of imdcon follows *** + end function imdcon +!----------------------------------------------------------------------------- + real(kind=8) function rmdcon(k) +! +! *** return machine dependent constants used by nl2sol *** +! +! +++ comments below contain data statements for various machines. +++ +! +++ to convert to another machine, place a c in column 1 of the +++ +! +++ data statement line(s) that correspond to the current machine +++ +! +++ and remove the c from column 1 of the data statement line(s) +++ +! +++ that correspond to the new machine. +++ +! + integer :: k +! +! *** the constant returned depends on k... +! +! *** k = 1... smallest pos. eta such that -eta exists. +! *** k = 2... square root of eta. +! *** k = 3... unit roundoff = smallest pos. no. machep such +! *** that 1 + machep .gt. 1 .and. 1 - machep .lt. 1. +! *** k = 4... square root of machep. +! *** k = 5... square root of big (see k = 6). +! *** k = 6... largest machine no. big such that -big exists. +! + real(kind=8) :: big, eta, machep + integer :: bigi(4), etai(4), machei(4) +!/+ +!el real(kind=8) :: dsqrt +!/ + equivalence (big,bigi(1)), (eta,etai(1)), (machep,machei(1)) +! +! +++ ibm 360, ibm 370, or xerox +++ +! +! data big/z7fffffffffffffff/, eta/z0010000000000000/, +! 1 machep/z3410000000000000/ +! +! +++ data general +++ +! +! data big/0.7237005577d+76/, eta/0.5397605347d-78/, +! 1 machep/2.22044605d-16/ +! +! +++ dec 11 +++ +! +! data big/1.7d+38/, eta/2.938735878d-39/, machep/2.775557562d-17/ +! +! +++ hp3000 +++ +! +! data big/1.157920892d+77/, eta/8.636168556d-78/, +! 1 machep/5.551115124d-17/ +! +! +++ honeywell +++ +! +! data big/1.69d+38/, eta/5.9d-39/, machep/2.1680435d-19/ +! +! +++ dec10 +++ +! +! data big/"377777100000000000000000/, +! 1 eta/"002400400000000000000000/, +! 2 machep/"104400000000000000000000/ +! +! +++ burroughs +++ +! +! data big/o0777777777777777,o7777777777777777/, +! 1 eta/o1771000000000000,o7770000000000000/, +! 2 machep/o1451000000000000,o0000000000000000/ +! +! +++ control data +++ +! +! data big/37767777777777777777b,37167777777777777777b/, +! 1 eta/00014000000000000000b,00000000000000000000b/, +! 2 machep/15614000000000000000b,15010000000000000000b/ +! +! +++ prime +++ +! +! data big/1.0d+9786/, eta/1.0d-9860/, machep/1.4210855d-14/ +! +! +++ univac +++ +! +! data big/8.988d+307/, eta/1.2d-308/, machep/1.734723476d-18/ +! +! +++ vax +++ +! + data big/1.7d+38/, eta/2.939d-39/, machep/1.3877788d-17/ +! +! +++ cray 1 +++ +! +! data bigi(1)/577767777777777777777b/, +! 1 bigi(2)/000007777777777777776b/, +! 2 etai(1)/200004000000000000000b/, +! 3 etai(2)/000000000000000000000b/, +! 4 machei(1)/377224000000000000000b/, +! 5 machei(2)/000000000000000000000b/ +! +! +++ port library -- requires more than just a data statement... +++ +! +! external d1mach +! double precision d1mach, zero +! data big/0.d+0/, eta/0.d+0/, machep/0.d+0/, zero/0.d+0/ +! if (big .gt. zero) go to 1 +! big = d1mach(2) +! eta = d1mach(1) +! machep = d1mach(4) +!1 continue +! +! +++ end of port +++ +! +!------------------------------- body -------------------------------- +! + go to (10, 20, 30, 40, 50, 60), k +! + 10 rmdcon = eta + go to 999 +! + 20 rmdcon = dsqrt(256.d+0*eta)/16.d+0 + go to 999 +! + 30 rmdcon = machep + go to 999 +! + 40 rmdcon = dsqrt(machep) + go to 999 +! + 50 rmdcon = dsqrt(big/256.d+0)*16.d+0 + go to 999 +! + 60 rmdcon = big +! + 999 return +! *** last card of rmdcon follows *** + end function rmdcon +!----------------------------------------------------------------------------- +! sc_move.F +!----------------------------------------------------------------------------- + subroutine sc_move(n_start,n_end,n_maxtry,e_drop,n_fun,etot) + + use control + use random, only: iran_num + use energy, only: esc +! Perform a quick search over side-chain arrangments (over +! residues n_start to n_end) for a given (frozen) CA trace +! Only side-chains are minimized (at most n_maxtry times each), +! not CA positions +! Stops if energy drops by e_drop, otherwise tries all residues +! in the given range +! If there is an energy drop, full minimization may be useful +! n_start, n_end CAN be modified by this routine, but only if +! out of bounds (n_start <= 1, n_end >= nres, n_start < n_end) +! NOTE: this move should never increase the energy +!rc implicit none + +! Includes +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' + include 'mpif.h' +! include 'COMMON.GEO' +! include 'COMMON.VAR' +! include 'COMMON.HEADER' +! include 'COMMON.IOUNITS' +! include 'COMMON.CHAIN' +! include 'COMMON.FFIELD' + +! External functions +!el integer iran_num +!el external iran_num + +! Input arguments + integer :: n_start,n_end,n_maxtry + real(kind=8) :: e_drop + +! Output arguments + integer :: n_fun + real(kind=8) :: etot + +! Local variables +! real(kind=8) :: energy(0:n_ene) + real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1) + real(kind=8) :: orig_e,cur_e + integer :: n,n_steps,n_first,n_cur,n_tot !,i + real(kind=8) :: orig_w(0:n_ene) + real(kind=8) :: wtime + +!elwrite(iout,*) "in sc_move etot= ", etot +! Set non side-chain weights to zero (minimization is faster) +! NOTE: e(2) does not actually depend on the side-chain, only CA + orig_w(2)=wscp + orig_w(3)=welec + orig_w(4)=wcorr + orig_w(5)=wcorr5 + orig_w(6)=wcorr6 + orig_w(7)=wel_loc + orig_w(8)=wturn3 + orig_w(9)=wturn4 + orig_w(10)=wturn6 + orig_w(11)=wang + orig_w(13)=wtor + orig_w(14)=wtor_d + orig_w(15)=wvdwpp + + wscp=0.D0 + welec=0.D0 + wcorr=0.D0 + wcorr5=0.D0 + wcorr6=0.D0 + wel_loc=0.D0 + wturn3=0.D0 + wturn4=0.D0 + wturn6=0.D0 + wang=0.D0 + wtor=0.D0 + wtor_d=0.D0 + wvdwpp=0.D0 + +! Make sure n_start, n_end are within proper range + if (n_start.lt.2) n_start=2 + if (n_end.gt.nres-1) n_end=nres-1 +!rc if (n_start.lt.n_end) then + if (n_start.gt.n_end) then + n_start=2 + n_end=nres-1 + endif + +! Save the initial values of energy and coordinates +!d call chainbuild +!d call etotal(energy) +!d write (iout,*) 'start sc ene',energy(0) +!d call enerprint(energy(0)) +!rc etot=energy(0) + n_fun=0 +!rc orig_e=etot +!rc cur_e=orig_e +!rc do i=2,nres-1 +!rc cur_alph(i)=alph(i) +!rc cur_omeg(i)=omeg(i) +!rc enddo + +!t wtime=MPI_WTIME() +! Try (one by one) all specified residues, starting from a +! random position in sequence +! Stop early if the energy has decreased by at least e_drop + n_tot=n_end-n_start+1 + n_first=iran_num(0,n_tot-1) + n_steps=0 + n=0 +!rc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop) + do while (n.lt.n_tot) + n_cur=n_start+mod(n_first+n,n_tot) + call single_sc_move(n_cur,n_maxtry,e_drop,& + n_steps,n_fun,etot) +!elwrite(iout,*) "after msingle sc_move etot= ", etot +! If a lower energy was found, update the current structure... +!rc if (etot.lt.cur_e) then +!rc cur_e=etot +!rc do i=2,nres-1 +!rc cur_alph(i)=alph(i) +!rc cur_omeg(i)=omeg(i) +!rc enddo +!rc else +! ...else revert to the previous one +!rc etot=cur_e +!rc do i=2,nres-1 +!rc alph(i)=cur_alph(i) +!rc omeg(i)=cur_omeg(i) +!rc enddo +!rc endif + n=n+1 +!d +!d call chainbuild +!d call etotal(energy) +!d print *,'running',n,energy(0) + enddo + +!d call chainbuild +!d call etotal(energy) +!d write (iout,*) 'end sc ene',energy(0) + +! Put the original weights back to calculate the full energy + wscp=orig_w(2) + welec=orig_w(3) + wcorr=orig_w(4) + wcorr5=orig_w(5) + wcorr6=orig_w(6) + wel_loc=orig_w(7) + wturn3=orig_w(8) + wturn4=orig_w(9) + wturn6=orig_w(10) + wang=orig_w(11) + wtor=orig_w(13) + wtor_d=orig_w(14) + wvdwpp=orig_w(15) + +!rc n_fun=n_fun+1 +!t write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime + return + end subroutine sc_move +!----------------------------------------------------------------------------- + subroutine single_sc_move(res_pick,n_maxtry,e_drop,n_steps,n_fun,e_sc) + +! Perturb one side-chain (res_pick) and minimize the +! neighbouring region, keeping all CA's and non-neighbouring +! side-chains fixed +! Try until e_drop energy improvement is achieved, or n_maxtry +! attempts have been made +! At the start, e_sc should contain the side-chain-only energy(0) +! nsteps and nfun for this move are ADDED to n_steps and n_fun +!rc implicit none + use energy, only: esc + use geometry, only:dist +! Includes +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.VAR' +! include 'COMMON.INTERACT' +! include 'COMMON.CHAIN' +! include 'COMMON.MINIM' +! include 'COMMON.FFIELD' +! include 'COMMON.IOUNITS' + +! External functions +!el double precision dist +!el external dist + +! Input arguments + integer :: res_pick,n_maxtry + real(kind=8) :: e_drop + +! Input/Output arguments + integer :: n_steps,n_fun + real(kind=8) :: e_sc + +! Local variables + logical :: fail + integer :: i,j + integer :: nres_moved + integer :: iretcode,loc_nfun,orig_maxfun,n_try + real(kind=8) :: sc_dist,sc_dist_cutoff +! real(kind=8) :: energy_(0:n_ene) + real(kind=8) :: evdw,escloc,orig_e,cur_e + real(kind=8) :: cur_alph(2:nres-1),cur_omeg(2:nres-1) + real(kind=8) :: var(6*nres) !(maxvar) (maxvar=6*maxres) + + real(kind=8) :: orig_theta(1:nres),orig_phi(1:nres),& + orig_alph(1:nres),orig_omeg(1:nres) + +!elwrite(iout,*) "in sinle etot/ e_sc",e_sc +! Define what is meant by "neighbouring side-chain" + sc_dist_cutoff=5.0D0 + +! Don't do glycine or ends + i=itype(res_pick) + if (i.eq.10 .or. i.eq.ntyp1) return + +! Freeze everything (later will relax only selected side-chains) + mask_r=.true. + do i=1,nres + mask_phi(i)=0 + mask_theta(i)=0 + mask_side(i)=0 + enddo + +! Find the neighbours of the side-chain to move +! and save initial variables +!rc orig_e=e_sc +!rc cur_e=orig_e + nres_moved=0 + do i=2,nres-1 +! Don't do glycine (itype(j)==10) + if (itype(i).ne.10) then + sc_dist=dist(nres+i,nres+res_pick) + else + sc_dist=sc_dist_cutoff + endif + if (sc_dist.lt.sc_dist_cutoff) then + nres_moved=nres_moved+1 + mask_side(i)=1 + cur_alph(i)=alph(i) + cur_omeg(i)=omeg(i) + endif + enddo + + call chainbuild + call egb1(evdw) + call esc(escloc) +!elwrite(iout,*) "in sinle etot/ e_sc",e_sc +!elwrite(iout,*) "in sinle wsc=",wsc,"evdw",evdw,"wscloc",wscloc,"escloc",escloc + e_sc=wsc*evdw+wscloc*escloc +!elwrite(iout,*) "in sinle etot/ e_sc",e_sc +!d call etotal(energy) +!d print *,'new ',(energy(k),k=0,n_ene) + orig_e=e_sc + cur_e=orig_e + + n_try=0 + do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop) +! Move the selected residue (don't worry if it fails) + call gen_side(iabs(itype(res_pick)),theta(res_pick+1),& + alph(res_pick),omeg(res_pick),fail) + +! Minimize the side-chains starting from the new arrangement + call geom_to_var(nvar,var) + orig_maxfun=maxfun + maxfun=7 + +!rc do i=1,nres +!rc orig_theta(i)=theta(i) +!rc orig_phi(i)=phi(i) +!rc orig_alph(i)=alph(i) +!rc orig_omeg(i)=omeg(i) +!rc enddo + +!elwrite(iout,*) "in sinle etot/ e_sc",e_sc + call minimize_sc1(e_sc,var,iretcode,loc_nfun) + +!elwrite(iout,*) "in sinle etot/ e_sc",e_sc +!v write(*,'(2i3,2f12.5,2i3)') +!v & res_pick,nres_moved,orig_e,e_sc-cur_e, +!v & iretcode,loc_nfun + +!$$$ if (iretcode.eq.8) then +!$$$ write(iout,*)'Coordinates just after code 8' +!$$$ call chainbuild +!$$$ call all_varout +!$$$ call flush(iout) +!$$$ do i=1,nres +!$$$ theta(i)=orig_theta(i) +!$$$ phi(i)=orig_phi(i) +!$$$ alph(i)=orig_alph(i) +!$$$ omeg(i)=orig_omeg(i) +!$$$ enddo +!$$$ write(iout,*)'Coordinates just before code 8' +!$$$ call chainbuild +!$$$ call all_varout +!$$$ call flush(iout) +!$$$ endif + + n_fun=n_fun+loc_nfun + maxfun=orig_maxfun + call var_to_geom(nvar,var) + +! If a lower energy was found, update the current structure... + if (e_sc.lt.cur_e) then +!v call chainbuild +!v call etotal(energy) +!d call egb1(evdw) +!d call esc(escloc) +!d e_sc1=wsc*evdw+wscloc*escloc +!d print *,' new',e_sc1,energy(0) +!v print *,'new ',energy(0) +!d call enerprint(energy(0)) + cur_e=e_sc + do i=2,nres-1 + if (mask_side(i).eq.1) then + cur_alph(i)=alph(i) + cur_omeg(i)=omeg(i) + endif + enddo + else +! ...else revert to the previous one + e_sc=cur_e + do i=2,nres-1 + if (mask_side(i).eq.1) then + alph(i)=cur_alph(i) + omeg(i)=cur_omeg(i) + endif + enddo + endif + n_try=n_try+1 + + enddo + n_steps=n_steps+n_try + +! Reset the minimization mask_r to false + mask_r=.false. + + return + end subroutine single_sc_move +!----------------------------------------------------------------------------- + subroutine sc_minimize(etot,iretcode,nfun) + +! Minimizes side-chains only, leaving backbone frozen +!rc implicit none + use energy, only: etotal +! Includes +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.VAR' +! include 'COMMON.CHAIN' +! include 'COMMON.FFIELD' + +! Output arguments + real(kind=8) :: etot + integer :: iretcode,nfun + +! Local variables + integer :: i + real(kind=8) :: orig_w(0:n_ene),energy_(0:n_ene) + real(kind=8) :: var(6*nres) !(maxvar)(maxvar=6*maxres) + + +! Set non side-chain weights to zero (minimization is faster) +! NOTE: e(2) does not actually depend on the side-chain, only CA + orig_w(2)=wscp + orig_w(3)=welec + orig_w(4)=wcorr + orig_w(5)=wcorr5 + orig_w(6)=wcorr6 + orig_w(7)=wel_loc + orig_w(8)=wturn3 + orig_w(9)=wturn4 + orig_w(10)=wturn6 + orig_w(11)=wang + orig_w(13)=wtor + orig_w(14)=wtor_d + + wscp=0.D0 + welec=0.D0 + wcorr=0.D0 + wcorr5=0.D0 + wcorr6=0.D0 + wel_loc=0.D0 + wturn3=0.D0 + wturn4=0.D0 + wturn6=0.D0 + wang=0.D0 + wtor=0.D0 + wtor_d=0.D0 + +! Prepare to freeze backbone + do i=1,nres + mask_phi(i)=0 + mask_theta(i)=0 + mask_side(i)=1 + enddo + +! Minimize the side-chains + mask_r=.true. + call geom_to_var(nvar,var) + call minimize(etot,var,iretcode,nfun) + call var_to_geom(nvar,var) + mask_r=.false. + +! Put the original weights back and calculate the full energy + wscp=orig_w(2) + welec=orig_w(3) + wcorr=orig_w(4) + wcorr5=orig_w(5) + wcorr6=orig_w(6) + wel_loc=orig_w(7) + wturn3=orig_w(8) + wturn4=orig_w(9) + wturn6=orig_w(10) + wang=orig_w(11) + wtor=orig_w(13) + wtor_d=orig_w(14) + + call chainbuild + call etotal(energy_) + etot=energy_(0) + + return + end subroutine sc_minimize +!----------------------------------------------------------------------------- + subroutine minimize_sc1(etot,x,iretcode,nfun) + + use energy, only: func,gradient,fdum,etotal,enerprint + use comm_srutu +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' + integer,parameter :: liv=60 + integer :: iretcode,nfun +! integer :: lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) +! include 'COMMON.IOUNITS' +! include 'COMMON.VAR' +! include 'COMMON.GEO' +! include 'COMMON.MINIM' +!el integer :: icall +!el common /srutu/ icall + integer,dimension(liv) :: iv + real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv) + real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) + real(kind=8) :: energia(0:n_ene) +!el real(kind=8) :: fdum +! external gradient,fdum !func, +! external func_restr1,grad_restr1 + logical :: not_done,change,reduce +!el common /przechowalnia/ v + + integer :: nvar_restr,lv,i,j + integer :: idum(1) + real(kind=8) :: rdum(1),etot !,fdum + + lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) + if (.not. allocated(v)) allocate(v(1:lv)) + + call deflt(2,iv,liv,lv,v) +! 12 means fresh start, dont call deflt + iv(1)=12 +! max num of fun calls + if (maxfun.eq.0) maxfun=500 + iv(17)=maxfun +! max num of iterations + if (maxmin.eq.0) maxmin=1000 + iv(18)=maxmin +! controls output + iv(19)=2 +! selects output unit +! iv(21)=iout + iv(21)=0 +! 1 means to print out result + iv(22)=0 +! 1 means to print out summary stats + iv(23)=0 +! 1 means to print initial x and d + iv(24)=0 +! min val for v(radfac) default is 0.1 + v(24)=0.1D0 +! max val for v(radfac) default is 4.0 + v(25)=2.0D0 +! v(25)=4.0D0 +! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +! the sumsl default is 0.1 + v(26)=0.1D0 +! false conv if (act fnctn decrease) .lt. v(34) +! the sumsl default is 100*machep + v(34)=v(34)/100.0D0 +! absolute convergence + if (tolf.eq.0.0D0) tolf=1.0D-4 + v(31)=tolf +! relative convergence + if (rtolf.eq.0.0D0) rtolf=1.0D-4 + v(32)=rtolf +! controls initial step size + v(35)=1.0D-1 +! large vals of d correspond to small components of step + do i=1,nphi + d(i)=1.0D-1 + enddo + do i=nphi+1,nvar + d(i)=1.0D-1 + enddo +!elmask_r=.false. + IF (mask_r) THEN + call x2xx(x,xx,nvar_restr) + call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,& + iv,liv,lv,v,idum,rdum,fdum) + call xx2x(x,xx) + ELSE + call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum) + ENDIF +!el--------------------- +! write (iout,'(/a)') & +! "Cartesian coordinates of the reference structure after SUMSL" +! write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & +! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" +! do i=1,nres +! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') & +! restyp(itype(i)),i,(c(j,i),j=1,3),& +! (c(j,i+nres),j=1,3) +! enddo +! call etotal(energia) +! call enerprint(energia) +!el---------------------------- + etot=v(10) + iretcode=iv(1) + nfun=iv(6) + + return + end subroutine minimize_sc1 +!----------------------------------------------------------------------------- + subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm) + + use comm_chu + use energy, only: zerograd,esc,sc_grad +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.DERIV' +! include 'COMMON.IOUNITS' +! include 'COMMON.GEO' +! include 'COMMON.FFIELD' +! include 'COMMON.INTERACT' +! include 'COMMON.TIME1' + integer :: n,nf,i,j +!el common /chuju/ jjj + real(kind=8) :: energia(0:n_ene),evdw,escloc + real(kind=8) :: e1,e2,f + real(kind=8),external :: ufparm + integer :: uiparm(1) + real(kind=8) :: urparm(1) + real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres) + nfl=nf + icg=mod(nf,2)+1 + +#ifdef OSF +! Intercept NaNs in the coordinates, before calling etotal + x_sum=0.D0 + do i=1,n + x_sum=x_sum+x(i) + enddo + FOUND_NAN=.false. + if (x_sum.ne.x_sum) then + write(iout,*)" *** func_restr1 : Found NaN in coordinates" + f=1.0D+73 + FOUND_NAN=.true. + return + endif +#endif + + call var_to_geom_restr(n,x) + call zerograd + call chainbuild +!d write (iout,*) 'ETOTAL called from FUNC' + call egb1(evdw) + call esc(escloc) + f=wsc*evdw+wscloc*escloc +!d call etotal(energia(0)) +!d f=wsc*energia(1)+wscloc*energia(12) +!d print *,f,evdw,escloc,energia(0) +! +! Sum up the components of the Cartesian gradient. +! + do i=1,nct + do j=1,3 + gradx(j,i,icg)=wsc*gvdwx(j,i) + enddo + enddo + + return + end subroutine func_restr1 +!----------------------------------------------------------------------------- + subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm) + + use energy, only: cartder,zerograd,esc,sc_grad +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.VAR' +! include 'COMMON.INTERACT' +! include 'COMMON.FFIELD' +! include 'COMMON.IOUNITS' +!el external ufparm + integer :: i,j,k,ind,n,nf,uiparm(1) + real(kind=8) :: f,urparm(1) + real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres) + integer :: ig,igall,ij + real(kind=8) :: gphii,gthetai,galphai,gomegai + real(kind=8),external :: ufparm + + icg=mod(nf,2)+1 + if (nf-nfl+1) 20,30,40 + 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm) +! write (iout,*) 'grad 20' + if (nf.eq.0) return + goto 40 + 30 call var_to_geom_restr(n,x) + call chainbuild +! +! Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +! + 40 call cartder +! +! Convert the Cartesian gradient into internal-coordinate gradient. +! + + ig=0 + ind=nres-2 + do i=2,nres-2 + IF (mask_phi(i+2).eq.1) THEN + gphii=0.0D0 + do j=i+1,nres-1 + ind=ind+1 + do k=1,3 + gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) + gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) + enddo + enddo + ig=ig+1 + g(ig)=gphii + ELSE + ind=ind+nres-1-i + ENDIF + enddo + + + ind=0 + do i=1,nres-2 + IF (mask_theta(i+2).eq.1) THEN + ig=ig+1 + gthetai=0.0D0 + do j=i+1,nres-1 + ind=ind+1 + do k=1,3 + gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) + gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) + enddo + enddo + g(ig)=gthetai + ELSE + ind=ind+nres-1-i + ENDIF + enddo + + do i=2,nres-1 + if (itype(i).ne.10) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + galphai=0.0D0 + do k=1,3 + galphai=galphai+dxds(k,i)*gradx(k,i,icg) + enddo + g(ig)=galphai + ENDIF + endif + enddo + + + do i=2,nres-1 + if (itype(i).ne.10) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + gomegai=0.0D0 + do k=1,3 + gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) + enddo + g(ig)=gomegai + ENDIF + endif + enddo + +! +! Add the components corresponding to local energy terms. +! + + ig=0 + igall=0 + do i=4,nres + igall=igall+1 + if (mask_phi(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do i=3,nres + igall=igall+1 + if (mask_theta(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do ij=1,2 + do i=2,nres-1 + if (itype(i).ne.10) then + igall=igall+1 + if (mask_side(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + endif + enddo + enddo + +!d do i=1,ig +!d write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) +!d enddo + return + end subroutine grad_restr1 +!----------------------------------------------------------------------------- + subroutine egb1(evdw) +! +! This subroutine calculates the interaction energy of nonbonded side chains +! assuming the Gay-Berne potential of interaction. +! + use calc_data + use energy, only: sc_grad +! use control, only:stopx +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.GEO' +! include 'COMMON.VAR' +! include 'COMMON.LOCAL' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.NAMES' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.CALC' +! include 'COMMON.CONTROL' + logical :: lprn + real(kind=8) :: evdw +!el local variables + integer :: iint,ind,itypi,itypi1,itypj + real(kind=8) :: xi,yi,zi,rrij,sig,sig0ij,rij_shift,fac,e1,e2,& + sigm,epsi +!elwrite(iout,*) "check evdw" +! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon + evdw=0.0D0 + lprn=.false. +! if (icall.eq.0) lprn=.true. + ind=0 + do i=iatsc_s,iatsc_e + + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + dsci_inv=dsc_inv(itypi) +!elwrite(iout,*) itypi,itypi1,xi,yi,zi,dxi,dyi,dzi,dsci_inv +! +! Calculate SC interaction energy. +! + do iint=1,nint_gr(i) + do j=istart(i,iint),iend(i,iint) + IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN + ind=ind+1 + itypj=iabs(itype(j)) + dscj_inv=dsc_inv(itypj) + sig0ij=sigma(itypi,itypj) + chi1=chi(itypi,itypj) + chi2=chi(itypj,itypi) + chi12=chi1*chi2 + chip1=chip(itypi) + chip2=chip(itypj) + chip12=chip1*chip2 + alf1=alp(itypi) + alf2=alp(itypj) + alf12=0.5D0*(alf1+alf2) +! For diagnostics only!!! +! chi1=0.0D0 +! chi2=0.0D0 +! chi12=0.0D0 +! chip1=0.0D0 +! chip2=0.0D0 +! chip12=0.0D0 +! alf1=0.0D0 +! alf2=0.0D0 +! alf12=0.0D0 + xj=c(1,nres+j)-xi + yj=c(2,nres+j)-yi + zj=c(3,nres+j)-zi + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + rij=dsqrt(rrij) +! Calculate angle-dependent terms of energy and contributions to their +! derivatives. + call sc_angular + sigsq=1.0D0/sigsq + sig=sig0ij*dsqrt(sigsq) + rij_shift=1.0D0/rij-sig+sig0ij +! I hate to put IF's in the loops, but here don't have another choice!!!! + if (rij_shift.le.0.0D0) then + evdw=1.0D20 +!d write (iout,'(2(a3,i3,2x),17(0pf7.3))') & +!d restyp(itypi),i,restyp(itypj),j, & +!d rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) + return + endif + sigder=-sig*sigsq +!--------------------------------------------------------------- + rij_shift=1.0D0/rij_shift + fac=rij_shift**expon + e1=fac*fac*aa(itypi,itypj) + e2=fac*bb(itypi,itypj) + evdwij=eps1*eps2rt*eps3rt*(e1+e2) + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt + evdwij=evdwij*eps2rt*eps3rt + evdw=evdw+evdwij + if (lprn) then + sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +!d write (iout,'(2(a3,i3,2x),17(0pf7.3))') & +!d restyp(itypi),i,restyp(itypj),j, & +!d epsi,sigm,chi1,chi2,chip1,chip2, & +!d eps1,eps2rt**2,eps3rt**2,sig,sig0ij, & +!d om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, & +!d evdwij + endif + + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & + 'evdw',i,j,evdwij + +! Calculate gradient components. + e1=e1*eps1*eps2rt**2*eps3rt**2 + fac=-expon*(e1+evdwij)*rij_shift + sigder=fac*sigder + fac=rij*fac +! Calculate the radial part of the gradient + gg(1)=xj*fac + gg(2)=yj*fac + gg(3)=zj*fac +! Calculate angular part of the gradient. + +!elwrite(iout,*) evdw + call sc_grad +!elwrite(iout,*) "evdw=",evdw,j,iint,i + ENDIF +!elwrite(iout,*) evdw + enddo ! j +!elwrite(iout,*) evdw + enddo ! iint +!elwrite(iout,*) evdw + enddo ! i +!elwrite(iout,*) evdw,i + end subroutine egb1 +!----------------------------------------------------------------------------- +! sumsld.f +!----------------------------------------------------------------------------- + subroutine sumsl(n,d,x,calcf,calcg,iv,liv,lv,v,uiparm,urparm,ufparm) +! +! *** minimize general unconstrained objective function using *** +! *** analytic gradient and hessian approx. from secant update *** +! +! use control + integer :: n, liv, lv + integer :: iv(liv), uiparm(1) + real(kind=8) :: d(n), x(n), v(lv), urparm(1) + real(kind=8),external :: ufparm !funtion name as an argument + +! dimension v(71 + n*(n+15)/2), uiparm(*), urparm(*) + external :: calcf, calcg !subroutine name as an argument +! +! *** purpose *** +! +! this routine interacts with subroutine sumit in an attempt +! to find an n-vector x* that minimizes the (unconstrained) +! objective function computed by calcf. (often the x* found is +! a local minimizer rather than a global one.) +! +!-------------------------- parameter usage -------------------------- +! +! n........ (input) the number of variables on which f depends, i.e., +! the number of components in x. +! d........ (input/output) a scale vector such that d(i)*x(i), +! i = 1,2,...,n, are all in comparable units. +! d can strongly affect the behavior of sumsl. +! finding the best choice of d is generally a trial- +! and-error process. choosing d so that d(i)*x(i) +! has about the same value for all i often works well. +! the defaults provided by subroutine deflt (see i +! below) require the caller to supply d. +! x........ (input/output) before (initially) calling sumsl, the call- +! er should set x to an initial guess at x*. when +! sumsl returns, x contains the best point so far +! found, i.e., the one that gives the least value so +! far seen for f(x). +! calcf.... (input) a subroutine that, given x, computes f(x). calcf +! must be declared external in the calling program. +! it is invoked by +! call calcf(n, x, nf, f, uiparm, urparm, ufparm) +! when calcf is called, nf is the invocation +! count for calcf. nf is included for possible use +! with calcg. if x is out of bounds (e.g., if it +! would cause overflow in computing f(x)), then calcf +! should set nf to 0. this will cause a shorter step +! to be attempted. (if x is in bounds, then calcf +! should not change nf.) the other parameters are as +! described above and below. calcf should not change +! n, p, or x. +! calcg.... (input) a subroutine that, given x, computes g(x), the gra- +! dient of f at x. calcg must be declared external in +! the calling program. it is invoked by +! call calcg(n, x, nf, g, uiparm, urparm, ufaprm) +! when calcg is called, nf is the invocation +! count for calcf at the time f(x) was evaluated. the +! x passed to calcg is usually the one passed to calcf +! on either its most recent invocation or the one +! prior to it. if calcf saves intermediate results +! for use by calcg, then it is possible to tell from +! nf whether they are valid for the current x (or +! which copy is valid if two copies are kept). if g +! cannot be computed at x, then calcg should set nf to +! 0. in this case, sumsl will return with iv(1) = 65. +! (if g can be computed at x, then calcg should not +! changed nf.) the other parameters to calcg are as +! described above and below. calcg should not change +! n or x. +! iv....... (input/output) an integer value array of length liv (see +! below) that helps control the sumsl algorithm and +! that is used to store various intermediate quanti- +! ties. of particular interest are the initialization/ +! return code iv(1) and the entries in iv that control +! printing and limit the number of iterations and func- +! tion evaluations. see the section on iv input +! values below. +! liv...... (input) length of iv array. must be at least 60. if li +! is too small, then sumsl returns with iv(1) = 15. +! when sumsl returns, the smallest allowed value of +! liv is stored in iv(lastiv) -- see the section on +! iv output values below. (this is intended for use +! with extensions of sumsl that handle constraints.) +! lv....... (input) length of v array. must be at least 71+n*(n+15)/2. +! (at least 77+n*(n+17)/2 for smsno, at least +! 78+n*(n+12) for humsl). if lv is too small, then +! sumsl returns with iv(1) = 16. when sumsl returns, +! the smallest allowed value of lv is stored in +! iv(lastv) -- see the section on iv output values +! below. +! v........ (input/output) a floating-point value array of length l +! (see below) that helps control the sumsl algorithm +! and that is used to store various intermediate +! quantities. of particular interest are the entries +! in v that limit the length of the first step +! attempted (lmax0) and specify convergence tolerances +! (afctol, lmaxs, rfctol, sctol, xctol, xftol). +! uiparm... (input) user integer parameter array passed without change +! to calcf and calcg. +! urparm... (input) user floating-point parameter array passed without +! change to calcf and calcg. +! ufparm... (input) user external subroutine or function passed without +! change to calcf and calcg. +! +! *** iv input values (from subroutine deflt) *** +! +! iv(1)... on input, iv(1) should have a value between 0 and 14...... +! 0 and 12 mean this is a fresh start. 0 means that +! deflt(2, iv, liv, lv, v) +! is to be called to provide all default values to iv and +! v. 12 (the value that deflt assigns to iv(1)) means the +! caller has already called deflt and has possibly changed +! some iv and/or v entries to non-default values. +! 13 means deflt has been called and that sumsl (and +! sumit) should only do their storage allocation. that is, +! they should set the output components of iv that tell +! where various subarrays arrays of v begin, such as iv(g) +! (and, for humsl and humit only, iv(dtol)), and return. +! 14 means that a storage has been allocated (by a call +! with iv(1) = 13) and that the algorithm should be +! started. when called with iv(1) = 13, sumsl returns +! iv(1) = 14 unless liv or lv is too small (or n is not +! positive). default = 12. +! iv(inith).... iv(25) tells whether the hessian approximation h should +! be initialized. 1 (the default) means sumit should +! initialize h to the diagonal matrix whose i-th diagonal +! element is d(i)**2. 0 means the caller has supplied a +! cholesky factor l of the initial hessian approximation +! h = l*(l**t) in v, starting at v(iv(lmat)) = v(iv(42)) +! (and stored compactly by rows). note that iv(lmat) may +! be initialized by calling sumsl with iv(1) = 13 (see +! the iv(1) discussion above). default = 1. +! iv(mxfcal)... iv(17) gives the maximum number of function evaluations +! (calls on calcf) allowed. if this number does not suf- +! fice, then sumsl returns with iv(1) = 9. default = 200. +! iv(mxiter)... iv(18) gives the maximum number of iterations allowed. +! it also indirectly limits the number of gradient evalua- +! tions (calls on calcg) to iv(mxiter) + 1. if iv(mxiter) +! iterations do not suffice, then sumsl returns with +! iv(1) = 10. default = 150. +! iv(outlev)... iv(19) controls the number and length of iteration sum- +! mary lines printed (by itsum). iv(outlev) = 0 means do +! not print any summary lines. otherwise, print a summary +! line after each abs(iv(outlev)) iterations. if iv(outlev) +! is positive, then summary lines of length 78 (plus carri- +! age control) are printed, including the following... the +! iteration and function evaluation counts, f = the current +! function value, relative difference in function values +! achieved by the latest step (i.e., reldf = (f0-v(f))/f01, +! where f01 is the maximum of abs(v(f)) and abs(v(f0)) and +! v(f0) is the function value from the previous itera- +! tion), the relative function reduction predicted for the +! step just taken (i.e., preldf = v(preduc) / f01, where +! v(preduc) is described below), the scaled relative change +! in x (see v(reldx) below), the step parameter for the +! step just taken (stppar = 0 means a full newton step, +! between 0 and 1 means a relaxed newton step, between 1 +! and 2 means a double dogleg step, greater than 2 means +! a scaled down cauchy step -- see subroutine dbldog), the +! 2-norm of the scale vector d times the step just taken +! (see v(dstnrm) below), and npreldf, i.e., +! v(nreduc)/f01, where v(nreduc) is described below -- if +! npreldf is positive, then it is the relative function +! reduction predicted for a newton step (one with +! stppar = 0). if npreldf is negative, then it is the +! negative of the relative function reduction predicted +! for a step computed with step bound v(lmaxs) for use in +! testing for singular convergence. +! if iv(outlev) is negative, then lines of length 50 +! are printed, including only the first 6 items listed +! above (through reldx). +! default = 1. +! iv(parprt)... iv(20) = 1 means print any nondefault v values on a +! fresh start or any changed v values on a restart. +! iv(parprt) = 0 means skip this printing. default = 1. +! iv(prunit)... iv(21) is the output unit number on which all printing +! is done. iv(prunit) = 0 means suppress all printing. +! default = standard output unit (unit 6 on most systems). +! iv(solprt)... iv(22) = 1 means print out the value of x returned (as +! well as the gradient and the scale vector d). +! iv(solprt) = 0 means skip this printing. default = 1. +! iv(statpr)... iv(23) = 1 means print summary statistics upon return- +! ing. these consist of the function value, the scaled +! relative change in x caused by the most recent step (see +! v(reldx) below), the number of function and gradient +! evaluations (calls on calcf and calcg), and the relative +! function reductions predicted for the last step taken and +! for a newton step (or perhaps a step bounded by v(lmaxs) +! -- see the descriptions of preldf and npreldf under +! iv(outlev) above). +! iv(statpr) = 0 means skip this printing. +! iv(statpr) = -1 means skip this printing as well as that +! of the one-line termination reason message. default = 1. +! iv(x0prt).... iv(24) = 1 means print the initial x and scale vector d +! (on a fresh start only). iv(x0prt) = 0 means skip this +! printing. default = 1. +! +! *** (selected) iv output values *** +! +! iv(1)........ on output, iv(1) is a return code.... +! 3 = x-convergence. the scaled relative difference (see +! v(reldx)) between the current parameter vector x and +! a locally optimal parameter vector is very likely at +! most v(xctol). +! 4 = relative function convergence. the relative differ- +! ence between the current function value and its lo- +! cally optimal value is very likely at most v(rfctol). +! 5 = both x- and relative function convergence (i.e., the +! conditions for iv(1) = 3 and iv(1) = 4 both hold). +! 6 = absolute function convergence. the current function +! value is at most v(afctol) in absolute value. +! 7 = singular convergence. the hessian near the current +! iterate appears to be singular or nearly so, and a +! step of length at most v(lmaxs) is unlikely to yield +! a relative function decrease of more than v(sctol). +! 8 = false convergence. the iterates appear to be converg- +! ing to a noncritical point. this may mean that the +! convergence tolerances (v(afctol), v(rfctol), +! v(xctol)) are too small for the accuracy to which +! the function and gradient are being computed, that +! there is an error in computing the gradient, or that +! the function or gradient is discontinuous near x. +! 9 = function evaluation limit reached without other con- +! vergence (see iv(mxfcal)). +! 10 = iteration limit reached without other convergence +! (see iv(mxiter)). +! 11 = stopx returned .true. (external interrupt). see the +! usage notes below. +! 14 = storage has been allocated (after a call with +! iv(1) = 13). +! 17 = restart attempted with n changed. +! 18 = d has a negative component and iv(dtype) .le. 0. +! 19...43 = v(iv(1)) is out of range. +! 63 = f(x) cannot be computed at the initial x. +! 64 = bad parameters passed to assess (which should not +! occur). +! 65 = the gradient could not be computed at x (see calcg +! above). +! 67 = bad first parameter to deflt. +! 80 = iv(1) was out of range. +! 81 = n is not positive. +! iv(g)........ iv(28) is the starting subscript in v of the current +! gradient vector (the one corresponding to x). +! iv(lastiv)... iv(44) is the least acceptable value of liv. (it is +! only set if liv is at least 44.) +! iv(lastv).... iv(45) is the least acceptable value of lv. (it is +! only set if liv is large enough, at least iv(lastiv).) +! iv(nfcall)... iv(6) is the number of calls so far made on calcf (i.e., +! function evaluations). +! iv(ngcall)... iv(30) is the number of gradient evaluations (calls on +! calcg). +! iv(niter).... iv(31) is the number of iterations performed. +! +! *** (selected) v input values (from subroutine deflt) *** +! +! v(bias)..... v(43) is the bias parameter used in subroutine dbldog -- +! see that subroutine for details. default = 0.8. +! v(afctol)... v(31) is the absolute function convergence tolerance. +! if sumsl finds a point where the function value is less +! than v(afctol) in absolute value, and if sumsl does not +! return with iv(1) = 3, 4, or 5, then it returns with +! iv(1) = 6. this test can be turned off by setting +! v(afctol) to zero. default = max(10**-20, machep**2), +! where machep is the unit roundoff. +! v(dinit).... v(38), if nonnegative, is the value to which the scale +! vector d is initialized. default = -1. +! v(lmax0).... v(35) gives the maximum 2-norm allowed for d times the +! very first step that sumsl attempts. this parameter can +! markedly affect the performance of sumsl. +! v(lmaxs).... v(36) is used in testing for singular convergence -- if +! the function reduction predicted for a step of length +! bounded by v(lmaxs) is at most v(sctol) * abs(f0), where +! f0 is the function value at the start of the current +! iteration, and if sumsl does not return with iv(1) = 3, +! 4, 5, or 6, then it returns with iv(1) = 7. default = 1. +! v(rfctol)... v(32) is the relative function convergence tolerance. +! if the current model predicts a maximum possible function +! reduction (see v(nreduc)) of at most v(rfctol)*abs(f0) +! at the start of the current iteration, where f0 is the +! then current function value, and if the last step attempt- +! ed achieved no more than twice the predicted function +! decrease, then sumsl returns with iv(1) = 4 (or 5). +! default = max(10**-10, machep**(2/3)), where machep is +! the unit roundoff. +! v(sctol).... v(37) is the singular convergence tolerance -- see the +! description of v(lmaxs) above. +! v(tuner1)... v(26) helps decide when to check for false convergence. +! this is done if the actual function decrease from the +! current step is no more than v(tuner1) times its predict- +! ed value. default = 0.1. +! v(xctol).... v(33) is the x-convergence tolerance. if a newton step +! (see v(nreduc)) is tried that has v(reldx) .le. v(xctol) +! and if this step yields at most twice the predicted func- +! tion decrease, then sumsl returns with iv(1) = 3 (or 5). +! (see the description of v(reldx) below.) +! default = machep**0.5, where machep is the unit roundoff. +! v(xftol).... v(34) is the false convergence tolerance. if a step is +! tried that gives no more than v(tuner1) times the predict- +! ed function decrease and that has v(reldx) .le. v(xftol), +! and if sumsl does not return with iv(1) = 3, 4, 5, 6, or +! 7, then it returns with iv(1) = 8. (see the description +! of v(reldx) below.) default = 100*machep, where +! machep is the unit roundoff. +! v(*)........ deflt supplies to v a number of tuning constants, with +! which it should ordinarily be unnecessary to tinker. see +! section 17 of version 2.2 of the nl2sol usage summary +! (i.e., the appendix to ref. 1) for details on v(i), +! i = decfac, incfac, phmnfc, phmxfc, rdfcmn, rdfcmx, +! tuner2, tuner3, tuner4, tuner5. +! +! *** (selected) v output values *** +! +! v(dgnorm)... v(1) is the 2-norm of (diag(d)**-1)*g, where g is the +! most recently computed gradient. +! v(dstnrm)... v(2) is the 2-norm of diag(d)*step, where step is the +! current step. +! v(f)........ v(10) is the current function value. +! v(f0)....... v(13) is the function value at the start of the current +! iteration. +! v(nreduc)... v(6), if positive, is the maximum function reduction +! possible according to the current model, i.e., the func- +! tion reduction predicted for a newton step (i.e., +! step = -h**-1 * g, where g is the current gradient and +! h is the current hessian approximation). +! if v(nreduc) is negative, then it is the negative of +! the function reduction predicted for a step computed with +! a step bound of v(lmaxs) for use in testing for singular +! convergence. +! v(preduc)... v(7) is the function reduction predicted (by the current +! quadratic model) for the current step. this (divided by +! v(f0)) is used in testing for relative function +! convergence. +! v(reldx).... v(17) is the scaled relative change in x caused by the +! current step, computed as +! max(abs(d(i)*(x(i)-x0(i)), 1 .le. i .le. p) / +! max(d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p), +! where x = x0 + step. +! +!------------------------------- notes ------------------------------- +! +! *** algorithm notes *** +! +! this routine uses a hessian approximation computed from the +! bfgs update (see ref 3). only a cholesky factor of the hessian +! approximation is stored, and this is updated using ideas from +! ref. 4. steps are computed by the double dogleg scheme described +! in ref. 2. the steps are assessed as in ref. 1. +! +! *** usage notes *** +! +! after a return with iv(1) .le. 11, it is possible to restart, +! i.e., to change some of the iv and v input values described above +! and continue the algorithm from the point where it was interrupt- +! ed. iv(1) should not be changed, nor should any entries of i +! and v other than the input values (those supplied by deflt). +! those who do not wish to write a calcg which computes the +! gradient analytically should call smsno rather than sumsl. +! smsno uses finite differences to compute an approximate gradient. +! those who would prefer to provide f and g (the function and +! gradient) by reverse communication rather than by writing subrou- +! tines calcf and calcg may call on sumit directly. see the com- +! ments at the beginning of sumit. +! those who use sumsl interactively may wish to supply their +! own stopx function, which should return .true. if the break key +! has been pressed since stopx was last invoked. this makes it +! possible to externally interrupt sumsl (which will return with +! iv(1) = 11 if stopx returns .true.). +! storage for g is allocated at the end of v. thus the caller +! may make v longer than specified above and may allow calcg to use +! elements of g beyond the first n as scratch storage. +! +! *** portability notes *** +! +! the sumsl distribution tape contains both single- and double- +! precision versions of the sumsl source code, so it should be un- +! necessary to change precisions. +! only the functions imdcon and rmdcon contain machine-dependent +! constants. to change from one machine to another, it should +! suffice to change the (few) relevant lines in these functions. +! intrinsic functions are explicitly declared. on certain com- +! puters (e.g. univac), it may be necessary to comment out these +! declarations. so that this may be done automatically by a simple +! program, such declarations are preceded by a comment having c/+ +! in columns 1-3 and blanks in columns 4-72 and are followed by +! a comment having c/ in columns 1 and 2 and blanks in columns 3-72. +! the sumsl source code is expressed in 1966 ansi standard +! fortran. it may be converted to fortran 77 by commenting out all +! lines that fall between a line having c/6 in columns 1-3 and a +! line having c/7 in columns 1-3 and by removing (i.e., replacing +! by a blank) the c in column 1 of the lines that follow the c/7 +! line and precede a line having c/ in columns 1-2 and blanks in +! columns 3-72. these changes convert some data statements into +! parameter statements, convert some variables from real to +! character*4, and make the data statements that initialize these +! variables use character strings delimited by primes instead +! of hollerith constants. (such variables and data statements +! appear only in modules itsum and parck. parameter statements +! appear nearly everywhere.) these changes also add save state- +! ments for variables given machine-dependent constants by rmdcon. +! +! *** references *** +! +! 1. dennis, j.e., gay, d.m., and welsch, r.e. (1981), algorithm 573 -- +! an adaptive nonlinear least-squares algorithm, acm trans. +! math. software 7, pp. 369-383. +! +! 2. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti- +! mization algorithms which use function and gradient +! values, j. optim. theory applic. 28, pp. 453-482. +! +! 3. dennis, j.e., and more, j.j. (1977), quasi-newton methods, motiva- +! tion and theory, siam rev. 19, pp. 46-89. +! +! 4. goldfarb, d. (1976), factorized variable metric methods for uncon- +! strained optimization, math. comput. 30, pp. 796-811. +! +! *** general *** +! +! coded by david m. gay (winter 1980). revised summer 1982. +! this subroutine was written in connection with research +! supported in part by the national science foundation under +! grants mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, +! and mcs-7906671. +!. +! +!---------------------------- declarations --------------------------- +! +!el external deflt, sumit +! +! deflt... supplies default iv and v input components. +! sumit... reverse-communication routine that carries out sumsl algo- +! rithm. +! + integer :: g1, iv1, nf + real(kind=8) :: f +! +! *** subscripts for iv *** +! +!el integer nextv, nfcall, nfgcal, g, toobig, vneed +! +!/6 +! data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, toobig/2/, vneed/4/ +!/7 + integer,parameter :: nextv=47, nfcall=6, nfgcal=7, g=28,& + toobig=2, vneed=4 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! +!elwrite(iout,*) "in sumsl" + if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v) + iv1 = iv(1) + if (iv1 .eq. 12 .or. iv1 .eq. 13) iv(vneed) = iv(vneed) + n + if (iv1 .eq. 14) go to 10 + if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10 + g1 = 1 + if (iv1 .eq. 12) iv(1) = 13 + go to 20 +! + 10 g1 = iv(g) +!elwrite(iout,*) "in sumsl go to 10" + +! +!elwrite(iout,*) "in sumsl" + 20 call sumit(d, f, v(g1), iv, liv, lv, n, v, x) +!elwrite(iout,*) "in sumsl, go to 20" + +!elwrite(iout,*) "in sumsl, go to 20, po sumit" +!elwrite(iout,*) "in sumsl iv()", iv(1)-2 + if (iv(1) - 2) 30, 40, 50 +! + 30 nf = iv(nfcall) +!elwrite(iout,*) "in sumsl iv",iv(nfcall) + call calcf(n, x, nf, f, uiparm, urparm, ufparm) +!elwrite(iout,*) "in sumsl" + if (nf .le. 0) iv(toobig) = 1 + go to 20 +! +!elwrite(iout,*) "in sumsl" + 40 call calcg(n, x, iv(nfgcal), v(g1), uiparm, urparm, ufparm) +!elwrite(iout,*) "in sumsl" + go to 20 +! + 50 if (iv(1) .ne. 14) go to 999 +! +! *** storage allocation +! + iv(g) = iv(nextv) + iv(nextv) = iv(g) + n + if (iv1 .ne. 13) go to 10 +!elwrite(iout,*) "in sumsl" +! + 999 return +! *** last card of sumsl follows *** + end subroutine sumsl +!----------------------------------------------------------------------------- + subroutine sumit(d,fx,g,iv,liv,lv,n,v,x) + + use control, only:stopx +! +! *** carry out sumsl (unconstrained minimization) iterations, using +! *** double-dogleg/bfgs steps. +! +! *** parameter declarations *** +! + integer :: liv, lv, n + integer :: iv(liv) + real(kind=8) :: d(n), fx, g(n), v(lv), x(n) +! +!-------------------------- parameter usage -------------------------- +! +! d.... scale vector. +! fx... function value. +! g.... gradient vector. +! iv... integer value array. +! liv.. length of iv (at least 60). +! lv... length of v (at least 71 + n*(n+13)/2). +! n.... number of variables (components in x and g). +! v.... floating-point value array. +! x.... vector of parameters to be optimized. +! +! *** discussion *** +! +! parameters iv, n, v, and x are the same as the corresponding +! ones to sumsl (which see), except that v can be shorter (since +! the part of v that sumsl uses for storing g is not needed). +! moreover, compared with sumsl, iv(1) may have the two additional +! output values 1 and 2, which are explained below, as is the use +! of iv(toobig) and iv(nfgcal). the value iv(g), which is an +! output value from sumsl (and smsno), is not referenced by +! sumit or the subroutines it calls. +! fx and g need not have been initialized when sumit is called +! with iv(1) = 12, 13, or 14. +! +! iv(1) = 1 means the caller should set fx to f(x), the function value +! at x, and call sumit again, having changed none of the +! other parameters. an exception occurs if f(x) cannot be +! (e.g. if overflow would occur), which may happen because +! of an oversized step. in this case the caller should set +! iv(toobig) = iv(2) to 1, which will cause sumit to ig- +! nore fx and try a smaller step. the parameter nf that +! sumsl passes to calcf (for possible use by calcg) is a +! copy of iv(nfcall) = iv(6). +! iv(1) = 2 means the caller should set g to g(x), the gradient vector +! of f at x, and call sumit again, having changed none of +! the other parameters except possibly the scale vector d +! when iv(dtype) = 0. the parameter nf that sumsl passes +! to calcg is iv(nfgcal) = iv(7). if g(x) cannot be +! evaluated, then the caller may set iv(nfgcal) to 0, in +! which case sumit will return with iv(1) = 65. +!. +! *** general *** +! +! coded by david m. gay (december 1979). revised sept. 1982. +! this subroutine was written in connection with research supported +! in part by the national science foundation under grants +! mcs-7600324 and mcs-7906671. +! +! (see sumsl for references.) +! +!+++++++++++++++++++++++++++ declarations ++++++++++++++++++++++++++++ +! +! *** local variables *** +! + integer :: dg1, dummy, g01, i, k, l, lstgst, nwtst1, step1,& + temp1, w, x01, z + real(kind=8) :: t +!el logical :: lstopx +! +! *** constants *** +! +!el real(kind=8) :: half, negone, one, onep2, zero +! +! *** no intrinsic functions *** +! +! *** external functions and subroutines *** +! +!el external assst, dbdog, deflt, dotprd, itsum, litvmu, livmul, +!el 1 ltvmul, lupdat, lvmul, parck, reldst, stopx, vaxpy, +!el 2 vcopy, vscopy, vvmulp, v2norm, wzbfgs +!el logical stopx +!el real(kind=8) :: dotprd, reldst, v2norm +! +! assst.... assesses candidate step. +! dbdog.... computes double-dogleg (candidate) step. +! deflt.... supplies default iv and v input components. +! dotprd... returns inner product of two vectors. +! itsum.... prints iteration summary and info on initial and final x. +! litvmu... multiplies inverse transpose of lower triangle times vector. +! livmul... multiplies inverse of lower triangle times vector. +! ltvmul... multiplies transpose of lower triangle times vector. +! lupdt.... updates cholesky factor of hessian approximation. +! lvmul.... multiplies lower triangle times vector. +! parck.... checks validity of input iv and v values. +! reldst... computes v(reldx) = relative step size. +! stopx.... returns .true. if the break key has been pressed. +! vaxpy.... computes scalar times one vector plus another. +! vcopy.... copies one vector to another. +! vscopy... sets all elements of a vector to a scalar. +! vvmulp... multiplies vector by vector raised to power (componentwise). +! v2norm... returns the 2-norm of a vector. +! wzbfgs... computes w and z for lupdat corresponding to bfgs update. +! +! *** subscripts for iv and v *** +! +!el integer afctol +!el integer cnvcod, dg, dgnorm, dinit, dstnrm, dst0, f, f0, fdif, +!el 1 gthg, gtstep, g0, incfac, inith, irc, kagqt, lmat, lmax0, +!el 2 lmaxs, mode, model, mxfcal, mxiter, nextv, nfcall, nfgcal, +!el 3 ngcall, niter, nreduc, nwtstp, preduc, radfac, radinc, +!el 4 radius, rad0, reldx, restor, step, stglim, stlstg, toobig, +!el 5 tuner4, tuner5, vneed, xirc, x0 +! +! *** iv subscript values *** +! +!/6 +! data cnvcod/55/, dg/37/, g0/48/, inith/25/, irc/29/, kagqt/33/, +! 1 mode/35/, model/5/, mxfcal/17/, mxiter/18/, nfcall/6/, +! 2 nfgcal/7/, ngcall/30/, niter/31/, nwtstp/34/, radinc/8/, +! 3 restor/9/, step/40/, stglim/11/, stlstg/41/, toobig/2/, +! 4 vneed/4/, xirc/13/, x0/43/ +!/7 + integer,parameter :: cnvcod=55, dg=37, g0=48, inith=25, irc=29, kagqt=33,& + mode=35, model=5, mxfcal=17, mxiter=18, nfcall=6,& + nfgcal=7, ngcall=30, niter=31, nwtstp=34, radinc=8,& + restor=9, step=40, stglim=11, stlstg=41, toobig=2,& + vneed=4, xirc=13, x0=43 +!/ +! +! *** v subscript values *** +! +!/6 +! data afctol/31/ +! data dgnorm/1/, dinit/38/, dstnrm/2/, dst0/3/, f/10/, f0/13/, +! 1 fdif/11/, gthg/44/, gtstep/4/, incfac/23/, lmat/42/, +! 2 lmax0/35/, lmaxs/36/, nextv/47/, nreduc/6/, preduc/7/, +! 3 radfac/16/, radius/8/, rad0/9/, reldx/17/, tuner4/29/, +! 4 tuner5/30/ +!/7 + integer,parameter :: afctol=31 + integer,parameter :: dgnorm=1, dinit=38, dstnrm=2, dst0=3, f=10, f0=13,& + fdif=11, gthg=44, gtstep=4, incfac=23, lmat=42,& + lmax0=35, lmaxs=36, nextv=47, nreduc=6, preduc=7,& + radfac=16, radius=8, rad0=9, reldx=17, tuner4=29,& + tuner5=30 +!/ +! +!/6 +! data half/0.5d+0/, negone/-1.d+0/, one/1.d+0/, onep2/1.2d+0/, +! 1 zero/0.d+0/ +!/7 + real(kind=8),parameter :: half=0.5d+0, negone=-1.d+0, one=1.d+0,& + onep2=1.2d+0,zero=0.d+0 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! +! Following SAVE statement inserted. + save l + i = iv(1) + if (i .eq. 1) go to 50 + if (i .eq. 2) go to 60 +! +! *** check validity of iv and v input values *** +! + if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v) + if (iv(1) .eq. 12 .or. iv(1) .eq. 13) & + iv(vneed) = iv(vneed) + n*(n+13)/2 + call parck(2, d, iv, liv, lv, n, v) + i = iv(1) - 2 + if (i .gt. 12) go to 999 + go to (180, 180, 180, 180, 180, 180, 120, 90, 120, 10, 10, 20), i +! +! *** storage allocation *** +! +10 l = iv(lmat) + iv(x0) = l + n*(n+1)/2 + iv(step) = iv(x0) + n + iv(stlstg) = iv(step) + n + iv(g0) = iv(stlstg) + n + iv(nwtstp) = iv(g0) + n + iv(dg) = iv(nwtstp) + n + iv(nextv) = iv(dg) + n + if (iv(1) .ne. 13) go to 20 + iv(1) = 14 + go to 999 +! +! *** initialization *** +! + 20 iv(niter) = 0 + iv(nfcall) = 1 + iv(ngcall) = 1 + iv(nfgcal) = 1 + iv(mode) = -1 + iv(model) = 1 + iv(stglim) = 1 + iv(toobig) = 0 + iv(cnvcod) = 0 + iv(radinc) = 0 + v(rad0) = zero + if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit)) + if (iv(inith) .ne. 1) go to 40 +! +! *** set the initial hessian approximation to diag(d)**-2 *** +! + l = iv(lmat) + call vscopy(n*(n+1)/2, v(l), zero) + k = l - 1 + do 30 i = 1, n + k = k + i + t = d(i) + if (t .le. zero) t = one + v(k) = t + 30 continue +! +! *** compute initial function value *** +! + 40 iv(1) = 1 + go to 999 +! + 50 v(f) = fx + if (iv(mode) .ge. 0) go to 180 + iv(1) = 2 + if (iv(toobig) .eq. 0) go to 999 + iv(1) = 63 + go to 300 +! +! *** make sure gradient could be computed *** +! + 60 if (iv(nfgcal) .ne. 0) go to 70 + iv(1) = 65 + go to 300 +! + 70 dg1 = iv(dg) + call vvmulp(n, v(dg1), g, d, -1) + v(dgnorm) = v2norm(n, v(dg1)) +! +! *** test norm of gradient *** +! + if (v(dgnorm) .gt. v(afctol)) go to 75 + iv(irc) = 10 + iv(cnvcod) = iv(irc) - 4 +! + 75 if (iv(cnvcod) .ne. 0) go to 290 + if (iv(mode) .eq. 0) go to 250 +! +! *** allow first step to have scaled 2-norm at most v(lmax0) *** +! + v(radius) = v(lmax0) +! + iv(mode) = 0 +! +! +!----------------------------- main loop ----------------------------- +! +! +! *** print iteration summary, check iteration limit *** +! + 80 call itsum(d, g, iv, liv, lv, n, v, x) + 90 k = iv(niter) + if (k .lt. iv(mxiter)) go to 100 + iv(1) = 10 + go to 300 +! +! *** update radius *** +! + 100 iv(niter) = k + 1 + if(k.gt.0)v(radius) = v(radfac) * v(dstnrm) +! +! *** initialize for start of next iteration *** +! + g01 = iv(g0) + x01 = iv(x0) + v(f0) = v(f) + iv(irc) = 4 + iv(kagqt) = -1 +! +! *** copy x to x0, g to g0 *** +! + call vcopy(n, v(x01), x) + call vcopy(n, v(g01), g) +! +! *** check stopx and function evaluation limit *** +! +! AL 4/30/95 + dummy=iv(nfcall) +!el lstopx = stopx(dummy) +!elwrite(iout,*) "lstopx",lstopx,dummy + 110 if (.not. stopx(dummy)) go to 130 + iv(1) = 11 +! write (iout,*) "iv(1)=11 !!!!" + go to 140 +! +! *** come here when restarting after func. eval. limit or stopx. +! + 120 if (v(f) .ge. v(f0)) go to 130 + v(radfac) = one + k = iv(niter) + go to 100 +! + 130 if (iv(nfcall) .lt. iv(mxfcal)) go to 150 + iv(1) = 9 + 140 if (v(f) .ge. v(f0)) go to 300 +! +! *** in case of stopx or function evaluation limit with +! *** improved v(f), evaluate the gradient at x. +! + iv(cnvcod) = iv(1) + go to 240 +! +!. . . . . . . . . . . . . compute candidate step . . . . . . . . . . +! + 150 step1 = iv(step) + dg1 = iv(dg) + nwtst1 = iv(nwtstp) + if (iv(kagqt) .ge. 0) go to 160 + l = iv(lmat) + call livmul(n, v(nwtst1), v(l), g) + v(nreduc) = half * dotprd(n, v(nwtst1), v(nwtst1)) + call litvmu(n, v(nwtst1), v(l), v(nwtst1)) + call vvmulp(n, v(step1), v(nwtst1), d, 1) + v(dst0) = v2norm(n, v(step1)) + call vvmulp(n, v(dg1), v(dg1), d, -1) + call ltvmul(n, v(step1), v(l), v(dg1)) + v(gthg) = v2norm(n, v(step1)) + iv(kagqt) = 0 + 160 call dbdog(v(dg1), lv, n, v(nwtst1), v(step1), v) + if (iv(irc) .eq. 6) go to 180 +! +! *** check whether evaluating f(x0 + step) looks worthwhile *** +! + if (v(dstnrm) .le. zero) go to 180 + if (iv(irc) .ne. 5) go to 170 + if (v(radfac) .le. one) go to 170 + if (v(preduc) .le. onep2 * v(fdif)) go to 180 +! +! *** compute f(x0 + step) *** +! + 170 x01 = iv(x0) + step1 = iv(step) + call vaxpy(n, x, one, v(step1), v(x01)) + iv(nfcall) = iv(nfcall) + 1 + iv(1) = 1 + iv(toobig) = 0 + go to 999 +! +!. . . . . . . . . . . . . assess candidate step . . . . . . . . . . . +! + 180 x01 = iv(x0) + v(reldx) = reldst(n, d, x, v(x01)) + call assst(iv, liv, lv, v) + step1 = iv(step) + lstgst = iv(stlstg) + if (iv(restor) .eq. 1) call vcopy(n, x, v(x01)) + if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1)) + if (iv(restor) .ne. 3) go to 190 + call vcopy(n, v(step1), v(lstgst)) + call vaxpy(n, x, one, v(step1), v(x01)) + v(reldx) = reldst(n, d, x, v(x01)) +! + 190 k = iv(irc) + go to (200,230,230,230,200,210,220,220,220,220,220,220,280,250), k +! +! *** recompute step with changed radius *** +! + 200 v(radius) = v(radfac) * v(dstnrm) + go to 110 +! +! *** compute step of length v(lmaxs) for singular convergence test. +! + 210 v(radius) = v(lmaxs) + go to 150 +! +! *** convergence or false convergence *** +! + 220 iv(cnvcod) = k - 4 + if (v(f) .ge. v(f0)) go to 290 + if (iv(xirc) .eq. 14) go to 290 + iv(xirc) = 14 +! +!. . . . . . . . . . . . process acceptable step . . . . . . . . . . . +! + 230 if (iv(irc) .ne. 3) go to 240 + step1 = iv(step) + temp1 = iv(stlstg) +! +! *** set temp1 = hessian * step for use in gradient tests *** +! + l = iv(lmat) + call ltvmul(n, v(temp1), v(l), v(step1)) + call lvmul(n, v(temp1), v(l), v(temp1)) +! +! *** compute gradient *** +! + 240 iv(ngcall) = iv(ngcall) + 1 + iv(1) = 2 + go to 999 +! +! *** initializations -- g0 = g - g0, etc. *** +! + 250 g01 = iv(g0) + call vaxpy(n, v(g01), negone, v(g01), g) + step1 = iv(step) + temp1 = iv(stlstg) + if (iv(irc) .ne. 3) go to 270 +! +! *** set v(radfac) by gradient tests *** +! +! *** set temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x))) *** +! + call vaxpy(n, v(temp1), negone, v(g01), v(temp1)) + call vvmulp(n, v(temp1), v(temp1), d, -1) +! +! *** do gradient tests *** +! + if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) & + go to 260 + if (dotprd(n, g, v(step1)) & + .ge. v(gtstep) * v(tuner5)) go to 270 + 260 v(radfac) = v(incfac) +! +! *** update h, loop *** +! + 270 w = iv(nwtstp) + z = iv(x0) + l = iv(lmat) + call wzbfgs(v(l), n, v(step1), v(w), v(g01), v(z)) +! +! ** use the n-vectors starting at v(step1) and v(g01) for scratch.. + call lupdat(v(temp1), v(step1), v(l), v(g01), v(l), n, v(w), v(z)) + iv(1) = 2 + go to 80 +! +!. . . . . . . . . . . . . . misc. details . . . . . . . . . . . . . . +! +! *** bad parameters to assess *** +! + 280 iv(1) = 64 + go to 300 +! +! *** print summary of final iteration and other requested items *** +! + 290 iv(1) = iv(cnvcod) + iv(cnvcod) = 0 + 300 call itsum(d, g, iv, liv, lv, n, v, x) +! + 999 return +! +! *** last line of sumit follows *** + end subroutine sumit +!----------------------------------------------------------------------------- + subroutine dbdog(dig,lv,n,nwtstp,step,v) +! +! *** compute double dogleg step *** +! +! *** parameter declarations *** +! + integer :: lv, n + real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv) +! +! *** purpose *** +! +! this subroutine computes a candidate step (for use in an uncon- +! strained minimization code) by the double dogleg algorithm of +! dennis and mei (ref. 1), which is a variation on powell*s dogleg +! scheme (ref. 2, p. 95). +! +!-------------------------- parameter usage -------------------------- +! +! dig (input) diag(d)**-2 * g -- see algorithm notes. +! g (input) the current gradient vector. +! lv (input) length of v. +! n (input) number of components in dig, g, nwtstp, and step. +! nwtstp (input) negative newton step -- see algorithm notes. +! step (output) the computed step. +! v (i/o) values array, the following components of which are +! used here... +! v(bias) (input) bias for relaxed newton step, which is v(bias) of +! the way from the full newton to the fully relaxed newton +! step. recommended value = 0.8 . +! v(dgnorm) (input) 2-norm of diag(d)**-1 * g -- see algorithm notes. +! v(dstnrm) (output) 2-norm of diag(d) * step, which is v(radius) +! unless v(stppar) = 0 -- see algorithm notes. +! v(dst0) (input) 2-norm of diag(d) * nwtstp -- see algorithm notes. +! v(grdfac) (output) the coefficient of dig in the step returned -- +! step(i) = v(grdfac)*dig(i) + v(nwtfac)*nwtstp(i). +! v(gthg) (input) square-root of (dig**t) * (hessian) * dig -- see +! algorithm notes. +! v(gtstep) (output) inner product between g and step. +! v(nreduc) (output) function reduction predicted for the full newton +! step. +! v(nwtfac) (output) the coefficient of nwtstp in the step returned -- +! see v(grdfac) above. +! v(preduc) (output) function reduction predicted for the step returned. +! v(radius) (input) the trust region radius. d times the step returned +! has 2-norm v(radius) unless v(stppar) = 0. +! v(stppar) (output) code telling how step was computed... 0 means a +! full newton step. between 0 and 1 means v(stppar) of the +! way from the newton to the relaxed newton step. between +! 1 and 2 means a true double dogleg step, v(stppar) - 1 of +! the way from the relaxed newton to the cauchy step. +! greater than 2 means 1 / (v(stppar) - 1) times the cauchy +! step. +! +!------------------------------- notes ------------------------------- +! +! *** algorithm notes *** +! +! let g and h be the current gradient and hessian approxima- +! tion respectively and let d be the current scale vector. this +! routine assumes dig = diag(d)**-2 * g and nwtstp = h**-1 * g. +! the step computed is the same one would get by replacing g and h +! by diag(d)**-1 * g and diag(d)**-1 * h * diag(d)**-1, +! computing step, and translating step back to the original +! variables, i.e., premultiplying it by diag(d)**-1. +! +! *** references *** +! +! 1. dennis, j.e., and mei, h.h.w. (1979), two new unconstrained opti- +! mization algorithms which use function and gradient +! values, j. optim. theory applic. 28, pp. 453-482. +! 2. powell, m.j.d. (1970), a hybrid method for non-linear equations, +! in numerical methods for non-linear equations, edited by +! p. rabinowitz, gordon and breach, london. +! +! *** general *** +! +! coded by david m. gay. +! this subroutine was written in connection with research supported +! by the national science foundation under grants mcs-7600324 and +! mcs-7906671. +! +!------------------------ external quantities ------------------------ +! +! *** functions and subroutines called *** +! +!el external dotprd, v2norm +!el real(kind=8) :: dotprd, v2norm +! +! dotprd... returns inner product of two vectors. +! v2norm... returns 2-norm of a vector. +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dsqrt +!/ +!-------------------------- local variables -------------------------- +! + integer :: i + real(kind=8) :: cfact, cnorm, ctrnwt, ghinvg, femnsq, gnorm,& + nwtnrm, relax, rlambd, t, t1, t2 +!el real(kind=8) :: half, one, two, zero +! +! *** v subscripts *** +! +!el integer bias, dgnorm, dstnrm, dst0, grdfac, gthg, gtstep, +!el 1 nreduc, nwtfac, preduc, radius, stppar +! +! *** data initializations *** +! +!/6 +! data half/0.5d+0/, one/1.d+0/, two/2.d+0/, zero/0.d+0/ +!/7 + real(kind=8),parameter :: half=0.5d+0, one=1.d+0, two=2.d+0, zero=0.d+0 +!/ +! +!/6 +! data bias/43/, dgnorm/1/, dstnrm/2/, dst0/3/, grdfac/45/, +! 1 gthg/44/, gtstep/4/, nreduc/6/, nwtfac/46/, preduc/7/, +! 2 radius/8/, stppar/5/ +!/7 + integer,parameter :: bias=43, dgnorm=1, dstnrm=2, dst0=3, grdfac=45,& + gthg=44, gtstep=4, nreduc=6, nwtfac=46, preduc=7,& + radius=8, stppar=5 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! + nwtnrm = v(dst0) + rlambd = one + if (nwtnrm .gt. zero) rlambd = v(radius) / nwtnrm + gnorm = v(dgnorm) + ghinvg = two * v(nreduc) + v(grdfac) = zero + v(nwtfac) = zero + if (rlambd .lt. one) go to 30 +! +! *** the newton step is inside the trust region *** +! + v(stppar) = zero + v(dstnrm) = nwtnrm + v(gtstep) = -ghinvg + v(preduc) = v(nreduc) + v(nwtfac) = -one + do 20 i = 1, n + 20 step(i) = -nwtstp(i) + go to 999 +! + 30 v(dstnrm) = v(radius) + cfact = (gnorm / v(gthg))**2 +! *** cauchy step = -cfact * g. + cnorm = gnorm * cfact + relax = one - v(bias) * (one - gnorm*cnorm/ghinvg) + if (rlambd .lt. relax) go to 50 +! +! *** step is between relaxed newton and full newton steps *** +! + v(stppar) = one - (rlambd - relax) / (one - relax) + t = -rlambd + v(gtstep) = t * ghinvg + v(preduc) = rlambd * (one - half*rlambd) * ghinvg + v(nwtfac) = t + do 40 i = 1, n + 40 step(i) = t * nwtstp(i) + go to 999 +! + 50 if (cnorm .lt. v(radius)) go to 70 +! +! *** the cauchy step lies outside the trust region -- +! *** step = scaled cauchy step *** +! + t = -v(radius) / gnorm + v(grdfac) = t + v(stppar) = one + cnorm / v(radius) + v(gtstep) = -v(radius) * gnorm + v(preduc) = v(radius)*(gnorm - half*v(radius)*(v(gthg)/gnorm)**2) + do 60 i = 1, n + 60 step(i) = t * dig(i) + go to 999 +! +! *** compute dogleg step between cauchy and relaxed newton *** +! *** femur = relaxed newton step minus cauchy step *** +! + 70 ctrnwt = cfact * relax * ghinvg / gnorm +! *** ctrnwt = inner prod. of cauchy and relaxed newton steps, +! *** scaled by gnorm**-1. + t1 = ctrnwt - gnorm*cfact**2 +! *** t1 = inner prod. of femur and cauchy step, scaled by +! *** gnorm**-1. + t2 = v(radius)*(v(radius)/gnorm) - gnorm*cfact**2 + t = relax * nwtnrm + femnsq = (t/gnorm)*t - ctrnwt - t1 +! *** femnsq = square of 2-norm of femur, scaled by gnorm**-1. + t = t2 / (t1 + dsqrt(t1**2 + femnsq*t2)) +! *** dogleg step = cauchy step + t * femur. + t1 = (t - one) * cfact + v(grdfac) = t1 + t2 = -t * relax + v(nwtfac) = t2 + v(stppar) = two - t + v(gtstep) = t1*gnorm**2 + t2*ghinvg + v(preduc) = -t1*gnorm * ((t2 + one)*gnorm) & + - t2 * (one + half*t2)*ghinvg & + - half * (v(gthg)*t1)**2 + do 80 i = 1, n + 80 step(i) = t1*dig(i) + t2*nwtstp(i) +! + 999 return +! *** last line of dbdog follows *** + end subroutine dbdog +!----------------------------------------------------------------------------- + subroutine ltvmul(n,x,l,y) +! +! *** compute x = (l**t)*y, where l is an n x n lower +! *** triangular matrix stored compactly by rows. x and y may +! *** occupy the same storage. *** +! + integer :: n +!al real(kind=8) :: x(n), l(1), y(n) + real(kind=8) :: x(n), l(n*(n+1)/2), y(n) +! dimension l(n*(n+1)/2) + integer :: i, ij, i0, j + real(kind=8) :: yi !el, zero +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +! + i0 = 0 + do 20 i = 1, n + yi = y(i) + x(i) = zero + do 10 j = 1, i + ij = i0 + j + x(j) = x(j) + yi*l(ij) + 10 continue + i0 = i0 + i + 20 continue + 999 return +! *** last card of ltvmul follows *** + end subroutine ltvmul +!----------------------------------------------------------------------------- + subroutine lupdat(beta,gamma,l,lambda,lplus,n,w,z) +! +! *** compute lplus = secant update of l *** +! +! *** parameter declarations *** +! + integer :: n +!al double precision beta(n), gamma(n), l(1), lambda(n), lplus(1), + real(kind=8) :: beta(n), gamma(n), l(n*(n+1)/2), lambda(n), & + lplus(n*(n+1)/2),w(n), z(n) +! dimension l(n*(n+1)/2), lplus(n*(n+1)/2) +! +!-------------------------- parameter usage -------------------------- +! +! beta = scratch vector. +! gamma = scratch vector. +! l (input) lower triangular matrix, stored rowwise. +! lambda = scratch vector. +! lplus (output) lower triangular matrix, stored rowwise, which may +! occupy the same storage as l. +! n (input) length of vector parameters and order of matrices. +! w (input, destroyed on output) right singular vector of rank 1 +! correction to l. +! z (input, destroyed on output) left singular vector of rank 1 +! correction to l. +! +!------------------------------- notes ------------------------------- +! +! *** application and usage restrictions *** +! +! this routine updates the cholesky factor l of a symmetric +! positive definite matrix to which a secant update is being +! applied -- it computes a cholesky factor lplus of +! l * (i + z*w**t) * (i + w*z**t) * l**t. it is assumed that w +! and z have been chosen so that the updated matrix is strictly +! positive definite. +! +! *** algorithm notes *** +! +! this code uses recurrence 3 of ref. 1 (with d(j) = 1 for all j) +! to compute lplus of the form l * (i + z*w**t) * q, where q +! is an orthogonal matrix that makes the result lower triangular. +! lplus may have some negative diagonal elements. +! +! *** references *** +! +! 1. goldfarb, d. (1976), factorized variable metric methods for uncon- +! strained optimization, math. comput. 30, pp. 796-811. +! +! *** general *** +! +! coded by david m. gay (fall 1979). +! this subroutine was written in connection with research supported +! by the national science foundation under grants mcs-7600324 and +! mcs-7906671. +! +!------------------------ external quantities ------------------------ +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dsqrt +!/ +!-------------------------- local variables -------------------------- +! + integer :: i, ij, j, jj, jp1, k, nm1, np1 + real(kind=8) :: a, b, bj, eta, gj, lj, lij, ljj, nu, s, theta,& + wj, zj +!el real(kind=8) :: one, zero +! +! *** data initializations *** +! +!/6 +! data one/1.d+0/, zero/0.d+0/ +!/7 + real(kind=8),parameter :: one=1.d+0, zero=0.d+0 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! + nu = one + eta = zero + if (n .le. 1) go to 30 + nm1 = n - 1 +! +! *** temporarily store s(j) = sum over k = j+1 to n of w(k)**2 in +! *** lambda(j). +! + s = zero + do 10 i = 1, nm1 + j = n - i + s = s + w(j+1)**2 + lambda(j) = s + 10 continue +! +! *** compute lambda, gamma, and beta by goldfarb*s recurrence 3. +! + do 20 j = 1, nm1 + wj = w(j) + a = nu*z(j) - eta*wj + theta = one + a*wj + s = a*lambda(j) + lj = dsqrt(theta**2 + a*s) + if (theta .gt. zero) lj = -lj + lambda(j) = lj + b = theta*wj + s + gamma(j) = b * nu / lj + beta(j) = (a - b*eta) / lj + nu = -nu / lj + eta = -(eta + (a**2)/(theta - lj)) / lj + 20 continue + 30 lambda(n) = one + (nu*z(n) - eta*w(n))*w(n) +! +! *** update l, gradually overwriting w and z with l*w and l*z. +! + np1 = n + 1 + jj = n * (n + 1) / 2 + do 60 k = 1, n + j = np1 - k + lj = lambda(j) + ljj = l(jj) + lplus(jj) = lj * ljj + wj = w(j) + w(j) = ljj * wj + zj = z(j) + z(j) = ljj * zj + if (k .eq. 1) go to 50 + bj = beta(j) + gj = gamma(j) + ij = jj + j + jp1 = j + 1 + do 40 i = jp1, n + lij = l(ij) + lplus(ij) = lj*lij + bj*w(i) + gj*z(i) + w(i) = w(i) + lij*wj + z(i) = z(i) + lij*zj + ij = ij + i + 40 continue + 50 jj = jj - j + 60 continue +! + 999 return +! *** last card of lupdat follows *** + end subroutine lupdat +!----------------------------------------------------------------------------- + subroutine lvmul(n,x,l,y) +! +! *** compute x = l*y, where l is an n x n lower triangular +! *** matrix stored compactly by rows. x and y may occupy the same +! *** storage. *** +! + integer :: n +!al double precision x(n), l(1), y(n) + real(kind=8) :: x(n), l(n*(n+1)/2), y(n) +! dimension l(n*(n+1)/2) + integer :: i, ii, ij, i0, j, np1 + real(kind=8) :: t !el, zero +!/6 +! data zero/0.d+0/ +!/7 + real(kind=8),parameter :: zero=0.d+0 +!/ +! + np1 = n + 1 + i0 = n*(n+1)/2 + do 20 ii = 1, n + i = np1 - ii + i0 = i0 - i + t = zero + do 10 j = 1, i + ij = i0 + j + t = t + l(ij)*y(j) + 10 continue + x(i) = t + 20 continue + 999 return +! *** last card of lvmul follows *** + end subroutine lvmul +!----------------------------------------------------------------------------- + subroutine vvmulp(n,x,y,z,k) +! +! *** set x(i) = y(i) * z(i)**k, 1 .le. i .le. n (for k = 1 or -1) *** +! + integer :: n, k + real(kind=8) :: x(n), y(n), z(n) + integer :: i +! + if (k .ge. 0) go to 20 + do 10 i = 1, n + 10 x(i) = y(i) / z(i) + go to 999 +! + 20 do 30 i = 1, n + 30 x(i) = y(i) * z(i) + 999 return +! *** last card of vvmulp follows *** + end subroutine vvmulp +!----------------------------------------------------------------------------- + subroutine wzbfgs(l,n,s,w,y,z) +! +! *** compute y and z for lupdat corresponding to bfgs update. +! + integer :: n +!al double precision l(1), s(n), w(n), y(n), z(n) + real(kind=8) :: l(n*(n+1)/2), s(n), w(n), y(n), z(n) +! dimension l(n*(n+1)/2) +! +!-------------------------- parameter usage -------------------------- +! +! l (i/o) cholesky factor of hessian, a lower triang. matrix stored +! compactly by rows. +! n (input) order of l and length of s, w, y, z. +! s (input) the step just taken. +! w (output) right singular vector of rank 1 correction to l. +! y (input) change in gradients corresponding to s. +! z (output) left singular vector of rank 1 correction to l. +! +!------------------------------- notes ------------------------------- +! +! *** algorithm notes *** +! +! when s is computed in certain ways, e.g. by gqtstp or +! dbldog, it is possible to save n**2/2 operations since (l**t)*s +! or l*(l**t)*s is then known. +! if the bfgs update to l*(l**t) would reduce its determinant to +! less than eps times its old value, then this routine in effect +! replaces y by theta*y + (1 - theta)*l*(l**t)*s, where theta +! (between 0 and 1) is chosen to make the reduction factor = eps. +! +! *** general *** +! +! coded by david m. gay (fall 1979). +! this subroutine was written in connection with research supported +! by the national science foundation under grants mcs-7600324 and +! mcs-7906671. +! +!------------------------ external quantities ------------------------ +! +! *** functions and subroutines called *** +! +!el external dotprd, livmul, ltvmul +!el real(kind=8) :: dotprd +! dotprd returns inner product of two vectors. +! livmul multiplies l**-1 times a vector. +! ltvmul multiplies l**t times a vector. +! +! *** intrinsic functions *** +!/+ +!el real(kind=8) :: dsqrt +!/ +!-------------------------- local variables -------------------------- +! + integer :: i + real(kind=8) :: cs, cy, epsrt, shs, ys, theta !el, eps, one +! +! *** data initializations *** +! +!/6 +! data eps/0.1d+0/, one/1.d+0/ +!/7 + real(kind=8),parameter :: eps=0.1d+0, one=1.d+0 +!/ +! +!+++++++++++++++++++++++++++++++ body ++++++++++++++++++++++++++++++++ +! + call ltvmul(n, w, l, s) + shs = dotprd(n, w, w) + ys = dotprd(n, y, s) + if (ys .ge. eps*shs) go to 10 + theta = (one - eps) * shs / (shs - ys) + epsrt = dsqrt(eps) + cy = theta / (shs * epsrt) + cs = (one + (theta-one)/epsrt) / shs + go to 20 + 10 cy = one / (dsqrt(ys) * dsqrt(shs)) + cs = one / shs + 20 call livmul(n, z, l, y) + do 30 i = 1, n + 30 z(i) = cy * z(i) - cs * w(i) +! + 999 return +! *** last card of wzbfgs follows *** + end subroutine wzbfgs +!----------------------------------------------------------------------------- +!----------------------------------------------------------------------------- + end module minimm