Changeset d5042e in git


Ignore:
Timestamp:
Jun 22, 2009, 5:17:37 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
Children:
587cc52fbea6eb31ce24d13c81748bbe2743326e
Parents:
08c2d636b6e8ba70134c33364e4e4a463aac2d13
Message:
Modified facet::facet() to set everything to 0 in the beginning
New method enqueueNewFacets
BUG in the way getConeNormals creates the list of normals. There is one facet to much created everytime!


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r08c2d6 rd5042e  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-06-19 11:22:06 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.63 2009-06-19 11:22:06 monerjan Exp $
    6 $Id: gfan.cc,v 1.63 2009-06-19 11:22:06 monerjan Exp $
     4$Date: 2009-06-22 15:17:37 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.64 2009-06-22 15:17:37 monerjan Exp $
     6$Id: gfan.cc,v 1.64 2009-06-22 15:17:37 monerjan Exp $
    77*/
    88
    … …  
    105105                        // Pointer to next facet.  */
    106106                        /* Defaults to NULL. This way there is no need to check explicitly */
    107                         this->next=NULL;
    108                         this->UCN=NULL;
     107                        this->fNormal=NULL;
     108                        this->interiorPoint=NULL;               
     109                        this->UCN=0;
    109110                        this->codim2Ptr=NULL;
    110111                        this->codim=1;          //default for (codim-1)-facets
    111                         this->numCodim2Facets=0;
     112                        this->flipGB=NULL;
     113                        this->isIncoming=FALSE;
     114                        this->next=NULL;
     115                        this->codim2Ptr=NULL;
     116                        this->numCodim2Facets=0;                                               
     117                        this->flipRing=NULL;
    112118                }
    113119               
    … …  
    117123                {
    118124                        this->next=NULL;
    119                         this->UCN=NULL;
     125                        this->UCN=0;
    120126                        this->codim2Ptr=NULL;
    121127                        if(n>1)
    … …  
    489495                        // End of var declaration
    490496#ifdef gfan_DEBUG
    491                         cout << "The Groebner basis has " << lengthGB << " elements" << endl;
    492                         cout << "The current ring has " << numvar << " variables" << endl;
     497//                      cout << "The Groebner basis has " << lengthGB << " elements" << endl;
     498//                      cout << "The current ring has " << numvar << " variables" << endl;
    493499#endif
    494500                        cols = numvar;
    … …  
    502508                        }
    503509#ifdef gfan_DEBUG
    504                         cout << "rows=" << rows << endl;
    505                         cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
     510//                      cout << "rows=" << rows << endl;
     511//                      cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
    506512#endif
    507513                        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
    … …  
    518524                                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
    519525#ifdef gfan_DEBUG
    520                                 cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
     526//                              cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
    521527#endif
    522528               
    … …  
    585591                        /*The pointer *fRoot should be the return value of this function*/
    586592                        facet *fRoot = new facet();     //instantiate new facet
     593                        fRoot->setUCN(this->getUCN());
    587594                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
    588595                        facet *fAct;                    //instantiate pointer to active facet
    … …  
    597604                                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
    598605#ifdef gfan_DEBUG
    599                                         std::cout << "fAct is " << foo << " at " << fAct << std::endl;
     606//                                      std::cout << "fAct is " << foo << " at " << fAct << std::endl;
    600607#endif
    601608                                        (*load)[jj-1] = (int)foo;               //store typecasted entry at pos jj-1 of load           
    … …  
    608615                                        intvec *ivCanonical = new intvec(load->length());
    609616                                        (*ivCanonical)[jj]=1;
    610                                         cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl;
     617//                                      cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl;
    611618                                        if (dotProduct(*load,*ivCanonical)<0)   
    612619                                        //if (ivMult(load,ivCanonical)<0)
    … …  
    624631                                else
    625632                                {       /*Now load should be full and we can call setFacetNormal*/
     633                                        //NOTE Check whether another facet is needed at all. Otherwise
     634                                        //this results in one facet too much being constructed!
    626635                                        fAct->setFacetNormal(load);
    627636                                        fAct->next = new facet();                                       
    628637                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
     638                                        fAct->setUCN(this->getUCN());
    629639                                        this->numFacets++;
    630640                                }//if (isFlippable==FALSE)
    … …  
    693703                                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    694704                                P=dd_CopyGenerators(ddpolyh);
    695                                 dd_WriteMatrix(stdout,P);
     705//                              dd_WriteMatrix(stdout,P);
    696706                                       
    697707                                /* We loop through each row of P
    … …  
    831841                        ina=idrCopyR(initialForm,srcRing);                     
    832842#ifdef gfan_DEBUG
    833                         cout << "ina=";
    834                         idShow(ina); cout << endl;
     843//                      cout << "ina=";
     844//                      idShow(ina); cout << endl;
    835845#endif
    836846                        ideal H;
    … …  
    15801590                        gcone *gcAct;
    15811591                        gcAct = &gcRoot;
    1582                         gcone *gcPtr;
     1592                        gcone *gcPtr;   //Pointer to end of linked list of cones
    15831593                        gcPtr = &gcRoot;
     1594                        gcone *gcNext;  //Pointer to next cone to be visited
     1595                        gcNext = &gcRoot;
     1596                        gcone *gcHead;
     1597                        gcHead = &gcRoot;
    15841598                       
    15851599                        facet *fAct;
    … …  
    16111625                        SearchListAct = SearchListRoot; //Set to beginning of list
    16121626                       
    1613                         fAct = gcAct->facetPtr;
    1614                         if(areEqual(fAct->getFacetNormal(),fAct->next->getFacetNormal()))
    1615                         {
    1616                                 cout <<"HI" << endl;
    1617                         }
    1618                        
    1619                         gcAct->writeConeToFile(*gcAct);
     1627                        fAct = gcAct->facetPtr;                 
     1628                        //gcAct->writeConeToFile(*gcAct);
    16201629                       
    16211630                        /*End of initialisation*/
    … …  
    16261635                        */
    16271636                        while(SearchListAct->next!=NULL)
    1628                         {//NOTE See to it that the cone is only changed after ALL facets have been flipped!
    1629                                 //As of now this is not the case!
    1630                                 int flag=1;
     1637                        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
    16311638                                fAct = SearchListAct;
    16321639                                //while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) )
    … …  
    16341641                                {
    16351642                                        gcAct->flip(gcAct->gcBasis,fAct);
    1636                                         ring rTmp=rCopy(SearchListAct->flipRing);
     1643                                        ring rTmp=rCopy(fAct->flipRing);
    16371644                                        rComplete(rTmp);
    16381645                                        rChangeCurrRing(rTmp);
    1639                                         gcone *gcTmp = new gcone::gcone(*gcAct,*SearchListAct);
     1646                                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);
    16401647                                        gcTmp->getConeNormals(gcTmp->gcBasis, FALSE);
    16411648                                        gcTmp->getCodim2Normals(*gcTmp);
    1642                                         /*add facets to SLA here*/                                     
     1649                                        gcTmp->normalize();
     1650                                        /*add facets to SLA here*/
     1651                                        gcTmp->enqueueNewFacets(*SearchListRoot);
    16431652                                        rChangeCurrRing(gcAct->baseRing);
    16441653                                        gcPtr->next=gcTmp;
    1645                                         if(flag==1)
    1646                                                 gcAct->next=gcTmp;
    16471654                                        gcPtr=gcPtr->next;
    16481655                                        if(fAct->getUCN() == fAct->next->getUCN())
    … …  
    16521659                                        else
    16531660                                                break;
    1654                                         flag++;
    1655                                 }while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );   
    1656                                 if (fAct->next!=NULL)
    1657                                 {
    1658                                         SearchListAct=fAct->next;
    1659                                         ring r=rCopy(SearchListAct->flipRing);
    1660                                         rComplete(r);
    1661                                         rChangeCurrRing(r);
    1662                                         gcAct = gcAct->next;
     1661                                }while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
     1662                                gcNext = gcHead;
     1663                                while(gcNext->next!=NULL)
     1664                                {
     1665                                        if(gcNext->getUCN() > gcNext->next->getUCN() )
     1666                                        {
     1667                                                gcAct = gcNext->next;
     1668                                        }
     1669                                        gcNext = gcNext->next;
    16631670                                }
    1664                                 else
    1665                                 {
    1666                                         SearchListAct=NULL;
    1667                                 }
     1671                                SearchListAct = SearchListAct->next;
     1672//                              if (fAct->next!=NULL)
     1673//                              {
     1674//                                      SearchListAct=fAct->next;
     1675//                                      ring r=rCopy(SearchListAct->flipRing);
     1676//                                      rComplete(r);
     1677//                                      rChangeCurrRing(r);
     1678//                                      gcAct = gcAct->next;
     1679//                              }
     1680//                              else
     1681//                              {
     1682//                                      SearchListAct=NULL;
     1683//                              }
    16681684                               
    1669                         }                       
     1685                        }
     1686                                               
     1687//                      while(SearchListAct->next!=NULL)
     1688//                      {
     1689//                              fAct = SearchListAct;
     1690//                              do
     1691//                              {
     1692//                                      gcAct->flip(gcAct->gcBasis,fAct);
     1693//                                      ring rTmp = rCopy(fAct->flipRing);
     1694//                                      rComplete(rTmp);
     1695//                                      rChangeCurrRing(rTmp);
     1696//                                      gcone *gcTmp = new gcone(*gcAct,*fAct);
     1697//                                      gcTmp->getConeNormals(gcTmp->gcBasis,FALSE);
     1698//                                      gcTmp->getCodim2Normals(*gcTmp);
     1699//                                      gcPtr->next = gcTmp;
     1700//                                      gcPtr = gcPtr->next;
     1701//                                      /*add facets to SLA here*/
     1702//                                      rChangeCurrRing(gcAct->baseRing);
     1703//                                      fAct = fAct->next;
     1704//                              }while( (fAct->next!=NULL) &&  (fAct->getUCN()==fAct->next->getUCN()));
     1705//                              SearchListAct = SearchListAct->next;//fAct;                             
     1706//                      }
     1707                       
     1708//                      while(SearchListAct->next!=NULL)
     1709//                      {
     1710//                              gcAct->flip(gcAct->gcBasis,SearchListAct);
     1711//                              ring rTmp=rCopy(fAct->flipRing);
     1712//                              rComplete(rTmp);
     1713//                              rChangeCurrRing(rTmp);
     1714//                              gcone *gcTmp = new gcone(*gcAct,*SearchListAct);
     1715//                              gcTmp->getConeNormals(gcTmp->gcBasis,FALSE);
     1716//                              gcTmp->getCodim2Normals(*gcTmp);
     1717//                              rChangeCurrRing(gcAct->baseRing);
     1718//                              SearchListAct = SearchListAct->next;
     1719//                             
     1720//                              //gcone *gcNew = new gcone(SearchListAct->flipRing,SearchListAct->flipGB);
     1721//                              //gcAct = gcNew;
     1722//                      }
    16701723               
    16711724                        //NOTE Hm, comment in and get a crash for free...
    … …  
    17761829                        }                       
    17771830                }
     1831                /**
     1832                * Takes ptr to search list root
     1833                */
     1834                void enqueueNewFacets(facet &f)
     1835                {
     1836                        facet *slAct;   //called with f=SearchListRoot
     1837                        slAct = &f;
     1838                        facet *slEnd;   //Pointer to end of SLA
     1839                        slEnd = &f;             
     1840                        facet *fAct;
     1841                        fAct = this->facetPtr;
     1842                        facet *codim2Act;
     1843                        codim2Act = this->facetPtr->codim2Ptr;
     1844                        while(slEnd->next!=NULL)
     1845                        {
     1846                                slEnd=slEnd->next;
     1847                        }
     1848                        /*1st step: compare facetNormals*/
     1849                        intvec *fNormal = new intvec(this->numVars);
     1850                        intvec *f2Normal = new intvec(this->numVars);
     1851                        intvec *slNormal = new intvec(this->numVars);
     1852                        intvec *sl2Normal = new intvec(this->numVars);
     1853                        while(slAct->next!=NULL)
     1854                        {
     1855                                slNormal = slAct->getFacetNormal();
     1856                                fAct = this->facetPtr;
     1857                                while(fAct->next!=NULL)
     1858                                {
     1859                                        fNormal = fAct->getFacetNormal();
     1860                                        if(!isParallel(fNormal,slNormal))
     1861                                        {                                               
     1862                                                slEnd->next = new facet();
     1863                                                slEnd = slEnd->next;
     1864                                                slEnd->setUCN(this->getUCN());
     1865                                                slEnd->setFacetNormal(fNormal);                                         
     1866                                                //memcpy(slEnd,fAct,sizeof(facet));
     1867                                        }
     1868                                        fAct = fAct->next;
     1869                                }
     1870                                slAct = slAct->next;
     1871                        }
     1872//                      delete sl2Normal;
     1873//                      delete slNormal;
     1874//                      delete f2Normal;
     1875//                      delete fNormal;
     1876                }//addC2N
    17781877               
    17791878                /** \brief Compute the gcd of two ints
Note: See TracChangeset for help on using the changeset viewer.