jacobi_rotation.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_JACOBI_ROTATION_H
#define pRC_ALGORITHMS_JACOBI_ROTATION_H
 
#include <prc/core/tensor/declarations.hpp>
 
namespace pRC
{
    template<class T>
        requires IsValue<T> || IsComplex<T>
    class JacobiRotation
    {
    public:
        template<IsConvertible<T> R1, IsConvertible<T> R2, IsConvertible<T> R3>
        static constexpr auto MakeJacobi(R1 const &x, R2 const &y, R3 const &z)
        {
            JacobiRotation<T> r;
            r.makeJacobi(x, y, z);
            return r;
        }
 
        template<IsConvertible<T> R1, IsConvertible<T> R2>
        static constexpr auto MakeGivens(R1 const &x, R2 const &y)
        {
            JacobiRotation<T> r;
            r.makeGivens(x, y);
            return r;
        }
 
    public:
        ~JacobiRotation() = default;
        constexpr JacobiRotation(JacobiRotation const &) = default;
        constexpr JacobiRotation(JacobiRotation &&) = default;
        constexpr JacobiRotation &operator=(JacobiRotation const &) & = default;
        constexpr JacobiRotation &operator=(JacobiRotation &&) & = default;
        constexpr JacobiRotation() = default;
 
        constexpr JacobiRotation(T const &c, T const &s)
            : mC(c)
            , mS(s)
        {
        }
 
        template<class X, IsTensorish R = RemoveReference<X>>
            requires IsConvertible<typename R::Type, T> && (R::Dimension == 2)
        constexpr JacobiRotation(X &&a, Index const p, Index const q)
        {
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                Tensor<T, 2, 2> m;
                m(0, 0) = a(p, p);
                m(1, 0) = a(q, p);
                m(0, 1) = a(p, q);
                m(1, 1) = a(q, q);
                if(!isSelfAdjoint(m))
                {
                    Logging::error(
                        "JacobiRotation input 2x2 matrix is not self-adjoint.");
                }
            }
 
            auto const &x = real(a(p, p));
            auto const y = mean(a(p, q), conj(a(q, p)));
            auto const &z = real(a(q, q));
 
            makeJacobi(x, y, z);
        }
 
        template<class X, IsTensorish R = RemoveReference<X>>
            requires IsConvertible<typename R::Type, T> && (R::Dimension == 1)
        constexpr JacobiRotation(X &&a, Index const p, Index const q)
        {
            auto const &x = a(p);
            auto const &y = a(q);
 
            makeGivens(x, y);
        }
 
        template<IsConvertible<T> R>
        constexpr JacobiRotation(JacobiRotation<R> const &other)
            : mC(other.c())
            , mS(other.s())
        {
        }
 
        template<IsConvertible<T> R>
        constexpr auto &operator=(JacobiRotation<R> const &rhs) &
        {
            mC = rhs.c();
            mS = rhs.s();
            return *this;
        }
 
        constexpr decltype(auto) c() &&
        {
            return move(mC);
        }
 
        constexpr decltype(auto) c() const &&
        {
            return move(mC);
        }
 
        constexpr auto &c() &
        {
            return mC;
        }
 
        constexpr auto &c() const &
        {
            return mC;
        }
 
        constexpr decltype(auto) s() &&
        {
            return move(mS);
        }
 
        constexpr decltype(auto) s() const &&
        {
            return move(mS);
        }
 
        constexpr auto &s() &
        {
            return mS;
        }
 
        constexpr auto &s() const &
        {
            return mS;
        }
 
    private:
        template<IsConvertible<T> R1, IsConvertible<T> R2, IsConvertible<T> R3>
        constexpr auto makeJacobi(R1 const &x, R2 const &y, R3 const &z)
        {
            auto denominator = identity<NonComplex<T>>(2) * abs(y);
            if(denominator <= NumericLimits<Value<T>>::min())
            {
                mC = identity();
                mS = zero();
            }
            else
            {
                auto const beta = (x - z) / denominator;
                auto const w = sqrt(norm<2, 1>(beta) + unit<NonComplex<T>>());
                auto const t = [&beta, &w]()
                {
                    if(beta > zero())
                    {
                        return rcp(beta + w);
                    }
                    else
                    {
                        return rcp(beta - w);
                    }
                }();
 
                mC = rcp(sqrt(norm<2, 1>(t) + unit<NonComplex<T>>()));
                mS = conj(y) / abs(y);
                if(t > zero())
                {
                    mS = -mS;
                }
                mS *= abs(t) * mC;
            }
        }
 
        template<IsConvertible<T> R1, IsConvertible<T> R2>
        constexpr auto makeGivens(R1 const &x, R2 const &y)
        {
            using R = Common<R1, R2>;
 
            if(x == zero() && y == zero())
            {
                mC = identity();
                mS = zero();
                return;
            }
 
#ifdef __FAST_MATH__
            auto const t = sqrt(norm<2, 1>(x) + norm<2, 1>(y));
            mC = conj(x) / t;
            mS = y / t;
#else
            if(y == zero())
            {
                mC = sign(conj(x));
                mS = zero();
                return;
            }
 
            if(abs(x) <= NumericLimits<Value<R>>::min())
            {
                mS = sign(y);
                mC = zero();
                return;
            }
 
            if constexpr(IsComplex<R>)
            {
                auto const x1 = norm<1>(x);
                auto const y1 = norm<1>(y);
 
                if(x1 >= y1)
                {
                    auto const xS = x / x1;
                    auto const x2 = norm<2, 1>(xS);
                    auto const yS = y / x1;
                    auto const y2 = norm<2, 1>(yS);
 
                    auto const u = sqrt(unit<NonComplex<T>>() + y2 / x2);
 
                    mC = rcp(u);
                    if(x.real() < zero())
                    {
                        mC = -mC;
                    }
                    mS = yS * conj(xS) * (mC / x2);
                }
                else
                {
                    auto const xS = x / y1;
                    auto const x2 = norm<2, 1>(xS);
                    auto const yS = y / y1;
                    auto const y2 = norm<2, 1>(yS);
 
                    auto const u = y1 * sqrt(x2 + y2);
 
                    mC = abs(x) / u;
                    mS = conj(x) / abs(x) * (y / u);
 
                    if(x.real() < zero())
                    {
                        mC = -mC;
                        mS = -mS;
                    }
                }
            }
            else
            {
                if(abs(x) > abs(y))
                {
                    auto const t = y / x;
                    auto const u = sign(x) * sqrt(unit<T>() + norm<2, 1>(t));
                    mC = rcp(u);
                    mS = t * mC;
                }
                else
                {
                    auto const t = x / y;
                    auto const u = sign(y) * sqrt(unit<T>() + norm<2, 1>(t));
                    mS = rcp(u);
                    mC = t * mS;
                }
            }
#endif // __FAST_MATH__
        }
 
    private:
        T mC;
        T mS;
    };
 
    template<class T, Size M, Size N>
    JacobiRotation(Tensor<T, M, N> const &a, Index const p, Index const q)
        -> JacobiRotation<T>;
 
    template<class V, class T, Size M, Size N>
    JacobiRotation(TensorViews::View<T, Sizes<M, N>, V> const &a, Index const p,
        Index const q) -> JacobiRotation<T>;
 
    template<class T, Size N>
    JacobiRotation(Tensor<T, N> const &a, Index const p, Index const q)
        -> JacobiRotation<T>;
 
    template<class V, class T, Size N>
    JacobiRotation(TensorViews::View<T, Sizes<N>, V> const &a, Index const p,
        Index const q) -> JacobiRotation<T>;
 
    template<class TA, class TB>
    static inline constexpr auto operator==(JacobiRotation<TA> const &a,
        JacobiRotation<TB> const &b)
    {
        return a.c() == b.c() && a.s() == b.s();
    }
 
    template<class TA, class TB>
    static inline constexpr auto operator!=(JacobiRotation<TA> const &a,
        JacobiRotation<TB> const &b)
    {
        return !(a == b);
    }
 
    template<class TA, class TB>
    static inline constexpr auto operator*(JacobiRotation<TA> const &a,
        JacobiRotation<TB> const &b)
    {
        return JacobiRotation(a.c() * b.c() - conj(a.s()) * b.s(),
            conj(a.c() * conj(b.s()) + conj(a.s()) * conj(b.c())));
    }
 
    template<class T>
    static inline constexpr auto transpose(JacobiRotation<T> const &a)
    {
        return JacobiRotation(a.c(), -conj(a.s()));
    }
 
    template<class T>
    static inline constexpr auto adjoint(JacobiRotation<T> const &a)
    {
        return JacobiRotation(conj(a.c()), -a.s());
    }
 
    template<class T, class X, IsTensorish R = RemoveReference<X>>
        requires(R::Dimension == 2) &&
        ((IsTensor<R> && !IsReference<X>) || IsAssignable<X>)
    static inline constexpr decltype(auto)
        apply(JacobiRotation<T> const &r, X &&m, Index const p, Index const q)
    {
        auto x = chip<0>(m, p);
        auto y = chip<0>(m, q);
 
        range<R::size(1)>(
            [&r, &x, &y](auto const i)
            {
                auto const tX = x(i);
                auto const tY = y(i);
                x(i) = r.c() * tX + conj(r.s()) * tY;
                y(i) = -r.s() * tX + conj(r.c()) * tY;
            });
 
        return forward<X>(m);
    }
 
    template<class T, class X, IsTensorish R = RemoveReference<X>>
        requires(R::Dimension == 2) &&
        ((IsTensor<R> && !IsReference<X>) || IsAssignable<X>)
    static inline constexpr decltype(auto)
        apply(X &&m, JacobiRotation<T> const &r, Index const p, Index const q)
    {
        auto x = chip<1>(m, p);
        auto y = chip<1>(m, q);
 
        range<R::size(0)>(
            [&r, &x, &y](auto const i)
            {
                auto const tX = x(i);
                auto const tY = y(i);
                x(i) = r.c() * tX - r.s() * tY;
                y(i) = conj(r.s()) * tX + conj(r.c()) * tY;
            });
 
        return forward<X>(m);
    }
 
    template<class T, class X, IsTensorish R = RemoveReference<X>>
        requires(R::Dimension == 1) &&
        ((IsTensor<R> && !IsReference<X>) || IsAssignable<X>)
    static inline constexpr decltype(auto)
        apply(JacobiRotation<T> const &r, X &&v, Index const p, Index const q)
    {
        auto const tX = v(p);
        auto const tY = v(q);
 
        v(p) = r.c() * tX + conj(r.s()) * tY;
        v(q) = -r.s() * tX + conj(r.c()) * tY;
 
        return forward<X>(v);
    }
 
    template<class T, class X, IsTensorish R = RemoveReference<X>>
        requires(R::Dimension == 1) &&
        ((IsTensor<R> && !IsReference<X>) || IsAssignable<X>)
    static inline constexpr decltype(auto)
        apply(X &&v, JacobiRotation<T> const &r, Index const p, Index const q)
    {
        auto const tX = v(p);
        auto const tY = v(q);
 
        v(p) = r.c() * tX - r.s() * tY;
        v(q) = conj(r.s()) * tX + conj(r.c()) * tY;
 
        return forward<X>(v);
    }
}
#endif // pRC_ALGORITHMS_JACOBI_ROTATION_H