generkeropt.hpp

Go to the documentation of this file.

#include <map>
#include <cstdlib>
#include <borderbasix/mdebug.hpp>
#include <borderbasix/mpoly.hpp>
#include <borderbasix/mpoly/matrix_of.hpp>
#include <borderbasix/solversdp_sdpa.h>
#include <borderbasix/solversdp_gmp.h>
#include <borderbasix/solversdp_csdp.h>
#include <borderbasix/solversdp_mosek.h>
//#include <mosek.h>


#include "tensor-decomposition/Decomposer.hpp"
//#include "../arithm/Scl.hpp"
#include "../arithm/Scl_mpf.hpp"

template<typename typP, typename typdump,typename Base, typename typserv,typename coeff>
void generickernelNEW(typP &newpol, int k,const typdump &dump, Base &b, typserv &serv, typP F,coeff *&sol2,typP rest, coeff *&solfinal,int solver, int &method, int &flat, string param, double &tresh)
{
    using namespace coinor;
    typedef typename typP::value_type::coeff_t Scalar;
    typedef typename typP::value_type::order_t ord;
    typedef Monom<double, dynamicexp<'x'> > Monomial;
    typedef Dlex<Monomial> Ordering;
    typedef MPoly<std::vector<Monomial>, Ordering > Polynomial;
    //typedef MPoly<std::vector<mon>,Dlex<mon> > mpol;
    //mdebug()<<"entra en mi generic kernel";
    //
    //Allocation du problem et des
    //contraintes

    struct sparseblockv<COEFF> *blockptr2;
    double *bb,*sol1;//*sol3;

    int i,aveccont=0,numsol=0,solst;

    struct constraintmatrixv<COEFF> *constraints2;

    mdebug()<<"tresh dedans de generker"<<tresh;
    //
    //On introduit un variable par monome de degre d en dehors de B
    //On va contruire l'operateur de de Hankel associé

    //
    //On commence par enumerer les monomes de degre plus petit que k
    //
    //
    typP stockprod,stockproj,stockproj2,stockproj2b,listmin,listminmon,listrest,listmon,listLambda;//dans stockprod on a des polynomes
    //
    //
    //
    int degmax,nbconstraint=0;
    list<typename typserv::monom_t> stockksur2,listmonLambda;
    list<Monomial> listmonLambda1;
    map<mon,int,ord> mapmonbase;
    map<Poly,Poly> mapmonproj;
    list<Monomial> monomesadd1;
    int grado=k;

    //
    //ICI on cree un index pour le monomes de degre k

    enumeratemon(stockksur2,k/2,b,dump);

    enumeratemoninv(listmonLambda,k/2,b,dump);
    enumeratemoninv(listmonLambda1,k,b,dump);


    //
    //pour le suite il nous faut les monomes de B de degre <=k/2
    //

    for(typename list<typename typserv::monom_t>::iterator
        iter=stockksur2.begin(); iter!=stockksur2.end();iter++)
    {
        if(!IsinB(*iter,b)) stockksur2.erase(iter--);
    }

    for(typename list<typename typserv::monom_t>::iterator
        iter=listmonLambda.begin(); iter!=listmonLambda.end();iter++)
    {

        if(!IsinB(*iter,b)) listmonLambda.erase(iter--);
    }
    for(typename list<Monomial>::iterator
        iter=listmonLambda1.begin(); iter!=listmonLambda1.end();iter++)
    {

        if(!IsinB(*iter,b))
        {
            mdebug()<<*iter;
            monomesadd1.push_back(*iter);
            listmonLambda1.erase(iter--);

        }
        if (*iter==1) listmonLambda1.erase(iter--);
    }
#ifdef DEBUGSDP
    mdebug()<<"monome en degre k/2 ";
    for(typename list<typename typserv::monom_t>::iterator
        iter=stockksur2.begin(); iter!=stockksur2.end();iter++)
        mdebug()<<*iter<<" ";

    mdebug();
#endif
    for(typename list<typename typserv::monom_t>::iterator
        iter=stockksur2.begin(); iter!=stockksur2.end();iter++)
        degmax=2*Degree(*iter);
    if (degmax<k) k=degmax;
    //mdebug()<<"paso3";

    int dim=stockksur2.size();
    int dimrest=rest.size();
    int numvar=b.nbvar();
    int r1=0,u,file,contmon,sizerest=0,relat,s,maxdegrest=0,nummonineq;
    //COEFF *listdegrest;
    int listsizerest[dimrest];
    int listdegrest[dimrest];
    int listnumeq[dimrest];


    for(typename Base::iterator itb=b.begin();itb!=b.end();itb++)
        for(int i=0;i<itb->taille2;i++)

            if (Degree(itb->refuse[i])==1 ) numvar--;

    for(typename typP::iterator iter=rest.begin();iter!=rest.end();
        iter++)
        if (Degree(*iter)>=maxdegrest) maxdegrest=Degree(*iter);

    mdebug()<<"maxdegrest"<<maxdegrest;
    //
    //  On  genere les formes normales des produits de monomes de
    // de tous les monomes de   degre <=k/2 dez B.
    //

    computeprod(stockprod,stockksur2,b,dump);

    computeprod(listLambda,listmonLambda,b,dump);

//    for(typename list<Monomial>::iterator i=listmonLambda1.begin();
//        i!=listmonLambda1.end();i++)
//    {
//        for(typename list<Monomial>::iterator j=i;j!=listmonLambda1.end();j++)
//            listLambda1.push_back(*i * *j);
//    }
//    listLambda1.unique();
    //listmonLambda1.sort(Dlex<Monomial>());

    mdebug()<<" MonNoproj";
    for(typename list<Monomial>::iterator
        iter=listmonLambda1.begin(); iter!=listmonLambda1.end();iter++)
        mdebug()<<*iter<<" ";

    mdebug()<<" Monproj";
    for(typename list<Monomial>::iterator
        iter1=monomesadd1.begin(); iter1!=monomesadd1.end();iter1++)
        mdebug()<<*iter1<<" ";


    //computeprodMonomial(listLambda1,listmonLambda,b,dump);



    stockproj=stockprod;

    //on calcule la forme normale de la liste de monomes de order k
//    stockproj2=stockprod;
//    proj(stockproj2,dump,b,serv);

    for(typename typP::iterator iter=F.begin();iter!=F.end();
        iter++)
        mdebug()<<*iter<<" ";
    //on calcule la forme normale de la function f
    proj(F,dump,b,serv);

    //on calcule la forme normale des inegalites
    proj(rest,dump,b,serv);




    //on calcule les polynomes qui sont parties de la matrices de localization associes a
    //les inegalites
    for(typename typP::iterator iter=rest.begin();iter!=rest.end();
        iter++,r1++)
    {


        listdegrest[r1]=Degree(*iter);
        u=floor(Degree(*iter)/2);
        if (Degree(*iter)%2!=0) u++;
        s=k/2-u;
        sizerest=factorial(s+numvar)/(factorial(s)*factorial(numvar));
        listsizerest[r1]=sizerest;
        nummonineq=0;

        if (s!=0)
        {

            relat=0;
            file=1;
            contmon=0;
            listmon.erase(listmon.begin(),listmon.end());
            for(typename list<typename typserv::monom_t>::iterator
                iters=stockksur2.begin(); iters!=stockksur2.end();iters++)

                if (contmon<sizerest)
                {

                    listrest.push_back(*iter * *iters);
                    listmon.push_back(*iters);
                    contmon++;
                    nummonineq++;

                }

            typename typP::iterator iterm=listmon.begin();

            for ( ;iterm!=listmon.end(); iterm++)
            {
                typename typP::iterator iterm2=listmon.begin();
                if (file!=sizerest)
                {
                    file++;
                    relat=file-1;
                    for(int j=1; j<file; j++)
                    {
                        if(file==2) iterm++;
                        iterm2++;
                    }


                    for (;iterm2!=listmon.end(); iterm2++)
                        if (relat!=sizerest)
                        {
                            listrest.push_back(*iter * *iterm * *iterm2);
                            relat++;
                            nummonineq++;
                        }

                }

            }



        }

        else{

            listrest.push_back(*iter);
            nummonineq++;
        }
        listnumeq[r1]=nummonineq;
    }

    //on calcule les formes normales de cette liste de polynomes (inegalites)
    proj(listrest,dump,b,serv);
    selectnozero(listrest);


#ifdef DEBUGSDP
    mdebug()<<"produits ";
    for(typename typP::iterator iter=stockprod.begin();iter!=stockprod.end();
        iter++)
        mdebug()<<*iter<<" ";
    mdebug();
    mdebug()<<"proj2 ";
    mdebug()<<"stockprod.size()"<<stockprod.size();
    for(typename typP::iterator iter=stockproj2.begin();iter!=stockproj2.end();
        iter++)
        mdebug()<<*iter<<" ";
    mdebug();
    mdebug()<<"Funcionproj";
    for(typename typP::iterator iter=F.begin();iter!=F.end();
        iter++)
        mdebug()<<*iter<<" ";
    mdebug();
    mdebug()<<"rest";
    for(typename typP::iterator iter=listrest.begin();iter!=listrest.end();
        iter++)
        mdebug()<<*iter<<" ";
    mdebug();

    mdebug()<<"monomios de Lambda";
    for(typename list<typename typserv::monom_t>::iterator
        iter=listmonLambda.begin(); iter!=listmonLambda.end();iter++)
    {

        mdebug()<<*iter<<"";
        mdebug();
    }

#endif


    for(typename typP::iterator iter=listLambda.begin();iter!=listLambda.end();
        iter++)
    {
        if(!IsinB(*(iter->begin()),b)) listLambda.erase(iter--);
        if (*(iter->begin())==1) listLambda.erase(iter--);
    }

//    for(typename list<Monomial>::iterator iter=listLambda1.begin();iter!=listLambda1.end();iter++)
//    {
//        if(!IsinB(*iter,b)) listLambda1.erase(iter--);
//        if (*iter==1) listLambda1.erase(iter--);
//    }
//    for(typename list<Monomial>::iterator
//        iter=listLambda1.begin(); iter!=listLambda1.end();iter++)
//        mdebug()<<*iter<<" ";

    //
    //

    //
    //
    //
    //
    // §Ici on va avoir besoin d'une structure contenant un monome de tete
    // hors de B et sa forme normal


    typP listnormalform,stockmonproj;


    int countmonb=0,countmont=0,countmonp=0;
    for(typename typP::iterator iter=stockprod.begin();
        iter!=stockprod.end();iter++)
    {
        if(IsinB(*(iter->begin()),b))
        {
            if (*iter !=1)
            {
                listmin.push_back(*iter);
            }
            listminmon.push_back(*iter);
            nbconstraint++;
        }
        else
        {
            countmonp++;
            mdebug()<<"monproj"<<countmonp ;
            stockmonproj.push_back(*iter);
            listnormalform.push_back(*iter);
        }
    }


    int count=0;
    for(typename typP::iterator itermonb=listminmon.begin();itermonb!=listminmon.end();
        itermonb++)
    {
        count++;
        mapmonbase[*itermonb->begin()]=count;
    }
    mdebug()<<"antes de definir las constraint";


    int iterac=0;
    if (!rest.empty())
    {
        aveccont=1;
        iterac=0;
        for(int k=0; k<dimrest; k++)
            iterac+=listsizerest[k];

    }

    int dimmonap=listminmon.size();
    mdebug()<<"dim mon en b"<<dimmonap;
    mdebug()<<"dim mon qu'on proj"<<countmonp;



    int nummatrix;

    mdebug()<<"dim"<<dim;





    mdebug()<<"avant proj";
    proj(listnormalform,dump,b,serv);
    typename typP::iterator iter=listnormalform.begin();
    for(typename typP::iterator iterproj=stockmonproj.begin();
        iterproj!=stockmonproj.end();iterproj++,iter++)
        mapmonproj[*iterproj]=*iter;

    mdebug()<<"apres proj";

    solversdp* slv;



    if (solver==1) slv= new solversdp_sdpa;
    else if (solver==2) slv=new solversdp_gmp;
       else if (solver==3) slv = new solversdp_csdp;
          else if (solver==4) slv = new solversdp_mosek;



    slv->dim()=dim;
    for(unsigned k=0;k<dimrest;k++)
        slv->dim(k+1)=listsizerest[k];


    int i2_counter=1,j2_counter;

    for(typename list<typename typserv::monom_t>::iterator i=stockksur2.begin();
        i!=stockksur2.end();i++,i2_counter++)
    {
        j2_counter=i2_counter;
        for(typename list<typename typserv::monom_t>::iterator j=i;
            j!=stockksur2.end();j++,j2_counter++)
        {


            if (!IsinB(*i * *j, b))
            {


                Poly polproj;
                polproj=mapmonproj[*i * *j];
                COEFF value;

                for(typename typP::value_type::iterator iterpol=polproj.begin();
                    iterpol!=polproj.end();iterpol++)
                {
                    countmont++;
                    nummatrix=mapmonbase[*iterpol];
                    //std::cout<<nummatrix<<std::endl;
                    value=(COEFF)iterpol->GetCoeff();
                    if (nummatrix==dimmonap)
                    {
                        slv->coeff(0,0,i2_counter-1,j2_counter-1)=Scalar::castDouble(value);
                        slv->nentries(0)++;
                    }
                    else
                    {
                        slv->coeff(nummatrix,0,i2_counter-1,j2_counter-1)=Scalar::castDouble(value);
                        slv->nentries(nummatrix)++;
                    }




                }




            }
            else
            {
                countmonb++;
                countmont++;
                nummatrix=mapmonbase[*i * *j];
                if (nummatrix==dimmonap)
                {
                    slv->coeff(0,0,i2_counter-1,j2_counter-1)=1;
                    slv->nentries(0)++;
                }
                else
                {
                 slv->coeff(nummatrix,0,i2_counter-1,j2_counter-1)=1;
                 slv->nentries(nummatrix)++;
                }

            }

        }
    }




    mdebug()<<"num mon total"<<countmont;
    int  indice=0,indicej=0,eqind,indcont,numine=1,numeqmatrix;
    int nummatrixcon=1;




   // cout<<"avant ine"<<endl;
    //numeqmatrix=listrest.size()/dimrest;
    int nineq=0;
    int kk=0;
    for(typename typP::iterator itera=listrest.begin();
        itera!=listrest.end();itera++,numine++)
    {
        eqind=0;
        mdebug()<<"numine"<<numine;
        mdebug()<<"listnumeq[nineq]"<<listnumeq[nineq];
        if (numine==listnumeq[nineq]+1)
        {
            nummatrixcon++;
            indice=0;
            indicej=0;
            numine=1;
            nineq++;
            kk++;
        }

        if (indice==0)
            indice++;
        indcont=0;

        indcont=listsizerest[kk];
        mdebug()<<"indicej"<<indicej;
        if (indicej==indcont)
        {
                    indice++;
                    indicej=indice;

                    eqind=1;

         }


         if (eqind!=1)
           indicej++;




        COEFF value3;
        for(typename typP::value_type::iterator itera1=itera->begin();
            itera1!=itera->end();itera1++)
        {


            nummatrix=mapmonbase[*itera1];

            value3=(COEFF)itera1->GetCoeff();

            if (nummatrix==dimmonap) nummatrix=0;
            slv->coeff(nummatrix,nummatrixcon,indice-1,indicej-1)=Scalar::castDouble(value3);
            slv->nentries(nummatrix)++;



        }

    }





    mdebug()<<"avant obj";

    struct blockmatrix<COEFF> X,Z;
    COEFF *y,*pp;
    //COEFF pobj,dobj;
    int tmpsize=(int)stockksur2.size(),funcionnonnull=0;

    //on define pp comme le coeff independ de la function a minimizer et on l'initialize a 0
    pp=new COEFF[1];
    pp[0]=0;

    //on define bb comme les coeffs de la function a minimizer et on l'initialize a 0
//#ifndef NONEW
//    bb=new COEFF[nbconstraint+1];
//#else
//    bb=(COEFF *)malloc((nbconstraint+1)*sizeof(COEFF));
//#endif


    for(i=1;i<=dimmonap-1;i++)
        slv->obj_coeff(i)=0.0;
      //  bb[i]=0.0;





    int count_lis;
    COEFF value2;


    //on rempli bb comme les coeffs de la function
    typename typP::iterator iterfuncion=F.begin();
    for(typename typP::value_type::iterator itermon=iterfuncion->begin();
        itermon!=iterfuncion->end();itermon++)
    {

        value2=(COEFF)itermon->GetCoeff();

        if (*itermon!=itermon->GetCoeff())
        {
            count_lis=1;
            for(typename typP::iterator iterlis=listmin.begin();iterlis!=listmin.end();
                iterlis++,count_lis++)

                if (isEgal(*itermon,*(iterlis->begin()),b))
                    //bb[count_lis]=itermon->GetCoeff();
                    slv->obj_coeff(count_lis)=Scalar::castDouble(value2);

        }
        else
            pp[0]=Scalar::castDouble(value2);
    }


//    for(i=0;i<=nbconstraint;i++)
//        if (bb[i]!=0) funcionnonnull=1;


    mdebug()<<"avant run";
    slv->run();
    mdebug()<<"apres run";




    //std::cout<<"Extended moments:"<<moments<<std::endl;


    char nom[50],command[1024],copy[1024],copy1[1024];
    static int d=0,t=0;
    int r;

    Polynomial p=1;

    int j=0;
    listmonLambda1.sort(Rlex<Monomial>());
    for(typename list<Monomial>::iterator iter=listmonLambda1.begin();iter!=listmonLambda1.end();
            iter++)
    {


            mdebug()<<"mon"<<(*iter)<<(COEFF)slv->output[j];
            p+=(*iter)*slv->output[j];

        j++;

    }
    mdebug()<<"p"<<p;

        //on appele a la methode de reconstruction (decomposer), d abord on construit la function entree avec les donnes de sortie
        //du solver SDP
        //if ((method==0 && solver!=4 && (abs(sol2[0])*1e+2<1)) || (method==0 && solver==4 && r==0 && (solst==1 ||solst==8 || solst==0)))
        //if ((method==0 && solver!=4 ) || (method==0 && solver==4 && r==0 && (solst==1 ||solst==8 || solst==0)))
        if (method==0 )
        {
            // il faut faire transforme le vecteur lambda en p
            //il faut definir le variable par defaut


            double  cof;
            double *listcoef;

            typP  projmonomesadd,prodmonomestot, monomesadd;

            list<typename typserv::monom_t> monomestot;




            enumeratemon(monomestot,k/2,b,dump);
            //enumeratemon(monomestot,1,b,dump);
            computeprod(prodmonomestot,monomestot,b,dump);
            //mdebug()<<"falla aqui";
            int nuevomon;
            for(typename typP::iterator itert=prodmonomestot.begin(); itert!=prodmonomestot.end();itert++)
            {
                nuevomon=0;
                for(typename typP::iterator iter=listLambda.begin();iter!=listLambda.end();
                    iter++)
                {
                    //               mdebug()<<"itert->begin()"<<*itert->begin();
                    //               mdebug()<<"iter->begin()"<<*iter->begin();
                    if (isEgal(*(itert->begin()),*(iter->begin()),b))
                    {
                        nuevomon=1;
                    }
                }
                if (nuevomon==0 && *itert!=1)
                {
                   // mdebug()<<"itert"<<*itert;
                    monomesadd.push_back(*itert);
                    projmonomesadd.push_back(*itert);
                }
            }


            //mdebug()<<"monomesadd"<<monomesadd1.size();
            mdebug()<<"monomesadd"<<monomesadd.size();
            int k=0;
            COEFF coeff2;
            //if (!monomesadd1.empty())
            if (!monomesadd.empty())
            {
                proj(projmonomesadd,dump,b,serv);
                listcoef=new double[projmonomesadd.size()];
                //listcoef=new COEFF[projmonomesadd.size()];
                for(typename typP::iterator itadd=projmonomesadd.begin(); itadd!=projmonomesadd.end();
                    itadd++)
                {
                    //mdebug()<<"itadd"<<*itadd;
                    cof=0;
                    for(typename typP::value_type::iterator itaddmon=itadd->begin();
                        itaddmon!=itadd->end();itaddmon++)
                    {
                        for(Polynomial::iterator itp=p.begin();itp!=p.end();
                            itp++)
//                        for(Poly::iterator itp=p.begin();itp!=p.end();
//                                itp++)
                           {
//                            mdebug()<<"itaddmon"<<*itaddmon;
//                            mdebug()<<"itp"<<*itp;
                            if (isEgal(*itaddmon,*itp,b))
                            {
                               // mdebug()<<"coeff(itp)"<<itp->GetCoeff();
                                coeff2=(COEFF)itaddmon->GetCoeff();
                                cof+=Scl<MPF>::castDouble(coeff2)*(itp->GetCoeff());
                            }
                           }

                    }
                    mdebug()<<"coef"<<cof;
                    listcoef[k]=cof;
                    k++;
                }

                int j=0;
                int s=projmonomesadd.size();
                for(typename list<Monomial>::iterator itp1=monomesadd1.begin(); itp1!=monomesadd1.end();
                    itp1++)
//                for(typename typP::iterator itp1=monomesadd.begin(); itp1!=monomesadd.end();
//                                    itp1++)
                {
                    if (*itp1!=1)
                    {
                     if (listcoef[s-j-1]==0) listcoef[s-j-1]=10e-20;
                     mdebug()<<"ipt"<<*itp1;
                     p+=(*itp1)*listcoef[s-j-1];
                     j++;
                    }

                }

            }

            int rank_stop=factorial(grado/2+(b.nbvar()))/(factorial(grado/2)*factorial(b.nbvar()));
            mdebug()<<"rank_stop"<<rank_stop;

            //double threshold=0.00003;
            mdebug()<<tresh;
            mmx::Decomposer<Polynomial> tens(p,tresh,0,-1,-1,rank_stop+1,0);

            tens.main_algo();
            mdebug()<<tens.lambdasopt.size();


            for (int i=0;i<tens.lambdasopt.size();i++)
                if (tens.lambdasopt[i].real()>tresh) numsol++;


            mdebug()<<"nombre de solutions"<<numsol;

            if (numsol!=0)
            {
                solfinal= new coeff[numsol*tens.nb_var+2];
                solfinal[0]=numsol;
                solfinal[1]=slv->fminprimal+pp[0];
                mdebug()<<"fmin"<<slv->fminprimal;
                mdebug()<<"p0"<<pp[0];


                int k=2;
                for (int j=0;j< tens.lambdasopt.size();j++)
                    if (tens.lambdasopt[j].real()>tresh)

                        for(unsigned i=0;i<tens.linear_forms.nbrow();i++)
                        {
                            solfinal[k]=tens.linear_forms(i,j).real();
                            k++;
                        }
                mdebug()<<"on a les bonnes solutions";
            }
            else
            {
                solfinal= new coeff[2];
                solfinal[0] = 0;
                solfinal[1]=slv->fmindual+pp[0];

            }


        }
        else{

            solfinal= new coeff[2];
            solfinal[0] = 0;
            solfinal[1]=slv->fminprimal+pp[0];
        }


   // }
//    else
//    { //si dimension est egal a 1
//        //list<COEFF> talist;
//        sol1=new coeff[dim];
//        sol1[0]=(coeff)1;
//        sol2=new coeff[6];
//        for(i=0;i<6;i++)
//            sol2[i]=0;

//        method=1;
//        solfinal=new coeff[1];
//        solfinal[0]=-1;

//    }
    mdebug()<<"method"<<method;
    //
    //On calcul le noyau
    //

    if (method==1 && solst!=5)
    {
        mdebug()<<"error";
        long n=stockksur2.size();
        double *V=new double[n*n];
        double *S=new double[n];
        sol1=new double[n*n];
        double *work=new double[6*n];
        long eps1=1e5;
        long eps2=1e7;//AVANT CT 8
        long eps3=1e17;
        long eps4=1e14;
        long lwork=6*n,info;
        int kk=0;//servira pour le rang 0
        int rank=0;
        double *listrank;
        listrank=new double[dim];
        std::map<Monomial,double> coeffm;
        Monomial m;

        std::map<pair<int,int>,double>  mat;


        for(typename Polynomial::const_iterator it=p.begin();it!= p.end();it++) {
            m =*it;
            m.SetCoeff(1);
            coeffm[m]=it->GetCoeff();
        }
        Polynomial lm = matrix_of<Polynomial>::choose(numvar,k/2);

        mdebug()<<"lm:"<<lm<<"lm.size()"<<lm.size();
        int i=0,j=0;
        for(typename Polynomial::const_iterator it= lm.begin();it!= lm.end();it++)
        {
            j=i;
            for(typename Polynomial::const_iterator jt= it; jt!= lm.end();jt++) {
                m=(*it)*(*jt);
                //std::cout<<"mon"<<m<<std::endl;


                    mat[make_pair(i,j)]=coeffm[m];
                    mat[make_pair(j,i)]=coeffm[m];


                    mdebug()<<",i:"<<i<<",j:"<<j<<",value"<<coeffm[m];
                j++;
            }
            i++;
        }
        int z=0;
        for (int s=0;s<n;s++)
            for (int t=0;t<n;t++)
            {
                sol1[z]=mat[make_pair(s,t)];
                z++;
            }

        r=dim-2;
//        int s=slv->output.size();
//        sol1[0]=1;
//        for (int k=0;k<n;k++)
//          //for (int j=0;j<n;j++)
//            for(int i=0;i<n;i++)
//        {
//            if (k*n+i==0) sol1[k*n+i]=0;
//            else
//            {
//                if (k>=1) sol1[k*n+i]= slv->output[s-i]
//             sol1[k*n+i]=slv->output[s-i];

//            }
//        }
        //   mdebug()<<sol2[0];
        //  if (abs(sol2[0])*1e+2<1  )
        // {
        list<double> sol1bis;

        for (int i=0;i<n*n;i++)
        {
            sol1bis.push_back(sol1[i]);
            // mdebug()<<"sol1"<<sol1[i];
        }


        dgesvd_((char*)"N",(char*)"A",&n,&n,sol1,
                &n,S,(double*)NULL,&n,V,&n,work,&lwork,&info);

        for (long m=dim; m!=1;m--)
        {
            double *sol1b;

            typename list<double>::iterator it1,it2;




            it1=sol1bis.begin();

            for (int j=1;j<=m;j++)
            {

                if (j!=m)
                {
                    advance(it1,m-1);

                    it1=sol1bis.erase(it1);
                }
                else
                {
                    it2=it1;
                    advance(it2,m);
                    sol1bis.erase(it1,it2);
                }

            }

            sol1b=new double[sol1bis.size()];
            it1=sol1bis.begin();
            for(int k=0;k<sol1bis.size();it1++,k++)
            {

                sol1b[k]=*it1;
            }

            n=m-1;
            double *V1=new double[n*n];
            double *S1=new double[n];
            double *work1=new double[6*n];
            long lwork1=6*n,info1;
            int k1=0;


            dgesvd_((char*)"N",(char*)"A",&n,&n,sol1b,
                    &n,S1,(double*)NULL,&n,V1,&n,work1,&lwork1,&info1);

            for(k1=1;k1<n;k1++)

            {
                if (k1==1)
                {

                    if (solver==1 || solver==3 ||solver==4)
                        if (abs(S1[0])>abs(S1[k1])*eps1) break;

                    if (solver==2)
                        if (abs(S1[0])>abs(S1[k1])*eps3) break;

                }
                else
                {

                    if (solver==1 || solver==3 ||solver==4)
                        if(abs(S1[0])>abs(S1[k1])*eps1 || abs(S1[k1-1])>abs(S1[k1])*eps1) break;

                    if (solver==2)
                    {

                        // if(abs(S[0])>abs(S[kk])*eps3 || abs(S[kk-1])>abs(S[kk])*eps2)break;
                        if(abs(S1[0])>abs(S1[k1])*eps4 )break;
                    }

                }

            }//MAGIC

            rank=k1;

            listrank[r]=k1;

            r--;
            mdebug()<<"rank"<<rank;
            delete[] sol1b;
            delete[] S1;
            delete[] V1;

        }
        // }

        if(method==1)
        {
            mdebug()<<"info svd "<<info;
            mdebug()<<"val sing X";
            for(int kk=0;kk<dim;kk++)
                mdebug()<<S[kk]<<" ";
            mdebug();





            // mdebug()<<"coiucou";
            for(kk=1;kk<dim;kk++)

            {   /*mdebug()<< abs(S[kk-1])<< endl;
    mdebug()<<abs(S[kk])*1e6;*/
                if (kk==1)
                {
                    //#if defined(SDPA) || defined(CSDP)
                    // mdebug()<<"svd0SDPA";
                    if (solver==1 || solver==3 ||solver==4)
                        if (abs(S[0])>abs(S[kk])*eps1) break;
                    //#endif
                    //#ifdef GMPSDPA
                    // mdebug()<<"svd0";
                    if (solver==2)
                        if (abs(S[0])>abs(S[kk])*eps3) break;
                    //#endif
                }
                else
                {
                    // #if defined(SDPA) || defined(CSDP)
                    if (solver==1 || solver==3 || solver==4)
                        if(abs(S[0])>abs(S[kk])*eps1 || abs(S[kk-1])>abs(S[kk])*eps1) break;
                    // #endif
                    // #ifdef GMPSDPA
                    if (solver==2)
                    {
                        mdebug()<<"svdvoisin"<<abs(S[kk])*eps3;
                        // if(abs(S[0])>abs(S[kk])*eps3 || abs(S[kk-1])>abs(S[kk])*eps2)break;
                        if(abs(S[0])>abs(S[kk])*eps4 )break;
                    }
                    // #endif
                }

            }//MAGIC
            //  if((abs(S[kk-1])>abs(S[kk])*eps1)) break;}//MAGIC
            mdebug()<<"kk ici "<<kk<<" et n "<<dim;
            //     mdebug()<<"V";

            //     for(int j=0;j<(int)dim;j++)
            //     {
            //       for(int jj=0;jj<dim;jj++)
            //        mdebug()<<V[dim*(jj)+j]<<" ";
            //       mdebug();
            //     }

            //sol2[4]=dim;
            //sol2[5]=kk;
            rank=kk;
            listrank[dim-1]=kk;

            mdebug()<<"grado"<<grado;
            mdebug()<<"rank final"<<rank;
            for(int l=0;l<dim;l++)
                mdebug()<<listrank[l];

            //flat extension condition gloptipoly
            int flatext=0;
            int j=0,j1,j2,c1,c2;
            int rankfinal;

            if (maxdegrest==0) j=1;
            else
            {
                if (maxdegrest%2==1) maxdegrest++;
                j=maxdegrest/2;
            }
            j1=Degree(*F.begin());
            if (j1%2==1)j1++;
            j1/=2;
            j2=max(j1,j);
            mdebug()<<"j"<<j;
            mdebug()<<"j2"<<j2;
            for (;j2<=grado/2;j2++)
            {

                c1=(factorial(numvar+j2)/(factorial(j2)*factorial(numvar)))-1;
                c2=(factorial(numvar+j2-j)/(factorial(j2-j)*factorial(numvar)))-1;
                mdebug()<<"c1 "<<c1<<"c2 "<<c2;

                if (listrank[c1]==listrank[c2])
                {
                    flatext=1;
                    rank=listrank[c1];
                     mdebug()<<"rank"<<rank;
                    break;
                }
            }
            mdebug()<<"flat extension"<<flatext;
            flat=flatext;


            //
            //On construit les nouvelles equations
            //

            int count_eq;
            int keri,kerf;
            kerf=c1;
            keri=c2+1;
            rankfinal=rank;

            for(;rank<dim;rank++)
            {
#if 1
                // mdebug()<<"equation dans le noyau";
#endif
                count_eq=1;
                typedef typename typP::value_type polyalp;
                polyalp tmp(0);
                for(typename list<typename typserv::monom_t>::iterator
                    iter=stockksur2.begin();
                    iter!=stockksur2.end();iter++,count_eq++)
                {


                    if(!Iszero(abs(V[rank+dim*(count_eq-1)])))

                    {

                        tmp+=polyalp(*iter*(COEFF)V[rank+dim*(count_eq-1)]);

                    }

                }

                mdebug()<<"g trouve ca comme nvlle eq "<<tmp;
                newpol.push_back(tmp);

            }

            solfinal[0]=rankfinal;
        }


        delete[] S;
        delete[] V;
        //delete[] work;

    }
    //delete[] sol1;
    // #ifdef CSDP
    //if (method==1 ) delete[] sol1;
//    if (solve3)
//        if (dim!=1)
//            delete[] sol3;
    //#endif





    delete [] constraints2;
    delete[] blockptr2;
    delete [] bb;

}

template<typename mon, typename mon2, typename Base>
int isEgal(const mon & m1,const mon2 & m2,const Base & b)
{
    const int nbvar=b.nbvar();
    int res=1;
    if((m1.GetCoeff())==0)
        res=0;
    else
        for(int i=0;i<nbvar;i++)
        {
            //mdebug()<<"    isdivisble  "<<m2.GetDegree(i)<<" "<<m1.GetDegree(i);
            if(m2.GetDegree(i)!=m1.GetDegree(i))
            { res=0;break;};
        }
    return res;
}

template<typename T>
int dgesvd_(char *jobu, char *jobvt, long *m, long *n,
            T *a, long *lda, T *s, T *u, long *
            ldu, T *vt, long *ldvt, T *work, long *lwork,
            long *info);