Changeset b7510d in git for kernel/gfan.cc


Ignore:
Timestamp:
Jun 23, 2009, 6:21:07 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
3449c9b5bf93acbaa0ef759a0a9f94f85ef0c032
Parents:
ae51772364a45a7d84beed113717c0c6f320de8c
Message:
BUGFIX: Creation of linked list of facets now working correctly
New method showSLA() - debugging purposes only
Fixed enumeration of cones via UCN -> works now
Enqueueing into SLA done
TODO: Find out why it crashes in flip(). Maybe the usual wrong ring issue


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    rae5177 rb7510d  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    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 $
     4$Date: 2009-06-23 16:21:07 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.65 2009-06-23 16:21:07 monerjan Exp $
     6$Id: gfan.cc,v 1.65 2009-06-23 16:21:07 monerjan Exp $
    77*/
    88
    … …  
    110110                        this->codim2Ptr=NULL;
    111111                        this->codim=1;          //default for (codim-1)-facets
     112                        this->numCodim2Facets=0;
    112113                        this->flipGB=NULL;
    113114                        this->isIncoming=FALSE;
    114                         this->next=NULL;
    115                         this->codim2Ptr=NULL;
    116                         this->numCodim2Facets=0;                                               
     115                        this->next=NULL;               
    117116                        this->flipRing=NULL;
    118117                }
    … …  
    122121                facet(int const &n)
    123122                {
    124                         this->next=NULL;
     123                        this->fNormal=NULL;
     124                        this->interiorPoint=NULL;                       
    125125                        this->UCN=0;
    126126                        this->codim2Ptr=NULL;
    … …  
    128128                        {
    129129                                this->codim=n;
    130                         }//NOTE Handle exception here!
     130                        }//NOTE Handle exception here!                 
    131131                        this->numCodim2Facets=0;
     132                        this->flipGB=NULL;
     133                        this->isIncoming=FALSE;
     134                        this->next=NULL;
     135                        this->flipRing=NULL;
    132136                }
    133137               
    … …  
    222226                int getUCN()
    223227                {
    224                         return this->UCN;
     228                        if(this!=NULL)
     229                                return this->UCN;
     230                        else
     231                                return -1;
    225232                }
    226233               
    … …  
    233240                {
    234241                        return this->interiorPoint;
    235                 }              
     242                }       
    236243               
    237244                /*bool isFlippable(intvec &load)
    … …  
    294301                /* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
    295302                intvec *ivIntPt;        //an interior point of the cone
    296                 static int UCN;         //unique number of the cone
     303                int UCN;                //unique number of the cone
     304                static int counter;
    297305               
    298306        public:
    … …  
    321329                        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
    322330                        this->baseRing=currRing;
    323                         this->UCN++;
     331                        this->counter++;
     332                        this->UCN=this->counter;                       
    324333                        this->numFacets=0;
    325334                }
    … …  
    339348                        this->inputIdeal=I;
    340349                        this->baseRing=currRing;
    341                         this->UCN++;
     350                        this->counter++;
     351                        this->UCN=this->counter;                       
    342352                        this->numFacets=0;
    343353                }
    … …  
    381391                {
    382392                        this->next=NULL;
    383                         this->numVars=gc.numVars;
    384                         //this->UCN=(gc.UCN)+1;
    385                         this->UCN++;
     393                        this->numVars=gc.numVars;                                               
     394                        this->counter++;
     395                        this->UCN=this->counter;                       
    386396                        this->facetPtr=NULL;
    387397                        this->gcBasis=idCopy(f.flipGB);
    … …  
    435445                        intvec *iv = new intvec(this->numVars);
    436446                       
    437                         while (f->next!=NULL)
     447                        //while (f->next!=NULL)
     448                        while(f!=NULL)
    438449                        {
    439450                                iv = f->getFacetNormal();
    … …  
    442453                        }
    443454                        //delete iv;                   
     455                }
     456               
     457                /** For debugging purposes only */
     458                void showSLA(facet &f)
     459                {
     460                        facet *fAct;
     461                        fAct = &f;
     462                        intvec *n = new intvec(this->numVars);
     463                        cout << endl;
     464                        while(fAct!=NULL)
     465                        {
     466                                n=fAct->getFacetNormal();
     467                                n->show();
     468                                fAct = fAct->next;
     469                                cout << endl;
     470                                cout << "---" << endl;
     471                               
     472                        }
    444473                }
    445474               
    … …  
    590619#endif
    591620                        /*The pointer *fRoot should be the return value of this function*/
    592                         facet *fRoot = new facet();     //instantiate new facet
    593                         fRoot->setUCN(this->getUCN());
    594                         this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
    595                         facet *fAct;                    //instantiate pointer to active facet
    596                         fAct = fRoot;                   //Seems to do the trick. fRoot and fAct have to point to the same adress!
     621                        //facet *fRoot = new facet();   //instantiate new facet
     622                        //fRoot->setUCN(this->getUCN());
     623                        //this->facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
     624                        facet *fAct;                    //pointer to active facet
     625                        //fAct = fRoot;                 //Seems to do the trick. fRoot and fAct have to point to the same adress!
    597626                        //std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
    598627                        for (int kk = 0; kk<ddrows; kk++)
    … …  
    625654                                if (isFlippable==FALSE)
    626655                                {
     656#ifdef gfan_DEBUG
    627657                                        cout << "Ignoring facet";
    628658                                        load->show();
    629659                                        //fAct->next=NULL;
     660#endif
    630661                                }
    631662                                else
    632                                 {       /*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!
     663                                {       /*Now load should be full and we can call setFacetNormal*/                                     
     664                                        this->numFacets++;
     665                                        if(this->numFacets==1)
     666                                        {
     667                                                facet *fRoot = new facet();
     668                                                this->facetPtr = fRoot;
     669                                                fAct = fRoot;
     670                                        }
     671                                        else
     672                                        {
     673                                                fAct->next = new facet();
     674                                                fAct = fAct->next;
     675                                        }
    635676                                        fAct->setFacetNormal(load);
    636                                         fAct->next = new facet();                                       
    637                                         fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
    638                                         fAct->setUCN(this->getUCN());
    639                                         this->numFacets++;
     677                                        fAct->setUCN(this->getUCN());                           
     678                                        //fAct->setFacetNormal(load);
     679                                        //fAct->next = new facet();                                     
     680                                        //fAct=fAct->next;      //this should definitely not be called in the above while-loop :D
     681                                        //fAct->setUCN(this->getUCN());
     682                                        //this->numFacets++;
    640683                                }//if (isFlippable==FALSE)
    641684                                //delete load;
    … …  
    678721                void getCodim2Normals(gcone const &gc)
    679722                {
    680                         this->facetPtr->codim2Ptr = new facet(2);       //instantiate a (codim-2)-facet                 
     723                        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet                 
    681724                        facet *codim2Act;
    682                         codim2Act = this->facetPtr->codim2Ptr;
     725                        //codim2Act = this->facetPtr->codim2Ptr;
    683726                       
    684727                        dd_set_global_constants();
    … …  
    700743                                dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                   
    701744                                       
    702                                         //dd_WriteMatrix(stdout,ddakt);
     745#ifdef gfan_DEBUG
     746//                              dd_WriteMatrix(stdout,ddakt);
     747#endif
    703748                                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    704749                                P=dd_CopyGenerators(ddpolyh);
    705 //                              dd_WriteMatrix(stdout,P);
     750#ifdef gfan_DEBUG
     751//                              dd_WriteMatrix(stdout,P);
     752#endif
    706753                                       
    707754                                /* We loop through each row of P
    … …  
    711758                                for (int jj=1;jj<=P->rowsize;jj++)
    712759                                {
     760                                        this->facetPtr->numCodim2Facets++;
     761                                        if(this->facetPtr->numCodim2Facets==1)
     762                                        {
     763                                                this->facetPtr->codim2Ptr = new facet(2);
     764                                                codim2Act = this->facetPtr->codim2Ptr;
     765                                        }
     766                                        else
     767                                        {
     768                                                codim2Act->next = new facet(2);
     769                                                codim2Act = codim2Act->next;
     770                                        }
    713771                                        intvec *n = new intvec(this->numVars);
     772                                        makeInt(P,jj,*n);
     773                                        codim2Act->setFacetNormal(n);
     774                                        delete n;                                       
     775                                        /*intvec *n = new intvec(this->numVars);
    714776                                        makeInt(P,jj,*n);                                               
    715777                                        codim2Act->setFacetNormal(n);
    … …  
    717779                                        codim2Act->next = new facet(2);
    718780                                        codim2Act = codim2Act->next;                                           
    719                                         delete n;                                                                       
     781                                        delete n;*/
    720782                                }                                                                               
    721783                                       
    … …  
    15861648                        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
    15871649                        facet *SearchListAct;
    1588                         SearchListAct = SearchListRoot;
     1650                        //SearchListAct = SearchListRoot;
    15891651                       
    15901652                        gcone *gcAct;
    … …  
    15991661                        facet *fAct;
    16001662                        fAct = gcAct->facetPtr;                 
     1663                       
     1664                        ring rAct;
     1665                        rAct = currRing;
    16011666                                               
    16021667                        int UCNcounter=gcAct->getUCN();
    … …  
    16161681                        we must not memcpy them before these ops!
    16171682                        */
    1618                         while(fAct->next!=NULL)
    1619                         {                               
    1620                                 memcpy(SearchListAct,fAct,sizeof(facet));
    1621                                 SearchListAct->next = new facet();
    1622                                 SearchListAct = SearchListAct->next;
     1683                       
     1684                        for(int ii=0;ii<this->numFacets;ii++)
     1685                        {
     1686                                if(ii==0)
     1687                                {
     1688                                        //facet *SearchListRoot = new facet();
     1689                                        SearchListAct = SearchListRoot;
     1690                                        memcpy(SearchListAct,fAct,sizeof(facet));                                       
     1691                                }
     1692                                else
     1693                                {
     1694                                        SearchListAct->next = new facet();
     1695                                        SearchListAct = SearchListAct->next;
     1696                                        memcpy(SearchListAct,fAct,sizeof(facet));
     1697                                }
    16231698                                fAct = fAct->next;
    16241699                        }
     1700                       
    16251701                        SearchListAct = SearchListRoot; //Set to beginning of list
    16261702                       
    … …  
    16341710                        We always choose the first facet from fListPtr as facet to be flipped
    16351711                        */
    1636                         while(SearchListAct->next!=NULL)
     1712                        while(SearchListAct!=NULL)
    16371713                        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
    16381714                                fAct = SearchListAct;
    16391715                                //while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) )
    1640                                 do
     1716                                //do
     1717                                while(fAct!=NULL)
    16411718                                {
    16421719                                        gcAct->flip(gcAct->gcBasis,fAct);
    … …  
    16501727                                        /*add facets to SLA here*/
    16511728                                        gcTmp->enqueueNewFacets(*SearchListRoot);
     1729                                        gcTmp->showSLA(*SearchListRoot);
    16521730                                        rChangeCurrRing(gcAct->baseRing);
    16531731                                        gcPtr->next=gcTmp;
    … …  
    16591737                                        else
    16601738                                                break;
    1661                                 }while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
     1739                                }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
     1740                                //Search for cone with smallest UCN
    16621741                                gcNext = gcHead;
    1663                                 while(gcNext->next!=NULL)
    1664                                 {
    1665                                         if(gcNext->getUCN() > gcNext->next->getUCN() )
    1666                                         {
    1667                                                 gcAct = gcNext->next;
     1742                                while(gcNext!=NULL)
     1743                                {
     1744                                        if( gcNext->getUCN() == UCNcounter+1 )
     1745                                        {//NOTE THIS IS BUGGY. Apparently changes to the wrong ring
     1746                                                gcAct = gcNext;
     1747                                                rAct=rCopy(gcAct->baseRing);
     1748                                                rComplete(rAct);
     1749                                                rChangeCurrRing(rAct);                                         
     1750                                                break;                                         
    16681751                                        }
    16691752                                        gcNext = gcNext->next;
    16701753                                }
    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 //                              }
    1684                                
     1754                                UCNcounter++;
     1755                                //SearchListAct = SearchListAct->next;
     1756                                SearchListAct = fAct->next;                             
    16851757                        }
    16861758                                               
    … …  
    18031875                        intvec *n = new intvec(this->numVars);
    18041876                       
    1805                         while(codim2Act->next!=NULL)
     1877                        //while(codim2Act->next!=NULL)
     1878                        while(codim2Act!=NULL)
    18061879                        {                               
    18071880                                n=codim2Act->getFacetNormal();
    … …  
    18141887                        //delete n;
    18151888                        codim2Act = this->facetPtr->codim2Ptr;  //reset to start of linked list
    1816                         while(codim2Act->next!=NULL)
     1889                        //while(codim2Act->next!=NULL)
     1890                        while(codim2Act!=NULL)
    18171891                        {
    18181892                                //intvec *n = new intvec(this->numVars);
    … …  
    18421916                        facet *codim2Act;
    18431917                        codim2Act = this->facetPtr->codim2Ptr;
     1918                        bool doNotAdd=FALSE;
    18441919                        while(slEnd->next!=NULL)
    18451920                        {
    … …  
    18501925                        intvec *f2Normal = new intvec(this->numVars);
    18511926                        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                                         {                                               
     1927                        intvec *sl2Normal = new intvec(this->numVars);                 
     1928                        while(fAct!=NULL)
     1929                        {
     1930                                doNotAdd=FALSE;
     1931                                fNormal = fAct->getFacetNormal();
     1932                                slAct = &f;     //return to start of list
     1933                                while(slAct!=NULL)
     1934                                {
     1935                                        slNormal = slAct->getFacetNormal();
     1936                                        /*if(!isParallel(fNormal,slNormal))
     1937                                        {
    18621938                                                slEnd->next = new facet();
    18631939                                                slEnd = slEnd->next;
    18641940                                                slEnd->setUCN(this->getUCN());
    1865                                                 slEnd->setFacetNormal(fNormal);                                         
    1866                                                 //memcpy(slEnd,fAct,sizeof(facet));
     1941                                                slEnd->setFacetNormal(fNormal);
     1942                                                break;
     1943                                                                                       
    18671944                                        }
    1868                                         fAct = fAct->next;
    1869                                 }
    1870                                 slAct = slAct->next;
     1945                                        else
     1946                                        {
     1947                                                //NOTE check codim2facets here
     1948                                                break;  //hopefully breaks the while loop
     1949                                        }*/
     1950                                        if(isParallel(fNormal,slNormal))
     1951                                        {
     1952                                                //NOTE check codim2facets here
     1953                                                doNotAdd=TRUE;
     1954                                                break;
     1955                                        }
     1956                                        slAct = slAct->next;
     1957                                }
     1958                                if(doNotAdd==FALSE)
     1959                                {
     1960                                        slEnd->next = new facet();
     1961                                        slEnd = slEnd->next;
     1962                                        slEnd->setUCN(this->getUCN());
     1963                                        slEnd->setFacetNormal(fNormal);
     1964                                }
     1965                                fAct = fAct->next;
    18711966                        }
    18721967//                      delete sl2Normal;
    … …  
    20072102};//class gcone
    20082103
    2009 int gcone::UCN=0;
     2104int gcone::counter=0;
    20102105
    20112106ideal gfan(ideal inputIdeal)
Note: See TracChangeset for help on using the changeset viewer.