MatrixProduct.h

#ifndef _theplu_yat_utility_expression_matrix_product
#define _theplu_yat_utility_expression_matrix_product

// $Id: MatrixProduct.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 "MatrixTraits.h"
#include "ScaledMatrix.h"
#include "TransposedMatrix.h"

#include "yat/utility/BasicMatrix.h"
#include "yat/utility/MatrixExpression.h"
#include "yat/utility/yat_assert.h"

#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>

namespace theplu {
namespace yat {
namespace utility {


namespace expression {


  template<class MATRIX1, class MATRIX2>
  class MatrixProduct : public MatrixExpression<MatrixProduct<MATRIX1,MATRIX2> >
  {
  public:
    MatrixProduct(const MATRIX1& A, const MATRIX2& B)
      : A_(A), B_(B) {}

    size_t rows(void) const { return A_.rows(); }
    size_t columns(void) const { return B_.columns(); }

    double operator()(size_t row, size_t column) const
    { return gsl_matrix_get(this->gsl_matrix_p(), row, column); }

    void calculate_matrix(gsl_matrix*& result) const
    {
      detail::reallocate(result, this->rows(), this->columns());
      YAT_ASSERT(detail::rows(result) == this->rows());
      YAT_ASSERT(detail::columns(result) == this->columns());
      YAT_ASSERT(A_.columns() == B_.rows());
      MatrixTraits<MATRIX1> traits1;
      MatrixTraits<MATRIX2> traits2;
      gsl_blas_dgemm(traits1.transpose_type(), traits2.transpose_type(),
                     traits1.factor(A_) * traits2.factor(B_),
                     traits1.get_gsl_matrix_p(A_),
                     traits2.get_gsl_matrix_p(B_),
                     0.0, result);
      YAT_ASSERT(result);
    }
  private:
    const MATRIX1& A_;
    const MATRIX2& B_;
  };
} // end namespace expression


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

#endif