synaps/msolve/sbdslv.h

Go to the documentation of this file.
00001 /********************************************************************
00002  *   This file is part of the source code of the SYNAPS library.
00003  *   Author(s): J.P. Pavone, GALAAD, INRIA
00004  *   $Id: sbdslv.h,v 1.1 2005/07/11 10:03:51 jppavone Exp $
00005  ********************************************************************/
00007 #ifndef SYNAPS_SOLVE_SBDSLV_H
00008 #define SYNAPS_SOLVE_SBDSLV_H
00009 /********************************************************************/
00010 #include <set>
00011 #include<vector>
00012 #include <synaps/init.h>
00013 #include <synaps/base/xml.h>
00014 #include <synaps/mpol.h>
00015 #include <synaps/linalg/VectDse.h>
00016 #include <synaps/linalg/MatrDse.h>
00017 #include <synaps/msolve/sbdslv/method.hpp>
00018 #include <synaps/msolve/sbdslv/support.hpp>
00019 #include <synaps/base/chrono.h>
00020 #include <synaps/mpol/MPolDse.h>
00021 //====================================================================
00022 __BEGIN_NAMESPACE_SYNAPS
00023 //====================================================================
00024 namespace sbdslv
00025 {
00026 #define SBDSLV_RD    0
00027 #define SBDSLV_RDL   1
00028 #define SBDSLV_RDS   2
00029 #define SBDSLV_RDRDL 3
00030 #define SBDSLV_RDRDLRDS 4
00031 #define SBDSLV_RDRDS    5
00032 #define SBDSLV_RDLRDS   6
00033 #define SBDSLV_RDSRDL   7
00034 #ifndef SBDSLV_MTH
00035 #define SBDSLV_MTH SBDSLV_RDRDL
00036 #endif
00037   
00038   template< class system, int mth > /* selection de methodes */
00039   struct select_mth
00040   { typedef sbdslv::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00041   template< class system > 
00042   struct select_mth<system,SBDSLV_RD>
00043   { typedef sbdslv::method<system,strgy::simple<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00044   template< class system > 
00045   struct select_mth<system,SBDSLV_RDL>
00046   { typedef sbdslv::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00047   template< class system>
00048   struct select_mth<system,SBDSLV_RDRDL>
00049   { typedef sbdslv::method<system,strgy::newton< system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> 
00050     result_t; };
00051   template< class system >
00052   struct select_mth<system,SBDSLV_RDRDLRDS>
00053   { typedef sbdslv::method<system,strgy::cvmatrix< system, strgy::newton<system,strgy::simple<system> > >,rdslv::parallel,sbdrl::parametric> result_t; };
00054   template< class system > 
00055     struct select_mth<system,SBDSLV_RDS>
00056   { typedef sbdslv::method<system,strgy::cvmatrix<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00057 
00058   template< class system > 
00059     struct select_mth<system,SBDSLV_RDRDS>
00060   { typedef sbdslv::method<system,strgy::cvmatrix<system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00061   
00062   template< class system > 
00063     struct select_mth<system,SBDSLV_RDLRDS>
00064   { typedef sbdslv::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00065 
00066   template< class system > 
00067     struct select_mth<system,SBDSLV_RDSRDL>
00068   { typedef sbdslv::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00069   
00070 };
00071 
00083 template< class T, int mth = SBDSLV_RDRDL>
00084 struct SBDSLV
00085 {
00086   typedef T coeff_t;
00087   T prec;
00089   SBDSLV( const T& precision = (T)1e-6 ): prec(precision){};
00090 };
00091 
00092 namespace SolveSBDSLV
00093 {
00094   template< class T, class eenv_t, class MPOL >
00095   void fill_data( T * data, eenv_t * env, const MPOL& p, binomials<T>& benv  )
00096   {
00097     typename MPOL::const_iterator mn; int i,a;
00098     std::fill(data,data+env->sz(),0); 
00099     mpolfill(data,p,*env);
00100     convertm2b (data,*env,benv);
00101   };
00102   template< class monom >
00103   void scan_monom_env ( std::set< int >& env, const monom& m )
00104   { for ( unsigned i = 0; i < (m.rep()).size(); i++ ) if ( (m.rep())[i] != 0 ) env.insert(i); };
00105 }
00106 
00107 
00108 
00115 template< class MPC, class C, int MTH, class Domains >
00116 Seq<VectDse<C> >  solve( const MPC& l, const SBDSLV< C, MTH >& mth, 
00117                          const Domains & dom, 
00118                          bool verbose = false  )
00119 {
00120   typedef C T;
00121   typedef typename MPC::value_type          MPoly;
00122   typedef MPoly                             mpoly_t;
00123   typedef MPC                               plist_t;
00124   typedef typename plist_t::const_iterator  pptr_t;
00125   typedef typename MPoly::const_iterator    mptr_t;
00126   typedef typename MPoly::monom_t           monom_t;
00127   typedef SBDSLV< C, MTH >         mth_t;
00128   typedef int            sz_t;
00129   sz_t nvars, neqs, * nvrs, ** vrs, ** szs; 
00130   pptr_t cpl;
00131   int e;
00132   nvars = 0;
00133   // recuperation de l'environement global
00134   for ( neqs = 0, cpl = l.begin(); cpl != l.end(); cpl++, neqs++  )
00135     for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ )
00136       nvars = std::max(nvars,mn->nvars()+1);
00137   // recuperation de l'environement de chaque equation
00138   std::set<int> * seenvs = new std::set< int > [ neqs ];
00139   for ( e = 0, cpl = l.begin(); e < neqs; e++, cpl++  )
00140     for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ ) 
00141       SolveSBDSLV::scan_monom_env(seenvs[e],*mn);
00142   
00143   //  assert(sizeof(sz_t)==sizeof(sz_t*));
00144   
00145   sz_t  s = 0;
00146   for ( e = 0; e < neqs; s += seenvs[e].size(), e++ );
00147   nvrs = new sz_t [neqs+2*s];
00148   vrs  = new sz_t*[2*neqs];
00149   szs  = vrs + neqs;
00150   sz_t * base = nvrs + neqs;
00151   std::set<int>::iterator seenvit;  
00152   for ( e = 0, cpl = l.begin(); e < neqs; base += 2*nvrs[e++], cpl++ )
00153     {
00154       nvrs[e] = seenvs[e].size();
00155       vrs [e] = base;
00156       szs [e] = base + nvrs[e];
00157       int k = 0;
00158       for ( seenvit = seenvs[e].begin(); k < nvrs[e]; seenvit++, k++  ) 
00159         {
00160           vrs[e][k]  = *seenvit;
00161           szs[e][k]  = MPOLDST::degree( *cpl, *seenvit ) + 1;
00162         };
00163     };
00164   
00165   std::vector<C> result;
00166   C * dmns = 0;
00167   if ( dom.size() >= 2*nvars )
00168     {
00169       using namespace let;
00170       dmns = new C[ 2*nvars ];
00171       for ( int i = 0; i < 2*nvars; i ++ ) assign(dmns[i],dom[i]);
00172     };
00173   
00174   chrono tm;
00175   std::vector<T> tmp;
00176   
00177   {  /* remplissage et appel */
00178     using namespace sbdslv;
00179     binomials<T> bznv;
00180     typedef  sbdslv::system<T> system_t;
00181     system_t sys( neqs, nvars, nvrs, (const int**)vrs, (const int**)szs, (T)(mth.prec) );
00182     for ( e = 0, cpl = l.begin(); e < neqs; e ++, cpl++  )
00183       SolveSBDSLV::fill_data( sys.data(e), sys.env(e), *cpl, bznv );
00184     sbdslv::system_ctrl< std::vector< T > >  intf(tmp);
00185     tm.start();
00186     typename sbdslv::select_mth<system_t,MTH>::result_t __mth;
00187     __mth.launch(&sys,intf,dmns);
00188     tm.stop();
00189     if ( verbose )
00190       {
00191         std::cout << " iterations  = " << __mth.m_niter << std::endl;
00192         std::cout << "subdivisions = "<< __mth.m_nsbd  << std::endl;
00193       };
00194   };
00195   
00196   if ( !verbose && dmns ) delete[] dmns;
00197   delete[] seenvs;
00198   delete[] nvrs;
00199   delete[] vrs;
00200   if ( verbose )
00201     {
00202       std::cout << "<domain>";
00203       if  ( dmns )
00204         {
00205           std::cout << std::endl;
00206           for ( int i = 0; i < nvars; i ++ )
00207             std::cout << "\t" 
00208                       << "x" << i//monom_t::var_of_index(i) 
00209                       << " = [ " << dmns[2*i] 
00210                       << "," 
00211                       << dmns[2*i+1] 
00212                       << " ]";
00213         }
00214       else std::cout << "[0,1]^" << nvars;
00215       std::cout << "</domain>\n";
00216       //        std::cout << "\treductions   = " << bsys.nrest() << std::endl;
00217       //        std::cout << "\titerations   = " << bsys.niter() << std::endl;
00218       std::cout << "<solver_time>"; tm.write(std::cout); 
00219       std::cout <<"</solver_time>"<< std::endl; 
00220       if ( dmns ) delete[] dmns;
00221     };
00222   
00223   int nsols = tmp.size() / (2*nvars);
00224   T boxes[nsols*2*nvars];
00225   std::copy(&tmp[0], &tmp[0]+2*nvars*nsols,boxes);
00226   int nsols2 = sbdslv::clean_result( &tmp[0], nvars, nsols,(T)1e-3 );
00227   if ( verbose ) { 
00228     xml::print(std::cout,nsols,"nb_box")<< std::endl;
00229     for ( unsigned i = 0; i < nsols; i ++ )
00230       {
00231         std::cout << "\t[ ";
00232         for ( unsigned j = 0; j < nvars; j ++ )
00233           std::cout << (boxes[2*(i*nvars+j)]+boxes[2*(i*nvars+j)+1])/2.0 << " ";
00234         std::cout << "]";
00235         std::cout << std::endl;
00236       };
00237     xml::print(std::cout,nsols2,"nb_sol")<< std::endl; 
00238     //    if ( nsols ) vctops::print(&tmp[0],nvars*nsols );
00239   };
00240 
00241 
00242   tmp.resize( nsols2 * nvars );
00243   //  cout <<"-----"<<endl;
00244   VectDse<T> vsol(nvars, AsSize());
00245   Seq<VectDse<T> > res;
00246   for(unsigned i=0;i<nsols2;++i)
00247     {
00248       for(unsigned j=0; j<nvars; j++)
00249         vsol[j]=tmp[i*nvars+j];
00250       res.push_back(vsol);
00251     }
00252   return res;
00253 };
00254 
00255 //----------------------------------------------------------------------
00256 template< class MPC, class C, int MTH>
00257 Seq<VectDse<C> >  solve( const MPC& l, const SBDSLV< C, MTH >& mth, 
00258                          bool verbose = false  )
00259 {
00260   return solve(l,mth,VectDse<C>(),verbose);
00261 }
00262 //----------------------------------------------------------------------
00263 
00264 template<class C, int METHOD>
00265 struct SbdSolver
00266 {
00267   typedef sbdslv::system<C> system_t;
00268   system_t  * sys;
00269   
00270 
00271   
00272   template< class MPDI >
00273   SbdSolver( MPDI  first, MPDI last, const C& eps )
00274   {
00275 
00276     int    ne = 0;
00277     MPDI eptr;
00278     for ( eptr = first; eptr != last; eptr ++ ) ne++;
00279     int nvrs[ne];
00280     const int *szs[ne];
00281     int e;
00282     int mxvr = 0;
00283     std::set<int> gvrs;
00284     for ( e = 0, eptr = first; eptr != last; eptr ++, e++  ) { 
00285       nvrs[e] = eptr->nvr(); 
00286       szs[e]  = eptr->szs();
00287       const int * tmp = eptr->vrs();
00288       for ( int v = 0; v < eptr->nvr(); v ++ ) gvrs.insert(tmp[v]);
00289     };
00290     int mxv = *(gvrs.rend());
00291     int lnames[ mxv+1 ];
00292     int nv = 0;
00293     for ( std::set<int>::const_iterator it = gvrs.begin(); it != gvrs.end(); lnames[*it++] = nv++ );
00294     int * vrdat = new int[ nv * ne ];
00295     int * vrptr = vrdat;
00296     const int *  vrs[ne];
00297     e = 0;
00298     for ( eptr = first; eptr != last; vrptr+= eptr->nvr(), e++, eptr ++ ) {
00299       vrs[e] = vrptr;
00300       const int * tmp = eptr->vrs();
00301       for ( int v = 0; v < eptr->nvr();  v ++ )
00302         vrptr[v] = lnames[ tmp[v] ];    
00303     };
00304     //    std::cout << ne << " " << nv << std::endl;
00305     assert(ne>=nv);
00306     sys = new sbdslv::system<C>( ne, nv, nvrs, (const int**)vrs, (const int**)szs, eps );
00307     delete[] vrdat;
00308   };
00309   
00310   void add( const mpoldse::bernstein<C>& data, int e )
00311   {
00312     sys->set_bernstein2( data.begin(), e );
00313   };
00314   
00315   void add( const mpoldse::monomial<C>& data, int e ) 
00316   {
00317     sys->set_monomial2( data.begin(), e );
00318   };
00319   
00320   template<class MPDI>
00321   void solve( std::vector<C> & R, int& ns, int& nv, MPDI first, MPDI last )
00322   {
00323     R.resize(0);
00324     typedef SBDSLV< C, METHOD > MTH;
00325     typename sbdslv::select_mth<system_t,METHOD>::result_t __mth;
00326     int e = 0;
00327     for ( MPDI eptr = first; eptr != last; add( eptr->rep(), e++  ), eptr++ );
00328     //std::vector<C> tmp;
00329     sbdslv::system_ctrl< std::vector< C > >  intf(R);
00330     int nvars = sys->nvr();
00331     __mth.launch(sys,intf,0);
00332 
00333     int nsols = R.size() / (2*nvars);
00334     int nsols2 = sbdslv::clean_result( &R[0], nvars, nsols,(C)1e-3 );
00335     R.resize(nsols2) ;
00336     ns = nsols2;
00337     nv = nvars;
00338   };
00339   
00340   ~SbdSolver() { delete sys; };
00341 };
00342 
00343 
00344 //--------------------------------------------------------------------
00345 __END_NAMESPACE_SYNAPS
00346 /********************************************************************/
00347 #endif //SYNAPS_SOLVE_SBDSLV_H

SYNAPS DOCUMENTATION
logo