NumCpp  2.1.0
A C++ implementation of the Python Numpy library
lu_decomposition.hpp
Go to the documentation of this file.
1 #pragma once
34 
37 #include "NumCpp/Core/Types.hpp"
39 #include "NumCpp/NdArray.hpp"
41 
42 #include <cmath>
43 #include <utility>
44 
45 namespace nc
46 {
47  namespace linalg
48  {
49  //============================================================================
50  // Method Description:
57  template<typename dtype>
58  std::pair<NdArray<double>, NdArray<double> > lu_decomposition(const NdArray<dtype>& inMatrix)
59  {
61 
62  const auto shape = inMatrix.shape();
63  if(!shape.issquare())
64  {
65  THROW_RUNTIME_ERROR("Input matrix should be square.");
66  }
67 
68  NdArray<double> lMatrix = zeros_like<double>(inMatrix);
69  NdArray<double> uMatrix = inMatrix.template astype<double>();
70 
71  for(uint32 col = 0; col < shape.cols; ++col)
72  {
73  lMatrix(col, col) = 1;
74 
75  for(uint32 row = col + 1; row < shape.rows; ++row)
76  {
77  const double& divisor = uMatrix(col, col);
78  if (utils::essentiallyEqual(divisor, double{0.0}))
79  {
80  THROW_RUNTIME_ERROR("Division by 0.");
81  }
82 
83  lMatrix(row, col) = uMatrix(row, col) / divisor;
84 
85  for(uint32 col2 = col; col2 < shape.cols; ++col2)
86  {
87  uMatrix(row, col2) -= lMatrix(row, col) * uMatrix(col, col2);
88  }
89  }
90  }
91 
92  return std::make_pair(lMatrix, uMatrix);
93  }
94  } // namespace linalg
95 } // namespace nc
StaticAsserts.hpp
nc::NdArray::shape
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4312
Error.hpp
STATIC_ASSERT_ARITHMETIC
#define STATIC_ASSERT_ARITHMETIC(dtype)
Definition: StaticAsserts.hpp:38
nc::shape
Shape shape(const NdArray< dtype > &inArray) noexcept
Definition: Functions/Shape.hpp:45
nc::utils::essentiallyEqual
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:53
nc::linalg::lu_decomposition
std::pair< NdArray< double >, NdArray< double > > lu_decomposition(const NdArray< dtype > &inMatrix)
Definition: lu_decomposition.hpp:58
nc::NdArray< double >
nc::Shape::issquare
constexpr bool issquare() const noexcept
Definition: Core/Shape.hpp:124
nc::uint32
std::uint32_t uint32
Definition: Types.hpp:41
zeros_like.hpp
NdArray.hpp
nc::Shape::cols
uint32 cols
Definition: Core/Shape.hpp:46
THROW_RUNTIME_ERROR
#define THROW_RUNTIME_ERROR(msg)
Definition: Error.hpp:38
nc
Definition: Coordinate.hpp:45
nc::Shape::rows
uint32 rows
Definition: Core/Shape.hpp:45
essentiallyEqual.hpp
Types.hpp