VVector.hpp

00001 //
00002 // $Id: VVector.hpp,v 1.8 2007/03/06 23:01:17 werner Exp $
00003 //
00004 // $Log: VVector.hpp,v $
00005 // Revision 1.8  2007/03/06 23:01:17  werner
00006 // Integration of the Vector library into the Eagle library, for the sake
00007 // of type registration.
00008 //
00009 // Revision 1.7  2007/02/22 17:08:43  werner
00010 // A major revision of the vector/fixed array interface, unifying three independent
00011 // libraries into one. The vecal (Vector Algebra), fish/vector (Multidimensional
00012 // vector arrays) and the Eagle library (Geometric Algebra) thus share a common
00013 // root now.
00014 //
00015 // Revision 1.6  2006/08/29 04:23:36  werner
00016 // 64 bit long long support
00017 //
00018 // Revision 1.5  2006/08/05 20:41:03  werner
00019 // better doc
00020 //
00021 // Revision 1.4  2006/07/02 22:07:30  werner
00022 // Revision and improvements of the vector/matrix interface, and improvment
00023 // of the examples and self-check.
00024 //
00025 // Revision 1.3  2006/05/26 03:03:12  werner
00026 // Iterator part finished
00027 //
00028 // Revision 1.2  2006/05/24 13:38:03  werner
00029 // better docu, Vector docu does not work...
00030 //
00031 // Revision 1.1  2006/05/20 18:04:12  werner
00032 // Matrix operations embedded and documented
00033 //
00034 // Revision 1.23  2004/11/22 14:54:20  werner
00035 // Added generic interpolator interface and activated HDF5 reading.
00036 //
00037 // Revision 1.22  2004/10/20 20:01:11  werner
00038 // moving towards windows
00039 //
00040 // Revision 1.21  2004/10/05 11:58:54  werner
00041 // Documentation improvements
00042 //
00043 // Revision 1.20  2004/08/25 16:57:33  werner
00044 // Adjustments for SGI
00045 //
00046 // Revision 1.19  2004/08/24 09:09:05  werner
00047 // Towards VC7
00048 //
00049 // Revision 1.18  2004/08/17 16:59:52  werner
00050 // Towards compilation with gcc 3.3.3
00051 //
00052 // Revision 1.17  2004/08/12 16:08:45  werner
00053 // Various improvements for implementing the tangential space
00054 //
00055 // Revision 1.16  2004/07/14 15:45:14  werner
00056 // Outsourcing of quadratic matrices into its own header file.
00057 // Added great gaus decomposition from Marcus Weber that allows
00058 // multiple right-hand sides.
00059 //
00060 // Revision 1.15  2004/07/12 10:30:44  werner
00061 // Minor docu adjustmens.
00062 //
00063 // Revision 1.14  2004/07/06 17:16:05  werner
00064 // Added docu.
00065 //
00066 // Revision 1.13  2004/05/13 16:40:06  werner
00067 // Adjusted iostream interface for IRIX.
00068 //
00069 // Revision 1.12  2004/05/13 10:48:44  werner
00070 // Bug fixes on template argument forwarding.
00071 //
00072 // Revision 1.11  2004/05/12 16:05:23  werner
00073 // Basic library ported to SGI C++.
00074 //
00075 // Revision 1.10  2004/05/11 17:56:10  werner
00076 // *** empty log message ***
00077 //
00078 // Revision 1.9  2004/05/03 13:33:33  werner
00079 // integration improved
00080 //
00081 // Revision 1.8  2004/04/26 13:44:05  werner
00082 // Added const member function on arrays, and minor documentation cleanup.
00083 //
00084 // Revision 1.7  2004/03/30 13:43:10  werner
00085 // Adjusted to Linux environment.
00086 //
00087 // Revision 1.6  2004/03/30 13:16:16  werner
00088 // Spin Interface!
00089 //
00090 // Revision 1.5  2004/03/29 15:23:58  werner
00091 // Adjusted to SGI.
00092 //
00093 // Revision 1.4  2004/03/26 13:35:35  werner
00094 // Improved documentation.
00095 //
00096 // Revision 1.3  2004/03/24 13:22:33  werner
00097 // Scalar constructor added, and inner/outer product of quaternions.
00098 //
00099 // Revision 1.2  2004/03/22 11:38:13  werner
00100 // Improvements and Bugfixes, better conversion functions and operators.
00101 //
00102 // Revision 1.1  2004/03/11 18:23:48  werner
00103 // Fixed the Intel CPU vectorization/alignment issue.
00104 // Vectorization is now done through the explicit VVector<> class, which splits up
00105 // a sequence of types into a vectorizable Vector<> and a remainder.
00106 // As the vectorizable component is based on native SIMD types, they get
00107 // correctly 16-byte aligned, which is all what we want.
00108 //
00110 #ifndef __FIBER_VVector_HPP
00111 #define __FIBER_VVector_HPP "Created 27.02.2001 21:42:27 by werner"
00112 
00113 #include <eagle/Vector.hpp>
00114 
00115 namespace Eagle
00116 {
00117 //using namespace std;
00118 typedef Eagle::Mult Mult;
00119 typedef Eagle::Div  Div;
00120 typedef Eagle::Add  Add;
00121 typedef Eagle::Sub  Sub;
00122 
00129 template <class FirstVector, class SecondVector, class value>
00130 class   VectorPair
00131 {
00132         FirstVector     first;
00133         SecondVector    second;
00134 public:
00136         VectorPair()
00137         {}
00138 
00140         VectorPair(const VectorPair&V)
00141         : first( V.first)
00142         , second(V.second)
00143         {}
00144 
00145 
00147 
00149         VectorPair(const VectorPair&A, const value&V, const Mult&)
00150         : first (A. first, V, Mult())
00151         , second(A.second, V, Mult())
00152         {}
00153 
00155         VectorPair(const VectorPair&A, const value&V, const Div&)
00156         : first (A. first, V, Div())
00157         , second(A.second, V, Div())
00158         {}
00159 
00161         VectorPair(const VectorPair&L, const VectorPair&R, const Add&)
00162         : first (L.first , R.first , Add())
00163         , second(L.second, R.second, Add())
00164         {}
00165 
00167         VectorPair(const VectorPair&L, const VectorPair&R, const Sub&)
00168         : first (L.first , R.first , Sub())
00169         , second(L.second, R.second, Sub())
00170         {}
00171 
00173         VectorPair(const VectorPair&L, const Sub&)
00174         : first (L.first , Sub())
00175         , second(L.second, Sub())
00176         {}
00178 
00179 
00181 
00183         VectorPair&operator+=(const VectorPair&A)
00184         {
00185                 first  += A.first;
00186                 second += A.second;
00187                 return *this;
00188         }
00189 
00191         VectorPair&operator-=(const VectorPair&A)
00192         {
00193                 first  -= A.first;
00194                 second -= A.second;
00195                 return *this;
00196         }
00197 
00199         VectorPair&operator*=(const value&v)
00200         {
00201                 first  *= v;
00202                 second *= v;
00203                 return *this;
00204         }
00205 
00207         VectorPair&operator/=(const value&v)
00208         {
00209                 first  /= v;
00210                 second /= v;
00211                 return *this;
00212         }
00214 
00216         friend value InnerProduct(const VectorPair&A, const VectorPair&B)
00217         {
00218                 return InnerProduct(A.first, B.first) +
00219                         InnerProduct(A.second, B.second);
00220         }
00221 };
00222 
00223 
00224 template <int VN, class Vvalue, int N, class value>
00225 class   MixedVector
00226 {
00227 public:
00228         typedef VectorPair< Eagle::Vector<Vvalue, VN>, Eagle::Vector<value, N>, value > result;
00229 };
00230 
00231 template <class Vvalue, int N, class value>
00232 class   MixedVector<0, Vvalue, N, value>
00233 {
00234 public:
00235         typedef Eagle::Vector<value,N>  result;
00236 };
00237 
00238 template <int VN, class Vvalue, class value>
00239 class   MixedVector<VN, Vvalue, 0, value>
00240 {
00241 public:
00242         typedef Eagle::Vector<Vvalue, VN> result;
00243 };
00244 
00245 //using Eagle::VectorizationTrait;
00246 
00255 template <int A, int B>
00256 struct GrandMacOsGCCDivisionCalculator
00257 {
00258         enum { result = A/B };
00259 };
00260 
00283 template <int N, class value>
00284         class ALIGN16 VVector : public MixedVector< GrandMacOsGCCDivisionCalculator<N,VectorizationTrait<value>::N>::result, 
00285                                               typename VectorizationTrait<value>::vectorized_t,
00286                                               N%VectorizationTrait<value>::N,
00287                                               value>::result
00288 {
00289 public:
00295         typedef typename MixedVector< GrandMacOsGCCDivisionCalculator<N,VectorizationTrait<value>::N>::result, 
00296                                      typename VectorizationTrait<value>::vectorized_t,
00297                                      N%VectorizationTrait<value>::N, 
00298                                      value >::result MixedVector_t;
00299 
00301         typedef Eagle::Vector<value, N> Vector_t;
00302 
00304         typedef VVector VVector_t;
00305 
00307         typedef value value_type;
00308 
00310         VVector()
00311         {}
00312 
00314         void set(const value&s)
00315         {
00316                 vec().set(s);
00317         }
00318 
00320         VVector(const VVector&V)
00321         : MixedVector_t(V)
00322         {}
00323 
00325         template <class T>
00326         explicit VVector(const Eagle::Vector<T, N>&V)
00327         {
00328         value*v = ptr(0);
00329                 for(int i=0; i<N; i++)
00330                         v[i] = value(V[i]);
00331         }
00332 
00336         template <class T>
00337         explicit VVector(const VVector<N, T>&V)
00338         {
00339         value*v = ptr(0);
00340                 for(int i=0; i<N; i++)
00341                         v[i] = value(V[i]);
00342         }
00343 
00345         explicit VVector(const MixedVector_t&MV)
00346         : MixedVector_t(MV)
00347         {}
00348 
00350 
00352         const Vector_t&vec() const 
00353         {
00354                 return *reinterpret_cast<const Vector_t*>(this);        
00355         }
00356 
00358               Vector_t&vec()
00359         {
00360                 return *reinterpret_cast<Vector_t*>(this);      
00361         }
00362 
00364         const VVector&vvec() const
00365         {
00366                 return *this;
00367         }
00368 
00370               VVector&vvec()
00371         {
00372                 return *this;
00373         }
00374 
00376         value&operator[](int i)
00377         {       return vec()[i];        }
00378 
00380         const value&operator[](int i) const
00381         {       return vec()[i];        }
00382 
00383 
00385         value*  ptr(int i=0)
00386         {
00387                 return vec().ptr(i);
00388         }
00389 
00391         const value*ptr(int i=0) const
00392         {
00393                 return vec().ptr(i);
00394         }
00395 
00397         const value*const_ptr(int i=0) const
00398         {
00399                 return vec().const_ptr(i);
00400         }
00402 
00404 
00405         VVector(const VVector&A, const value&V, const Mult&)
00406         : MixedVector_t(A,V, Mult())
00407         {}
00408 
00409         VVector(const VVector&A, const value&V, const Div&)
00410         : MixedVector_t(A,V, Div())
00411         {}
00412 
00413         VVector(const VVector&L, const VVector&R, const Add&)
00414         : MixedVector_t(L,R, Add())
00415         {}
00416 
00417         VVector(const VVector&L, const VVector&R, const Sub&)
00418         : MixedVector_t(L,R, Sub())
00419         {}
00420 
00421         VVector(const VVector&L, const Sub&)
00422         : MixedVector_t(L, Sub())
00423         {}
00425 
00426 
00428 
00430         VVector&operator+=(const VVector&V)
00431         {
00432                 MixedVector_t::operator+=(V);
00433                 return *this;
00434         }
00435 
00437         VVector&operator-=(const VVector&V)
00438         {
00439                 MixedVector_t::operator-=(V);
00440                 return *this;
00441         }
00442 
00444         VVector&operator*=(const value&v)
00445         {
00446                 MixedVector_t::operator*=(v);
00447                 return *this;
00448         }
00449 
00451         VVector&operator/=(const value&v)
00452         {
00453                 MixedVector_t::operator/=(v);
00454                 return *this;
00455         }
00457 
00458 
00460 
00462         friend VVector operator*(const VVector&A, const value&V)
00463         {
00464                 return VVector(A, V, Mult() );
00465         }
00466 
00468         friend VVector operator*(const value&V, const VVector&A)
00469         {
00470                 return VVector(A, V, Mult() );
00471         }
00472 
00474         friend VVector operator/(const VVector&A, const value&V)
00475         {
00476                 return VVector(A, V, Div() );
00477         }
00478 
00480         friend VVector operator+(const VVector&A, const VVector&B)
00481         {
00482                 return VVector(A, B, Add() );
00483         }
00484 
00486         friend VVector operator-(const VVector&A, const VVector&B)
00487         {
00488                 return VVector(A, B, Sub() );
00489         }
00490 
00492 
00493 
00495         Eagle::Assignment<Vector_t, 0> operator=(const value&x)
00496         {
00497                 // we need to perform this operation on the elementary type,
00498                 // since the vectorization cannot be used here
00499                 return Eagle::Assignment<Vector_t, 0>( vec(), x);
00500         }
00501 };
00502 
00503 //}  __attribute__ ((aligned (16)));;
00504 
00505 template <int N, class T>
00506         struct MetaInfo<VVector<N, T> > : MetaInfo<Vector<T, N> >
00507 {
00508         using MetaInfo<Vector<T, N> >::MULTIPLICITY;
00509         using MetaInfo<Vector<T, N> >::element_t;
00510         using MetaInfo<Vector<T, N> >::Chart_t;
00511 
00512 static  const T&getComponent(const FixedArray<T,N>&F, int i)
00513         {
00514                 return F[i];
00515         }
00516 };
00517 
00518 
00520 template <int N, class value>
00521 inline VVector<N, value> operator*(const VVector<N, value>&A, double V)
00522 {
00523         return VVector<N, value>(A, V, Mult() );
00524 }
00525 
00527 template <int N, class value>
00528 inline VVector<N, value> operator/(const VVector<N, value>&A, double V)
00529 {
00530         return VVector<N, value>(A, V, Div() );
00531 }
00532 
00533 
00534 
00536 template <int N, class value>
00537 inline  value norm2(const VVector<N,value>&A)
00538 {
00539         return InnerProduct(A, A);
00540 }
00541 
00543 template <int N, class value>
00544 inline value norm(const VVector<N,value>&A)
00545 {
00546         return sqrt(norm2(A) );
00547 }
00548 
00549 
00550 //               MSVC                           GCC             G++                   SGI C++
00551 #if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)
00552 
00553 template <int N, class value>
00554 std::ostream&operator<<(std::ostream&os, const VVector<N, value>&VV)
00555 {
00556 const Eagle::Vector<value,N>&A = VV.vec();
00557         os << '(';
00558         for(int i=0; i<N-1; i++)
00559                 os << A[i] << ',';
00560         os << A[N-1] << ')';
00561         return os;
00562 }
00563 
00564 #endif
00565 
00566 } /* namespace Eagle */ 
00567 
00569 template <int N, class value>
00570 inline Eagle::VVector<N, value> operator/(const Eagle::VVector<N, value>&A, double V)
00571 {
00572         return Eagle::VVector<N, value>(A, V, Eagle::Div() );
00573 }
00574 
00575 
00576 #endif /* __FIBER_VVector_HPP */
00577 
00578