source: trunk/yat/utility/Matrix.h @ 1887

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: Matrix.h 1887 2009-03-31 17:38:16Z peter $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2007, 2008, 2009 Jari Häkkinen, Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "Exception.h"
#include "StrideIterator.h"
#include "Vector.h"
#include "VectorConstView.h"
#include "VectorView.h"

#include <gsl/gsl_matrix.h>

#include <cstddef> // size_t
#include <iosfwd>

namespace theplu {
namespace yat {
namespace utility {

  class VectorBase;

  ///
  /// @brief Interface to GSL matrix.
  ///
  /// For the time being 'double' is the only type supported.
  ///
  /// \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.
  ///
  /// @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,
  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
  /// superdiagonals.
  ///
  class Matrix
  {
  public:
    /**
       value_type is double

       \since New in yat 0.5
     */
    typedef double value_type;

    /**
       reference type is double&

       \since New in yat 0.5
     */
    typedef double& reference;

    /**
       const_reference type is const double&

       \since New in yat 0.5
     */
    typedef const double& const_reference;

    /**
       Mutable iterator that iterates over all elements
     */
    typedef StrideIterator<double*> iterator;

    /**
       Read-only iterator that iterates over all elements
     */
    typedef StrideIterator<const double*> const_iterator;

    /**
       Mutable iterator that iterates over one column
     */
    typedef StrideIterator<double*> column_iterator;

    /**
       Read-only iterator that iterates over one column
     */
    typedef StrideIterator<const double*> const_column_iterator;

    /**
       Mutable iterator that iterates over one row
     */
    typedef StrideIterator<double*> row_iterator;

    /**
       Read-only iterator that iterates over one row
     */
    typedef StrideIterator<const double*> const_row_iterator;

    /**
       @brief The default constructor.
    
       This constructor does not initialize underlying (essential)
       structures.
    */
    Matrix(void);

    /**
       \brief Constructor allocating memory space for \a r times \a c
       elements, and sets all elements to \a init_value.

       \throw GSL_error if memory allocation fails.
    */
    Matrix(const size_t& r, const size_t& c, double init_value=0);

    /**
       \brief The copy constructor.

       \throw A GSL_error is indirectly thrown if memory allocation
       fails.
    */
    Matrix(const Matrix&);

    /**
       \brief The istream constructor.

       The std::istream will be interpreted as outlined here:

       Missing values, i.e. empty elements, are treated as NaN values
       (std::numeric_limits<double>::quiet_NaN() to be specific).

       Matrix rows are separated with the new line character.

       Column element separation has two modes depending on the value
       of \a sep.

       - If \a sep is the default '\\0' value then column elements are
       separated with white space characters except the new line
       character. Multiple sequential white space characters are
       treated as one separator.

       - Setting \a sep to something else than the default value will
       change the behaviour to use the \a sep character as the
       separator between column elements. Multiple sequential \a sep
       characters will be treated as separating elements with missing
       values.

       End of input is the end of file marker and this treatment
       cannot be redefined using the provided API.

       \throw GSL_error if memory allocation fails, IO_error if
       unexpected input is found in the input stream.
    */
    explicit Matrix(std::istream &, char sep='\0') 
      throw(utility::IO_error, std::exception);

    /**
       \brief The destructor.
    */
    ~Matrix(void);

    ///
    /// Set all elements to \a value.
    ///
    void all(const double value);

    /**
       Iterator iterates along a row. When end of row is reached it
       jumps to beginning of next row.

       \return iterator pointing to upper-left element.
     */
    iterator begin(void);

    /**
       Iterator iterates along a row. When end of row is reached it
       jumps to beginning of next row.

       \return const_iterator pointing to upper-left element.
     */
    const_iterator begin(void) const;

    /**
       Iterator iterates along a column.

       \return iterator pointing to first element of column \a i.
     */
    iterator begin_column(size_t i);

    /**
       Iterator iterates along a column.

       \return const_iterator pointing to first element of column \a i.
     */
    const_iterator begin_column(size_t i) const;

    /**
       Iterator iterates along a row.

       \return iterator pointing to first element of row \a i.
     */
    iterator begin_row(size_t i);

    /**
       Iterator iterates along a row.

       \return const_iterator pointing to first element of row \a i.
     */
    const_iterator begin_row(size_t i) const;

    /**
       \return Vector view of column \a i
     */
    VectorView column_view(size_t i);

    /**
       \return const vector view of column \a i
     */
    const VectorConstView column_const_view(size_t) const;

    ///
    /// @return The number of columns in the matrix.
    ///
    size_t columns(void) const;

    /**
       Elementwise division of the elements 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.

       \throw GSL_error if dimensions mismatch.
    */
    void div(const Matrix& b);

    /**
       \return iterator pointing to end of matrix
     */
    iterator end(void);

    /**
       \return const_iterator pointing to end of matrix
     */
    const_iterator end(void) const;

    /**
       \return iterator pointing to end of column \a i
     */
    iterator end_column(size_t i);

    /**
       \return const_iterator pointing to end of column \a i
     */
    const_iterator end_column(size_t i) const;

    /**
       \return iterator pointing to end of row \a i
     */
    iterator end_row(size_t i);

    /**
       \return const_iterator pointing to end of row \a i
     */
    const_iterator end_row(size_t i) const;

    /**
       \brief 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.

       \see operator== and operator!=
    */
    bool equal(const Matrix&, const double precision=0) const;

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

    ///
    /// @return A pointer to the internal GSL matrix.
    ///
    gsl_matrix* gsl_matrix_p(void);

    /**
       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.

       \throw GSL_error if dimensions mismatch.
    */
    void mul(const Matrix& b);

    /**
       \brief Resize Matrix
       
       All elements are set to @a init_value.

       \note underlying GSL matrix is destroyed and views into this
       Matrix becomes invalid.
    */
    void resize(size_t, size_t, double init_value=0);

    ///
    /// @return The number of rows in the matrix.
    ///
    size_t rows(void) const;

    /**
       \return Vector view of \a i
     */
    VectorView row_view(size_t);

    /**
       \return const vector view of \a i
     */
    const VectorConstView row_const_view(size_t) const;

    /**
       \brief Swap columns \a i and \a j.

       \throw GSL_error if either index is out of bounds.
    */
    void swap_columns(const size_t i,const size_t j);

    /**
       \brief Swap row \a i and column \a j.

       The Matrix must be square.

       \throw GSL_error if either index is out of bounds, or if matrix
       is not square.
    */
    void swap_rowcol(const size_t i,const size_t j);

    /**
       \brief Swap rows \a i and \a j.

       \throw GSL_error if either index is out of bounds.
    */
    void swap_rows(const size_t i, const size_t j);

    /**
       \brief Transpose the matrix.

       \throw GSL_error if memory allocation fails for the new
       transposed matrix.
    */
    void transpose(void);

    /**
       \brief Element access operator.

       \return Reference to the element position (\a row, \a column).

       \throw If GSL range checks are enabled in the underlying GSL
       library a GSL_error exception is thrown if either index is out
       of range.
    */
    double& operator()(size_t row,size_t column);

    /**
       \brief Element access operator.

       \return Const reference to the element position (\a row, \a
       column).

       \throw If GSL range checks are enabled in the underlying GSL
       library a GSL_error exception is thrown if either index is out
       of range.
    */
    const double& operator()(size_t row,size_t column) const;

    /**
       \brief Comparison operator. Takes squared time.

       Checks are performed with exact matching, i.e., rounding off
       effects may destroy comparison. Use the equal function for
       comparing elements within a user defined precision.

       \return True if all elements are equal otherwise false.

       \see equal
    */
    bool operator==(const Matrix& other) const;

    /**
       \brief Comparison operator. Takes squared time.

       Checks are performed with exact matching, i.e., rounding off
       effects may destroy comparison. Use the equal function for
       comparing elements within a user defined precision.

       \return False if all elements are equal otherwise true.

       \see equal
    */
    bool operator!=(const Matrix& other) const;

    /**
       \brief The assignment operator.

       \return A const reference to the resulting Matrix.
    */
    const Matrix& operator=(const Matrix& other);

    /**
       \brief Add and assign operator.

       Elementwise addition of the elements of Matrix \a b to the left
       hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
       \f$.

       \return A const reference to the resulting Matrix.

       \throw GSL_error if dimensions mismatch.
    */
    const Matrix& operator+=(const Matrix& b);

    /**
       \brief Add and assign operator

       Add the scalar value \a d to the left hand side Matrix, \f$
       a_{ij} = a_{ij} + d \; \forall i,j \f$.
    */
    const Matrix& operator+=(const double d);

    /**
       \brief Subtract and assign operator.

       Elementwise subtraction of the elements of Matrix \a b to the
       left hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
       i,j \f$.

       \return A const reference to the resulting Matrix.

       \throw GSL_error if dimensions mismatch.
    */
    const Matrix& operator-=(const Matrix&);

    /**
       \brief Subtract and assign operator

       Subtract the scalar value \a d to the left hand side Matrix,
       \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
    */
    const Matrix& operator-=(const double d);

    /**
       \brief Multiply and assignment operator.

       \return Const reference to the resulting Matrix.

       \throw GSL_error if memory allocation fails.
    */
    const Matrix& operator*=(const Matrix&);

    /**
       \brief Multiply and assignment operator

       Multiply the elements of the left hand side Matrix with a
       scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.

       \throw GSL_error if memory allocation fails.
    */
    const Matrix& operator*=(double d);

  private:

    /**
       \brief 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.

       \throw GSL_error if memory cannot be allocated for the new
       copy, or if dimensions mismatch.
    */
    gsl_matrix* create_gsl_matrix_copy(void) const;

    /**
       \brief Clear all dynamically allocated memory.

       Internal utility function.
    */
    void delete_allocated_memory(void);

    // blas_result_ is used to temporarily store result in BLAS calls
    // and when not NULL it should always have the same size as
    // m_. Memory is not allocated for blas_result_ until it is
    // needed.
    gsl_matrix* blas_result_;
    gsl_matrix* m_;
  };

  /**
     \brief Check if all elements of the Matrix are zero.

     \return True if all elements in the Matrix is zero, false
     otherwise.

     \relates Matrix
  */
  bool isnull(const Matrix&);

  /**
     \brief Get the maximum value of the Matrix.

     \return The maximum value of the Matrix.

     \relates Matrix
  */
  double max(const Matrix&);

  /**
     \brief Get the minimum value of the Matrix.

     \return The minimum value of the Matrix.

     \relates Matrix
  */
  double min(const Matrix&);

  /**
     \brief Locate the maximum and minimum element in the Matrix.

     \return The indices to the element with the minimum and maximum
     values of the Matrix, respectively.

     \note Lower index has precedence (searching in row-major order).

     \relates Matrix
  */
  void minmax_index(const Matrix&,
                    std::pair<size_t,size_t>& min,
                    std::pair<size_t,size_t>& max);

  /**
     \brief Create a Matrix \a flag indicating NaN's in another Matrix
     \a templat.

     The \a flag Matrix is changed to contain 1's and 0's only. A 1
     means that the corresponding element in the \a templat Matrix is
     valid and a zero means that the corresponding element is a NaN.

     \note Space for Matrix \a flag is reallocated to fit the size of
     Matrix \a templat if sizes mismatch.

     \return True if the \a templat Matrix contains at least one NaN.

     \relates Matrix
  */
  bool nan(const Matrix& templat, Matrix& flag);

  /**
     \brief Exchange all elements between the matrices by copying.

     The two matrices must have the same size.

     \throw GSL_error if sizes are not equal.

     \relates Matrix
  */
  void swap(Matrix&, Matrix&);

  /**
     \brief The output operator for the Matrix class.

     \relates Matrix
  */
  std::ostream& operator<< (std::ostream& s, const Matrix&);

  /**
     \brief Vector Matrix multiplication

     \relates Matrix
   */
  Vector operator*(const Matrix&, const VectorBase&);

  /**
     \brief Matrix Vector multiplication

     \relatesalso Matrix

     \relatesalso VectorBase
   */
  Vector operator*(const VectorBase&, const Matrix&);

}}} // of namespace utility, yat, and theplu

#endif