
//	module		: conjugate gradient driver for 1-D poisson equation 
//	file name	: main.cpp
//	version		: 1.0
//	author		: Lung Sheng Chien
//	date		: 2004 / 1 / 24
//  
//	description : this is for CG test
//			Laplace U = f  on [0,1]
//			Dirichlet B.C
//  (1) parameters are hot coded
//	(2) output file 
//		fort.41 : u - u(exact) 
//		fort.42 : u 
//
//	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
//
//
//  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
//	3. implement FFT-based PCG
//


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "def.h"
#include "cg_1d_array.h"
#include "cg_1d_list.h"
#include "pcg_1d_list.h"
#include "math.h"
#include "testcase_1.h"



#define		_OUTFILE_DU		"fort.41"
#define		_OUTFILE_U		"fort.42"


linear_eqn *A ;
double *x ;
double *b ;
double *r ;  
double *v ;
double *z ;
double *u ;
double *f ;

double *r_hat ;

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

int imax ;
double eps ;
double delta ;



void  alloc_linear_sys(int nmax) ;
void  free_linear_sys(void) ;

void  display_error(int nx , double dx , FILE* fout1,FILE* fout2 ) ;

int main(int argc, char* argv[])
{
	int k ;
	int nx ;
	double dx ;
	FILE*  fout1 ;
	FILE*  fout2 ;
	char	filename[32] ;

//	strcpy(filename,_DEBUG_PATH);
	strcpy(filename,_OUTFILE_DU);
	fout1 = fopen(filename,"w+");
	assert(fout1);

//	strcpy(filename,_DEBUG_PATH);
	strcpy(filename,_OUTFILE_U);
	fout2 = fopen(filename,"w+");
	assert(fout2);

	alloc_linear_sys( COL_DIM ) ;

	nx = 256 ;
	assert( nx < COL_DIM ) ;
	dx = 1.0 / nx ;

	setupfft(COL_DIM,nx,dx,sx,wx) ;

	testcase_1_gen_operator( nx , dx , A ) ;
	testcase_1_gen_u_f(nx , dx , u, f ) ;
	testcase_1_setup_init( nx , dx , u , f, x, b) ;
	testcase_1_setup_bc( nx , dx , u , f , x , b) ;

	eps		= 1.0e-8 ;
	delta	= 1.0e-8 ;
	imax	= nx ;

	k = pcg_1d_list( COL_DIM , nx , A , x , b , r , r_hat, v , z, imax , eps , delta ) ;
//	k = cg_1d_list( COL_DIM , nx , A , x , b , r , v , z, imax , eps , delta ) ;

	display_error(nx , dx , fout1 , fout2 );
	printf("nx = %d , dx = %g \n",nx,dx);
	printf("eps = %g, delta = %g , imax = %d\n",eps,delta,imax);

	fclose(fout1);
	fclose(fout2);

	free_linear_sys() ;

	return 0 ;
}


void  display_error(int nx , double dx , FILE* fout1 , FILE* fout2 ) 
{
	int i ;
	double  errinf, errL2 , errL1 ;
	double  err ;

	errinf = 0 ;
	errL2 = 0 ;
	errL1 = 0 ;
	
	for(i=0 ; i <= nx ; i++){
		err = x[i] - u[i] ;
		fprintf(fout2,"%g\n",x[i]);
		fprintf(fout1,"%g\n",err);
		if ( 0.0 > err){
			err = -err ;
		}
		if ( err > errinf ){
			errinf = err ;
		}
//		printf("inf=%g\t",errinf);
		errL1 += err ;
		errL2 += err*err ;
	}
	errL2 = sqrt(errL2*dx) ;
	errL1 = errL1*dx ;
	printf("errinf = %g , errL1 = %g , errL2 = %g\n",errinf,errL1,errL2);
}

void  alloc_linear_sys(int nmax)
{
	A = (linear_eqn*) malloc(sizeof(linear_eqn)*nmax) ;
	assert(A) ;

	x = (double*) malloc(sizeof(double)*nmax);
	assert(x) ;
	b = (double*) malloc(sizeof(double)*nmax);
	assert(b) ;
	r = (double*) malloc(sizeof(double)*nmax);
	assert(r) ;
	v = (double*) malloc(sizeof(double)*nmax);
	assert(v) ;
	z = (double*) malloc(sizeof(double)*nmax);
	assert(z) ;
	u = (double*) malloc(sizeof(double)*nmax);
	assert(u) ;
	f = (double*) malloc(sizeof(double)*nmax);
	assert(f) ;

	r_hat = (double*) malloc(sizeof(double)*nmax);
	assert(r_hat) ;

	wx = (double*) malloc(sizeof(double)*nmax*3);
	assert(wx) ;

	sx = (double*) malloc(sizeof(double)*nmax);
	assert(sx) ;

	work1 = (double*) malloc(sizeof(double)*nmax);
	assert(work1) ;

	work2 = (double*) malloc(sizeof(double)*nmax);
	assert(work2) ;

}

void  free_linear_sys(void) 
{
	free(A) ;
	free(x) ;
	free(b) ;
	free(r) ;
	free(v) ;
	free(z) ;
	free(u) ;
	free(f) ;

	free(r_hat) ;

	free(wx) ;
	free(sx) ;
	free(work1) ;
	free(work2) ;


}