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

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: Matrix.h 3999 2020-10-08 23:22:32Z 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
  Copyright (C) 2012, 2017, 2018, 2019, 2020 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 <cassert>

#include "config_public.h"

// include these operations as they are expected when using Matrix
#include "BLAS_level2.h"
#include "BLAS_level3.h"

#include "Exception.h"
#include "MatrixExpression.h"
#include "StrideIterator.h"
#include "Vector.h"
#include "VectorConstView.h"
#include "VectorExpression.h"
#include "VectorView.h"
#include "yat_assert.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 BasicMatrix<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.

       If \a r is zero \a c must be zero and vice versa.

       \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&);

    /**
       Constructor from a matrix expression. A matrix expression is
       result \c operator+, \c operator-, and operator*, or
       combinations of them.

       A typical usage looks like

       \code
       Matrix A( B * C + D)
       \endcode

       where B, C, and D are all instances of class Matrix.

       Typically this constructor is not used, if rvalues are
       enabled. \see Matrix(MatrixExpression<T>&&).

       \since new in yat 0.15
     */
    template<class T>
    Matrix(const MatrixExpression<T>& other);

    /**
       \brief The move constructor.

       \since new in yat 0.15
    */
    Matrix(Matrix&&) noexcept;

    /**
       Move constructor from a matrix expression. A matrix expression is
       result \c operator+, \c operator-, and operator*, or
       combinations of them.

       A typical usage looks like

       \code
       Matrix A( B * C + D)
       \endcode

       where B, C, and D are all instances of class Matrix.

       \since new in yat 0.15
    */
    template<class T>
    Matrix(MatrixExpression<T>&& other)
      : blas_result_(NULL), m_(NULL)
    {
      other.move(m_);
    }


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

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

       If \a r is zero \a c must be zero and vice versa.

       \note underlying GSL matrix is destroyed and views and
       iterators are invalidated
    */
    void resize(size_t r, size_t c, 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.

       \note Unless matrix is square, views and iterators are
       invalidated.

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

       \note If lhs needs to be resized, views and iterators are
       invalidated.

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

    /**
       Assignment from a matrix expression. A matrix expression is the
       result of \c operator+, \c operator-, and operator*, or
       combinations of them.

       A typical usage looks like

       \code
       A = B * C + D;
       \endcode

       where A, B, C, and D are all instances of class Matrix.

       Typically this constructor is not used, if rvalues are
       enabled. \see Matrix(MatrixExpression<T>&&).

       \since new in yat 0.15
     */
    template<class T>
    Matrix& operator=(const MatrixExpression<T>& rhs);


    /**
       \brief Move assignment operator

       \since new in yat 0.15
     */
    Matrix& operator=(Matrix&& other);

    /**
       Move assignment from a matrix expression. A matrix expression
       is the result of \c operator+, \c operator-, and operator*, or
       combinations of them.

       A typical usage looks like

       \code
       A = B * C + D;
       \endcode

       where B, C, and D are all instances of class Matrix.

       Typically this constructor is not used, if rvalues are
       enabled. \see Matrix(MatrixExpression<T>&&).

       \since new in yat 0.15
     */
    template<class T>
    Matrix& operator=(MatrixExpression<T>&& rhs)
    {
      // This function steals data from rhs into this, and give the old
      // m_ in return to rhs so destructor of rhs takes care of
      // deallocation.

      rhs.move(m_);
      return *this;
    }


    /**
       \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 Addition and assign operator.

       \return A reference to resulting Matrix
    */
    template<class T>
    Matrix& operator+=(const MatrixExpression<T>& rhs)
    {
      *this = *this + rhs;
      return *this;
    }

    /**
       \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 Subtraction and assign operator.

       \return A reference to the resulting VectorMutable.
    */
    template<class T>
    Matrix& operator-=(const MatrixExpression<T>& rhs)
    {
      *this = *this - rhs;
      return *this;
    }

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

       \note invalidates views and iterators

       \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_;
  };

  /**
     Calculate inverse of a (square) Matrix using SVD

     \since New in yat 0.15

     \relates Matrix
   */
  void inverse_svd(const Matrix& input, Matrix& result);

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

     This function swaps element by element and matrices must have the
     same size. References and iterators are not invalidated.

     If you use c++11, do not care about iterators and references,
     std::swap is typically faster (only swapping a couple of pointers)

     \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 Trace of matrix

     \relates Matrix

     \return sum of diagonal elements

     \since New in yat 0.9
  */
  double trace(const Matrix&);

  /////  template implemantions
  template<class T>
  Matrix::Matrix(const MatrixExpression<T>& other)
    : blas_result_(NULL), m_(NULL)
  {
    // In principle, we could implement this as a move constructor
    // and steal the data (rather than copy), but it's a bit hackish
    // to workaround the constness, so users who are eager for speed
    // should enable rvalues and get access to the faster move
    // constructor below.

    try {
      detail::copy(m_, other.gsl_matrix_p());
    }
    catch (GSL_error& e) {
      detail::deallocate(m_);
      throw e;
    }
  }


  template<class T>
  Matrix& Matrix::operator=(const MatrixExpression<T>& rhs)
  {
    // Be careful because *this might be embedded in rhs.
    rhs.get(blas_result_);
    std::swap(m_, blas_result_);
    return *this;
  }


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

#endif