| 48 | // 'sign'. |
| 49 | template <class Scalar> |
| 50 | static typename Eigen::NumTraits<Scalar>::Real SLogDet( |
| 51 | const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>& inputs, |
| 52 | Scalar* sign) { |
| 53 | using RealScalar = typename Eigen::NumTraits<Scalar>::Real; |
| 54 | RealScalar log_abs_det = 0; |
| 55 | *sign = 1; |
| 56 | // An empty matrix' determinant is defined to be 1. |
| 57 | // (https://en.wikipedia.org/wiki/Determinant) |
| 58 | if (inputs.size() > 0) { |
| 59 | // Compute the log determinant through a Partially Pivoted LU decomposition |
| 60 | using Eigen::Dynamic; |
| 61 | Eigen::PartialPivLU<Eigen::Matrix<Scalar, Dynamic, Dynamic>> lu(inputs); |
| 62 | Eigen::Matrix<Scalar, Dynamic, Dynamic> LU = lu.matrixLU(); |
| 63 | *sign = lu.permutationP().determinant(); |
| 64 | auto diag = LU.diagonal().array().eval(); |
| 65 | auto abs_diag = diag.cwiseAbs().eval(); |
| 66 | log_abs_det += abs_diag.log().sum(); |
| 67 | *sign *= (diag / abs_diag).prod(); |
| 68 | } |
| 69 | if (!Eigen::numext::isfinite(log_abs_det)) { |
| 70 | *sign = 0; |
| 71 | log_abs_det = |
| 72 | log_abs_det > 0 ? -std::log(RealScalar(0)) : std::log(RealScalar(0)); |
| 73 | } |
| 74 | return log_abs_det; |
| 75 | } |
| 76 | |
| 77 | template <class Scalar> |
| 78 | class LogDeterminantOp : public LinearAlgebraOp<Scalar> { |