yat/utility/VectorExpression.h

Code
Comments
Other
Rev Date Author Line
3605 27 Jan 17 peter 1 #ifndef _theplu_yat_utility_vector_expression
3605 27 Jan 17 peter 2 #define _theplu_yat_utility_vector_expression
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 /*
4359 23 Aug 23 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 "config_public.h"
3605 27 Jan 17 peter 26
3605 27 Jan 17 peter 27 #include "BasicVector.h"
3605 27 Jan 17 peter 28 #include "MatrixExpression.h"
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_blas.h>
3605 27 Jan 17 peter 32 #include <gsl/gsl_matrix.h>
3605 27 Jan 17 peter 33 #include <gsl/gsl_vector.h>
3605 27 Jan 17 peter 34
3605 27 Jan 17 peter 35 #include <utility>
3605 27 Jan 17 peter 36
3605 27 Jan 17 peter 37 namespace theplu {
3605 27 Jan 17 peter 38 namespace yat {
3605 27 Jan 17 peter 39 namespace utility {
3605 27 Jan 17 peter 40
3605 27 Jan 17 peter 41   /**
3605 27 Jan 17 peter 42      \brief An expression that can be converted to a Vector.
3605 27 Jan 17 peter 43
3605 27 Jan 17 peter 44      yat uses template expression to implement fast linear algebra
3605 27 Jan 17 peter 45      operations. This class defines the interface of expression that
3605 27 Jan 17 peter 46      can be converted to a vector.
3605 27 Jan 17 peter 47
3605 27 Jan 17 peter 48      The class is used in two ways: 1) Classes that implements an
3605 27 Jan 17 peter 49      expression that translates to a vector, for example, A * x,
3605 27 Jan 17 peter 50      should inherit from this class. 2) Functions that takes
3605 27 Jan 17 peter 51      expression that looks like a vector, but do not accept
3605 27 Jan 17 peter 52      VectorBase, can pass this class.
3605 27 Jan 17 peter 53
3605 27 Jan 17 peter 54      To inherit from this class do the following:
3605 27 Jan 17 peter 55
3605 27 Jan 17 peter 56      1) declare you class
3605 27 Jan 17 peter 57      class MyClass : public VectorExpression<MyClass>
3605 27 Jan 17 peter 58
3605 27 Jan 17 peter 59      Implement functions
3605 27 Jan 17 peter 60      - size_t size(void) const
3605 27 Jan 17 peter 61      - double operator()(size_t) const
3605 27 Jan 17 peter 62      - void calculate_gsl_vector_p(void) const
3605 27 Jan 17 peter 63
3605 27 Jan 17 peter 64      \since new in yat 0.15
3605 27 Jan 17 peter 65   */
3605 27 Jan 17 peter 66   template<class Derived>
3605 27 Jan 17 peter 67   class VectorExpression
3605 27 Jan 17 peter 68     : public BasicVector<VectorExpression<Derived> >
3605 27 Jan 17 peter 69   {
3605 27 Jan 17 peter 70   public:
3605 27 Jan 17 peter 71     /// Derived class
3605 27 Jan 17 peter 72     typedef Derived derived_type;
3605 27 Jan 17 peter 73
3605 27 Jan 17 peter 74     /// default constructor
3605 27 Jan 17 peter 75     VectorExpression(void) : v_(NULL) {}
3605 27 Jan 17 peter 76
3605 27 Jan 17 peter 77     /// destructor, free allocated memory
3605 27 Jan 17 peter 78     ~VectorExpression(void)
3605 27 Jan 17 peter 79     {
3605 27 Jan 17 peter 80       delete_allocated_memory();
3605 27 Jan 17 peter 81     }
3605 27 Jan 17 peter 82
3605 27 Jan 17 peter 83     /// Cop constructor
3605 27 Jan 17 peter 84     VectorExpression(const VectorExpression& other)
3605 27 Jan 17 peter 85       : v_(NULL)
3605 27 Jan 17 peter 86     {
3605 27 Jan 17 peter 87       if (other.v_) {
3605 27 Jan 17 peter 88         allocate_memory(other.size());
3605 27 Jan 17 peter 89         if (gsl_vector_memcpy(v_, other.v_))
3605 27 Jan 17 peter 90           throw GSL_error("VectorExpression: gsl_vector_memcpy failed");
3605 27 Jan 17 peter 91       }
3605 27 Jan 17 peter 92     }
3605 27 Jan 17 peter 93
3605 27 Jan 17 peter 94     /// Move constructor
3605 27 Jan 17 peter 95     VectorExpression(VectorExpression&& other)
3605 27 Jan 17 peter 96       : v_(NULL)
3605 27 Jan 17 peter 97     {
3605 27 Jan 17 peter 98       std::swap(v_, other.v_);
3605 27 Jan 17 peter 99     }
3605 27 Jan 17 peter 100
3605 27 Jan 17 peter 101     /// \return number of elements in vector
3605 27 Jan 17 peter 102     size_t size(void) const
3605 27 Jan 17 peter 103     {  return static_cast<const Derived*>(this)->size();  }
3605 27 Jan 17 peter 104
3605 27 Jan 17 peter 105     /// access an element
3605 27 Jan 17 peter 106     double operator()(size_t i) const
3605 27 Jan 17 peter 107     {  return (*static_cast<const Derived*>(this))(i);  }
3605 27 Jan 17 peter 108
3605 27 Jan 17 peter 109     /// return a GSL vector pointer
3605 27 Jan 17 peter 110     const gsl_vector* gsl_vector_p(void) const
3605 27 Jan 17 peter 111     {
3605 27 Jan 17 peter 112       if (!v_)
3605 27 Jan 17 peter 113         static_cast<const Derived*>(this)->calculate_gsl_vector_p();
3605 27 Jan 17 peter 114       return v_;
3605 27 Jan 17 peter 115     }
3605 27 Jan 17 peter 116
3605 27 Jan 17 peter 117     /// return a mutable GSL vector pointer
3605 27 Jan 17 peter 118     gsl_vector* gsl_vector_p(void)
3605 27 Jan 17 peter 119     {
3605 27 Jan 17 peter 120       if (!v_)
3605 27 Jan 17 peter 121         static_cast<const Derived*>(this)->calculate_gsl_vector_p();
3605 27 Jan 17 peter 122       return v_;
3605 27 Jan 17 peter 123     }
3605 27 Jan 17 peter 124
3605 27 Jan 17 peter 125
3605 27 Jan 17 peter 126     /// set underlying GSL vector
3605 27 Jan 17 peter 127     void gsl_vector_p(gsl_vector* v)
3605 27 Jan 17 peter 128     {
3605 27 Jan 17 peter 129       v_ = v;
3605 27 Jan 17 peter 130     }
3605 27 Jan 17 peter 131
3605 27 Jan 17 peter 132   protected:
3605 27 Jan 17 peter 133     /// cache variable that is calculated in derived classes via
3605 27 Jan 17 peter 134     /// calculate_gsl_vector_p(void);
3605 27 Jan 17 peter 135     mutable gsl_vector* v_; // lazy eval
3605 27 Jan 17 peter 136
3605 27 Jan 17 peter 137     /// \brief allocate memory and throw if fails
3605 27 Jan 17 peter 138     void allocate_memory(size_t n) const
3605 27 Jan 17 peter 139     {
3605 27 Jan 17 peter 140       v_ = gsl_vector_alloc(n);
3605 27 Jan 17 peter 141       if (!v_)
3605 27 Jan 17 peter 142         throw GSL_error("Vector: failed to allocate memory");
3605 27 Jan 17 peter 143     }
3605 27 Jan 17 peter 144
3605 27 Jan 17 peter 145     /// \brief clean up
3605 27 Jan 17 peter 146     void delete_allocated_memory(void)
3605 27 Jan 17 peter 147     {
3907 09 May 20 peter 148       gsl_vector_free(v_);
3605 27 Jan 17 peter 149       v_ = NULL;
3605 27 Jan 17 peter 150     }
3605 27 Jan 17 peter 151
3606 27 Jan 17 peter 152   private:
3606 27 Jan 17 peter 153     // assignment not allowed
3606 27 Jan 17 peter 154     VectorExpression& operator=(const VectorExpression& rhs);
3605 27 Jan 17 peter 155   };
3605 27 Jan 17 peter 156
3605 27 Jan 17 peter 157
3605 27 Jan 17 peter 158 }}} // of namespace utility, yat, and theplu
3605 27 Jan 17 peter 159
3605 27 Jan 17 peter 160 #endif