| # Gaudi Framework, version v20r4
Generated: 8 Jan 2009

# 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   : GaudiTool       ( 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   declareInterface <IEqSolver>(this);
    00088   declareProperty ( "Algorithm", m_algType       );
    00090   declareProperty ( "Iteration", m_max_iter      );
    00092   declareProperty ( "Residual" , m_norm_residual );
    00093 }
    00094
    00095 namespace
    00096 {
    00097   using namespace Genfun;
    00098
    00101   int fun_gsl  ( const gsl_vector* v      ,
    00102                  void*             params ,
    00103                  gsl_vector*       f      )
    00104   {
    00105     EqSolver::EqSolverMisc* local =
    00106       static_cast <EqSolver::EqSolverMisc*> (params);
    00107     const EqSolver::Equations& eqs = *(local->equations()) ;
    00108     Argument&                  arg = local->argument()  ;
    00109
    00110     {for ( unsigned int i = 0; i < v->size; ++i)
    00111       {
    00112         arg[i] = gsl_vector_get (v, i);
    00113       }
    00114     }
    00115
    00116     {for ( unsigned int i = 0; i < v->size; ++i)
    00117       {
    00118         gsl_vector_set(f, i, (*eqs[i])(arg));
    00119       }
    00120     }
    00121     return GSL_SUCCESS;
    00122   }
    00123
    00124
    00127   int dfun_gsl  ( const gsl_vector* v      ,
    00128                   void*             params ,
    00129                   gsl_matrix*       df     )
    00130   {
    00131     EqSolver::EqSolverMisc* local =
    00132       static_cast <EqSolver::EqSolverMisc*> (params);
    00133     const EqSolver::Jacobi&     jac = local -> jacobi()    ;
    00134     Argument&                   arg = local -> argument()  ;
    00135
    00136     {for ( unsigned int i = 0; i < v->size; ++i)
    00137       {
    00138         arg[i] = gsl_vector_get (v, i);
    00139       }
    00140     }
    00141
    00142     {for( unsigned int i = 0 ; i < v->size ; ++i )
    00143       {
    00144         for (unsigned int j = 0; j < v->size; ++j)
    00145           {
    00146             Genfun::GENFUNCTION  f = *(jac[i][j]) ;
    00147             gsl_matrix_set      ( df , i , j , f ( arg ) ) ;
    00148           }
    00149       }
    00150     }
    00151     return GSL_SUCCESS;
    00152   }
    00153
    00157   int fdfun_gsl ( const gsl_vector* v ,
    00158                   void*       params  ,
    00159                   gsl_vector* f       ,
    00160                   gsl_matrix* df      )
    00161   {
    00162     EqSolver::EqSolverMisc* local =
    00163       static_cast <EqSolver::EqSolverMisc*> (params);
    00164     const EqSolver::Equations&  eqs = *(local->equations()) ;
    00165     const EqSolver::Jacobi&     jac = local->jacobi()    ;
    00166     Argument&                   arg = local->argument()  ;
    00167
    00168     {for ( unsigned int i = 0; i < v->size; ++i)
    00169       {
    00170         arg[i] = gsl_vector_get (v, i);
    00171       }
    00172     }
    00173
    00174     {for( unsigned int i = 0 ; i < v->size ; ++i )
    00175       {
    00176         gsl_vector_set(f, i, (*eqs[i])(arg ) ) ;
    00177         for (unsigned int j = 0; j < v->size; ++j)
    00178           {
    00179             Genfun::GENFUNCTION  f = *(jac[i][j]) ;
    00180             gsl_matrix_set      ( df , i , j , f(arg) ) ;
    00181           }
    00182       }
    00183     }
    00184     return GSL_SUCCESS;
    00185   }
    00186 }
    00187
    00188 //=============================================================================
    00193 StatusCode EqSolver::solver (const Equations&  funcs ,
    00194                              Arg&              arg   ) const
    00195 //=============================================================================
    00196 {
    00197   using namespace Genfun;
    00198
    00199   gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
    00200                                                  arg.dimension() );
    00201   MsgStream log( msgSvc(), name() );
    00202
    00203   EqSolverMisc local (funcs, arg);
    00204
    00205   const gsl_multiroot_fdfsolver_type *T = m_type;
    00206   gsl_multiroot_fdfsolver *s;
    00207   int    status = 0 ;
    00208   size_t iter   = 0 ;
    00209
    00210   gsl_multiroot_function_fdf function;
    00211
    00212   function.f = &fun_gsl;
    00213   function.df = &dfun_gsl;
    00214   function.fdf = &fdfun_gsl;
    00215   function.n = vect.vector.size;
    00216   function.params = (void*) &local;
    00217
    00218   s = gsl_multiroot_fdfsolver_alloc(T, vect.vector.size);
    00219   gsl_multiroot_fdfsolver_set (s, &function, &vect.vector);
    00220
    00221   for (iter = 0; iter < m_max_iter; ++iter)
    00222     {
    00223       status = gsl_multiroot_fdfsolver_iterate (s);
    00224       if (status)
    00225         {
    00226           return Error
    00227             ("Error from gsl_gsl_multiroot_fdfsolver_iterate '"
    00228              + std::string(gsl_strerror(status)) + "'") ;
    00229           break;
    00230         }
    00231
    00232       status = gsl_multiroot_test_residual (s->f,
    00233                                             m_norm_residual);
    00234
    00235       if ( status != GSL_CONTINUE ) { break; }
    00236     }
    00237
    00238   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00239     {
    00240       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00241     }
    00242
    00243   if (status == GSL_SUCCESS)
    00244     {
    00245       log << MSG::DEBUG
    00246           << "We stopped in the method on the " << iter
    00247           << " iteration (we have maximum "     << m_max_iter
    00248           << " iterations)"                     << endreq;
    00249     }
    00250   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00251     {
    00252       return Error ( "Method finished with '"
    00253                      + std::string(gsl_strerror(status))
    00254                      +  "' error" );
    00255     }
    00256   else
    00257     {
    00258       return Error ( "Method finished with '" +
    00259                      std::string(gsl_strerror(status))
    00260                      +  "' error" );
    00261     }
    00262
    00263   gsl_multiroot_fdfsolver_free (s);
    00264
    00265   if (status)
    00266     {
    00267       return Error ( "Method finished with '"
    00268                      + std::string(gsl_strerror(status))
    00269                      +  "' error" );
    00270     }
    00271
    00272   return GSL_SUCCESS;
    00273 }
    00274
    00275 //=============================================================================
    00276 StatusCode  EqSolver::initialize()
    00277
    00278 {
    00279   StatusCode sc = GaudiTool::initialize() ;
    00280
    00281   MsgStream log( msgSvc() , name() ) ;
    00282
    00283   if( sc.isFailure() )
    00284     {
    00285       return Error ("Could not initiliaze base class GaudiTool", sc);
    00286     }
    00287   /* The algorithm for multiional root-finding
    00288      (solving nonlinear systems with n equations in n unknowns)*/
    00289
    00290   if(       "fdfsolver_hybridsj" == m_algType )
    00291     {
    00292       m_type = gsl_multiroot_fdfsolver_hybridsj  ;
    00293       log << MSG::DEBUG
    00294           << "Root findind algoritms to be used: "
    00295           << "'gsl_multiroot_fdfsolver_hybridsj'"
    00296           << endreq;
    00297     }
    00298   else if ( "fdfsolver_hybridj" == m_algType )
    00299     {
    00300       m_type = gsl_multiroot_fdfsolver_hybridj   ;
    00301       log << MSG::DEBUG
    00302           << "Root findind algoritms to be used: "
    00303           << "'gsl_multiroot_fdfsolver_hybridj'"
    00304           << endreq;
    00305     }
    00306   else if ( "fdfsolver_newton" == m_algType )
    00307     {
    00308       m_type = gsl_multiroot_fdfsolver_newton    ;
    00309       log << MSG::DEBUG
    00310           << "Root findind algoritms to be used: "
    00311           << "'gsl_multiroot_fdfsolver_newton'"
    00312           << endreq;
    00313     }
    00314   else if ( "fdfsolver_gnewton" == m_algType )
    00315     {
    00316       m_type = gsl_multiroot_fdfsolver_gnewton   ;
    00317       log << MSG::DEBUG
    00318           << "Root findind algoritms to be used: "
    00319           << "'gsl_multiroot_fdfsolver_gnewton'"
    00320           << endreq;
    00321     }
    00322   else
    00323     {
    00324       return Error(" Unknown algorith type '"
    00325                    + std::string(m_algType) + "'");
    00326      }
    00327
    00328   return StatusCode::SUCCESS;
    00329 }
    00330 //=============================================================================
    00331 StatusCode EqSolver::finalize   ()
    00332 {
    00333   StatusCode sc = GaudiTool::finalize() ;
    00334
    00335   MsgStream log( msgSvc() , name() ) ;
    00336
    00337   if ( sc.isFailure() )
    00338     {
    00339       return Error("Could not finaliaze base class GaudiTool", sc);
    00340     }
    00341   return StatusCode::SUCCESS;
    00342 };
    00343 //=============================================================================
    00344 EqSolver::~EqSolver( )
    00345 {}
    00346 //=============================================================================