Matrix.h

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: Matrix.h 4207 2022-08-26 04:36:28Z 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, 2021, 2022 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 "Exception.h"
#include "MatrixExpression.h"
#include "MatrixMutable.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;


  namespace detail {

    /*
      Class involved in move assignment of MatrixMutable to determine
      when both lhs and rhs are a Matrix. When the behavoiur only
      depends on the type of one object, dynamic polymorphism with a
      virtual function suffices; when the behaviour depends on two
      type of two objects, we need to simulate double dispatch and
      this class is helping the Matrix and MatrixMutable to implement
      that.

      In brief, it's constructed inside lhs (being a Matrix), taking
      and storing a reference to lhs. It then calls operator() taking
      rhs. It then uses the Visitor pattern calling rhs::visit(Mover),
      where visit is a virtual function having a special implemntation
      for Matrix where it's using fast move assignment, whereas for
      other rhs types move assignment can't be used (because the rhs
      doesn't own its data) and we use copy assignment instead.
     */
    class Mover
    {
    public:
      Mover(Matrix& lhs);
      void operator()(MatrixMutable& rhs);
      Matrix& lhs(void);
      void copy_assign(const MatrixMutable& rhs);
    private:
      Matrix& lhs_;
    };
  }

  class Matrix : public MatrixMutable
  {
  public:
    Matrix(void);

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

    Matrix(const Matrix&);

    explicit Matrix(const MatrixBase&);

    virtual ~Matrix(void);

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

    Matrix(Matrix&&) noexcept;

    Matrix(MatrixMutable&&);

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


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

    gsl_matrix* gsl_matrix_p(void);

    const gsl_matrix* gsl_matrix_p(void) const;

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

    void transpose(void);

    const Matrix& operator=(const Matrix& other);
    using MatrixMutable::operator=; // access to =(const MatrixBase&)

    Matrix& operator=(Matrix&& other);

    // to get access to operator*=(double)
    using MatrixMutable::operator*=;

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

  protected:
    void copy_assign(const gsl_matrix* rhs);
    using MatrixMutable::copy_assign;

    void move_assign(MatrixMutable&& rhs);

    void move_assign(gsl_matrix*&& rhs);

  private:
    gsl_matrix* m_;

    gsl_matrix* create_gsl_matrix_copy(void) const;

    // Used in move_assignment, see docs for Mover class.
    friend class detail::Mover;
    void visit(detail::Mover& mover);
  };

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

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

  template<class T>
  Matrix::Matrix(const MatrixExpression<T>& other)
  {
    // 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;
    }
  }

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

#endif