#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 件のコメント:
コメントを投稿