51 #include <visp3/core/vpCPUFeatures.h>
52 #include <visp3/core/vpColVector.h>
53 #include <visp3/core/vpConfig.h>
54 #include <visp3/core/vpException.h>
55 #include <visp3/core/vpMath.h>
56 #include <visp3/core/vpMatrix.h>
57 #include <visp3/core/vpTranslationVector.h>
59 #if defined(VISP_HAVE_SIMDLIB)
60 #include <Simd/SimdLib.h>
63 #ifdef VISP_HAVE_LAPACK
65 #include <gsl/gsl_eigen.h>
66 #include <gsl/gsl_linalg.h>
67 #include <gsl/gsl_math.h>
68 #elif defined(VISP_HAVE_MKL)
75 #ifdef VISP_HAVE_LAPACK
77 #elif defined(VISP_HAVE_MKL)
78 typedef MKL_INT integer;
80 void vpMatrix::blas_dsyev(
char jobz,
char uplo,
unsigned int n_,
double *a_data,
unsigned int lda_,
double *w_data,
81 double *work_data,
int lwork_,
int &info_)
83 MKL_INT n =
static_cast<MKL_INT
>(n_);
84 MKL_INT lda =
static_cast<MKL_INT
>(lda_);
85 MKL_INT lwork =
static_cast<MKL_INT
>(lwork_);
86 MKL_INT info =
static_cast<MKL_INT
>(info_);
88 dsyev(&jobz, &uplo, &n, a_data, &lda, w_data, work_data, &lwork, &info);
92 #if defined(VISP_HAVE_LAPACK_BUILT_IN)
93 typedef long int integer;
97 extern "C" integer dsyev_(
char *jobz,
char *uplo, integer *n,
double *a, integer *lda,
double *w,
double *WORK,
98 integer *lwork, integer *info);
100 void vpMatrix::blas_dsyev(
char jobz,
char uplo,
unsigned int n_,
double *a_data,
unsigned int lda_,
double *w_data,
101 double *work_data,
int lwork_,
int &info_)
103 integer n =
static_cast<integer
>(n_);
104 integer lda =
static_cast<integer
>(lda_);
105 integer lwork =
static_cast<integer
>(lwork_);
106 integer info =
static_cast<integer
>(info_);
108 dsyev_(&jobz, &uplo, &n, a_data, &lda, w_data, work_data, &lwork, &info);
110 lwork_ =
static_cast<int>(lwork);
111 info_ =
static_cast<int>(info);
116 #if !defined(VISP_USE_MSVC) || (defined(VISP_USE_MSVC) && !defined(VISP_BUILD_SHARED_LIBS))
117 const unsigned int vpMatrix::m_lapack_min_size_default = 0;
118 unsigned int vpMatrix::m_lapack_min_size = vpMatrix::m_lapack_min_size_default;
131 if (((r + nrows) > M.
rowNum) || ((c + ncols) > M.
colNum)) {
133 "Cannot construct a sub matrix (%dx%d) starting at "
134 "position (%d,%d) that is not contained in the "
135 "original matrix (%dx%d)",
139 init(M, r, c, nrows, ncols);
205 #if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11)
276 vpMatrix::vpMatrix(
unsigned int nrows,
unsigned int ncols,
const std::initializer_list<double> &list)
359 unsigned int rnrows = r + nrows;
360 unsigned int cncols = c + ncols;
370 resize(nrows, ncols,
false,
false);
372 if (this->
rowPtrs ==
nullptr) {
375 for (
unsigned int i = 0; i < nrows; ++i) {
376 memcpy((*
this)[i], &M[i + r][c], ncols *
sizeof(
double));
426 unsigned int rnrows = r + nrows;
427 unsigned int cncols = c + ncols;
439 M.
resize(nrows, ncols,
false,
false);
440 for (
unsigned int i = 0; i < nrows; ++i) {
441 memcpy(M[i], &(*
this)[i + r][c], ncols *
sizeof(
double));
499 for (
unsigned int i = 0; i < column_size; ++i) {
500 c[i] = (*this)[i_begin + i][j];
641 if ((r.
data !=
nullptr) && (
data !=
nullptr)) {
642 memcpy(r.
data, (*
this)[i] + j_begin, row_size *
sizeof(
double));
689 unsigned int min_size = std::min<unsigned int>(
rowNum,
colNum);
693 diag.resize(min_size,
false);
695 for (
unsigned int i = 0; i < min_size; ++i) {
696 diag[i] = (*this)[i][i];
780 unsigned int nca = A.
getCols();
781 unsigned int ncb = B.
getCols();
790 if ((B.
getRows() == 0) || ((nca + ncb) == 0)) {
791 std::cerr <<
"B.getRows() == 0 || nca+ncb == 0" << std::endl;
824 int vpMatrix::print(std::ostream &s,
unsigned int length,
const std::string &intro)
const
826 typedef std::string::size_type size_type;
831 std::vector<std::string> values(m * n);
832 std::ostringstream oss;
833 std::ostringstream ossFixed;
834 std::ios_base::fmtflags original_flags = oss.flags();
836 ossFixed.setf(std::ios::fixed, std::ios::floatfield);
838 size_type maxBefore = 0;
839 size_type maxAfter = 0;
841 for (
unsigned int i = 0; i < m; ++i) {
842 for (
unsigned int j = 0; j < n; ++j) {
844 oss << (*this)[i][j];
845 if (oss.str().find(
"e") != std::string::npos) {
847 ossFixed << (*this)[i][j];
848 oss.str(ossFixed.str());
851 values[(i * n) + j] = oss.str();
852 size_type thislen = values[(i * n) + j].
size();
853 size_type p = values[(i * n) + j].find(
'.');
855 if (p == std::string::npos) {
866 size_type totalLength = length;
870 maxAfter = std::min<size_type>(maxAfter, totalLength - maxBefore);
872 if (!intro.empty()) {
875 s <<
"[" << m <<
"," << n <<
"]=\n";
877 for (
unsigned int i = 0; i < m; ++i) {
879 for (
unsigned int j = 0; j < n; ++j) {
880 size_type p = values[(i * n) + j].find(
'.');
881 s.setf(std::ios::right, std::ios::adjustfield);
882 s.width(
static_cast<std::streamsize
>(maxBefore));
883 s << values[(i * n) + j].substr(0, p).c_str();
886 s.setf(std::ios::left, std::ios::adjustfield);
887 if (p != std::string::npos) {
888 s.width(
static_cast<std::streamsize
>(maxAfter));
889 s << values[(i * n) + j].substr(p, maxAfter).c_str();
892 s.width(
static_cast<std::streamsize
>(maxAfter));
902 s.flags(original_flags);
904 return static_cast<int>(maxBefore + maxAfter);
949 unsigned int this_rows = this->
getRows();
950 unsigned int this_col = this->
getCols();
952 for (
unsigned int i = 0; i < this_rows; ++i) {
953 for (
unsigned int j = 0; j < this_col; ++j) {
954 os << (*this)[i][j] <<
", ";
956 if (this->
getRows() != (i + 1)) {
957 os <<
";" << std::endl;
960 os <<
"]" << std::endl;
1000 unsigned int this_rows = this->
getRows();
1001 unsigned int this_col = this->
getCols();
1002 os <<
"([ " << std::endl;
1003 for (
unsigned int i = 0; i < this_rows; ++i) {
1005 for (
unsigned int j = 0; j < this_col; ++j) {
1006 os << (*this)[i][j] <<
", ";
1008 os <<
"]," << std::endl;
1010 os <<
"])" << std::endl;
1047 unsigned int this_rows = this->
getRows();
1048 unsigned int this_col = this->
getCols();
1049 for (
unsigned int i = 0; i < this_rows; ++i) {
1050 for (
unsigned int j = 0; j < this_col; ++j) {
1051 os << (*this)[i][j];
1052 if (!(j == (this->
getCols() - 1))) {
1102 os <<
"vpMatrix " << matrixName <<
" (" << this->
getRows() <<
", " << this->
getCols() <<
"); " << std::endl;
1104 unsigned int this_rows = this->
getRows();
1105 unsigned int this_col = this->
getCols();
1106 for (
unsigned int i = 0; i < this_rows; ++i) {
1107 for (
unsigned int j = 0; j < this_col; ++j) {
1109 os << matrixName <<
"[" << i <<
"][" << j <<
"] = " << (*this)[i][j] <<
"; " << std::endl;
1112 for (
unsigned int k = 0; k <
sizeof(double); ++k) {
1113 os <<
"((unsigned char*)&(" << matrixName <<
"[" << i <<
"][" << j <<
"]) )[" << k <<
"] = 0x" << std::hex
1114 <<
static_cast<unsigned int>(((
unsigned char *)&((*
this)[i][j]))[k]) <<
"; " << std::endl;
1140 unsigned int a_rows = A.
getRows();
1141 for (
unsigned int i = r; i < (r + a_rows); ++i) {
1203 for (
unsigned int i = 0; i <
rowNum; ++i) {
1204 for (
unsigned int j = 0; j <
rowNum; ++j) {
1206 if (std::fabs(At_A[i][j]) > std::numeric_limits<double>::epsilon()) {
1212 #if defined(VISP_HAVE_LAPACK)
1213 #if defined(VISP_HAVE_GSL)
1215 gsl_vector *eval = gsl_vector_alloc(
rowNum);
1218 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(
rowNum);
1221 unsigned int Atda = (
unsigned int)m->tda;
1222 for (
unsigned int i = 0; i <
rowNum; i++) {
1223 unsigned int k = i * Atda;
1224 for (
unsigned int j = 0; j <
colNum; j++)
1225 m->data[k + j] = (*
this)[i][j];
1227 gsl_eigen_symmv(m, eval, evec, w);
1229 gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);
1231 for (
unsigned int i = 0; i <
rowNum; i++) {
1232 evalue[i] = gsl_vector_get(eval, i);
1235 gsl_eigen_symmv_free(w);
1236 gsl_vector_free(eval);
1238 gsl_matrix_free(evec);
1242 const char jobz =
'N';
1243 const char uplo =
'U';
1250 lwork =
static_cast<int>(wkopt);
1258 "You should install Lapack 3rd party"));
1326 for (
unsigned int i = 0; i <
rowNum; ++i) {
1327 for (
unsigned int j = 0; j <
rowNum; ++j) {
1329 if (std::fabs(At_A[i][j]) > std::numeric_limits<double>::epsilon()) {
1339 #if defined(VISP_HAVE_LAPACK)
1340 #if defined(VISP_HAVE_GSL)
1342 gsl_vector *eval = gsl_vector_alloc(
rowNum);
1345 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(
rowNum);
1348 unsigned int Atda = (
unsigned int)m->tda;
1349 for (
unsigned int i = 0; i <
rowNum; i++) {
1350 unsigned int k = i * Atda;
1351 for (
unsigned int j = 0; j <
colNum; j++)
1352 m->data[k + j] = (*
this)[i][j];
1354 gsl_eigen_symmv(m, eval, evec, w);
1356 gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);
1358 for (
unsigned int i = 0; i <
rowNum; i++) {
1359 evalue[i] = gsl_vector_get(eval, i);
1361 Atda = (
unsigned int)evec->tda;
1362 for (
unsigned int i = 0; i <
rowNum; i++) {
1363 unsigned int k = i * Atda;
1364 for (
unsigned int j = 0; j <
rowNum; j++) {
1365 evector[i][j] = evec->
data[k + j];
1369 gsl_eigen_symmv_free(w);
1370 gsl_vector_free(eval);
1372 gsl_matrix_free(evec);
1376 const char jobz =
'V';
1377 const char uplo =
'U';
1384 lwork =
static_cast<int>(wkopt);
1393 "You should install Lapack 3rd party"));
1418 unsigned int nbline =
getRows();
1419 unsigned int nbcol =
getCols();
1424 V.
resize(nbcol, nbcol,
false);
1430 if (nbline < nbcol) {
1431 U.
resize(nbcol, nbcol,
true);
1434 U.
resize(nbline, nbcol,
false);
1443 for (
unsigned int i = 0; i < nbcol; ++i) {
1444 if (sv[i] > maxsv) {
1449 unsigned int rank = 0;
1450 for (
unsigned int i = 0; i < nbcol; ++i) {
1451 if (sv[i] >(maxsv * svThreshold)) {
1456 kerAt.
resize(nbcol - rank, nbcol);
1457 if (rank != nbcol) {
1459 for (
unsigned int j = 0; j < nbcol; ++j) {
1461 if ((sv[j] <= (maxsv * svThreshold)) &&
1462 (std::fabs(V.
getCol(j).
sumSquare()) > std::numeric_limits<double>::epsilon())) {
1463 unsigned int v_rows = V.
getRows();
1464 for (
unsigned int i = 0; i < v_rows; ++i) {
1465 kerAt[k][i] = V[i][j];
1493 unsigned int nbrow =
getRows();
1494 unsigned int nbcol =
getCols();
1499 V.
resize(nbcol, nbcol,
false);
1505 if (nbrow < nbcol) {
1506 U.
resize(nbcol, nbcol,
true);
1509 U.
resize(nbrow, nbcol,
false);
1517 double maxsv = sv[0];
1519 unsigned int rank = 0;
1520 for (
unsigned int i = 0; i < nbcol; ++i) {
1521 if (sv[i] >(maxsv * svThreshold)) {
1526 kerA.
resize(nbcol, nbcol - rank);
1527 if (rank != nbcol) {
1529 for (
unsigned int j = 0; j < nbcol; ++j) {
1531 if (sv[j] <= (maxsv * svThreshold)) {
1532 for (
unsigned int i = 0; i < nbcol; ++i) {
1533 kerA[i][k] = V[i][j];
1540 return (nbcol - rank);
1561 unsigned int nbrow =
getRows();
1562 unsigned int nbcol =
getCols();
1563 unsigned int dim_ =
static_cast<unsigned int>(dim);
1568 V.
resize(nbcol, nbcol,
false);
1574 if (nbrow < nbcol) {
1575 U.
resize(nbcol, nbcol,
true);
1578 U.
resize(nbrow, nbcol,
false);
1585 kerA.
resize(nbcol, dim_);
1587 unsigned int rank = nbcol - dim_;
1588 for (
unsigned int k = 0; k < dim_; ++k) {
1589 unsigned int j = k + rank;
1590 for (
unsigned int i = 0; i < nbcol; ++i) {
1591 kerA[i][k] = V[i][j];
1596 double maxsv = sv[0];
1597 unsigned int rank = 0;
1598 for (
unsigned int i = 0; i < nbcol; ++i) {
1599 if (sv[i] >(maxsv * 1e-6)) {
1603 return (nbcol - rank);
1654 #if defined(VISP_HAVE_EIGEN3)
1656 #elif defined(VISP_HAVE_LAPACK)
1658 #elif defined(VISP_HAVE_OPENCV)
1662 "Install Lapack, Eigen3 or OpenCV 3rd party"));
1680 #ifdef VISP_HAVE_GSL
1682 double *b =
new double[size_];
1683 for (
size_t i = 0; i < size_; i++)
1685 gsl_matrix_view m = gsl_matrix_view_array(this->
data,
rowNum, colNum);
1686 gsl_matrix_view em = gsl_matrix_view_array(b,
rowNum,
colNum);
1687 gsl_linalg_exponential_ss(&m.matrix, &em.matrix, 0);
1691 memcpy(expA.
data, b, size_ *
sizeof(
double));
1712 for (
unsigned int i = 0; i <
rowNum; ++i) {
1714 for (
unsigned int j = 0; j <
colNum; ++j) {
1715 sum += fabs((*
this)[i][j]);
1717 if ((
sum > nA) || (i == 0)) {
1726 double sca = 1.0 / pow(2.0, s);
1729 v_expE = (c * exp) + v_eye;
1730 v_expD = (-c * exp) + v_eye;
1731 for (
int k = 2; k <= q; ++k) {
1732 c = (c * (
static_cast<double>((q - k) + 1))) / (
static_cast<double>(k * (((2 * q) - k) + 1)));
1733 v_expcX = exp * v_expX;
1735 v_expcX = c * v_expX;
1736 v_expE = v_expE + v_expcX;
1738 v_expD = v_expD + v_expcX;
1741 v_expD = v_expD - v_expcX;
1746 exp = v_expX * v_expE;
1747 for (
int k = 1; k <= s; ++k) {
1780 unsigned int m_rows = M.
getRows();
1781 unsigned int m_col = M.
getCols();
1782 for (
unsigned int i = 0; i < col; ++i) {
1783 for (
unsigned int j = 0; j < row; ++j) {
1784 M_comp[i][j] = M[i][j];
1786 for (
unsigned int j = row + 1; j < m_rows; ++j) {
1787 M_comp[i][j - 1] = M[i][j];
1790 for (
unsigned int i = col + 1; i < m_col; ++i) {
1791 for (
unsigned int j = 0; j < row; ++j) {
1792 M_comp[i - 1][j] = M[i][j];
1794 for (
unsigned int j = row + 1; j < m_rows; ++j) {
1795 M_comp[i - 1][j - 1] = M[i][j];
1812 unsigned int nbline =
getRows();
1813 unsigned int nbcol =
getCols();
1818 V.
resize(nbcol, nbcol,
false);
1824 if (nbline < nbcol) {
1825 U.
resize(nbcol, nbcol,
true);
1828 U.
resize(nbline, nbcol,
false);
1837 for (
unsigned int i = 0; i < nbcol; ++i) {
1838 if (sv[i] > maxsv) {
1844 unsigned int rank = 0;
1845 for (
unsigned int i = 0; i < nbcol; ++i) {
1846 if (sv[i] >(maxsv * svThreshold)) {
1852 double minsv = maxsv;
1853 for (
unsigned int i = 0; i < rank; ++i) {
1854 if (sv[i] < minsv) {
1859 if (std::fabs(minsv) > std::numeric_limits<double>::epsilon()) {
1860 return maxsv / minsv;
1863 return std::numeric_limits<double>::infinity();
1881 unsigned int h_col = H.
getCols();
1882 for (
unsigned int i = 0; i < h_col; ++i) {
1883 HLM[i][i] += alpha * H[i][i];
1897 for (
unsigned int i = 0; i <
dsize; ++i) {
1898 double x = *(
data + i);
1915 if (this->
dsize != 0) {
1924 unsigned int maxRank = std::min<unsigned int>(this->
getCols(), this->
getRows());
1927 unsigned int boundary = std::min<unsigned int>(maxRank, w.
size());
1932 for (
unsigned int i = 0; i < boundary; ++i) {
1957 for (
unsigned int i = 0; i <
rowNum; ++i) {
1959 for (
unsigned int j = 0; j <
colNum; ++j) {
1960 x += fabs(*(*(
rowPtrs + i) + j));
1977 double sum_square = 0.0;
1980 for (
unsigned int i = 0; i <
rowNum; ++i) {
1981 for (
unsigned int j = 0; j <
colNum; ++j) {
1983 sum_square += x * x;
1989 #if defined(VISP_BUILD_DEPRECATED_FUNCTIONS)
2023 for (
unsigned int j = 0; j <
getCols(); ++j) {
2024 c[j] = (*this)[i - 1][j];
2050 for (
unsigned int i = 0; i <
getRows(); ++i) {
2051 c[i] = (*this)[i][j - 1];
2064 for (
unsigned int i = 0; i <
rowNum; ++i) {
2065 for (
unsigned int j = 0; j <
colNum; ++j) {
2067 (*this)[i][j] = val;
Implementation of a generic 2D array used as base class for matrices and vectors.
unsigned int getCols() const
double * data
Address of the first element of the data array.
double ** rowPtrs
Address of the first element of each rows.
void insert(const vpArray2D< Type > &A, unsigned int r, unsigned int c)
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 dsize
Current array size (rowNum * colNum)
unsigned int size() const
Return the number of elements of the 2D array.
unsigned int getRows() const
unsigned int colNum
Number of columns in the array.
Implementation of column vector and the associated operations.
void resize(unsigned int i, bool flagNullify=true)
error that can be emitted by ViSP classes.
@ functionNotImplementedError
Function not implemented.
@ dimensionError
Bad dimension.
Implementation of an homogeneous matrix and operations on such kind of matrices.
static Type maximum(const Type &a, const Type &b)
Implementation of a matrix and operations on matrices.
vpColVector eigenValues() const
std::ostream & cppPrint(std::ostream &os, const std::string &matrixName="A", bool octet=false) const
double cond(double svThreshold=1e-6) const
VP_DEPRECATED void setIdentity(const double &val=1.0)
int print(std::ostream &s, unsigned int length, const std::string &intro="") const
std::ostream & maplePrint(std::ostream &os) const
vpMatrix cholesky() const
unsigned int kernel(vpMatrix &kerAt, double svThreshold=1e-6) const
vpMatrix inverseByLU() const
static vpMatrix juxtaposeMatrices(const vpMatrix &A, const vpMatrix &B)
VP_DEPRECATED double euclideanNorm() const
VP_DEPRECATED vpColVector column(unsigned int j)
void svd(vpColVector &w, vpMatrix &V)
void diag(const double &val=1.0)
vpRowVector getRow(unsigned int i) const
vpMatrix choleskyByOpenCV() const
Compute the Cholesky decomposition of a Hermitian positive-definite matrix using OpenCV library.
vpMatrix choleskyByEigen3() const
Compute the Cholesky decomposition of a Hermitian positive-definite matrix using Eigen3 library.
double inducedL2Norm() const
double infinityNorm() const
double frobeniusNorm() const
vpColVector getDiag() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
VP_DEPRECATED vpRowVector row(unsigned int i)
vpColVector getCol(unsigned int j) const
vpMatrix choleskyByLapack() const
Compute the Cholesky decomposition of a Hermitian positive-definite matrix using Lapack library.
double det(vpDetMethod method=LU_DECOMPOSITION) const
void insert(const vpMatrix &A, unsigned int r, unsigned int c)
VP_DEPRECATED void init()
std::ostream & csvPrint(std::ostream &os) const
unsigned int nullSpace(vpMatrix &kerA, double svThreshold=1e-6) const
std::ostream & matlabPrint(std::ostream &os) const
vpMatrix extract(unsigned int r, unsigned int c, unsigned int nrows, unsigned int ncols) const
Implementation of a rotation matrix and operations on such kind of matrices.
Implementation of row vector and the associated operations.
Class that consider the case of a translation vector.