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()) {
144 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
149 inv[1][1] = (*this)[0][0] * d;
150 inv[0][0] = (*this)[1][1] * d;
151 inv[0][1] = -(*this)[0][1] * d;
152 inv[1][0] = -(*this)[1][0] * d;
159 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
164 inv[0][0] = (((*this)[1][1] * (*this)[2][2]) - ((*
this)[1][2] * (*this)[2][1])) * d;
165 inv[0][1] = (((*this)[0][2] * (*this)[2][1]) - ((*
this)[0][1] * (*this)[2][2])) * d;
166 inv[0][2] = (((*this)[0][1] * (*this)[1][2]) - ((*
this)[0][2] * (*this)[1][1])) * d;
168 inv[1][0] = (((*this)[1][2] * (*this)[2][0]) - ((*
this)[1][0] * (*this)[2][2])) * d;
169 inv[1][1] = (((*this)[0][0] * (*this)[2][2]) - ((*
this)[0][2] * (*this)[2][0])) * d;
170 inv[1][2] = (((*this)[0][2] * (*this)[1][0]) - ((*
this)[0][0] * (*this)[1][2])) * d;
172 inv[2][0] = (((*this)[1][0] * (*this)[2][1]) - ((*
this)[1][1] * (*this)[2][0])) * d;
173 inv[2][1] = (((*this)[0][1] * (*this)[2][0]) - ((*
this)[0][0] * (*this)[2][1])) * d;
174 inv[2][2] = (((*this)[0][0] * (*this)[1][1]) - ((*
this)[0][1] * (*this)[1][0])) * d;
178 #if defined(VISP_HAVE_LAPACK)
180 #elif defined(VISP_HAVE_EIGEN3)
182 #elif defined(VISP_HAVE_OPENCV)
186 "Install Lapack, Eigen3 or OpenCV 3rd party"));
227 return (*
this)[0][0];
230 return (((*
this)[0][0] * (*
this)[1][1]) - ((*
this)[0][1] * (*
this)[1][0]));
233 return ((((*
this)[0][0] * (((*
this)[1][1] * (*
this)[2][2]) - ((*
this)[1][2] * (*
this)[2][1]))) -
234 ((*
this)[0][1] * (((*
this)[1][0] * (*
this)[2][2]) - ((*
this)[1][2] * (*
this)[2][0])))) +
235 ((*
this)[0][2] * (((*
this)[1][0] * (*
this)[2][1]) - ((*
this)[1][1] * (*
this)[2][0]))));
238 #if defined(VISP_HAVE_LAPACK)
240 #elif defined(VISP_HAVE_EIGEN3)
242 #elif defined(VISP_HAVE_OPENCV)
246 "Install Lapack, Eigen3 or OpenCV 3rd party"));
251 #if defined(VISP_HAVE_LAPACK)
284 #if defined(VISP_HAVE_GSL)
293 unsigned int tda =
static_cast<unsigned int>(A->tda);
294 for (
unsigned int i = 0; i <
rowNum; ++i) {
295 unsigned int k = i * tda;
296 for (
unsigned int j = 0; j <
colNum; ++j)
297 A->data[k + j] = (*
this)[i][j];
305 inverse.tda = inverse.size2;
306 inverse.data = Ainv.
data;
310 gsl_permutation *p = gsl_permutation_alloc(
rowNum);
314 gsl_linalg_LU_decomp(A, p, &s);
315 gsl_linalg_LU_invert(A, p, &inverse);
317 gsl_permutation_free(p);
328 integer dim = (integer)
rowNum;
331 integer lwork = dim * dim;
332 integer *ipiv =
new integer[dim + 1];
333 double *work =
new double[lwork];
337 dgetrf_(&dim, &dim, A.
data, &lda, &ipiv[1], &info);
344 dgetri_(&dim, A.
data, &dim, &ipiv[1], work, &lwork, &info);
381 #if defined(VISP_HAVE_GSL)
393 unsigned int tda = (
unsigned int)A->tda;
394 for (
unsigned int i = 0; i <
rowNum; i++) {
395 unsigned int k = i * tda;
396 for (
unsigned int j = 0; j <
colNum; j++)
397 A->data[k + j] = (*
this)[i][j];
400 gsl_permutation *p = gsl_permutation_alloc(
rowNum);
404 gsl_linalg_LU_decomp(A, p, &s);
405 det = gsl_linalg_LU_det(A, s);
407 gsl_permutation_free(p);
419 integer dim = (integer)
rowNum;
422 integer *ipiv =
new integer[dim + 1];
426 dgetrf_(&dim, &dim, A.
data, &lda, &ipiv[1], &info);
432 double det = A[0][0];
433 for (
unsigned int i = 1; i <
rowNum; ++i) {
438 for (
int i = 1; i <= dim; ++i) {
453 #if defined(VISP_HAVE_OPENCV)
492 cv::Mat Minv = M.inv(cv::DECOMP_LU);
495 memcpy(A.
data, Minv.data, (
size_t)(8 * Minv.rows * Minv.cols));
535 det = cv::determinant(M);
541 #if defined(VISP_HAVE_EIGEN3)
580 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->
data, this->
getRows(),
582 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > A_(A.
data, this->getRows(),
622 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->
data, this->
getRows(),
625 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