NumCpp  2.9.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
43 {
44  namespace linalg
45  {
46  //============================================================================
47  // Method Description:
55  template<typename dtype>
57  {
59 
60  const Shape inShape = inArray.shape();
61  if (inShape.rows != inShape.cols)
62  {
63  THROW_INVALID_ARGUMENT_ERROR("input array must be square.");
64  }
65 
66  NdArray<double> inArrayDouble = inArray.template astype<double>();
67  NdArray<int> incidence = nc::zeros<int>(inShape);
68 
69  for (uint32 k = 0; k < inShape.rows - 1; ++k)
70  {
71  if (utils::essentiallyEqual(inArrayDouble(k, k), 0.))
72  {
73  uint32 l = k;
74  while (l < inShape.cols && utils::essentiallyEqual(inArrayDouble(k, l), 0.))
75  {
76  ++l;
77  }
78 
79  inArrayDouble.swapRows(k, l);
80  incidence(k, k) = 1;
81  incidence(k, l) = 1;
82  }
83  }
84 
85  NdArray<double> result(inShape);
86 
87  for (uint32 k = 0; k < inShape.rows; ++k)
88  {
89  result(k, k) = -1. / inArrayDouble(k, k);
90  for (uint32 i = 0; i < inShape.rows; ++i)
91  {
92  for (uint32 j = 0; j < inShape.cols; ++j)
93  {
94  if ((i - k) && (j - k))
95  {
96  result(i, j) =
97  inArrayDouble(i, j) + inArrayDouble(k, j) * inArrayDouble(i, k) * result(k, k);
98  }
99  else if ((i - k) && !(j - k))
100  {
101  result(i, k) = inArrayDouble(i, k) * result(k, k);
102  }
103  else if (!(i - k) && (j - k))
104  {
105  result(k, j) = inArrayDouble(k, j) * result(k, k);
106  }
107  }
108  }
109 
110  inArrayDouble = result;
111  }
112 
113  result *= -1.;
114 
115  for (int i = static_cast<int>(inShape.rows) - 1; i >= 0; --i)
116  {
117  if (incidence(i, i) != 1)
118  {
119  continue;
120  }
121 
122  int k = 0;
123  for (; k < static_cast<int>(inShape.cols); ++k)
124  {
125  if ((k - i) && incidence(i, k) != 0)
126  {
127  result.swapCols(i, k);
128  break;
129  }
130  }
131  }
132 
133  return result;
134  }
135  } // namespace linalg
136 } // namespace nc
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:36
#define STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype)
Definition: StaticAsserts.hpp:49
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4092
void swapCols(int32 colIdx1, int32 colIdx2) noexcept
Definition: NdArrayCore.hpp:4251
void swapRows(int32 rowIdx1, int32 rowIdx2) noexcept
Definition: NdArrayCore.hpp:4266
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: Constants.hpp:45
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:56
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:51
Definition: Coordinate.hpp:45
std::uint32_t uint32
Definition: Types.hpp:40