15

LAPACK に特定の部分行列の要素を与える関数はありますか? もしそうなら、C ++の構文はどうですか?

それともコード化する必要がありますか？

4

1 に答える 1

26

サブマトリックスにアクセスするための関数はありません。ただし、マトリックスデータがLAPACKルーチンに格納される方法のため、マトリックスデータは必要ありません。これにより多くのコピーが節約され、データレイアウトは（部分的に）次の理由で選択されました。

LAPACKの密な（つまり、縞模様ではない、三角形、エルミートなど）行列は、次の4つの値で定義されることを思い出してください。

  • マトリックスの左上隅へのポインタ
  • マトリックスの行数
  • 行列の列数
  • マトリックスの「主要な次元」。通常、これは行の隣接する要素間のメモリ内の距離です。

ほとんどの場合、ほとんどの人は行数に等しい先頭のディメンションのみを使用します。3x3行列は通常、次のように格納されます。

a[0] a[3] a[6] 
a[1] a[4] a[7]
a[2] a[5] a[8]

代わりに、先頭の次元を持つ巨大な行列の3x3の部分行列が必要だとしldaます。左上隅がa（15,42）にある3x3部分行列が特に必要だとします。

         .             .            .
         .             .            .
... a[15+42*lda] a[15+43*lda] a[15+44*lda] ...
... a[16+42*lda] a[16+43*lda] a[16+44*lda] ...
... a[17+42*lda] a[17+43*lda] a[17+44*lda] ...
         .             .            .
         .             .            .

この3x3マトリックスを連続したストレージにコピーすることもできますが、入力（または出力）マトリックスとしてLAPACKルーチンに渡す場合は、その必要はありません。パラメータを適切に定義するだけで済みます。この部分行列と呼びましょうb; 次に、次のように定義します。

// pointer to the top-left corner of b:
float *b = &a[15 + 42*lda];
// number of rows in b:
const int nb = 3;
// number of columns in b:
const int mb = 3;
// leading dimension of b:
const int ldb = lda;

驚くべきことは、ldb;の値だけです。lda「大きな行列」の値を使用することで、コピーせずに部分行列をアドレス指定し、その場で操作できます。

しかし 私は嘘をついた（一種の）。場合によっては、サブマトリックスを適切に操作できず、本当にコピーする必要があります。それはめったにないので、話したくありませんでした。可能な限りインプレース操作を使用する必要がありますが、それが可能であると言わないのは残念です。ルーチン：

SLACPY(UPLO,M,N,A,LDA,B,LDB)

左上隅が先行次元で格納されているx行列を、左上隅が先行次元であるMx行列にコピーします。このパラメーターは、上三角、下三角、または行列全体のいずれをコピーするかを示します。NALDAMNBLDBUPLO

上記の例では、次のように使用します（clapackバインディングを想定）。

...
const int m = 3;
const int n = 3;
float b[9];
const int ldb = 3;
slacpy("A",  // anything except "U" or "L" means "copy everything"
       &m,   // number of rows to copy
       &n,   // number of columns to copy
       &a[15 + 42*lda], // pointer to top-left element to copy
       lda,  // leading dimension of a (something huge)
       b,    // pointer to top-left element of destination
       ldb); // leading dimension of b (== m, so storage is dense)
...
于 2011-02-17T16:43:26.197 に答える