MetaMat.hpp

Go to the documentation of this file.

/*******************************************************************************
 * Copyright (C) 2017-2024 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/>.
 ******************************************************************************/
#ifndef METAMAT_HPP
#define METAMAT_HPP
 
#include "triplet_form.hpp"
#include "IterativeSolver.hpp"
#include "ILU.hpp"
#include "Jacobi.hpp"
 
template<typename T, typename U> concept ArmaContainer = std::is_floating_point_v<U> && (std::is_convertible_v<T, Mat<U>> || std::is_convertible_v<T, SpMat<U>>) ;
 
template<sp_d T> class MetaMat;
 
template<sp_d T> class op_add {
    friend MetaMat<T>;
 
    shared_ptr<MetaMat<T>> mat_a, mat_b;
 
public:
    explicit op_add(const shared_ptr<MetaMat<T>>& A)
        : mat_a(A)
        , mat_b(nullptr) {}
 
    op_add(const shared_ptr<MetaMat<T>>& A, const shared_ptr<MetaMat<T>>& B)
        : mat_a(A)
        , mat_b(B) {}
};
 
template<sp_d T> class op_scale {
    friend MetaMat<T>;
 
    T scalar;
    op_add<T> bracket;
 
public:
    op_scale(const T A, const shared_ptr<MetaMat<T>>& B)
        : scalar(A)
        , bracket(B) {}
 
    op_scale(const T A, op_add<T>&& B)
        : scalar(A)
        , bracket(std::forward<op_add<T>>(B)) {}
};
 
template<sp_d T> class MetaMat {
protected:
    bool factored = false;
 
    SolverSetting<T> setting{};
 
    virtual int direct_solve(Mat<T>&, const Mat<T>&) = 0;
 
    virtual int direct_solve(Mat<T>&, Mat<T>&&) = 0;
 
    int direct_solve(Mat<T>& X, const SpMat<T>& B) { return this->direct_solve(X, Mat<T>(B)); }
 
    int direct_solve(Mat<T>& X, SpMat<T>&& B) { return this->direct_solve(X, B); }
 
    int iterative_solve(Mat<T>&, const Mat<T>&);
 
    int iterative_solve(Mat<T>& X, const SpMat<T>& B) { return this->iterative_solve(X, Mat<T>(B)); }
 
    template<std::invocable<fmat&> F> int mixed_trs(mat& X, mat&& B, F trs) {
        auto INFO = 0;
 
        X = arma::zeros(size(B));
 
        auto multiplier = norm(B);
 
        auto counter = 0u;
        while(counter++ < this->setting.iterative_refinement) {
            if(multiplier < this->setting.tolerance) break;
 
            auto residual = conv_to<fmat>::from(B / multiplier);
 
            if(0 != (INFO = trs(residual))) break;
 
            const mat incre = multiplier * conv_to<mat>::from(residual);
 
            X += incre;
 
            suanpan_debug("Mixed precision algorithm multiplier: {:.5E}.\n", multiplier = arma::norm(B -= this->operator*(incre)));
        }
 
        return INFO;
    }
 
public:
    triplet_form<T, uword> triplet_mat;
 
    const uword n_rows;
    const uword n_cols;
    const uword n_elem;
 
    MetaMat(const uword in_rows, const uword in_cols, const uword in_elem)
        : triplet_mat(in_rows, in_cols)
        , n_rows(in_rows)
        , n_cols(in_cols)
        , n_elem(in_elem) {}
 
    MetaMat(const MetaMat&) = default;
    MetaMat(MetaMat&&) noexcept = delete;
    MetaMat& operator=(const MetaMat&) = delete;
    MetaMat& operator=(MetaMat&&) noexcept = delete;
    virtual ~MetaMat() = default;
 
    void set_solver_setting(const SolverSetting<T>& SS) { setting = SS; }
 
    [[nodiscard]] SolverSetting<T>& get_solver_setting() { return setting; }
 
    void set_factored(const bool F) { factored = F; }
 
    [[nodiscard]] virtual bool is_empty() const = 0;
    virtual void zeros() = 0;
 
    virtual unique_ptr<MetaMat> make_copy() = 0;
 
    void unify(const uword K) {
        this->nullify(K);
        this->at(K, K) = T(1);
    }
 
    virtual void nullify(uword) = 0;
 
    [[nodiscard]] virtual T max() const = 0;
    [[nodiscard]] virtual Col<T> diag() const = 0;
 
    virtual T operator()(uword, uword) const = 0;
    virtual T& unsafe_at(const uword I, const uword J) { return this->at(I, J); }
 
    virtual T& at(uword, uword) = 0;
 
    [[nodiscard]] virtual const T* memptr() const = 0;
    virtual T* memptr() = 0;
 
    virtual void scale_accu(T, const shared_ptr<MetaMat>&) = 0;
    virtual void scale_accu(T, const triplet_form<T, uword>&) = 0;
 
    void operator+=(const shared_ptr<MetaMat>& M) { return this->scale_accu(1., M); }
 
    void operator-=(const shared_ptr<MetaMat>& M) { return this->scale_accu(-1., M); }
 
    void operator+=(const op_scale<T>& M) {
        const auto& bracket = M.bracket;
        if(nullptr != bracket.mat_a) this->scale_accu(M.scalar, bracket.mat_a);
        if(nullptr != bracket.mat_b) this->scale_accu(M.scalar, bracket.mat_b);
    }
 
    void operator-=(const op_scale<T>& M) {
        const auto& bracket = M.bracket;
        if(nullptr != bracket.mat_a) this->scale_accu(-M.scalar, bracket.mat_a);
        if(nullptr != bracket.mat_b) this->scale_accu(-M.scalar, bracket.mat_b);
    }
 
    void operator+=(const triplet_form<T, uword>& M) { return this->scale_accu(1., M); }
 
    void operator-=(const triplet_form<T, uword>& M) { return this->scale_accu(-1., M); }
 
    virtual Mat<T> operator*(const Mat<T>&) const = 0;
 
    virtual void operator*=(T) = 0;
 
    template<ArmaContainer<T> C> int solve(Mat<T>& X, C&& B) { return IterativeSolver::NONE == this->setting.iterative_solver ? this->direct_solve(X, std::forward<C>(B)) : this->iterative_solve(X, std::forward<C>(B)); }
 
    template<ArmaContainer<T> C> Mat<T> solve(C&& B) {
        Mat<T> X;
 
        if(SUANPAN_SUCCESS != this->solve(X, std::forward<C>(B))) throw std::runtime_error("fail to solve the system");
 
        return X;
    }
 
    [[nodiscard]] virtual int sign_det() const = 0;
 
    void save(const char* name) {
        if(!to_mat(*this).save(name, raw_ascii))
            suanpan_error("Cannot save to file \"{}\".\n", name);
    }
 
    virtual void csc_condense() {}
 
    virtual void csr_condense() {}
 
    [[nodiscard]] Col<T> evaluate(const Col<T>& X) const { return this->operator*(X); }
};
 
template<sp_d T> int MetaMat<T>::iterative_solve(Mat<T>& X, const Mat<T>& B) {
    this->csc_condense();
 
    X.zeros(arma::size(B));
 
    unique_ptr<Preconditioner<T>> preconditioner;
    if(PreconditionerType::JACOBI == this->setting.preconditioner_type) preconditioner = std::make_unique<Jacobi<T>>(this->diag());
#ifndef SUANPAN_SUPERLUMT
    else if(PreconditionerType::ILU == this->setting.preconditioner_type) {
        if(this->triplet_mat.is_empty()) preconditioner = std::make_unique<ILU<T>>(to_triplet_form<T, int>(this));
        else preconditioner = std::make_unique<ILU<T>>(this->triplet_mat);
    }
#endif
    else if(PreconditionerType::NONE == this->setting.preconditioner_type) preconditioner = std::make_unique<UnityPreconditioner<T>>();
 
    if(SUANPAN_SUCCESS != preconditioner->init()) return SUANPAN_FAIL;
 
    this->setting.preconditioner = preconditioner.get();
 
    std::atomic_int code = 0;
 
    if(IterativeSolver::GMRES == setting.iterative_solver)
        suanpan::for_each(B.n_cols, [&](const uword I) {
            Col<T> sub_x(X.colptr(I), X.n_rows, false, true);
            const Col<T> sub_b(B.colptr(I), B.n_rows);
            auto col_setting = setting;
            code += GMRES(this, sub_x, sub_b, col_setting);
        });
    else if(IterativeSolver::BICGSTAB == setting.iterative_solver)
        suanpan::for_each(B.n_cols, [&](const uword I) {
            Col<T> sub_x(X.colptr(I), X.n_rows, false, true);
            const Col<T> sub_b(B.colptr(I), B.n_rows);
            auto col_setting = setting;
            code += BiCGSTAB(this, sub_x, sub_b, col_setting);
        });
    else throw invalid_argument("no proper iterative solver assigned but somehow iterative solving is called");
 
    return 0 == code ? SUANPAN_SUCCESS : SUANPAN_FAIL;
}
 
template<sp_d T> Mat<T> to_mat(const MetaMat<T>& in_mat) {
    Mat<T> out_mat(in_mat.n_rows, in_mat.n_cols);
    for(uword J = 0; J < in_mat.n_cols; ++J) for(uword I = 0; I < in_mat.n_rows; ++I) out_mat(I, J) = in_mat(I, J);
    return out_mat;
}
 
template<sp_d T> Mat<T> to_mat(const shared_ptr<MetaMat<T>>& in_mat) { return to_mat(*in_mat); }
 
template<sp_d data_t, sp_i index_t> Mat<data_t> to_mat(const triplet_form<data_t, index_t>& in_mat) {
    Mat<data_t> out_mat(in_mat.n_rows, in_mat.n_cols, fill::zeros);
    for(index_t I = 0; I < in_mat.n_elem; ++I) out_mat(in_mat.row(I), in_mat.col(I)) += in_mat.val(I);
    return out_mat;
}
 
template<sp_d data_t, sp_i index_t> Mat<data_t> to_mat(const csr_form<data_t, index_t>& in_mat) {
    Mat<data_t> out_mat(in_mat.n_rows, in_mat.n_cols, fill::zeros);
 
    index_t c_idx = 1;
    for(index_t I = 0; I < in_mat.n_elem; ++I) {
        if(I >= in_mat.row_mem()[c_idx]) ++c_idx;
        out_mat(c_idx - 1, in_mat.col_mem()[I]) += in_mat.val_mem()[I];
    }
 
    return out_mat;
}
 
template<sp_d data_t, sp_i index_t> Mat<data_t> to_mat(const csc_form<data_t, index_t>& in_mat) {
    Mat<data_t> out_mat(in_mat.n_rows, in_mat.n_cols, fill::zeros);
 
    index_t c_idx = 1;
    for(index_t I = 0; I < in_mat.n_elem; ++I) {
        if(I >= in_mat.col_mem()[c_idx]) ++c_idx;
        out_mat(in_mat.row_mem()[I], c_idx - 1) += in_mat.val_mem()[I];
    }
 
    return out_mat;
}
 
template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t> to_triplet_form(MetaMat<data_t>* in_mat) {
    if(!in_mat->triplet_mat.is_empty()) return triplet_form<data_t, index_t>(in_mat->triplet_mat);
 
    const sp_i auto n_rows = index_t(in_mat->n_rows);
    const sp_i auto n_cols = index_t(in_mat->n_cols);
    const sp_i auto n_elem = index_t(in_mat->n_elem);
 
    triplet_form<data_t, index_t> out_mat(n_rows, n_cols, n_elem);
    for(index_t J = 0; J < n_cols; ++J) for(index_t I = 0; I < n_rows; ++I) out_mat.at(I, J) = in_mat->operator()(I, J);
 
    return out_mat;
}
 
template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t> to_triplet_form(const shared_ptr<MetaMat<data_t>>& in_mat) { return to_triplet_form<data_t, index_t>(in_mat.get()); }
 
#endif