Matrix.h

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: Matrix.h 3855 2020-01-02 01:11:34Z 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 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;

  class Matrix : public BasicMatrix<Matrix>
  {
  public:
    typedef double value_type;

    typedef double& reference;

    typedef const double& const_reference;

    typedef StrideIterator<double*> iterator;

    typedef StrideIterator<const double*> const_iterator;

    typedef StrideIterator<double*> column_iterator;

    typedef StrideIterator<const double*> const_column_iterator;

    typedef StrideIterator<double*> row_iterator;

    typedef StrideIterator<const double*> const_row_iterator;

    Matrix(void);

    Matrix(const size_t& r, const size_t& c, double init_value=0);

    Matrix(const Matrix&);

    template<class T>
    Matrix(const MatrixExpression<T>& other);

#ifdef YAT_HAVE_RVALUE

    Matrix(Matrix&&) noexcept;

    template<class T>
    Matrix(MatrixExpression<T>&& other)
      : blas_result_(NULL), m_(NULL)
    {
      other.move(m_);
    }

#endif

    explicit Matrix(std::istream &, char sep='\0');

    ~Matrix(void);

    void all(const double value);

    iterator begin(void);

    const_iterator begin(void) const;

    iterator begin_column(size_t i);

    const_iterator begin_column(size_t i) const;

    iterator begin_row(size_t i);

    const_iterator begin_row(size_t i) const;

    VectorView column_view(size_t i);

    const VectorConstView column_const_view(size_t) const;

    size_t columns(void) const;

    void div(const Matrix& b);

    iterator end(void);

    const_iterator end(void) const;

    iterator end_column(size_t i);

    const_iterator end_column(size_t i) const;

    iterator end_row(size_t i);

    const_iterator end_row(size_t i) const;

    bool equal(const Matrix&, const double precision=0) const;

    const gsl_matrix* gsl_matrix_p(void) const;

    gsl_matrix* gsl_matrix_p(void);

    void mul(const Matrix& b);

    void resize(size_t r, size_t c, double init_value=0);

    size_t rows(void) const;

    VectorView row_view(size_t);

    const VectorConstView row_const_view(size_t) const;

    void swap_columns(const size_t i,const size_t j);

    void swap_rowcol(const size_t i,const size_t j);

    void swap_rows(const size_t i, const size_t j);

    void transpose(void);

    double& operator()(size_t row,size_t column);

    const double& operator()(size_t row,size_t column) const;

    bool operator==(const Matrix& other) const;

    bool operator!=(const Matrix& other) const;

    const Matrix& operator=(const Matrix& other);

    template<class T>
    Matrix& operator=(const MatrixExpression<T>& rhs);

#ifdef YAT_HAVE_RVALUE

    Matrix& operator=(Matrix&& other);

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

#endif

    const Matrix& operator+=(const Matrix& b);

    template<class T>
    Matrix& operator+=(const MatrixExpression<T>& rhs)
    {
      *this = *this + rhs;
      return *this;
    }

    const Matrix& operator+=(const double d);

    const Matrix& operator-=(const Matrix&);

    template<class T>
    Matrix& operator-=(const MatrixExpression<T>& rhs)
    {
      *this = *this - rhs;
      return *this;
    }

    const Matrix& operator-=(const double d);

    const Matrix& operator*=(const Matrix&);

    const Matrix& operator*=(double d);

  private:

    gsl_matrix* create_gsl_matrix_copy(void) const;

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

  void inverse_svd(const Matrix& input, Matrix& result);

  bool isnull(const Matrix&);

  double max(const Matrix&);

  double min(const Matrix&);

  void minmax_index(const Matrix&,
                    std::pair<size_t,size_t>& min,
                    std::pair<size_t,size_t>& max);

  bool nan(const Matrix& templat, Matrix& flag);

  void swap(Matrix&, Matrix&);

  std::ostream& operator<< (std::ostream& s, const Matrix&);

  double trace(const Matrix&);

  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