Skip to content

Custom Matrix Class

ezp supports customised matrix objects. Besides the provided wrappers full_mat, band_mat, and band_symm_mat, one may define or reuse custom matrix classes to call the solvers. All custom matrix classes must expose public members .n_rows and .n_cols. A band matrix must additionally expose .kl and .ku, the numbers of sub-diagonals (lower bandwidth) and super-diagonals (upper bandwidth). A band symmetric matrix must expose the additional member .klu, the bandwidth.

Every custom class must define at least one of the following public methods, which return a pointer to the first element. The memory layout is assumed to be contiguous. std::vector<T>, std::array<T>, and std::unique_ptr<T[]> are therefore all valid storages.

  1. .mem() -> T*
  2. .memptr() -> T*
  3. .data() -> T*
  4. a contiguous iterator pair .begin() and .end(), so that &(*begin()) -> T*

The possibilities are unlimited. The simplest approach is to subclass std::vector<T> and add the required public members. The posv example below demonstrates a custom matrix class.

example.pposv.cpp:27:54
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#include <ezp/pposv.hpp>
#include <iomanip>
#include <iostream>

using namespace ezp;

class mat {
public:
    int_t n_rows, n_cols;

private:
    std::vector<double> storage;

public:
    mat(const int_t rows, const int cols)
        : n_rows(rows)
        , n_cols(cols) {}

    auto init() { storage.resize(n_rows * n_cols); }

    auto& operator()(const int_t i, const int_t j) { return storage[i + j * n_cols]; }

    auto& operator[](const int_t i) { return storage[i]; }

    auto begin() { return storage.begin(); }

    auto end() { return storage.end(); }
};

All processes must provide the same dimensional inputs. Thus .n_rows and .n_cols must be valid on every process. The actual storage, however, is only initialised on the root process.

example.pposv.cpp:56:91
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
int main() {
    // get the current blacs environment
    const auto& env = get_env<>();

    constexpr auto N = 6, NRHS = 2;

    // storage for the matrices A and B using the custom class
    // acceptable objects shall have members `.n_rows` and `.n_cols`
    // and have any of the following methods `.mem()`, `.memptr()`, `.data()`
    // or contiguous iterators
    // all above methods shall return a pointer to the first element
    mat A(N, N), B(N, NRHS);

    if(0 == env.rank()) {
        // the matrices are only initialized on the root process
        A.init();
        B.init();

        for(auto I = 0; I < N; ++I) {
            B[I] = A(I, I) = I + 1;
            B[I + N] = I + 4;
        }
    }

    // create a parallel solver
    // and send custom matrix objects to the solver
    const auto info = par_dposv<int_t>().solve(A, B);

    if(0 == env.rank() && 0 == info) {
        std::cout << std::setprecision(10) << "Info: " << info << '\n';
        std::cout << "Solution:\n";
        for(const auto i : B) std::cout << i << '\n';
    }

    return info;
}