Matrix.h

#ifndef _theplu_yat_utility_matrix_
#define _theplu_yat_utility_matrix_

// $Id: Matrix.h 2735 2012-05-29 10:03:15Z 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 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 "Exception.h"
#include "StrideIterator.h"
#include "Vector.h"
#include "VectorConstView.h"
#include "VectorView.h"

#include <gsl/gsl_matrix.h>

#include <cstddef> // size_t
#include <iosfwd>

namespace theplu {
namespace yat {
namespace utility {

  class VectorBase;

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

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

    ~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);

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

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

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

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

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

  Vector operator*(const Matrix&, const VectorBase&);

  Vector operator*(const VectorBase&, const Matrix&);

  double trace(const Matrix&);

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

#endif