#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 " ) ; }
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
] ) ; }
}
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
] ) ; }
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
] ) ; }
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
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KCiNkZWZpbmUgTiAxMAojZGVmaW5lIFRNQVggMTAwCiNkZWZpbmUgRVBTICgxMGUtNikKCnZvaWQgaW5pdF92ZWN0b3IoKTsKZG91YmxlIGNnX21ldGhvZChkb3VibGUgKiptYXRyaXhfYSwgZG91YmxlICp2ZWN0b3JfeCwgZG91YmxlICp2ZWN0b3JfYik7CnZvaWQgaW5pdF92ZWN0b3IoZG91YmxlICp2ZWN0b3JfeCk7CnZvaWQgdmVjdG9yX3hfbWF0cml4KGRvdWJsZSAqdmVjdG9yX3ksIGRvdWJsZSAqKm1hdHJpeF9hLGRvdWJsZSAqdmVjdG9yX3gpOwpkb3VibGUgZG90X3Byb2R1Y3QoZG91YmxlICp2ZWN0b3JfeCwgZG91YmxlICp2ZWN0b3JfeSk7CmRvdWJsZSB2ZWN0b3Jfbm9ybShkb3VibGUgKnZlY3Rvcl94KTsKCmludCBtYWluICh2b2lkKXsKRklMRSAqZmluLCAqZm91dDsKc3RhdGljIGRvdWJsZSAqKm1hdHJpeF9hOwpzdGF0aWMgZG91YmxlICp2ZWN0b3JfYjsKc3RhdGljIGRvdWJsZSAqdmVjdG9yX3g7CmludCBpLGo7CgptYXRyaXhfYSA9IChkb3VibGUqKiltYWxsb2Moc2l6ZW9mKGRvdWJsZSopKk4pOwoJZm9yKGk9MDtpPE47aSsrKXsKCQltYXRyaXhfYVtpXSA9IChkb3VibGUqKW1hbGxvYyhzaXplb2YoZG91YmxlKSpOKTsKCX0KCnZlY3Rvcl9iID0gKGRvdWJsZSopbWFsbG9jKHNpemVvZihkb3VibGUpKk4pOwp2ZWN0b3JfeCA9IChkb3VibGUqKW1hbGxvYyhzaXplb2YoZG91YmxlKSpOKTsKCglpZigoZmluID0gZm9wZW4oImlucHV0X3NwLmRhdCIsInIiKSk9PU5VTEwpCgl7CgkJcHJpbnRmKCJubyBmaWxlIDppbnB1dF9zcC5kYXRcbiIpOwoJfQkKCWlmKChmb3V0ID0gZm9wZW4oIm91dHB1dF9zcC5kYXQiLCJ3IikpPT1OVUxMKQoJewoJCXByaW50ZigiZG8gbm90IG1haWtlIGZpbGUgOiBvdXRwdXRfc3AuZGF0XG4iKTsKCQlleGl0KDEpOwoJfQoKZnByaW50ZiggZm91dCwgIkFcbiIpOwpwcmludGYoIm1hdHJpeF9hXG4iKTsKCWZvciggaSA9IDAgOyBpIDwgTiA7IGkrKykKCQl7CgkJCWZvciggaiA9IDAgOyBqIDwgTiA7IGorKykKCQkJCXsKCQkJCQkJZnNjYW5mKGZpbiwgIiVsZiIsICZtYXRyaXhfYVtpXVtqXSk7CgkJCQkJCWZwcmludGYoZm91dCwgIiU1LjJmXHQiLCBtYXRyaXhfYVtpXVtqXSk7CgkJCQkJCXByaW50ZigiJWxmXG4iLG1hdHJpeF9hW2ldW2pdKTsKCQkJfQoJCQlwcmludGYoIlxuIik7CgkJCWZwcmludGYoIGZvdXQsICJcbiIpOwoJfQoKZnByaW50ZiggZm91dCwgImJcbiIpOwpwcmludGYoInZlY3Rvcl9iXG4iKTsKCWZvciggaSA9IDAgOyBpIDwgTiA7IGkrKykKCXsKCQlmc2NhbmYoZmluLCAiJWxmIiwgJnZlY3Rvcl9iW2ldKTsKCQlmcHJpbnRmKGZvdXQsICIlNS4yZlx0IiwgdmVjdG9yX2JbaV0pOwoJCWZwcmludGYoIGZvdXQsICJcbiIpOwoJCXByaW50ZigiJWxmXG4iLHZlY3Rvcl9iW2ldKTsKCX0KCnByaW50ZigiXG4iKTsKCmZwcmludGYoIGZvdXQsICJ4XG4iKTsKcHJpbnRmKCJ2ZWN0b3JfeFxuIik7Cglmb3IoIGkgPSAwIDsgaSA8IE4gOyBpKyspCgl7CgkJZnNjYW5mKGZpbiwgIiVsZiIsICZ2ZWN0b3JfeFtpXSk7CgkJZnByaW50Zihmb3V0LCAiJTUuMmZcdCIsIHZlY3Rvcl94W2ldKTsKCQlmcHJpbnRmKCBmb3V0LCAiXG4iKTsKCQlwcmludGYoIiVsZlxuIix2ZWN0b3JfeFtpXSk7Cgl9CgpwcmludGYoIlxuIik7Cgp2ZWN0b3JfeCA9IGNnX21ldGhvZChkb3VibGUgKiptYXRyaXhfYSxkb3VibGUgKnZlY3Rvcl94LGRvdWJsZSAqdmVjdG9yX2IpOwoKZnByaW50ZiggZm91dCwgIkF4PWIgaXMgc29sdmVkXG4iKTsKZm9yKCBpID0gMSA7IGkgPD0gTiA7IGkrKyApCgl7CgkJCWZwcmludGYoZm91dCwgIiVmXG4iLCB2ZWN0b3JfeFtpXSk7Cgl9CmZjbG9zZShmaW4pOyAKZmNsb3NlKGZvdXQpOwoKCnJldHVybiAwOwp9Cgp2b2lkIGluaXRfdmVjdG9yKGRvdWJsZSAqdmVjdG9yX3gpCnsKICBpbnQgaTsKICBmb3IoaSA9IDA7IGkgPCBOOyBpKyspewogICAgdmVjdG9yX3hbaV0gPSAwOwogIH0KfQoKdm9pZCB2ZWN0b3JfeF9tYXRyaXgoZG91YmxlICp2ZWN0b3JfeSxkb3VibGUgKm1hdHJpeF9hLGRvdWJsZSAqdmVjdG9yX3gpCnsKICBpbnQgICAgaSwgajsKICBkb3VibGUgdnhtOwogIGZvcihpID0gMDsgaSA8IE47IGkrKyl7CiAgICB2eG0gPSAwOwogICAgZm9yKGogPSAwOyBqIDwgTjsgaisrKXsKICAgICAgdnhtICs9IG1hdHJpeF9hW2ldW2pdICogdmVjdG9yX3hbaV07CiAgICB9CiAgICB2ZWN0b3JfeVtpXSA9IHZ4bTsKICB9Cn0KCmRvdWJsZSBkb3RfcHJvZHVjdChkb3VibGUgKnZlY3Rvcl94LGRvdWJsZSAqdmVjdG9yX3kpCnsKICBpbnQgICAgaTsKICBkb3VibGUgZG90X3AgPSAwOwogIGZvcihpID0gMDsgaSA8IE47IGkrKyl7CiAgICBkb3RfcCArPSB2ZWN0b3JfeFtpXSAqIHZlY3Rvcl95W2ldOwogIH0KICByZXR1cm4gZG90X3A7Cn0KCmRvdWJsZSB2ZWN0b3Jfbm9ybSh2ZWN0b3JfeCkKewogIGludCBpOwogIGRvdWJsZSBub3JtID0gMDsKICBmb3IoaSA9IDA7IGkgPCBOOyBpKyspewogICAgbm9ybSArPSBmYWJzKHZlY3Rvcl94W2ldKTsKICB9CiAgcmV0dXJuIG5vcm07Cn0KCmRvdWJsZSBjZ19tZXRob2QoZG91YmxlICoqbWF0cml4X2EsIGRvdWJsZSAqdmVjdG9yX3gsIGRvdWJsZSAqdmVjdG9yX2IpCnsKICBkb3VibGUgKnZlY3Rvcl9wOwogIGRvdWJsZSAqdmVjdG9yX3I7CiAgZG91YmxlICp2ZWN0b3JfYXg7CiAgZG91YmxlICp2ZWN0b3JfYXA7CiAgaW50IGksIGl0ZXI7CgogIHZlY3Rvcl9wID0gKGRvdWJsZSopbWFsbG9jKHNpemVvZihkb3VibGUpKk4pOyAKICB2ZWN0b3JfciA9IChkb3VibGUqKW1hbGxvYyhzaXplb2YoZG91YmxlKSpOKTsKICB2ZWN0b3JfYXggPSAoZG91YmxlKiltYWxsb2Moc2l6ZW9mKGRvdWJsZSkqTik7CiAgdmVjdG9yX2FwID0gKGRvdWJsZSopbWFsbG9jKHNpemVvZihkb3VibGUpKk4pOwogCiAgdmVjdG9yX3hfbWF0cml4KGRvdWJsZSAqdmVjdG9yX2F4LCBkb3VibGUgKm1hdHJpeF9hLCBkb3VibGUgKnZlY3Rvcl94KTsKCiAgZm9yKGkgPSAwOyBpIDwgTjsgaSsrKXsKICAgIHZlY3Rvcl9wW2ldID0gdmVjdG9yX2JbaV0gLSB2ZWN0b3JfYXhbaV07CiAgICB2ZWN0b3JfcltpXSA9IHZlY3Rvcl9wW2ldOwogIH0KCiAgZm9yKGl0ZXIgPSAxOyBpdGVyIDwgVE1BWDsgaXRlcisrKXsKICAgIGRvdWJsZSBhbHBoYSwgYmV0YSwgZXJyID0gMDsKCiAgICB2ZWN0b3JfeF9tYXRyaXgoZG91YmxlICp2ZWN0b3JfYXAsIGRvdWJsZSAqbWF0cml4X2EsIGRvdWJsZSAqdmVjdG9yX3ApOwogICAgYWxwaGEgPSArZG90X3Byb2R1Y3QoZG91YmxlICp2ZWN0b3JfcCwgKnZlY3Rvcl9yKS9kb3RfcHJvZHVjdChkb3VibGUgKnZlY3Rvcl9wLCBkb3VibGUgKnZlY3Rvcl9hcCk7CiAgICAKICAgIGZvcihpID0gMDsgaSA8IE47IGkrKyl7CiAgICAgIHZlY3Rvcl94W2ldICs9ICthbHBoYSAqIHZlY3Rvcl9wW2ldOwogICAgICB2ZWN0b3JfcltpXSArPSAtYWxwaGEgKiB2ZWN0b3JfYXBbaV07CiAgICB9CiAgICAKICAgIGVyciA9IHZlY3Rvcl9ub3JtKHIpOyAgCglwcmludGYoIkxPT1AgOiAlZFx0IEVycm9yIDogJWdcbiIsIGl0ZXIsIGVycik7CiAgICBpZihFUFMgPiBlcnIpIGJyZWFrOwoJCgliZXRhID0gLWRvdF9wcm9kdWN0KGRvdWJsZSAqdmVjdG9yX3IsIGRvdWJsZSAqdmVjdG9yX2FwKS9kb3RfcHJvZHVjdChkb3VibGUgKnZlY3Rvcl9wLCBkb3VibGUgKnZlY3Rvcl9hcCk7CiAgICBmb3IoaSA9IDA7IGkgPCBOOyBpKyspewogICAgICB2ZWN0b3JfcFtpXSA9IHZlY3Rvcl9yW2ldICsgYmV0YSAqIHZlY3Rvcl9wW2ldOwogICAgfSAKICB9Cn0KCgotLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQppbnB1dF9zcC5kYXQKCjEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIDEuMDAwMDAwIAoxLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAzLjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgNS4wMDAwMDAgCjEuMDAwMDAwIDMuMDAwMDAwIDUuMDAwMDAwIDcuMDAwMDAwIDcuMDAwMDAwIDcuMDAwMDAwIDcuMDAwMDAwIDcuMDAwMDAwIDcuMDAwMDAwIDcuMDAwMDAwIAoxLjAwMDAwMCAzLjAwMDAwMCA1LjAwMDAwMCA3LjAwMDAwMCA5LjAwMDAwMCA5LjAwMDAwMCA5LjAwMDAwMCA5LjAwMDAwMCA5LjAwMDAwMCA5LjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNy4wMDAwMDAgOS4wMDAwMDAgMTEuMDAwMDAwIDExLjAwMDAwMCAxMS4wMDAwMDAgMTEuMDAwMDAwIDExLjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNy4wMDAwMDAgOS4wMDAwMDAgMTEuMDAwMDAwIDEzLjAwMDAwMCAxMy4wMDAwMDAgMTMuMDAwMDAwIDEzLjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNy4wMDAwMDAgOS4wMDAwMDAgMTEuMDAwMDAwIDEzLjAwMDAwMCAxNS4wMDAwMDAgMTUuMDAwMDAwIDE1LjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNy4wMDAwMDAgOS4wMDAwMDAgMTEuMDAwMDAwIDEzLjAwMDAwMCAxNS4wMDAwMDAgMTcuMDAwMDAwIDE3LjAwMDAwMCAKMS4wMDAwMDAgMy4wMDAwMDAgNS4wMDAwMDAgNy4wMDAwMDAgOS4wMDAwMDAgMTEuMDAwMDAwIDEzLjAwMDAwMCAxNS4wMDAwMDAgMTcuMDAwMDAwIDE5LjAwMDAwMCAKCjEwLjAwMDAwMCAyOC4wMDAwMDAgNDQuMDAwMDAwIDU4LjAwMDAwMCA3MC4wMDAwMDAgODAuMDAwMDAwIDg4LjAwMDAwMCA5NC4wMDAwMDAgOTguMDAwMDAwIDEwMC4wMDAwMDAgCgowLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAwLjAwMDAwMCAK