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