source: trunk/lib/gslapi/vector.cc @ 314

// $Id: vector.cc 310 2005-05-03 19:17:18Z peter $

#include <c++_tools/gslapi/vector.h>
#include <c++_tools/gslapi/matrix.h>

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <utility>





namespace theplu {
namespace gslapi {



  vector::vector(void)
    : v_(NULL), view_(NULL)
  {
  }



  vector::vector(const size_t n,const double init_value)
    : view_(NULL)
  {
    v_ = gsl_vector_alloc(n);
    set_all(init_value);
  }



  vector::vector(const vector& other)
    : view_(NULL)
  {
    v_ = other.TEMP_gsl_vector_copy();
  }



  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
  {
    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
                                                                 stride,n));
    v_ = &(view_->vector);
  }



  vector::vector(gsl_vector* vec)
    : v_(vec), view_(NULL)
  {
  }



  vector::vector(std::istream& is)
    : view_(NULL)
  {
    matrix m(is);
    
    // convert the data to a gsl vector and check that data file is a
    // column vector or a row vector
    if(m.columns()==1) {
      v_ = gsl_vector_alloc ( m.rows() );
      for(u_int i=0;i<m.rows();i++) 
        gsl_vector_set( v_, i, m(i,0) );
    }
    else if(m.rows()==1){
      v_ = gsl_vector_alloc ( m.columns() );
      for(u_int i=0;i<m.columns();i++) 
        gsl_vector_set( v_, i, m(0,i) );
    }
    else {
      // Jari. change to throw exception, remove unnecessary includes
      std::cerr << "vector::vector(std::istream&) data file error: "
                << "file has " << m.rows() << " rows and " << m.columns() 
                << " columns; expected a row vector or a column vector"
                << std::endl;
      exit(1);
    }
  }



  vector::~vector(void)
  {
    if (view_)
      delete view_;
    else if (v_)
      gsl_vector_free(v_);
    v_=NULL;
  }



  gsl_vector* vector::TEMP_gsl_vector_copy(void) const
  {
    gsl_vector* vec = gsl_vector_alloc(size());
    gsl_vector_memcpy(vec,v_);
    return vec;
  }



  std::pair<double,double> vector::minmax(void) const
  {
    double min, max;
    gsl_vector_minmax(v_,&min,&max);
    return std::pair<double,double>(min,max);
  }



  std::pair<size_t,size_t> vector::minmax_index(void) const
  {
    size_t min_index, max_index;
    gsl_vector_minmax_index(v_,&min_index,&max_index);
    return std::pair<size_t,size_t>(min_index,max_index);
  }



  double vector::sum(void) const
  {
    double sum = 0;
    for (size_t i=0; i<size(); i++)
      sum += gsl_vector_get( v_, i );
    return( sum );
  }  



  vector vector::operator-(void) const
  {
    vector res( *this );
    gsl_vector_scale (res.v_,-1.);
    return res;
  }



  double vector::operator*( const vector &other ) const
  {
    // Jari, check for gsl_support
    double res = 0.0;;
    for ( size_t i = 0; i < size(); ++i ) 
      res += other[i] * (*this)[i];
    return res;
  }



  vector vector::operator+(const vector &other) const
  {
    vector res(*this);
    gsl_vector_add(res.v_,other.v_);
    return res;
  }



  vector vector::operator-( const vector &other ) const
  {
    vector res( *this );
    gsl_vector_sub(res.v_,other.v_);
    return res;
  }



  vector& vector::operator=( const vector& other )
  {
    if( this != &other ) {
      if ( v_ )
        gsl_vector_free( v_ );
      v_ = other.TEMP_gsl_vector_copy();
    }
    return *this;
  } 



  std::ostream& operator<<(std::ostream& s, const vector& a)
  {
    s.setf(std::ios::fixed);

    for (size_t j = 0; j < a.size(); ++j) {
      s << a[j];
      if ( (j+1)<a.size() )
        s << " ";
    }

    return s;
  }


}} // of namespace gslapi and namespace thep