BLAS_level2.h

#ifndef _theplu_yat_utility_blas_level2
#define _theplu_yat_utility_blas_level2

// $Id: BLAS_level2.h 3999 2020-10-08 23:22:32Z peter $

/*
  Copyright (C) 2017, 2019, 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 "BasicVector.h"
#include "BLAS_level3.h"
#include "VectorExpression.h"
#include "yat_assert.h"

#include "expression/MatrixTraits.h"

#include "gsl/gsl_blas.h"

#include <cstddef>

namespace theplu {
namespace yat {
namespace utility {

  // This file defines operations using both Vector and Matrix



  namespace expression {

    template<typename MATRIX, typename VECTOR>
    class MatrixVector : public VectorExpression<MatrixVector<MATRIX, VECTOR> >
    {
    public:

      MatrixVector(const BasicMatrix<MATRIX>& lhs,
                   const BasicVector<VECTOR>& rhs)
      {
        YAT_ASSERT(lhs.columns() == rhs.size());
        YAT_ASSERT(lhs.rows());
        this->allocate_memory(lhs.rows());

        MatrixTraits<BasicMatrix<MATRIX> > traits;

        gsl_blas_dgemv(traits.transpose_type(), traits.factor(lhs),
                       traits.get_gsl_matrix_p(lhs),
                       rhs.gsl_vector_p(), 0.0, this->v_);
      }


      // access an element
      double operator()(size_t i) const
      {
        YAT_ASSERT(this->v_);
        return gsl_vector_get(this->v_, i);
      }


      size_t size(void) const
      {
        YAT_ASSERT(this->v_);
        return this->v_->size;
      }


      void calculate_gsl_vector_p(void) const
      {
        // This should never be called since v_ is constructed in
        // constructor
        YAT_ASSERT(this->v_);
        YAT_ASSERT(0);
      }
    };

  } // end namespace expression



  template<class MATRIX, class VECTOR>
  expression::MatrixVector<MATRIX, VECTOR>
  operator*(const BasicMatrix<MATRIX>& lhs, const BasicVector<VECTOR>& rhs)
  {
    YAT_ASSERT(lhs.columns() == rhs.size());
    return expression::MatrixVector<MATRIX, VECTOR>(lhs, rhs);
  }


  template<class MATRIX, class VECTOR>
  expression::MatrixVector<MatrixExpression<
                             expression::TransposedMatrix<MATRIX> >,
                           VECTOR>
  operator*(const BasicVector<VECTOR>& lhs, const BasicMatrix<MATRIX>& rhs)
  {
    YAT_ASSERT(lhs.size() == rhs.rows());
    return transpose(rhs) * lhs;
  }


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

#endif