pivotLU_decomposition.hpp

Go to the documentation of this file.

#pragma once
 
#include <cmath>
#include <tuple>
 
#include "NumCpp/Core/Internal/Error.hpp"
#include "NumCpp/Core/Internal/StaticAsserts.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Functions/eye.hpp"
#include "NumCpp/Functions/zeros_like.hpp"
#include "NumCpp/NdArray.hpp"
#include "NumCpp/Utils/essentiallyEqual.hpp"
 
namespace nc::linalg
{
    //============================================================================
    // Method Description:
    template<typename dtype>
    std::tuple<NdArray<double>, NdArray<double>, NdArray<double>> pivotLU_decomposition(const NdArray<dtype>& inMatrix)
    {
        STATIC_ASSERT_ARITHMETIC(dtype);
 
        const auto& shape = inMatrix.shape();
 
        if (!shape.issquare())
        {
            THROW_RUNTIME_ERROR("Input matrix should be square.");
        }
 
        NdArray<double> lMatrix = zeros_like<double>(inMatrix);
        NdArray<double> uMatrix = inMatrix.template astype<double>();
        NdArray<double> pMatrix = eye<double>(shape.rows);
 
        for (uint32 k = 0; k < shape.rows; ++k)
        {
            double max = 0.;
            uint32 pk  = 0;
            for (uint32 i = k; i < shape.rows; ++i)
            {
                double s = 0.;
                for (uint32 j = k; j < shape.cols; ++j)
                {
                    s += std::fabs(uMatrix(i, j));
                }
 
                const double q = std::fabs(uMatrix(i, k)) / s;
                if (q > max)
                {
                    max = q;
                    pk  = i;
                }
            }
 
            if (utils::essentiallyEqual(max, double{ 0. }))
            {
                THROW_RUNTIME_ERROR("Division by 0.");
            }
 
            if (pk != k)
            {
                for (uint32 j = 0; j < shape.cols; ++j)
                {
                    std::swap(pMatrix(k, j), pMatrix(pk, j));
                    std::swap(lMatrix(k, j), lMatrix(pk, j));
                    std::swap(uMatrix(k, j), uMatrix(pk, j));
                }
            }
 
            for (uint32 i = k + 1; i < shape.rows; ++i)
            {
                lMatrix(i, k) = uMatrix(i, k) / uMatrix(k, k);
 
                for (uint32 j = k; j < shape.cols; ++j)
                {
                    uMatrix(i, j) = uMatrix(i, j) - lMatrix(i, k) * uMatrix(k, j);
                }
            }
        }
 
        for (uint32 k = 0; k < shape.rows; ++k)
        {
            lMatrix(k, k) = 1.;
        }
 
        return std::make_tuple(lMatrix, uMatrix, pMatrix);
    }
} // namespace nc::linalg