scalapack入門 (in C++)

はじめに

scalapackとはマルチプロセス環境で線形代数計算を実行するためのライブラリである。lapackの並列版といえるオープンソースのライブラリで、多くのスパコンにプリインストールされている。

この記事ではscalapackを導入し、線型方程式を解くところまでを解説する。
行列を複数のプロセスで分散して持つことになるが、どのように分散するかということも含めて解説していく。また、lapackと同様にscalapackはFortranで実装されているため使い方にクセがあり、特にC,C++から使う時はネット上に情報が不足している。
この記事ではC++からscalapackを使うための手順を示す。

まずはscalapackをインストールしてサンプルコードを動かすまでの手順を示し、その後サンプルコードの解説を行っていく。

インストール

ここではmac上でscalapackを利用するための手順を示す。homebrewで簡単にインストールすることができる。

$ brew install open-mpi         # MPI環境をインストールしていない場合はインストールする
$ brew install scalapack

サンプルコード

以下のサンプルコードを準備する。
ちなみにこちらのサンプルコードでは変数を出力するのに便利なように icecream-cppというライブラリを使用している。
こちらの記事で、このライブラリの使用方法を紹介しているので参考にして欲しい。

pdgesv_example.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 [e-mail]: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プロセスが２次元の格子状に配置されるということである。
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の１次元配列として持ち、先頭から(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の領域を確保し、さらに５行目に相当する箇所は空けておく必要がある。つまり(4,0)に相当する箇所は空白にしておき、(0,1)に対応する部分は配列の６番目(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となっているはず。
• プログラムの終了時には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の列のサイズ。ここでは１だが、複数列を持つ行列を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ルーチンが使えるわけではない）