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>
53 typedef MKL_INT integer;
55 #ifdef VISP_HAVE_LAPACK_BUILT_IN
56 typedef long int integer;
60 extern "C" int dgetrf_(integer *m, integer *n,
double *a, integer *lda, integer *ipiv, integer *info);
61 extern "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 rowNum
Number of rows in the array.
unsigned int getRows() const
unsigned int colNum
Number of columns in the array.
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