jacobi.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef cMHN_NONTT_JACOBI_H
#define cMHN_NONTT_JACOBI_H
 
#include <cmhn/nontt/mhn_operator.hpp>
 
#include <prc.hpp>
 
namespace cMHN::nonTT
{
    template<pRC::Operator::Transform OT = pRC::Operator::Transform::None,
        class T, class X, pRC::Size D, class RX = pRC::RemoveConstReference<X>>
    X jacobi(nonTT::MHNOperator<T, D> const &op, X const &b)
    {
        // starting value for pTheta
        auto pTheta = eval(pRC::unit<X>(T(1) / pow(T(D), T(2))));
 
        // get Tensor with diagonal entries of (1-Q)
        auto dg = apply_diag(op, eval(pRC::unit<X>()));
 
        for(pRC::Index i = 0; i <= D; ++i)
        {
            pTheta = b + apply_offdiag<OT>(op, pTheta);
            pTheta /= dg;
        }
 
        return pTheta;
    }
}
 
#endif // cMHN_NONTT_JACOBI_H