#! /bin/csh
# conjgrad.sh  Conjugate-Gradient Iteration on Radon Inversion
# 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 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 up Zero "Previous Step" for 1st Iteration
cp temp.ds step.ds
# Set up Zero "Previous Conjugate Step" for 1st Iteration
fac 22500 0 <test.dt >step.dt

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

	# Find Residual Size R.R
	set RR=`dot 22500 resid.dt resid.dt `
	echo $i $RR >>conjgrad.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 Conjugate Gradient G = (A)g
	radon par=slpar i=grad.ds o=grad.dt >&/dev/null

	# Compute Dot Products
	set GG=`dot 22500 grad.dt grad.dt `
	set RG=`dot 22500 resid.dt grad.dt `
	if ( $i == 0 ) then
		# First Iteration uses Steepest Descent
		set alpha=`divide $RG $GG `
		set beta=0
	else
		# Next Iterations use Conjugate Gradients
		# Compute Dot Products
		set RS=`dot 22500 resid.dt step.dt `
		set GS=`dot 22500 grad.dt step.dt `
		set SS=`dot 22500 step.dt step.dt `
		# Compute alpha = 1/det((S.S)(G.R) - (S.G)(S.R))
		set alpha=`getalpha $GG $SS $GS $RG $RS `
		# Compute beta = 1/det( -(S.G)(G.R) + (G.G)(S.R))
		set beta=`getbeta $GG $SS $GS $RG $RS `
	endif
	echo Iteration $i : residual $RR ,  alpha $alpha , beta $beta

	# Compute alpha.g
	fac 22500 $alpha <grad.ds >agrad.ds
	rm grad.ds
	# Compute beta.s
	fac 22500 $beta <step.ds >bstep.ds
	# Compute Next Step s = alpha.g + beta.s
	add 22500 agrad.ds bstep.ds >step.ds
	rm agrad.ds bstep.ds
	# Step s Saved for Next Iteration
	# Estimate Next Model x <- x + s
	add 22500 step.ds temp.ds >new.ds
	mv new.ds temp.ds

	# Compute alpha.G
	fac 22500 $alpha <grad.dt >agrad.dt
	rm grad.dt
	# Compute beta.S
	fac 22500 $beta <step.dt >bstep.dt
	# Compute Next Conjugate Step S = alpha.G + beta.S
	add 22500 agrad.dt bstep.dt >step.dt
	rm agrad.dt bstep.dt
	# Conjugate Step S Saved for Next Iteration
	# Estimate Next Residual R <- R - S
	neg 22500 <step.dt >mstep.dt
	add 22500 resid.dt mstep.dt >new.dt
	rm mstep.dt
	mv new.dt resid.dt

	#cat temp.ds >>all.ds
	@ i = ( $i + 1)
end
mv temp.ds conjgrad.ds