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
dirscalingmatrix.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 __DirScalingMatrix_H
30 #define __DirScalingMatrix_H 1
31 
32 #include <algorithm>
33 
34 
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/banded.hpp>
37 #include <boost/numeric/ublas/io.hpp>
38 
39 namespace Feel
40 {
41 using namespace boost::numeric::ublas;
42 
50 template<typename NumType>
52 {
53 public:
54 
55 
59 
60  typedef enum mode_type { NO_JACOBIAN, WITH_JACOBIAN };
61 
62  typedef NumType value_type;
63 
64  typedef vector<value_type> vector_type;
65  typedef banded_matrix<value_type> matrix_type;
66 
68 
72 
74  :
75  M_lb(),
76  M_ub(),
77  M_lb_ub(),
78  M_value(),
79  M_jacobian(),
80  M_trust_region_active( false )
81  {}
82 
83  DirScalingMatrix( vector_type const& __lb, vector_type const& __ub )
84  :
85  M_lb( __lb ),
86  M_ub( __ub ),
87  M_lb_ub( M_ub - M_lb ),
88  M_value( __lb.size(), __lb.size(), 0, 0 ),
89  M_jacobian( __lb.size(), __lb.size(), 0, 0 ),
90  M_trust_region_active( false )
91  {
92  GST_SMART_ASSERT( __lb.size() == __ub.size() )( __lb )( __ub )( "inconsistent bounds definition" );
93  GST_SMART_ASSERT( *std::min_element( M_lb_ub.begin(), M_lb_ub.end() ) >= 0 )
94  ( M_lb )( M_ub )( "lower and upper bounds are not properly defined" );
95  }
97  {}
99  {}
100 
102 
106 
107 
109 
113 
114  value_type zeta( vector_type const& __x ) const;
115  value_type zeta() const
116  {
117  return M_zeta;
118  }
119 
120  matrix_type const& operator()() const
121  {
122  return M_value;
123  }
124 
125  matrix_type const& jacobian() const
126  {
127  return M_jacobian;
128  }
129 
130  bool isTrustRegionActive() const
131  {
132  return M_trust_region_active;
133  }
135 
139 
140  void update( value_type const&, vector_type const&, vector_type const&, mode_type = WITH_JACOBIAN );
141 
142  void setBounds( vector_type const& __lb, vector_type const& __up )
143  {
144  GST_SMART_ASSERT( __lb.size() == __up.size() )( __lb )( __up )( "inconsistent bounds definition" );
145  M_lb = __lb;
146  M_ub = __up;
147  M_lb_ub = __up - __lb;
148 
149  GST_SMART_ASSERT( *std::min_element( M_lb_ub.begin(), M_lb_ub.end() ) >= 0 )
150  ( M_lb )( M_ub )( "lower and upper bounds are not properly defined" );
151  }
153 
157 
159 
160 
161 
162 protected:
163 
164  vector_type distanceToLB( vector_type const& __x ) const
165  {
166  GST_SMART_ASSERT( __x.size() == M_lb.size() )( __x )( M_lb )( "inconsistent bounds definition" );
167  return element_div( __x - M_lb, M_lb_ub );
168  }
169 
170  vector_type distanceToUB( vector_type const& __x ) const
171  {
172  GST_SMART_ASSERT( __x.size() == M_ub.size() )( __x )( M_ub )( "inconsistent bounds definition" );
173  return element_div( __x - M_ub, M_lb_ub );
174  }
175 
176 private:
177 
178  vector_type M_lb;
179  vector_type M_ub;
180  vector_type M_lb_ub;
181 
182  matrix_type M_value;
183  matrix_type M_jacobian;
184 
185  mutable value_type M_zeta;
186 
187  bool M_trust_region_active;
188 };
189 
190 template<typename NumType>
191 typename DirScalingMatrix<NumType>::value_type
193 {
194  M_zeta = 0.9;//M_options.zeta_min;
195 
196  M_zeta = std::max( M_zeta, std::max( ublas::norm_inf( distanceToLB( __x ) ),
197  ublas::norm_inf( distanceToUB( __x ) ) ) );
198  return M_zeta;
199 }
200 template<typename NumType>
201 void
202 DirScalingMatrix<NumType>::update( value_type const& __Delta,
203  vector_type const& __x,
204  vector_type const& __s,
205  mode_type __mode )
206 {
207  GST_SMART_ASSERT( __x.size() == M_lb.size() )( __x )( M_lb )( "inconsistent bounds definition" );
208  GST_SMART_ASSERT( __x.size() == M_ub.size() )( __x )( M_ub )( "inconsistent bounds definition" );
209 
210  M_value.resize( __x.size(), __x.size(), 0, 0 );
211  M_jacobian.resize( __x.size(), __x.size(), 0, 0 );
212 
213 
214  M_zeta = zeta( __x );
215 
216  vector_type __dl = distanceToLB( __x );
217  vector_type __du = distanceToUB( __x );
218 
219  if ( M_zeta * std::min( norm_inf( __dl ), norm_inf( __du ) ) > __Delta )
220  {
221  // we are in the trust region
222  M_value = identity_matrix<value_type>( M_value.size1(), M_value.size2() );
223  }
224 
225  else
226  {
227  M_trust_region_active = true;
228 
229  for ( size_t __i = 0; __i < __x.size(); ++__i )
230  {
231  if ( __s ( __i ) < 0 )
232  M_value ( __i, __i ) = M_zeta * std::min( 1. , __dl( __i ) ) / __Delta;
233 
234  else
235  M_value ( __i, __i ) = M_zeta * std::min( 1. , __du( __i ) ) / __Delta;
236  }
237  }
238 
239  if ( __mode == WITH_JACOBIAN )
240  {
241  for ( size_t i = 0; i < __x.size(); i++ )
242  {
243  if ( ( __s( i ) < 0 ) && ( __x( i ) < M_lb( i ) + __Delta ) )
244  {
245  M_jacobian ( i, i ) = M_zeta / __Delta;
246  }
247 
248  else if ( ( __s ( i ) > 0 ) && ( __x ( i ) > M_ub( i ) - __Delta ) )
249  {
250  M_jacobian ( i, i ) = -M_zeta / __Delta;
251  }
252 
253  else
254  {
255  M_jacobian ( i, i ) = 0;
256  }
257  }
258  }
259 }
260 
261 }
262 #endif /* __DirScalingMatrix_H */
263 

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