34 #include <visp3/core/vpConfig.h>
36 #include <visp3/core/vpColVector.h>
37 #include <visp3/core/vpMath.h>
38 #include <visp3/core/vpMatrix.h>
40 #ifdef VISP_HAVE_EIGEN3
44 #ifdef VISP_HAVE_LAPACK
46 #include <gsl/gsl_linalg.h>
47 #include <gsl/gsl_permutation.h>
51 typedef MKL_INT integer;
53 #ifdef VISP_HAVE_LAPACK_BUILT_IN
54 typedef long int integer;
58 extern "C" int dgetrf_(integer *m, integer *n,
double *a, integer *lda, integer *ipiv, integer *info);
59 extern "C" void dgetri_(integer *n,
double *a, integer *lda, integer *ipiv,
double *work, integer *lwork,
64 #if defined(VISP_HAVE_OPENCV)
65 #include <opencv2/core/core.hpp>
69 #include <visp3/core/vpException.h>
70 #include <visp3/core/vpMatrixException.h>
136 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
147 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
152 inv[1][1] = (*this)[0][0] * d;
153 inv[0][0] = (*this)[1][1] * d;
154 inv[0][1] = -(*this)[0][1] * d;
155 inv[1][0] = -(*this)[1][0] * d;
162 if (std::fabs(d) < std::numeric_limits<double>::epsilon()) {
167 const unsigned int index_0 = 0;
168 const unsigned int index_1 = 1;
169 const unsigned int index_2 = 2;
170 inv[index_0][index_0] = (((*this)[index_1][index_1] * (*this)[index_2][index_2]) - ((*
this)[index_1][index_2] * (*this)[index_2][index_1])) * d;
171 inv[index_0][index_1] = (((*this)[index_0][index_2] * (*this)[index_2][index_1]) - ((*
this)[index_0][index_1] * (*this)[index_2][index_2])) * d;
172 inv[index_0][index_2] = (((*this)[index_0][index_1] * (*this)[index_1][index_2]) - ((*
this)[index_0][index_2] * (*this)[index_1][index_1])) * d;
174 inv[index_1][index_0] = (((*this)[index_1][index_2] * (*this)[index_2][index_0]) - ((*
this)[index_1][index_0] * (*this)[index_2][index_2])) * d;
175 inv[index_1][index_1] = (((*this)[index_0][index_0] * (*this)[index_2][index_2]) - ((*
this)[index_0][index_2] * (*this)[index_2][index_0])) * d;
176 inv[index_1][index_2] = (((*this)[index_0][index_2] * (*this)[index_1][index_0]) - ((*
this)[index_0][index_0] * (*this)[index_1][index_2])) * d;
178 inv[index_2][index_0] = (((*this)[index_1][index_0] * (*this)[index_2][index_1]) - ((*
this)[index_1][index_1] * (*this)[index_2][index_0])) * d;
179 inv[index_2][index_1] = (((*this)[index_0][index_1] * (*this)[index_2][index_0]) - ((*
this)[index_0][index_0] * (*this)[index_2][index_1])) * d;
180 inv[index_2][index_2] = (((*this)[index_0][index_0] * (*this)[index_1][index_1]) - ((*
this)[index_0][index_1] * (*this)[index_1][index_0])) * d;
184 #if defined(VISP_HAVE_LAPACK)
186 #elif defined(VISP_HAVE_EIGEN3)
188 #elif defined(VISP_HAVE_OPENCV)
192 "Install Lapack, Eigen3 or OpenCV 3rd party"));
237 return (*
this)[0][0];
240 return (((*
this)[0][0] * (*
this)[1][1]) - ((*
this)[0][1] * (*
this)[1][0]));
243 const unsigned int index_0 = 0;
244 const unsigned int index_1 = 1;
245 const unsigned int index_2 = 2;
246 return ((((*
this)[index_0][index_0] * (((*
this)[index_1][index_1] * (*
this)[index_2][index_2]) - ((*
this)[index_1][index_2] * (*
this)[index_2][index_1]))) -
247 ((*
this)[index_0][index_1] * (((*
this)[index_1][index_0] * (*
this)[index_2][index_2]) - ((*
this)[index_1][index_2] * (*
this)[index_2][index_0])))) +
248 ((*
this)[index_0][index_2] * (((*
this)[index_1][index_0] * (*
this)[index_2][index_1]) - ((*
this)[index_1][index_1] * (*
this)[index_2][index_0]))));
251 #if defined(VISP_HAVE_LAPACK)
253 #elif defined(VISP_HAVE_EIGEN3)
255 #elif defined(VISP_HAVE_OPENCV)
259 "Install Lapack, Eigen3 or OpenCV 3rd party"));
264 #if defined(VISP_HAVE_LAPACK)
301 #if defined(VISP_HAVE_GSL)
310 unsigned int tda =
static_cast<unsigned int>(A->tda);
311 for (
unsigned int i = 0; i <
rowNum; ++i) {
312 unsigned int k = i * tda;
313 for (
unsigned int j = 0; j <
colNum; ++j)
314 A->data[k + j] = (*
this)[i][j];
322 inverse.tda = inverse.size2;
323 inverse.data = Ainv.
data;
327 gsl_permutation *p = gsl_permutation_alloc(
rowNum);
331 gsl_linalg_LU_decomp(A, p, &s);
332 gsl_linalg_LU_invert(A, p, &inverse);
334 gsl_permutation_free(p);
345 integer dim = (integer)
rowNum;
348 integer lwork = dim * dim;
349 integer *ipiv =
new integer[dim + 1];
350 double *work =
new double[lwork];
354 dgetrf_(&dim, &dim, A.
data, &lda, &ipiv[1], &info);
361 dgetri_(&dim, A.
data, &dim, &ipiv[1], work, &lwork, &info);
402 #if defined(VISP_HAVE_GSL)
414 unsigned int tda = (
unsigned int)A->tda;
415 for (
unsigned int i = 0; i <
rowNum; i++) {
416 unsigned int k = i * tda;
417 for (
unsigned int j = 0; j <
colNum; j++)
418 A->data[k + j] = (*
this)[i][j];
421 gsl_permutation *p = gsl_permutation_alloc(
rowNum);
425 gsl_linalg_LU_decomp(A, p, &s);
426 det = gsl_linalg_LU_det(A, s);
428 gsl_permutation_free(p);
440 integer dim = (integer)
rowNum;
443 integer *ipiv =
new integer[dim + 1];
447 dgetrf_(&dim, &dim, A.
data, &lda, &ipiv[1], &info);
453 double det = A[0][0];
454 for (
unsigned int i = 1; i <
rowNum; ++i) {
459 for (
int i = 1; i <= dim; ++i) {
474 #if defined(VISP_HAVE_OPENCV)
517 cv::Mat Minv = M.inv(cv::DECOMP_LU);
520 memcpy(A.
data, Minv.data, (
size_t)(8 * Minv.rows * Minv.cols));
564 det = cv::determinant(M);
570 #if defined(VISP_HAVE_EIGEN3)
613 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->
data, this->
getRows(),
615 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > A_(A.
data, this->getRows(),
659 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->
data, this->
getRows(),
662 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