LU分解について

こんにちは。
行列のＬＵ分解のｃによるプログラムを勉強している者です。
ＬＵ分解のｃによるプログラムについてお尋ねしたいことがございます。


ニューメリカルレシピに、ＬＵ分解の方法として
以下のようなコードが記載せれておりました。


void ludcmp(double **a, int n, int *indx, double *d)
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;// 各行の暗黙のスケーリングを記録する．

vv=vector(n);
*d=1.0;// まだ行交換していない．
for (i=0;i<n;i++) {// 行についてループし，暗黙のスケーリングの情報を得る．
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
if (big == 0.0) printf("Singular matrix in routine ludcmp\n");// 最大要素が０なら特異行列である．
vv[i]=1.0/big;// スケーリングを記録する．
}
for (j=0;j<n;j++) {// Crout法，列についてのループ
for (i=0;i<j;i++) {// 方程式(2.3.12)のi=j以外
sum=a[i][j];
for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++) {
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big) {
big=dum;
imax=i;
}
}
if (j != imax) {
for (k=0;k<n;k++) {
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d = -(*d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != n) {
dum=1.0/(a[j][j]);
for (i=j+1;i<n;i++) a[i][j] *= dum;
}
}
free_vector(vv);
}


例えば、ＬＵ分解させたい行列３×３の（｛4,7,6｝｛2,5,9｝｛3,1,8｝）だとしたら、
メイン関数にどのように書いたらＬＵ分解を実行できるでしょうか？

ニューメリカルレシピを買って読んでみたものの、よく理解できず、質問させていただきます。
稚拙な質問かもしれませんが、どうぞよろしくお願いいたします。

No.2 ベストアンサー 回答者： AsanoNagi 回答日時： 2011/06/07 22:47

あ、もしかして、

double a[3][3];
じゃなくて、

double a1[] = {4, 7, 6};
double a2[] = {2, 5, 9};
doulle a3[] = {3, 1, 8};
double *a[] = {a1, a2, a3};

なら、とりあえずはうまくいきますね。

でも、まだ、ludcmp(); をじっくり眺めるとおかしなところが残っています。
特異行列の時の処理とか。
あと、vector() はたぶん、どこかで独自に定義しているでしょうね。

この回答への補足

丁寧にお答えくださいましてありがとうございます。
さっそくお答えいただいた通りに、以下のようにして動かしてみたのですが、動きませんでした・・・。
どうしてでしょうか？？（質問ばかりで申し訳ありません）

#include <stdio.h>
#include <math.h>

#include "nr.h"


#define TINY 1.0e-20;


int main()
{
double d;
int indx;
double a1[] = {4, 7, 6};
double a2[] = {2, 5, 9};
double a3[] = {3, 1, 8};
double *a[] = {a1, a2, a3};


ludcmp(a, 3, &indx, &d);


return(0);

}

void ludcmp(double **a, int n, int *indx, double *d)
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;// 各行の暗黙のスケーリングを記録する．

vv=vector(n);
*d=1.0;// まだ行交換していない．
for (i=0;i<n;i++) {// 行についてループし，暗黙のスケーリングの情報を得る．
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
if (big == 0.0) printf("Singular matrix in routine ludcmp\n");// 最大要素が０なら特異行列である．
vv[i]=1.0/big;// スケーリングを記録する．
}
for (j=0;j<n;j++) {// Crout法，列についてのループ
for (i=0;i<j;i++) {// 方程式(2.3.12)のi=j以外
sum=a[i][j];
for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++) {
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big) {
big=dum;
imax=i;
}
}
if (j != imax) {
for (k=0;k<n;k++) {
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d = -(*d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != n) {
dum=1.0/(a[j][j]);
for (i=j+1;i<n;i++) a[i][j] *= dum;
}
}
free_vector(vv);
}


ちなみに

>あと、vector() はたぶん、どこかで独自に定義しているでしょうね。

その通りです。ヘッダファイルで定義されております。
記載せず、申し訳ありませんでした。

/*nr.h*/

#define DIM3

void ludcmp(double **a, int n, int *indx, double *d);
double *vector(long size);
void free_vector(double *v);
double **matrix(long rows, long cols);
void free_matrix(double **m);
int *ivector(long nh);
void free_ivector(int *v);

補足日時：2011/06/08 11:45

この回答へのお礼

お礼欄に書くのはおかしいですが、

でたエラーメッセージは

1>------ ビルド開始: プロジェクト: LU, 構成: Debug Win32 ------
1> LU.cpp
1>LU.obj : error LNK2019: 未解決の外部シンボル "void __cdecl free_vector(double *)" (?free_vector@@YAXPAN@Z) が関数 "void __cdecl ludcmp(double * *,int,int *,double *)" (?ludcmp@@YAXPAPANHPAHPAN@Z) で参照されました。
1>LU.obj : error LNK2019: 未解決の外部シンボル "double * __cdecl vector(long)" (?vector@@YAPANJ@Z) が関数 "void __cdecl ludcmp(double * *,int,int *,double *)" (?ludcmp@@YAXPAPANHPAHPAN@Z) で参照されました。
1>C:\Users\njdata\Documents\Visual Studio 2010\Projects\LU\Debug\LU.exe : fatal error LNK1120: 外部参照 2 が未解決です。
========== ビルド: 0 正常終了、1 失敗、0 更新不要、0 スキップ ==========

でした。お伝えし忘れていました、すみません。

お礼日時：2011/06/08 12:05

No.3 回答者： AsanoNagi 回答日時： 2011/06/08 13:01

単に、vector の実体がリンクの指定にないだけです。

void ludcmp(double **a, int n, int *indx, double *d)
と同じように、vector の中身が書いてあるファイルをリンクしてください。

この回答へのお礼

ご丁寧にありがとうございます。
どうしてもうまくいかないのですが、
方向性は示して頂けたので、もう少しじっくり自分で勉強してみたいと思います。
何度もご丁寧にありがとうございました。大変参考になりました！

お礼日時：2011/06/10 10:28

No.1 回答者： AsanoNagi 回答日時： 2011/06/07 22:34

気持ちとしては、こんな感じなんだろうなと思います。

でも、これでは、うまくいきません。
肝心の ludcmp(); 自体にいくつかの間違いがあるからです。
勉強中なのだし、まずは、ludcmp(); の間違いを修正するというのはいいかもしれません。

int main()
{
double d;
int indx;
double a[3][3] =
{
{4, 7, 6},
{2, 5, 9},
{3, 1, 8}
};

ludcmp(a, 3, &indx, &d);
// その他の処理
}