Program Listing for File PartialPivLU.hpp

↰ Return to documentation for file (include/eigenpy/decompositions/PartialPivLU.hpp)

/*
 * Copyright 2025 INRIA
 */

#ifndef __eigenpy_decompositions_partialpivlu_hpp__
#define __eigenpy_decompositions_partialpivlu_hpp__

#include <Eigen/LU>
#include <Eigen/Core>

#include "eigenpy/eigenpy.hpp"
#include "eigenpy/utils/scalar-name.hpp"
#include "eigenpy/eigen/EigenBase.hpp"

namespace eigenpy {

template <typename _MatrixType>
struct PartialPivLUSolverVisitor : public boost::python::def_visitor<
                                       PartialPivLUSolverVisitor<_MatrixType>> {
  typedef _MatrixType MatrixType;
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
      VectorXs;
  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
                        MatrixType::Options>
      MatrixXs;
  typedef Eigen::PartialPivLU<MatrixType> Solver;

  template <class PyClass>
  void visit(PyClass &cl) const {
    cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
        .def(bp::init<Eigen::DenseIndex>(
            bp::args("self", "size"),
            "Default constructor with memory preallocation"))
        .def(bp::init<MatrixType>(
            bp::args("self", "matrix"),
            "Constructs a LU factorization from a given matrix."))

        .def(EigenBaseVisitor<Solver>())

        .def("determinant", &Solver::determinant, bp::arg("self"),
             "Returns the determinant of the matrix of which *this is the LU "
             "decomposition.")
        .def(
            "compute",
            (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
                Solver::compute,
            bp::args("self", "matrix"),
            "Computes the LU factorization of given matrix.",
            bp::return_self<>())
        .def(
            "inverse",
            +[](const Solver &self) -> MatrixType { return self.inverse(); },
            bp::arg("self"),
            "Returns the inverse of the matrix of which *this is the LU "
            "decomposition.")
        .def("matrixLU", &Solver::matrixLU, bp::arg("self"),
             "Returns the LU decomposition matrix.",
             bp::return_internal_reference<>())

        .def("permutationP", &Solver::permutationP, bp::arg("self"),
             "Returns the permutation P.",
             bp::return_value_policy<bp::copy_const_reference>())

        .def("rcond", &Solver::rcond, bp::arg("self"),
             "Returns an estimate of the reciprocal condition number of the "
             "matrix.")
        .def("reconstructedMatrix", &Solver::reconstructedMatrix,
             bp::arg("self"),
             "Returns the matrix represented by the decomposition, i.e., it "
             "returns the product: P-1LUQ-1. This function is provided for "
             "debug "
             "purpose.")
        .def("solve", &solve<VectorXs>, bp::args("self", "b"),
             "Returns the solution x of A x = b using the current "
             "decomposition of A.")
        .def("solve", &solve<MatrixXs>, bp::args("self", "B"),
             "Returns the solution X of A X = B using the current "
             "decomposition of A where B is a right hand side matrix.");
  }

  static void expose() {
    static const std::string classname =
        "FullPivLU" + scalar_name<Scalar>::shortname();
    expose(classname);
  }

  static void expose(const std::string &name) {
    bp::class_<Solver>(
        name.c_str(),
        "LU decomposition of a matrix with partial pivoting, "
        "and related features. \n\n"
        "This class represents a LU decomposition of a square "
        "invertible matrix, "
        "with partial pivoting: the matrix A is decomposed as A "
        "= PLU where L is "
        "unit-lower-triangular, U is upper-triangular, and P is "
        "a permutation matrix.\n\n"
        "Typically, partial pivoting LU decomposition is only "
        "considered numerically "
        "stable for square invertible matrices. Thus LAPACK's "
        "dgesv and dgesvx require "
        "the matrix to be square and invertible. The present "
        "class does the same. It "
        "will assert that the matrix is square, but it won't "
        "(actually it can't) check "
        "that the matrix is invertible: it is your task to "
        "check that you only use this "
        "decomposition on invertible matrices. \n\n"
        "The guaranteed safe alternative, working for all matrices, "
        "is the full pivoting LU decomposition, provided by class "
        "FullPivLU. \n\n"
        "This is not a rank-revealing LU decomposition. Many features "
        "are intentionally absent from this class, such as "
        "rank computation. If you need these features, use class "
        "FullPivLU. \n\n"
        "This LU decomposition is suitable to invert invertible "
        "matrices. It is what MatrixBase::inverse() uses in the "
        "general case. On the other hand, it is not suitable to "
        "determine whether a given matrix is invertible. \n\n"
        "The data of the LU decomposition can be directly accessed "
        "through the methods matrixLU(), permutationP().",
        bp::no_init)
        .def(IdVisitor<Solver>())
        .def(PartialPivLUSolverVisitor());
  }

 private:
  template <typename MatrixOrVector>
  static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
    return self.solve(vec);
  }
};

}  // namespace eigenpy

#endif  // ifndef __eigenpy_decompositions_partialpivlu_hpp__