NGSolve
4.9
|
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