strgy_rdslv_parallel.hpp

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: parallel.h,v 1.1 2005/07/11 10:03:56 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_PARALLEL_HPP
#define realroot_SOLVE_SBDSLV_PARALLEL_HPP
//--------------------------------------------------------------------
//#include <realroot/debug_stream.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace realroot
{
  /*                         solveur de reduction parallèle                                 */
  /*                         -------------------------------                                */
  /* améliorations:                                                                         */
  /*   (1) rejet acceleré: subdiviser au niveau des intersections entre bornes              */
  /*   (2) si (1) est implementé on devrait assez vite arriver a la situation ou il n'y a   */
  /*       qu'une borne à traiter => solveur de descartes                                   */
  /* oops: la version utilisée dans la thèse ne permettait pas les réductions multiples !!! */
  /*       je m'en suis rendu compte un peu tard, à une époque je prétendait que certains   */
  /*       systèmes pouvaient être résolus sans subdivision, cela a dû être vrai avec le    */
  /*       solveur de descartes que je n'ai plus utilisé par la suite ...                   */

  namespace rdslv
  {
    template< class system >
    struct parallel :  system::vstack_t
    {
      typedef typename system::creal_t    creal_t;
      typedef typename system::interval_t interval_t;
      typedef typename system::vstack_t   vstack_t;
      typedef typename vstack_t::sz_t sz_t;
      
      /* incrementation et decrementation de la pile */
      sz_t pinc ( sz_t iptr ) 
      {
    assert(iptr%4 == 0);
    sz_t i,s;
    for ( s = i = 0; i < 4*this->nbp(); s += this->pmsz(iptr+i)+this->pMsz(iptr+i), i += 4 ) ;
    s = this->allocc(s+2)+2;
    sz_t inext = this->alloci(this->nbp()*4);
    for ( i = 0; i < 4*this->nbp(); s += this->pmsz(iptr+i)+this->pMsz(iptr+i), i += 4 )
      {
        this->pmsz(inext+i) = this->pmsz(iptr+i);
        this->pmad(inext+i) = s;
        this->pMad(inext+i) = s + this->pmsz(iptr+i);
        this->pMsz(inext+i) = this->pMsz(iptr+i);
      };
    return inext;
      };
      
      sz_t pdec( sz_t iptr )
      { 
    sz_t s,i;
    for ( s = i = 0; i < 4*this->m_nbp; s += this->pmsz(iptr+i)+this->pMsz(iptr+i), i += 4 ) ;
    this->m_cpchnk.dec( s + 2 );
    this->m_ipchnk.dec(4*this->nbp());
    iptr = iptr - 4*this->nbp();
    return iptr;
      }; 
      

      bool process( std::vector<interval_t>& reductions,  const vstack_t& vp, const creal_t& eps = 1e-3, int k = 12  )
      {
    this->copy(vp); /* duplication  de la pile de projections */
    sz_t  iptr, inext, mxsgn, sgn, i, nmult;
    int   itk = 0;
    nmult = iptr = 0;

    for ( ; itk >= 0; )
      {
        sgn = mxsgn = 0;

        /* projections basses */
        for ( mxsgn = i = 0; i < 4*this->nbp(); i += 4 )
          {
        if ( this->pmsz(iptr+i) && !(sgn = vctops::sgnchg( this->pmdata(iptr+i),this->pmsz(iptr+i))) )
            { 
              this->pmsz(iptr+i) = 0; 
              if ( * this->pmdata(iptr+i) > 0.0 ) break;
            };
        if ( mxsgn < sgn ) mxsgn = sgn;
          };
        if ( i < 4*this->nbp() ) { iptr = pdec(iptr); itk --; continue; };
        
        /* projection hautes */
        for ( i = 0; i < 4*this->nbp(); i += 4 )
          {
        if ( this->pMsz(iptr+i) && !(sgn = vctops::sgnchg( this->pMdata(iptr+i), this->pMsz(iptr+i))) )
            {
              this->pMsz(iptr+i) = 0;
              if ( * this->pMdata(iptr+i) < 0.0 ) break;
            };
        if ( mxsgn < sgn ) mxsgn = sgn;
          };
        if ( i < 4*this->nbp() ) { iptr = pdec(iptr); itk --; continue; };
        
        /* si le domaine a été accepté, si la précision ou le quota de multiplication est dépassé. */
        if ( !mxsgn || (this->pupper(iptr)-this->plower(iptr)) < eps )/* || nmult > m_ssz * m_esz ) */
          {
        if ( reductions.size() == 0  || reductions.back().upper() != this->plower(iptr) )
          reductions.push_back(interval_t(this->plower(iptr),this->pupper(iptr)));
        else 
          reductions.back().define(reductions.back().lower(),this->pupper( iptr ));
        iptr = pdec( iptr ); itk --;
        continue;
          };

        /* incrémentation de la pile */
        inext = pinc(iptr);

        //      sz_t nbcstl = 0;
        /* subdivision des projections basses */
        {
          numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
          for ( i = 0; i < 4*this->nbp(); i += 4 )
        if ( this->pmsz(iptr+i) ) 
          {
            nmult += this->pmsz(iptr+i)*(this->pmsz(iptr+i)-1);
            brnops::decasteljau(this->pmdata(iptr+i),this->pmdata(inext+i),this->pmsz(iptr+i));
            //    nbcstl++;
          };
        };
        
        /* subdivision des projections hautes */
        {
          numerics::rounding<creal_t> rnd( numerics::rnd_up() );
          for ( i = 0; i < 4*this->nbp(); i += 4 )
        if ( this->pMsz(iptr+i) ) 
          {
            nmult += this->pMsz(iptr+i)*(this->pMsz(iptr+i)-1);
            brnops::decasteljau(this->pMdata(iptr+i),this->pMdata(inext+i),this->pMsz(iptr+i) );
            //          nbcstl++;
          };
        };
        
        /*      if ( nbcstl == 0 )
          {
        print(iptr);
        }; */
        
        // assert(nbcstl);
        this->plower(inext) =  this->plower(iptr);
        this->pupper(inext) = (this->plower(iptr)+this->pupper(iptr))/2;
        this->plower(iptr)  =  this->pupper(inext);
        iptr = inext;
        itk ++;
      };

    return reductions.size() != 0;
    assert(k>0);
      }

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