Changeset 2126c0d in git for kernel/gfan.cc


Ignore:
Timestamp:
Jun 29, 2009, 4:46:58 PM (15 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
7379a771e910d618071fc586a19837c93dd9e01e
Parents:
832b70058672520c5e80a49ba96723f3a316ad73
Message:
Fixed bug in ring construction! Works on most 2D-examples except for some esoteric ideals like
<x-1,y-1> whose restricted cone is just the pos orth. Easy to fix.
However code hangs on 3D-examples like <x+y+z,x3z+x+y2> which is due to the enqueueNewFacets not
working as it should.
Some minor cleanup


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r832b70 r2126c0d  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-06-25 08:52:03 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.67 2009-06-25 08:52:03 monerjan Exp $
    6 $Id: gfan.cc,v 1.67 2009-06-25 08:52:03 monerjan Exp $
     4$Date: 2009-06-29 14:46:58 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.68 2009-06-29 14:46:58 monerjan Exp $
     6$Id: gfan.cc,v 1.68 2009-06-29 14:46:58 monerjan Exp $
    77*/
    88
    … …  
    397397                        this->gcBasis=idCopy(f.flipGB);
    398398                        this->inputIdeal=idCopy(this->gcBasis);
    399                         this->baseRing=rCopy0(f.flipRing);
     399                        this->baseRing=rCopy(f.flipRing);
    400400                        this->numFacets=0;
    401401                        //rComplete(this->baseRing);
    … …  
    587587
    588588#ifdef gfan_DEBUG
    589                         cout << "The inequality matrix is" << endl;
    590                         dd_WriteMatrix(stdout, ddineq);
     589                        cout << "The inequality matrix is" << endl;
     590                        dd_WriteMatrix(stdout, ddineq);
    591591#endif
    592592
    … …  
    608608                        ddcols = ddineq->colsize;
    609609#ifdef gfan_DEBUG
    610                         cout << "Having removed redundancies, the normals now read:" << endl;
    611                         dd_WriteMatrix(stdout,ddineq);
    612                         cout << "Rows = " << ddrows << endl;
    613                         cout << "Cols = " << ddcols << endl;
     610                        cout << "Having removed redundancies, the normals now read:" << endl;
     611                        dd_WriteMatrix(stdout,ddineq);
     612                        cout << "Rows = " << ddrows << endl;
     613                        cout << "Cols = " << ddcols << endl;
    614614#endif
    615615                       
    … …  
    744744                                       
    745745#ifdef gfan_DEBUG
    746 //                              dd_WriteMatrix(stdout,ddakt);
     746//                              cout << "Codim2 matrix"<<endl;
     747//                              dd_WriteMatrix(stdout,ddakt);
    747748#endif
    748749                                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    749750                                P=dd_CopyGenerators(ddpolyh);
    750751#ifdef gfan_DEBUG
    751 //                              dd_WriteMatrix(stdout,P);
     752//                              dd_WriteMatrix(stdout,P);
    752753#endif
    753754                                       
    … …  
    807808                        std::cout << "===" << std::endl;
    808809                        std::cout << "running gcone::flip" << std::endl;
    809 //                      std::cout << "fNormal=";
    810 //                      fNormal->show();
    811 //                      std::cout << std::endl;
     810                        std::cout << "flipping" << endl;
     811                        for(int ii=0;ii<IDELEMS(gb);ii++)
     812                        {
     813                                pWrite((poly)gb->m[ii]);
     814                        }
     815                        cout << "over facet" << endl;
     816                        fNormal->show();
     817                        std::cout << std::endl;
    812818#endif                         
    813819                        /*1st step: Compute the initial ideal*/
    … …  
    902908                        ina=idrCopyR(initialForm,srcRing);                     
    903909#ifdef gfan_DEBUG
    904                         cout << "ina=";
    905                         idShow(ina); cout << endl;
     910//                      cout << "ina=";
     911//                      idShow(ina); cout << endl;
    906912#endif
    907913                        ideal H;
    … …  
    910916                        idSkipZeroes(H);
    911917#ifdef gfan_DEBUG
    912                         cout << "H="; idShow(H); cout << endl;
     918//                      cout << "H="; idShow(H); cout << endl;
    913919#endif
    914920                        /*Substep 2.2
    … …  
    920926                        srcRing_H=idrCopyR(H,tmpRing);
    921927#ifdef gfan_DEBUG
    922                         cout << "srcRing_H = ";
    923                         idShow(srcRing_H); cout << endl;
     928//                      cout << "srcRing_H = ";
     929//                      idShow(srcRing_H); cout << endl;
    924930#endif
    925931                        srcRing_HH=ffG(srcRing_H,this->gcBasis);               
    926932#ifdef gfan_DEBUG
    927                         cout << "srcRing_HH = ";
    928                         idShow(srcRing_HH); cout << endl;
     933//                      cout << "srcRing_HH = ";
     934//                      idShow(srcRing_HH); cout << endl;
    929935#endif
    930936                        /*Substep 2.2.1
    … …  
    10331039                                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
    10341040                        }
    1035                         dd_WriteMatrix(stdout,intPointMatrix);
     1041#ifdef gfan_DEBUG
     1042//                      dd_WriteMatrix(stdout,intPointMatrix);
     1043#endif
    10361044                        intvec *iv_weight = new intvec(this->numVars);
    10371045                        interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
    … …  
    10811089                        ring dstRing=rCopy0(tmpRing);
    10821090                        int length=iv_weight->length();
    1083                         int *A=(int *)omAlloc(length*sizeof(int));
     1091                        int *A=(int *)omAlloc0(length*sizeof(int));
    10841092                        for(int jj=0;jj<length;jj++)
    10851093                        {
    … …  
    11161124                       
    11171125                        f->setFlipGB(dstRing_I);//store the flipped GB
    1118                         f->flipRing=rCopy0(dstRing);    //store the ring on the other side
     1126                        f->flipRing=rCopy(dstRing);     //store the ring on the other side
    11191127#ifdef gfan_DEBUG
    11201128                        cout << "Flipped GB is: " << endl;
    … …  
    11341142                poly restOfDiv(poly const &f, ideal const &I)
    11351143                {
    1136                         cout << "Entering restOfDiv" << endl;
     1144//                      cout << "Entering restOfDiv" << endl;
    11371145                        poly p=f;
    11381146                        //pWrite(p);
    … …  
    11961204                ideal ffG(ideal const &H, ideal const &G)
    11971205                {
    1198                         cout << "Entering ffG" << endl;
     1206//                      cout << "Entering ffG" << endl;
    11991207                        int size=IDELEMS(H);
    12001208                        ideal res=idInit(size,1);
    … …  
    17251733                                SearchListAct = fAct->next;                             
    17261734                        }
    1727                                                
    1728 //                      while(SearchListAct->next!=NULL)
    1729 //                      {
    1730 //                              fAct = SearchListAct;
    1731 //                              do
    1732 //                              {
    1733 //                                      gcAct->flip(gcAct->gcBasis,fAct);
    1734 //                                      ring rTmp = rCopy(fAct->flipRing);
    1735 //                                      rComplete(rTmp);
    1736 //                                      rChangeCurrRing(rTmp);
    1737 //                                      gcone *gcTmp = new gcone(*gcAct,*fAct);
    1738 //                                      gcTmp->getConeNormals(gcTmp->gcBasis,FALSE);
    1739 //                                      gcTmp->getCodim2Normals(*gcTmp);
    1740 //                                      gcPtr->next = gcTmp;
    1741 //                                      gcPtr = gcPtr->next;
    1742 //                                      /*add facets to SLA here*/
    1743 //                                      rChangeCurrRing(gcAct->baseRing);
    1744 //                                      fAct = fAct->next;
    1745 //                              }while( (fAct->next!=NULL) &&  (fAct->getUCN()==fAct->next->getUCN()));
    1746 //                              SearchListAct = SearchListAct->next;//fAct;                             
    1747 //                      }
    1748                        
    1749 //                      while(SearchListAct->next!=NULL)
    1750 //                      {
    1751 //                              gcAct->flip(gcAct->gcBasis,SearchListAct);
    1752 //                              ring rTmp=rCopy(fAct->flipRing);
    1753 //                              rComplete(rTmp);
    1754 //                              rChangeCurrRing(rTmp);
    1755 //                              gcone *gcTmp = new gcone(*gcAct,*SearchListAct);
    1756 //                              gcTmp->getConeNormals(gcTmp->gcBasis,FALSE);
    1757 //                              gcTmp->getCodim2Normals(*gcTmp);
    1758 //                              rChangeCurrRing(gcAct->baseRing);
    1759 //                              SearchListAct = SearchListAct->next;
    1760 //                             
    1761 //                              //gcone *gcNew = new gcone(SearchListAct->flipRing,SearchListAct->flipGB);
    1762 //                              //gcAct = gcNew;
    1763 //                      }
    17641735               
    17651736                        //NOTE Hm, comment in and get a crash for free...
    … …  
    18021773                        }
    18031774                       
     1775                        //Check whether denom is all ones, in which case we will divide out the gcd of the nominators
     1776//                      mpz_t checksum; mpz_t rop;
     1777//                      mpz_init(checksum);
     1778//                      mpz_init(rop);
     1779//                      bool divideOutGcd=FALSE;
     1780//                      for(int ii=0;ii<this->numVars;ii++)
     1781//                      {
     1782//                              mpz_add(rop, checksum, denom[ii]);
     1783//                              mpz_set(checksum, rop);
     1784//                      }
     1785//                      if( (int)mpz_get_ui(checksum)==this->numVars)
     1786//                      {
     1787//                              divideOutGcd=TRUE;
     1788//                      }
     1789                       
    18041790                        /*Compute lcm of the denominators*/
    18051791                        mpz_set(tmp,denom[0]);
    … …  
    18331819                }
    18341820                /**
    1835                  * We compute the gcd of the components of the codim-2-facets and
    1836                  * multiply the each codim-2facet by it. Thus we get a normalized representation of each
    1837                  * (codim-2)-facet normal.
     1821                 * For all codim-2-facets we compute the gcd of the components of the facet normal and
     1822                 * divide it out. Thus we get a normalized representation of each
     1823                 * (codim-2)-facet normal, i.e. a primitive vector
    18381824                */
    18391825                void normalize()
    … …  
    18521838                                        ggT = intgcd(ggT,(*n)[ii]);
    18531839                                }
     1840                                for(int ii=0;ii<this->numVars;ii++)
     1841                                {
     1842                                        (*n)[ii] = ((*n)[ii])/ggT;
     1843                                }
    18541844                                codim2Act = codim2Act->next;                           
    18551845                        }
    18561846                        //delete n;
    1857                         codim2Act = this->facetPtr->codim2Ptr;  //reset to start of linked list
    1858                         //while(codim2Act->next!=NULL)
     1847                        /*codim2Act = this->facetPtr->codim2Ptr;        //reset to start of linked list                 
    18591848                        while(codim2Act!=NULL)
    1860                         {
    1861                                 //intvec *n = new intvec(this->numVars);
     1849                        {                               
    18621850                                n=codim2Act->getFacetNormal();
    18631851                                intvec *new_n = new intvec(this->numVars);
    … …  
    18701858                                //delete n;
    18711859                                //delete new_n;
    1872                         }                      
     1860                        }       */             
    18731861                }
    18741862                /**
    … …  
    18851873                        facet *codim2Act;
    18861874                        codim2Act = this->facetPtr->codim2Ptr;
     1875                        facet *sl2Act;
     1876                        sl2Act = f.codim2Ptr;
    18871877                        bool doNotAdd=FALSE;
    18881878                        while(slEnd->next!=NULL)
    … …  
    19021892                                while(slAct!=NULL)
    19031893                                {
    1904                                         slNormal = slAct->getFacetNormal();
    1905                                         /*if(!isParallel(fNormal,slNormal))
    1906                                         {
    1907                                                 slEnd->next = new facet();
    1908                                                 slEnd = slEnd->next;
    1909                                                 slEnd->setUCN(this->getUCN());
    1910                                                 slEnd->setFacetNormal(fNormal);
    1911                                                 break;
    1912                                                                                        
    1913                                         }
    1914                                         else
    1915                                         {
    1916                                                 //NOTE check codim2facets here
    1917                                                 break;  //hopefully breaks the while loop
    1918                                         }*/
     1894                                        slNormal = slAct->getFacetNormal();                                     
     1895                                        /*If the normals are parallel we check whether the
     1896                                        codim-2-normals coincide as well*/
    19191897                                        if(isParallel(fNormal,slNormal))
    19201898                                        {
    19211899                                                //NOTE check codim2facets here
    1922                                                 doNotAdd=TRUE;
     1900//                                              while(codim2Act!=NULL)
     1901//                                              {
     1902//                                                      f2Normal = codim2Act->getFacetNormal();
     1903//                                                      sl2Act = f.codim2Ptr;
     1904//                                                      while(sl2Act!=NULL)
     1905//                                                      {
     1906//                                                              sl2Normal = sl2Act->getFacetNormal();
     1907//                                                              if(!isParallel(f2Normal,sl2Normal))
     1908//                                                              {
     1909//                                                                      doNotAdd=FALSE;
     1910//                                                                      break;
     1911//                                                              }
     1912//                                                              sl2Act = sl2Act->next;
     1913//                                                      }
     1914//                                                      if(doNotAdd==FALSE)
     1915//                                                              break;
     1916//                                                      codim2Act = codim2Act->next;
     1917//                                                     
     1918//                                              }
     1919                                                doNotAdd=FALSE;
    19231920                                                break;
    19241921                                        }
    … …  
    19751972                       
    19761973                        int jj=0;
    1977                         while(fAct->next!=NULL)
     1974                        //while(fAct->next!=NULL)
     1975                        while(fAct!=NULL)
    19781976                        {
    19791977                                intvec *comp;
    … …  
    20372035                               
    20382036                                gcOutputFile << "FACETS" << endl;                                                               
    2039                                 while(fAct->next!=NULL)
     2037                                //while(fAct->next!=NULL)
     2038                                while(fAct!=NULL)
    20402039                                {       
    20412040                                        intvec *iv = new intvec(gc.numVars);
Note: See TracChangeset for help on using the changeset viewer.