lu.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_LU_H
#define pRC_ALGORITHMS_LU_H
 
#include <prc/core/tensor/declarations.hpp>
#include <prc/core/tensor/operator/hint.hpp>
 
namespace pRC
{
    template<Operator::Hint H = Operator::Hint::None, class X,
        IsTensorish R = RemoveReference<X>>
        requires(R::Dimension == 2) && (R::size(0) == R::size(1))
    static inline constexpr auto lu(X &&a)
    {
        constexpr auto N = R::size(0);
 
        decltype(auto) u = copy<!(!IsReference<X> && !IsConst<R>)>(eval(a));
 
        if constexpr(H != Operator::Hint::UpperTriangular)
        {
            if(u(0, 0) != pRC::zero())
            {
                for(Index j = 0; j < N; ++j)
                {
                    for(Index i = 0; i <= j; ++i)
                    {
                        for(Index k = 0; k < i; ++k)
                        {
                            u(i, j) -= innerProduct(u(i, k), u(k, j));
                        }
                    }
                    for(Index i = j + 1; i < N; ++i)
                    {
                        for(Index k = 0; k < j; ++k)
                        {
                            u(i, j) -= innerProduct(u(i, k), u(k, j));
                        }
                        if(u(i, j) != pRC::zero())
                        {
                            u(i, j) /= u(j, j);
                        }
                    }
                }
            }
            else
            {
                Logging::error(
                    "LU decomposition failed: Input is has 0 at entry (0,0)");
            }
        }
        else
        {
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                if(!isLowerTriangular(a))
                {
                    Logging::error(
                        "LU decomposition failed: Input is not lower "
                        "triangular part of matrix.");
                }
            }
        }
 
        return u;
    }
}
#endif // pRC_ALGORITHMS_LU_H