-
プログラムを書くときは、すべてを１つの main(){ に書かないで、
-
LU分解するモジュール
( L )(y) = (b)を解くモジュール
( U )(x) = (y)を解くモジュール
を関数とかで別に作っておいて、（子会社化）
-
main() (親会社）はそれを呼び出して使うだけ、
という感じに作りましょう。
モジュールごとに検査が出てきて、どこが悪いか調べやすいし、
別のモジュールに差し替えたりできるので、
汎用性も保守性も上がります。
main()
が短くなって、全体の流れがわかりやすくなることは、
将来、他の人と仕事をしたり、仕事を引き継いだりするときに役立ちます。
2次元配列の受け渡しが苦手な方は、
2次元配列だけmain(){ の前に型宣言して
グローバル変数にしてしまうという手もあります。
検算は手動でやってもＯＫですが
自動でやってくれるようにプログラムしておくと後々楽です。
なんかこんな感じ
-
#include < stdio.h >
#define N 4 .... 何元連立方程式か、問題によって変更できるように
void lu_decomposition(int, double a[][N+1],double l[][N+1], double u[][N+1]);
void solve_l(int, double l[][N+1], double b[], double x[]);
void solve_u(int, double u[][N+1], double b[], double x[]); ....プロトタイプ宣言
int main(){
-
int n,i,j; ....型宣言 このほかにも必要なら自分で宣言
double a[N+1][N+1], l[N+1][N+1], u[N+1][N+1]; ....a[N][N]まで使うなら配列のサイズは[N+1][N+1]
double b[N+1], x[N+1], y[N+1]; ....b[N]まで使うなら配列のサイズは[N+1]、使わなくてもb[0]があるため
n = N;
-
/*---- 初期化 (0番は使わないですが)---*/
for(i=0; i<=n; i++){
-
for(j=0; j<=n; j++){
-
$l$[i][j] = 0.0;
u[i][j] = 0.0;
}
$l$[i][i] = 1.0;
}
/*---- 係数a[][] と右辺 b[] を入力したり代入したりする ---*/
：
自作してね
：
-
/*---AをLUに分解する ---*/
lu_decomposition(n,a,l,u);
/*---ここでLとUを出力して確認するとよい---*/
/*--- Ly=b を解く (Lとbからyを求める）---*/
solve_l(n,l,b,y);
/*---ここでyを出力して確認するとよい---*/
/*--- Ux=y を解く (Uとyからxを求める）---*/
solve_u(n,u,y,x);
/*---ここで 解xを出力 ---*/
：
自作してね
：
/*---検算する(オプション）--*/
：
自作してね
：
-
}
void lu_decomposition(int n, double a[][N+1], double $l$[][N+1],
double u[][N+1] ){
/*--- LU分解する関数---*/
-
int i,j,k,m;
:
自作してね
：
-
}
void solve_l(int n, double l[][N+1], double b[], double x[]){
/*--- Lx=b を解く関数---*/
-
int i,k;
:
自作してね
：
-
}
void solve_u(int n, double u[][N+1], double b[], double x[]){
/*--- Ux=b を解く関数---*/
-
int i,k;
:
自作してね
：
-
}
void product_matrix(int n, double a[][N+1], double x[], double z[]){
/*---検算のためz= Ux を計算する関数---*/
：
(オプション）手計算で検算してもいいけど作っておくと楽だよね
：
-
}
a11から始まる例
a00から始まる例
|