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