| # Gaudi Framework, version v25r0 |

Home |  | Generated: Mon Feb 17 2014

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
EqSolver.cpp
Go to the documentation of this file.
1 // $Id: EqSolver.cpp,v 1.4 2006/01/10 20:00:05 hmd Exp $
2 // ============================================================================
3 // Include files
4
5 #include <stdlib.h>
6 #include <stdio.h>
7 // from Gaudi
8 #include "GaudiKernel/ToolFactory.h"
9 #include "GaudiKernel/MsgStream.h"
10 //from GSL
11 #include "gsl/gsl_vector.h"
12 #include "gsl/gsl_multimin.h"
13 #include "gsl/gsl_math.h"
14 #include "gsl/gsl_matrix.h"
15 #include "gsl/gsl_linalg.h"
16 #include "gsl/gsl_blas.h"
17 #include "gsl/gsl_errno.h"
18
19 // local
20 #include "EqSolver.h"
21
29 // Declaration of the Tool Factory
30 DECLARE_COMPONENT(EqSolver)
31
32
33 // ============================================================================
34 // ============================================================================
35 EqSolver::EqSolverMisc::EqSolverMisc
36 ( const EqSolver::Equations& funcs ,
37 EqSolver::Arg& arg )
38 : m_argum ( arg )
39 , m_eqs ( &funcs )
40 , m_jac ()
41 {
42 const size_t N = funcs.size () ;
43 for( size_t i = 0 ; i < N ; ++i )
44 {
45 Equations last;
46 for( size_t j = 0 ; j < N ; ++j )
47 {
48 Genfun::GENFUNCTION fij = funcs[i]->partial(j);
49 last.push_back( fij.clone() ) ;
50 }
51 m_jac.push_back( last );
52 }
53 }
54 // ============================================================================
55
56 // ============================================================================
57 EqSolver::EqSolverMisc::~EqSolverMisc()
58 {
59 while( !m_jac.empty() )
60 {
61 Equations& last = m_jac.back() ;
62 while( !last.empty() )
63 {
64 delete last.back() ;
65 last.pop_back () ;
66 }
67 m_jac.pop_back();
68 }
69 }
70 // ============================================================================
71
72 //=============================================================================
73 // Standard constructor, initializes variables
74 //=============================================================================
75 EqSolver::EqSolver( const std::string& type,
76 const std::string& name,
77 const IInterface* parent )
78 : base_class ( type, name , parent )
79 , m_algType ( "fdfsolver_hybridsj" )
80 , m_max_iter ( 1000 )
81 , m_norm_residual ( 1.0e-7 )
82 , m_type ( 0 )
83
84 {
86 declareProperty ( "Algorithm", m_algType );
88 declareProperty ( "Iteration", m_max_iter );
90 declareProperty ( "Residual" , m_norm_residual );
91 }
92
93 namespace
94 {
95 using namespace Genfun;
96
99 int fun_gsl ( const gsl_vector* v ,
100 void* params ,
101 gsl_vector* f )
102 {
103 EqSolver::EqSolverMisc* local =
104 static_cast <EqSolver::EqSolverMisc*> (params);
105 const EqSolver::Equations& eqs = *(local->equations()) ;
106 Argument& arg = local->argument() ;
107
108 {for ( unsigned int i = 0; i < v->size; ++i)
109 {
110 arg[i] = gsl_vector_get (v, i);
111 }
112 }
113
114 {for ( unsigned int i = 0; i < v->size; ++i)
115 {
116 gsl_vector_set(f, i, (*eqs[i])(arg));
117 }
118 }
119 return GSL_SUCCESS;
120 }
121
122
125 int dfun_gsl ( const gsl_vector* v ,
126 void* params ,
127 gsl_matrix* df )
128 {
129 EqSolver::EqSolverMisc* local =
130 static_cast <EqSolver::EqSolverMisc*> (params);
131 const EqSolver::Jacobi& jac = local -> jacobi() ;
132 Argument& arg = local -> argument() ;
133
134 {for ( unsigned int i = 0; i < v->size; ++i)
135 {
136 arg[i] = gsl_vector_get (v, i);
137 }
138 }
139
140 {for( unsigned int i = 0 ; i < v->size ; ++i )
141 {
142 for (unsigned int j = 0; j < v->size; ++j)
143 {
144 Genfun::GENFUNCTION f = *(jac[i][j]) ;
145 gsl_matrix_set ( df , i , j , f ( arg ) ) ;
146 }
147 }
148 }
149 return GSL_SUCCESS;
150 }
151
155 int fdfun_gsl ( const gsl_vector* v ,
156 void* params ,
157 gsl_vector* f ,
158 gsl_matrix* df )
159 {
160 EqSolver::EqSolverMisc* local =
161 static_cast <EqSolver::EqSolverMisc*> (params);
162 const EqSolver::Equations& eqs = *(local->equations()) ;
163 const EqSolver::Jacobi& jac = local->jacobi() ;
164 Argument& arg = local->argument() ;
165
166 {for ( unsigned int i = 0; i < v->size; ++i)
167 {
168 arg[i] = gsl_vector_get (v, i);
169 }
170 }
171
172 {for( unsigned int i = 0 ; i < v->size ; ++i )
173 {
174 gsl_vector_set(f, i, (*eqs[i])(arg ) ) ;
175 for (unsigned int j = 0; j < v->size; ++j)
176 {
177 Genfun::GENFUNCTION f1 = *(jac[i][j]) ;
178 gsl_matrix_set ( df , i , j , f1(arg) ) ;
179 }
180 }
181 }
182 return GSL_SUCCESS;
183 }
184 }
185
186 //=============================================================================
191 StatusCode EqSolver::solver (const Equations& funcs ,
192 Arg& arg ) const
193 //=============================================================================
194 {
195 using namespace Genfun;
196
197 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
198 arg.dimension() );
199 MsgStream log( msgSvc(), name() );
200
201 EqSolverMisc local (funcs, arg);
202
203 const gsl_multiroot_fdfsolver_type *T = m_type;
204 gsl_multiroot_fdfsolver *s;
205 int status = 0 ;
206 size_t iter = 0 ;
207
208 gsl_multiroot_function_fdf function;
209
210 function.f = &fun_gsl;
211 function.df = &dfun_gsl;
212 function.fdf = &fdfun_gsl;
213 function.n = vect.vector.size;
214 function.params = (void*) &local;
215
216 s = gsl_multiroot_fdfsolver_alloc(T, vect.vector.size);
217 gsl_multiroot_fdfsolver_set (s, &function, &vect.vector);
218
219 for (iter = 0; iter < m_max_iter; ++iter)
220 {
221 status = gsl_multiroot_fdfsolver_iterate (s);
222 if (status)
223 {
224 return Error
225 ("Error from gsl_gsl_multiroot_fdfsolver_iterate '"
226 + std::string(gsl_strerror(status)) + "'") ;
227 }
228
229 status = gsl_multiroot_test_residual (s->f,
230 m_norm_residual);
231
232 if ( status != GSL_CONTINUE ) { break; }
233 }
234
235 for (unsigned int i = 0; i < vect.vector.size; ++i)
236 {
237 gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
238 }
239
240 if (status == GSL_SUCCESS)
241 {
242 log << MSG::DEBUG
243 << "We stopped in the method on the " << iter
244 << " iteration (we have maximum " << m_max_iter
245 << " iterations)" << endmsg;
246 }
247 else if (status == GSL_CONTINUE && iter <= m_max_iter )
248 {
249 return Error ( "Method finished with '"
250 + std::string(gsl_strerror(status))
251 + "' error" );
252 }
253 else
254 {
255 return Error ( "Method finished with '" +
256 std::string(gsl_strerror(status))
257 + "' error" );
258 }
259
260 gsl_multiroot_fdfsolver_free (s);
261
262 if (status)
263 {
264 return Error ( "Method finished with '"
265 + std::string(gsl_strerror(status))
266 + "' error" );
267 }
268
269 return GSL_SUCCESS;
270 }
271
272 //=============================================================================
273 StatusCode EqSolver::initialize()
274
275 {
276 StatusCode sc = GaudiTool::initialize() ;
277
278 MsgStream log( msgSvc() , name() ) ;
279
280 if( sc.isFailure() )
281 {
282 return Error ("Could not initiliaze base class GaudiTool", sc);
283 }
284 /* The algorithm for multiional root-finding
285 (solving nonlinear systems with n equations in n unknowns)*/
286
287 if( "fdfsolver_hybridsj" == m_algType )
288 {
289 m_type = gsl_multiroot_fdfsolver_hybridsj ;
290 log << MSG::DEBUG
291 << "Root findind algoritms to be used: "
292 << "'gsl_multiroot_fdfsolver_hybridsj'"
293 << endmsg;
294 }
295 else if ( "fdfsolver_hybridj" == m_algType )
296 {
297 m_type = gsl_multiroot_fdfsolver_hybridj ;
298 log << MSG::DEBUG
299 << "Root findind algoritms to be used: "
300 << "'gsl_multiroot_fdfsolver_hybridj'"
301 << endmsg;
302 }
303 else if ( "fdfsolver_newton" == m_algType )
304 {
305 m_type = gsl_multiroot_fdfsolver_newton ;
306 log << MSG::DEBUG
307 << "Root findind algoritms to be used: "
308 << "'gsl_multiroot_fdfsolver_newton'"
309 << endmsg;
310 }
311 else if ( "fdfsolver_gnewton" == m_algType )
312 {
313 m_type = gsl_multiroot_fdfsolver_gnewton ;
314 log << MSG::DEBUG
315 << "Root findind algoritms to be used: "
316 << "'gsl_multiroot_fdfsolver_gnewton'"
317 << endmsg;
318 }
319 else
320 {
321 return Error(" Unknown algorith type '"
322 + std::string(m_algType) + "'");
323 }
324
325 return StatusCode::SUCCESS;
326 }
327 //=============================================================================
328 StatusCode EqSolver::finalize ()
329 {
330 StatusCode sc = GaudiTool::finalize() ;
331
332 MsgStream log( msgSvc() , name() ) ;
333
334 if ( sc.isFailure() )
335 {
336 return Error("Could not finaliaze base class GaudiTool", sc);
337 }
338 return StatusCode::SUCCESS;
339 }
340 //=============================================================================
341 EqSolver::~EqSolver( )
342 {}
343 //=============================================================================