2010年9月7日火曜日

LU分解

色々合って自分でLU分解しなければならなくなりました.

行列方程式の解を求めたいのですが,諸事情により,よく使ってた便利なOpenCVの関数を使えないということなのです.

いい勉強の機会だと思い,まずはLU分解の勉強.

線形代数とその応用という本で勉強したのですが,こいつがとんでもなく良い本ですね.

こういう語り口調でわかりやすい本は大好きです.

とりあえず仕組みは分かった.じゃあ分かりやすいソースコードがどこかに無いかなーと探すw

読みやすいソースコードを発見した.

神奈川大学の内田先生の公開コードを参照させていただきました.当該コードの著作権は内田先生に帰属します.参照元: http://www.inf.ie.kanagawa-u.ac.jp/c_learn/index.html#ouyou
よって私が改変した以下のコードはその二次的著作物となります.

上記コードではピボット選択をしない仕様でしたので,改変しできるようにしました.これで特異でない行列を含む行列方程式なら解を導出できることになります.

pivotが0の時のよさげなpivotを探して行を交換.それに伴いL行列も更新する.という操作を加えただけなんですが.

改変した主要な部分はLUpivot()とLUsolve()を加えた点です.うごかねーぞこら!ということに関して一切責任を負いかねます.もし参照するなら自己責任でお願いします.

いやはやエラー処理があまりに適当すぎるのが難点…

main.c
#include 
#include 
#include 
#include 
#include "LUutil.h"

int main( void )
{
int      n ; /* 次数 n : カウンタ i         */
double *A, *b ; /* 係数行列 A : 定数ベクトル b */
double *x; /* ベクトル x      */

n = get_dimension( MAX_DIM, MIN_DIM );           /* 次数を取得する */
if( ( A = get_matrix( n ) ) == (double *) NULL || /* 行列の A 領域   */
( b = get_vector( n ) ) == (double *) NULL || /* ベクトル b の領域 */
( x = get_vector( n ) ) == (double *) NULL )  /* ベクトル x の領域 */
{
fprintf( stderr, "メモリ取得エラーです!!\n" );
}
else
{
set_matrix( n, A );    /* 行列要素の入力      */
set_vector( n, b );    /* ベクトル要素の入力 */

// LU分解法によって行列方程式を解く
if( LUsolve(n, x, A, b) == -1 ) {
printf("LUsolve failure\n");
return -1;
}
MulMatandVec(b, A, x, n); /* 検算 */
printf( "LU2の解: " );
pr_vector( n, x ) ; /* 計算結果の出力 */
printf( "LU2の検算: " );
pr_vector( n, b ) ; /* 計算結果の出力 */
}

system("pause");
return EXIT_SUCCESS;
}


LUutil.h
#define MAX_DIM     10 /* 次数の最大値 */
#define MIN_DIM      1 /* 次数の最小値 */
#ifndef TRUE
#define TRUE      1
#endif
#ifndef FALSE
#define FALSE      0
#endif

void   MulMatandVec(double result[], double matA[], double vecB[], int size);
extern int      get_dimension( int max, int min );
double *get_matrix( int n );
double *get_vector( int n );
void   set_matrix( int n, double *a );
void   set_vector( int n, double *b );
void   mswap( double *a, double *b );
void   LU( int n, double *a, double *l, double *u );
int   LUsolve( int n, double *x, double *A, double *b );
int   LUpivot( int n, double *a, double *b, double *l, double *u );
void   L_equ( int n, double *l, double *c, double *b );
void   U_equ( int n, double *u, double *x, double *c );
void   pr_vector( int n, double *b );
void   matrixcpy( int n, double *u, double *a );
void   make_I( int n, double *l );


LUutil.c
#pragma once
#include 
#include 
#include 
#include 
#include "LUutil.h"

/* ----------------------------------------------------- *
* 行列とベクトルの積
*/
void MulMatandVec(double result[], double matA[], double vecB[], int size)
{
double sum;
int i, k;
for(i = 0; i < size; i++) {
sum = 0;
for(k = 0; k < size; k++)
sum += matA[i*size + k]*vecB[k];
result[i] = sum;
}
}
/* ------------------------------------------------------ *
*   行列・ベクトルの次数を取得する。
*   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;
}
}
}

/* ------------------------------------------------------ *
*    与えられたサイズの行列を格納する領域を確保する。
*    確保された領域の先頭アドレスを返す。
*    エラーならば、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 );
}


/* ------------------------------------------------------ *
*    与えられた行列に要素を格納する。
*/
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 );
scanf( "%lf", &a[ i*n + j ] );
}
}
}


/* ------------------------------------------------------ *
*    与えられたベクトルに要素を格納する。
*/
void set_vector( int n, double *b )
{
int i ; /* カウンタ */

for ( i = 0 ; i < n ; i ++ )
{
printf( "ベクトルの %d行の要素は ? >>>", i+1 );
scanf( "%lf", &b[ i ] );
}
}

/* ------------------------------------------------------ *
*    与えられた値の入れ替えを行う。
*/
void mswap( double *a, double *b )
{
double temp;

temp = *a;
*a = *b;
*b = temp;
}

/* ------------------------------------------------------ *
*   LU分解を行う(ピボットの選択は行わない)。
*   エラー時の動作は不定。
*/
void LU( int n, double *a, double *l, double *u )
{
int p, i, j;
double    c;

matrixcpy( n, (double *)u, (double *)a ); /* 行列 Uに 行列 Aをコピーする */
make_I( n, (double *)l );              /* 行列 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 ];
l[ i*n + p ] = c;
for ( j = p ; j < n ; j++ )
{
u[ i*n + j ] -= c * u[ p*n + j ];
}
}
}
}


int LUsolve( int n, double *x, double *A, double *b)
{
double   *L, *U ; /* 下三角行列 L : 上三角行列 U */
double   *c;

if(   ( L = get_matrix( n ) ) == (double *) NULL || /* 行列の L 領域   */
( U = get_matrix( n ) ) == (double *) NULL || /* 行列の U 領域   */
( c = get_vector( n ) ) == (double *) NULL ) /* ベクトル c の領域 */
{
fprintf( stderr, "メモリ取得エラーです!!\n" );
}

if( LUpivot( n, A, b, L, U) == -1) {    /* LU分解を行う */
system("pause");
return -1;
}
L_equ( n, L, c, b ); /* Lに関する連立方程式を解く */
U_equ( n, U, x, c ); /* Uに関する連立方程式を解く */
//   printf( "LU2の解は," );
//   pr_vector( n, x ); /* 計算結果の出力 */
return 0;
}

/* ------------------------------------------------------ *
*   LU分解を行う(ピボットの選択を行う)。
*/
int LUpivot( int n, double *a, double *b, double *l, double *u )
{
int      p, i, j, maxPos, ii;
int      pflg;
double   max = 0;
double   c;

matrixcpy( n, (double *)u, (double *)a ); /* 行列 Uに 行列 Aをコピーする */
make_I( n, (double *)l );              /* 行列 Lを単位行列にする       */

for ( p = 0 ; p < n ; p++ ) /* 掃き出し操作を行いながら */
{                     /* 行列を L,Uに分解する      */
maxPos = p;
pflg = FALSE;
// u[p][p]が0のとき,ピボット選択が必要
if( u[ p*n + p ] == 0) {
// permulation flug を立ち上げる
pflg = TRUE;
// 同一列(第p列)中の最大値(絶対値)を検索
for ( j = p + 1 ; j < n ; j++ ) {
max = 0;
if(fabs(u[ j*n + p ]) > max) {
max = u[ j*n + p ];
maxPos = j;
}
}
// もし最大値が0のままなら,Aは特異行列
if(max == 0) {
printf("Matrix A is singular! (in a column)\n");
return -1;
}
// 行列uの第p行と第maxPos行を交換する
for( j = p; j < n; j++)
mswap( &u[ p*n + j ], &u[ maxPos*n + j ]);

// 定数ベクトルの第i行と第maxPos行を交換する
mswap( &b[ p ], &b[ maxPos ]);

// L行列を更新する
// 行を交換したので,Lについても交換した行を更新する必要がある
for(ii = 0; ii < p; ii++) {
mswap( &l[maxPos*n + ii], &l[p*n + ii]);
}
}
for ( i = p + 1 ; i < n ; i++ )
{
c = u[ i*n + p ] / u[ p*n + p ];
l[ i*n + p ] = c;

for ( j = p ; j < n ; j++ )
{
u[ i*n + j ] -= c * u[ p*n + j ];
}
}
}
// もしu[n][n]が0なら,それは特異行列
if(u[ n*n - 1 ]==0) {
printf("Matrix A is singular! (in a row)\n");
return -1;
}
return 0;
}

/* ------------------------------------------------------ *
*   下三角行列・ベクトルで与えられた連立一次方程式を解く。
*/
void L_equ( int n, double *l, 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 ] * l[ 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 pr_vector( int n, double *b )
{
int i ; /* カウンタ */

for( i = 0 ; i < n ; i ++ )
{
printf ( " %g ", b[ i ] );
}
printf ( "\n" );
}

/* ------------------------------------------------------ *
*   行列のコピーを行う。
*/
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 *l )
{
int i, j ; /* カウンタ */

for ( i = 0 ; i < n ; i++ )      /* 行列の要素を 0にする */
{
for ( j = 0 ; j < n ; j++ )
{
l[ i*n + j ] = 0;
}
}
for ( i = 0 ; i < n ; i++ )      /* 対角要素を 1にする   */
{
l[ i*n + i ] = 1;
}
}


あぁ,何とかできた.てかこれくらいのコードもっとちゃちゃっとかけるようになりたいわ.

0 件のコメント:

コメントを投稿