source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/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) = mi
 575       iv(lastv) = m
 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 c      dummy=iv(niter)
2028 c      write (iout,*) "dummy",dummy
2029  150  if (.not. stopx(iv(niter))) go to 170
2030          iv(1) = 11
2031          go to 180
2032 c
2033 c     ***  come here when restarting after func. eval. limit or stopx.
2034 c
2035  160  if (v(f) .ge. v(f0)) go to 170
2036          v(radfac) = one
2037          k = iv(niter)
2038          go to 130
2039 c
2040  170  if (iv(nfcall) .lt. iv(mxfcal)) go to 190
2041          iv(1) = 9
2042  180     if (v(f) .ge. v(f0)) go to 350
2043 c
2044 c        ***  in case of stopx or function evaluation limit with
2045 c        ***  improved v(f), evaluate the gradient at x.
2046 c
2047               iv(cnvcod) = iv(1)
2048               go to 290
2049 c
2050 c. . . . . . . . . . . . .  compute candidate step  . . . . . . . . . .
2051 c
2052  190  step1 = iv(step)
2053       dg1 = iv(dg)
2054       l = iv(lmat)
2055       w1 = iv(w)
2056       call gqtst(d, v(dg1), h, iv(kagqt), v(l), n, v(step1), v, v(w1))
2057       if (iv(irc) .eq. 6) go to 210
2058 c
2059 c  ***  check whether evaluating f(x0 + step) looks worthwhile  ***
2060 c
2061       if (v(dstnrm) .le. zero) go to 210
2062       if (iv(irc) .ne. 5) go to 200
2063       if (v(radfac) .le. one) go to 200
2064       if (v(preduc) .le. onep2 * v(fdif)) go to 210
2065 c
2066 c  ***  compute f(x0 + step)  ***
2067 c
2068  200  x01 = iv(x0)
2069       step1 = iv(step)
2070       call vaxpy(n, x, one, v(step1), v(x01))
2071       iv(nfcall) = iv(nfcall) + 1
2072       iv(1) = 1
2073       iv(toobig) = 0
2074       go to 999
2075 c
2076 c. . . . . . . . . . . . .  assess candidate step  . . . . . . . . . . .
2077 c
2078  210  x01 = iv(x0)
2079       v(reldx) = reldst(n, d, x, v(x01))
2080       call assst(iv, liv, lv, v)
2081       step1 = iv(step)
2082       lstgst = iv(stlstg)
2083       if (iv(restor) .eq. 1) call vcopy(n, x, v(x01))
2084       if (iv(restor) .eq. 2) call vcopy(n, v(lstgst), v(step1))
2085       if (iv(restor) .ne. 3) go to 220
2086          call vcopy(n, v(step1), v(lstgst))
2087          call vaxpy(n, x, one, v(step1), v(x01))
2088          v(reldx) = reldst(n, d, x, v(x01))
2089 c
2090  220  k = iv(irc)
2091       go to (230,260,260,260,230,240,250,250,250,250,250,250,330,300), k
2092 c
2093 c     ***  recompute step with new radius  ***
2094 c
2095  230     v(radius) = v(radfac) * v(dstnrm)
2096          go to 150
2097 c
2098 c  ***  compute step of length v(lmaxs) for singular convergence test.
2099 c
2100  240  v(radius) = v(lmaxs)
2101       go to 190
2102 c
2103 c  ***  convergence or false convergence  ***
2104 c
2105  250  iv(cnvcod) = k - 4
2106       if (v(f) .ge. v(f0)) go to 340
2107          if (iv(xirc) .eq. 14) go to 340
2108               iv(xirc) = 14
2109 c
2110 c. . . . . . . . . . . .  process acceptable step  . . . . . . . . . . .
2111 c
2112  260  if (iv(irc) .ne. 3) go to 290
2113          temp1 = lstgst
2114 c
2115 c     ***  prepare for gradient tests  ***
2116 c     ***  set  temp1 = hessian * step + g(x0)
2117 c     ***             = diag(d) * (h * step + g(x0))
2118 c
2119 c        use x0 vector as temporary.
2120          k = x01
2121          do 270 i = 1, n
2122               v(k) = d(i) * v(step1)
2123               k = k + 1
2124               step1 = step1 + 1
2125  270          continue
2126          call slvmul(n, v(temp1), h, v(x01))
2127          do 280 i = 1, n
2128               v(temp1) = d(i) * v(temp1) + g(i)
2129               temp1 = temp1 + 1
2130  280          continue
2131 c
2132 c  ***  compute gradient and hessian  ***
2133 c
2134  290  iv(ngcall) = iv(ngcall) + 1
2135       iv(1) = 2
2136       go to 999
2137 c
2138  300  iv(1) = 2
2139       if (iv(irc) .ne. 3) go to 110
2140 c
2141 c  ***  set v(radfac) by gradient tests  ***
2142 c
2143       temp1 = iv(stlstg)
2144       step1 = iv(step)
2145 c
2146 c     ***  set  temp1 = diag(d)**-1 * (hessian*step + (g(x0)-g(x)))  ***
2147 c
2148       k = temp1
2149       do 310 i = 1, n
2150          v(k) = (v(k) - g(i)) / d(i)
2151          k = k + 1
2152  310     continue
2153 c
2154 c     ***  do gradient tests  ***
2155 c
2156       if (v2norm(n, v(temp1)) .le. v(dgnorm) * v(tuner4)) go to 320
2157            if (dotprd(n, g, v(step1))
2158      1               .ge. v(gtstep) * v(tuner5))  go to 110
2159  320            v(radfac) = v(incfac)
2160                 go to 110
2161 c
2162 c. . . . . . . . . . . . . .  misc. details  . . . . . . . . . . . . . .
2163 c
2164 c  ***  bad parameters to assess  ***
2165 c
2166  330  iv(1) = 64
2167       go to 350
2168 c
2169 c  ***  print summary of final iteration and other requested items  ***
2170 c
2171  340  iv(1) = iv(cnvcod)
2172       iv(cnvcod) = 0
2173  350  call itsum(d, g, iv, liv, lv, n, v, x)
2174 c
2175  999  return
2176 c
2177 c  ***  last card of humit follows  ***
2178       end
2179       subroutine dupdu(d, hdiag, iv, liv, lv, n, v)
2180 c
2181 c  ***  update scale vector d for humsl  ***
2182 c
2183 c  ***  parameter declarations  ***
2184 c
2185       integer liv, lv, n
2186       integer iv(liv)
2187       double precision d(n), hdiag(n), v(lv)
2188 c
2189 c  ***  local variables  ***
2190 c
2191       integer dtoli, d0i, i
2192       double precision t, vdfac
2193 c
2194 c  ***  intrinsic functions  ***
2195 c/+
2196       double precision dabs, dmax1, dsqrt
2197 c/
2198 c  ***  subscripts for iv and v  ***
2199 c
2200       integer dfac, dtol, dtype, niter
2201 c/6
2202 c     data dfac/41/, dtol/59/, dtype/16/, niter/31/
2203 c/7
2204       parameter (dfac=41, dtol=59, dtype=16, niter=31)
2205 c/
2206 c
2207 c-------------------------------  body  --------------------------------
2208 c
2209       i = iv(dtype)
2210       if (i .eq. 1) go to 10
2211          if (iv(niter) .gt. 0) go to 999
2212 c
2213  10   dtoli = iv(dtol)
2214       d0i = dtoli + n
2215       vdfac = v(dfac)
2216       do 20 i = 1, n
2217          t = dmax1(dsqrt(dabs(hdiag(i))), vdfac*d(i))
2218          if (t .lt. v(dtoli)) t = dmax1(v(dtoli), v(d0i))
2219          d(i) = t
2220          dtoli = dtoli + 1
2221          d0i = d0i + 1
2222  20      continue
2223 c
2224  999  return
2225 c  ***  last card of dupdu follows  ***
2226       end
2227       subroutine gqtst(d, dig, dihdi, ka, l, p, step, v, w)
2228 c
2229 c  *** compute goldfeld-quandt-trotter step by more-hebden technique ***
2230 c  ***  (nl2sol version 2.2), modified a la more and sorensen  ***
2231 c
2232 c  ***  parameter declarations  ***
2233 c
2234       integer ka, p
2235 cal   double precision d(p), dig(p), dihdi(1), l(1), v(21), step(p),
2236 cal  1                 w(1)
2237       double precision d(p), dig(p), dihdi(p*(p+1)/2), l(p*(p+1)/2),
2238      1    v(21), step(p),w(4*p+7)
2239 c     dimension dihdi(p*(p+1)/2), l(p*(p+1)/2), w(4*p+7)
2240 c
2241 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2242 c
2243 c  ***  purpose  ***
2244 c
2245 c        given the (compactly stored) lower triangle of a scaled
2246 c     hessian (approximation) and a nonzero scaled gradient vector,
2247 c     this subroutine computes a goldfeld-quandt-trotter step of
2248 c     approximate length v(radius) by the more-hebden technique.  in
2249 c     other words, step is computed to (approximately) minimize
2250 c     psi(step) = (g**t)*step + 0.5*(step**t)*h*step  such that the
2251 c     2-norm of d*step is at most (approximately) v(radius), where
2252 c     g  is the gradient,  h  is the hessian, and  d  is a diagonal
2253 c     scale matrix whose diagonal is stored in the parameter d.
2254 c     (gqtst assumes  dig = d**-1 * g  and  dihdi = d**-1 * h * d**-1.)
2255 c
2256 c  ***  parameter description  ***
2257 c
2258 c     d (in)  = the scale vector, i.e. the diagonal of the scale
2259 c              matrix  d  mentioned above under purpose.
2260 c   dig (in)  = the scaled gradient vector, d**-1 * g.  if g = 0, then
2261 c              step = 0  and  v(stppar) = 0  are returned.
2262 c dihdi (in)  = lower triangle of the scaled hessian (approximation),
2263 c              i.e., d**-1 * h * d**-1, stored compactly by rows., i.e.,
2264 c              in the order (1,1), (2,1), (2,2), (3,1), (3,2), etc.
2265 c    ka (i/o) = the number of hebden iterations (so far) taken to deter-
2266 c              mine step.  ka .lt. 0 on input means this is the first
2267 c              attempt to determine step (for the present dig and dihdi)
2268 c              -- ka is initialized to 0 in this case.  output with
2269 c              ka = 0  (or v(stppar) = 0)  means  step = -(h**-1)*g.
2270 c     l (i/o) = workspace of length p*(p+1)/2 for cholesky factors.
2271 c     p (in)  = number of parameters -- the hessian is a  p x p  matrix.
2272 c  step (i/o) = the step computed.
2273 c     v (i/o) contains various constants and variables described below.
2274 c     w (i/o) = workspace of length 4*p + 6.
2275 c
2276 c  ***  entries in v  ***
2277 c
2278 c v(dgnorm) (i/o) = 2-norm of (d**-1)*g.
2279 c v(dstnrm) (output) = 2-norm of d*step.
2280 c v(dst0)   (i/o) = 2-norm of d*(h**-1)*g (for pos. def. h only), or
2281 c             overestimate of smallest eigenvalue of (d**-1)*h*(d**-1).
2282 c v(epslon) (in)  = max. rel. error allowed for psi(step).  for the
2283 c             step returned, psi(step) will exceed its optimal value
2284 c             by less than -v(epslon)*psi(step).  suggested value = 0.1.
2285 c v(gtstep) (out) = inner product between g and step.
2286 c v(nreduc) (out) = psi(-(h**-1)*g) = psi(newton step)  (for pos. def.
2287 c             h only -- v(nreduc) is set to zero otherwise).
2288 c v(phmnfc) (in)  = tol. (together with v(phmxfc)) for accepting step
2289 c             (more*s sigma).  the error v(dstnrm) - v(radius) must lie
2290 c             between v(phmnfc)*v(radius) and v(phmxfc)*v(radius).
2291 c v(phmxfc) (in)  (see v(phmnfc).)
2292 c             suggested values -- v(phmnfc) = -0.25, v(phmxfc) = 0.5.
2293 c v(preduc) (out) = psi(step) = predicted obj. func. reduction for step.
2294 c v(radius) (in)  = radius of current (scaled) trust region.
2295 c v(rad0)   (i/o) = value of v(radius) from previous call.
2296 c v(stppar) (i/o) is normally the marquardt parameter, i.e. the alpha
2297 c             described below under algorithm notes.  if h + alpha*d**2
2298 c             (see algorithm notes) is (nearly) singular, however,
2299 c             then v(stppar) = -alpha.
2300 c
2301 c  ***  usage notes  ***
2302 c
2303 c     if it is desired to recompute step using a different value of
2304 c     v(radius), then this routine may be restarted by calling it
2305 c     with all parameters unchanged except v(radius).  (this explains
2306 c     why step and w are listed as i/o).  on an initial call (one with
2307 c     ka .lt. 0), step and w need not be initialized and only compo-
2308 c     nents v(epslon), v(stppar), v(phmnfc), v(phmxfc), v(radius), and
2309 c     v(rad0) of v must be initialized.
2310 c
2311 c  ***  algorithm notes  ***
2312 c
2313 c        the desired g-q-t step (ref. 2, 3, 4, 6) satisfies
2314 c     (h + alpha*d**2)*step = -g  for some nonnegative alpha such that
2315 c     h + alpha*d**2 is positive semidefinite.  alpha and step are
2316 c     computed by a scheme analogous to the one described in ref. 5.
2317 c     estimates of the smallest and largest eigenvalues of the hessian
2318 c     are obtained from the gerschgorin circle theorem enhanced by a
2319 c     simple form of the scaling described in ref. 7.  cases in which
2320 c     h + alpha*d**2 is nearly (or exactly) singular are handled by
2321 c     the technique discussed in ref. 2.  in these cases, a step of
2322 c     (exact) length v(radius) is returned for which psi(step) exceeds
2323 c     its optimal value by less than -v(epslon)*psi(step).  the test
2324 c     suggested in ref. 6 for detecting the special case is performed
2325 c     once two matrix factorizations have been done -- doing so sooner
2326 c     seems to degrade the performance of optimization routines that
2327 c     call this routine.
2328 c
2329 c  ***  functions and subroutines called  ***
2330 c
2331 c dotprd - returns inner product of two vectors.
2332 c litvmu - applies inverse-transpose of compact lower triang. matrix.
2333 c livmul - applies inverse of compact lower triang. matrix.
2334 c lsqrt  - finds cholesky factor (of compactly stored lower triang.).
2335 c lsvmin - returns approx. to min. sing. value of lower triang. matrix.
2336 c rmdcon - returns machine-dependent constants.
2337 c v2norm - returns 2-norm of a vector.
2338 c
2339 c  ***  references  ***
2340 c
2341 c 1.  dennis, j.e., gay, d.m., and welsch, r.e. (1981), an adaptive
2342 c             nonlinear least-squares algorithm, acm trans. math.
2343 c             software, vol. 7, no. 3.
2344 c 2.  gay, d.m. (1981), computing optimal locally constrained steps,
2345 c             siam j. sci. statist. computing, vol. 2, no. 2, pp.
2346 c             186-197.
2347 c 3.  goldfeld, s.m., quandt, r.e., and trotter, h.f. (1966),
2348 c             maximization by quadratic hill-climbing, econometrica 34,
2349 c             pp. 541-551.
2350 c 4.  hebden, m.d. (1973), an algorithm for minimization using exact
2351 c             second derivatives, report t.p. 515, theoretical physics
2352 c             div., a.e.r.e. harwell, oxon., england.
2353 c 5.  more, j.j. (1978), the levenberg-marquardt algorithm, implemen-
2354 c             tation and theory, pp.105-116 of springer lecture notes
2355 c             in mathematics no. 630, edited by g.a. watson, springer-
2356 c             verlag, berlin and new york.
2357 c 6.  more, j.j., and sorensen, d.c. (1981), computing a trust region
2358 c             step, technical report anl-81-83, argonne national lab.
2359 c 7.  varga, r.s. (1965), minimal gerschgorin sets, pacific j. math. 15,
2360 c             pp. 719-729.
2361 c
2362 c  ***  general  ***
2363 c
2364 c     coded by david m. gay.
2365 c     this subroutine was written in connection with research
2366 c     supported by the national science foundation under grants
2367 c     mcs-7600324, dcr75-10143, 76-14311dss, mcs76-11989, and
2368 c     mcs-7906671.
2369 c
2370 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2371 c
2372 c  ***  local variables  ***
2373 c
2374       logical restrt
2375       integer dggdmx, diag, diag0, dstsav, emax, emin, i, im1, inc, irc,
2376      1        j, k, kalim, kamin, k1, lk0, phipin, q, q0, uk0, x
2377       double precision alphak, aki, akk, delta, dst, eps, gtsta, lk,
2378      1                 oldphi, phi, phimax, phimin, psifac, rad, radsq,
2379      2                 root, si, sk, sw, t, twopsi, t1, t2, uk, wi
2380 c
2381 c     ***  constants  ***
2382       double precision big, dgxfac, epsfac, four, half, kappa, negone,
2383      1                 one, p001, six, three, two, zero
2384 c
2385 c  ***  intrinsic functions  ***
2386 c/+
2387       double precision dabs, dmax1, dmin1, dsqrt
2388 c/
2389 c  ***  external functions and subroutines  ***
2390 c
2391       external dotprd, litvmu, livmul, lsqrt, lsvmin, rmdcon, v2norm
2392       double precision dotprd, lsvmin, rmdcon, v2norm
2393 c
2394 c  ***  subscripts for v  ***
2395 c
2396       integer dgnorm, dstnrm, dst0, epslon, gtstep, stppar, nreduc,
2397      1        phmnfc, phmxfc, preduc, radius, rad0
2398 c/6
2399 c     data dgnorm/1/, dstnrm/2/, dst0/3/, epslon/19/, gtstep/4/,
2400 c    1     nreduc/6/, phmnfc/20/, phmxfc/21/, preduc/7/, radius/8/,
2401 c    2     rad0/9/, stppar/5/
2402 c/7
2403       parameter (dgnorm=1, dstnrm=2, dst0=3, epslon=19, gtstep=4,
2404      1           nreduc=6, phmnfc=20, phmxfc=21, preduc=7, radius=8,
2405      2           rad0=9, stppar=5)
2406 c/
2407 c
2408 c/6
2409 c     data epsfac/50.0d+0/, four/4.0d+0/, half/0.5d+0/,
2410 c    1     kappa/2.0d+0/, negone/-1.0d+0/, one/1.0d+0/, p001/1.0d-3/,
2411 c    2     six/6.0d+0/, three/3.0d+0/, two/2.0d+0/, zero/0.0d+0/
2412 c/7
2413       parameter (epsfac=50.0d+0, four=4.0d+0, half=0.5d+0,
2414      1     kappa=2.0d+0, negone=-1.0d+0, one=1.0d+0, p001=1.0d-3,
2415      2     six=6.0d+0, three=3.0d+0, two=2.0d+0, zero=0.0d+0)
2416       save dgxfac
2417 c/
2418       data big/0.d+0/, dgxfac/0.d+0/
2419 c
2420 c  ***  body  ***
2421 c
2422 c     ***  store largest abs. entry in (d**-1)*h*(d**-1) at w(dggdmx).
2423       dggdmx = p + 1
2424 c     ***  store gerschgorin over- and underestimates of the largest
2425 c     ***  and smallest eigenvalues of (d**-1)*h*(d**-1) at w(emax)
2426 c     ***  and w(emin) respectively.
2427       emax = dggdmx + 1
2428       emin = emax + 1
2429 c     ***  for use in recomputing step, the final values of lk, uk, dst,
2430 c     ***  and the inverse derivative of more*s phi at 0 (for pos. def.
2431 c     ***  h) are stored in w(lk0), w(uk0), w(dstsav), and w(phipin)
2432 c     ***  respectively.
2433       lk0 = emin + 1
2434       phipin = lk0 + 1
2435       uk0 = phipin + 1
2436       dstsav = uk0 + 1
2437 c     ***  store diag of (d**-1)*h*(d**-1) in w(diag),...,w(diag0+p).
2438       diag0 = dstsa
2439       diag = diag0 + 1
2440 c     ***  store -d*step in w(q),...,w(q0+p).
2441       q0 = diag0 + p
2442       q = q0 + 1
2443 c     ***  allocate storage for scratch vector x  ***
2444       x = q + p
2445       rad = v(radius)
2446       radsq = rad**2
2447 c     ***  phitol = max. error allowed in dst = v(dstnrm) = 2-norm of
2448 c     ***  d*step.
2449       phimax = v(phmxfc) * rad
2450       phimin = v(phmnfc) * rad
2451       psifac = two * v(epslon) / (three * (four * (v(phmnfc) + one) *
2452      1                       (kappa + one)  +  kappa  +  two) * rad**2)
2453 c     ***  oldphi is used to detect limits of numerical accuracy.  if
2454 c     ***  we recompute step and it does not change, then we accept it.
2455       oldphi = zero
2456       eps = v(epslon)
2457       irc = 0
2458       restrt = .false.
2459       kalim = ka + 50
2460 c
2461 c  ***  start or restart, depending on ka  ***
2462 c
2463       if (ka .ge. 0) go to 290
2464 c
2465 c  ***  fresh start  ***
2466 c
2467       k = 0
2468       uk = negone
2469       ka = 0
2470       kalim = 50
2471       v(dgnorm) = v2norm(p, dig)
2472       v(nreduc) = zero
2473       v(dst0) = zero
2474       kamin = 3
2475       if (v(dgnorm) .eq. zero) kamin = 0
2476 c
2477 c     ***  store diag(dihdi) in w(diag0+1),...,w(diag0+p)  ***
2478 c
2479       j = 0
2480       do 10 i = 1, p
2481          j = j + i
2482          k1 = diag0 + i
2483          w(k1) = dihdi(j)
2484  10      continue
2485 c
2486 c     ***  determine w(dggdmx), the largest element of dihdi  ***
2487 c
2488       t1 = zero
2489       j = p * (p + 1) / 2
2490       do 20 i = 1, j
2491          t = dabs(dihdi(i))
2492          if (t1 .lt. t) t1 = t
2493  20      continue
2494       w(dggdmx) = t1
2495 c
2496 c  ***  try alpha = 0  ***
2497 c
2498  30   call lsqrt(1, p, l, dihdi, irc)
2499       if (irc .eq. 0) go to 50
2500 c        ***  indef. h -- underestimate smallest eigenvalue, use this
2501 c        ***  estimate to initialize lower bound lk on alpha.
2502          j = irc*(irc+1)/2
2503          t = l(j)
2504          l(j) = one
2505          do 40 i = 1, irc
2506  40           w(i) = zero
2507          w(irc) = one
2508          call litvmu(irc, w, l, w)
2509          t1 = v2norm(irc, w)
2510          lk = -t / t1 / t1
2511          v(dst0) = -lk
2512          if (restrt) go to 210
2513          go to 70
2514 c
2515 c     ***  positive definite h -- compute unmodified newton step.  ***
2516  50   lk = zero
2517       t = lsvmin(p, l, w(q), w(q))
2518       if (t .ge. one) go to 60
2519          if (big .le. zero) big = rmdcon(6)
2520          if (v(dgnorm) .ge. t*t*big) go to 70
2521  60   call livmul(p, w(q), l, dig)
2522       gtsta = dotprd(p, w(q), w(q))
2523       v(nreduc) = half * gtsta
2524       call litvmu(p, w(q), l, w(q))
2525       dst = v2norm(p, w(q))
2526       v(dst0) = dst
2527       phi = dst - rad
2528       if (phi .le. phimax) go to 260
2529       if (restrt) go to 210
2530 c
2531 c  ***  prepare to compute gerschgorin estimates of largest (and
2532 c  ***  smallest) eigenvalues.  ***
2533 c
2534  70   k = 0
2535       do 100 i = 1, p
2536          wi = zero
2537          if (i .eq. 1) go to 90
2538          im1 = i - 1
2539          do 80 j = 1, im1
2540               k = k + 1
2541               t = dabs(dihdi(k))
2542               wi = wi + t
2543               w(j) = w(j) + t
2544  80           continue
2545  90      w(i) = wi
2546          k = k + 1
2547  100     continue
2548 c
2549 c  ***  (under-)estimate smallest eigenvalue of (d**-1)*h*(d**-1)  ***
2550 c
2551       k = 1
2552       t1 = w(diag) - w(1)
2553       if (p .le. 1) go to 120
2554       do 110 i = 2, p
2555          j = diag0 + i
2556          t = w(j) - w(i)
2557          if (t .ge. t1) go to 110
2558               t1 = t
2559               k = i
2560  110     continue
2561 c
2562  120  sk = w(k)
2563       j = diag0 + k
2564       akk = w(j)
2565       k1 = k*(k-1)/2 + 1
2566       inc = 1
2567       t = zero
2568       do 150 i = 1, p
2569          if (i .eq. k) go to 130
2570          aki = dabs(dihdi(k1))
2571          si = w(i)
2572          j = diag0 + i
2573          t1 = half * (akk - w(j) + si - aki)
2574          t1 = t1 + dsqrt(t1*t1 + sk*aki)
2575          if (t .lt. t1) t = t1
2576          if (i .lt. k) go to 140
2577  130     inc = i
2578  140     k1 = k1 + inc
2579  150     continue
2580 c
2581       w(emin) = akk - t
2582       uk = v(dgnorm)/rad - w(emin)
2583       if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2584       if (uk .le. zero) uk = p001
2585 c
2586 c  ***  compute gerschgorin (over-)estimate of largest eigenvalue  ***
2587 c
2588       k = 1
2589       t1 = w(diag) + w(1)
2590       if (p .le. 1) go to 170
2591       do 160 i = 2, p
2592          j = diag0 + i
2593          t = w(j) + w(i)
2594          if (t .le. t1) go to 160
2595               t1 = t
2596               k = i
2597  160     continue
2598 c
2599  170  sk = w(k)
2600       j = diag0 + k
2601       akk = w(j)
2602       k1 = k*(k-1)/2 + 1
2603       inc = 1
2604       t = zero
2605       do 200 i = 1, p
2606          if (i .eq. k) go to 180
2607          aki = dabs(dihdi(k1))
2608          si = w(i)
2609          j = diag0 + i
2610          t1 = half * (w(j) + si - aki - akk)
2611          t1 = t1 + dsqrt(t1*t1 + sk*aki)
2612          if (t .lt. t1) t = t1
2613          if (i .lt. k) go to 190
2614  180     inc = i
2615  190     k1 = k1 + inc
2616  200     continue
2617 c
2618       w(emax) = akk + t
2619       lk = dmax1(lk, v(dgnorm)/rad - w(emax))
2620 c
2621 c     ***  alphak = current value of alpha (see alg. notes above).  we
2622 c     ***  use more*s scheme for initializing it.
2623       alphak = dabs(v(stppar)) * v(rad0)/rad
2624 c
2625       if (irc .ne. 0) go to 210
2626 c
2627 c  ***  compute l0 for positive definite h  ***
2628 c
2629       call livmul(p, w, l, w(q))
2630       t = v2norm(p, w)
2631       w(phipin) = dst / t / t
2632       lk = dmax1(lk, phi*w(phipin))
2633 c
2634 c  ***  safeguard alphak and add alphak*i to (d**-1)*h*(d**-1)  ***
2635 c
2636  210  ka = ka + 1
2637       if (-v(dst0) .ge. alphak .or. alphak .lt. lk .or. alphak .ge. uk)
2638      1                      alphak = uk * dmax1(p001, dsqrt(lk/uk))
2639       if (alphak .le. zero) alphak = half * uk
2640       if (alphak .le. zero) alphak = uk
2641       k = 0
2642       do 220 i = 1, p
2643          k = k + i
2644          j = diag0 + i
2645          dihdi(k) = w(j) + alphak
2646  220     continue
2647 c
2648 c  ***  try computing cholesky decomposition  ***
2649 c
2650       call lsqrt(1, p, l, dihdi, irc)
2651       if (irc .eq. 0) go to 240
2652 c
2653 c  ***  (d**-1)*h*(d**-1) + alphak*i  is indefinite -- overestimate
2654 c  ***  smallest eigenvalue for use in updating lk  ***
2655 c
2656       j = (irc*(irc+1))/2
2657       t = l(j)
2658       l(j) = one
2659       do 230 i = 1, irc
2660  230     w(i) = zero
2661       w(irc) = one
2662       call litvmu(irc, w, l, w)
2663       t1 = v2norm(irc, w)
2664       lk = alphak - t/t1/t1
2665       v(dst0) = -lk
2666       go to 210
2667 c
2668 c  ***  alphak makes (d**-1)*h*(d**-1) positive definite.
2669 c  ***  compute q = -d*step, check for convergence.  ***
2670 c
2671  240  call livmul(p, w(q), l, dig)
2672       gtsta = dotprd(p, w(q), w(q))
2673       call litvmu(p, w(q), l, w(q))
2674       dst = v2norm(p, w(q))
2675       phi = dst - rad
2676       if (phi .le. phimax .and. phi .ge. phimin) go to 270
2677       if (phi .eq. oldphi) go to 270
2678       oldphi = phi
2679       if (phi .lt. zero) go to 330
2680 c
2681 c  ***  unacceptable alphak -- update lk, uk, alphak  ***
2682 c
2683  250  if (ka .ge. kalim) go to 270
2684 c     ***  the following dmin1 is necessary because of restarts  ***
2685       if (phi .lt. zero) uk = dmin1(uk, alphak)
2686 c     *** kamin = 0 only iff the gradient vanishes  ***
2687       if (kamin .eq. 0) go to 210
2688       call livmul(p, w, l, w(q))
2689       t1 = v2norm(p, w)
2690       alphak = alphak  +  (phi/t1) * (dst/t1) * (dst/rad)
2691       lk = dmax1(lk, alphak)
2692       go to 210
2693 c
2694 c  ***  acceptable step on first try  ***
2695 c
2696  260  alphak = zero
2697 c
2698 c  ***  successful step in general.  compute step = -(d**-1)*q  ***
2699 c
2700  270  do 280 i = 1, p
2701          j = q0 + i
2702          step(i) = -w(j)/d(i)
2703  280     continue
2704       v(gtstep) = -gtsta
2705       v(preduc) = half * (dabs(alphak)*dst*dst + gtsta)
2706       go to 410
2707 c
2708 c
2709 c  ***  restart with new radius  ***
2710 c
2711  290  if (v(dst0) .le. zero .or. v(dst0) - rad .gt. phimax) go to 310
2712 c
2713 c     ***  prepare to return newton step  ***
2714 c
2715          restrt = .true.
2716          ka = ka + 1
2717          k = 0
2718          do 300 i = 1, p
2719               k = k + i
2720               j = diag0 + i
2721               dihdi(k) = w(j)
2722  300          continue
2723          uk = negone
2724          go to 30
2725 c
2726  310  kamin = ka + 3
2727       if (v(dgnorm) .eq. zero) kamin = 0
2728       if (ka .eq. 0) go to 50
2729 c
2730       dst = w(dstsav)
2731       alphak = dabs(v(stppar))
2732       phi = dst - rad
2733       t = v(dgnorm)/rad
2734       uk = t - w(emin)
2735       if (v(dgnorm) .eq. zero) uk = uk + p001 + p001*uk
2736       if (uk .le. zero) uk = p001
2737       if (rad .gt. v(rad0)) go to 320
2738 c
2739 c        ***  smaller radius  ***
2740          lk = zero
2741          if (alphak .gt. zero) lk = w(lk0)
2742          lk = dmax1(lk, t - w(emax))
2743          if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2744          go to 250
2745 c
2746 c     ***  bigger radius  ***
2747  320  if (alphak .gt. zero) uk = dmin1(uk, w(uk0))
2748       lk = dmax1(zero, -v(dst0), t - w(emax))
2749       if (v(dst0) .gt. zero) lk = dmax1(lk, (v(dst0)-rad)*w(phipin))
2750       go to 250
2751 c
2752 c  ***  decide whether to check for special case... in practice (from
2753 c  ***  the standpoint of the calling optimization code) it seems best
2754 c  ***  not to check until a few iterations have failed -- hence the
2755 c  ***  test on kamin below.
2756 c
2757  330  delta = alphak + dmin1(zero, v(dst0))
2758       twopsi = alphak*dst*dst + gtsta
2759       if (ka .ge. kamin) go to 340
2760 c     *** if the test in ref. 2 is satisfied, fall through to handle
2761 c     *** the special case (as soon as the more-sorensen test detects
2762 c     *** it).
2763       if (delta .ge. psifac*twopsi) go to 370
2764 c
2765 c  ***  check for the special case of  h + alpha*d**2  (nearly)
2766 c  ***  singular.  use one step of inverse power method with start
2767 c  ***  from lsvmin to obtain approximate eigenvector corresponding
2768 c  ***  to smallest eigenvalue of (d**-1)*h*(d**-1).  lsvmin returns
2769 c  ***  x and w with  l*w = x.
2770 c
2771  340  t = lsvmin(p, l, w(x), w)
2772 c
2773 c     ***  normalize w  ***
2774       do 350 i = 1, p
2775  350     w(i) = t*w(i)
2776 c     ***  complete current inv. power iter. -- replace w by (l**-t)*w.
2777       call litvmu(p, w, l, w)
2778       t2 = one/v2norm(p, w)
2779       do 360 i = 1, p
2780  360     w(i) = t2*w(i)
2781       t = t2 * t
2782 c
2783 c  ***  now w is the desired approximate (unit) eigenvector and
2784 c  ***  t*x = ((d**-1)*h*(d**-1) + alphak*i)*w.
2785 c
2786       sw = dotprd(p, w(q), w)
2787       t1 = (rad + dst) * (rad - dst)
2788       root = dsqrt(sw*sw + t1)
2789       if (sw .lt. zero) root = -root
2790       si = t1 / (sw + root)
2791 c
2792 c  ***  the actual test for the special case...
2793 c
2794       if ((t2*si)**2 .le. eps*(dst**2 + alphak*radsq)) go to 380
2795 c
2796 c  ***  update upper bound on smallest eigenvalue (when not positive)
2797 c  ***  (as recommended by more and sorensen) and continue...
2798 c
2799       if (v(dst0) .le. zero) v(dst0) = dmin1(v(dst0), t2**2 - alphak)
2800       lk = dmax1(lk, -v(dst0))
2801 c
2802 c  ***  check whether we can hope to detect the special case in
2803 c  ***  the available arithmetic.  accept step as it is if not.
2804 c
2805 c     ***  if not yet available, obtain machine dependent value dgxfac.
2806  370  if (dgxfac .eq. zero) dgxfac = epsfac * rmdcon(3)
2807 c
2808       if (delta .gt. dgxfac*w(dggdmx)) go to 250
2809          go to 270
2810 c
2811 c  ***  special case detected... negate alphak to indicate special case
2812 c
2813  380  alphak = -alphak
2814       v(preduc) = half * twopsi
2815 c
2816 c  ***  accept current step if adding si*w would lead to a
2817 c  ***  further relative reduction in psi of less than v(epslon)/3.
2818 c
2819       t1 = zero
2820       t = si*(alphak*sw - half*si*(alphak + t*dotprd(p,w(x),w)))
2821       if (t .lt. eps*twopsi/six) go to 390
2822          v(preduc) = v(preduc) + t
2823          dst = rad
2824          t1 = -si
2825  390  do 400 i = 1, p
2826          j = q0 + i
2827          w(j) = t1*w(i) - w(j)
2828          step(i) = w(j) / d(i)
2829  400     continue
2830       v(gtstep) = dotprd(p, dig, w(q))
2831 c
2832 c  ***  save values for use in a possible restart  ***
2833 c
2834  410  v(dstnrm) = dst
2835       v(stppar) = alphak
2836       w(lk0) = lk
2837       w(uk0) = uk
2838       v(rad0) = rad
2839       w(dstsav) = dst
2840 c
2841 c     ***  restore diagonal of dihdi  ***
2842 c
2843       j = 0
2844       do 420 i = 1, p
2845          j = j + i
2846          k = diag0 + i
2847          dihdi(j) = w(k)
2848  420     continue
2849 c
2850  999  return
2851 c
2852 c  ***  last card of gqtst follows  ***
2853       end
2854       subroutine lsqrt(n1, n, l, a, irc)
2855 c
2856 c  ***  compute rows n1 through n of the cholesky factor  l  of
2857 c  ***  a = l*(l**t),  where  l  and the lower triangle of  a  are both
2858 c  ***  stored compactly by rows (and may occupy the same storage).
2859 c  ***  irc = 0 means all went well.  irc = j means the leading
2860 c  ***  principal  j x j  submatrix of  a  is not positive definite --
2861 c  ***  and  l(j*(j+1)/2)  contains the (nonpos.) reduced j-th diagonal.
2862 c
2863 c  ***  parameters  ***
2864 c
2865       integer n1, n, irc
2866 cal   double precision l(1), a(1)
2867       double precision l(n*(n+1)/2), a(n*(n+1)/2)
2868 c     dimension l(n*(n+1)/2), a(n*(n+1)/2)
2869 c
2870 c  ***  local variables  ***
2871 c
2872       integer i, ij, ik, im1, i0, j, jk, jm1, j0, k
2873       double precision t, td, zero
2874 c
2875 c  ***  intrinsic functions  ***
2876 c/+
2877       double precision dsqrt
2878 c/
2879 c/6
2880 c     data zero/0.d+0/
2881 c/7
2882       parameter (zero=0.d+0)
2883 c/
2884 c
2885 c  ***  body  ***
2886 c
2887       i0 = n1 * (n1 - 1) / 2
2888       do 50 i = n1, n
2889          td = zero
2890          if (i .eq. 1) go to 40
2891          j0 = 0
2892          im1 = i - 1
2893          do 30 j = 1, im1
2894               t = zero
2895               if (j .eq. 1) go to 20
2896               jm1 = j - 1
2897               do 10 k = 1, jm1
2898                    ik = i0 + k
2899                    jk = j0 + k
2900                    t = t + l(ik)*l(jk)
2901  10                continue
2902  20           ij = i0 + j
2903               j0 = j0 + j
2904               t = (a(ij) - t) / l(j0)
2905               l(ij) = t
2906               td = td + t*t
2907  30           continue
2908  40      i0 = i0 + i
2909          t = a(i0) - td
2910          if (t .le. zero) go to 60
2911          l(i0) = dsqrt(t)
2912  50      continue
2913 c
2914       irc = 0
2915       go to 999
2916 c
2917  60   l(i0) = t
2918       irc = i
2919 c
2920  999  return
2921 c
2922 c  ***  last card of lsqrt  ***
2923       end
2924       double precision function lsvmin(p, l, x, y)
2925 c
2926 c  ***  estimate smallest sing. value of packed lower triang. matrix l
2927 c
2928 c  ***  parameter declarations  ***
2929 c
2930       integer p
2931 cal   double precision l(1), x(p), y(p)
2932       double precision l(p*(p+1)/2), x(p), y(p)
2933 c     dimension l(p*(p+1)/2)
2934 c
2935 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2936 c
2937 c  ***  purpose  ***
2938 c
2939 c     this function returns a good over-estimate of the smallest
2940 c     singular value of the packed lower triangular matrix l.
2941 c
2942 c  ***  parameter description  ***
2943 c
2944 c  p (in)  = the order of l.  l is a  p x p  lower triangular matrix.
2945 c  l (in)  = array holding the elements of  l  in row order, i.e.
2946 c             l(1,1), l(2,1), l(2,2), l(3,1), l(3,2), l(3,3), etc.
2947 c  x (out) if lsvmin returns a positive value, then x is a normalized
2948 c             approximate left singular vector corresponding to the
2949 c             smallest singular value.  this approximation may be very
2950 c             crude.  if lsvmin returns zero, then some components of x
2951 c             are zero and the rest retain their input values.
2952 c  y (out) if lsvmin returns a positive value, then y = (l**-1)*x is an
2953 c             unnormalized approximate right singular vector correspond-
2954 c             ing to the smallest singular value.  this approximation
2955 c             may be crude.  if lsvmin returns zero, then y retains its
2956 c             input value.  the caller may pass the same vector for x
2957 c             and y (nonstandard fortran usage), in which case y over-
2958 c             writes x (for nonzero lsvmin returns).
2959 c
2960 c  ***  algorithm notes  ***
2961 c
2962 c     the algorithm is based on (1), with the additional provision that
2963 c     lsvmin = 0 is returned if the smallest diagonal element of l
2964 c     (in magnitude) is not more than the unit roundoff times the
2965 c     largest.  the algorithm uses a random number generator proposed
2966 c     in (4), which passes the spectral test with flying colors -- see
2967 c     (2) and (3).
2968 c
2969 c  ***  subroutines and functions called  ***
2970 c
2971 c        v2norm - function, returns the 2-norm of a vector.
2972 c
2973 c  ***  references  ***
2974 c
2975 c     (1) cline, a., moler, c., stewart, g., and wilkinson, j.h.(1977),
2976 c         an estimate for the condition number of a matrix, report
2977 c         tm-310, applied math. div., argonne national laboratory.
2978 c
2979 c     (2) hoaglin, d.c. (1976), theoretical properties of congruential
2980 c         random-number generators --  an empirical view,
2981 c         memorandum ns-340, dept. of statistics, harvard univ.
2982 c
2983 c     (3) knuth, d.e. (1969), the art of computer programming, vol. 2
2984 c         (seminumerical algorithms), addison-wesley, reading, mass.
2985 c
2986 c     (4) smith, c.s. (1971), multiplicative pseudo-random number
2987 c         generators with prime modulus, j. assoc. comput. mach. 18,
2988 c         pp. 586-593.
2989 c
2990 c  ***  history  ***
2991 c
2992 c     designed and coded by david m. gay (winter 1977/summer 1978).
2993 c
2994 c  ***  general  ***
2995 c
2996 c     this subroutine was written in connection with research
2997 c     supported by the national science foundation under grants
2998 c     mcs-7600324, dcr75-10143, 76-14311dss, and mcs76-11989.
2999 c
3000 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3001 c
3002 c  ***  local variables  ***
3003 c
3004       integer i, ii, ix, j, ji, jj, jjj, jm1, j0, pm1
3005       double precision b, sminus, splus, t, xminus, xplus
3006 c
3007 c  ***  constants  ***
3008 c
3009       double precision half, one, r9973, zero
3010 c
3011 c  ***  intrinsic functions  ***
3012 c/+
3013       integer mod
3014       real float
3015       double precision dabs
3016 c/
3017 c  ***  external functions and subroutines  ***
3018 c
3019       external dotprd, v2norm, vaxpy
3020       double precision dotprd, v2norm
3021 c
3022 c/6
3023 c     data half/0.5d+0/, one/1.d+0/, r9973/9973.d+0/, zero/0.d+0/
3024 c/7
3025       parameter (half=0.5d+0, one=1.d+0, r9973=9973.d+0, zero=0.d+0)
3026 c/
3027 c
3028 c  ***  body  ***
3029 c
3030       ix = 2
3031       pm1 = p - 1
3032 c
3033 c  ***  first check whether to return lsvmin = 0 and initialize x  ***
3034 c
3035       ii = 0
3036       j0 = p*pm1/2
3037       jj = j0 + p
3038       if (l(jj) .eq. zero) go to 110
3039       ix = mod(3432*ix, 9973)
3040       b = half*(one + float(ix)/r9973)
3041       xplus = b / l(jj)
3042       x(p) = xplus
3043       if (p .le. 1) go to 60
3044       do 10 i = 1, pm1
3045          ii = ii + i
3046          if (l(ii) .eq. zero) go to 110
3047          ji = j0 + i
3048          x(i) = xplus * l(ji)
3049  10      continue
3050 c
3051 c  ***  solve (l**t)*x = b, where the components of b have randomly
3052 c  ***  chosen magnitudes in (.5,1) with signs chosen to make x large.
3053 c
3054 c     do j = p-1 to 1 by -1...
3055       do 50 jjj = 1, pm1
3056          j = p - jjj
3057 c       ***  determine x(j) in this iteration. note for i = 1,2,...,j
3058 c       ***  that x(i) holds the current partial sum for row i.
3059          ix = mod(3432*ix, 9973)
3060          b = half*(one + float(ix)/r9973)
3061          xplus = (b - x(j))
3062          xminus = (-b - x(j))
3063          splus = dabs(xplus)
3064          sminus = dabs(xminus)
3065          jm1 = j - 1
3066          j0 = j*jm1/2
3067          jj = j0 + j
3068          xplus = xplus/l(jj)
3069          xminus = xminus/l(jj)
3070          if (jm1 .eq. 0) go to 30
3071          do 20 i = 1, jm1
3072               ji = j0 + i
3073               splus = splus + dabs(x(i) + l(ji)*xplus)
3074               sminus = sminus + dabs(x(i) + l(ji)*xminus)
3075  20           continue
3076  30      if (sminus .gt. splus) xplus = xminus
3077          x(j) = xplus
3078 c       ***  update partial sums  ***
3079          if (jm1 .gt. 0) call vaxpy(jm1, x, xplus, l(j0+1), x)
3080  50      continue
3081 c
3082 c  ***  normalize x  ***
3083 c
3084  60   t = one/v2norm(p, x)
3085       do 70 i = 1, p
3086  70      x(i) = t*x(i)
3087 c
3088 c  ***  solve l*y = x and return lsvmin = 1/twonorm(y)  ***
3089 c
3090       do 100 j = 1, p
3091          jm1 = j - 1
3092          j0 = j*jm1/2
3093          jj = j0 + j
3094          t = zero
3095          if (jm1 .gt. 0) t = dotprd(jm1, l(j0+1), y)
3096          y(j) = (x(j) - t) / l(jj)
3097  100     continue
3098 c
3099       lsvmin = one/v2norm(p, y)
3100       go to 999
3101 c
3102  110  lsvmin = zero
3103  999  return
3104 c  ***  last card of lsvmin follows  ***
3105       end
3106       subroutine slvmul(p, y, s, x)
3107 c
3108 c  ***  set  y = s * x,  s = p x p symmetric matrix.  ***
3109 c  ***  lower triangle of  s  stored rowwise.         ***
3110 c
3111 c  ***  parameter declarations  ***
3112 c
3113       integer p
3114 cal   double precision s(1), x(p), y(p)
3115       double precision s(p*(p+1)/2), x(p), y(p)
3116 c     dimension s(p*(p+1)/2)
3117 c
3118 c  ***  local variables  ***
3119 c
3120       integer i, im1, j, k
3121       double precision xi
3122 c
3123 c  ***  no intrinsic functions  ***
3124 c
3125 c  ***  external function  ***
3126 c
3127       external dotprd
3128       double precision dotprd
3129 c
3130 c-----------------------------------------------------------------------
3131 c
3132       j = 1
3133       do 10 i = 1, p
3134          y(i) = dotprd(i, s(j), x)
3135          j = j + i
3136  10      continue
3137 c
3138       if (p .le. 1) go to 999
3139       j = 1
3140       do 40 i = 2, p
3141          xi = x(i)
3142          im1 = i - 1
3143          j = j + 1
3144          do 30 k = 1, im1
3145               y(k) = y(k) + s(j)*xi
3146               j = j + 1
3147  30           continue
3148  40      continue
3149 c
3150  999  return
3151 c  ***  last card of slvmul follows  ***
3152       end