
//	module		: implementation of linear equation 
//	file name	: cg_1d_list.cpp
//	version		: 1.0
//	author		: Lung Sheng Chien
//	date		: 2004 / 1 / 26
//  
//	description : 
//			
//
#include "def.h"
#include <stdio.h>
#include <math.h>
#include "cg_1d_array.h"
#include "cg_1d_list.h"

int cg_1d_list(int nmax , int nx 
				 , linear_eqn *A , double* x
				 , double* b     , double* r   
				 , double* v     , double* z
				 , int imax , double eps , double delta )
{
	int i , k ;
	double c , d , t  ;

// 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( nmax , nx , r , r ) ;

	for( k = 1 ; k < imax ; k++){
		if  ( delta > sqrt(inner_product(nmax , nx , 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( nmax , nx, A, v, z ) ;
		t = c / inner_product( nmax , nx , 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( nmax , nx , 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(int nmax , int nx , double* x , double* y )
{
	int i ;
	double sum = 0; 
	
	for(i=1 ; i <= nx-1 ; i++){
		sum = sum + x[i]*y[i] ;
	}
	return sum ; 
}

void matrix_vector(int nmax , int nx , linear_eqn* A , double* x , 
				   double* y )
{
	int i ;
	
	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] ; 
	}

}