NGSolve
4.9
|
00001 #ifndef FILE_NGBLA_EXPR 00002 #define FILE_NGBLA_EXPR 00003 00004 /**************************************************************************/ 00005 /* File: expr.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 01. Jan. 02 */ 00008 /**************************************************************************/ 00009 00010 00011 00012 namespace ngbla 00013 { 00014 00015 template <int H, int W, typename T> class Mat; 00016 template <int S, typename T> class Vec; 00017 00018 00019 /* 00020 Matrix expression templates 00021 */ 00022 00023 00024 template <typename TVEC> 00025 inline typename TVEC::TELEM & Access (TVEC & vec, int nr) 00026 { 00027 return vec(nr); 00028 } 00029 00030 template <typename TVEC> 00031 inline typename TVEC::TELEM Access (const TVEC & vec, int nr) 00032 { 00033 return vec(nr); 00034 } 00035 00036 inline double & Access (double & vec, int nr) 00037 { 00038 return vec; 00039 } 00040 00041 inline double Access (const double & vec, int nr) 00042 { 00043 return vec; 00044 } 00045 00046 inline Complex & Access (Complex & vec, int nr) 00047 { 00048 return vec; 00049 } 00050 00051 inline Complex Access (const Complex & vec, int nr) 00052 { 00053 return vec; 00054 } 00055 00056 00057 00058 00059 00060 00061 template <typename TM> 00062 inline typename TM::TELEM & Access (TM & mat, int i, int j) 00063 { 00064 return mat(i,j); 00065 } 00066 00067 template <typename TM> 00068 inline typename TM::TELEM Access (const TM & mat, int i, int j) 00069 { 00070 return mat(i,j); 00071 } 00072 00073 inline double & Access (double & mat, int i, int j) 00074 { 00075 return mat; 00076 } 00077 00078 inline double Access (const double & mat, int i, int j) 00079 { 00080 return mat; 00081 } 00082 00083 inline Complex & Access (Complex & mat, int i, int j) 00084 { 00085 return mat; 00086 } 00087 00088 inline Complex Access (const Complex & mat, int i, int j) 00089 { 00090 return mat; 00091 } 00092 00093 00094 00095 00096 00097 00098 00099 00100 00101 00106 template <class T> 00107 class mat_traits 00108 { 00109 public: 00111 typedef typename T::TELEM TELEM; 00113 typedef typename T::TSCAL TSCAL; 00115 typedef typename T::TV_COL TV_COL; 00117 typedef typename T::TV_ROW TV_ROW; 00119 enum { HEIGHT = T::HEIGHT }; 00121 enum { WIDTH = T::WIDTH }; 00123 enum { IS_COMPLEX = mat_traits<TSCAL>::IS_COMPLEX }; 00124 }; 00125 00126 00127 // some compiler thinks it needs Mat<int>... 00128 template <> 00129 class mat_traits<Mat<1,1,int> > 00130 { 00131 public: 00132 typedef int TELEM; 00133 typedef int TSCAL; 00134 typedef int TV_COL; 00135 typedef int TV_ROW; 00136 enum { HEIGHT = 1 }; 00137 enum { WIDTH = 1 }; 00138 enum { IS_COMPLEX = 0 }; 00139 }; 00140 00141 template <> 00142 class mat_traits<int> 00143 { 00144 public: 00145 typedef int TELEM; 00146 typedef int TSCAL; 00147 typedef int TV_COL; 00148 typedef int TV_ROW; 00149 enum { HEIGHT = 1 }; 00150 enum { WIDTH = 1 }; 00151 enum { IS_COMPLEX = 0 }; 00152 }; 00153 00154 00155 template <> 00156 class mat_traits<double> 00157 { 00158 public: 00159 typedef double TELEM; 00160 typedef double TSCAL; 00161 typedef double TV_COL; 00162 typedef double TV_ROW; 00163 enum { HEIGHT = 1 }; 00164 enum { WIDTH = 1 }; 00165 enum { IS_COMPLEX = 0 }; 00166 }; 00167 00168 template <> 00169 class mat_traits<float> 00170 { 00171 public: 00172 typedef float TELEM; 00173 typedef float TSCAL; 00174 typedef float TV_COL; 00175 typedef float TV_ROW; 00176 enum { HEIGHT = 1 }; 00177 enum { WIDTH = 1 }; 00178 enum { IS_COMPLEX = 0 }; 00179 }; 00180 00181 template <> 00182 class mat_traits<Complex> 00183 { 00184 public: 00185 typedef Complex TELEM; 00186 typedef Complex TSCAL; 00187 typedef Complex TV_COL; 00188 typedef Complex TV_ROW; 00189 enum { HEIGHT = 1 }; 00190 enum { WIDTH = 1 }; 00191 enum { IS_COMPLEX = 1 }; 00192 }; 00193 00194 template <int D, typename SCAL> 00195 class mat_traits<AutoDiff<D,SCAL> > 00196 { 00197 public: 00198 typedef AutoDiff<D,SCAL> TELEM; 00199 typedef AutoDiff<D,SCAL> TSCAL; 00200 typedef AutoDiff<D,SCAL> TV_COL; 00201 typedef AutoDiff<D,SCAL> TV_ROW; 00202 enum { HEIGHT = 1 }; 00203 enum { WIDTH = 1 }; 00204 enum { IS_COMPLEX = mat_traits<SCAL>::IS_COMPLEX }; 00205 }; 00206 00207 template <int D> 00208 class mat_traits<AutoDiffDiff<D> > 00209 { 00210 public: 00211 typedef AutoDiffDiff<D> TELEM; 00212 typedef AutoDiffDiff<D> TSCAL; 00213 typedef AutoDiffDiff<D> TV_COL; 00214 typedef AutoDiffDiff<D> TV_ROW; 00215 enum { HEIGHT = 1 }; 00216 enum { WIDTH = 1 }; 00217 enum { IS_COMPLEX = false }; 00218 }; 00219 00220 00221 template <> 00222 class mat_traits<const int> 00223 { 00224 public: 00225 typedef int TELEM; 00226 typedef int TSCAL; 00227 typedef int TV_COL; 00228 typedef int TV_ROW; 00229 enum { HEIGHT = 1 }; 00230 enum { WIDTH = 1 }; 00231 enum { IS_COMPLEX = 0 }; 00232 }; 00233 00234 template <> 00235 class mat_traits<const double> 00236 { 00237 public: 00238 typedef double TELEM; 00239 typedef double TSCAL; 00240 typedef double TV_COL; 00241 typedef double TV_ROW; 00242 enum { HEIGHT = 1 }; 00243 enum { WIDTH = 1 }; 00244 enum { IS_COMPLEX = 0 }; 00245 }; 00246 00247 template <> 00248 class mat_traits<const Complex> 00249 { 00250 public: 00251 typedef Complex TELEM; 00252 typedef Complex TSCAL; 00253 typedef Complex TV_COL; 00254 typedef Complex TV_ROW; 00255 enum { HEIGHT = 1 }; 00256 enum { WIDTH = 1 }; 00257 enum { IS_COMPLEX = 1 }; 00258 }; 00259 00260 00262 template <typename TV_COL, typename TV_ROW> 00263 class mat_from_vecs 00264 { 00265 enum { HEIGHT = mat_traits<TV_COL>::HEIGHT }; 00266 enum { WIDTH = mat_traits<TV_ROW>::HEIGHT }; 00267 typedef typename mat_from_vecs<typename TV_COL::TELEM, typename TV_ROW::TELEM>::TMAT TELEM; 00268 typedef Mat<HEIGHT,WIDTH,TELEM> TMAT; 00269 }; 00270 00271 template <> class mat_from_vecs<double,double> { typedef double TMAT; }; 00272 template <> class mat_from_vecs<double,Complex> { typedef Complex TMAT; }; 00273 template <> class mat_from_vecs<Complex,double> { typedef Complex TMAT; }; 00274 template <> class mat_from_vecs<Complex,Complex> { typedef Complex TMAT; }; 00275 00276 00277 00279 template <typename TA, typename TB> 00280 class mat_prod_type 00281 { 00282 public: 00283 enum { HEIGHT = mat_traits<TA>::HEIGHT }; 00284 enum { WIDTH = mat_traits<TB>::WIDTH }; 00285 00286 typedef typename mat_prod_type<typename mat_traits<TA>::TELEM, 00287 typename mat_traits<TB>::TELEM>::TMAT TELEM; 00288 00289 typedef Mat<HEIGHT,WIDTH,TELEM> TMAT; 00290 }; 00291 00292 /* 00293 template <int S, typename T> class mat_prod_type<double, Vec<S,T> > { public: typedef Vec<S,T> TMAT; }; 00294 template <int S, typename T> class mat_prod_type<Complex, Vec<S,T> > { public: typedef Vec<S,T> TMAT; }; 00295 template <int S, typename T> class mat_prod_type<double, const FlatVec<S,T> > { public: typedef Vec<S,T> TMAT; }; 00296 template <int S, typename T> class mat_prod_type<Complex, const FlatVec<S,T> > { public: typedef Vec<S,T> TMAT; }; 00297 */ 00298 00299 template <> class mat_prod_type<double,double> { public: typedef double TMAT; }; 00300 template <> class mat_prod_type<double,Complex> { public: typedef Complex TMAT; }; 00301 template <> class mat_prod_type<Complex,double> { public: typedef Complex TMAT; }; 00302 template <> class mat_prod_type<Complex,Complex> { public: typedef Complex TMAT; }; 00303 00304 template <int D, typename TA, typename TB> 00305 class mat_prod_type<AutoDiff<D, TA>, TB> 00306 { public: typedef AutoDiff<D, typename mat_prod_type<TA,TB>::TMAT> TMAT; }; 00307 00308 template <int D, typename TA, typename TB> 00309 class mat_prod_type<TB, AutoDiff<D, TA> > 00310 { public: typedef AutoDiff<D, typename mat_prod_type<TA,TB>::TMAT> TMAT; }; 00311 00312 template <int D, typename TA, int E, typename TB> 00313 class mat_prod_type<AutoDiff<D, TA>, AutoDiff<E,TB> > 00314 { public: typedef AutoDiff<D, typename mat_prod_type<TA,TB>::TMAT> TMAT; }; 00315 00316 00317 00319 template <typename TA, typename TB> 00320 class mat_sum_type 00321 { 00322 public: 00323 enum { HEIGHT = mat_traits<TA>::HEIGHT }; 00324 enum { WIDTH = mat_traits<TA>::WIDTH }; 00325 typedef typename mat_sum_type<typename TA::TELEM, 00326 typename TB::TELEM>::TMAT TELEM; 00327 typedef Mat<HEIGHT,WIDTH,TELEM> TMAT; 00328 }; 00329 00330 template <> class mat_sum_type<float,float> { public: typedef float TMAT; }; 00331 template <> class mat_sum_type<double,double> { public: typedef double TMAT; }; 00332 template <> class mat_sum_type<double,Complex> { public: typedef Complex TMAT; }; 00333 template <> class mat_sum_type<Complex,double> { public: typedef Complex TMAT; }; 00334 template <> class mat_sum_type<Complex,Complex> { public: typedef Complex TMAT; }; 00335 00336 template <int D, typename TA, typename TB> 00337 class mat_sum_type<AutoDiff<D, TA>, TB> 00338 { public: typedef AutoDiff<D, typename mat_sum_type<TA,TB>::TMAT> TMAT; }; 00339 00340 00341 template <int D, typename TA, typename TB> 00342 class mat_sum_type<TB, AutoDiff<D, TA> > 00343 { public: typedef AutoDiff<D, typename mat_sum_type<TA,TB>::TMAT> TMAT; }; 00344 00345 00346 template <int D, typename TA, int E, typename TB> 00347 class mat_sum_type<AutoDiff<D, TA>, AutoDiff<E,TB> > 00348 { public: typedef AutoDiff<D, typename mat_sum_type<TA,TB>::TMAT> TMAT; }; 00349 00350 00351 00352 00353 00354 00355 00357 template <typename TM, typename TS> 00358 class mat_scale_type 00359 { 00360 public: 00361 enum { HEIGHT = mat_traits<TM>::HEIGHT }; 00362 enum { WIDTH = mat_traits<TM>::WIDTH }; 00363 typedef typename mat_scale_type<typename TM::TELEM, TS>::TMAT TELEM; 00364 typedef Mat<HEIGHT,WIDTH,TELEM> TMAT; 00365 }; 00366 00367 template <> class mat_scale_type<double,double> { public: typedef double TMAT; }; 00368 template <> class mat_scale_type<double,Complex> { public: typedef Complex TMAT; }; 00369 template <> class mat_scale_type<Complex,double> { public: typedef Complex TMAT; }; 00370 template <> class mat_scale_type<Complex,Complex> { public: typedef Complex TMAT; }; 00371 00372 template <int D, typename TA, typename TB> 00373 class mat_scale_type<TB, AutoDiff<D, TA> > 00374 { public: typedef AutoDiff<D, TA> TMAT; }; 00375 00376 template <int D, typename TA, typename TB> 00377 class mat_scale_type<AutoDiff<D, TA>, TB> 00378 { public: typedef AutoDiff<D, TA> TMAT; }; 00379 00380 template <int D, typename TA, int E, typename TB> 00381 class mat_scale_type<AutoDiff<D, TA>, AutoDiff<E,TB> > 00382 { public: typedef AutoDiff<D, typename mat_scale_type<TA,TB>::TMAT> TMAT; }; 00383 00384 00385 00386 00388 template <class TM> 00389 inline int Height (const TM & m) 00390 { 00391 return m.Height(); 00392 } 00393 00395 template <class TM> 00396 inline int Width (const TM & m) 00397 { 00398 return m.Width(); 00399 } 00400 00401 template <> inline int Height<double> (const double&) { return 1; } 00402 template <> inline int Height<Complex> (const Complex&) { return 1; } 00403 template <> inline int Width<double> (const double&) { return 1; } 00404 template <> inline int Width<Complex> (const Complex&) { return 1; } 00405 00406 00408 class Complex2RealException : public Exception 00409 { 00410 public: 00411 Complex2RealException() 00412 : Exception("Assignment of Complex 2 Real") { ; } 00413 }; 00414 00415 00416 00417 template <typename TO> 00418 inline TO ConvertTo (double f) 00419 { 00420 return TO(f); 00421 } 00422 00423 template <typename TO> 00424 inline TO ConvertTo (Complex f) 00425 { 00426 return TO(f); 00427 } 00428 00429 00430 template <typename TO> 00431 inline TO ConvertTo (const AutoDiff<1, Complex> & f) 00432 { 00433 return TO(f); 00434 } 00435 00436 00437 template <> 00438 inline double ConvertTo (Complex f) 00439 { 00440 throw Complex2RealException(); 00441 } 00442 00443 00444 template <class TA> class RowsArrayExpr; 00445 template <class TA> class ColsArrayExpr; 00446 template <class TA> class SubMatrixExpr; 00447 00448 00449 00450 00451 00452 #ifdef CHECK_RANGE 00453 00454 class MatrixNotFittingException : public Exception 00455 { 00456 public: 00457 MatrixNotFittingException (const string & where, int h1, int w1, 00458 int h2, int w2) 00459 : Exception ("") 00460 { 00461 stringstream str; 00462 str << where << " mat1 = " << h1 << " x " << w1 << ", mat2 = " << h2 << " x " << w2 << endl; 00463 Append (str.str()); 00464 } 00465 }; 00466 #endif 00467 00468 00477 template <typename T> 00478 class Expr 00479 { 00480 public: 00481 Expr () { ; } 00482 00484 T & Spec() { return static_cast<T&> (*this); } 00485 00487 const T & Spec() const { return static_cast<const T&> (*this); } 00488 00489 00491 int Height() const { return Spec().T::Height(); } 00492 int Width() const { return Spec().T::Width(); } 00493 00494 void Dump (ostream & ost) const { Spec().T::Dump(ost); } 00495 00496 SubMatrixExpr<T> 00497 Rows (int first, int next) 00498 { 00499 return SubMatrixExpr<T> (static_cast<T&> (*this), first, 0, next-first, Width()); 00500 } 00501 00502 SubMatrixExpr<T> 00503 Cols (int first, int next) 00504 { 00505 return SubMatrixExpr<T> (static_cast<T&> (*this), 0, first, Height(), next-first); 00506 } 00507 00508 SubMatrixExpr<T> 00509 Rows (IntRange range) 00510 { 00511 return Rows (range.First(), range.Next()); 00512 } 00513 00514 SubMatrixExpr<T> 00515 Cols (IntRange range) 00516 { 00517 return Cols (range.First(), range.Next()); 00518 } 00519 00520 00521 00522 RowsArrayExpr<T> 00523 Rows (FlatArray<int> rows) 00524 { 00525 return RowsArrayExpr<T> (static_cast<const T&> (*this), rows); 00526 } 00527 00528 ColsArrayExpr<T> 00529 Cols (FlatArray<int> cols) 00530 { 00531 return ColsArrayExpr<T> (static_cast<const T&> (*this), cols); 00532 } 00533 00534 }; 00535 00536 00537 00542 template <typename T> 00543 class SymExpr : public Expr<SymExpr<T> > 00544 { 00545 const T & a; 00546 public: 00547 typedef typename T::TELEM TELEM; 00548 typedef typename mat_traits<TELEM>::TSCAL TSCAL; 00549 00550 SymExpr (const T & aa) : a(aa) { ; } 00551 00552 TELEM operator() (int i) const { return a(i); } 00553 TELEM operator() (int i, int j) const { return a(i,j); } 00554 int Height() const { return a.Height(); } 00555 int Width() const { return a.Width(); } 00556 enum { IS_LINEAR = T::IS_LINEAR }; 00557 void Dump (ostream & ost) const 00558 { ost << "Sym ("; a.Dump(ost); ost << ")"; } 00559 }; 00560 00561 00565 template <typename T> 00566 inline SymExpr<T> Symmetric (const Expr<T> & a) 00567 { 00568 return SymExpr<T> (static_cast<const T&> (a)); 00569 } 00570 00571 00572 00573 template <typename TA> 00574 class LapackExpr : public Expr<LapackExpr<TA> > 00575 { 00576 const TA & a; 00577 public: 00578 LapackExpr (const TA & aa) : a(aa) { ; } 00579 const TA & A() const { return a; } 00580 int Height() const { return a.Height(); } 00581 int Width() const { return a.Width(); } 00582 }; 00583 00584 typedef int T_Lapack; 00585 static T_Lapack Lapack; 00586 00587 template <typename TA> 00588 LapackExpr<TA> operator| (const Expr<TA> & a, T_Lapack /* tl */) 00589 { 00590 return LapackExpr<TA> (a.Spec()); 00591 } 00592 00593 00594 00595 00596 00597 template <typename TA> 00598 class LocalHeapExpr : public Expr<LocalHeapExpr<TA> > 00599 { 00600 const TA & a; 00601 LocalHeap * lh; 00602 public: 00603 LocalHeapExpr (const TA & aa, LocalHeap & alh) : a(aa), lh(&alh) { ; } 00604 const TA & A() const { return a; } 00605 int Height() const { return a.Height(); } 00606 int Width() const { return a.Width(); } 00607 LocalHeap & GetLocalHeap() const { return *lh; } 00608 }; 00609 00610 template <typename TA> 00611 LocalHeapExpr<TA> operator| (const Expr<TA> & a, LocalHeap & lh) 00612 { 00613 return LocalHeapExpr<TA> (a.Spec(), lh); 00614 } 00615 00616 00617 00618 00619 00620 00621 template <class TA, class TB> class MultExpr; 00622 00623 00627 template <class T> 00628 class MatExpr : public Expr<T> 00629 { 00630 public: 00631 00632 MatExpr () { ; } 00633 00634 using Expr<T>::Spec; 00635 using Expr<T>::Height; 00636 using Expr<T>::Width; 00637 00638 enum { IS_LINEAR = 1 }; 00639 00640 void Dump (ostream & ost) const { ost << "Matrix"; } 00641 00642 template<typename TB> 00643 T & operator= (const Expr<TB> & v) 00644 { 00645 #ifdef CHECK_RANGE 00646 if (Height() != v.Height() || Width() != v.Width()) 00647 { 00648 throw MatrixNotFittingException ("operator=", 00649 Height(), Width(), 00650 v.Height(), v.Width()); 00651 } 00652 #endif 00653 00654 if (TB::IS_LINEAR) 00655 { 00656 if (T::IS_LINEAR) 00657 { 00658 int hw = Expr<T>::Height() * Expr<T>::Width(); 00659 for (int i = 0; i < hw; i++) 00660 Spec()(i) = v.Spec()(i); 00661 } 00662 else 00663 { 00664 int h = Expr<T>::Height(); 00665 int w = Expr<T>::Width(); 00666 for (int i = 0, k = 0; i < h; i++) 00667 for (int j = 0; j < w; j++, k++) 00668 Spec()(i,j) = v.Spec()(k); 00669 } 00670 } 00671 else 00672 { 00673 int h = Expr<T>::Height(); 00674 int w = Expr<T>::Width(); 00675 00676 if (T::IS_LINEAR) 00677 for (int i = 0, k = 0; i < h; i++) 00678 for (int j = 0; j < w; j++, k++) 00679 Spec()(k) = v.Spec()(i,j); 00680 else 00681 for (int i = 0; i < h; i++) 00682 for (int j = 0; j < w; j++) 00683 Spec()(i,j) = v.Spec()(i,j); 00684 } 00685 return Spec(); 00686 } 00687 00688 00689 template<typename TB> 00690 T & operator+= (const Expr<TB> & v) 00691 { 00692 #ifdef CHECK_RANGE 00693 if (Height() != v.Height() || Width() != v.Width()) 00694 throw MatrixNotFittingException ("operator+=", 00695 Height(), Width(), 00696 v.Height(), v.Width()); 00697 #endif 00698 00699 if (TB::IS_LINEAR) 00700 { 00701 if (T::IS_LINEAR) 00702 { 00703 int hw = Height() * Width(); 00704 for (int i = 0; i < hw; i++) 00705 Spec()(i) += v.Spec()(i); 00706 } 00707 else 00708 { 00709 int h = Height(); 00710 int w = Width(); 00711 for (int i = 0, k = 0; i < h; i++) 00712 for (int j = 0; j < w; j++, k++) 00713 Spec()(i,j) += v.Spec()(k); 00714 } 00715 } 00716 else 00717 { 00718 int h = Height(); 00719 int w = Width(); 00720 00721 if (T::IS_LINEAR) 00722 for (int i = 0, k = 0; i < h; i++) 00723 for (int j = 0; j < w; j++, k++) 00724 Spec()(k) += v.Spec()(i,j); 00725 else 00726 for (int i = 0; i < h; i++) 00727 for (int j = 0; j < w; j++) 00728 Spec()(i,j) += v.Spec()(i,j); 00729 } 00730 return Spec(); 00731 } 00732 00733 template <typename TA, typename TB> 00734 T & operator= (const Expr<LapackExpr<MultExpr<TA, TB> > > & prod) 00735 { 00736 LapackMultAdd (prod.Spec().A().A(), prod.Spec().A().B(), 1.0, Spec(), 0.0); 00737 return Spec(); 00738 } 00739 00740 template <typename TA, typename TB> 00741 T & operator+= (const Expr<LapackExpr<MultExpr<TA, TB> > > & prod) 00742 { 00743 LapackMultAdd (prod.Spec().A().A(), prod.Spec().A().B(), 1.0, Spec(), 1.0); 00744 return Spec(); 00745 } 00746 00747 template <typename TA, typename TB> 00748 T & operator-= (const Expr<LapackExpr<MultExpr<TA, TB> > > & prod) 00749 { 00750 LapackMultAdd (prod.Spec().A().A(), prod.Spec().A().B(), -1.0, Spec(), 1.0); 00751 return Spec(); 00752 } 00753 00754 00755 00756 template<typename TB> 00757 T & operator+= (const Expr<SymExpr<TB> > & v) 00758 { 00759 #ifdef CHECK_RANGE 00760 if (Height() != v.Height() || Width() != v.Width()) 00761 throw MatrixNotFittingException ("operator+=", 00762 Height(), Width(), 00763 v.Height(), v.Width()); 00764 #endif 00765 int h = Height(); 00766 for (int i = 0; i < h; i++) 00767 { 00768 for (int j = 0; j < i; j++) 00769 { 00770 double val = v.Spec()(i,j); 00771 Spec()(i,j) += val; 00772 Spec()(j,i) += val; 00773 } 00774 Spec()(i,i) += v.Spec()(i,i); 00775 } 00776 return Spec(); 00777 } 00778 00779 00780 /* 00781 MatExpr<T> & operator+= (double scal) 00782 { 00783 int hw = Height() * Width(); 00784 for (int i = 0; i < hw; i++) 00785 Spec()(i) += scal; 00786 return *this; 00787 } 00788 */ 00789 00790 template<typename TB> 00791 MatExpr<T> & operator-= (const Expr<TB> & v) 00792 { 00793 #ifdef CHECK_RANGE 00794 if (Height() != v.Height() || Width() != v.Width()) 00795 throw MatrixNotFittingException ("operator-=", 00796 Height(), Width(), 00797 v.Height(), v.Width()); 00798 #endif 00799 if (TB::IS_LINEAR) 00800 { 00801 if (T::IS_LINEAR) 00802 { 00803 int hw = Height() * Width(); 00804 for (int i = 0; i < hw; i++) 00805 Spec()(i) -= v.Spec()(i); 00806 } 00807 else 00808 { 00809 int h = Height(); 00810 int w = Width(); 00811 for (int i = 0, k = 0; i < h; i++) 00812 for (int j = 0; j < w; j++, k++) 00813 Spec()(i,j) -= v.Spec()(k); 00814 } 00815 } 00816 else 00817 { 00818 int h = Height(); 00819 int w = Width(); 00820 00821 if (T::IS_LINEAR) 00822 for (int i = 0, k = 0; i < h; i++) 00823 for (int j = 0; j < w; j++, k++) 00824 Spec()(k) -= v.Spec()(i,j); 00825 else 00826 for (int i = 0; i < h; i++) 00827 for (int j = 0; j < w; j++) 00828 Spec()(i,j) -= v.Spec()(i,j); 00829 } 00830 return Spec(); 00831 } 00832 00833 00834 template <class SCAL2> 00835 T & operator*= (const SCAL2 & s) 00836 { 00837 if (T::IS_LINEAR) 00838 { 00839 int hw = Height() * Width(); 00840 for (int i = 0; i < hw; i++) 00841 Spec()(i) *= s; 00842 } 00843 else 00844 for (int i = 0; i < Height(); i++) 00845 for (int j = 0; j < Width(); j++) 00846 Spec()(i,j) *= s; 00847 00848 return Spec(); 00849 } 00850 00851 template <class SCAL2> 00852 T & operator/= (const SCAL2 & s) 00853 { 00854 return (*this) *= (1./s); 00855 } 00856 }; 00857 00858 00859 00860 00861 00862 00863 00864 00865 00871 template <class T> 00872 class CMCPMatExpr : public MatExpr<T> 00873 { 00874 public: 00875 // int Height() const { return Spec().T::Height(); } 00876 // int Width() const { return Spec().T::Width(); } 00877 00878 CMCPMatExpr () { ; } 00879 00880 using MatExpr<T>::Spec; 00881 using MatExpr<T>::Height; 00882 using MatExpr<T>::Width; 00883 00884 // T & Spec() { return static_cast<T&> (*this); } 00885 // const T & Spec() const { return static_cast<const T&> (*this); } 00886 // enum { IS_LINEAR = 1 }; 00887 00888 template<typename TB> 00889 const T & operator= (const Expr<TB> & v) const 00890 { 00891 const_cast<CMCPMatExpr*> (this) -> MatExpr<T>::operator= (v); 00892 return Spec(); 00893 } 00894 00895 template<typename TB> 00896 const T & operator+= (const Expr<TB> & v) const 00897 { 00898 const_cast<CMCPMatExpr*> (this) -> MatExpr<T>::operator+= (v); 00899 return Spec(); 00900 } 00901 00902 template<typename TB> 00903 const T & operator+= (const Expr<SymExpr<TB> > & v) const 00904 { 00905 const_cast<CMCPMatExpr*> (this) -> MatExpr<T>::operator+= (v); 00906 return Spec(); 00907 } 00908 00909 template<typename TB> 00910 const T & operator-= (const Expr<TB> & v) const 00911 { 00912 const_cast<CMCPMatExpr*> (this) -> MatExpr<T>::operator-= (v); 00913 return Spec(); 00914 } 00915 00916 template <class SCAL2> 00917 const T & operator*= (const SCAL2 & s) const 00918 { 00919 const_cast<CMCPMatExpr*> (this) -> MatExpr<T>::operator*= (s); 00920 return Spec(); 00921 } 00922 00923 template <class SCAL2> 00924 const T & operator/= (const SCAL2 & s) const 00925 { 00926 return (*this) *= (1.0/s); 00927 } 00928 00929 SubMatrixExpr<T> 00930 Rows (int first, int next) const 00931 { 00932 return SubMatrixExpr<T> (static_cast<const T&> (*this), first, 0, next-first, Width()); 00933 } 00934 00935 SubMatrixExpr<T> 00936 Cols (int first, int next) const 00937 { 00938 return SubMatrixExpr<T> (static_cast<const T&> (*this), 0, first, Height(), next-first); 00939 } 00940 00941 SubMatrixExpr<T> 00942 Rows (IntRange range) const 00943 { 00944 return Rows (range.First(), range.Next()); 00945 } 00946 00947 SubMatrixExpr<T> 00948 Cols (IntRange range) const 00949 { 00950 return Cols (range.First(), range.Next()); 00951 } 00952 00953 00954 00955 RowsArrayExpr<T> 00956 Rows (FlatArray<int> rows) const 00957 { 00958 return RowsArrayExpr<T> (static_cast<const T&> (*this), rows); 00959 } 00960 00961 ColsArrayExpr<T> 00962 Cols (FlatArray<int> cols) const 00963 { 00964 return ColsArrayExpr<T> (static_cast<const T&> (*this), cols); 00965 } 00966 }; 00967 00968 00969 00970 00971 00972 00973 00974 00975 00976 00977 00978 00979 00980 00981 00982 00983 00984 /* *************************** SumExpr **************************** */ 00985 00990 template <class TA, class TB> 00991 class SumExpr : public Expr<SumExpr<TA,TB> > 00992 { 00993 const TA & a; 00994 const TB & b; 00995 public: 00996 typedef typename mat_sum_type<typename TA::TELEM, 00997 typename TB::TELEM>::TMAT TELEM; 00998 // typedef typename mat_traits<TELEM>::TSCAL TSCAL; 00999 enum { IS_LINEAR = TA::IS_LINEAR && TB::IS_LINEAR }; 01000 01001 SumExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; } 01002 01003 TELEM operator() (int i) const { return a(i)+b(i); } 01004 TELEM operator() (int i, int j) const { return a(i,j)+b(i,j); } 01005 int Height() const { return a.Height(); } 01006 int Width() const { return a.Width(); } 01007 01008 void Dump (ostream & ost) const 01009 { ost << "("; a.Dump(ost); ost << ") + ("; b.Dump(ost); ost << ")"; } 01010 }; 01011 01012 template <typename TA, typename TB> 01013 inline SumExpr<TA, TB> 01014 operator+ (const Expr<TA> & a, const Expr<TB> & b) 01015 { 01016 return SumExpr<TA, TB> (a.Spec(), b.Spec()); 01017 } 01018 01019 01020 01021 01022 /* *************************** SubExpr **************************** */ 01023 01024 01029 template <class TA, class TB> 01030 class SubExpr : public Expr<SubExpr<TA,TB> > 01031 { 01032 const TA & a; 01033 const TB & b; 01034 public: 01035 typedef typename mat_sum_type<typename TA::TELEM, 01036 typename TB::TELEM>::TMAT TELEM; 01037 // typedef typename mat_traits<TELEM>::TSCAL TSCAL; 01038 enum { IS_LINEAR = TA::IS_LINEAR && TB::IS_LINEAR }; 01039 01040 SubExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; } 01041 01042 TELEM operator() (int i) const { return a(i)-b(i); } 01043 TELEM operator() (int i, int j) const { return a(i,j)-b(i,j); } 01044 int Height() const { return a.Height(); } 01045 int Width() const { return a.Width(); } 01046 }; 01047 01048 01049 template <typename TA, typename TB> 01050 inline SubExpr<TA, TB> 01051 operator- (const Expr<TA> & a, const Expr<TB> & b) 01052 { 01053 return SubExpr<TA, TB> (a.Spec(), b.Spec()); 01054 } 01055 01056 01057 01058 01059 01060 01061 01062 /* *************************** MinusExpr **************************** */ 01063 01064 01069 template <class TA> 01070 class MinusExpr : public Expr<MinusExpr<TA> > 01071 { 01072 const TA & a; 01073 public: 01074 typedef typename TA::TELEM TELEM; 01075 01076 MinusExpr (const TA & aa) : a(aa) { ; } 01077 01078 TELEM operator() (int i) const { return -a(i); } 01079 TELEM operator() (int i, int j) const { return -a(i,j); } 01080 int Height() const { return a.Height(); } 01081 int Width() const { return a.Width(); } 01082 01083 enum { IS_LINEAR = TA::IS_LINEAR }; 01084 }; 01085 01086 01087 template <typename TA> 01088 inline MinusExpr<TA> 01089 operator- (const Expr<TA> & a) 01090 { 01091 return MinusExpr<TA> (a.Spec()); 01092 } 01093 01094 01095 /* *************************** ScaleExpr **************************** */ 01096 01097 01101 template <class TA, class TS> 01102 class ScaleExpr : public Expr<ScaleExpr<TA,TS> > 01103 { 01104 const TA & a; 01105 TS s; 01106 public: 01107 typedef typename mat_scale_type<typename TA::TELEM, TS>::TMAT TELEM; 01108 typedef typename mat_traits<TELEM>::TSCAL TSCAL; 01109 enum { IS_LINEAR = TA::IS_LINEAR }; 01110 01111 ScaleExpr (const TA & aa, TS as) : a(aa), s(as) { ; } 01112 01113 TELEM operator() (int i) const { return s * a(i); } 01114 TELEM operator() (int i, int j) const { return s * a(i,j); } 01115 01116 int Height() const { return a.Height(); } 01117 int Width() const { return a.Width(); } 01118 void Dump (ostream & ost) const 01119 { ost << "Scale, s=" << s << " * "; a.Dump(ost); } 01120 }; 01121 01122 template <typename TA> 01123 inline ScaleExpr<TA, double> 01124 operator* (double b, const Expr<TA> & a) 01125 { 01126 return ScaleExpr<TA, double> (a.Spec(), b); 01127 } 01128 01129 template <typename TA> 01130 inline ScaleExpr<TA, Complex> 01131 operator* (Complex b, const Expr<TA> & a) 01132 { 01133 return ScaleExpr<TA, Complex> (a.Spec(), b); 01134 } 01135 01136 template <int D, typename TAD, typename TA> 01137 inline ScaleExpr<TA, AutoDiff<D,TAD> > 01138 operator* (const AutoDiff<D,TAD> & b, const Expr<TA> & a) 01139 { 01140 return ScaleExpr<TA, AutoDiff<D,TAD> > (a.Spec(), b ); 01141 } 01142 01143 01144 01145 01146 /* ************************* MultExpr ************************* */ 01147 01148 01152 template <class TA, class TB> class MultExpr : public Expr<MultExpr<TA,TB> > 01153 { 01154 const TA & a; 01155 const TB & b; 01156 public: 01157 typedef typename mat_prod_type<typename TA::TELEM, 01158 typename TB::TELEM>::TMAT TELEM; 01159 typedef typename mat_traits<TELEM>::TSCAL TSCAL; 01160 01161 MultExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; } 01162 01163 TELEM operator() (int i) const { return operator()(i,0); } // TELEM(TSCAL(0)); } JS, 310508 01164 TELEM operator() (int i, int j) const 01165 { 01166 int wa = a.Width(); 01167 TELEM sum; 01168 if (wa >= 1) 01169 { 01170 sum = a(i,0) * b(0,j); 01171 for (int k = 1; k < wa; k++) 01172 sum += a(i,k) * b(k,j); 01173 } 01174 else 01175 sum = TSCAL(0); 01176 01177 return sum; 01178 } 01179 01180 const TA & A() const { return a; } 01181 const TB & B() const { return b; } 01182 int Height() const { return a.Height(); } 01183 int Width() const { return b.Width(); } 01184 enum { IS_LINEAR = 0 }; 01185 }; 01186 01187 01188 template <typename TA, typename TB> 01189 inline MultExpr<TA, TB> 01190 operator* (const Expr<TA> & a, const Expr<TB> & b) 01191 { 01192 return MultExpr<TA, TB> (a.Spec(), b.Spec()); 01193 // (static_cast <const TA&> (a), static_cast <const TB&> (b)); 01194 } 01195 01196 01197 /* ************************** Trans *************************** */ 01198 01199 01200 inline double Trans (double a) { return a; } 01201 inline Complex Trans (Complex a) { return a; } 01202 template<int D, typename TAD> 01203 inline AutoDiff<D,TAD> Trans (const AutoDiff<D,TAD> & a) { return a; } 01204 01205 01209 template <class TA> class TransExpr : public Expr<TransExpr<TA> > 01210 { 01211 const TA & a; 01212 public: 01213 typedef typename TA::TELEM TELEM; 01214 01215 TransExpr (const TA & aa) : a(aa) { ; } 01216 01217 int Height() const { return a.Width(); } 01218 int Width() const { return a.Height(); } 01219 01220 TELEM operator() (int i, int j) const { return Trans (a(j,i)); } 01221 TELEM operator() (int i) const { return 0; } 01222 01223 enum { IS_LINEAR = 0 }; 01224 01225 const TA & A() { return a; } 01226 }; 01227 01228 01230 template <typename TA> 01231 inline TransExpr<TA> 01232 Trans (const Expr<TA> & a) 01233 { 01234 return TransExpr<TA> (a.Spec()); 01235 } 01236 01237 01238 /* ************************* SubMatrix ************************ */ 01239 01243 template <class TA> class SubMatrixExpr : public MatExpr<SubMatrixExpr<TA> > 01244 { 01245 TA & a; 01246 int first_row, first_col; 01247 int height, width; 01248 public: 01249 typedef typename TA::TELEM TELEM; 01250 typedef typename TA::TSCAL TSCAL; 01251 01252 SubMatrixExpr (TA & aa, int fr, int fc, int ah, int aw) 01253 : a(aa), first_row(fr), first_col(fc), height(ah), width(aw) { ; } 01254 01255 int Height() const { return height; } 01256 int Width() const { return width; } 01257 01258 TELEM & operator() (int i, int j) { return a(i+first_row, j+first_col); } 01259 TELEM & operator() (int i) { return a(i+first_row); } 01260 const TELEM & operator() (int i, int j) const { return a(i+first_row, j+first_col); } 01261 const TELEM & operator() (int i) const { return a(i+first_row); } 01262 01263 enum { IS_LINEAR = 0 }; 01264 01265 template<typename TB> 01266 const SubMatrixExpr & operator= (const Expr<TB> & m) 01267 { 01268 MatExpr<SubMatrixExpr<TA> >::operator= (m); 01269 return *this; 01270 } 01271 }; 01272 01273 01274 01275 /* ************************* RowsArray ************************ */ 01276 01280 template <class TA> class RowsArrayExpr : public MatExpr<RowsArrayExpr<TA> > 01281 { 01282 const TA & a; 01283 FlatArray<int> rows; 01284 public: 01285 typedef typename TA::TELEM TELEM; 01286 typedef typename TA::TREF TREF; 01287 typedef typename TA::TSCAL TSCAL; 01288 01289 RowsArrayExpr (const TA & aa, FlatArray<int> arows) : a(aa), rows(arows) { ; } 01290 01291 int Height() const { return rows.Size(); } 01292 int Width() const { return a.Width(); } 01293 01294 /* 01295 TELEM & operator() (int i, int j) const { return a(rows[i], j); } 01296 TELEM & operator() (int i) const { return a(rows[i]); } 01297 */ 01298 01299 TREF operator() (int i, int j) const { return a(rows[i], j); } 01300 TREF operator() (int i) const { return a(rows[i]); } 01301 01302 01303 enum { IS_LINEAR = 0 }; 01304 01305 template<typename TB> 01306 const RowsArrayExpr & operator= (const Expr<TB> & m) 01307 { 01308 MatExpr<RowsArrayExpr<TA> >::operator= (m); 01309 return *this; 01310 } 01311 01312 const RowsArrayExpr & operator= (const RowsArrayExpr & m) 01313 { 01314 MatExpr<RowsArrayExpr<TA> >::operator= (m); 01315 return *this; 01316 } 01317 }; 01318 01319 01323 template <class TA> class ColsArrayExpr : public MatExpr<ColsArrayExpr<TA> > 01324 { 01325 const TA & a; 01326 FlatArray<int> cols; 01327 public: 01328 typedef typename TA::TELEM TELEM; 01329 typedef typename TA::TSCAL TSCAL; 01330 01331 ColsArrayExpr (const TA & aa, FlatArray<int> acols) : a(aa), cols(acols) { ; } 01332 01333 int Height() const { return a.Height(); } 01334 int Width() const { return cols.Size(); } 01335 01336 TELEM & operator() (int i, int j) const { return a(i, cols[j]); } 01337 TELEM & operator() (int i) const { return a(i, cols[0]); } 01338 01339 enum { IS_LINEAR = 0 }; 01340 01341 template<typename TB> 01342 const ColsArrayExpr & operator= (const Expr<TB> & m) 01343 { 01344 MatExpr<ColsArrayExpr<TA> >::operator= (m); 01345 return *this; 01346 } 01347 01348 }; 01349 01350 01351 01352 01353 /* ************************* Conjugate *********************** */ 01354 01355 01356 inline double Conj (double a) 01357 { 01358 return a; 01359 } 01360 01361 inline Complex Conj (Complex a) 01362 { 01363 return conj(a); 01364 } 01365 01369 // Attention NOT transpose !!! elemwise conjugate !!! 01370 template <class TA> class ConjExpr : public Expr<ConjExpr<TA> > 01371 { 01372 const TA & a; 01373 public: 01374 typedef typename TA::TELEM TELEM; 01375 typedef typename TA::TSCAL TSCAL; 01376 01377 ConjExpr (const TA & aa) : a(aa) { ; } 01378 01379 int Height() const { return a.Height(); } 01380 int Width() const { return a.Width(); } 01381 01382 TELEM operator() (int i, int j) const { return Conj(a(i,j)); } 01383 TELEM operator() (int i) const { return Conj(a(i)); } 01384 01385 enum { IS_LINEAR = 0 }; 01386 }; 01387 01388 01390 template <typename TA> 01391 inline ConjExpr<TA> 01392 Conj (const Expr<TA> & a) 01393 { 01394 return ConjExpr<TA> (static_cast <const TA&> (a)); 01395 } 01396 01397 template<int D, typename TAD> 01398 inline AutoDiff<D,TAD> Conj (const AutoDiff<D,TAD> & a) 01399 { 01400 AutoDiff<D,TAD> b; 01401 b.Value() = conj(a.Value()); 01402 01403 for(int i=0;i<D;i++) 01404 b.DValue(i) = conj(a.DValue(i)); 01405 01406 return b; 01407 } 01408 01409 01410 01411 01412 01413 /* ************************* InnerProduct ********************** */ 01414 01415 01416 inline double InnerProduct (double a, double b) {return a * b;} 01417 inline Complex InnerProduct (Complex a, Complex b) {return a * b;} 01418 inline Complex InnerProduct (double a, Complex b) {return a * b;} 01419 inline Complex InnerProduct (Complex a, double b) {return a * b;} 01420 01421 01425 template <class TA, class TB> 01426 inline 01427 // typename TA::TSCAL 01428 typename mat_prod_type<typename TA::TSCAL, typename TB::TSCAL>::TMAT 01429 InnerProduct (const Expr<TA> & a, const Expr<TB> & b) 01430 { 01431 typedef typename mat_prod_type<typename TA::TSCAL, typename TB::TSCAL>::TMAT TSCAL; 01432 // typedef typename TA::TSCAL TSCAL; 01433 01434 if (a.Height() == 0) return TSCAL(0); 01435 01436 TSCAL sum = InnerProduct (a.Spec()(0), b.Spec()(0)); 01437 for (int i = 1; i < a.Height(); i++) 01438 sum += InnerProduct (a.Spec()(i), b.Spec()(i)); 01439 01440 return sum; 01441 } 01442 01443 01444 /* **************************** Trace **************************** */ 01445 01446 01450 template <class TA> 01451 inline typename TA::TELEM Trace (const Expr<TA> & a) 01452 { 01453 typedef typename TA::TELEM TELEM; 01454 TELEM sum = 0; 01455 for (int i = 0; i < Height(a); i++) 01456 sum += a.Spec()(i,i); 01457 return sum; 01458 } 01459 01460 /* **************************** L2Norm **************************** */ 01461 01463 inline double L2Norm2 (double v) 01464 { 01465 return v*v; 01466 } 01467 01468 inline double L2Norm2 (Complex v) 01469 { 01470 return v.real()*v.real()+v.imag()*v.imag(); 01471 } 01472 01473 template<int D, typename SCAL> 01474 inline double L2Norm2 (const AutoDiff<D,SCAL> & x) throw() 01475 { 01476 return L2Norm2(x.Value()); 01477 } 01478 01479 template <class TA> 01480 inline double L2Norm2 (const Expr<TA> & v) 01481 { 01482 // typedef typename TA::TSCAL TSCAL; 01483 double sum = 0; 01484 if (TA::IS_LINEAR) 01485 for (int i = 0; i < v.Height()*v.Width(); i++) 01486 sum += L2Norm2 (v.Spec()(i)); 01487 else 01488 for (int i = 0; i < v.Height(); i++) 01489 for (int j = 0; j < v.Width(); j++) 01490 sum += L2Norm2 (v.Spec()(i,j)); 01491 01492 return sum; 01493 } 01494 01498 template <class TA> 01499 inline double L2Norm (const Expr<TA> & v) 01500 { 01501 return sqrt (L2Norm2(v)); 01502 } 01503 01504 01505 01506 /* **************************** MaxNorm **************************** */ 01507 01509 inline double MaxNorm (double v) 01510 { 01511 return fabs(v); 01512 } 01513 01514 inline double MaxNorm (Complex v) 01515 { 01516 return fabs(v); 01517 } 01518 01519 template<int D, typename SCAL> 01520 inline double L2Norm (const AutoDiff<D,SCAL> & x) throw() 01521 { 01522 return MaxNorm(x.Value()); 01523 } 01524 01525 template <class TA> 01526 inline double MaxNorm (const Expr<TA> & v) 01527 { 01528 double sum = 0; 01529 01530 if (TA::IS_LINEAR) 01531 for (int i = 0; i < v.Height()*v.Width(); i++) 01532 sum = max(sum, MaxNorm( v.Spec()(i)) ); 01533 else 01534 for (int i = 0; i < v.Height(); i++) 01535 for (int j = 0; j < v.Width(); j++) 01536 sum = max(sum, MaxNorm ( v.Spec()(i,j)) ); 01537 01538 return sum; 01539 } 01540 01541 01542 /* *************************** Output ****************************** */ 01543 01544 01545 01547 template<typename T> 01548 std::ostream & operator<< (std::ostream & s, const Expr<T> & v) 01549 { 01550 for (int i = 0; i < v.Height(); i++) 01551 { 01552 for (int j = 0 ; j < v.Width(); j++) 01553 s << " " << setw(7) << v.Spec()(i,j); 01554 s << endl; 01555 } 01556 return s; 01557 } 01558 01559 01560 01561 /* **************************** Inverse *************************** */ 01562 01563 template <typename T> class FlatMatrix; 01564 template <typename T> class Matrix; 01565 template <int H, int W, typename T> class Mat; 01566 01567 01569 template <class T2> 01570 extern NGS_DLL_HEADER void CalcInverse (FlatMatrix<T2> inv); 01571 01572 extern NGS_DLL_HEADER void CalcInverse (FlatMatrix<double> inv); 01573 01574 template <class T, class T2> 01575 inline void CalcInverse (const FlatMatrix<T> m, FlatMatrix<T2> inv) 01576 { 01577 inv = m; 01578 CalcInverse (inv); 01579 } 01580 01581 template <class T, class T2> 01582 inline void CalcInverse (const FlatMatrix<T> m, Matrix<T2> & inv) 01583 { 01584 inv = m; 01585 CalcInverse (inv); 01586 } 01587 01588 01589 01590 01594 template <typename T> 01595 inline Matrix<T> Inv (const FlatMatrix<T> & m) 01596 { 01597 Matrix<T> inv(m.Height(),m.Height()); 01598 CalcInverse (m, inv); 01599 return inv; 01600 } 01601 01602 inline void CalcInverse (double & m) 01603 { 01604 m = 1 / m; 01605 } 01606 01607 inline void CalcInverse (Complex & m) 01608 { 01609 m = 1.0 / m; 01610 } 01611 01612 template <int H, int W, typename T> 01613 inline void CalcInverse (Mat<H,W,T> & m) 01614 { 01615 FlatMatrix<T> fm(m); 01616 CalcInverse (fm); 01617 } 01618 01619 01620 inline void CalcInverse (const double & m, double & inv) 01621 { 01622 inv = 1 / m; 01623 } 01624 01625 inline void CalcInverse (const Complex & m, Complex & inv) 01626 { 01627 inv = 1.0 / m; 01628 } 01629 01630 template <int H, int W, typename T, typename TINV> 01631 inline void CalcInverse (const Mat<H,W,T> & m, TINV & inv) 01632 { 01633 /* 01634 switch (H) 01635 { 01636 case 1: 01637 { 01638 cout << "mat 1 inv, should not be called" << endl; 01639 inv(0,0) = 1.0 / m(0,0); 01640 return; 01641 } 01642 case 2: 01643 { 01644 cout << "mat 2 inv, should not be called" << endl; 01645 T idet = 1.0 / (m(0,0) * m(1,1) - m(0,1) * m(1,0)); 01646 inv(0,0) = idet * m(1,1); 01647 inv(0,1) = -idet * m(0,1); 01648 inv(1,0) = -idet * m(1,0); 01649 inv(1,1) = idet * m(0,0); 01650 return; 01651 } 01652 case 3: 01653 { 01654 cout << "mat 3 inv, should not be called" << endl; 01655 T h0 = m(4)*m(8)-m(5)*m(7); 01656 T h1 = m(5)*m(6)-m(3)*m(8); 01657 T h2 = m(3)*m(7)-m(4)*m(6); 01658 T det = m(0) * h0 + m(1) * h1 + m(2) * h2; 01659 T idet = 1.0 / det; 01660 01661 inv(0,0) = idet * h0; 01662 inv(0,1) = -idet * (m(1) * m(8) - m(2) * m(7)); 01663 inv(0,2) = idet * (m(1) * m(5) - m(2) * m(4)); 01664 01665 inv(1,0) = idet * h1; 01666 inv(1,1) = idet * (m(0) * m(8) - m(2) * m(6)); 01667 inv(1,2) = -idet * (m(0) * m(5) - m(2) * m(3)); 01668 01669 inv(2,0) = idet * h2; 01670 inv(2,1) = -idet * (m(0) * m(7) - m(1) * m(6)); 01671 inv(2,2) = idet * (m(0) * m(4) - m(1) * m(3)); 01672 return; 01673 } 01674 default: 01675 { 01676 */ 01677 FlatMatrix<T> fm(m); 01678 FlatMatrix<T> finv(inv); 01679 CalcInverse (fm, finv); 01680 } 01681 01682 template <typename T, typename TINV> 01683 inline void CalcInverse (const Mat<1,1,T> & m, TINV & inv) 01684 { 01685 inv(0,0) = 1.0 / m(0,0); 01686 } 01687 01688 template <typename T, typename TINV> 01689 inline void CalcInverse (const Mat<2,2,T> & m, TINV & inv) 01690 { 01691 T idet = 1.0 / (m(0,0) * m(1,1) - m(0,1) * m(1,0)); 01692 inv(0,0) = idet * m(1,1); 01693 inv(0,1) = -idet * m(0,1); 01694 inv(1,0) = -idet * m(1,0); 01695 inv(1,1) = idet * m(0,0); 01696 } 01697 01698 template <typename T, typename TINV> 01699 inline void CalcInverse (const Mat<3,3,T> & m, TINV & inv) 01700 { 01701 T h0 = m(4)*m(8)-m(5)*m(7); 01702 T h1 = m(5)*m(6)-m(3)*m(8); 01703 T h2 = m(3)*m(7)-m(4)*m(6); 01704 T det = m(0) * h0 + m(1) * h1 + m(2) * h2; 01705 T idet = 1.0 / det; 01706 01707 inv(0,0) = idet * h0; 01708 inv(0,1) = -idet * (m(1) * m(8) - m(2) * m(7)); 01709 inv(0,2) = idet * (m(1) * m(5) - m(2) * m(4)); 01710 01711 inv(1,0) = idet * h1; 01712 inv(1,1) = idet * (m(0) * m(8) - m(2) * m(6)); 01713 inv(1,2) = -idet * (m(0) * m(5) - m(2) * m(3)); 01714 01715 inv(2,0) = idet * h2; 01716 inv(2,1) = -idet * (m(0) * m(7) - m(1) * m(6)); 01717 inv(2,2) = idet * (m(0) * m(4) - m(1) * m(3)); 01718 return; 01719 } 01720 01721 01722 01723 template <int H, int W, typename T> 01724 inline Mat<H,W,T> Inv (const Mat<H,W,T> & m) 01725 { 01726 Mat<H,W,T> inv; 01727 CalcInverse (m, inv); 01728 return inv; 01729 } 01730 01731 01732 01733 /* ********************** Determinant *********************** */ 01734 01735 01739 template <class T> 01740 // inline typename T::TSCAL Det (const MatExpr<T> & m) 01741 // inline typename mat_traits<typename T::TELEM>::TSCAL Det (const MatExpr<T> & m) 01742 inline typename T::TELEM Det (const MatExpr<T> & m) 01743 { 01744 const T & sm = m.Spec(); 01745 switch (sm.Height()) 01746 { 01747 case 1: 01748 { 01749 return sm(0,0); 01750 } 01751 case 2: 01752 { 01753 return ( sm(0,0)*sm(1,1) - sm(0,1)*sm(1,0) ); 01754 } 01755 case 3: 01756 { 01757 return 01758 sm(0) * (sm(4) * sm(8) - sm(5) * sm(7)) + 01759 sm(1) * (sm(5) * sm(6) - sm(3) * sm(8)) + 01760 sm(2) * (sm(3) * sm(7) - sm(4) * sm(6)); 01761 } 01762 default: 01763 { 01764 cerr << "general det not implemented" << endl; 01765 } 01766 } 01767 // return typename mat_traits<typename T::TELEM>::TSCAL (0); // typename T::TSCAL(0); 01768 return typename T::TELEM (0); // typename T::TSCAL(0); 01769 } 01770 } 01771 01772 #endif