#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を入力
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 " ) ; fin
= fopen ( inputfile
, "r" ) ; if ( fin== NULL)
{
}
printf ( "please type output file name.\n " ) ; fout
= fopen ( outputfile
, "w" ) ; if ( fout== NULL)
{
printf ( "do not maike file : output_sp.dat\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
] ) ; }
}
for ( i = 0 ; i < n ; i++ )
{
fscanf ( fin
, "%lf" , & vector_b
[ i
] ) ; fprintf ( fout
, "%5.2f\t " , vector_b
[ i
] ) ; }
for ( i = 0 ; i < n ; i++ )
{
fscanf ( fin
, "%lf" , & vector_x
[ i
] ) ; fprintf ( fout
, "%5.2f\t " , vector_x
[ i
] ) ; }
//ターミナルへ出力
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
] ) ; }
}
for ( i = 0 ; i < n; i ++ )
{
printf ( "matrix_b[%d] = %lf\n " , i
, vector_b
[ i
] ) ; }
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 ;
}
}
//ターミナルへ出力
for ( i= 0 ; i< c; i++ ) {
}
/*
printf("\nrow_ind =");
for(i=0;i<c;i++){
printf(" %d", row_ind[i]);
}
*/
for ( i= 0 ; i< c; i++ ) {
}
for ( i= 0 ; i< n; i++ ) {
}
dot_crs_cpu( n, vector_b, vector_x, val, col_ind, row_ptr) ;
for ( i = 0 ; i < n; i ++ )
{
printf ( "matrix_b[%d] = %lf\n " , i
, vector_b
[ i
] ) ; }
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=