vector.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
/*
 * File:   vector.hh
 * Author: ngo
 *
 * Created on April 14th, 2011
 */
 
#ifndef _VECTOR_HH
#define _VECTOR_HH
 
#include <assert.h>
 
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>
 
#include "exceptions.hh"
 
namespace hdnum {
 
  template<typename REAL>
  class Vector : public std::vector<REAL>  // inherit from the STL vector
  {
  public:
    typedef std::size_t size_type;
 
  private:
    static bool bScientific;
    static std::size_t nIndexWidth;
    static std::size_t nValueWidth;
    static std::size_t nValuePrecision;
 
  public:
 
    Vector() : std::vector<REAL>()
    {
    }
 
    Vector( const size_t size,                 // user must specify the size
            const REAL defaultvalue_ = 0    // if not specified, the value 0 will take effect
            )
      : std::vector<REAL>( size, defaultvalue_ )
    {
    }
 
    Vector (const std::initializer_list<REAL> &v)
    {
      for (auto elem : v) this->push_back(elem);
    }
 
    // Methods:
 
    Vector& operator=( const REAL value )
    {
      const size_t s = this->size();
      Vector & self = *this;
      for(size_t i=0; i<s; ++i)
        self[i] = value;
      return *this;
    }
 
    Vector sub (size_type i, size_type m)
    {
      Vector v(m);
      Vector &self = *this;
      size_type k=0;
      for (size_type j=i; j<i+m; j++){
        v[k]=self[j];
        k++;
      }
      return v;
    }
 
 
 
#ifdef DOXYGEN
    Vector& operator=( const Vector& y )
    {
      // It is already implemented in the STL vector class itself!
    }
#endif
 
 
 
    Vector& operator*=( const REAL value )
    {
      Vector &self = *this;
      for (size_t i = 0; i < this->size(); ++i)
        self[i] *= value;
      return *this;
    }
 
 
    Vector& operator/=( const REAL value )
    {
      Vector &self = *this;
      for (size_t i = 0; i < this->size(); ++i)
        self[i] /= value;
      return *this;
    }
 
 
    Vector& operator+=( const Vector & y )
    {
      assert( this->size() == y.size());
      Vector &self = *this;
      for (size_t i = 0; i < this->size(); ++i)
        self[i] += y[i];
      return *this;
    }
 
 
    Vector& operator-=( const Vector & y )
    {
      assert( this->size() == y.size());
      Vector &self = *this;
      for (size_t i = 0; i < this->size(); ++i)
        self[i] -= y[i];
      return *this;
    }
 
 
    Vector & update(const REAL alpha, const Vector & y)
    {
      assert( this->size() == y.size());
      Vector &self = *this;
      for (size_t i = 0; i < this->size(); ++i)
        self[i] += alpha * y[i];
      return *this;
    }
 
 
    REAL operator*(Vector & x) const
    {
      assert( x.size() == this->size() );   // checks if the dimensions of the two vectors are equal
      REAL sum( 0 );
      const Vector & self = *this;
      for( size_t i = 0; i < this->size(); ++i )
        sum += self[i] * x[i];
      return sum;
    }
 
 
 
 
    Vector operator+(Vector & x) const
    {
      assert( x.size() == this->size() );   // checks if the dimensions of the two vectors are equal
      Vector sum( *this );
      sum += x;
      return sum;
    }
 
 
 
    Vector operator-(Vector & x) const
    {
      assert( x.size() == this->size() );   // checks if the dimensions of the two vectors are equal
      Vector sum( *this );
      sum -= x;
      return sum;
    }
 
 
 
    REAL two_norm_2() const
    {
      REAL sum( 0 );
      const Vector & self = *this;
      for (size_t i = 0; i < (size_t) this->size(); ++i)
        sum += self[i] * self[i];
      return sum;
    }
 
    REAL two_norm() const
    {
      return sqrt(two_norm_2());
    }
 
    bool scientific() const
    {
      return bScientific;
    }
 
    void scientific(bool b) const
    {
      bScientific=b;
    }
 
    std::size_t iwidth () const
    {
      return nIndexWidth;
    }
 
    std::size_t width () const
    {
      return nValueWidth;
    }
 
    std::size_t precision () const
    {
      return nValuePrecision;
    }
 
    void iwidth (std::size_t i) const
    {
      nIndexWidth=i;
    }
 
    void width (std::size_t i) const
    {
      nValueWidth=i;
    }
 
    void precision (std::size_t i) const
    {
      nValuePrecision=i;
    }
 
  };
 
 
 
  template<typename REAL>
  bool Vector<REAL>::bScientific = true;
 
  template<typename REAL>
  std::size_t Vector<REAL>::nIndexWidth = 2;
 
  template<typename REAL>
  std::size_t Vector<REAL>::nValueWidth = 15;
 
  template<typename REAL>
  std::size_t Vector<REAL>::nValuePrecision = 7;
 
 
  template <typename REAL>
  inline std::ostream & operator <<(std::ostream & os, const Vector<REAL> & x)
  {
    os << std::endl;
 
    for (size_t r = 0; r < x.size(); ++r)
      {
        if( x.scientific() )
          {
            os << "["
               << std::setw(x.iwidth())
               << r
               << "]"
               << std::scientific
               << std::showpoint
               << std::setw( x.width() )
               << std::setprecision( x.precision() )
              << x[r]
               << std::endl;
          }
        else
          {
            os << "["
               << std::setw(x.iwidth())
               << r
               << "]"
               << std::fixed
               << std::showpoint
               << std::setw( x.width() )
               << std::setprecision( x.precision() )
               << x[r]
               << std::endl;
          }
      }
    return os;
  }
 
 
 
  template<typename REAL>
  inline void gnuplot(
                      const std::string& fname,
                      const Vector<REAL> x
                      )
  {
    std::fstream f(fname.c_str(),std::ios::out);
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      {
        if( x.scientific() )
          {
            f << std::setw(x.width())
              << i
              << std::scientific
              << std::showpoint
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << x[i]
              << std::endl;
          }
        else
          {
            f << std::setw(x.width())
              << i
              << std::fixed
              << std::showpoint
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << x[i]
              << std::endl;
          }
      }
    f.close();
  }
 
  template<typename REAL>
  inline void gnuplot(
                      const std::string& fname,
                      const std::vector<std::string>& t,
                      const Vector<REAL> x
                      )
  {
    std::fstream f(fname.c_str(),std::ios::out);
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      {
        if( x.scientific() )
          {
            f << t[i] << " "
              << std::scientific
              << std::showpoint
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << x[i]
              << std::endl;
          }
        else
          {
            f << t[i] << " "
              << std::fixed
              << std::showpoint
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << x[i]
              << std::endl;
          }
      }
    f.close();
  }
 
  template<typename REAL>
  inline void gnuplot(
                      const std::string& fname,
                      const Vector<REAL> x,
                      const Vector<REAL> y
                      )
  {
    std::fstream f(fname.c_str(),std::ios::out);
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      {
        if( x.scientific() )
          {
            f << std::setw(x.width())
              << i
              << std::scientific
              << std::showpoint
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << x[i]
              << " "
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << y[i]
              << std::endl;
          }
        else
          {
            f << std::setw(x.width())
              << i
              << std::fixed
              << std::showpoint
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << x[i]
              << " "
              << std::setw( x.width() )
              << std::setprecision( x.precision() )
              << y[i]
              << std::endl;
          }
      }
 
    f.close();
  }
 
 
 
  template<typename REAL>
  inline void readVectorFromFile (const std::string& filename, Vector<REAL> &vector)
  {
    std::string buffer;
    std::ifstream fin( filename.c_str() );
    if( fin.is_open() ){
      while( fin ){
        std::string sub;
        fin >> sub;
        //std::cout << " sub = " << sub.c_str() << ": ";
        if( sub.length()>0 ){
          REAL a = atof(sub.c_str());
          //std::cout << std::fixed << std::setw(10) << std::setprecision(5) << a;
          vector.push_back(a);
        }
      }
      fin.close();
    }
    else{
      HDNUM_ERROR("Could not open file!");
    }
  }
 
 
  template<class REAL>
  inline void zero (Vector<REAL>& x)
  {
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      x[i] = REAL(0);
  }
 
  template<class REAL>
  inline REAL norm (Vector<REAL> x)
  {
    REAL sum(0.0);
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      sum += x[i]*x[i];
    return sqrt(sum);
  }
 
  template<class REAL>
  inline void fill (Vector<REAL>& x, const REAL t)
  {
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      x[i] = t;
  }
 
  template<class REAL>
  inline void fill (Vector<REAL>& x, const REAL& t, const REAL& dt)
  {
    REAL myt(t);
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      {
        x[i] = myt;
        myt += dt;
      }
  }
 
 
  template<class REAL>
  inline void unitvector (Vector<REAL> & x, std::size_t j)
  {
    for (typename Vector<REAL>::size_type i=0; i<x.size(); i++)
      if (i==j)
        x[i] = REAL(1);
      else
        x[i] = REAL(0);
  }
 
 
} // end of namespace hdnum
 
#endif  /* _VECTOR_HH */