基本は2重対角化と同じです。
と、その前に。
前回 LAPACK の関数を紹介しましたが、LAPACK や BLAS の関数には、LDA だとか INCX だとか良く分からない引数が出てきます。
LDA については、マニュアルなどでは Leading Dimension の事で、行列の行数を M として
LDA >= MAX (1, M)
などとなっています。
LAPACK や BLAS では、本来2次元であるべき行列も 1次元の連続したメモリの中に格納されます。
因みにその格納方法は col major で、例えば
という行列は
double data[] = { 1、 4、 7、 2、 5、 8、 3、 6、 9 }
と、列方向の順に行列の要素を格納しなければいけません。
で、上の様な1次元配列だとこれがどんなディメンジョンの行列なのか分からなくなるため、1列にいくつの要素が含まれるかを LDA で指定しなければなりません。
上の行列の場合だと LDA = 3 ですね。
<- LDA -> <- LDA -> <- LDA ->
double data[] = { 1、 4、 7、| 2、 5、 8、| 3、 6、 9 }
でも大抵の場合、LAPACK では LDA とは別に行数 M も指定します。
上の行列だと M = 3 で LDA と同じです。
上の行列なら LDA も M も、どちらも行数の 3 になるんですが、なんでこんな風におんなじパラメータを名前を変えて指定しなければいけないんでしょう？？
でも実は LDA は結構便利なんです。
LDA が重要な働きをするのが部分行列です。
例えば行列 A の 2 行 2 列目から 2 行 2 列部分の部分行列を取りだすと
となります。
で、B の先頭要素は data[ ] の何番目かというと・・・
数えてみれば B の先頭要素 （5） は double data[ ] の 5 番目になります。
<- LDA -> <- LDA -> <- LDA ->
double data[] = { 1、 4、 7、| 2、 5、 8、| 3、 6、 9 }
この 5 番目、というのをもう少し詳しく見てみると・・・
B の1 行 1 列成分は、A の2 行 2 列になります。なので col major で順番を数えると、
A の 1列分の 3 個の要素を飛ばして、さらに 2 行目なので 2 個目、
これで data[ ] の 5番目、となる訳ですね。
一般に、行数 MA である行列 A の i 行 j 列目の要素のポインタは、 A の先頭要素のポインタを *data として
data + i + MA * j
となります。
この様に、元の行列 A とその i 行 j 列目から n1 行 n2 列部分を抜き出した部分行列 B では、両者の行数が異なるので
LDA = MA (A の行数)
M = n1 (B の行数)
N = n2 (B の列数)
先頭要素のポインタ = data + i + j * LDA
とすれば部分行列 B が指定出来ます。
GSL の Householder 変換の関数では、入力行列を変換後の行列に逐次書き換えて返します。
LAPACK で相当するのは
void
dlarfx_ (char *side, int *m, int *n, double *v, double *tau, double *a, int *lda, double *work)
ですが、これも入力行列 double *a を逐次書き換えます。
こんな場合、上の様に部分行列を指定すれば、部分行列部分をコピーして変換し元の行列にコピーし直して・・・などという手間もメモリも節約できます。
さらに BLAS では、例えば（CBLASですが・・・）
void cblas_dscal (int n, double alpha, double *x, int incX)
などと incX なるものが出てきます。
この関数は、ベクトル（とは言っても単なる1次元配列） *x を、 incX 置きに alpha 倍したもので置き換える、つまり
for (i = 0; i < n; i++) x[ i * incX ] *= alpha;
に相当する操作をします。
この incX もうまく使えば結構便利です。
例えば行数 M、 Leading Dimension が LDA である行列 double *a の i 行を alpha 倍したければ
cblas_dscal (M, alpha, a + i, LDA)
で出来ます。
GSL では row major なので上に書いた事と少し違ってきて混乱しますが・・・
LDA に当たる物が、 gsl_matrix *a のメンバ
size_t a->tda;
incX に当たる物が gsl_vector *v のメンバ
size_t v->stride;
になります。
そもそも、gsl_matrix、 gsl_vector とは（1部省略してますが）
typedef struct {
size_t size1; // 行数
size_t size2; // 列数
size_t tda; // LDAに相当。但しrow major なので列数
double *data; // 行列先頭要素のポインタ
int owner; // オーナー情報
} gsl_matrix;
typedef struct {
size_t size; // 要素数
size_t stride; // incX に相当
double *data; // ベクトル先頭要素のポインタ
int owner; // オーナー情報
} gsl_vector;
と定義されています。
で、例えば、行列 gsl_matrix *a の i 行 j 列目から n1 行 n2 列部分の部分行列（の像）を求める関数
gsl_matrix_view m = gsl_matrix_submatrix (a, i, j, n1, n2);
では
&m.matrix.size1 = n1;
&m.matrix.size2 = n2;
&m.matrix.tda = a->tda;
&m.matrix.data = a->data + j + i * a->tda; // row major なので
&m.matrix.owner = false;
と、また行列 gsl_matrix *a の j 列ベクトル（の像）を求める関数
gsl_vector_view c = gsl_matrix_column ( a, j );
では
&c.vector.size = a->size1;
&c.vector.stride = a->tda; // row major なので 列ベクトルの時にこう指定される
&c.vector.data = a->data + j; // row major なので （しつこい・・・）
&c.vector.owner = false;
となります。
owner は data が自前のものか余所を参照しているかを区別する者で、これが false の場合には行列、ベクトルを gsl_matrix_free や gsl_vector_free しても参照先の data を解放しません。
おっと、3重対角化でしたっけ・・・
こちらをお読みください（汗