Changeset be7ab3f in git


Ignore:
Timestamp:
Mar 3, 2010, 4:36:29 PM (13 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
9f857dbf446399fd8232a210c134547a6b54d8a2
Parents:
591e530bebccea49fd635f208dd10c42a3b894d7
Message:
Some asserts
Facet copy constructor: no ivCopy of f.interiorPoint, just referencing.
Bad style, should go into some kind of shallow copy
Skeleton for shallow copy
getExtremalRays: For hom input, compute rays, add them up and add a
suitable mulitple of (1,_1); for non hom we add the constraints for the
pos orthant. Both ways we get strictly pos int points.
flip2: Different approach to flipping using an ordering (omega,-N,dp);
requires no marking and is appreciably faster than the old way
interiorPoint2: Legacy code, no longer used since int points are now
computed in getExtremalRays
rCopyAndAddWeight2 (a(),a(),dp)
Earlier computation of linspace and ext rays for 1st cone
interiorPoint goes into the searchlist
replaceDouble_ringorder_a_ByASingleOne()
changes in readConeFromFile - to be undone


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r591e530 rbe7ab3f  
    1313#include "kstd1.h"
    1414#include "kutil.h"      //ksCreateSpoly
    15 // #include "intvec.h"
     15// #include "int64vec.h"
    1616#include "polys.h"
    1717#include "ideals.h"
    … …  
    3131#include <string>
    3232#include <sstream>
    33 #include <time.h>
     33// #include <time.h>
    3434#include <stdlib.h>
    3535//#include <gmpxx.h>
    3636// #include <vector>
     37#include <assert.h>
    3738
    3839
    … …  
    141142        this->UCN=f.UCN;
    142143        this->isFlippable=f.isFlippable;
     144        //Needed for flip2
     145        //NOTE ONLY REFERENCE
     146        this->interiorPoint=/*ivCopy(*/f.interiorPoint/*)*/;//only referencing is prolly not the best thing to do in a copy constructor
    143147        facet* f2Copy;
    144148        f2Copy=f.codim2Ptr;
    … …  
    169173        }       
    170174}
     175
     176/** \brief Shallow copy constructor for facets
     177*/
     178facet::facet(const facet& f, bool shallow)
     179{}
    171180               
    172181/** The default destructor */
    … …  
    207216{
    208217        return(this->fNormal);
    209 }               
     218}       
     219       
     220bool gcone::areEqual2(facet* f, facet *g)
     221{
     222#ifdef gfanp
     223        gcone::numberOfFacetChecks++;
     224        timeval start, end;
     225        gettimeofday(&start, 0);
     226#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        {
     236                for(int ii=0;ii<this->numVars;ii++)
     237                {
     238                        if( (*fNormal)[ii]!=(*gNormal)[ii] )
     239                        {
     240                                res=TRUE;
     241                                break;
     242                        }
     243                }
     244        }
     245        return res;
     246#ifdef gfanp
     247        gettimeofday(&end, 0);
     248        t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
     249#endif 
     250}
     251
    210252/** \brief Comparison of facets
    211253 * called from enqueueNewFacets
    … …  
    419461* prints the facet normal an all (codim-2)-facets that belong to it
    420462*/
    421 inline void facet::fDebugPrint()
     463volatile void facet::fDebugPrint()
    422464{                       
    423465        facet *codim2Act;                       
    … …  
    12791321        //Add lineality space - dd_LinealitySpace
    12801322        dd_MatrixPtr ddineq;
    1281         dd_ErrorType err;
     1323        dd_ErrorType err;       
    12821324        if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
    12831325                ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
    12841326        else
    12851327                ddineq = dd_CopyMatrix(gc.ddFacets);
     1328        /* In case the input is non-homogeneous we add constrains for the positive orthant.
     1329        * This is justified by the fact that for non-homog ideals we only consider the
     1330        * restricted fan. This way we can be sure to find strictly positive interior points.
     1331        * This in turn makes life easy when checking for flippability!
     1332        * Drawback: Makes the LP larger so probably slows down computations a wee bit.
     1333        */
     1334        dd_MatrixPtr ddPosRestr;
     1335        if(hasHomInput==FALSE)
     1336        {
     1337                dd_MatrixPtr tmp;
     1338                ddPosRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
     1339                for(int ii=0;ii<this->numVars;ii++)
     1340                        dd_set_si(ddPosRestr->matrix[ii][ii+1],1);
     1341                dd_MatrixAppendTo(&ddineq,ddPosRestr);
     1342                assert(ddineq);
     1343                dd_FreeMatrix(ddPosRestr);
     1344        }       
    12861345        dd_PolyhedraPtr ddPolyh;
    12871346        ddPolyh = dd_DDMatrix2Poly(ddineq, &err);
    12881347        dd_MatrixPtr P;
    12891348        P=dd_CopyGenerators(ddPolyh);
    1290         dd_FreePolyhedra(ddPolyh);     
     1349        /*dd_rowset ddlinset,ddredrows; dd_rowindex ddnewpos;
     1350        dd_MatrixCanonicalize(&P, &ddlinset, &ddredrows, &ddnewpos, &err);*/
     1351        dd_FreePolyhedra(ddPolyh);
     1352        /* Compute interior point on the fly*/
     1353        intvec *ivIntPointOfCone = new intvec(this->numVars);
     1354        mpq_t *colSum = new mpq_t[this->numVars];
     1355        int denom[this->numVars];//denominators of colSum
     1356        for(int jj=0;jj<this->numVars;jj++)
     1357        {
     1358                mpq_init(colSum[jj]);           
     1359                for(int ii=0;ii<P->rowsize;ii++)
     1360                {
     1361                        mpq_t tmp; mpq_init(tmp);
     1362                        mpq_t sum; mpq_init(sum);
     1363                        mpq_set(sum,colSum[jj]);
     1364                        mpq_add(tmp,sum,P->matrix[ii][jj+1]);
     1365                        mpq_set(colSum[jj],tmp);
     1366                        mpq_clear(tmp);
     1367                }
     1368                mpz_t den; mpz_init(den);
     1369                mpq_get_den(den,colSum[jj]);
     1370                denom[jj]=(int)mpz_get_si(den);
     1371                mpz_clear(den);
     1372        }
     1373        //Now compute lcm of denominators of colSum[jj]
     1374        mpz_t kgV; mpz_init(kgV);
     1375        mpz_set_str(kgV,"1",10);
     1376        mpz_t den; mpz_init(den);
     1377        mpz_t tmp; mpz_init(tmp);
     1378        mpq_get_den(tmp,colSum[0]);
     1379        for (int ii=0;ii<(this->numVars)-1;ii++)
     1380        {
     1381                mpq_get_den(den,colSum[ii+1]);
     1382                mpz_lcm(kgV,tmp,den);
     1383                mpz_set(tmp, kgV);
     1384        }
     1385        mpq_t qkgV;
     1386        mpq_init(qkgV);
     1387        mpq_set_z(qkgV,kgV);
     1388        mpz_clear(kgV);
     1389        mpz_clear(den);
     1390        mpz_clear(tmp);
     1391        for (int ii=0;ii<(this->numVars);ii++)
     1392        {
     1393                mpq_t product;
     1394                mpq_init(product);
     1395                mpq_mul(product,qkgV,colSum[ii]);
     1396                (*ivIntPointOfCone)[ii]=(int)mpz_get_d(mpq_numref(product));
     1397                mpq_clear(product);
     1398        }
     1399        mpq_clear(qkgV);
     1400        /*For homogeneous input (like Det3,3,5) the int points may be negative. So add a suitable multiple of (1,_,1)*/
     1401        if(hasHomInput==TRUE && iv64isStrictlyPositive(ivIntPointOfCone)==FALSE)
     1402        {
     1403                intvec *ivOne = new intvec(this->numVars);
     1404                for(int ii=0;ii<this->numVars;ii++)
     1405                        (*ivOne)[ii]=1;
     1406                while( !iv64isStrictlyPositive(ivIntPointOfCone) )
     1407                {
     1408                        intvec *tmp = ivIntPointOfCone;
     1409                        for(int jj=0;jj<this->numVars;jj++)
     1410                                (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
     1411                        ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
     1412                        delete tmp;                             
     1413                }
     1414                delete ivOne;
     1415        }
     1416        assert(iv64isStrictlyPositive(ivIntPointOfCone));
     1417        this->setIntPoint(ivIntPointOfCone);
     1418        delete(ivIntPointOfCone);
     1419        /* end of interior point computation*/
     1420       
    12911421        //Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
    12921422        int rows=P->rowsize;
    … …  
    13031433                        if(dotProduct(*fNormal,*rowvec)==0)
    13041434                        {
    1305                                 intvec *tmp = ivIntPointOfFacet;
     1435                                intvec *tmp = ivIntPointOfFacet;//Prevent memleak
    13061436                                fAct->numCodim2Facets++;
    13071437                                facet *codim2Act;
    … …  
    13181448                                codim2Act->setFacetNormal(rowvec);
    13191449                                fAct->numRays++;
    1320                                 //TODO Add up to interior point here!
     1450                                //TODO Add up to interior point here!
     1451                                //Memleak sinve ivAdd returns a new intvec
    13211452                                ivIntPointOfFacet=ivAdd(ivIntPointOfFacet,rowvec);
     1453                                //Now tmp still points to the OLD address of ivIntPt
    13221454                                delete(tmp);
    13231455                                       
    … …  
    13251457                        delete(rowvec);
    13261458                }
    1327                 fAct->setInteriorPoint(ivIntPointOfFacet);
     1459                //if we have no strictly pos ray but the input is homogeneous
     1460                //then add a suitable multiple of (1,...,1)
     1461                if( !iv64isStrictlyPositive(ivIntPointOfFacet) && hasHomInput==TRUE)
     1462                {
     1463                        intvec *ivOne = new intvec(this->numVars);
     1464                        for(int ii=0;ii<this->numVars;ii++)
     1465                                (*ivOne)[ii]=1;
     1466//                      intvec *diff = new intvec(this->numVars);
     1467//                      diff=ivSub(ivIntPointOfFacet,ivOne);
     1468                        while( !iv64isStrictlyPositive(ivIntPointOfFacet) )
     1469                        {
     1470                                intvec *tmp = ivIntPointOfFacet;
     1471                                for(int jj=0;jj<this->numVars;jj++)
     1472                                {
     1473                                        (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
     1474//                                      (*diff)[jj]= (*diff)[jj] << 1;
     1475                                }
     1476                                ivIntPointOfFacet = ivAdd(ivIntPointOfFacet/*diff*/,ivOne);
     1477                                delete tmp;                             
     1478                        }
     1479                        fAct->setInteriorPoint(ivIntPointOfFacet);
     1480                        delete ivOne;
     1481                }
     1482                else
     1483                        fAct->setInteriorPoint(ivIntPointOfFacet);
     1484               
    13281485                delete(ivIntPointOfFacet);
    13291486                //Now (if we have at least 3 variables) do a bubblesort on the rays
    … …  
    13681525                fAct = fAct->next;
    13691526        }
     1527        dd_FreeMatrix(P);
    13701528        //Now all extremal rays should be set w.r.t their respective fNormal
    13711529        //TODO Not sufficient -> vol2 II/125&127
    13721530        //NOTE Sufficient according to cddlibs doc. These ARE rays
     1531        //What the hell... let's just take interior points
    13731532        if(gcone::hasHomInput==FALSE)
    13741533        {
    … …  
    13761535                while(fAct!=NULL)
    13771536                {
    1378                         bool containsStrictlyPosRay=FALSE;
    1379                         facet *codim2Act;
    1380                         codim2Act = fAct->codim2Ptr;
    1381                         while(codim2Act!=NULL)
    1382                         {
    1383         //                      containsStrictlyPosRay=TRUE;
    1384                                 intvec *rayvec;
    1385                                 rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
    1386                                 //int negCtr=0;
    1387                                 if(iv64isStrictlyPositive(rayvec))
    1388                                 {
    1389                                         containsStrictlyPosRay=TRUE;
    1390         //                              delete(rayvec);
    1391                                         break;
    1392                                 }
    1393                                 /*for(int ii=0;ii<rayvec->length();ii++)
    1394                                 {
    1395                                         if( (*rayvec)[ii] < 0 )
    1396                                         {
    1397                                                 containsStrictlyPosRay=FALSE;                                   
    1398                                                 break;
    1399                                         }
    1400                                 }
    1401                                 if(containsStrictlyPosRay==TRUE)
    1402                                 {
    1403                                         delete(rayvec);
    1404                                         break;
    1405                                 }*/                     
    1406         //                      delete(rayvec);
    1407                                 codim2Act = codim2Act->next;
    1408                         }
    1409                         if(containsStrictlyPosRay==FALSE)
     1537//                      bool containsStrictlyPosRay=FALSE;
     1538//                      facet *codim2Act;
     1539//                      codim2Act = fAct->codim2Ptr;
     1540//                      while(codim2Act!=NULL)
     1541//                      {
     1542//                              intvec *rayvec;
     1543//                              rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
     1544//                              //int negCtr=0;
     1545//                              if(iv64isStrictlyPositive(rayvec))
     1546//                              {
     1547//                                      containsStrictlyPosRay=TRUE;
     1548//                                      delete(rayvec);
     1549//                                      break;
     1550//                              }                               
     1551//                              delete(rayvec);
     1552//                              codim2Act = codim2Act->next;
     1553//                      }
     1554//                      if(containsStrictlyPosRay==FALSE)
     1555//                              fAct->isFlippable=FALSE;
     1556                        if(!iv64isStrictlyPositive(fAct->interiorPoint))
    14101557                                fAct->isFlippable=FALSE;
    14111558                        fAct = fAct->next;
    … …  
    14181565}               
    14191566       
    1420 inline bool gcone::iv64isStrictlyPositive(intvec * iv64)
     1567inline bool gcone::iv64isStrictlyPositive(const intvec * iv64)
    14211568{
    14221569        bool res=TRUE;
    … …  
    14271574                        res=FALSE;
    14281575                        break;
    1429                 }
     1576                }               
    14301577        }
    14311578        return res;
    … …  
    14531600        intvec *fNormal;// = new intvec(this->numVars); //facet normal, check for parallelity                   
    14541601        fNormal = f->getFacetNormal();  //read this->fNormal;
    1455 
     1602#ifdef gfan_DEBUG
    14561603//      std::cout << "running gcone::flip" << std::endl;
    1457 //      std::cout << "flipping UCN " << this->getUCN() << endl;
    1458 //      cout << "over facet (";
    1459 //      fNormal->show(1,0);
    1460 //      cout << ") with UCN " << f->getUCN();
    1461 //      std::cout << std::endl;
     1604        std::cout << "flipping UCN " << this->getUCN();
     1605        cout << " over facet (";
     1606        fNormal->show(1,0);
     1607        cout << ") with UCN " << f->getUCN();
     1608        std::cout << std::endl;
     1609#endif
    14621610        if(this->getUCN() != f->getUCN())
    14631611        {
    … …  
    15881736                        {
    15891737#ifdef gfan_DEBUG
    1590 //                                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
     1738//                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
    15911739#endif
    15921740                                if (src_ExpV[kk]!=dst_ExpV[kk])
    … …  
    15941742                                        expVAreEqual=FALSE;
    15951743                                }
    1596                         }                                       
    1597                                         //if (*src_ExpV == *dst_ExpV)
     1744                        }                                                               
    15981745                        if (expVAreEqual==TRUE)
    15991746                        {
    … …  
    16491796                        omFree(v);
    16501797                }
    1651 //              omFree(v);
    16521798                omFree(leadExpV);
    16531799        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
    … …  
    17101856        //NOTE Here we should remove interiorPoint and instead
    17111857        // create and ordering like (a(omega),a(fNormal),dp)
    1712         interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
     1858//      if(this->ivIntPt==NULL)
     1859                interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
     1860//      else
     1861//              iv_weight=this->getIntPoint();
    17131862        dd_FreeMatrix(intPointMatrix);
    17141863        /*Crude attempt for interior point */
    … …  
    17341883        gettimeofday(&t_dd_end, 0);
    17351884        t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec));
    1736 #endif                 
     1885#endif 
     1886       
    17371887        /*Step 3
    17381888        * turn the minimal basis into a reduced one */                 
    … …  
    18031953}//void flip(ideal gb, facet *f)
    18041954
     1955/** \brief A slightly different approach to flipping
     1956* Here we use the fact that in_v(in_u(I))=in_(u+eps*v)(I). Therefore, we do no longer
     1957* need to compute an interior point and run BBA on the minimal basis but we can rather
     1958* use the ordering (a(omega),a(fNormal),dp)
     1959* The second parameter facet *f must not be const since we need to store f->flipGB
     1960* Problem: Assume we start in a cone with ordering (dp,C). Then \f$ in_\omega(I) \f$
     1961* will be from a ring with (a(),dp,C) and our resulting cone from (a(),a(),dp,C). Hence a way
     1962* must be found to circumvent the sequence of a()'s growing to a ridiculous size.
     1963* Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we
     1964* do have an interior point of the cone by adding the extremal rays. So we replace
     1965* the latter cone by a cone with (a(sum_of_rays),dp,C).
     1966* Con: It's incredibly ugly
     1967* Pro: No messing around with readConeFromFile()
     1968* Is there a way to construct a vector from \f$ \omega \f$ and the facet normal?
     1969*/
     1970inline void gcone::flip2(const ideal gb, facet *f)
     1971{
     1972#ifdef gfanp
     1973        timeval start, end;
     1974        gettimeofday(&start, 0);
     1975#endif
     1976        /*const*/ intvec *fNormal;
     1977        fNormal = f/*->getRef2FacetNormal();*/->getFacetNormal();       //read this->fNormal;
     1978#ifdef gfan_DEBUG
     1979        std::cout << "flipping UCN " << this->getUCN();
     1980        cout << " over facet (";
     1981        fNormal->show(1,0);
     1982        cout << ") with UCN " << f->getUCN();
     1983        std::cout << std::endl;
     1984#endif
     1985        if(this->getUCN() != f->getUCN())
     1986        {
     1987                WerrorS("Uh oh... Trying to flip over facet with incompatible UCN");
     1988                exit(-1);
     1989        }
     1990        /*1st step: Compute the initial ideal*/
     1991        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);     
     1992        computeInv( gb, initialForm, *fNormal );
     1993        ring srcRing=currRing;
     1994        ring tmpRing;
     1995       
     1996        const intvec *intPointOfFacet;
     1997        intPointOfFacet=f->getInteriorPoint();
     1998        //Now we need two blocks of ringorder_a!
     1999        //May assume the same situation as in flip() here                       
     2000        if( (srcRing->order[0]!=ringorder_a) && (srcRing->order[1]!=ringorder_a) )
     2001        {
     2002                intvec *iv = new intvec(this->numVars);//init with 1s, since we do not need a 2nd block here but later
     2003//              intvec *iv_foo = new intvec(this->numVars,1);//placeholder
     2004                intvec *ivw = ivNeg(fNormal);           
     2005                tmpRing=rCopyAndAddWeight2(srcRing,ivw/*intPointOfFacet*/,iv);
     2006                delete iv;delete ivw;
     2007//              delete iv_foo;
     2008        }
     2009        else
     2010        {
     2011                intvec *iv=new intvec(this->numVars);
     2012                intvec *ivw=ivNeg(fNormal);
     2013                tmpRing=rCopyAndAddWeight2(srcRing,ivw,iv);
     2014                delete iv; delete ivw;
     2015                /*tmpRing=rCopy0(srcRing);
     2016                int length=fNormal->length();
     2017                int *A1=(int *)omAlloc0(length*sizeof(int));
     2018                int *A2=(int *)omAlloc0(length*sizeof(int));
     2019                for(int jj=0;jj<length;jj++)
     2020                {
     2021                        A1[jj] = -(*fNormal)[jj];
     2022                        A2[jj] = 1;//-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
     2023                }
     2024                omFree(tmpRing->wvhdl[0]);
     2025                if(tmpRing->wvhdl[1]!=NULL)
     2026                        omFree(tmpRing->wvhdl[1]);
     2027                tmpRing->wvhdl[0]=(int*)A1;             
     2028                tmpRing->block1[0]=length;
     2029                tmpRing->wvhdl[1]=(int*)A2;             
     2030                tmpRing->block1[1]=length;
     2031                rComplete(tmpRing);*/
     2032        }
     2033        delete fNormal;
     2034        rChangeCurrRing(tmpRing);       
     2035        //Now currRing should have (a(),a(),dp,C)               
     2036        ideal ina;                     
     2037        ina=idrCopyR(initialForm,srcRing);
     2038        idDelete(&initialForm);
     2039        ideal H;
     2040#ifdef gfanp
     2041        timeval t_kStd_start, t_kStd_end;
     2042        gettimeofday(&t_kStd_start,0);
     2043#endif
     2044        BITSET save=test;
     2045        test|=Sy_bit(OPT_REDSB);
     2046        test|=Sy_bit(OPT_REDTAIL);
     2047//      if(gcone::hasHomInput==TRUE)
     2048                H=kStd(ina,NULL,testHomog/*isHomog*/,NULL/*,gcone::hilbertFunction*/);
     2049//      else
     2050//              H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
     2051        test=save;
     2052#ifdef gfanp
     2053        gettimeofday(&t_kStd_end, 0);
     2054        t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
     2055#endif
     2056        idSkipZeroes(H);
     2057        idDelete(&ina);
     2058       
     2059        rChangeCurrRing(srcRing);
     2060        ideal srcRing_H;
     2061        ideal srcRing_HH;                       
     2062        srcRing_H=idrCopyR(H,tmpRing);
     2063        //H is needed further below, so don't idDelete here
     2064        srcRing_HH=ffG(srcRing_H,this->gcBasis);
     2065        idDelete(&srcRing_H);
     2066        //Now BBA(srcRing_HH) with (a(),a(),dp)
     2067        /* Evil modification of currRing */
     2068        ring dstRing=rCopy0(tmpRing);
     2069        int length=this->numVars;
     2070        int *A1=(int *)omAlloc0(length*sizeof(int));
     2071        int *A2=(int *)omAlloc0(length*sizeof(int));
     2072        const intvec *ivw=f->getRef2FacetNormal();
     2073        for(int jj=0;jj<length;jj++)
     2074        {
     2075                A1[jj] = (*intPointOfFacet)[jj];
     2076                A2[jj] = -(*ivw)[jj];//TODO Or minus (*ivw)[ii] ??? NOTE minus
     2077        }
     2078        omFree(dstRing->wvhdl[0]);
     2079        if(dstRing->wvhdl[1]!=NULL)
     2080                omFree(dstRing->wvhdl[1]);
     2081        dstRing->wvhdl[0]=(int*)A1;             
     2082        dstRing->block1[0]=length;
     2083        dstRing->wvhdl[1]=(int*)A2;             
     2084        dstRing->block1[1]=length;
     2085        rComplete(dstRing);
     2086        rChangeCurrRing(dstRing);
     2087        ideal dstRing_I;                       
     2088        dstRing_I=idrCopyR(srcRing_HH,srcRing);
     2089        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more       
     2090        save=test;
     2091        test|=Sy_bit(OPT_REDSB);
     2092        test|=Sy_bit(OPT_REDTAIL);
     2093        ideal tmpI;
     2094        tmpI = dstRing_I;
     2095#ifdef gfanp
     2096//      timeval t_kStd_start, t_kStd_end;
     2097        gettimeofday(&t_kStd_start,0);
     2098#endif
     2099//      if(gcone::hasHomInput==TRUE)
     2100//              dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
     2101//      else
     2102                dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
     2103#ifdef gfanp
     2104        gettimeofday(&t_kStd_end, 0);
     2105        t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
     2106#endif
     2107        idDelete(&tmpI);
     2108        idNorm(dstRing_I);                     
     2109        idSkipZeroes(dstRing_I);
     2110        test=save;
     2111        /*End of step 3 - reduction*/
     2112       
     2113        f->setFlipGB(dstRing_I);
     2114        f->flipRing=rCopy(dstRing);
     2115        rDelete(tmpRing);
     2116        rDelete(dstRing);
     2117        //Now we should have dstRing with (a(),a(),dp,C)
     2118        //This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
     2119        //of cones in noRevS
     2120        rChangeCurrRing(srcRing);
     2121#ifdef gfanp
     2122        gettimeofday(&end, 0);
     2123        time_flip2 += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
     2124#endif
     2125}//flip2
     2126
    18052127/** \brief Compute initial ideal
    18062128 * Compute the initial ideal in_v(G) wrt a (possible) facet normal
    … …  
    18082130 * and in gcone::flip for obvious reasons.
    18092131*/
    1810 inline void gcone::computeInv(ideal &gb, ideal &initialForm, intvec &fNormal)
     2132inline void gcone::computeInv(const ideal &gb, ideal &initialForm, const intvec &fNormal)
    18112133{
    18122134#ifdef gfanp
    … …  
    20002322/** \brief Compute the negative of a given intvec
    20012323        */             
    2002 inline intvec *gcone::ivNeg(intvec *iv)
     2324inline intvec *gcone::ivNeg(const intvec *iv)
    20032325{       //Hm, switching to intvec const intvec does no longer work
    20042326        intvec *res;// = new intvec(iv->length());
    … …  
    20122334*
    20132335*/
    2014 inline int gcone::dotProduct(intvec &iva, intvec &ivb)                         
    2015 {                       
    2016         int res=0;
    2017         for (int i=0;i<this->numVars;i++)
    2018         {
    2019                 res = res+(iva[i]*ivb[i]);
    2020         }
    2021         return res;
    2022 }//int dotProduct
     2336// inline int gcone::dotProduct(intvec &iva, intvec &ivb)                               
     2337// {                   
     2338//      int res=0;
     2339//      for (int i=0;i<this->numVars;i++)
     2340//      {
     2341//              res = res+(iva[i]*ivb[i]);
     2342//      }
     2343//      return res;
     2344// }//int dotProduct
    20232345inline int gcone::dotProduct(const intvec &iva, const intvec &ivb)                             
    20242346{                       
    … …  
    20262348        for (int i=0;i<this->numVars;i++)
    20272349        {
     2350// #ifndef NDEBUG
     2351//      (const_cast<intvec*>(&iva))->show(1,0); (const_cast<intvec*>(&ivb))->show(1,0);
     2352// #endif
    20282353                res = res+(iva[i]*ivb[i]);
    20292354        }
    … …  
    20342359 * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
    20352360 */
    2036 inline bool gcone::isParallel(intvec &a, intvec &b)
     2361inline bool gcone::isParallel(const intvec &a,const intvec &b)
    20372362{                       
    20382363        int lhs,rhs;
    … …  
    21802505}//void interiorPoint(dd_MatrixPtr const &M)
    21812506
    2182 
     2507/** Computes an interior point of a cone by taking two interior points a,b from two different facets
     2508* and then computing b+(a-b)/2
     2509* Of course this only works for flippable facets
     2510* Two cases may occur:
     2511* 1st: There are only two facets who share the only strictly positive ray
     2512* 2nd: There are at least two facets which have a distinct positive ray
     2513* In the former case we use linear algebra to determine an interior point,
     2514* in the latter case we simply add up the two rays
     2515*
     2516* Way too bad! The case may occur that the cone is spanned by three rays, of which only two are strictly
     2517* positive => these lie in a plane and thus their sum is not from relative interior.
     2518* So let's just sum up all rays, find one strictly positive and shift the point along that ray
     2519*
     2520* Used by noRevS
     2521*/
     2522inline void gcone::interiorPoint2()
     2523{//idPrint(this->gcBasis);
     2524#ifdef gfan_DEBUG
     2525        if(this->ivIntPt!=NULL)
     2526                WarnS("Interior point already exists - ovrewriting!");
     2527#endif
     2528        facet *f1 = this->facetPtr;
     2529        facet *f2 = NULL;
     2530        intvec *intF1=NULL;
     2531        while(f1!=NULL)
     2532        {
     2533                if(f1->isFlippable)
     2534                {
     2535                        facet *f1Ray = f1->codim2Ptr;
     2536                        while(f1Ray!=NULL)
     2537                        {
     2538                                const intvec *check=f1Ray->getRef2FacetNormal();
     2539                                if(iv64isStrictlyPositive(check))
     2540                                {
     2541                                        intF1=ivCopy(check);
     2542                                        break;
     2543                                }                               
     2544                                f1Ray=f1Ray->next;
     2545                        }
     2546                }
     2547                if(intF1!=NULL)
     2548                        break;
     2549                f1=f1->next;
     2550        }
     2551        if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
     2552                f2=f1->next;
     2553        else
     2554                f2=this->facetPtr;
     2555        if(intF1==NULL && hasHomInput==TRUE)
     2556        {
     2557                intF1 = new intvec(this->numVars);
     2558                for(int ii=0;ii<this->numVars;ii++)
     2559                        (*intF1)[ii]=1;
     2560        }
     2561        assert(f1); assert(f2);
     2562        intvec *intF2=f2->getInteriorPoint();
     2563        mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
     2564        mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
     2565        mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
     2566        for(int ii=0;ii<this->numVars;ii++)
     2567        {
     2568                mpq_init(qPosRay[ii]);
     2569                mpq_init(qIntPt[ii]);
     2570                mpq_init(qPosIntPt[ii]);
     2571        }       
     2572        //Compute a+((b-a)/2) && Convert intF1 to mpq
     2573        for(int ii=0;ii<this->numVars;ii++)
     2574        {
     2575                mpq_t a,b;
     2576                mpq_init(a); mpq_init(b);
     2577                mpq_set_si(a,(*intF1)[ii],1);
     2578                mpq_set_si(b,(*intF2)[ii],1);
     2579                mpq_t diff;
     2580                mpq_init(diff);
     2581                mpq_sub(diff,b,a);      //diff=b-a
     2582                mpq_t quot;
     2583                mpq_init(quot);
     2584                mpq_div_2exp(quot,diff,1);      //quot=diff/2=(b-a)/2
     2585                mpq_clear(diff);
     2586                //Don't be clever and reuse diff here
     2587                mpq_t sum; mpq_init(sum);
     2588                mpq_add(sum,b,quot);    //sum=b+quot=a+(b-a)/2
     2589                mpq_set(qIntPt[ii],sum);
     2590                mpq_clear(sum);
     2591                mpq_clear(quot);
     2592                mpq_clear(a); mpq_clear(b);
     2593                //Now for intF1
     2594                mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
     2595        }
     2596        //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
     2597        while(TRUE)
     2598        {       
     2599                bool success=FALSE;
     2600                int posCtr=0;   
     2601                for(int ii=0;ii<this->numVars;ii++)
     2602                {
     2603                        mpq_t sum; mpq_init(sum);
     2604                        mpq_add(sum,qPosRay[ii],qIntPt[ii]);
     2605                        mpq_set(qPosIntPt[ii],sum);
     2606                        mpq_clear(sum);
     2607                        if(mpq_sgn(qPosIntPt[ii])==1)
     2608                                posCtr++;
     2609                }
     2610                if(posCtr==this->numVars)//qPosIntPt > 0
     2611                        break;
     2612                else
     2613                {
     2614                        mpq_t qTwo; mpq_init(qTwo);
     2615                        mpq_set_ui(qTwo,2,1);
     2616                        for(int jj=0;jj<this->numVars;jj++)
     2617                        {
     2618                                mpq_t tmp; mpq_init(tmp);
     2619                                mpq_mul(tmp,qPosRay[jj],qTwo);
     2620                                mpq_set( qPosRay[jj], tmp);
     2621                                mpq_clear(tmp);
     2622                        }
     2623                        mpq_clear(qTwo);
     2624                }
     2625        }//while
     2626        //Now qPosIntPt ought to be >0, so convert back to int :D
     2627        /*Compute lcm of the denominators*/
     2628        mpz_t *denom = new mpz_t[this->numVars];
     2629        mpz_t tmp,kgV;
     2630        mpz_init(tmp); mpz_init(kgV);
     2631        for (int ii=0;ii<this->numVars;ii++)
     2632        {
     2633                mpz_t z;
     2634                mpz_init(z);
     2635                mpq_get_den(z,qPosIntPt[ii]);
     2636                mpz_init(denom[ii]);
     2637                mpz_set( denom[ii], z);
     2638                mpz_clear(z);                           
     2639        }
     2640               
     2641        mpz_set(tmp,denom[0]);
     2642        for (int ii=0;ii<this->numVars;ii++)
     2643        {
     2644                mpz_lcm(kgV,tmp,denom[ii]);
     2645                mpz_set(tmp,kgV);                               
     2646        }
     2647        mpz_clear(tmp);
     2648        /*Multiply the nominators by kgV*/
     2649        mpq_t qkgV,res;
     2650        mpq_init(qkgV);
     2651        mpq_canonicalize(qkgV);         
     2652        mpq_init(res);
     2653        mpq_canonicalize(res);
     2654                               
     2655        mpq_set_num(qkgV,kgV);
     2656        intvec *n=new intvec(this->numVars);
     2657        for (int ii=0;ii<this->numVars;ii++)
     2658        {
     2659                mpq_canonicalize(qPosIntPt[ii]);
     2660                mpq_mul(res,qkgV,qPosIntPt[ii]);
     2661                (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
     2662        }
     2663        this->setIntPoint(n);
     2664        delete n;
     2665        delete [] qPosIntPt;
     2666        delete [] denom;
     2667        delete [] qPosRay;
     2668        delete [] qIntPt;
     2669        mpz_clear(kgV);
     2670        mpq_clear(qkgV); mpq_clear(res);
     2671       
     2672/***************************************************/
     2673//      //We need to find out if there are at least two pos rays
     2674//      facet *f1 = this->facetPtr;
     2675//      facet *f2 = NULL;
     2676//      intvec *intF1=NULL;
     2677//      intvec *intF2=NULL;
     2678//      while(f1!=NULL)
     2679//      {
     2680//              if(f1->isFlippable)
     2681//              {
     2682//                      facet *f1Ray=f1->codim2Ptr;
     2683//                      while(f1Ray!=NULL)
     2684//                      {
     2685//                              intvec *check = f1Ray->getFacetNormal();
     2686//                              if(iv64isStrictlyPositive(check))
     2687//                              {
     2688//                                      intF1=ivCopy(check);
     2689//                                      delete check;
     2690//                                      break;
     2691//                              }
     2692//                              delete check;
     2693//                              f1Ray = f1Ray->next;
     2694//                      }       
     2695//              }
     2696//              if(intF1!=NULL)//We have found the first strictly positive ray
     2697//                      break;
     2698//              f1 = f1->next;
     2699//      }
     2700//      if(f1->next!=NULL)
     2701//              f2=f1->next;
     2702//      while(f2!=NULL)
     2703//      {
     2704//              if(f2->isFlippable)
     2705//              {
     2706//                      facet *f2Ray=f2->codim2Ptr;
     2707//                      while(f2Ray!=NULL)
     2708//                      {
     2709//                              intvec *check = f2Ray->getFacetNormal();
     2710//                              if(iv64isStrictlyPositive(check))
     2711//                              {
     2712//                                      //and not equal to intF1
     2713//                                      if(check->compare(intF1)!=0)//we found a distinct 2nd ray
     2714//                                      {
     2715//                                              intF2=ivCopy(check);
     2716//                                              delete check;
     2717//                                      }
     2718//                                      break;
     2719//                              }
     2720//                              delete check;
     2721//                              f2Ray = f2Ray->next;
     2722//                      }
     2723//              }
     2724//              if(intF2!=NULL)
     2725//                      break;
     2726//              f2=f2->next;
     2727//      }
     2728//      if(intF2==NULL)//do some linear algebra
     2729//      {
     2730//              f2=f1->next;
     2731//              intvec *arbitraryRay = f2->codim2Ptr->getFacetNormal();
     2732//              //Now add
     2733//              intvec *ivSum = ivAdd(intF1, arbitraryRay);
     2734//              while(iv64isStrictlyPositive(ivSum)==FALSE)
     2735//              {
     2736//                      intvec *tmp = intF1;
     2737//                      for(int ii=0;ii<this->numVars;ii++)
     2738//                              (*intF1)[ii] = (*intF1)[ii] << 1;       //times 2
     2739//                      intF1 = ivAdd(intF1,ivSum);
     2740//                      delete tmp;
     2741//                     
     2742//              }
     2743//              this->setIntPoint(ivSum);
     2744//
     2745//      }
     2746//      else    //just add intF1+intF2
     2747//      {
     2748//              intvec *sum = ivAdd(intF1, intF2);
     2749//              this->setIntPoint(sum);
     2750//              delete sum;
     2751//      }
     2752/***************************************************/
     2753//      //Find 1st flippable facet
     2754//      facet *f1 = this->facetPtr;
     2755//      facet *f2 = NULL; //idPrint(this->gcBasis);
     2756//      intvec *intF1=NULL;
     2757//      intvec *intF2=NULL;
     2758//      while(f1->isFlippable==FALSE && f1!=NULL)
     2759//              f1=f1->next;
     2760//      if(f1!=NULL)
     2761//      {
     2762//              intF1 = f1->getInteriorPoint();
     2763//              f2 = f1->next;
     2764//      }
     2765//      while(f2!=NULL && f2->isFlippable==FALSE)
     2766//              f2=f2->next;
     2767//      if(f1==NULL || f2==NULL) //Is there only one flippable facet?
     2768//      {       //f1==NULL would mean there ain't any flippable facet...
     2769//              /*If f2==NULL we ignore f2 and try to find an int point
     2770//              * using f1 and linear algebra
     2771//              */
     2772//              bool foundIntPoint=FALSE;
     2773//              intvec *f1Ray=f1->codim2Ptr->getFacetNormal();
     2774//              const intvec *fNormal=f1->getRef2FacetNormal();
     2775//              while(foundIntPoint==FALSE)
     2776//              {                       
     2777//                      intvec *res;// = new intvec(this->numVars);
     2778//                      res = ivAdd(const_cast<intvec*>(fNormal), f1Ray);
     2779//                      int posCtr=0;   //Count the pos entries of res
     2780//                      for(int ii=0;ii<this->numVars;ii++)
     2781//                      {
     2782//                              if( (*res)[ii]>0)
     2783//                                      posCtr++;                       
     2784//                      }                       
     2785//                      if(posCtr==this->numVars)
     2786//                      {
     2787//                              foundIntPoint=TRUE;
     2788//                              this->setIntPoint(res);
     2789//                      }
     2790//                      else
     2791//                      {
     2792//                              for(int ii=0;ii<this->numVars;ii++)
     2793//                                      (*f1Ray)[ii] = ((*f1Ray)[ii]) << 1;
     2794//                      }
     2795//                      delete res;
     2796//              }
     2797//              //Now we should have an int point
     2798//              delete f1Ray;
     2799//      }
     2800//      else    //ok, we have two flippable facets
     2801//      {               
     2802//              intF2=f2->getInteriorPoint();
     2803//              mpq_t *ratIntPt = new mpq_t[this->numVars];
     2804//              for(int ii=0;ii<this->numVars;ii++)
     2805//                      mpq_init(ratIntPt[ii]);
     2806//              for(int ii=0;ii<this->numVars;ii++)
     2807//              {
     2808//                      mpq_t a,b;
     2809//                      mpq_init(a); mpq_init(b);
     2810//                      mpq_set_si(a,(*intF1)[ii],1);
     2811//                      mpq_set_si(b,(*intF2)[ii],1);
     2812//                      mpq_t diff;
     2813//                      mpq_init(diff);
     2814//                      mpq_sub(diff,a,b);      //diff=a-b
     2815//                      mpq_t quot;
     2816//                      mpq_init(quot);
     2817//                      mpq_div_2exp(quot,diff,1);      //quot=diff/2=(a-b)/2
     2818//                      mpq_clear(diff);
     2819//                      //Don't be clever and reuse diff here
     2820//                      mpq_t sum; mpq_init(sum);
     2821//                      mpq_add(sum,b,quot);    //sum=b+quot=b+(a-b)/2
     2822//                      mpq_set(ratIntPt[ii],sum);
     2823//                      mpq_clear(sum);
     2824//                      mpq_clear(quot);
     2825//                      mpq_clear(a); mpq_clear(b);
     2826//              }       
     2827//              /*Compute lcm of the denominators*/
     2828//              mpz_t *denom = new mpz_t[this->numVars];
     2829//              mpz_t tmp,kgV;
     2830//              mpz_init(tmp); mpz_init(kgV);
     2831//              for (int ii=0;ii<this->numVars;ii++)
     2832//              {
     2833//                      mpz_t z;
     2834//                      mpz_init(z);
     2835//                      mpq_get_den(z,ratIntPt[ii]);
     2836//                      mpz_init(denom[ii]);
     2837//                      mpz_set( denom[ii], z);
     2838//                      mpz_clear(z);                           
     2839//              }
     2840//             
     2841//              mpz_set(tmp,denom[0]);
     2842//              for (int ii=0;ii<this->numVars;ii++)
     2843//              {
     2844//                      mpz_lcm(kgV,tmp,denom[ii]);
     2845//                      mpz_set(tmp,kgV);                               
     2846//              }
     2847//              mpz_clear(tmp);
     2848//              /*Multiply the nominators by kgV*/
     2849//              mpq_t qkgV,res;
     2850//              mpq_init(qkgV);
     2851//              /*mpq_set_str(qkgV,"1",10);*/
     2852//              mpq_canonicalize(qkgV);
     2853//                             
     2854//              mpq_init(res);
     2855//              /*mpq_set_str(res,"1",10);*/
     2856//              mpq_canonicalize(res);
     2857//                             
     2858//              mpq_set_num(qkgV,kgV);
     2859//              intvec *n=new intvec(this->numVars);
     2860//              for (int ii=0;ii<this->numVars;ii++)
     2861//              {
     2862//                      mpq_canonicalize(ratIntPt[ii]);mpq_mul(res,qkgV,ratIntPt[ii]);
     2863//                      (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
     2864//              }
     2865//              this->setIntPoint(n);
     2866//              delete n;
     2867//              delete [] ratIntPt;
     2868//              delete [] denom;
     2869//              mpz_clear(kgV);
     2870//              mpq_clear(qkgV); mpq_clear(res);
     2871//      }
     2872}
    21832873       
    21842874/** \brief Copy a ring and add a weighted ordering in first place
    … …  
    22192909        return res;
    22202910}//rCopyAndAdd
     2911               
     2912ring gcone::rCopyAndAddWeight2(const ring &r,const intvec *ivw, const intvec *fNormal)                         
     2913{
     2914        ring res=rCopy0(r);
     2915                       
     2916        omFree(res->order);
     2917        res->order =(int *)omAlloc0(5*sizeof(int));
     2918        omFree(res->block0);
     2919        res->block0=(int *)omAlloc0(5*sizeof(int));
     2920        omFree(res->block1);
     2921        res->block1=(int *)omAlloc0(5*sizeof(int));
     2922        omfree(res->wvhdl);
     2923        res->wvhdl =(int **)omAlloc0(5*sizeof(int**));
     2924                       
     2925        res->order[0]=ringorder_a;
     2926        res->block0[0]=1;
     2927        res->block1[0]=res->N;
     2928        res->order[1]=ringorder_a;
     2929        res->block0[1]=1;
     2930        res->block1[1]=res->N;
     2931       
     2932        res->order[2]=ringorder_dp;
     2933        res->block0[2]=1;
     2934        res->block1[2]=res->N;
     2935       
     2936        res->order[3]=ringorder_C;
     2937
     2938        int length=ivw->length();
     2939        int *A1=(int *)omAlloc0(length*sizeof(int));
     2940        int *A2=(int *)omAlloc0(length*sizeof(int));
     2941        for (int jj=0;jj<length;jj++)
     2942        {                               
     2943                A1[jj]=(*ivw)[jj];
     2944                A2[jj]=-(*fNormal)[jj];
     2945        }                       
     2946        res->wvhdl[0]=(int *)A1;
     2947        res->block1[0]=length;
     2948        res->wvhdl[1]=(int *)A2;
     2949        res->block1[1]=length;
     2950        rComplete(res);
     2951        return res;
     2952}
    22212953               
    22222954ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
    … …  
    24143146        /*Compute lineality space here
    24153147        * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/
    2416         dd_LinealitySpace = gcAct->computeLinealitySpace();
     3148        if(dd_LinealitySpace==NULL)
     3149                dd_LinealitySpace = gcAct->computeLinealitySpace();
    24173150        /*End of lineality space computation*/         
    24183151//      gcAct->getCodim2Normals(*gcAct);
    2419         gcAct->getExtremalRays(*gcAct);
     3152        if(fAct->codim2Ptr==NULL)
     3153                gcAct->getExtremalRays(*gcAct);
    24203154                               
    24213155        //Compute unique representation of Facets and rays, i.e. primitive vectors
    … …  
    24483182                        SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
    24493183                        SearchListAct->isFlippable=TRUE;
     3184                        //Copy int point as well
     3185                        intvec *ivIntPt;// = new intvec(this->numVars);
     3186                        ivIntPt = fAct->getInteriorPoint();
     3187                        SearchListAct->setInteriorPoint(ivIntPt);
     3188                        delete(ivIntPt);
    24503189                        //Copy codim2-facets
    24513190                        facet *codim2Act;
    … …  
    25043243          We always choose the first facet from SearchList as facet to be flipped
    25053244        */                     
    2506         while((SearchListAct!=NULL))// && gcone::counter<10)
     3245        while((SearchListAct!=NULL))// && counter<363)
    25073246        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
    25083247                fAct = SearchListAct;
    … …  
    25113250//              while( (fAct->getUCN() == fAct->next->getUCN()) )               
    25123251                {       //Since SLA should only contain flippables there should be no need to check for that                   
    2513                         gcAct->flip(gcAct->gcBasis,fAct);                       
     3252                        gcAct->flip2(gcAct->gcBasis,fAct);
    25143253                        //NOTE rCopy needed?
    25153254                        ring rTmp=rCopy(fAct->flipRing);
    … …  
    25343273                        gcTmp->normalize();// will be done in g2c
    25353274                        //gcTmp->ddFacets should not be needed anymore, so
     3275//                      //NOTE If flip2 is used we need to get an interior point of gcTmp
     3276//                      // and replace gcTmp->baseRing with an appropriate ring with only
     3277//                      // one weight
     3278//                      gcTmp->interiorPoint2();
     3279                        gcTmp->replaceDouble_ringorder_a_ByASingleOne();
    25363280                        dd_FreeMatrix(gcTmp->ddFacets);
    25373281#ifdef gfan_DEBUG
    … …  
    28263570                                                if(slAct==NULL)
    28273571                                                {
    2828                                                         slAct = new facet(*fCopy);
     3572                                                        slAct = new facet(*fCopy);//copy constructor
    28293573                                                        slHead = slAct;                                                         
    28303574                                                }
    … …  
    29153659                                slEnd->isFlippable = TRUE;
    29163660                                slEnd->setFacetNormal(fNormal);
     3661                                //NOTE Add interior point here
     3662                                //This is ugly but needed for flip2
     3663                                //Better: have slEnd->interiorPoint point to fAct->interiorPoint
     3664                                //NOTE Only reference -> c.f. copy constructor
     3665//                              slEnd->setInteriorPoint(fAct->getInteriorPoint());
     3666                                slEnd->interiorPoint = fAct->interiorPoint;
    29173667                                slEnd->prev = marker;
    29183668                                //Copy codim2-facets
    … …  
    29573707        return slHead;
    29583708}//addC2N
    2959 
     3709/**
     3710* During flip2 every gc->baseRing gets two ringorder_a
     3711* To avoid having an awful lot of those in the end we endow
     3712* gc->baseRing by a suitable ring with (a,dp,C) and copy all
     3713* necessary stuff over
     3714* But first we will try to just do an inplace exchange and copying only the
     3715* gc->gcBasis
     3716*/
     3717inline void gcone::replaceDouble_ringorder_a_ByASingleOne()
     3718{
     3719        ring srcRing=currRing;
     3720        ring replacementRing=rCopy0((ring)this->baseRing);
     3721        /*We assume we have (a(),a(),dp) here*/
     3722        omFree(replacementRing->order);
     3723        replacementRing->order =(int *)omAlloc0(4*sizeof(int));
     3724        omFree(replacementRing->block0);
     3725        replacementRing->block0=(int *)omAlloc0(4*sizeof(int));
     3726        omFree(replacementRing->block1);
     3727        replacementRing->block1=(int *)omAlloc0(4*sizeof(int));
     3728        omfree(replacementRing->wvhdl);
     3729        replacementRing->wvhdl =(int **)omAlloc0(4*sizeof(int**));
     3730                       
     3731        replacementRing->order[0]=ringorder_a;
     3732        replacementRing->block0[0]=1;
     3733        replacementRing->block1[0]=replacementRing->N;
     3734               
     3735        replacementRing->order[1]=ringorder_dp;
     3736        replacementRing->block0[1]=1;
     3737        replacementRing->block1[1]=replacementRing->N;
     3738       
     3739        replacementRing->order[2]=ringorder_C;
     3740
     3741        intvec *ivw = this->getIntPoint();
     3742//      assert(this->ivIntPt);
     3743        int length=ivw->length();       
     3744        int *A=(int *)omAlloc0(length*sizeof(int));
     3745        for (int jj=0;jj<length;jj++)
     3746                A[jj]=(*ivw)[jj];
     3747        delete ivw;
     3748        replacementRing->wvhdl[0]=(int *)A;
     3749        replacementRing->block1[0]=length;
     3750        /*Finish*/
     3751        rComplete(replacementRing);
     3752        rChangeCurrRing(replacementRing);
     3753        ideal temporaryGroebnerBasis=idrCopyR(this->gcBasis,this->baseRing);
     3754        rDelete(this->baseRing);
     3755        this->baseRing=rCopy(replacementRing);
     3756        this->gcBasis=idCopy(temporaryGroebnerBasis);
     3757        /*And back to where we came from*/
     3758        rChangeCurrRing(srcRing);                       
     3759}
    29603760//A try with the STL
    29613761facet * gcone::enqueue2(facet *f)
    … …  
    30743874        f2Act=fAct->codim2Ptr;
    30753875       
    3076         char *ringString = rString(currRing);
     3876        char *ringString = rString(gc.baseRing);
    30773877                       
    30783878        if (!gcOutputFile)
    … …  
    31413941               
    31423942/** \brief Reads a cone from a file identified by its number
    3143  */
     3943* ||depending on whether flip or flip2 is used, switch the flag flipFlag
     3944* ||defaults to 0 => flip
     3945* ||1 => flip2
     3946*/
    31443947inline void gcone::readConeFromFile(int UCN, gcone *gc)
    31453948{
    … …  
    31753978                        string strweight;
    31763979                        strweight=line.substr(0,line.find_first_of(")"));
    3177                         intvec *iv=new intvec(this->numVars);
     3980                                               
     3981                        intvec *iv=new intvec(this->numVars);//                         
    31783982                        for(int ii=0;ii<this->numVars;ii++)
    31793983                        {
    … …  
    31833987                                line.erase(0,line.find_first_of(",)")+1);
    31843988                        }
     3989                        found = line.find("a(");
     3990//                      intvec *iv2=new intvec(this->numVars);
     3991//                      if(found!=string::npos)
     3992//                      {
     3993//                              line.erase(0,found+2);
     3994//                              string strweight2=line.substr(0,line.find_first_of(")"));       
     3995//                              for(int ii=0;ii<this->numVars;ii++)
     3996//                              {
     3997//                                      string weight;
     3998//                                      weight=line.substr(0,line.find_first_of(",)"));
     3999//                                      (*iv2)[ii]=atol(weight.c_str());
     4000//                                      line.erase(0,line.find_first_of(",)")+1);
     4001//                              }
     4002//                      }
     4003                       
    31854004                        ring newRing;
    31864005                        if(currRing->order[0]!=ringorder_a)
    31874006                        {
    3188                                 newRing=rCopyAndAddWeight(currRing,iv);
     4007//                              if(iv2!=NULL)
     4008//                                      newRing=rCopyAndAddWeight2(currRing,iv,iv2);
     4009//                              else
     4010                                        newRing=rCopyAndAddWeight(currRing,iv);
    31894011                        }
    31904012                        else
    3191                         {                       
    3192                                 newRing=rCopy0(currRing);
    3193                                 int length=this->numVars;
    3194                                 int *A=(int *)omAlloc0(length*sizeof(int));
    3195                                 for(int jj=0;jj<length;jj++)
    3196                                 {
    3197                                         A[jj]=(*iv)[jj];
    3198                                 }
    3199                                 omFree(newRing->wvhdl[0]);
    3200                                 newRing->wvhdl[0]=(int*)A;
    3201                                 newRing->block1[0]=length;
     4013                        {       
     4014//                              if(iv2==NULL)
     4015//                              {
     4016                                        newRing=rCopy0(currRing);
     4017                                        int length=this->numVars;
     4018                                        int *A=(int *)omAlloc0(length*sizeof(int));
     4019                                        for(int jj=0;jj<length;jj++)
     4020                                        {
     4021                                                A[jj]=(*iv)[jj];
     4022                                        }
     4023                                        omFree(newRing->wvhdl[0]);
     4024                                        newRing->wvhdl[0]=(int*)A;
     4025                                        newRing->block1[0]=length;
     4026//                              }
     4027//                              else
     4028//                              {
     4029//                                      newRing=rCopy0(currRing);
     4030//                                      int length=this->numVars;
     4031//                                      int *A1=(int *)omAlloc0(length*sizeof(int));
     4032//                                      int *A2=(int *)omAlloc0(length*sizeof(int));
     4033//                                      for(int jj=0;jj<length;jj++)
     4034//                                      {
     4035//                                              A1[jj]=(*iv2)[jj];
     4036//                                              A2[jj]=(*iv)[jj];                                                               
     4037//                                      }
     4038//                                      omFree(newRing->wvhdl[0]);
     4039//                                      newRing->wvhdl[0]=(int*)A1;
     4040//                                      newRing->block1[0]=length;
     4041//                                      if(newRing->wvhdl[1]!=NULL)
     4042//                                              omFree(newRing->wvhdl[0]);
     4043//                                      newRing->block1[1]=length;
     4044//                              }
    32024045                        }
    32034046                        delete iv;
     4047//                      delete iv2;
    32044048                        rComplete(newRing);
    32054049                        gc->baseRing=rCopy(newRing);
    … …  
    34684312float gcone::t_getExtremalRays;
    34694313float gcone::time_flip;
     4314float gcone::time_flip2;
    34704315float gcone::t_areEqual;
    34714316float gcone::t_markings;
    … …  
    35354380                        }                       
    35364381                        gcAct->getConeNormals(gcAct->gcBasis);
     4382                        gcone::dd_LinealitySpace = gcAct->computeLinealitySpace();
     4383                        gcAct->getExtremalRays(*gcAct);
    35374384                        gcAct->noRevS(*gcAct);  //Here we go!
    35384385                        //Switch back to the ring the computation was started in
    … …  
    35784425        cout << "  t_dd:" << gcone::t_dd << endl;
    35794426        cout << "  t_kStd:" << gcone::t_kStd << endl;
     4427        cout << "t_Flip2:" << gcone::time_flip2 << endl;
     4428        cout << "  t_dd:" << gcone::t_dd << endl;
     4429        cout << "  t_kStd:" << gcone::t_kStd << endl;
    35804430        cout << "t_computeInv:" << gcone::time_computeInv << endl;
    35814431        cout << "t_enqueue:" << gcone::time_enqueue << endl;

  • kernel/gfan.h

    r591e530 rbe7ab3f  
    8383                /**  The copy constructor */
    8484                facet(const facet& f);
    85                
     85                facet(const facet& f, bool shallow);
    8686                /** The default destructor */
    8787                ~facet();
    88                                
     88                /** Comparison operator*/
     89//              inline bool operator==(const facet *f,const facet *g);                 
    8990                /** \brief Comparison of facets*/
    9091                inline bool areEqual(facet *f, facet *g);
    … …  
    115116                 * prints the facet normal an all (codim-2)-facets that belong to it
    116117                 */
    117                 inline void fDebugPrint();
     118                volatile void fDebugPrint();
    118119                friend class gcone;             
    119120};
    … …  
    145146                static float t_getExtremalRays;
    146147                static float time_flip;
     148                static float time_flip2;
    147149                static float t_areEqual;
    148150                static float t_ffG;
    … …  
    205207                inline void invPrint(const ideal &I);
    206208                inline bool isMonomial(const ideal &I);
    207                 inline intvec *ivNeg(intvec *iv);
    208                 inline int dotProduct(intvec &iva, intvec &ivb);
     209                inline intvec *ivNeg(const intvec *iv);
     210//              inline int dotProduct(intvec &iva, intvec &ivb);
    209211                inline int dotProduct(const intvec &iva, const intvec &ivb);
    210                 inline bool isParallel(intvec &a, intvec &b);
     212                inline bool isParallel(const intvec &a, const intvec &b);
    211213//              inline int dotProduct(const intvec* a, const intvec *b);
    212214//              inline bool isParallel(const intvec* a, const intvec* b);
    213215                inline bool areEqual(intvec &a, intvec &b);
    214216                inline bool areEqual( facet *f,  facet *g);
     217                inline bool areEqual2(facet* f, facet *g);
    215218                inline int intgcd(int a, int b);
    216219                inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE);
    … …  
    223226                inline void getExtremalRays(const gcone &gc);
    224227                inline void flip(ideal gb, facet *f);
    225                 inline void computeInv(ideal &gb, ideal &inv, intvec &f);
    226 //              poly restOfDiv(poly const &f, ideal const &I); removed with r12286
     228                inline void flip2(const ideal gb, facet *f);
     229                inline void computeInv(const ideal &gb, ideal &inv, const intvec &f);//                 poly restOfDiv(poly const &f, ideal const &I); removed with r12286
    227230                inline ideal ffG(const ideal &H, const ideal &G);
    228231                inline void getGB(ideal const &inputIdeal);             
    229232                inline void interiorPoint( dd_MatrixPtr &M, intvec &iv);
    230 //              inline void interiorPoint2(const dd_MatrixPtr &M, intvec &iv);//removed Feb 8th, 2010
     233                inline void interiorPoint2(); //removed Feb 8th, 2010, new method Feb 19th, 2010
    231234                inline void preprocessInequalities(dd_MatrixPtr &M);
    232235                ring rCopyAndAddWeight(const ring &r, intvec *ivw);
     236                ring rCopyAndAddWeight2(const ring &, const intvec *, const intvec *);
    233237                ring rCopyAndChangeWeight(const ring &r, intvec *ivw);         
    234238//              void reverseSearch(gcone *gcAct); //NOTE both removed from r12286
    … …  
    242246                /** Compute the lineality space Ax=0 and return it as dd_MatrixPtr dd_LinealitySpace*/
    243247                inline dd_MatrixPtr computeLinealitySpace();
    244                 inline bool iv64isStrictlyPositive(intvec *);
     248                inline bool iv64isStrictlyPositive(const intvec *);
     249                /** Exchange 2 ordertype_a by just 1 */
     250                inline void replaceDouble_ringorder_a_ByASingleOne();
    245251//              static void gcone::idPrint(ideal &I);           
    246252                friend class facet;     
Note: See TracChangeset for help on using the changeset viewer.