NGSolve
4.9
|
00001 #ifndef FILE_VECTOR_EXPR 00002 #define FILE_VECTOR_EXPR 00003 00004 /**************************************************************************/ 00005 /* File: vector.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 01. Jan. 02 */ 00008 /**************************************************************************/ 00009 00010 00011 namespace ngbla 00012 { 00013 00014 00015 template <int S, class T> class Vec; 00016 template <int S, typename T> class FlatVec; 00017 template <class T> class SysVector; 00018 template <class T> class FlatVector; 00019 template <class T> class Vector; 00020 template <class T> class SliceVector; 00021 00022 00023 extern void CheckVecRange(int s, int i); 00024 extern void CheckVecRange(int s, int i, int j); 00025 00026 00032 template <typename T = double> 00033 class FlatVector : public CMCPMatExpr<FlatVector<T> > 00034 { 00035 protected: 00037 int s; 00039 T * data; 00040 public: 00042 typedef T TELEM; 00043 typedef T& TREF; 00045 typedef typename mat_traits<T>::TSCAL TSCAL; 00046 00048 FlatVector () { ; } 00050 FlatVector (unsigned int as, T * adata) : s(as), data(adata) { ; } 00051 00053 FlatVector (unsigned int as, void * adata) : s(as), data(static_cast<TELEM*> (adata)) { ; } 00054 00056 template <int S> 00057 FlatVector (const Vec<S,TSCAL> & v) 00058 : s(v.Size()), data(const_cast<T*>(&v(0))) 00059 { ; } 00060 00061 template <int S> 00062 FlatVector (const FlatVec<S,TSCAL> & v) 00063 : s(v.Size()), data(const_cast<T*>(&v(0))) 00064 { ; } 00065 00067 FlatVector (int as, LocalHeap & lh) 00068 : s(as), data((T*)lh.Alloc(s*sizeof(T))) { ; } 00069 00071 FlatVector (const SysVector<TSCAL> & sv) 00072 : s(sv.Size()*sv.BlockSize() / mat_traits<T>::VDIM), 00073 data (sv(0)) 00074 { 00075 ; 00076 } 00077 00079 void AssignMemory (int as, LocalHeap & lh) 00080 { 00081 s = as; 00082 data = (T*)lh.Alloc(s*sizeof(T)); 00083 } 00084 00086 void AssignMemory (int as, T * mem) 00087 { 00088 s = as; 00089 data = mem; 00090 } 00091 00093 const FlatVector & operator= (const FlatVector & v) const 00094 { 00095 for (int i = 0; i < s; i++) 00096 data[i] = v(i); 00097 return *this; 00098 } 00099 00101 template<typename TB> 00102 const FlatVector & operator= (const Expr<TB> & v) const 00103 { 00104 return CMCPMatExpr<FlatVector>::operator= (v); 00105 } 00106 00108 const FlatVector & operator= (TSCAL scal) const 00109 { 00110 for (int i = 0; i < s; i++) 00111 data[i] = scal; 00112 return *this; 00113 } 00114 00115 template<typename TB> 00116 const FlatVector & operator+= (const Expr<TB> & v) const 00117 { 00118 if (TB::IS_LINEAR) 00119 for (int i = 0; i < s; i++) 00120 data[i] += v.Spec()(i); 00121 else 00122 for (int i = 0; i < s; i++) 00123 data[i] += v.Spec()(i,0); 00124 return *this; 00125 } 00126 00128 TELEM & operator() (int i) const 00129 { 00130 #ifdef CHECK_RANGE 00131 CheckVecRange(s,i); 00132 #endif 00133 return data[i]; 00134 } 00135 00136 RowsArrayExpr<FlatVector> operator() (FlatArray<int> rows) const 00137 { 00138 return RowsArrayExpr<FlatVector> (*this, rows); 00139 } 00140 00141 00143 TELEM & operator() (int i, int j) const 00144 { 00145 #ifdef CHECK_RANGE 00146 CheckVecRange(s,i); 00147 #endif 00148 return data[i]; 00149 } 00150 00152 TELEM & operator[] (int i) const 00153 { 00154 #ifdef CHECK_RANGE 00155 CheckVecRange(s,i); 00156 #endif 00157 return data[i]; 00158 } 00159 00160 // shape functions had a problem with icc v9.1 00161 const CArray<T> Addr(int i) const 00162 { 00163 return CArray<T> (data+i); 00164 } 00165 00166 /* 00167 T * const Addr (int i) const // const not respected by icc ??? 00168 { 00169 return data+i; 00170 } 00171 */ 00172 00174 const FlatVector<T> Range (int first, int next) const 00175 { return FlatVector<T> (next-first, data+first); } 00176 00178 const FlatVector<T> Range (IntRange range) const 00179 { return FlatVector<T> (range.Next()-range.First(), data+range.First()); } 00180 00181 00183 int Size () const { return s; } 00184 00186 int Height () const { return s; } 00187 00189 int Width () const { return 1; } 00190 00192 SliceVector<T> Slice (int first, int dist) 00193 { 00194 return SliceVector<T> (s/dist, dist, data+first); 00195 } 00196 00198 const SliceVector<T> Slice (int first, int dist) const 00199 { 00200 return SliceVector<T> (s/dist, dist, data+first); 00201 } 00202 00204 const void * Data () const { return static_cast<const void*>(data); } 00206 void * Data () { return static_cast<void*>(data); } 00207 00208 00209 // new for SysVectors: 00210 typedef FlatVector<T> TV_COL; 00211 typedef double TV_ROW; 00212 enum { HEIGHT = 1 }; 00213 enum { WIDTH = 1 }; 00214 }; 00215 00216 00217 00218 00219 00220 00221 00222 00223 template <int S, typename T> 00224 class FlatVector<Vec<S, T> > : public CMCPMatExpr<FlatVector<Vec<S, T> > > 00225 { 00226 protected: 00228 int s; 00230 T * data; 00231 public: 00233 typedef Vec<S,T> TELEM; 00234 typedef FlatVec<S,T> TREF; 00236 typedef typename mat_traits<T>::TSCAL TSCAL; 00237 00239 FlatVector () { ; } 00241 FlatVector (const FlatVector & v2) : s(v2.s), data(v2.data) { ; } 00243 FlatVector (unsigned int as, T * adata) : s(as), data(adata) { ; } 00244 00245 // set size and mem 00246 // FlatVector (unsigned int as, TELEM * adata) : s(as), data(&(*adata)(0)) { ; } 00247 00249 FlatVector (unsigned int as, void * adata) : s(as), data(static_cast<T*> (adata)) { ; } 00250 00251 /* 00253 template <int S> 00254 FlatVector (const Vec<S,TSCAL> & v) 00255 : s(v.Size()), data(const_cast<T*>(&v(0))) 00256 { ; } 00257 */ 00258 00260 FlatVector (int as, LocalHeap & lh) 00261 : s(as), data((T*)lh.Alloc(s*S*sizeof(T))) { ; } 00262 00263 /* 00265 FlatVector (const SysVector<TSCAL> & sv) 00266 : s(sv.Size()*sv.BlockSize() / mat_traits<T>::VDIM), 00267 data (sv(0)) 00268 { 00269 ; 00270 } 00271 */ 00272 00274 void AssignMemory (int as, LocalHeap & lh) 00275 { 00276 s = as; 00277 data = (T*)lh.Alloc(s*S*sizeof(T)); 00278 } 00279 00281 void AssignMemory (int as, T * mem) 00282 { 00283 s = as; 00284 data = mem; 00285 } 00286 00288 const FlatVector & operator= (const FlatVector & v) const 00289 { 00290 for (int i = 0; i < s; i++) 00291 (*this)(i) = v(i); 00292 return *this; 00293 } 00294 00296 template<typename TB> 00297 const FlatVector & operator= (const Expr<TB> & v) const 00298 { 00299 return CMCPMatExpr<FlatVector>::operator= (v); 00300 } 00301 00303 const FlatVector & operator= (TSCAL scal) const 00304 { 00305 for (int i = 0; i < s; i++) 00306 (*this)(i) = scal; 00307 return *this; 00308 } 00309 00310 template<typename TB> 00311 const FlatVector & operator+= (const Expr<TB> & v) const 00312 { 00313 if (TB::IS_LINEAR) 00314 for (int i = 0; i < s; i++) 00315 (*this)(i) += v.Spec()(i); 00316 else 00317 for (int i = 0; i < s; i++) 00318 (*this)(i) += v.Spec()(i,0); 00319 return *this; 00320 } 00321 00323 const FlatVec<S,T> operator() (int i) const 00324 { 00325 return FlatVec<S,T> (data+i*S); 00326 } 00327 00329 const FlatVec<S,T> operator() (int i, int j) const 00330 { 00331 return FlatVec<S,T> (data+i*S); 00332 } 00333 00335 const FlatVec<S,T> operator[] (int i) const 00336 { 00337 return FlatVec<S,T> (data+i*S); 00338 } 00339 00340 RowsArrayExpr<FlatVector> operator() (FlatArray<int> rows) const 00341 { 00342 return RowsArrayExpr<FlatVector> (*this, rows); 00343 } 00344 00345 00346 // shape functions had a problem with icc v9.1 00347 00348 const CArray<Vec<S,T> > Addr(int i) const 00349 { 00350 return CArray<Vec<S,T> > (static_cast<Vec<S,T>*> ((void*) (data+i*S))); 00351 } 00352 /* 00353 const CArray<T> Addr(int i) const 00354 { 00355 return CArray<T> (data+i*S); 00356 } 00357 */ 00358 00360 const FlatVector<Vec<S,T> > Range (int first, int next) const 00361 { return FlatVector<Vec<S,T> > (next-first, data+S*first); } 00362 00364 const FlatVector<Vec<S,T> > Range (IntRange range) const 00365 { return FlatVector<Vec<S,T> > (range.Next()-range.First(), data+S*range.First()); } 00366 00368 int Size () const { return s; } 00369 00371 int Height () const { return s; } 00372 00374 int Width () const { return 1; } 00375 00376 /* 00378 SliceVector<T> Slice (int first, int dist) 00379 { 00380 return SliceVector<T> (s/dist, dist, data+first); 00381 } 00382 00384 const SliceVector<T> Slice (int first, int dist) const 00385 { 00386 return SliceVector<T> (s/dist, dist, data+first); 00387 } 00388 */ 00389 00391 const void * Data () const { return static_cast<const void*>(data); } 00393 void * Data () { return static_cast<void*>(data); } 00394 00395 // new for SysVectors: 00396 // typedef FlatVector<T> TV_COL; 00397 // typedef double TV_ROW; 00398 enum { HEIGHT = 1 }; 00399 enum { WIDTH = 1 }; 00400 }; 00401 00402 00403 00404 00405 00406 00407 00411 template <typename T = double> 00412 class Vector : public FlatVector<T> 00413 { 00414 public: 00415 typedef typename mat_traits<T>::TSCAL TSCAL; 00416 00418 Vector () : FlatVector<T> (0, (T*)0) { ; } 00419 00421 explicit Vector (int as) : FlatVector<T> (as, new T[as]) { ; } 00422 00423 00425 Vector (const Vector & v2) 00426 : FlatVector<T> (v2.Size(), new T[v2.Size()]) 00427 { 00428 FlatVector<T>::operator= (v2); 00429 } 00430 00432 template<typename TB> 00433 Vector (const Expr<TB> & v2) 00434 : FlatVector<T> (v2.Height(), new T[v2.Height()]) 00435 { 00436 FlatVector<T>::operator= (v2); 00437 } 00438 00439 00441 ~Vector() { delete [] this->data; } 00442 00444 Vector & operator= (TSCAL scal) 00445 { 00446 FlatVector<T>::operator= (scal); 00447 return *this; 00448 } 00449 00451 void SetSize(int as) 00452 { 00453 if (this->s == as) return; 00454 delete [] this->data; 00455 this->s = as; 00456 this->data = new T[this->s]; 00457 } 00458 00460 template<typename TB> 00461 Vector & operator= (const Expr<TB> & v) 00462 { 00463 MatExpr<FlatVector<T> >::operator= (v); 00464 return *this; 00465 } 00466 00467 Vector & operator= (const Vector & v2) 00468 { 00469 FlatVector<T>::operator= (static_cast<FlatVector<T> >(v2)); 00470 // MatExpr<FlatVector<T> >::operator= (v2); // does not work, we don't know why 00471 return *this; 00472 } 00473 00474 }; 00475 00476 00477 00478 00479 00480 template <int S, typename T> 00481 class Vector<Vec<S,T> > : public FlatVector<Vec<S,T> > 00482 { 00483 public: 00484 typedef typename mat_traits<T>::TSCAL TSCAL; 00485 00487 Vector () : FlatVector<Vec<S,T> > (0, (T*)0) { ; } 00488 00490 explicit Vector (int as) : FlatVector<Vec<S,T> > (as, new T[as*S]) { ; } 00491 00492 00494 Vector (const Vector & v2) 00495 : FlatVector<Vec<S,T> > (v2.Size(), new T[S*v2.Size()]) 00496 { 00497 FlatVector<Vec<S,T> >::operator= (v2); 00498 } 00499 00501 template<typename TB> 00502 Vector (const Expr<TB> & v2) 00503 : FlatVector<Vec<S,T> > (v2.Height(), new T[S*v2.Height()]) 00504 { 00505 FlatVector<Vec<S,T> >::operator= (v2); 00506 } 00507 00509 ~Vector() { delete [] this->data; } 00510 00512 Vector & operator= (TSCAL scal) 00513 { 00514 FlatVector<Vec<S,T> >::operator= (scal); 00515 return *this; 00516 } 00517 00519 void SetSize(int as) 00520 { 00521 if (this->s == as) return; 00522 delete [] this->data; 00523 this->s = as; 00524 this->data = new T[S*this->s]; 00525 } 00526 00528 template<typename TB> 00529 Vector & operator= (const Expr<TB> & v) 00530 { 00531 MatExpr<FlatVector<Vec<S,T> > >::operator= (v); 00532 return *this; 00533 } 00534 00535 Vector & operator= (const Vector & v2) 00536 { 00537 FlatVector<Vec<S,T> >::operator= (static_cast<FlatVector<Vec<S,T> > >(v2)); 00538 // MatExpr<FlatVector<Vec<S,T> > >::operator= (v2); // does not work, we don't know why 00539 return *this; 00540 } 00541 00542 }; 00543 00544 00545 00546 00547 00548 00549 00556 template <int S, typename T = double> 00557 class VectorMem : public FlatVector<T> 00558 { 00560 double mem[(S*sizeof(T)+7) / 8]; // alignment (on ia64 machines) 00561 // T mem[S]; // should be best, but calls trivial default constructor 00562 public: 00564 typedef typename mat_traits<T>::TSCAL TSCAL; 00565 00569 explicit VectorMem (int as) : FlatVector<T> (as, (as <= S) ? 00570 static_cast<void*> (&mem[0]) : 00571 static_cast<void*> (new T[as])) { ; } 00572 00574 ~VectorMem() { if (this->s > S) delete [] this->data; } 00575 00577 VectorMem & operator= (TSCAL scal) 00578 { 00579 FlatVector<T>::operator= (scal); 00580 return *this; 00581 } 00582 00584 template<typename TB> 00585 VectorMem & operator= (const Expr<TB> & v) 00586 { 00587 MatExpr<FlatVector<T> >::operator= (v); 00588 return *this; 00589 } 00590 }; 00591 00592 00593 00594 00595 00596 00597 00598 00599 00600 00601 00602 00603 00604 00605 00606 // A system vector (not completely functional) 00607 template <typename T = double> 00608 class FlatSysVector : public MatExpr<FlatSysVector<T> > 00609 { 00610 protected: 00611 int s; 00612 int blocksize; 00613 T * data; 00614 public: 00615 typedef FlatVector<T> TELEM; 00616 typedef typename mat_traits<T>::TSCAL TSCAL; 00617 00618 FlatSysVector (int as, int bs, T * adata) 00619 : s(as), blocksize(bs), data(adata) { ; } 00620 00621 FlatSysVector (int as, int bs, LocalHeap & lh) 00622 : s(as), blocksize(bs), data (new (lh) T[as*bs]) 00623 { ; } 00624 00625 /* 00626 FlatSysVector (const FlatSysVector<TSCAL> & sv) 00627 : s(sv.Size()), 00628 blocksize(sv.BlockSize()), 00629 data (sv(0)) 00630 { 00631 ; 00632 } 00633 */ 00634 00635 FlatSysVector & operator= (const FlatSysVector & v) 00636 { 00637 for (int i = 0; i < this->s * this->bs; i++) 00638 data[i] = v.data[i]; 00639 return *this; 00640 } 00641 00642 template<typename TB> 00643 FlatSysVector & operator= (const Expr<TB> & v) 00644 { 00645 return MatExpr<FlatSysVector>::operator= (v); 00646 } 00647 00648 FlatSysVector & operator= (TSCAL s) 00649 { 00650 for (int i = 0; i < this->s*this->bs; i++) 00651 data[i] = s; 00652 return *this; 00653 } 00654 00655 00656 00657 TELEM operator() (int i) 00658 { 00659 return FlatVector<T> (blocksize, &data[i*blocksize]); 00660 } 00661 00662 const TELEM operator() (int i) const 00663 { 00664 return FlatVector<T> (blocksize, &data[i*blocksize]); 00665 } 00666 00667 const TELEM operator() (int i, int j) const 00668 { 00669 return FlatVector<T> (blocksize, &data[i*blocksize]); 00670 } 00671 00672 00673 FlatSysVector<T> Range(int first, int last) 00674 { return FlatSysVector<T> (last-first+1, blocksize, data+(first*blocksize)); } 00675 const FlatSysVector<T> Range(int first, int last) const 00676 { return FlatSysVector<T> (last-first+1, blocksize, data+(first*blocksize)); } 00677 00678 int Size () const { return s; } 00679 int Height () const { return s; } 00680 int Width () const { return 1; } 00681 }; 00682 00683 00684 00689 template <int S, typename T = double> 00690 class Vec : public MatExpr<Vec<S,T> > // , protected BaseVec 00691 { 00693 T data[(S>0)?S:1]; 00694 public: 00696 typedef T TELEM; 00698 typedef typename mat_traits<T>::TSCAL TSCAL; 00700 typedef Vec<S, typename mat_traits<T>::TV_COL> TV_COL; 00702 typedef Vec<1, typename mat_traits<T>::TV_ROW> TV_ROW; 00703 00704 enum { SIZE = S }; 00706 enum { HEIGHT = S }; 00708 enum { WIDTH = 1 }; 00709 00711 Vec () { ; } 00713 Vec (const Vec & v) : MatExpr<Vec> () 00714 { 00715 for (int i = 0; i < S; i++) 00716 data[i] = v.data[i]; 00717 } 00718 00720 Vec (const TSCAL & scal) 00721 { 00722 for (int i = 0; i < S; i++) 00723 data[i] = scal; 00724 } 00725 00727 template<typename TB> 00728 Vec (const Expr<TB> & v) 00729 { 00730 (*this) = v; 00731 } 00732 00734 Vec (const T & s1, const T & s2) 00735 { 00736 data[0] = s1; 00737 data[1] = s2; 00738 } 00739 00741 Vec (const T & s1, const T & s2, const T & s3) 00742 { 00743 data[0] = s1; 00744 data[1] = s2; 00745 data[2] = s3; 00746 } 00747 00749 Vec (const T & s1, const T & s2, const T & s3, const T & s4) 00750 { 00751 data[0] = s1; 00752 data[1] = s2; 00753 data[2] = s3; 00754 data[3] = s4; 00755 } 00756 00757 00758 00760 Vec & operator= (const Vec & v) 00761 { 00762 for (int i = 0; i < S; i++) 00763 data[i] = v.data[i]; 00764 return *this; 00765 } 00766 00767 00769 Vec & operator= (TSCAL scal) 00770 { 00771 for (int i = 0; i < S; i++) 00772 data[i] = scal; 00773 return *this; 00774 } 00775 00777 template<typename TB> 00778 Vec & operator= (const Expr<TB> & v) 00779 { 00780 for (int i = 0; i < S; i++) 00781 data[i] = v.Spec()(i,0); 00782 return *this; 00783 } 00784 00786 TELEM & operator() (int i) 00787 { 00788 #ifdef CHECK_RANGE 00789 CheckVecRange(S,i); 00790 #endif 00791 return data[i]; 00792 } 00793 00795 const TELEM & operator() (int i) const 00796 { 00797 #ifdef CHECK_RANGE 00798 CheckVecRange(S,i); 00799 #endif 00800 return data[i]; 00801 } 00802 00803 00805 TELEM & operator[] (int i) 00806 { 00807 #ifdef CHECK_RANGE 00808 CheckVecRange(S,i); 00809 #endif 00810 return data[i]; 00811 } 00812 00814 const TELEM & operator[] (int i) const 00815 { 00816 #ifdef CHECK_RANGE 00817 CheckVecRange(S,i); 00818 #endif 00819 return data[i]; 00820 } 00821 00823 TELEM & operator() (int i, int j) 00824 { 00825 #ifdef CHECK_RANGE 00826 CheckVecRange(S,i); 00827 #endif 00828 return data[i]; 00829 } 00830 00832 const TELEM & operator() (int i, int j) const 00833 { 00834 #ifdef CHECK_RANGE 00835 CheckVecRange(S,i); 00836 #endif 00837 return data[i]; 00838 } 00839 00841 int Size () const { return S; } 00843 int Height () const { return S; } 00845 int Width () const { return 1; } 00846 00847 const FlatVector<const T> Range(int first, int next) const 00848 { return FlatVector<const T> (next-first, data+first); } 00849 00850 const FlatVector<T> Range(int first, int next) 00851 { return FlatVector<T> (next-first, data+first); } 00852 }; 00853 00854 00855 00856 00857 00859 template <typename S> 00860 inline Vec<3,S> Cross (const Vec<3,S> & a, const Vec<3,S> & b) 00861 { 00862 Vec<3,S> prod; 00863 prod(0) = a(1) * b(2) - a(2) * b(1); 00864 prod(1) = a(2) * b(0) - a(0) * b(2); 00865 prod(2) = a(0) * b(1) - a(1) * b(0); 00866 return prod; 00867 } 00868 00870 template<int S, typename T> 00871 inline ostream & operator<< (ostream & ost, const Vec<S,T> & v) 00872 { 00873 for (int i = 0; i < S; i++) 00874 ost << " " << setw(7) << v(i); 00875 return ost; 00876 } 00877 00878 00879 template<int S, typename TB> 00880 Vec<S> & operator+= (Vec<S> & v, const Expr<TB> & v2) 00881 { 00882 for (int i = 0; i < S; i++) 00883 v(i) += v2.Spec()(i,0); 00884 return v; 00885 } 00886 00887 00888 00889 00890 00891 00892 00893 00897 template <int S, typename T = double> 00898 class FlatVec : public CMCPMatExpr<FlatVec<S,T> > 00899 { 00901 T * data; 00902 public: 00904 typedef T TELEM; 00906 typedef typename mat_traits<T>::TSCAL TSCAL; 00908 typedef Vec<S, typename mat_traits<T>::TV_COL> TV_COL; 00910 typedef Vec<1, typename mat_traits<T>::TV_ROW> TV_ROW; 00911 00912 enum { SIZE = S }; 00914 enum { HEIGHT = S }; 00916 enum { WIDTH = 1 }; 00917 00919 FlatVec (T * adata) : data(adata) { ; } 00920 00922 FlatVec (Vec<S,T> & v2) : data(&v2(0)) { ; } 00923 00924 00926 const FlatVec & operator= (const FlatVec & v) const 00927 { 00928 for (int i = 0; i < S; i++) 00929 data[i] = v.data[i]; 00930 return *this; 00931 } 00932 00934 const FlatVec & operator= (TSCAL scal) const 00935 { 00936 for (int i = 0; i < S; i++) 00937 data[i] = scal; 00938 return *this; 00939 } 00940 00942 template<typename TB> 00943 const FlatVec & operator= (const Expr<TB> & v) const 00944 { 00945 for (int i = 0; i < S; i++) 00946 data[i] = v.Spec()(i,0); 00947 return *this; 00948 } 00949 00950 template<typename TB> 00951 const FlatVec & operator+= (const Expr<TB> & v) const 00952 { 00953 for (int i = 0; i < S; i++) 00954 data[i] += v.Spec()(i,0); 00955 return *this; 00956 } 00957 00959 TELEM & operator() (int i) const 00960 { 00961 #ifdef CHECK_RANGE 00962 CheckVecRange(S,i); 00963 #endif 00964 return data[i]; 00965 } 00966 00968 TELEM & operator[] (int i) const 00969 { 00970 #ifdef CHECK_RANGE 00971 CheckVecRange(S,i); 00972 #endif 00973 return data[i]; 00974 } 00975 00977 TELEM & operator() (int i, int j) const 00978 { 00979 #ifdef CHECK_RANGE 00980 CheckVecRange(S,i); 00981 #endif 00982 return data[i]; 00983 } 00984 00985 const FlatVector<T> Range(int first, int next) const 00986 { return FlatVector<T> (next-first, data+first); } 00987 00989 int Size () const { return S; } 00991 int Height () const { return S; } 00993 int Width () const { return 1; } 00994 }; 00995 00997 template<int S, typename T> 00998 inline ostream & operator<< (ostream & ost, const FlatVec<S,T> & v) 00999 { 01000 for (int i = 0; i < S; i++) 01001 ost << " " << setw(7) << v(i); 01002 return ost; 01003 } 01004 01005 01006 01007 01008 01009 01010 01011 01012 01013 01014 01015 01016 01017 01018 01024 template <typename T = double> 01025 class SliceVector : public CMCPMatExpr<SliceVector<T> > 01026 { 01027 protected: 01029 int s; 01031 int dist; 01033 T * data; 01034 public: 01036 typedef T TELEM; 01038 typedef typename mat_traits<T>::TSCAL TSCAL; 01039 01041 enum { IS_LINEAR = 0 }; 01042 01044 SliceVector (unsigned int as, unsigned int ad, T * adata) 01045 : s(as), dist(ad), data(adata) { ; } 01046 01048 template<typename TB> 01049 const SliceVector & operator= (const Expr<TB> & v) const 01050 { 01051 return CMCPMatExpr<SliceVector>::operator= (v); 01052 } 01053 01055 const SliceVector & operator= (TSCAL scal) const 01056 { 01057 for (int i = 0; i < s; i++) 01058 data[i*dist] = scal; 01059 return *this; 01060 } 01061 01063 const SliceVector & operator= (const SliceVector & v2) const 01064 { 01065 for (int i = 0; i < s; i++) 01066 data[i*dist] = v2(i); 01067 return *this; 01068 } 01069 01070 /* 01071 template<typename TB> 01072 const SliceVector & operator= (const Expr<TB> & v) const 01073 { 01074 if (TB::IS_LINEAR) 01075 for (int i = 0; i < s; i++) 01076 data[i*dist] = v.Spec()(i); 01077 else 01078 for (int i = 0; i < s; i++) 01079 data[i*dist] = v.Spec()(i,0); 01080 return *this; 01081 } 01082 */ 01083 01084 template<typename TB> 01085 const SliceVector & operator+= (const Expr<TB> & v) const 01086 { 01087 if (TB::IS_LINEAR) 01088 for (int i = 0; i < s; i++) 01089 data[i*dist] += v.Spec()(i); 01090 else 01091 for (int i = 0; i < s; i++) 01092 data[i*dist] += v.Spec()(i,0); 01093 return *this; 01094 } 01095 01096 01097 01099 TELEM & operator() (int i) 01100 { 01101 #ifdef CHECK_RANGE 01102 CheckVecRange(s,i); 01103 #endif 01104 return data[i*dist]; 01105 } 01106 01108 TELEM & operator() (int i) const 01109 { 01110 #ifdef CHECK_RANGE 01111 CheckVecRange(s,i); 01112 #endif 01113 return data[i*dist]; 01114 } 01115 01117 TELEM & operator() (int i, int j) const 01118 { 01119 #ifdef CHECK_RANGE 01120 CheckVecRange(s,i); 01121 #endif 01122 return data[i*dist]; 01123 } 01124 01126 TELEM & operator() (int i, int j) 01127 { 01128 #ifdef CHECK_RANGE 01129 CheckVecRange(s,i); 01130 #endif 01131 return data[i*dist]; 01132 } 01133 01135 TELEM & operator[] (int i) 01136 { 01137 #ifdef CHECK_RANGE 01138 CheckVecRange(s,i); 01139 #endif 01140 return data[i*dist]; 01141 } 01142 01144 TELEM & operator[] (int i) const 01145 { 01146 #ifdef CHECK_RANGE 01147 CheckVecRange(s,i); 01148 #endif 01149 return data[i*dist]; 01150 } 01151 01152 TELEM * Addr (int i) const 01153 { 01154 return data+i*dist; 01155 } 01156 01157 01159 int Size () const { return s; } 01160 01162 int Height () const { return s; } 01164 int Width () const { return 1; } 01165 01166 const SliceVector<T> Range (int first, int next) const 01167 { 01168 return SliceVector<T> (next-first, dist, data+first*dist); 01169 } 01170 01171 const SliceVector<T> Range (IntRange range) const 01172 { 01173 return Range (range.First(), range.Next()); 01174 } 01175 01176 01177 const SliceVector<T> Slice (int first, int adist) const 01178 { 01179 return SliceVector<T> (s/adist, dist*adist, data+first*dist); 01180 } 01181 }; 01182 01183 01184 01185 01186 01187 01188 01189 01190 01191 01192 01193 01194 01195 01196 01197 01198 01199 01200 01206 template <int DIST, typename T = double> 01207 class FixSliceVector : public CMCPMatExpr<FixSliceVector<DIST, T> > 01208 { 01209 protected: 01211 int s; 01213 T * data; 01214 public: 01216 typedef T TELEM; 01218 typedef typename mat_traits<T>::TSCAL TSCAL; 01219 01221 enum { IS_LINEAR = 0 }; 01222 01224 FixSliceVector (unsigned int as, T * adata) 01225 : s(as), data(adata) { ; } 01226 01228 template<typename TB> 01229 FixSliceVector & operator= (const Expr<TB> & v) 01230 { 01231 return MatExpr<FixSliceVector>::operator= (v); 01232 } 01233 01235 const FixSliceVector & operator= (TSCAL scal) const 01236 { 01237 for (int i = 0; i < s; i++) 01238 data[i*DIST] = scal; 01239 return *this; 01240 } 01241 01243 const FixSliceVector & operator= (const FixSliceVector & v2) const 01244 { 01245 for (int i = 0; i < s; i++) 01246 data[i*DIST] = v2(i); 01247 return *this; 01248 } 01249 01250 template<typename TB> 01251 const FixSliceVector & operator= (const Expr<TB> & v) const 01252 { 01253 if (TB::IS_LINEAR) 01254 for (int i = 0; i < s; i++) 01255 data[i*DIST] = v.Spec()(i); 01256 else 01257 for (int i = 0; i < s; i++) 01258 data[i*DIST] = v.Spec()(i,0); 01259 return *this; 01260 } 01261 01262 01263 template<typename TB> 01264 const FixSliceVector & operator+= (const Expr<TB> & v) const 01265 { 01266 if (TB::IS_LINEAR) 01267 for (int i = 0; i < s; i++) 01268 data[i*DIST] += v.Spec()(i); 01269 else 01270 for (int i = 0; i < s; i++) 01271 data[i*DIST] += v.Spec()(i,0); 01272 return *this; 01273 } 01274 01275 01276 01278 TELEM & operator() (int i) 01279 { 01280 #ifdef CHECK_RANGE 01281 CheckVecRange(s,i); 01282 #endif 01283 return data[i*DIST]; 01284 } 01285 01287 TELEM & operator() (int i) const 01288 { 01289 #ifdef CHECK_RANGE 01290 CheckVecRange(s,i); 01291 #endif 01292 return data[i*DIST]; 01293 } 01294 01296 TELEM & operator() (int i, int j) const 01297 { 01298 #ifdef CHECK_RANGE 01299 CheckVecRange(s,i); 01300 #endif 01301 return data[i*DIST]; 01302 } 01303 01305 TELEM & operator() (int i, int j) 01306 { 01307 #ifdef CHECK_RANGE 01308 CheckVecRange(s,i); 01309 #endif 01310 return data[i*DIST]; 01311 } 01312 01314 TELEM & operator[] (int i) 01315 { 01316 #ifdef CHECK_RANGE 01317 CheckVecRange(s,i); 01318 #endif 01319 return data[i*DIST]; 01320 } 01321 01323 TELEM & operator[] (int i) const 01324 { 01325 #ifdef CHECK_RANGE 01326 CheckVecRange(s,i); 01327 #endif 01328 return data[i*DIST]; 01329 } 01330 01331 TELEM * Addr (int i) const 01332 { 01333 return data+i*DIST; 01334 } 01335 01336 01338 int Size () const { return s; } 01339 01341 int Height () const { return s; } 01343 int Width () const { return 1; } 01344 01345 const FixSliceVector Range (int first, int next) const 01346 { 01347 return FixSliceVector (next-first, data+first*DIST); 01348 } 01349 01350 const FixSliceVector Range (IntRange range) const 01351 { 01352 return Range (range.First(), range.Next()); 01353 } 01354 01355 01356 const SliceVector<T> Slice (int first, int adist) const 01357 { 01358 return SliceVector<T> (s/adist, DIST*adist, data+first*DIST); 01359 } 01360 }; 01361 01362 01363 01364 01365 01366 01367 01368 01369 01370 01371 template <class TV, class TSCAL> class Scalar2ElemVector 01372 { 01373 public: 01374 const FlatVector<TSCAL> & vec; 01375 Scalar2ElemVector (const FlatVector<TSCAL> & avec) : vec(avec) { ; } 01376 01377 enum { H = mat_traits<TV>::HEIGHT }; 01378 01379 FlatVec<H,TSCAL> operator() (int i) const 01380 { 01381 return FlatVec<H,TSCAL> (&vec(i*H)); 01382 } 01383 01384 }; 01385 01386 01387 template <class TSCAL> class Scalar2ElemVector<TSCAL,TSCAL> 01388 { 01389 public: 01390 const FlatVector<TSCAL> & vec; 01391 Scalar2ElemVector (const FlatVector<TSCAL> & avec) : vec(avec) { ; } 01392 01393 01394 const TSCAL & operator() (int i) const 01395 { 01396 return vec(i); 01397 } 01398 01399 TSCAL & operator() (int i) 01400 { 01401 return vec(i); 01402 } 01403 01404 }; 01405 01406 01407 01408 01409 01410 01411 template <class T> 01412 class mat_traits<FlatVector<T> > 01413 { 01414 public: 01415 typedef T TELEM; 01416 typedef T TSCAL; 01417 }; 01418 01419 01420 } 01421 01422 #ifdef PARALLEL 01423 namespace ngstd 01424 { 01425 template<int S, typename T> 01426 class MPI_Traits<ngbla::Vec<S, T> > 01427 { 01428 public: 01429 static MPI_Datatype MPIType () 01430 { 01431 static MPI_Datatype MPI_T = 0; 01432 if (!MPI_T) 01433 { 01434 MPI_Type_contiguous ( S, MPI_Traits<T>::MPIType(), &MPI_T); 01435 MPI_Type_commit ( &MPI_T ); 01436 } 01437 return MPI_T; 01438 } 01439 }; 01440 } 01441 #endif 01442 01443 #endif