42 #include <visp3/core/vpConfig.h> 44 #include <visp3/core/vpColVector.h> 45 #include <visp3/core/vpMath.h> 46 #include <visp3/core/vpMatrix.h> 49 #include <visp3/core/vpException.h> 50 #include <visp3/core/vpMatrixException.h> 53 #include <visp3/core/vpDebug.h> 55 #ifdef VISP_HAVE_LAPACK 56 #ifdef VISP_HAVE_LAPACK_BUILT_IN 57 typedef long int integer;
62 extern "C" int dgeqrf_(integer *m, integer *n,
double *a, integer *lda,
double *tau,
double *work, integer *lwork,
64 extern "C" int dormqr_(
char *side,
char *trans, integer *m, integer *n, integer *k,
double *a, integer *lda,
65 double *tau,
double *c__, integer *ldc,
double *work, integer *lwork, integer *info);
66 extern "C" int dorgqr_(integer *, integer *, integer *,
double *, integer *,
double *,
double *, integer *, integer *);
67 extern "C" int dtrtri_(
char *uplo,
char *diag, integer *n,
double *a, integer *lda, integer *info);
68 extern "C" int dgeqp3_(integer *m, integer *n,
double*a, integer *lda, integer *p,
69 double *tau,
double *work, integer* lwork, integer *info);
71 int allocate_work(
double **work);
73 int allocate_work(
double **work)
75 unsigned int dimWork = (
unsigned int)((*work)[0]);
77 *work =
new double[dimWork];
82 #ifdef VISP_HAVE_LAPACK 119 integer rowNum_ = (integer)this->
getRows();
120 integer colNum_ = (integer)this->
getCols();
121 integer lda = (integer)rowNum_;
122 integer dimTau = (std::min)(rowNum_, colNum_);
123 integer dimWork = -1;
124 double *tau =
new double[dimTau];
125 double *work =
new double[1];
152 std::cout <<
"dgeqrf_:Preparation:" << -info <<
"th element had an illegal value" << std::endl;
155 dimWork = allocate_work(&work);
177 std::cout <<
"dgeqrf_:" << -info <<
"th element had an illegal value" << std::endl;
186 dtrtri_((
char *)
"U", (
char *)
"N", &dimTau, C.
data, &lda, &info);
189 std::cout <<
"dtrtri_:" << -info <<
"th element had an illegal value" << std::endl;
191 std::cout <<
"dtrtri_:R(" << info <<
"," << info <<
")" 192 <<
" is exactly zero. The triangular matrix is singular " 193 "and its inverse can not be computed." 195 std::cout <<
"R=" << std::endl << C << std::endl;
204 for (
unsigned int i = 0; i < C.
getRows(); i++)
205 for (
unsigned int j = 0; j < C.
getRows(); j++)
216 dormqr_((
char *)
"R", (
char *)
"T", &rowNum_, &colNum_, &dimTau, A.
data, &lda, tau, C.
data, &ldc, work, &dimWork,
219 std::cout <<
"dormqr_:Preparation" << -info <<
"th element had an illegal value" << std::endl;
222 dimWork = allocate_work(&work);
224 dormqr_((
char *)
"R", (
char *)
"T", &rowNum_, &colNum_, &dimTau, A.
data, &lda, tau, C.
data, &ldc, work, &dimWork,
228 std::cout <<
"dormqr_:" << -info <<
"th element had an illegal value" << std::endl;
275 #ifdef VISP_HAVE_LAPACK 339 #ifdef VISP_HAVE_LAPACK 340 integer m = (integer)
rowNum;
341 integer n = (integer)
colNum;
342 integer r = std::min(n,m);
354 Q.
resize(static_cast<unsigned int>(m), static_cast<unsigned int>(q));
356 R.
resize(static_cast<unsigned int>(r), static_cast<unsigned int>(r));
358 R.
resize(static_cast<unsigned int>(r), static_cast<unsigned int>(n));
362 integer dimWork = -1;
363 double * qrdata =
new double[m*na];
364 double *tau =
new double[std::min(m,q)];
365 double *work =
new double[1];
369 for(integer i = 0; i < m; ++i)
371 for(integer j = 0; j < n; ++j)
372 qrdata[i+m*j] =
data[j + n*i];
373 for(integer j = n; j < na; ++j)
391 std::cout <<
"dgeqrf_:Preparation:" << -info <<
"th element had an illegal value" << std::endl;
397 dimWork = allocate_work(&work);
411 std::cout <<
"dgeqrf_:" << -info <<
"th element had an illegal value" << std::endl;
424 for(
int i=0;i<na;i++)
426 for(
int j=i;j<na;j++)
427 R[i][j] = qrdata[i+m*j];
428 if(std::abs(qrdata[i+m*i]) < tol)
434 for(
int i=0;i<na;i++)
437 R[i][j] = qrdata[i+m*j];
438 if(std::abs(qrdata[i+m*i]) < tol)
456 for(
int i = 0; i < m; ++i)
457 for(
int j = 0; j < q; ++j)
458 Q[i][j] = qrdata[i+m*j];
463 return (
unsigned int) r;
539 #ifdef VISP_HAVE_LAPACK 540 integer m = (integer)
rowNum;
541 integer n = (integer)
colNum;
542 integer r = std::min(n,m);
554 Q.
resize(static_cast<unsigned int>(m), static_cast<unsigned int>(q));
560 P.
resize(0, static_cast<unsigned int>(n));
564 R.
resize(static_cast<unsigned int>(r), static_cast<unsigned int>(n));
565 P.
resize(static_cast<unsigned int>(n), static_cast<unsigned int>(n));
570 integer dimWork = -1;
571 double* qrdata =
new double[m*na];
572 double* tau =
new double[std::min(q,m)];
573 double* work =
new double[1];
574 integer* p =
new integer[na];
575 for(
int i = 0; i < na; ++i)
581 for(integer i = 0; i < m; ++i)
583 for(integer j = 0; j < n; ++j)
584 qrdata[i+m*j] =
data[j + n*i];
585 for(integer j = n; j < na; ++j)
604 std::cout <<
"dgeqp3_:Preparation:" << -info <<
"th element had an illegal value" << std::endl;
612 dimWork = allocate_work(&work);
628 std::cout <<
"dgeqp3_:" << -info <<
" th element had an illegal value" << std::endl;
640 for(
int i = 0; i < na; ++i)
641 if(std::abs(qrdata[i+m*i]) < tol)
647 R.
resize(static_cast<unsigned int>(r), static_cast<unsigned int>(r));
650 R[i][j] = qrdata[i+m*j];
653 P.
resize(static_cast<unsigned int>(r), static_cast<unsigned int>(n));
654 for(
int i = 0; i < r; ++i)
659 R.
resize(static_cast<unsigned int>(na), static_cast<unsigned int>(n));
660 for(
int i=0;i<na;i++)
662 R[i][j] = qrdata[i+m*j];
664 P.
resize(static_cast<unsigned int>(n), static_cast<unsigned int>(n));
665 for(
int i = 0; i < n; ++i)
682 for(
int i = 0; i < m; ++i)
683 for(
int j = 0; j < q; ++j)
684 Q[i][j] = qrdata[i+m*j];
690 return (
unsigned int) r;
715 #ifdef VISP_HAVE_LAPACK 719 integer n = (integer)
rowNum;
726 dtrtri_((
char *)
"L", (
char *)
"N", &n, R.
data, &n, &info);
728 dtrtri_((
char *)
"U", (
char *)
"N", &n, R.
data, &n, &info);
732 std::cout <<
"dtrtri_:" << -info <<
"th element had an illegal value" << std::endl;
734 std::cout <<
"dtrtri_:R(" << info <<
"," << info <<
")" 735 <<
" is exactly zero. The triangular matrix is singular " 736 "and its inverse can not be computed." 795 unsigned int r =
t().
qrPivot(Q, R, P,
false,
true);
vpMatrix extract(unsigned int r, unsigned int c, unsigned int nrows, unsigned int ncols) const
Used to indicate that a value is not in the allowed range.
Implementation of a matrix and operations on matrices.
vpMatrix inverseByQR() const
unsigned int qrPivot(vpMatrix &Q, vpMatrix &R, vpMatrix &P, bool full=false, bool squareR=false, double tol=1e-6) const
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
error that can be emited by ViSP classes.
unsigned int getRows() const
Type * data
Address of the first element of the data array.
vpMatrix inverseTriangular(bool upper=true) const
vpMatrix inverseByQRLapack() const
unsigned int getCols() const
unsigned int rowNum
Number of rows in the array.
unsigned int qr(vpMatrix &Q, vpMatrix &R, bool full=false, bool squareR=false, double tol=1e-6) const
void solveByQR(const vpColVector &b, vpColVector &x) const
unsigned int colNum
Number of columns in the array.
Implementation of column vector and the associated operations.
error that can be emited by the vpMatrix class and its derivates