source: trunk/lib/gslapi/matrix.h @ 295

// $Id: matrix.h 295 2005-04-29 09:15:58Z peter $

#ifndef _theplu_gslapi_matrix_ 
#define _theplu_gslapi_matrix_ 

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

#include <gsl/gsl_matrix.h>
#include <iostream>


namespace theplu {
namespace gslapi {

  class vector;

  ///
  /// This is the C++ tools interface to GSL matrix. 'double' is the
  /// only type supported, maybe we should add a 'complex' type as
  /// well in the future.
  ///
  class matrix
  {
  public:

    ///
    /// The default constructor.
    ///
    matrix(void);

    ///
    /// Constructor. Allocates memory space for \a r times \a c
    /// elements, and sets all elements to \a init_value.
    ///
    matrix(const size_t& r, const size_t& c, const double init_value=0);
 
    ///
    /// The copy constructor.
    ///
    matrix(const matrix&);

    ///
    /// The istream constructor.
    ///
    /// Matrix rows are sepearated with the new line character, and
    /// column elements in a row must be separated with white space
    /// characters except the new line (\\n) character. Rows, and
    /// columns, are read consequently and only rectangular matrices
    /// are supported. Empty lines are ignored. End of input is at end
    /// of file marker.
    ///
    explicit matrix(std::istream &);

    ///
    /// The destructor.
    ///
    ~matrix(void);

    ///
    /// @return The number of columns in the matrix.
    ///
    inline size_t columns(void) const { return m_->size2; }

    ///
    /// Swap column \a i with \a j.
    ///
    inline void column_swap(const size_t& i,const size_t& j)
      { gsl_matrix_swap_columns(m_,i,j); }

    ///
    /// @return true if each element deviates less or equal than \a
    /// d. If any matrix contain a nan false is always returned.
    ///
    bool equal(const matrix&, const double d=0) const;

    ///
    /// Create a new copy of the internal GSL matrix.
    ///
    /// Necessary memory for the new GSL matrix is allocated and the
    /// caller is responsible for freeing the allocated memory.
    ///
    /// @return A pointer to a copy of the internal GSL matrix.
    ///
    gsl_matrix* gsl_matrix_copy(void) const;

    ///
    /// @return A const pointer to the internal GSL matrix.
    ///
    inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; }

    ///
    /// @return A pointer to the internal GSL matrix.
    ///
    inline gsl_matrix* gsl_matrix_pointer(void) { return m_; };

    ///
    /// This function multiplies all elements between two matrices and
    /// the result is stored back into the calling object, \f$a_{ij} =
    /// a_{ij} * b_{ij} \; \forall i,j\f$.
    ///
    inline void mul_elements(const matrix& b)
      { gsl_matrix_mul_elements(m_,b.m_); }

    ///
    /// Swap row \a i with \a j.
    ///
    inline void row_swap( const size_t& i, const size_t& j)
      { gsl_matrix_swap_rows(m_,i,j); }

    ///
    /// @return The number of rows in the matrix.
    ///
    inline size_t rows(void) const { return m_->size1; }

    ///
    /// Set all elements to \a value.
    ///
    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }

    ///
    /// Set column \a i to \a vec. 
    ///
    inline void set_column(const size_t i, const vector vec) 
    {gsl_matrix_set_col(m_, i, vec.gsl_vector_pointer());}

    ///
    /// Set row \a i to \a vec. 
    ///
    inline void set_row(const size_t i, const vector vec) 
    {gsl_matrix_set_row(m_, i, vec.gsl_vector_pointer());}

    ///
    /// Set column \a i to \a vec
    ///
    //void matrix set_column(const size_t i, const vector vec) const;

    ///
    /// Transpose the matrix.
    ///
    /// @return Transpose of the matrix.
    ///
    matrix transpose(void) const;
    
    ///
    /// @return Reference to the element position (\a row, \a column).
    ///
    inline double& operator()(size_t row,size_t column)
      { return (*gsl_matrix_ptr(m_,row,column)); }

    ///
    /// @return Const reference to the element position (\a row, \a
    /// column).
    ///
    inline const double& operator()(size_t row,size_t column) const
      { return (*gsl_matrix_const_ptr(m_,row,column)); }

    ///
    /// Jari, To be implemented. Should return a reference to a vector
    /// object. Underlying GSL supports only GSL_vector views, thus
    /// some redesign of the vector class is needed.
    ///
    /// inline vector& operator[](size_t i);
    /// So for now, stupid copying is done, and actually a matrix
    /// row is returned using this function.
    vector operator[](size_t row) const;

    ///
    /// Jari, to be properly implemented. Should return a reference to
    /// a vector object (matrix row). Underlying GSL supports only GSL_vector
    /// views, thus some redesign of the vector class is needed.
    ///
    /// inline const vector& operator[](size_t i) const;
    ///
    /// So for now, stupid copying is done, and actually a matrix
    /// column is returned using this function.
    vector TEMP_col_return(size_t column) const;

    ///
    /// Matrix multiplication.
    ///
    /// @return The resulting matrix.
    ///
    matrix operator*(const matrix&) const;

    ///
    /// Matrix-vector multiplication.
    ///
    /// @return The resulting vector.
    ///
    vector operator*(const vector&) const;

    ///
    /// Matrix addition.
    ///
    /// @return The resulting matrix.
    ///
    matrix operator+(const matrix&) const;

    ///
    /// Matrix subtraction.
    ///
    /// @return The resulting matrix.
    ///
    matrix operator-(const matrix&) const;

    ///
    /// Comparison operator.
    ///
    /// @return True if all elements are equal otherwise False.
    ///
    /// @see equal
    ///
    inline bool operator==( const matrix& other ) const { return equal(other);}

    ///
    /// Comparison operator.
    ///
    /// @return False if all elements are equal otherwise True.
    ///
    /// @see equal
    ///
    inline bool operator!=( const matrix& other ) const {return !equal(other);}

    ///
    /// The assignment operator. There is no requirements on
    /// dimensions.
    ///
    matrix& operator=(const matrix& other);

    ///
    /// Subtract and assign operator.
    ///
    inline matrix& operator-=(const matrix& m)
      { gsl_matrix_sub(m_,m.m_); return *this; }

    ///
    /// Multiply and assign operator.
    ///
    inline matrix& operator*=(const double d)
      { gsl_matrix_scale(m_,d); return *this; }


  private:

    gsl_matrix* m_;

  };

  ///
  /// The output operator for the matrix class.
  ///
  std::ostream& operator<< (std::ostream& s, const matrix&);



}} // of namespace gslapi and namespace theplu

#endif