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

Generated on Wed May 15 2019 03:27:28 for yat by  doxygen 1.8.11