00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 public:
00033
00034
00035 typedef Scalar value_type;
00036
00037 enum { static_size = RowsAtCompileTime*ColsAtCompileTime };
00038
00039
00040
00041
00042
00043 typedef Scalar* iterator;
00044 typedef const Scalar* const_iterator;
00045
00046 EIGEN_STRONG_INLINE iterator begin() { return derived().data(); }
00047 EIGEN_STRONG_INLINE iterator end() { return &(derived().data()[size()-1]); }
00048 EIGEN_STRONG_INLINE const_iterator begin() const { return derived().data(); }
00049 EIGEN_STRONG_INLINE const_iterator end() const { return &(derived().data()[size()-1]); }
00050
00051
00052
00053
00054
00055
00056
00057
00058 EIGEN_STRONG_INLINE void fill(const Scalar v) { derived().setConstant(v); }
00059
00060
00061 EIGEN_STRONG_INLINE void assign(const Scalar v) { derived().setConstant(v); }
00062
00063 EIGEN_STRONG_INLINE void assign(size_t N, const Scalar v) { derived().resize(N); derived().setConstant(v); }
00064
00065
00066 EIGEN_STRONG_INLINE size_t getRowCount() const { return rows(); }
00067
00068 EIGEN_STRONG_INLINE size_t getColCount() const { return cols(); }
00069
00070
00071 EIGEN_STRONG_INLINE void unit(const size_t nRows, const Scalar diag_vals) {
00072 if (diag_vals==1)
00073 derived().setIdentity(nRows,nRows);
00074 else {
00075 derived().setZero(nRows,nRows);
00076 derived().diagonal().setConstant(diag_vals);
00077 }
00078 }
00079
00080
00081 EIGEN_STRONG_INLINE void unit() { derived().setIdentity(); }
00082
00083 EIGEN_STRONG_INLINE void eye() { derived().setIdentity(); }
00084
00085
00086 EIGEN_STRONG_INLINE void zeros() { derived().setZero(); }
00087
00088 EIGEN_STRONG_INLINE void zeros(const size_t row,const size_t col) { derived().setZero(row,col); }
00089
00090
00091 EIGEN_STRONG_INLINE void ones(const size_t row, const size_t col) { derived().setOnes(row,col); }
00092
00093 EIGEN_STRONG_INLINE void ones() { derived().setOnes(); }
00094
00095
00096
00097
00098 EIGEN_STRONG_INLINE Scalar * get_unsafe_row(size_t row) { return &derived().coeffRef(row,0); }
00099 EIGEN_STRONG_INLINE const Scalar* get_unsafe_row(size_t row) const { return &derived().coeff(row,0); }
00100
00101
00102 EIGEN_STRONG_INLINE Scalar get_unsafe(const size_t row, const size_t col) const {
00103 #ifdef _DEBUG
00104 return derived()(row,col);
00105 #else
00106 return derived().coeff(row,col);
00107 #endif
00108 }
00109
00110 EIGEN_STRONG_INLINE Scalar& get_unsafe(const size_t row, const size_t col) {
00111 #ifdef _DEBUG
00112 return derived()(row,col);
00113 #endif
00114 return derived().coeffRef(row,col);
00115 }
00116
00117 EIGEN_STRONG_INLINE void set_unsafe(const size_t row, const size_t col, const Scalar val) {
00118 #ifdef _DEBUG
00119 derived()(row,col) = val;
00120 #endif
00121 derived().coeffRef(row,col) = val;
00122 }
00123
00124
00125 EIGEN_STRONG_INLINE void push_back(Scalar val)
00126 {
00127 const Index N = size();
00128 derived().conservativeResize(N+1);
00129 coeffRef(N) = val;
00130 }
00131
00132 EIGEN_STRONG_INLINE bool isSquare() const { return cols()==rows(); }
00133 EIGEN_STRONG_INLINE bool isSingular(const Scalar absThreshold = 0) const { return std::abs(derived().determinant())<absThreshold; }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 bool fromMatlabStringFormat(const std::string &s, bool dumpErrorMsgToStdErr = true);
00147
00148
00149
00150
00151
00152
00153 std::string inMatlabFormat(const size_t decimal_digits = 6) const;
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164 void saveToTextFile(
00165 const std::string &file,
00166 mrpt::math::TMatrixTextFileFormat fileFormat = mrpt::math::MATRIX_FORMAT_ENG,
00167 bool appendMRPTHeader = false,
00168 const std::string &userHeader = std::string()
00169 ) const;
00170
00171
00172
00173
00174
00175
00176 void loadFromTextFile(const std::string &file);
00177
00178
00179
00180 void loadFromTextFile(std::istream &_input_text_stream);
00181
00182
00183 EIGEN_STRONG_INLINE void multiplyColumnByScalar(size_t c, Scalar s) { derived().col(c)*=s; }
00184 EIGEN_STRONG_INLINE void multiplyRowByScalar(size_t r, Scalar s) { derived().row(r)*=s; }
00185
00186 EIGEN_STRONG_INLINE void swapCols(size_t i1,size_t i2) { derived().col(i1).swap( derived().col(i2) ); }
00187 EIGEN_STRONG_INLINE void swapRows(size_t i1,size_t i2) { derived().row(i1).swap( derived().row(i2) ); }
00188
00189 EIGEN_STRONG_INLINE size_t countNonZero() const { return ((*static_cast<const Derived*>(this))!= 0).count(); }
00190
00191
00192
00193
00194 EIGEN_STRONG_INLINE Scalar maximum() const
00195 {
00196 if (size()==0) throw std::runtime_error("maximum: container is empty");
00197 return derived().maxCoeff();
00198 }
00199
00200
00201
00202
00203 EIGEN_STRONG_INLINE Scalar minimum() const
00204 {
00205 if (size()==0) throw std::runtime_error("minimum: container is empty");
00206 return derived().minCoeff();
00207 }
00208
00209
00210
00211
00212 EIGEN_STRONG_INLINE void minimum_maximum(
00213 Scalar & out_min,
00214 Scalar & out_max) const
00215 {
00216 out_min = minimum();
00217 out_max = maximum();
00218 }
00219
00220
00221
00222
00223
00224 EIGEN_STRONG_INLINE Scalar maximum(size_t *maxIndex) const
00225 {
00226 if (size()==0) throw std::runtime_error("maximum: container is empty");
00227 Index idx;
00228 const Scalar m = derived().maxCoeff(&idx);
00229 if (maxIndex) *maxIndex = idx;
00230 return m;
00231 }
00232
00233
00234
00235
00236 void find_index_max_value(size_t &u,size_t &v,Scalar &valMax) const
00237 {
00238 if (cols()==0 || rows()==0) throw std::runtime_error("find_index_max_value: container is empty");
00239 Index idx1,idx2;
00240 valMax = derived().maxCoeff(&idx1,&idx2);
00241 u = idx1; v = idx2;
00242 }
00243
00244
00245
00246
00247
00248 EIGEN_STRONG_INLINE Scalar minimum(size_t *minIndex) const
00249 {
00250 if (size()==0) throw std::runtime_error("minimum: container is empty");
00251 Index idx;
00252 const Scalar m =derived().minCoeff(&idx);
00253 if (minIndex) *minIndex = idx;
00254 return m;
00255 }
00256
00257
00258
00259
00260 EIGEN_STRONG_INLINE void minimum_maximum(
00261 Scalar & out_min,
00262 Scalar & out_max,
00263 size_t *minIndex,
00264 size_t *maxIndex) const
00265 {
00266 out_min = minimum(minIndex);
00267 out_max = maximum(maxIndex);
00268 }
00269
00270
00271 EIGEN_STRONG_INLINE Scalar norm_inf() const { return lpNorm<Eigen::Infinity>(); }
00272
00273
00274 EIGEN_STRONG_INLINE Scalar squareNorm() const { return squaredNorm(); }
00275
00276
00277 EIGEN_STRONG_INLINE Scalar sumAll() const { return derived().sum(); }
00278
00279
00280
00281
00282 template<typename OtherDerived>
00283 EIGEN_STRONG_INLINE void laplacian(Eigen::MatrixBase <OtherDerived>& ret) const
00284 {
00285 if (rows()!=cols()) throw std::runtime_error("laplacian: Defined for square matrixes only");
00286 const Index N = rows();
00287 ret = -(*this);
00288 for (Index i=0;i<N;i++)
00289 {
00290 Scalar deg = 0;
00291 for (Index j=0;j<N;j++) deg+= derived().coeff(j,i);
00292 ret.coeffRef(i,i)+=deg;
00293 }
00294 }
00295
00296
00297
00298
00299
00300 EIGEN_STRONG_INLINE void setSize(size_t row, size_t col)
00301 {
00302 #ifdef _DEBUG
00303 if ((Derived::RowsAtCompileTime!=Eigen::Dynamic && Derived::RowsAtCompileTime!=int(row)) || (Derived::ColsAtCompileTime!=Eigen::Dynamic && Derived::ColsAtCompileTime!=int(col))) {
00304 std::stringstream ss; ss << "setSize: Trying to change a fixed sized matrix from " << rows() << "x" << cols() << " to " << row << "x" << col;
00305 throw std::runtime_error(ss.str());
00306 }
00307 #endif
00308 const Index oldCols = cols();
00309 const Index oldRows = rows();
00310 const int nNewCols = int(col) - int(cols());
00311 const int nNewRows = int(row) - int(rows());
00312 ::mrpt::math::detail::TAuxResizer<Eigen::MatrixBase<Derived>,SizeAtCompileTime>::internal_resize(*this,row,col);
00313 if (nNewCols>0) derived().block(0,oldCols,row,nNewCols).setZero();
00314 if (nNewRows>0) derived().block(oldRows,0,nNewRows,col).setZero();
00315 }
00316
00317
00318 template <class OUTVECT>
00319 void largestEigenvector(
00320 OUTVECT &x,
00321 Scalar resolution = Scalar(0.01),
00322 size_t maxIterations = 6,
00323 int *out_Iterations = NULL,
00324 float *out_estimatedResolution = NULL ) const
00325 {
00326
00327 size_t iter=0;
00328 const Index n = rows();
00329 x.resize(n);
00330 x.setConstant(1);
00331 Scalar dif;
00332 do
00333 {
00334 Eigen::Matrix<Scalar,Derived::RowsAtCompileTime,1> xx = (*this) * x;
00335 xx *= Scalar(1.0/xx.norm());
00336 dif = (x-xx).array().abs().sum();
00337 x = xx;
00338 iter++;
00339 } while (iter<maxIterations && dif>resolution);
00340 if (out_Iterations) *out_Iterations=static_cast<int>(iter);
00341 if (out_estimatedResolution) *out_estimatedResolution=dif;
00342 }
00343
00344
00345 MatrixBase<Derived>& operator ^= (const unsigned int pow)
00346 {
00347 if (pow==0)
00348 derived().setIdentity();
00349 else
00350 for (unsigned int i=1;i<pow;i++)
00351 derived() *= derived();
00352 return *this;
00353 }
00354
00355
00356 EIGEN_STRONG_INLINE void scalarPow(const Scalar s) { (*this)=array().pow(s); }
00357
00358
00359 EIGEN_STRONG_INLINE bool isDiagonal() const
00360 {
00361 for (Index c=0;c<cols();c++)
00362 for (Index r=0;r<rows();r++)
00363 if (r!=c && coeff(r,c)!=0)
00364 return false;
00365 return true;
00366 }
00367
00368
00369 EIGEN_STRONG_INLINE Scalar maximumDiagonal() const { return diagonal().maximum(); }
00370
00371
00372
00373 EIGEN_STRONG_INLINE double mean() const
00374 {
00375 if ( size()==0) throw std::runtime_error("mean: Empty container.");
00376 return derived().sum()/static_cast<double>(size());
00377 }
00378
00379
00380
00381
00382
00383 template <class VEC>
00384 void meanAndStd(
00385 VEC &outMeanVector,
00386 VEC &outStdVector,
00387 const bool unbiased_variance = true ) const
00388 {
00389 const double N = rows();
00390 if (N==0) throw std::runtime_error("meanAndStd: Empty container.");
00391 const double N_inv = 1.0/N;
00392 const double N_ = unbiased_variance ? (N>1 ? 1.0/(N-1) : 1.0) : 1.0/N;
00393 outMeanVector.resize(cols());
00394 outStdVector.resize(cols());
00395 for (Index i=0;i<cols();i++)
00396 {
00397 outMeanVector[i]= col(i).array().sum() * N_inv;
00398 outStdVector[i] = std::sqrt( (col(i).array()-outMeanVector[i]).square().sum() * N_ );
00399 }
00400 }
00401
00402
00403
00404
00405
00406 void meanAndStdAll(
00407 double &outMean,
00408 double &outStd,
00409 const bool unbiased_variance = true ) const
00410 {
00411 const double N = size();
00412 if (N==0) throw std::runtime_error("meanAndStdAll: Empty container.");
00413 const double N_ = unbiased_variance ? (N>1 ? 1.0/(N-1) : 1.0) : 1.0/N;
00414 outMean = derived().array().sum()/static_cast<double>(size());
00415 outStd = std::sqrt( (this->array() - outMean).square().sum()*N_);
00416 }
00417
00418
00419 template <typename MAT>
00420 EIGEN_STRONG_INLINE void insertMatrix(size_t r,size_t c, const MAT &m) { derived().block(r,c,m.rows(),m.cols())=m; }
00421
00422 template <typename MAT>
00423 EIGEN_STRONG_INLINE void insertMatrixTranspose(size_t r,size_t c, const MAT &m) { derived().block(r,c,m.cols(),m.rows())=m.adjoint(); }
00424
00425 template <typename MAT> EIGEN_STRONG_INLINE void insertRow(size_t nRow, const MAT & aRow) { this->row(nRow) = aRow; }
00426 template <typename MAT> EIGEN_STRONG_INLINE void insertCol(size_t nCol, const MAT & aCol) { this->col(nCol) = aCol; }
00427
00428
00429 EIGEN_STRONG_INLINE void removeColumns(const std::vector<size_t> &idxsToRemove)
00430 {
00431 std::vector<size_t> idxs = idxsToRemove;
00432 std::sort( idxs.begin(), idxs.end() );
00433 std::vector<size_t>::iterator itEnd = std::unique( idxs.begin(), idxs.end() );
00434 idxs.resize( itEnd - idxs.begin() );
00435
00436 unsafeRemoveColumns( idxs );
00437 }
00438
00439
00440 EIGEN_STRONG_INLINE void unsafeRemoveColumns(const std::vector<size_t> &idxs)
00441 {
00442 size_t k = 1;
00443 for (std::vector<size_t>::const_reverse_iterator it = idxs.rbegin(); it != idxs.rend(); it++, k++)
00444 {
00445 const size_t nC = cols() - *it - k;
00446 if( nC > 0 )
00447 derived().block(0,*it,rows(),nC) = derived().block(0,*it+1,rows(),nC).eval();
00448 }
00449 derived().conservativeResize(NoChange,cols()-idxs.size());
00450 }
00451
00452
00453 EIGEN_STRONG_INLINE void removeRows(const std::vector<size_t> &idxsToRemove)
00454 {
00455 std::vector<size_t> idxs = idxsToRemove;
00456 std::sort( idxs.begin(), idxs.end() );
00457 std::vector<size_t>::iterator itEnd = std::unique( idxs.begin(), idxs.end() );
00458 idxs.resize( itEnd - idxs.begin() );
00459
00460 unsafeRemoveRows( idxs );
00461 }
00462
00463
00464 EIGEN_STRONG_INLINE void unsafeRemoveRows(const std::vector<size_t> &idxs)
00465 {
00466 size_t k = 1;
00467 for (std::vector<size_t>::reverse_iterator it = idxs.rbegin(); it != idxs.rend(); it++, k++)
00468 {
00469 const size_t nR = rows() - *it - k;
00470 if( nR > 0 )
00471 derived().block(*it,0,nR,cols()) = derived().block(*it+1,0,nR,cols()).eval();
00472 }
00473 derived().conservativeResize(rows()-idxs.size(),NoChange);
00474 }
00475
00476
00477 EIGEN_STRONG_INLINE const AdjointReturnType t() const { return derived().adjoint(); }
00478
00479 EIGEN_STRONG_INLINE PlainObject inv() const { PlainObject outMat = derived().inverse().eval(); return outMat; }
00480 template <class MATRIX> EIGEN_STRONG_INLINE void inv(MATRIX &outMat) const { outMat = derived().inverse().eval(); }
00481 template <class MATRIX> EIGEN_STRONG_INLINE void inv_fast(MATRIX &outMat) const { outMat = derived().inverse().eval(); }
00482 EIGEN_STRONG_INLINE Scalar det() const { return derived().determinant(); }
00483
00484
00485
00486
00487
00488
00489
00490 EIGEN_STRONG_INLINE bool empty() const { return this->getColCount()==0 || this->getRowCount()==0; }
00491
00492
00493 template<typename OTHERMATRIX> EIGEN_STRONG_INLINE void add_Ac(const OTHERMATRIX &m,const Scalar c) { (*this)+=c*m; }
00494
00495 template<typename OTHERMATRIX> EIGEN_STRONG_INLINE void substract_Ac(const OTHERMATRIX &m,const Scalar c) { (*this) -= c*m; }
00496
00497
00498 template<typename OTHERMATRIX> EIGEN_STRONG_INLINE void substract_At(const OTHERMATRIX &m) { (*this) -= m.adjoint(); }
00499
00500
00501 template<typename OTHERMATRIX> EIGEN_STRONG_INLINE void substract_An(const OTHERMATRIX& m, const size_t n) { this->noalias() -= n * m; }
00502
00503
00504 template<typename OTHERMATRIX> EIGEN_STRONG_INLINE void add_AAt(const OTHERMATRIX &A) { this->noalias() += A; this->noalias() += A.adjoint(); }
00505
00506 \
00507 template<typename OTHERMATRIX> EIGEN_STRONG_INLINE void substract_AAt(const OTHERMATRIX &A) { this->noalias() -= A; this->noalias() -= A.adjoint(); }
00508
00509
00510 template <class MATRIX1,class MATRIX2> EIGEN_STRONG_INLINE void multiply( const MATRIX1& A, const MATRIX2 &B ) { (*this)= A*B; }
00511
00512 template <class MATRIX1,class MATRIX2>
00513 EIGEN_STRONG_INLINE void multiply_AB( const MATRIX1& A, const MATRIX2 &B ) {
00514 (*this)= A*B;
00515 }
00516
00517 template <typename MATRIX1,typename MATRIX2>
00518 EIGEN_STRONG_INLINE void multiply_AtB(const MATRIX1 &A,const MATRIX2 &B) {
00519 *this = A.adjoint() * B;
00520 }
00521
00522
00523 template<typename OTHERVECTOR1,typename OTHERVECTOR2>
00524 EIGEN_STRONG_INLINE void multiply_Ab(const OTHERVECTOR1 &vIn,OTHERVECTOR2 &vOut,bool accumToOutput = false) const {
00525 if (accumToOutput) vOut.noalias() += (*this) * vIn;
00526 else vOut = (*this) * vIn;
00527 }
00528
00529 \
00530 template<typename OTHERVECTOR1,typename OTHERVECTOR2>
00531 EIGEN_STRONG_INLINE void multiply_Atb(const OTHERVECTOR1 &vIn,OTHERVECTOR2 &vOut,bool accumToOutput = false) const {
00532 if (accumToOutput) vOut.noalias() += this->adjoint() * vIn;
00533 else vOut = this->adjoint() * vIn;
00534 }
00535
00536 template <typename MAT_C, typename MAT_R>
00537 EIGEN_STRONG_INLINE void multiply_HCHt(const MAT_C &C,MAT_R &R,bool accumResultInOutput=false) const {
00538 if (accumResultInOutput)
00539 R.noalias() += (*this) * C * this->adjoint();
00540 else R.noalias() = (*this) * C * this->adjoint();
00541 }
00542
00543 template <typename MAT_C, typename MAT_R>
00544 EIGEN_STRONG_INLINE void multiply_HtCH(const MAT_C &C,MAT_R &R,bool accumResultInOutput=false) const {
00545 if (accumResultInOutput)
00546 R.noalias() += this->adjoint() * C * (*this);
00547 else R.noalias() = this->adjoint() * C * (*this);
00548 }
00549
00550
00551 template <typename MAT_C>
00552 EIGEN_STRONG_INLINE Scalar multiply_HCHt_scalar(const MAT_C &C) const {
00553 return ( (*this) * C * this->adjoint() ).eval()(0,0);
00554 }
00555
00556
00557 template <typename MAT_C>
00558 EIGEN_STRONG_INLINE Scalar multiply_HtCH_scalar(const MAT_C &C) const {
00559 return ( this->adjoint() * C * (*this) ).eval()(0,0);
00560 }
00561
00562
00563 template<typename MAT_A>
00564 EIGEN_STRONG_INLINE void multiply_AAt_scalar(const MAT_A &A,typename MAT_A::value_type f) {
00565 *this = (A * A.adjoint()) * f;
00566 }
00567
00568
00569 template<typename MAT_A> EIGEN_STRONG_INLINE void multiply_AtA_scalar(const MAT_A &A,typename MAT_A::value_type f) {
00570 *this = (A.adjoint() * A) * f;
00571 }
00572
00573
00574 template <class MAT_A,class SKEW_3VECTOR> void multiply_A_skew3(const MAT_A &A,const SKEW_3VECTOR &v) {
00575 mrpt::math::multiply_A_skew3(A,v,*this); }
00576
00577
00578 template <class SKEW_3VECTOR,class MAT_A> void multiply_skew3_A(const SKEW_3VECTOR &v,const MAT_A &A) {
00579 mrpt::math::multiply_skew3_A(v,A,*this); }
00580
00581
00582
00583 template <class MAT_A,class MAT_OUT>
00584 EIGEN_STRONG_INLINE void multiply_subMatrix(const MAT_A &A,MAT_OUT &outResult,const size_t A_cols_offset,const size_t A_rows_offset,const size_t A_col_count) const {
00585 outResult = derived() * A.block(A_rows_offset,A_cols_offset,derived().cols(),A_col_count);
00586 }
00587
00588 template <class MAT_A,class MAT_B,class MAT_C>
00589 void multiply_ABC(const MAT_A &A, const MAT_B &B, const MAT_C &C) {
00590 *this = A*B*C;
00591 }
00592
00593 template <class MAT_A,class MAT_B,class MAT_C>
00594 void multiply_ABCt(const MAT_A &A, const MAT_B &B, const MAT_C &C) {
00595 *this = A*B*C.adjoint();
00596 }
00597
00598 template <class MAT_A,class MAT_B,class MAT_C>
00599 void multiply_AtBC(const MAT_A &A, const MAT_B &B, const MAT_C &C) {
00600 *this = A.adjoint()*B*C;
00601 }
00602
00603 template <class MAT_A,class MAT_B>
00604 EIGEN_STRONG_INLINE void multiply_ABt(const MAT_A &A,const MAT_B &B) {
00605 *this = A*B.adjoint();
00606 }
00607
00608 template <class MAT_A>
00609 EIGEN_STRONG_INLINE void multiply_AAt(const MAT_A &A) {
00610 *this = A*A.adjoint();
00611 }
00612
00613 template <class MAT_A>
00614 EIGEN_STRONG_INLINE void multiply_AtA(const MAT_A &A) {
00615 *this = A.adjoint()*A;
00616 }
00617
00618 template <class MAT_A,class MAT_B>
00619 EIGEN_STRONG_INLINE void multiply_result_is_symmetric(const MAT_A &A,const MAT_B &B) {
00620 *this = A*B;
00621 }
00622
00623
00624
00625 template<class MAT2,class MAT3 >
00626 EIGEN_STRONG_INLINE void leftDivideSquare(const MAT2 &A, MAT3 &RES) const
00627 {
00628 Eigen::ColPivHouseholderQR<PlainObject> QR = A.colPivHouseholderQr();
00629 if (!QR.isInvertible()) throw std::runtime_error("leftDivideSquare: Matrix A is not invertible");
00630 RES = QR.inverse() * (*this);
00631 }
00632
00633
00634 template<class MAT2,class MAT3 >
00635 EIGEN_STRONG_INLINE void rightDivideSquare(const MAT2 &B, MAT3 &RES) const
00636 {
00637 Eigen::ColPivHouseholderQR<PlainObject> QR = B.colPivHouseholderQr();
00638 if (!QR.isInvertible()) throw std::runtime_error("rightDivideSquare: Matrix B is not invertible");
00639 RES = (*this) * QR.inverse();
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652 template <class MATRIX1,class MATRIX2>
00653 EIGEN_STRONG_INLINE void eigenVectors( MATRIX1 & eVecs, MATRIX2 & eVals ) const;
00654
00655
00656
00657
00658
00659
00660 template <class MATRIX1,class VECTOR1>
00661 EIGEN_STRONG_INLINE void eigenVectorsVec( MATRIX1 & eVecs, VECTOR1 & eVals ) const;
00662
00663
00664
00665
00666
00667
00668 template <class VECTOR>
00669 EIGEN_STRONG_INLINE void eigenValues( VECTOR & eVals ) const
00670 {
00671 PlainObject eVecs;
00672 eVecs.resizeLike(*this);
00673 this->eigenVectorsVec(eVecs,eVals);
00674 }
00675
00676
00677
00678 template <class MATRIX1,class MATRIX2>
00679 EIGEN_STRONG_INLINE void eigenVectorsSymmetric( MATRIX1 & eVecs, MATRIX2 & eVals ) const;
00680
00681
00682
00683
00684 template <class MATRIX1,class VECTOR1>
00685 EIGEN_STRONG_INLINE void eigenVectorsSymmetricVec( MATRIX1 & eVecs, VECTOR1 & eVals ) const;
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697 template <class MATRIX> EIGEN_STRONG_INLINE bool chol(MATRIX &U) const
00698 {
00699 Eigen::LLT<PlainObject> Chol = derived().selfadjointView<Eigen::Lower>().llt();
00700 if (Chol.info()==Eigen::NoConvergence)
00701 return false;
00702 U = PlainObject(Chol.matrixU());
00703 return true;
00704 }
00705
00706
00707
00708
00709 EIGEN_STRONG_INLINE size_t rank(double threshold=0) const
00710 {
00711 Eigen::ColPivHouseholderQR<PlainObject> QR = this->colPivHouseholderQr();
00712 if (threshold>0) QR.setThreshold(threshold);
00713 return QR.rank();
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723 EIGEN_STRONG_INLINE MatrixBase<Derived>& Sqrt() { (*this) = this->array().sqrt(); return *this; }
00724 EIGEN_STRONG_INLINE PlainObject Sqrt() const { PlainObject res = this->array().sqrt(); return res; }
00725
00726 EIGEN_STRONG_INLINE MatrixBase<Derived>& Abs() { (*this) = this->array().abs(); return *this; }
00727 EIGEN_STRONG_INLINE PlainObject Abs() const { PlainObject res = this->array().abs(); return res; }
00728
00729 EIGEN_STRONG_INLINE MatrixBase<Derived>& Log() { (*this) = this->array().log(); return *this; }
00730 EIGEN_STRONG_INLINE PlainObject Log() const { PlainObject res = this->array().log(); return res; }
00731
00732 EIGEN_STRONG_INLINE MatrixBase<Derived>& Exp() { (*this) = this->array().exp(); return *this; }
00733 EIGEN_STRONG_INLINE PlainObject Exp() const { PlainObject res = this->array().exp(); return res; }
00734
00735 EIGEN_STRONG_INLINE MatrixBase<Derived>& Square() { (*this) = this->array().square(); return *this; }
00736 EIGEN_STRONG_INLINE PlainObject Square() const { PlainObject res = this->array().square(); return res; }
00737
00738
00739 void normalize(Scalar valMin, Scalar valMax)
00740 {
00741 if (size()==0) return;
00742 Scalar curMin,curMax;
00743 minimum_maximum(curMin,curMax);
00744 Scalar minMaxDelta = curMax - curMin;
00745 if (minMaxDelta==0) minMaxDelta = 1;
00746 const Scalar minMaxDelta_ = (valMax-valMin)/minMaxDelta;
00747 this->array() = (this->array()-curMin)*minMaxDelta_+valMin;
00748 }
00749
00750 inline void adjustRange(Scalar valMin, Scalar valMax) { normalize(valMin,valMax); }
00751
00752
00753
00754
00755
00756 template <class OtherDerived> EIGEN_STRONG_INLINE void extractRow(size_t nRow, Eigen::EigenBase<OtherDerived> &v, size_t startingCol = 0) const {
00757 v = derived().block(nRow,startingCol,1,cols()-startingCol);
00758 }
00759
00760 template <typename T> inline void extractRow(size_t nRow, std::vector<T> &v, size_t startingCol = 0) const {
00761 const size_t N = cols()-startingCol;
00762 v.resize(N);
00763 for (size_t i=0;i<N;i++) v[i]=(*this)(nRow,startingCol+i);
00764 }
00765
00766 template <class VECTOR> EIGEN_STRONG_INLINE void extractRowAsCol(size_t nRow, VECTOR &v, size_t startingCol = 0) const
00767 {
00768 v = derived().adjoint().block(startingCol,nRow,cols()-startingCol,1);
00769 }
00770
00771
00772
00773 template <class VECTOR> EIGEN_STRONG_INLINE void extractCol(size_t nCol, VECTOR &v, size_t startingRow = 0) const {
00774 v = derived().block(startingRow,nCol,rows()-startingRow,1);
00775 }
00776
00777 template <typename T> inline void extractCol(size_t nCol, std::vector<T> &v, size_t startingRow = 0) const {
00778 const size_t N = rows()-startingRow;
00779 v.resize(N);
00780 for (size_t i=0;i<N;i++) v[i]=(*this)(startingRow+i,nCol);
00781 }
00782
00783 template <class MATRIX> EIGEN_STRONG_INLINE void extractMatrix(const size_t firstRow, const size_t firstCol, MATRIX &m) const
00784 {
00785 m = derived().block(firstRow,firstCol,m.rows(),m.cols());
00786 }
00787 template <class MATRIX> EIGEN_STRONG_INLINE void extractMatrix(const size_t firstRow, const size_t firstCol, const size_t nRows, const size_t nCols, MATRIX &m) const
00788 {
00789 m.resize(nRows,nCols);
00790 m = derived().block(firstRow,firstCol,nRows,nCols);
00791 }
00792
00793
00794 template <class MATRIX>
00795 EIGEN_STRONG_INLINE void extractSubmatrix(const size_t row_first,const size_t row_last,const size_t col_first,const size_t col_last,MATRIX &out) const
00796 {
00797 out.resize(row_last-row_first+1,col_last-col_first+1);
00798 out = derived().block(row_first,col_first,row_last-row_first+1,col_last-col_first+1);
00799 }
00800
00801
00802
00803
00804
00805 template <class MATRIX>
00806 void extractSubmatrixSymmetricalBlocks(
00807 const size_t block_size,
00808 const std::vector<size_t> &block_indices,
00809 MATRIX& out) const
00810 {
00811 if (block_size<1) throw std::runtime_error("extractSubmatrixSymmetricalBlocks: block_size must be >=1");
00812 if (cols()!=rows()) throw std::runtime_error("extractSubmatrixSymmetricalBlocks: Matrix is not square.");
00813
00814 const size_t N = block_indices.size();
00815 const size_t nrows_out=N*block_size;
00816 out.resize(nrows_out,nrows_out);
00817 if (!N) return;
00818 for (size_t dst_row_blk=0;dst_row_blk<N; ++dst_row_blk )
00819 {
00820 for (size_t dst_col_blk=0;dst_col_blk<N; ++dst_col_blk )
00821 {
00822 #if defined(_DEBUG)
00823 if (block_indices[dst_col_blk]*block_size + block_size-1>=size_t(cols())) throw std::runtime_error("extractSubmatrixSymmetricalBlocks: Indices out of range!");
00824 #endif
00825 out.block(dst_row_blk * block_size,dst_col_blk * block_size, block_size,block_size)
00826 =
00827 derived().block(block_indices[dst_row_blk] * block_size, block_indices[dst_col_blk] * block_size, block_size,block_size);
00828 }
00829 }
00830 }
00831
00832
00833
00834
00835
00836
00837 template <class MATRIX>
00838 void extractSubmatrixSymmetrical(
00839 const std::vector<size_t> &indices,
00840 MATRIX& out) const
00841 {
00842 if (cols()!=rows()) throw std::runtime_error("extractSubmatrixSymmetrical: Matrix is not square.");
00843
00844 const size_t N = indices.size();
00845 const size_t nrows_out=N;
00846 out.resize(nrows_out,nrows_out);
00847 if (!N) return;
00848 for (size_t dst_row_blk=0;dst_row_blk<N; ++dst_row_blk )
00849 for (size_t dst_col_blk=0;dst_col_blk<N; ++dst_col_blk )
00850 out.coeffRef(dst_row_blk,dst_col_blk) = this->coeff(indices[dst_row_blk],indices[dst_col_blk]);
00851 }
00852