BLAS_level3.h

#ifndef _theplu_yat_utility_blas_level3
#define _theplu_yat_utility_blas_level3

// $Id: BLAS_level3.h 3654 2017-07-10 05:36:55Z peter $

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

#include <gsl/gsl_matrix.h>

namespace theplu {
namespace yat {
namespace utility {

  // This file defines operations using both Matrix (but not Vector)


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

      void calculate_matrix(gsl_matrix*& result, Multiplies) const
      {
        YAT_ASSERT(detail::rows(result) == this->rows());
        YAT_ASSERT(detail::columns(result) == this->columns());
        YAT_ASSERT(lhs_.columns() == rhs_.rows());
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
                       lhs_.gsl_matrix_p(), rhs_.gsl_matrix_p(),0.0, result);
        YAT_ASSERT(result);
      }

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


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


    template<class T>
    class ScaledMatrix : public MatrixExpression<ScaledMatrix<T> >
    {
    public:
      ScaledMatrix(double factor, const BasicMatrix<T>& A)
        : A_(A), factor_(factor) {}

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


      double operator()(size_t i, size_t j) const
      {
        return factor_ * A_(i, j);
      }


      void calculate_matrix(gsl_matrix*& result) const
      {
        detail::copy(result, A_.gsl_matrix_p());
        gsl_matrix_scale(result, factor_);
      }

    private:
      const BasicMatrix<T>& A_;
      double factor_;
    };


    template<class T>
    class TransposedMatrix : public MatrixExpression<TransposedMatrix<T> >
    {
    public:
      TransposedMatrix(const BasicMatrix<T>& A)
        : A_(A) {}

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

      double operator()(size_t i, size_t j) const
      {
        return A_(j, i);
      }

      void calculate_matrix(gsl_matrix*& result) const
      {
        detail::reallocate(result, rows(), columns());
        gsl_matrix_transpose_memcpy(result, A_.gsl_matrix_p());
      }

    private:
      const BasicMatrix<T>& A_;
    };


  } // end namespace expression


  template<class Derived1, class Derived2>
  expression::MatrixBinary<Derived1, Derived2, expression::Plus>
  operator+(const BasicMatrix<Derived1>& lhs, const BasicMatrix<Derived2>& rhs)
  {
    return expression::MatrixBinary<Derived1, Derived2,
                                    expression::Plus>(lhs, rhs);
  }


  template<class Derived1, class Derived2>
  expression::MatrixBinary<Derived1, Derived2, expression::Minus>
  operator-(const BasicMatrix<Derived1>& lhs, const BasicMatrix<Derived2>& rhs)
  {
    return expression::MatrixBinary<Derived1, Derived2,
                                    expression::Minus>(lhs, rhs);
  }


  template<class Derived1, class Derived2>
  expression::MatrixBinary<Derived1, Derived2, expression::Multiplies>
  operator*(const BasicMatrix<Derived1>& lhs, const BasicMatrix<Derived2>& rhs)
  {
    YAT_ASSERT(lhs.columns() == rhs.rows());
    return expression::MatrixBinary<Derived1, Derived2,
                                    expression::Multiplies>(lhs, rhs);
  }


  template<class T>
  expression::ScaledMatrix<T>
  operator*(const BasicMatrix<T>& A, double k)
  {
    return expression::ScaledMatrix<T>(k, A);
  }


  template<class T>
  expression::ScaledMatrix<T>
  operator*(double k, const BasicMatrix<T>& A)
  {
    return expression::ScaledMatrix<T>(k, A);
  }


  template<typename T>
  expression::ScaledMatrix<T>
  operator-(const BasicMatrix<T>& m)
  {
    return expression::ScaledMatrix<T>(-1.0, m);
  }


  template<typename T>
  expression::TransposedMatrix<T> transpose(const BasicMatrix<T>& A)
  {
    return expression::TransposedMatrix<T>(A);
  }

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

#endif