tnt/tnt_vector.h

/*

*

* Template Numerical Toolkit (TNT)

*

* Mathematical and Computational Sciences Division

* National Institute of Technology,

* Gaithersburg, MD USA

*

*

* This software was developed at the National Institute of Standards and

* Technology (NIST) by employees of the Federal Government in the course

* of their official duties. Pursuant to title 17 Section 105 of the

* United States Code, this software is not subject to copyright protection

* and is in the public domain. NIST assumes no responsibility whatsoever for

* its use by other parties, and makes no guarantees, expressed or implied,

* about its quality, reliability, or any other characteristic.

*

*/







#ifndef TNT_VECTOR_H

#define TNT_VECTOR_H



#include "tnt_subscript.h"

#include <cstdlib>

#include <cassert>

#include <iostream>

#include <sstream>

#include <cmath>



using namespace std;



namespace TNT

{



//namespace Linear_Algebra

// {





template <class T>

class Vector 

{





  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;}

 

  private:

    T* v_;                  

    T* vm1_;        // pointer adjustment for optimzied 1-offset indexing

    Subscript n_;



    // internal helper function to create the array

    // of row pointers



    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);

        v_ = new T[N];

        assert(v_  != NULL);

        vm1_ = v_-1;

        n_ = N;

    }

   

    void copy(const T*  v)

    {

        Subscript N = n_;

        Subscript i;



#ifdef TNT_UNROLL_LOOPS

        Subscript Nmod4 = N & 3;

        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_LOOPS

        Subscript Nmod4 = N & 3;

        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      

    }

    





    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_);     



        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_; }



                operator const T* const() { return v_; } 

                operator T*() { return v_; }



    // destructor



    ~Vector() 

    {

        destroy();

    }



    // constructors



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



    Vector(const Vector<T> &A) : v_(0), vm1_(0), n_(0)

    {

        initialize(A.n_);

        copy(A.v_);

    }



    Vector(Subscript N, const T& value = T()) :  v_(0), vm1_(0), n_(0)

    {

        initialize(N);

        set(value);

    }



    Vector(Subscript N, const T* v) :  v_(0), vm1_(0), n_(0)

    {

        initialize(N);

        copy(v);

    }



    Vector(Subscript N, char *s) :  v_(0), vm1_(0), n_(0)

    {

        initialize(N);

        std::istringstream ins(s);



        Subscript i;



        for (i=0; i<N; i++)

                ins >> v_[i];

    }





    // methods

    // 

    Vector<T>& newsize(Subscript N)

    {

        if (n_ == N) return *this;



        destroy();

        initialize(N);



        return *this;

    }





    // assignments

    //

    Vector<T>& operator=(const Vector<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;

    }

        

    Vector<T>& operator=(const T& scalar)

    { 

        set(scalar);  

        return *this;

    }



    inline Subscript dim() const 

    {

        return  n_; 

    }



    inline Subscript size() const 

    {

        return  n_; 

    }





    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]; 

    }







};





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



template <class T>

std::ostream& operator<<(std::ostream &s, const Vector<T> &A)

{

    Subscript N=A.dim();



    s <<  N << "\n";



    for (Subscript i=0; i<N; i++)

        s   << A[i] << " " << "\n";

    s << "\n";



    return s;

}



template <class T>

std::istream & operator>>(std::istream &s, Vector<T> &A)

{



    Subscript N;



    s >> N;



    if ( !(N == A.size() ))

    {

        A.newsize(N);

    }





    for (Subscript i=0; i<N; i++)

            s >>  A[i];





    return s;

}



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





template <class T>

Vector<T> operator+(const Vector<T> &A, 

    const Vector<T> &B)

{

    Subscript N = A.dim();



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



    Vector<T> tmp(N);

    Subscript i;



    for (i=0; i<N; i++)

            tmp[i] = A[i] + B[i];



    return tmp;

}



template <class T>

Vector<T> operator+=(Vector<T> &A, 

    const Vector<T> &B)

{

    Subscript N = A.dim();



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



    Subscript i;



    for (i=0; i<N; i++)

            A[i] += B[i];



    return A;

}



template <class T>

Vector<T> operator-(const Vector<T> &A, 

    const Vector<T> &B)

{

    Subscript N = A.dim();



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



    Vector<T> tmp(N);

    Subscript i;



    for (i=0; i<N; i++)

            tmp[i] = A[i] - B[i];



    return tmp;

}



template <class T>

Vector<T> operator-=(Vector<T> &A, 

    const Vector<T> &B)

{

    Subscript N = A.dim();



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



    Subscript i;



    for (i=0; i<N; i++)

            A[i] -= B[i];



    return A;

}









template <class T>

Vector<T> elementwise_mult(const Vector<T> &A, const Vector<T> &B)

{

    Subscript N = A.dim();



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



    Vector<T> tmp(N);

    Subscript i;



    for (i=0; i<N; i++)

            tmp[i] = A[i] * B[i];



    return tmp;

}





template <class T>

double norm(const Vector<T> &A)

{

        Subscript N = A.dim();



        double sum = 0.0;

        for (int i=0; i<N; i++)

                sum +=  abs(A[i])*abs(A[i]);

        return sqrt(sum);

}







template <class T>

T dot_prod(const Vector<T> &A, const Vector<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;

}



template <class T>

inline T dot_product(const Vector<T> &A, const Vector<T> &B)

{

        return dot_prod(A, B);

}





template <class T>

inline T operator*(const Vector<T> &A, 

    const Vector<T> &B)

{

    return dot_prod(A,B);

}





template <class T>

Vector<T> operator*(const T &a, const Vector<T> &A)

{

        Subscript N = A.dim();

        Vector<T> r(N);



        for (int i=0; i<N; i++)

                r[i] = A[i] * a;



        return r;

}



template <class T>

inline Vector<T> operator*(const Vector<T> &A, const T& a)

{

        return a * A;

}



//}   /* namspace TNT::Linear_Algebra */

}   /* namespace TNT */





#endif

// TNT_VECTOR_H

---