| # Gaudi Framework, version v21r11 |

Home |  | Generated: 30 Sep 2010

# EqSolver.cpp

Go to the documentation of this file.

    00001 // $Id: EqSolver.cpp,v 1.4 2006/01/10 20:00:05 hmd Exp $
    00002 // ============================================================================
    00003 // Include files
    00004
    00005 #include <stdlib.h>
    00006 #include <stdio.h>
    00007 // from Gaudi
    00008 #include "GaudiKernel/ToolFactory.h"
    00009 #include "GaudiKernel/MsgStream.h"
    00010 //from GSL
    00011 #include "gsl/gsl_vector.h"
    00012 #include "gsl/gsl_multimin.h"
    00013 #include "gsl/gsl_math.h"
    00014 #include "gsl/gsl_matrix.h"
    00015 #include "gsl/gsl_linalg.h"
    00016 #include "gsl/gsl_blas.h"
    00017 #include "gsl/gsl_errno.h"
    00018
    00019 // local
    00020 #include "EqSolver.h"
    00021
    00029 // Declaration of the Tool Factory
    00030 DECLARE_TOOL_FACTORY(EqSolver)
    00031
    00032
    00033 // ============================================================================
    00034 // ============================================================================
    00035 EqSolver::EqSolverMisc::EqSolverMisc
    00036 ( const EqSolver::Equations& funcs ,
    00037   EqSolver::Arg&             arg   )
    00038   : m_argum ( arg    )
    00039   , m_eqs   ( &funcs )
    00040   , m_jac   ()
    00041 {
    00042   const size_t N = funcs.size () ;
    00043   for( size_t i = 0 ; i < N ; ++i )
    00044     {
    00045       Equations last;
    00046       for( size_t j = 0 ; j < N ; ++j )
    00047         {
    00048           Genfun::GENFUNCTION fij = funcs[i]->partial(j);
    00049           last.push_back( fij.clone() ) ;
    00050         }
    00051       m_jac.push_back( last );
    00052     }
    00053 }
    00054 // ============================================================================
    00055
    00056 // ============================================================================
    00057 EqSolver::EqSolverMisc::~EqSolverMisc()
    00058 {
    00059   while( !m_jac.empty() )
    00060     {
    00061       Equations& last = m_jac.back() ;
    00062       while( !last.empty() )
    00063         {
    00064           delete last.back() ;
    00065           last.pop_back () ;
    00066         }
    00067       m_jac.pop_back();
    00068     }
    00069 }
    00070 // ============================================================================
    00071
    00072 //=============================================================================
    00073 // Standard constructor, initializes variables
    00074 //=============================================================================
    00075 EqSolver::EqSolver( const std::string& type,
    00076                     const std::string& name,
    00077                     const IInterface* parent )
    00078   : base_class      ( type, name , parent )
    00079   , m_algType       ( "fdfsolver_hybridsj" )
    00080   , m_max_iter      ( 1000                 )
    00081   , m_norm_residual ( 1.0e-7               )
    00082   , m_type          ( 0                    )
    00083
    00084 {
    00086   declareProperty ( "Algorithm", m_algType       );
    00088   declareProperty ( "Iteration", m_max_iter      );
    00090   declareProperty ( "Residual" , m_norm_residual );
    00091 }
    00092
    00093 namespace
    00094 {
    00095   using namespace Genfun;
    00096
    00099   int fun_gsl  ( const gsl_vector* v      ,
    00100                  void*             params ,
    00101                  gsl_vector*       f      )
    00102   {
    00103     EqSolver::EqSolverMisc* local =
    00104       static_cast <EqSolver::EqSolverMisc*> (params);
    00105     const EqSolver::Equations& eqs = *(local->equations()) ;
    00106     Argument&                  arg = local->argument()  ;
    00107
    00108     {for ( unsigned int i = 0; i < v->size; ++i)
    00109       {
    00110         arg[i] = gsl_vector_get (v, i);
    00111       }
    00112     }
    00113
    00114     {for ( unsigned int i = 0; i < v->size; ++i)
    00115       {
    00116         gsl_vector_set(f, i, (*eqs[i])(arg));
    00117       }
    00118     }
    00119     return GSL_SUCCESS;
    00120   }
    00121
    00122
    00125   int dfun_gsl  ( const gsl_vector* v      ,
    00126                   void*             params ,
    00127                   gsl_matrix*       df     )
    00128   {
    00129     EqSolver::EqSolverMisc* local =
    00130       static_cast <EqSolver::EqSolverMisc*> (params);
    00131     const EqSolver::Jacobi&     jac = local -> jacobi()    ;
    00132     Argument&                   arg = local -> argument()  ;
    00133
    00134     {for ( unsigned int i = 0; i < v->size; ++i)
    00135       {
    00136         arg[i] = gsl_vector_get (v, i);
    00137       }
    00138     }
    00139
    00140     {for( unsigned int i = 0 ; i < v->size ; ++i )
    00141       {
    00142         for (unsigned int j = 0; j < v->size; ++j)
    00143           {
    00144             Genfun::GENFUNCTION  f = *(jac[i][j]) ;
    00145             gsl_matrix_set      ( df , i , j , f ( arg ) ) ;
    00146           }
    00147       }
    00148     }
    00149     return GSL_SUCCESS;
    00150   }
    00151
    00155   int fdfun_gsl ( const gsl_vector* v ,
    00156                   void*       params  ,
    00157                   gsl_vector* f       ,
    00158                   gsl_matrix* df      )
    00159   {
    00160     EqSolver::EqSolverMisc* local =
    00161       static_cast <EqSolver::EqSolverMisc*> (params);
    00162     const EqSolver::Equations&  eqs = *(local->equations()) ;
    00163     const EqSolver::Jacobi&     jac = local->jacobi()    ;
    00164     Argument&                   arg = local->argument()  ;
    00165
    00166     {for ( unsigned int i = 0; i < v->size; ++i)
    00167       {
    00168         arg[i] = gsl_vector_get (v, i);
    00169       }
    00170     }
    00171
    00172     {for( unsigned int i = 0 ; i < v->size ; ++i )
    00173       {
    00174         gsl_vector_set(f, i, (*eqs[i])(arg ) ) ;
    00175         for (unsigned int j = 0; j < v->size; ++j)
    00176           {
    00177             Genfun::GENFUNCTION  f1 = *(jac[i][j]) ;
    00178             gsl_matrix_set      ( df , i , j , f1(arg) ) ;
    00179           }
    00180       }
    00181     }
    00182     return GSL_SUCCESS;
    00183   }
    00184 }
    00185
    00186 //=============================================================================
    00191 StatusCode EqSolver::solver (const Equations&  funcs ,
    00192                              Arg&              arg   ) const
    00193 //=============================================================================
    00194 {
    00195   using namespace Genfun;
    00196
    00197   gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
    00198                                                  arg.dimension() );
    00199   MsgStream log( msgSvc(), name() );
    00200
    00201   EqSolverMisc local (funcs, arg);
    00202
    00203   const gsl_multiroot_fdfsolver_type *T = m_type;
    00204   gsl_multiroot_fdfsolver *s;
    00205   int    status = 0 ;
    00206   size_t iter   = 0 ;
    00207
    00208   gsl_multiroot_function_fdf function;
    00209
    00210   function.f = &fun_gsl;
    00211   function.df = &dfun_gsl;
    00212   function.fdf = &fdfun_gsl;
    00213   function.n = vect.vector.size;
    00214   function.params = (void*) &local;
    00215
    00216   s = gsl_multiroot_fdfsolver_alloc(T, vect.vector.size);
    00217   gsl_multiroot_fdfsolver_set (s, &function, &vect.vector);
    00218
    00219   for (iter = 0; iter < m_max_iter; ++iter)
    00220     {
    00221       status = gsl_multiroot_fdfsolver_iterate (s);
    00222       if (status)
    00223         {
    00224           return Error
    00225             ("Error from gsl_gsl_multiroot_fdfsolver_iterate '"
    00226              + std::string(gsl_strerror(status)) + "'") ;
    00227         }
    00228
    00229       status = gsl_multiroot_test_residual (s->f,
    00230                                             m_norm_residual);
    00231
    00232       if ( status != GSL_CONTINUE ) { break; }
    00233     }
    00234
    00235   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00236     {
    00237       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00238     }
    00239
    00240   if (status == GSL_SUCCESS)
    00241     {
    00242       log << MSG::DEBUG
    00243           << "We stopped in the method on the " << iter
    00244           << " iteration (we have maximum "     << m_max_iter
    00245           << " iterations)"                     << endmsg;
    00246     }
    00247   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00248     {
    00249       return Error ( "Method finished with '"
    00250                      + std::string(gsl_strerror(status))
    00251                      +  "' error" );
    00252     }
    00253   else
    00254     {
    00255       return Error ( "Method finished with '" +
    00256                      std::string(gsl_strerror(status))
    00257                      +  "' error" );
    00258     }
    00259
    00260   gsl_multiroot_fdfsolver_free (s);
    00261
    00262   if (status)
    00263     {
    00264       return Error ( "Method finished with '"
    00265                      + std::string(gsl_strerror(status))
    00266                      +  "' error" );
    00267     }
    00268
    00269   return GSL_SUCCESS;
    00270 }
    00271
    00272 //=============================================================================
    00273 StatusCode  EqSolver::initialize()
    00274
    00275 {
    00276   StatusCode sc = GaudiTool::initialize() ;
    00277
    00278   MsgStream log( msgSvc() , name() ) ;
    00279
    00280   if( sc.isFailure() )
    00281     {
    00282       return Error ("Could not initiliaze base class GaudiTool", sc);
    00283     }
    00284   /* The algorithm for multiional root-finding
    00285      (solving nonlinear systems with n equations in n unknowns)*/
    00286
    00287   if(       "fdfsolver_hybridsj" == m_algType )
    00288     {
    00289       m_type = gsl_multiroot_fdfsolver_hybridsj  ;
    00290       log << MSG::DEBUG
    00291           << "Root findind algoritms to be used: "
    00292           << "'gsl_multiroot_fdfsolver_hybridsj'"
    00293           << endmsg;
    00294     }
    00295   else if ( "fdfsolver_hybridj" == m_algType )
    00296     {
    00297       m_type = gsl_multiroot_fdfsolver_hybridj   ;
    00298       log << MSG::DEBUG
    00299           << "Root findind algoritms to be used: "
    00300           << "'gsl_multiroot_fdfsolver_hybridj'"
    00301           << endmsg;
    00302     }
    00303   else if ( "fdfsolver_newton" == m_algType )
    00304     {
    00305       m_type = gsl_multiroot_fdfsolver_newton    ;
    00306       log << MSG::DEBUG
    00307           << "Root findind algoritms to be used: "
    00308           << "'gsl_multiroot_fdfsolver_newton'"
    00309           << endmsg;
    00310     }
    00311   else if ( "fdfsolver_gnewton" == m_algType )
    00312     {
    00313       m_type = gsl_multiroot_fdfsolver_gnewton   ;
    00314       log << MSG::DEBUG
    00315           << "Root findind algoritms to be used: "
    00316           << "'gsl_multiroot_fdfsolver_gnewton'"
    00317           << endmsg;
    00318     }
    00319   else
    00320     {
    00321       return Error(" Unknown algorith type '"
    00322                    + std::string(m_algType) + "'");
    00323      }
    00324
    00325   return StatusCode::SUCCESS;
    00326 }
    00327 //=============================================================================
    00328 StatusCode EqSolver::finalize   ()
    00329 {
    00330   StatusCode sc = GaudiTool::finalize() ;
    00331
    00332   MsgStream log( msgSvc() , name() ) ;
    00333
    00334   if ( sc.isFailure() )
    00335     {
    00336       return Error("Could not finaliaze base class GaudiTool", sc);
    00337     }
    00338   return StatusCode::SUCCESS;
    00339 }
    00340 //=============================================================================
    00341 EqSolver::~EqSolver( )
    00342 {}
    00343 //=============================================================================