| The Gaudi Framework v29r1 (5cd95fa2)

- 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 <stdio.h>
4 #include <stdlib.h>
5 // from Gaudi
6 #include "GaudiKernel/MsgStream.h"
7 // from GSL
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"
15
16 // local
17 #include "EqSolver.h"
18
26 // Declaration of the Tool Factory
27 DECLARE_COMPONENT( EqSolver )
28
29 // ============================================================================
30 // ============================================================================
31 EqSolver::EqSolverMisc::EqSolverMisc( const EqSolver::Equations& funcs, EqSolver::Arg& arg )
32 : m_argum( arg ), m_eqs( &funcs )
33 {
34 const size_t N = funcs.size();
35 for ( size_t i = 0; i < N; ++i ) {
36 Equations last;
37 for ( size_t j = 0; j < N; ++j ) {
38 Genfun::GENFUNCTION fij = funcs[i]->partial( j );
39 last.push_back( fij.clone() );
40 }
41 m_jac.push_back( last );
42 }
43 }
44 // ============================================================================
45
46 // ============================================================================
47 EqSolver::EqSolverMisc::~EqSolverMisc()
48 {
49 while ( !m_jac.empty() ) {
50 Equations& last = m_jac.back();
51 while ( !last.empty() ) {
52 delete last.back();
53 last.pop_back();
54 }
55 m_jac.pop_back();
56 }
57 }
58 // ============================================================================
59
60 namespace
61 {
62 using namespace Genfun;
63
66 int fun_gsl( const gsl_vector* v, void* params, gsl_vector* f )
67 {
68 EqSolver::EqSolverMisc* local = static_cast<EqSolver::EqSolverMisc*>( params );
69 const EqSolver::Equations& eqs = *( local->equations() );
70 Argument& arg = local->argument();
71
72 for ( unsigned int i = 0; i < v->size; ++i ) {
73 arg[i] = gsl_vector_get( v, i );
74 }
75
76 for ( unsigned int i = 0; i < v->size; ++i ) {
77 gsl_vector_set( f, i, ( *eqs[i] )( arg ) );
78 }
79 return GSL_SUCCESS;
80 }
81
84 int dfun_gsl( const gsl_vector* v, void* params, gsl_matrix* df )
85 {
86 EqSolver::EqSolverMisc* local = static_cast<EqSolver::EqSolverMisc*>( params );
87 const EqSolver::Jacobi& jac = local->jacobi();
88 Argument& arg = local->argument();
89
90 for ( unsigned int i = 0; i < v->size; ++i ) {
91 arg[i] = gsl_vector_get( v, i );
92 }
93
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 ) );
98 }
99 }
100 return GSL_SUCCESS;
101 }
102
106 int fdfun_gsl( const gsl_vector* v, void* params, gsl_vector* f, gsl_matrix* df )
107 {
108 EqSolver::EqSolverMisc* local = static_cast<EqSolver::EqSolverMisc*>( params );
109 const EqSolver::Equations& eqs = *( local->equations() );
110 const EqSolver::Jacobi& jac = local->jacobi();
111 Argument& arg = local->argument();
112
113 for ( unsigned int i = 0; i < v->size; ++i ) {
114 arg[i] = gsl_vector_get( v, i );
115 }
116
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 ) );
122 }
123 }
124 return GSL_SUCCESS;
125 }
126 }
127
128 //=============================================================================
133 StatusCode EqSolver::solver( const Equations& funcs, Arg& arg ) const
134 //=============================================================================
135 {
136 using namespace Genfun;
137
138 gsl_vector_view vect = gsl_vector_view_array( &arg[0], arg.dimension() );
139 EqSolverMisc local( funcs, arg );
140
141 const gsl_multiroot_fdfsolver_type* T = m_type;
142 gsl_multiroot_fdfsolver* s = nullptr;
143 int status = 0;
144 size_t iter = 0;
145
146 gsl_multiroot_function_fdf function;
147
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;
153
154 s = gsl_multiroot_fdfsolver_alloc( T, vect.vector.size );
155 gsl_multiroot_fdfsolver_set( s, &function, &vect.vector );
156
157 for ( iter = 0; iter < m_max_iter; ++iter ) {
158 status = gsl_multiroot_fdfsolver_iterate( s );
159 if ( status ) {
160 return Error( "Error from gsl_gsl_multiroot_fdfsolver_iterate '" + std::string( gsl_strerror( status ) ) + "'" );
161 }
162
163 status = gsl_multiroot_test_residual( s->f, m_norm_residual );
164
165 if ( status != GSL_CONTINUE ) {
166 break;
167 }
168 }
169
170 for ( unsigned int i = 0; i < vect.vector.size; ++i ) {
171 gsl_vector_set( &vect.vector, i, gsl_vector_get( s->x, i ) );
172 }
173
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" );
179 } else {
180 return Error( "Method finished with '" + std::string( gsl_strerror( status ) ) + "' error" );
181 }
182
183 gsl_multiroot_fdfsolver_free( s );
184
185 if ( status ) {
186 return Error( "Method finished with '" + std::string( gsl_strerror( status ) ) + "' error" );
187 }
188
189 return GSL_SUCCESS;
190 }
191
192 //=============================================================================
193 StatusCode EqSolver::initialize()
194 {
195 StatusCode sc = GaudiTool::initialize();
196 if ( sc.isFailure() ) {
197 return Error( "Could not initiliaze base class GaudiTool", sc );
198 }
199
200 /* The algorithm for multiional root-finding
201 (solving nonlinear systems with n equations in n unknowns)*/
202
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;
219 } else {
220 return Error( " Unknown algorithm type '" + std::string( m_algType ) + "'" );
221 }
222
223 return StatusCode::SUCCESS;
224 }
225 //=============================================================================
StatusCode::SUCCESS
Definition: StatusCode.h:31
std::vector::empty
T empty(T...args)
MsgStream.h
EqSolver::EqSolverMisc::jacobi
const Jacobi & jacobi() const
Definition: EqSolver.h:53
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:33
std::string
STL class.
EqSolver::m_algType
Gaudi::Property< std::string > m_algType
Definition: EqSolver.h:81
IOTest.N
int N
Definition: IOTest.py:101
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:83
EqSolver::initialize
StatusCode initialize() override
Overriding initialize.
Definition: EqSolver.cpp:193
EqSolver::EqSolverMisc::argument
const Arg & argument() const
Definition: EqSolver.h:50
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:28
EqSolver::EqSolverMisc::~EqSolverMisc
~EqSolverMisc()
Definition: EqSolver.cpp:47
EqSolver::EqSolverMisc
Definition: EqSolver.h:39
EqSolver::EqSolverMisc::equations
const Equations * equations() const
Definition: EqSolver.h:52
EqSolver.h
std::vector::pop_back
T pop_back(T...args)
EqSolver::EqSolverMisc::m_jac
Jacobi m_jac
Definition: EqSolver.h:66
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:133
std::vector< Equations >
CommonMessagingBase::debug
MsgStream & debug() const
shortcut for the method msgStream(MSG::DEBUG)
Definition: CommonMessaging.h:115
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
std::vector::back
T back(T...args)
gaudirun.s
string s
Definition: gaudirun.py:253
EqSolver::m_type
const gsl_multiroot_fdfsolver_type * m_type
Definition: EqSolver.h:86
Genfun
CLHEP.
Definition: IEqSolver.h:13
endmsg
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Definition: MsgStream.h:209
EqSolver::m_norm_residual
Gaudi::Property< double > m_norm_residual
Definition: EqSolver.h:84
GaudiTool::initialize
StatusCode initialize() override
standard initialization method
Definition: GaudiTool.cpp:151