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｝）だとしたら、
メイン関数にどのように書いたらＬＵ分解を実行できるでしょうか？
ニューメリカルレシピを買って読んでみたものの、よく理解できず、質問させていただきます。
稚拙な質問かもしれませんが、どうぞよろしくお願いいたします。
補足
こういう方の解答は求めてません。 面白と思いですか？(笑)