#ifndef __MYVECTOR_H
#define __MYVECTOR_H
/*!
  \file    vector.hh
  \brief   dynamic numerical vector
  \date    30 April 2010
  \author  Clemens Pechstein

  Institute of Computational Mathematics
  Johannes Kepler University Linz, Austria

  Tutorials -- Numerical methods for elliptic PDEs
  summer semester 2010
*/
#include <cassert>
#include <valarray>

class Vector
{
public:
  /// constructor
  Vector(int size=0) : data_(size) {}
  
  /// copy constructor
  Vector(const Vector& v) : data_(v.data_) {}

  /// assignment
  Vector& operator= (const Vector& v)
  {
    if (&v != this)
      {
        data_.resize (v.size());
        data_ = v.data_;
      }
    return *this;
  }

  /// destructor
  ~Vector () {}
  
  /// size/length/dimension of vector
  int size () const { return (int)(data_.size()); }

  /// change size and set all entries to value
  void resize (int newsize, const double value=0) { data_.resize ((size_t)newsize, value); }

  /// const element access
  const double operator[] (int i) const { assert (0 <= i && i < size()); return data_[i]; }
  
  /// element access
  double& operator[] (int i) { assert (0 <= i && i < size()); return data_[i]; }
  
  /// set all elements to value
  Vector& operator= (const double value) { for (int i=0; i<size(); ++i) { data_[i] = value; } return *this; }

  /// increment
  Vector& operator+= (const Vector& v)
  {
    for (int i=0; i<size(); ++i) data_[i] += v.data_[i];
    return *this;
  }
  
  /// decrement
  Vector& operator-= (const Vector& v)
  {
    for (int i=0; i<size(); ++i) data_[i] -= v.data_[i];
    return *this;
  }
  
  /// scale
  Vector& operator*= (const double factor)
  {
    for (int i=0; i<size(); ++i) data_[i] *= factor;
    return *this;
  }

  /// output
  void print (std::ostream& os) const
  {
    os << "[";
    for (int i=0; i<size(); ++i) os << ' ' << data_[i];
    os << " ]";
  }
  
private:
  std::valarray<double> data_;

}; // class Vector



/// vector addition
inline
Vector operator+ (const Vector& v, const Vector& w)
{
  Vector res(v); // copy
  res += w;
  return res;
}

 
/// vector subtraction
inline
Vector operator- (const Vector& v, const Vector& w)
{
  Vector res(v); // copy
  res -= w;
  return res;
}


/// scale vector
inline
Vector operator* (const double factor, const Vector& v)
{
  Vector res(v);
  res *= factor;
  return res;
}


/// Euclidean scalar product
inline
double inner_product (const Vector& v, const Vector& w)
{
  double res = 0;
  for (int i=0; i<(int)v.size (); ++i) res += v[i]*w[i];
  return res;
}



/// Euclidean norm
inline
double l2norm (const Vector& v)
{
  double res = 0.;
  for (int i=0; i<(int)v.size (); ++i) res += v[i]*v[i];
  return sqrt(res);
}



/// output
inline
std::ostream& operator<< (std::ostream& os, const Vector& v)
{
  v.print (os);
  return os;
}

#endif // __MYVECTOR_H