subroutine dlmmv(n,m,s,y,rho,scale,mark,v,wa)
integer n, m, mark
double precision scale
double precision s(n,m), y(n,m), rho(m), v(n), wa(m)
c **********
c
c This subroutine computes the matrix-vector product H*v
c where H is the inverse BFGS approximation.
c
c The matrix H depends on an initial matrix H0, m steps
c s(1),...,s(m), and m gradient differences y(1),...,y(m).
c These vectors are stored in the arrays s and y.
c The most recent step and gradient difference are stored
c in columns s(1,mark) and y(1,mark), respectively.
c The initial matrix H0 is assumed to be scale*I.
c
c The subroutine statement is
c
c subroutine dlmmv(n,m,s,y,rho,scale,mark,v,wa)
c
c where
c
c n is an integer variable.
c On entry n is the number of variables.
c On exit n is unchanged.
c
c m is an integer variable.
c On entry m specifies the number of steps and gradient
c differences that are stored.
c On exit m is unchanged.
c
c s is a double precision array of dimension (n,m)
c On entry s contains the m steps.
c On exit s is unchanged.
c
c y is a double precision array of dimension (n,m)
c On entry y contains the m gradient differences.
c On exit y is unchanged.
c
c rho is a double precision array of dimension m
c On entry rho contains the m innerproducts (s(i),y(i)).
c On exit rho is unchanged.
c
c scale is a double precision variable
c On entry scale specifies the initial matrix H0 = scale*I.
c On exit scale is unchanged.
c
c mark is an integer variable.
c On entry mark points to the current s(i) and y(i).
c On exit mark is unchanged.
c
c v is a double precision array of dimension n.
c On entry v contains the vector v.
c On exit v contains the matrix-vector product H*v.
c
c wa is a double precision work array of dimension m.
c
c Subprograms called
c
c Level 1 BLAS ... daxpy, ddot, dscal
c
c MINPACK-2 project. November 1993.
c Argonne National Laboratory and University of Minnesota.
c Brett M. Averick, Richard G. Carter, and Jorge J. More'.
c
c **********
integer i, k
double precision beta
double precision ddot
external daxpy, ddot, dscal
k = mark + 1
do 10 i = 1, m
k = k - 1
if (k .eq. 0) k = m
wa(k) = ddot(n,s(1,k),1,v,1)/rho(k)
call daxpy(n,-wa(k),y(1,k),1,v,1)
10 continue
call dscal(n,scale,v,1)
do 20 i = 1, m
beta = wa(k) - ddot(n,y(1,k),1,v,1)/rho(k)
call daxpy(n,beta,s(1,k),1,v,1)
k = k + 1
if (k .eq. m+1) k = 1
20 continue
end