iterative_solver_base.h

Go to the documentation of this file.

/*
 * Copyright 2023 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include <Eigen/Core>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/SparseCore>
 
#include "num_collect/base/concepts/dense_vector_of.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/precondition.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/util/assert.h"
 
namespace num_collect::linear {
 
namespace impl {
 
template <typename Solver, typename /*reserved for SFINAE*/ = void>
struct iterative_solver_traits;  // IWYU pragma: keep
 
}  // namespace impl
 
template <typename Derived>
class iterative_solver_base : public Eigen::SparseSolverBase<Derived> {
protected:
    using Base = Eigen::SparseSolverBase<Derived>;
    using Base::m_isInitialized;
 
public:
    using MatrixType =
        typename impl::iterative_solver_traits<Derived>::matrix_type;
    using Scalar = typename MatrixType::Scalar;
    using StorageIndex = typename MatrixType::StorageIndex;
    using RealScalar = typename MatrixType::RealScalar;
 
    using matrix_type = MatrixType;
    using scalar_type = Scalar;
    using storage_index_type = StorageIndex;
    using real_scalar_type = RealScalar;
 
    enum {  // NOLINT(performance-enum-size): Preserve the same implementation as Eigen library.
        ColsAtCompileTime = MatrixType::ColsAtCompileTime,  // NOLINT
        MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime  // NOLINT
    };
 
    using Base::derived;
 
    iterative_solver_base() { m_isInitialized = false; }
 
    auto compute(const matrix_type& coeff) -> Derived& {
        coeff_ = &coeff;
        m_isInitialized = true;
        return derived();
    }
 
    [[nodiscard]] auto rows() const noexcept -> Eigen::Index {
        return coeff().rows();
    }
 
    [[nodiscard]] auto cols() const noexcept -> Eigen::Index {
        return coeff().cols();
    }
 
    [[nodiscard]] auto tolerance() const noexcept -> real_scalar_type {
        return tolerance_;
    }
 
    auto tolerance(const real_scalar_type& val) -> Derived& {
        NUM_COLLECT_PRECONDITION(val > static_cast<real_scalar_type>(0),
            "Tolerance of rate of residual must be positive value.");
        tolerance_ = val;
        return derived();
    }
 
    [[nodiscard]] auto max_iterations() const noexcept -> index_type {
        return max_iterations_;
    }
 
    auto max_iterations(index_type val) -> Derived& {
        NUM_COLLECT_PRECONDITION(val > 0,
            "Maximum number of iterations must be a positive integer.");
        max_iterations_ = val;
        return derived();
    }
 
    template <base::concepts::dense_vector_of<scalar_type> Right>
    [[nodiscard]] auto solve(const Right& right) const
        -> Eigen::Solve<Derived, Right> {
        // TODO: Version for matrices.
        return Eigen::Solve<Derived, Right>(derived(), right);
    }
 
    template <base::concepts::dense_vector_of<scalar_type> Right,
        base::concepts::dense_vector_of<scalar_type> Solution>
    [[nodiscard]] auto solve_with_guess(
        const Right& right, const Solution& solution) const
        -> Eigen::SolveWithGuess<Derived, Right, Solution> {
        // TODO: Version for matrices.
        return Eigen::SolveWithGuess<Derived, Right, Solution>(
            derived(), right, solution);
    }
 
    template <base::concepts::dense_vector_of<scalar_type> Right,
        base::concepts::dense_vector_of<scalar_type> Solution>
    void _solve_impl(  // NOLINT(readability-identifier-naming
                       // name required by Eigen.
        const Right& right, Solution& solution) const {
        // TODO: Version for matrices.
        solution.setZero();
        derived()._solve_with_guess_impl(right, solution);
    }
 
    template <base::concepts::dense_vector_of<scalar_type> Right,
        base::concepts::dense_vector_of<scalar_type> Solution>
    void _solve_with_guess_impl(  // NOLINT(readability-identifier-naming
                                  // name required by Eigen.
        const Right& right, Solution& solution) const {
        // TODO: Version for matrices.
        derived().solve_vector_in_place(right, solution);
    }
 
protected:
    [[nodiscard]] auto coeff() const noexcept -> const matrix_type& {
        NUM_COLLECT_ASSERT(static_cast<const void*>(coeff_) != nullptr);
        return *coeff_;
    }
 
private:
    static constexpr index_type default_max_iterations = 10000;
 
    index_type max_iterations_{default_max_iterations};
 
    static constexpr auto default_tolerance =
        Eigen::NumTraits<real_scalar_type>::dummy_precision();
 
    real_scalar_type tolerance_{default_tolerance};
 
    const matrix_type* coeff_{nullptr};
};
 
}  // namespace num_collect::linear