include/realroot/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 //

---