Changeset 73809b in git


Ignore:
Timestamp:
Mar 6, 2010, 2:36:13 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b52fc4b2495505785981d640dcf7eb3e456778ef')
Children:
e6ddbdacf047de338a179b12ccb98096d4f02871
Parents:
c57a83edd49526195f04ddf602fca071cafed3b7
Message:
#defined SHALLOW
facet::shallowCopy, facet::shallowDelete. Needed for gcone::enqueue2
Fixes in gcone::areEqual2(facet*,facet*), now the standard for
comparisons of facet
facet::getRef2InteriorPoint
getExtremalRays: Fixed dividing out of gcd
makeInt: now returns primitve vector
renamed gcone::areEqual(intvec&, intvec&) to gcone::ivAreEqual
gcone::normalize no longer needed



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

Legend:

Unmodified
Added
Removed
  • kernel/gfan.cc

    rc57a83e r73809b  
    7575#endif
    7676
     77//NOTE Defining this will slow things down!
     78//Only good for very coarse profiling
    7779#define gfanp
    7880#ifdef gfanp
    7981#include <sys/time.h>
     82#endif
     83
     84#ifndef SHALLOW
     85#define SHALLOW
    8086#endif
    8187
     
    114120* could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
    115121*/
    116 facet::facet(int const &n)
     122facet::facet(const int &n)
    117123{
    118124        this->fNormal=NULL;
     
    144150        //Needed for flip2
    145151        //NOTE ONLY REFERENCE
    146         this->interiorPoint=/*ivCopy(*/f.interiorPoint/*)*/;//only referencing is prolly not the best thing to do in a copy constructor
     152        this->interiorPoint=ivCopy(f.interiorPoint);//only referencing is prolly not the best thing to do in a copy constructor
    147153        facet* f2Copy;
    148154        f2Copy=f.codim2Ptr;
     
    151157        while(f2Copy!=NULL)
    152158        {
    153                 if(f2Act==NULL)
     159                if(f2Act==NULL || f2Act==(facet*)0xfefefefefefefefe)
    154160                {
    155161                        f2Act=new facet(2);
     
    175181
    176182/** \brief Shallow copy constructor for facets
     183* We only need the interior point for equality testing
    177184*/
    178 facet::facet(const facet& f, bool shallow)
    179 {}
     185facet* facet::shallowCopy(const facet& f)
     186{
     187        facet *res = new facet();
     188        res->fNormal=(intvec* const)f.fNormal;
     189        res->UCN=f.UCN;
     190        res->isFlippable=f.isFlippable;
     191        res->interiorPoint=(intvec * const)f.interiorPoint;
     192        res->codim2Ptr=(facet * const)f.codim2Ptr;
     193        res->prev=NULL;
     194        res->next=NULL;
     195        res->flipGB=NULL;
     196        return res;
     197}
     198
     199void facet::shallowDelete()
     200{
     201        this->fNormal=NULL;
     202        this->UCN=0;
     203        this->interiorPoint=NULL;
     204        this->codim2Ptr=NULL;
     205        this->prev=NULL;
     206        this->next=NULL;
     207        this->flipGB=NULL;
     208//      delete this;
     209}
    180210               
    181211/** The default destructor */
     
    217247        return(this->fNormal);
    218248}       
    219        
    220 bool gcone::areEqual2(facet* f, facet *g)
     249
     250/** Equality check for facets based on unique interior points*/
     251inline bool gcone::areEqual2(facet* f, facet *g)
    221252{
    222253#ifdef gfanp
     
    225256        gettimeofday(&start, 0);
    226257#endif
    227         bool res = FALSE;
    228         const intvec* fNormal;
    229         const intvec* gNormal;
    230         fNormal = f->getRef2FacetNormal();
    231         gNormal = g->getRef2FacetNormal();
    232         intvec *fNRef=const_cast<intvec*>(fNormal);
    233         intvec *gNRef=const_cast<intvec*>(gNormal);
    234         if(isParallel(*fNRef,*gNRef))
    235         {
     258        bool res = TRUE;
     259//      const intvec *fNormal;
     260//      const intvec *gNormal;
     261//      fNormal = f->getRef2FacetNormal();
     262//      gNormal = g->getRef2FacetNormal();
     263//      intvec *fNRef=const_cast<intvec*>(fNormal);
     264//      intvec *gNRef=const_cast<intvec*>(gNormal);
     265        /*if(isParallel(fNRef,*gNRef))
     266        {
     267                const intvec *fIntP = f->getRef2InteriorPoint();
     268                const intvec *gIntP = g->getRef2InteriorPoint();
    236269                for(int ii=0;ii<this->numVars;ii++)
    237270                {
    238                         if( (*fNormal)[ii]!=(*gNormal)[ii] )
    239                         {
    240                                 res=TRUE;
     271                        if( (*fIntP)[ii] != (*gIntP)[ii] )
     272                        {
     273                                res=FALSE;
    241274                                break;
    242275                        }
    243276                }
    244277        }
    245         return res;
     278        else
     279                res=FALSE;*/
     280        const intvec *fIntP = f->getRef2InteriorPoint();
     281        const intvec *gIntP = g->getRef2InteriorPoint();
     282        for(int ii=0;ii<this->numVars;ii++)
     283        {
     284                if( (*fIntP)[ii] != (*gIntP)[ii] )
     285                {
     286                        res=FALSE;
     287                        break;
     288                }
     289        }
    246290#ifdef gfanp
    247291        gettimeofday(&end, 0);
    248292        t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    249293#endif 
     294        return res;
    250295}
    251296
     
    385430//      return this->fNormal;
    386431}
    387                
     432
    388433/** Method to print the facet normal*/
    389434inline void facet::printNormal()
     
    456501{
    457502        return ivCopy(this->interiorPoint);
    458 }       
     503}
     504
     505inline const intvec *facet::getRef2InteriorPoint()
     506{
     507        return (this->interiorPoint);
     508}
    459509               
    460510/** \brief Debugging function
     
    565615gcone::~gcone()
    566616{
     617#ifndef NDEBUG
     618  #if SIZEOF_LONG==8
     619        if(this->gcBasis!=(ideal)(0xfbfbfbfbfbfbfbfb))
     620                idDelete((ideal*)&this->gcBasis);
     621  #elif SIZEOF_LONG!=8
     622        if(this->gcBasis!=(ideal)0xfbfbfbfb)
     623                idDelete((ideal *)&this->gcBasis);
     624  #endif
     625#else
    567626        if(this->gcBasis!=NULL)
    568627                idDelete((ideal *)&this->gcBasis);
     628#endif
     629//      idDelete((ideal *)&this->gcBasis);
    569630//      if(this->inputIdeal!=NULL)
    570631//              idDelete((ideal *)&this->inputIdeal);
     
    589650        //dd_FreeMatrix(this->ddFacets);
    590651}                       
    591        
     652
     653/** Returns the number of cones existing at the time*/
    592654inline int gcone::getCounter()
    593655{
     
    603665}
    604666               
    605 /** \brief Return the interior point */
    606 inline intvec *gcone::getIntPoint()
    607 {
    608         return ivCopy(this->ivIntPt);
     667/** \brief Returns either a physical copy the interior point of a cone or just a reference to it.*/
     668inline intvec *gcone::getIntPoint(bool shallow)
     669{
     670        if(shallow==TRUE)
     671                return this->ivIntPt;
     672        else
     673                return ivCopy(this->ivIntPt);
    609674}
    610675               
     
    762827{
    763828        return this->pred;
     829}
     830/** Returns a copy of the this->baseRing */
     831inline ring gcone::getBaseRing()
     832{
     833        return rCopy(this->baseRing);
    764834}
    765835
     
    796866        for (int ii=0;ii<IDELEMS(I);ii++)
    797867        {
    798                 aktpoly=(poly)I->m[ii];
    799                 rows=rows+pLength(aktpoly)-1;
     868//              aktpoly=(poly)I->m[ii];
     869//              rows=rows+pLength(aktpoly)-1;
     870                rows=rows+pLength((poly)I->m[ii])-1;
    800871        }
    801872
     
    10051076        for (int kk = 0; kk<ddrows; kk++)
    10061077        {
     1078                int ggT=1;//NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
    10071079                intvec *load = new intvec(this->numVars);//intvec to store a single facet normal that will then be stored via setFacetNormal
    10081080                for (int jj = 1; jj <ddcols; jj++)
     
    10101082                        double foo;
    10111083                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
    1012                         (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load           
     1084                        (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load
     1085                        ggT = intgcd(ggT,(int&)foo);
    10131086                }//for (int jj = 1; jj <ddcols; jj++)
    1014                                
     1087                if(ggT>1)
     1088                {
     1089                        for(int ll=0;ll<this->numVars;ll++)
     1090                                (*load)[ll] /= ggT;//make primitive vector
     1091                }
    10151092                /*Quick'n'dirty hack for flippability. Executed only if gcone::hasHomInput==FALSE
    10161093                * Otherwise every facet intersects the positive orthant
     
    13121389* So we keep pointer tmp to iv and delete(tmp), so there should not occur a
    13131390* memleak
     1391* TODO normalization
    13141392*/
    13151393void gcone::getExtremalRays(const gcone &gc)
     
    13181396        timeval start, end;
    13191397        gettimeofday(&start, 0);
     1398        timeval poly_start, poly_end;
     1399        gettimeofday(&poly_start,0);
    13201400#endif
    13211401        //Add lineality space - dd_LinealitySpace
     
    13471427        dd_MatrixPtr P;
    13481428        P=dd_CopyGenerators(ddPolyh);
    1349         /*dd_rowset ddlinset,ddredrows; dd_rowindex ddnewpos;
    1350         dd_MatrixCanonicalize(&P, &ddlinset, &ddredrows, &ddnewpos, &err);*/
    13511429        dd_FreePolyhedra(ddPolyh);
     1430        dd_FreeMatrix(ddineq);
     1431#ifdef gfanp
     1432        gettimeofday(&poly_end,0);
     1433        t_ddPolyh += (poly_end.tv_sec - poly_start.tv_sec + 1e-6*(poly_end.tv_usec - poly_start.tv_usec));
     1434#endif
    13521435        /* Compute interior point on the fly*/
    13531436        intvec *ivIntPointOfCone = new intvec(this->numVars);
     1437//      for(int ii=0;ii<P->rowsize;ii++)
     1438//      {}
    13541439        mpq_t *colSum = new mpq_t[this->numVars];
    13551440        int denom[this->numVars];//denominators of colSum
     1441        //NOTE TODO need to gcd of rows and factor out! -> makeInt
    13561442        for(int jj=0;jj<this->numVars;jj++)
    13571443        {
     
    13721458        }
    13731459        //Now compute lcm of denominators of colSum[jj]
     1460        //NOTE gcd as well and normalise instantly?
    13741461        mpz_t kgV; mpz_init(kgV);
    13751462        mpz_set_str(kgV,"1",10);
     
    13891476        mpz_clear(den);
    13901477        mpz_clear(tmp);
     1478        int ggT=(*ivIntPointOfCone)[0];
    13911479        for (int ii=0;ii<(this->numVars);ii++)
    13921480        {
     
    13961484                (*ivIntPointOfCone)[ii]=(int)mpz_get_d(mpq_numref(product));
    13971485                mpq_clear(product);
     1486                //Compute intgcd
     1487                ggT=intgcd(ggT,(*ivIntPointOfCone)[ii]);
     1488        }
     1489        //Divide out a common gcd > 1
     1490        if(ggT>1)
     1491        {
     1492                for(int ii=0;ii<this->numVars;ii++)
     1493                        (*ivIntPointOfCone)[ii] /= ggT;
    13981494        }
    13991495        mpq_clear(qkgV);
     1496        delete [] colSum;
    14001497        /*For homogeneous input (like Det3,3,5) the int points may be negative. So add a suitable multiple of (1,_,1)*/
    14011498        if(hasHomInput==TRUE && iv64isStrictlyPositive(ivIntPointOfCone)==FALSE)
     
    14131510                }
    14141511                delete ivOne;
     1512                int ggT=(*ivIntPointOfCone)[0];
     1513                for(int ii=0;ii<this->numVars;ii++)
     1514                        ggT=intgcd( ggT, (*ivIntPointOfCone)[ii]);
     1515                if(ggT>1)
     1516                {
     1517                        for(int jj=0;jj<this->numVars;jj++)
     1518                                (*ivIntPointOfCone)[jj] /= ggT;
     1519                }
    14151520        }
    14161521        assert(iv64isStrictlyPositive(ivIntPointOfCone));
     1522       
    14171523        this->setIntPoint(ivIntPointOfCone);
    14181524        delete(ivIntPointOfCone);
     
    14301536                {                       
    14311537                        intvec *rowvec = new intvec(this->numVars);
    1432                         makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational                     
     1538                        makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve                 
    14331539                        if(dotProduct(*fNormal,*rowvec)==0)
    14341540                        {
     
    14561562                        }
    14571563                        delete(rowvec);
    1458                 }
     1564                }//For non-homog input ivIntPointOfFacet should already be >0 here
     1565                if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
    14591566                //if we have no strictly pos ray but the input is homogeneous
    14601567                //then add a suitable multiple of (1,...,1)
     
    14771584                                delete tmp;                             
    14781585                        }
    1479                         fAct->setInteriorPoint(ivIntPointOfFacet);
     1586                        //fAct->setInteriorPoint(ivIntPointOfFacet);
    14801587                        delete ivOne;
    14811588                }
    1482                 else
    1483                         fAct->setInteriorPoint(ivIntPointOfFacet);
     1589                //else
     1590                int ggT=(*ivIntPointOfFacet)[0];
     1591                for(int ii=0;ii<this->numVars;ii++)
     1592                        ggT=intgcd(ggT,(*ivIntPointOfFacet)[ii]);
     1593                if(ggT>1)
     1594                {
     1595                        for(int ii=0;ii<this->numVars;ii++)
     1596                                (*ivIntPointOfFacet)[ii] /= ggT;
     1597                }                       
     1598                fAct->setInteriorPoint(ivIntPointOfFacet);
    14841599               
    14851600                delete(ivIntPointOfFacet);
     
    16371752        if( (srcRing->order[0]!=ringorder_a))
    16381753        {
    1639                 intvec *iv = new intvec(this->numVars);
    1640                 iv = ivNeg(fNormal);
     1754                intvec *iv;// = new intvec(this->numVars);
     1755                iv = ivNeg(fNormal);//ivNeg uses ivCopy -> new
    16411756//              tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
    16421757                tmpRing=rCopyAndAddWeight(srcRing,iv);
     
    23452460inline int gcone::dotProduct(const intvec &iva, const intvec &ivb)                             
    23462461{                       
    2347         int res=0;
     2462        int res=0;     
    23482463        for (int i=0;i<this->numVars;i++)
    23492464        {
     
    23602475 */
    23612476inline bool gcone::isParallel(const intvec &a,const intvec &b)
    2362 {                       
     2477{       
     2478/*#ifdef gfanp 
     2479        timeval start, end;
     2480        gettimeofday(&start, 0);
     2481#endif*/               
    23632482        int lhs,rhs;
    23642483        bool res;
     
    23742493                res = FALSE;
    23752494        }
     2495// #ifdef gfanp
     2496//      gettimeofday(&end, 0);
     2497//      t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
     2498// #endif       
    23762499        return res;
    23772500}//bool isParallel
     
    29523075}
    29533076               
     3077//NOTE not needed anywhere
    29543078ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
    29553079{
     
    29733097/** \brief Check for equality of two intvecs
    29743098 */
    2975 inline bool gcone::areEqual(intvec &a, intvec &b)
     3099inline bool gcone::ivAreEqual(intvec &a, intvec &b)
    29763100{
    29773101        bool res=TRUE;
     
    31023226}       
    31033227
    3104 /** Compute the kernel of a given mxn-matrix */
    3105 // dd_MatrixPtr gcone::kernel(const dd_MatrixPtr &M)
    3106 // {
    3107 // }
    31083228
    31093229/** \brief The new method of Markwig and Jensen
     
    31543274                               
    31553275        //Compute unique representation of Facets and rays, i.e. primitive vectors
    3156         gcAct->normalize();
     3276//      gcAct->normalize();
    31573277                       
    31583278        /* Make a copy of the facet list of first cone
     
    31673287                if(fAct->isFlippable==TRUE)
    31683288                {
     3289                        /*Using shallow copy here*/
     3290#ifdef SHALLOW
     3291                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
     3292                        {
     3293                                if(SearchListRoot!=NULL)
     3294                                        delete(SearchListRoot);
     3295                                SearchListRoot = fAct->shallowCopy(*fAct);
     3296                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
     3297                        }
     3298                        else
     3299                        {                       
     3300                                SearchListAct->next = fAct->shallowCopy(*fAct);
     3301                                SearchListAct = SearchListAct->next;                                           
     3302                        }
     3303                        fAct=fAct->next;
     3304#else
     3305                        /*End of shallow copy*/
    31693306                        intvec *fNormal;
    31703307                        fNormal = fAct->getFacetNormal();
     
    32173354                        }
    32183355                        fAct = fAct->next;
    3219                         delete fNormal;                 
     3356                        delete fNormal;
     3357#endif
    32203358                }//if(fAct->isFlippable==TRUE)                 
    32213359                else {fAct = fAct->next;}
     
    32243362        SearchListAct = SearchListRoot; //Set to beginning of list
    32253363        /*Make SearchList doubly linked*/
     3364#ifndef NDEBUG
     3365  #if SIZEOF_LONG==8
     3366        while(SearchListAct!=(facet*)0xfefefefefefefefe && SearchListAct!=NULL)
     3367        {
     3368                if(SearchListAct->next!=(facet*)0xfefefefefefefefe && SearchListAct->next!=NULL){
     3369  #elif SIZEOF_LONG!=8
     3370        while(SearchListAct!=(facet*)0xfefefefe)
     3371        {
     3372                if(SearchListAct->next!=(facet*)0xfefefefe){
     3373  #endif
     3374#else
    32263375        while(SearchListAct!=NULL)
    32273376        {
    3228                 if(SearchListAct->next!=NULL)
    3229                 {
     3377                if(SearchListAct->next!=NULL){
     3378#endif
     3379//              if(SearchListAct->next!=NULL)
     3380//              {
    32303381                        SearchListAct->next->prev = SearchListAct;                                     
    32313382                }
     
    32433394          We always choose the first facet from SearchList as facet to be flipped
    32443395        */                     
    3245         while((SearchListAct!=NULL))// && counter<363)
     3396        while( (SearchListAct!=NULL))//&& counter<490)
    32463397        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
    32473398                fAct = SearchListAct;
     
    32713422//                      gcTmp->getCodim2Normals(*gcTmp);
    32723423                        gcTmp->getExtremalRays(*gcTmp);
    3273                         gcTmp->normalize();// will be done in g2c
     3424//                      gcTmp->normalize();// will be done in g2c
    32743425                        //gcTmp->ddFacets should not be needed anymore, so
    32753426//                      //NOTE If flip2 is used we need to get an interior point of gcTmp
     
    32833434#endif
    32843435                        /*add facets to SLA here*/
     3436#ifdef SHALLOW
     3437                        SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
     3438#else
    32853439                        SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
     3440#endif
    32863441                        if(gfanHeuristic==1)
    32873442                        {
     
    33773532/** \brief Make a set of rational vectors into integers
    33783533 *
    3379  * We compute the lcm of the denominators and multiply with this to get integer values.         
     3534 * We compute the lcm of the denominators and multiply with this to get integer values.
     3535 * If the gcd of the nominators > 1 we divide by the gcd => primitive vector
    33803536 * \param dd_MatrixPtr,intvec
    33813537 */
     
    34233579                       
    34243580//                      mpq_canonicalize(qkgV);
     3581        int ggT=1;
    34253582        for (int ii=0;ii<(M->colsize)-1;ii++)
    34263583        {
    34273584                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
    34283585                n[ii]=(int)mpz_get_d(mpq_numref(res));
     3586                ggT=intgcd(ggT,n[ii]);
     3587        }
     3588        //Normalization
     3589        if(ggT>1)
     3590        {
     3591                for(int ii=0;ii<this->numVars;ii++)
     3592                        n[ii] /= ggT;
    34293593        }
    34303594        delete [] denom;
     
    35413705         */
    35423706        volatile bool removalOccured=FALSE;
    3543         int ctr=0;      //encountered equalities in SLA
    3544         int notParallelCtr=0;
     3707//      int ctr=0;      //encountered equalities in SLA
     3708//      int notParallelCtr=0;
    35453709//      gcone::lengthOfSearchList=1;
    35463710        while(slEnd->next!=NULL)
     
    35553719                {
    35563720                        intvec *fNormal=NULL;
    3557                         fNormal=fAct->getFacetNormal();
     3721                        fNormal=fAct->getFacetNormal();                 
    35583722                        slAct = slHead;
    3559                         notParallelCtr=0;
     3723                        //notParallelCtr=0;
    35603724                        /*If slAct==NULL and fAct!=NULL
    35613725                        we just copy all remaining facets into SLA*/
     
    35703734                                                if(slAct==NULL)
    35713735                                                {
    3572                                                         slAct = new facet(*fCopy);//copy constructor
     3736                                                        slAct = new facet(*fCopy/*,TRUE*/);//copy constructor
    35733737                                                        slHead = slAct;                                                         
    35743738                                                }
    35753739                                                else
    35763740                                                {
    3577                                                         slAct->next = new facet(*fCopy);
     3741                                                        slAct->next = new facet(*fCopy/*,TRUE*/);
    35783742                                                        slAct = slAct->next;
    35793743                                                }
     
    35923756                                cout << "Checking facet (";fNormal->show(1,1);cout << ") against (";slNormal->show(1,1);cout << ")" << endl;
    35933757#endif                         
    3594                                 if(areEqual(fAct,slAct))
    3595                                 {
     3758//                              if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
     3759//                                      exit(-1);
     3760                                if(areEqual2(fAct,slAct))
     3761                                {                                       
    35963762                                        deleteMarker = slAct;
    35973763                                        if(slAct==slHead)
     
    36483814                                facet *f2Act;
    36493815                                f2Act = fAct->codim2Ptr;
    3650                                                
     3816
    36513817                                slEnd->next = new facet();
    36523818                                slEnd = slEnd->next;//Update slEnd
     
    36633829                                //Better: have slEnd->interiorPoint point to fAct->interiorPoint
    36643830                                //NOTE Only reference -> c.f. copy constructor
    3665 //                              slEnd->setInteriorPoint(fAct->getInteriorPoint());
     3831                                //slEnd->setInteriorPoint(fAct->getInteriorPoint());
    36663832                                slEnd->interiorPoint = fAct->interiorPoint;
    36673833                                slEnd->prev = marker;
    36683834                                //Copy codim2-facets
    3669 //                              intvec *f2Normal=new intvec(this->numVars);
     3835                                //intvec *f2Normal=new intvec(this->numVars);
    36703836                                while(f2Act!=NULL)
    36713837                                {
     
    37073873        return slHead;
    37083874}//addC2N
     3875
     3876/** Try using shallow copies*/
     3877facet * gcone::enqueue2(facet *f)
     3878{
     3879        assert(f!=NULL);
     3880#ifdef gfanp
     3881        timeval start, end;
     3882        gettimeofday(&start, 0);
     3883#endif
     3884        facet *slHead;
     3885        slHead = f;
     3886        facet *slAct;   //called with f=SearchListRoot
     3887        slAct = f;
     3888        facet *slEnd;   //Pointer to end of SLA
     3889        slEnd = f;
     3890       
     3891        facet *fAct;
     3892        fAct = this->facetPtr;//New facets to compare
     3893        facet *codim2Act;
     3894        codim2Act = this->facetPtr->codim2Ptr;
     3895//      facet *sl2Act;
     3896//      sl2Act = f->codim2Ptr;
     3897        volatile bool removalOccured=FALSE;
     3898        while(slEnd->next!=NULL)
     3899        {
     3900                slEnd=slEnd->next;
     3901        }       
     3902        while(fAct!=NULL)
     3903        {
     3904                if(fAct->isFlippable)
     3905                {
     3906                        slAct = slHead;
     3907                        if(slAct==NULL)
     3908                        {
     3909                                facet *fCopy;
     3910                                fCopy = fAct;                           
     3911                                while(fCopy!=NULL)
     3912                                {
     3913                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
     3914                                        {
     3915                                                if(slAct==NULL)
     3916                                                {
     3917                                                        slAct = slAct->shallowCopy(*fCopy);//shallow copy constructor
     3918                                                        slHead = slAct;                                                         
     3919                                                }
     3920                                                else
     3921                                                {
     3922                                                        slAct->next = slAct->shallowCopy(*fCopy);
     3923                                                        slAct = slAct->next;
     3924                                                }
     3925                                        }
     3926                                        fCopy = fCopy->next;
     3927                                }
     3928                                break; 
     3929                        }
     3930                        /*Comparison starts here*/
     3931                        while(slAct!=NULL)
     3932                        {
     3933                                removalOccured=FALSE;
     3934#ifdef gfan_DEBUG
     3935cout << "Checking facet (";fAct->fNormal->show(1,1);cout << ") against (";slAct->fNormal->show(1,1);cout << ")" << endl;
     3936#endif 
     3937                                if(areEqual2(fAct,slAct))
     3938                                {
     3939                                        facet *marker=slAct;
     3940                                        if(slAct==slHead)
     3941                                        {                                               
     3942                                                slHead = slAct->next;                                           
     3943                                                if(slHead!=NULL)
     3944                                                        slHead->prev = NULL;
     3945                                        }
     3946                                        else if (slAct==slEnd)
     3947                                        {
     3948                                                slEnd=slEnd->prev;
     3949                                                slEnd->next = NULL;
     3950                                        }                                                               
     3951                                        else
     3952                                        {
     3953                                                slAct->prev->next = slAct->next;
     3954                                                slAct->next->prev = slAct->prev;
     3955                                        }
     3956                                        removalOccured=TRUE;
     3957                                        gcone::lengthOfSearchList--;
     3958#ifdef gfan_DEBUG
     3959                                        cout << "Removing (";fAct->fNormal->show(1,1);cout << ") from list" << endl;
     3960#endif
     3961                                        marker->shallowDelete();
     3962//                                      delete(marker);
     3963                                        break;
     3964                                }
     3965                                slAct = slAct->next;
     3966                        }//while(slAct!=NULL)
     3967                        if(removalOccured==FALSE)
     3968                        {
     3969                                facet *marker=slEnd;
     3970                                slEnd->next = fAct->shallowCopy(*fAct);
     3971                                slEnd = slEnd->next;
     3972                                slEnd->prev=marker;
     3973                                gcone::lengthOfSearchList++;
     3974                        }
     3975                        fAct = fAct->next;
     3976                }
     3977                else
     3978                        fAct = fAct->next;
     3979        }//while(fAct!=NULL)
     3980       
     3981#ifdef gfanp
     3982        gettimeofday(&end, 0);
     3983        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
     3984#endif 
     3985        return slHead;
     3986}
     3987
    37093988/**
    37103989* During flip2 every gc->baseRing gets two ringorder_a
     
    37584037        rChangeCurrRing(srcRing);                       
    37594038}
    3760 //A try with the STL
    3761 facet * gcone::enqueue2(facet *f)
    3762 {
    3763 //      vector<facet*> searchList;
    3764 //      facet *fAct;
    3765 //      fAct=this->facetPtr;
    3766 //      facet *slAct;
    3767 //      slAct = f;
    3768 //      facet *slHead;
    3769 //      while(fAct!=NULL)
    3770 //      {
    3771 //              if(fAct->isFlippable==TRUE)
    3772 //              {
    3773 //                      if(searchList.size()==0)
    3774 //                      {
    3775 //                              facet *fCopy;
    3776 //                              fCopy=fAct;
    3777 //                              while(fCopy!=NULL)
    3778 //                              {
    3779 //                                      if(fCopy->isFlippable==TRUE)
    3780 //                                              searchList.push_back(fCopy);
    3781 //                                      fCopy=fCopy->next;
    3782 //                              }
    3783 //                      }
    3784 //              }
    3785 //      }
    3786 }
    37874039
    37884040/** \brief Compute the gcd of two ints
    37894041 */
    3790 inline int gcone::intgcd(const int a, const int b)
     4042inline int gcone::intgcd(const int &a, const int &b)
    37914043{
    37924044        int r, p=a, q=b;
     
    38694121                                               
    38704122        ofstream gcOutputFile(filename.c_str());
     4123        assert(gcOutputFile);
    38714124        facet *fAct;
    38724125        fAct = gc.facetPtr;                     
     
    41674420}
    41684421       
    4169 ring gcone::getBaseRing()
    4170 {
    4171         return rCopy(this->baseRing);
    4172 }
     4422
    41734423/** \brief Sort the rays of a facet lexicographically
    41744424*/
     
    41914441*
    41924442*/
    4193 lists lprepareResult(gcone *gc, int n)
     4443lists lprepareResult(gcone *gc, const int n)
    41944444{
    41954445        gcone *gcAct;
     
    42054455                {
    42064456                        gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
    4207                         rChangeCurrRing(gcAct->getBaseRing());
     4457                        rChangeCurrRing(gcAct->getBaseRing());//NOTE memleak?
    42084458                }               
    42094459                res->m[ii].rtyp=LIST_CMD;
     
    42244474//              rChangeCurrRing(tmpRing);
    42254475//              l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
    4226                 l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());
     4476                l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());//NOTE memleak?
    42274477//              rChangeCurrRing(saveRing);
    42284478
    42294479                l->m[2].rtyp=INTVEC_CMD;
    4230                 intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
     4480                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));//NOTE memleak?
    42314481                l->m[2].data=(void*)ivCopy(&iv);
    42324482               
     
    43114561float gcone::time_getCodim2Normals;
    43124562float gcone::t_getExtremalRays;
     4563float gcone::t_ddPolyh;
    43134564float gcone::time_flip;
    43144565float gcone::time_flip2;
     
    43224573float gcone::t_mI;
    43234574float gcone::t_iP;
     4575float gcone::t_isParallel;
    43244576unsigned gcone::parallelButNotEqual=0;
    43254577unsigned gcone::numberOfFacetChecks=0;
     
    44214673//      cout << "  t_iP:" << gcone::t_iP << endl;
    44224674        cout << "t_getExtremalRays:" << gcone::t_getExtremalRays << endl;
     4675        cout << "  t_ddPolyh:" << gcone::t_ddPolyh << endl;
    44234676        cout << "t_Flip:" << gcone::time_flip << endl;
    44244677        cout << "  t_markings:" << gcone::t_markings << endl;
     
    44314684        cout << "t_enqueue:" << gcone::time_enqueue << endl;
    44324685        cout << "  t_areEqual:" <<gcone::t_areEqual << endl;
     4686        cout << "t_isParallel:" <<gcone::t_isParallel << endl;
    44334687        cout << endl;
    44344688        cout << "Checked " << gcone::numberOfFacetChecks << " Facets" << endl;
  • kernel/gfan.h

    rc57a83e r73809b  
    2323extern int gfanHeuristic;
    2424// extern dd_MatrixPtr ddLinealitySpace;
    25 #define gfanp
     25
    2626// #ifdef gfanp
    2727// extern       static float time_getConeNormals;
     
    7575                facet *codim2Ptr;       //Pointer to (codim-2)-facet. Bit of recursion here ;-)
    7676                int numCodim2Facets;    //#of (codim-2)-facets of this facet. Set in getCodim2Normals()
     77                unsigned numRays;       //Number of spanning rays of the cone
    7778                ring flipRing;          //the ring on the other side of the facet
    78                 unsigned numRays;       //Number of spanning rays of the cone                   
     79                                       
    7980                /** The default constructor. */
    8081                facet();
    8182                /** Constructor for lower dimensional faces*/
    82                 facet(int const &n);
     83                facet(const int &n);
    8384                /**  The copy constructor */
    8485                facet(const facet& f);
    85                 facet(const facet& f, bool shallow);
     86                facet* shallowCopy(const facet& f);
     87                void shallowDelete();
    8688                /** The default destructor */
    8789                ~facet();
     
    113115                inline void setInteriorPoint(intvec *iv);
    114116                inline intvec *getInteriorPoint();
     117                inline const intvec *getRef2InteriorPoint();
    115118                /** \brief Debugging function
    116119                 * prints the facet normal an all (codim-2)-facets that belong to it
     
    145148                static float time_getCodim2Normals;
    146149                static float t_getExtremalRays;
     150                static float t_ddPolyh;
    147151                static float time_flip;
    148152                static float time_flip2;
     
    157161                static float t_mI;
    158162                static float t_iP;
     163                static float t_isParallel;
    159164                static unsigned parallelButNotEqual;
    160165                static unsigned numberOfFacetChecks;
     
    196201                inline ring getBaseRing();
    197202                inline void setIntPoint(intvec *iv);
    198                 inline intvec *getIntPoint();
     203                inline intvec *getIntPoint(bool shallow=FALSE);
    199204                inline void showIntPoint();
    200205                inline void setNumFacets();
     
    213218//              inline int dotProduct(const intvec* a, const intvec *b);
    214219//              inline bool isParallel(const intvec* a, const intvec* b);
    215                 inline bool areEqual(intvec &a, intvec &b);
     220                inline bool ivAreEqual(intvec &a, intvec &b);
    216221                inline bool areEqual( facet *f,  facet *g);
    217222                inline bool areEqual2(facet* f, facet *g);
    218                 inline int intgcd(int a, int b);
     223                inline int intgcd(const int &a, const int &b);
    219224                inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE);
    220225                inline void readConeFromFile(int gcNum, gcone *gc);
     
    252257                friend class facet;     
    253258};
    254 lists lprepareResult(gcone *gc, int n);
     259lists lprepareResult(gcone *gc, const int n);
    255260// bool iv64isStrictlyPositive(intvec *);
    256261#endif
Note: See TracChangeset for help on using the changeset viewer.