* ************************************************************************ subroutine scrk1( n, h, s, y, epsmch, ynoise, upd, work ) * ************************************************************************ * Purpose : * --------- * This routine applies the rank one 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 (y-hs)Ts is too small, the update is skipped and UPD is returned * .false.. * 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 rk1 * updating formula, if its denominator has been found * to be large 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. * upd ( log ) * input : arbitrary. * output : .true. iff the update has been performed on H, and * .false. otherwise. * work ( dble ) * input : a vector of length n whose content is arbitrary. * output : meaningless. * Routines used : * --------------- * dcopy, dnrm2, ddot (blasd) * dtscal, ltcmvp * Programming : * ------------- * Ph.L. Toint * ======================================================================== * Routine parameters integer n logical upd double precision h(*), s(*), y(*), work(*), epsmch, ynoise * Internal variables integer i, j, l double precision rts double precision dnrm2, ddot double precision zero, one, two, three, half, tenm1, tenm2, tenm4 common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4 call dcopy( n, y, 1, work, 1 ) call ltcmvp( n, h, s, work, .false. ) rts = ddot( n, s, 1, work, 1 ) upd = abs(rts) .gt. sqrt(epsmch)*dnrm2(n,s,1)*dnrm2(n,work,1) if( upd ) then do 100 i = 1 , n if( abs(work(i)).gt.ynoise ) go to 200 100 continue upd = .false. endif 200 continue if( upd ) then call dtscal( n, (one/rts), 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 endif return end