Changeset 276932 in git


Ignore:
Timestamp:
Feb 9, 2010, 4:18:18 PM (13 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
d02c1bdf6308fb7d278cd3aa9dfdfb014a07d639
Parents:
8e77eb857f444e4cf338e91ad0e4dc348c8b8a96
Message:
options.h
isParallel vs. intvec::compare in gcone::areEqual(facet, facet). The
latter is supposed to be better but miraculously causes a segfault in
noRevS while trying to change into a apparently rDeleted ring...
conditional #elifs
conditional appending of linearity space (might be 0)
interiorPoint2 removed
Ring changes in lPrepareResult
gcone::f2M casts int64vec to intvec due to LIST_CMD issues


git-svn-id: file:///usr/local/Singular/svn/trunk@12536 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
kernel
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r8e77eb r276932  
    1010
    1111#ifdef HAVE_GFAN
    12 
     12#include "options.h"
    1313#include "kstd1.h"
    1414#include "kutil.h"      //ksCreateSpoly
    … …  
    6868
    6969#ifndef gfan_DEBUG
    70 // #define gfan_DEBUG
     70#define gfan_DEBUG
    7171#ifndef gfan_DEBUGLEVEL
    7272#define gfan_DEBUGLEVEL 1
    … …  
    193193                idDelete((ideal *)&this->flipGB);
    194194        if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
    195                 rDelete(this->flipRing);
     195//              rDelete(this->flipRing); //See vol II/134
    196196//      this->flipRing=NULL;
    197197        this->prev=NULL;
    198198        this->next=NULL;
    199         //this=NULL;
    200 }
    201                
    202                
    203 /** \brief Comparison of facets*/
     199}
     200               
     201               
     202/** \brief Comparison of facets
     203 * called from enqueueNewFacets
     204* The facet normals are primitve vectors since we call gcone::normalize() on all cones.
     205* Hence it should suffice to check whether facet normal f equals minus facet normal s.
     206* If so we check the extremal rays
     207*/
    204208inline bool gcone::areEqual(facet *f, facet *s)
    205209{
    … …  
    216220        //Do not need parallelity. Too time consuming
    217221        if(!isParallel(*fNormal,*sNormal))
     222//      if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
    218223                notParallelCtr++;
    219 //      int64vec *foo=ivNeg(sNormal);
    220 //      if(fNormal->compare(foo)!=0) //facet normals
    221 //              return TRUE;   
    222         else//parallelity, so we check the codim2-facets
     224//      else//parallelity, so we check the codim2-facets
     225//      if(isParallel(*fNormal,*sNormal))
     226//      if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
    223227        {
    224228                facet* f2Act;
    … …  
    258262        }
    259263        return res;
     264//      int64vec *foo=ivNeg(sNormal);
     265//      if(fNormal->compare(foo)!=0) //facet normals
     266//      {
     267//              delete foo;
     268//              res=FALSE;
     269//      }
     270//      else
     271//      {
     272//              facet* f2Act;
     273//              facet* s2Act;
     274//              f2Act = f->codim2Ptr;           
     275//              ctr=0;
     276//              while(f2Act!=NULL)
     277//              {
     278//                      int64vec* f2Normal;
     279//                      f2Normal = f2Act->getFacetNormal();
     280//                      s2Act = s->codim2Ptr;
     281//                      while(s2Act!=NULL)
     282//                      {
     283//                              int64vec* s2Normal;
     284//                              s2Normal = s2Act->getFacetNormal();
     285// //                           bool foo=areEqual(f2Normal,s2Normal);
     286//                              int foo=f2Normal->compare(s2Normal);
     287//                              if(foo==0)
     288//                                      ctr++;
     289//                              s2Act = s2Act->next;
     290//                              delete s2Normal;
     291//                      }
     292//                      delete f2Normal;
     293//                      f2Act = f2Act->next;
     294//              }
     295//      }               
     296//      delete fNormal;
     297//      delete sNormal;
     298//      if(ctr==f->numCodim2Facets)
     299//              res=TRUE;
     300//      return res;
    260301}       
    261302               
    … …  
    624665                        res = FALSE;
    625666                        break;
    626                 }                                               
     667                }
    627668        }
    628669        return res;
    … …  
    11991240        dd_MatrixPtr ddineq;
    12001241        dd_ErrorType err;
    1201         ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
     1242        if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
     1243                ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
     1244        else
     1245                ddineq = dd_CopyMatrix(gc.ddFacets);
    12021246        dd_PolyhedraPtr ddPolyh;
    12031247        ddPolyh = dd_DDMatrix2Poly(ddineq, &err);
    … …  
    16651709                       
    16661710        f->setFlipGB(dstRing_I);//store the flipped GB
    1667         idDelete(&dstRing_I);
     1711//      idDelete(&dstRing_I);
    16681712        f->flipRing=rCopy(dstRing);     //store the ring on the other side
    1669 //#ifdef gfan_DEBUG
    1670 //      cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
    1671 //      this->idDebugPrint(dstRing_I);
    1672 //      cout << endl;
    1673 //#endif                       
     1713#ifdef gfan_DEBUG
     1714        cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
     1715        this->idDebugPrint(dstRing_I);
     1716//      idPrint(dstRing_I);
     1717        cout << endl;
     1718#endif 
     1719        idDelete(&dstRing_I);   
    16741720        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
    16751721        rDelete(dstRing);
    … …  
    20232069}//void interiorPoint(dd_MatrixPtr const &M)
    20242070
    2025 inline void gcone::interiorPoint2(const dd_MatrixPtr &M, int64vec &iv)
    2026 {
    2027         KMatrix<Rational> mMat(M->rowsize+1,M->colsize);
    2028         for(int ii=0;ii<M->rowsize;ii++)
    2029         {
    2030                 for(int jj=0;jj<M->colsize-1;jj++)
    2031                 {
    2032                         if(mpq_sgn(M->matrix[ii][jj+1])<-1)
    2033                         {
    2034                                 mMat.set(ii,jj,-(Rational)mpq_get_d(M->matrix[ii][jj+1]));
    2035                         }
    2036                         else
    2037                                 mMat.set(ii,jj,(Rational)mpq_get_d(M->matrix[ii][jj+1]));
    2038                                
    2039 //                      mMat.set(ii,jj,&(M->matrix[ii][jj+1]) );
    2040                         cout << mpq_get_d(M->matrix[ii][jj+1]) << " ";
    2041 //                      int val=(int)mMat.get(ii,jj);
    2042 //                      cout << ii << "," << jj << endl;;
    2043 //                      mpq_out_str (NULL, 10, (__mpq_struct)mMat.get(ii,jj));
    2044                 }
    2045                 cout << endl;
    2046                 mMat.set(ii,M->colsize-1,1);
    2047         }
    2048         dd_WriteMatrix(stdout,M);
    2049 //      for(int ii=0;ii<M->rowsize;ii++)
    2050 //      {
    2051 //              cout << mMat.get(ii,ii+M->colsize) << " ";
    2052 //              if((ii+M->colsize)%M->colsize==0)
    2053 //                      cout << endl;
    2054 //      }
    2055        
    2056         Rational* mSol;
    2057         int rank;
    2058         int c;
    2059 //      dd_WriteMatrix(stdout,M);
    2060         rank=mMat.solve(&mSol,&c);
    2061 //      for(int ii=0;ii<c;ii++)
    2062 //              iv[ii]=mSol[ii];
    2063 //              cout << mSol[ii].get_den() << "/" << mSol[ii].get_num() << endl;
    2064         int gcd=1;
    2065         for(int ii=0;ii<c-1;ii++)
    2066                 gcd += intgcd(mSol[ii].get_den_si(),mSol[ii+1].get_den_si());
    2067         cout << gcd << endl;
    2068         for(int ii=0;ii<iv.length();ii++)
    2069                 iv[ii]=(int)((mSol[ii].get_num_si()*gcd)/mSol[ii].get_den_si());
    2070 
    2071 }
    2072        
     2071
    20732072       
    20742073/** \brief Copy a ring and add a weighted ordering in first place
    … …  
    21502149/** \brief Compute the lineality/homogeneity space
    21512150* It is the kernel of the inequality matrix Ax=0
     2151* As a result gcone::dd_LinealitySpace is set
    21522152*/
    21532153dd_MatrixPtr gcone::computeLinealitySpace()
    … …  
    22962296        int UCNcounter=gcAct->getUCN();
    22972297       
    2298         /* Use pure SLA or writeCone2File */
    2299 //      enum heuristic {
    2300 //              ram,
    2301 //              hdd
    2302 //      };
    2303 //      heuristic h;
    2304 //      h=hdd;
    2305        
    23062298#ifdef gfan_DEBUG
    23072299        cout << "NoRevs" << endl;
    … …  
    23102302#endif                 
    23112303        /*Compute lineality space here
    2312         * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set
    2313         */
     2304        * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/
    23142305        dd_LinealitySpace = gcAct->computeLinealitySpace();
    23152306        /*End of lineality space computation*/         
    … …  
    23172308        gcAct->getExtremalRays(*gcAct);
    23182309                               
    2319         //Compute unique representation of Facets and rays
     2310        //Compute unique representation of Facets and rays, i.e. primitive vectors
    23202311        gcAct->normalize();
    23212312                       
    23222313        /* Make a copy of the facet list of first cone
    23232314           Since the operations getCodim2Normals and normalize affect the facets
    2324            we must not memcpy them before these ops!
    2325         */
    2326        
    2327         facet *codim2Act; codim2Act = NULL;                     
    2328         facet *sl2Root; //sl2Root = new facet(2);
    2329         facet *sl2Act;  //sl2Act = sl2Root;
    2330                        
     2315           we must not memcpy them before these ops! */
     2316        /*facet *codim2Act; codim2Act = NULL;                   
     2317        facet *sl2Root;
     2318        facet *sl2Act;*/                       
    23312319        for(int ii=0;ii<this->numFacets;ii++)
    23322320        {
    … …  
    23382326                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
    23392327                        {                                               
    2340                                 SearchListAct = SearchListRoot;                        
     2328                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
    23412329                        }
    23422330                        else
    … …  
    23492337                        SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
    23502338                        SearchListAct->isFlippable=TRUE;
    2351                         //Copy codim2-facets                           
     2339                        //Copy codim2-facets
     2340                        facet *codim2Act;
     2341                        codim2Act = NULL;                       
     2342                        facet *sl2Root;
     2343                        facet *sl2Act;                 
    23522344                        codim2Act=fAct->codim2Ptr;
    23532345                        SearchListAct->codim2Ptr = new facet(2);
    … …  
    23922384        SearchListAct = SearchListRoot; //Set to beginning of List
    23932385       
    2394         fAct = gcAct->facetPtr;
    2395         //NOTE Disabled until code works as expected. Reenabled 2.11.2009
     2386//      fAct = gcAct->facetPtr;//???
    23962387        gcAct->writeConeToFile(*gcAct);                 
    23972388        /*End of initialisation*/
    … …  
    23992390        fAct = SearchListAct;
    24002391        /*2nd step
    2401           Choose a facet from fListPtr, flip it and forget the previous cone
    2402           We always choose the first facet from fListPtr as facet to be flipped
     2392          Choose a facet from SearchList, flip it and forget the previous cone
     2393          We always choose the first facet from SearchList as facet to be flipped
    24032394        */                     
    2404         while((SearchListAct!=NULL))// && counter<10)
     2395        while((SearchListAct!=NULL) && gcone::counter<5)
    24052396        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
    24062397                fAct = SearchListAct;
    … …  
    24422433                                gcTmp->writeConeToFile(*gcTmp);
    24432434                                idDelete((ideal*)&gcTmp->gcBasis);//Whonder why?
    2444                                 //If you use the following make sure it is uncommented in readConeFromFile
    24452435                                rDelete(gcTmp->baseRing);
    24462436                        }                       
    … …  
    24552445                        gcPtr->next=gcTmp;
    24562446                        gcPtr=gcPtr->next;
     2447                        //Cleverly disguised exit condition follows
    24572448                        if(fAct->getUCN() == fAct->next->getUCN())
    24582449                        {
    … …  
    24682459  #if SIZEOF_LONG==8    //64 bit
    24692460                while(gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && SearchListRoot!=NULL)
    2470   #elif SIZEOF_LONG==4
     2461  #elif SIZEOF_LONG == 4
    24712462                while(gcNext!=(gcone * const)0xfbfbfbfb && SearchListRoot!=NULL)
    24722463  #endif
    … …  
    24972488                                        readConeFromFile(gcAct->getUCN(), gcAct);
    24982489                                        //Kill the baseRing but ONLY if it is not the ring the computation started in!
    2499                                         if(gcDel->getUCN()!=1)
    2500                                                 rDelete(gcDel->baseRing);
     2490//                                      if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
     2491//                                              rDelete(gcDel->baseRing);
    25012492//                                      rAct=rCopy(gcAct->baseRing);
    25022493                                        /*The ring change occurs in readConeFromFile, so as to
    … …  
    25072498                                        break;
    25082499                                }                               
    2509                         }                       
     2500                        }
     2501                        /*else
     2502                        {
     2503                                if(gfanHeuristic==1)
     2504                                {
     2505                                        gcone *gcDel;
     2506                                        gcDel = gcNext;
     2507                                        if(gcDel->getUCN()!=1)
     2508                                                rDelete(gcDel->baseRing);
     2509                                }                                       
     2510                        }*/                     
    25102511                        gcNext = gcNext->next;
    25112512                }
    … …  
    27092710                                while(fCopy!=NULL)
    27102711                                {
    2711                                         if(slAct==NULL)
     2712                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
    27122713                                        {
    2713                                                 slAct = new facet(*fCopy);
    2714                                                 slHead = slAct;                                                         
    2715                                         }
    2716                                         else
    2717                                         {
    2718                                                 slAct->next = new facet(*fCopy);
    2719                                                 slAct = slAct->next;
     2714                                                if(slAct==NULL)
     2715                                                {
     2716                                                        slAct = new facet(*fCopy);
     2717                                                        slHead = slAct;                                                         
     2718                                                }
     2719                                                else
     2720                                                {
     2721                                                        slAct->next = new facet(*fCopy);
     2722                                                        slAct = slAct->next;
     2723                                                }
    27202724                                        }
    27212725                                        fCopy = fCopy->next;
    27222726                                }
    2723                                 break;
     2727                                break;//Where does this lead to?
    27242728                        }
    27252729                        /*End of dumping into SLA*/
    … …  
    28452849facet * gcone::enqueue2(facet *f)
    28462850{
    2847 //      vector<facet> searchList;
    2848 //      facet fAct;
     2851//      vector<facet*> searchList;
     2852//      facet *fAct;
    28492853//      fAct=this->facetPtr;
    28502854//      facet *slAct;
    28512855//      slAct = f;
     2856//      facet *slHead;
    28522857//      while(fAct!=NULL)
    28532858//      {
    28542859//              if(fAct->isFlippable==TRUE)
    28552860//              {
    2856 //                      if(slAct==NULL)
     2861//                      if(searchList.size()==0)
     2862//                      {
     2863//                              facet *fCopy;
     2864//                              fCopy=fAct;
     2865//                              while(fCopy!=NULL)
     2866//                              {
     2867//                                      if(fCopy->isFlippable==TRUE)
     2868//                                              searchList.push_back(fCopy);
     2869//                                      fCopy=fCopy->next;
     2870//                              }
     2871//                      }
    28572872//              }
    28582873//      }
    … …  
    32173232* allows us to rDelete(gcDel->baseRing) and the such in gcone::noRevS.
    32183233*\param *gc Pointer to gcone, preferably gcRoot ;-)
    3219 *\param n the number of cones
     3234*\param n the number of cones as determined by gcRoot->getCounter()
    32203235*
    32213236*/
    … …  
    32313246        for(int ii=0;ii<n;ii++)
    32323247        {
     3248                if(gfanHeuristic==1)// && gcAct->getUCN()>1)
     3249                {
     3250                        gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
     3251                        rChangeCurrRing(gcAct->getBaseRing());
     3252                }               
    32333253                res->m[ii].rtyp=LIST_CMD;
    32343254                lists l=(lists)omAllocBin(slists_bin);
    … …  
    32423262                */
    32433263//              if(gcAct->gcBasis->m[0]==(poly)NULL)
    3244                 if(gfanHeuristic==1)
    3245                         gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
     3264//              if(gfanHeuristic==1 && gcAct->getUCN()>1)
     3265//                      gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
    32463266//              ring saveRing=currRing;
    32473267//              ring tmpRing=gcAct->getBaseRing;
    … …  
    32523272
    32533273                l->m[2].rtyp=INTVEC_CMD;
    3254                 int64vec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
    3255                 l->m[2].data=(void*)iv64Copy(&iv);
     3274                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
     3275                l->m[2].data=(void*)ivCopy(&iv);
    32563276               
    32573277                l->m[3].rtyp=LIST_CMD;
    … …  
    32633283                        {
    32643284                                lCodim2List->m[jj].rtyp=INTVEC_CMD;
    3265                                 int64vec ivC2=(gcAct->f2M(gcAct,fAct,2));
    3266                                 lCodim2List->m[jj].data=(void*)iv64Copy(&ivC2);
     3285                                intvec ivC2=(gcAct->f2M(gcAct,fAct,2));
     3286                                lCodim2List->m[jj].data=(void*)ivCopy(&ivC2);
    32673287                                jj++;
    32683288                                fAct = fAct->next;
    … …  
    32823302* f should always point to gc->facetPtr
    32833303* param n is used to determine whether it operates in codim 1 or 2
    3284 *
     3304* We have to cast the int64vecs to intvec due to issues with list structure
    32853305*/
    3286 inline int64vec gcone::f2M(gcone *gc, facet *f, int n)
     3306inline intvec gcone::f2M(gcone *gc, facet *f, int n)
    32873307{
    32883308        facet *fAct;
    3289         int64vec *res;//=new int64vec(this->numVars);   
     3309        intvec *res;//=new int64vec(this->numVars);     
    32903310//      int codim=n;
    32913311//      int bound;
    … …  
    32933313        if(n==1)
    32943314        {
    3295                 int64vec *m1Res=new int64vec(gc->numFacets,gc->numVars,0);
    3296                 res = iv64Copy(m1Res);
     3315                intvec *m1Res=new intvec(gc->numFacets,gc->numVars,0);
     3316                res = ivCopy(m1Res);
    32973317                fAct = gc->facetPtr;
    32983318                delete m1Res;
    … …  
    33023322        {
    33033323                fAct = f->codim2Ptr;
    3304                 int64vec *m2Res = new int64vec(f->numCodim2Facets,gc->numVars,0);
    3305                 res = iv64Copy(m2Res);
     3324                intvec *m2Res = new intvec(f->numCodim2Facets,gc->numVars,0);
     3325                res = ivCopy(m2Res);
    33063326                delete m2Res;   
    33073327//              bound = fAct->numCodim2Facets*(this->numVars);
    … …  
    33153335                for(int jj=0;jj<this->numVars;jj++)
    33163336                {
    3317                         (*res)[ii]=(*fNormal)[jj];
     3337                        (*res)[ii]=(int)(*fNormal)[jj];//This is ugly and prone to overflow
    33183338                        ii++;
    33193339                }
    … …  
    33673387                ring inputRing=currRing;        // The ring the user entered
    33683388                ring rootRing;                  // The ring associated to the target ordering
    3369 //              ideal res;
     3389
    33703390                dd_set_global_constants();
    33713391                if(method==noRevS)
    … …  
    33893409                        {
    33903410                                WerrorS("Monomial input - terminating");
    3391                                 exit(-1);
    3392                                 gcAct->getConeNormals(gcAct->gcBasis);
    3393                                 lResList=lprepareResult(gcAct,1);
     3411                                lResList->Init(1);
     3412//                              exit(-1);
     3413//                              gcAct->getConeNormals(gcAct->gcBasis);
     3414//                              lResList=lprepareResult(gcAct,1);
    33943415                                dd_free_global_constants;
    33953416                                //This is filthy
    … …  
    33993420                        gcAct->noRevS(*gcAct);  //Here we go!
    34003421                        //Switch back to the ring the computation was started in
    3401                         rChangeCurrRing(inputRing);
     3422//                      rChangeCurrRing(inputRing);
    34023423                        //res=gcAct->gcBasis;
    34033424                        //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
    … …  
    34203441        else
    34213442        {
     3443                //Simply return an empty list
    34223444                WerrorS("Ring has non-global ordering - terminating");
    34233445                lResList->Init(1);
    3424                 lResList->m[0].rtyp=INT_CMD;
    3425                 int ires=0;
    3426                 lResList->m[0].data=(void*)&ires;
     3446//              lResList->m[0].rtyp=INT_CMD;
     3447//              int ires=0;
     3448//              lResList->m[0].data=(void*)&ires;
    34273449        }
    34283450        //gcone::counter=0;

  • kernel/gfan.h

    r8e77eb r276932  
    209209                inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE);
    210210                inline void readConeFromFile(int gcNum, gcone *gc);
    211                 inline int64vec f2M(gcone *gc, facet *f, int n=1);
     211                inline intvec f2M(gcone *gc, facet *f, int n=1);
    212212                inline void sortRays(gcone *gc);
    213213                //The real stuff
    … …  
    221221                inline void getGB(ideal const &inputIdeal);             
    222222                inline void interiorPoint( dd_MatrixPtr &M, int64vec &iv);
    223                 inline void interiorPoint2(const dd_MatrixPtr &M, int64vec &iv);
     223//              inline void interiorPoint2(const dd_MatrixPtr &M, int64vec &iv);//removed Feb 8th, 2010
    224224                inline void preprocessInequalities(dd_MatrixPtr &M);
    225225                ring rCopyAndAddWeight(const ring &r, int64vec *ivw);
Note: See TracChangeset for help on using the changeset viewer.