strgy_newton.h

Go to the documentation of this file.

/**********************************************************************
 *    This file is part of the source code of the realroot library.     
 *    Author(s): J.P. Pavone, GALAAD, INRIA                           
 *    $Id: newton.h,v 1.1 2005/07/11 10:03:57 jppavone Exp $
 **********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_STRGY_NEWTON_HPP
#define realroot_SOLVE_SBDSLV_STRGY_NEWTON_HPP
//--------------------------------------------------------------------
#include <realroot/linear_doolittle.hpp>
#include <realroot/loops_vctops.hpp>
#include <realroot/strgy_empty.h>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace realroot
{
  namespace strgy
  {
    template< class system, class _previous_ = empty<system> >
    struct newton : _previous_
    {
      typedef typename system::creal_t creal_t;
      typedef typename system::sz_t    sz_t;
      typedef typename system::interval_t interval_t;
      
      bool process( system * slv )
      {
    /* application de la strategie de la classe de base */
    if ( !_previous_::process(slv) ) return false;
    /* estimation de la jacobienne en utilisant une différence finie sur le polygone  */
    /* de contrôle, la position correspond à l'abscisse de greville la plus proche    */
    /* du centre du domaine                                                           */
    sz_t   nvars    = slv->nvr();
    sz_t   neqs     = slv->neqs();
    sz_t * selected = slv->m_eqssel;
    assert(nvars<=neqs);
    slv->restrict(selected,neqs); /*std::min(neqs,nvars));*/
    slv->elevate(selected,neqs);  /*selected,std::min(neqs,nvars));*/
    creal_t jac[nvars*neqs];
    sz_t add = 0;
    sz_t * str = slv->m_elenv.str();
    sz_t * szs = slv->m_elenv.szs();
    for ( sz_t v = 0; v < nvars; v ++ ) 
      add += (szs[v]-1)/2*str[v];
    for ( sz_t s = 0; s < nvars; s ++ ) 
      {
        sz_t e = selected[s];
        creal_t * src = slv->elldata(e);
        /* calcul du coefficient               */
        for ( sz_t v = 0; v < nvars; v ++ ) 
        jac[s*nvars+v] = (src[add+str[v]]-src[add])/(slv->top()[v].size()*szs[v]);
      };
    
    //  vctops::print(jac,(unsigned)nvars*nvars);
    //  std::cout << std::endl;
    //  vctops::scale(jac,(unsigned)nvars*nvars);
    //  vctops::print(jac,(unsigned)nvars*nvars);
    //  std::cout << std::endl;
    /* inversion de la jacobienne */
    int c = 0;
    creal_t ijac[nvars*nvars];
    /* tant que la matrice est singuliere */
    while( !linear::LUinverse(ijac, jac, nvars) ) 
      {
        /* on bruite ses coefficients */
        //    for ( sz_t v = 0; v < nvars; v ++ )
        for ( sz_t e = 0; e < nvars; e ++ )
          jac[e*nvars+e] += 1;
        c ++ ;
      };
    //  std::cout << std::endl;
    //  std::cout << " try =  " << c << std::endl;
    //  for ( unsigned i = 0; i < nvars; i ++ )
    //    {
    //      vctops::print(ijac + i*nvars, nvars );
    //      std::cout << std::endl;
        //    };
    //  std::cout << std::endl;

    slv->transform( ijac,selected, nvars );
    slv->project(selected,nvars);
    /* faire tourner la selection (devrait être cacher dans select(n) du solver */
    sz_t sv = selected[0];
    for ( sz_t e = 0;  e<  neqs-1; e ++ )
    selected[e] = selected[e+1];
    selected[neqs-1] = sv;
    return true;
      };
    };
  };
};

//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif