# A step of conjugate-gradient descent.
#
subroutine cgstep( iter,   n, x, g, s,   m, rr, gg, ss)
integer i,         iter,   n,            m
real	x(n), rr(m),	# solution, residual
	g(n), gg(m),	# gradient, conjugate gradient
	s(n), ss(m)	# step,     conjugate step
real dot, sds, gdg, gds, determ, gdr, sdr, alfa, beta
if( iter == 0 ) {
	do i= 1, n
		s(i) = 0.
	do i= 1, m
		ss(i) = 0.
	if( dot(m,gg,gg)==0 )	call erexit('cgstep: grad vanishes identically')
	alfa =  dot(m,gg,rr) / dot(m,gg,gg)
	beta = 0.
	}
else {				# search plane by solving 2-by-2
	gdg = dot(m,gg,gg)	#  G . (R - G*alfa - S*beta) = 0
	sds = dot(m,ss,ss)	#  S . (R - G*alfa - S*beta) = 0
	gds = dot(m,gg,ss)
	determ = gdg * sds - gds * gds + 1.e-15
	gdr = dot(m,gg,rr)
	sdr = dot(m,ss,rr)
	alfa = ( sds * gdr - gds * sdr ) / determ
	beta = (-gds * gdr + gdg * sdr ) / determ
	}
do i= 1, n					# s = model step
	s(i)  = alfa * g(i)  + beta * s(i)
do i= 1, m					# ss = conjugate 
	ss(i) = alfa * gg(i) + beta * ss(i)
do i= 1, n					# update solution
	x(i)  =  x(i) +  s(i)
do i= 1, m					# update residual
	rr(i) = rr(i) - ss(i)
return; end

real function dot( n, x, y )
integer i, n;		real val, x(n), y(n)
val = 0.;	do i=1,n  { val = val + x(i) * y(i) }
dot = val;	return;		end