source: trunk/yat/utility/matrix.h @ 1098

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: matrix.h 1098 2008-02-18 03:13:53Z jari $

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

  This file is part of the yat library, http://trac.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 2 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 this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

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

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

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:
    /**
     */
    typedef StrideIterator<double*> iterator;

    /**
     */
    typedef StrideIterator<const double*> const_iterator;

    /**
       @brief The default constructor.
    
       This contructor 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.

       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.

       \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 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.

       \throw GSL_error if dimensions mis-match.
    */
    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 mis-match.
    */
    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 mis-match.
    */
    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 mis-match.
    */
    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 assigment 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 mis-match.
    */
    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
     othwerwise.
  */
  bool isnull(const matrix&);

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

     \return The maximum value of the matrix.
  */
  double max(const matrix&);

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

     \return The minimum value of the matrix.
  */
  double min(const matrix&);

  /**
     \brief Locate the maximum and minumum element in the matrix.

     \return The indecies to the element with the minimum and maximum
     values of the matrix, respectively.

     \note Lower index has precedence (searching in row-major order).
  */
  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.
  */
  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 either index is out of bounds.
  */
  void swap(matrix&, matrix&);

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

  /**
     \brief vector matrix multiplication
   */
  vector operator*(const matrix&, const VectorBase&);

  /**
     \brief matrix vector multiplication
   */
  vector operator*(const VectorBase&, const matrix&);

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

#endif