vec.h

Go to the documentation of this file.

#ifndef VEC_H
#define VEC_H

#include "haconst.h"
#include <stdlib.h>
#include <assert.h>

template <class T>
class NumVector
{

  public:

    typedef Subscript   size_type;
    typedef         T   value_type;
    typedef         T   element_type;
    typedef         T*  pointer;
    typedef         T*  iterator;
    typedef         T&  reference;
    typedef const   T*  const_iterator;
    typedef const   T&  const_reference;

    Subscript lbound() const { return 1;}
        
  protected:
    T* v_;                  
    T* vm1_;        // pointer adjustment for optimzied 1-offset indexing
    Subscript n_;

    // internal helper function to create the array
    // of row pointers

  public:

    void initialize(Subscript N)
    {
        // adjust pointers so that they are 1-offset:
        // v_[] is the internal contiguous array, it is still 0-offset
        //
        assert(v_ == NULL);
        if(N==0) 
                        return;
//        v_ = new T[N];
        v_=(pointer)malloc(N*sizeof(T));
        assert(v_  != NULL);
        vm1_ = v_-1;
        n_ = N;
    }
   
    void copy(const T*  v)
    {
        Subscript N = n_;
        Subscript i;

#ifdef TNT_UNROLL
        Subscript Nmod4 = N & 4;
        Subscript N4 = N - Nmod4;

        for (i=0; i<N4; i+=4)
        {
            v_[i] = v[i];
            v_[i+1] = v[i+1];
            v_[i+2] = v[i+2];
            v_[i+3] = v[i+3];
        }

        for (i=N4; i< N; i++)
            v_[i] = v[i];
#else

        for (i=0; i< N; i++)
            v_[i] = v[i];
#endif      
    }

    void set(const T& val)
    {
        Subscript N = n_;
        Subscript i;

#ifdef TNT_UNROLL
        Subscript Nmod4 = N & 4;
        Subscript N4 = N - Nmod4;

        for (i=0; i<N4; i+=4)
        {
            v_[i] = val;
            v_[i+1] = val;
            v_[i+2] = val;
            v_[i+3] = val; 
        }

        for (i=N4; i< N; i++)
            v_[i] = val;
#else

        for (i=0; i< N; i++)
            v_[i] = val;
        
#endif      
    }

    
        bool operator == (const NumVector<T>& rhs ) const       
        {
                if(n_ != rhs.n_) return false;
                int i;
                for(i = 0; i < n_; i++)
                {
                        if(v_[i] != rhs.v_[i]) return false;
                }
                return true;
        }

        bool operator < (const NumVector<T>& rhs) const
        {
                if(n_ > rhs.n_) return false;
                if(n_ < rhs.n_) return true;
                int i;
                for(i = 0; i < n_; i++)
                {
                        if(v_[i] > rhs.v_[i]) return false;
                        if(v_[i] < rhs.v_[i]) return true;
                }
                return false;
        }

    void destroy()
    {     
        /* do nothing, if no memory has been previously allocated */
        if (v_ == NULL) return ;

        /* if we are here, then matrix was previously allocated */
//        delete [] (v_);     
        free(v_);
            n_ = 0;
        v_ = NULL;
        vm1_ = NULL;
      }


  public:

    // access

    iterator begin() { return v_;}
    iterator end()   { return v_ + n_; }
    const iterator begin() const { return v_;}
    const iterator end() const  { return v_ + n_; }

    // destructor

    ~NumVector() 
    {
        destroy();
    }

    // constructors

    NumVector() : v_(0), vm1_(0), n_(0)  {};

    NumVector(const NumVector<T> &A) : v_(0), vm1_(0), n_(0)
    {
        initialize(A.n_);
        copy(A.v_);
    }

    NumVector(Subscript N, const T& value = T(0)) :  v_(0), vm1_(0), n_(0)
    {
        initialize(N);
        set(value);
    }

    NumVector(Subscript N, const T* v) :  v_(0), vm1_(0), n_(0)
    {
        initialize(N);
        copy(v);
    }

    NumVector(Subscript N, char *s) :  v_(0), vm1_(0), n_(0)
    {
        initialize(N);
        istrstream ins(s);

        Subscript i;

        for (i=0; i<N; i++)
                ins >> v_[i];
    }


    // methods
    // 
    void newsize(Subscript N)
    {
        if (n_ == N) return;

        destroy();
        initialize(N);

        return;
    }


    // assignments
    //
    NumVector<T>& operator=(const NumVector<T> &A)
    {
        if (v_ == A.v_)
            return *this;

        if (n_ == A.n_)         // no need to re-alloc
            copy(A.v_);

        else
        {
            destroy();
            initialize(A.n_);
            copy(A.v_);
        }

        return *this;
    }
        
    NumVector<T>& operator=(const T& scalar)
    { 
        set(scalar);  
        return *this;
    }

    Subscript dim() const 
    {
        return  n_; 
    }

    Subscript size() const 
    {
        return  n_; 
    }

        inline value_type GetVal(Subscript i) const
        {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= n_) ;
#endif
        return vm1_[i]; 
    }

        inline void SetVal(Subscript i, const value_type new_val)
        {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= n_) ;
#endif
                vm1_[i] = new_val;
        }



    inline reference operator()(Subscript i)
    { 
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= n_) ;
#endif
        return vm1_[i]; 
    }

    inline const_reference operator() (Subscript i) const
    {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= n_) ;
#endif
        return vm1_[i]; 
    }

    inline reference operator[](Subscript i)
    { 
#ifdef TNT_BOUNDS_CHECK
        assert(0<=i);
        assert(i < n_) ;
#endif
        return v_[i]; 
    }

    inline const_reference operator[](Subscript i) const
    {
#ifdef TNT_BOUNDS_CHECK
        assert(0<=i);
        assert(i < n_) ;
#endif
        return v_[i]; 
    }

#if defined(_MSC_VER) || defined(__DECCXX) || defined(__xlC__) 
    friend istream& operator>>(istream &s, NumVector<T> &A);
#else
    friend istream& operator>><>(istream &s, NumVector<T> &A);
#endif

};


/* ***************************  I/O  ********************************/

template <class T>
ostream& operator<<(ostream &s, const NumVector<T> &A)
{
    Subscript N=A.dim();

    s <<  N << endl;

    for (Subscript i=0; i<N; i++)
        s   << A[i] << " " << endl;
    s << endl;

    return s;
}

template <class T>
istream& operator>>(istream &s, NumVector<T> &A)
{

    Subscript N;

    s >> N;

    if ( !(N == A.n_) )
    {
        A.destroy();
        A.initialize(N);
    }


    for (Subscript i=0; i<N; i++)
            s >>  A[i];


    return s;
}

//*******************[ basic matrix algorithms ]***************************


template <class T>
NumVector<T> operator+(const NumVector<T> &A, 
    const NumVector<T> &B)
{
    Subscript N = A.dim();

    assert(N==B.dim());

    NumVector<T> tmp(N);
    Subscript i;

    for (i=0; i<N; i++)
            tmp[i] = A[i] + B[i];

    return tmp;
}

template <class T>
NumVector<T> operator-(const NumVector<T> &A, 
    const NumVector<T> &B)
{
    Subscript N = A.dim();

    assert(N==B.dim());

    NumVector<T> tmp(N);
    Subscript i;

    for (i=0; i<N; i++)
            tmp[i] = A[i] - B[i];

    return tmp;
}

template <class T>
NumVector<T> operator*(const NumVector<T> &A, 
    const NumVector<T> &B)
{
    Subscript N = A.dim();

    assert(N==B.dim());

    NumVector<T> tmp(N);
    Subscript i;

    for (i=0; i<N; i++)
            tmp[i] = A[i] * B[i];

    return tmp;
}

template <class T>
NumVector<T> operator/(const NumVector<T> &A, 
    const NumVector<T> &B)
{
    Subscript N = A.dim();

    assert(N==B.dim());

    NumVector<T> tmp(N);
    Subscript i;

    for (i=0; i<N; i++)
            tmp[i] = A[i] / B[i];

    return tmp;
}

/*
template <class T>
int operator < (const NumVector<T> &A, 
    const NumVector<T> &B)
{
   Subscript N = A.dim();

    if(N > B.dim()) return 0;
        if(N < B.dim()) return 1;

    Subscript i;

    for (i=0; i<N; i++)
        {
       if( (A[i] - B[i]) > 1.0E-10)  return 0;
           if( (A[i] - B[i]) < -1.0E-10) return 1;
        }

    return 1;
}

template <class T>
int operator == (const NumVector<T> &A, 
    const NumVector<T> &B)
{
    Subscript N = A.dim();

    if(N != B.dim()) return 0;

    Subscript i;

    for (i=0; i<N; i++)
        {
       if( fabs(A[i] - B[i]) > 1.0E-10) return 0;
        }

    return 1;
}
*/

template <class T>
int operator == (const NumVector<T> &A, 
    const T & t)
{
    Subscript N = A.dim();

    Subscript i;

    for (i=0; i<N; i++)
        {
       if( A[i] != t) return 0;
        }

    return 1;
}



template <class T>
T dot_prod(const NumVector<T> &A, const NumVector<T> &B)
{
    Subscript N = A.dim();
    assert(N == B.dim());

    Subscript i;
    T sum = 0;

    for (i=0; i<N; i++)
        sum += A[i] * B[i];

    return sum;
}

#endif /* if(!Defined(VEC_H))  */
 

---