はじめに
scalapackとはマルチプロセス環境で線形代数計算を実行するためのライブラリである。lapackの並列版といえるオープンソースのライブラリで、多くのスパコンにプリインストールされている。
この記事ではscalapackを導入し、線型方程式を解くところまでを解説する。
行列を複数のプロセスで分散して持つことになるが、どのように分散するかということも含めて解説していく。また、lapackと同様にscalapackはFortranで実装されているため使い方にクセがあり、特にC,C++から使う時はネット上に情報が不足している。
この記事ではC++からscalapackを使うための手順を示す。
まずはscalapackをインストールしてサンプルコードを動かすまでの手順を示し、その後サンプルコードの解説を行っていく。
インストール
ここではmac上でscalapackを利用するための手順を示す。homebrewで簡単にインストールすることができる。
$ brew install open-mpi # MPI環境をインストールしていない場合はインストールする
$ brew install scalapack
サンプルコード
以下のサンプルコードを準備する。
ちなみにこちらのサンプルコードでは変数を出力するのに便利なように icecream-cppというライブラリを使用している。
こちらの記事で、このライブラリの使用方法を紹介しているので参考にして欲しい。
# include <cstdlib>
# include <iostream>
# include <mpi.h>
# include <random>
# include <vector>
# include "icecream-cpp/icecream.hpp"
// Declaration of Fortran subroutines
extern "C" {
void sl_init_(int *icontext, int *nprow, int *npcolumn);
// SL_INIT initializes an NPROW x NPCOL process grid using a row-major ordering
// of the processes. This routine retrieves a default system context which will
// include all available processes. (out) ictxt, (in) nprow, npcolumn
void blacs_gridinfo_(int *icontext, int *nprow, int *npcolumn, int *myrow,
int *mycolumn);
// (in) icontext: BLACS context
// (out) nprow, npcolumn: the numbers of rows and columns in this process grid
// (out) myrow, mycolumn: the process grid row- and column-index
void blacs_exit_(int *cont);
// (in) continue: if 0, all the resources are released. If nonzero, MPI
// resources are not released.
void blacs_gridexit_(int *icontext);
// (in) icontext: BLACS context
void descinit_(int *desc, int *m, int *n, int *mb, int *nb, int *irsrc,
int *icsrc, int *icontext, int *lld, int *info);
// (out) descriptor for the global matrix. `desc` must be an array of int of
// length 9. int[9]
// (in) m, n: rows and columns of the matrix (in) mb, nb: row,
// column block sizes
// (in) irsrc, icsrc: the process row (column) over which the
// first row of the global matrix is distributed.
// (in) icontext: BLACS context
// (in) lld: leading dimension of the local array
// (out) info: 0 => completed successfully
void dgesd2d_(int *icontext, int *m, int *n, double *A, int *lda, int *r_dest,
int *c_dest);
// Takes a general rectangular matrix and sends it to the destination process.
// (in) icontext: BLACS context
// (in) m,n: matrix sizes
// (in) A: matrix
// (in) lda: leading dimension (m)
// (in) r_dest, c_dest: the process corrdinate of the process to send the
// message to
void dgerv2d_(int *icontext, int *m, int *n, double *A, int *lda, int *row_src,
int *col_src);
// Receives a message from the process into the general rectangular matrix.
// (in) icontext: BLACS context
// (in) m,n,lda: sizes of the matrix
// (out) A: matrix
// (in) row_src, col_src: the process coordinate of the source of the message
void pdgesv_(int *n, int *nrhs, double *A, int *ia, int *ja, int desc_a[9],
int *ipvt, double *B, int *ib, int *jb, int desc_b[9], int *info);
// These subroutines solve the following systems of equations for multiple
// right-hand sides: AX = B
// (in) n: order of the submatrix = the number of rows of B
// (in) nrhs: the number of columns of B
// (in/out) A: the local part of the global general matrix A.
// (in) ia, ja: the row and the column indices of the
// global matrix A, identifying the first row and column of the submatrix A.
// (in) desc_a: descriptor of A matrix
// (out) ipvt: the local part of the global vector ipvt, containing the pivot
// indices.
// (in/out) B: the local part of the global general matrix B,
// containing the right-hand sides of the system.
// (in) ib, jb: the row and the column indices of the global matrix B,
// identifying the first row and column of the submatrix B.
// (in) desc_b: descriptor of B matrix (out) info: error code
}
int main(int argc, char *argv[]) {
icecream::ic.disable();
// Get command line arguments
if (argc != 5) {
exit(1);
}
// Number of columns in the global matrix A
int N = std::strtol(argv[1], nullptr, 10);
int M = N; // Number of rows in the global matrix A
// block size for the matrix A
int NB = std::strtol(argv[2], nullptr, 10);
int MB = NB;
// Number of rows and columns in the process grid
int NPROW = std::strtol(argv[3], nullptr, 10);
int NPCOL = std::strtol(argv[4], nullptr, 10);
// Initialize the process grid
int ICTXT, MYROW, MYCOL;
sl_init_(&ICTXT, &NPROW, &NPCOL);
blacs_gridinfo_(&ICTXT, &NPROW, &NPCOL, &MYROW, &MYCOL);
int my_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
icecream::ic.prefix(
"ic (", [my_rank] { return my_rank; }, "): ");
IC(ICTXT, NPROW, NPCOL, MYROW, MYCOL);
// If I'm not in the process grid, go to the end of the program
if (MYROW == -1) {
int blacs_exitcode = 0;
blacs_exit_(&blacs_exitcode);
exit(0);
}
// global matrix A and B
std::vector<std::vector<double>> A;
std::vector<double> B, X;
// Generate matrices A and B
if (MYROW == 0 && MYCOL == 0) {
A.resize(M, std::vector<double>(N, 0.0));
B.resize(M);
X.resize(M, 0.0);
// Do not ever do like this in real world. Use high performance random array
// generators
std::mt19937_64 rnd(123456789ull);
std::uniform_real_distribution<double> uni;
std::cout << "Matrix A:\n";
// A is a row-major matrix (as in C,C++) but SUB_A is
// the column major matrix as in Fortran
for (int i = 0; i < M; i++) {
for (int j = 0; j < N; j++) {
// ROW-MAJOR matrix
A[i][j] = uni(rnd);
std::cout << A[i][j] << ' ';
}
std::cout << std::endl;
}
std::cout << "Matrix B:\n";
for (int i = 0; i < M; i++) {
B[i] = uni(rnd);
std::cout << B[i] << ' ';
}
std::cout << std::endl;
}
//*****
// Distribute matrices to the process grid
//*****
int CSRC = 0, RSRC = 0; // root process coordinate
int NRHS = 1; // Number of columns in the global solution matrix B
// size of SUB_A
int SUB_A_COLS = (N / (NB * NPCOL)) * NB + std::min(N % (NB * NPCOL), NB);
int SUB_A_ROWS = (M / (MB * NPROW)) * MB + std::min(M % (MB * NPROW), MB);
// size of SUB_B
int SUB_B_ROWS = SUB_A_ROWS;
int SUB_B_COLS = 1;
IC(RSRC, CSRC, SUB_A_COLS, SUB_A_ROWS, SUB_A_ROWS, SUB_B_ROWS);
// Initialize array descriptors
int DESCA[9], DESCB[9];
int INFO;
descinit_(DESCA, &M, &N, &MB, &NB, &RSRC, &CSRC, &ICTXT, &SUB_A_ROWS, &INFO);
descinit_(DESCB, &N, &NRHS, &MB, &NRHS, &RSRC, &CSRC, &ICTXT, &SUB_B_ROWS,
&INFO);
std::vector<int> IPIV(SUB_A_ROWS + NB);
// COLUMN-MAJOR matrices
std::vector<double> SUB_A(SUB_A_ROWS * SUB_A_COLS, 0.0);
std::vector<double> SUB_B(SUB_B_ROWS * SUB_B_COLS, 0.0);
// convert submatrix index (i,j) at process (p_row, p_col) into global
// coordinate (I,J)
auto to_global_coordinate = [SUB_A_ROWS, SUB_A_COLS, MB, NB, NPROW,
NPCOL](int i, int j, int p_row,
int p_col) -> std::pair<int, int> {
// block coordinate (bi, bj)
int bi = i / MB;
int bj = j / NB;
// local coordinate inside the block
int ii = i % MB;
int jj = j % NB;
// calculate global coordinate
int I = bi * (MB * NPROW) + p_row * MB + ii;
int J = bj * (NB * NPCOL) + p_col * NB + jj;
return std::make_pair(I, J);
// set value to local matrix
};
auto copy_A_to_SUB_A = [&A, &SUB_A, SUB_A_ROWS, SUB_A_COLS, M, N,
&to_global_coordinate](int p_row, int p_col) {
for (int i = 0; i < SUB_A_ROWS; i++) {
for (int j = 0; j < SUB_A_COLS; j++) {
auto IJ = to_global_coordinate(i, j, p_row, p_col);
// set value to local matrix. Note SUB_A is column major
int I = IJ.first, J = IJ.second;
if (I >= M || J >= N)
continue;
// IC(p_row, p_col, i, j, I, J);
SUB_A[i + j * SUB_A_ROWS] = A[I][J];
}
}
};
auto copy_B_to_SUB_B = [&B, &SUB_B, SUB_B_ROWS, SUB_B_COLS, M, NRHS,
&to_global_coordinate](int p_row, int p_col) {
for (int i = 0; i < SUB_B_ROWS; i++) {
for (int j = 0; j < SUB_B_COLS; j++) {
auto IJ = to_global_coordinate(i, j, p_row, p_col);
// set value to local matrix. Note SUB_A is column major
int I = IJ.first, J = IJ.second;
if (I >= M || J >= NRHS)
continue;
SUB_B[i + j * SUB_B_ROWS] = B[I * NRHS + J];
}
}
};
// distribute SUB_A
int LDA = 1; // The distance between two elements in matrix row
if (MYROW == 0 && MYCOL == 0) {
for (int p_row = 0; p_row < NPROW; p_row++) {
for (int p_col = 0; p_col < NPCOL; p_col++) {
if (p_row == 0 && p_col == 0)
continue;
copy_A_to_SUB_A(p_row, p_col);
dgesd2d_(&ICTXT, &SUB_A_ROWS, &SUB_A_COLS, SUB_A.data(), &LDA, &p_row,
&p_col);
}
}
copy_A_to_SUB_A(0, 0); // set local submatrix for rank 0
} else {
dgerv2d_(&ICTXT, &SUB_A_ROWS, &SUB_A_COLS, SUB_A.data(), &LDA, &RSRC,
&CSRC);
}
// distribute SUB_B
if (MYROW == 0 && MYCOL == 0) {
for (int p_row = 1; p_row < NPROW; p_row++) {
int p_col = 0;
copy_B_to_SUB_B(p_row, p_col);
dgesd2d_(&ICTXT, &SUB_B_ROWS, &SUB_B_COLS, SUB_B.data(), &LDA, &p_row,
&p_col);
}
copy_B_to_SUB_B(0, 0);
} else if (MYCOL == 0) {
dgerv2d_(&ICTXT, &SUB_B_ROWS, &SUB_B_COLS, SUB_B.data(), &LDA, &RSRC,
&CSRC);
}
if (my_rank == 1) {
IC(SUB_A, SUB_B);
}
//*****
// Call the SCALAPACK routine
//*****
int IA = 1;
int JA = 1;
int IB = 1;
int JB = 1;
pdgesv_(&N, &NRHS, SUB_A.data(), &IA, &JA, DESCA, IPIV.data(), SUB_B.data(),
&IB, &JB, DESCB, &INFO);
// collect SUB_B
auto copy_SUB_B_to_X = [&X, &SUB_B, SUB_B_ROWS, SUB_B_COLS, M, NRHS,
&to_global_coordinate](int p_row, int p_col) {
for (int i = 0; i < SUB_B_ROWS; i++) {
for (int j = 0; j < SUB_B_COLS; j++) {
auto IJ = to_global_coordinate(i, j, p_row, p_col);
// set value to local matrix. Note SUB_A is column major
int I = IJ.first, J = IJ.second;
if (I >= M || J >= NRHS)
continue;
X[I] = SUB_B[i + j * SUB_B_ROWS];
}
}
};
if (MYROW == 0 && MYCOL == 0) {
copy_SUB_B_to_X(0, 0);
for (int p_row = 1; p_row < NPROW; p_row++) {
int p_col = 0;
dgerv2d_(&ICTXT, &SUB_B_ROWS, &SUB_B_COLS, SUB_B.data(), &LDA, &p_row,
&p_col);
copy_SUB_B_to_X(p_row, p_col);
}
} else if (MYCOL == 0) {
dgesd2d_(&ICTXT, &SUB_B_ROWS, &SUB_B_COLS, SUB_B.data(), &LDA, &RSRC,
&CSRC);
}
if ((MYROW == 0) && (MYCOL == 0)) {
std::cout << "Matrix X:\n";
for (int i = 0; i < M; i++) {
std::cout << X[i] << ' ';
}
std::cout << std::endl;
std::cout << "Matrix AX:\n";
for (int i = 0; i < M; i++) {
double ans = 0.0;
for (int j = 0; j < N; j++) {
ans += A[i][j] * X[j];
}
std::cout << ans << ' ';
}
std::cout << std::endl;
}
// Release the process grid
blacs_gridexit_(&ICTXT);
int blacs_exitcode = 0;
blacs_exit_(&blacs_exitcode);
return 0;
}
このファイルを"pdgesv_example.cpp"というファイルに保存する。
同じディレクトリにicecream-cppを次の手順で導入する。ヘッダオンリーのライブラリなのでコードをクローンしてくるだけで利用できる。
$ git clone git@github.com:renatoGarcia/icecream-cpp.git
以上の準備が整ったら、ビルドと実行は次のように行う。難しいことはなく、scalapackのライブラリをリンクしてビルドしたら、あとは普通にmpiexec
するだけである。
$ mpicxx -std=c++11 pdgesv_example.cpp -lscalapack
$ mpirun -np 4 ./a.out 8 3 2 2
(実行結果)
Matrix A:
0.348872 0.266886 0.136646 0.0285569 0.868933 0.350855 0.0711051 0.323368
0.555103 0.875991 0.204709 0.892759 0.584466 0.369779 0.850631 0.391382
0.119661 0.754243 0.695023 0.686615 0.931935 0.454888 0.0674011 0.337989
0.974885 0.726438 0.0454151 0.745967 0.496126 0.716716 0.859742 0.134076
0.488442 0.871219 0.766468 0.251256 0.166365 0.743796 0.980511 0.729577
0.901105 0.264365 0.885651 0.882112 0.748933 0.919626 0.693453 0.215403
0.828589 0.0442154 0.863038 0.352605 0.77204 0.58612 0.322777 0.172931
0.805364 0.306002 0.2191 0.724731 0.696487 0.911934 0.679563 0.354942
Matrix B:
0.73897 0.187402 0.314613 0.137569 0.653774 0.270132 0.899839 0.573423
Matrix X:
2.8599 -0.0194602 -0.0624903 0.561362 -1.11953 -0.847793 -2.77944 3.73202
Matrix AX:
0.73897 0.187402 0.314613 0.137569 0.653774 0.270132 0.899839 0.573423
$Ax = B$を解いて、実際に$Ax$の値を計算して出力している。ここで$A$は8x8の行列、$B$は8x1の行列(列ベクトル)で、各要素は$[0,1]$の一様分布からランダムにサンプリングしている。
解いて得られた$X$を用いて、$Ax$の値も計算し、実際に$B$と同じ結果が得られていることを確かめている。
行列の分割方法
scalapackでは行列を複数のプロセス(別々のメモリー空間)に分散して保持する。
その行列の持ち方をここでは解説する。
まず理解したいのは、MPIプロセスが2次元の格子状に配置されるということである。
rank0のプロセスから順番に(NPROW x NPCOL)
の格子状に配置される。このプロセスが並んだ格子を「プロセスグリッド」と呼ぶ。
上記のサンプルコードでは実行時の引数としてプロセッサーグリッドのサイズを指定するようにしている。
例えば、NPROW=2, NPCOL=3 の場合、各ランクのプロセスは次のように配置される。
(上の図では今後わかりやすく説明するためプロセッサごとに異なる色をつけた。)
続いて行列の分割について解説する。scalapackでは"ブロックサイクリック分割"という方法でデータを分散して保持する。
まず、$M \times N$の行列$A$があり、それらをブロックサイズ$MB \times NB$のブロックで分割する。
例えば$M=N=9$, $MB=NB=2$の場合、こんな感じで分割される。
これらのブロックに対してプロセスグリッドのプロセッサを割り当てていく。
各色が、どのプロセスがどのブロックを担当しているかを表しており、それぞれのプロセスが保持している局所的な行列SUB_A
はそれぞれの色のブロックを詰め合わせることによって得られる。
たとえば青色のSUB_Aは、以下のように作られ、5x5の行列になる。このように行列のとびとびの部分をそれぞれのプロセスで保持することになる。
以上をふまえると、プロセスグリッドの座標(p_row,p_col)のプロセスにおけるSUB_Aの(i,j)要素が、グローバルな行列Aのどの座標(I,J)と対応するかを計算するコードは次のように書ける。
auto to_global_coordinate = [SUB_A_ROWS, SUB_A_COLS, MB, NB, NPROW, NPCOL]
(int i, int j, int p_row, int p_col) -> std::pair<int, int> {
// block coordinate (bi, bj)
int bi = i / MB;
int bj = j / NB;
// local coordinate inside the block
int ii = i % MB;
int jj = j % NB;
// calculate global coordinate
int I = bi * (MB * NPROW) + p_row * MB + ii;
int J = bj * (NB * NPCOL) + p_col * NB + jj;
return std::make_pair(I, J);
};
SUB_Aをメモリ上で保持する方法についても注意が必要である。scalapackはFortranで実装されているので、column-majorで要素を持つ必要がある。
具体的にはSUB_Aは5x5=25の長さのdoubleの1次元配列として持ち、先頭から(0,0)->(1,0)->(2,0)-> ... (3,4)->(4,4)
の順番で格納する必要がある。
つまり、以下のように保存する。
SUB_A[i + j * SUB_A_ROW] = SUB_Aの(i,j)要素
上の図を見てもらうと気づくと思うが、SUB_Aの大きさはプロセスごとに異なる。しかし、どのプロセスでもSUB_A
は同じ大きさの配列として確保されなくてはならないので注意が必要である。
例えば上図の黄色のプロセスのSUB_Aの行数は4であるが、5x5の領域を確保し、さらに5行目に相当する箇所は空けておく必要がある。つまり(4,0)に相当する箇所は空白にしておき、(0,1)に対応する部分は配列の6番目(SUB_A[5]
)に格納する。(実際にはSUB_A[4]
は参照されないので任意の数字を入れておいても問題なく動作する)
コードの解説
では上のコードがどのように動作しているのかより詳細に解説していく。
scalapack関数の宣言
まず、前提としてscalapackの関数をCやC++から実行しようとするときに、scalapackの関数の宣言が必要になる。
一般にC,C++ライブラリではこのような宣言を含んだヘッダファイルが提供されるのが一般的だが、scalapackの場合そのようなヘッダファイルは(私が調べた限り)存在しない。
よって、自分で宣言を書く必要がある。"scalapack pdgesv" などのように関数名で検索するとその関数のドキュメントがヒットすると思うので、そのドキュメントを参考に宣言を書いていく。
サンプルコードでは
extern "C" {
void sl_init_(int *icontext, int *nprow, int *npcolumn);
....
となっている箇所が該当する部分である。Fortranの関数を呼ぶことになるので、extern "C"
としてCリンケージを指定する必要がある。
どのように関数宣言を準備すれば良いかFortranを理解できていないと自明ではないのでややこしいが、コツコツとやるしかない。しかも、引数のどれが入力でどれが出力なのかわかりづらいのもツラい。(C言語のヘッダファイルを用意してくれればいいのに...)
初期化、終了処理
scalapackを使った処理を記述するときに、プログラムの最初と最後におまじない的に書く処理がある。
初期化処理
// Initialize the process grid
int ICTXT, MYROW, MYCOL;
sl_init_(&ICTXT, &NPROW, &NPCOL);
blacs_gridinfo_(&ICTXT, &NPROW, &NPCOL, &MYROW, &MYCOL);
終端処理
// Release the process grid
blacs_gridexit_(&ICTXT);
int blacs_exitcode = 0;
blacs_exit_(&blacs_exitcode);
- 初期化処理の
sl_init_
は初期化処理を行、ICTXT
を返す。ICTXTはMPIのcommunicatorのようなもの。通常、ICTXT=0となるはず。入力としてプロセスグリッドのサイズNPROW
,NPCOL
を与える。-
sl_init_
の中ではMPI_Init
も呼び出されるらしく、これとは別にMPI_Init
を呼ぶ必要はない。
-
-
blacs_gridinfo
はプロセスグリッド上の自分のプロセスの座標(MYROW
,MYCOL
) を取得する関数。- MPIのランクを
my_rank
とすると、MYROW == my_rank / NPCOL
,MYCOL == my_rank % NPCOL
となっているはず。
- MPIのランクを
- プログラムの終了時には
blacs_gridexit_
とblacs_exit_
を呼ぶ必要がある。- この中で
MPI_Finalize
も呼ばれているようだ。
- この中で
行列A,Bの準備
行列A,Bについて、このサンプルコードではランク0のプロセスで構築し、サブ行列を各プロセスに転送するという処理を行っている。
これはサンプルコードのわかりやすさと、解の検証のためにやっているだけで、より大きな行列の処理を行うときには各プロセスで分散して構築する必要がある。
このような処理を行っている箇所はサンプルコードの以下の箇所である。
ランク0ではAからプロセス(p_row,p_col)用のSUB_Aを構築し、該当プロセスに転送している。データ転送にはdgesd2d_
, dgerv2d_
というルーチンを利用しているが、単純にMPI_Send
, MPI_Recv
を使っても良いはず。
// distribute SUB_A
int LDA = 1; // The distance between two elements in matrix row
if (MYROW == 0 && MYCOL == 0) {
for (int p_row = 0; p_row < NPROW; p_row++) {
for (int p_col = 0; p_col < NPCOL; p_col++) {
if (p_row == 0 && p_col == 0)
continue;
copy_A_to_SUB_A(p_row, p_col);
dgesd2d_(&ICTXT, &SUB_A_ROWS, &SUB_A_COLS, SUB_A.data(), &LDA, &p_row,
&p_col);
}
}
copy_A_to_SUB_A(0, 0); // set local submatrix for rank 0
} else {
dgerv2d_(&ICTXT, &SUB_A_ROWS, &SUB_A_COLS, SUB_A.data(), &LDA, &RSRC,
&CSRC);
}
-
copy_A_to_SUB_A
はA行列のうち、他のプロセスに転送する部分をSUB_A
にコピーする処理。 - 同様のことをBに対しても行っている。
また、これより以前の場所でdescinit_
も読んでおく必要がある。DESCA
,DESCB
は行列A,Bのdescriptorと呼ばれるもので、どのように行列が分散して保持されているかという情報を有する長さ9のintの配列である。
// Initialize array descriptors
int DESCA[9], DESCB[9];
int INFO;
descinit_(DESCA, &M, &N, &MB, &NB, &RSRC, &CSRC, &ICTXT, &SUB_A_ROWS, &INFO);
descinit_(DESCB, &N, &NRHS, &MB, &NRHS, &RSRC, &CSRC, &ICTXT, &SUB_B_ROWS, &INFO);
求解
ここまでの準備ができてしまえば、実際に解を求める部分はとてもシンプルである。pdgesv_
を呼ぶだけでよい。
//*****
// Call the SCALAPACK routine
//*****
int IA = 1;
int JA = 1;
int IB = 1;
int JB = 1;
pdgesv_(&N, &NRHS, SUB_A.data(), &IA, &JA, DESCA, IPIV.data(), SUB_B.data(),
&IB, &JB, DESCB, &INFO);
- scalapackの関数名の命名則はlapackとほぼ同等である。lapackの関数名に接頭辞
p
をつけるだけで該当する関数が見つけられる。 -
N
は行列Aのサイズ。NRHS
は行列Bの列のサイズ。ここでは1だが、複数列を持つ行列をBに指定することで、複数のBベクトルに対して同時に解を求めることができる。 -
IA
,JA
,IB
,JB
は行列A,Bの最初の要素のインデックス。通常1を指定すれば良い。 -
IPIV
は解を求める際に行や列がどのように入れ替えられたかという情報が格納される配列。 長さSUB_A_ROWS + NB
のintの配列を事前に用意しておく必要がある。 -
INFO
はリターンコードが入る。INFO == 0
ならば正常終了している。
解の出力と検証
求解したら、解$X$は配列SUB_B
に保存されている。このコードでは結果をランク0に集めて出力している。
if (MYROW == 0 && MYCOL == 0) {
copy_SUB_B_to_X(0, 0);
for (int p_row = 1; p_row < NPROW; p_row++) {
int p_col = 0;
dgerv2d_(&ICTXT, &SUB_B_ROWS, &SUB_B_COLS, SUB_B.data(), &LDA, &p_row,
&p_col);
copy_SUB_B_to_X(p_row, p_col);
}
} else if (MYCOL == 0) {
dgesd2d_(&ICTXT, &SUB_B_ROWS, &SUB_B_COLS, SUB_B.data(), &LDA, &RSRC,
&CSRC);
}
各プロセスが持っているSUB_B
をランク0のプロセスのSUB_B
にコピーして、最終的な解の配列X
にコピーしている。
最後に本当に$Ax = B$となっているか確認のため、$Ax$を計算して出力している。
std::cout << "Matrix AX:\n";
for (int i = 0; i < M; i++) {
double ans = 0.0;
for (int j = 0; j < N; j++) {
ans += A[i][j] * X[j];
}
std::cout << ans << ' ';
}
std::cout << std::endl;
C++ラッパー
(2021/12/5 追記)
C++のクラスとしてラップしたコードを以下で公開した。PDGESVを比較的わかりやすいインターフェースで使えるようになっている。(すべてのlapackルーチンが使えるわけではない)