program main

c*********************************************************************72
c
cc MAIN is the main program for TOMS611_PRB.
c
c  Discussion:
c
c    TOMS611_PRB tests the TOMS611 library.
c
c     ***  simple test program for smsno, sumsl, and humsl  ***
c
c  in these examples,  n = 4  and  f(x) = (1.0 + 0.5*(x1**t)*a*x1)**0.5,
c  where  x1(i) = d1(i)*x(i) - i,  **t denotes transpose, and  a  is a
c  matrix having fives on the main diagonal and ones everywhere else.
c  the scale vector  d1  is passed to qdrtf, the subroutine that
c  evaluates  f,  as part of urparm.  specifically, the matrix  urp
c  declared below is passed for ufparm, and  d1  is urp(*,1), the first
c  column of urp.  this main program repeatedly minimizes f, starting
c  from  x = 0,  by calling smsno, sumsl, and humsl.  we actually
c  use two different objective functions, since we change  d1  after
c  the first call on sumsl.  all runs but the last use  d = d1.
c
c     f(x) is minimized at  x1 = 0  (a vector of zeros), i.e., at
c  x(i) = i/d1(i).
c
      external deflt, humsl, imdcon, qdrtf, qdrtg, qdrtgh, smsno,
     1         sumsl
      integer imdcon
c
c imdcon - supplies nout, the output unit number.
c qdrtf  - passed for calcf to smsno, sumsl, and humsl.
c qdrtg  - passed for calcg to sumsl.
c qdrtgh - passed for calcgh to humsl.
c
      integer i, iv(60), liv, lv, nout, uip(1)
      double precision d(4), urp(4,3), v(150), x(4)
c
c  the length of  v  is dictated by humsl.
c
      data liv/60/, lv/150/
c
c  initialize nout, d, d1, and x.
c
      nout = imdcon(1)
c
      do 10 i = 1, 4
         d(i) = 1.d+0
         urp(i,1) = d(i)
         x(i) = 0.d+0
 10      continue
      write(nout,20)
 20   format(/16h smsno on qdrtf )
c
c  before this first call, we set iv(1) to 0 so that all input
c  components of iv and v will be given default values.  before
c  subsequent calls, we set iv(1) to 12 so that the old input values
c  of iv and v are used.
c
c  qdrtf does not make use of ufparm.  in calling smsno, we
c  arbitrarily pass qdrtf for ufparm to satifsy the calling sequence.
c
      iv(1) = 0
      call smsno(4, d, x, qdrtf, iv, liv, lv, v, uip, urp, qdrtf)
c
c  we reinitialize x and minimize  f  again, this time using sumsl.
c  qdrtg, the subroutine passed for calcg, assumes that ufparm is qdrtf.
c
      do 30 i = 1, 4
 30      x(i) = 0.d+0
      write(nout,40)
 40   format(/16h sumsl on qdrtf )
c
      iv(1) = 12
      call sumsl(4, d, x, qdrtf, qdrtg, iv, liv, lv, v, uip,urp,qdrtf)
c
c  now we modify  f  by using a different choice of d1.  we still use
c  d = d1, so the performance of sumsl should stay the same -- only d
c  and the final x and gradient should be affected.
c
      do 50 i = 1, 4
         x(i) = 0.d+0
         d(i) = 1.d2 ** i
         urp(i,1) = d(i)
 50      continue
      write(nout,40)
c
      iv(1) = 12
      call sumsl(4, d, x, qdrtf, qdrtg, iv, liv, lv, v, uip,urp,qdrtf)
c
c  using the last choice of d and d1, we now use humsl to minimize f.
c  like qdrtg, qdrtgh assumes that ufparm is qdrtf.
c
      do 60 i = 1, 4
 60      x(i) = 0.d+0
      write(nout,70)
 70   format(/16h humsl on qdrtf )
c
      iv(1) = 12
      call humsl(4, d, x, qdrtf, qdrtgh, iv, liv, lv, v, uip,urp,qdrtf)
c
c  we repeat the last run with iv(dtype) = 1 and v(dinit) = 0.0, so
c  that humsl will determine  d  from the diagonal of the hessian.
c  this run also demonstrates the use of subroutine deflt and the
c  passing of nondefault parameters.  (since the iv and v input
c  components still have their default values at his point, it is not
c  really necessary to call deflt.  it is necessary to reset iv(1) to
c  12, however, and deflt does this for us.)
c
      write(nout,80)
 80   format(/18h humsl updating d )
c
      do 90 i = 1, 4
 90      x(i) = 0.d+0
c
      call deflt(2, iv, liv, lv, v)
      iv(16) = 1
      v(38) = 0.d+0
      call humsl(4, d, x, qdrtf, qdrtgh, iv, liv, lv, v, uip,urp,qdrtf)
c
 999  stop
      end
c
c***********************************************************************
c
c     q d r t f
c
c***********************************************************************
c
      subroutine qdrtf(n, x, nf, f, uip, urp, ufp)
c
c  this routine evaluates the objective function f(x) described in the
c  main program above.  it stores in urp(*,2) and urp(*,3) some
c  information useful in evaluating the gradient and hessian of f, and
c  it stores  nf  in uip(1) to identify the x corresponding to this
c  information.  f(x) has the form  f(x) = phi(q(x)),  where  q  is a
c  quadratic form and  phi(y) = y**0.5.  the gradient of  f  is
c  g(x) = phiprm(q(x))*gq(x),  where  phiprm  is the derivative of phi
c  and  gq  is the gradient of q.  this routine stores phiprm(q(x)) in
c  urp(1,3) and gq(x) in urp(*,2).  the hessian of f is
c  h(x) = phi2prm(q(x))*gq(x)*gq(x)**t + phiprm(q(x))*hq(x),  where
c  phi2prm  is the second derivative of phi, **t denotes transpose,
c  and  hq  is the hessian of q.  this routine stores phi2prm(q(x)) in
c  urp(2,3).  the subroutines qdrtg and qdrtgh given below would work
c  without change on any other choice of phi.  qdrtg would also work
c  with any other differentiable function q.  qdrtgh, on the other
c  hand, assumes that  hq(x)  is the matrix  a  described in the main
c  program above.
c
      integer n, nf, uip(1)
      double precision x(n), f, urp(n,3)
      external ufp
c/+
      real float
      double precision dsqrt
c/
      integer i
      double precision dn, f2, t, t1
c
c
      uip(1) = nf
      dn = n
      t = 0.d+0
      do 10 i = 1, n
         urp(i,2) = urp(i,1)*x(i) - float(i)
         t = t + urp(i,2)
 10      continue
      f2 = 0.d+0
      do 20 i = 1, n
         t1 = dn*urp(i,2) + t
         f2 = f2 + t1*urp(i,2)
         urp(i,2) = urp(i,1) * t1
 20      continue
      f2 = 1.d+0  +   0.5d+0 * f2
      f = dsqrt(f2)
      urp(1,3) = 0.5d+0 / f
      urp(2,3) = -0.5d+0 / (f * f2)
 999  return
      end
c
c***********************************************************************
c
c     q d r t g
c
c***********************************************************************
c
      subroutine qdrtg(n, x, nf, g, uip, urp, qdrtf)
c
c  this routine evaluates the gradient of the objective function f(x)
c  described in the main program above.  see the comments there and in
c  subroutine qdrtf above.
c
      integer n, nf, uip(1)
      double precision x(n), g(n), urp(n,3)
      external qdrtf
c
      integer i
      double precision f
c
      if (nf .ne. uip(1)) call qdrtf(n, x, nf, f, uip, urp, qdrtf)
      do 10 i = 1, n
 10      g(i) = urp(1,3) * urp(i,2)
 999  return
      end
c
c***********************************************************************
c
c     q d r t g h
c
c***********************************************************************
c
      subroutine qdrtgh(n, x, nf, g, h, uip, urp, qdrtf)
c
c  this routine evaluates the gradient and hessian of the objective
c  function f(x) described in the main program above.  see the comments
c  there and in subroutine qdrtf above.  note that the  h  returned is
c  the lower triangle of the hessian, stored row-wise.
c
      integer n, nf, uip(1)
      double precision x(n), g(n), h(1), urp(n,3)
c     dimension h(n*(n+1)/2)
      external qdrtf
c
      integer i, j, k
      double precision dn, f, t1, t2
c
      if (nf .ne. uip(1)) call qdrtf(n, x, nf, f, uip, urp, qdrtf)
      k = 0
      dn = n
      do 20 i = 1, n
         g(i) = urp(1,3) * urp(i,2)
         t1 = urp(1,3) * urp(i,1)
         t2 = urp(i,2) * urp(2,3)
         do 10 j = 1, i
              k = k + 1
              h(k) = t2*urp(j,2) + t1*urp(j,1)
 10           continue
         h(k) = h(k) + dn*urp(i,1)*t1
 20      continue
 999  return
      end