BLAS_level1.h

#ifndef _theplu_yat_utility_blas_level1
#define _theplu_yat_utility_blas_level1

// $Id: BLAS_level1.h 3678 2017-08-17 05:07:11Z 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 "BasicVector.h"
#include "BLAS_utility.h"
#include "VectorExpression.h"

#include "yat_assert.h"

#include "gsl/gsl_vector.h"

#include <cstddef>

namespace theplu {
namespace yat {
namespace utility {

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


  namespace expression {

    template<typename LHS, typename RHS, class OP>
    class VectorBinary
      : public VectorExpression<VectorBinary<LHS, RHS, OP> >
    {
    public:
      VectorBinary(const BasicVector<LHS>& lhs, const BasicVector<RHS>& rhs)
        : lhs_(lhs), rhs_(rhs)
      {
        YAT_ASSERT(lhs.size() == rhs.size());
      }

      size_t size(void) const { return lhs_.size(); }

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

      void calculate_gsl_vector_p(void) const
      {
        this->allocate_memory(size());
        for (size_t i=0; i<size(); ++i)
          gsl_vector_set(this->v_, i, get(i, op_));
      }

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

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


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


    template<class T>
    class ScaledVector : public VectorExpression<ScaledVector<T> >
    {
    public:
      ScaledVector(double factor, const BasicVector<T>& vec)
        : vec_(vec), factor_(factor) {}

      size_t size(void) const
      {
        return vec_.size();
      }


      double operator()(size_t i) const
      {
        return factor_ * vec_(i);
      }


      void calculate_gsl_vector_p(void) const
      {
        YAT_ASSERT(this->v_ == NULL);
        this->allocate_memory(size());
        gsl_vector_memcpy(this->v_, vec_.gsl_vector_p());
        gsl_vector_scale(this->v_, factor_);
      }

    private:
      const BasicVector<T>& vec_;
      double factor_;
    };


  } // end namespace expression



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


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


  template<class Derived1, class Derived2>
  double operator*(const BasicVector<Derived1>& lhs,
                   const BasicVector<Derived2>& rhs)
  {
    double res = 0.0;;
    gsl_blas_ddot(lhs.gsl_vector_p(), rhs.gsl_vector_p(), &res);
    return res;
  }


  // TODO specialization for when vec is ScaledVector
  template<typename T>
  expression::ScaledVector<T>
  operator*(double k, const BasicVector<T>& vec)
  {
    return expression::ScaledVector<T>(k, vec);
  }


  template<typename T>
  expression::ScaledVector<T>
  operator*(const BasicVector<T>& vec, double k)
  {
    return expression::ScaledVector<T>(k, vec);
  }


  template<typename T>
  expression::ScaledVector<T>
  operator-(const BasicVector<T>& vec)
  {
    return expression::ScaledVector<T>(-1.0, vec);
  }

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

#endif