101 requires(is_matrix_or_view_v<M>)
103 using value_t = get_value_t<M>;
104 static_assert(std::is_convertible_v<value_t,
double>
or std::is_convertible_v<value_t, std::complex<
double>>,
105 "Error in nda::determinant_in_place: Value type needs to be convertible to double or std::complex<double>");
106 static_assert(
not std::is_const_v<M>,
"Error in nda::determinant_in_place: Value type cannot be const");
109 if (m.empty())
return value_t{1};
112 if (m.extent(0) != m.extent(1))
NDA_RUNTIME_ERROR <<
"Error in nda::determinant_in_place: Matrix is not square: " << m.shape();
115 const int dim = m.extent(0);
117 int info = lapack::getrf(m, ipiv);
118 if (info < 0)
NDA_RUNTIME_ERROR <<
"Error in nda::determinant_in_place: info = " << info;
121 auto det = value_t{1};
123 for (
int i = 0; i < dim; i++) {
126 if (ipiv(i) != i + 1) ++n_flips;
129 return ((n_flips % 2 == 1) ? -det : det);
212 auto b00 = +m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1);
213 auto b10 = -m(1, 0) * m(2, 2) + m(1, 2) * m(2, 0);
214 auto b20 = +m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0);
215 auto b01 = -m(0, 1) * m(2, 2) + m(0, 2) * m(2, 1);
216 auto b11 = +m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0);
217 auto b21 = -m(0, 0) * m(2, 1) + m(0, 1) * m(2, 0);
218 auto b02 = +m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1);
219 auto b12 = -m(0, 0) * m(1, 2) + m(0, 2) * m(1, 0);
220 auto b22 = +m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
223 auto det = m(0, 0) * b00 + m(0, 1) * b10 + m(0, 2) * b20;
224 if (det == 0.0)
NDA_RUNTIME_ERROR <<
"Error in nda::inverse3_in_place: Matrix is not invertible";
225 auto detinv = 1.0 / det;
228 m(0, 0) = detinv * b00;
229 m(0, 1) = detinv * b01;
230 m(0, 2) = detinv * b02;
231 m(1, 0) = detinv * b10;
232 m(1, 1) = detinv * b11;
233 m(1, 2) = detinv * b12;
234 m(2, 0) = detinv * b20;
235 m(2, 1) = detinv * b21;
236 m(2, 2) = detinv * b22;