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
!#ifdef LBFGS
!      double precision funcgrad_dc
!#endif
!-----------------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------------
      contains
!-----------------------------------------------------------------------------
! cored.f
!-----------------------------------------------------------------------------
      subroutine assst(iv, liv, lv, v)
!
!  ***  assess candidate step (***sol version 2.3)  ***
!
      integer :: liv, l
      integer ::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
      integer::lv
      integer :: iv(liv)
      integer :: alg 
      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, p
      integer::lv
      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, n
      integer :: lv
      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
      integer::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, n
      integer:: lv
      integer:: iv(liv)
      integer :: 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, n
      integer::lv
      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, n
      integer ::lv
      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)
#ifdef LBFGS
      use minima
      use inform
      use output
      use iounit
      use scales
      use energy, only: funcgrad,etotal,enerprint
#else
      use energy, only: func,gradient,fdum!,etotal,enerprint
#endif
      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),grdmin
!      external func,gradient,fdum
!      external func_restr,grad_restr
      logical :: not_done,change,reduce 
!el      common /przechowalnia/ v
!el local variables
      integer :: iretcode,nfun,nvar_restr,idum(1),j
      integer :: lv
      real(kind=8) :: etot,rdum(1)	!,fdum
#ifdef LBFGS
      external optsave
#endif
      lv=(77+(6*nres)*(6*nres+17)/2)	!(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
#ifndef LBFGS
      if (.not.allocated(v)) allocate(v(1:lv))
#endif
#ifdef LBFGS
      maxiter=maxmin
      coordtype='RIGIDBODY'
      grdmin=tolf
      jout=iout
      jprint=print_min_stat
      iwrite=0
      if (.not. allocated(scale))  allocate (scale(nvar))
!c
!c     set scaling parameter for function and derivative values;
!c     use square root of median eigenvalue of typical Hessian
!c
      set_scale = .true.
!c      nvar = 0
      do i = 1, nvar
!c         if (use(i)) then
!c            do j = 1, 3
!c               nvar = nvar + 1
               scale(i) = 12.0d0
!c            end do
!c         end if
      end do
!c      write (iout,*) "Calling lbfgs"
      write (iout,*) 'Calling LBFGS, minimization in angles'
      call var_to_geom(nvar,x)
      call chainbuild
      call etotal(energia(0))
      call enerprint(energia(0))
      call lbfgs (nvar,x,etot,grdmin,funcgrad,optsave)
      deallocate(scale)
      write (iout,*) "Minimized energy",etot
#else

      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
      write(iout,*) 'Warning calling chainbuild'
      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,1)),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
#endif
      return
      end subroutine minimize
!-----------------------------------------------------------------------------
! gradient_p.F
!-----------------------------------------------------------------------------
#ifndef LBFGS
      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)
      write(iout,*) 'Warning calling chainbuild'
      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,1).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,1).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,1).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
#endif
#ifndef LBFGS
!-----------------------------------------------------------------------------
      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
      write(iout,*) 'Warning calling chainbuild'
      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
#endif
!-----------------------------------------------------------------------------
!      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,1).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)
#ifdef LBFGS
      use minima
      use inform
      use output
      use iounit
      use scales
#endif

      use MPI_data
      use energy, only: fdum,check_ecartint,etotal,enerprint
!      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,idum(1)
      integer :: lv,nvarx,nf
      integer,dimension(liv) :: iv                                               
      real(kind=8) :: minval	!,v(1:77+(6*nres)*(6*nres+17)/2)	!(1:lv)
      real(kind=8),dimension(:),allocatable :: x,d,xx	!(maxvar) (maxvar=6*maxres)
!el      common /przechowalnia/ v

      real(kind=8) :: energia(0:n_ene),grdmin
!      external func_dc,grad_dc	,fdum
      logical :: not_done,change,reduce 
      real(kind=8) :: g(6*nres),f1,etot,rdum(1)	!,fdum
#ifdef LBFGS
      external optsave
      maxiter=maxmin
      coordtype='CARTESIAN'
      grdmin=tolf
      jout=iout
      jprint=print_min_stat
      iwrite=0
      write(iout,*) "in minim_dc LBFGS"
#else

      lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
      print *,"lv value",lv
      if (.not. allocated(v)) allocate(v(1:lv))

      call deflt(2,iv,liv,lv,v)                                         
      print *,"after delft"
! 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
#endif
      if (.not.allocated(x)) then
       allocate(x(6*nres))
       allocate(xx(6*nres))
       allocate(d(6*nres))
      endif
      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
      nvarx=k
      call chainbuild_cart
      call etotal(energia(0))
      call enerprint(energia(0))
#ifdef LBFGS
!c
!c     From tinker
!c
!c     perform dynamic allocation of some global arrays
!c
      if (.not. allocated(scale))  allocate (scale(nvarx))
!c
!c     set scaling parameter for function and derivative values;
!c     use square root of median eigenvalue of typical Hessian
!c
      set_scale = .true.
!c      nvar = 0
      do i = 1, nvarx
!c         if (use(i)) then
!c            do j = 1, 3
!c               nvar = nvar + 1
               scale(i) = 12.0d0
!c            end do
!c         end if
      end do
      write (iout,*) "minim_dc Calling lbfgs"
      call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
      deallocate(scale)
#else
!      call check_ecartint
      call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)      
!      call check_ecartint
#endif
      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,1)),i,(c(j,i),j=1,3),&
!          (c(j,i+nres),j=1,3)
!      enddo
!el----------------------------
#ifndef LBFGS
      etot=v(10)                                                      
      iretcode=iv(1)
      nfun=iv(6)
#endif
      return
      end subroutine  minim_dc


#ifdef LBFGS
      real(kind=8) function funcgrad_dc(x,g)
      use energy, only: etotal, zerograd,cartgrad
!      implicit real*8 (a-h,o-z)
#ifdef MPI
      include 'mpif.h'
#endif
      integer k,nf,i,j
      real(kind=8),dimension(nres*6):: x,g
      double precision energia(0:n_ene)
!      common /gacia/ nf
!c
      nf=nf+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(0))
!c      write (iout,*) "energia",energia(0)
      funcgrad_dc=energia(0)
      call cartgrad
      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

      return
      end function
#else
!-----------------------------------------------------------------------------
      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
#endif
!-----------------------------------------------------------------------------
! minim_mcmf.F
!-----------------------------------------------------------------------------
#ifdef MPI
      subroutine minim_mcmf

      use MPI_data
      use csa_data
#ifndef LBFGS
      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 :: idum(1),nf	!
      integer :: lv
      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
#endif
      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,1)
      if (i.eq.10 .or. i.eq.ntyp1 .or. molnum(res_pick).eq.5) 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,1)==10)
        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1) &
        .and.(molnum(i).ne.5)) 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
!      write(iout,*) 'Warning calling chainbuild'
      call chainbuild
!      write(iout,*) "before egb1"
      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,molnum(res_pick))),theta(res_pick+1),&
             alph(res_pick),omeg(res_pick),fail,molnum(res_pick))

!     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)
      write(iout,*) 'Warning calling chainbuild'
      call chainbuild
      call etotal(energy_)
      etot=energy_(0)

      return
      end subroutine sc_minimize
!-----------------------------------------------------------------------------
      subroutine minimize_sc1(etot,x,iretcode,nfun)
#ifndef LBFGS
      use energy, only: func,gradient,fdum,etotal,enerprint
#endif
      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:: lv
      integer :: nvar_restr,i,j
      integer :: idum(1)
      real(kind=8) :: rdum(1),etot	!,fdum
#ifndef LBFGS
      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,1)),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)
#endif
      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
      write(iout,*) 'Warning calling chainbuild'
      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)
      write(iout,*) 'Warning calling chainbuild'
      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,1).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,1).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,1).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
        if ((itype(i,1).eq.ntyp1).or.(molnum(i).gt.1)) cycle
        itypi=iabs(itype(i,1))
        itypi1=iabs(itype(i+1,1))
!         print *,"ebg1",i,itypi,itypi1
        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,1))
            if ((itype(j,1).eq.ntyp1).or.(molnum(j).gt.1)) cycle
!            print *,"ebg1",j,itypj

            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_aq(itypi,itypj)
            e2=fac*bb_aq(itypi,itypj)
            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
            eps2der=evdwij*eps3rt
            eps3der=evdwij*eps2rt
            evdwij=evdwij*eps2rt*eps3rt
            evdw=evdw+evdwij
!            if (wliptran.gt.0.0) print *,"WARNING eps_aq used!"
            if (lprn) then
            sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
            epsi=bb_aq(itypi,itypj)**2/aa_aq(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
      integer::lv
      integer:: iv(liv)
      integer ::  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 ::  iv1, nf
      real(kind=8) :: f
      integer::g1
!
!  ***  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) 
      print *,"my new g1",g1
!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
      use MD_data,  only: itime_mat
!
!  ***  carry out sumsl (unconstrained minimization) iterations, using
!  ***  double-dogleg/bfgs steps.
!
!  ***  parameter declarations  ***
!
      integer :: liv, n
      integer :: lv
      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)
      itime_mat=k
!      imat_update=k
      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 ::  n
      integer :: lv
      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
!-----------------------------------------------------------------------------

!
!
!     ###################################################
!     ##  COPYRIGHT (C)  1999  by  Jay William Ponder  ##
!     ##              All Rights Reserved              ##
!     ###################################################
!
!     ##############################################################
!     ##                                                          ##
!     ##  subroutine lbfgs  --  limited memory BFGS optimization  ##
!     ##                                                          ##
!     ##############################################################
!
!
!     "lbfgs" is a limited memory BFGS quasi-newton nonlinear
!     optimization routine
!
!     literature references:
!
!     J. Nocedal, "Updating Quasi-Newton Matrices with Limited
!     Storage", Mathematics of Computation, 35, 773-782 (1980)
!
!     D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method
!     for Large Scale Optimization", Mathematical Programming,
!     45, 503-528 (1989)
!
!     J. Nocedal and S. J. Wright, "Numerical Optimization",
!     Springer-Verlag, New York, 1999, Section 9.1
!
!     variables and parameters:
!
!     nvar      number of parameters in the objective function
!     x0        contains starting point upon input, upon return
!                 contains the best point found
!     minimum   during optimization contains best current function
!                 value; returns final best function value
!     grdmin    normal exit if rms gradient gets below this value
!     ncalls    total number of function/gradient evaluations
!
!     required external routines:
!
!     fgvalue    function to evaluate function and gradient values
!     optsave    subroutine to write out info about current status
!
!
      subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave)
      use inform
      use iounit
      use keys
      use linmin
      use math2
      use minima
      use output
      use scales
      use MD_data,  only: itime_mat
#ifdef LBFGS
      use energy_data,only:statusbf,niter,ncalls
#endif
      implicit none
      integer i,j,k,m
      integer nvar,next
      integer msav,muse
#ifndef LBFGS
      integer niter,ncalls
      character*9 statusbf
#endif
      integer nerr,maxerr
      real*8 f,f_old,fgvalue
      real*8 f_move,x_move
      real*8 g_norm,g_rms
      real*8 minimum,grdmin
      real*8 angle,rms,beta
      real*8 ys,yy,gamma
      real*8 x0(*)
      real*8, allocatable :: rho(:)
      real*8, allocatable :: alpha(:)
      real*8, allocatable :: x_old(:)
      real*8, allocatable :: g(:)
      real*8, allocatable :: g_old(:)
      real*8, allocatable :: p(:)
      real*8, allocatable :: q(:)
      real*8, allocatable :: r(:)
      real*8, allocatable :: h0(:)
      real*8, allocatable :: s(:,:)
      real*8, allocatable :: y(:,:)
      logical done
      character*9 blank
      character*20 keyword
      character*240 record
      character*240 string
      external fgvalue,optsave
!      common /lbfgstat/ status,niter,ncalls
!c
!c
!c     initialize some values to be used below
!c
      ncalls = 0
      rms = sqrt(dble(nvar))
      if (coordtype .eq. 'CARTESIAN') then
         rms = rms / sqrt(3.0d0)
      else if (coordtype .eq. 'RIGIDBODY') then
         rms = rms / sqrt(6.0d0)
      end if
      blank = '         '
      done = .false.
      nerr = 0
      maxerr = 2
!c
!c     perform dynamic allocation of some global arrays
!c
      if (.not. allocated(scale))  allocate (scale(nvar))
!c
!c     set default values for variable scale factors
!c
      if (.not. set_scale) then
         do i = 1, nvar
            if (scale(i) .eq. 0.0d0)  scale(i) = 1.0d0
         end do
      end if
!c
!c     set default parameters for the optimization
!c
      msav = min(nvar,20)
      if (fctmin .eq. 0.0d0)  fctmin = -100000000.0d0
      if (maxiter .eq. 0)  maxiter = 1000000
      if (nextiter .eq. 0)  nextiter = 1
      if (jprint .lt. 0)  jprint = 1
      if (iwrite .lt. 0)  iwrite = 1
      write(iout,*) "maxiter",maxiter
!c
!c     set default parameters for the line search
!c
      if (stpmax .eq. 0.0d0)  stpmax = 5.0d0
      stpmin = 1.0d-16
      cappa = 0.9d0
      slpmax = 10000.0d0
      angmax = 180.0d0
      intmax = 5
!c
!c     search the keywords for optimization parameters
!c
#ifdef LBFGSREAD
      do i = 1, nkey
         next = 1
         record = keyline(i)
         call gettext (record,keyword,next)
         call upcase (keyword)
         string = record(next:240)
         if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then
            read (string,*,err=10,end=10)  msav
            msav = max(0,min(msav,nvar))
         else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then
            msav = 0
         else if (keyword(1:7) .eq. 'FCTMIN ') then
            read (string,*,err=10,end=10)  fctmin
         else if (keyword(1:8) .eq. 'MAXITER ') then
            read (string,*,err=10,end=10)  maxiter
         else if (keyword(1:8) .eq. 'STEPMAX ') then
            read (string,*,err=10,end=10)  stpmax
         else if (keyword(1:8) .eq. 'STEPMIN ') then
            read (string,*,err=10,end=10)  stpmin
         else if (keyword(1:6) .eq. 'CAPPA ') then
            read (string,*,err=10,end=10)  cappa
         else if (keyword(1:9) .eq. 'SLOPEMAX ') then
            read (string,*,err=10,end=10)  slpmax
         else if (keyword(1:7) .eq. 'ANGMAX ') then
            read (string,*,err=10,end=10)  angmax
         else if (keyword(1:7) .eq. 'INTMAX ') then
            read (string,*,err=10,end=10)  intmax
         end if
   10    continue
      end do
#endif
!c
!c     print header information about the optimization method
!c
      if (jprint .gt. 0) then
         if (msav .eq. 0) then
            write (jout,20)
   20       format (/,' Steepest Descent Gradient Optimization :')
            write (jout,30)
   30       format (/,' SD Iter     F Value      G RMS      F Move',&
                      '   X Move   Angle  FG Call  Comment',/)
         else
            write (jout,40)
   40       format (/,' Limited Memory BFGS Quasi-Newton',&
                      ' Optimization :')
            write (jout,50)
   50       format (/,' QN Iter     F Value      G RMS      F Move',&
                      '   X Move   Angle  FG Call  Comment',/)
         end if
         flush (jout)
      end if
!c
!c     perform dynamic allocation of some local arrays
!c
      allocate (x_old(nvar))
      allocate (g(nvar))
      allocate (g_old(nvar))
      allocate (p(nvar))
      allocate (q(nvar))
      allocate (r(nvar))
      allocate (h0(nvar))
      if (msav .ne. 0) then
         allocate (rho(msav))
         allocate (alpha(msav))
         allocate (s(nvar,msav))
         allocate (y(nvar,msav))
      end if
!c
!c     evaluate the function and get the initial gradient
!c
      niter = nextiter - 1
      maxiter = niter + maxiter
      ncalls = ncalls + 1
!      itime_mat=niter
      f = fgvalue (x0,g)
      f_old = f
      m = 0
      gamma = 1.0d0
      g_norm = 0.0d0
      g_rms = 0.0d0
      do i = 1, nvar
         g_norm = g_norm + g(i)*g(i)
         g_rms = g_rms + (g(i)*scale(i))**2
      end do
      g_norm = sqrt(g_norm)
      g_rms = sqrt(g_rms) / rms
      f_move = 0.5d0 * stpmax * g_norm
!c
!c     print initial information prior to first iteration
!c
      write(jout,*) "before first"
      if (jprint .gt. 0) then
         if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then
            write (jout,60)  niter,f,g_rms,ncalls
   60       format (i6,f14.4,f11.4,2x,i7)
         else
            write (jout,70)  niter,f,g_rms,ncalls
   70       format (i6,d14.4,d11.4,2x,i7)
         end if
         flush (jout)
      end if
!c
!c     write initial intermediate prior to first iteration
!c
      if (iwrite .gt. 0)  call optsave (niter,f,x0)
!c
!c     tests of the various termination criteria
!c
      if (niter .ge. maxiter) then
         statusbf = 'IterLimit'
         done = .true.
      end if
      if (f .le. fctmin) then
         statusbf = 'SmallFct '
         done = .true.
      end if
      if (g_rms .le. grdmin) then
         statusbf = 'SmallGrad'
         done = .true.
      end if
!c
!c     start of a new limited memory BFGS iteration
!c
      do while (.not. done)
         niter = niter + 1
!c         write (jout,*) "LBFGS niter",niter
         muse = min(niter-1,msav)
         m = m + 1
         if (m .gt. msav)  m = 1
!c
!c     estimate Hessian diagonal and compute the Hg product
!c
         do i = 1, nvar
            h0(i) = gamma
            q(i) = g(i)
         end do
         k = m
         do j = 1, muse
            k = k - 1
            if (k .eq. 0)  k = msav
            alpha(k) = 0.0d0
            do i = 1, nvar
               alpha(k) = alpha(k) + s(i,k)*q(i)
            end do
            alpha(k) = alpha(k) * rho(k)
            do i = 1, nvar
               q(i) = q(i) - alpha(k)*y(i,k)
            end do
         end do
         do i = 1, nvar
            r(i) = h0(i) * q(i)
         end do
         do j = 1, muse
            beta = 0.0d0
            do i = 1, nvar
               beta = beta + y(i,k)*r(i)
            end do
            beta = beta * rho(k)
            do i = 1, nvar
               r(i) = r(i) + s(i,k)*(alpha(k)-beta)
            end do
            k = k + 1
            if (k .gt. msav)  k = 1
         end do
!c
!c     set search direction and store current point and gradient
!c
         do i = 1, nvar
            p(i) = -r(i)
            x_old(i) = x0(i)
            g_old(i) = g(i)
         end do
!c
!c     perform line search along the new conjugate direction
!c
         statusbf = blank
!c         write (jout,*) "Before search"
         call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,statusbf)
!c         write (jout,*) "After search"
!c
!c     update variables based on results of this iteration
!c
         if (msav .ne. 0) then
            ys = 0.0d0
            yy = 0.0d0
            do i = 1, nvar
               s(i,m) = x0(i) - x_old(i)
               y(i,m) = g(i) - g_old(i)
               ys = ys + y(i,m)*s(i,m)
               yy = yy + y(i,m)*y(i,m)
            end do
            gamma = abs(ys/yy)
            rho(m) = 1.0d0 / ys
         end if
!c
!c     get the sizes of the moves made during this iteration
!c
         f_move = f_old - f
!         if (f_move.eq.0) f_move=1.0d0
         f_old = f
         x_move = 0.0d0
         do i = 1, nvar
            x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
         end do
         x_move = sqrt(x_move) / rms
         if (coordtype .eq. 'INTERNAL') then
            x_move = radian * x_move
         end if
!c
!c     compute the rms gradient per optimization parameter
!c
         g_rms = 0.0d0
         do i = 1, nvar
            g_rms = g_rms + (g(i)*scale(i))**2
         end do
         g_rms = sqrt(g_rms) / rms
!c
!c     test for error due to line search problems
!c
         if (statusbf.eq.'BadIntpln' .or. statusbf.eq.'IntplnErr') then
            nerr = nerr + 1
            if (nerr .ge. maxerr)  done = .true.
         else
            nerr = 0
         end if
!c
!c     test for too many total iterations
!c
         if (niter .ge. maxiter) then
            statusbf = 'IterLimit'
            done = .true.
         end if
!c
!c     test the normal termination criteria
!c
         if (f .le. fctmin) then
            statusbf = 'SmallFct '
            done = .true.
         end if
         if (g_rms .le. grdmin) then
            statusbf = 'SmallGrad'
            done = .true.
         end if
!c
!c     print intermediate results for the current iteration
!c
!         write(iout,*) "niter1",niter
         itime_mat=niter
         if (jprint .gt. 0) then
            if (done .or. mod(niter,jprint).eq.0) then
               if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.&
                  g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.&
                  f_move.gt.-1.0d5) then 
                  write (jout,80)  niter,f,g_rms,f_move,x_move,&
                                  angle,ncalls,statusbf
   80             format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9)
               else
                  write (jout,90)  niter,f,g_rms,f_move,x_move,&
                                  angle,ncalls,statusbf
   90             format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9)
               end if
            end if
            flush (jout)
         end if
!c
!c     write intermediate results for the current iteration
!c
         if (iwrite .gt. 0) then
            if (done .or. mod(niter,iwrite).eq.0) then
               call optsave (niter,f,x0)
            end if
         end if
      end do
!c
!c     perform deallocation of some local arrays
!c
      deallocate (x_old)
      deallocate (g)
      deallocate (g_old)
      deallocate (p)
      deallocate (q)
      deallocate (r)
      deallocate (h0)
      if (msav .ne. 0) then
         deallocate (rho)
         deallocate (alpha)
         deallocate (s)
         deallocate (y)
      end if
!c
!c     set final value of the objective function
!c
      minimum = f
      if (jprint .gt. 0) then
         if (statusbf.eq.'SmallGrad' .or. statusbf.eq.'SmallFct ') then
            write (jout,100)  statusbf
  100       format (/,' LBFGS  --  Normal Termination due to ',a9)
         else
            write (jout,110)  statusbf
  110       format (/,' LBFGS  --  Incomplete Convergence due to ',a9)
         end if
         flush (jout)
      end if
      return
      end subroutine
!ifdef LBFGS
!     double precision function funcgrad(x,g)
!     implicit none
!     include 'DIMENSIONS'
!     include 'COMMON.CONTROL'
!     include 'COMMON.CHAIN'
!     include 'COMMON.DERIV'
!     include 'COMMON.VAR'
!     include 'COMMON.INTERACT'
!     include 'COMMON.FFIELD'
!     include 'COMMON.MD'
!     include 'COMMON.QRESTR'
!     include 'COMMON.IOUNITS'
!     include 'COMMON.GEO'
!     double precision energia(0:n_ene)
!     double precision x(nvar),g(nvar)
!     integer i
!     if (jjj.gt.0) then
!      write (iout,*) "in func x"
!      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
!     endif
!     call var_to_geom(nvar,x)
!     call zerograd
!     call chainbuild_extconf
!     call etotal(energia(0))
!     call sum_gradient
!     funcgrad=energia(0)
!     call cart2intgrad(nvar,g)
!
! Add the components corresponding to local energy terms.
!
! Add the usampl contributions
!     if (usampl) then
!        do i=1,nres-3
!          gloc(i,icg)=gloc(i,icg)+dugamma(i)
!        enddo
!        do i=1,nres-2
!          gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
!        enddo
!     endif
!     do i=1,nvar
!d      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
!       g(i)=g(i)+gloc(i,icg)
!     enddo
!     return
!     end
!endif
!-----------------------------------------------------------------------------
      end module minimm