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