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
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 // ============================================================================
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 //=============================================================================
271 
272 {
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 //=============================================================================
326 {
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 //=============================================================================
339 {}
340 //=============================================================================
std::vector< Equations > Jacobi
Definition: EqSolver.h:25
Definition of the MsgStream class used to transmit messages.
Definition: MsgStream.h:24
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Definition: MsgStream.h:244
EqSolver()
default constructor is private
const Jacobi & jacobi() const
Definition: EqSolver.h:56
bool isFailure() const
Test for a status code of FAILURE.
Definition: StatusCode.h:86
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
int N
Definition: IOTest.py:90
std::string m_algType
Definition: EqSolver.h:92
The simplest concrete implementation of IEqSolver interface.
Definition: EqSolver.h:23
StatusCode initialize() override
Overriding initialize.
Definition: EqSolver.cpp:270
const Arg & argument() const
Definition: EqSolver.h:53
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
double m_norm_residual
Definition: EqSolver.h:94
#define DECLARE_COMPONENT(type)
Definition: PluginService.h:36
const Equations * equations() const
Definition: EqSolver.h:55
Definition of the basic interface.
Definition: IInterface.h:234
double m_max_iter
Definition: EqSolver.h:93
Base class used to extend a class implementing other interfaces.
Definition: extends.h:10
StatusCode finalize() override
standard finalization method
Definition: GaudiTool.cpp:179
StatusCode finalize() override
Definition: EqSolver.cpp:325
string s
Definition: gaudirun.py:245
const gsl_multiroot_fdfsolver_type * m_type
Definition: EqSolver.h:95
CLHEP.
Definition: IEqSolver.h:13
~EqSolver() override
Destructor.
Definition: EqSolver.cpp:338
list i
Definition: ana.py:128
StatusCode initialize() override
standard initialization method
Definition: GaudiTool.cpp:161
string type
Definition: gaudirun.py:151