Changeset b1bed35 in git for kernel/gfan.cc


Ignore:
Timestamp:
Apr 17, 2009, 3:44:27 PM (15 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
694257b820eac75a8100417ebf8f059879a873fd
Parents:
ee5ace2b2bf7d30c23b3fce50222dcc3f12e88dc
Message:
Markings work fine, though the code is pretty ugly. Needs some clean-up.
Newly introduced method void gcone::interiorPoint(dd_MatrixPtr &M) works well. Drawback: As of now no rational or integer is returned. Intended workaround: Continued fractions


git-svn-id: file:///usr/local/Singular/svn/trunk@11730 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    ree5ace rb1bed35  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-04-16 09:59:16 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.34 2009-04-16 09:59:16 monerjan Exp $
    6 $Id: gfan.cc,v 1.34 2009-04-16 09:59:16 monerjan Exp $
     4$Date: 2009-04-17 13:44:27 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.35 2009-04-17 13:44:27 monerjan Exp $
     6$Id: gfan.cc,v 1.35 2009-04-17 13:44:27 monerjan Exp $
    77*/
    88
    … …  
    463463                        srcRing_H=idrCopyR(H,dstRing);
    464464                        idShow(srcRing_H);                     
    465                         srcRing_HH=ffG(srcRing_H,this->gcBasis);                       
     465                        srcRing_HH=ffG(srcRing_H,this->gcBasis);               
    466466                        idShow(srcRing_HH);
    467467                        /*Substep 2.2.1
    … …  
    483483                                iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
    484484                        }
    485                         intPointMatrix = dd_CreateMatrix(iPMatrixRows,this->numVars+1);
     485                        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
     486                        construction of the differences
     487                        */
     488                        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
    486489                       
    487490                        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
    488491                        {
     492                                markingsAreCorrect=FALSE;       //crucial to initialise here
    489493                                poly aktpoly=srcRing_HH->m[ii];
    490494                                /*Comparison of leading monomials is done via exponent vectors*/
    … …  
    494498                                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
    495499                                        pGetExpV(aktpoly,src_ExpV);
     500                                        rChangeCurrRing(dstRing);       //this ring change is crucial!
    496501                                        pGetExpV(pCopy(H->m[ii]),dst_ExpV);
    497                                         cout << *src_ExpV << endl;
    498                                         cout << *dst_ExpV << endl;
    499                                         if (src_ExpV == dst_ExpV)
     502                                        rChangeCurrRing(srcRing);
     503                                        bool expVAreEqual=TRUE;
     504                                        for (int kk=1;kk<=this->numVars;kk++)
     505                                        {
     506                                                cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
     507                                                if (src_ExpV[kk]!=dst_ExpV[kk])
     508                                                {
     509                                                        expVAreEqual=FALSE;
     510                                                }
     511                                        }                                       
     512                                        //if (*src_ExpV == *dst_ExpV)
     513                                        if (expVAreEqual==TRUE)
    500514                                        {
    501515                                                markingsAreCorrect=TRUE; //everything is fine
    … …  
    514528                                else
    515529                                {
    516                                         pGetExpV(pCopy(pHead(H->m[ii])),leadExpV); //We use H->m[ii] as leading monomial
     530                                        rChangeCurrRing(dstRing);
     531                                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
     532                                        rChangeCurrRing(srcRing);
    517533                                }
    518                                 /*compute differences of the expvects*/
     534                                /*compute differences of the expvects*/                         
    519535                                while (pNext(aktpoly)!=NULL)
    520536                                {
    521                                         aktpoly=pNext(aktpoly);
    522                                         pGetExpV(aktpoly,v);                                           
     537                                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
     538                                        is not omitted when computing the differences*/
     539                                        if(markingsAreCorrect==TRUE)
     540                                        {
     541                                                aktpoly=pNext(aktpoly);
     542                                                pGetExpV(aktpoly,v);
     543                                        }
     544                                        else
     545                                        {
     546                                                pGetExpV(pHead(aktpoly),v);
     547                                                markingsAreCorrect=TRUE;
     548                                        }
     549
    523550                                        for (int jj=0;jj<this->numVars;jj++)
    524551                                        {                               
    525                                                 /*Store into ddMatrix*/                 
    526                                                 /*FIXME Wrong values*/
    527                                                 dd_set_si(intPointMatrix->matrix[aktrow][jj+1],v[jj+1]-leadExpV[jj+1]);
     552                                                /*Store into ddMatrix*/                                                         
     553                                                dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
    528554                                        }
    529555                                        aktrow +=1;
    … …  
    532558                                delete leadExpV;
    533559                        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
     560                        /*Now we add the constraint for the standard simplex*/
     561                        /*NOTE:Might actually work without the standard simplex*/
     562                        //dd_set_si(intPointMatrix->matrix[aktrow][0],100);
     563                        for (int jj=1;jj<=this->numVars;jj++)
     564                        {
     565                                //dd_set_si(intPointMatrix->matrix[aktrow][jj],-1);
     566                        }
    534567                        dd_WriteMatrix(stdout,intPointMatrix);
     568                        interiorPoint(intPointMatrix);  //TODO This should finally return an intvec
    535569                        dd_FreeMatrix(intPointMatrix);
    536570                       
    … …  
    548582                * compute the factors \f$ a_i \f$
    549583                */
     584                //NOTE: Should be replaced by kNF or kNF2
    550585                poly restOfDiv(poly const &f, ideal const &I)
    551586                {
    … …  
    614649                /** \brief Compute \f$ f-f^{\mathcal{G}} \f$
    615650                */
     651                //NOTE: use kNF or kNF2 instead of restOfDivision
    616652                ideal ffG(ideal const &H, ideal const &G)
    617653                {
    … …  
    623659                        {
    624660                                res->m[ii]=restOfDiv(H->m[ii],G);
     661                                //res->m[ii]=kNF(H->m[ii],G);
    625662                                temp1=H->m[ii];
    626663                                temp2=res->m[ii];                               
    … …  
    631668                                //pSetm(res->m[ii]);
    632669                                cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                         
    633                         }
     670                        }                       
    634671                        return res;
    635672                }
    … …  
    686723                        }
    687724                }//bool isParallel
     725               
     726                void interiorPoint(dd_MatrixPtr &M) //no const &M here since we want to remove redundant rows
     727                {
     728                        dd_LPPtr lp,lpInt;
     729                        dd_ErrorType err=dd_NoError;
     730                        dd_LPSolverType solver=dd_DualSimplex;
     731                        dd_LPSolutionPtr lpSol;
     732                        dd_rowset ddlinset,ddredrows;
     733                        dd_rowindex ddnewpos;
     734                        dd_NumberType numb;
     735                        numb=dd_Integer;
     736                       
     737                        M->numbtype=numb;
     738                        dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
     739                        dd_WriteMatrix(stdout,M);
     740                       
     741                        lp=dd_Matrix2LP(M, &err);
     742                        lpInt=dd_MakeLPforInteriorFinding(lp);
     743                        dd_LPSolve(lpInt,solver,&err);
     744                        lpSol=dd_CopyLPSolution(lpInt);
     745                        lpSol->numbtype=numb;
     746                        cout << "Interior point: ";
     747                        for (int ii=1; ii<(lpSol->d)-1;ii++)
     748                        {
     749                                dd_WriteNumber(stdout,lpSol->sol[ii]);
     750                        }
     751                        dd_FreeLPData(lp);
     752                        dd_FreeLPSolution(lpSol);
     753                        dd_FreeLPData(lpInt);                   
     754                }//void interiorPoint(dd_MatrixPtr const &M)
    688755};//class gcone
    689756
Note: See TracChangeset for help on using the changeset viewer.