8 #include "gsl/gsl_blas.h" 9 #include "gsl/gsl_errno.h" 10 #include "gsl/gsl_linalg.h" 11 #include "gsl/gsl_math.h" 12 #include "gsl/gsl_matrix.h" 13 #include "gsl/gsl_multimin.h" 14 #include "gsl/gsl_vector.h" 32 : m_argum( arg ), m_eqs( &funcs )
34 const size_t N = funcs.size();
35 for (
size_t i = 0; i <
N; ++i ) {
37 for (
size_t j = 0; j <
N; ++j ) {
38 Genfun::GENFUNCTION fij = funcs[i]->partial( j );
39 last.push_back( fij.clone() );
41 m_jac.push_back( last );
51 while ( !last.empty() ) {
66 int fun_gsl(
const gsl_vector* v,
void* params, gsl_vector* f )
69 const EqSolver::Equations& eqs = *( local->
equations() );
72 for (
unsigned int i = 0; i < v->size; ++i ) {
73 arg[i] = gsl_vector_get( v, i );
76 for (
unsigned int i = 0; i < v->size; ++i ) {
77 gsl_vector_set( f, i, ( *eqs[i] )( arg ) );
84 int dfun_gsl(
const gsl_vector* v,
void* params, gsl_matrix* df )
90 for (
unsigned int i = 0; i < v->size; ++i ) {
91 arg[i] = gsl_vector_get( v, i );
94 for (
unsigned int i = 0; i < v->size; ++i ) {
95 for (
unsigned int j = 0; j < v->size; ++j ) {
96 Genfun::GENFUNCTION f = *( jac[i][j] );
97 gsl_matrix_set( df, i, j, f( arg ) );
106 int fdfun_gsl(
const gsl_vector* v,
void* params, gsl_vector* f, gsl_matrix* df )
109 const EqSolver::Equations& eqs = *( local->
equations() );
113 for (
unsigned int i = 0; i < v->size; ++i ) {
114 arg[i] = gsl_vector_get( v, i );
117 for (
unsigned int i = 0; i < v->size; ++i ) {
118 gsl_vector_set( f, i, ( *eqs[i] )( arg ) );
119 for (
unsigned int j = 0; j < v->size; ++j ) {
120 Genfun::GENFUNCTION f1 = *( jac[i][j] );
121 gsl_matrix_set( df, i, j, f1( arg ) );
138 gsl_vector_view vect = gsl_vector_view_array( &arg[0], arg.dimension() );
141 const gsl_multiroot_fdfsolver_type* T =
m_type;
142 gsl_multiroot_fdfsolver*
s =
nullptr;
146 gsl_multiroot_function_fdf
function;
148 function.f = &fun_gsl;
149 function.df = &dfun_gsl;
150 function.fdf = &fdfun_gsl;
151 function.n = vect.vector.size;
152 function.params = &local;
154 s = gsl_multiroot_fdfsolver_alloc( T, vect.vector.size );
155 gsl_multiroot_fdfsolver_set( s, &
function, &vect.vector );
158 status = gsl_multiroot_fdfsolver_iterate( s );
160 return Error(
"Error from gsl_gsl_multiroot_fdfsolver_iterate '" +
std::string( gsl_strerror( status ) ) +
"'" );
165 if ( status != GSL_CONTINUE ) {
170 for (
unsigned int i = 0; i < vect.vector.size; ++i ) {
171 gsl_vector_set( &vect.vector, i, gsl_vector_get( s->x, i ) );
174 if ( status == GSL_SUCCESS ) {
175 debug() <<
"We stopped in the method on the " << iter <<
" iteration (we have maximum " << m_max_iter
176 <<
" iterations)" <<
endmsg;
177 }
else if ( status == GSL_CONTINUE && iter <= m_max_iter ) {
178 return Error(
"Method finished with '" +
std::string( gsl_strerror( status ) ) +
"' error" );
180 return Error(
"Method finished with '" +
std::string( gsl_strerror( status ) ) +
"' error" );
183 gsl_multiroot_fdfsolver_free( s );
186 return Error(
"Method finished with '" +
std::string( gsl_strerror( status ) ) +
"' error" );
197 return Error(
"Could not initiliaze base class GaudiTool", sc );
203 if (
"fdfsolver_hybridsj" ==
m_algType ) {
204 m_type = gsl_multiroot_fdfsolver_hybridsj;
205 debug() <<
"Root finding algorithm to be used: " 206 <<
"'gsl_multiroot_fdfsolver_hybridsj'" <<
endmsg;
207 }
else if (
"fdfsolver_hybridj" ==
m_algType ) {
208 m_type = gsl_multiroot_fdfsolver_hybridj;
209 debug() <<
"Root finding algorithm to be used: " 210 <<
"'gsl_multiroot_fdfsolver_hybridj'" <<
endmsg;
211 }
else if (
"fdfsolver_newton" ==
m_algType ) {
212 m_type = gsl_multiroot_fdfsolver_newton;
213 debug() <<
"Root findind algorithm to be used: " 214 <<
"'gsl_multiroot_fdfsolver_newton'" <<
endmsg;
215 }
else if (
"fdfsolver_gnewton" ==
m_algType ) {
216 m_type = gsl_multiroot_fdfsolver_gnewton;
217 debug() <<
"Root findind algorithm to be used: " 218 <<
"'gsl_multiroot_fdfsolver_gnewton'" <<
endmsg;
const Jacobi & jacobi() const
#define DECLARE_COMPONENT(type)
Gaudi::Property< std::string > m_algType
The simplest concrete implementation of IEqSolver interface.
Gaudi::Property< double > m_max_iter
StatusCode initialize() override
Overriding initialize.
const Arg & argument() const
This class is used for returning status codes from appropriate routines.
const Equations * equations() const
constexpr static const auto SUCCESS
StatusCode solver(const Equations &funcs, Arg &arg) const override
Solving nonlinear system with N equations in N unknowns of the function "GenFunc".
MsgStream & debug() const
shortcut for the method msgStream(MSG::DEBUG)
const gsl_multiroot_fdfsolver_type * m_type
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Gaudi::Property< double > m_norm_residual