Changeset be7ab3f in git
 Timestamp:
 Mar 3, 2010, 4:36:29 PM (13 years ago)
 Branches:
 (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
 Children:
 9f857dbf446399fd8232a210c134547a6b54d8a2
 Parents:
 591e530bebccea49fd635f208dd10c42a3b894d7
 Location:
 kernel
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

kernel/gfan.cc
r591e530 rbe7ab3f 13 13 #include "kstd1.h" 14 14 #include "kutil.h" //ksCreateSpoly 15 // #include "int vec.h"15 // #include "int64vec.h" 16 16 #include "polys.h" 17 17 #include "ideals.h" … … 31 31 #include <string> 32 32 #include <sstream> 33 #include <time.h>33 // #include <time.h> 34 34 #include <stdlib.h> 35 35 //#include <gmpxx.h> 36 36 // #include <vector> 37 #include <assert.h> 37 38 38 39 … … 141 142 this>UCN=f.UCN; 142 143 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 143 147 facet* f2Copy; 144 148 f2Copy=f.codim2Ptr; … … 169 173 } 170 174 } 175 176 /** \brief Shallow copy constructor for facets 177 */ 178 facet::facet(const facet& f, bool shallow) 179 {} 171 180 172 181 /** The default destructor */ … … 207 216 { 208 217 return(this>fNormal); 209 } 218 } 219 220 bool 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 + 1e6*(end.tv_usec  start.tv_usec)); 249 #endif 250 } 251 210 252 /** \brief Comparison of facets 211 253 * called from enqueueNewFacets … … 419 461 * prints the facet normal an all (codim2)facets that belong to it 420 462 */ 421 inline void facet::fDebugPrint()463 volatile void facet::fDebugPrint() 422 464 { 423 465 facet *codim2Act; … … 1279 1321 //Add lineality space  dd_LinealitySpace 1280 1322 dd_MatrixPtr ddineq; 1281 dd_ErrorType err; 1323 dd_ErrorType err; 1282 1324 if(dd_LinealitySpace>rowsize>0)//The linspace might be 0 1283 1325 ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace); 1284 1326 else 1285 1327 ddineq = dd_CopyMatrix(gc.ddFacets); 1328 /* In case the input is nonhomogeneous we add constrains for the positive orthant. 1329 * This is justified by the fact that for nonhomog 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 } 1286 1345 dd_PolyhedraPtr ddPolyh; 1287 1346 ddPolyh = dd_DDMatrix2Poly(ddineq, &err); 1288 1347 dd_MatrixPtr P; 1289 1348 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 1291 1421 //Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal 1292 1422 int rows=P>rowsize; … … 1303 1433 if(dotProduct(*fNormal,*rowvec)==0) 1304 1434 { 1305 intvec *tmp = ivIntPointOfFacet; 1435 intvec *tmp = ivIntPointOfFacet;//Prevent memleak 1306 1436 fAct>numCodim2Facets++; 1307 1437 facet *codim2Act; … … 1318 1448 codim2Act>setFacetNormal(rowvec); 1319 1449 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 1321 1452 ivIntPointOfFacet=ivAdd(ivIntPointOfFacet,rowvec); 1453 //Now tmp still points to the OLD address of ivIntPt 1322 1454 delete(tmp); 1323 1455 … … 1325 1457 delete(rowvec); 1326 1458 } 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 1328 1485 delete(ivIntPointOfFacet); 1329 1486 //Now (if we have at least 3 variables) do a bubblesort on the rays … … 1368 1525 fAct = fAct>next; 1369 1526 } 1527 dd_FreeMatrix(P); 1370 1528 //Now all extremal rays should be set w.r.t their respective fNormal 1371 1529 //TODO Not sufficient > vol2 II/125&127 1372 1530 //NOTE Sufficient according to cddlibs doc. These ARE rays 1531 //What the hell... let's just take interior points 1373 1532 if(gcone::hasHomInput==FALSE) 1374 1533 { … … 1376 1535 while(fAct!=NULL) 1377 1536 { 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)) 1410 1557 fAct>isFlippable=FALSE; 1411 1558 fAct = fAct>next; … … 1418 1565 } 1419 1566 1420 inline bool gcone::iv64isStrictlyPositive( intvec * iv64)1567 inline bool gcone::iv64isStrictlyPositive(const intvec * iv64) 1421 1568 { 1422 1569 bool res=TRUE; … … 1427 1574 res=FALSE; 1428 1575 break; 1429 } 1576 } 1430 1577 } 1431 1578 return res; … … 1453 1600 intvec *fNormal;// = new intvec(this>numVars); //facet normal, check for parallelity 1454 1601 fNormal = f>getFacetNormal(); //read this>fNormal; 1455 1602 #ifdef gfan_DEBUG 1456 1603 // 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 1462 1610 if(this>getUCN() != f>getUCN()) 1463 1611 { … … 1588 1736 { 1589 1737 #ifdef gfan_DEBUG 1590 // 1738 // cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl; 1591 1739 #endif 1592 1740 if (src_ExpV[kk]!=dst_ExpV[kk]) … … 1594 1742 expVAreEqual=FALSE; 1595 1743 } 1596 } 1597 //if (*src_ExpV == *dst_ExpV) 1744 } 1598 1745 if (expVAreEqual==TRUE) 1599 1746 { … … 1649 1796 omFree(v); 1650 1797 } 1651 // omFree(v);1652 1798 omFree(leadExpV); 1653 1799 }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++) … … 1710 1856 //NOTE Here we should remove interiorPoint and instead 1711 1857 // 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(); 1713 1862 dd_FreeMatrix(intPointMatrix); 1714 1863 /*Crude attempt for interior point */ … … 1734 1883 gettimeofday(&t_dd_end, 0); 1735 1884 t_dd += (t_dd_end.tv_sec  t_dd_start.tv_sec + 1e6*(t_dd_end.tv_usec  t_dd_start.tv_usec)); 1736 #endif 1885 #endif 1886 1737 1887 /*Step 3 1738 1888 * turn the minimal basis into a reduced one */ … … 1803 1953 }//void flip(ideal gb, facet *f) 1804 1954 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 */ 1970 inline 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 + 1e6*(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 + 1e6*(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 + 1e6*(end.tv_usec  start.tv_usec)); 2124 #endif 2125 }//flip2 2126 1805 2127 /** \brief Compute initial ideal 1806 2128 * Compute the initial ideal in_v(G) wrt a (possible) facet normal … … 1808 2130 * and in gcone::flip for obvious reasons. 1809 2131 */ 1810 inline void gcone::computeInv( ideal &gb, ideal &initialForm,intvec &fNormal)2132 inline void gcone::computeInv(const ideal &gb, ideal &initialForm, const intvec &fNormal) 1811 2133 { 1812 2134 #ifdef gfanp … … 2000 2322 /** \brief Compute the negative of a given intvec 2001 2323 */ 2002 inline intvec *gcone::ivNeg( intvec *iv)2324 inline intvec *gcone::ivNeg(const intvec *iv) 2003 2325 { //Hm, switching to intvec const intvec does no longer work 2004 2326 intvec *res;// = new intvec(iv>length()); … … 2012 2334 * 2013 2335 */ 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 dotProduct2336 // 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 2023 2345 inline int gcone::dotProduct(const intvec &iva, const intvec &ivb) 2024 2346 { … … 2026 2348 for (int i=0;i<this>numVars;i++) 2027 2349 { 2350 // #ifndef NDEBUG 2351 // (const_cast<intvec*>(&iva))>show(1,0); (const_cast<intvec*>(&ivb))>show(1,0); 2352 // #endif 2028 2353 res = res+(iva[i]*ivb[i]); 2029 2354 } … … 2034 2359 * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$ 2035 2360 */ 2036 inline bool gcone::isParallel( intvec &a,intvec &b)2361 inline bool gcone::isParallel(const intvec &a,const intvec &b) 2037 2362 { 2038 2363 int lhs,rhs; … … 2180 2505 }//void interiorPoint(dd_MatrixPtr const &M) 2181 2506 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+(ab)/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 */ 2522 inline 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+((ba)/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+((ba)/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=ba 2582 mpq_t quot; 2583 mpq_init(quot); 2584 mpq_div_2exp(quot,diff,1); //quot=diff/2=(ba)/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+(ba)/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=ab 2815 // mpq_t quot; 2816 // mpq_init(quot); 2817 // mpq_div_2exp(quot,diff,1); //quot=diff/2=(ab)/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+(ab)/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 } 2183 2873 2184 2874 /** \brief Copy a ring and add a weighted ordering in first place … … 2219 2909 return res; 2220 2910 }//rCopyAndAdd 2911 2912 ring 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 } 2221 2953 2222 2954 ring rCopyAndChangeWeight(ring const &r, intvec *ivw) … … 2414 3146 /*Compute lineality space here 2415 3147 * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/ 2416 dd_LinealitySpace = gcAct>computeLinealitySpace(); 3148 if(dd_LinealitySpace==NULL) 3149 dd_LinealitySpace = gcAct>computeLinealitySpace(); 2417 3150 /*End of lineality space computation*/ 2418 3151 // gcAct>getCodim2Normals(*gcAct); 2419 gcAct>getExtremalRays(*gcAct); 3152 if(fAct>codim2Ptr==NULL) 3153 gcAct>getExtremalRays(*gcAct); 2420 3154 2421 3155 //Compute unique representation of Facets and rays, i.e. primitive vectors … … 2448 3182 SearchListAct>numCodim2Facets=fAct>numCodim2Facets; 2449 3183 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); 2450 3189 //Copy codim2facets 2451 3190 facet *codim2Act; … … 2504 3243 We always choose the first facet from SearchList as facet to be flipped 2505 3244 */ 2506 while((SearchListAct!=NULL))// && gcone::counter<10)3245 while((SearchListAct!=NULL))// && counter<363) 2507 3246 {//NOTE See to it that the cone is only changed after ALL facets have been flipped! 2508 3247 fAct = SearchListAct; … … 2511 3250 // while( (fAct>getUCN() == fAct>next>getUCN()) ) 2512 3251 { //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); 2514 3253 //NOTE rCopy needed? 2515 3254 ring rTmp=rCopy(fAct>flipRing); … … 2534 3273 gcTmp>normalize();// will be done in g2c 2535 3274 //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(); 2536 3280 dd_FreeMatrix(gcTmp>ddFacets); 2537 3281 #ifdef gfan_DEBUG … … 2826 3570 if(slAct==NULL) 2827 3571 { 2828 slAct = new facet(*fCopy); 3572 slAct = new facet(*fCopy);//copy constructor 2829 3573 slHead = slAct; 2830 3574 } … … 2915 3659 slEnd>isFlippable = TRUE; 2916 3660 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; 2917 3667 slEnd>prev = marker; 2918 3668 //Copy codim2facets … … 2957 3707 return slHead; 2958 3708 }//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 */ 3717 inline 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 } 2960 3760 //A try with the STL 2961 3761 facet * gcone::enqueue2(facet *f) … … 3074 3874 f2Act=fAct>codim2Ptr; 3075 3875 3076 char *ringString = rString( currRing);3876 char *ringString = rString(gc.baseRing); 3077 3877 3078 3878 if (!gcOutputFile) … … 3141 3941 3142 3942 /** \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 */ 3144 3947 inline void gcone::readConeFromFile(int UCN, gcone *gc) 3145 3948 { … … 3175 3978 string strweight; 3176 3979 strweight=line.substr(0,line.find_first_of(")")); 3177 intvec *iv=new intvec(this>numVars); 3980 3981 intvec *iv=new intvec(this>numVars);// 3178 3982 for(int ii=0;ii<this>numVars;ii++) 3179 3983 { … … 3183 3987 line.erase(0,line.find_first_of(",)")+1); 3184 3988 } 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 3185 4004 ring newRing; 3186 4005 if(currRing>order[0]!=ringorder_a) 3187 4006 { 3188 newRing=rCopyAndAddWeight(currRing,iv); 4007 // if(iv2!=NULL) 4008 // newRing=rCopyAndAddWeight2(currRing,iv,iv2); 4009 // else 4010 newRing=rCopyAndAddWeight(currRing,iv); 3189 4011 } 3190 4012 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 // } 3202 4045 } 3203 4046 delete iv; 4047 // delete iv2; 3204 4048 rComplete(newRing); 3205 4049 gc>baseRing=rCopy(newRing); … … 3468 4312 float gcone::t_getExtremalRays; 3469 4313 float gcone::time_flip; 4314 float gcone::time_flip2; 3470 4315 float gcone::t_areEqual; 3471 4316 float gcone::t_markings; … … 3535 4380 } 3536 4381 gcAct>getConeNormals(gcAct>gcBasis); 4382 gcone::dd_LinealitySpace = gcAct>computeLinealitySpace(); 4383 gcAct>getExtremalRays(*gcAct); 3537 4384 gcAct>noRevS(*gcAct); //Here we go! 3538 4385 //Switch back to the ring the computation was started in … … 3578 4425 cout << " t_dd:" << gcone::t_dd << endl; 3579 4426 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; 3580 4430 cout << "t_computeInv:" << gcone::time_computeInv << endl; 3581 4431 cout << "t_enqueue:" << gcone::time_enqueue << endl; 
kernel/gfan.h
r591e530 rbe7ab3f 83 83 /** The copy constructor */ 84 84 facet(const facet& f); 85 85 facet(const facet& f, bool shallow); 86 86 /** The default destructor */ 87 87 ~facet(); 88 88 /** Comparison operator*/ 89 // inline bool operator==(const facet *f,const facet *g); 89 90 /** \brief Comparison of facets*/ 90 91 inline bool areEqual(facet *f, facet *g); … … 115 116 * prints the facet normal an all (codim2)facets that belong to it 116 117 */ 117 inline void fDebugPrint();118 volatile void fDebugPrint(); 118 119 friend class gcone; 119 120 }; … … 145 146 static float t_getExtremalRays; 146 147 static float time_flip; 148 static float time_flip2; 147 149 static float t_areEqual; 148 150 static float t_ffG; … … 205 207 inline void invPrint(const ideal &I); 206 208 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); 209 211 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); 211 213 // inline int dotProduct(const intvec* a, const intvec *b); 212 214 // inline bool isParallel(const intvec* a, const intvec* b); 213 215 inline bool areEqual(intvec &a, intvec &b); 214 216 inline bool areEqual( facet *f, facet *g); 217 inline bool areEqual2(facet* f, facet *g); 215 218 inline int intgcd(int a, int b); 216 219 inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE); … … 223 226 inline void getExtremalRays(const gcone &gc); 224 227 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 r12286228 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 227 230 inline ideal ffG(const ideal &H, const ideal &G); 228 231 inline void getGB(ideal const &inputIdeal); 229 232 inline void interiorPoint( dd_MatrixPtr &M, intvec &iv); 230 // inline void interiorPoint2(const dd_MatrixPtr &M, intvec &iv);//removed Feb 8th, 2010233 inline void interiorPoint2(); //removed Feb 8th, 2010, new method Feb 19th, 2010 231 234 inline void preprocessInequalities(dd_MatrixPtr &M); 232 235 ring rCopyAndAddWeight(const ring &r, intvec *ivw); 236 ring rCopyAndAddWeight2(const ring &, const intvec *, const intvec *); 233 237 ring rCopyAndChangeWeight(const ring &r, intvec *ivw); 234 238 // void reverseSearch(gcone *gcAct); //NOTE both removed from r12286 … … 242 246 /** Compute the lineality space Ax=0 and return it as dd_MatrixPtr dd_LinealitySpace*/ 243 247 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(); 245 251 // static void gcone::idPrint(ideal &I); 246 252 friend class facet;
Note: See TracChangeset
for help on using the changeset viewer.