NumCpp  2.11.0
A Templatized Header Only C++ Implementation of the Python NumPy Library
inv.hpp
Go to the documentation of this file.
1 #pragma once
29 
30 #include <algorithm>
31 #include <string>
32 
35 #include "NumCpp/Core/Shape.hpp"
36 #include "NumCpp/Core/Types.hpp"
38 #include "NumCpp/Linalg/det.hpp"
39 #include "NumCpp/NdArray.hpp"
41 
42 namespace nc::linalg
43 {
44  //============================================================================
45  // Method Description:
53  template<typename dtype>
55  {
57 
58  const Shape inShape = inArray.shape();
59  if (inShape.rows != inShape.cols)
60  {
61  THROW_INVALID_ARGUMENT_ERROR("input array must be square.");
62  }
63 
64  NdArray<double> inArrayDouble = inArray.template astype<double>();
65  NdArray<int> incidence = nc::zeros<int>(inShape);
66 
67  for (uint32 k = 0; k < inShape.rows - 1; ++k)
68  {
69  if (utils::essentiallyEqual(inArrayDouble(k, k), 0.))
70  {
71  uint32 l = k;
72  while (l < inShape.cols && utils::essentiallyEqual(inArrayDouble(k, l), 0.))
73  {
74  ++l;
75  }
76 
77  inArrayDouble.swapRows(k, l);
78  incidence(k, k) = 1;
79  incidence(k, l) = 1;
80  }
81  }
82 
83  NdArray<double> result(inShape);
84 
85  for (uint32 k = 0; k < inShape.rows; ++k)
86  {
87  result(k, k) = -1. / inArrayDouble(k, k);
88  for (uint32 i = 0; i < inShape.rows; ++i)
89  {
90  for (uint32 j = 0; j < inShape.cols; ++j)
91  {
92  if ((i - k) && (j - k))
93  {
94  result(i, j) = inArrayDouble(i, j) + inArrayDouble(k, j) * inArrayDouble(i, k) * result(k, k);
95  }
96  else if ((i - k) && !(j - k))
97  {
98  result(i, k) = inArrayDouble(i, k) * result(k, k);
99  }
100  else if (!(i - k) && (j - k))
101  {
102  result(k, j) = inArrayDouble(k, j) * result(k, k);
103  }
104  }
105  }
106 
107  inArrayDouble = result;
108  }
109 
110  result *= -1.;
111 
112  for (int i = static_cast<int>(inShape.rows) - 1; i >= 0; --i)
113  {
114  if (incidence(i, i) != 1)
115  {
116  continue;
117  }
118 
119  int k = 0;
120  for (; k < static_cast<int>(inShape.cols); ++k)
121  {
122  if ((k - i) && incidence(i, k) != 0)
123  {
124  result.swapCols(i, k);
125  break;
126  }
127  }
128  }
129 
130  return result;
131  }
132 } // namespace nc::linalg
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:37
#define STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype)
Definition: StaticAsserts.hpp:56
self_type & swapCols(index_type colIdx1, index_type colIdx2) noexcept
Definition: NdArrayCore.hpp:4625
self_type & swapRows(index_type rowIdx1, index_type rowIdx2) noexcept
Definition: NdArrayCore.hpp:4643
const Shape & shape() const noexcept
Definition: NdArrayCore.hpp:4464
A Shape Class for NdArrays.
Definition: Core/Shape.hpp:41
uint32 rows
Definition: Core/Shape.hpp:44
uint32 cols
Definition: Core/Shape.hpp:45
constexpr auto j
Definition: Core/Constants.hpp:42
Definition: cholesky.hpp:41
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:54
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:48
std::uint32_t uint32
Definition: Types.hpp:40