symmetric_successive_over_relaxation.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 <cmath>
 
#include <Eigen/Core>
 
#include "num_collect/base/concepts/dense_vector_of.h"
#include "num_collect/base/concepts/real_scalar.h"
#include "num_collect/base/concepts/sparse_matrix.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/precondition.h"
#include "num_collect/linear/iterative_solver_base.h"
#include "num_collect/logging/logging_macros.h"
 
namespace num_collect::linear {
 
template <base::concepts::sparse_matrix Matrix>
class symmetric_successive_over_relaxation;
 
namespace impl {
 
template <base::concepts::sparse_matrix Matrix>
struct iterative_solver_traits<symmetric_successive_over_relaxation<Matrix>> {
    using matrix_type = Matrix;
};
 
}  // namespace impl
 
template <base::concepts::sparse_matrix Matrix>
class symmetric_successive_over_relaxation
    : public iterative_solver_base<
          symmetric_successive_over_relaxation<Matrix>> {
    static_assert(Matrix::IsRowMajor == 1, "Row major matrix is required.");
    static_assert(base::concepts::real_scalar<typename Matrix::Scalar>,
        "Complex matrices are not supported.");
 
public:
    using base_type =
        iterative_solver_base<symmetric_successive_over_relaxation<Matrix>>;
 
    using typename base_type::matrix_type;
    using typename base_type::real_scalar_type;
    using typename base_type::scalar_type;
    using typename base_type::storage_index_type;
 
protected:
    using base_type::coeff;
 
public:
    using vector_type = Eigen::VectorX<scalar_type>;
 
    symmetric_successive_over_relaxation() = default;
 
    void compute(const matrix_type& coeff) {
        base_type::compute(coeff);
        diag_ = coeff.diagonal();
        inv_diag_ = diag_.cwiseInverse();
        intermidiate_solution_.resize(coeff.cols());
        NUM_COLLECT_PRECONDITION(inv_diag_.array().isFinite().all(),
            "All diagonal elements of the coefficient matrix must not be "
            "zero.");
    }
 
    template <base::concepts::dense_vector_of<scalar_type> Right,
        base::concepts::dense_vector_of<scalar_type> Solution>
    void solve_vector_in_place(const Right& right, Solution& solution) const {
        const auto& coeff_ref = coeff();
 
        NUM_COLLECT_PRECONDITION(coeff_ref.rows() == coeff_ref.cols(),
            "Coefficient matrix must be a square matrix.");
        NUM_COLLECT_PRECONDITION(right.rows() == coeff_ref.cols(),
            "Right-hand-side vector must have the number of elements same as "
            "the size of the coefficient matrix.");
        NUM_COLLECT_PRECONDITION(solution.rows() == coeff_ref.cols(),
            "Solution vector must have the number of elements same as the size "
            "of the coefficient matrix.");
 
        iterations_ = 0;
        const scalar_type right_norm = right.squaredNorm();
        const index_type max_iterations = base_type::max_iterations();
        while (iterations_ < max_iterations) {
            iterate(coeff_ref, right, solution);
            if (!std::isfinite(residual_)) {
                NUM_COLLECT_LOG_AND_THROW(algorithm_failure,
                    "Failure in "
                    "symmetric_successive_over_relaxation.");
            }
            ++iterations_;
            using std::sqrt;
            residual_rate_ = sqrt(residual_ / right_norm);
            if (residual_rate_ < base_type::tolerance()) {
                break;
            }
        }
    }
 
    [[nodiscard]] auto iterations() const noexcept -> index_type {
        return iterations_;
    }
 
    [[nodiscard]] auto residual_rate() const noexcept -> index_type {
        return residual_rate_;
    }
 
    auto relaxation_coeff(const scalar_type& val)
        -> symmetric_successive_over_relaxation& {
        NUM_COLLECT_PRECONDITION(static_cast<scalar_type>(0) < val &&
                val < static_cast<scalar_type>(2),
            "Relaxation coefficient must be in the range (0, 2).");
        relaxation_coeff_ = val;
        return *this;
    }
 
private:
    template <base::concepts::dense_vector_of<scalar_type> Right,
        base::concepts::dense_vector_of<scalar_type> Solution>
    void iterate(const matrix_type& coeff_ref, const Right& right,
        Solution& solution) const {
        const index_type size = coeff_ref.rows();
        const scalar_type prev_sol_coeff =
            static_cast<scalar_type>(1) - relaxation_coeff_;
        residual_ = static_cast<scalar_type>(0);
 
        // Forward update.
        for (index_type i = 0; i < size; ++i) {
            scalar_type numerator = right(i);
            for (typename matrix_type::InnerIterator iter(coeff_ref, i); iter;
                ++iter) {
                if (iter.index() < i) {
                    numerator -=
                        iter.value() * intermidiate_solution_(iter.index());
                } else if (iter.index() > i) {
                    numerator -= iter.value() * solution(iter.index());
                }
            }
            const scalar_type row_residual = numerator - diag_(i) * solution(i);
            intermidiate_solution_(i) =
                relaxation_coeff_ * numerator * inv_diag_(i) +
                prev_sol_coeff * solution(i);
            residual_ += row_residual * row_residual;
        }
 
        // Backward update.
        for (index_type i = size - 1; i >= 0; --i) {
            scalar_type numerator = right(i);
            for (typename matrix_type::InnerIterator iter(coeff_ref, i); iter;
                ++iter) {
                if (iter.index() < i) {
                    numerator -=
                        iter.value() * intermidiate_solution_(iter.index());
                } else if (iter.index() > i) {
                    numerator -= iter.value() * solution(iter.index());
                }
            }
            solution(i) = relaxation_coeff_ * numerator * inv_diag_(i) +
                prev_sol_coeff * intermidiate_solution_(i);
        }
    }
 
    mutable index_type iterations_{};
 
    mutable scalar_type residual_{};
 
    mutable scalar_type residual_rate_{};
 
    scalar_type relaxation_coeff_{static_cast<scalar_type>(1)};
 
    vector_type diag_{};
 
    vector_type inv_diag_{};
 
    mutable vector_type intermidiate_solution_{};
};
 
}  // namespace num_collect::linear