source/unres/src_CSA/cored.f

   1       subroutine assst(iv, liv, lv, v)
   2 c
   3 c  ***  assess candidate step (***sol version 2.3)  ***
   4 c
   5       integer liv, l
   6       integer iv(liv)
   7       double precision v(lv)
   8 c
   9 c  ***  purpose  ***
  10 c
  11 c        this subroutine is called by an unconstrained minimization
  12 c     routine to assess the next candidate step.  it may recommend one
  13 c     of several courses of action, such as accepting the step, recom-
  14 c     puting it using the same or a new quadratic model, or halting due
  15 c     to convergence or false convergence.  see the return code listing
  16 c     below.
  17 c
  18 c--------------------------  parameter usage  --------------------------
  19 c
  20 c  iv (i/o) integer parameter and scratch vector -- see description
  21 c             below of iv values referenced.
  22 c liv (in)  length of iv array.
  23 c  lv (in)  length of v array.
  24 c   v (i/o) real parameter and scratch vector -- see description
  25 c             below of v values referenced.
  26 c
  27 c  ***  iv values referenced  ***
  28 c
  29 c    iv(irc) (i/o) on input for the first step tried in a new iteration,
  30 c             iv(irc) should be set to 3 or 4 (the value to which it is
  31 c             set when step is definitely to be accepted).  on input
  32 c             after step has been recomputed, iv(irc) should be
  33 c             unchanged since the previous return of assst.
  34 c                on output, iv(irc) is a return code having one of the
  35 c             following values...
  36 c                  1 = switch models or try smaller step.
  37 c                  2 = switch models or accept step.
  38 c                  3 = accept step and determine v(radfac) by gradient
  39 c                       tests.
  40 c                  4 = accept step, v(radfac) has been determined.
  41 c                  5 = recompute step (using the same model).
  42 c                  6 = recompute step with radius = v(lmaxs) but do not
  43 c                       evaulate the objective function.
  44 c                  7 = x-convergence (see v(xctol)).
  45 c                  8 = relative function convergence (see v(rfctol)).
  46 c                  9 = both x- and relative function convergence.
  47 c                 10 = absolute function convergence (see v(afctol)).
  48 c                 11 = singular convergence (see v(lmaxs)).
  49 c                 12 = false convergence (see v(xftol)).
  50 c                 13 = iv(irc) was out of range on input.
  51 c             return code i has precdence over i+1 for i = 9, 10, 11.
  52 c iv(mlstgd) (i/o) saved value of iv(model).
  53 c  iv(model) (i/o) on input, iv(model) should be an integer identifying
  54 c             the current quadratic model of the objective function.
  55 c             if a previous step yielded a better function reduction,
  56 c             then iv(model) will be set to iv(mlstgd) on output.
  57 c iv(nfcall) (in)  invocation count for the objective function.
  58 c iv(nfgcal) (i/o) value of iv(nfcall) at step that gave the biggest
  59 c             function reduction this iteration.  iv(nfgcal) remains
  60 c             unchanged until a function reduction is obtained.
  61 c iv(radinc) (i/o) the number of radius increases (or minus the number
  62 c             of decreases) so far this iteration.
  63 c iv(restor) (out) set to 1 if v(f) has been restored and x should be
  64 c             restored to its initial value, to 2 if x should be saved,
  65 c             to 3 if x should be restored from the saved value, and to
  66 c             0 otherwise.
  67 c  iv(stage) (i/o) count of the number of models tried so far in the
  68 c             current iteration.
  69 c iv(stglim) (in)  maximum number of models to consider.
  70 c iv(switch) (out) set to 0 unless a new model is being tried and it
  71 c             gives a smaller function value than the previous model,
  72 c             in which case assst sets iv(switch) = 1.
  73 c iv(toobig) (in)  is nonzero if step was too big (e.g. if it caused
  74 c             overflow).
  75 c   iv(xirc) (i/o) value that iv(irc) would have in the absence of
  76 c             convergence, false convergence, and oversized steps.
  77 c
  78 c  ***  v values referenced  ***
  79 c
  80 c v(afctol) (in)  absolute function convergence tolerance.  if the
  81 c             absolute value of the current function value v(f) is less
  82 c             than v(afctol), then assst returns with iv(irc) = 10.
  83 c v(decfac) (in)  factor by which to decrease radius when iv(toobig) is
  84 c             nonzero.
  85 c v(dstnrm) (in)  the 2-norm of d*step.
  86 c v(dstsav) (i/o) value of v(dstnrm) on saved step.
  87 c   v(dst0) (in)  the 2-norm of d times the newton step (when defined,
  88 c             i.e., for v(nreduc) .ge. 0).
  89 c      v(f) (i/o) on both input and output, v(f) is the objective func-
  90 c             tion value at x.  if x is restored to a previous value,
  91 c             then v(f) is restored to the corresponding value.
  92 c   v(fdif) (out) the function reduction v(f0) - v(f) (for the output
  93 c             value of v(f) if an earlier step gave a bigger function
  94 c             decrease, and for the input value of v(f) otherwise).
  95 c v(flstgd) (i/o) saved value of v(f).
  96 c     v(f0) (in)  objective function value at start of iteration.
  97 c v(gtslst) (i/o) value of v(gtstep) on saved step.
  98 c v(gtstep) (in)  inner product between step and gradient.
  99 c v(incfac) (in)  minimum factor by which to increase radius.
 100 c  v(lmaxs) (in)  maximum reasonable step size (and initial step bound).
 101 c             if the actual function decrease is no more than twice
 102 c             what was predicted, if a return with iv(irc) = 7, 8, 9,
 103 c             or 10 does not occur, if v(dstnrm) .gt. v(lmaxs), and if
 104 c             v(preduc) .le. v(sctol) * abs(v(f0)), then assst re-
 105 c             turns with iv(irc) = 11.  if so doing appears worthwhile,
 106 c             then assst repeats this test with v(preduc) computed for
 107 c             a step of length v(lmaxs) (by a return with iv(irc) = 6).
 108 c v(nreduc) (i/o)  function reduction predicted by quadratic model for
 109 c             newton step.  if assst is called with iv(irc) = 6, i.e.,
 110 c             if v(preduc) has been computed with radius = v(lmaxs) for
 111 c             use in the singular convervence test, then v(nreduc) is
 112 c             set to -v(preduc) before the latter is restored.
 113 c v(plstgd) (i/o) value of v(preduc) on saved step.
 114 c v(preduc) (i/o) function reduction predicted by quadratic model for
 115 c             current step.
 116 c v(radfac) (out) factor to be used in determining the new radius,
 117 c             which should be v(radfac)*dst, where  dst  is either the
 118 c             output value of v(dstnrm) or the 2-norm of
 119 c             diag(newd)*step  for the output value of step and the
 120 c             updated version, newd, of the scale vector d.  for
 121 c             iv(irc) = 3, v(radfac) = 1.0 is returned.
 122 c v(rdfcmn) (in)  minimum value for v(radfac) in terms of the input
 123 c             value of v(dstnrm) -- suggested value = 0.1.
 124 c v(rdfcmx) (in)  maximum value for v(radfac) -- suggested value = 4.0.
 125 c  v(reldx) (in) scaled relative change in x caused by step, computed
 126 c             (e.g.) by function  reldst  as
 127 c                 max (d(i)*abs(x(i)-x0(i)), 1 .le. i .le. p) /
 128 c                    max (d(i)*(abs(x(i))+abs(x0(i))), 1 .le. i .le. p).
 129 c v(rfctol) (in)  relative function convergence tolerance.  if the
 130 c             actual function reduction is at most twice what was pre-
 131 c             dicted and  v(nreduc) .le. v(rfctol)*abs(v(f0)),  then
 132 c             assst returns with iv(irc) = 8 or 9.
 133 c v(stppar) (in)  marquardt parameter -- 0 means full newton step.
 134 c v(tuner1) (in)  tuning constant used to decide if the function
 135 c             reduction was much less than expected.  suggested
 136 c             value = 0.1.
 137 c v(tuner2) (in)  tuning constant used to decide if the function
 138 c             reduction was large enough to accept step.  suggested
 139 c             value = 10**-4.
 140 c v(tuner3) (in)  tuning constant used to decide if the radius
 141 c             should be increased.  suggested value = 0.75.
 142 c  v(xctol) (in)  x-convergence criterion.  if step is a newton step
 143 c             (v(stppar) = 0) having v(reldx) .le. v(xctol) and giving
 144 c             at most twice the predicted function decrease, then
 145 c             assst returns iv(irc) = 7 or 9.
 146 c  v(xftol) (in)  false convergence tolerance.  if step gave no or only
 147 c             a small function decrease and v(reldx) .le. v(xftol),
 148 c             then assst returns with iv(irc) = 12.
 149 c
 150 c-------------------------------  notes  -------------------------------
 151 c
 152 c  ***  application and usage restrictions  ***
 153 c
 154 c        this routine is called as part of the nl2sol (nonlinear
 155 c     least-squares) package.  it may be used in any unconstrained
 156 c     minimization solver that uses dogleg, goldfeld-quandt-trotter,
 157 c     or levenberg-marquardt steps.
 158 c
 159 c  ***  algorithm notes  ***
 160 c
 161 c        see (1) for further discussion of the assessing and model
 162 c     switching strategies.  while nl2sol considers only two models,
 163 c     assst is designed to handle any number of models.
 164 c
 165 c  ***  usage notes  ***
 166 c
 167 c        on the first call of an iteration, only the i/o variables
 168 c     step, x, iv(irc), iv(model), v(f), v(dstnrm), v(gtstep), and
 169 c     v(preduc) need have been initialized.  between calls, no i/o
 170 c     values execpt step, x, iv(model), v(f) and the stopping toler-
 171 c     ances should be changed.
 172 c        after a return for convergence or false convergence, one can
 173 c     change the stopping tolerances and call assst again, in which
 174 c     case the stopping tests will be repeated.
 175 c
 176 c  ***  references  ***
 177 c
 178 c     (1) dennis, j.e., jr., gay, d.m., and welsch, r.e. (1981),
 179 c        an adaptive nonlinear least-squares algorithm,
 180 c        acm trans. math. software, vol. 7, no. 3.
 181 c
 182 c     (2) powell, m.j.d. (1970)  a fortran subroutine for solving
 183 c        systems of nonlinear algebraic equations, in numerical
 184 c        methods for nonlinear algebraic equations, edited by
 185 c        p. rabinowitz, gordon and breach, london.
 186 c
 187 c  ***  history  ***
 188 c
 189 c        john dennis designed much of this routine, starting with
 190 c     ideas in (2). roy welsch suggested the model switching strategy.
 191 c        david gay and stephen peters cast this subroutine into a more
 192 c     portable form (winter 1977), and david gay cast it into its
 193 c     present form (fall 1978).
 194 c
 195 c  ***  general  ***
 196 c
 197 c     this subroutine was written in connection with research
 198 c     supported by the national science foundation under grants
 199 c     mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
 200 c     mcs-7906671.
 201 c
 202 c------------------------  external quantities  ------------------------
 203 c
 204 c  ***  no external functions and subroutines  ***
 205 c
 206 c  ***  intrinsic functions  ***
 207 c/+
 208       double precision dabs, dmax1
 209 c/
 210 c  ***  no common blocks  ***
 211 c
 212 c--------------------------  local variables  --------------------------
 213 c
 214       logical goodx
 215       integer i, nfc
 216       double precision emax, emaxs, gts, rfac1, xmax
 217       double precision half, one, onep2, two, zero
 218 c
 219 c  ***  subscripts for iv and v  ***
 220 c
 221       integer afctol, decfac, dstnrm, dstsav, dst0, f, fdif, flstgd, f0,
 222      1        gtslst, gtstep, incfac, irc, lmaxs, mlstgd, model, nfcall,
 223      2        nfgcal, nreduc, plstgd, preduc, radfac, radinc, rdfcmn,
 224      3        rdfcmx, reldx, restor, rfctol, sctol, stage, stglim,
 225      4        stppar, switch, toobig, tuner1, tuner2, tuner3, xctol,
 226      5        xftol, xirc
 227 c
 228 c  ***  data initializations  ***
 229 c
 230 c/6
 231 c     data half/0.5d+0/, one/1.d+0/, onep2/1.2d+0/, two/2.d+0/,
 232 c    1     zero/0.d+0/
 233 c/7
 234       parameter (half=0.5d+0, one=1.d+0, onep2=1.2d+0, two=2.d+0,
 235      1           zero=0.d+0)
 236 c/
 237 c
 238 c/6
 239 c     data irc/29/, mlstgd/32/, model/5/, nfcall/6/, nfgcal/7/,
 240 c    1     radinc/8/, restor/9/, stage/10/, stglim/11/, switch/12/,
 241 c    2     toobig/2/, xirc/13/
 242 c/7
 243       parameter (irc=29, mlstgd=32, model=5, nfcall=6, nfgcal=7,
 244      1           radinc=8, restor=9, stage=10, stglim=11, switch=12,
 245      2           toobig=2, xirc=13)
 246 c/
 247 c/6
 248 c     data afctol/31/, decfac/22/, dstnrm/2/, dst0/3/, dstsav/18/,
 249 c    1     f/10/, fdif/11/, flstgd/12/, f0/13/, gtslst/14/, gtstep/4/,
 250 c    2     incfac/23/, lmaxs/36/, nreduc/6/, plstgd/15/, preduc/7/,
 251 c    3     radfac/16/, rdfcmn/24/, rdfcmx/25/, reldx/17/, rfctol/32/,
 252 c    4     sctol/37/, stppar/5/, tuner1/26/, tuner2/27/, tuner3/28/,
 253 c    5     xctol/33/, xftol/34/
 254 c/7
 255       parameter (afctol=31, decfac=22, dstnrm=2, dst0=3, dstsav=18,
 256      1           f=10, fdif=11, flstgd=12, f0=13, gtslst=14, gtstep=4,
 257      2           incfac=23, lmaxs=36, nreduc=6, plstgd=15, preduc=7,
 258      3           radfac=16, rdfcmn=24, rdfcmx=25, reldx=17, rfctol=32,
 259      4           sctol=37, stppar=5, tuner1=26, tuner2=27, tuner3=28,
 260      5           xctol=33, xftol=34)
 261 c/
 262 c
 263 c+++++++++++++++++++++++++++++++  body  ++++++++++++++++++++++++++++++++
 264 c
 265       nfc = iv(nfcall)
 266       iv(switch) = 0
 267       iv(restor) = 0
 268       rfac1 = one
 269       goodx = .true.
 270       i = iv(irc)
 271       if (i .ge. 1 .and. i .le. 12)
 272      1             go to (20,30,10,10,40,280,220,220,220,220,220,170), i
 273          iv(irc) = 13
 274          go to 999
 275 c
 276 c  ***  initialize for new iteration  ***
 277 c
 278  10   iv(stage) = 1
 279       iv(radinc) = 0
 280       v(flstgd) = v(f0)
 281       if (iv(toobig) .eq. 0) go to 110
 282          iv(stage) = -1
 283          iv(xirc) = i
 284          go to 60
 285 c
 286 c  ***  step was recomputed with new model or smaller radius  ***
 287 c  ***  first decide which  ***
 288 c
 289  20   if (iv(model) .ne. iv(mlstgd)) go to 30
 290 c        ***  old model retained, smaller radius tried  ***
 291 c        ***  do not consider any more new models this iteration  ***
 292          iv(stage) = iv(stglim)
 293          iv(radinc) = -1
 294          go to 110
 295 c
 296 c  ***  a new model is being tried.  decide whether to keep it.  ***
 297 c
 298  30   iv(stage) = iv(stage) + 1
 299 c
 300 c     ***  now we add the possibiltiy that step was recomputed with  ***
 301 c     ***  the same model, perhaps because of an oversized step.     ***
 302 c
 303  40   if (iv(stage) .gt. 0) go to 50
 304 c
 305 c        ***  step was recomputed because it was too big.  ***
 306 c
 307          if (iv(toobig) .ne. 0) go to 60
 308 c
 309 c        ***  restore iv(stage) and pick up where we left off.  ***
 310 c
 311          iv(stage) = -iv(stage)
 312          i = iv(xirc)
 313          go to (20, 30, 110, 110, 70), i
 314 c
 315  50   if (iv(toobig) .eq. 0) go to 70
 316 c
 317 c  ***  handle oversize step  ***
 318 c
 319       if (iv(radinc) .gt. 0) go to 80
 320          iv(stage) = -iv(stage)
 321          iv(xirc) = iv(irc)
 322 c
 323  60      v(radfac) = v(decfac)
 324          iv(radinc) = iv(radinc) - 1
 325          iv(irc) = 5
 326          iv(restor) = 1
 327          go to 999
 328 c
 329  70   if (v(f) .lt. v(flstgd)) go to 110
 330 c
 331 c     *** the new step is a loser.  restore old model.  ***
 332 c
 333       if (iv(model) .eq. iv(mlstgd)) go to 80
 334          iv(model) = iv(mlstgd)
 335          iv(switch) = 1
 336 c
 337 c     ***  restore step, etc. only if a previous step decreased v(f).
 338 c
 339  80   if (v(flstgd) .ge. v(f0)) go to 110
 340          iv(restor) = 1
 341          v(f) = v(flstgd)
 342          v(preduc) = v(plstgd)
 343          v(gtstep) = v(gtslst)
 344          if (iv(switch) .eq. 0) rfac1 = v(dstnrm) / v(dstsav)
 345          v(dstnrm) = v(dstsav)
 346          nfc = iv(nfgcal)
 347          goodx = .false.
 348 c
 349  110  v(fdif) = v(f0) - v(f)
 350       if (v(fdif) .gt. v(tuner2) * v(preduc)) go to 140
 351       if(iv(radinc).gt.0) go to 140
 352 c
 353 c        ***  no (or only a trivial) function decrease
 354 c        ***  -- so try new model or smaller radius
 355 c
 356          if (v(f) .lt. v(f0)) go to 120
 357               iv(mlstgd) = iv(model)
 358               v(flstgd) = v(f)
 359               v(f) = v(f0)
 360               iv(restor) = 1
 361               go to 130
 362  120     iv(nfgcal) = nfc
 363  130     iv(irc) = 1
 364          if (iv(stage) .lt. iv(stglim)) go to 160
 365               iv(irc) = 5
 366               iv(radinc) = iv(radinc) - 1
 367               go to 160
 368 c
 369 c  ***  nontrivial function decrease achieved  ***
 370 c
 371  140  iv(nfgcal) = nfc
 372       rfac1 = one
 373       v(dstsav) = v(dstnrm)
 374       if (v(fdif) .gt. v(preduc)*v(tuner1)) go to 190
 375 c
 376 c  ***  decrease was much less than predicted -- either change models
 377 c  ***  or accept step with decreased radius.
 378 c
 379       if (iv(stage) .ge. iv(stglim)) go to 150
 380 c        ***  consider switching models  ***
 381          iv(irc) = 2
 382          go to 160
 383 c
 384 c     ***  accept step with decreased radius  ***
 385 c
 386  150  iv(irc) = 4
 387 c
 388 c  ***  set v(radfac) to fletcher*s decrease factor  ***
 389 c
 390  160  iv(xirc) = iv(irc)
 391       emax = v(gtstep) + v(fdif)
 392       v(radfac) = half * rfac1
 393       if (emax .lt. v(gtstep)) v(radfac) = rfac1 * dmax1(v(rdfcmn),
 394      1                                           half * v(gtstep)/emax)
 395 c
 396 c  ***  do false convergence test  ***
 397 c
 398  170  if (v(reldx) .le. v(xftol)) go to 180
 399          iv(irc) = iv(xirc)
 400          if (v(f) .lt. v(f0)) go to 200
 401               go to 230
 402 c
 403  180  iv(irc) = 12
 404       go to 240
 405 c
 406 c  ***  handle good function decrease  ***
 407 c
 408  190  if (v(fdif) .lt. (-v(tuner3) * v(gtstep))) go to 210
 409 c
 410 c     ***  increasing radius looks worthwhile.  see if we just
 411 c     ***  recomputed step with a decreased radius or restored step
 412 c     ***  after recomputing it with a larger radius.
 413 c
 414       if (iv(radinc) .lt. 0) go to 210
 415       if (iv(restor) .eq. 1) go to 210
 416 c
 417 c        ***  we did not.  try a longer step unless this was a newton
 418 c        ***  step.
 419 c
 420          v(radfac) = v(rdfcmx)
 421          gts = v(gtstep)
 422          if (v(fdif) .lt. (half/v(radfac) - one) * gts)
 423      1            v(radfac) = dmax1(v(incfac), half*gts/(gts + v(fdif)))
 424          iv(irc) = 4
 425          if (v(stppar) .eq. zero) go to 230
 426          if (v(dst0) .ge. zero .and. (v(dst0) .lt. two*v(dstnrm)
 427      1             .or. v(nreduc) .lt. onep2*v(fdif)))  go to 230
 428 c             ***  step was not a newton step.  recompute it with
 429 c             ***  a larger radius.
 430               iv(irc) = 5
 431               iv(radinc) = iv(radinc) + 1
 432 c
 433 c  ***  save values corresponding to good step  ***
 434 c
 435  200  v(flstgd) = v(f)
 436       iv(mlstgd) = iv(model)
 437       if (iv(restor) .ne. 1) iv(restor) = 2
 438       v(dstsav) = v(dstnrm)
 439       iv(nfgcal) = nfc
 440       v(plstgd) = v(preduc)
 441       v(gtslst) = v(gtstep)
 442       go to 230
 443 c
 444 c  ***  accept step with radius unchanged  ***
 445 c
 446  210  v(radfac) = one
 447       iv(irc) = 3
 448       go to 230
 449 c
 450 c  ***  come here for a restart after convergence  ***
 451 c
 452  220  iv(irc) = iv(xirc)
 453       if (v(dstsav) .ge. zero) go to 240
 454          iv(irc) = 12
 455          go to 240
 456 c
 457 c  ***  perform convergence tests  ***
 458 c
 459  230  iv(xirc) = iv(irc)
 460  240  if (iv(restor) .eq. 1 .and. v(flstgd) .lt. v(f0)) iv(restor) = 3
 461       if (half * v(fdif) .gt. v(preduc)) go to 999
 462       emax = v(rfctol) * dabs(v(f0))
 463       emaxs = v(sctol) * dabs(v(f0))
 464       if (v(dstnrm) .gt. v(lmaxs) .and. v(preduc) .le. emaxs)
 465      1                       iv(irc) = 11
 466       if (v(dst0) .lt. zero) go to 250
 467       i = 0
 468       if ((v(nreduc) .gt. zero .and. v(nreduc) .le. emax) .or.
 469      1    (v(nreduc) .eq. zero. and. v(preduc) .eq. zero))  i = 2
 470       if (v(stppar) .eq. zero .and. v(reldx) .le. v(xctol)
 471      1                        .and. goodx)                  i = i + 1
 472       if (i .gt. 0) iv(irc) = i + 6
 473 c
 474 c  ***  consider recomputing step of length v(lmaxs) for singular
 475 c  ***  convergence test.
 476 c
 477  250  if (iv(irc) .gt. 5 .and. iv(irc) .ne. 12) go to 999
 478       if (v(dstnrm) .gt. v(lmaxs)) go to 260
 479          if (v(preduc) .ge. emaxs) go to 999
 480               if (v(dst0) .le. zero) go to 270
 481                    if (half * v(dst0) .le. v(lmaxs)) go to 999
 482                         go to 270
 483  260  if (half * v(dstnrm) .le. v(lmaxs)) go to 999
 484       xmax = v(lmaxs) / v(dstnrm)
 485       if (xmax * (two - xmax) * v(preduc) .ge. emaxs) go to 999
 486  270  if (v(nreduc) .lt. zero) go to 290
 487 c
 488 c  ***  recompute v(preduc) for use in singular convergence test  ***
 489 c
 490       v(gtslst) = v(gtstep)
 491       v(dstsav) = v(dstnrm)
 492       if (iv(irc) .eq. 12) v(dstsav) = -v(dstsav)
 493       v(plstgd) = v(preduc)
 494       i = iv(restor)
 495       iv(restor) = 2
 496       if (i .eq. 3) iv(restor) = 0
 497       iv(irc) = 6
 498       go to 999
 499 c
 500 c  ***  perform singular convergence test with recomputed v(preduc)  ***
 501 c
 502  280  v(gtstep) = v(gtslst)
 503       v(dstnrm) = dabs(v(dstsav))
 504       iv(irc) = iv(xirc)
 505       if (v(dstsav) .le. zero) iv(irc) = 12
 506       v(nreduc) = -v(preduc)
 507       v(preduc) = v(plstgd)
 508       iv(restor) = 3
 509  290  if (-v(nreduc) .le. v(sctol) * dabs(v(f0))) iv(irc) = 11
 510 c
 511  999  return
 512 c
 513 c  ***  last card of assst follows  ***
 514       end
 515       subroutine deflt(alg, iv, liv, lv, v)
 516 c
 517 c  ***  supply ***sol (version 2.3) default values to iv and v  ***
 518 c
 519 c  ***  alg = 1 means regression constants.
 520 c  ***  alg = 2 means general unconstrained optimization constants.
 521 c
 522       integer liv, l
 523       integer alg, iv(liv)
 524       double precision v(lv)
 525 c
 526       external imdcon, vdflt
 527       integer imdcon
 528 c imdcon... returns machine-dependent integer constants.
 529 c vdflt.... provides default values to v.
 530 c
 531       integer miv, m
 532       integer miniv(2), minv(2)
 533 c
 534 c  ***  subscripts for iv  ***
 535 c
 536       integer algsav, covprt, covreq, dtype, hc, ierr, inith, inits,
 537      1        ipivot, ivneed, lastiv, lastv, lmat, mxfcal, mxiter,
 538      2        nfcov, ngcov, nvdflt, outlev, parprt, parsav, perm,
 539      3        prunit, qrtyp, rdreq, rmat, solprt, statpr, vneed,
 540      4        vsave, x0prt
 541 c
 542 c  ***  iv subscript values  ***
 543 c
 544 c/6
 545 c     data algsav/51/, covprt/14/, covreq/15/, dtype/16/, hc/71/,
 546 c    1     ierr/75/, inith/25/, inits/25/, ipivot/76/, ivneed/3/,
 547 c    2     lastiv/44/, lastv/45/, lmat/42/, mxfcal/17/, mxiter/18/,
 548 c    3     nfcov/52/, ngcov/53/, nvdflt/50/, outlev/19/, parprt/20/,
 549 c    4     parsav/49/, perm/58/, prunit/21/, qrtyp/80/, rdreq/57/,
 550 c    5     rmat/78/, solprt/22/, statpr/23/, vneed/4/, vsave/60/,
 551 c    6     x0prt/24/
 552 c/7
 553       parameter (algsav=51, covprt=14, covreq=15, dtype=16, hc=71,
 554      1           ierr=75, inith=25, inits=25, ipivot=76, ivneed=3,
 555      2           lastiv=44, lastv=45, lmat=42, mxfcal=17, mxiter=18,
 556      3           nfcov=52, ngcov=53, nvdflt=50, outlev=19, parprt=20,
 557      4           parsav=49, perm=58, prunit=21, qrtyp=80, rdreq=57,
 558      5           rmat=78, solprt=22, statpr=23, vneed=4, vsave=60,
 559      6           x0prt=24)
 560 c/
 561       data miniv(1)/80/, miniv(2)/59/, minv(1)/98/, minv(2)/71/
 562 c
 563 c-------------------------------  body  --------------------------------
 564 c
 565       if (alg .lt. 1 .or. alg .gt. 2) go to 40
 566       miv = miniv(alg)
 567       if (liv .lt. miv) go to 20
 568       mv = minv(alg)
 569       if (lv .lt. mv) go to 30
 570       call vdflt(alg, lv, v)
 571       iv(1) = 12
 572       iv(algsav) = alg
 573       iv(ivneed) = 0
 574       iv(lastiv) = miv
 575       iv(lastv) = mv
 576       iv(lmat) = mv + 1
 577       iv(mxfcal) = 200
 578       iv(mxiter) = 150
 579       iv(outlev) = 1
 580       iv(parprt) = 1
 581       iv(perm) = miv + 1
 582       iv(prunit) = imdcon(1)
 583       iv(solprt) = 1
 584       iv(statpr) = 1
 585       iv(vneed) = 0
 586       iv(x0prt) = 1
 587 c
 588       if (alg .ge. 2) go to 10
 589 c
 590 c  ***  regression  values
 591 c
 592       iv(covprt) = 3
 593       iv(covreq) = 1
 594       iv(dtype) = 1
 595       iv(hc) = 0
 596       iv(ierr) = 0
 597       iv(inits) = 0
 598       iv(ipivot) = 0
 599       iv(nvdflt) = 32
 600       iv(parsav) = 67
 601       iv(qrtyp) = 1
 602       iv(rdreq) = 3
 603       iv(rmat) = 0
 604       iv(vsave) = 58
 605       go to 999
 606 c
 607 c  ***  general optimization values
 608 c
 609  10   iv(dtype) = 0
 610       iv(inith) = 1
 611       iv(nfcov) = 0
 612       iv(ngcov) = 0
 613       iv(nvdflt) = 25
 614       iv(parsav) = 47
 615       go to 999
 616 c
 617  20   iv(1) = 15
 618       go to 999
 619 c
 620  30   iv(1) = 16
 621       go to 999
 622 c
 623  40   iv(1) = 67
 624 c
 625  999  return
 626 c  ***  last card of deflt follows  ***
 627       end
 628       double precision function dotprd(p, x, y)
 629 c
 630 c  ***  return the inner product of the p-vectors x and y.  ***
 631 c
 632       integer p
 633       double precision x(p), y(p)
 634 c
 635       integer i
 636       double precision one, sqteta, t, zero
 637 c/+
 638       double precision dmax1, dabs
 639 c/
 640       external rmdcon
 641       double precision rmdcon
 642 c
 643 c  ***  rmdcon(2) returns a machine-dependent constant, sqteta, which
 644 c  ***  is slightly larger than the smallest positive number that
 645 c  ***  can be squared without underflowing.
 646 c
 647 c/6
 648 c     data one/1.d+0/, sqteta/0.d+0/, zero/0.d+0/
 649 c/7
 650       parameter (one=1.d+0, zero=0.d+0)
 651       data sqteta/0.d+0/
 652 c/
 653 c
 654       dotprd = zero
 655       if (p .le. 0) go to 999
 656 crc      if (sqteta .eq. zero) sqteta = rmdcon(2)
 657       do 20 i = 1, p
 658 crc         t = dmax1(dabs(x(i)), dabs(y(i)))
 659 crc         if (t .gt. one) go to 10
 660 crc         if (t .lt. sqteta) go to 20
 661 crc         t = (x(i)/sqteta)*y(i)
 662 crc         if (dabs(t) .lt. sqteta) go to 20
 663  10      dotprd = dotprd + x(i)*y(i)
 664  20   continue
 665 c
 666  999  return
 667 c  ***  last card of dotprd follows  ***
 668       end
 669       subroutine itsum(d, g, iv, liv, lv, p, v, x)
 670 c
 671 c  ***  print iteration summary for ***sol (version 2.3)  ***
 672 c
 673 c  ***  parameter declarations  ***
 674 c
 675       integer liv, lv, p
 676       integer iv(liv)
 677       double precision d(p), g(p), v(lv), x(p)
 678 c
 679 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 680 c
 681 c  ***  local variables  ***
 682 c
 683       integer alg, i, iv1, m, nf, ng, ol, pu
 684 c/6
 685 c     real model1(6), model2(6)
 686 c/7
 687       character*4 model1(6), model2(6)
 688 c/
 689       double precision nreldf, oldf, preldf, reldf, zero
 690 c
 691 c  ***  intrinsic functions  ***
 692 c/+
 693       integer iabs
 694       double precision dabs, dmax1
 695 c/
 696 c  ***  no external functions or subroutines  ***
 697 c
 698 c  ***  subscripts for iv and v  ***
 699 c
 700       integer algsav, dstnrm, f, fdif, f0, needhd, nfcall, nfcov, ngcov,
 701      1        ngcall, niter, nreduc, outlev, preduc, prntit, prunit,
 702      2        reldx, solprt, statpr, stppar, sused, x0prt
 703 c
 704 c  ***  iv subscript values  ***
 705 c
 706 c/6
 707 c     data algsav/51/, needhd/36/, nfcall/6/, nfcov/52/, ngcall/30/,
 708 c    1     ngcov/53/, niter/31/, outlev/19/, prntit/39/, prunit/21/,
 709 c    2     solprt/22/, statpr/23/, sused/64/, x0prt/24/
 710 c/7
 711       parameter (algsav=51, needhd=36, nfcall=6, nfcov=52, ngcall=30,
 712      1           ngcov=53, niter=31, outlev=19, prntit=39, prunit=21,
 713      2           solprt=22, statpr=23, sused=64, x0prt=24)
 714 c/
 715 c
 716 c  ***  v subscript values  ***
 717 c
 718 c/6
 719 c     data dstnrm/2/, f/10/, f0/13/, fdif/11/, nreduc/6/, preduc/7/,
 720 c    1     reldx/17/, stppar/5/
 721 c/7
 722       parameter (dstnrm=2, f=10, f0=13, fdif=11, nreduc=6, preduc=7,
 723      1           reldx=17, stppar=5)
 724 c/
 725 c
 726 c/6
 727 c     data zero/0.d+0/
 728 c/7
 729       parameter (zero=0.d+0)
 730 c/
 731 c/6
 732 c     data model1(1)/4h    /, model1(2)/4h    /, model1(3)/4h    /,
 733 c    1     model1(4)/4h    /, model1(5)/4h  g /, model1(6)/4h  s /,
 734 c    2     model2(1)/4h g  /, model2(2)/4h s  /, model2(3)/4hg-s /,
 735 c    3     model2(4)/4hs-g /, model2(5)/4h-s-g/, model2(6)/4h-g-s/
 736 c/7
 737       data model1/'    ','    ','    ','    ','  g ','  s '/,
 738      1     model2/' g  ',' s  ','g-s ','s-g ','-s-g','-g-s'/
 739 c/
 740 c
 741 c-------------------------------  body  --------------------------------
 742 c
 743       pu = iv(prunit)
 744       if (pu .eq. 0) go to 999
 745       iv1 = iv(1)
 746       if (iv1 .gt. 62) iv1 = iv1 - 51
 747       ol = iv(outlev)
 748       alg = iv(algsav)
 749       if (iv1 .lt. 2 .or. iv1 .gt. 15) go to 370
 750       if (iv1 .ge. 12) go to 120
 751       if (iv1 .eq. 2 .and. iv(niter) .eq. 0) go to 390
 752       if (ol .eq. 0) go to 120
 753       if (iv1 .ge. 10 .and. iv(prntit) .eq. 0) go to 120
 754       if (iv1 .gt. 2) go to 10
 755          iv(prntit) = iv(prntit) + 1
 756          if (iv(prntit) .lt. iabs(ol)) go to 999
 757  10   nf = iv(nfcall) - iabs(iv(nfcov))
 758       iv(prntit) = 0
 759       reldf = zero
 760       preldf = zero
 761       oldf = dmax1(dabs(v(f0)), dabs(v(f)))
 762       if (oldf .le. zero) go to 20
 763          reldf = v(fdif) / oldf
 764          preldf = v(preduc) / oldf
 765  20   if (ol .gt. 0) go to 60
 766 c
 767 c        ***  print short summary line  ***
 768 c
 769          if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,30)
 770  30   format(/10h   it   nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,
 771      1       2x,13hmodel  stppar)
 772          if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,40)
 773  40   format(/11h    it   nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,
 774      1       3x,6hstppar)
 775          iv(needhd) = 0
 776          if (alg .eq. 2) go to 50
 777          m = iv(sused)
 778          write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),
 779      1                 model1(m), model2(m), v(stppar)
 780          go to 120
 781 c
 782  50      write(pu,110) iv(niter), nf, v(f), reldf, preldf, v(reldx),
 783      1                 v(stppar)
 784          go to 120
 785 c
 786 c     ***  print long summary line  ***
 787 c
 788  60   if (iv(needhd) .eq. 1 .and. alg .eq. 1) write(pu,70)
 789  70   format(/11h    it   nf,6x,1hf,7x,5hreldf,3x,6hpreldf,3x,5hreldx,
 790      1       2x,13hmodel  stppar,2x,6hd*step,2x,7hnpreldf)
 791       if (iv(needhd) .eq. 1 .and. alg .eq. 2) write(pu,80)
 792  80   format(/11h    it   nf,7x,1hf,8x,5hreldf,4x,6hpreldf,4x,5hreldx,
 793      1       3x,6hstppar,3x,6hd*step,3x,7hnpreldf)
 794       iv(needhd) = 0
 795       nreldf = zero
 796       if (oldf .gt. zero) nreldf = v(nreduc) / oldf
 797       if (alg .eq. 2) go to 90
 798       m = iv(sused)
 799       write(pu,100) iv(niter), nf, v(f), reldf, preldf, v(reldx),
 800      1             model1(m), model2(m), v(stppar), v(dstnrm), nreldf
 801       go to 120
 802 c
 803  90   write(pu,110) iv(niter), nf, v(f), reldf, preldf,
 804      1             v(reldx), v(stppar), v(dstnrm), nreldf
 805  100  format(i6,i5,d10.3,2d9.2,d8.1,a3,a4,2d8.1,d9.2)
 806  110  format(i6,i5,d11.3,2d10.2,3d9.1,d10.2)
 807 c
 808  120  if (iv(statpr) .lt. 0) go to 430
 809       go to (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,
 810      1       330, 350, 520), iv1
 811 c
 812  130  write(pu,140)
 813  140  format(/26h ***** x-convergence *****)
 814       go to 430
 815 c
 816  150  write(pu,160)
 817  160  format(/42h ***** relative function convergence *****)
 818       go to 430
 819 c
 820  170  write(pu,180)
 821  180  format(/49h ***** x- and relative function convergence *****)
 822       go to 430
 823 c
 824  190  write(pu,200)
 825  200  format(/42h ***** absolute function convergence *****)
 826       go to 430
 827 c
 828  210  write(pu,220)
 829  220  format(/33h ***** singular convergence *****)
 830       go to 430
 831 c
 832  230  write(pu,240)
 833  240  format(/30h ***** false convergence *****)
 834       go to 430
 835 c
 836  250  write(pu,260)
 837  260  format(/38h ***** function evaluation limit *****)
 838       go to 430
 839 c
 840  270  write(pu,280)
 841  280  format(/28h ***** iteration limit *****)
 842       go to 430
 843 c
 844  290  write(pu,300)
 845  300  format(/18h ***** stopx *****)
 846       go to 430
 847 c
 848  310  write(pu,320)
 849  320  format(/44h ***** initial f(x) cannot be computed *****)
 850 c
 851       go to 390
 852 c
 853  330  write(pu,340)
 854  340  format(/37h ***** bad parameters to assess *****)
 855       go to 999
 856 c
 857  350  write(pu,360)
 858  360  format(/43h ***** gradient could not be computed *****)
 859       if (iv(niter) .gt. 0) go to 480
 860       go to 390
 861 c
 862  370  write(pu,380) iv(1)
 863  380  format(/14h ***** iv(1) =,i5,6h *****)
 864       go to 999
 865 c
 866 c  ***  initial call on itsum  ***
 867 c
 868  390  if (iv(x0prt) .ne. 0) write(pu,400) (i, x(i), d(i), i = 1, p)
 869  400  format(/23h     i     initial x(i),8x,4hd(i)//(1x,i5,d17.6,d14.3))
 870 c     *** the following are to avoid undefined variables when the
 871 c     *** function evaluation limit is 1...
 872       v(dstnrm) = zero
 873       v(fdif) = zero
 874       v(nreduc) = zero
 875       v(preduc) = zero
 876       v(reldx)  = zero
 877       if (iv1 .ge. 12) go to 999
 878       iv(needhd) = 0
 879       iv(prntit) = 0
 880       if (ol .eq. 0) go to 999
 881       if (ol .lt. 0 .and. alg .eq. 1) write(pu,30)
 882       if (ol .lt. 0 .and. alg .eq. 2) write(pu,40)
 883       if (ol .gt. 0 .and. alg .eq. 1) write(pu,70)
 884       if (ol .gt. 0 .and. alg .eq. 2) write(pu,80)
 885       if (alg .eq. 1) write(pu,410) v(f)
 886       if (alg .eq. 2) write(pu,420) v(f)
 887  410  format(/11h     0    1,d10.3)
 888 c365  format(/11h     0    1,e11.3)
 889  420  format(/11h     0    1,d11.3)
 890       go to 999
 891 c
 892 c  ***  print various information requested on solution  ***
 893 c
 894  430  iv(needhd) = 1
 895       if (iv(statpr) .eq. 0) go to 480
 896          oldf = dmax1(dabs(v(f0)), dabs(v(f)))
 897          preldf = zero
 898          nreldf = zero
 899          if (oldf .le. zero) go to 440
 900               preldf = v(preduc) / oldf
 901               nreldf = v(nreduc) / oldf
 902  440     nf = iv(nfcall) - iv(nfcov)
 903          ng = iv(ngcall) - iv(ngcov)
 904          write(pu,450) v(f), v(reldx), nf, ng, preldf, nreldf
 905  450  format(/9h function,d17.6,8h   reldx,d17.3/12h func. evals,
 906      1   i8,9x,11hgrad. evals,i8/7h preldf,d16.3,6x,7hnpreldf,d15.3)
 907 c
 908          if (iv(nfcov) .gt. 0) write(pu,460) iv(nfcov)
 909  460     format(/1x,i4,50h extra func. evals for covariance and diagnost
 910      1ics.)
 911          if (iv(ngcov) .gt. 0) write(pu,470) iv(ngcov)
 912  470     format(1x,i4,50h extra grad. evals for covariance and diagnosti
 913      1cs.)
 914 c
 915  480  if (iv(solprt) .eq. 0) go to 999
 916          iv(needhd) = 1
 917          write(pu,490)
 918  490  format(/22h     i      final x(i),8x,4hd(i),10x,4hg(i)/)
 919          do 500 i = 1, p
 920               write(pu,510) i, x(i), d(i), g(i)
 921  500          continue
 922  510     format(1x,i5,d16.6,2d14.3)
 923       go to 999
 924 c
 925  520  write(pu,530)
 926  530  format(/24h inconsistent dimensions)
 927  999  return
 928 c  ***  last card of itsum follows  ***
 929       end
 930       subroutine litvmu(n, x, l, y)
 931 c
 932 c  ***  solve  (l**t)*x = y,  where  l  is an  n x n  lower triangular
 933 c  ***  matrix stored compactly by rows.  x and y may occupy the same
 934 c  ***  storage.  ***
 935 c
 936       integer n
 937 cal   double precision x(n), l(1), y(n)
 938       double precision x(n), l(n*(n+1)/2), y(n)
 939       integer i, ii, ij, im1, i0, j, np1
 940       double precision xi, zero
 941 c/6
 942 c     data zero/0.d+0/
 943 c/7
 944       parameter (zero=0.d+0)
 945 c/
 946 c
 947       do 10 i = 1, n
 948  10      x(i) = y(i)
 949       np1 = n + 1
 950       i0 = n*(n+1)/2
 951       do 30 ii = 1, n
 952          i = np1 - ii
 953          xi = x(i)/l(i0)
 954          x(i) = xi
 955          if (i .le. 1) go to 999
 956          i0 = i0 - i
 957          if (xi .eq. zero) go to 30
 958          im1 = i - 1
 959          do 20 j = 1, im1
 960               ij = i0 + j
 961               x(j) = x(j) - xi*l(ij)
 962  20           continue
 963  30      continue
 964  999  return
 965 c  ***  last card of litvmu follows  ***
 966       end
 967       subroutine livmul(n, x, l, y)
 968 c
 969 c  ***  solve  l*x = y, where  l  is an  n x n  lower triangular
 970 c  ***  matrix stored compactly by rows.  x and y may occupy the same
 971 c  ***  storage.  ***
 972 c
 973       integer n
 974 cal   double precision x(n), l(1), y(n)
 975       double precision x(n), l(n*(n+1)/2), y(n)
 976       external dotprd
 977       double precision dotprd
 978       integer i, j, k
 979       double precision t, zero
 980 c/6
 981 c     data zero/0.d+0/
 982 c/7
 983       parameter (zero=0.d+0)
 984 c/
 985 c
 986       do 10 k = 1, n
 987          if (y(k) .ne. zero) go to 20
 988          x(k) = zero
 989  10      continue
 990       go to 999
 991  20   j = k*(k+1)/2
 992       x(k) = y(k) / l(j)
 993       if (k .ge. n) go to 999
 994       k = k + 1
 995       do 30 i = k, n
 996          t = dotprd(i-1, l(j+1), x)
 997          j = j + i
 998          x(i) = (y(i) - t)/l(j)
 999  30      continue
1000  999  return
1001 c  ***  last card of livmul follows  ***
1002       end
1003       subroutine parck(alg, d, iv, liv, lv, n, v)
1004 c
1005 c  ***  check ***sol (version 2.3) parameters, print changed values  ***
1006 c
1007 c  ***  alg = 1 for regression, alg = 2 for general unconstrained opt.
1008 c
1009       integer alg, liv, lv, n
1010       integer iv(liv)
1011       double precision d(n), v(lv)
1012 c
1013       external rmdcon, vcopy, vdflt
1014       double precision rmdcon
1015 c rmdcon -- returns machine-dependent constants.
1016 c vcopy  -- copies one vector to another.
1017 c vdflt  -- supplies default parameter values to v alone.
1018 c/+
1019       integer max0
1020 c/
1021 c
1022 c  ***  local variables  ***
1023 c
1024       integer i, ii, iv1, j, k, l, m, miv1, miv2, ndfalt, parsv1, pu
1025       integer ijmp, jlim(2), miniv(2), ndflt(2)
1026 c/6
1027 c     integer varnm(2), sh(2)
1028 c     real cngd(3), dflt(3), vn(2,34), which(3)
1029 c/7
1030       character*1 varnm(2), sh(2)
1031       character*4 cngd(3), dflt(3), vn(2,34), which(3)
1032 c/
1033       double precision big, machep, tiny, vk, vm(34), vx(34), zero
1034 c
1035 c  ***  iv and v subscripts  ***
1036 c
1037       integer algsav, dinit, dtype, dtype0, epslon, inits, ivneed,
1038      1        lastiv, lastv, lmat, nextiv, nextv, nvdflt, oldn,
1039      2        parprt, parsav, perm, prunit, vneed
1040 c
1041 c
1042 c/6
1043 c     data algsav/51/, dinit/38/, dtype/16/, dtype0/54/, epslon/19/,
1044 c    1     inits/25/, ivneed/3/, lastiv/44/, lastv/45/, lmat/42/,
1045 c    2     nextiv/46/, nextv/47/, nvdflt/50/, oldn/38/, parprt/20/,
1046 c    3     parsav/49/, perm/58/, prunit/21/, vneed/4/
1047 c/7
1048       parameter (algsav=51, dinit=38, dtype=16, dtype0=54, epslon=19,
1049      1           inits=25, ivneed=3, lastiv=44, lastv=45, lmat=42,
1050      2           nextiv=46, nextv=47, nvdflt=50, oldn=38, parprt=20,
1051      3           parsav=49, perm=58, prunit=21, vneed=4)
1052       save big, machep, tiny
1053 c/
1054 c
1055       data big/0.d+0/, machep/-1.d+0/, tiny/1.d+0/, zero/0.d+0/
1056 c/6
1057 c     data vn(1,1),vn(2,1)/4hepsl,4hon../
1058 c     data vn(1,2),vn(2,2)/4hphmn,4hfc../
1059 c     data vn(1,3),vn(2,3)/4hphmx,4hfc../
1060 c     data vn(1,4),vn(2,4)/4hdecf,4hac../
1061 c     data vn(1,5),vn(2,5)/4hincf,4hac../
1062 c     data vn(1,6),vn(2,6)/4hrdfc,4hmn../
1063 c     data vn(1,7),vn(2,7)/4hrdfc,4hmx../
1064 c     data vn(1,8),vn(2,8)/4htune,4hr1../
1065 c     data vn(1,9),vn(2,9)/4htune,4hr2../
1066 c     data vn(1,10),vn(2,10)/4htune,4hr3../
1067 c     data vn(1,11),vn(2,11)/4htune,4hr4../
1068 c     data vn(1,12),vn(2,12)/4htune,4hr5../
1069 c     data vn(1,13),vn(2,13)/4hafct,4hol../
1070 c     data vn(1,14),vn(2,14)/4hrfct,4hol../
1071 c     data vn(1,15),vn(2,15)/4hxcto,4hl.../
1072 c     data vn(1,16),vn(2,16)/4hxfto,4hl.../
1073 c     data vn(1,17),vn(2,17)/4hlmax,4h0.../
1074 c     data vn(1,18),vn(2,18)/4hlmax,4hs.../
1075 c     data vn(1,19),vn(2,19)/4hscto,4hl.../
1076 c     data vn(1,20),vn(2,20)/4hdini,4ht.../
1077 c     data vn(1,21),vn(2,21)/4hdtin,4hit../
1078 c     data vn(1,22),vn(2,22)/4hd0in,4hit../
1079 c     data vn(1,23),vn(2,23)/4hdfac,4h..../
1080 c     data vn(1,24),vn(2,24)/4hdltf,4hdc../
1081 c     data vn(1,25),vn(2,25)/4hdltf,4hdj../
1082 c     data vn(1,26),vn(2,26)/4hdelt,4ha0../
1083 c     data vn(1,27),vn(2,27)/4hfuzz,4h..../
1084 c     data vn(1,28),vn(2,28)/4hrlim,4hit../
1085 c     data vn(1,29),vn(2,29)/4hcosm,4hin../
1086 c     data vn(1,30),vn(2,30)/4hhube,4hrc../
1087 c     data vn(1,31),vn(2,31)/4hrspt,4hol../
1088 c     data vn(1,32),vn(2,32)/4hsigm,4hin../
1089 c     data vn(1,33),vn(2,33)/4heta0,4h..../
1090 c     data vn(1,34),vn(2,34)/4hbias,4h..../
1091 c/7
1092       data vn(1,1),vn(2,1)/'epsl','on..'/
1093       data vn(1,2),vn(2,2)/'phmn','fc..'/
1094       data vn(1,3),vn(2,3)/'phmx','fc..'/
1095       data vn(1,4),vn(2,4)/'decf','ac..'/
1096       data vn(1,5),vn(2,5)/'incf','ac..'/
1097       data vn(1,6),vn(2,6)/'rdfc','mn..'/
1098       data vn(1,7),vn(2,7)/'rdfc','mx..'/
1099       data vn(1,8),vn(2,8)/'tune','r1..'/
1100       data vn(1,9),vn(2,9)/'tune','r2..'/
1101       data vn(1,10),vn(2,10)/'tune','r3..'/
1102       data vn(1,11),vn(2,11)/'tune','r4..'/
1103       data vn(1,12),vn(2,12)/'tune','r5..'/
1104       data vn(1,13),vn(2,13)/'afct','ol..'/
1105       data vn(1,14),vn(2,14)/'rfct','ol..'/
1106       data vn(1,15),vn(2,15)/'xcto','l...'/
1107       data vn(1,16),vn(2,16)/'xfto','l...'/
1108       data vn(1,17),vn(2,17)/'lmax','0...'/
1109       data vn(1,18),vn(2,18)/'lmax','s...'/
1110       data vn(1,19),vn(2,19)/'scto','l...'/
1111       data vn(1,20),vn(2,20)/'dini','t...'/
1112       data vn(1,21),vn(2,21)/'dtin','it..'/
1113       data vn(1,22),vn(2,22)/'d0in','it..'/
1114       data vn(1,23),vn(2,23)/'dfac','....'/
1115       data vn(1,24),vn(2,24)/'dltf','dc..'/
1116       data vn(1,25),vn(2,25)/'dltf','dj..'/
1117       data vn(1,26),vn(2,26)/'delt','a0..'/
1118       data vn(1,27),vn(2,27)/'fuzz','....'/
1119       data vn(1,28),vn(2,28)/'rlim','it..'/
1120       data vn(1,29),vn(2,29)/'cosm','in..'/
1121       data vn(1,30),vn(2,30)/'hube','rc..'/
1122       data vn(1,31),vn(2,31)/'rspt','ol..'/
1123       data vn(1,32),vn(2,32)/'sigm','in..'/
1124       data vn(1,33),vn(2,33)/'eta0','....'/
1125       data vn(1,34),vn(2,34)/'bias','....'/
1126 c/
1127 c
1128       data vm(1)/1.0d-3/, vm(2)/-0.99d+0/, vm(3)/1.0d-3/, vm(4)/1.0d-2/,
1129      1     vm(5)/1.2d+0/, vm(6)/1.d-2/, vm(7)/1.2d+0/, vm(8)/0.d+0/,
1130      2     vm(9)/0.d+0/, vm(10)/1.d-3/, vm(11)/-1.d+0/, vm(13)/0.d+0/,
1131      3     vm(15)/0.d+0/, vm(16)/0.d+0/, vm(19)/0.d+0/, vm(20)/-10.d+0/,
1132      4     vm(21)/0.d+0/, vm(22)/0.d+0/, vm(23)/0.d+0/, vm(27)/1.01d+0/,
1133      5     vm(28)/1.d+10/, vm(30)/0.d+0/, vm(31)/0.d+0/, vm(32)/0.d+0/,
1134      6     vm(34)/0.d+0/
1135       data vx(1)/0.9d+0/, vx(2)/-1.d-3/, vx(3)/1.d+1/, vx(4)/0.8d+0/,
1136      1     vx(5)/1.d+2/, vx(6)/0.8d+0/, vx(7)/1.d+2/, vx(8)/0.5d+0/,
1137      2     vx(9)/0.5d+0/, vx(10)/1.d+0/, vx(11)/1.d+0/, vx(14)/0.1d+0/,
1138      3     vx(15)/1.d+0/, vx(16)/1.d+0/, vx(19)/1.d+0/, vx(23)/1.d+0/,
1139      4     vx(24)/1.d+0/, vx(25)/1.d+0/, vx(26)/1.d+0/, vx(27)/1.d+10/,
1140      5     vx(29)/1.d+0/, vx(31)/1.d+0/, vx(32)/1.d+0/, vx(33)/1.d+0/,
1141      6     vx(34)/1.d+0/
1142 c
1143 c/6
1144 c     data varnm(1)/1hp/, varnm(2)/1hn/, sh(1)/1hs/, sh(2)/1hh/
1145 c     data cngd(1),cngd(2),cngd(3)/4h---c,4hhang,4hed v/,
1146 c    1     dflt(1),dflt(2),dflt(3)/4hnond,4hefau,4hlt v/
1147 c/7
1148       data varnm(1)/'p'/, varnm(2)/'n'/, sh(1)/'s'/, sh(2)/'h'/
1149       data cngd(1),cngd(2),cngd(3)/'---c','hang','ed v'/,
1150      1     dflt(1),dflt(2),dflt(3)/'nond','efau','lt v'/
1151 c/
1152       data ijmp/33/, jlim(1)/0/, jlim(2)/24/, ndflt(1)/32/, ndflt(2)/25/
1153       data miniv(1)/80/, miniv(2)/59/
1154 c
1155 c...............................  body  ................................
1156 c
1157       pu = 0
1158       if (prunit .le. liv) pu = iv(prunit)
1159       if (alg .lt. 1 .or. alg .gt. 2) go to 340
1160       if (iv(1) .eq. 0) call deflt(alg, iv, liv, lv, v)
1161       iv1 = iv(1)
1162       if (iv1 .ne. 13 .and. iv1 .ne. 12) go to 10
1163       miv1 = miniv(alg)
1164       if (perm .le. liv) miv1 = max0(miv1, iv(perm) - 1)
1165       if (ivneed .le. liv) miv2 = miv1 + max0(iv(ivneed), 0)
1166       if (lastiv .le. liv) iv(lastiv) = miv2
1167       if (liv .lt. miv1) go to 300
1168       iv(ivneed) = 0
1169       iv(lastv) = max0(iv(vneed), 0) + iv(lmat) - 1
1170       iv(vneed) = 0
1171       if (liv .lt. miv2) go to 300
1172       if (lv .lt. iv(lastv)) go to 320
1173  10   if (alg .eq. iv(algsav)) go to 30
1174          if (pu .ne. 0) write(pu,20) alg, iv(algsav)
1175  20      format(/39h the first parameter to deflt should be,i3,
1176      1          12h rather than,i3)
1177          iv(1) = 82
1178          go to 999
1179  30   if (iv1 .lt. 12 .or. iv1 .gt. 14) go to 60
1180          if (n .ge. 1) go to 50
1181               iv(1) = 81
1182               if (pu .eq. 0) go to 999
1183               write(pu,40) varnm(alg), n
1184  40           format(/8h /// bad,a1,2h =,i5)
1185               go to 999
1186  50      if (iv1 .ne. 14) iv(nextiv) = iv(perm)
1187          if (iv1 .ne. 14) iv(nextv) = iv(lmat)
1188          if (iv1 .eq. 13) go to 999
1189          k = iv(parsav) - epslon
1190          call vdflt(alg, lv-k, v(k+1))
1191          iv(dtype0) = 2 - alg
1192          iv(oldn) = n
1193          which(1) = dflt(1)
1194          which(2) = dflt(2)
1195          which(3) = dflt(3)
1196          go to 110
1197  60   if (n .eq. iv(oldn)) go to 80
1198          iv(1) = 17
1199          if (pu .eq. 0) go to 999
1200          write(pu,70) varnm(alg), iv(oldn), n
1201  70      format(/5h /// ,1a1,14h changed from ,i5,4h to ,i5)
1202          go to 999
1203 c
1204  80   if (iv1 .le. 11 .and. iv1 .ge. 1) go to 100
1205          iv(1) = 80
1206          if (pu .ne. 0) write(pu,90) iv1
1207  90      format(/13h ///  iv(1) =,i5,28h should be between 0 and 14.)
1208          go to 999
1209 c
1210  100  which(1) = cngd(1)
1211       which(2) = cngd(2)
1212       which(3) = cngd(3)
1213 c
1214  110  if (iv1 .eq. 14) iv1 = 12
1215       if (big .gt. tiny) go to 120
1216          tiny = rmdcon(1)
1217          machep = rmdcon(3)
1218          big = rmdcon(6)
1219          vm(12) = machep
1220          vx(12) = big
1221          vx(13) = big
1222          vm(14) = machep
1223          vm(17) = tiny
1224          vx(17) = big
1225          vm(18) = tiny
1226          vx(18) = big
1227          vx(20) = big
1228          vx(21) = big
1229          vx(22) = big
1230          vm(24) = machep
1231          vm(25) = machep
1232          vm(26) = machep
1233          vx(28) = rmdcon(5)
1234          vm(29) = machep
1235          vx(30) = big
1236          vm(33) = machep
1237  120  m = 0
1238       i = 1
1239       j = jlim(alg)
1240       k = epslon
1241       ndfalt = ndflt(alg)
1242       do 150 l = 1, ndfalt
1243          vk = v(k)
1244          if (vk .ge. vm(i) .and. vk .le. vx(i)) go to 140
1245               m = k
1246               if (pu .ne. 0) write(pu,130) vn(1,i), vn(2,i), k, vk,
1247      1                                    vm(i), vx(i)
1248  130          format(/6h ///  ,2a4,5h.. v(,i2,3h) =,d11.3,7h should,
1249      1               11h be between,d11.3,4h and,d11.3)
1250  140     k = k + 1
1251          i = i + 1
1252          if (i .eq. j) i = ijmp
1253  150     continue
1254 c
1255       if (iv(nvdflt) .eq. ndfalt) go to 170
1256          iv(1) = 51
1257          if (pu .eq. 0) go to 999
1258          write(pu,160) iv(nvdflt), ndfalt
1259  160     format(/13h iv(nvdflt) =,i5,13h rather than ,i5)
1260          go to 999
1261  170  if ((iv(dtype) .gt. 0 .or. v(dinit) .gt. zero) .and. iv1 .eq. 12)
1262      1                  go to 200
1263       do 190 i = 1, n
1264          if (d(i) .gt. zero) go to 190
1265               m = 18
1266               if (pu .ne. 0) write(pu,180) i, d(i)
1267  180     format(/8h ///  d(,i3,3h) =,d11.3,19h should be positive)
1268  190     continue
1269  200  if (m .eq. 0) go to 210
1270          iv(1) = m
1271          go to 999
1272 c
1273  210  if (pu .eq. 0 .or. iv(parprt) .eq. 0) go to 999
1274       if (iv1 .ne. 12 .or. iv(inits) .eq. alg-1) go to 230
1275          m = 1
1276          write(pu,220) sh(alg), iv(inits)
1277  220     format(/22h nondefault values..../5h init,a1,14h..... iv(25) =,
1278      1          i3)
1279  230  if (iv(dtype) .eq. iv(dtype0)) go to 250
1280          if (m .eq. 0) write(pu,260) which
1281          m = 1
1282          write(pu,240) iv(dtype)
1283  240     format(20h dtype..... iv(16) =,i3)
1284  250  i = 1
1285       j = jlim(alg)
1286       k = epslon
1287       l = iv(parsav)
1288       ndfalt = ndflt(alg)
1289       do 290 ii = 1, ndfalt
1290          if (v(k) .eq. v(l)) go to 280
1291               if (m .eq. 0) write(pu,260) which
1292  260          format(/1h ,3a4,9halues..../)
1293               m = 1
1294               write(pu,270) vn(1,i), vn(2,i), k, v(k)
1295  270          format(1x,2a4,5h.. v(,i2,3h) =,d15.7)
1296  280     k = k + 1
1297          l = l + 1
1298          i = i + 1
1299          if (i .eq. j) i = ijmp
1300  290     continue
1301 c
1302       iv(dtype0) = iv(dtype)
1303       parsv1 = iv(parsav)
1304       call vcopy(iv(nvdflt), v(parsv1), v(epslon))
1305       go to 999
1306 c
1307  300  iv(1) = 15
1308       if (pu .eq. 0) go to 999
1309       write(pu,310) liv, miv2
1310  310  format(/10h /// liv =,i5,17h must be at least,i5)
1311       if (liv .lt. miv1) go to 999
1312       if (lv .lt. iv(lastv)) go to 320
1313       go to 999
1314 c
1315  320  iv(1) = 16
1316       if (pu .eq. 0) go to 999
1317       write(pu,330) lv, iv(lastv)
1318  330  format(/9h /// lv =,i5,17h must be at least,i5)
1319       go to 999
1320 c
1321  340  iv(1) = 67
1322       if (pu .eq. 0) go to 999
1323       write(pu,350) alg
1324  350  format(/10h /// alg =,i5,15h must be 1 or 2)
1325 c
1326  999  return
1327 c  ***  last card of parck follows  ***
1328       end
1329       double precision function reldst(p, d, x, x0)
1330 c
1331 c  ***  compute and return relative difference between x and x0  ***
1332 c  ***  nl2sol version 2.2  ***
1333 c
1334       integer p
1335       double precision d(p), x(p), x0(p)
1336 c/+
1337       double precision dabs
1338 c/
1339       integer i
1340       double precision emax, t, xmax, zero
1341 c/6
1342 c     data zero/0.d+0/
1343 c/7
1344       parameter (zero=0.d+0)
1345 c/
1346 c
1347       emax = zero
1348       xmax = zero
1349       do 10 i = 1, p
1350          t = dabs(d(i) * (x(i) - x0(i)))
1351          if (emax .lt. t) emax = t
1352          t = d(i) * (dabs(x(i)) + dabs(x0(i)))
1353          if (xmax .lt. t) xmax = t
1354  10      continue
1355       reldst = zero
1356       if (xmax .gt. zero) reldst = emax / xmax
1357  999  return
1358 c  ***  last card of reldst follows  ***
1359       end
1360 c     logical function stopx(idummy)
1361 c     *****parameters...
1362 c     integer idummy
1363 c
1364 c     ..................................................................
1365 c
1366 c     *****purpose...
1367 c     this function may serve as the stopx (asynchronous interruption)
1368 c     function for the nl2sol (nonlinear least-squares) package at
1369 c     those installations which do not wish to implement a
1370 c     dynamic stopx.
1371 c
1372 c     *****algorithm notes...
1373 c     at installations where the nl2sol system is used
1374 c     interactively, this dummy stopx should be replaced by a
1375 c     function that returns .true. if and only if the interrupt
1376 c     (break) key has been pressed since the last call on stopx.
1377 c
1378 c     ..................................................................
1379 c
1380 c     stopx = .false.
1381 c     return
1382 c     end
1383       subroutine vaxpy(p, w, a, x, y)
1384 c
1385 c  ***  set w = a*x + y  --  w, x, y = p-vectors, a = scalar  ***
1386 c
1387       integer p
1388       double precision a, w(p), x(p), y(p)
1389 c
1390       integer i
1391 c
1392       do 10 i = 1, p
1393  10      w(i) = a*x(i) + y(i)
1394       return
1395       end
1396       subroutine vcopy(p, y, x)
1397 c
1398 c  ***  set y = x, where x and y are p-vectors  ***
1399 c
1400       integer p
1401       double precision x(p), y(p)
1402 c
1403       integer i
1404 c
1405       do 10 i = 1, p
1406  10      y(i) = x(i)
1407       return
1408       end
1409       subroutine vdflt(alg, lv, v)
1410 c
1411 c  ***  supply ***sol (version 2.3) default values to v  ***
1412 c
1413 c  ***  alg = 1 means regression constants.
1414 c  ***  alg = 2 means general unconstrained optimization constants.
1415 c
1416       integer alg, l
1417       double precision v(lv)
1418 c/+
1419       double precision dmax1
1420 c/
1421       external rmdcon
1422       double precision rmdcon
1423 c rmdcon... returns machine-dependent constants
1424 c
1425       double precision machep, mepcrt, one, sqteps, three
1426 c
1427 c  ***  subscripts for v  ***
1428 c
1429       integer afctol, bias, cosmin, decfac, delta0, dfac, dinit, dltfdc,
1430      1        dltfdj, dtinit, d0init, epslon, eta0, fuzz, huberc,
1431      2        incfac, lmax0, lmaxs, phmnfc, phmxfc, rdfcmn, rdfcmx,
1432      3        rfctol, rlimit, rsptol, sctol, sigmin, tuner1, tuner2,
1433      4        tuner3, tuner4, tuner5, xctol, xftol
1434 c
1435 c/6
1436 c     data one/1.d+0/, three/3.d+0/
1437 c/7
1438       parameter (one=1.d+0, three=3.d+0)
1439 c/
1440 c
1441 c  ***  v subscript values  ***
1442 c
1443 c/6
1444 c     data afctol/31/, bias/43/, cosmin/47/, decfac/22/, delta0/44/,
1445 c    1     dfac/41/, dinit/38/, dltfdc/42/, dltfdj/43/, dtinit/39/,
1446 c    2     d0init/40/, epslon/19/, eta0/42/, fuzz/45/, huberc/48/,
1447 c    3     incfac/23/, lmax0/35/, lmaxs/36/, phmnfc/20/, phmxfc/21/,
1448 c    4     rdfcmn/24/, rdfcmx/25/, rfctol/32/, rlimit/46/, rsptol/49/,
1449 c    5     sctol/37/, sigmin/50/, tuner1/26/, tuner2/27/, tuner3/28/,
1450 c    6     tuner4/29/, tuner5/30/, xctol/33/, xftol/34/
1451 c/7
1452       parameter (afctol=31, bias=43, cosmin=47, decfac=22, delta0=44,
1453      1           dfac=41, dinit=38, dltfdc=42, dltfdj=43, dtinit=39,
1454      2           d0init=40, epslon=19, eta0=42, fuzz=45, huberc=48,
1455      3           incfac=23, lmax0=35, lmaxs=36, phmnfc=20, phmxfc=21,
1456      4           rdfcmn=24, rdfcmx=25, rfctol=32, rlimit=46, rsptol=49,
1457      5           sctol=37, sigmin=50, tuner1=26, tuner2=27, tuner3=28,
1458      6           tuner4=29, tuner5=30, xctol=33, xftol=34)
1459 c/
1460 c
1461 c-------------------------------  body  --------------------------------
1462 c
1463       machep = rmdcon(3)
1464       v(afctol) = 1.d-20
1465       if (machep .gt. 1.d-10) v(afctol) = machep**2
1466       v(decfac) = 0.5d+0
1467       sqteps = rmdcon(4)
1468       v(dfac) = 0.6d+0
1469       v(delta0) = sqteps
1470       v(dtinit) = 1.d-6
1471       mepcrt = machep ** (one/three)
1472       v(d0init) = 1.d+0
1473       v(epslon) = 0.1d+0
1474       v(incfac) = 2.d+0
1475       v(lmax0) = 1.d+0
1476       v(lmaxs) = 1.d+0
1477       v(phmnfc) = -0.1d+0
1478       v(phmxfc) = 0.1d+0
1479       v(rdfcmn) = 0.1d+0
1480       v(rdfcmx) = 4.d+0
1481       v(rfctol) = dmax1(1.d-10, mepcrt**2)
1482       v(sctol) = v(rfctol)
1483       v(tuner1) = 0.1d+0
1484       v(tuner2) = 1.d-4
1485       v(tuner3) = 0.75d+0
1486       v(tuner4) = 0.5d+0
1487       v(tuner5) = 0.75d+0
1488       v(xctol) = sqteps
1489       v(xftol) = 1.d+2 * machep
1490 c
1491       if (alg .ge. 2) go to 10
1492 c
1493 c  ***  regression  values
1494 c
1495       v(cosmin) = dmax1(1.d-6, 1.d+2 * machep)
1496       v(dinit) = 0.d+0
1497       v(dltfdc) = mepcrt
1498       v(dltfdj) = sqteps
1499       v(fuzz) = 1.5d+0
1500       v(huberc) = 0.7d+0
1501       v(rlimit) = rmdcon(5)
1502       v(rsptol) = 1.d-3
1503       v(sigmin) = 1.d-4
1504       go to 999
1505 c
1506 c  ***  general optimization values
1507 c
1508  10   v(bias) = 0.8d+0
1509       v(dinit) = -1.0d+0
1510       v(eta0) = 1.0d+3 * machep
1511 c
1512  999  return
1513 c  ***  last card of vdflt follows  ***
1514       end
1515       subroutine vscopy(p, y, s)
1516 c
1517 c  ***  set p-vector y to scalar s  ***
1518 c
1519       integer p
1520       double precision s, y(p)
1521 c
1522       integer i
1523 c
1524       do 10 i = 1, p
1525  10      y(i) = s
1526       return
1527       end
1528       double precision function v2norm(p, x)
1529 c
1530 c  ***  return the 2-norm of the p-vector x, taking  ***
1531 c  ***  care to avoid the most likely underflows.    ***
1532 c
1533       integer p
1534       double precision x(p)
1535 c
1536       integer i, j
1537       double precision one, r, scale, sqteta, t, xi, zero
1538 c/+
1539       double precision dabs, dsqrt
1540 c/
1541       external rmdcon
1542       double precision rmdcon
1543 c
1544 c/6
1545 c     data one/1.d+0/, zero/0.d+0/
1546 c/7
1547       parameter (one=1.d+0, zero=0.d+0)
1548       save sqteta
1549 c/
1550       data sqteta/0.d+0/
1551 c
1552       if (p .gt. 0) go to 10
1553          v2norm = zero
1554          go to 999
1555  10   do 20 i = 1, p
1556          if (x(i) .ne. zero) go to 30
1557  20      continue
1558       v2norm = zero
1559       go to 999
1560 c
1561  30   scale = dabs(x(i))
1562       if (i .lt. p) go to 40
1563          v2norm = scale
1564          go to 999
1565  40   t = one
1566       if (sqteta .eq. zero) sqteta = rmdcon(2)
1567 c
1568 c     ***  sqteta is (slightly larger than) the square root of the
1569 c     ***  smallest positive floating point number on the machine.
1570 c     ***  the tests involving sqteta are done to prevent underflows.
1571 c
1572       j = i + 1
1573       do 60 i = j, p
1574          xi = dabs(x(i))
1575          if (xi .gt. scale) go to 50
1576               r = xi / scale
1577               if (r .gt. sqteta) t = t + r*r
1578               go to 60
1579  50           r = scale / xi
1580               if (r .le. sqteta) r = zero
1581               t = one  +  t * r*r
1582               scale = xi
1583  60      continue
1584 c
1585       v2norm = scale * dsqrt(t)
1586  999  return
1587 c  ***  last card of v2norm follows  ***
1588       end
1589       subroutine humsl(n, d, x, calcf, calcgh, iv, liv, lv, v,
1590      1                  uiparm, urparm, ufparm)
1591 c
1592 c  ***  minimize general unconstrained objective function using   ***
1593 c  ***  (analytic) gradient and hessian provided by the caller.   ***
1594 c
1595       integer liv, lv, n
1596       integer iv(liv), uiparm(1)
1597       double precision d(n), x(n), v(lv), urparm(1)
1598 c     dimension v(78 + n*(n+12)), uiparm(*), urparm(*)
1599       external calcf, calcgh, ufparm
1600 c
1601 c------------------------------  discussion  ---------------------------
1602 c
1603 c        this routine is like sumsl, except that the subroutine para-
1604 c     meter calcg of sumsl (which computes the gradient of the objec-
1605 c     tive function) is replaced by the subroutine parameter calcgh,
1606 c     which computes both the gradient and (lower triangle of the)
1607 c     hessian of the objective function.  the calling sequence is...
1608 c             call calcgh(n, x, nf, g, h, uiparm, urparm, ufparm)
1609 c     parameters n, x, nf, g, uiparm, urparm, and ufparm are the same
1610 c     as for sumsl, while h is an array of length n*(n+1)/2 in which
1611 c     calcgh must store the lower triangle of the hessian at x.  start-
1612 c     ing at h(1), calcgh must store the hessian entries in the order
1613 c     (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ...
1614 c        the value printed (by itsum) in the column labelled stppar
1615 c     is the levenberg-marquardt used in computing the current step.
1616 c     zero means a full newton step.  if the special case described in
1617 c     ref. 1 is detected, then stppar is negated.  the value printed
1618 c     in the column labelled npreldf is zero if the current hessian
1619 c     is not positive definite.
1620 c        it sometimes proves worthwhile to let d be determined from the
1621 c     diagonal of the hessian matrix by setting iv(dtype) = 1 and
1622 c     v(dinit) = 0.  the following iv and v components are relevant...
1623 c
1624 c iv(dtol)..... iv(59) gives the starting subscript in v of the dtol
1625 c             array used when d is updated.  (iv(dtol) can be
1626 c             initialized by calling humsl with iv(1) = 13.)
1627 c iv(dtype).... iv(16) tells how the scale vector d should be chosen.
1628 c             iv(dtype) .le. 0 means that d should not be updated, and
1629 c             iv(dtype) .ge. 1 means that d should be updated as
1630 c             described below with v(dfac).  default = 0.
1631 c v(dfac)..... v(41) and the dtol and d0 arrays (see v(dtinit) and
1632 c             v(d0init)) are used in updating the scale vector d when
1633 c             iv(dtype) .gt. 0.  (d is initialized according to
1634 c             v(dinit), described in sumsl.)  let
1635 c                  d1(i) = max(sqrt(abs(h(i,i))), v(dfac)*d(i)),
1636 c             where h(i,i) is the i-th diagonal element of the current
1637 c             hessian.  if iv(dtype) = 1, then d(i) is set to d1(i)
1638 c             unless d1(i) .lt. dtol(i), in which case d(i) is set to
1639 c                  max(d0(i), dtol(i)).
1640 c             if iv(dtype) .ge. 2, then d is updated during the first
1641 c             iteration as for iv(dtype) = 1 (after any initialization
1642 c             due to v(dinit)) and is left unchanged thereafter.
1643 c             default = 0.6.
1644 c v(dtinit)... v(39), if positive, is the value to which all components
1645 c             of the dtol array (see v(dfac)) are initialized.  if
1646 c             v(dtinit) = 0, then it is assumed that the caller has
1647 c             stored dtol in v starting at v(iv(dtol)).
1648 c             default = 10**-6.
1649 c v(d0init)... v(40), if positive, is the value to which all components
1650 c             of the d0 vector (see v(dfac)) are initialized.  if
1651 c             v(dfac) = 0, then it is assumed that the caller has
1652 c             stored d0 in v starting at v(iv(dtol)+n).  default = 1.0.
1653 c
1654 c  ***  reference  ***
1655 c
1656 c 1. gay, d.m. (1981), computing optimal locally constrained steps,
1657 c         siam j. sci. statist. comput. 2, pp. 186-197.
1658 c.
1659 c  ***  general  ***
1660 c
1661 c     coded by david m. gay (winter 1980).  revised sept. 1982.
1662 c     this subroutine was written in connection with research supported
1663 c     in part by the national science foundation under grants
1664 c     mcs-7600324 and mcs-7906671.
1665 c
1666 c----------------------------  declarations  ---------------------------
1667 c
1668       external deflt, humit
1669 c
1670 c deflt... provides default input values for iv and v.
1671 c humit... reverse-communication routine that does humsl algorithm.
1672 c
1673       integer g1, h1, iv1, lh, nf
1674       double precision f
1675 c
1676 c  ***  subscripts for iv   ***
1677 c
1678       integer g, h, nextv, nfcall, nfgcal, toobig, vneed
1679 c
1680 c/6
1681 c     data nextv/47/, nfcall/6/, nfgcal/7/, g/28/, h/56/, toobig/2/,
1682 c    1     vneed/4/
1683 c/7
1684       parameter (nextv=47, nfcall=6, nfgcal=7, g=28, h=56, toobig=2,
1685      1           vneed=4)
1686 c/
1687 c
1688 c+++++++++++++++++++++++++++++++  body  ++++++++++++++++++++++++++++++++
1689 c
1690       lh = n * (n + 1) / 2
1691       if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1692       if (iv(1) .eq. 12 .or. iv(1) .eq. 13)
1693      1     iv(vneed) = iv(vneed) + n*(n+3)/2
1694       iv1 = iv(1)
1695       if (iv1 .eq. 14) go to 10
1696       if (iv1 .gt. 2 .and. iv1 .lt. 12) go to 10
1697       g1 = 1
1698       h1 = 1
1699       if (iv1 .eq. 12) iv(1) = 13
1700       go to 20
1701 c
1702  10   g1 = iv(g)
1703       h1 = iv(h)
1704 c
1705  20   call humit(d, f, v(g1), v(h1), iv, lh, liv, lv, n, v, x)
1706       if (iv(1) - 2) 30, 40, 50
1707 c
1708  30   nf = iv(nfcall)
1709       call calcf(n, x, nf, f, uiparm, urparm, ufparm)
1710       if (nf .le. 0) iv(toobig) = 1
1711       go to 20
1712 c
1713  40   call calcgh(n, x, iv(nfgcal), v(g1), v(h1), uiparm, urparm,
1714      1            ufparm)
1715       go to 20
1716 c
1717  50   if (iv(1) .ne. 14) go to 999
1718 c
1719 c  ***  storage allocation
1720 c
1721       iv(g) = iv(nextv)
1722       iv(h) = iv(g) + n
1723       iv(nextv) = iv(h) + n*(n+1)/2
1724       if (iv1 .ne. 13) go to 10
1725 c
1726  999  return
1727 c  ***  last card of humsl follows  ***
1728       end
1729       subroutine humit(d, fx, g, h, iv, lh, liv, lv, n, v, x)
1730 c
1731 c  ***  carry out humsl (unconstrained minimization) iterations, using
1732 c  ***  hessian matrix provided by the caller.
1733 c
1734 c  ***  parameter declarations  ***
1735 c
1736       integer lh, liv, lv, n
1737       integer iv(liv)
1738       double precision d(n), fx, g(n), h(lh), v(lv), x(n)
1739 c
1740 c--------------------------  parameter usage  --------------------------
1741 c
1742 c d.... scale vector.
1743 c fx... function value.
1744 c g.... gradient vector.
1745 c h.... lower triangle of the hessian, stored rowwise.
1746 c iv... integer value array.
1747 c lh... length of h = p*(p+1)/2.
1748 c liv.. length of iv (at least 60).
1749 c lv... length of v (at least 78 + n*(n+21)/2).
1750 c n.... number of variables (components in x and g).
1751 c v.... floating-point value array.
1752 c x.... parameter vector.
1753 c
1754 c  ***  discussion  ***
1755 c
1756 c        parameters iv, n, v, and x are the same as the corresponding
1757 c     ones to humsl (which see), except that v can be shorter (since
1758 c     the part of v that humsl uses for storing g and h is not needed).
1759 c     moreover, compared with humsl, iv(1) may have the two additional
1760 c     output values 1 and 2, which are explained below, as is the use
1761 c     of iv(toobig) and iv(nfgcal).  the value iv(g), which is an
1762 c     output value from humsl, is not referenced by humit or the
1763 c     subroutines it calls.
1764 c
1765 c iv(1) = 1 means the caller should set fx to f(x), the function value
1766 c             at x, and call humit again, having changed none of the
1767 c             other parameters.  an exception occurs if f(x) cannot be
1768 c             computed (e.g. if overflow would occur), which may happen
1769 c             because of an oversized step.  in this case the caller
1770 c             should set iv(toobig) = iv(2) to 1, which will cause
1771 c             humit to ignore fx and try a smaller step.  the para-
1772 c             meter nf that humsl passes to calcf (for possible use by
1773 c             calcgh) is a copy of iv(nfcall) = iv(6).
1774 c iv(1) = 2 means the caller should set g to g(x), the gradient of f at
1775 c             x, and h to the lower triangle of h(x), the hessian of f
1776 c             at x, and call humit again, having changed none of the
1777 c             other parameters except perhaps the scale vector d.
1778 c                  the parameter nf that humsl passes to calcg is
1779 c             iv(nfgcal) = iv(7).  if g(x) and h(x) cannot be evaluated,
1780 c             then the caller may set iv(nfgcal) to 0, in which case
1781 c             humit will return with iv(1) = 65.
1782 c                  note -- humit overwrites h with the lower triangle
1783 c             of  diag(d)**-1 * h(x) * diag(d)**-1.
1784 c.
1785 c  ***  general  ***
1786 c
1787 c     coded by david m. gay (winter 1980).  revised sept. 1982.
1788 c     this subroutine was written in connection with research supported
1789 c     in part by the national science foundation under grants
1790 c     mcs-7600324 and mcs-7906671.
1791 c
1792 c        (see sumsl and humsl for references.)
1793 c
1794 c+++++++++++++++++++++++++++  declarations  ++++++++++++++++++++++++++++
1795 c
1796 c  ***  local variables  ***
1797 c
1798       integer dg1, dummy, i, j, k, l, lstgst, nn1o2, step1,
1799      1        temp1, w1, x01
1800       double precision t
1801 c
1802 c     ***  constants  ***
1803 c
1804       double precision one, onep2, zero
1805 c
1806 c  ***  no intrinsic functions  ***
1807 c
1808 c  ***  external functions and subroutines  ***
1809 c
1810       external assst, deflt, dotprd, dupdu, gqtst, itsum, parck,
1811      1         reldst, slvmul, stopx, vaxpy, vcopy, vscopy, v2norm
1812       logical stopx
1813       double precision dotprd, reldst, v2norm
1814 c
1815 c assst.... assesses candidate step.
1816 c deflt.... provides default iv and v input values.
1817 c dotprd... returns inner product of two vectors.
1818 c dupdu.... updates scale vector d.
1819 c gqtst.... computes optimally locally constrained step.
1820 c itsum.... prints iteration summary and info on initial and final x.
1821 c parck.... checks validity of input iv and v values.
1822 c reldst... computes v(reldx) = relative step size.
1823 c slvmul... multiplies symmetric matrix times vector, given the lower
1824 c             triangle of the matrix.
1825 c stopx.... returns .true. if the break key has been pressed.
1826 c vaxpy.... computes scalar times one vector plus another.
1827 c vcopy.... copies one vector to another.
1828 c vscopy... sets all elements of a vector to a scalar.
1829 c v2norm... returns the 2-norm of a vector.
1830 c
1831 c  ***  subscripts for iv and v  ***
1832 c
1833       integer cnvcod, dg, dgnorm, dinit, dstnrm, dtinit, dtol,
1834      1        dtype, d0init, f, f0, fdif, gtstep, incfac, irc, kagqt,
1835      2        lmat, lmax0, lmaxs, mode, model, mxfcal, mxiter, nextv,
1836      3        nfcall, nfgcal, ngcall, niter, preduc, radfac, radinc,
1837      4        radius, rad0, reldx, restor, step, stglim, stlstg, stppar,
1838      5        toobig, tuner4, tuner5, vneed, w, xirc, x0
1839 c
1840 c  ***  iv subscript values  ***
1841 c
1842 c/6
1843 c     data cnvcod/55/, dg/37/, dtol/59/, dtype/16/, irc/29/, kagqt/33/,
1844 c    1     lmat/42/, mode/35/, model/5/, mxfcal/17/, mxiter/18/,
1845 c    2     nextv/47/, nfcall/6/, nfgcal/7/, ngcall/30/, niter/31/,
1846 c    3     radinc/8/, restor/9/, step/40/, stglim/11/, stlstg/41/,
1847 c    4     toobig/2/, vneed/4/, w/34/, xirc/13/, x0/43/
1848 c/7
1849       parameter (cnvcod=55, dg=37, dtol=59, dtype=16, irc=29, kagqt=33,
1850      1           lmat=42, mode=35, model=5, mxfcal=17, mxiter=18,
1851      2           nextv=47, nfcall=6, nfgcal=7, ngcall=30, niter=31,
1852      3           radinc=8, restor=9, step=40, stglim=11, stlstg=41,
1853      4           toobig=2, vneed=4, w=34, xirc=13, x0=43)
1854 c/
1855 c
1856 c  ***  v subscript values  ***
1857 c
1858 c/6
1859 c     data dgnorm/1/, dinit/38/, dstnrm/2/, dtinit/39/, d0init/40/,
1860 c    1     f/10/, f0/13/, fdif/11/, gtstep/4/, incfac/23/, lmax0/35/,
1861 c    2     lmaxs/36/, preduc/7/, radfac/16/, radius/8/, rad0/9/,
1862 c    3     reldx/17/, stppar/5/, tuner4/29/, tuner5/30/
1863 c/7
1864       parameter (dgnorm=1, dinit=38, dstnrm=2, dtinit=39, d0init=40,
1865      1           f=10, f0=13, fdif=11, gtstep=4, incfac=23, lmax0=35,
1866      2           lmaxs=36, preduc=7, radfac=16, radius=8, rad0=9,
1867      3           reldx=17, stppar=5, tuner4=29, tuner5=30)
1868 c/
1869 c
1870 c/6
1871 c     data one/1.d+0/, onep2/1.2d+0/, zero/0.d+0/
1872 c/7
1873       parameter (one=1.d+0, onep2=1.2d+0, zero=0.d+0)
1874 c/
1875 c
1876 c+++++++++++++++++++++++++++++++  body  ++++++++++++++++++++++++++++++++
1877 c
1878       i = iv(1)
1879       if (i .eq. 1) go to 30
1880       if (i .eq. 2) go to 40
1881 c
1882 c  ***  check validity of iv and v input values  ***
1883 c
1884       if (iv(1) .eq. 0) call deflt(2, iv, liv, lv, v)
1885       if (iv(1) .eq. 12 .or. iv(1) .eq. 13)
1886      1     iv(vneed) = iv(vneed) + n*(n+21)/2 + 7
1887       call parck(2, d, iv, liv, lv, n, v)
1888       i = iv(1) - 2
1889       if (i .gt. 12) go to 999
1890       nn1o2 = n * (n + 1) / 2
1891       if (lh .ge. nn1o2) go to (210,210,210,210,210,210,160,120,160,
1892      1                          10,10,20), i
1893          iv(1) = 66
1894          go to 350
1895 c
1896 c  ***  storage allocation  ***
1897 c
1898  10   iv(dtol) = iv(lmat) + nn1o2
1899       iv(x0) = iv(dtol) + 2*n
1900       iv(step) = iv(x0) + n
1901       iv(stlstg) = iv(step) + n
1902       iv(dg) = iv(stlstg) + n
1903       iv(w) = iv(dg) + n
1904       iv(nextv) = iv(w) + 4*n + 7
1905       if (iv(1) .ne. 13) go to 20
1906          iv(1) = 14
1907          go to 999
1908 c
1909 c  ***  initialization  ***
1910 c
1911  20   iv(niter) = 0
1912       iv(nfcall) = 1
1913       iv(ngcall) = 1
1914       iv(nfgcal) = 1
1915       iv(mode) = -1
1916       iv(model) = 1
1917       iv(stglim) = 1
1918       iv(toobig) = 0
1919       iv(cnvcod) = 0
1920       iv(radinc) = 0
1921       v(rad0) = zero
1922       v(stppar) = zero
1923       if (v(dinit) .ge. zero) call vscopy(n, d, v(dinit))
1924       k = iv(dtol)
1925       if (v(dtinit) .gt. zero) call vscopy(n, v(k), v(dtinit))
1926       k = k + n
1927       if (v(d0init) .gt. zero) call vscopy(n, v(k), v(d0init))
1928       iv(1) = 1
1929       go to 999
1930 c
1931  30   v(f) = fx
1932       if (iv(mode) .ge. 0) go to 210
1933       iv(1) = 2
1934       if (iv(toobig) .eq. 0) go to 999
1935          iv(1) = 63
1936          go to 350
1937 c
1938 c  ***  make sure gradient could be computed  ***
1939 c
1940  40   if (iv(nfgcal) .ne. 0) go to 50
1941          iv(1) = 65
1942          go to 350
1943 c
1944 c  ***  update the scale vector d  ***
1945 c
1946  50   dg1 = iv(dg)
1947       if (iv(dtype) .le. 0) go to 70
1948       k = dg1
1949       j = 0
1950       do 60 i = 1, n
1951          j = j + i
1952          v(k) = h(j)
1953          k = k + 1
1954  60      continue
1955       call dupdu(d, v(dg1), iv, liv, lv, n, v)
1956 c
1957 c  ***  compute scaled gradient and its norm  ***
1958 c
1959  70   dg1 = iv(dg)
1960       k = dg1
1961       do 80 i = 1, n
1962          v(k) = g(i) / d(i)
1963          k = k + 1
1964  80      continue
1965       v(dgnorm) = v2norm(n, v(dg1))
1966 c
1967 c  ***  compute scaled hessian  ***
1968 c
1969       k = 1
1970       do 100 i = 1, n
1971          t = one / d(i)
1972          do 90 j = 1, i
1973               h(k) = t * h(k) / d(j)
1974               k = k + 1
1975  90           continue
1976  100     continue
1977 c
1978       if (iv(cnvcod) .ne. 0) go to 340
1979       if (iv(mode) .eq. 0) go to 300
1980 c
1981 c  ***  allow first step to have scaled 2-norm at most v(lmax0)  ***
1982 c
1983       v(radius) = v(lmax0)
1984 c
1985       iv(mode) = 0
1986 c
1987 c
1988 c-----------------------------  main loop  -----------------------------
1989 c
1990 c
1991 c  ***  print iteration summary, check iteration limit  ***
1992 c
1993  110  call itsum(d, g, iv, liv, lv, n, v, x)
1994  120  k = iv(niter)
1995       if (k .lt. iv(mxiter)) go to 130
1996          iv(1) = 10
1997          go to 350
1998 c
1999  130  iv(niter) = k + 1
2000 c
2001 c  ***  initialize for start of next iteration  ***
2002 c
2003       dg1 = iv(dg)
2004       x01 = iv(x0)
2005       v(f0) = v(f)
2006       iv(irc) = 4
2007       iv(kagqt) = -1
2008 c
2009 c     ***  copy x to x0  ***
2010 c
2011       call vcopy(n, v(x01), x)
2012 c
2013 c  ***  update radius  ***
2014 c
2015       if (k .eq. 0) go to 150
2016       step1 = iv(step)
2017       k = step1
2018       do 140 i = 1, n
2019          v(k) = d(i) * v(k)
2020          k = k + 1
2021  140     continue
2022       v(radius) = v(radfac) * v2norm(n, v(step1))
2023 c
2024 c  ***  check stopx and function evaluation limit  ***
2025 c
2026 C AL 4/30/95
2027       dummy=iv(nfcall)
2028  150  if (.not. stopx(dummy)) go to 170
2029          iv(1) = 11
2030          go to 180
2031 c
2032 c     ***  come here when restarting after func. eval. limit or stopx.
2033 c
2034  160  if (v(f) .ge. v(f0)) go to 170
2035          v(radfac) = one
2036          k = iv(niter)
2037          go to 130
2038 c
2039  170  if (iv(nfcall) .lt. iv(mxfcal)) go to 190
2040          iv(1) = 9
2041  180     if (v(f) .ge. v(f0)) go to 350
2042 c
2043 c        ***  in case of stopx or function evaluation limit with
2044 c        ***  improved v(f), evaluate the gradient at x.
2045 c
2046               iv(cnvcod) = iv(1)
2047               go to 290
2048 c
2049 c. . . . . . . . . . . . .  compute candidate step  . . . . . . . . . .
2050 c
2051  190  step1 = iv(step)
2052       dg1 = iv(dg)
2053       l = iv(lmat)
2054       w1 = iv(w)
2055       call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1))
2056       if (iv(irc) .eq. 6) go to 210
2057 c
2058 c  ***  check whether evaluating f(x0 + step) looks worthwhile  ***
2059 c
2060       if (v(dstnrm) .le. zero) go to 210
2061       if (iv(irc) .ne. 5) go to 200
2062       if (v(radfac) .le. one) go to 200
2063       if (v(preduc) .le. onep2 * v(fdif)) go to 210
2064 c
2065 c  ***  compute f(x0 + step)  ***
2066 c
2067  200  x01 = iv(x0)
2068       step1 = iv(step)
2069       call vaxpy(n, x, one, v(step1), v(x01))
2070       iv(nfcall) = iv(nfcall) + 1
2071       iv(1) = 1
2072       iv(toobig) = 0
2073       go to 999
2074 c
2075 c. . . . . . . . . . . . .  assess candidate step  . . . . . . . . . . .
2076 c
2077  210  x01 = iv(x0)
2078       v(reldx) = reldst(n, d, x, v(x01))
2079       call assst(iv, liv, lv, v)
2080       step1 = iv(step)
2081       lstgst = iv(stlstg)
2082       if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
2083       if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
2084       if (iv(restor) .ne. 3) go to 220
2085          call vcopy(n, v(step1), v(lstgst))
2086          call vaxpy(n, x, one, v(step1), v(x01))
2087          v(reldx) = reldst(n, d, x, v(x01))
2088 c
2089  220  k = iv(irc)
2090       go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k
2091 c
2092 c     ***  recompute step with new radius  ***
2093 c
2094  230     v(radius) = v(radfac) * v(dstnrm)
2095          go to 150
2096 c
2097 c  ***  compute step of length v(lmaxs) for singular convergence test.
2098 c
2099  240  v(radius) = v(lmaxs)
2100       go to 190
2101 c
2102 c  ***  convergence or false convergence  ***
2103 c
2104  250  iv(cnvcod) = k - 4
2105       if (v(f) .ge. v(f0)) go to 340
2106          if (iv(xirc) .eq. 14) go to 340
2107               iv(xirc) = 14
2108 c
2109 c. . . . . . . . . . . .  process acceptable step  . . . . . . . . . . .
2110 c
2111  260  if (iv(irc) .ne. 3) go to 290
2112          temp1 = lstgst
2113 c
2114 c     ***  prepare for gradient tests  ***
2115 c     ***  set  temp1 = hessian * step + g(x0)
2116 c     ***             = diag(d) * (h * step + g(x0))
2117 c
2118 c        use x0 vector as temporary.
2119          k = x01
2120          do 270 i = 1, n
2121               v(k) = d(i) * v(step1)
2122               k = k + 1
2123               step1 = step1 + 1
2124  270          continue
2125          call slvmul(n, v(temp1), h, v(x01))
2126          do 280 i = 1, n
2127               v(temp1) = d(i) * v(temp1) + g(i)
2128               temp1 = temp1 + 1
2129  280          continue
2130 c
2131 c  ***  compute gradient and hessian  ***
2132 c
2133  290  iv(ngcall) = iv(ngcall) + 1
2134       iv(1) = 2
2135       go to 999
2136 c
2137  300  iv(1) = 2
2138       if (iv(irc) .ne. 3) go to 110
2139 c
2140 c  ***  set v(radfac) by gradient tests  ***
2141 c
2142       temp1 = iv(stlstg)
2143       step1 = iv(step)
2144 c
2145 c     ***  set  temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x)))  ***
2146 c
2147       k = temp1
2148       do 310 i = 1, n
2149          v(k) = (v(k) - g(i)) / d(i)
2150          k = k + 1
2151  310     continue
2152 c
2153 c     ***  do gradient tests  ***
2154 c
2155       if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320
2156            if (dotprd(n, g, v(step1))
2157      1               .ge. v(gtstep) * v(tuner5))  go to 110
2158  320            v(radfac) = v(incfac)
2159                 go to 110
2160 c
2161 c. . . . . . . . . . . . . .  misc. details  . . . . . . . . . . . . . .
2162 c
2163 c  ***  bad parameters to assess  ***
2164 c
2165  330  iv(1) = 64
2166       go to 350
2167 c
2168 c  ***  print summary of final iteration and other requested items  ***
2169 c
2170  340  iv(1) = iv(cnvcod)
2171       iv(cnvcod) = 0
2172  350  call itsum(d, g, iv, liv, lv, n, v, x)
2173 c
2174  999  return
2175 c
2176 c  ***  last card of humit follows  ***
2177       end
2178       subroutine dupdu(d, hdiag, iv, liv, lv, n, v)
2179 c
2180 c  ***  update scale vector d for humsl  ***
2181 c
2182 c  ***  parameter declarations  ***
2183 c
2184       integer liv, lv, n
2185       integer iv(liv)
2186       double precision d(n), hdiag(n), v(lv)
2187 c
2188 c  ***  local variables  ***
2189 c
2190       integer dtoli, d0i, i
2191       double precision t, vdfac
2192 c
2193 c  ***  intrinsic functions  ***
2194 c/+
2195       double precision dabs, dmax1, dsqrt
2196 c/
2197 c  ***  subscripts for iv and v  ***
2198 c
2199       integer dfac, dtol, dtype, niter
2200 c/6
2201 c     data dfac/41/, dtol/59/, dtype/16/, niter/31/
2202 c/7
2203       parameter (dfac=41, dtol=59, dtype=16, niter=31)
2204 c/
2205 c
2206 c-------------------------------  body  --------------------------------
2207 c
2208       i = iv(dtype)
2209       if (i .eq. 1) go to 10
2210          if (iv(niter) .gt. 0) go to 999
2211 c
2212  10   dtoli = iv(dtol)
2213       d0i = dtoli + n
2214       vdfac = v(dfac)
2215       do 20 i = 1, n
2216          t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i))
2217          if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i))
2218          d(i) = t
2219          dtoli = dtoli + 1
2220          d0i = d0i + 1
2221  20      continue
2222 c
2223  999  return
2224 c  ***  last card of dupdu follows  ***
2225       end
2226       subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w)
2227 c
2228 c  *** compute goldfeld-quandt-trotter step by more-hebden technique ***
2229 c  ***  (nl2sol version 2.2), modified a la more and sorensen  ***
2230 c
2231 c  ***  parameter declarations  ***
2232 c
2233       integer ka, p
2234 cal   double precision d(p), dig(p), dihdi(1), l(1), v(21), step(p),
2235 cal  1                 w(1)
2236       double precision d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),
2237      1    v(21), step(p),w(4*p+7)
2238 c     dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7)
2239 c
2240 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2241 c
2242 c  ***  purpose  ***
2243 c
2244 c        given the (compactly stored) lower triangle of a scaled
2245 c     hessian (approximation) and a nonzero scaled gradient vector,
2246 c     this subroutine computes a goldfeld-quandt-trotter step of
2247 c     approximate length v(radius) by the more-hebden technique.  in
2248 c     other words, step is computed to (approximately) minimize
2249 c     psi(step) = (g**t)*step + 0.5*(step**t)*h*step  such that the
2250 c     2-norm of d*step is at most (approximately) v(radius), where
2251 c     g  is the gradient,  h  is the hessian, and  d  is a diagonal
2252 c     scale matrix whose diagonal is stored in the parameter d.
2253 c     (gqtst assumes  dig = d**-1 * g  and  dihdi = d**-1 * h * d**-1.)
2254 c
2255 c  ***  parameter description  ***
2256 c
2257 c     d (in)  = the scale vector, i.e. the diagonal of the scale
2258 c              matrix  d  mentioned above under purpose.
2259 c   dig (in)  = the scaled gradient vector, d**-1 * g.  if g = 0, then
2260 c              step = 0  and  v(stppar) = 0  are returned.
2261 c dihdi (in)  = lower triangle of the scaled hessian (approximation),
2262 c              i.e., d**-1 * h * d**-1, stored compactly by rows., i.e.,
2263 c              in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc.
2264 c    ka (i/o) = the number of hebden iterations (so far) taken to deter-
2265 c              mine step.  ka .lt. 0 on input means this is the first
2266 c              attempt to determine step (for the present dig and dihdi)
2267 c              -- ka is initialized to 0 in this case.  output with
2268 c              ka = 0  (or v(stppar) = 0)  means  step = -(h**-1)*g.
2269 c     l (i/o) = workspace of length p*(p+1)/2 for cholesky factors.
2270 c     p (in)  = number of parameters -- the hessian is a  p x p  matrix.
2271 c  step (i/o) = the step computed.
2272 c     v (i/o) contains various constants and variables described below.
2273 c     w (i/o) = workspace of length 4*p + 6.
2274 c
2275 c  ***  entries in v  ***
2276 c
2277 c v(dgnorm) (i/o) = 2-norm of (d**-1)*g.
2278 c v(dstnrm) (output) = 2-norm of d*step.
2279 c v(dst0)   (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or
2280 c             overestimate of smallest eigenvalue of (d**-1)*h*(d**-1).
2281 c v(epslon) (in)  = max. rel. error allowed for psi(step).  for the
2282 c             step returned, psi(step) will exceed its optimal value
2283 c             by less than -v(epslon)*psi(step).  suggested value = 0.1.
2284 c v(gtstep) (out) = inner product between g and step.
2285 c v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step)  (for pos. def.
2286 c             h only -- v(nreduc) is set to zero otherwise).
2287 c v(phmnfc) (in)  = tol. (together with v(phmxfc)) for accepting step
2288 c             (more*s sigma).  the error v(dstnrm) - v(radius) must lie
2289 c             between v(phmnfc)*v(radius) and v(phmxfc)*v(radius).
2290 c v(phmxfc) (in)  (see v(phmnfc).)
2291 c             suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5.
2292 c v(preduc) (out) = psi(step) = predicted obj. func. reduction for step.
2293 c v(radius) (in)  = radius of current (scaled) trust region.
2294 c v(rad0)   (i/o) = value of v(radius) from previous call.
2295 c v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha
2296 c             described below under algorithm notes.  if h + alpha*d**2
2297 c             (see algorithm notes) is (nearly) singular, however,
2298 c             then v(stppar) = -alpha.
2299 c
2300 c  ***  usage notes  ***
2301 c
2302 c     if it is desired to recompute step using a different value of
2303 c     v(radius), then this routine may be restarted by calling it
2304 c     with all parameters unchanged except v(radius).  (this explains
2305 c     why step and w are listed as i/o).  on an initial call (one with
2306 c     ka .lt. 0), step and w need not be initialized and only compo-
2307 c     nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and
2308 c     v(rad0) of v must be initialized.
2309 c
2310 c  ***  algorithm notes  ***
2311 c
2312 c        the desired g-q-t step (ref. 2, 3, 4, 6) satisfies
2313 c     (h + alpha*d**2)*step = -g  for some nonnegative alpha such that
2314 c     h + alpha*d**2 is positive semidefinite.  alpha and step are
2315 c     computed by a scheme analogous to the one described in ref. 5.
2316 c     estimates of the smallest and largest eigenvalues of the hessian
2317 c     are obtained from the gerschgorin circle theorem enhanced by a
2318 c     simple form of the scaling described in ref. 7.  cases in which
2319 c     h + alpha*d**2 is nearly (or exactly) singular are handled by
2320 c     the technique discussed in ref. 2.  in these cases, a step of
2321 c     (exact) length v(radius) is returned for which psi(step) exceeds
2322 c     its optimal value by less than -v(epslon)*psi(step).  the test
2323 c     suggested in ref. 6 for detecting the special case is performed
2324 c     once two matrix factorizations have been done -- doing so sooner
2325 c     seems to degrade the performance of optimization routines that
2326 c     call this routine.
2327 c
2328 c  ***  functions and subroutines called  ***
2329 c
2330 c dotprd - returns inner product of two vectors.
2331 c litvmu - applies inverse-transpose of compact lower triang. matrix.
2332 c livmul - applies inverse of compact lower triang. matrix.
2333 c lsqrt  - finds cholesky factor (of compactly stored lower triang.).
2334 c lsvmin - returns approx. to min. sing. value of lower triang. matrix.
2335 c rmdcon - returns machine-dependent constants.
2336 c v2norm - returns 2-norm of a vector.
2337 c
2338 c  ***  references  ***
2339 c
2340 c 1.  dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive
2341 c             nonlinear least-squares algorithm, acm trans. math.
2342 c             software, vol. 7, no. 3.
2343 c 2.  gay, d.m. (1981), computing optimal locally constrained steps,
2344 c             siam j. sci. statist. computing, vol. 2, no. 2, pp.
2345 c             186-197.
2346 c 3.  goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966),
2347 c             maximization by quadratic hill-climbing, econometrica 34,
2348 c             pp. 541-551.
2349 c 4.  hebden, m.d. (1973), an algorithm for minimization using exact
2350 c             second derivatives, report t.p. 515, theoretical physics
2351 c             div., a.e.r.e. harwell, oxon., england.
2352 c 5.  more, j.j. (1978), the levenberg-marquardt algorithm, implemen-
2353 c             tation and theory, pp.105-116 of springer lecture notes
2354 c             in mathematics no. 630, edited by g.a. watson, springer-
2355 c             verlag, berlin and new york.
2356 c 6.  more, j.j., and sorensen, d.c. (1981), computing a trust region
2357 c             step, technical report anl-81-83, argonne national lab.
2358 c 7.  varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15,
2359 c             pp. 719-729.
2360 c
2361 c  ***  general  ***
2362 c
2363 c     coded by david m. gay.
2364 c     this subroutine was written in connection with research
2365 c     supported by the national science foundation under grants
2366 c     mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
2367 c     mcs-7906671.
2368 c
2369 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2370 c
2371 c  ***  local variables  ***
2372 c
2373       logical restrt
2374       integer dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,
2375      1        j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x
2376       double precision alphak, aki, akk, delta, dst, eps, gtsta, lk,
2377      1                 oldphi, phi, phimax, phimin, psifac, rad, radsq,
2378      2                 root, si, sk, sw, t, twopsi, t1, t2, uk, wi
2379 c
2380 c     ***  constants  ***
2381       double precision big, dgxfac, epsfac, four, half, kappa, negone,
2382      1                 one, p001, six, three, two, zero
2383 c
2384 c  ***  intrinsic functions  ***
2385 c/+
2386       double precision dabs, dmax1, dmin1, dsqrt
2387 c/
2388 c  ***  external functions and subroutines  ***
2389 c
2390       external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm
2391       double precision dotprd, lsvmin, rmdcon, v2norm
2392 c
2393 c  ***  subscripts for v  ***
2394 c
2395       integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc,
2396      1        phmnfc, phmxfc, preduc, radius, rad0
2397 c/6
2398 c     data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/,
2399 c    1     nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/,
2400 c    2     rad0/9/, stppar/5/
2401 c/7
2402       parameter (dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,
2403      1           nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,
2404      2           rad0=9, stppar=5)
2405 c/
2406 c
2407 c/6
2408 c     data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/,
2409 c    1     kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/,
2410 c    2     six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/
2411 c/7
2412       parameter (epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,
2413      1     kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,
2414      2     six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0)
2415       save dgxfac
2416 c/
2417       data big/0.d+0/, dgxfac/0.d+0/
2418 c
2419 c  ***  body  ***
2420 c
2421 c     ***  store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx).
2422       dggdmx = p + 1
2423 c     ***  store gerschgorin over- and underestimates of the largest
2424 c     ***  and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax)
2425 c     ***  and w(emin) respectively.
2426       emax = dggdmx + 1
2427       emin = emax + 1
2428 c     ***  for use in recomputing step, the final values of lk, uk, dst,
2429 c     ***  and the inverse derivative of more*s phi at 0 (for pos. def.
2430 c     ***  h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin)
2431 c     ***  respectively.
2432       lk0 = emin + 1
2433       phipin = lk0 + 1
2434       uk0 = phipin + 1
2435       dstsav = uk0 + 1
2436 c     ***  store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p).
2437       diag0 = dstsav
2438       diag = diag0 + 1
2439 c     ***  store -d*step in w(q),...,w(q0+p).
2440       q0 = diag0 + p
2441       q = q0 + 1
2442 c     ***  allocate storage for scratch vector x  ***
2443       x = q + p
2444       rad = v(radius)
2445       radsq = rad**2
2446 c     ***  phitol = max. error allowed in dst = v(dstnrm) = 2-norm of
2447 c     ***  d*step.
2448       phimax = v(phmxfc) * rad
2449       phimin = v(phmnfc) * rad
2450       psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) *
2451      1                       (kappa + one)  +  kappa  +  two) * rad**2)
2452 c     ***  oldphi is used to detect limits of numerical accuracy.  if
2453 c     ***  we recompute step and it does not change, then we accept it.
2454       oldphi = zero
2455       eps = v(epslon)
2456       irc = 0
2457       restrt = .false.
2458       kalim = ka + 50
2459 c
2460 c  ***  start or restart, depending on ka  ***
2461 c
2462       if (ka .ge. 0) go to 290
2463 c
2464 c  ***  fresh start  ***
2465 c
2466       k = 0
2467       uk = negone
2468       ka = 0
2469       kalim = 50
2470       v(dgnorm) = v2norm(p, dig)
2471       v(nreduc) = zero
2472       v(dst0) = zero
2473       kamin = 3
2474       if (v(dgnorm) .eq. zero) kamin = 0
2475 c
2476 c     ***  store diag(dihdi) in w(diag0+1),...,w(diag0+p)  ***
2477 c
2478       j = 0
2479       do 10 i = 1, p
2480          j = j + i
2481          k1 = diag0 + i
2482          w(k1) = dihdi(j)
2483  10      continue
2484 c
2485 c     ***  determine w(dggdmx), the largest element of dihdi  ***
2486 c
2487       t1 = zero
2488       j = p * (p + 1) / 2
2489       do 20 i = 1, j
2490          t = dabs(dihdi(i))
2491          if (t1 .lt. t) t1 = t
2492  20      continue
2493       w(dggdmx) = t1
2494 c
2495 c  ***  try alpha = 0  ***
2496 c
2497  30   call lsqrt(1, p, l, dihdi, irc)
2498       if (irc .eq. 0) go to 50
2499 c        ***  indef. h -- underestimate smallest eigenvalue, use this
2500 c        ***  estimate to initialize lower bound lk on alpha.
2501          j = irc*(irc+1)/2
2502          t = l(j)
2503          l(j) = one
2504          do 40 i = 1, irc
2505  40           w(i) = zero
2506          w(irc) = one
2507          call litvmu(irc, w, l, w)
2508          t1 = v2norm(irc, w)
2509          lk = -t / t1 / t1
2510          v(dst0) = -lk
2511          if (restrt) go to 210
2512          go to 70
2513 c
2514 c     ***  positive definite h -- compute unmodified newton step.  ***
2515  50   lk = zero
2516       t = lsvmin(p, l, w(q), w(q))
2517       if (t .ge. one) go to 60
2518          if (big .le. zero) big = rmdcon(6)
2519          if (v(dgnorm) .ge. t*t*big) go to 70
2520  60   call livmul(p, w(q), l, dig)
2521       gtsta = dotprd(p, w(q), w(q))
2522       v(nreduc) = half * gtsta
2523       call litvmu(p, w(q), l, w(q))
2524       dst = v2norm(p, w(q))
2525       v(dst0) = dst
2526       phi = dst - rad
2527       if (phi .le. phimax) go to 260
2528       if (restrt) go to 210
2529 c
2530 c  ***  prepare to compute gerschgorin estimates of largest (and
2531 c  ***  smallest) eigenvalues.  ***
2532 c
2533  70   k = 0
2534       do 100 i = 1, p
2535          wi = zero
2536          if (i .eq. 1) go to 90
2537          im1 = i - 1
2538          do 80 j = 1, im1
2539               k = k + 1
2540               t = dabs(dihdi(k))
2541               wi = wi + t
2542               w(j) = w(j) + t
2543  80           continue
2544  90      w(i) = wi
2545          k = k + 1
2546  100     continue
2547 c
2548 c  ***  (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1)  ***
2549 c
2550       k = 1
2551       t1 = w(diag) - w(1)
2552       if (p .le. 1) go to 120
2553       do 110 i = 2, p
2554          j = diag0 + i
2555          t = w(j) - w(i)
2556          if (t .ge. t1) go to 110
2557               t1 = t
2558               k = i
2559  110     continue
2560 c
2561  120  sk = w(k)
2562       j = diag0 + k
2563       akk = w(j)
2564       k1 = k*(k-1)/2 + 1
2565       inc = 1
2566       t = zero
2567       do 150 i = 1, p
2568          if (i .eq. k) go to 130
2569          aki = dabs(dihdi(k1))
2570          si = w(i)
2571          j = diag0 + i
2572          t1 = half * (akk - w(j) + si - aki)
2573          t1 = t1 + dsqrt(t1*t1 + sk*aki)
2574          if (t .lt. t1) t = t1
2575          if (i .lt. k) go to 140
2576  130     inc = i
2577  140     k1 = k1 + inc
2578  150     continue
2579 c
2580       w(emin) = akk - t
2581       uk = v(dgnorm)/rad - w(emin)
2582       if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2583       if (uk .le. zero) uk = p001
2584 c
2585 c  ***  compute gerschgorin (over-)estimate of largest eigenvalue  ***
2586 c
2587       k = 1
2588       t1 = w(diag) + w(1)
2589       if (p .le. 1) go to 170
2590       do 160 i = 2, p
2591          j = diag0 + i
2592          t = w(j) + w(i)
2593          if (t .le. t1) go to 160
2594               t1 = t
2595               k = i
2596  160     continue
2597 c
2598  170  sk = w(k)
2599       j = diag0 + k
2600       akk = w(j)
2601       k1 = k*(k-1)/2 + 1
2602       inc = 1
2603       t = zero
2604       do 200 i = 1, p
2605          if (i .eq. k) go to 180
2606          aki = dabs(dihdi(k1))
2607          si = w(i)
2608          j = diag0 + i
2609          t1 = half * (w(j) + si - aki - akk)
2610          t1 = t1 + dsqrt(t1*t1 + sk*aki)
2611          if (t .lt. t1) t = t1
2612          if (i .lt. k) go to 190
2613  180     inc = i
2614  190     k1 = k1 + inc
2615  200     continue
2616 c
2617       w(emax) = akk + t
2618       lk = dmax1(lk, v(dgnorm)/rad - w(emax))
2619 c
2620 c     ***  alphak = current value of alpha (see alg. notes above).  we
2621 c     ***  use more*s scheme for initializing it.
2622       alphak = dabs(v(stppar)) * v(rad0)/rad
2623 c
2624       if (irc .ne. 0) go to 210
2625 c
2626 c  ***  compute l0 for positive definite h  ***
2627 c
2628       call livmul(p, w, l, w(q))
2629       t = v2norm(p, w)
2630       w(phipin) = dst / t / t
2631       lk = dmax1(lk, phi*w(phipin))
2632 c
2633 c  ***  safeguard alphak and add alphak*i to (d**-1)*h*(d**-1)  ***
2634 c
2635  210  ka = ka + 1
2636       if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk)
2637      1                      alphak = uk * dmax1(p001, dsqrt(lk/uk))
2638       if (alphak .le. zero) alphak = half * uk
2639       if (alphak .le. zero) alphak = uk
2640       k = 0
2641       do 220 i = 1, p
2642          k = k + i
2643          j = diag0 + i
2644          dihdi(k) = w(j) + alphak
2645  220     continue
2646 c
2647 c  ***  try computing cholesky decomposition  ***
2648 c
2649       call lsqrt(1, p, l, dihdi, irc)
2650       if (irc .eq. 0) go to 240
2651 c
2652 c  ***  (d**-1)*h*(d**-1) + alphak*i  is indefinite -- overestimate
2653 c  ***  smallest eigenvalue for use in updating lk  ***
2654 c
2655       j = (irc*(irc+1))/2
2656       t = l(j)
2657       l(j) = one
2658       do 230 i = 1, irc
2659  230     w(i) = zero
2660       w(irc) = one
2661       call litvmu(irc, w, l, w)
2662       t1 = v2norm(irc, w)
2663       lk = alphak - t/t1/t1
2664       v(dst0) = -lk
2665       go to 210
2666 c
2667 c  ***  alphak makes (d**-1)*h*(d**-1) positive definite.
2668 c  ***  compute q = -d*step, check for convergence.  ***
2669 c
2670  240  call livmul(p, w(q), l, dig)
2671       gtsta = dotprd(p, w(q), w(q))
2672       call litvmu(p, w(q), l, w(q))
2673       dst = v2norm(p, w(q))
2674       phi = dst - rad
2675       if (phi .le. phimax .and. phi .ge. phimin) go to 270
2676       if (phi .eq. oldphi) go to 270
2677       oldphi = phi
2678       if (phi .lt. zero) go to 330
2679 c
2680 c  ***  unacceptable alphak -- update lk, uk, alphak  ***
2681 c
2682  250  if (ka .ge. kalim) go to 270
2683 c     ***  the following dmin1 is necessary because of restarts  ***
2684       if (phi .lt. zero) uk = dmin1(uk, alphak)
2685 c     *** kamin = 0 only iff the gradient vanishes  ***
2686       if (kamin .eq. 0) go to 210
2687       call livmul(p, w, l, w(q))
2688       t1 = v2norm(p, w)
2689       alphak = alphak  +  (phi/t1) * (dst/t1) * (dst/rad)
2690       lk = dmax1(lk, alphak)
2691       go to 210
2692 c
2693 c  ***  acceptable step on first try  ***
2694 c
2695  260  alphak = zero
2696 c
2697 c  ***  successful step in general.  compute step = -(d**-1)*q  ***
2698 c
2699  270  do 280 i = 1, p
2700          j = q0 + i
2701          step(i) = -w(j)/d(i)
2702  280     continue
2703       v(gtstep) = -gtsta
2704       v(preduc) = half * (dabs(alphak)*dst*dst + gtsta)
2705       go to 410
2706 c
2707 c
2708 c  ***  restart with new radius  ***
2709 c
2710  290  if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310
2711 c
2712 c     ***  prepare to return newton step  ***
2713 c
2714          restrt = .true.
2715          ka = ka + 1
2716          k = 0
2717          do 300 i = 1, p
2718               k = k + i
2719               j = diag0 + i
2720               dihdi(k) = w(j)
2721  300          continue
2722          uk = negone
2723          go to 30
2724 c
2725  310  kamin = ka + 3
2726       if (v(dgnorm) .eq. zero) kamin = 0
2727       if (ka .eq. 0) go to 50
2728 c
2729       dst = w(dstsav)
2730       alphak = dabs(v(stppar))
2731       phi = dst - rad
2732       t = v(dgnorm)/rad
2733       uk = t - w(emin)
2734       if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2735       if (uk .le. zero) uk = p001
2736       if (rad .gt. v(rad0)) go to 320
2737 c
2738 c        ***  smaller radius  ***
2739          lk = zero
2740          if (alphak .gt. zero) lk = w(lk0)
2741          lk = dmax1(lk, t - w(emax))
2742          if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2743          go to 250
2744 c
2745 c     ***  bigger radius  ***
2746  320  if (alphak .gt. zero) uk = dmin1(uk, w(uk0))
2747       lk = dmax1(zero, -v(dst0), t - w(emax))
2748       if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2749       go to 250
2750 c
2751 c  ***  decide whether to check for special case... in practice (from
2752 c  ***  the standpoint of the calling optimization code) it seems best
2753 c  ***  not to check until a few iterations have failed -- hence the
2754 c  ***  test on kamin below.
2755 c
2756  330  delta = alphak + dmin1(zero, v(dst0))
2757       twopsi = alphak*dst*dst + gtsta
2758       if (ka .ge. kamin) go to 340
2759 c     *** if the test in ref. 2 is satisfied, fall through to handle
2760 c     *** the special case (as soon as the more-sorensen test detects
2761 c     *** it).
2762       if (delta .ge. psifac*twopsi) go to 370
2763 c
2764 c  ***  check for the special case of  h + alpha*d**2  (nearly)
2765 c  ***  singular.  use one step of inverse power method with start
2766 c  ***  from lsvmin to obtain approximate eigenvector corresponding
2767 c  ***  to smallest eigenvalue of (d**-1)*h*(d**-1).  lsvmin returns
2768 c  ***  x and w with  l*w = x.
2769 c
2770  340  t = lsvmin(p, l, w(x), w)
2771 c
2772 c     ***  normalize w  ***
2773       do 350 i = 1, p
2774  350     w(i) = t*w(i)
2775 c     ***  complete current inv. power iter. -- replace w by (l**-t)*w.
2776       call litvmu(p, w, l, w)
2777       t2 = one/v2norm(p, w)
2778       do 360 i = 1, p
2779  360     w(i) = t2*w(i)
2780       t = t2 * t
2781 c
2782 c  ***  now w is the desired approximate (unit) eigenvector and
2783 c  ***  t*x = ((d**-1)*h*(d**-1) + alphak*i)*w.
2784 c
2785       sw = dotprd(p, w(q), w)
2786       t1 = (rad + dst) * (rad - dst)
2787       root = dsqrt(sw*sw + t1)
2788       if (sw .lt. zero) root = -root
2789       si = t1 / (sw + root)
2790 c
2791 c  ***  the actual test for the special case...
2792 c
2793       if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380
2794 c
2795 c  ***  update upper bound on smallest eigenvalue (when not positive)
2796 c  ***  (as recommended by more and sorensen) and continue...
2797 c
2798       if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak)
2799       lk = dmax1(lk, -v(dst0))
2800 c
2801 c  ***  check whether we can hope to detect the special case in
2802 c  ***  the available arithmetic.  accept step as it is if not.
2803 c
2804 c     ***  if not yet available, obtain machine dependent value dgxfac.
2805  370  if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3)
2806 c
2807       if (delta .gt. dgxfac*w(dggdmx)) go to 250
2808          go to 270
2809 c
2810 c  ***  special case detected... negate alphak to indicate special case
2811 c
2812  380  alphak = -alphak
2813       v(preduc) = half * twopsi
2814 c
2815 c  ***  accept current step if adding si*w would lead to a
2816 c  ***  further relative reduction in psi of less than v(epslon)/3.
2817 c
2818       t1 = zero
2819       t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w)))
2820       if (t .lt. eps*twopsi/six) go to 390
2821          v(preduc) = v(preduc) + t
2822          dst = rad
2823          t1 = -si
2824  390  do 400 i = 1, p
2825          j = q0 + i
2826          w(j) = t1*w(i) - w(j)
2827          step(i) = w(j) / d(i)
2828  400     continue
2829       v(gtstep) = dotprd(p, dig, w(q))
2830 c
2831 c  ***  save values for use in a possible restart  ***
2832 c
2833  410  v(dstnrm) = dst
2834       v(stppar) = alphak
2835       w(lk0) = lk
2836       w(uk0) = uk
2837       v(rad0) = rad
2838       w(dstsav) = dst
2839 c
2840 c     ***  restore diagonal of dihdi  ***
2841 c
2842       j = 0
2843       do 420 i = 1, p
2844          j = j + i
2845          k = diag0 + i
2846          dihdi(j) = w(k)
2847  420     continue
2848 c
2849  999  return
2850 c
2851 c  ***  last card of gqtst follows  ***
2852       end
2853       subroutine lsqrt(n1, n, l, a, irc)
2854 c
2855 c  ***  compute rows n1 through n of the cholesky factor  l  of
2856 c  ***  a = l*(l**t),  where  l  and the lower triangle of  a  are both
2857 c  ***  stored compactly by rows (and may occupy the same storage).
2858 c  ***  irc = 0 means all went well.  irc = j means the leading
2859 c  ***  principal  j x j  submatrix of  a  is not positive definite --
2860 c  ***  and  l(j*(j+1)/2)  contains the (nonpos.) reduced j-th diagonal.
2861 c
2862 c  ***  parameters  ***
2863 c
2864       integer n1, n, irc
2865 cal   double precision l(1), a(1)
2866       double precision l(n*(n+1)/2), a(n*(n+1)/2)
2867 c     dimension l(n*(n+1)/2), a(n*(n+1)/2)
2868 c
2869 c  ***  local variables  ***
2870 c
2871       integer i, ij, ik, im1, i0, j, jk, jm1, j0, k
2872       double precision t, td, zero
2873 c
2874 c  ***  intrinsic functions  ***
2875 c/+
2876       double precision dsqrt
2877 c/
2878 c/6
2879 c     data zero/0.d+0/
2880 c/7
2881       parameter (zero=0.d+0)
2882 c/
2883 c
2884 c  ***  body  ***
2885 c
2886       i0 = n1 * (n1 - 1) / 2
2887       do 50 i = n1, n
2888          td = zero
2889          if (i .eq. 1) go to 40
2890          j0 = 0
2891          im1 = i - 1
2892          do 30 j = 1, im1
2893               t = zero
2894               if (j .eq. 1) go to 20
2895               jm1 = j - 1
2896               do 10 k = 1, jm1
2897                    ik = i0 + k
2898                    jk = j0 + k
2899                    t = t + l(ik)*l(jk)
2900  10                continue
2901  20           ij = i0 + j
2902               j0 = j0 + j
2903               t = (a(ij) - t) / l(j0)
2904               l(ij) = t
2905               td = td + t*t
2906  30           continue
2907  40      i0 = i0 + i
2908          t = a(i0) - td
2909          if (t .le. zero) go to 60
2910          l(i0) = dsqrt(t)
2911  50      continue
2912 c
2913       irc = 0
2914       go to 999
2915 c
2916  60   l(i0) = t
2917       irc = i
2918 c
2919  999  return
2920 c
2921 c  ***  last card of lsqrt  ***
2922       end
2923       double precision function lsvmin(p, l, x, y)
2924 c
2925 c  ***  estimate smallest sing. value of packed lower triang. matrix l
2926 c
2927 c  ***  parameter declarations  ***
2928 c
2929       integer p
2930 cal   double precision l(1), x(p), y(p)
2931       double precision l(p*(p+1)/2), x(p), y(p)
2932 c     dimension l(p*(p+1)/2)
2933 c
2934 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2935 c
2936 c  ***  purpose  ***
2937 c
2938 c     this function returns a good over-estimate of the smallest
2939 c     singular value of the packed lower triangular matrix l.
2940 c
2941 c  ***  parameter description  ***
2942 c
2943 c  p (in)  = the order of l.  l is a  p x p  lower triangular matrix.
2944 c  l (in)  = array holding the elements of  l  in row order, i.e.
2945 c             l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc.
2946 c  x (out) if lsvmin returns a positive value, then x is a normalized
2947 c             approximate left singular vector corresponding to the
2948 c             smallest singular value.  this approximation may be very
2949 c             crude.  if lsvmin returns zero, then some components of x
2950 c             are zero and the rest retain their input values.
2951 c  y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an
2952 c             unnormalized approximate right singular vector correspond-
2953 c             ing to the smallest singular value.  this approximation
2954 c             may be crude.  if lsvmin returns zero, then y retains its
2955 c             input value.  the caller may pass the same vector for x
2956 c             and y (nonstandard fortran usage), in which case y over-
2957 c             writes x (for nonzero lsvmin returns).
2958 c
2959 c  ***  algorithm notes  ***
2960 c
2961 c     the algorithm is based on (1), with the additional provision that
2962 c     lsvmin = 0 is returned if the smallest diagonal element of l
2963 c     (in magnitude) is not more than the unit roundoff times the
2964 c     largest.  the algorithm uses a random number generator proposed
2965 c     in (4), which passes the spectral test with flying colors -- see
2966 c     (2) and (3).
2967 c
2968 c  ***  subroutines and functions called  ***
2969 c
2970 c        v2norm - function, returns the 2-norm of a vector.
2971 c
2972 c  ***  references  ***
2973 c
2974 c     (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977),
2975 c         an estimate for the condition number of a matrix, report
2976 c         tm-310, applied math. div., argonne national laboratory.
2977 c
2978 c     (2) hoaglin, d.c. (1976), theoretical properties of congruential
2979 c         random-number generators --  an empirical view,
2980 c         memorandum ns-340, dept. of statistics, harvard univ.
2981 c
2982 c     (3) knuth, d.e. (1969), the art of computer programming, vol. 2
2983 c         (seminumerical algorithms), addison-wesley, reading, mass.
2984 c
2985 c     (4) smith, c.s. (1971), multiplicative pseudo-random number
2986 c         generators with prime modulus, j. assoc. comput. mach. 18,
2987 c         pp. 586-593.
2988 c
2989 c  ***  history  ***
2990 c
2991 c     designed and coded by david m. gay (winter 1977/summer 1978).
2992 c
2993 c  ***  general  ***
2994 c
2995 c     this subroutine was written in connection with research
2996 c     supported by the national science foundation under grants
2997 c     mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989.
2998 c
2999 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3000 c
3001 c  ***  local variables  ***
3002 c
3003       integer i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1
3004       double precision b, sminus, splus, t, xminus, xplus
3005 c
3006 c  ***  constants  ***
3007 c
3008       double precision half, one, r9973, zero
3009 c
3010 c  ***  intrinsic functions  ***
3011 c/+
3012       integer mod
3013       real float
3014       double precision dabs
3015 c/
3016 c  ***  external functions and subroutines  ***
3017 c
3018       external dotprd, v2norm, vaxpy
3019       double precision dotprd, v2norm
3020 c
3021 c/6
3022 c     data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/
3023 c/7
3024       parameter (half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0)
3025 c/
3026 c
3027 c  ***  body  ***
3028 c
3029       ix = 2
3030       pm1 = p - 1
3031 c
3032 c  ***  first check whether to return lsvmin = 0 and initialize x  ***
3033 c
3034       ii = 0
3035       j0 = p*pm1/2
3036       jj = j0 + p
3037       if (l(jj) .eq. zero) go to 110
3038       ix = mod(3432*ix, 9973)
3039       b = half*(one + float(ix)/r9973)
3040       xplus = b / l(jj)
3041       x(p) = xplus
3042       if (p .le. 1) go to 60
3043       do 10 i = 1, pm1
3044          ii = ii + i
3045          if (l(ii) .eq. zero) go to 110
3046          ji = j0 + i
3047          x(i) = xplus * l(ji)
3048  10      continue
3049 c
3050 c  ***  solve (l**t)*x = b, where the components of b have randomly
3051 c  ***  chosen magnitudes in (.5,1) with signs chosen to make x large.
3052 c
3053 c     do j = p-1 to 1 by -1...
3054       do 50 jjj = 1, pm1
3055          j = p - jjj
3056 c       ***  determine x(j) in this iteration. note for i = 1,2,...,j
3057 c       ***  that x(i) holds the current partial sum for row i.
3058          ix = mod(3432*ix, 9973)
3059          b = half*(one + float(ix)/r9973)
3060          xplus = (b - x(j))
3061          xminus = (-b - x(j))
3062          splus = dabs(xplus)
3063          sminus = dabs(xminus)
3064          jm1 = j - 1
3065          j0 = j*jm1/2
3066          jj = j0 + j
3067          xplus = xplus/l(jj)
3068          xminus = xminus/l(jj)
3069          if (jm1 .eq. 0) go to 30
3070          do 20 i = 1, jm1
3071               ji = j0 + i
3072               splus = splus + dabs(x(i) + l(ji)*xplus)
3073               sminus = sminus + dabs(x(i) + l(ji)*xminus)
3074  20           continue
3075  30      if (sminus .gt. splus) xplus = xminus
3076          x(j) = xplus
3077 c       ***  update partial sums  ***
3078          if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x)
3079  50      continue
3080 c
3081 c  ***  normalize x  ***
3082 c
3083  60   t = one/v2norm(p, x)
3084       do 70 i = 1, p
3085  70      x(i) = t*x(i)
3086 c
3087 c  ***  solve l*y = x and return lsvmin = 1/twonorm(y)  ***
3088 c
3089       do 100 j = 1, p
3090          jm1 = j - 1
3091          j0 = j*jm1/2
3092          jj = j0 + j
3093          t = zero
3094          if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y)
3095          y(j) = (x(j) - t) / l(jj)
3096  100     continue
3097 c
3098       lsvmin = one/v2norm(p, y)
3099       go to 999
3100 c
3101  110  lsvmin = zero
3102  999  return
3103 c  ***  last card of lsvmin follows  ***
3104       end
3105       subroutine slvmul(p, y, s, x)
3106 c
3107 c  ***  set  y = s * x,  s = p x p symmetric matrix.  ***
3108 c  ***  lower triangle of  s  stored rowwise.         ***
3109 c
3110 c  ***  parameter declarations  ***
3111 c
3112       integer p
3113 cal   double precision s(1), x(p), y(p)
3114       double precision s(p*(p+1)/2), x(p), y(p)
3115 c     dimension s(p*(p+1)/2)
3116 c
3117 c  ***  local variables  ***
3118 c
3119       integer i, im1, j, k
3120       double precision xi
3121 c
3122 c  ***  no intrinsic functions  ***
3123 c
3124 c  ***  external function  ***
3125 c
3126       external dotprd
3127       double precision dotprd
3128 c
3129 c-----------------------------------------------------------------------
3130 c
3131       j = 1
3132       do 10 i = 1, p
3133          y(i) = dotprd(i, s(j), x)
3134          j = j + i
3135  10      continue
3136 c
3137       if (p .le. 1) go to 999
3138       j = 1
3139       do 40 i = 2, p
3140          xi = x(i)
3141          im1 = i - 1
3142          j = j + 1
3143          do 30 k = 1, im1
3144               y(k) = y(k) + s(j)*xi
3145               j = j + 1
3146  30           continue
3147  40      continue
3148 c
3149  999  return
3150 c  ***  last card of slvmul follows  ***
3151       end