NGSolve  4.9
basiclinalg/expr.hpp
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