Solvestacke2ter.hpp

Go to the documentation of this file.

#ifndef ALREADY_Solvestack
#define ALREADY_Solvestack

//#include"auxfun3.hpp"
#define NUMBITS 8
#define MASKBIT 255
#define LEVELSTACK 40
//#define NBTHREAD 7
#define NBTHREAD 7
template<>
void INIT_ELEM<stack_elem>(void * ptr, int size)
{
  int nbelem=size/sizeof(stack_elem);
  for(int k=0;k<nbelem;k++)
    {
     
      INIT_ELEM<COEFF>(&((stack_elem*)ptr)[k].coeff,sizeof(COEFF));
    }
}

template<>
void CLEAR_ELEM<stack_elem>(void * ptr, int size)
{
  int nbelem=size/sizeof(stack_elem);
  for(int k=0;k<nbelem;k++)
    {
     
      CLEAR_ELEM<COEFF>(&((stack_elem*)ptr)[k].coeff,sizeof(COEFF));
    }
}
#ifdef ASYNC2
template<typename compo, typename workspace, typename mon, typename T>
void vider_pile( stack_elem * stack, int stacksize, 
         pol<mon,T> *p, mon * &exceed_p, int &sizeexceed_p,workspace *w);
template<typename compo, typename workspace, typename mon, typename T>
void wrap_vider_pile( stack_elem * stack, int stacksize, 
              pol<mon,T> *p, mon * &exceed_p, int &sizeexceed_p, 
              //              int * nbstack,recursive_mutex &mymutex, workspace *w)
              int * nbstack,mutex &mymutex, workspace *w)
{
  //  workspace<COEFF> w;
  //  cout<<"je vide la pile "<<endl;
  //w.sizenf=0;
  //w.size=0;
  //w.tabnf=NULL;
  //w.tabcoeff=NULL;
  vider_pile<compo,workspace,mon,T>(stack,stacksize,p,exceed_p
                       ,sizeexceed_p,w);
  *nbstack=0;
  mymutex.unlock();
  //cout<<"j'ai vidé la pile "<<endl;
}

#endif
template<typename compo, typename workspace, typename mon, typename T>
void vider_pile( stack_elem * stack, int stacksize, 
         pol<mon,T> *p, mon * &exceed_p, int &sizeexceed_p,workspace *w)
{
  //le but de cette fonction est de realiser p-=la stack
  //en respectant aussi les exceed
  //TODO peut etre mettre les exceed sous forme de map (unoirdered cxx11) 
  //directement dans mult
  typedef typename pol<mon,T>::coeff_t coeff_t;
  typedef typename mon::rep_t rep_t;
  map<rep_t,coeff_t> mapexceed;
  int i,maxnfsize,j,compt,k;
  unsigned long *tmpptr;
  unsigned long mask=1;
  // cout<<"ds vider pile pour "<<p->ind<<endl;
  if (stacksize==0) return;
  for(i=0;i<stacksize;i++)
    for(;w->sizenf<stack[i].p->sizenf;increase_space(w));
  for(;w->sizenf<p->sizenf;increase_space(w));
  //
  for(i=0;i<w->size;i++) w->tabcoeff[i]=0;
  for(i=0;i<w->sizenf;i++) w->tabnf[i]=0;
   //
  compt=0;
  tmpptr=(unsigned long *)p->nfind;
  //#pragma omp parallel for shared(compt)
  for(i=0;i<p->sizenf/sizeof(unsigned long);i++)
    if(tmpptr[i])
      {
    ((unsigned long*)w->tabnf)[i]=tmpptr[i];
    mask=tmpptr[i];
    for(j=0;j<8*sizeof(unsigned long);j+=NUMBITS,mask>>=NUMBITS)
      if(mask&MASKBIT)
        {
          tabfun1[mask&MASKBIT](w->tabcoeff+(8*sizeof(unsigned long)*i+j),p->nf+compt);
          compt+=nbbits[mask&MASKBIT];
        }
      }
  //au tour des exceed
  for(i=0;i<sizeexceed_p;i++)
    mapexceed[exceed_p[i].rep]=exceed_p[i].GetCoeff();
  //
  for(i=p->sizenf-(p->sizenf % sizeof(unsigned long));i<p->sizenf;i++)
    //  for(i=0;i<p->sizenf;i++)
    if(p->nfind[i])
      {
    w->tabnf[i]=p->nfind[i];
    tabfun1[w->tabnf[i]](w->tabcoeff+(8*i),p->nf+compt);
    compt+=nbbits[p->nfind[i]];
      }
  //  cout<<"stacksize "<<stacksize<<endl;
  for(k=0;k<stacksize;k++)
    {
      compt=0;
      pol<mon,T> *tmpp=stack[k].p;
      T tmpcoeff=stack[k].coeff;
      //  cout<<"tmpcoeff "<<tmpcoeff<<endl;
      tmpptr=(unsigned long*)tmpp->nfind;
      //ici race condition omp sur comp
      //#pragma omp parallel for shared(compt)
      for(i=0;i<tmpp->sizenf/sizeof(unsigned long);i++)
    if(tmpptr[i])
      {
        ((unsigned long*)w->tabnf)[i]|=tmpptr[i];
        mask=tmpptr[i];
        for(j=0;j<8*sizeof(unsigned long);j+=NUMBITS,mask>>=NUMBITS)
          {
        if(mask&MASKBIT)
          {
            tabfun2[mask&MASKBIT]
              (w->tabcoeff+(8*sizeof(unsigned long)*i+j),
                  &tmpcoeff,tmpp->nf+compt);
            compt+=nbbits[mask&MASKBIT];
          }
#if 0   
        for(int k=0;k<8;k++)
          if(Iszero(w->tabcoeff[8*i*sizeof(unsigned long)+j+k]))
            ((unsigned long*)w->tabnf)[i]&=~(1<<(j+k));//mask;
#endif
          }
      }
      for(i=tmpp->sizenf-(tmpp->sizenf%sizeof(unsigned long));i<tmpp->sizenf;i++)
    if(tmpp->nfind[i])
      {
        w->tabnf[i]|=tmpp->nfind[i];
        mask=1;
        for(j=0;j<8;j++,mask<<=1)
          {
        if(tmpp->nfind[i]&mask)
          {
            //          cout<<"coeff dans w avant "<<w->tabcoeff[8*i+j]<<endl;
            //  cout<<"tmpcoeff "<<tmpcoeff<<endl;
            // cout<<"tmpp->nf[compt]; "<<tmpp->nf[compt]<<endl;
            w->tabcoeff[8*i+j]-=tmpcoeff*tmpp->nf[compt++];
            //cout<<"coeff dans w apres "<<w->tabcoeff[8*i+j]<<endl;
          }
        //if(Iszero(w->tabcoeff[8*i+j]))
        //  w->tabnf[i]&=~mask;
          }
      }
      //
      //ici on fait la meme operation sur les exceed!
      //
      //a voir//
      //
      //  cout<<"k ici "<<k<<" stack[k].sizeexceed "<<stack[k].sizeexceed<<endl;
      for(int i=0;i<stack[k].sizeexceed;i++)
    {
      //  cout<<"stack[k].exceed[i].rep "<<stack[k].exceed[i].rep<<" "<<stack[k].exceed[i].GetCoeff()<<" tmpcoeff "<<tmpcoeff<<endl;
      mapexceed[stack[k].exceed[i].rep]-=stack[k].exceed[i].GetCoeff()*tmpcoeff;
    }
    }
  compt=0;j=-1;
  for(i=0;i<w->sizenf;i++)
    if(w->tabnf[i])
      {
    for(int k=0;k<8;k++)
      if(Iszero(w->tabcoeff[8*i+k]))
        w->tabnf[i]&=~(1<<k);
    if(w->tabnf[i])
      {
        j=i;
        compt+=nbbits[w->tabnf[i]];
      }
      }
  //  cout<<" ici la size "<<compt<<endl;
  p->nf=(T*)MAC_REV_REALLOC<T>(p->nf,p->size*sizeof(T),
                   compt*sizeof(T));
  p->size=compt;
  p->nfind=(unsigned char*)realloc(p->nfind,j+1);
  p->sizenf=j+1;
  compt=0;
  for(i=0;i<=j;i++)
    {
      
      p->nfind[i]=w->tabnf[i];
      w->tabnf[i]=0;
      if(p->nfind[i])
    {
      unsigned char mask=1;
      for(k=0;k<8;k++,mask<<=1)
        {
          if(p->nfind[i]&mask)
        {
          p->nf[compt++]=w->tabcoeff[8*i+k];
          w->tabcoeff[8*i+k]=0;
        }
        }
    }
    }
  //
  //Il reste a reconstruire un exceed à partir de la mapexceed
  //
#ifdef DEBUG
 {
  cout<<"dans vider pile exceed a la fin"<<endl;
  auto itermap=mapexceed.begin();
  for(int i=0;itermap!=mapexceed.end();i++,itermap++)
    cout<<mon(itermap->second,itermap->first)<<" ---- ";
  cout<<endl;
  } 
#endif
  MAC_REV_FREE<mon>(exceed_p,sizeexceed_p*sizeof(mon));
  exceed_p=(mon*)MAC_REV_MALLOC<mon>(mapexceed.size()*sizeof(mon));
  auto itermap=mapexceed.begin();
  for(int i=0;itermap!=mapexceed.end();i++,itermap++)
    exceed_p[i]=mon(itermap->second,itermap->first);
  sizeexceed_p=mapexceed.size();
}

template<typename coeff, typename workspace, typename component>
void SolveU(const harewell<coeff> & u, component &compo, workspace &w)
{
  // typedef list<pol<mon,COEFF> > typPk;
  // typedef pol<mon,COEFF> elem_typ;
  typedef typename component::pol pol;
  int t;
  const  int taille=compo.poly.size();
  auto &tmpPk=compo.poly;//c le vecteur des polynomes de la component
  stack_elem *stack=(stack_elem*)MAC_REV_MALLOC<stack_elem>
    (LEVELSTACK*(taille+1)*sizeof(stack_elem));
  int *nbstack=(int*)malloc((taille+1)*sizeof(int));
#ifdef ASYNC2
  list<std::thread > listthread;
  mutex *tabmut=new mutex[taille];
  vector<workspace> vw(taille);
  for(auto w :vw)
    {
      w.sizenf=0;
      w.size=0;
      w.tabnf=NULL;
      w.tabcoeff=NULL;
    }
#endif
  cout<<"ds solveU"<<endl;
  for(int i=0;i<taille;i++)
    nbstack[i]=0;

  //copy(tmpPk.begin(),Pk.begin(),Pk.end());
  //  for(typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
  //  tmpPk.push_back(*iter);
  //Pk.erase(Pk.begin(),Pk.end());
  //cout<<"u.ncol "<<u.ncol<<" u.nrow "<<u.nrow<<endl; 
  for(int i=u.ncol-1;i>=0;i--)
      {
    coeff mult_coeff=(coeff)1;
    int j;
    //
#ifdef ASYNC2
    tabmut[i].lock();
#endif
      vider_pile<component,workspace,mon,COEFF>
    (stack+i*LEVELSTACK,nbstack[i],&(tmpPk[i]),compo.exceed[i],
           compo.sizeexceed[i],&w);

#ifdef ASYNC2
      tabmut[i].unlock();
#endif
    //
    nbstack[i]=0;
    for (j=((i==(u.ncol-1))?(u.nnz-1):(u.colptr[i+1]-1))
           ;(j>=(u.colptr)[i]);j--)
      if (i==u.rowind[j]) 
        {
          mult_coeff=((COEFF*)u.nzval)[j];
          break;
        }
    if(Iszero(mult_coeff)) continue;
    mult_coeff=((coeff)1)/mult_coeff;
    tmpPk[i]*mult_coeff;
    //on reporte aussi sur l'exceed les modifs
    for(j=0;j<compo.sizeexceed[i];j++)
      compo.exceed[i][j]*=mult_coeff;
    //fin du report
    for (j=((i==(u.ncol-1))?(u.nnz-1):(u.colptr[i+1]-1))
           ;(j>=(u.colptr)[i]);j--)
      {
        if (i!=u.rowind[j])//t==j 
          {
#ifdef ASYNC2
        tabmut[u.rowind[j]].lock();
#endif

        if(((COEFF*)u.nzval)[j]!=(coeff)0)
          {
          stack[u.rowind[j]*LEVELSTACK+nbstack[u.rowind[j]]].coeff
            =((COEFF*)u.nzval)[j];
          stack[u.rowind[j]*LEVELSTACK+nbstack[u.rowind[j]]].p
            =&tmpPk[i];
          stack[u.rowind[j]*LEVELSTACK+nbstack[u.rowind[j]]].exceed=compo.exceed[i];
          stack[u.rowind[j]*LEVELSTACK+nbstack[u.rowind[j]]].sizeexceed=
            compo.sizeexceed[i];
          nbstack[u.rowind[j]]++;     
      if(nbstack[u.rowind[j]]==LEVELSTACK)
            {
#ifndef ASYNC2
              vider_pile<component,workspace,mon,COEFF>
            (stack+u.rowind[j]*LEVELSTACK,LEVELSTACK
             ,&tmpPk[u.rowind[j]],compo.exceed[u.rowind[j]],
             compo.sizeexceed[u.rowind[j]],&w);
    
         // ,&tmpPk[u.rowind[j]],compo.exceed[u.rowind[i]],
             // compo.sizeexceed[u.rowind[i]],&w);
              nbstack[u.rowind[j]]=0;
#else
             //       tabmut[u.rowind[j]].lock();//on lock le mutex ici
              //et c le thread qui delockera qd il aura fini
              listthread.push_back(std::thread(
                wrap_vider_pile<component,workspace,mon,COEFF>,
            stack+u.rowind[j]*LEVELSTACK,LEVELSTACK
            ,&tmpPk[u.rowind[j]],std::ref(compo.exceed[u.rowind[j]]),
            std::ref(compo.sizeexceed[u.rowind[j]]),
            nbstack+u.rowind[j],
            std::ref(tabmut[u.rowind[j]]),&vw[u.rowind[j]]));
              if(listthread.size()==NBTHREAD)
            {   
              for(auto iterfuture=listthread.begin();
                   iterfuture!=listthread.end();iterfuture++)
                iterfuture->join();
              listthread.erase(listthread.begin(),listthread.end());
            }      
#endif
            }
#ifdef ASYNC2
      else
        {
         tabmut[u.rowind[j]].unlock(); 
        }
#endif
          }
#ifdef ASYNC2
        else 
          {
            tabmut[u.rowind[j]].unlock();
          }
#endif


          }
      }
      }
#ifdef ASYNC2
for(auto iterfuture=listthread.begin();
    iterfuture!=listthread.end();iterfuture++)
  iterfuture->join();
 listthread.erase(listthread.begin(),listthread.end());
 for(auto w : vw)
   destroy_space(&w);
#endif
  //  for(int i=0;i<u.ncol;i++)
  //   cout<<"u nbstack "<<i<<" "<<nbstack[i]<<endl;
#ifdef PASCOMPRIS
  typename list<pol >::iterator iterPk=Pk.begin();
  for(typename vector<pol >::iterator iter=tmpPk.begin()
    ;iter!=tmpPk.end();iter++)
        Pk.push_back(*iter);
#endif
  MAC_REV_FREE<stack_elem>(stack,(LEVELSTACK*taille*sizeof(stack_elem)));
  free(nbstack);
#ifdef DEBUG
  cout<<"apres SolveU"<<endl;
  for(int i=0;i<tmpPk.size();i++)
    {
      cout<<"compo.sizeexceed[i]<<endl; "<<compo.sizeexceed[i]<<endl;
      copy(compo.exceed[i],compo.exceed[i]+compo.sizeexceed[i],ostream_iterator<mon>(cout," "));
      cout<<endl;
      copy(compo.poly[i].nf,compo.poly[i].nf+compo.poly[i].size,ostream_iterator<coeff>(cout," "));
      cout<<endl;
      cout<<"compo "<<compo.poly[i].ind<<endl;
    }
#endif
  //return Pk;
}

template<typename workspace, typename component, typename coeff>
void SolveL(const harewell<coeff> & l, component &compo, workspace &w)
//list<pol<mon,COEFF> >& SolveL(const harewell<COEFF> &l
//                ,list<pol<mon,COEFF> > &Pk
//                ,workspace &w)
{
  //  typedef list<pol<mon,COEFF> > typPk;
  //typedef pol<mon,COEFF> elem_typ;
  //vector<pol<mon,COEFF> > tmpPk;
  
 typedef typename component::pol pol_t;
  int t;
  const  int taille=l.nrow;//compo.poly.size();
  auto &tmpPk=compo.poly;//c le vecteur des polynomes de la component 
  
  stack_elem *stack=(stack_elem*)MAC_REV_MALLOC<stack_elem>
    (LEVELSTACK*(taille)*sizeof(stack_elem));
  int *nbstack=(int*)malloc((taille)*sizeof(int));
#ifdef ASYNC2
  list<std::thread > listthread;
  mutex *tabmut=new mutex[taille];
  vector<workspace> vw(taille);
  for(auto w :vw)
    {
      w.sizenf=0;
      w.size=0;
      w.tabnf=NULL;
      w.tabcoeff=NULL;
    }

#endif
  
  for(int i=0;i<taille;i++)
    nbstack[i]=0;

#ifdef DEBUG
cout<<"avant SolveL"<<endl;
  for(int i=0;i<tmpPk.size();i++)
    {
      cout<<"compo "<<compo.poly[i].ind<<endl;
      copy(compo.exceed[i],compo.exceed[i]+compo.sizeexceed[i],ostream_iterator<mon>(cout," "));
      copy(compo.poly[i].nf,compo.poly[i].nf+compo.poly[i].size,ostream_iterator<mon>(cout," "));
      cout<<endl;

    }
#endif
  //ALGO A REVOIR

  //  for(typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
  //  tmpPk.push_back(*iter);
  //  Pk.erase(Pk.begin(),Pk.end());
    // cout<<"l.ncol "<<l.ncol<<" l.nrow "<<l.nrow<<endl;
  //  auto itertmppk=tmpPk.begin();
  for(int i=0;i<l.ncol;i++)
    {
      int j=l.colptr[i];
#ifdef DEBUG      
      cout<<"vider pile"<<endl;
      copy(compo.exceed[i],compo.exceed[i]+compo.sizeexceed[i],ostream_iterator<mon>(cout," "));
        cout<<endl;
    cout<<"compo "<<compo.poly[i].ind<<endl;
    copy(compo.poly[i].nf,compo.poly[i].nf+compo.poly[i].size,ostream_iterator<coeff>(cout," --- "));
#endif   
#ifdef ASYNC2
    tabmut[i].lock();
#endif

      vider_pile<component,workspace,mon,COEFF>
    (stack+i*LEVELSTACK,nbstack[i],&(tmpPk[i]),compo.exceed[i],
           compo.sizeexceed[i],&w);
#ifdef ASYNC2
    tabmut[i].unlock();
#endif

      nbstack[i]=0;
      for (;(j<l.nnz)&&(j<(l.colptr)[i+1]);j++)
      {
    if (i!=l.rowind[j]) 
      {
#ifdef ASYNC2
            tabmut[l.rowind[j]].lock();
#endif

        if((((COEFF*)l.nzval)[j])!=(COEFF)0)
          {
        //avant#ifdef ASYNC2
        //    tabmut[l.rowind[j]].lock();
        //#endif

        //      cout<<"L:rowind "<<l.rowind[j]<<endl;
        //  cout<<"nbstack "<<nbstack[l.rowind[j]]<<endl;
        //cout<<"taille*LEVELSTACK "<<taille*LEVELSTACK<<endl;
            stack[l.rowind[j]*LEVELSTACK+nbstack[l.rowind[j]]].coeff
              =((COEFF*)l.nzval)[j];
            stack[l.rowind[j]*LEVELSTACK+nbstack[l.rowind[j]]].p
              =&(tmpPk[i]);
            stack[l.rowind[j]*LEVELSTACK+nbstack[l.rowind[j]]].exceed
              =compo.exceed[i];
            stack[l.rowind[j]*LEVELSTACK+nbstack[l.rowind[j]]].sizeexceed
              =compo.sizeexceed[i];
            nbstack[l.rowind[j]]++;
            //AVANT #ifdef ASYNC2
              //avnat tabmut[l.rowind[j]].unlock();
            //avant#endif
            //cout<<"nbstack ici"<<nbstack[l.rowind[j]]<<endl;
            if(nbstack[l.rowind[j]]==LEVELSTACK)
              {
#ifndef ASYNC2
            vider_pile<component,workspace,mon,COEFF>
              (stack+l.rowind[j]*LEVELSTACK,LEVELSTACK
               ,&tmpPk[l.rowind[j]],compo.exceed[l.rowind[j]],
               compo.sizeexceed[l.rowind[j]],&w);
            //ICI ON A BESOIN de l'acces
            nbstack[l.rowind[j]]=0;
#else
            //cout<<"je declenche une vidange"<<endl;
            //tabmut[l.rowind[j]].lock();//on lock le mutex tout de suite
            listthread.push_back(std::thread(
             wrap_vider_pile<component,workspace,mon,COEFF>,
             stack+l.rowind[j]*LEVELSTACK,LEVELSTACK
             ,&tmpPk[l.rowind[j]],std::ref(compo.exceed[l.rowind[j]]),
             std::ref(compo.sizeexceed[l.rowind[j]]),
             nbstack+l.rowind[j],std::ref(tabmut[l.rowind[j]]),
             &vw[l.rowind[j]]));
              if(listthread.size()==NBTHREAD)
            {   
              for(auto iterfuture=listthread.begin();
                   iterfuture!=listthread.end();iterfuture++)
                iterfuture->join();
              listthread.erase(listthread.begin(),listthread.end());
            }      
#endif
        
              
              }
#ifdef ASYNC2
            else//on est non nul et on vide pas la pile
              {
            tabmut[l.rowind[j]].unlock();
              }
#endif
          }
#ifdef ASYNC2   
        else//coeffciient nul
          {
        tabmut[l.rowind[j]].unlock();
          }
#endif
      } 
      }
    }
  // for(int i=0;i<l.ncol;i++)
  //   cout<<"nbstack "<<i<<" "<<nbstack[i]<<endl;
     //vider_pile(stack+i*LEVELSTACK,nbstack[i],&tmpPk[i],&w);
  //  list<pol<mon, COEFF> >::iterator iterPk=Pk.begin();
  //for(vector<pol<mon , COEFF> >::iterator iter=tmpPk.begin()
  //    ;iter!=tmpPk.end();iter++,iterPk++)
  // Pk.push_back(*iter);
#ifdef ASYNC2
for(auto iterfuture=listthread.begin();
    iterfuture!=listthread.end();iterfuture++)
  iterfuture->join();
 listthread.erase(listthread.begin(),listthread.end());
 for(auto &w : vw)
     destroy_space(&w);
#endif

  MAC_REV_FREE<stack_elem>(stack,(LEVELSTACK*taille*sizeof(stack_elem)));
  free(nbstack);
#ifdef DEBUG 
cout<<"apres SolveL"<<endl;
  for(int i=0;i<tmpPk.size();i++)
    {
      copy(compo.exceed[i],compo.exceed[i]+compo.sizeexceed[i],ostream_iterator<mon>(cout," "));
      cout<<endl;
      copy(compo.poly[i].nf,compo.poly[i].nf+compo.poly[i].size,ostream_iterator<coeff>(cout," "));
      cout<<endl;
    
      cout<<"compo "<<compo.poly[i].ind<<endl;
    }
#endif
  //return Pk;
}


template<typename typmat ,typename typPk, typename mon, typename typMk, typename workspace,
typename typserv>
//typserv a virer ne sert que pour le debug 
void
Solve(typmat &S, typmat &L,
      component<typPk,mon> &compo,int **pr,int **pc
      ,typMk &tmpMk, workspace &w, typserv & serv)
{
 typmat U;
 int *Pr,*Pc,maxsize=0;
 int *tmppr;
 int *tmppc;
 struct timeval t1,t2;
 //#ifdef DEBUG
 // vector<mon> tmpMk(compo.poly.size());
 //#endif
 gettimeofday(&t1,NULL);
 LUdecomp(S,&L,&U,&Pr,&Pc);
 gettimeofday(&t2,NULL);
 cout<<"\e[32m Temps LU "<<t2.tv_sec-t1.tv_sec+1e-6*(t2.tv_usec-t1.tv_usec);
 cout<<"\e(B\e[m"<<endl;
 t1=t2;
#if 0
 cout<<"L"<<endl;
affdebug(L);
 cout<<"R"<<endl;
 affdebug(S); 
cout<<"PR "<<endl;
 copy(Pr,Pr+L.nrow,ostream_iterator<int>(cout," "));cout<<endl; 
#endif
 tmppr=(int*)malloc(sizeof(int)*L.nrow);
 tmppc=(int*)malloc(sizeof(int)*L.nrow);

#if 0
for(int i=0; i<compo.poly.size();i++)
   {
     tmpMk[i]=compo.poly[i].ind;
     tmpMk[i].SetCoeff(1);
   }
#endif


 inverse(Pr,L.nrow);
 copy(Pr,Pr+L.nrow,tmppr); 
 Dopermu(compo.poly,tmppr);
 copy(Pr,Pr+L.nrow,tmppr); 
 Dopermu(compo.exceed,tmppr);
 copy(Pr,Pr+L.nrow,tmppr) ;
 Dopermu(compo.sizeexceed,tmppr);
 free(tmppr);


 copy(Pc,Pc+S.ncol,tmppc); 
 inverse(tmppc,S.ncol);
 Dopermu(tmpMk,tmppc);

 free(tmppc);



  
 auto tmpPk=compo.poly; 
#ifdef DEBUG
 cout<<"avant SolveL"<<endl;
 for(int i=0;i<tmpPk.size();i++)
   {
     cout<<"compo "<<compo.poly[i].ind<<endl;
     copy(compo.exceed[i],compo.exceed[i]+compo.sizeexceed[i],ostream_iterator<mon>(cout," "));
     cout<<endl;
     copy(compo.poly[i].nf,compo.poly[i].nf+compo.poly[i].size,ostream_iterator<typename typPk::value_type::coeff_t>(cout," "));
     cout<<endl;
     cout<<invconv<Poly>(compo.poly[i],serv)<<endl;
     cout<<"----------------------------------------------"<<endl;
   }
#endif
 SolveL(L,compo,w);
 //
 reorder_ind(S,tmpMk,compo,serv);
 //

#ifdef DEBUG
cout<<"apres SolveLLL"<<endl;
 for(int i=0;i<tmpPk.size();i++)
   {
     cout<<"l'ind "<<compo.poly[i].ind<<endl;    
     for(int j=0;j<S.ncol;j++)
       if(S(i,j)!=0) 
     cout<<tmpMk[j]*S(i,j)<<"+ ";
     copy(compo.exceed[i],compo.exceed[i]+compo.sizeexceed[i],ostream_iterator<mon>(cout," "));
     cout<<endl;
     cout<<invconv<Poly>(compo.poly[i],serv)<<endl;
   }
#endif
 SolveU(S,compo,w);//SolveU(S,TildePk);
 *pc=Pc;*pr=Pr;
 gettimeofday(&t2,NULL);
 cout<<"\e[32m Temps RESTE "<<t2.tv_sec-t1.tv_sec+1e-6*(t2.tv_usec-t1.tv_usec);
 cout<<"\e(B\e[m"<<endl;
 t1=t2;
}

template<typename typmat, typename typMk, typename component, typename typserv>
void reorder_ind(const typmat &tmpmat, const typMk &verytmpMk, component &compo, typserv &serv) 
{
  //ici on remet les ind a la bonne valeur 
  //on s'occuppe aussi de rajouter en exceed les monomes des colonnes nulles
  //on suppose verytmpMk deja permute par Pc
  typedef typename typMk::value_type::coeff_t coeff;
  int j=0;
  typename typMk::const_iterator itermon=verytmpMk.begin();
  
  for(auto iter_in_poly=compo.poly.begin();
      itermon!=verytmpMk.end();iter_in_poly++,itermon++,j++)
    {
      //      cout<<" dans verytmp "<<*itermon<<
      //    invconv<Poly>(*iter_in_poly,serv)<<endl;
      if(!Iszero(tmpmat(j,j)))
    {
      iter_in_poly->ind=*itermon;
      iter_in_poly->ind.SetCoeff(1);
    }
      else
    {
      //On rajoute aux exceed des autre polyn\çomes de la composante les
      //les monomes correspondant a cette colonne 
      for(int i=tmpmat.colptr[j];i<tmpmat.colptr[j+1];i++)
        if(!(((coeff*)tmpmat.nzval)[i]==0))
          {
        //il y a eu permutation des colonnes il
        // faut donc se fier a itermon plutot qu'aux ind
        mon tmpmon=*itermon;
        coeff tmpcoeff=((coeff*)tmpmat.nzval)[i];
        
        tmpmon.SetCoeff(tmpcoeff);
        compo.exceed[tmpmat.rowind[i]]=(mon*)
          MAC_REV_REALLOC<mon>(compo.exceed[tmpmat.rowind[i]],
                       compo.sizeexceed[tmpmat.rowind[i]]*sizeof(mon),
                        (compo.sizeexceed[tmpmat.rowind[i]]+1)*sizeof(mon));
        
        compo.exceed[tmpmat.rowind[i]][compo.sizeexceed[tmpmat.rowind[i]]]=tmpmon;
        compo.sizeexceed[tmpmat.rowind[i]]++;
          }
      iter_in_poly->ind=mon(0);
      
    }
    }
}

#endif //ALREADY_Solvestack