+++ /dev/null
- 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
projects / unres4.git / blobdiff