block_precond.h

#ifndef BLOCK_PRECOND_H
#define BLOCK_PRECOND_H
 
#include <stdexcept>
#include <vector>
 
#include <eigen3/Eigen/Dense>
 
#include "sparseMat.h"
 
namespace algebra
{
template <int BS>
class BlockDiagPrecond
    {
public:
    BlockDiagPrecond(const SparseMatrix &A, const int nbBlocks, const std::vector<int> &ld)
        : inv(nbBlocks)
        {
        std::vector<char> dir(BS * nbBlocks, 0);
        for (const int i : ld)
            { dir[i] = 1; }
        for (int n = 0; n < nbBlocks; n++)
            {
            Eigen::Matrix<double, BS, BS> a;
            for (int r = 0; r < BS; r++)
                for (int c = 0; c < BS; c++)
                    {
                    const int I = BS * n + r, J = BS * n + c;
                    a(r, c) = (dir[I] || dir[J]) ? (r == c ? 1.0 : 0.0) : A(I, J);
                    }
            if (a.determinant() == 0.0)
                { throw std::runtime_error("BlockDiagPrecond: singular diagonal block"); }
            inv[n] = a.inverse();
            }
        }
 
    void apply(const std::vector<double> &x, std::vector<double> &y) const
        {
        for (size_t n = 0; n < inv.size(); n++)
            {
            Eigen::Map<const Eigen::Matrix<double, BS, 1>> xb(&x[BS * n]);
            Eigen::Map<Eigen::Matrix<double, BS, 1>> yb(&y[BS * n]);
            yb = inv[n] * xb;
            }
        }
 
private:
    std::vector<Eigen::Matrix<double, BS, BS>> inv;
    };
 
}  // namespace algebra
#endif