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

// $Id: matrix.h 536 2006-03-03 21:19:29Z peter $

#ifndef _theplu_gslapi_matrix_
#define _theplu_gslapi_matrix_

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

#include <cassert>
#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.
  ///
  /// \par[File streams] Reading and writing vectors to file streams
  /// are of course supported. These are implemented without using GSL
  /// functionality, and thus binary read and write to streams are not
  /// supported.
  ///
  /// \par[Matrix views] GSL matrix views are supported and these are
  /// disguised as ordinary gslapi::matrix'es. A support function is
  /// added, gslapi::matrix::isview(), that can be used to check if a
  /// matrix object is a view. Note that view matricess do not own the
  /// undelying data, and a view is not valid if the matrix owning the
  /// data is deallocated.
  ///
  /// @note All GSL matrix related functions are not implement but
  /// most functionality defined for GSL matrices can be achieved with
  /// this interface class. Most notable GSL functionality not
  /// supported are no binary file support and views on arrays,
  /// gslapi::vectors, gsl_vectors, diagonals, subdiagonals, and
  /// superdiagonals. If there is an interest from the user community,
  /// the omitted functionality could be included.
  ///
  /// @todo Maybe it would be smart to create temporary objects need
  /// for BLAS in constructors. As it is now temporary objects are
  /// called before BLAS functionality iss used, cf. const matrix&
  /// matrix::operator*=(const matrix& other)
  ///
  class matrix
  {
  public:

    ///
    /// The default constructor.
    ///
    /// This contructor does not initialize underlying (essential)
    /// structures.
    ///
    inline matrix(void) : m_(NULL), view_(NULL) {}

    ///
    /// Constructor. Allocates memory space for \a r times \a c
    /// elements, and sets all elements to \a init_value.
    ///
    inline matrix(const size_t& r, const size_t& c, double init_value=0)
      : view_(NULL) { m_ = gsl_matrix_alloc(r,c); set_all(init_value); }

    ///
    /// The copy constructor.
    ///
    /// @note If the object to be copied is a matrix view, the values
    /// of the view will be copied, i.e. the view is not copied.
    ///
    inline matrix(const matrix& o) : view_(NULL) {m_=o.create_gsl_matrix_copy();}

    ///
    /// The matrix view constructor.
    ///
    /// Create a view of matrix \a m, with starting row \a offset_row,
    /// starting column \a offset_column, row size \a n_row, and
    /// column size \a n_column.
    ///
    /// A matrix view can be used as any matrix with the difference
    /// that changes made to the view will also change the object that
    /// is viewed. Also, using the copy constructor will create a new
    /// matrix object that is a copy of whatever is viewed. If a copy
    /// of the view is needed then you should use this constructor to
    /// obtain a copy.
    ///
    /// @note If the object viewed by the view goes out of scope or is
    /// deleted, the view becomes invalid and the result of further
    /// use is undefined.
    ///
    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
           size_t n_column);

    ///
    /// The istream constructor.
    ///
    /// Matrix rows are sepearated with the new line character, and
    /// column elements in a row must be separated either with white space
    /// characters except the new line (\\n) character (default), or 
    /// by the delimiter \a sep.
    ///. Rows, and
    /// columns, are read consecutively and only rectangular matrices
    /// are supported. Empty lines are ignored. End of input is at end
    /// of file marker.
    ///
    //    explicit matrix(std::istream &) throw (utility::IO_error,std::exception);
    explicit matrix(std::istream &, char sep='\0') 
      throw (utility::IO_error,std::exception);


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

    ///
    /// Elementwise addition of the elements of matrix \a b to the
    /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \;
    /// \forall i,j\f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }

    ///
    /// Add the scalar value \a d to the elements of the calling
    /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result
    /// is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int
    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }

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

    ///
    /// Elementwise division of the elemnts of the calling matrix by
    /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \;
    /// \forall i,j\f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int
    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }

    ///
    /// Check whether matrices are equal within a user defined
    /// precision, set by \a precision.
    ///
    /// @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 precision=0) const;

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

    ///
    /// @return A pointer to the internal GSL matrix.
    ///
    // Jari, is this needed?
    inline gsl_matrix* gsl_matrix_p(void) { return m_; };

    ///
    /// @return True if all elements in the matrix is zero, false
    /// othwerwise;
    ///
    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }

    ///
    /// Check if the matrix object is a view (sub-matrix) to another
    /// matrix.
    ///
    /// @return True if the object is a view, false othwerwise.
    ///
    inline bool isview(void) const { return view_; }

    ///
    /// @return The maximum value of the matrix.
    ///
    inline double max(void) const { return gsl_matrix_max(m_); }

    ///
    /// @return The index to the element with the maximum value of the
    /// matrix. The lowest index has precedence (searching in
    /// row-major order).
    ///
    inline std::pair<size_t,size_t> max_index(void) const;

    ///
    /// @return The minimum value of the matrix.
    ///
    inline double min(void) const { return gsl_matrix_min(m_); }

    ///
    /// @return The index to the element with the minimum value of the
    /// matrix. The lowest index has precedence (searching in
    /// row-major order).
    ///
    inline std::pair<size_t,size_t> min_index(void) const;

    ///
    /// @return The indecies to the element with the minimum and
    /// maximum values of the matrix, respectively. The lowest index
    /// has precedence (searching in row-major order).
    ///
    inline void minmax_index(std::pair<size_t,size_t>& min,
                             std::pair<size_t,size_t>& max) const
      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
                                &max.first,&max.second); }

    ///
    /// @return The minimum and maximum values of the matrix, as the
    /// \a first and \a second member of the returned \a pair,
    /// respectively.
    ///
    std::pair<double,double> minmax(void) const;

    ///
    /// @return matrix containing 1 and 0 only. 1 means element in
    /// matrix is valid and a zero means element is a NaN.
    ///
    matrix nan(void) const;

    ///
    /// Multiply the elements of matrix \a b with the elements of the
    /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall
    /// i,j\f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int
    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }

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

    ///
    /// Multiply the elements of the calling matrix with a scalar \a
    /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is
    /// stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }

    ///
    /// Set element values to values in \a mat. This function is
    /// needed for assignment of viewed elements.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    /// @see const matrix& operator=(const matrix&)
    ///
    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }

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

    ///
    /// Set \a column values to values in \a vec.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    inline int set_column(const size_t column, const vector& vec)
      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }

    ///
    /// Set \a row values to values in \a vec.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    inline int set_row(const size_t row, const vector& vec)
      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }

    ///
    /// Subtract the elements of matrix \a b from the elements of the
    /// calling matrix ,\f$a_{ij} = a_{ij} - b_{ij} \; \forall
    /// i,j\f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }

    ///
    /// Exchange the elements of the this and \a other by copying. The
    /// two matrices must have the same size.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }

    ///
    /// Swap columns \a i and \a j.
    ///
    inline int swap_columns(const size_t i,const size_t j)
      { return gsl_matrix_swap_columns(m_,i,j); }

    ///
    /// Swap row \a i and column \a j.
    ///
    inline int swap_rowcol(const size_t i,const size_t j)
      { return gsl_matrix_swap_rowcol(m_,i,j); }

    ///
    /// Swap rows \a i and \a j.
    ///
    inline int swap_rows(const size_t i, const size_t j)
      { return gsl_matrix_swap_rows(m_,i,j); }

    ///
    /// Transpose the matrix.
    ///
    void transpose(void);

    ///
    /// @return Reference to the element position (\a row, \a column).
    ///
    inline double& operator()(size_t row,size_t column)
    { 
      assert(row<rows());
      assert(column<columns());
      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
    {       
      assert(row<rows());
      assert(column<columns());
      return (*gsl_matrix_const_ptr(m_,row,column)); 
    }

    ///
    /// Matrix-vector multiplication.
    ///
    /// @return The resulting vector.
    ///
    // Jari, where should this go?
    const vector operator*(const vector&) 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, i.e. the matrix is remapped in memory if
    /// necessary. This implies that in general views cannot be
    /// assigned using this operator. Views will be mutated into
    /// normal matrices. The only exception to this behaviour on views
    /// is when self-assignemnt is done, since self-assignment is
    /// ignored.
    ///
    /// @return A const reference to the resulting matrix.
    ///
    /// @see int set(const matrix&)
    ///
    const matrix& operator=(const matrix& other);

    ///
    /// Add and assign operator.
    ///
    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }

    ///
    /// Add and assign operator.
    ///
    inline const matrix&
    operator+=(const double d) { add_constant(d); return *this; }

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

    ///
    /// Multiply and assigment operator.
    ///
    /// @return Const reference to the resulting matrix.
    ///
    const matrix& operator*=(const matrix&);

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


  private:

    ///
    /// 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* create_gsl_matrix_copy(void) const;

    gsl_matrix* m_;
    gsl_matrix_view* view_;
  };

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


}} // of namespace gslapi and namespace theplu

#endif