36#include <visp3/core/vpConfig.h>
38#include <visp3/core/vpColVector.h>
39#include <visp3/core/vpMath.h>
40#include <visp3/core/vpMatrix.h>
42#ifdef VISP_HAVE_EIGEN3
46#ifdef VISP_HAVE_LAPACK
48#include <gsl/gsl_linalg.h>
49#include <gsl/gsl_permutation.h>
53typedef MKL_INT integer;
55#ifdef VISP_HAVE_LAPACK_BUILT_IN
56typedef long int integer;
60extern "C" int dgetrf_(integer *m, integer *n,
double *a, integer *lda, integer *ipiv, integer *info);
61extern "C" void dgetri_(integer *n,
double *a, integer *lda, integer *ipiv,
double *work, integer *lwork,
66#if defined(VISP_HAVE_OPENCV)
67#include <opencv2/core/core.hpp>
71#include <visp3/core/vpException.h>
72#include <visp3/core/vpMatrixException.h>
133 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
143 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
148 inv[1][1] = (*this)[0][0] * d;
149 inv[0][0] = (*this)[1][1] * d;
150 inv[0][1] = -(*this)[0][1] * d;
151 inv[1][0] = -(*this)[1][0] * d;
157 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
162 inv[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1]) * d;
163 inv[0][1] = ((*this)[0][2] * (*this)[2][1] - (*this)[0][1] * (*this)[2][2]) * d;
164 inv[0][2] = ((*this)[0][1] * (*this)[1][2] - (*this)[0][2] * (*this)[1][1]) * d;
166 inv[1][0] = ((*this)[1][2] * (*this)[2][0] - (*this)[1][0] * (*this)[2][2]) * d;
167 inv[1][1] = ((*this)[0][0] * (*this)[2][2] - (*this)[0][2] * (*this)[2][0]) * d;
168 inv[1][2] = ((*this)[0][2] * (*this)[1][0] - (*this)[0][0] * (*this)[1][2]) * d;
170 inv[2][0] = ((*this)[1][0] * (*this)[2][1] - (*this)[1][1] * (*this)[2][0]) * d;
171 inv[2][1] = ((*this)[0][1] * (*this)[2][0] - (*this)[0][0] * (*this)[2][1]) * d;
172 inv[2][2] = ((*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]) * d;
175#if defined(VISP_HAVE_LAPACK)
177#elif defined(VISP_HAVE_EIGEN3)
179#elif defined(VISP_HAVE_OPENCV)
183 "Install Lapack, Eigen3 or OpenCV 3rd party"));
224 return (*
this)[0][0];
226 return ((*
this)[0][0] * (*
this)[1][1] - (*
this)[0][1] * (*
this)[1][0]);
228 return ((*
this)[0][0] * ((*
this)[1][1] * (*
this)[2][2] - (*
this)[1][2] * (*
this)[2][1]) -
229 (*
this)[0][1] * ((*
this)[1][0] * (*
this)[2][2] - (*
this)[1][2] * (*
this)[2][0]) +
230 (*
this)[0][2] * ((*
this)[1][0] * (*
this)[2][1] - (*
this)[1][1] * (*
this)[2][0]));
232#if defined(VISP_HAVE_LAPACK)
234#elif defined(VISP_HAVE_EIGEN3)
236#elif defined(VISP_HAVE_OPENCV)
240 "Install Lapack, Eigen3 or OpenCV 3rd party"));
245#if defined(VISP_HAVE_LAPACK)
278#if defined(VISP_HAVE_GSL)
287 unsigned int tda = (
unsigned int)A->tda;
288 for (
unsigned int i = 0; i <
rowNum; i++) {
289 unsigned int k = i * tda;
290 for (
unsigned int j = 0; j <
colNum; j++)
291 A->data[k + j] = (*
this)[i][j];
299 inverse.tda = inverse.size2;
300 inverse.data = Ainv.
data;
304 gsl_permutation *p = gsl_permutation_alloc(
rowNum);
308 gsl_linalg_LU_decomp(A, p, &s);
309 gsl_linalg_LU_invert(A, p, &inverse);
311 gsl_permutation_free(p);
322 integer dim = (integer)
rowNum;
325 integer lwork = dim * dim;
326 integer *ipiv =
new integer[dim + 1];
327 double *work =
new double[lwork];
331 dgetrf_(&dim, &dim, A.
data, &lda, &ipiv[1], &info);
338 dgetri_(&dim, A.
data, &dim, &ipiv[1], work, &lwork, &info);
375#if defined(VISP_HAVE_GSL)
387 unsigned int tda = (
unsigned int)A->tda;
388 for (
unsigned int i = 0; i <
rowNum; i++) {
389 unsigned int k = i * tda;
390 for (
unsigned int j = 0; j <
colNum; j++)
391 A->data[k + j] = (*
this)[i][j];
394 gsl_permutation *p = gsl_permutation_alloc(
rowNum);
398 gsl_linalg_LU_decomp(A, p, &s);
399 det = gsl_linalg_LU_det(A, s);
401 gsl_permutation_free(p);
413 integer dim = (integer)
rowNum;
416 integer *ipiv =
new integer[dim + 1];
420 dgetrf_(&dim, &dim, A.
data, &lda, &ipiv[1], &info);
426 double det = A[0][0];
427 for (
unsigned int i = 1; i <
rowNum; i++) {
432 for (
int i = 1; i <= dim; i++) {
447#if defined(VISP_HAVE_OPENCV)
486 cv::Mat Minv = M.inv(cv::DECOMP_LU);
489 memcpy(A.
data, Minv.data, (
size_t)(8 * Minv.rows * Minv.cols));
529 det = cv::determinant(M);
535#if defined(VISP_HAVE_EIGEN3)
574 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->
data, this->
getRows(),
576 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > A_(A.
data, this->getRows(),
616 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->
data, this->
getRows(),
619 return M.determinant();
unsigned int getCols() const
Type * data
Address of the first element of the data array.
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
unsigned int getRows() const
error that can be emitted by ViSP classes.
Implementation of a matrix and operations on matrices.
vpMatrix inverseByLU() const
vpMatrix inverseByLUEigen3() const
vpMatrix inverseByLUOpenCV() const
double detByLUEigen3() const
double detByLUOpenCV() const
vpMatrix inverseByLULapack() const
double det(vpDetMethod method=LU_DECOMPOSITION) const
double detByLULapack() const