$search
00001 #include "../config.h" 00002 #ifdef __MATRIXWRAPPER_EIGEN__ 00003 00004 #include "matrix_EIGEN.h" 00005 #include "vector_EIGEN.h" 00006 00007 #include <Eigen/LU> 00008 00009 using namespace std; 00010 00011 // Passing the constructor arguments... 00012 MyMatrix::Matrix() : EigenMatrix() {} 00013 MyMatrix::Matrix(int num_rows, int num_cols) : EigenMatrix(num_rows, 00014 num_cols){} 00015 00016 // Destructor 00017 MyMatrix::~Matrix(){} 00018 00019 // Copy constructor 00020 MyMatrix::Matrix(const MyMatrix& a) : EigenMatrix(a){} 00021 00022 // ill-designed 00023 MyMatrix::Matrix(const EigenMatrix & a) : EigenMatrix(a){} 00024 00025 MyMatrix::Matrix(int num_rows,const RowVector& v):EigenMatrix(num_rows,v.columns()){ 00026 EigenMatrix & m = *this; 00027 const EigenRowVector & r = v; 00028 for(unsigned int i=0;i<num_rows;i++) 00029 m.row(i) = r; 00030 } 00031 00032 MyRowVector MyMatrix::operator[](unsigned int i) const{ 00033 return this->rowCopy(i); 00034 } 00035 00036 // Size/Capacity 00037 unsigned int MyMatrix::size() const { return this->rows();} 00038 unsigned int MyMatrix::capacity() const { return this->rows();} 00039 00040 // Number of Rows/Cols 00041 unsigned int MyMatrix::rows() const { return ((const EigenMatrix *)this)->rows();} 00042 unsigned int MyMatrix::columns() const { return ((const EigenMatrix *)this)->cols();} 00043 00044 // MATRIX - SCALAR operators 00045 MyMatrix& MyMatrix::operator+= (double a) 00046 { 00047 EigenMatrix & op1 = *this; 00048 op1 += EigenMatrix::Constant(op1.rows(), op1.cols(), a); 00049 return (MyMatrix&)op1; 00050 } 00051 00052 MyMatrix& MyMatrix::operator-= (double a) 00053 { 00054 EigenMatrix & op1 = (*this); 00055 op1 -= EigenMatrix::Constant(op1.rows(), op1.cols(), a); 00056 return (MyMatrix&) op1; 00057 } 00058 00059 MyMatrix& MyMatrix::operator*= (double a) 00060 { 00061 EigenMatrix & op1 = (*this); 00062 op1 *= a; 00063 return *this; 00064 } 00065 00066 MyMatrix& MyMatrix::operator/= (double a) 00067 { 00068 EigenMatrix & op1 = (*this); 00069 op1 /= a; 00070 return (MyMatrix&) op1; 00071 } 00072 00073 MyMatrix MyMatrix::operator+ (double a) const 00074 { 00075 return (MyMatrix)(((EigenMatrix)(*this)) + EigenMatrix::Constant(rows(), cols(), a)); 00076 } 00077 00078 MyMatrix MyMatrix::operator- (double a) const 00079 { 00080 return (MyMatrix)(((EigenMatrix)(*this)) - EigenMatrix::Constant(rows(), cols(), a)); 00081 } 00082 00083 MyMatrix MyMatrix::operator* (double a) const 00084 { 00085 const EigenMatrix& op1 = (*this); 00086 return (MyMatrix) (op1 * a); 00087 } 00088 00089 MyMatrix MyMatrix::operator/ (double a) const 00090 { 00091 const EigenMatrix& op1 = (*this); 00092 return (MyMatrix) (op1 / a); 00093 } 00094 00095 MyMatrix& 00096 MyMatrix::operator =(const MySymmetricMatrix& a) 00097 { 00098 *this =(MyMatrix) a; 00099 00100 return *this; 00101 } 00102 00103 // MATRIX - MATRIX Operators 00104 MyMatrix MyMatrix::operator- (const MyMatrix& a) const 00105 { 00106 const EigenMatrix& op1 = *this; 00107 const EigenMatrix& op2 = a; 00108 00109 return (MyMatrix)(op1 - op2); 00110 } 00111 00112 MyMatrix MyMatrix::operator+ (const MyMatrix& a) const 00113 { 00114 const EigenMatrix& op1 = *this; 00115 const EigenMatrix& op2 = a; 00116 00117 return (MyMatrix)(op1 + op2); 00118 } 00119 00120 MyMatrix MyMatrix::operator* (const MyMatrix& a) const 00121 { 00122 const EigenMatrix& op1 = *this; 00123 const EigenMatrix& op2 = a; 00124 00125 return (MyMatrix)(op1 * op2); 00126 } 00127 00128 MyMatrix & MyMatrix::operator+= (const MyMatrix& a) 00129 { 00130 EigenMatrix & op1 = (*this); 00131 const EigenMatrix & op2 = a; 00132 op1 += op2; 00133 return (MyMatrix &) op1; 00134 } 00135 00136 MyMatrix & MyMatrix::operator-= (const MyMatrix& a) 00137 { 00138 EigenMatrix & op1 = (*this); 00139 const EigenMatrix & op2 = a; 00140 op1 -= op2; 00141 return (MyMatrix &) op1; 00142 } 00143 00144 00145 // MATRIX - VECTOR Operators 00146 MyColumnVector MyMatrix::operator* (const MyColumnVector &b) const 00147 { 00148 const EigenMatrix& op1 = (*this); 00149 return (MyColumnVector) (op1 * ((const EigenColumnVector&)b)); 00150 } 00151 00152 00153 00154 double& MyMatrix::operator()(unsigned int a, unsigned int b) 00155 { 00156 EigenMatrix & op1 = (*this); 00157 return op1(a-1,b-1); 00158 } 00159 00160 double MyMatrix::operator()(unsigned int a, unsigned int b) const 00161 { 00162 const EigenMatrix & op1(*this); 00163 return op1(a-1,b-1); 00164 } 00165 00166 bool MyMatrix::operator==(const MyMatrix& a) const 00167 { 00168 if (this->rows() != a.rows()) return false; 00169 if (this->columns() != a.columns()) return false; 00170 return(((EigenMatrix)(*this)-(EigenMatrix)a).isApproxToConstant(0.0)); 00171 } 00172 00173 00174 // Set all elements equal to a 00175 MyMatrix& 00176 MyMatrix::operator=(double a) 00177 { 00178 ((EigenMatrix&)(*this)).setConstant(a); 00179 return *this; 00180 } 00181 00182 00183 MyRowVector MyMatrix::rowCopy(unsigned int r) const 00184 { 00185 return (MyRowVector) (*this).row(r); 00186 } 00187 00188 MyColumnVector MyMatrix::columnCopy(unsigned int c) const 00189 { 00190 return (MyColumnVector) (*this).col(c); 00191 } 00192 00193 00194 00195 00196 MyMatrix MyMatrix::transpose() const 00197 { 00198 const EigenMatrix &op1 = (*this); 00199 return (MyMatrix) op1.transpose(); 00200 } 00201 00202 double MyMatrix::determinant() const 00203 { 00204 unsigned int r = this->rows(); 00205 assert(r == this->columns()); 00206 const EigenMatrix& A = (*this); 00207 return A.determinant(); 00208 } 00209 00210 00211 MyMatrix MyMatrix::inverse() const 00212 { 00213 unsigned int r = this->rows(); 00214 assert(r == this->columns()); 00215 const EigenMatrix& A = (*this); 00216 return (MyMatrix) A.inverse(); 00217 } 00218 00219 00220 int 00221 MyMatrix::convertToSymmetricMatrix(MySymmetricMatrix& sym) 00222 { 00223 // test if matrix is square matrix 00224 assert(this->rows() == this->columns()); 00225 00226 const EigenMatrix & A = (EigenMatrix &) (*this); 00227 sym = MySymmetricMatrix(A.selfadjointView<Eigen::Upper>()); 00228 return 0; 00229 } 00230 00231 void 00232 MyMatrix::resize(unsigned int i, unsigned int j, bool copy, bool initialize) 00233 { 00234 EigenMatrix & temp = (EigenMatrix &) (*this); 00235 temp.resize(i,j); 00236 } 00237 00238 // get sub matrix 00239 MyMatrix MyMatrix::sub(int i_start, int i_end, int j_start , int j_end) const 00240 { 00241 const EigenMatrix & A = (EigenMatrix &) (*this); 00242 MyMatrix submatrix(A.block(i_start-1,j_start-1,i_end-i_start+1,j_end-j_start+1)); 00243 return submatrix; 00244 } 00245 00247 // CLASS SYMMETRIC MATRIX // 00249 00250 MySymmetricMatrix::SymmetricMatrix() : EigenSymmetricMatrix() {} 00251 MySymmetricMatrix::SymmetricMatrix(int n) : EigenSymmetricMatrix(n,n) {} 00252 MySymmetricMatrix::SymmetricMatrix(int num_rows,const RowVector& v):EigenSymmetricMatrix(num_rows,v.columns()){ 00253 EigenSymmetricMatrix & m = *this; 00254 const EigenRowVector & r = v; 00255 for(unsigned int i=0;i<num_rows;i++) 00256 m.row(i) = r; 00257 } 00258 00259 MyRowVector MySymmetricMatrix::operator[](unsigned int i) const{ 00260 return this->rowCopy(i); 00261 } 00262 00263 00264 00265 // Copy constructor 00266 MySymmetricMatrix::SymmetricMatrix(const SymmetricMatrix& a) : EigenSymmetricMatrix(a){} 00267 MySymmetricMatrix::SymmetricMatrix(const EigenSymmetricMatrix& a) : EigenSymmetricMatrix(a){} 00268 MySymmetricMatrix::SymmetricMatrix(const EigenSymmetricView & a) : EigenSymmetricMatrix(a){} 00269 00270 // Destructor 00271 MySymmetricMatrix::~SymmetricMatrix(){} 00272 00273 // Size/Capacity 00274 unsigned int MySymmetricMatrix::size() const { return this->rows();} 00275 unsigned int MySymmetricMatrix::capacity() const { return this->rows();} 00276 00277 // Ask Number of Rows and Columns 00278 unsigned int MySymmetricMatrix::rows() const { return ((const EigenSymmetricMatrix *)this)->rows();} 00279 unsigned int MySymmetricMatrix::columns() const { return ((const EigenSymmetricMatrix *)this)->cols();} 00280 00281 00282 MyRowVector MySymmetricMatrix::rowCopy(unsigned int r) const 00283 { 00284 00285 unsigned int cols = columns(); 00286 EigenRowVector temp(cols); 00287 for (unsigned int i=0; i<cols; i++) 00288 temp(i) = (*this)(r,i+1); 00289 return (MyRowVector) temp; 00290 } 00291 00292 MySymmetricMatrix MySymmetricMatrix::transpose() const {return (*this);} 00293 00294 MySymmetricMatrix MySymmetricMatrix::inverse() const 00295 { 00296 unsigned int r = this->rows(); 00297 assert(r == this->columns()); 00298 const EigenSymmetricMatrix& A = (*this); 00299 // EigenSymmetricView A = ((const EigenSymmetricMatrix *)(this))->selfadjointView<Eigen::Upper>(); 00300 return MySymmetricMatrix(A.inverse()); 00301 } 00302 00303 double MySymmetricMatrix::determinant() const 00304 { 00305 unsigned int r = this->rows(); 00306 assert(r == this->columns()); 00307 const EigenSymmetricMatrix& A = (*this); 00308 // EigenSymmetricView A = ((const EigenSymmetricMatrix *)(this))->selfadjointView<Eigen::Upper>(); 00309 return A.determinant(); 00310 } 00311 00312 00313 // Set all elements equal to a 00314 MySymmetricMatrix& MySymmetricMatrix::operator=(const double a) 00315 { 00316 ((EigenSymmetricMatrix&)(*this)).setConstant(a); 00317 return *this; 00318 } 00319 00320 00321 // SYMMETRICMATRIX - SCALAR operators 00322 MySymmetricMatrix& MySymmetricMatrix::operator +=(double a) 00323 { 00324 EigenSymmetricMatrix & op1 = *this; 00325 op1 += EigenSymmetricMatrix::Constant(op1.rows(), op1.cols(), a); 00326 return (MySymmetricMatrix&)op1; 00327 } 00328 00329 MySymmetricMatrix& MySymmetricMatrix::operator -=(double a) 00330 { 00331 EigenSymmetricMatrix & op1 = *this; 00332 op1 -= EigenSymmetricMatrix::Constant(op1.rows(), op1.cols(), a); 00333 return (MySymmetricMatrix&)op1; 00334 } 00335 00336 MySymmetricMatrix& MySymmetricMatrix::operator *=(double b) 00337 { 00338 EigenSymmetricMatrix & op1 = (*this); 00339 op1 *= b; 00340 return (MySymmetricMatrix&) op1; 00341 } 00342 00343 MySymmetricMatrix& MySymmetricMatrix::operator /=(double b) 00344 { 00345 EigenSymmetricMatrix & op1 = (*this); 00346 op1 /= b; 00347 return (MySymmetricMatrix&) op1; 00348 } 00349 00350 MySymmetricMatrix MySymmetricMatrix::operator +(double a) const 00351 { 00352 return (MySymmetricMatrix)(((EigenSymmetricMatrix)(*this)) + EigenSymmetricMatrix::Constant(rows(), cols(), a)); 00353 } 00354 00355 MySymmetricMatrix MySymmetricMatrix::operator -(double a) const 00356 { 00357 return (MySymmetricMatrix)(((EigenSymmetricMatrix)(*this)) - EigenSymmetricMatrix::Constant(rows(), cols(), a)); 00358 } 00359 00360 MySymmetricMatrix MySymmetricMatrix::operator *(double b) const 00361 { 00362 const EigenSymmetricMatrix& op1 = (*this); 00363 return (MySymmetricMatrix) (op1 * b); 00364 } 00365 00366 MySymmetricMatrix MySymmetricMatrix::operator /(double b) const 00367 { 00368 const EigenSymmetricMatrix& op1 = (*this); 00369 return (MySymmetricMatrix) (op1 / b); 00370 } 00371 00372 00373 00374 00375 // SYMMETRICMATRIX - MATRIX operators 00376 MyMatrix& MySymmetricMatrix::operator +=(const MyMatrix& a) 00377 { 00378 EigenSymmetricMatrix & op1 = (*this); 00379 op1 += a; 00380 return (MyMatrix &) op1; 00381 } 00382 00383 MyMatrix& MySymmetricMatrix::operator -=(const MyMatrix& a) 00384 { 00385 EigenSymmetricMatrix & op1 = (*this); 00386 op1 -= a; 00387 return (MyMatrix &) op1; 00388 } 00389 00390 00391 MyMatrix MySymmetricMatrix::operator+ (const MyMatrix &a) const 00392 { 00393 const EigenSymmetricMatrix& op1 = *this; 00394 const EigenMatrix& op2 = a; 00395 00396 return (MyMatrix) (op1 + op2); 00397 } 00398 00399 MyMatrix MySymmetricMatrix::operator- (const MyMatrix &a) const 00400 { 00401 const EigenSymmetricMatrix& op1 = *this; 00402 const EigenMatrix& op2 = a; 00403 00404 return (MyMatrix) (op1 - op2); 00405 } 00406 00407 MyMatrix MySymmetricMatrix::operator* (const MyMatrix &a) const 00408 { 00409 const EigenSymmetricMatrix& op1 = *this; 00410 const EigenMatrix& op2 = a; 00411 00412 return (MyMatrix) (op1 * op2); 00413 } 00414 00415 00416 00417 // SYMMETRICMATRIX - SYMMETRICMATRIX operators 00418 MySymmetricMatrix& MySymmetricMatrix::operator +=(const MySymmetricMatrix& a) 00419 { 00420 EigenSymmetricMatrix & op1 = (*this); 00421 const EigenSymmetricMatrix & op2 = a; 00422 op1 += op2; 00423 return (MySymmetricMatrix &) op1; 00424 } 00425 00426 MySymmetricMatrix& MySymmetricMatrix::operator -=(const MySymmetricMatrix& a) 00427 { 00428 EigenSymmetricMatrix & op1 = (*this); 00429 const EigenSymmetricMatrix & op2 = a; 00430 op1 -= op2; 00431 return (MySymmetricMatrix &) op1; 00432 } 00433 00434 MySymmetricMatrix MySymmetricMatrix::operator+ (const MySymmetricMatrix &a) const 00435 { 00436 const EigenSymmetricMatrix& op1 = *this; 00437 const EigenSymmetricMatrix& op2 = a; 00438 00439 return (MySymmetricMatrix) (op1 + op2); 00440 } 00441 00442 MySymmetricMatrix MySymmetricMatrix::operator- (const MySymmetricMatrix &a) const 00443 { 00444 const EigenSymmetricMatrix& op1 = *this; 00445 const EigenSymmetricMatrix& op2 = a; 00446 00447 return (MySymmetricMatrix) (op1 - op2); 00448 } 00449 00450 MyMatrix MySymmetricMatrix::operator* (const MySymmetricMatrix &a) const 00451 { 00452 const EigenSymmetricMatrix& op1 = *this; 00453 const EigenSymmetricMatrix& op2 = a; 00454 00455 return (MyMatrix) (op1 * op2); 00456 } 00457 00458 00459 00460 00461 MyColumnVector MySymmetricMatrix::operator* (const MyColumnVector &b) const 00462 { 00463 const EigenSymmetricMatrix& op1 = (EigenSymmetricMatrix) *this; 00464 return (MyColumnVector) (op1 * ((const EigenColumnVector&)b)); 00465 } 00466 00467 void MySymmetricMatrix::multiply (const MyColumnVector &b, MyColumnVector &result) const 00468 { 00469 const EigenSymmetricMatrix& op1 = (EigenSymmetricMatrix) *this; 00470 result = (MyColumnVector) (op1 * ((const EigenColumnVector&)b)); 00471 } 00472 00473 MyMatrix MySymmetricMatrix::sub(int i_start, int i_end, int j_start , int j_end) const 00474 { 00475 MyMatrix submatrix(i_end-i_start+1, j_end-j_start+1); 00476 for (int i=i_start; i<=i_end; i++) 00477 for (int j=j_start; j<=j_end; j++) 00478 submatrix(i-i_start+1,j-j_start+1) = (*this)(i,j); 00479 00480 return submatrix; 00481 } 00482 00483 00484 00485 double& MySymmetricMatrix::operator()(unsigned int a, unsigned int b) 00486 { 00487 EigenSymmetricMatrix & op1 = (*this); 00488 return op1(a-1,b-1); 00489 } 00490 00491 double MySymmetricMatrix::operator()(unsigned int a, unsigned int b) const 00492 { 00493 const EigenSymmetricMatrix & op1(*this); 00494 return op1(a-1,b-1); 00495 } 00496 00497 bool MySymmetricMatrix::operator==(const MySymmetricMatrix& a) const 00498 { 00499 if (this->rows() != a.rows()) return false; 00500 if (this->columns() != a.columns()) return false; 00501 return(((EigenSymmetricMatrix)(*this)-(EigenSymmetricMatrix)a).isApproxToConstant(0.0)); 00502 } 00503 00504 void 00505 MySymmetricMatrix::resize(unsigned int i, bool copy, bool initialize) 00506 { 00507 EigenSymmetricMatrix & temp = (EigenSymmetricMatrix &) (*this); 00508 temp.resize(i,i); 00509 } 00510 00511 00512 #endif