Cazoo! Japan CLAPACK DGETRF/ZGETRF

DGETRF/ZGETRF (CLAPACK)

DGETRF/ZGETRF

m行n列の一般行列[A]を、LU分解する。

LU分解後行列[A]は [A]=[P]*[L]*[U]の形体をとっている。

[U]に上三角行列を、[L]に下三角行列を、[P]に置換行列を表す。

ただし、m<nのとき[U]は上台形行列、m>nのとき[P]は下台形行列を表す。

＜関数＞

int dgetrf_( integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info);

＜引数＞

long int m - (input) 行列[A]の列数。(m≧0)

long int n - (input) 行列[A]の行数。(n≧0)

double a - (input/output) 配列形式はa[lda×n]。

(input)LU分解をするm行n列の行列[A]。

(output)LU分解後の行列因子[L]と[U]。ただし、行列[L]の単位対角要素は格納されない。

A=[P]*[L]*[U]であり、[P]は行と列を入れ替える操作に対応し、置換行列と呼ばれ、0か1が格納される。

long int lda - (input) 行列Aの第一次元(のメモリ格納数)。lda≧max(1,m) つまり、1～mの中で一番大きい数よりも、ldaが大きくなければならない。

通常はlda=mで良い。

long int ipiv - (output) 次元数min(m,n)の配列。軸選択用添え字。つまり、ipivはmとnのうちで小さい方の数だけメモリを確保する必要がある。

1 <= i <= min(m,n)なる行列Aのi行は、IPIV[i]と交換された。

long int info - (output)

info = 0: 正常終了

info < 0: info = -i ならば、i番目の引数の値が間違えていることを示す。

info > 0: info = -i ならば、U[i][i]は0である。LU分解は行われたが、上三角因子Uが特異になっており、

この結果を利用して方程式を解くとゼロ割が発生する。

を解く。

#include <stdio.h>

#define M 4

#define N 3

double A[M*N];

int main(void)

{

static int i;

static long int m=M,n=N,lda=M,info,piv[N];

A[0]=1.;A[1]=3.;A[2]=1.; A[3]=1.;

A[4]=1.;A[5]=1.;A[6]=2.;A[7]=4.;

A[8]=1.;A[9]=3.;A[10]=-5.;A[11]=3.;

dgetrf_( &m, &n, A, &lda, piv, &info);

for(i=0;i<(N*M);++i) printf("%lf\n",A[i]);

return 0;

}