BLAS_level3.h

#ifndef _theplu_yat_utility_blas_level3
#define _theplu_yat_utility_blas_level3

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

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

#include "expression/MatrixBinary.h"
#include "expression/MatrixProduct.h"
#include "expression/ScaledMatrix.h"
#include "expression/TransposedMatrix.h"

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

namespace theplu {
namespace yat {
namespace utility {

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

  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::MatrixProduct<BasicMatrix<Derived1>, BasicMatrix<Derived2> >
  operator*(const BasicMatrix<Derived1>& lhs, const BasicMatrix<Derived2>& rhs)
  {
    YAT_ASSERT(lhs.columns() == rhs.rows());
    return expression::MatrixProduct<BasicMatrix<Derived1>,
                                     BasicMatrix<Derived2>
                                     >(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<class T>
  expression::ScaledMatrix<T>
  operator*(const BasicMatrix<MatrixExpression<expression::ScaledMatrix<T> > >& A,
            double k)
  {
    const expression::ScaledMatrix<T>& sm =
      static_cast<const expression::ScaledMatrix<T>&>(A);
    return expression::ScaledMatrix<T>(k*sm.factor(), sm.base());
  }


  template<class T>
  expression::ScaledMatrix<T>
  operator*(double k,
            const BasicMatrix<MatrixExpression<expression::ScaledMatrix<T> > >& A)
  {
    const expression::ScaledMatrix<T>& sm =
      static_cast<const expression::ScaledMatrix<T>&>(A);
    return expression::ScaledMatrix<T>(k*sm.factor(), sm.base());
  }


  template<class T>
  expression::ScaledMatrix<T>
  operator-(const BasicMatrix<MatrixExpression<expression::ScaledMatrix<T> > >& A)
  {
    const expression::ScaledMatrix<T>& sm =
      static_cast<const expression::ScaledMatrix<T>&>(A);
    return expression::ScaledMatrix<T>(-sm.factor(), sm.base());
  }


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

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

#endif