#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 10
#define TMAX 100
#define EPS (10e-6)

void init_vector();
double cg_method(double **matrix_a, double *vector_x, double *vector_b);
void init_vector(double *vector_x);
void vector_x_matrix(double *vector_y, double **matrix_a,double *vector_x);
double dot_product(double *vector_x, double *vector_y);
double vector_norm(double *vector_x);

int main (void){
FILE *fin, *fout;
static double **matrix_a;
static double *vector_b;
static double *vector_x;
int i,j;

matrix_a = (double**)malloc(sizeof(double*)*N);
	for(i=0;i<N;i++){
		matrix_a[i] = (double*)malloc(sizeof(double)*N);
	}

vector_b = (double*)malloc(sizeof(double)*N);
vector_x = (double*)malloc(sizeof(double)*N);

	if((fin = fopen("input_sp.dat","r"))==NULL)
	{
		printf("no file :input_sp.dat\n");
	}	
	if((fout = fopen("output_sp.dat","w"))==NULL)
	{
		printf("do not maike file : output_sp.dat\n");
		exit(1);
	}

fprintf( fout, "A\n");
printf("matrix_a\n");
	for( i = 0 ; i < N ; i++)
		{
			for( j = 0 ; j < N ; j++)
				{
						fscanf(fin, "%lf", &matrix_a[i][j]);
						fprintf(fout, "%5.2f\t", matrix_a[i][j]);
						printf("%lf\n",matrix_a[i][j]);
			}
			printf("\n");
			fprintf( fout, "\n");
	}

fprintf( fout, "b\n");
printf("vector_b\n");
	for( i = 0 ; i < N ; i++)
	{
		fscanf(fin, "%lf", &vector_b[i]);
		fprintf(fout, "%5.2f\t", vector_b[i]);
		fprintf( fout, "\n");
		printf("%lf\n",vector_b[i]);
	}

printf("\n");

fprintf( fout, "x\n");
printf("vector_x\n");
	for( i = 0 ; i < N ; i++)
	{
		fscanf(fin, "%lf", &vector_x[i]);
		fprintf(fout, "%5.2f\t", vector_x[i]);
		fprintf( fout, "\n");
		printf("%lf\n",vector_x[i]);
	}

printf("\n");

vector_x = cg_method(double **matrix_a,double *vector_x,double *vector_b);

fprintf( fout, "Ax=b is solved\n");
for( i = 1 ; i <= N ; i++ )
	{
			fprintf(fout, "%f\n", vector_x[i]);
	}
fclose(fin); 
fclose(fout);


return 0;
}

void init_vector(double *vector_x)
{
  int i;
  for(i = 0; i < N; i++){
    vector_x[i] = 0;
  }
}

void vector_x_matrix(double *vector_y,double *matrix_a,double *vector_x)
{
  int    i, j;
  double vxm;
  for(i = 0; i < N; i++){
    vxm = 0;
    for(j = 0; j < N; j++){
      vxm += matrix_a[i][j] * vector_x[i];
    }
    vector_y[i] = vxm;
  }
}

double dot_product(double *vector_x,double *vector_y)
{
  int    i;
  double dot_p = 0;
  for(i = 0; i < N; i++){
    dot_p += vector_x[i] * vector_y[i];
  }
  return dot_p;
}

double vector_norm(vector_x)
{
  int i;
  double norm = 0;
  for(i = 0; i < N; i++){
    norm += fabs(vector_x[i]);
  }
  return norm;
}

double cg_method(double **matrix_a, double *vector_x, double *vector_b)
{
  double *vector_p;
  double *vector_r;
  double *vector_ax;
  double *vector_ap;
  int i, iter;

  vector_p = (double*)malloc(sizeof(double)*N); 
  vector_r = (double*)malloc(sizeof(double)*N);
  vector_ax = (double*)malloc(sizeof(double)*N);
  vector_ap = (double*)malloc(sizeof(double)*N);
 
  vector_x_matrix(double *vector_ax, double *matrix_a, double *vector_x);

  for(i = 0; i < N; i++){
    vector_p[i] = vector_b[i] - vector_ax[i];
    vector_r[i] = vector_p[i];
  }

  for(iter = 1; iter < TMAX; iter++){
    double alpha, beta, err = 0;

    vector_x_matrix(double *vector_ap, double *matrix_a, double *vector_p);
    alpha = +dot_product(double *vector_p, *vector_r)/dot_product(double *vector_p, double *vector_ap);
    
    for(i = 0; i < N; i++){
      vector_x[i] += +alpha * vector_p[i];
      vector_r[i] += -alpha * vector_ap[i];
    }
    
    err = vector_norm(r);  
	printf("LOOP : %d\t Error : %g\n", iter, err);
    if(EPS > err) break;
	
	beta = -dot_product(double *vector_r, double *vector_ap)/dot_product(double *vector_p, double *vector_ap);
    for(i = 0; i < N; i++){
      vector_p[i] = vector_r[i] + beta * vector_p[i];
    } 
  }
}


-----------------------------------
input_sp.dat

1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 
1.000000 3.000000 3.000000 3.000000 3.000000 3.000000 3.000000 3.000000 3.000000 3.000000 
1.000000 3.000000 5.000000 5.000000 5.000000 5.000000 5.000000 5.000000 5.000000 5.000000 
1.000000 3.000000 5.000000 7.000000 7.000000 7.000000 7.000000 7.000000 7.000000 7.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 9.000000 9.000000 9.000000 9.000000 9.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 11.000000 11.000000 11.000000 11.000000 11.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 11.000000 13.000000 13.000000 13.000000 13.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 11.000000 13.000000 15.000000 15.000000 15.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 11.000000 13.000000 15.000000 17.000000 17.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 11.000000 13.000000 15.000000 17.000000 19.000000 

10.000000 28.000000 44.000000 58.000000 70.000000 80.000000 88.000000 94.000000 98.000000 100.000000 

0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000