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



void dot_crs_cpu(int n,double *b,double *val,double *x,int *row_ind,int *col_ind);

int main(void){
int n;
printf("please type dimention n.\n");
//今回のテストデータではn=5を入力
scanf("%d",&n);

double *matrix_a;
double *vector_b;
double *vector_x;

matrix_a = (double*)malloc(sizeof(double)*(n*n));
vector_b = (double*)malloc(sizeof(double)*n);
vector_x = (double*)malloc(sizeof(double)*n);

double *val;
val = (double*)malloc(sizeof(double)*(n*n));

static int *row_ind, *col_ind, *row_ptr;
row_ind = (int *)malloc( sizeof( int ) * (n*n) );
col_ind = (int *)malloc( sizeof( int ) * (n*n) );
row_ptr = (int *)malloc( sizeof( int ) * (n*n) );

double *p_val=val;
int *p_row_ind = row_ind,*p_col_ind = col_ind;
/*
p_row_ind = (int *)malloc( sizeof( int ) * (n*n) );
p_col_ind = (int *)malloc( sizeof( int ) * (n*n) );
*p_row_ind=row_ind,*p_col_ind=col_ind;
*/
int i, j, c;
FILE *fin,*fout;
char inputfile[50],outputfile[50];

    printf("please type input file name.\n");
	scanf("%s",inputfile);
	fin = fopen(inputfile ,"r");
    if(fin==NULL)
    {
        printf("no file\n");
		exit(1);
    }    

	printf("please type output file name.\n");
	scanf("%s",outputfile);
	fout = fopen(outputfile,"w");
    if(fout==NULL)
    {
        printf("do not maike file : output_sp.dat\n");
        exit(1);
    }

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

	fprintf( fout, "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");
    }
 
 
    fprintf( fout, "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("\nmatrix_a\n");
	for(i = 0; i < n; i ++)
    {
        for(j = 0; j < n; j ++)
        {
          printf("matrix_a[%d][%d] = %lf\n",i,j,matrix_a[n * i + j]);
		}
	}

	printf("\nvector_b\n");
	for(i = 0; i < n; i ++)
    {
          printf("matrix_b[%d] = %lf\n",i,vector_b[i]);
	}

	printf("\nvector_x\n");
	for(i = 0; i < n; i ++)
    {
          printf("matrix_x[%d] = %lf\n",i,vector_x[i]);
	}


	//matrix_aをCRS形式へ変換
	for( i = 0, c = 0 ; i < n ; i ++){
		for( j = 0; j < n; j ++){
			if(matrix_a[n * i + j]){
				*p_val++ = matrix_a[n * i + j];
				*p_row_ind++ = i+1;
				*p_col_ind++ = j+1;
				c++;
			}
		}
	}

	for(i=0;i<c;i++) 
	{
		if(row_ptr[row_ind[i]-1]==0) 
		{
				row_ptr[row_ind[i]-1]=i+1;
		}
	}

	//ターミナルへ出力
	printf("\nval ="); 
	for(i=0;i<c;i++){
			printf(" %lf", val[i]);
	}
	/*
	printf("\nrow_ind ="); 
	for(i=0;i<c;i++){
		printf(" %d", row_ind[i]);
	}
	*/

	printf("\ncol_ind ="); 
	for(i=0;i<c;i++){
			printf(" %d", col_ind[i]);
	}

	printf("\nrow_ptr ="); 
	for(i=0;i<n;i++){
		printf(" %d", row_ptr[i]);
	}


	dot_crs_cpu(n,vector_b,vector_x,val,col_ind,row_ptr);
	
	printf("\nvector_b\n");
	for(i = 0; i < n; i ++)
    {
          printf("matrix_b[%d] = %lf\n",i,vector_b[i]);
	}

	
	fclose(fin); 
    fclose(fout);

	printf("\n");

return 0;
}

/*test.dat
//matrix_a
//vector_b
//vector_x
1.000000 1.000000 1.000000 1.000000 1.000000 
1.000000 3.000000 3.000000 3.000000 3.000000 
1.000000 3.000000 5.000000 5.000000 5.000000 
1.000000 3.000000 5.000000 7.000000 7.000000 
1.000000 3.000000 5.000000 7.000000 9.000000 

5.000000 13.000000 19.000000 23.000000 25.000000 

0.000000 0.000000 0.000000 0.000000 0.000000 
 */


void dot_crs_cpu(int n,double *b,double *x,double *val,int *col_ind,int *row_ptr){
// Ax = b
int i,j,n;
for(i = 0 ; i < n ; i++) {
b[i] = 0;
  for(j = row_ptr[i] ; j < row_ptr[i + 1] ; j++) {
    b[i] += val[j] * x[ col_ind[j] ];
	printf("vector_b[%d] = %lf",n,b[i]);
  }
}
}