
//	module		: implementation of linear equation 
//	file name	: cg_1d_array.cpp
//	version		: 1.0
//	author		: Lung Sheng Chien
//	date		: 2004 / 1 / 24
//  
//	description : implement with array ( call by reference )
//				re-implement with dynamic allocation later 	
//

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "cg_1d_array.h"

int cg_1d_array(linear_eqn A[COL_DIM] , double x[COL_DIM] 
				 , double b[COL_DIM]  , double r[COL_DIM]  
				 , double v[COL_DIM]  , double z[COL_DIM]
				 , int imax , double eps , double delta )
{
	int i , k ;
	double c , d , t  ;
	int nx ;

	nx = COL_DIM - 1 ;
// step 1 : r := b - Ax		
	for(i=1 ; i <= nx-1 ; i++){
		r[i] = A[i].coeff[0]*x[i  ] 
			 + A[i].coeff[1]*x[i-1]
			 + A[i].coeff[2]*x[i+1] ;
		r[i] = b[i] - r[i] ;
	}
// step 2 : v := r 
	for(i=1 ; i <= nx-1 ; i++){
		v[i] = r[i] ;
	}
// step 3 : c := <r,r>
	c = inner_product( r , r ) ;

	for( k = 1 ; k < imax ; k++){
		if  ( delta > sqrt(inner_product(v,v)) ){
			printf("k=%d , <r,r> = %g\n",k,d);
			return k ;
		}
// step 4 : update  x(k+1) = x(k) + t(k)v(k)
//					r(k+1) = r(k) - t(k)Av(k)
		matrix_vector( A, v, z ) ;
		t = c / inner_product( v , z ) ;
		for(i=1 ; i <= nx-1 ; i++){
			x[i] = x[i] + t*v[i] ;
			r[i] = r[i] - t*z[i] ;
		}

		d = inner_product( r , r ) ;
		if ( eps > d  ) {
			printf("k=%d , <r,r> = %g\n",k,d);
			return k ;
		}
		for(i=1 ; i <= nx-1 ; i++){
			v[i] = r[i] + (d/c)*v[i] ;
		}
		c = d ;
		printf("k=%d , <r,r> = %g\n",k,d);
	}// end for(k)
	return k ;
}

double inner_product(double x[COL_DIM], double y[COL_DIM] )
{
	int i ;
	double sum = 0; 
	int nx ;

	nx = COL_DIM - 1 ;
	for(i=1 ; i <= nx-1 ; i++){
		sum = sum + x[i]*y[i] ;
	}
	return sum ; 
}


void matrix_vector(linear_eqn A[COL_DIM], double x[COL_DIM], 
				   double y[COL_DIM] )
{
	int i ;
	int nx ;

	nx = COL_DIM - 1 ;
	for(i=1 ; i <= nx-1 ; i++){
		y[i] = A[i].coeff[0]*x[i  ] 
			 + A[i].coeff[1]*x[i-1]
			 + A[i].coeff[2]*x[i+1] ; 
	}

}
