resolve.hpp

Go to the documentation of this file.

//#include<g2c.h>
//struct doublecomplex 
//{
//  double r;
//  double i;
//};
typedef double doublereal;
typedef long int logical;

#include "linalg/clapack.h"
#include<complex>
#include"borderbasix/linalg/Lapack.H"
#include"borderbasix/linalg/Eigen.H"
#include<stdlib.h>
//extern "C" {
  //bon je sais la decalration ci apres est un peu crade mais bon 
  //entre f2c et g2c tout n'est pas tjrs compatible alors comme ca 
  //meme avec le -Wall on cloue le bec a gcc.. :-0
//  void ztrsen_(...);
//  void zgebal_(...);
//  void zgees_(...);
//  void dgees_(...);
  //  void dgees_(char *JOBVS,char *SORT,bool (*SELECT)(double,double),int *N,double *A,int *LDA,int *SDIM,double *WR,double *WI,double *VS, int *LDVS,double *WORK,int *LWORK,double *BWORK,int *info);
  //  void zgees_(char *JOBVS,char *SORT,logical (*SELECT)(doublecomplex *), tla_int *N,doublecomplex *A,tla_int *LDA,tla_int *SDIM,doublecomplex *W,doublecomplex *VS,tla_int *LDVS,doublecomplex *WORK,tla_int *LWORK,doublereal * RWORK,tla_int *BWORK,tla_int *info);
//}

//wrapper au dessus de lapack pour calculer 
//une decomposition d schur
//de la matrice qui doit etre un lapack<double>
template<typename Mat>
Mat wraplapack_dgees(const Mat &M)
{
    typedef typename Mat::rep_t rept;
    rept m=M.rep,res;
    bool (*SELECT)(double,double)=(bool (*)(double,double))1;
    double *bwork=(double*)malloc(m.nbcol()*sizeof(void (*)(void)));
    double *wr=(double*)malloc(m.nbrow()*sizeof(double));
    double *wi=(double*)malloc(m.nbrow()*sizeof(double));
    Mat Mres(m.nbrow(),m.nbrow());
    double *work=(double*)malloc(10*sizeof(double)*m.nbrow());
    int info,n,i,l,lres,lwork,sdim;
    double* tab =(double*)malloc(m.nbcol()*m.nbrow()*sizeof(double));
    double* tabres =(double*)malloc(m.nbcol_*m.nbrow_*sizeof(double));
    n=m.nbrow();l=m.nbrow();lres=m.nbrow();lwork=10*m.nbrow();
    for(i=0;i<m.nbcol_*m.nbrow_;i++) tab[i]=m.tab_[i];
    dgees_("V","N",SELECT,&n,tab,&l,&sdim,wr,wi,tabres,&lres,work,
            &lwork,bwork,&info);
    (Mres.rep).reserve(m.nbcol_,m.nbcol_);
    for(i=0;i<m.nbcol_*m.nbrow_;i++) (Mres.rep).tab_[i]=tabres[i];
    free(wr);
    free(wi);
    free(work);
    free(tab);
    free(tabres);
    return Mres;
}
//fonction qui sert a rien.. :-)

extern "C"{
long int select_true(doublecomplex *toto)
{
  static long int rep=1;
  return rep;
}
}

//la meme que wraplapack_dgees 
//sauf qu'elle marche sur le matrice du type
//lapack<complex<double> >
//Si PbS sous linux verifie la fonction dznrm2!!!
//Elle est pas compliee de maniere adequate dans la distrib semir
template<typename Mat>
Mat wraplapack_zgees(Mat M)
{
      typedef typename Mat::rep_t rept;
      float toto;
      rept m=M.rep;
      logical (*SELECT)(doublecomplex *)=select_true;
      doublecomplex *w=(doublecomplex *)malloc(m.nbrow_*sizeof(doublecomplex));
      double *rwork=(double*)malloc(m.nbrow_*sizeof(double));
      logical *bwork=(logical*)malloc(m.nbrow_*sizeof(logical));
      Mat Mres(m.nbrow_,m.nbrow_);
      doublecomplex* work=
    (doublecomplex*)malloc(10*sizeof(doublecomplex)*m.nbrow_);
      int info,n,l,lres,lwork,sdim;
      unsigned int i;
      doublecomplex* tab =
    (doublecomplex*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex));
      doublecomplex* tabres =
    (doublecomplex*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex));
      char *JOBS="V",*SORT="N";
      
      n=m.nbrow_;l=m.nbrow_;lres=m.nbrow_;lwork=5*m.nbrow_;
      
      for(i=0;i<m.nbcol_*m.nbrow_;i++) {
    tab[i].r=m.tab_[i].real();
    tab[i].i=m.tab_[i].imag();
    toto=tab[i].r;
      }
      zgees_(JOBS,SORT,SELECT,&n,tab,&l,&sdim,
         w,tabres,&lres,work,&lwork,rwork,bwork,&info);
      (Mres.rep).reserve(m.nbcol_,m.nbcol_);//matrice carree exclusivement!!
      for(i=0;i<m.nbcol_*m.nbrow_;i++) 
        {
            (Mres.rep).tab_[i]=complex<double>(tab[i].r,tab[i].i);
        toto=*(((float*)tabres)+i+2);
        }
      cout<<"info "<<info<<endl;
      //cout<<Mres<<endl<<endl;
      for(i=0;i<m.nbcol_*m.nbrow_;i++) 
            {
        (Mres.rep).tab_[i]=complex<double>(tabres[i].r,tabres[i].i);
            }
      free(w);
      free(rwork);
      free(work);
      free(bwork);
      free(tab);
      free(tabres);
      return Mres;
}

//debut de l'implementation de l'algo de Corless dans 
//ISSAC 97

//trigo simultanee sur des doubles
template<typename Mat>
void trigsim(list<Mat> &lmat)
{
    //on supose que la list n'est pas vide!!!
    Mat Shurvect=wraplapack_dgees(lmat.front());
    Mat invShurvect;
    Shurvect.transpose();//on est reel jusque la
    invShurvect=Shurvect;
    Shurvect.transpose();
    typename list<Mat>::iterator iter;
    for(iter=lmat.begin();iter!=lmat.end();iter++)
    {
        Mat tmp=invShurvect*(*iter)*Shurvect;
        *iter=tmp;
    }
}
//un wrap de lapack rien d'important
template<typename Mat>
void wraplapack_ztrsen(Mat &M,logical * select,doublereal *res)
{
  char *JOBS="E",*COMPQ="N";
  int N=M.nbcol(),LDT=M.nbcol();
  doublecomplex *T=(doublecomplex*)
    malloc(M.nbcol()*M.nbcol()*sizeof(doublecomplex));
  doublecomplex *Q=NULL;
  doublecomplex *W=(doublecomplex*)malloc(M.nbcol()*sizeof(doublecomplex));
  int m=1,LDQ=1;
  doublereal S=1,SEP=2;
  doublecomplex *WORK=(doublecomplex*)malloc(2*N*N*sizeof(doublecomplex));
  int LWORK=2*N*N,info;
  //maintenant on initialise le tableau T
  for(int i=0;i<N*N;i++)
        {
      T[i].r=(M.rep).tab_[i].real();
      T[i].i=(M.rep).tab_[i].imag();
    }
  ztrsen_(JOBS,COMPQ,select,&N,T,&LDT,Q,&LDQ,W,&m,&S,&SEP,WORK,&LWORK,&info);
  *res=S;
  free(WORK);
  free(W);
  free(T);
}
//evalue le conditionnement non clsuetrise de 
//chacunes des valeurs propres de la matrice
  
template<typename Mat>
doublereal cond_number(Mat &Combilin,int i)
{
  doublereal res;
  logical *select;
  select=(logical*)calloc(Combilin.nbcol(),sizeof(logical));
  select[i]=1;
  wraplapack_ztrsen(Combilin,select,&res);
  return res;
}
//
//creation des clusters 
//utilisation de l'heuritique de Stetter!!
//

template<typename Mat,typename typcond>
int getcluster(Mat &M,typcond *cond)
{
  double EPSILON=1e-16;
  doublecomplex *truncnum=(doublecomplex*)
    malloc(M.nbcol()*sizeof(doublecomplex));
  int *digits=(int*)malloc(M.nbcol()*sizeof(int));
  double normM=0;
  char tampon1[60],tampon2[60],*rien;
  int comp=0;
  for(int i=0;i<M.nbcol();i++)
    normM=max(sqrt((M(i,i)*conj(M(i,i))).real()),normM);
  for(int i=0;i<M.nbcol();i++)
    {
      // cout<<"cond[i]"<<cond[i]<<endl;
      double lambda=1.0/cond[i]*EPSILON*normM;
      double nbdigits=(1/lambda);
      int k;
      for(k=1;nbdigits>1;k++,nbdigits/=10);//
      //cout<<"nbdigits "<<nbdigits<<endl;;
      nbdigits=k;
      digits[i]=k;
      //rajout des tests de taille apres
      snprintf(tampon1,60,"%%.%dg\0",(k>16)?16:k);
      //cout<<tampon1<<endl;
      snprintf(tampon2,60,tampon1,M(i,i).real());
      //cout<<tampon2<<endl;
      truncnum[i].r=strtod(tampon2,&rien);
      snprintf(tampon2,60,tampon1,M(i,i).imag());
      //cout<<tampon2<<endl;
      truncnum[i].i=strtod(tampon2,&rien);
    }
  for(int i=0;i<M.nbcol();i++)
    {
      //if(cond[i]>=0)
      {
    //  cout<<"toto"<<endl;
      cond[i]=-i;
      for(int j=0;j<M.nbcol();j++)
        {
          double val=1.0;
          //cout<<"digitis_i "<<digits[i]<<" "<<digits[j]<<endl;
          for(int k=0;k<min(digits[i],digits[j]);k++,val*=10);
          //cout<<"val "<<val<<endl;
          //cout<<"la diference "<<abs(truncnum[j].r-truncnum[i].r)<<endl;
          //cout<<" en complex "<<abs(truncnum[j].i-truncnum[i].i)<<endl;
          // cout<<"prod "<<abs(truncnum[j].r-truncnum[i].r)*val<<endl;
          if ((abs(truncnum[j].r-truncnum[i].r)*val<1.0)
          &&(abs(truncnum[j].i-truncnum[i].i)*val<1.0))
        {
          typcond stock=cond[j];
          for(int k=0;k<M.nbcol();k++)
            if(cond[k]==stock) cond[k]=-i;
          //on merge les cluster
          //truncnum[i].r=-1;
          //      cout<<"clustering ok plus d'un"<<endl;
        }
        }
    }
    }
  //cout<<"cond"<<endl;
  //for(int k=0;k<M.nbcol();k++)
  //  cout<<cond[k]<<endl;
  free(truncnum);
  return M.nbcol();//cluster numerotes de -1 a ...
}

//
//appliquation des clusters crees avec la fonction d'avant 
//

template<typename Mat>
void docluster(Mat &M,Mat &Shurvect,doublereal * cond,int i)
{
  //0n recherche tous les -i dans le tableau cond et apres on
  //fait un cluster avec eux......
  logical *select=(logical*)malloc(M.nbcol()*sizeof(logical));
  int N=M.nbcol(),LDT=M.nbcol(),LDQ=M.nbcol();
  int m,info,LWORK=N*N*2;
  int flag=0;
  double *stockcond=(double*)malloc(N*sizeof(double));
  double *WORK=(double*)malloc(2*sizeof(double)*LWORK);
  double S,SEP;
  doublecomplex *stockW=(doublecomplex*)malloc(N*sizeof(doublecomplex));
  doublecomplex *W=(doublecomplex*)malloc(N*sizeof(doublecomplex));
  doublecomplex *tmpM=(doublecomplex*)malloc(N*N*sizeof(doublecomplex));
  doublecomplex *tmpQ=(doublecomplex*)malloc(N*N*sizeof(doublecomplex));
  for(int k=0;k<N*N;k++)
    {
      tmpM[k].r=(M.rep).tab_[k].real();
      tmpM[k].i=(M.rep).tab_[k].imag();
    }
  for(int k=0;k<N*N;k++)
    {
      tmpQ[k].r=(Shurvect.rep).tab_[k].real();
      tmpQ[k].i=(Shurvect.rep).tab_[k].imag();
    }
  
  for(int k=0;k<M.nbcol();k++)
    {
      stockW[k].r=M(k,k).real();
      stockW[k].i=M(k,k).imag();
    }
  //cout<<"icisel "<<i<<endl;
  for(int k=1;k<=M.nbcol();k++)
    {
      flag+=(select[k]=(cond[k]==-i)?1l:0l);
      //  cout<<"select "<<k<<" "<<select[k]<<endl;
    }
      //on fait la selection
  if (flag)
    {
      //  cout<<"avant reorde"<<M<<endl;
      ztrsen_("B","V",select,&N,tmpM,&LDT,tmpQ,
          &LDQ,W,&m,&S,&SEP,WORK,&LWORK,&info);
      for(int k=0;k<N*N;k++)
    (M.rep).tab_[k]=complex<double>(tmpM[k].r,tmpM[k].i);
      for(int k=0;k<N*N;k++)
    (Shurvect.rep).tab_[k]=complex<double>(tmpQ[k].r,tmpQ[k].i);
      //cout<<"apres reord"<<M<<endl;
      //donc maintenant on a reordonne les vp du cluster i 
      //il faut maintenant reconnaitre qui est qui pour 
      //remettre a jour cond 
      for(int k=0;k<N;k++)
    {
      //on prends pour valeur celle qui minimise la difference
      double mindiff=abs(
                 (W[0].r-stockW[k].r)*(W[0].r-stockW[k].r)
                 +(W[0].i-stockW[k].i)*(W[0].i-stockW[k].i)
                 );
      int stockpos=0;
      for(int l=0;l<N;l++)
        if(abs((W[l].r-stockW[k].r)*(W[l].r-stockW[k].r)+
           (W[l].i-stockW[k].i)*(W[l].i-stockW[k].i)<mindiff)
           )
          {
        mindiff=abs((W[l].r-stockW[k].r)*(W[l].r-stockW[k].r)
                +(W[l].i-stockW[k].i)*(W[l].i-stockW[k].i)
                );
        stockpos=l;
          }
      stockcond[stockpos]=cond[k];
    }
      for(int k=0;k<N;k++)
    cond[k]=stockcond[k];
    }
  free(select);
  free(stockcond);
  free(WORK);
  free(stockW);
  free(W);
      
}

//
//fonction de reordenement en enbtier            
//

template<typename Mat>
void reorder(Mat &Combilin,Mat &Shurvect,Mat &invShurvect)
{
  doublereal *cond;
  int i,nbcluster;
  cond=(doublereal*)malloc(sizeof(doublereal)*Combilin.nbcol());
  for(i=0;i<Combilin.nbcol();i++)
    {
      cond[i]=cond_number(Combilin,i);
      //cout<<"coond "<<i<<" "<<cond[i]<<endl;
    }
  //cout<<"nbcluster"<<endl;
  nbcluster=getcluster(Combilin,cond);
  //cout<<"docluster"<<endl;
  for(i=0;i<nbcluster;i++)
    docluster(Combilin,Shurvect,cond,i);
  //Shurvect.transpose();
  for(i=0;i<invShurvect.nbrow();i++)
    for(int j=0;j<invShurvect.nbcol();j++)
      invShurvect(i,j)=conj(Shurvect(j,i));
  //Shurvect.transpose();
  free(cond);
}

//
//trigonlisation simultanne d'une liste de matrices
//complexe de type lapack<complex<double>>
//

template<typename Mat>
void compltrigsim(list<Mat> &lmat)
{
    //on supose que la list n'est pas vide!!!
        int sum=0;
        Mat Shurvect,Combilin,tmp;
        Mat invShurvect;
    //<ugly>
    double *tab=(double*)malloc(lmat.size()*sizeof(double));
      //<\ugly>
    typename list<Mat>::iterator iterComb=lmat.begin();
    for(int i=0;i<lmat.size();i++)
      tab[i]=rand()% 32003;
    for(int i=0;i<lmat.size();i++)
      {
        sum+=tab[i];
        //  cout<<"tab["<<i<<"] "<<tab[i]<<endl;
      }
    iterComb++;
    //cout<<"completrigsim"<<endl;
    Combilin=lmat.front();
    for(int i=0;iterComb!=lmat.end();iterComb++,i++)
      {
        Combilin+=tab[i]*(*iterComb);
      }
    Combilin/=sum;
    //cout<<Combilin<<endl;
    Shurvect=wraplapack_zgees(Combilin);
    Shurvect.transpose();
        invShurvect=Shurvect;
        Shurvect.transpose();
    //tmp=invShurvect*Combilin*Shurvect;
    //Combilin=tmp;
    //cout<<"apres trigo "<<endl<<Combilin<<endl;
    for(int i=0;i<invShurvect.nbrow();i++)
        for(int j=0;j<invShurvect.nbcol();j++)
            invShurvect(i,j)=conj(invShurvect(i,j));
    //cout<<"reorder"<<endl;
    tmp=invShurvect*Combilin*Shurvect;
    Combilin=tmp;
    //cout<<"apres trigo "<<endl<<Combilin<<endl;
    reorder(Combilin,Shurvect,invShurvect);
    
    typename list<Mat>::iterator iter;
        for(iter=lmat.begin();iter!=lmat.end();iter++)
    {
        Mat tmp=invShurvect;
        tmp*=(*iter);
        tmp*=Shurvect;
        *iter=tmp;
    }
}

//
//resolution une fois le boulot fait!!
//on a plus qu'a lire sur la diagonale
//

template<typename Mat,typename Complex_Mat>
void solvefromcompltrigsim(const list<Mat> &lmat,Complex_Mat &res)
{
        int increment,i;
        typename list<Mat>::const_iterator iter,deb=lmat.begin();
        for(i=0;i<deb->nbrow();i++)
        {
          int comp=0;
           for(iter=lmat.begin();iter!=lmat.end();iter++)
              {
                res(i,comp,(*iter)(i,i));
                comp++;
              }
          }
}
//
//la meme mais en reel                             
//c plus dur parce que il peut ya vavoir de blocs 2x2   
//
template<typename Mat,typename Complex_Mat>
void solvefromtrigsim(const list<Mat> &lmat,Complex_Mat &res)
{
    int increment,i;
    typename list<Mat>::const_iterator iter,deb=lmat.begin();
    for(i=0;i<deb->nbrow();)
    {
        if ((*deb)(i+1,i)+(*deb)(i,i)==(*deb)(i,i))//cas reel
        {
            int comp=0;
            increment=1;
            for(iter=lmat.begin();iter!=lmat.end();iter++)
            {
                res(i,comp)=(*iter)(i,i);
                comp++;
            }
        }
        else //cas complex
        {
            int comp=0;
            increment=2;
            Complex_Mat tmp(2,2),tmpShu(2,2),invtmpShu(2,2);
            tmp(0,0)=(*deb)(i,i);
            tmp(1,0)=(*deb)(i+1,i);
            tmp(1,1)=(*deb)(i+1,i+1);
            tmp(0,1)=(*deb)(i,i+1);
            tmpShu=wraplapack_zgees(tmp);
            tmpShu.transpose();
            invtmpShu=tmpShu;
            tmpShu.transpose();
            invtmpShu(0,0)=conj(invtmpShu(0,0));
            invtmpShu(1,0)=conj(invtmpShu(1,0));
            invtmpShu(1,1)=conj(invtmpShu(1,1));
            invtmpShu(0,1)=conj(invtmpShu(0,1));
            for(iter=lmat.begin();iter!=lmat.end();iter++)
            {
                 tmp(0,0)=(*iter)(i,i); 
                 tmp(1,0)=(*iter)(i+1,i); 
                 tmp(1,1)=(*iter)(i+1,i+1);
                 tmp(0,1)=(*iter)(i,i+1);
                 tmp=invtmpShu*tmp*tmpShu;
                 res(i,comp)=tmp(0,0);
                 res(i+1,comp)=tmp(1,1);
                 comp++;
            }
        }
        i+=increment;
    }
}

template<typename Mat>
Mat wraplapack_zgebal(Mat M)
{
      typedef typename Mat::rep_t rept;
      float toto;
      rept m=M.rep;
      Mat Mres(m.nbrow(),m.nbrow());
      int info,n,l,ilo,ihi;
      unsigned int i;
      doublecomplex* tab =
    (doublecomplex*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex));
      doublereal *scale =
    (doublereal*)calloc(m.nbcol_,sizeof(doublereal));
      char *JOBS="B";
      
      n=m.nbrow_;l=m.nbrow_;
      
      for(i=0;i<m.nbcol_*m.nbrow_;i++) {
    tab[i].r=m.tab_[i].real();
    tab[i].i=m.tab_[i].imag();
    
    //cout<<"tab["<<2*i<<"]";
    //cout<<tab[i].r<<"    "<<tab[i].i<<endl;
    toto=tab[i].r;
    //cout<<"tentative toto"<<toto<<endl;
      }
      //cout<<"avant lapack"<<endl;
      //cout<<"int "<<sizeof(int)<<endl<<"int *"<<sizeof(int*)<<endl;
      //cout<<"fun "<<sizeof(int (*)(void *))<<endl<<"DOUBLE "<<sizeof(double)<<endl;
      zgebal_(JOBS,&n,tab,&l,&ilo,&ihi,scale,&info);
      // cout<<"apres lapack"<<tab<<" info "<<info<<" ilo "<<ilo<<" ihi "<<ihi<<endl;
      (Mres.rep).reserve(n,n);//matrice carree exclusivement!!
      for(unsigned int i=0;i<m.nbcol_*m.nbrow_;i++) 
        {
            (Mres.rep).tab_[i]=complex<double>(tab[i].r,tab[i].i);
        //cout<<"tab["<<2*i<<"]";
        cout<<tab[i].r<<"    "<<tab[i].i<<endl;
        //  cout<<"scale "<<scale[i/n]<<endl;
        }
      //      cout<<"info "<<info<<endl;
      //cout<<Mres<<endl<<endl;
      

      free(tab);
      free(scale);
      return Mres;
}

//
//resoud le system en effectuant une decomposition de Schgur simultannee 
//telle que decrite dans Corless Trager issac 97 
//

template<typename lmat,typename typres>
void GetRootMat(lmat &lres,typres &res)
{
  compltrigsim(lres);
  {
    typres tmpres((lres.front()).nbrow(),lres.size());
    solvefromcompltrigsim(lres,tmpres);
    res=tmpres;
  }
}

template<typename typres,typename Base>
void realsolve2(const list<harewell<QQ> > &lmat,const Base &b,typres &res)
{
  int i;
  list<MatDense<lapack<complex<double > > > >lres;
  
  for(list<harewell<QQ> >::const_iterator iter=lmat.begin()
    ;iter!=lmat.end();iter++)
    {
       MatDense<lapack<complex<double> > > tmp(iter->nrow,iter->ncol);
      (tmp.rep).reserve(iter->nrow,iter->ncol);
      for(int j=0;j<iter->nrow;j++)
    for(int k=0;k<iter->ncol;k++)
      {
        tmp(j,k)=mpq_get_d(&(*iter)(j,k).rep());
        //tmp(i,k).imag()=0;
        //fortement sous optimal
      }
      lres.push_back(tmp);
    }
  GetRootMat(lres,res);
}

template<typename typres, typename Base>
void realsolve2(const list<harewell<RR> > &lmat,const Base &b,typres &res)
{
  int i;
  list<MatDense<lapack<complex<double > > > >lres;
  
  for(list<harewell<RR> >::const_iterator iter=lmat.begin()
    ;iter!=lmat.end();iter++)
    {
       MatDense<lapack<complex<double> > > tmp(iter->nrow,iter->ncol);
      (tmp.rep).reserve(iter->nrow,iter->ncol);
      for(int j=0;j<iter->nrow;j++)
    for(int k=0;k<iter->ncol;k++)
      {
        tmp(j,k)=mpf_get_d((*iter)(j,k).rep());
        //tmp(i,k).imag()=0;
        //fortement sous optimal
      }
      lres.push_back(tmp);
    }
  GetRootMat(lres,res);
}

template<typename typres, 
  typename typlmat,typename Base>
void realsolve3(const typlmat &lmat,const Base &b,typres &res)
{
  int i;
  list<MatDense<lapack<complex<double > > > >lres;
  
  for(typename typlmat::const_iterator iter=lmat.begin()
    ;iter!=lmat.end();iter++)
    {
       MatDense<lapack<complex<double> > > tmp(iter->nrow,iter->ncol);
      (tmp.rep).reserve(iter->nrow,iter->ncol);
      for(int j=0;j<iter->nrow;j++)
    for(int k=0;k<iter->ncol;k++)
      {
        tmp(j,k)=mpq_get_d(&(*iter)(j,k).rep());
        //tmp(i,k).imag()=0;
        //fortement sous optimal
      }
      lres.push_back(tmp);
    }
  GetRootMat(lres,res);
}