resolve2.hpp

Go to the documentation of this file.

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

#include "borderbasix/arithm/Scl_mpf.hpp"
#include "linalg/clapack.h"
#include<complex>
#include"linpkg/Lapack.H"
#include"linpkg/Eigen.H"
#include<stdlib.h>
#include<cstdlib>
//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_(...);
int 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);

template<typename T>
int zgees_(const char *jobvs,const char *sort, 
       logical (*select)(doublecomplex<T>*),tla_int *n,
        doublecomplex<T> *a, tla_int *lda, tla_int *sdim, doublecomplex<T> *w,
        doublecomplex<T> *vs, tla_int *ldvs, doublecomplex<T> *work,
       tla_int *lwork, T *rwork, logical *bwork, tla_int *info);
template<typename T>
int ztrsen_(const char *job, const char *compq, logical *select, tla_int 
        *n, doublecomplex<T> *t, tla_int *ldt, doublecomplex<T> *q, 
        tla_int *ldq, doublecomplex<T> *w, tla_int *m, T *s, T *sep, 
        doublecomplex<T> *work, tla_int *lwork, 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"{
template<typename T>
long int select_true(doublecomplex<T> *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 verifier la fonction dznrm2!!!
//Elle est pas compilee de maniere adequate dans la distrib semir
template<typename Mat>
Mat wraplapack_zgees(Mat M)
{
      typedef typename Mat::rep_t rept;
      typedef typename Mat::coeff_t::value_type T;
      rept m=M.rep;
      logical (*SELECT)(doublecomplex<T> *)=select_true;
      doublecomplex<T> *w=(doublecomplex<T> *)
    malloc(m.nbrow_*sizeof(doublecomplex<T>));
      T *rwork=(T*)malloc(m.nbrow_*sizeof(T));
      logical *bwork=(logical*)malloc(m.nbrow_*sizeof(logical));
      Mat Mres(m.nbrow_,m.nbrow_);
      doublecomplex<T>* work=
    (doublecomplex<T>*)malloc(10*sizeof(doublecomplex<T>)*m.nbrow_);
      tla_int info,n,l,lres,lwork,sdim;
      unsigned int i;
      doublecomplex<T>* tab =
    (doublecomplex<T>*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex<T>));
      doublecomplex<T>* tabres =
    (doublecomplex<T>*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex<T>));
      const 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<T>(tab[i].r,tab[i].i);
        //toto=*(((T*)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<T>(tabres[i].r,tabres[i].i);
            }
      free(w);
      free(rwork);
      free(work);
      free(bwork);
      free(tab);
      free(tabres);
      return Mres;
}

#if 0
template<>
MatDense<lapack<complex<mpf_class> > >
 wraplapack_zgees(MatDense<lapack<complex<mpf_class> > >M)
{
  typedef MatDense<lapack<complex<mpf_class> > > Mat;
  typedef lapack<complex<mpf_class> > rept;
  typedef mpf_class T;
  rept m=M.rep;
  logical (*SELECT)(doublecomplex<T> *)=select_true;
  doublecomplex<T> *w=(doublecomplex<T> *)
    malloc(m.nbrow_*sizeof(doublecomplex<T>));
  T *rwork=(T*)malloc(m.nbrow_*sizeof(T));
  logical *bwork=(logical*)malloc(m.nbrow_*sizeof(logical));
  Mat Mres(m.nbrow_,m.nbrow_);
  doublecomplex<T>* work=
    (doublecomplex<T>*)new doublecomplex<T>[10*m.nbrow_];
  tla_int info,n,l,lres,lwork,sdim;
  unsigned int i;
  doublecomplex<T>* tab =
    (doublecomplex<T>*) new doublecomplex<T>[m.nbcol_*m.nbrow_];
  doublecomplex<T>* tabres = new doublecomplex<T>[m.nbcol_*m.nbrow_];
  //    (doublecomplex<T>*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex<T>));
  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<T>(tab[i].r,tab[i].i);
      //toto=*(((T*)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<T>(tabres[i].r,tabres[i].i);
    }
  free(w);
  free(rwork);
  delete[] work;
  free(bwork);
  delete[] tab;
  delete[] tabres;
  return Mres;
}
#endif

//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,typename coeff>
void wraplapack_ztrsen(Mat &M,logical * select,coeff *res)
{
  //typedef typename Mat::coeff_t::value_type coeff;
  const char *JOBS="E",*COMPQ="N";
  tla_int N=M.nbcol(),LDT=M.nbcol();
  doublecomplex<coeff> *T=(doublecomplex<coeff>*)
    malloc(M.nbcol()*M.nbcol()*sizeof(doublecomplex<coeff>));
  doublecomplex<coeff> *Q=NULL;
  doublecomplex<coeff> *W=(doublecomplex<coeff>*)
    malloc(M.nbcol()*sizeof(doublecomplex<coeff>));
  tla_int m=1,LDQ=1;
  coeff S=1,SEP=2;
  doublecomplex<coeff> *WORK=(doublecomplex<coeff>*)
    malloc(2*N*N*sizeof(doublecomplex<coeff>));
  tla_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 clusterise de 
//chacunes des valeurs propres de la matrice
  
template<typename T,typename Mat>
T cond_number(Mat &Combilin,int i)
{
  //typedef typename Mat::coeff_t::value_type T;
  T res;
  logical *select;
  select=(logical*)calloc(Combilin.nbcol(),sizeof(logical));
  select[i]=1;
  wraplapack_ztrsen(Combilin,select,&res);
  free(select);
 return res;
}
//
//creation des clusters 
//utilisation de l'heuritique de Stetter!!
//

template<typename Mat,typename typcond>
int getcluster(Mat &M,typcond *cond)
{
  typedef typename Mat::coeff_t::value_type T;
  // typedef MYTYPE T;
  T EPSILON=1e-16;
  doublecomplex<T> *truncnum=(doublecomplex<T>*)
    malloc(M.nbcol()*sizeof(doublecomplex<T>));
  int *digits=(int*)malloc(M.nbcol()*sizeof(int));
  T normM=0;
  char tampon1[60],tampon2[60],*rien;
  int comp=0;
  for(unsigned i=0;i<M.nbcol();i++)
    normM=max(sqrt((M(i,i)*conj(M(i,i))).real()),normM);
  if(normM<EPSILON)
    {
      for(unsigned i=0;i<M.nbcol();i++)
        cond[i]=0;
      //cout<<"ca foire ici"<<endl;
        return M.nbcol();
    }
  //cout<<"M.nbcol()"<<M.nbcol()<<endl;
  for(unsigned i=0;i<M.nbcol();i++)
    {
      //cout<<"cond[i]"<<cond[i]<<endl;
      //cout<<"normM"<<normM<<endl;
      //cout<<EPSILON*normM<<endl;
      //cout<<1.0/cond[i]<<endl;
      T lambda(1.0/cond[i]*EPSILON*normM);
      //cout<<"lada"<<lambda<<endl;
      T nbdigits=(1.0/lambda);
      int k;
      for(k=1;nbdigits>1;k++,nbdigits/=10);//
      //cout<<"nbdigits "<<nbdigits<<" K "<<k<<endl;;
      nbdigits=k;
      digits[i]=k;
      //rajout des tests de taille apres
      sprintf(tampon1,"%%.%dLg\0",(k>16)?16:k);
      //et pas snprintf parce que ca compile pas sur alpha.
      //      cout<<tampon1<<endl;
      sprintf(tampon2,tampon1,M(i,i).real());
      //et pas snprintf parce que ca compile pas sur alpha.
      //cout<<tampon2<<endl;
      truncnum[i].r=(T)(strtod(tampon2,&rien));
      sprintf(tampon2,tampon1,M(i,i).imag());
      //et pas snprintf parce que ca compile pas sur alpha.
      //cout<<tampon2<<endl;
      truncnum[i].i=(T)strtod(tampon2,&rien);
    }
  for(int i=0;i<M.nbcol();i++)
    {
      //if(cond[i]>=0)
      //{
    //  cout<<"toto"<<endl;
      cond[i]=-i;
          cout<<"i"<<i<<"condi"<<cond[i]<<endl;
      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)*(T)val<(T)1.0)
          &&(abs(truncnum[j].i-truncnum[i].i)*(T)val<(T)1.0))
        {
          typcond stock=cond[j];
                  cout<<"stock"<<stock<<endl;
                  for(int k=0;k<M.nbcol();k++)
                  {
                      cout<<"condk"<< cond[k]<<endl;
            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);
  free(digits);
  return M.nbcol();//cluster numerotes de -1 a ...
}

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

template<typename Mat,typename typcond>
int docluster(Mat &M,Mat &Shurvect,typcond * cond,int i)
{
  //0n recherche tous les -i dans le tableau cond et apres on
  //fait un cluster avec eux......
  typedef typename Mat::coeff_t::value_type T;
  //typedef MYTYPE T;
  logical *select=(logical*)malloc(M.nbcol()*sizeof(logical));
  tla_int N=M.nbcol(),LDT=M.nbcol(),LDQ=M.nbcol();
  tla_int m,info,LWORK=N*N*2;
  long flag=0;
  T *stockcond=(T*)malloc(N*sizeof(T));
  doublecomplex<T> *WORK=(doublecomplex<T>*)
    malloc(sizeof(doublecomplex<T>)*LWORK);
  T S,SEP;
  doublecomplex<T> *stockW=(doublecomplex<T>*)
    malloc(N*sizeof(doublecomplex<T>));
  doublecomplex<T> *W=(doublecomplex<T>*)malloc(N*sizeof(doublecomplex<T>));
  doublecomplex<T> *tmpM=(doublecomplex<T>*)
    malloc(N*N*sizeof(doublecomplex<T>));
  doublecomplex<T> *tmpQ=(doublecomplex<T>*)
    malloc(N*N*sizeof(doublecomplex<T>));
  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();
    }
  //cout<<"1"<<endl;
  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=0;k<M.nbcol();k++)
    {
     // cout<<"i"<<i<<endl;
     // cout<<"M"<<M<<endl;
    //  cout<<"select "<<k<<" "<<select[k]<<endl;
     // cout<<"cond "<<k<<" "<<cond[k]<<endl;
      select[k]=((cond[k]==-i)?1:0);
      flag+=select[k];
       //cout<<"select "<<k<<" "<<select[k]<<endl;
   }
      //on fait la selection
  //cout<<"flag"<<flag<<endl;
//  cout<<"2"<<endl;
  if (flag>1)
    {
      //      cout<<"on fait quelquechose reorde"<<M<<endl;
      const char *JOBS="B",*COMPQ="V";
      ztrsen_(JOBS,COMPQ,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<T>(tmpM[k].r,tmpM[k].i);
      for(int k=0;k<N*N;k++)
    (Shurvect.rep).tab_[k]=complex<T>(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
      char * posbusy=(char*)malloc(N);
      for(int i=0;i<N;i++)
    posbusy[i]=0;
      for(int k=0;k<N;k++)
    {
      //on prends pour valeur celle qui minimise la difference
      int stockpos=0;
      for(;posbusy[stockpos];stockpos++);
      T mindiff=abs(
            (W[stockpos].r-stockW[k].r)*(W[stockpos].r-stockW[k].r)
            +(W[stockpos].i-stockW[k].i)*(W[stockpos].i-stockW[k].i)
                 );
      for(int l=0;l<N;l++)
        if((!posbusy[l])&&(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;
          }
      posbusy[stockpos]=1;
      stockcond[stockpos]=(T)cond[k];
    }
      free(posbusy);
      for(int k=0;k<N;k++)
    cond[k]=stockcond[k];
    }
 // cout<<"cool"<<endl;
  free(select);
  free(stockcond);
  free(WORK);
  free(stockW);
  free(W);
  free(tmpM);
  free(tmpQ);
  return flag;
}

//
//fonction de reordenement en enbtier            
//

template<typename T,typename Mat>
void reorder(Mat &Combilin,Mat &Shurvect,Mat &invShurvect,
         list<int> &size_clust)
{
  T *cond;
  int i,nbcluster;
  cond=new T[Combilin.nbcol()];
    //(doublereal*)malloc(sizeof(doublereal)*Combilin.nbcol());
  for(i=0;i<(int)Combilin.nbcol();i++)
    {
      cond[i]=cond_number<T>(Combilin,i);
      //cout<<"coond "<<i<<" "<<cond[i]<<endl;
    }

//  cout<<"cond"<< cond<<endl;
  nbcluster=getcluster(Combilin,cond);
 // cout<<"nbcluster"<< nbcluster<<endl;
 // cout<<"combilin"<<Combilin<<endl;

  for(i=0;i<nbcluster;i++)
    {
//      cout<<"i"<<i<<endl;
//      cout<<"nbcluster"<<nbcluster<<endl;
      int tmp_size=docluster(Combilin,Shurvect,cond,i);
     // cout<<"tmp_size"<<tmp_size<<endl;
      if (tmp_size!=0)
    {
      if(tmp_size>1)
        size_clust.push_back(tmp_size);
      else
        size_clust.push_front(tmp_size);
    }
    //  cout<<"final"<<endl;
    }
//  cout<<"sortie"<<endl;
//  cout<<"d1clu"<<endl;
  //cout<<"nb cluster "<<size_clust.size()<<endl;
 //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);
  delete[] cond;
}

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

template<typename Mat>
void compltrigsim(list<Mat> &lmat,list<int> &size_clust)
{
    //on supose que la list n'est pas vide!!!
  typedef typename Mat::coeff_t::value_type T;      
  double 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(unsigned i=0;i<lmat.size();i++)
    tab[i]=rand()% 32003+1;
  for(unsigned i=0;i<lmat.size();i++)
    {
      sum+=tab[i];
      //  cout<<"tab["<<i<<"] "<<tab[i]<<endl;
    }
  iterComb++;
//  cout<<"completrigsim"<<endl;
  Combilin=lmat.front();
  for(unsigned i=0;iterComb!=lmat.end();iterComb++,i++)
    {
      Combilin+=T(tab[i])*(*iterComb);
    }
  Combilin/=T(sum);
  //cout<<"Combilin"<<endl;
  Shurvect=wraplapack_zgees(Combilin);
//  cout<<"apres Combilin"<<endl;
  Shurvect.transpose();
  invShurvect=Shurvect;
  Shurvect.transpose();
  //tmp=invShurvect*Combilin*Shurvect;
  //Combilin=tmp;
 // cout<<"apres trigo "<<Combilin<<endl;
  for(unsigned i=0;i<invShurvect.nbrow();i++)
    for(unsigned j=0;j<invShurvect.nbcol();j++)
      invShurvect(i,j)=conj(invShurvect(i,j));
//  cout<<"reorder"<<endl;
  tmp=invShurvect*Combilin*Shurvect;
  Combilin=tmp;
  //cout<<"apres trigo2 "<<Combilin<<endl;
  reorder<T>(Combilin,Shurvect,invShurvect,size_clust);
 // cout<<"apres trigo3 "<<endl;
  typename list<Mat>::iterator iter;
  for(iter=lmat.begin();iter!=lmat.end();iter++)
    {
      Mat tmp=invShurvect;
      tmp*=(*iter);
      tmp*=Shurvect;
      *iter=tmp;
    }
 // cout<<"apres trigo4 "<<endl;
}

//
//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,const list<int> &size_clust,
               Complex_Mat &res)
{
  int increment,i;
    typename list<Mat>::const_iterator iter,deb=lmat.begin();
        typename list<int>::const_reverse_iterator
          iter_size=size_clust.rbegin();
//               for( iter=lmat.begin();iter!=lmat.end();iter++)
//          cout<<*iter<<endl<<"--------------"<<endl;
//        cout<<"deb->nbrow()"<<deb->nbrow()<<endl;
//        for( typename list<int>::const_iterator i=size_clust.begin();i!=size_clust.end();i++)
//            cout<<*i<<" toto "<<endl;
        for(int i=0;i<deb->nbrow();)
        {
          int comp=0;
          cout<<"Solution: ";
          for(iter=lmat.begin();iter!=lmat.end();iter++)
              {
             //  cout<<"deux fois"<<endl;
               //cout<<"*iter_size"<<*iter_size<<endl;
        typename list<Mat>::value_type::coeff_t value(0) ,tmp(*iter_size);
                // cout<<"bug"<<endl;
        //on calcule la valeur

        // Segmentation fault dans cette boucle pour des systemes de solutions 0
        for(int j=0;j<*iter_size;j++) {
//                    cout<<"no puedes"<<endl;
          value+=(*iter)(i+j,i+j);
                  // cout << "i " << i << " j " << j << " *iter.size() " << (*iter).nbcol() << endl;
        }
                //cout << "bug 2" << endl;
                //cout<<"value"<<value<<endl;
        value/=tmp;
        for(int j=i;j<i+*iter_size;j++) {
          res(j,comp,value);
                  cout <<" "<< res(j,comp,value);
        }
                comp++;
                increment=*iter_size;
              }
           std::cout<<std::endl;
          // cout<<"taille cluster "<<*iter_size<<endl;
           //i+=increment;
           i+=*iter_size++;
          }
}
//
//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;
      typedef typename Mat::coeff_t::value_type T;
      //      typedef long double T;
      float toto;
      rept m=M.rep;
      Mat Mres(m.nbrow(),m.nbrow());
      int info,n,l,ilo,ihi;
      unsigned int i;
      doublecomplex<T>* tab =
    (doublecomplex<T>*)calloc(m.nbcol_*m.nbrow_,sizeof(doublecomplex<T>));
      T *scale =
    (T*)calloc(m.nbcol_,sizeof(T));
      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)
{
  list<int> size_clust;
  compltrigsim(lres,size_clust);
 // cout<<"fincomple--------------------------------"<<endl;
 for(typename lmat::iterator iter=lres.begin();iter!=lres.end();iter++)

  {
     cout<<*iter<<endl<<endl;
    typres tmpres((lres.front()).nbrow(),lres.size());
   // cout<<"si"<<endl;
    solvefromcompltrigsim(lres,size_clust,tmpres);
   // cout<<"no"<<endl;

    res=tmpres;

  }
}

template<typename typres,typename Base>
void realsolve2(const list<harewell<QQ> > &lmat,const Base &b,typres &res)
{
  //int i;
  list<MatDense<lapack<typename typres::coeff_t > > >lres;
  
  cout<<"RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR"<<endl;
  for(list<harewell<QQ> >::const_iterator iter=lmat.begin()
    ;iter!=lmat.end();iter++)
    {
      MatDense<lapack<typename typres::coeff_t > > 
    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(j,k)=&(*iter)(j,k).rep();
        //tmp(i,k).imag()=0;
        //fortement sous optimal
      }
      lres.push_back(tmp);
    }
  cout<<"RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR"<<endl;
  GetRootMat(lres,res);
}

template<typename typres, typename Base>
void realsolve2(const list<harewell<RR> > &lmat
        ,const Base &b,typres &res)
{
    //  int i;
    // cout<<"entra aqui"<<endl;
    typedef typename typres::coeff_t T;//Ca c deja complex
    list<MatDense<lapack<typename typres::coeff_t > > > lres;

    for(list<harewell<RR> >::const_iterator iter=lmat.begin()
        ;iter!=lmat.end();iter++)
    {
        MatDense<lapack<typename typres::coeff_t> >
                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++)
            {
                //char *toto=NULL;
                char _toto [256]; char* toto= & (_toto[0]);
                mp_exp_t expon;
                int sign=1;
                toto= mpf_get_str(toto,&expon,10,50,(*iter)(j,k).rep());
                tmp(j,k)=(T)(0);
                if (toto[0]=='-')
                {
                    sign=-1;
                    toto[0]='0';
                }
                int l;
                for(l=0;toto[l]!='\0';l++)
                {
                    tmp(j,k)=(T)(10)*tmp(j,k)+(T)(toto[l]-'0');
                    //cout<<toto<<endl;
                    //cout<<tmp(j,k)<<endl;
                }
                // cout<<-l+1+expon-(sign+1)/2<<endl;
                tmp(j,k)*=(T)(sign*pow(10.0,(int)(-l+1+expon-((sign+1)/2))));

                //             cout<<"toto "<<toto<<endl;
                //             cout<<"tmp(j,k) a la fin "<<tmp(j,k)<<" "
                //             <<mpf_get_d((*iter)(j,k).rep())<<endl;;
                //free(toto);
                //tmp(i,k).imag()=0;
                //fortement sous optimal
            }
        lres.push_back(tmp);
    }
    GetRootMat(lres,res);
}
template<typename typcoeff,typename typres, typename Base>
void realsolve2(const list<harewell<typcoeff> > &lmat
        ,const Base &b,typres &res)
{

   int i;
  list<MatDense<lapack<typename typres::coeff_t > > >lres;
  
  for(typename list<harewell<typcoeff> >::const_iterator iter=lmat.begin()
    ;iter!=lmat.end();iter++)
    {
      MatDense<lapack<typename typres::coeff_t > > 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)= (*iter)(j,k);
        //fortement sous optimal
      }
      lres.push_back(tmp);
    }
  GetRootMat(lres,res);
}

template<int p,typename typrep,typename typres, typename Base>
void realsolve2(const list<harewell<Z<p,typrep> > > &lmat
        ,const Base &b,typres &res)
{
  res.ncol=0;
  res.nrow=0;
}
template<typename typres, 
  typename typlmat,typename Base>
void realsolve3(const typlmat &lmat,const Base &b,typres &res)
{
  int i;
  typedef typename typlmat::value_type::coeff_t::value_type T;
  list<MatDense<lapack<complex<T > > > >lres;
  
  for(typename typlmat::const_iterator iter=lmat.begin()
    ;iter!=lmat.end();iter++)
    {
       MatDense<lapack<complex<T> > > 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);
}
#undef MYTYPE