#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");

    cg_method(matrix_a, vector_x, 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(double *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(vector_ax, matrix_a, 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(vector_ap, matrix_a, vector_p);
        alpha = dot_product(vector_p, vector_r)/dot_product(vector_p, 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(vector_r);
        printf("LOOP : %d\t Error : %g\n", iter, err);
        if(EPS > err) break;

        beta = -dot_product(vector_r, vector_ap)/dot_product(vector_p, vector_ap);
        for(i = 0; i < N; i++){
            vector_p[i] = vector_r[i] + beta * vector_p[i];
        }
    }
}