mmx::realroot Namespace Reference

Namespaces

Classes

Functions

Variables

Function Documentation

unsigned mmx::realroot::clean_result ( real_t *  isols,
int  nvars,
int  nsols,
const real_t &  prec 
) [inline]

Definition at line 34 of file system_support.hpp.

References mmx::abs(), mmx::sparse::copy(), and mmx::max().

Referenced by solver< C, ProjRd< MTH > >::run(), and solver< C, ProjRd< MTH > >::solve_monomial().

00035   {
00036     int a,c,n,i;
00037     real_t cd;
00038     if ( nsols == 0 ) return 0;
00039     for ( a = 0; a < nsols*nvars; isols[a] = (isols[2*a]+isols[2*a+1])/2.0, a ++ ) ;
00040 
00041     for ( n = 1, c = 0, a = 0; a < nsols*nvars; a += nvars )
00042       {
00043         for ( cd = 0, i = 0; i < nvars; cd = std::max(cd,std::abs((isols+c)[i]-(isols+a)[i])), i ++ ) ;
00044         if ( cd < prec )
00045         {
00046           double fa = ((double)n)/(n+1);
00047           double fb = ((double)1.0)/(n+1);
00048           for ( int i = 0; i < nvars; i ++ ) (isols+c)[i] = fa*(isols+c)[i]+fb*(isols+a)[i];
00049           n ++;
00050         }
00051         else
00052           {
00053             c += nvars;
00054             std::copy(isols+a,isols+a+nvars,isols+c);
00055             n  = 1;
00056           };
00057       };
00058     return (c+nvars)/nvars;
00059   };

T mmx::realroot::det ( T *  mat,
int  order 
) [inline]

Definition at line 142 of file solver_mv_monomial_tests.hpp.

References submatrix().

Referenced by exclude2(), include1(), include2(), and signof_det().

00143     {
00144       if(order == 1)
00145         return mat[0]; // Matrix is of order one
00146       int i;
00147       T ev(0);
00148 
00149       for(i = 0; i < order; i++)
00150       { 
00151         //std::cout<<"entering"<<order<<",i="<<i <<std::endl;
00152         ev +=  T(i%2==0?1:-1)* mat[i*order] *  det(submatrix(mat, order, i, 0), order - 1);
00153 
00154         //std::cout<<"OK, ("<<order<<" "<<i<<")"<<std::endl;
00155       }
00156       return ev;
00157     }

T* mmx::realroot::eval_poly_matrix ( Seq< T >  val,
polynomial< double, with< MonomialTensor > > *  mat,
int  n 
) [inline]

Definition at line 90 of file solver_mv_monomial_tests.hpp.

References mmx::eval(), mmx::rep(), and Seq< C, R >::size().

Referenced by exclude2(), include1(), include2(), and include3().

00091     {
00092       int m(val.size());
00093       T *ev;//[m*n];
00094       ev= new T [m*n];
00095 
00096       for( int i=0;i<m;i++)
00097       {
00098         for( int k=0;k<n;k++)
00099         {
00100           tensor::eval( ev[i*m+k] , mat[i*m+k].rep(),
00101                         val , m  );
00102         }
00103       }
00104       return ev;
00105     }

bool mmx::realroot::exclude1 ( box_rep< POL > *  b,
Seq< POL * > &  S0 
) [inline]

Exclusion criteria (inteval arithmetic).

Definition at line 350 of file solver_mv_monomial_tests.hpp.

References mmx::assign(), DPOL, mmx::eval(), Interval< T, r >::M, and Interval< T, r >::m.

00351     {
00352       Interval<double> ev;
00353       Seq<Interval<double> > dom;
00354       dom= b->template domain<double>();
00355       DPOL tmp;
00356 
00357       for (unsigned i=0; i!=b->nbpol(); ++i)
00358       {
00359         let::assign(tmp, *S0[i] );
00360         tensor::eval( ev , tmp.rep(), 
00361                       dom , b->nbvar() );       
00362         if ( ev.m*ev.M > 0 )
00363         {   
00364           //   std::cout<<i<<" ev"<< 
00365           //    dom<<"\nf= "<< *(S0[1])<<
00366           //   "\n f([])= " <<ev<<std::endl;
00367           //  if (td!=0)
00368           // std::cout<<"!!!!!!----td: "<<td <<std::endl;
00369           return true;
00370         }       
00371         
00372       }
00373       return false;
00374     };

bool mmx::realroot::exclude2 ( box_rep< POL > *  b,
polynomial< double, with< MonomialTensor > > *  J 
) [inline]

Exclusion criteria (topological degree).

Definition at line 379 of file solver_mv_monomial_tests.hpp.

References det(), eval_poly_matrix(), Interval< T, r >::M, and Interval< T, r >::m.

Referenced by solver_mv_monomial< FT, POL >::solve_system().

00380   {
00381     Interval<double> ev;
00382     Interval<double> * m;
00383     
00384     //Seq<POL> SS= b->system();
00385     unsigned dim= b->nbvar();
00386     
00387     // Check Jacobian's sign
00388     m= eval_poly_matrix( b->template domain<double>() ,J,dim);
00389     ev= det(m,dim);
00390     if ( ev.m*ev.M < 0 ) 
00391       return false;
00392 
00393     int td= topological_degree_2d<POL, double>(b);
00394     
00395     return ( td==0 );
00396   }

bool mmx::realroot::exclude3 ( box_rep< POL > *  b  )  [inline]

Definition at line 399 of file solver_mv_monomial_tests.hpp.

References no_variation().

00400     {
00401       Seq<POL> S= b->system();
00402       for (unsigned i=0; i!=b->nbpol(); ++i)
00403         if ( no_variation( S[i]) ) 
00404           return true;
00405       return false;
00406     }

void mmx::realroot::fill_data ( T *  data,
eenv_t *  env,
const MPOL &  p,
binomials< T > &  benv 
) [inline]

Definition at line 120 of file solver_mv_bernstein.hpp.

References mmx::tensor::convertm2b(), mmx::vct::fill(), and mmx::tensor::mpolfill().

Referenced by solver< C, ProjRd< MTH > >::solve_monomial().

00121   {
00122     std::fill(data,data+env->sz(),0); 
00123     mpolfill(data,p.rep(),*env);
00124     convertm2b (data,*env,benv);
00125   };

bool mmx::realroot::include1 ( box_rep< POL > *  b,
polynomial< double, with< MonomialTensor > > *  J 
) [inline]

Inclusion criteria (Miranda Test).

Definition at line 209 of file solver_mv_monomial_tests.hpp.

References det(), eval_poly_matrix(), Interval< T, r >::M, Interval< T, r >::m, miranda_test(), and POL.

00210     {
00211       Interval<double> ev(0);
00212       Interval<double> * m;
00213       unsigned i,j,c; //,r ;
00214       POL u,l;
00215       Seq<POL> SS= b->system();
00216       unsigned dim= b->nbvar();
00217 
00218       bool t[dim][dim];
00219 
00220       m= eval_poly_matrix( b->template domain<double>() ,J,dim);
00221 
00222 // std::cout<<"Ok " <<std::endl;     
00223 // double dm[9];
00224 // for(i = 0; i < 9; i++)
00225 // {  dm[i]= m[i].M;
00226 //   std::cout<<"dm["<<i<<"]="<<dm[i] <<std::endl;
00227 // }
00228 // ev= det(dm,dim);
00229 // std::cout<<"Ok ="<<ev <<std::endl;     
00230 
00231 // for(i = 0; i < 9; i++)
00232 //   std::cout<<"m["<<i<<"]="<<m[i] <<std::endl;
00233 
00234  //   std::cout<<"SUB "<<std::endl;
00235  // ev= det(submatrix(m, 3, 2, 0), 2);
00236  //   std::cout<<"FULL "<<std::endl;
00237 
00238       ev=det(m,dim);
00239 
00240       if ( ev.m*ev.M < 0 ) 
00241         return false;
00242       
00243       for (i=0; i!=dim;++i)
00244         for (j=0; j!=dim;++j){
00245           t[i][j]= miranda_test(SS,i,j);
00246           //std::cout<<"mir("<<i<<","<<j<<")="<< t[i][j] <<std::endl;
00247         }
00248 
00249       c=0;
00250       for (i=0; i!=dim;++i)
00251         for (j=0; j!=dim;++j)
00252           if (t[i][j]==true) 
00253           {     
00254             c++; break;
00255           }
00256       if (c<dim) return false;
00257       
00258       c=0;
00259       for (i=0; i!=dim;++i)
00260         for (j=0; j!=dim;++j)
00261           if (t[j][i]==true) 
00262           {     
00263             c++; break;
00264           }
00265       if (c<dim) return false;
00266       
00267       //std::cout<<"INCLUDE. ev="<<ev<<", c="<<c <<"*dim:"<<dim <<std::endl;
00268       return true;
00269     };

bool mmx::realroot::include2 ( box_rep< POL > *  b,
polynomial< double, with< MonomialTensor > > *  J 
) [inline]

Inclusion criteria (Jacobian+Topological Degree).

Definition at line 273 of file solver_mv_monomial_tests.hpp.

References det(), eval_poly_matrix(), Interval< T, r >::M, and Interval< T, r >::m.

00274   {
00275     Interval<double> ev;
00276     Interval<double> * m;
00277     
00278     //Seq<POL> SS= b->system();
00279     unsigned dim= b->nbvar();
00280     
00281     // Check Jacobian's sign
00282     m= eval_poly_matrix( b->template domain<double>() ,J,dim);
00283     ev= det(m,dim);
00284     if ( ev.m*ev.M < 0 ) 
00285       return false;
00286 
00287     int td= topological_degree_2d<POL, double>(b);
00288 
00289     //if ( (td==-1 || td==1) )
00290     //{ std::cout<<"INCLUDE. ev="<<ev<<", td="<<td <<std::endl;
00291     //this->print();  }
00292     
00293     return ( (td==-1 || td==1) );
00294   }

bool mmx::realroot::include3 ( box_rep< POL > *  b,
polynomial< double, with< MonomialTensor > > *  J,
Seq< POL * > &  S0 
) [inline]

Inclusion criteria based on Rump's test.

Definition at line 298 of file solver_mv_monomial_tests.hpp.

References mmx::assign(), DPOL, mmx::eval(), eval_poly_matrix(), mmx::linear::LUinverse(), Interval< T, r >::M, and Interval< T, r >::m.

00299         {
00300           double *fc, *jc, *ijc;
00301           Interval<double> ev, *jb, *m;
00302           Seq<double> c= b->point( b->middle() ) ;
00303           Seq<Interval<double> > bx;
00304           unsigned i,j,n= b->nbvar();
00305           DPOL tmp;
00306 
00307           fc= new double[n];                // system S0 evaluated on c
00308           for(i=0; i< n;i++ )
00309           {
00310             let::assign(tmp, *S0[i]) ;
00311             tensor::eval( fc[i] , tmp.rep(), c , n  );
00312           }
00313           bx= b->template domain<double>(); // box 
00314           jc= eval_poly_matrix( c, J, n);   // Jacobian evaluated on c
00315           ijc= new double[n*n];
00316           linear::LUinverse(ijc, jc, n);     // Inverse of above
00317           jb= eval_poly_matrix( bx, J, n);  // Jacobian evaluated on box
00318           
00319           m= new Interval<double>[n*n];
00320           for(i=0; i< n;i++ )
00321             for(j=0; j< n;j++ )
00322               m[i*n+j] += ijc[i*n+j]*jb[j*n+i];
00323           free(jb);
00324 
00325           // check if -im*fc + (I-im*jb)*bx is inside bx.
00326           
00327           for(i=0; i< n;i++ )
00328           {
00329             bx[i]= bx[i]- c[i];
00330             ev= bx[i];
00331             for(j=0; j<n; j++ )
00332             {
00333               ev+= -ijc[i*n+j]*fc[j] - m[i*n+j] *bx[i];
00334             }
00335 //            std::cout<<ev<<" --- "<<bx[i]<<std::endl;
00336             if (  ev.m < bx[i].m ||
00337                   ev.M > bx[i].M )
00338             {
00339 //              std::cout<<ev<<" not in "<<bx[i]<<std::endl;
00340                   return false;
00341             }
00342           }
00343 //          std::cout<<"found root in "<<bx<< "(ev="<<ev<<")"<<std::endl;
00344             return ( true );
00345         }

polynomial<double,with<MonomialTensor> >* mmx::realroot::jacobian ( Seq< POL * > &  S0  )  [inline]

Definition at line 48 of file solver_mv_monomial_tests.hpp.

References mmx::assign(), mmx::diff(), DPOL, POL, and Seq< C, R >::size().

Referenced by solver_mv_monomial< FT, POL >::approximate(), and solver_mv_monomial< FT, POL >::isolate().

00049     {
00050       int n=S0.size();
00051       DPOL * jac;//[n*n]; 
00052       POL df;
00053       jac= new DPOL [n*n];
00054       
00055       for( int i=0;i<n;i++)
00056       {
00057         for( int k=0;k<n;k++)
00058         {
00059           df= diff( *(S0[i]),k);
00060           let::assign( jac[i*n+k] , df);
00061         }
00062       
00063       }
00064       return jac;
00065     }

bool mmx::realroot::miranda_test ( Seq< POL > &  S,
const int  i,
const int  j 
) [inline]

Definition at line 193 of file solver_mv_monomial_tests.hpp.

References mmx::tensor::face(), no_variation(), and POL.

Referenced by include1().

00194   {
00195     POL u,l;
00196     
00197     tensor::face(l, S[i], j, 0);
00198     tensor::face(u, S[i], j, 1);
00199     //std::cout<<"l="<<l<<",u="<<u<<std::endl;
00200 
00201     return  ( no_variation(l)    && 
00202               no_variation(u)    &&
00203               (l[0]>0)!=(u[0]>0) );
00204   };

bool mmx::realroot::no_variation ( POL &  p  )  [inline]

True iff all coefficients of p have the same sign.

Definition at line 39 of file solver_mv_monomial_tests.hpp.

References mmx::has_sign_variation().

Referenced by exclude3(), miranda_test(), and box_rep< POL >::miranda_test().

00040   {
00041     return !has_sign_variation(p.begin(),p.end()) ;
00042   };

void mmx::realroot::precondition ( Seq< T >  val,
polynomial< double, with< MonomialTensor > > *  mat,
int  n 
) [inline]

Definition at line 69 of file solver_mv_monomial_tests.hpp.

References mmx::eval(), mmx::rep(), and Seq< C, R >::size().

00070     {
00071       int m(val.size());
00072       T *ev;//[m*n];
00073       ev= new T [m*n];
00074 
00075       for( int i=0;i<m;i++)
00076       {
00077         for( int k=0;k<n;k++)
00078         {
00079           tensor::eval( ev[i*m+k] , mat[i*m+k].rep(),
00080                         val , m  );
00081         }
00082       }
00083       return ev;
00084     }

void mmx::realroot::scan_monom_env ( std::set< int > &  env,
const monom &  m 
) [inline]

Definition at line 127 of file solver_mv_bernstein.hpp.

References monom< C, E >::rep(), and mmx::size().

Referenced by solver< C, ProjRd< MTH > >::solve_monomial().

00128   { for ( unsigned i = 0; i < (m.rep()).size(); i++ ) if ( (m.rep())[i] != 0 ) env.insert(i); };

int mmx::realroot::sgn ( POL &  f,
FT  a 
) [inline]

Sign of f at point a>=0 (-1=inf).

Definition at line 423 of file solver_mv_monomial_tests.hpp.

References sgn().

00424   {
00425     if (a==FT(-1) ) 
00426       return sgn ( f[ f.rep().szs()[0]-1 ]);
00427     else 
00428       return sgn ( f(a) );
00429   };

int mmx::realroot::sgn ( FT  a  )  [inline]

Sign of a.

Definition at line 415 of file solver_mv_monomial_tests.hpp.

Referenced by parallel< system >::process(), sgn(), mmx::sign(), descartes_solver< real_t, local_method >::solve(), and topological_degree_2d().

00416   {
00417     if (a==FT(0) ) return 0;
00418     return ( (a>0) ? 1 : -1 );
00419   }

int mmx::realroot::signof ( polynomial< double, with< MonomialTensor > > *  p,
Seq< FT >  t 
) [inline]

Definition at line 162 of file solver_mv_monomial_tests.hpp.

References mmx::eval(), and Seq< C, R >::size().

00163     {
00164       
00165       Interval<double> ev;
00166       //int dim= p->nbvar();
00167       
00168       tensor::eval( ev , p->rep() , 
00169                     t , t.size() );
00170       
00171       if (ev< .1e-10) return (0);
00172       return (ev>0?1:-1);
00173     }

int mmx::realroot::signof_det ( T *  mat,
const int &  n 
) [inline]

Definition at line 108 of file solver_mv_monomial_tests.hpp.

References det().

00109     { // mat= double[n*n]
00110       double ev;
00111       ev= det(mat,n);
00112 
00113       if (ev< .1e-15) return (0);
00114       return (ev>0?1:-1);
00115     }

T* mmx::realroot::submatrix ( T *  matrix,
int  order,
int  i,
int  j 
) [inline]

Definition at line 118 of file solver_mv_monomial_tests.hpp.

Referenced by det().

00119     {
00120       int n=order-1;
00121       T * subm;//[n*n];
00122       int p, q;         // Indexes for matrix
00123       int a = 0, b;     // Indexes for subm
00124       subm = new T [n*n];
00125       
00126       for(p = 0; p < order; p++) {
00127         if(p==i) continue;      //Skip ith row
00128         
00129         b = 0;
00130         
00131         for(q = 0; q < order; q++) {
00132           if(q==j) continue;    //Skip jth column
00133           subm[a*order + b++] = matrix[p*order+q];
00134         }
00135         a++; //Increment row index
00136       }
00137       return subm;
00138     }

int topological_degree_2d ( box_rep< POL > *  b  )  [inline]

Definition at line 433 of file solver_mv_monomial_tests.hpp.

References assert, C, POL, sgn(), and Seq< C, R >::size().

00434   {
00435     /*-1=inf*/
00436     typedef typename POL::scalar_t C;
00437     
00438     assert( b->nbvar() ==2 );
00439     Seq<FT> p;
00440     Seq<FT> tmp;
00441     unsigned ev1(0), ev2(0), ev3(0), ev4(0);
00442     FT p1[2],p2[2];
00443     int tdeg;
00444     unsigned i;
00445     
00446       POL fl0, fr0,fl1,fr1, s,
00447         gl0, gr0,gl1,gr1  ;
00448 
00449       fl0=b->lface(0,1); gl0=b->lface(1,1);
00450       fr0=b->rface(0,1); gr0=b->rface(1,1);
00451       fl1=b->lface(0,0); gl1=b->lface(1,0);
00452       fr1=b->rface(0,0); gr1=b->rface(1,0);
00453       
00454       s= fl0 * gl0 ;
00455       p<< FT(0);p << FT(0) ;
00456       tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00457       for (i=0;i< tmp.size(); ++i )
00458         p<< tmp[i] << FT(0) ;
00459      
00460       s= fr1 * gr1; 
00461       p<< FT(-1);p << FT(0);
00462       tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00463       for (i=0;i< tmp.size(); ++i )
00464         p<< FT(-1)  << tmp[i]; 
00465 
00466       s=  fr0 * gr0 ;
00467       p<<  FT(-1);p << FT(-1) ;
00468       tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00469       for (i=0;i< tmp.size(); ++i )
00470         p<<  tmp[i] << FT(-1)  ;
00471 
00472       s= fl1 * gl1 ; 
00473       p<< FT(0);p << FT(-1) ;
00474       tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00475       for (i=0;i< tmp.size(); ++i )
00476         p<<  FT(0) << tmp[i]  ;
00477       
00478       p<< FT(0);p << FT(0) ;
00479 
00480       tdeg= 0;
00481       for (i=0; i<p.size()-3; i+=2)
00482       {
00483         p1[0]= p[i]  ;p1[1]=p[i+1];
00484         p2[0]= p[i+2];p2[1]=p[i+3];
00485         
00486         if (p1[1]==FT(0) && p2[1]==FT(0)){
00487           ev1= sgn( fl0, p1[0] );
00488           ev3= sgn( fl0, p2[0] );
00489           ev4= sgn( gl0, p1[0] );
00490           ev2= sgn( gl0, p2[0] );
00491         }else
00492           if (p1[0]==FT(-1) && p2[0]==FT(-1)){
00493             ev1= sgn( fr1, p1[1] );
00494             ev3= sgn( fr1, p2[1] );
00495             ev4= sgn( gr1, p1[1] );
00496             ev2= sgn( gr1, p2[1] );
00497           }else
00498             if (p1[1]==FT(-1) && p2[1]==FT(-1)){
00499               ev1= sgn( fr0, p1[0] );
00500               ev3= sgn( fr0, p2[0] );        
00501               ev4= sgn( gr0, p1[0] );
00502               ev2= sgn( gr0, p2[0] );
00503             }else
00504               if (p1[0]==FT(0)  && p2[0]==FT(0)){
00505                 ev1= sgn( fl1, p1[1] );
00506                 ev3= sgn( fl1, p2[1] );
00507                 ev4= sgn( gl1, p1[1] );
00508                 ev2= sgn( gl1, p2[1] );
00509               }
00510               else 
00511                 std::cout<< "TROUBLE!"<<std::endl;
00512         
00513         tdeg+=  ( (ev1)*(ev2) - (ev3)*(ev4) );
00514       }
00515       return (tdeg/8);
00516     }

Variable Documentation

const int C_ACCEPT = 5

Definition at line 32 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::accept(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::launch().

const int D_REJECT = 11

Definition at line 42 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::launch().

const int E_CTRL = 0

Definition at line 25 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::launch().

const int E_INIT = 4

Definition at line 29 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::launch().

const int E_RDSLV = 2

Definition at line 27 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::reduction().

const int E_SBDRL = 3

Definition at line 28 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), method< system, _strgy_, _rdslv_, _sbdrl_ >::launch(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::subdivision().

const int E_STRGY = 1

Definition at line 26 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::launch().

const int R_ERROR = 10

Definition at line 41 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error().

const int R_FAIL = 9

Definition at line 40 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::error(), method< system, _strgy_, _rdslv_, _sbdrl_ >::launch(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::reduction().

const int R_ISOK = 7

Definition at line 38 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::launch(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::reduction().

const int R_REJECT = 6

Definition at line 36 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::launch(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::reduction().

const int R_WEAK = 8

Definition at line 39 of file system_method.hpp.

Referenced by method< system, _strgy_, _rdslv_, _sbdrl_ >::launch(), and method< system, _strgy_, _rdslv_, _sbdrl_ >::reduction().

Generated on 6 Dec 2012 for realroot by  doxygen 1.6.1