codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
#include<stdio.h> #include<stdlib.h> #include<float.h> #define MAX_DIM 10 /* 次数の最大値 */ #define MIN_DIM 1 /* 次数の最小値 */ int get_dimension( int, int ); double *get_matrix( int ); double *get_vector( int ); void set_matrix( int, double* ); void set_vector( int, double* ); void LU( int, double*, double*, double* ); void L_equ( int, double*, double*, double* ); void U_equ( int, double*, double*, double* ); void pr_vector( int, double* ); void matrixcpy( int, double*, double* ); void make_I( int, double* ); void main( void ) { int n , i ; /* 次数 n : カウンタ i */ double *A, *b ; /* 係数行列 A : 定数ベクトル b */ double *L, *U ; /* 下三角行列 L : 上三角行列 U */ double *c, *x ; /* ベクトルc : ベクトル x */ n = get_dimension( MAX_DIM, MIN_DIM ); /* 次数を取得する */ if( ( A = get_matrix( n ) ) == (double *) NULL || /* 行列の A 領域 */ ( L = get_matrix( n ) ) == (double *) NULL || /* 行列の L 領域 */ ( U = get_matrix( n ) ) == (double *) NULL || /* 行列の U 領域 */ ( b = get_vector( n ) ) == (double *) NULL || /* ベクトル b 領域 */ ( c = get_vector( n ) ) == (double *) NULL || /* ベクトル c 領域 */ ( x = get_vector( n ) ) == (double *) NULL ) /* ベクトル x 領域 */ { fprintf( stderr, "メモリ取得エラーです!!\n" ); } else { set_matrix( n, A ); /* 行列要素の入力 */ set_vector( n, b ); /* ベクトル要素の入力 */ LU( n, A, L, U ); /* LU分解を行う */ L_equ( n, L, c, b ); /* Lに関する連立方程式を解く */ U_equ( n, U, x, c ); /* Uに関する連立方程式を解く */ printf( "解は, " ); pr_vector( n, x ) ; /* 計算結果の出力 */ } } /* ----------------------------------------------------------- * * LU分解を行う(ピボットの選択は行わない)。 * エラー時の動作は不定。 */ void LU( int n, double *a, double *I, double *u ) { int max, p , i, j; double c; matrixcpy( n, (double *)u, (double *)a ); /* 行列 Uに 行列 Aをコピーする */ make_I( n, (double *)I ); /* 行列 Lを単位行列にする */ for ( p = 0 ; p < n ; p++ ) /* 掃き出し操作を行いながら */ { /* 行列を L , Uに分解する */ for ( i = p + 1 ; i < n ; i++ ) { c = u[ i*n + p ] / u[ p*n + p ]; I[ i*n +p ] = c; for ( j = p ; j < n ; j++ ) { u[ i*n + j ] -= c * u[ p*n + j ]; } } } } /*----------------------------------------* * 下三角行列・ベクトルで与えられた連立一次方程式を解く。 */ void L_equ( int n, double *I, double *c, double *b ) { int i, j ; for ( i = 0 ; i < n ; i++ ) { c[ i ] = b[ i ]; for ( j= 0 ; j< i ; j++ ) { c[ i ] -= c[ j ] * I[ i*n + j ]; } } } /*------------------------------------------* * 上三角行列・ベクトルで与えられた連立一次方程式を解く * (後退代入) */ void U_equ( int n, double *u, double *x, double *c ) { int i, j; for ( i = n -1 ; i >= 0 ; i-- ) { x[ i ] = c[ i ]; for ( j = i + 1 ; j < n ; j++) { x[ i ] -= x[ j ] * u[ i*n + j ]; } x[ i ] /= u[ i*n + i ]; } } /*--------------------------------------------* * 行列のコピーを行う。 */ void matrixcpy( int n, double *u, double *a ) { int max, i; max = n * n ; for ( i= 0 ; i < max ; i++ ) /* 一次元配列としてコピー */ { u[ i ] = a[ i ]; } } /*---------------------------------------------* * 行列を単位行列化する。 */ void make_I( int n, double *I ) { int i, j ; /* カウンタ */ for ( i = 0 ; i < n ; i++ ) /* 行列の要素を 0にする */ { for ( j = 0 ; j < n ; j++ ) { I[ i*n + j] = 0; } } for ( i = 0 ; i < n ; i++ ) /* 対角要素を 1にする */ { I[ i*n + i ] = 1; } } /* -----------------------------------------------* * 与えられたサイズの行列を格納する領域を確保する。 * 確保された領域の先頭アドレスを返す。 * エラーならば、NULL を返す。 */ double *get_matrix( int n ) { return( ( double * ) malloc( sizeof( double ) * n * n ) ); } /* -----------------------------------------------* * 与えられたサイズのベクトルを格納する領域を確保する。 * 確保された領域の先頭アドレスを返す。 * エラーならば、NULL を返す。 */ double *get_vector( int n ) { return( ( double * ) malloc( sizeof ( double ) * n ) ); } /*--------------------------------------------------* * 行列・ベクトルの次元を取得する。 * max min で指定された範囲が入力されるまで、繰り返す。 */ int get_dimension( int max, int min ) { int d ; /* Dimension */ for ( ; ; ) { printf( "行列・ベクトルの次数を入力してください >>>" ); scanf( "%d", &d ); if ( d > max || d < min ) { fprintf(stderr, "その値は,許されていません。 \n" ); } else { return ( d ); } } } /* ---------------------------------------------------* * 与えられた行列に要素を格納する */ void set_matrix( int n, double *a ) { int i, j ; /* カウンタ */ for ( i = 0 ; i < n ; i++ ) { for ( j = 0 ; j < n ; j ++ ) { printf( "行列の %d行 %d列の要素は ? >>>", i+1 , j+1 ); } } } /* ----------------------------------------------------* * 与えられたベクトルに要素を格納する */ void set_vector( int n, double *b ) { int i ; /* カウンタ */ for ( i = 0 ; i < n ; i ++ ) { printf( "ベクトルの %d行の要素は ? >>>", i+1 ); scanf( "%1f", &b[ i ] ); } }
Private
[
?
]
Run code
Submit