Skip to content

Custom Matrix Class

ezp supports customised matrix objects.

Apart from using the provided wrappers such as full_mat, band_mat and band_symm_mat, one can also define or use customised matrix classes to call the solver.

All custom matrix classes shall have the public member .n_rows and .n_cols. The band matrix further requires the public member .kl and .ku which are the numbers of sub-diagonals (lower bandwidth) and super-diagonals (upper bandwidth). The band symmetric matrix requires the additional public member .klu which is the bandwidth.

All custom classes must define at least one of the following public methods, that return a pointer to the first element. It is assumed the memory layout is contiguous, thus, one can use std::vector<T>, std::array<T>, std::unique_ptr<T[]> as the storage.

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

The possibility is unlimited. The simplest is to subclass std::vector<T> with additional public members. In the example of posv solver, a custom matrix class is used.

example.pposv.cpp
 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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#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:
    auto init(const int_t rows, const int cols) {
        n_rows = rows;
        n_cols = cols;
        storage.resize(rows * cols);
    }

    explicit mat(const int_t rows = 0, const int cols = 0) { init(rows, cols); }

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

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

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

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

int main() {
    // get the current blacs environment
    const auto& env = get_env<int_t>();

    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 of above methods shall return a pointer to the first element
    mat A, B;

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

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

    // 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;
}

//! @}