VVector.hpp

//
// $Id: VVector.hpp,v 1.8 2007/03/06 23:01:17 werner Exp $
//
// $Log: VVector.hpp,v $
// Revision 1.8  2007/03/06 23:01:17  werner
// Integration of the Vector library into the Eagle library, for the sake
// of type registration.
//
// Revision 1.7  2007/02/22 17:08:43  werner
// A major revision of the vector/fixed array interface, unifying three independent
// libraries into one. The vecal (Vector Algebra), fish/vector (Multidimensional
// vector arrays) and the Eagle library (Geometric Algebra) thus share a common
// root now.
//
// Revision 1.6  2006/08/29 04:23:36  werner
// 64 bit long long support
//
// Revision 1.5  2006/08/05 20:41:03  werner
// better doc
//
// Revision 1.4  2006/07/02 22:07:30  werner
// Revision and improvements of the vector/matrix interface, and improvment
// of the examples and self-check.
//
// Revision 1.3  2006/05/26 03:03:12  werner
// Iterator part finished
//
// Revision 1.2  2006/05/24 13:38:03  werner
// better docu, Vector docu does not work...
//
// Revision 1.1  2006/05/20 18:04:12  werner
// Matrix operations embedded and documented
//
// Revision 1.23  2004/11/22 14:54:20  werner
// Added generic interpolator interface and activated HDF5 reading.
//
// Revision 1.22  2004/10/20 20:01:11  werner
// moving towards windows
//
// Revision 1.21  2004/10/05 11:58:54  werner
// Documentation improvements
//
// Revision 1.20  2004/08/25 16:57:33  werner
// Adjustments for SGI
//
// Revision 1.19  2004/08/24 09:09:05  werner
// Towards VC7
//
// Revision 1.18  2004/08/17 16:59:52  werner
// Towards compilation with gcc 3.3.3
//
// Revision 1.17  2004/08/12 16:08:45  werner
// Various improvements for implementing the tangential space
//
// Revision 1.16  2004/07/14 15:45:14  werner
// Outsourcing of quadratic matrices into its own header file.
// Added great gaus decomposition from Marcus Weber that allows
// multiple right-hand sides.
//
// Revision 1.15  2004/07/12 10:30:44  werner
// Minor docu adjustmens.
//
// Revision 1.14  2004/07/06 17:16:05  werner
// Added docu.
//
// Revision 1.13  2004/05/13 16:40:06  werner
// Adjusted iostream interface for IRIX.
//
// Revision 1.12  2004/05/13 10:48:44  werner
// Bug fixes on template argument forwarding.
//
// Revision 1.11  2004/05/12 16:05:23  werner
// Basic library ported to SGI C++.
//
// Revision 1.10  2004/05/11 17:56:10  werner
// *** empty log message ***
//
// Revision 1.9  2004/05/03 13:33:33  werner
// integration improved
//
// Revision 1.8  2004/04/26 13:44:05  werner
// Added const member function on arrays, and minor documentation cleanup.
//
// Revision 1.7  2004/03/30 13:43:10  werner
// Adjusted to Linux environment.
//
// Revision 1.6  2004/03/30 13:16:16  werner
// Spin Interface!
//
// Revision 1.5  2004/03/29 15:23:58  werner
// Adjusted to SGI.
//
// Revision 1.4  2004/03/26 13:35:35  werner
// Improved documentation.
//
// Revision 1.3  2004/03/24 13:22:33  werner
// Scalar constructor added, and inner/outer product of quaternions.
//
// Revision 1.2  2004/03/22 11:38:13  werner
// Improvements and Bugfixes, better conversion functions and operators.
//
// Revision 1.1  2004/03/11 18:23:48  werner
// Fixed the Intel CPU vectorization/alignment issue.
// Vectorization is now done through the explicit VVector<> class, which splits up
// a sequence of types into a vectorizable Vector<> and a remainder.
// As the vectorizable component is based on native SIMD types, they get
// correctly 16-byte aligned, which is all what we want.
//
#ifndef __FIBER_VVector_HPP
#define __FIBER_VVector_HPP "Created 27.02.2001 21:42:27 by werner"

#include <eagle/Vector.hpp>

namespace Eagle
{
//using namespace std;
typedef Eagle::Mult Mult;
typedef Eagle::Div  Div;
typedef Eagle::Add  Add;
typedef Eagle::Sub  Sub;

template <class FirstVector, class SecondVector, class value>
class   VectorPair
{
        FirstVector     first;
        SecondVector    second;
public:
        VectorPair()
        {}

        VectorPair(const VectorPair&V)
        : first( V.first)
        , second(V.second)
        {}



        VectorPair(const VectorPair&A, const value&V, const Mult&)
        : first (A. first, V, Mult())
        , second(A.second, V, Mult())
        {}

        VectorPair(const VectorPair&A, const value&V, const Div&)
        : first (A. first, V, Div())
        , second(A.second, V, Div())
        {}

        VectorPair(const VectorPair&L, const VectorPair&R, const Add&)
        : first (L.first , R.first , Add())
        , second(L.second, R.second, Add())
        {}

        VectorPair(const VectorPair&L, const VectorPair&R, const Sub&)
        : first (L.first , R.first , Sub())
        , second(L.second, R.second, Sub())
        {}

        VectorPair(const VectorPair&L, const Sub&)
        : first (L.first , Sub())
        , second(L.second, Sub())
        {}



        VectorPair&operator+=(const VectorPair&A)
        {
                first  += A.first;
                second += A.second;
                return *this;
        }

        VectorPair&operator-=(const VectorPair&A)
        {
                first  -= A.first;
                second -= A.second;
                return *this;
        }

        VectorPair&operator*=(const value&v)
        {
                first  *= v;
                second *= v;
                return *this;
        }

        VectorPair&operator/=(const value&v)
        {
                first  /= v;
                second /= v;
                return *this;
        }

        friend value InnerProduct(const VectorPair&A, const VectorPair&B)
        {
                return InnerProduct(A.first, B.first) +
                        InnerProduct(A.second, B.second);
        }
};


template <int VN, class Vvalue, int N, class value>
class   MixedVector
{
public:
        typedef VectorPair< Eagle::Vector<Vvalue, VN>, Eagle::Vector<value, N>, value > result;
};

template <class Vvalue, int N, class value>
class   MixedVector<0, Vvalue, N, value>
{
public:
        typedef Eagle::Vector<value,N>  result;
};

template <int VN, class Vvalue, class value>
class   MixedVector<VN, Vvalue, 0, value>
{
public:
        typedef Eagle::Vector<Vvalue, VN> result;
};

//using Eagle::VectorizationTrait;

template <int A, int B>
struct GrandMacOsGCCDivisionCalculator
{
        enum { result = A/B };
};

template <int N, class value>
        class ALIGN16 VVector : public MixedVector< GrandMacOsGCCDivisionCalculator<N,VectorizationTrait<value>::N>::result, 
                                              typename VectorizationTrait<value>::vectorized_t,
                                              N%VectorizationTrait<value>::N,
                                              value>::result
{
public:
        typedef typename MixedVector< GrandMacOsGCCDivisionCalculator<N,VectorizationTrait<value>::N>::result, 
                                     typename VectorizationTrait<value>::vectorized_t,
                                     N%VectorizationTrait<value>::N, 
                                     value >::result MixedVector_t;

        typedef Eagle::Vector<value, N> Vector_t;

        typedef VVector VVector_t;

        typedef value value_type;

        VVector()
        {}

        void set(const value&s)
        {
                vec().set(s);
        }

        VVector(const VVector&V)
        : MixedVector_t(V)
        {}

        template <class T>
        explicit VVector(const Eagle::Vector<T, N>&V)
        {
        value*v = ptr(0);
                for(int i=0; i<N; i++)
                        v[i] = value(V[i]);
        }

        template <class T>
        explicit VVector(const VVector<N, T>&V)
        {
        value*v = ptr(0);
                for(int i=0; i<N; i++)
                        v[i] = value(V[i]);
        }

        explicit VVector(const MixedVector_t&MV)
        : MixedVector_t(MV)
        {}


        const Vector_t&vec() const 
        {
                return *reinterpret_cast<const Vector_t*>(this);        
        }

              Vector_t&vec()
        {
                return *reinterpret_cast<Vector_t*>(this);      
        }

        const VVector&vvec() const
        {
                return *this;
        }

              VVector&vvec()
        {
                return *this;
        }

        value&operator[](int i)
        {       return vec()[i];        }

        const value&operator[](int i) const
        {       return vec()[i];        }


        value*  ptr(int i=0)
        {
                return vec().ptr(i);
        }

        const value*ptr(int i=0) const
        {
                return vec().ptr(i);
        }

        const value*const_ptr(int i=0) const
        {
                return vec().const_ptr(i);
        }


        VVector(const VVector&A, const value&V, const Mult&)
        : MixedVector_t(A,V, Mult())
        {}

        VVector(const VVector&A, const value&V, const Div&)
        : MixedVector_t(A,V, Div())
        {}

        VVector(const VVector&L, const VVector&R, const Add&)
        : MixedVector_t(L,R, Add())
        {}

        VVector(const VVector&L, const VVector&R, const Sub&)
        : MixedVector_t(L,R, Sub())
        {}

        VVector(const VVector&L, const Sub&)
        : MixedVector_t(L, Sub())
        {}



        VVector&operator+=(const VVector&V)
        {
                MixedVector_t::operator+=(V);
                return *this;
        }

        VVector&operator-=(const VVector&V)
        {
                MixedVector_t::operator-=(V);
                return *this;
        }

        VVector&operator*=(const value&v)
        {
                MixedVector_t::operator*=(v);
                return *this;
        }

        VVector&operator/=(const value&v)
        {
                MixedVector_t::operator/=(v);
                return *this;
        }



        friend VVector operator*(const VVector&A, const value&V)
        {
                return VVector(A, V, Mult() );
        }

        friend VVector operator*(const value&V, const VVector&A)
        {
                return VVector(A, V, Mult() );
        }

        friend VVector operator/(const VVector&A, const value&V)
        {
                return VVector(A, V, Div() );
        }

        friend VVector operator+(const VVector&A, const VVector&B)
        {
                return VVector(A, B, Add() );
        }

        friend VVector operator-(const VVector&A, const VVector&B)
        {
                return VVector(A, B, Sub() );
        }



        Eagle::Assignment<Vector_t, 0> operator=(const value&x)
        {
                // we need to perform this operation on the elementary type,
                // since the vectorization cannot be used here
                return Eagle::Assignment<Vector_t, 0>( vec(), x);
        }
};

//}  __attribute__ ((aligned (16)));;

template <int N, class T>
        struct MetaInfo<VVector<N, T> > : MetaInfo<Vector<T, N> >
{
        using MetaInfo<Vector<T, N> >::MULTIPLICITY;
        using MetaInfo<Vector<T, N> >::element_t;
        using MetaInfo<Vector<T, N> >::Chart_t;

static  const T&getComponent(const FixedArray<T,N>&F, int i)
        {
                return F[i];
        }
};


template <int N, class value>
inline VVector<N, value> operator*(const VVector<N, value>&A, double V)
{
        return VVector<N, value>(A, V, Mult() );
}

template <int N, class value>
inline VVector<N, value> operator/(const VVector<N, value>&A, double V)
{
        return VVector<N, value>(A, V, Div() );
}



template <int N, class value>
inline  value norm2(const VVector<N,value>&A)
{
        return InnerProduct(A, A);
}

template <int N, class value>
inline value norm(const VVector<N,value>&A)
{
        return sqrt(norm2(A) );
}


//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)

template <int N, class value>
std::ostream&operator<<(std::ostream&os, const VVector<N, value>&VV)
{
const Eagle::Vector<value,N>&A = VV.vec();
        os << '(';
        for(int i=0; i<N-1; i++)
                os << A[i] << ',';
        os << A[N-1] << ')';
        return os;
}

#endif

} /* namespace Eagle */ 

template <int N, class value>
inline Eagle::VVector<N, value> operator/(const Eagle::VVector<N, value>&A, double V)
{
        return Eagle::VVector<N, value>(A, V, Eagle::Div() );
}


#endif /* __FIBER_VVector_HPP */