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

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: matrix.h 717 2006-12-25 20:30:00Z 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

  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.
    ///
    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.
    ///
    matrix(const matrix&);

    ///
    /// 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 &, 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.
    ///
    /// @return Whatever GSL returns.
    ///
    int 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.
    ///
    /// @return Whatever GSL returns.
    ///
    int add_constant(const 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.
    ///
    /// @return Whatever GSL returns.
    ///
    int 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.
    ///
    // Jari, is this needed?
    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.
    ///
    /// @return Whatever GSL returns.
    ///
    int
    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.
    ///
    /// @return Whatever GSL returns.
    ///
    int scale(const double 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&)
    ///
    int set(const matrix& mat);

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

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

    ///
    /// @brief Set \a row values to values in \a vec.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    int 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.
    ///
    /// @return Whatever GSL returns.
    ///
    int sub(const matrix& b);

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

    ///
    /// @brief Swap columns \a i and \a j.
    ///
    int swap_columns(const size_t i,const size_t j);

    ///
    /// @brief Swap row \a i and column \a j.
    ///
    int swap_rowcol(const size_t i,const size_t j);

    ///
    /// @brief Swap rows \a i and \a j.
    ///
    int swap_rows(const size_t i, const size_t j);

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

    ///
    /// @return Reference to the element position (\a row, \a column).
    ///
    double& operator()(size_t row,size_t column);

    ///
    /// @return Const reference to the element position (\a row, \a
    /// column).
    ///
    const double& operator()(size_t row,size_t column) const;

    ///
    /// @brief Matrix-vector multiplication.
    ///
    /// @return The resulting vector.
    ///
    // Jari, where should this go?
    const vector operator*(const vector&) 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 int set(const matrix&)
    ///
    const matrix& operator=(const matrix& other);

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

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

    ///
    /// @brief Subtract and assign operator.
    ///
    const matrix& operator-=(const matrix&);

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

    ///
    /// @brief Multiply and assignment operator
    ///
    const matrix& operator*=(const double);

  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 utility, yat, and theplu

#endif