26 #ifndef __AitkenExtrapolation
27 #define __AitkenExtrapolation 1
33 #include <boost/assign/list_of.hpp>
34 #include <boost/assert.hpp>
35 #include <boost/timer.hpp>
36 #include <boost/tuple/tuple.hpp>
37 #include <feel/feeldiscr/functionspace.hpp>
38 #include <feel/feelalg/vector.hpp>
82 template<
typename fs_type >
90 typedef fs_type functionspace_type;
91 typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
93 typedef typename functionspace_type::element_type element_type;
95 typedef typename functionspace_type::template Element<
typename functionspace_type::value_type,
105 typedef std::map<int, convergence_iteration_type> convergence_type;
111 AitkenType _aitkenType = AITKEN_STANDARD,
112 double _failsafeParameter = 1.0,
114 double _minParam = 1e-4 )
117 M_failsafeParameter( _failsafeParameter ),
118 M_previousParameter( _failsafeParameter ),
119 M_maxParameter( _failsafeParameter ),
120 M_minParameter( _minParam ),
121 M_previousResidual( M_Xh,
"previous residual" ),
122 M_previousElement( M_Xh,
"previous element" ),
123 M_currentResidual( M_Xh,
"current residual" ),
124 M_currentElement( M_Xh,
"current element" ),
126 M_aitkenType( _aitkenType ),
128 M_residualConvergence( 1. ),
129 M_hasConverged( false )
139 M_failsafeParameter( tc.M_failsafeParameter ),
140 M_previousParameter( tc.M_previousParameter ),
141 M_maxParameter( tc.M_maxParameter ),
142 M_minParameter( tc.M_minParameter ),
143 M_previousResidual( tc.M_previousResidual ),
144 M_previousElement( tc.M_previousElement ),
145 M_currentResidual( tc.M_currentResidual ),
146 M_currentElement( tc.M_currentElement ),
147 M_cptIteration( tc.M_cptIteration ),
148 M_aitkenType( tc.M_aitkenType ),
149 M_tolerance( tc.M_tolerance ),
150 M_residualConvergence( tc.M_residualConvergence ),
151 M_hasConverged( tc.M_hasConverged )
163 BOOST_PARAMETER_MEMBER_FUNCTION(
171 initializeimpl( residual,currentElt );
177 BOOST_PARAMETER_MEMBER_FUNCTION(
186 ( forceRelaxation, (
bool ),
false )
190 element_type newElt( M_Xh );
191 applyimpl( newElt,residual,currentElt,forceRelaxation );
198 template<
typename eltType >
199 element_type
operator()( element_type
const& residual, eltType
const& elem,
bool _forceRelax=
false )
201 element_type newElt( M_Xh );
202 applyimpl( newElt,residual,elem,_forceRelax );
209 BOOST_PARAMETER_MEMBER_FUNCTION(
219 ( forceRelaxation, (
bool ),
false )
223 applyimpl( newElt,residual,currentElt,forceRelaxation );
258 return M_previousParameter;
266 M_previousParameter=v;
275 return M_cptIteration;
280 return M_hasConverged;
283 double residualNorm()
285 return M_residualConvergence;
296 void forceConvergence(
bool b )
310 void setElement( element_type
const& residual, element_type
const& elem );
315 void setElement( element_type
const& residual, element_range_type
const& elem );
326 template<
typename eltType >
332 return M_convergence;
339 void initializeimpl( element_type
const& residual, element_type
const& elem );
344 void initializeimpl( element_type
const& residual, element_range_type
const& elem );
349 void applyimpl( element_type & new_elem,element_type
const& residual, element_type
const& elem,
bool _forceRelax=
false );
354 void applyimpl( element_range_type new_elem,element_type
const& residual, element_range_type
const& elem,
bool _forceRelax=
false );
372 functionspace_ptrtype M_Xh;
374 double M_failsafeParameter, M_previousParameter, M_maxParameter;
375 double M_minParameter;
377 element_type M_previousResidual, M_previousElement, M_currentResidual, M_currentElement;
380 AitkenType M_aitkenType;
382 double M_residualConvergence;
384 convergence_type M_convergence;
387 mutable boost::timer M_timer;
393 template<
typename fs_type >
395 Aitken<fs_type>::initializeimpl( element_type
const& residual, element_type
const& elem )
397 M_previousResidual = residual;
398 M_previousElement = elem;
403 template<
typename fs_type >
405 Aitken<fs_type>::initializeimpl( element_type
const& residual, element_range_type
const& elem )
407 M_previousResidual = residual;
408 M_previousElement.zero();
409 M_previousElement.add( 1.,elem );
417 template<
typename fs_type >
421 auto oldEltL2Norm = M_previousElement.l2Norm();
423 if ( oldEltL2Norm > 1e-8 )
424 M_residualConvergence = M_currentResidual.l2Norm()/oldEltL2Norm;
427 M_residualConvergence = M_currentResidual.l2Norm();
429 if ( M_residualConvergence < M_tolerance ) M_hasConverged=
true;
435 template<
typename fs_type >
439 M_currentResidual = residual;
440 M_currentElement = elem;
445 template<
typename fs_type >
449 M_currentResidual = residual;
450 M_currentElement.zero();
451 M_currentElement.add( 1.,elem );
459 template<
typename fs_type >
464 setElement( residual,elem );
466 if ( M_cptIteration>=2 )
468 computeResidualNorm();
470 if ( !M_hasConverged || _forceRelax )
472 calculateParameter();
476 if ( !M_hasConverged || _forceRelax )
477 relaxationStep( new_elem );
483 template<
typename fs_type >
485 Aitken<fs_type>::applyimpl( element_range_type new_elem,element_type
const& residual, element_range_type
const& elem,
bool _forceRelax )
488 setElement( residual,elem );
490 if ( M_cptIteration>=2 )
492 computeResidualNorm();
494 if ( !M_hasConverged || _forceRelax )
496 calculateParameter();
500 if ( !M_hasConverged || _forceRelax )
501 relaxationStep( new_elem );
506 template<
typename fs_type >
507 template<
typename eltType >
512 new_elem = M_currentResidual;
514 new_elem.scale( -( 1-M_previousParameter ) );
517 new_elem.add( 1.,M_currentElement );
519 new_elem = M_previousElement;
520 new_elem.add(M_previousParameter,M_currentResidual);
522 M_currentElement = new_elem;
528 template<
typename fs_type >
533 double rel_residual = 1;
535 if ( M_cptIteration > 1 )
536 rel_residual = M_currentResidual.l2Norm()/M_convergence.find( 1 )->second.find(
"absolute_residual" )->second;
539 boost::assign::map_list_of
540 (
"absolute_residual", M_currentResidual.l2Norm() )
541 (
"relative_residual", rel_residual )
542 (
"time", M_timer.elapsed() )
543 (
"relaxation_parameter", M_previousParameter );
544 M_convergence.insert( std::make_pair( M_cptIteration, cit ) );
547 M_previousResidual = M_currentResidual;
548 M_previousElement = M_currentElement;
556 template<
typename fs_type >
567 template<
typename fs_type >
571 if ( !M_convergence.empty() )
573 std::string entry =
"relaxation_parameter";
574 auto it = std::max_element( M_convergence.begin(), M_convergence.end(),
575 [entry] ( std::pair<int, convergence_iteration_type>
const& x,
576 std::pair<int, convergence_iteration_type>
const& y )
578 return x.second.find( entry )->second < y.second.find( entry )->second;
581 if ( it != M_convergence.end() )
582 M_maxParameter = std::max( M_maxParameter,it->second.find( entry )->second );
584 M_previousParameter = M_maxParameter;
588 M_hasConverged=
false;
589 M_convergence.clear();
595 template<
typename fs_type >
599 std::cout <<
"[Aitken] iteration : "<< M_cptIteration
600 <<
" theta=" << M_previousParameter
601 <<
" residualNorm : " << M_residualConvergence
604 LOG(INFO) <<
"[Aitken] iteration : "<< M_cptIteration
605 <<
" theta=" << M_previousParameter
606 <<
" residualNorm : " << M_residualConvergence
610 template<
typename fs_type >
614 std::ofstream ofs( fname.c_str() );
618 ofs.setf( std::ios::scientific );
619 BOOST_FOREACH(
auto cit, M_convergence )
622 ofs << std::setw( 6 ) << cit.first <<
" ";
623 BOOST_FOREACH(
auto it, cit.second )
625 ofs << std::setw( 20 ) << std::setprecision( 16 ) << it.second <<
" ";
633 std::cerr <<
"[Aitken] convergence history filename " << fname <<
" could not be opened"
639 template<
typename fs_type >
643 if ( M_aitkenType==AITKEN_STANDARD )
644 calculateParameter( mpl::int_<AITKEN_STANDARD>() );
646 if ( M_aitkenType==AITKEN_METHOD_1 )
647 calculateParameter( mpl::int_<AITKEN_METHOD_1>() );
649 if ( M_aitkenType==FIXED_RELAXATION_METHOD )
650 M_previousParameter = M_failsafeParameter;
655 template<
typename fs_type >
660 element_type aux( M_Xh,
"aux" );
662 aux = M_currentResidual;
663 aux -= M_previousResidual;
667 aux.scale( 1.0/scalar );
669 element_type aux2( M_Xh,
"aux2" );
671 aux2 = M_currentElement;
672 aux2 -= M_previousElement;
677 scalar = M_failsafeParameter;
680 scalar = M_failsafeParameter;
682 M_previousParameter = 1-scalar;
687 template<
typename fs_type >
692 element_type aux( M_Xh,
"aux" );
694 aux = M_currentResidual;
695 aux -= M_previousResidual;
699 aux.scale( 1.0/scalar );
707 scalar = -M_previousParameter*scalar;
711 scalar = M_failsafeParameter;
713 if ( scalar < M_minParameter )
714 scalar = M_failsafeParameter;
716 M_previousParameter = scalar;
725 template <
typename fs_type >
727 operator++( boost::shared_ptr<Aitken<fs_type> > & aitk )
736 template<
typename SpaceType>
737 boost::shared_ptr<Aitken<SpaceType> >
738 aitkenNew( boost::shared_ptr<SpaceType>
const& _space,
745 boost::shared_ptr<Aitken<SpaceType> > Aitk(
new Aitken<SpaceType>( _space,_type,_init_theta,_tol,_minParam ) );
752 template<
typename Args>
753 struct compute_aitken_return
755 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::space>::type>::type::element_type space_type;
757 typedef Aitken<space_type> type;
758 typedef boost::shared_ptr<type> ptrtype;
762 BOOST_PARAMETER_FUNCTION(
763 (
typename compute_aitken_return<Args>::type ),
767 ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
770 ( type, ( AitkenType ), AITKEN_STANDARD )
771 ( initial_theta, *( boost::is_arithmetic<mpl::_> ), 1.0 )
772 ( min_theta, *( boost::is_arithmetic<mpl::_> ), 1e-4 )
773 ( tolerance, *( boost::is_arithmetic<mpl::_> ), 1.0e-6 )
777 return *aitkenNew( space,type,initial_theta,tolerance,min_theta );
780 BOOST_PARAMETER_FUNCTION(
781 (
typename compute_aitken_return<Args>::ptrtype ),
785 ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
788 ( type, ( AitkenType ), AITKEN_STANDARD )
789 ( initial_theta, *( boost::is_arithmetic<mpl::_> ), 1.0 )
790 ( min_theta, *( boost::is_arithmetic<mpl::_> ), 1e-4 )
791 ( tolerance, *( boost::is_arithmetic<mpl::_> ), 1.0e-6 )
795 return aitkenNew( space,type,initial_theta,tolerance,min_theta );