NGSolve  4.9
linalg/basevector.hpp
00001 #ifndef FILE_BASEVECTOR
00002 #define FILE_BASEVECTOR
00003 
00004 /*********************************************************************/
00005 /* File:   basevector.hpp                                            */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   7. Feb. 2003                                              */
00008 /*********************************************************************/
00009 
00010 
00011 
00012 
00013 namespace ngla
00014 {
00015 
00016   class BaseVector;
00017   typedef auto_ptr<BaseVector> TempVector;
00018   template <class SCAL> class S_BaseVector;
00019 
00020 
00021   class ComplexConjugate;
00022   class ComplexConjugate2;
00023 
00024   template<class IPTYPE>
00025   class SCAL_TRAIT
00026   {
00027   public:
00028     typedef double SCAL;
00029   };
00030 
00031   template<> class SCAL_TRAIT<Complex>
00032   {
00033   public:
00034     typedef Complex SCAL;
00035   };
00036 
00037   template<> class SCAL_TRAIT<ComplexConjugate>
00038   {
00039   public:
00040     typedef Complex SCAL;
00041   };
00042 
00043   template<> class SCAL_TRAIT<ComplexConjugate2>
00044   {
00045   public:
00046     typedef Complex SCAL;
00047   };
00048 
00049 
00050 
00054   template <class T>
00055   class VVecExpr 
00056   {
00057     const T data;
00058   public:
00060     VVecExpr (const T & d) : data(d) { ; }
00061 
00063     template <class TS>
00064     void AssignTo (TS s , BaseVector & v) const { data.AssignTo(s, v); }
00065 
00067     template <class TS>
00068     void AddTo (TS s, BaseVector & v) const { data.AddTo(s, v); }
00069   };
00070 
00071 
00072   enum PARALLEL_STATUS { DISTRIBUTED, CUMULATED, NOT_PARALLEL };
00073 
00074 
00075 
00079   class NGS_DLL_HEADER BaseVector
00080   {
00081   protected:
00083     int size;
00085     int entrysize;
00087     const ParallelDofs * paralleldofs;
00088   public:
00090     BaseVector () throw ()  
00091       : paralleldofs (NULL) 
00092     { ; }
00093 
00095     virtual ~BaseVector ();
00096 
00098     template <typename T> 
00099     BaseVector & operator= (const VVecExpr<T> & v)
00100     {
00101       v.AssignTo (1.0, *this);
00102       return *this;
00103     }
00104 
00106     BaseVector & operator= (const BaseVector & v);
00108     BaseVector & operator= (double s);
00110     BaseVector & operator= (Complex s);
00111 
00113     template <typename T>
00114     BaseVector & operator+= (const VVecExpr<T> & v)
00115     {
00116       v.AddTo (1.0, *this);
00117       return *this;
00118     }
00119 
00121     BaseVector & operator+= (const BaseVector & v)
00122     {
00123       Add (1.0, v);
00124       return *this;
00125     }
00126 
00128     template <typename T>
00129     BaseVector & operator-= (const VVecExpr<T> & v)
00130     {
00131       v.AddTo (-1.0, *this);
00132       return *this;
00133     }
00134 
00136     BaseVector & operator-= (const BaseVector & v)
00137     {
00138       Add (-1.0, v);
00139       return *this;
00140     }
00141 
00143     BaseVector & operator*= (double s)
00144     {
00145       return Scale (s);
00146     }
00147 
00149     BaseVector & operator*= (Complex s)
00150     {
00151       return Scale (s);
00152     }
00153 
00155     BaseVector & operator/= (double s)
00156     {
00157       if (s == 0)
00158         throw Exception ("BaseVector::operator/=: division by zero");
00159       return Scale (1/s);
00160     }
00161 
00163     BaseVector & operator/= (Complex s)
00164     {
00165       if (s == 0.0)
00166         throw Exception ("BaseVector::operator/=: division by zero");
00167       return Scale (1.0/s);
00168     }
00169 
00170     template <class SCAL>
00171     S_BaseVector<SCAL> & Spec()
00172     {
00173       return dynamic_cast<S_BaseVector<SCAL>&> (*this);
00174     }
00175   
00176     template <class SCAL>
00177     const S_BaseVector<SCAL> & Spec() const
00178     {
00179       return dynamic_cast<const S_BaseVector<SCAL>&> (*this);
00180     }
00181 
00182     int Size() const throw () { return size; }
00183     int EntrySize() const throw () { return entrysize; }
00184     virtual void * Memory () const throw () = 0;
00185 
00186     virtual FlatVector<double> FVDouble () const = 0;
00187     virtual FlatVector<Complex> FVComplex () const = 0;
00188 
00189     template <typename T>
00190     FlatVector<T> FV () const;
00191 
00192     template <class TSCAL>
00193     TSCAL InnerProduct (const BaseVector & v2) const 
00194     {
00195       return dynamic_cast<const S_BaseVector<TSCAL>&> (*this) . 
00196         InnerProduct (v2);
00197     }
00198 
00199     virtual double L2Norm () const
00200     {
00201       return ngbla::L2Norm (FVDouble());
00202     }
00203 
00204     virtual BaseVector & Scale (double scal);
00205     virtual BaseVector & Scale (Complex scal);
00206 
00207     virtual BaseVector & SetScalar (double scal);
00208     virtual BaseVector & SetScalar (Complex scal);
00209 
00210     virtual BaseVector & Set (double scal, const BaseVector & v);
00211     virtual BaseVector & Set (Complex scal, const BaseVector & v);
00212 
00213     virtual BaseVector & Add (double scal, const BaseVector & v);
00214     virtual BaseVector & Add (Complex scal, const BaseVector & v);
00215 
00216     virtual ostream & Print (ostream & ost) const;
00217     virtual void Save(ostream & ost) const;
00218     virtual void Load(istream & ist);
00219     virtual void SaveText(ostream & ost) const;
00220     virtual void LoadText(istream & ist);
00221 
00222     virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const;
00223 
00224     // create vector, procs is the set of processors on which the vector exists
00225     // default 0 pointer means all procs
00226     virtual BaseVector * CreateVector ( const Array<int> * procs = 0) const;
00227 
00228 
00229     virtual void SetRandom ();
00230 
00231     virtual BaseVector * Range (int begin, int end) const;
00232     virtual BaseVector * Range (IntRange range) const;
00233 
00234     void GetIndirect (const FlatArray<int> & ind, 
00235                       const FlatVector<double> & v) const;
00236     void GetIndirect (const FlatArray<int> & ind, 
00237                       const FlatVector<Complex> & v) const;
00238     void SetIndirect (const FlatArray<int> & ind, 
00239                       const FlatVector<double> & v);
00240     void SetIndirect (const FlatArray<int> & ind, 
00241                       const FlatVector<Complex> & v);
00242     void AddIndirect (const FlatArray<int> & ind, 
00243                       const FlatVector<double> & v);
00244     void AddIndirect (const FlatArray<int> & ind, 
00245                       const FlatVector<Complex> & v);
00246 
00247 
00248 
00249     template<int S>
00250     void GetIndirect (const Array<int> & ind, 
00251                       FlatVector< Vec<S,double> > & v) const
00252     { 
00253       FlatVector<double> fv = FVDouble();
00254       // int es = EntrySize();
00255       for (int i = 0, ii = 0; i < ind.Size(); i++)
00256         if (ind[i] != -1)
00257           {
00258             int base = S * ind[i];
00259             for (int j = 0; j < S; j++)
00260               v[ii++] = fv[base++];
00261           }
00262         else
00263           {
00264             for (int j = 0; j < S; j++)
00265               v[ii++] = 0;
00266           }
00267     }
00268 
00269     
00270     template<int S>
00271     void GetIndirect (const Array<int> & ind, 
00272                       FlatVector< Vec<S,Complex> > & v) const
00273     { 
00274       FlatVector<Complex> fv = FVComplex();
00275       // int es = EntrySize() / 2;
00276       for (int i = 0, ii = 0; i < ind.Size(); i++)
00277         if (ind[i] != -1)
00278           {
00279             int base = S * ind[i];
00280             for (int j = 0; j < S; j++)
00281               v[ii++] = fv[base++];
00282           }
00283         else
00284           {
00285             for (int j = 0; j < S; j++)
00286               v[ii++] = 0.0;
00287           }
00288     }
00289 
00290 
00291     template<int S>
00292     void AddIndirect (const Array<int> & ind, 
00293                       const FlatVector< Vec<S,double> > & v)
00294     { 
00295       FlatVector<double> fv = FVDouble();
00296       // int es = EntrySize();
00297     
00298       for (int i = 0; i < ind.Size(); i++)
00299         if (ind[i] != -1)
00300           {
00301             int base = S * ind[i];
00302             for (int j = 0; j < S; j++)
00303               fv[base++] += v[i](j);
00304           }
00305     }
00306 
00307     template<int S>
00308     void AddIndirect (const Array<int> & ind, 
00309                       const FlatVector< Vec<S,Complex> > & v)
00310     { 
00311       FlatVector<Complex> fv = FVComplex();
00312       // if(EntrySize() != 2*S)
00313       //       throw Exception("BaseVector::AddIndirect() wrong dimensions");
00314 
00315       for (int i = 0; i < ind.Size(); i++)
00316         if (ind[i] != -1)
00317           {
00318             int base = S * ind[i];
00319             for (int j = 0; j < S; j++)
00320               fv[base++] += v[i](j);
00321           }
00322     }
00323   
00324   
00325     virtual void Cumulate () const { ; }
00326     virtual void Distribute() const { ; }
00327     virtual PARALLEL_STATUS GetParallelStatus () const { return NOT_PARALLEL; }
00328     virtual void SetParallelStatus (PARALLEL_STATUS stat) const { ; }
00329   };
00330 
00331 
00332 
00333   template <>
00334   inline FlatVector<double> BaseVector::FV<double> () const
00335   {
00336     return FVDouble();
00337   }
00338 
00339   template <>
00340   inline FlatVector<Complex> BaseVector::FV<Complex> () const
00341   {
00342     return FVComplex();
00343   }
00344 
00345   template <typename T>
00346   inline FlatVector<T> BaseVector::FV () const
00347   {
00348     typedef typename mat_traits<T>::TSCAL TSCAL;
00349     return FlatVector<T> (Size(), FV<TSCAL>().Addr(0));
00350   }
00351 
00352 
00353 
00354 
00355 
00356 
00360   template <class SCAL>
00361   class NGS_DLL_HEADER S_BaseVector : virtual public BaseVector
00362   {
00363   public:
00364     S_BaseVector () throw () { ; }
00365     virtual ~S_BaseVector() { ; }
00366 
00367     S_BaseVector & operator= (double s);
00368 
00369     virtual SCAL InnerProduct (const BaseVector & v2) const
00370     {
00371       return ngbla::InnerProduct (FVScal(), 
00372                                   dynamic_cast<const S_BaseVector&>(v2).FVScal());
00373     }
00374 
00375     virtual FlatVector<double> FVDouble () const;
00376     virtual FlatVector<Complex> FVComplex () const;
00377 
00378     virtual FlatVector<SCAL> FVScal () const throw() 
00379     {
00380       return FlatVector<SCAL> (size * entrysize, Memory());
00381     }
00382   };
00383 
00384 
00385 
00386 
00387   template <>
00388   class NGS_DLL_HEADER S_BaseVector<Complex> : virtual public BaseVector
00389   {
00390   public:
00391     S_BaseVector () throw() { ; }
00392     ~S_BaseVector () { ; }
00393 
00394     virtual Complex InnerProduct (const BaseVector & v2) const
00395     {
00396       return ngbla::InnerProduct (FVScal(), 
00397                                   dynamic_cast<const S_BaseVector&>(v2).FVScal());
00398     }
00399 
00400     virtual FlatVector<double> FVDouble () const throw();
00401     virtual FlatVector<Complex> FVComplex () const throw();
00402     virtual FlatVector<Complex> FVScal () const throw() 
00403     {
00404       return FlatVector<Complex> (size * entrysize/2, Memory());
00405     }
00406   };
00407 
00408 
00409 
00410 
00411 
00412 
00413 
00414 
00415 
00416   /* ********************* Expression templates ***************** */
00417 
00418 
00419 
00420   template <> class VVecExpr<BaseVector>
00421   {
00422     const BaseVector & v;
00423   public:
00424     VVecExpr (const BaseVector & av) : v(av) { ; }
00425 
00426     template <class TS>
00427     void AssignTo (TS s, BaseVector & v2) const { v2.Set (s, v); }
00428     template <class TS>
00429     void AddTo (TS s, BaseVector & v2) const { v2.Add (s,  v); }
00430   };
00431 
00432 
00433 
00434   /* ***************************** VSumExpr ************************** */
00435 
00437   template <class TA, class TB>
00438   class VSumExpr
00439   {
00440     const TA a;
00441     const TB b;
00442   
00443   public:
00444     VSumExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; }
00445 
00446     template <class TS>
00447     void AssignTo (TS s, BaseVector & v) const
00448     { 
00449       a.AssignTo (s, v);
00450       b.AddTo (s, v);
00451     }
00452     template <class TS>
00453     void AddTo (TS s, BaseVector & v) const
00454     { 
00455       a.AddTo (s, v);
00456       b.AddTo (s, v);
00457     }
00458   };
00459 
00460 
00461 
00462   inline VVecExpr<VSumExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > >
00463   operator+ (const BaseVector & a, const BaseVector & b)
00464   {
00465     typedef VSumExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > TRES;
00466     return TRES (a, b);
00467   }
00468 
00469   template <class TA>
00470   inline VVecExpr<VSumExpr<VVecExpr<TA>, VVecExpr<BaseVector> > >
00471   operator+ (const VVecExpr<TA> & a, const BaseVector & b)
00472   {
00473     typedef VSumExpr<VVecExpr<TA>, VVecExpr<BaseVector> > TRES;
00474     return TRES (a, b);
00475   }
00476 
00477   template <class TB>
00478   inline VVecExpr<VSumExpr<VVecExpr<BaseVector>, VVecExpr<TB> > >
00479   operator+ (const BaseVector & a, const VVecExpr<TB> & b)
00480   {
00481     typedef VSumExpr<VVecExpr<BaseVector>, VVecExpr<TB> > TRES;
00482     return TRES (a, b);
00483   }
00484 
00485   template <class TA, class TB>
00486   inline VVecExpr<VSumExpr<VVecExpr<TA>, VVecExpr<TB> > >
00487   operator+ (const VVecExpr<TA> & a, const VVecExpr<TB> & b)
00488   {
00489     typedef VSumExpr<VVecExpr<TA>, VVecExpr<TB> > TRES;
00490     return TRES (a, b);
00491   }
00492 
00493 
00494 
00495 
00496 
00497 
00498 
00499 
00500   /* ***************************** VSubExpr ************************** */
00501 
00503   template <class TA, class TB>
00504   class VSubExpr
00505   {
00506     const TA a;
00507     const TB b;
00508   
00509   public:
00510     VSubExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; }
00511 
00512 
00513     template <class TS>
00514     void AssignTo (TS s, BaseVector & v) const
00515     { 
00516       a.AssignTo (s, v);
00517       b.AddTo (-s, v);
00518     }
00519     template <class TS>
00520     void AddTo (TS s, BaseVector & v) const
00521     { 
00522       a.AddTo (s, v);
00523       b.AddTo (-s, v);
00524     }
00525   };
00526 
00527 
00528 
00529   inline VVecExpr<VSubExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > >
00530   operator- (const BaseVector & a, const BaseVector & b)
00531   {
00532     typedef VSubExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > TRES;
00533     return TRES (a, b);
00534   }
00535 
00536   template <class TA>
00537   inline VVecExpr<VSubExpr<VVecExpr<TA>, VVecExpr<BaseVector> > >
00538   operator- (const VVecExpr<TA> & a, const BaseVector & b)
00539   {
00540     typedef VSubExpr<VVecExpr<TA>, VVecExpr<BaseVector> > TRES;
00541     return TRES (a, b);
00542   }
00543 
00544   template <class TB>
00545   inline VVecExpr<VSubExpr<VVecExpr<BaseVector>, VVecExpr<TB> > >
00546   operator- (const BaseVector & a, const VVecExpr<TB> & b)
00547   {
00548     typedef VSubExpr<VVecExpr<BaseVector>, VVecExpr<TB> > TRES;
00549     return TRES (a, b);
00550   }
00551 
00552   template <class TA, class TB>
00553   inline VVecExpr<VSubExpr<VVecExpr<TA>, VVecExpr<TB> > >
00554   operator- (const VVecExpr<TA> & a, const VVecExpr<TB> & b)
00555   {
00556     typedef VSubExpr<VVecExpr<TA>, VVecExpr<TB> > TRES;
00557     return TRES (a, b);
00558   }
00559 
00560 
00561 
00562 
00563 
00564   /* ************************* Scal * Vec ******************** */
00565 
00566 
00568   template <class TA, class TSCAL>
00569   class VScaleExpr
00570   {
00571     const TA a;
00572     const TSCAL scal;
00573   
00574   public:
00575     VScaleExpr (const TA & aa, const TSCAL & as) : a(aa), scal(as) { ; }
00576 
00577 
00578     template <class TS>
00579     void AssignTo (TS s, BaseVector & v) const
00580     { 
00581       a.AssignTo (scal * s, v);
00582     }
00583     template <class TS>
00584     void AddTo (TS s, BaseVector & v) const
00585     { 
00586       a.AddTo (scal * s, v);
00587     }
00588   };
00589 
00590 
00591   inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, double> >
00592   operator* (const BaseVector & a, const double & b)
00593   {
00594     typedef VScaleExpr<VVecExpr<BaseVector>, double> TRES;
00595     return TRES (a, b);
00596   }
00597 
00598   template <class TA>
00599   inline VVecExpr<VScaleExpr<VVecExpr<TA>, double> >
00600   operator* (const VVecExpr<TA> & a, const double & b)
00601   {
00602     typedef VScaleExpr<VVecExpr<TA>, double> TRES;
00603     return TRES (a, b);
00604   }
00605 
00606 
00607 
00608   inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, Complex> >
00609   operator* (const BaseVector & a, const Complex & b)
00610   {
00611     typedef VScaleExpr<VVecExpr<BaseVector>, Complex> TRES;
00612     return TRES (a, b);
00613   }
00614 
00615   template <class TA>
00616   inline VVecExpr<VScaleExpr<VVecExpr<TA>, Complex> >
00617   operator* (const VVecExpr<TA> & a, const Complex & b)
00618   {
00619     typedef VScaleExpr<VVecExpr<TA>, Complex> TRES;
00620     return TRES (a, b);
00621   }
00622 
00623 
00624 
00625 
00626 
00627   inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, double> >
00628   operator* (const double & b, const BaseVector & a)
00629   {
00630     typedef VScaleExpr<VVecExpr<BaseVector>, double> TRES;
00631     return TRES (a, b);
00632   }
00633 
00634   template <class TA>
00635   inline VVecExpr<VScaleExpr<VVecExpr<TA>, double> >
00636   operator* (const double & b, const VVecExpr<TA> & a)
00637   {
00638     typedef VScaleExpr<VVecExpr<TA>, double> TRES;
00639     return TRES (a, b);
00640   }
00641 
00642 
00643 
00644   inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, Complex> >
00645   operator* (const Complex & b, const BaseVector & a)
00646   {
00647     typedef VScaleExpr<VVecExpr<BaseVector>, Complex> TRES;
00648     return TRES (a, b);
00649   }
00650 
00651   template <class TA>
00652   inline VVecExpr<VScaleExpr<VVecExpr<TA>, Complex> >
00653   operator* (const Complex & b, const VVecExpr<TA> & a)
00654   {
00655     typedef VScaleExpr<VVecExpr<TA>, Complex> TRES;
00656     return TRES (a, b);
00657   }
00658 
00659 
00660   template <class TA>
00661   inline VVecExpr<VScaleExpr<VVecExpr<TA>,double > >
00662   operator- (const VVecExpr<TA> & a)
00663   {
00664     typedef VScaleExpr<VVecExpr<TA>, double> TRES;
00665     return TRES (a, -1);
00666   }
00667 
00668 
00669 
00670 
00671   /* *********************** operator<< ********************** */
00672 
00674   inline ostream & operator<< (ostream & ost, const BaseVector & v)
00675   {
00676     return v.Print(ost);
00677   }
00678 
00680   inline double InnerProduct (const BaseVector & v1, const BaseVector & v2)
00681   {
00682     return dynamic_cast<const S_BaseVector<double>&>(v1).InnerProduct(v2); 
00683   }
00684 
00685 
00687   template <class IPTYPE>
00688   inline typename SCAL_TRAIT<IPTYPE>::SCAL S_InnerProduct (const BaseVector & v1, const BaseVector & v2)
00689   {
00690     return dynamic_cast<const S_BaseVector<typename SCAL_TRAIT<IPTYPE>::SCAL>&>(v1).InnerProduct(v2); 
00691   }
00692 
00693   template <>
00694   inline Complex 
00695   S_InnerProduct<ComplexConjugate> (const BaseVector & v1, const BaseVector & v2)
00696   {
00697     return InnerProduct( v1.FVComplex(), Conj(v2.FVComplex()) );
00698   }
00699 
00700   template <>
00701   inline Complex S_InnerProduct<ComplexConjugate2> (const BaseVector & v1, const BaseVector & v2)
00702   {
00703     return InnerProduct( v2.FVComplex(), Conj(v1.FVComplex()) );
00704   }
00705 
00707   inline double L2Norm (const BaseVector & v)
00708   {
00709     return v.L2Norm();
00710   }
00711 
00712 }
00713 
00714 
00715 
00716 #endif