Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
adtypeorder1.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2008-02-14
7 
8  Copyright (C) 2008 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __ADTypeOrder1_H
30 #define __ADTypeOrder1_H 1
31 
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 
35 #include <boost/numeric/ublas/io.hpp>
36 
37 namespace Feel
38 {
55 template<typename T, int Nvar, int Var>
56 class ADType<T, Nvar, 1, Var>
57 {
58 public:
59 
60  enum { nvar = Nvar };
61  enum { order = 1 };
62 
63  typedef ADVariable<Var> variable_type;
64 
65  typedef T value_type;
66 
67  typedef boost::numeric::ublas::vector<value_type> gradient_type;
68  //typedef typename STinyVector<value_type, Nvar>::type gradient_type;
69  //typedef typename STinyVector<bool, Nvar>::type deps_type;
70  typedef std::vector<int> deps_type;
71  typedef ADType<T,Nvar, 1, Var> This;
72  typedef std::set<int> dependency_type;
73 
74  template<typename NumT, int NumVar, int Order, int VarNum> friend class ADType;
75 
78  //@ {
79  ADType()
80  :
81  M_val( 0. ),
82  M_grad( nvar ),
83  __dep( nvar, false ),
84  __deps()
85  {}
86 
87  ADType( value_type val )
88  :
89  M_val( val ),
90  M_grad( nvar ),
91  __dep( nvar, false ),
92  __deps()
93  {
94  M_grad = boost::numeric::ublas::zero_vector<value_type>( nvar );
95 
96  if ( Var != -1 )
97  {
98  __dep[Var]=true;
99  M_grad( Var ) = 1.;
100  }
101  }
102 
103  template<int VarNum>
104  ADType( ADType<T,Nvar,1,VarNum> const& sad )
105  :
106  M_val( sad.M_val ),
107  M_grad( sad.M_grad ),
108  __dep( sad.__dep ),
109  __deps()
110  {
111  }
112 
113  template<typename ExprT>
114  ADType ( const ADExpr<ExprT>& expr )
115  :
116  M_val( 0 ),
117  M_grad( nvar ),
118  __dep( nvar, false ),
119  __deps()
120  {
121  *this = expr;
122  }
124 
128  value_type value() const
129  {
130  return M_val;
131  }
132  deps_type deps() const
133  {
134  return __dep;
135  }
136 
138  bool deps( int __i ) const
139  {
140  ENSURE( __i >= 0 && __i < nvar );
141  return __dep[ __i ];
142  }
143 
144  value_type grad( int __i ) const
145  {
146  return M_grad( __i );
147  }
148 
149  gradient_type grad() const
150  {
151  return M_grad;
152  }
153 
155 
159  value_type& value()
160  {
161  return M_val;
162  }
164  bool& deps( int __i )
165  {
166  ENSURE( __i >= 0 && __i < nvar );
167  return __dep[ __i ];
168  }
169  value_type& grad( int __i )
170  {
171  return M_grad( __i );
172  }
173 
174  gradient_type& grad()
175  {
176  return M_grad;
177  }
178 
180 
184 
185  This& operator=( T const& );
186  This& operator=( This const& );
187  template <class ExprT> This& operator=( const ADExpr<ExprT>& expr );
188 
190  ADExpr< ADUnaryPlus< This > > operator+ () const
191  {
192  typedef ADUnaryPlus<This> expr_t;
193  return ADExpr<expr_t> ( expr_t ( *this ) );
194  }
195 
197  ADExpr< ADUnaryMinus< This > > operator- () const
198  {
199  typedef ADUnaryMinus<This> expr_t;
200  return ADExpr<expr_t> ( expr_t ( *this ) );
201  }
202 
203 #define AD_UNARY_OP( op ) \
204  This& operator op ( value_type val ) \
205  { \
206  M_val op val; \
207  return *this; \
208  }
209 
210  AD_UNARY_OP( += );
211  AD_UNARY_OP( -= );
212  AD_UNARY_OP( *= );
213  AD_UNARY_OP( /= );
214 
215 #undef AD_UNARY_OP
216 
217  This& operator += ( This const& sad )
218  {
219  M_val += sad.M_val;
220  M_grad += sad.M_grad;
221  return *this;
222  }
223  This& operator -= ( This const& sad )
224  {
225  M_val -= sad.M_val;
226  M_grad -= sad.M_grad;
227  return *this;
228  }
229 
230 
231  This& operator *= ( This const& sad )
232  {
233  M_val *= sad.M_val;
234  M_grad = M_grad * sad.M_val + M_val * sad.M_grad;
235  return *this;
236  }
237 
238  This& operator /= ( This const& sad )
239  {
240  M_val /= sad.M_val;
241  M_grad = ( M_grad * sad.M_val - M_val * sad.M_grad ) / ( sad.M_val * sad.M_val );
242  return *this;
243  }
244 
245  template<typename Expr>
246  This& operator += ( ADExpr<Expr> const& sad )
247  {
248  *this = *this + sad;
249  return *this;
250  }
251 
252  template<typename Expr>
253  This& operator -= ( ADExpr<Expr> const& sad )
254  {
255  *this = *this - sad;
256  return *this;
257  }
258 
259  template<typename Expr>
260  This& operator *= ( ADExpr<Expr> const& sad )
261  {
262  *this = *this * sad;
263  return *this;
264  }
265 
266  template<typename Expr>
267  This& operator /= ( ADExpr<Expr> const& sad )
268  {
269  *this = *this / sad;
270  return *this;
271  }
273 protected:
274 
275 private:
276  value_type M_val;
277  gradient_type M_grad;
278 
279  deps_type __dep;
280  dependency_type __deps;
281 };
282 
283 template<typename T,int Nvar, int Var>
284 ADType<T, Nvar, 1, Var>&
285 ADType<T, Nvar, 1, Var>::operator=( value_type const& val )
286 {
287  M_val = val;
288  M_grad = boost::numeric::ublas::zero_vector<value_type>( nvar );
289  __dep = false;
290  return *this;
291 }
292 
293 template<typename T,int Nvar, int Var>
294 ADType<T, Nvar, 1, Var>&
295 ADType<T, Nvar, 1, Var>::operator=( This const& sad )
296 {
297  M_val = sad.M_val;
298  M_grad = sad.M_grad;
299  __dep = sad.__dep;
300  return *this;
301 }
302 template<typename T,int Nvar, int Var>
303 template <class ExprT>
304 ADType<T,Nvar, 1, Var> &
305 ADType<T,Nvar, 1, Var>::operator=( const ADExpr<ExprT>& expr )
306 {
307  M_val = expr.value();
308 #if 0
309  __deps.clear();
310 
311  for ( int __i = 0; __i < nvar; ++__i )
312  {
313 
314  __dep( __i ) = expr.deps( __i );
315 
316  if ( __dep( __i ) )
317  {
318  //std::cerr << "depends on " << __i << "\n";
319  __deps.insert( __i );
320  }
321 
322  }
323 
324  if ( ! __deps.empty() && __deps.size() < nvar )
325  {
326  //M_grad = 0;
327  //__hessian = M_grad;
328 
329  dependency_type::iterator __begin = __deps.begin();
330  dependency_type::iterator __end = __deps.end();
331 
332  dependency_type::iterator __it = __begin;
333 
334  while ( __it != __end )
335  {
336  M_grad( *__it )= expr.grad( *__it );
337  ++__it;
338  }
339  }
340 
341  else
342  {
343 #endif
344 
345  for ( int __i = 0; __i < nvar; ++__i )
346  {
347  M_grad( __i )= expr.grad( __i );
348  }
349 
350 #if 0
351  }
352 
353 #endif
354  return *this;
355 }
356 
357 
358 
359 }
360 
361 
362 //------------------------------- AD ostream operator ------------------------------------------
363 template <class T, int Nvar, int Var>
364 std::ostream&
365 operator << ( std::ostream& os, const Feel::ADType<T, Nvar, 1, Var>& a )
366 {
367  os.setf( std::ios::fixed,std::ios::floatfield );
368  os.width( 12 );
369  os << "value = " << a.value() << " \n";
370  os << "gradient = " << a.grad() << "\n";
371  return os;
372 }
373 
374 
375 #endif /* __ADTypeOrder1_H */

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4