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