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

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: matrix.h 759 2007-02-19 19:41:25Z jari $

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

  This file is part of the yat library, http://lev.thep.lu.se/trac/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 "vector.h"
#include "Exception.h"

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

namespace theplu {
namespace yat {
namespace utility {

  class vector;

  ///
  /// This is the yat 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 utility::matrix'es. A support function is
  /// added, utility::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,
  /// utility::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:

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

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

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

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

       \throw GSL_error if a view cannot be set up.
    */
    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
           size_t n_column);

    /**
       \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);

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

       \throw GSL_error if dimensions mis-match.
    */
    void add(const matrix& b);

    /**
       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.
    */
    void add(double d);

    ///
    /// @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_elements(const matrix& b);

    ///
    /// 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.
    ///
    const gsl_matrix* gsl_matrix_p(void) const;

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

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

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

    ///
    /// @return The maximum value of the matrix.
    ///
    double max(void) const;

    ///
    /// @return The minimum value of the matrix.
    ///
    double min(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).
    ///
    void minmax_index(std::pair<size_t,size_t>& min,
                      std::pair<size_t,size_t>& max) const;

    ///
    /// The argument is changed into a matrix containing 1's and 0's
    /// only. 1 means element in matrix is valid and a zero means
    /// element is a NaN.
    ///
    /// @return true if matrix contains at least one NaN.
    ///
    bool nan(matrix&) 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.

       \throw GSL_error if dimensions mis-match.
    */
    void mul_elements(const matrix& b);

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

    /**
       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.
    */
    void scale(const double d);

    /**
       \brief Set element values to values in \a mat.

       This function is needed for assignment of viewed elements.

       \see const matrix& operator=(const matrix&)

       \throw GSL_error if dimensions mis-match.
    */
    void set(const matrix& mat);

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

    /**
       \brief Set \a column values to values in \a vec.

       \note No check on size is done.

       \throw GSL_error if index is out of range or mis-match in
       sizes.
    */
    void set_column(const size_t column, const vector& vec);

    /**
       \brief Set \a row values to values in \a vec.

       \note No check on size is done.

       \throw GSL_error if index is out of range or mis-match in
       sizes.
    */
    void set_row(const size_t row, const vector& vec);

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

       \throw GSL_error if dimensions mis-match.
    */
    void sub(const matrix& b);

    /**
       \brief Exchange the elements of the this and \a other by copying.

       The two matrices must have the same size.

       \throw GSL_error if either index is out of bounds.
    */
    void swap(matrix& other);

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

       \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.
    ///
    /// @return True if all elements are equal otherwise False.
    ///
    /// @see equal
    ///
    bool operator==(const matrix& other) const;

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

    /**
       \brief 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 void set(const matrix&).

       \throw A GSL_error is indirectly thrown if memory allocation
       fails, or if dimensions mis-match.
    */
    const matrix& operator=(const matrix& other);

    /**
       \brief Add and assign operator.

       \return A const reference to the resulting matrix.

       \throw GSL_error if dimensions mis-match.
    */
    const matrix& operator+=(const matrix&);

    ///
    /// @brief Add and assign operator.
    ///
    const matrix& operator+=(const double);

    /**
       \brief Subtract and assign operator.

       \return A const reference to the resulting matrix.

       \throw GSL_error if dimensions mis-match.
    */
    const matrix& operator-=(const matrix&);

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

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

  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;

    gsl_matrix* m_;
    gsl_matrix_view* view_;
  };

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

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

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

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

#endif