rename

[unres4.git] / source / unres / minim.F90

      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