| The Gaudi Framework v28r2p1 (f1a77ff4)

- 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 namespace
70 {
71 using namespace Genfun;
72
75 int fun_gsl ( const gsl_vector* v ,
76 void* params ,
77 gsl_vector* f )
78 {
79 EqSolver::EqSolverMisc* local =
80 static_cast <EqSolver::EqSolverMisc*> (params);
81 const EqSolver::Equations& eqs = *(local->equations()) ;
82 Argument& arg = local->argument() ;
83
84 {for ( unsigned int i = 0; i < v->size; ++i)
85 {
86 arg[i] = gsl_vector_get (v, i);
87 }
88 }
89
90 {for ( unsigned int i = 0; i < v->size; ++i)
91 {
92 gsl_vector_set(f, i, (*eqs[i])(arg));
93 }
94 }
95 return GSL_SUCCESS;
96 }
97
98
101 int dfun_gsl ( const gsl_vector* v ,
102 void* params ,
103 gsl_matrix* df )
104 {
105 EqSolver::EqSolverMisc* local =
106 static_cast <EqSolver::EqSolverMisc*> (params);
107 const EqSolver::Jacobi& jac = local -> jacobi() ;
108 Argument& arg = local -> argument() ;
109
110 {for ( unsigned int i = 0; i < v->size; ++i)
111 {
112 arg[i] = gsl_vector_get (v, i);
113 }
114 }
115
116 {for( unsigned int i = 0 ; i < v->size ; ++i )
117 {
118 for (unsigned int j = 0; j < v->size; ++j)
119 {
120 Genfun::GENFUNCTION f = *(jac[i][j]) ;
121 gsl_matrix_set ( df , i , j , f ( arg ) ) ;
122 }
123 }
124 }
125 return GSL_SUCCESS;
126 }
127
131 int fdfun_gsl ( const gsl_vector* v ,
132 void* params ,
133 gsl_vector* f ,
134 gsl_matrix* df )
135 {
136 EqSolver::EqSolverMisc* local =
137 static_cast <EqSolver::EqSolverMisc*> (params);
138 const EqSolver::Equations& eqs = *(local->equations()) ;
139 const EqSolver::Jacobi& jac = local->jacobi() ;
140 Argument& arg = local->argument() ;
141
142 {for ( unsigned int i = 0; i < v->size; ++i)
143 {
144 arg[i] = gsl_vector_get (v, i);
145 }
146 }
147
148 {for( unsigned int i = 0 ; i < v->size ; ++i )
149 {
150 gsl_vector_set(f, i, (*eqs[i])(arg ) ) ;
151 for (unsigned int j = 0; j < v->size; ++j)
152 {
153 Genfun::GENFUNCTION f1 = *(jac[i][j]) ;
154 gsl_matrix_set ( df , i , j , f1(arg) ) ;
155 }
156 }
157 }
158 return GSL_SUCCESS;
159 }
160 }
161
162 //=============================================================================
167 StatusCode EqSolver::solver (const Equations& funcs ,
168 Arg& arg ) const
169 //=============================================================================
170 {
171 using namespace Genfun;
172
173 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
174 arg.dimension() );
175 MsgStream log( msgSvc(), name() );
176
177 EqSolverMisc local (funcs, arg);
178
179 const gsl_multiroot_fdfsolver_type *T = m_type;
180 gsl_multiroot_fdfsolver *s;
181 int status = 0 ;
182 size_t iter = 0 ;
183
184 gsl_multiroot_function_fdf function;
185
186 function.f = &fun_gsl;
187 function.df = &dfun_gsl;
188 function.fdf = &fdfun_gsl;
189 function.n = vect.vector.size;
190 function.params = (void*) &local;
191
192 s = gsl_multiroot_fdfsolver_alloc(T, vect.vector.size);
193 gsl_multiroot_fdfsolver_set (s, &function, &vect.vector);
194
195 for (iter = 0; iter < m_max_iter; ++iter)
196 {
197 status = gsl_multiroot_fdfsolver_iterate (s);
198 if (status)
199 {
200 return Error
201 ("Error from gsl_gsl_multiroot_fdfsolver_iterate '"
202 + std::string(gsl_strerror(status)) + "'") ;
203 }
204
205 status = gsl_multiroot_test_residual (s->f,
206 m_norm_residual);
207
208 if ( status != GSL_CONTINUE ) { break; }
209 }
210
211 for (unsigned int i = 0; i < vect.vector.size; ++i)
212 {
213 gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
214 }
215
216 if (status == GSL_SUCCESS)
217 {
218 log << MSG::DEBUG
219 << "We stopped in the method on the " << iter
220 << " iteration (we have maximum " << m_max_iter
221 << " iterations)" << endmsg;
222 }
223 else if (status == GSL_CONTINUE && iter <= m_max_iter )
224 {
225 return Error ( "Method finished with '"
226 + std::string(gsl_strerror(status))
227 + "' error" );
228 }
229 else
230 {
231 return Error ( "Method finished with '" +
232 std::string(gsl_strerror(status))
233 + "' error" );
234 }
235
236 gsl_multiroot_fdfsolver_free (s);
237
238 if (status)
239 {
240 return Error ( "Method finished with '"
241 + std::string(gsl_strerror(status))
242 + "' error" );
243 }
244
245 return GSL_SUCCESS;
246 }
247
248 //=============================================================================
249 StatusCode EqSolver::initialize()
250
251 {
252 StatusCode sc = GaudiTool::initialize() ;
253
254 MsgStream log( msgSvc() , name() ) ;
255
256 if( sc.isFailure() )
257 {
258 return Error ("Could not initiliaze base class GaudiTool", sc);
259 }
260 /* The algorithm for multiional root-finding
261 (solving nonlinear systems with n equations in n unknowns)*/
262
263 if( "fdfsolver_hybridsj" == m_algType )
264 {
265 m_type = gsl_multiroot_fdfsolver_hybridsj ;
266 log << MSG::DEBUG
267 << "Root findind algoritms to be used: "
268 << "'gsl_multiroot_fdfsolver_hybridsj'"
269 << endmsg;
270 }
271 else if ( "fdfsolver_hybridj" == m_algType )
272 {
273 m_type = gsl_multiroot_fdfsolver_hybridj ;
274 log << MSG::DEBUG
275 << "Root findind algoritms to be used: "
276 << "'gsl_multiroot_fdfsolver_hybridj'"
277 << endmsg;
278 }
279 else if ( "fdfsolver_newton" == m_algType )
280 {
281 m_type = gsl_multiroot_fdfsolver_newton ;
282 log << MSG::DEBUG
283 << "Root findind algoritms to be used: "
284 << "'gsl_multiroot_fdfsolver_newton'"
285 << endmsg;
286 }
287 else if ( "fdfsolver_gnewton" == m_algType )
288 {
289 m_type = gsl_multiroot_fdfsolver_gnewton ;
290 log << MSG::DEBUG
291 << "Root findind algoritms to be used: "
292 << "'gsl_multiroot_fdfsolver_gnewton'"
293 << endmsg;
294 }
295 else
296 {
297 return Error(" Unknown algorith type '"
298 + std::string(m_algType) + "'");
299 }
300
301 return StatusCode::SUCCESS;
302 }
303 //=============================================================================
304 StatusCode EqSolver::finalize ()
305 {
306 StatusCode sc = GaudiTool::finalize() ;
307
308 MsgStream log( msgSvc() , name() ) ;
309
310 if ( sc.isFailure() )
311 {
312 return Error("Could not finaliaze base class GaudiTool", sc);
313 }
314 return StatusCode::SUCCESS;
315 }
316 //=============================================================================
317 EqSolver::~EqSolver( )
318 {}
319 //=============================================================================
MsgStream
Definition of the MsgStream class used to transmit messages.
Definition: MsgStream.h:24
StatusCode::SUCCESS
Definition: StatusCode.h:30
MsgStream.h
MSG::DEBUG
Definition: IMessageSvc.h:18
EqSolver::EqSolverMisc::jacobi
const Jacobi & jacobi() const
Definition: EqSolver.h:56
StatusCode::isFailure
bool isFailure() const
Test for a status code of FAILURE.
Definition: StatusCode.h:84
DECLARE_COMPONENT
#define DECLARE_COMPONENT(type)
Definition: PluginService.h:36
std::string
STL class.
EqSolver::m_algType
Gaudi::Property< std::string > m_algType
Definition: EqSolver.h:83
IOTest.N
int N
Definition: IOTest.py:90
EqSolver
The simplest concrete implementation of IEqSolver interface.
Definition: EqSolver.h:23
EqSolver::m_max_iter
Gaudi::Property< double > m_max_iter
Definition: EqSolver.h:85
EqSolver::initialize
StatusCode initialize() override
Overriding initialize.
Definition: EqSolver.cpp:249
EqSolver::EqSolverMisc::argument
const Arg & argument() const
Definition: EqSolver.h:53
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:42
EqSolver::EqSolverMisc::equations
const Equations * equations() const
Definition: EqSolver.h:55
EqSolver.h
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:167
AlgTool::name
const std::string & name() const override
Retrieve full identifying name of the concrete tool object.
Definition: AlgTool.cpp:63
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:671
GaudiTool::finalize
StatusCode finalize() override
standard finalization method
Definition: GaudiTool.cpp:175
EqSolver::finalize
StatusCode finalize() override
Definition: EqSolver.cpp:304
gaudirun.s
string s
Definition: gaudirun.py:245
CommonMessagingBase::msgSvc
SmartIF< IMessageSvc > & msgSvc() const
The standard message service.
Definition: CommonMessaging.h:70
EqSolver::m_type
const gsl_multiroot_fdfsolver_type * m_type
Definition: EqSolver.h:88
Genfun
CLHEP.
Definition: IEqSolver.h:13
EqSolver::~EqSolver
~EqSolver() override
Destructor.
Definition: EqSolver.cpp:317
endmsg
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Definition: MsgStream.h:244
EqSolver::m_norm_residual
Gaudi::Property< double > m_norm_residual
Definition: EqSolver.h:86
GaudiTool::initialize
StatusCode initialize() override
standard initialization method
Definition: GaudiTool.cpp:157
Gaudi.Configuration.log
log
Definition: Configuration.py:13