source: trunk/lib/gslapi/matrix.cc @ 357

// $Id: matrix.cc 357 2005-08-03 14:53:06Z peter $

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

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

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

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>


namespace theplu {
namespace gslapi {



  matrix::matrix(void)
    : m_(NULL)
  {
  }



  matrix::matrix(const size_t& r, const size_t& c, double init_value)
  {
    m_ = gsl_matrix_alloc(r,c);
    set_all(init_value);
  }



  matrix::matrix(const matrix& other)
  {
    m_ = other.gsl_matrix_copy();
  }



  // Constructor that gets data from istream
  matrix::matrix(std::istream& is)
  {
    // read the data file and store in stl vectors (dynamically
    // expandable)
    std::vector<std::vector<double> > data_matrix;
    u_int nof_columns=0;
    u_int nof_rows = 0;
    std::vector<double> v;
    for (nof_rows = 0; utility::read_to_double(is, v); nof_rows++) {
    
      // Ignoring empty lines
      if(!v.size()){
        nof_rows--;
        continue;
      }
      
      if(nof_columns==0)
        nof_columns=v.size();
      else if(v.size()!=nof_columns) {
        std::cerr << "matrix::matrix(std::istream&) data file error: "
                  << "line" << nof_rows+1 << " has " << v.size() 
                  << " columns; expected " << nof_columns
                  << " columns"
                  << std::endl; 
        exit(1);
      }
      data_matrix.push_back(v);
    }

    // manipulate the state of the stream to be good
    is.clear(std::ios::goodbit);
    // convert the data to a gsl matrix
    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    for(u_int i=0;i<nof_rows;i++) 
      for(u_int j=0;j<nof_columns;j++) 
        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
  }



  matrix::~matrix(void)
  {
    if (m_) {
      gsl_matrix_free(m_);
      m_=NULL;
    }
  }



  gsl_matrix* matrix::gsl_matrix_copy(void) const
  {
    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
    gsl_matrix_memcpy(m,m_);
    return m;
  }

  bool matrix::equal(const matrix& other, const double d) const
  {
    if (columns()!=other.columns() || rows()!=other.rows())
      return false;
    for (size_t i=0; i<rows(); i++) 
      for (size_t j=0; j<columns(); j++) 
        // two last cond. are to check for nans
        if (fabs( (*this)(i,j)-other(i,j) ) > d || 
            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
          return false;

    return true;
  }

  vector matrix::TEMP_col_return(size_t column) const
  {
    vector vec(rows());
    for (size_t i=0; i<rows(); ++i)
      vec(i)=operator()(i,column);
    return vec;
  }



  vector matrix::operator[](size_t row) const
  {
    vector vec(columns());
    for (size_t i=0; i<columns(); ++i)
      vec(i)=operator()(row,i);
    return vec;
  }



  matrix matrix::operator*( const matrix &other ) const
  {
    matrix res(rows(), other.columns());
    gsl_linalg_matmult(m_,other.m_,res.m_);
    return res;
  }



  vector matrix::operator*( const vector &other ) const
  {
    gsl_vector* gslvec=gsl_vector_alloc(rows());
    gsl_blas_dgemv(CblasNoTrans, 1.0, m_, other.gsl_vector_pointer(), 0.0,
                   gslvec );
    vector res(gslvec);
    return res;
  }



  matrix matrix::operator+( const matrix &other ) const
  {
    matrix res( *this );
    gsl_matrix_add( res.m_, other.m_ );
    return res;
  }



  matrix matrix::operator-( const matrix &other ) const
  {
    matrix res( *this );
    gsl_matrix_sub( res.m_, other.m_ );
    return res;
  }


  
  matrix& matrix::operator=( const matrix& other )
  {
    if ( this != &other ) {
      if ( m_ )
        gsl_matrix_free( m_ );
      m_ = other.gsl_matrix_copy();
    }
    return *this;
  }



  matrix matrix::transpose(void) const
  {
    matrix res(columns(),rows());
    gsl_matrix_transpose_memcpy(res.m_,m_);
    return res;
  }



  std::ostream& operator<<(std::ostream& s, const matrix& m)
  {
    using namespace std;
    s.setf( ios::dec );
    s.precision(12);
    for(size_t i=0, j=0; i<m.rows(); i++) 
      for (j=0; j<m.columns(); j++) {
        s << m(i,j);
        if (j<m.columns()-1)
          s << "\t";
        else if (i<m.rows()-1)
          s << "\n";
      }
    return s;
  }


}} // of namespace gslapi and namespace thep