#! /bin/csh
# hestenes.sh  Conjugate-Gradient Iteration on Radon Inversion, after
#		Hestenes and Stiefel (1952), Claerbout (1992).
# All programs here in /quake/data/rg
# Set Number of Iterations
set niter=200
# Data test.dt already computed from True Model test.ds

# Set Initial Model x to Zero
fac 22500 0 <test.ds >temp.ds
#cat temp.ds >all.ds
# Set Initial Residual to Data test.dt , R = Y
cp test.dt resid.dt

# Set Gradient g = (A transpose)Y
radon par=islpar i=test.dt o=grad.ds >&/dev/null
# Set Previous Step to Gradient
cp grad.ds step.ds
# Find gamma1 = g.g
set gamma1=`dot 22500 grad.ds grad.ds `

# Iterate, Two Transforms per Iteration
set i=0
while ($i < $niter )

	# Find Conjugate Step S = (A)s
	radon par=slpar i=step.ds o=step.dt >&/dev/null

	# Compute S.S
	set SS=`dot 22500 step.dt step.dt `
	# Compute alpha = gamma1/SS
	set alpha=`divide $gamma1 $SS `
	# Compute alpha.s
	fac 22500 $alpha <step.ds >astep.ds
	# Compute Next Model x <- x + alpha.s
	add 22500 astep.ds temp.ds >new.ds
	rm astep.ds
	mv new.ds temp.ds

	# Compute -alpha.S
	fac 22500 $alpha <step.dt >astep.dt
	neg 22500 <astep.dt >mastep.dt
	rm astep.dt
	# Find Residual R <- R - alpha.S
	add 22500 resid.dt mastep.dt >new.dt
	rm mastep.dt
	mv new.dt resid.dt

	# Find Residual Size R.R
	set RR=`dot 22500 resid.dt resid.dt `
	echo $i $RR >>hestenes.res

	# Find Gradient g = (A transpose)R
	radon par=islpar i=resid.dt o=grad.ds >&/dev/null

	# Condition Gradient by Smoothing
	smooth i=grad.ds o=grad.ds nx=150 nt=150 xw=3 tw=3

	# Find gamma = g.g
	set gamma=`dot 22500 grad.ds grad.ds `
	# Compute beta = gamma/gamma1
	set beta=`divide $gamma $gamma1 `
	# Set gamma1 = gamma
	set gamma1=$gamma

	# Compute beta.s
	fac 22500 $beta <step.ds >bstep.ds
	# Compute Next Step s <- g + beta.s
	add 22500 grad.ds bstep.ds >new.ds
	rm bstep.ds
	mv new.ds step.ds

	echo Iteration $i : Residual $RR  alpha $alpha  beta $beta
	#cat temp.ds >>all.ds
	@ i = ( $i + 1)
end
mv temp.ds hestenes.ds