| The Gaudi Framework v27r1

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

- File List
- File Members

EqSolver.cpp
Go to the documentation of this file.
1 // Include files
2
3 #include <stdlib.h>
4 #include <stdio.h>
5 // from Gaudi
6 #include "GaudiKernel/MsgStream.h"
7 //from GSL
8 #include "gsl/gsl_vector.h"
9 #include "gsl/gsl_multimin.h"
10 #include "gsl/gsl_math.h"
11 #include "gsl/gsl_matrix.h"
12 #include "gsl/gsl_linalg.h"
13 #include "gsl/gsl_blas.h"
14 #include "gsl/gsl_errno.h"
15
16 // local
17 #include "EqSolver.h"
18
26 // Declaration of the Tool Factory
27 DECLARE_COMPONENT(EqSolver)
28
29
30 // ============================================================================
31 // ============================================================================
32 EqSolver::EqSolverMisc::EqSolverMisc
33 ( const EqSolver::Equations& funcs ,
34 EqSolver::Arg& arg )
35 : m_argum ( arg )
36 , m_eqs ( &funcs )
37 , m_jac ()
38 {
39 const size_t N = funcs.size () ;
40 for( size_t i = 0 ; i < N ; ++i )
41 {
42 Equations last;
43 for( size_t j = 0 ; j < N ; ++j )
44 {
45 Genfun::GENFUNCTION fij = funcs[i]->partial(j);
46 last.push_back( fij.clone() ) ;
47 }
48 m_jac.push_back( last );
49 }
50 }
51 // ============================================================================
52
53 // ============================================================================
54 EqSolver::EqSolverMisc::~EqSolverMisc()
55 {
56 while( !m_jac.empty() )
57 {
58 Equations& last = m_jac.back() ;
59 while( !last.empty() )
60 {
61 delete last.back() ;
62 last.pop_back () ;
63 }
64 m_jac.pop_back();
65 }
66 }
67 // ============================================================================
68
69 //=============================================================================
70 // Standard constructor, initializes variables
71 //=============================================================================
72 EqSolver::EqSolver( const std::string& type,
73 const std::string& name,
74 const IInterface* parent )
75 : base_class ( type, name , parent )
76 , m_algType ( "fdfsolver_hybridsj" )
77 , m_max_iter ( 1000 )
78 , m_norm_residual ( 1.0e-7 )
79 , m_type ( nullptr )
80
81 {
83 declareProperty ( "Algorithm", m_algType );
85 declareProperty ( "Iteration", m_max_iter );
87 declareProperty ( "Residual" , m_norm_residual );
88 }
89
90 namespace
91 {
92 using namespace Genfun;
93
96 int fun_gsl ( const gsl_vector* v ,
97 void* params ,
98 gsl_vector* f )
99 {
100 EqSolver::EqSolverMisc* local =
101 static_cast <EqSolver::EqSolverMisc*> (params);
102 const EqSolver::Equations& eqs = *(local->equations()) ;
103 Argument& arg = local->argument() ;
104
105 {for ( unsigned int i = 0; i < v->size; ++i)
106 {
107 arg[i] = gsl_vector_get (v, i);
108 }
109 }
110
111 {for ( unsigned int i = 0; i < v->size; ++i)
112 {
113 gsl_vector_set(f, i, (*eqs[i])(arg));
114 }
115 }
116 return GSL_SUCCESS;
117 }
118
119
122 int dfun_gsl ( const gsl_vector* v ,
123 void* params ,
124 gsl_matrix* df )
125 {
126 EqSolver::EqSolverMisc* local =
127 static_cast <EqSolver::EqSolverMisc*> (params);
128 const EqSolver::Jacobi& jac = local -> jacobi() ;
129 Argument& arg = local -> argument() ;
130
131 {for ( unsigned int i = 0; i < v->size; ++i)
132 {
133 arg[i] = gsl_vector_get (v, i);
134 }
135 }
136
137 {for( unsigned int i = 0 ; i < v->size ; ++i )
138 {
139 for (unsigned int j = 0; j < v->size; ++j)
140 {
141 Genfun::GENFUNCTION f = *(jac[i][j]) ;
142 gsl_matrix_set ( df , i , j , f ( arg ) ) ;
143 }
144 }
145 }
146 return GSL_SUCCESS;
147 }
148
152 int fdfun_gsl ( const gsl_vector* v ,
153 void* params ,
154 gsl_vector* f ,
155 gsl_matrix* df )
156 {
157 EqSolver::EqSolverMisc* local =
158 static_cast <EqSolver::EqSolverMisc*> (params);
159 const EqSolver::Equations& eqs = *(local->equations()) ;
160 const EqSolver::Jacobi& jac = local->jacobi() ;
161 Argument& arg = local->argument() ;
162
163 {for ( unsigned int i = 0; i < v->size; ++i)
164 {
165 arg[i] = gsl_vector_get (v, i);
166 }
167 }
168
169 {for( unsigned int i = 0 ; i < v->size ; ++i )
170 {
171 gsl_vector_set(f, i, (*eqs[i])(arg ) ) ;
172 for (unsigned int j = 0; j < v->size; ++j)
173 {
174 Genfun::GENFUNCTION f1 = *(jac[i][j]) ;
175 gsl_matrix_set ( df , i , j , f1(arg) ) ;
176 }
177 }
178 }
179 return GSL_SUCCESS;
180 }
181 }
182
183 //=============================================================================
188 StatusCode EqSolver::solver (const Equations& funcs ,
189 Arg& arg ) const
190 //=============================================================================
191 {
192 using namespace Genfun;
193
194 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
195 arg.dimension() );
196 MsgStream log( msgSvc(), name() );
197
198 EqSolverMisc local (funcs, arg);
199
200 const gsl_multiroot_fdfsolver_type *T = m_type;
201 gsl_multiroot_fdfsolver *s;
202 int status = 0 ;
203 size_t iter = 0 ;
204
205 gsl_multiroot_function_fdf function;
206
207 function.f = &fun_gsl;
208 function.df = &dfun_gsl;
209 function.fdf = &fdfun_gsl;
210 function.n = vect.vector.size;
211 function.params = (void*) &local;
212
213 s = gsl_multiroot_fdfsolver_alloc(T, vect.vector.size);
214 gsl_multiroot_fdfsolver_set (s, &function, &vect.vector);
215
216 for (iter = 0; iter < m_max_iter; ++iter)
217 {
218 status = gsl_multiroot_fdfsolver_iterate (s);
219 if (status)
220 {
221 return Error
222 ("Error from gsl_gsl_multiroot_fdfsolver_iterate '"
223 + std::string(gsl_strerror(status)) + "'") ;
224 }
225
226 status = gsl_multiroot_test_residual (s->f,
227 m_norm_residual);
228
229 if ( status != GSL_CONTINUE ) { break; }
230 }
231
232 for (unsigned int i = 0; i < vect.vector.size; ++i)
233 {
234 gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
235 }
236
237 if (status == GSL_SUCCESS)
238 {
239 log << MSG::DEBUG
240 << "We stopped in the method on the " << iter
241 << " iteration (we have maximum " << m_max_iter
242 << " iterations)" << endmsg;
243 }
244 else if (status == GSL_CONTINUE && iter <= m_max_iter )
245 {
246 return Error ( "Method finished with '"
247 + std::string(gsl_strerror(status))
248 + "' error" );
249 }
250 else
251 {
252 return Error ( "Method finished with '" +
253 std::string(gsl_strerror(status))
254 + "' error" );
255 }
256
257 gsl_multiroot_fdfsolver_free (s);
258
259 if (status)
260 {
261 return Error ( "Method finished with '"
262 + std::string(gsl_strerror(status))
263 + "' error" );
264 }
265
266 return GSL_SUCCESS;
267 }
268
269 //=============================================================================
270 StatusCode EqSolver::initialize()
271
272 {
273 StatusCode sc = GaudiTool::initialize() ;
274
275 MsgStream log( msgSvc() , name() ) ;
276
277 if( sc.isFailure() )
278 {
279 return Error ("Could not initiliaze base class GaudiTool", sc);
280 }
281 /* The algorithm for multiional root-finding
282 (solving nonlinear systems with n equations in n unknowns)*/
283
284 if( "fdfsolver_hybridsj" == m_algType )
285 {
286 m_type = gsl_multiroot_fdfsolver_hybridsj ;
287 log << MSG::DEBUG
288 << "Root findind algoritms to be used: "
289 << "'gsl_multiroot_fdfsolver_hybridsj'"
290 << endmsg;
291 }
292 else if ( "fdfsolver_hybridj" == m_algType )
293 {
294 m_type = gsl_multiroot_fdfsolver_hybridj ;
295 log << MSG::DEBUG
296 << "Root findind algoritms to be used: "
297 << "'gsl_multiroot_fdfsolver_hybridj'"
298 << endmsg;
299 }
300 else if ( "fdfsolver_newton" == m_algType )
301 {
302 m_type = gsl_multiroot_fdfsolver_newton ;
303 log << MSG::DEBUG
304 << "Root findind algoritms to be used: "
305 << "'gsl_multiroot_fdfsolver_newton'"
306 << endmsg;
307 }
308 else if ( "fdfsolver_gnewton" == m_algType )
309 {
310 m_type = gsl_multiroot_fdfsolver_gnewton ;
311 log << MSG::DEBUG
312 << "Root findind algoritms to be used: "
313 << "'gsl_multiroot_fdfsolver_gnewton'"
314 << endmsg;
315 }
316 else
317 {
318 return Error(" Unknown algorith type '"
319 + std::string(m_algType) + "'");
320 }
321
322 return StatusCode::SUCCESS;
323 }
324 //=============================================================================
325 StatusCode EqSolver::finalize ()
326 {
327 StatusCode sc = GaudiTool::finalize() ;
328
329 MsgStream log( msgSvc() , name() ) ;
330
331 if ( sc.isFailure() )
332 {
333 return Error("Could not finaliaze base class GaudiTool", sc);
334 }
335 return StatusCode::SUCCESS;
336 }
337 //=============================================================================
338 EqSolver::~EqSolver( )
339 {}
340 //=============================================================================
MsgStream
Definition of the MsgStream class used to transmit messages.
Definition: MsgStream.h:24
gaudiComponentHelp.name
name
Definition: gaudiComponentHelp.py:28
MsgStream.h
Gaudi.Configuration.log
tuple log
Definition: Configuration.py:13
EqSolver::EqSolver
EqSolver()
default constructor is private
MSG::DEBUG
Definition: IMessageSvc.h:18
merge_files.sc
int sc
Definition: merge_files.py:167
EqSolver::EqSolverMisc::jacobi
const Jacobi & jacobi() const
Definition: EqSolver.h:57
StatusCode::isFailure
bool isFailure() const
Test for a status code of FAILURE.
Definition: StatusCode.h:86
DECLARE_COMPONENT
#define DECLARE_COMPONENT(type)
Definition: PluginService.h:36
EqSolver::solver
StatusCode solver(const Equations &funcs, Arg &arg) const override
Solving nonlinear system with N equations in N unknowns of the function "GenFunc".
Definition: EqSolver.cpp:188
std::string
STL class.
IOTest.N
int N
Definition: IOTest.py:90
EqSolver::m_algType
std::string m_algType
Definition: EqSolver.h:93
EqSolver
The simplest concrete implementation of IEqSolver interface.
Definition: EqSolver.h:23
EqSolver::initialize
StatusCode initialize() override
Overriding initialize.
Definition: EqSolver.cpp:270
EqSolver::EqSolverMisc::argument
const Arg & argument() const
Definition: EqSolver.h:54
StatusCode::SUCCESS
Definition: StatusCode.h:30
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
EqSolver::EqSolverMisc::~EqSolverMisc
~EqSolverMisc()
Definition: EqSolver.cpp:54
EqSolver::EqSolverMisc
Definition: EqSolver.h:43
EqSolver::m_norm_residual
double m_norm_residual
Definition: EqSolver.h:95
EqSolver::EqSolverMisc::equations
const Equations * equations() const
Definition: EqSolver.h:56
IInterface
Definition of the basic interface.
Definition: IInterface.h:234
EqSolver.h
GaudiTool::declareProperty
Property * declareProperty(const std::string &name, DataObjectHandle< T > &hndl, const std::string &doc="none") const
Definition: GaudiTool.h:834
EqSolver::m_max_iter
double m_max_iter
Definition: EqSolver.h:94
std::vector< Equations >
GaudiTool::Error
StatusCode Error(const std::string &msg, const StatusCode st=StatusCode::FAILURE, const size_t mx=10) const
Print the error message and return with the given StatusCode.
Definition: GaudiTool.h:689
AlgTool
Base class from which all the concrete tool classes should be derived.
Definition: AlgTool.h:45
GaudiTool::finalize
StatusCode finalize() override
standard finalization method
Definition: GaudiTool.cpp:179
EqSolver::finalize
StatusCode finalize() override
Definition: EqSolver.cpp:325
gaudirun.s
string s
Definition: gaudirun.py:245
CommonMessagingBase::msgSvc
SmartIF< IMessageSvc > & msgSvc() const
The standard message service.
Definition: CommonMessaging.h:73
EqSolver::m_type
const gsl_multiroot_fdfsolver_type * m_type
Definition: EqSolver.h:96
AlgTool::type
const std::string & type() const override
Retrieve type (concrete class) of the sub-algtool.
Definition: AlgTool.cpp:72
Genfun
CLHEP.
Definition: IEqSolver.h:13
EqSolver::~EqSolver
~EqSolver() override
Destructor.
Definition: EqSolver.cpp:338
ana.i
list i
Definition: ana.py:128
AlgTool::name
const std::string & name() const override
Retrieve full identifying name of the concrete tool object.
Definition: AlgTool.cpp:65
endmsg
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Definition: MsgStream.h:244
GaudiTool::initialize
StatusCode initialize() override
standard initialization method
Definition: GaudiTool.cpp:161
ForwardSchedulerSimpleTest.msgSvc
tuple msgSvc
Definition: ForwardSchedulerSimpleTest.py:112
AlgTool::parent
const IInterface * parent() const override
Retrieve parent of the sub-algtool.
Definition: AlgTool.cpp:79