MatrixBinary.h

#ifndef _theplu_yat_utility_expression_matrix_binary
#define _theplu_yat_utility_expression_matrix_binary

// $Id: MatrixBinary.h 3908 2020-05-13 06:58:38Z peter $

/*
  Copyright (C) 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 "yat/utility/BasicMatrix.h"
#include "yat/utility/BLAS_utility.h"
#include "yat/utility/MatrixExpression.h"
#include "yat/utility/yat_assert.h"

#include <gsl/gsl_matrix.h>

namespace theplu {
namespace yat {
namespace utility {


namespace expression {



  template<typename LHS, typename RHS, class OP>
  class MatrixBinary
    : public MatrixExpression<MatrixBinary<LHS, RHS, OP> >
  {
  public:
    MatrixBinary(const BasicMatrix<LHS>& lhs, const BasicMatrix<RHS>& rhs)
      : lhs_(lhs), rhs_(rhs)
    {
    }

    size_t rows(void) const { return lhs_.rows(); }
    size_t columns(void) const { return rhs_.columns(); }

    double operator()(size_t row, size_t column) const
    { return get(row, column, op_); }

    void calculate_matrix(gsl_matrix*& result) const
    {
      detail::reallocate(result, this->rows(), this->columns());
      calculate_matrix(result, op_);
    }

  private:
    const BasicMatrix<LHS>& lhs_;
    const BasicMatrix<RHS>& rhs_;
    OP op_;

    template<class T>
    void calculate_matrix(gsl_matrix*& result, T) const
    {
      YAT_ASSERT(detail::rows(result) == this->rows());
      YAT_ASSERT(detail::columns(result) == this->columns());
      for (size_t i=0; i<rows(); ++i)
        for (size_t j=0; j<columns(); ++j)
          gsl_matrix_set(result, i, j, (*this)(i, j));
    }


    double get(size_t row, size_t column, Plus) const
    {
      return lhs_(row, column) + rhs_(row, column);
    }


    double get(size_t row, size_t column, Minus) const
    {
      return lhs_(row, column) - rhs_(row, column);
    }
  };


} // end namespace expression


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

#endif