
//	module		: conjugate gradient driver for 1-D poisson equation 
//	file name	: pcg_1d_list.cpp
//	version		: 1.0
//	author		: Lung Sheng Chien
//	date		: 2004 / 1 / 28
//  
//	description : this is for PCG test
//			Laplace U = f  on [0,1]
//			Dirichlet B.C
//
//	how to generate vfftpack.c
//  (1) D:\fft>f2c -A  -r8 -P vfftpack.f
//	(2) use ulta-edit to merge F2c.h and vfftpack.P and store as f2c.h
//	(3) insert 
//		#include <math.h> 
//		to vfftpck.c and store as vfftpack.cpp
//	(4) add f2c.h, vfftpack.cpp, LIBF77.lib, LIBI77.lib to project
//  remark 
//	(1) in Fortran 77, array index start at 1 
//		in c		 , array index start at 0
//
//		in sin1d, subroutine call
//				vsint_( &m, &n, work2 , work1 , &m , wx);
//		it needs to setup RHS work2, remember index starts at 0 
//			for(i=1 ; i <= nx-1 ; i++){
//				work2[i-1]=u[i] ;
//			}
//	(2) all parameter pass by pointer 
//
//  result :
//	1. CG converges slowly, error would be reduced in the last few iteration
//	2. 2nd order convergence
//
//	future 
//	1. check 2-D version
//	2. connect to matlab to show result
//	2. implement FFT-based PCG
//

#include "def.h"
#include <stdio.h>
#include <math.h>
#include "cg_1d_array.h"
#include "cg_1d_list.h"
#include "f2c.h"
#include "pcg_1d_list.h"

extern double *wx ;
extern double *sx ;
extern double *work1 ;
extern double *work2 ;

void gauss_seidel(int nmax , int nx 
				 , linear_eqn *A , double* x , double* b  )
{
	int i  ;
	
	for(i=1; i <= nx - 1 ; i++ ){
//		x[i] = (- A[i].coeff[1]*x[i-1]  + b[i] ) / A[i].coeff[0] ;
		x[i] = b[i] / A[i].coeff[0] ;
//		x[i] = b[i] ;
	}
	

}

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

// 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 : solve Qr_hat = r ;
#ifdef	_GAUSS_SEIDEL
	gauss_seidel( nmax , nx , A , r_hat , r   ) ;
#endif

#ifdef	_FFT
	for(i=1 ; i <= nx-1 ; i++){
		r_hat[i] = r[i] ;
	}
	poisson1d(r_hat,COL_DIM,nx,sx,wx,work1,work2);
#endif

// step 3 : v := r_hat 
	for(i=1 ; i <= nx-1 ; i++){
		v[i] = r_hat[i] ;
	}

// step 4 : c := <r_hat,r>
	c = inner_product( nmax , nx , r_hat , r ) ;

	for( k = 1 ; k < imax ; k++){
		if  ( delta > sqrt(inner_product(nmax , nx , v,v)) ){
			printf("k=%d , <r_hat,r> = %g\n",k,d);
			return k ;
		}
// step 5 : update  x(k+1) = x(k) + t_hat(k)v(k)
//					r(k+1) = r(k) - t_hat(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] ;
		}

#ifdef	_GAUSS_SEIDEL
		gauss_seidel( nmax , nx , A , r_hat , r  ) ;
#endif

#ifdef	_FFT
		for(i=1 ; i <= nx-1 ; i++){
			r_hat[i] = r[i] ;
		}
		poisson1d(r_hat,COL_DIM,nx,sx,wx,work1,work2);
#endif

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


void  setupfft(int n1,int nx,double dx,double *sx,double *wx)
{
	double pi ;
	double si ;
	int i ;
	long int	n ;

	pi=3.14159265358979 ; 

	n = nx-1 ;
    vsinti_(&n, wx) ;

    for(i=0 ; i <= nx ; i++ ){
      si=pi*((double)i)/((double) nx)/2. ;
      sx[i]=-4.*sin(si)*sin(si)/dx/dx ;
    }
}

void sin1d(double *u,int n1,int nx,double *wx,double *work1,double *work2 )
{
	int i ;
	integer  m ;
	integer  n ;

// transform in the y-direction
	for(i=1 ; i <= nx-1 ; i++){
		work2[i-1]=u[i] ;
	}
	m = 1 ;
	n = nx-1 ;
	vsint_( &m, &n, work2 , work1 , &m , wx);
	for(i=1 ; i <= nx-1 ; i++){
		u[i]=work2[i-1] ;
	}
    
}

void  poisson1d(double *u,int n1,int nx,double *sx,double *wx,double *work1,double *work2)
{
	int i ;

    sin1d(u,n1,nx,wx,work1,work2 ) ;
    
	for(i=1 ; i <= nx-1 ; i++){
		u[i] = u[i] / sx[i] ;
	}

    sin1d(u,n1,nx,wx,work1,work2 ) ;

}