Detailed Description

Solver for general band matrices.

Note: Although the pgbsv solver supports KL=0 and/or KU=0, a zero (half) bandwidth would lead to unwanted warning message from ScaLAPACK.; See: https://github.com/Reference-ScaLAPACK/scalapack/issues/116

It solves the system of linear equations A * X = B with a general band matrix A. The band matrix A has KL sub-diagonals and KU super-diagonals. It shall be stored in the following format. The band storage scheme is illustrated by the following example, when M=N=6, KL=2, KU=1.

 .    .    .    .    .    .
 .    .    .    .    .    .
 .    .    .    .    .    .
 .   a12  a23  a34  a45  a56
a11  a22  a33  a44  a55  a66
a21  a32  a43  a54  a65   .
a31  a42  a53  a64   .    .

The lead dimension should be 2*(KL+KU)+1.

With zero based indexing, for a general band matrix A, the element at row i and column j is stored at A[IDX(i, j)].

const auto IDX = [&](const int i, const int j) {
    if(i - j > KL || j - i > KU) return -1;
    return 2 * KU + KL + i + 2 * j * (KL + KU);
};

The example usage can be seen as follows.

/*******************************************************************************
 * Copyright (C) 2025-2026 Theodore Chang
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ******************************************************************************/
#include <ezp/pgbsv.hpp>
#include <iomanip>
#include <iostream>
 
using namespace ezp;
 
int main() {
    // get the current blacs environment
    const auto& env = get_env<int_t>();
 
    constexpr auto N = 10, NRHS = 1, KL = 2, KU = 2;
    constexpr auto LDA = 2 * (KL + KU) + 1;
 
    // storage for the matrices A and B
    std::vector<double> A, B;
 
    // helper function to convert 2D indices to 1D indices
    // the band symmetric matrix used for gbsv subroutine requires the matrix to be stored with a leading dimension of (2 * (KL + KU) + 1)
    // see Fig. 4.10 https://netlib.org/scalapack/slug/node84.html
    const auto IDX = par_dgbsv<int_t>::indexer{N, KL, KU};
 
    if(0 == env.rank()) {
        // the matrices are only initialized on the root process
        A.resize(N * LDA, 0.);
        B.resize(N * NRHS, 1.);
 
        for(auto I = 0; I < N; ++I) A[IDX(I, I)] = I + 1;
    }
 
    // create a parallel solver
    // it uses a one-dimensional process grid
    // it takes the number of processes as arguments
    auto solver = par_dgbsv(env.size());
 
    // need to wrap the data in full_mat objects
    // it requires the number of rows and columns of the matrix, and a pointer to the data
    // on non-root processes, the data pointer is nullptr as the vector is empty
    // solver.solve(band_mat{N, N, KL, KU, A.data()}, full_mat{N, NRHS, B.data()});
    const auto info = solver.solve({N, N, KL, KU, A.data()}, {N, NRHS, B.data()});
 
    if(0 == env.rank() && 0 == info) {
        std::cout << std::setprecision(6) << std::fixed << "Info: " << info << '\n';
        std::cout << "Solution:\n";
        for(const double i : B) std::cout << i << '\n';
    }
 
    return info;
}
 

Author: tlc

Date: 07/03/2025

Version: 1.0.0

The documentation for this class was generated from the following file:

ezp/pgbsv.hpp