*     ************************************************************************

      subroutine scbfgs( n, h, s, y, epsmch, ynoise, negcur, work )

*     ************************************************************************

*     Purpose :
*     ---------

*     This routine applies the bfgs quasi-Newton update to the matrix H
*     using the step S and the difference in gradient Y.  The matrix H is
*     supposed to be stored columnwise and only its lower triangular part
*     is stored.
*     If yTs is not positive, NEGCUR is set to .true., indicating that
*     the curvature along the step is negative and another update should
*     be applied.

*     Parameters :
*     ------------

*     n       ( int )

*       input  : the dimension of the matrix h.
*       output : unmodified.

*     h       ( dble )

*       input  : a vector that contains the successive columns of the
*                lower triangular part of the n*n matrix to be updated.
*       output : the matrix on input has been updated using the bfgs
*                updating formula, if the curvature was found to be
*                positive enough and the vector Y-H*S sufficiently
*                different from zero. In all other cases, the update
*                is skipped.

*     s       ( dble )

*       input  : a vector of dimension n containing the step.
*       output : unmodified.

*     y       ( dble )

*       input  : a vector of dimension n containing the change in gradients
*                corresponding to the step S.
*       output : meaningless.

*     epsmch  ( dble )

*       input  : the relative precision of the machine.
*       output : unmodified.

*     ynoise  ( dble )

*       input  : an estimation of the noise present in the vector Y.
*       output : unmodified.

*     negcur  ( log )

*       input  : arbitrary.
*       output : .true. iff yTs is clearly negative.  In this case,
*                the curvature of the objective function along the step
*                S is negative and bfgs cannot and has not been applied.

*     work    ( dble )

*       input  : a vector of length n whose content is arbitrary.
*       output : meaningless.


*     Routines used :
*     ---------------

*     dnrm2,  ddot, dscal, dcopy (blasd)
*     dsetvl, dtscal, ltcmvp

*     Programming :
*     -------------

*     Ph.L. Toint 

*     ========================================================================

*     Routine parameters

      integer          n
      logical          negcur
      double precision h(*), s(*), y(*), work(*), epsmch, ynoise

*     Internal variables

      integer          i, j, l
      double precision yts, shs

      double precision dnrm2, ddot

      double precision  zero, one, two, three, half, tenm1, tenm2, tenm4
      common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4

      yts    = ddot( n, y, 1, s, 1 )
      negcur = yts .le. sqrt(epsmch)*dnrm2(n,s,1)*dnrm2(n,y,1)

      if( negcur ) return

      call dsetvl( n, work, 1, zero )
      call ltcmvp( n, h, s, work, .true. )

      do 100 i = 1 , n
        if( abs(work(i)).gt.ynoise ) go to 200
  100 continue
      return

  200 continue
      shs    = ddot( n, work, 1, s, 1 )
      negcur = shs.le.zero
      if( negcur ) return
      shs = one/sqrt(ddot(n,work,1,s,1))
      call dscal( n, shs, work, 1 )
      l = 1
      do 400 i = 1 , n
        do 500 j = i , n
          h(l) = h(l) - work(i)*work(j)
          l    = l + 1
  500   continue
  400 continue

      call dcopy( n,  y, 1, work, 1 )
      call dtscal( n, (one/yts), work, 1, y, 1 )
      l = 1
      do 700 i = 1 , n
        do 800 j = i , n
          h(l) = h(l) + y(i)*work(j)
          l    = l + 1
  800   continue
  700 continue

      return
      end