/Users/mourrain/Devel/mmx/realroot/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 += pmsz(iptr+i)+pMsz(iptr+i), i += 4 ) ;
        s = allocc(s+2)+2;
        sz_t inext = alloci(this->nbp()*4);
        for ( i = 0; i < 4*this->nbp(); s += pmsz(iptr+i)+pMsz(iptr+i), i += 4 )
          {
            pmsz(inext+i) = pmsz(iptr+i);
            pmad(inext+i) = s;
            pMad(inext+i) = s + pmsz(iptr+i);
            pMsz(inext+i) = 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 += pmsz(iptr+i)+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 ( pmsz(iptr+i) && !(sgn = vctops::sgnchg(pmdata(iptr+i),pmsz(iptr+i))) )
                    { 
                      pmsz(iptr+i) = 0; 
                      if ( *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 ( pMsz(iptr+i) && !(sgn = vctops::sgnchg( pMdata(iptr+i), pMsz(iptr+i))) )
                    {
                      pMsz(iptr+i) = 0;
                      if ( *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 || (pupper(iptr)-plower(iptr)) < eps )/* || nmult > m_ssz * m_esz ) */
              {
                if ( reductions.size() == 0  || reductions.back().upper() != plower(iptr) )
                  reductions.push_back(interval_t(plower(iptr),pupper(iptr)));
                else 
                  reductions.back().define(reductions.back().lower(),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 ( pmsz(iptr+i) ) 
                  {
                    nmult += pmsz(iptr+i)*(pmsz(iptr+i)-1);
                    brnops::decasteljau(pmdata(iptr+i),pmdata(inext+i),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 ( pMsz(iptr+i) ) 
                  {
                    nmult += pMsz(iptr+i)*(pMsz(iptr+i)-1);
                    brnops::decasteljau(pMdata(iptr+i),pMdata(inext+i),pMsz(iptr+i) );
                    //              nbcstl++;
                  };
            };
            
            /*      if ( nbcstl == 0 )
              {
                print(iptr);
                }; */
            
            // assert(nbcstl);
            plower(inext) =  plower(iptr);
            pupper(inext) = (plower(iptr)+pupper(iptr))/2;
            plower(iptr)  =  pupper(inext);
            iptr = inext;
            itk ++;
          };

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

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