ppbsv Class Reference

Solver for symmetric band positive definite matrices. More...

#include <ppbsv.hpp>

Detailed Description

Solver for symmetric band positive definite matrices.

Note

Although the ppbsv solver supports KLU=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 symmetric band positive definite matrix A. The band matrix A has KLU sub-diagonals. It shall be stored in the following format. The band storage scheme is illustrated by the following example, when M=N=6, KLU=2.

For ‘UPLO='L’`, the lower half is stored.

a11  a22  a33  a44  a55  a66
a21  a32  a43  a54  a65   .
a31  a42  a53  a64   .    .

The lead dimension should be KLU+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 = [&](int i, int j) {
    if(i < j) std::swap(i, j);
    if(i - j > KLU) return -1;
    return i + j * KLU;
};

For ‘UPLO='U’`, the upper half is stored.

 .    .   a13  a24  a35  a46
 .   a12  a23  a34  a45  a56
a11  a22  a33  a44  a55  a66

The lead dimension should be KLU+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 = [&](int i, int j) {
    if(i > j) std::swap(i, j);
    if(j - i > KLU) return -1;
    return 2 * j - i + (j + 1) * KLU;
};

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/ppbsv.hpp>
#include <iomanip>
#include <iostream>
 
// get the current blacs environment
const auto& env = ezp::get_env<>();
 
class general_mat {
public:
    int_t n_rows{}, n_cols{};
 
private:
    std::vector<double> storage;
 
public:
    auto init(const int_t rows, const int cols) {
        n_rows = rows;
        n_cols = cols;
        storage.resize(rows * cols);
    }
 
    explicit general_mat(const int_t rows = 0, const int cols = 0) { init(rows, cols); }
 
    auto& operator[](const int_t i) { return storage[i]; }
 
    auto begin() { return storage.begin(); }
 
    auto end() { return storage.end(); }
};
 
class bandsymm_mat {
public:
    int_t n_rows{}, n_cols{}, klu{};
 
private:
    std::vector<double> storage;
 
    ezp::ppbsv<double, int_t>::indexer indexer{0, 0};
 
public:
    auto init(const int_t n, const int bandwidth) {
        n_rows = n;
        n_cols = n;
        klu = bandwidth;
        if(0 == env.rank()) storage.resize(n * (klu + 1));
        indexer = {n, klu};
    }
 
    explicit bandsymm_mat(const int_t n = 0, const int bandwidth = 0) { init(n, bandwidth); }
 
    auto& operator()(const int_t i, const int_t j) { return storage[indexer(i, j)]; }
 
    auto data() { return storage.data(); }
};
 
int main() {
    // storage for the matrices A and B
    bandsymm_mat A;
    general_mat B;
 
    constexpr auto N = 6, NRHS = 2, KLU = 1;
    // the matrices are only initialized on the root process
    A.init(N, KLU);
    B.init(N, NRHS);
 
    if(0 == env.rank()) {
        static constexpr auto M = 5.10156648;
 
        for(auto I = 0; I < N; ++I) {
            B[I] = A(I, I) = I + 1;
            B[N + I] = (I + 1) * M;
        }
    }
 
    // create a parallel solver
    // it uses a one-dimensional process grid
    // it takes the number of processes as arguments
    auto solver = ezp::par_dpbsv(env.size());
 
    // use custom matrix objects
    const auto info = solver.solve(A, B);
 
    if(0 == env.rank() && 0 == info) {
        std::cout << std::setprecision(10) << "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/ppbsv.hpp