#include "header.h"

int CG(gsl_matrix_float *A,gsl_vector_float *b,gsl_vector_float *x,int n)
{	

	float temp1;
	float  rw,wz,rz,rw1,wz1;
	double alpha0= -1.0,alpha1=1.0;
	double beta=0.0;
	float a=0.0,B=0.0;
	int i;
	gsl_vector_float *temp,*w,*z,*temp2,*temp3,*temp4,*temp5;
	temp=gsl_vector_float_calloc((*b).size);
	w=gsl_vector_float_calloc((*b).size);
	z=gsl_vector_float_calloc((*b).size);
	temp4=gsl_vector_float_calloc((*b).size);
	temp5=gsl_vector_float_calloc((*b).size);
	
	temp2=gsl_vector_float_calloc((*b).size);
	temp3=gsl_vector_float_calloc((*b).size);	
	gsl_blas_sgemv(CblasNoTrans,alpha0,A,x,beta,temp);//compute -Ax, store it value in temp.
	
	gsl_blas_saxpy(alpha1,b,temp); //get b-Ax,store this value in temp(r);
	gsl_blas_scopy(temp,w); //copy r to w.
	gsl_blas_sscal(alpha0,w); //get w=-r, store the value in w.
	gsl_blas_sgemv(CblasNoTrans,alpha1,A,w,beta,z); //get z=Aw, store it in z.
	gsl_blas_sdot(temp,w,&rw); // get r'w. store it in rw.
	gsl_blas_sdot(w,z,&wz);	// get w'z. store it in wz.
        a=(rw)/(wz);
	gsl_blas_scopy(w,temp2); //copy w to temp2,temp2=w;
	gsl_blas_sscal(a,temp2); //temp2=a*wgsl_blas_sscal(a,temp3); //get a*z, store in temp3=a*z;
	
	gsl_vector_float_add(x,temp2); //x=x+a*w
	//gsl_vector_float_fprintf(stdout,z,"%lf");
	//printf("a=%lf",a);
	for (i=0;i<=n;i++) 
	{
		gsl_blas_scopy(z,temp3); //copy z to temp3,temp3=z;
		
		gsl_blas_sscal(a,temp3); //get a*z, store in temp3=a*z;
		
		
		gsl_vector_float_sub(temp,temp3);// r=r-a*z, store in temp
		
		if (gsl_blas_snrm2(temp)<qloop)
			break;
		gsl_blas_sdot(temp,z,&rz); //get r'*z,store in rz
		B=(rz)/(wz);
		
		gsl_blas_scopy(w,temp4); //copy w to temp4,temp4=w;
		
		gsl_blas_sscal(B,temp4); //temp4=B*w
		
		gsl_blas_scopy(temp,temp5);//copy temp(r) to temp5. temp5=r
		gsl_blas_sscal(alpha0,temp5); //get temp5=-r, store the value in temp5.
		
		gsl_vector_float_add(temp5,temp4); //temp5=temp5+temp4, i.e temp5=-r+B*w
		
		gsl_blas_scopy(temp5,w); //copy temp5 to w,w=-r+B*w;
		
		
		gsl_blas_sgemv(CblasNoTrans,alpha1,A,w,beta,z); //get z=Aw, store it in z.
		
		gsl_blas_sdot(temp,w,&rw1); // get r'w. store it in rw1.
		
		gsl_blas_sdot(w,z,&wz1); // get r'w. store it in wz1.
		
		a=(rw1)/(wz1);
		
		//gsl_vector_float_fprintf(stdout,temp5,"%lf");
		gsl_blas_sscal(a,temp5); //get a*z, store in temp5=a*w;
		
		gsl_vector_float_add(x,temp5); //x=x+a*w
    	}
	return 0;
}