Fabulous installation procedure

The matrix class will be simple, we will use a wrapper of the dense matrices of fabulous

class Matrix
{
private:
    Block _data;

To simplify, we add a function to load a matrix from a file

void load_matrix(const std::string &filename)
{
    fa::MatrixMarketLoader mml;
    mml.LoadMatrix(_data, filename);
}

In our case we only have one constructor using a file.

public:
    using value_type = std::complex<double>;
    using primary_type = double;

    explicit Matrix(const std::string &filename)
    {
        load_matrix(filename);
    }

The size of your matrix (for display mostly).

int size() const { return _data.get_nb_row(); }

The matrix vector product (required by fabulous), in our case it is the wrapper to the lapack gemm.

// B := \alpha * A * X + \beta * B
int64_t operator()(const Block &X, Block &B, S alpha=S{1.0}, S beta=S{0.0}) const
{
    int M = B.get_nb_row();
    int N = B.get_nb_col();
    int K = X.get_nb_row();
    fa::lapacke::gemm(
        M, N, K,
        _data.get_ptr(), _data.get_leading_dim(),
        X.get_ptr(), X.get_leading_dim(),
        B.get_ptr(), B.get_leading_dim(),
        alpha, beta
    );
    return fa::lapacke::flops::gemm<S>(M, N, K);
}

Fist let's set the solver accuracy, for instance, how small do we want \(\frac{||B-AX||}{||B||}\) ? To do so in fabulous you need to specify an array with either, \(1\) element if the accuracy is the same for all right-hand-sides, or a specific one for each right-hand-side.

In our let's take the same one for all of the right-hand-sides.

std::vector<double> epsilon{1e-2};

Iterative methods in fabulous (but also in general) may no converge and running indefinitely. In that regard we need a fail-safe.

Fabulous counts the number of matrix vector product (mvp) and the number of iterations separately. The mvp accounts for the number of right-hand-sides of \(X\) when you do a \(Y = AX\) The iteration describe the total number of step the algorithm needed to do to converge. If the number of right-hand-side is \(1\) then the mvp is basically the number of iterations.

The fail-safe over the mvp is a better indicator for multiple right-hand-side algorithm.

const int max_mvp = 200;

Fabulous implement Krylov space methods, those methods need a maximum space size for searching the solution. This maximum basically represent the memory cost of the algorithm, it needs to be big enough to be able to find the solution. After the Krylov space is at it's maximum size a restart occur, that is we compute the partial solution and empty the Krylov space.

const int max_krylov_space = 50;

Some algorithm have a feature that is called deflated restart (-DR variant of algorithm). The difference with a usual restart (also called cold restart) is that some spectral information of the Krylov space is keep. This is translated by choosing the number of eigen vector you want to kept.

Note: with real matrices the eigen values may be complex, thus the number of eigen vector might be implicitly increase by \(1\)

const int nb_eigen_pair = 5;

Krylov space methods also require a way to orthogonalize a base of the Krylov space. There are \(2\) types of orthogonality in fabulous, BLOCK for orthogonalizing all the right-hand-side at the same time and RUHE for orthogonalizing one right-hand-side at a time. Also there are \(4\) schemes of orthogonality in fabulous, Classical Gram-Schmidt (CGS), Modified Gram-Schmidt (MGS) and their iterative extension (ICGS, IMGS).

Let's pick BLOCK and CGS.

auto ortho = fa::OrthoType::BLOCK + fa::OrthoScheme::CGS;

Now that we have our matrices we want to solve the equation \(AX = B\).

To do so we need first to declare a placeholder, like

fa::noop_placeholder _;

This placeholder allows to declare an equation without needing to specify all the parameters.

auto eq = fa::equation(A, _, _, X, B, epsilon);

To be more precise, here are all the parameters

fa::equation(A, M, CB, X, B, epsilon)

• M: the preconditionner of the form AMy = b with My = x.
• CB: a callback function of prototype void callback(int, Logger&, Vector, Vector) called after each iteration.

All solver share the same parameters structure to pass all parameters at ones

auto param = fa::parameters(max_mvp, max_krylov_space);
//param.set_quiet(true); // turn off the solvers' verbose

At first let's choose the usual BGMRES algorithm (Block-GMRES)

auto gmres = fa::bgmres::std();
auto log = fa::bgmres::solve(eq, gmres, param, ortho);
//log.print("gmres"); // this prints comma-separated columns with information about convergence

Now let's change the orthogonalization method.

ortho = fa::OrthoType::RUHE + fa::OrthoScheme::MGS;

Let's do a new run to see

X.zero(); // Reset X
log = fa::bgmres::solve(eq, gmres, param, ortho);
//log.print("gmres_ruhe_mgs");

The IB prefix signify that, in the case where some right-hand-sides have already converged, they will be detected and treated as such.

X.zero(); // Reset X
auto ib_gmres = fa::bgmres::ib();
log = fa::bgmres::solve(eq, ib_gmres, param, ortho);
//log.print("ib_gmres");

As explain in Section Krylov space method, IB-BGMRES-DR is the IB-BGMRES with deflated restart

X.zero(); // Reset X
auto restart = fa::deflated_restart(nb_eigen_pair, S{0.0});
auto ib_gmres_dr = fa::bgmres::ibdr();
log = fa::bgmres::solve(eq, ib_gmres_dr, param, ortho, restart);
//log.print("ib_gmres_dr");

IB-BGCRO-DR is an algorithm that goes faster after each solve by saving spectral information in a block vector, \(U_k\) in our case

X.zero(); // Reset X

Block Uk{}; // delation space storage, /!\ start empty /!\

auto ib_gcro_dr = fa::bgcro::ibdr();
log = fa::bgcro::solve(eq, ib_gcro_dr, param, Uk, ortho, restart);
//log.print("ib_gcro_dr_1");

X.zero(); // Reset X
rnl.LoadMatrix(B); // new right-hand-side
log = fa::bgcro::solve(eq, ib_gcro_dr, param, Uk, ortho, restart);
//log.print("ib_gcro_dr_2");