yat/utility/BLAS_level1.h

Code
Comments
Other
Rev Date Author Line
3605 27 Jan 17 peter 1 #ifndef _theplu_yat_utility_blas_level1
3605 27 Jan 17 peter 2 #define _theplu_yat_utility_blas_level1
3605 27 Jan 17 peter 3
3605 27 Jan 17 peter 4 // $Id$
3605 27 Jan 17 peter 5
3605 27 Jan 17 peter 6 /*
3605 27 Jan 17 peter 7   Copyright (C) 2017 Peter Johansson
3605 27 Jan 17 peter 8
3605 27 Jan 17 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3605 27 Jan 17 peter 10
3605 27 Jan 17 peter 11   The yat library is free software; you can redistribute it and/or
3605 27 Jan 17 peter 12   modify it under the terms of the GNU General Public License as
3605 27 Jan 17 peter 13   published by the Free Software Foundation; either version 3 of the
3605 27 Jan 17 peter 14   License, or (at your option) any later version.
3605 27 Jan 17 peter 15
3605 27 Jan 17 peter 16   The yat library is distributed in the hope that it will be useful,
3605 27 Jan 17 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3605 27 Jan 17 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3605 27 Jan 17 peter 19   General Public License for more details.
3605 27 Jan 17 peter 20
3605 27 Jan 17 peter 21   You should have received a copy of the GNU General Public License
3605 27 Jan 17 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3605 27 Jan 17 peter 23 */
3605 27 Jan 17 peter 24
3605 27 Jan 17 peter 25 #include "BasicVector.h"
3605 27 Jan 17 peter 26 #include "BLAS_utility.h"
3605 27 Jan 17 peter 27 #include "VectorExpression.h"
3605 27 Jan 17 peter 28
3605 27 Jan 17 peter 29 #include "yat_assert.h"
3605 27 Jan 17 peter 30
3605 27 Jan 17 peter 31 #include "gsl/gsl_vector.h"
3605 27 Jan 17 peter 32
3605 27 Jan 17 peter 33 #include <cstddef>
3605 27 Jan 17 peter 34
3605 27 Jan 17 peter 35 namespace theplu {
3605 27 Jan 17 peter 36 namespace yat {
3605 27 Jan 17 peter 37 namespace utility {
3605 27 Jan 17 peter 38
3605 27 Jan 17 peter 39   // This file defines operations using both Vector (but not Matrix)
3605 27 Jan 17 peter 40
3605 27 Jan 17 peter 41   /// \cond IGNORE_DOXYGEN
3605 27 Jan 17 peter 42
3605 27 Jan 17 peter 43   namespace expression {
3605 27 Jan 17 peter 44
3605 27 Jan 17 peter 45     template<typename LHS, typename RHS, class OP>
3605 27 Jan 17 peter 46     class VectorBinary
3605 27 Jan 17 peter 47       : public VectorExpression<VectorBinary<LHS, RHS, OP> >
3605 27 Jan 17 peter 48     {
3605 27 Jan 17 peter 49     public:
3605 27 Jan 17 peter 50       VectorBinary(const BasicVector<LHS>& lhs, const BasicVector<RHS>& rhs)
3605 27 Jan 17 peter 51         : lhs_(lhs), rhs_(rhs)
3605 27 Jan 17 peter 52       {
3605 27 Jan 17 peter 53         YAT_ASSERT(lhs.size() == rhs.size());
3605 27 Jan 17 peter 54       }
3605 27 Jan 17 peter 55
3605 27 Jan 17 peter 56       size_t size(void) const { return lhs_.size(); }
3605 27 Jan 17 peter 57
3605 27 Jan 17 peter 58       double operator()(size_t i) const
3605 27 Jan 17 peter 59       { return get(i, op_); }
3605 27 Jan 17 peter 60
3605 27 Jan 17 peter 61       void calculate_gsl_vector_p(void) const
3605 27 Jan 17 peter 62       {
3605 27 Jan 17 peter 63         this->allocate_memory(size());
3605 27 Jan 17 peter 64         for (size_t i=0; i<size(); ++i)
3605 27 Jan 17 peter 65           gsl_vector_set(this->v_, i, get(i, op_));
3605 27 Jan 17 peter 66       }
3605 27 Jan 17 peter 67
3605 27 Jan 17 peter 68     private:
3605 27 Jan 17 peter 69       const BasicVector<LHS>& lhs_;
3605 27 Jan 17 peter 70       const BasicVector<RHS>& rhs_;
3605 27 Jan 17 peter 71       OP op_;
3605 27 Jan 17 peter 72
3605 27 Jan 17 peter 73       double get(size_t i, Plus) const
3605 27 Jan 17 peter 74       {
3605 27 Jan 17 peter 75         return lhs_(i) + rhs_(i);
3605 27 Jan 17 peter 76       }
3605 27 Jan 17 peter 77
3605 27 Jan 17 peter 78
3605 27 Jan 17 peter 79       double get(size_t i, Minus) const
3605 27 Jan 17 peter 80       {
3605 27 Jan 17 peter 81         return lhs_(i) - rhs_(i);
3605 27 Jan 17 peter 82       }
3605 27 Jan 17 peter 83     };
3605 27 Jan 17 peter 84
3605 27 Jan 17 peter 85
3605 27 Jan 17 peter 86     template<class T>
3605 27 Jan 17 peter 87     class ScaledVector : public VectorExpression<ScaledVector<T> >
3605 27 Jan 17 peter 88     {
3605 27 Jan 17 peter 89     public:
3605 27 Jan 17 peter 90       ScaledVector(double factor, const BasicVector<T>& vec)
3605 27 Jan 17 peter 91         : vec_(vec), factor_(factor) {}
3605 27 Jan 17 peter 92
3605 27 Jan 17 peter 93       size_t size(void) const
3605 27 Jan 17 peter 94       {
3605 27 Jan 17 peter 95         return vec_.size();
3605 27 Jan 17 peter 96       }
3605 27 Jan 17 peter 97
3605 27 Jan 17 peter 98
3605 27 Jan 17 peter 99       double operator()(size_t i) const
3605 27 Jan 17 peter 100       {
3605 27 Jan 17 peter 101         return factor_ * vec_(i);
3605 27 Jan 17 peter 102       }
3605 27 Jan 17 peter 103
3605 27 Jan 17 peter 104
3605 27 Jan 17 peter 105       void calculate_gsl_vector_p(void) const
3605 27 Jan 17 peter 106       {
3605 27 Jan 17 peter 107         YAT_ASSERT(this->v_ == NULL);
3605 27 Jan 17 peter 108         this->allocate_memory(size());
3605 27 Jan 17 peter 109         gsl_vector_memcpy(this->v_, vec_.gsl_vector_p());
3605 27 Jan 17 peter 110         gsl_vector_scale(this->v_, factor_);
3605 27 Jan 17 peter 111       }
3605 27 Jan 17 peter 112
3605 27 Jan 17 peter 113     private:
3605 27 Jan 17 peter 114       const BasicVector<T>& vec_;
3605 27 Jan 17 peter 115       double factor_;
3605 27 Jan 17 peter 116     };
3605 27 Jan 17 peter 117
3605 27 Jan 17 peter 118
3605 27 Jan 17 peter 119   } // end namespace expression
3605 27 Jan 17 peter 120
3605 27 Jan 17 peter 121   /// \endcond
3605 27 Jan 17 peter 122
3605 27 Jan 17 peter 123
3605 27 Jan 17 peter 124   /**
3678 17 Aug 17 peter 125      \relates VectorBase
3678 17 Aug 17 peter 126
3605 27 Jan 17 peter 127      Element-wise addition of vectors
3605 27 Jan 17 peter 128
3605 27 Jan 17 peter 129      \since new in yat 0.15
3605 27 Jan 17 peter 130    */
3605 27 Jan 17 peter 131   template<class Derived1, class Derived2>
3605 27 Jan 17 peter 132   expression::VectorBinary<Derived1, Derived2, expression::Plus>
3605 27 Jan 17 peter 133   operator+(const BasicVector<Derived1>& lhs, const BasicVector<Derived2>& rhs)
3605 27 Jan 17 peter 134   {
3605 27 Jan 17 peter 135     return expression::VectorBinary<Derived1, Derived2,
3605 27 Jan 17 peter 136                                     expression::Plus>(lhs, rhs);
3605 27 Jan 17 peter 137   }
3605 27 Jan 17 peter 138
3605 27 Jan 17 peter 139
3605 27 Jan 17 peter 140   /**
3624 10 Feb 17 peter 141      \relates VectorBase
3624 10 Feb 17 peter 142
3605 27 Jan 17 peter 143      Element-wise subtraction of vectors
3605 27 Jan 17 peter 144
3605 27 Jan 17 peter 145      \since new in yat 0.15
3605 27 Jan 17 peter 146    */
3605 27 Jan 17 peter 147   template<class Derived1, class Derived2>
3605 27 Jan 17 peter 148   expression::VectorBinary<Derived1, Derived2, expression::Minus>
3605 27 Jan 17 peter 149   operator-(const BasicVector<Derived1>& lhs,
3605 27 Jan 17 peter 150             const BasicVector<Derived2>& rhs)
3605 27 Jan 17 peter 151   {
3605 27 Jan 17 peter 152     return expression::VectorBinary<Derived1, Derived2,
3605 27 Jan 17 peter 153                                     expression::Minus>(lhs, rhs);
3605 27 Jan 17 peter 154   }
3605 27 Jan 17 peter 155
3605 27 Jan 17 peter 156
3605 27 Jan 17 peter 157   /**
3624 10 Feb 17 peter 158      \relates VectorBase
3624 10 Feb 17 peter 159
3605 27 Jan 17 peter 160      \since New in yat 0.15
3605 27 Jan 17 peter 161    */
3605 27 Jan 17 peter 162   template<class Derived1, class Derived2>
3605 27 Jan 17 peter 163   double operator*(const BasicVector<Derived1>& lhs,
3605 27 Jan 17 peter 164                    const BasicVector<Derived2>& rhs)
3605 27 Jan 17 peter 165   {
3605 27 Jan 17 peter 166     double res = 0.0;;
3605 27 Jan 17 peter 167     gsl_blas_ddot(lhs.gsl_vector_p(), rhs.gsl_vector_p(), &res);
3605 27 Jan 17 peter 168     return res;
3605 27 Jan 17 peter 169   }
3605 27 Jan 17 peter 170
3605 27 Jan 17 peter 171
3605 27 Jan 17 peter 172   /**
3605 27 Jan 17 peter 173      \brief scale a vector expression
3605 27 Jan 17 peter 174
3624 10 Feb 17 peter 175      \relates VectorBase
3624 10 Feb 17 peter 176
3605 27 Jan 17 peter 177      \since new in yat 0.15
3605 27 Jan 17 peter 178    */
3605 27 Jan 17 peter 179   // TODO specialization for when vec is ScaledVector
3605 27 Jan 17 peter 180   template<typename T>
3605 27 Jan 17 peter 181   expression::ScaledVector<T>
3605 27 Jan 17 peter 182   operator*(double k, const BasicVector<T>& vec)
3605 27 Jan 17 peter 183   {
3605 27 Jan 17 peter 184     return expression::ScaledVector<T>(k, vec);
3605 27 Jan 17 peter 185   }
3605 27 Jan 17 peter 186
3605 27 Jan 17 peter 187
3605 27 Jan 17 peter 188   /**
3605 27 Jan 17 peter 189      \brief scale a vector expression
3605 27 Jan 17 peter 190
3624 10 Feb 17 peter 191      \relates VectorBase
3624 10 Feb 17 peter 192
3605 27 Jan 17 peter 193      \since new in yat 0.15
3605 27 Jan 17 peter 194    */
3605 27 Jan 17 peter 195   template<typename T>
3605 27 Jan 17 peter 196   expression::ScaledVector<T>
3605 27 Jan 17 peter 197   operator*(const BasicVector<T>& vec, double k)
3605 27 Jan 17 peter 198   {
3605 27 Jan 17 peter 199     return expression::ScaledVector<T>(k, vec);
3605 27 Jan 17 peter 200   }
3605 27 Jan 17 peter 201
3610 27 Jan 17 peter 202
3610 27 Jan 17 peter 203   /**
3610 27 Jan 17 peter 204      \brief negation operator
3610 27 Jan 17 peter 205
3624 10 Feb 17 peter 206      \relates VectorBase
3624 10 Feb 17 peter 207
3610 27 Jan 17 peter 208      \since new in yat 0.15
3610 27 Jan 17 peter 209    */
3610 27 Jan 17 peter 210   template<typename T>
3610 27 Jan 17 peter 211   expression::ScaledVector<T>
3610 27 Jan 17 peter 212   operator-(const BasicVector<T>& vec)
3610 27 Jan 17 peter 213   {
3610 27 Jan 17 peter 214     return expression::ScaledVector<T>(-1.0, vec);
3610 27 Jan 17 peter 215   }
3610 27 Jan 17 peter 216
3605 27 Jan 17 peter 217 }}} // of namespace utility, yat, and theplu
3605 27 Jan 17 peter 218
3605 27 Jan 17 peter 219 #endif

Generated on Sat Oct 14 02:24:22 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)