source: branches/mtl-test/yat/utility/mtlmatrix.h @ 987

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: mtlmatrix.h 987 2007-10-26 20:31:39Z 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

  This file is part of the yat library, http://trac.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;

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

    /**
       \brief Make a copy of \a other.

       This function will make a deep copy of \a other. Memory is
       resized and view state is changed to same as in other.
    */
    const matrix& clone(const matrix& other);

    ///
    /// @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);

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

    ///
    /// @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;

    /**
       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 a view 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;

    ///
    /// Set all elements to \a value.
    ///
    void 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 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 row(const size_t row, const vector& vec);

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

       Dimensions of the matrices must match. If the LHS matrix is a
       view, the underlying data will be changed.

       \return A const reference to the resulting matrix.

       \see void set(const matrix&).

       \throw GSL_error if dimensions mis-match.
    */
    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_;
    gsl_matrix_view* view_;
    const gsl_matrix_const_view* view_const_;
    // proxy_m_ is used to access the proper underlying gsl_matrix in
    // all const member functions. It is not used by design for
    // non-const vector functions and operators. This is to make sure
    // that runtime errors occur if a const matrix is used in an
    // inappropriate manner such as on left hand side in assignment
    // (remember, m_ is null for const matrix views).
    const gsl_matrix* proxy_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 vector&);

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

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

#endif