6 #include "GaudiKernel/MsgStream.h"
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"
39 const size_t N = funcs.size () ;
40 for(
size_t i = 0 ;
i <
N ; ++
i )
43 for(
size_t j = 0 ; j <
N ; ++j )
45 Genfun::GENFUNCTION fij = funcs[
i]->partial(j);
46 last.push_back( fij.clone() ) ;
48 m_jac.push_back( last );
56 while( !m_jac.empty() )
58 Equations& last = m_jac.back() ;
59 while( !last.empty() )
73 const std::string&
name,
83 declareProperty (
"Algorithm",
m_algType );
96 int fun_gsl (
const gsl_vector* v ,
102 const EqSolver::Equations& eqs = *(local->
equations()) ;
105 {
for (
unsigned int i = 0;
i < v->size; ++
i)
107 arg[
i] = gsl_vector_get (v,
i);
111 {
for (
unsigned int i = 0;
i < v->size; ++
i)
113 gsl_vector_set(f,
i, (*eqs[
i])(arg));
122 int dfun_gsl (
const gsl_vector* v ,
129 Argument& arg = local -> argument() ;
131 {
for (
unsigned int i = 0;
i < v->size; ++
i)
133 arg[
i] = gsl_vector_get (v,
i);
137 {
for(
unsigned int i = 0 ;
i < v->size ; ++
i )
139 for (
unsigned int j = 0; j < v->size; ++j)
141 Genfun::GENFUNCTION f = *(jac[
i][j]) ;
142 gsl_matrix_set ( df ,
i , j , f ( arg ) ) ;
152 int fdfun_gsl (
const gsl_vector* v ,
159 const EqSolver::Equations& eqs = *(local->
equations()) ;
163 {
for (
unsigned int i = 0;
i < v->size; ++
i)
165 arg[
i] = gsl_vector_get (v,
i);
169 {
for(
unsigned int i = 0 ;
i < v->size ; ++
i )
171 gsl_vector_set(f,
i, (*eqs[
i])(arg ) ) ;
172 for (
unsigned int j = 0; j < v->size; ++j)
174 Genfun::GENFUNCTION f1 = *(jac[
i][j]) ;
175 gsl_matrix_set ( df , i , j , f1(arg) ) ;
194 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
200 const gsl_multiroot_fdfsolver_type *T = m_type;
201 gsl_multiroot_fdfsolver *
s;
205 gsl_multiroot_function_fdf
function;
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;
213 s = gsl_multiroot_fdfsolver_alloc(T, vect.vector.size);
214 gsl_multiroot_fdfsolver_set (s, &
function, &vect.vector);
216 for (iter = 0; iter < m_max_iter; ++iter)
218 status = gsl_multiroot_fdfsolver_iterate (s);
222 (
"Error from gsl_gsl_multiroot_fdfsolver_iterate '"
223 + std::string(gsl_strerror(status)) +
"'") ;
226 status = gsl_multiroot_test_residual (s->f,
229 if ( status != GSL_CONTINUE ) {
break; }
232 for (
unsigned int i = 0;
i < vect.vector.size; ++
i)
234 gsl_vector_set (&vect.vector,
i, gsl_vector_get (s->x,
i));
237 if (status == GSL_SUCCESS)
240 <<
"We stopped in the method on the " << iter
241 <<
" iteration (we have maximum " << m_max_iter
242 <<
" iterations)" <<
endmsg;
244 else if (status == GSL_CONTINUE && iter <= m_max_iter )
246 return Error (
"Method finished with '"
247 + std::string(gsl_strerror(status))
252 return Error (
"Method finished with '" +
253 std::string(gsl_strerror(status))
257 gsl_multiroot_fdfsolver_free (s);
261 return Error (
"Method finished with '"
262 + std::string(gsl_strerror(status))
279 return Error (
"Could not initiliaze base class GaudiTool", sc);
286 m_type = gsl_multiroot_fdfsolver_hybridsj ;
288 <<
"Root findind algoritms to be used: "
289 <<
"'gsl_multiroot_fdfsolver_hybridsj'"
292 else if (
"fdfsolver_hybridj" ==
m_algType )
294 m_type = gsl_multiroot_fdfsolver_hybridj ;
296 <<
"Root findind algoritms to be used: "
297 <<
"'gsl_multiroot_fdfsolver_hybridj'"
300 else if (
"fdfsolver_newton" ==
m_algType )
302 m_type = gsl_multiroot_fdfsolver_newton ;
304 <<
"Root findind algoritms to be used: "
305 <<
"'gsl_multiroot_fdfsolver_newton'"
308 else if (
"fdfsolver_gnewton" ==
m_algType )
310 m_type = gsl_multiroot_fdfsolver_gnewton ;
312 <<
"Root findind algoritms to be used: "
313 <<
"'gsl_multiroot_fdfsolver_gnewton'"
318 return Error(
" Unknown algorith type '"
333 return Error(
"Could not finaliaze base class GaudiTool", sc);
std::vector< Equations > Jacobi
Definition of the MsgStream class used to transmit messages.
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
EqSolver()
default constructor is private
const Jacobi & jacobi() const
bool isFailure() const
Test for a status code of FAILURE.
StatusCode solver(const Equations &funcs, Arg &arg) const override
Solving nonlinear system with N equations in N unknowns of the function "GenFunc".
The simplest concrete implementation of IEqSolver interface.
StatusCode initialize() override
Overriding initialize.
const Arg & argument() const
This class is used for returning status codes from appropriate routines.
#define DECLARE_COMPONENT(type)
const Equations * equations() const
Definition of the basic interface.
Base class used to extend a class implementing other interfaces.
StatusCode finalize() override
const gsl_multiroot_fdfsolver_type * m_type
~EqSolver() override
Destructor.