qr.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_QR_H
#define pRC_ALGORITHMS_QR_H
 
#include <prc/config.hpp>
#include <prc/core/basic/functions/ceil_div.hpp>
#include <prc/core/basic/functions/copy.hpp>
#include <prc/core/basic/functions/min.hpp>
#include <prc/core/basic/range.hpp>
#include <prc/core/tensor/functions/chip.hpp>
#include <prc/core/tensor/functions/norm.hpp>
#include <prc/core/tensor/functions/slice.hpp>
#include <prc/core/tensor/operator/functions/is_unitary.hpp>
#include <prc/core/tensor/type_traits.hpp>
#include <prc/core/value/functions/is_approx.hpp>
#include <prc/core/value/limits.hpp>
 
namespace pRC
{
    template<Size B = 32, class X, class R = RemoveReference<X>,
        If<IsTensorish<R>> = 0, If<IsFloat<typename R::Value>> = 0,
        If<IsSatisfied<(typename R::Dimension() == 2)>> = 0>
    static inline constexpr auto qr(X &&input)
    {
        using T = typename R::Type;
        constexpr auto M = R::size(0);
        constexpr auto N = R::size(1);
        constexpr auto D = min(M, N);
 
        decltype(auto) a = copy<!(!IsReference<X>() && !IsConst<R>() &&
            cDebugLevel < DebugLevel::High)>(eval(input));
 
        Tensor<T, M, M> h = identity();
 
        range<Context::CompileTime, ceilDiv(D, B)>(
            [&](auto const b)
            {
                constexpr auto s = b * B;
                range<s, min(D, B + s)>(
                    [&](auto const k)
                    {
                        auto bA = slice<M - s, N - s>(a, s, s);
                        Tensor v = chip<1>(bA, k - s);
                        for(Index i = 0; i <= k - s; ++i)
                        {
                            v(i) = zero();
                        }
 
                        auto const tau = [&a, &v, &k]()
                        {
                            auto const c = a(k, k);
                            auto const sqNorm = norm<2, 1>(v)();
                            if(max(sqNorm, norm<2, 1>(imag(c))) <=
                                NumericLimits<typename T::Value>::min())
                            {
                                v = zero();
                                return zero<T>();
                            }
                            else
                            {
                                auto beta = sqrt(norm<2, 1>(c) + sqNorm);
                                if(real(c) < zero<typename T::NonComplex>())
                                {
                                    beta = -beta;
                                }
                                v /= c + beta;
                                return conj((beta + c) / beta);
                            }
                        }();
 
                        v(k - s) = unit();
                        Tensor const vDa = conj(v) * bA;
                        bA -= tau * tensorProduct(v, vDa);
 
                        auto bH = slice<M - s, M>(h, s, 0);
                        Tensor const vDh = conj(v) * bH;
                        bH -= tau * tensorProduct(v, vDh);
                    });
            });
 
        auto const result = [&a, &h]()
        {
            using Tp = tuple<Tensor<T, M, D>, Tensor<T, D, N>>;
            if constexpr(M > N)
            {
                return Tp(slice<M, N>(adjoint(h), 0, 0), slice<N, N>(a, 0, 0));
            }
            else
            {
                return Tp(adjoint(h), move(a));
            }
        }();
 
        if constexpr(cDebugLevel >= DebugLevel::High)
        {
            auto const &q = get<0>(result);
            auto const &r = get<1>(result);
 
            if(!isApprox(input, q * r))
            {
                Logging::error("QR decomposition failed.");
            }
 
            if(!isUnitary(q))
            {
                Logging::error("QR decomposition failed: Q is not unitary.");
            }
        }
 
        return result;
    }
}
#endif // pRC_ALGORITHMS_QR_H