operator_conjugate_gradient.h

Go to the documentation of this file.

/*
 * Copyright 2024 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 <limits>
 
#include <Eigen/Core>
 
#include "num_collect/base/concepts/dense_vector.h"
#include "num_collect/base/concepts/invocable.h"
#include "num_collect/base/index_type.h"
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/logging/logging_mixin.h"
 
namespace num_collect::linear::impl {
 
constexpr auto operator_conjugate_gradient_tag = logging::log_tag_view(
    "num_collect::linear::impl::operator_conjugate_gradient");
 
template <base::concepts::dense_vector Vector>
class operator_conjugate_gradient : public logging::logging_mixin {
public:
    using vector_type = Vector;
 
    using scalar_type = typename vector_type::Scalar;
 
    operator_conjugate_gradient()
        : logging::logging_mixin(operator_conjugate_gradient_tag) {}
 
    template <base::concepts::invocable<const vector_type&, vector_type&>
            CoeffFunction>
    void solve(CoeffFunction&& coeff_function, const vector_type& rhs,
        vector_type& solution) {
        iterations_ = 0;
 
        coeff_function(solution, residual_);
        residual_ = rhs - residual_;
        const scalar_type tolerance =
            tolerance_rate_ * tolerance_rate_ * rhs.squaredNorm();
        NUM_COLLECT_LOG_TRACE(this->logger(), "tolerance={}", tolerance);
        scalar_type residual_norm = residual_.squaredNorm();
        if (residual_norm <= tolerance) {
            NUM_COLLECT_LOG_SUMMARY(this->logger(),
                "No iteration needed. residual_norm={}", residual_norm);
            return;
        }
        p_ = residual_;
        while (iterations_ < max_iterations_) {
            coeff_function(p_, coeff_p_);
            const scalar_type alpha = residual_.squaredNorm() /
                (p_.dot(coeff_p_) +
                    std::numeric_limits<scalar_type>::epsilon());
            solution += alpha * p_;
            previous_residual_ = residual_;
            residual_ -= alpha * coeff_p_;
            residual_norm = residual_.squaredNorm();
            if (residual_norm <= tolerance) {
                break;
            }
            const scalar_type beta = residual_norm /
                (previous_residual_.squaredNorm() +
                    std::numeric_limits<scalar_type>::epsilon());
            p_ = residual_ + beta * p_;
            ++iterations_;
        }
        NUM_COLLECT_LOG_SUMMARY(this->logger(),
            "Finished iterations: iterations={}, residual_norm={}", iterations_,
            residual_norm);
    }
 
    auto tolerance_rate(const scalar_type& val)
        -> operator_conjugate_gradient& {
        tolerance_rate_ = val;
        return *this;
    }
 
    auto iterations() -> index_type { return iterations_; }
 
private:
    index_type iterations_{};
 
    static constexpr index_type default_max_iterations = 1000;
 
    index_type max_iterations_{default_max_iterations};
 
    static constexpr scalar_type default_tolerance_rate =
        Eigen::NumTraits<scalar_type>::dummy_precision();
 
    scalar_type tolerance_rate_{default_tolerance_rate};
 
    vector_type residual_{};
 
    vector_type previous_residual_{};
    vector_type p_{};
    vector_type coeff_p_{};
};
 
}  // namespace num_collect::linear::impl