#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]);
  }
}
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KCgoKdm9pZCBkb3RfY3JzX2NwdShpbnQgbixkb3VibGUgKmIsZG91YmxlICp2YWwsZG91YmxlICp4LGludCAqcm93X2luZCxpbnQgKmNvbF9pbmQpOwoKaW50IG1haW4odm9pZCl7CmludCBuOwpwcmludGYoInBsZWFzZSB0eXBlIGRpbWVudGlvbiBuLlxuIik7Ci8v5LuK5Zue44Gu44OG44K544OI44OH44O844K/44Gn44Gvbj0144KS5YWl5YqbCnNjYW5mKCIlZCIsJm4pOwoKZG91YmxlICptYXRyaXhfYTsKZG91YmxlICp2ZWN0b3JfYjsKZG91YmxlICp2ZWN0b3JfeDsKCm1hdHJpeF9hID0gKGRvdWJsZSopbWFsbG9jKHNpemVvZihkb3VibGUpKihuKm4pKTsKdmVjdG9yX2IgPSAoZG91YmxlKiltYWxsb2Moc2l6ZW9mKGRvdWJsZSkqbik7CnZlY3Rvcl94ID0gKGRvdWJsZSopbWFsbG9jKHNpemVvZihkb3VibGUpKm4pOwoKZG91YmxlICp2YWw7CnZhbCA9IChkb3VibGUqKW1hbGxvYyhzaXplb2YoZG91YmxlKSoobipuKSk7CgpzdGF0aWMgaW50ICpyb3dfaW5kLCAqY29sX2luZCwgKnJvd19wdHI7CnJvd19pbmQgPSAoaW50ICopbWFsbG9jKCBzaXplb2YoIGludCApICogKG4qbikgKTsKY29sX2luZCA9IChpbnQgKiltYWxsb2MoIHNpemVvZiggaW50ICkgKiAobipuKSApOwpyb3dfcHRyID0gKGludCAqKW1hbGxvYyggc2l6ZW9mKCBpbnQgKSAqIChuKm4pICk7Cgpkb3VibGUgKnBfdmFsPXZhbDsKaW50ICpwX3Jvd19pbmQgPSByb3dfaW5kLCpwX2NvbF9pbmQgPSBjb2xfaW5kOwovKgpwX3Jvd19pbmQgPSAoaW50ICopbWFsbG9jKCBzaXplb2YoIGludCApICogKG4qbikgKTsKcF9jb2xfaW5kID0gKGludCAqKW1hbGxvYyggc2l6ZW9mKCBpbnQgKSAqIChuKm4pICk7CipwX3Jvd19pbmQ9cm93X2luZCwqcF9jb2xfaW5kPWNvbF9pbmQ7CiovCmludCBpLCBqLCBjOwpGSUxFICpmaW4sKmZvdXQ7CmNoYXIgaW5wdXRmaWxlWzUwXSxvdXRwdXRmaWxlWzUwXTsKCiAgICBwcmludGYoInBsZWFzZSB0eXBlIGlucHV0IGZpbGUgbmFtZS5cbiIpOwoJc2NhbmYoIiVzIixpbnB1dGZpbGUpOwoJZmluID0gZm9wZW4oaW5wdXRmaWxlICwiciIpOwogICAgaWYoZmluPT1OVUxMKQogICAgewogICAgICAgIHByaW50Zigibm8gZmlsZVxuIik7CgkJZXhpdCgxKTsKICAgIH0gICAgCgoJcHJpbnRmKCJwbGVhc2UgdHlwZSBvdXRwdXQgZmlsZSBuYW1lLlxuIik7CglzY2FuZigiJXMiLG91dHB1dGZpbGUpOwoJZm91dCA9IGZvcGVuKG91dHB1dGZpbGUsInciKTsKICAgIGlmKGZvdXQ9PU5VTEwpCiAgICB7CiAgICAgICAgcHJpbnRmKCJkbyBub3QgbWFpa2UgZmlsZSA6IG91dHB1dF9zcC5kYXRcbiIpOwogICAgICAgIGV4aXQoMSk7CiAgICB9CgoJZnByaW50ZiggZm91dCwgIkFcbiIpOwogICAgZm9yKCBpID0gMCA7IGkgPCBuIDsgaSsrKQogICAgewogICAgICAgIGZvciggaiA9IDAgOyBqIDwgbiA7IGorKykKICAgICAgICB7CiAgICAgICAgICAgIGZzY2FuZihmaW4sICIlbGYiLCAmbWF0cml4X2FbbippK2pdKTsKICAgICAgICAgICAgZnByaW50Zihmb3V0LCAiJTUuMmZcdCIsIG1hdHJpeF9hW24qaStqXSk7CiAgICAgICAgfQogICAgICAgIGZwcmludGYoIGZvdXQsICJcbiIpOwogICAgfQoKCWZwcmludGYoIGZvdXQsICJiXG4iKTsKICAgIGZvciggaSA9IDAgOyBpIDwgbiA7IGkrKykKICAgIHsKICAgICAgICBmc2NhbmYoZmluLCAiJWxmIiwgJnZlY3Rvcl9iW2ldKTsKICAgICAgICBmcHJpbnRmKGZvdXQsICIlNS4yZlx0IiwgdmVjdG9yX2JbaV0pOwogICAgICAgIGZwcmludGYoIGZvdXQsICJcbiIpOwogICAgfQogCiAKICAgIGZwcmludGYoIGZvdXQsICJ4XG4iKTsKICAgIGZvciggaSA9IDAgOyBpIDwgbiA7IGkrKykKICAgIHsKICAgICAgICBmc2NhbmYoZmluLCAiJWxmIiwgJnZlY3Rvcl94W2ldKTsKICAgICAgICBmcHJpbnRmKGZvdXQsICIlNS4yZlx0IiwgdmVjdG9yX3hbaV0pOwogICAgICAgIGZwcmludGYoIGZvdXQsICJcbiIpOwogICAgfQoKCS8v44K/44O844Of44OK44Or44G45Ye65YqbCglwcmludGYoIlxubWF0cml4X2FcbiIpOwoJZm9yKGkgPSAwOyBpIDwgbjsgaSArKykKICAgIHsKICAgICAgICBmb3IoaiA9IDA7IGogPCBuOyBqICsrKQogICAgICAgIHsKICAgICAgICAgIHByaW50ZigibWF0cml4X2FbJWRdWyVkXSA9ICVsZlxuIixpLGosbWF0cml4X2FbbiAqIGkgKyBqXSk7CgkJfQoJfQoKCXByaW50ZigiXG52ZWN0b3JfYlxuIik7Cglmb3IoaSA9IDA7IGkgPCBuOyBpICsrKQogICAgewogICAgICAgICAgcHJpbnRmKCJtYXRyaXhfYlslZF0gPSAlbGZcbiIsaSx2ZWN0b3JfYltpXSk7Cgl9CgoJcHJpbnRmKCJcbnZlY3Rvcl94XG4iKTsKCWZvcihpID0gMDsgaSA8IG47IGkgKyspCiAgICB7CiAgICAgICAgICBwcmludGYoIm1hdHJpeF94WyVkXSA9ICVsZlxuIixpLHZlY3Rvcl94W2ldKTsKCX0KCgoJLy9tYXRyaXhfYeOCkkNSU+W9ouW8j+OBuOWkieaPmwoJZm9yKCBpID0gMCwgYyA9IDAgOyBpIDwgbiA7IGkgKyspewoJCWZvciggaiA9IDA7IGogPCBuOyBqICsrKXsKCQkJaWYobWF0cml4X2FbbiAqIGkgKyBqXSl7CgkJCQkqcF92YWwrKyA9IG1hdHJpeF9hW24gKiBpICsgal07CgkJCQkqcF9yb3dfaW5kKysgPSBpKzE7CgkJCQkqcF9jb2xfaW5kKysgPSBqKzE7CgkJCQljKys7CgkJCX0KCQl9Cgl9CgoJZm9yKGk9MDtpPGM7aSsrKSAKCXsKCQlpZihyb3dfcHRyW3Jvd19pbmRbaV0tMV09PTApIAoJCXsKCQkJCXJvd19wdHJbcm93X2luZFtpXS0xXT1pKzE7CgkJfQoJfQoKCS8v44K/44O844Of44OK44Or44G45Ye65YqbCglwcmludGYoIlxudmFsID0iKTsgCglmb3IoaT0wO2k8YztpKyspewoJCQlwcmludGYoIiAlbGYiLCB2YWxbaV0pOwoJfQoJLyoKCXByaW50ZigiXG5yb3dfaW5kID0iKTsgCglmb3IoaT0wO2k8YztpKyspewoJCXByaW50ZigiICVkIiwgcm93X2luZFtpXSk7Cgl9CgkqLwoKCXByaW50ZigiXG5jb2xfaW5kID0iKTsgCglmb3IoaT0wO2k8YztpKyspewoJCQlwcmludGYoIiAlZCIsIGNvbF9pbmRbaV0pOwoJfQoKCXByaW50ZigiXG5yb3dfcHRyID0iKTsgCglmb3IoaT0wO2k8bjtpKyspewoJCXByaW50ZigiICVkIiwgcm93X3B0cltpXSk7Cgl9CgoKCWRvdF9jcnNfY3B1KG4sdmVjdG9yX2IsdmVjdG9yX3gsdmFsLGNvbF9pbmQscm93X3B0cik7CgkKCXByaW50ZigiXG52ZWN0b3JfYlxuIik7Cglmb3IoaSA9IDA7IGkgPCBuOyBpICsrKQogICAgewogICAgICAgICAgcHJpbnRmKCJtYXRyaXhfYlslZF0gPSAlbGZcbiIsaSx2ZWN0b3JfYltpXSk7Cgl9CgoJCglmY2xvc2UoZmluKTsgCiAgICBmY2xvc2UoZm91dCk7CgoJcHJpbnRmKCJcbiIpOwoKcmV0dXJuIDA7Cn0KCi8qdGVzdC5kYXQKLy9tYXRyaXhfYQovL3ZlY3Rvcl9iCi8vdmVjdG9yX3gKMS4wMDAwMDAgMS4wMDAwMDAgMS4wMDAwMDAgMS4wMDAwMDAgMS4wMDAwMDAgCjEuMDAwMDAwIDMuMDAwMDAwIDMuMDAwMDAwIDMuMDAwMDAwIDMuMDAwMDAwIAoxLjAwMDAwMCAzLjAwMDAwMCA1LjAwMDAwMCA1LjAwMDAwMCA1LjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNy4wMDAwMDAgNy4wMDAwMDAgCjEuMDAwMDAwIDMuMDAwMDAwIDUuMDAwMDAwIDcuMDAwMDAwIDkuMDAwMDAwIAoKNS4wMDAwMDAgMTMuMDAwMDAwIDE5LjAwMDAwMCAyMy4wMDAwMDAgMjUuMDAwMDAwIAoKMC4wMDAwMDAgMC4wMDAwMDAgMC4wMDAwMDAgMC4wMDAwMDAgMC4wMDAwMDAgCiAqLwoKCnZvaWQgZG90X2Nyc19jcHUoaW50IG4sZG91YmxlICpiLGRvdWJsZSAqeCxkb3VibGUgKnZhbCxpbnQgKmNvbF9pbmQsaW50ICpyb3dfcHRyKXsKLy8gQXggPSBiCmludCBpLGosbjsKZm9yKGkgPSAwIDsgaSA8IG4gOyBpKyspIHsKYltpXSA9IDA7CiAgZm9yKGogPSByb3dfcHRyW2ldIDsgaiA8IHJvd19wdHJbaSArIDFdIDsgaisrKSB7CiAgICBiW2ldICs9IHZhbFtqXSAqIHhbIGNvbF9pbmRbal0gXTsKCXByaW50ZigidmVjdG9yX2JbJWRdID0gJWxmIixuLGJbaV0pOwogIH0KfQp9Cgo=