| # Gaudi Framework, version v21r7p1 |

Home |  | Generated: 15 Feb 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  f = *(jac[i][j]) ;
    00178             gsl_matrix_set      ( df , i , j , f(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           break;
    00228         }
    00229
    00230       status = gsl_multiroot_test_residual (s->f,
    00231                                             m_norm_residual);
    00232
    00233       if ( status != GSL_CONTINUE ) { break; }
    00234     }
    00235
    00236   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00237     {
    00238       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00239     }
    00240
    00241   if (status == GSL_SUCCESS)
    00242     {
    00243       log << MSG::DEBUG
    00244           << "We stopped in the method on the " << iter
    00245           << " iteration (we have maximum "     << m_max_iter
    00246           << " iterations)"                     << endmsg;
    00247     }
    00248   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00249     {
    00250       return Error ( "Method finished with '"
    00251                      + std::string(gsl_strerror(status))
    00252                      +  "' error" );
    00253     }
    00254   else
    00255     {
    00256       return Error ( "Method finished with '" +
    00257                      std::string(gsl_strerror(status))
    00258                      +  "' error" );
    00259     }
    00260
    00261   gsl_multiroot_fdfsolver_free (s);
    00262
    00263   if (status)
    00264     {
    00265       return Error ( "Method finished with '"
    00266                      + std::string(gsl_strerror(status))
    00267                      +  "' error" );
    00268     }
    00269
    00270   return GSL_SUCCESS;
    00271 }
    00272
    00273 //=============================================================================
    00274 StatusCode  EqSolver::initialize()
    00275
    00276 {
    00277   StatusCode sc = GaudiTool::initialize() ;
    00278
    00279   MsgStream log( msgSvc() , name() ) ;
    00280
    00281   if( sc.isFailure() )
    00282     {
    00283       return Error ("Could not initiliaze base class GaudiTool", sc);
    00284     }
    00285   /* The algorithm for multiional root-finding
    00286      (solving nonlinear systems with n equations in n unknowns)*/
    00287
    00288   if(       "fdfsolver_hybridsj" == m_algType )
    00289     {
    00290       m_type = gsl_multiroot_fdfsolver_hybridsj  ;
    00291       log << MSG::DEBUG
    00292           << "Root findind algoritms to be used: "
    00293           << "'gsl_multiroot_fdfsolver_hybridsj'"
    00294           << endmsg;
    00295     }
    00296   else if ( "fdfsolver_hybridj" == m_algType )
    00297     {
    00298       m_type = gsl_multiroot_fdfsolver_hybridj   ;
    00299       log << MSG::DEBUG
    00300           << "Root findind algoritms to be used: "
    00301           << "'gsl_multiroot_fdfsolver_hybridj'"
    00302           << endmsg;
    00303     }
    00304   else if ( "fdfsolver_newton" == m_algType )
    00305     {
    00306       m_type = gsl_multiroot_fdfsolver_newton    ;
    00307       log << MSG::DEBUG
    00308           << "Root findind algoritms to be used: "
    00309           << "'gsl_multiroot_fdfsolver_newton'"
    00310           << endmsg;
    00311     }
    00312   else if ( "fdfsolver_gnewton" == m_algType )
    00313     {
    00314       m_type = gsl_multiroot_fdfsolver_gnewton   ;
    00315       log << MSG::DEBUG
    00316           << "Root findind algoritms to be used: "
    00317           << "'gsl_multiroot_fdfsolver_gnewton'"
    00318           << endmsg;
    00319     }
    00320   else
    00321     {
    00322       return Error(" Unknown algorith type '"
    00323                    + std::string(m_algType) + "'");
    00324      }
    00325
    00326   return StatusCode::SUCCESS;
    00327 }
    00328 //=============================================================================
    00329 StatusCode EqSolver::finalize   ()
    00330 {
    00331   StatusCode sc = GaudiTool::finalize() ;
    00332
    00333   MsgStream log( msgSvc() , name() ) ;
    00334
    00335   if ( sc.isFailure() )
    00336     {
    00337       return Error("Could not finaliaze base class GaudiTool", sc);
    00338     }
    00339   return StatusCode::SUCCESS;
    00340 };
    00341 //=============================================================================
    00342 EqSolver::~EqSolver( )
    00343 {}
    00344 //=============================================================================