NGSolve  4.9
basiclinalg/matrix.hpp
00001 #ifndef FILE_MATRIX_EXPR
00002 #define FILE_MATRIX_EXPR
00003 
00004 /**************************************************************************/
00005 /* File:   matrix.hpp                                                     */
00006 /* Author: Joachim Schoeberl                                              */
00007 /* Date:   01. Jan. 02                                                    */
00008 /**************************************************************************/
00009 
00010 namespace ngbla
00011 {
00012 
00013   template <int H, int W, typename T> class Mat;
00014   template <typename T> class SliceMatrix;
00015   template <typename T> class DoubleSliceMatrix;
00016 
00017 
00019   extern void CheckMatRange(int h, int w, int i);
00021   extern void CheckMatRange(int h, int w, int i, int j);
00022 
00023 
00024   // compatibility (for a while)
00025   // #define VRange Rows
00026   // #define HRange Cols
00027 
00033   template <typename T = double>
00034   class FlatMatrix : public CMCPMatExpr<FlatMatrix<T> >
00035   {
00036   protected:
00038     int h;
00040     int w;
00042     T * data;
00043   public:
00044 
00046     typedef T TELEM;
00047     typedef T& TREF;
00049     typedef typename mat_traits<T>::TSCAL TSCAL;
00050 
00052     FlatMatrix () throw () { ; }
00053   
00055     FlatMatrix (int ah, int aw, T * adata) throw ()
00056       : h(ah), w(aw), data(adata) { ; }
00057   
00059     FlatMatrix (int ah, T * adata) throw()
00060       : h(ah), w(ah), data(adata) { ; }
00061 
00063     FlatMatrix (int ah, int aw, LocalHeap & lh) // throw (LocalHeapOverflow)
00064     // : h(ah), w(aw), data((T*)lh.Alloc(ah*aw*sizeof(T))) { ; }
00065       : h(ah), w(aw), data (lh.Alloc<T>(ah*aw)) { ; }
00066   
00068     FlatMatrix (int ah, LocalHeap & lh) // throw (LocalHeapOverflow)
00069     // : h(ah), w(ah), data((T*)lh.Alloc(ah*ah*sizeof(T))) { ; }
00070       : h(ah), w(ah), data(lh.Alloc<T>(ah*ah)) { ; }
00071   
00073     FlatMatrix (const FlatMatrix & m) throw () 
00074       : CMCPMatExpr<FlatMatrix> (), h(m.h), w(m.w) , data(m.data)
00075     { ; }
00076   
00078     template<typename TB>
00079     FlatMatrix (const LocalHeapExpr<TB> & m2) 
00080     {
00081       h = m2.A().Height();
00082       w = m2.A().Width();
00083       LocalHeap & lh = m2.GetLocalHeap();
00084       data = lh.Alloc<T> (h*w);
00085       CMCPMatExpr<FlatMatrix<T> >::operator= (m2.A());
00086     }
00087 
00089     template <typename T2>
00090     explicit FlatMatrix (const MatExpr<T2> & m)
00091       : h(m.Height()), w(m.Width()),
00092         data(const_cast<T*>(&m.Spec()(0,0)))  
00093     { ; }
00094   
00096     template <int H, int W>
00097     FlatMatrix (const Mat<H,W,TSCAL> & m) throw()
00098       : h(H), w(W), data(const_cast<T*>(&m(0,0)))
00099     { ; }
00100 
00102     ~FlatMatrix () throw() { ; }
00103 
00105     void AssignMemory (int ah, int aw, LocalHeap & lh)  // throw (LocalHeapOverflow)
00106     {
00107       h = ah;
00108       w = aw;
00109       data = (T*)lh.Alloc(h*w*sizeof(T));
00110     }
00111   
00113     void AssignMemory (int ah, int aw, T * mem) throw()
00114     {
00115       h = ah;
00116       w = aw;
00117       data = mem;
00118     }
00119   
00120 
00122     template<typename TBxx>
00123     const FlatMatrix & operator= (const Expr<TBxx> & m) const
00124     {
00125       return CMCPMatExpr<FlatMatrix>::operator= (m);
00126     }
00127 
00128     /*
00129     template <typename TA, typename TB>
00130     const FlatMatrix & operator= (const LapackProduct<TA, TB> & prod) const
00131     {
00132       LapackMult (prod.A(), prod.B(), *this);
00133       return *this;
00134     }
00135     */
00136 
00138     const FlatMatrix & operator= (const FlatMatrix & m) throw()
00139     {
00140       int hw = h*w;
00141       for (int i = 0; i < hw; i++)
00142         data[i] = m(i);
00143       return *this;
00144     }
00145 
00147     FlatMatrix & operator= (TSCAL s) throw()
00148     {
00149       int hw = h*w;
00150       for (int i = 0; i < hw; i++)
00151         data[i] = s; 
00152       return *this;
00153     }
00154 
00156     FlatMatrix & Assign (const FlatMatrix & m) throw()
00157     {
00158       h = m.h;
00159       w = m.w;
00160       data = m.data;
00161       return *this;
00162     }
00163 
00164 
00166     TELEM & operator() (int i) const 
00167     { 
00168 #ifdef CHECK_RANGE
00169       CheckMatRange(h,w,i);
00170 #endif
00171       return data[i]; 
00172     }
00173 
00175     TELEM & operator() (int i, int j) const
00176     {
00177 #ifdef CHECK_RANGE
00178       CheckMatRange(h,w,i,j);
00179 #endif
00180       return data[i*size_t(w)+j]; 
00181     }
00182 
00184     int Height () const throw() { return h; }
00185 
00187     int Width () const throw() { return w; }
00188 
00189     const FlatVector<T> Row (int i) const
00190     {
00191       return FlatVector<T> (w, &data[i*size_t(w)]);
00192     }
00193 
00194     const SliceVector<T> Col (int i) const
00195     {
00196       return SliceVector<T> (h, w, &data[i]);
00197     }
00198 
00199     const SliceVector<T> Diag () const
00200     {
00201       return SliceVector<T> (h, w+1, &data[0]);
00202     }
00203 
00204     using CMCPMatExpr<FlatMatrix<T> >::Rows;
00205     using CMCPMatExpr<FlatMatrix<T> >::Cols;
00206 
00207     const FlatMatrix Rows (int first, int next) const
00208     {
00209       return FlatMatrix (next-first, w, data+first*w);
00210     }
00211 
00212     const SliceMatrix<T> Cols (int first, int next) const
00213     {
00214       return SliceMatrix<T> (h, next-first, w, data+first);
00215     }
00216 
00217     const FlatMatrix Rows (IntRange range) const
00218     {
00219       return FlatMatrix (range.Next()-range.First(), w, data+range.First()*w);
00220     }
00221 
00222     const SliceMatrix<T> Cols (IntRange range) const
00223     {
00224       return SliceMatrix<T> (h, range.Next()-range.First(), w, data+range.First());
00225     }
00226   };
00227 
00228 
00229 
00230 
00231 
00235   template <typename T = double>
00236   class Matrix : public FlatMatrix<T>
00237   {
00238   public:
00239 
00241     typedef T TELEM;
00243     typedef typename mat_traits<T>::TSCAL TSCAL;
00244 
00246     Matrix () throw () : FlatMatrix<T> (0, 0) { ; }
00247 
00249     Matrix (int ah) : FlatMatrix<T> (ah, new T[size_t(ah)*size_t(ah)]) { ; }
00250     
00252     Matrix (int ah, int aw) : FlatMatrix<T> (ah, aw, new T[size_t(ah)*size_t(aw)]) { ; }
00253 
00255     Matrix (const Matrix & m2) 
00256       : FlatMatrix<T> (m2.Height(), m2.Width(), new T[m2.Height()*m2.Width()]) 
00257     {
00258       FlatMatrix<T>::operator= (m2);
00259     }
00260 
00262     template<typename TB>
00263     Matrix (const Expr<TB> & m2) 
00264       : FlatMatrix<T> (m2.Height(), m2.Width(), new T[m2.Height()*m2.Width()]) 
00265     {
00266       CMCPMatExpr<FlatMatrix<T> >::operator= (m2);
00267     }
00268 
00269 
00271     ~Matrix() { delete [] this->data; }
00272 
00274     void SetSize(int ah, int aw)
00275     {
00276       if (this->h == ah && this->w == aw) return;
00277       delete [] this->data;
00278       this->h = ah;
00279       this->w = aw;
00280       this->data = new T[this->h * this->w];
00281     }
00282 
00284     void SetSize(int ah)
00285     {
00286       if (this->h == ah && this->w == ah) return;
00287       delete [] this->data;
00288       this->h = ah;
00289       this->w = ah;
00290       this->data = new T[this->h * this->w];
00291     }
00292 
00294     template<typename TB>
00295     Matrix & operator= (const Expr<TB> & m) 
00296     { 
00297       SetSize (m.Height(), m.Width());
00298       CMCPMatExpr<FlatMatrix<T> >::operator= (m);
00299       return *this;
00300     }
00301 
00302     /*
00303     template <typename TA, typename TB>
00304     Matrix & operator= (const LapackProduct<TA, TB> & prod) 
00305     {
00306       // SetSize (m.Height(), m.Width());
00307       LapackMult (prod.A(), prod.B(), *this);
00308       return *this;
00309     }
00310     */
00311 
00313     Matrix & operator= (TSCAL s) 
00314     {
00315       FlatMatrix<T>::operator= (s);
00316       return *this;
00317     }
00318 
00319   };
00320 
00321 
00322 
00323 
00328   template <int H, int W = H, typename T = double>
00329   class Mat : public MatExpr<Mat<H,W,T> >
00330   {
00331     T data[(H*W>0) ? H*W : 1];
00332   public:
00333     typedef T TELEM;
00334     typedef typename mat_traits<T>::TSCAL TSCAL;
00335     typedef Vec<H, typename mat_traits<T>::TV_COL> TV_COL;
00336     typedef Vec<W, typename mat_traits<T>::TV_ROW> TV_ROW;
00337     enum { HEIGHT = H };
00338     enum { WIDTH  = W };
00339 
00341     Mat () throw () { ; }
00342 
00344     Mat (const Mat & m) throw()
00345       : MatExpr<Mat> ()
00346     {
00347       (*this) = m;
00348     }
00349 
00351     template<typename TB>
00352     Mat (const Expr<TB> & m)
00353     {
00354       MatExpr<Mat>::operator= (m);
00355     }
00356 
00358     Mat (TSCAL s) throw()     // explicit removed JS, Aug '07
00359     {
00360       for (int i = 0; i < H*W; i++)
00361         data[i] = s;
00362     }
00363 
00365     template<typename TB>
00366     Mat & operator= (const Expr<TB> & m)
00367     {
00368       MatExpr<Mat>::operator= (m);
00369       return *this;
00370     }
00371 
00373     Mat & operator= (const Mat & m) throw()
00374     {
00375       for (int i = 0; i < H*W; i++)
00376         data[i] = m.data[i];
00377       return *this;
00378     }
00379  
00381     Mat & operator= (TSCAL s) throw()
00382     {
00383       for (int i = 0; i < H*W; i++)
00384         data[i] = s;
00385       return *this;
00386     }
00387 
00389     TELEM & operator() (int i) { return data[i]; }
00391     TELEM & operator() (int i, int j) { return data[i*W+j]; }
00393     const TELEM & operator() (int i) const { return data[i]; }
00395     const TELEM & operator() (int i, int j) const { return data[i*W+j]; }
00396 
00398     int Height () const throw() { return H; }
00400     int Width () const throw() { return W; }
00401 
00403     const FlatVec<W,T> Row (int i) 
00404     {
00405       return FlatVec<W,T> (&(*this)(i,0));
00406     }
00407 
00409     const FlatVec<W,const T> Row (int i) const
00410     {
00411       return FlatVec<W,const T> (&(*this)(i,0));
00412     }
00413 
00415     const FixSliceVector<W, T> Col (int i) 
00416     {
00417       return FixSliceVector<W,T> (H, &(*this)(0,i));
00418     }
00419 
00420   };
00421 
00422 
00423   /*
00424   template <class TB> ALWAYS_INLINE
00425   Mat<3,3,double> operator+= (class Mat<3,3,double> & m, const Expr<TB> & exp)
00426   {
00427     m(0,0) += exp.Spec()(0,0);
00428     m(0,1) += exp.Spec()(0,1);
00429     m(0,2) += exp.Spec()(0,2);
00430     m(1,0) += exp.Spec()(1,0);
00431     m(1,1) += exp.Spec()(1,1);
00432     m(1,2) += exp.Spec()(1,2);
00433     m(2,0) += exp.Spec()(2,0);
00434     m(2,1) += exp.Spec()(2,1);
00435     m(2,2) += exp.Spec()(2,2);
00436     return m;
00437   }
00438 
00439   template <class TB> ALWAYS_INLINE
00440   Mat<2,2,double> operator+= (class Mat<2,2,double> & m, const Expr<TB> & exp)
00441   {
00442     m(0,0) += exp.Spec()(0,0);
00443     m(0,1) += exp.Spec()(0,1);
00444     m(1,0) += exp.Spec()(1,0);
00445     m(1,1) += exp.Spec()(1,1);
00446     return m;
00447   }
00448   */
00449 
00450 
00451 
00452 
00453 
00454 
00455 
00456 
00457 
00458 
00459 
00464   template <int W, typename T = double>
00465   class FlatMatrixFixWidth : public CMCPMatExpr<FlatMatrixFixWidth<W,T> >
00466   {
00467   protected:
00469     T * data;
00471     int h;
00472   public:
00474     typedef T TELEM;
00475     typedef T& TREF;
00477     typedef typename mat_traits<T>::TSCAL TSCAL;
00478 
00480     FlatMatrixFixWidth () throw() { ; }
00481   
00483     FlatMatrixFixWidth (int ah, T * adata) throw()
00484       : data(adata), h(ah) { ; }
00485 
00487     FlatMatrixFixWidth (int ah, LocalHeap & lh) // throw (LocalHeapOverflow)
00488       : data((T*)lh.Alloc(ah*W*sizeof(T))), h(ah) { ; }
00489   
00491     FlatMatrixFixWidth (const FlatMatrixFixWidth & m) throw()
00492       : CMCPMatExpr<FlatMatrixFixWidth> (), data(m.data), h(m.h)
00493     { ; }
00494 
00496     FlatMatrixFixWidth (const FlatMatrix<TSCAL> & m) throw()
00497       : data(const_cast<T*> (&m(0))), h(m.Height())
00498     { ; }
00499 
00500     template <int H>
00501     FlatMatrixFixWidth (const Mat<H,W,TSCAL> & m) throw()
00502       : data(const_cast<T*>(&m(0,0))), h(H)
00503     { ; }  
00504 
00506     ~FlatMatrixFixWidth () throw() { ; }
00507 
00509     void AssignMemory (int ah, LocalHeap & lh) // throw (LocalHeapOverflow)
00510     {
00511       h = ah;
00512       data = (T*)lh.Alloc(h*W*sizeof(T));
00513     }
00514   
00516     void AssignMemory (int ah, T * mem) throw()
00517     {
00518       h = ah;
00519       data = mem;
00520     }
00521   
00522 
00524     template<typename TB>
00525     const FlatMatrixFixWidth & operator= (const Expr<TB> & m) const
00526     {
00527       return CMCPMatExpr<FlatMatrixFixWidth>::operator= (m);
00528     }
00529 
00531     FlatMatrixFixWidth & operator= (const FlatMatrixFixWidth & m) throw()
00532     {
00533       for (int i = 0; i < h*W; i++)
00534         data[i] = m(i);
00535       return *this;
00536     }
00537 
00539     FlatMatrixFixWidth & operator= (TSCAL s) throw()
00540     {
00541       for (int i = 0; i < h*W; i++)
00542         data[i] = s; 
00543       return *this;
00544     }
00545 
00547     FlatMatrixFixWidth & Assign (const FlatMatrixFixWidth & m) throw()
00548     {
00549       h = m.h;
00550       data = m.data;
00551       return *this;
00552     }
00553 
00555     TELEM & operator() (int i) const
00556     {
00557 #ifdef CHECK_RANGE
00558       CheckMatRange(h,W,i);
00559 #endif
00560       return data[i]; 
00561     }
00562 
00564     TELEM & operator() (int i, int j) const
00565     {
00566 #ifdef CHECK_RANGE
00567       CheckMatRange(h,W,i,j);
00568 #endif
00569       return data[i*W+j]; 
00570     }
00571 
00573     int Height () const throw() { return h; }
00574 
00576     int Width () const throw() { return W; }
00577 
00579     operator const FlatMatrix<T>() const { return FlatMatrix<T> (h, W, data); }
00580 
00582     const FlatVec<W,T> Row (int i) const
00583     {
00584       return FlatVec<W,T> (&(*this)(i,0));
00585     }
00586 
00587     const FixSliceVector<W,T> Col (int i)
00588     {
00589       return FixSliceVector<W,T> (h, &data[i]);
00590     }
00591 
00592     using CMCPMatExpr<FlatMatrixFixWidth<W,T> >::Rows;
00593     using CMCPMatExpr<FlatMatrixFixWidth<W,T> >::Cols;
00594 
00595     const FlatMatrixFixWidth Rows (int first, int next) const
00596     {
00597       return FlatMatrixFixWidth (next-first, data+first*W);
00598     }
00599 
00600     const SliceMatrix<T> Cols (int first, int next) const
00601     {
00602       return SliceMatrix<T> (h, next-first, W, data+first);
00603     }
00604 
00605     const FlatMatrixFixWidth Rows (IntRange range) const
00606     {
00607       return Rows (range.First(), range.Next());
00608     }
00609 
00610     const SliceMatrix<T> Cols (IntRange range) const
00611     {
00612       return Cols (range.First(), range.Next());
00613     }
00614   };
00615 
00616 
00620   template <int W, typename T = double>
00621   class MatrixFixWidth : public FlatMatrixFixWidth<W,T>
00622   {
00623   public:
00624     typedef typename mat_traits<T>::TSCAL TSCAL;
00625 
00626     MatrixFixWidth () { ; }
00627 
00629     MatrixFixWidth (int ah) : FlatMatrixFixWidth<W,T> (ah, new T[ah*W]) { ; }
00630 
00632     ~MatrixFixWidth() { delete [] this->data; }
00633 
00635     void SetSize(int ah)
00636     {
00637       if (this->h == ah) return;
00638       delete [] this->data;
00639       this->h = ah;
00640       this->data = new T[this->h * W];
00641     }
00642 
00644     template<typename TB>
00645     MatrixFixWidth & operator= (const Expr<TB> & m)
00646     {
00647       MatExpr<FlatMatrixFixWidth<W,T> >::operator= (m);
00648       return *this;
00649     }
00650 
00652     MatrixFixWidth & operator= (TSCAL s)
00653     {
00654       FlatMatrixFixWidth<W,T>::operator= (s);
00655       return *this;
00656     }
00657   };
00658 
00659 
00660 
00661 
00662 
00663 
00664 
00665 
00666 
00667 
00668 
00669 
00670 
00671 
00672 
00673 
00674 
00675 
00676 
00677 
00678 
00679 
00685   template <int H, typename T = double>
00686   class FlatMatrixFixHeight : public MatExpr<FlatMatrixFixHeight<H,T> >
00687   {
00688   protected:
00690     T * data;
00692     int w;
00693   public:
00695     typedef T TELEM;
00697     typedef typename mat_traits<T>::TSCAL TSCAL;
00698 
00700     FlatMatrixFixHeight () throw()
00701       : data(0), w(0) { ; }
00702   
00704     FlatMatrixFixHeight (int aw, T * adata) throw()
00705       : data(adata), w(aw) { ; }
00706 
00708     FlatMatrixFixHeight (int aw, LocalHeap & lh) // throw (LocalHeapOverflow)
00709       : data((T*)lh.Alloc(aw*H*sizeof(T))), w(aw) { ; }
00710   
00712     FlatMatrixFixHeight (const FlatMatrixFixHeight & m)
00713       : data(m.data), w(m.w)
00714     { ; }
00715   
00716 
00718     ~FlatMatrixFixHeight () throw() { ; }
00719 
00721     void AssignMemory (int aw, LocalHeap & lh) // throw (LocalHeapOverflow)
00722     {
00723       w = aw;
00724       data = (T*)lh.Alloc(w*H*sizeof(T));
00725     }
00726   
00728     void AssignMemory (int aw, T * mem) 
00729     {
00730       w = aw;
00731       data = mem;
00732     }
00733   
00734 
00736     template<typename TB>
00737     FlatMatrixFixHeight & operator= (const Expr<TB> & m)
00738     {
00739       for (int j = 0; j < w; j++)
00740         for (int i = 0; i < H; i++)
00741           (*this)(i,j) = m.Spec()(i,j);
00742       return *this;
00743     }
00744 
00746     FlatMatrixFixHeight & operator= (const FlatMatrixFixHeight & m)
00747     {
00748       for (int i = 0; i < w*H; i++)
00749         data[i] = m(i);
00750       return *this;
00751     }
00752 
00754     FlatMatrixFixHeight & operator= (TSCAL s)
00755     {
00756       for (int i = 0; i < w*H; i++)
00757         data[i] = s; 
00758       return *this;
00759     }
00760 
00762     FlatMatrixFixHeight & Assign (const FlatMatrixFixHeight & m)
00763     {
00764       w = m.w;
00765       data = m.data;
00766       return *this;
00767     }
00768 
00769 
00771     TELEM & operator() (int i)
00772     { 
00773 #ifdef CHECK_RANGE
00774       CheckMatRange(H,w,i);
00775 #endif
00776       return data[i]; 
00777     }
00778 
00780     TELEM & operator() (int i, int j) 
00781     {
00782 #ifdef CHECK_RANGE
00783       CheckMatRange(H,w,i,j);
00784 #endif
00785       return data[i+j*H]; 
00786     }
00787 
00788 
00790     const TELEM & operator() (int i) const
00791     {
00792 #ifdef CHECK_RANGE
00793       CheckMatRange(H,w,i);
00794 #endif
00795       return data[i]; 
00796     }
00797 
00799     const TELEM & operator() (int i, int j) const
00800     {
00801 #ifdef CHECK_RANGE
00802       CheckMatRange(H,w,i,j);
00803 #endif
00804       return data[i+j*H]; 
00805     }
00806 
00807 
00809     int Height () const { return H; }
00810 
00812     int Width () const { return w; }
00813 
00814 
00815     const FlatVec<H,T> Col (int i) const
00816     {
00817       return FlatVec<H,T> ( data+i*H ); 
00818     }
00819 
00820     const SliceVector<T> Row (int i) const
00821     {
00822       return SliceVector<T> (w, H, &data[i]);
00823     }
00824     
00825     /*
00826     const DoubleSliceMatrix<T> Rows (int first, int next) const
00827     {
00828       return DoubleSliceMatrix<T> (next-first, w, 1, H, data+first);
00829     }
00830 
00831     const FlatMatrixFixHeight<H,T> Cols (int first, int next) const
00832     {
00833       return FlatMatrixFixHeight<H,T> (next-first, data+first*H);
00834     }
00835 
00836     const DoubleSliceMatrix<T> Rows (IntRange range) const
00837     {
00838       return Rows (range.First(), range.Next());
00839     }
00840 
00841     const FlatMatrixFixHeight<H,T> Cols (IntRange range) const
00842     {
00843       return Cols (range.First(), range.Next());
00844     }
00845     */
00846 
00847 
00848     const FlatMatrixFixWidth<H,T> Trans () const
00849     {
00850       return FlatMatrixFixWidth<H,T> (w, data);
00851     }
00852 
00853 
00854 
00855     enum { IS_LINEAR = 0 };
00856     //  static bool IsLinear() { return false; }
00857   };
00858 
00859 
00863   template <int H, typename T = double>
00864   class MatrixFixHeight : public FlatMatrixFixHeight<H,T>
00865   {
00866   public:
00867     typedef typename mat_traits<T>::TSCAL TSCAL;
00868 
00870     MatrixFixHeight (int aw) : FlatMatrixFixHeight<H,T> (aw, new T[aw*H]) { ; }
00871 
00873     ~MatrixFixHeight() { delete [] this->data; }
00874 
00876     void SetSize(int aw)
00877     {
00878       if (this->w == aw) return;
00879       delete [] this->data;
00880       this->w = aw;
00881       this->data = new T[H*this->w];
00882     }
00883 
00885     template<typename TB>
00886     MatrixFixHeight & operator= (const Expr<TB> & m)
00887     {
00888       int k = 0;
00889       for (int j = 0; j < this->w; j++)
00890         for (int i = 0; i < H; i++)
00891           {
00892             this->data[k] = m.Spec()(i,j);
00893             k++;
00894           }
00895       return *this;
00896     }
00897 
00899     MatrixFixHeight & operator= (TSCAL s)
00900     {
00901       FlatMatrixFixHeight<H,T>::operator= (s);
00902       return *this;
00903     }
00904   };
00905 
00906 
00907 
00908 
00909 
00910 
00911 
00912 
00913 
00914 
00915 
00916 
00917   template <typename T = double>
00918   class SliceMatrix : public CMCPMatExpr<SliceMatrix<T> >
00919   {
00920   protected:
00922     int h;
00924     int w;
00926     int dist;
00928     T * data;
00929   public:
00930 
00932     typedef T TELEM;
00934     typedef typename mat_traits<T>::TSCAL TSCAL;
00935     enum { IS_LINEAR = 0 };
00936 
00938     SliceMatrix (int ah, int aw, int adist, T * adata) throw ()
00939       : h(ah), w(aw), dist(adist), data(adata) { ; }
00940   
00941     SliceMatrix (const FlatMatrix<T> & mat)
00942       : h(mat.Height()), w(mat.Width()), dist(mat.Width()), data(&mat(0,0))
00943     { ; }
00944 
00946     template<typename TB>
00947     const SliceMatrix & operator= (const Expr<TB> & m) const
00948     {
00949       return CMCPMatExpr<SliceMatrix>::operator= (m);
00950     }
00951 
00952     /*
00953     template <typename TA, typename TB>
00954     const SliceMatrix & operator= (const LapackProduct<TA, TB> & prod) const
00955     {
00956       LapackMult (prod.A(), prod.B(), *this);
00957       return *this;
00958     }
00959     */
00960 
00962     SliceMatrix & operator= (TSCAL s) throw()
00963     {
00964       for (int i = 0; i < h; i++)
00965         for (int j = 0; j < w; j++)
00966           data[i*dist+j] = s;
00967       return *this;
00968     }
00969 
00971     TELEM & operator() (int i, int j) const
00972     {
00973 #ifdef CHECK_RANGE
00974       CheckMatRange(h,w,i,j);
00975 #endif
00976       return data[i*dist+j]; 
00977     }
00978 
00980     TELEM & operator() (int i) const
00981     {
00982 #ifdef CHECK_RANGE
00983       CheckMatRange(h,w,i);
00984 #endif
00985       return data[i]; 
00986     }
00987 
00989     int Height () const throw() { return h; }
00990 
00992     int Width () const throw() { return w; }
00993 
00995     int Dist () const throw() { return dist; }
00996 
00997     const SliceMatrix Rows (int first, int next) const
00998     {
00999       return SliceMatrix (next-first, w, dist, data+first*dist);
01000     }
01001 
01002     const FlatVector<T> Row (int i) const
01003     {
01004       return FlatVector<T> (w, &data[i*size_t(dist)]);
01005     }
01006 
01007     const SliceVector<T> Col (int i) const
01008     {
01009       return SliceVector<T> (h, dist, &data[i]);
01010     }
01011 
01012     const SliceMatrix<T> Cols (int first, int next) const
01013     {
01014       return SliceMatrix<T> (h, next-first, dist, data+first);
01015     }
01016 
01017     const SliceMatrix Rows (IntRange range) const
01018     {
01019       return Rows (range.First(), range.Next());
01020     }
01021 
01022     const SliceMatrix<T> Cols (IntRange range) const
01023     {
01024       return Cols (range.First(), range.Next());
01025     }
01026 
01027   };
01028 
01029 
01030 
01031 
01032 
01033   template <typename T = double>
01034   class DoubleSliceMatrix : public CMCPMatExpr<DoubleSliceMatrix<T> >
01035   {
01036   protected:
01038     int h;
01040     int w;
01042     int distr, distc;
01044     T * data;
01045   public:
01046 
01048     typedef T TELEM;
01050     typedef typename mat_traits<T>::TSCAL TSCAL;
01051     enum { IS_LINEAR = 0 };
01052 
01054     DoubleSliceMatrix (int ah, int aw, int adistr, int adistc, T * adata) throw ()
01055       : h(ah), w(aw), distr(adistr), distc(adistc), data(adata) { ; }
01056   
01058     template<typename TB>
01059     const DoubleSliceMatrix & operator= (const Expr<TB> & m) const
01060     {
01061       return CMCPMatExpr<DoubleSliceMatrix>::operator= (m);
01062     }
01063 
01065     DoubleSliceMatrix & operator= (TSCAL s) throw()
01066     {
01067       for (int i = 0; i < h; i++)
01068         for (int j = 0; j < w; j++)
01069           data[i*distr+j*distc] = s;
01070       return *this;
01071     }
01072 
01074     TELEM & operator() (int i, int j) const
01075     {
01076 #ifdef CHECK_RANGE
01077       CheckMatRange(h,w,i,j);
01078 #endif
01079       return data[i*distr+j*distc]; 
01080     }
01081 
01083     TELEM & operator() (int i) const
01084     {
01085 #ifdef CHECK_RANGE
01086       CheckMatRange(h,w,i);
01087 #endif
01088       return data[i]; 
01089     }
01090 
01092     int Height () const throw() { return h; }
01093 
01095     int Width () const throw() { return w; }
01096 
01097     /*
01098     const DoubleSliceMatrix Rows (int first, int next) const
01099     {
01100       return DoubleSliceMatrix (next-first, w, distr, distc, data+first*distr);
01101     }
01102 
01103     const DoubleSliceMatrix<T> Cols (int first, int next) const
01104     {
01105       return DoubleSliceMatrix<T> (h, next-first, distr, distc, data+first*distc);
01106     }
01107 
01108     const DoubleSliceMatrix Rows (IntRange range) const
01109     {
01110       return Rows(range.First(), range.Next());
01111     }
01112 
01113     const DoubleSliceMatrix<T> Cols (IntRange range) const
01114     {
01115       return Cols(range.First(), range.Next());
01116     }
01117     */
01118 
01119   };
01120 
01121 
01122 
01123 
01124 
01125 
01126 
01127 
01128 
01129 
01130 
01131 
01132 
01133 
01134 
01135 
01136   template <class TM, class TSCAL>
01137   class Scalar2ElemMatrix 
01138   {
01139   public:
01140     const FlatMatrix<TSCAL> & mat;
01141     Scalar2ElemMatrix (const FlatMatrix<TSCAL> & amat) : mat(amat) { ; }
01142 
01143     enum { H = mat_traits<TM>::HEIGHT };
01144     enum { W = mat_traits<TM>::WIDTH };
01145 
01146     TM operator() (int i, int j) const
01147     {
01148       TM ret;
01149       for (int k = 0; k < H; k++)
01150         for (int l = 0; l < W; l++)
01151           Access(ret, k,l) = mat(i*H+k, j*W+l);
01152       return ret;
01153     }
01154   };
01155 
01156 
01157 
01158 
01159 
01160 
01164   template <int H>
01165   class Id : public MatExpr<Id<H> >
01166   {
01167   public:
01168     typedef double TELEM;
01169     typedef double TSCAL;
01170     typedef Vec<H, double> TV_COL;
01171     typedef Vec<H, double> TV_ROW;
01172     enum { HEIGHT = H };
01173     enum { WIDTH  = H };
01174     // static bool IsLinear() { return 0; }
01175     enum { IS_LINEAR = 0 };
01176 
01178     Id () { ; }
01179 
01181     double operator() (int i) const
01182     { cerr << "id, linear access" << endl; return 0; }
01184     double operator() (int i, int j) const { return (i == j) ? 1 : 0; }
01185 
01187     int Height () const { return H; }
01189     int Width () const { return H; }
01190   };
01191 
01192 
01193 
01194 
01196   class Identity : public MatExpr<Identity >
01197   {
01198     int size;
01199   public:
01200     typedef double TELEM;
01201     typedef double TSCAL;
01202     enum { IS_LINEAR = 0 };
01203 
01204     Identity (int s) : size(s) { ; }
01205 
01206     double operator() (int i) const
01207     { cerr << "Identity, linear access" << endl; return 0; }
01208 
01209     double operator() (int i, int j) const { return (i == j) ? 1 : 0; }
01210 
01211     int Height () const { return size; }
01212     int Width () const { return size; }
01213   };
01214 
01215 
01216 
01217 
01218 
01219   template<int H, int W, typename T>
01220   inline std::ostream & operator<< (std::ostream & s, const Mat<H,W,T> & m)
01221   {
01222     for (int i = 0; i < H*W; i++)
01223       s << " " << setw(7) << m(i);
01224     return s;
01225   }
01226 
01227 
01228 }
01229 
01230 
01231 #ifdef PARALLEL
01232 namespace ngstd
01233 {
01234   template<int N, int M, typename T>
01235   class MPI_Traits<ngbla::Mat<N, M, T> >
01236   {
01237   public:
01238     static MPI_Datatype MPIType () 
01239     { 
01240       static MPI_Datatype MPI_T = 0;
01241       if (!MPI_T)
01242         {
01243           int size = N * M;
01244           MPI_Type_contiguous ( size, MPI_Traits<T>::MPIType(), &MPI_T);
01245           MPI_Type_commit ( &MPI_T );
01246         }
01247       return MPI_T;
01248     }
01249   };
01250   
01251 
01252 }
01253 #endif
01254 
01255 #endif