callocで確保したメモリ領域を2次元配列として扱う

まず，変数nxの数だけrandを使って未知数を作成する． 次に，やはりrandを使って連立方程式を作成し，callocで確保した領域に格納する． 後は，それをGauseの消去法で解く．

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

int ncol, nrow;

#define val(x,i,j)        x[ nrow * i + j ]     /* メモリ領域を2次元配列のように扱うマクロ */

/////////////////////////////////////////////////////////////////////////
int pivot( double *a, int i ){   /* ピボットの部分選択を行う関数 */

      double p = val( a, i, i ), t;
      int ii, j, k = i;

      for ( ii = i + 1; ii < ncol; ii++ ){
          t = fabs( val( a, ii, i ) );
          if ( t > p ){
              k = ii;
              p = t;
          }
      }
      if ( k != i ){
         for ( j = i; j < nrow; j++ ){
             t = val( a, k, j );
             val( a, k, j ) = val( a, i, j );
             val( a, i, j ) = t;
         }
      }
      return k;
}

/////////////////////////////////////////////////////////////////////////
int main(){

    int nx, i, ii, j;
    double *p, *x, t, sum;


    nx = 5;
    ncol = nx;
    nrow = nx + 1;

    srand( ( unsigned )time( NULL ) );  /* rand関数を初期化する */ 

    x = calloc( nx, sizeof( double ) ); /* nx * sizeof( double ) のメモリ領域を確保 */
    printf( "x values generated from function rand()\n" );
    for( i = 0; i < nx; i++ ) {
        x[ i ] = rand() / ( double )RAND_MAX; /* 確保した領域に適当な値を格納 */
        printf( "%12.7f", x[ i ] );
    }
    putchar( '\n' );

    p = calloc( ncol * nrow, sizeof( double ) ); /* 行数 * 列数 * sizeof( double ) のメモリ領域を確保 */
    for ( i = 0; i < ncol; i++ ){  /* 確保した領域に連立方程式の係数を代入 */
        for ( j = 0, sum = 0; j < nrow; j++ ){
            if ( j != nrow - 1 ){
               t = rand();
               val( p, i, j ) = t;
               sum += t * x[ j ];
            }
            else {
               val( p, i, j ) = sum;            
            }
       }
    }

    printf( "\ninitial matrix\n" );
    for ( i = 0; i < ncol; i++ ){
        for ( j = 0; j < nrow; j++ )
            printf( "%12.4e", val( p, i, j ) );
        putchar( '\n' );
    }

    for ( i = 0; i < ncol; i++ ){ /* ガウスの消去法の前進消去 */
       pivot( p, i );
       for ( ii = i + 1; ii < ncol; ii++ ){
           t = val( p, ii, i ) / val( p, i, i );
           for ( j = i; j < nrow; j++ ) 
               val( p, ii, j ) -= t * val( p, i, j );			
       }
    }

    printf( "\nmatrix after Gausian elimination\n" );
    for ( i = 0; i < ncol; i++ ){
        for ( j = 0; j < nrow; j++ )
            printf( "%12.4e", val( p, i, j ) );
        putchar( '\n' );
    }

    for ( i = ncol - 1; i >= 0; i-- ){ /* ガウスの消去法の後退代入 */
       for ( j = i + 1, sum = 0; j < ncol; j++ ) 
           sum += val( p, i, j ) * x[ j ];
       x[ i ] = ( val( p, i, nrow - 1 ) - sum ) / val( p, i, i );
    }

    printf( "\nx values obtained using Gausian elimination method\n" );
    for( i = 0; i < nx; i++ ) {
        printf( "%12.7f", x[ i ] );
    }
    putchar( '\n' );

    free( p );
    free( x );


}

戻る

SEO		[PR] 爆速!無料ブログ 無料ホームページ開設 無料ライブ放送