Changeset dd2b206 in git
- Timestamp:
- Jul 13, 2009, 8:50:42 AM (14 years ago)
- Branches:
- (u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
- Children:
- 35c996104b9fc1ab5d679c8201904b6bb15f508e
- Parents:
- a9ca4ab64e5d10dd0a5613e146923d5b5c027a1b
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
kernel/gfan.cc
ra9ca4a rdd2b206 2 2 Compute the Groebner fan of an ideal 3 3 $Author: monerjan $ 4 $Date: 2009-07- 09 09:59:10$5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.7 4 2009-07-09 09:59:10monerjan Exp $6 $Id: gfan.cc,v 1.7 4 2009-07-09 09:59:10monerjan Exp $4 $Date: 2009-07-13 06:50:42 $ 5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.75 2009-07-13 06:50:42 monerjan Exp $ 6 $Id: gfan.cc,v 1.75 2009-07-13 06:50:42 monerjan Exp $ 7 7 */ 8 8 … … 75 75 76 76 /** \brief Universal Cone Number 77 * The number of the cone the facet belongs to 77 * The number of the cone the facet belongs to, Set in getConeNormals() 78 78 */ 79 79 int UCN; … … 93 93 94 94 public: 95 //bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i]96 bool isIncoming; //Is the facet incoming or outgoing ?95 bool isFlippable; //flippable facet? 96 bool isIncoming; //Is the facet incoming or outgoing in the reverse search? 97 97 facet *next; //Pointer to next facet 98 98 facet *prev; //Pointer to predecessor. Needed for the SearchList in noRevS … … 117 117 this->prev=NULL; 118 118 this->flipRing=NULL; 119 this->isFlippable=FALSE; 119 120 } 120 121 … … 137 138 this->prev=NULL; 138 139 this->flipRing=NULL; 140 this->isFlippable=FALSE; 139 141 } 140 142 … … 244 246 return this->interiorPoint; 245 247 } 248 249 void fDebugPrint() 250 { 251 //facet *f; 252 facet *codim2Act; 253 //f = this; 254 codim2Act = this->codim2Ptr; 255 intvec *fNormal; 256 fNormal = this->getFacetNormal(); 257 intvec *f2Normal; 258 cout << "=======================" << endl; 259 cout << "Facet normal = ("; 260 for(int ii=0;ii<pVariables;ii++) 261 { 262 cout << (*fNormal)[ii] << ","; 263 } 264 if(this->isFlippable==TRUE) 265 cout << ")" << endl; 266 else 267 cout << ")*" << endl; //This case should never happen! 268 cout << "-----------------------" << endl; 269 cout << "Codim2 facets:" << endl; 270 while(codim2Act!=NULL) 271 { 272 f2Normal = codim2Act->getFacetNormal(); 273 cout << "("; 274 for(int jj=0;jj<pVariables;jj++) 275 { 276 cout << (*f2Normal)[jj] << ","; 277 } 278 cout << ")" << endl; 279 codim2Act = codim2Act->next; 280 } 281 cout << "=======================" << endl; 282 } 246 283 247 284 /*bool isFlippable(intvec &load) … … 456 493 } 457 494 458 intvec *iv;// = new intvec(this->numVars); 459 460 //while (f->next!=NULL) 495 intvec *iv; 461 496 while(f!=NULL) 462 497 { 463 498 iv = f->getFacetNormal(); 464 iv->show(); 499 cout << "("; 500 iv->show(1,0); 501 if(f->isFlippable==FALSE) 502 cout << ")* "; 503 else 504 cout << ") "; 465 505 f=f->next; 466 506 } 467 //delete iv; 507 cout << endl; 468 508 } 469 509 … … 473 513 facet *fAct; 474 514 fAct = &f; 475 intvec *n;// = new intvec(this->numVars); 515 facet *codim2Act; 516 codim2Act = fAct->codim2Ptr; 517 intvec *fNormal;// = new intvec(this->numVars); 518 intvec *f2Normal; 476 519 cout << endl; 477 520 while(fAct!=NULL) 478 521 { 479 n=fAct->getFacetNormal(); 480 n->show(); 522 fNormal=fAct->getFacetNormal(); 523 cout << "("; 524 fNormal->show(1,0); 525 cout << ") "; 526 codim2Act = fAct->codim2Ptr; 527 cout << " Codim2: "; 528 while(codim2Act!=NULL) 529 { 530 f2Normal = codim2Act->getFacetNormal(); 531 cout << "("; 532 f2Normal->show(1,0); 533 cout << ") "; 534 codim2Act = codim2Act->next; 535 } 536 cout << "UCN = " << fAct->getUCN() << endl; 481 537 fAct = fAct->next; 482 cout << endl;483 cout << "---" << endl;484 485 538 } 486 539 } … … 600 653 601 654 #ifdef gfan_DEBUG 602 //cout << "The inequality matrix is" << endl;603 //dd_WriteMatrix(stdout, ddineq);655 cout << "The inequality matrix is" << endl; 656 dd_WriteMatrix(stdout, ddineq); 604 657 #endif 605 658 … … 624 677 ddFacets = dd_CopyMatrix(ddineq); 625 678 #ifdef gfan_DEBUG 626 //cout << "Having removed redundancies, the normals now read:" << endl;627 //dd_WriteMatrix(stdout,ddineq);679 cout << "Having removed redundancies, the normals now read:" << endl; 680 dd_WriteMatrix(stdout,ddineq); 628 681 // cout << "Rows = " << ddrows << endl; 629 682 // cout << "Cols = " << ddcols << endl; … … 655 708 656 709 /*Quick'n'dirty hack for flippability*/ 657 bool isFlip pable=FALSE;710 bool isFlip=FALSE; 658 711 for (int jj = 0; jj<load->length(); jj++) 659 712 { … … 664 717 //if (ivMult(load,ivCanonical)<0) 665 718 { 666 isFlip pable=TRUE;719 isFlip=TRUE; 667 720 break; //URGHS 668 721 } 669 722 } 670 if (isFlippable==FALSE) 671 { 672 #ifdef gfan_DEBUG 673 cout << "Ignoring facet" << endl; 674 load->show(); 675 cout << endl; 676 //fAct->next=NULL; 723 if (isFlip==FALSE) 724 { 725 this->numFacets++; 726 if(this->numFacets==1) 727 { 728 facet *fRoot = new facet(); 729 this->facetPtr = fRoot; 730 fAct = fRoot; 731 732 } 733 else 734 { 735 fAct->next = new facet(); 736 fAct = fAct->next; 737 } 738 fAct->isFlippable=FALSE; 739 fAct->setFacetNormal(load); 740 fAct->setUCN(this->getUCN()); 741 #ifdef gfan_DEBUG 742 cout << "Marking facet ("; 743 load->show(1,0); 744 cout << ") as non flippable" << endl; 677 745 #endif 678 746 } … … 691 759 fAct = fAct->next; 692 760 } 761 fAct->isFlippable=TRUE; 693 762 fAct->setFacetNormal(load); 694 763 fAct->setUCN(this->getUCN()); … … 760 829 { 761 830 ddakt = dd_CopyMatrix(ddineq); 831 ddakt->representation=dd_Inequality; 762 832 set_addelem(ddakt->linset,ii+1); 763 764 dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);833 //dd_WriteMatrix(stdout,ddakt); 834 //dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err); 765 835 766 836 #ifdef gfan_DEBUG 767 // cout << "Codim2 matrix"<<endl;768 // dd_WriteMatrix(stdout,ddakt);837 // cout << "Codim2 matrix"<<endl; 838 // dd_WriteMatrix(stdout,ddakt); 769 839 #endif 770 840 ddpolyh=dd_DDMatrix2Poly(ddakt, &err); 771 P=dd_CopyGenerators(ddpolyh); 772 #ifdef gfan_DEBUG 773 // cout << "Codim2 facet:" << endl;774 // dd_WriteMatrix(stdout,P);775 // cout << endl;841 P=dd_CopyGenerators(ddpolyh); 842 #ifdef gfan_DEBUG 843 // cout << "Codim2 facet:" << endl; 844 // dd_WriteMatrix(stdout,P); 845 // cout << endl; 776 846 #endif 777 847 … … 833 903 fNormal = f->getFacetNormal(); //read this->fNormal; 834 904 #ifdef gfan_DEBUG 835 std::cout << "===" << std::endl;905 //std::cout << "===" << std::endl; 836 906 std::cout << "running gcone::flip" << std::endl; 837 907 std::cout << "flipping" << endl; … … 840 910 pWrite((poly)gb->m[ii]); 841 911 } 842 cout << "over facet " << endl;843 fNormal->show( );912 cout << "over facet "; 913 fNormal->show(1,0); 844 914 std::cout << std::endl; 845 915 #endif … … 883 953 }//while 884 954 #ifdef gfan_DEBUG 885 cout << "Initial Form=";886 pWrite(initialFormElement[ii]);887 cout << "---" << endl;955 // cout << "Initial Form="; 956 // pWrite(initialFormElement[ii]); 957 // cout << "---" << endl; 888 958 #endif 889 959 /*Now initialFormElement must be added to (ideal)initialForm */ … … 891 961 }//for 892 962 #ifdef gfan_DEBUG 893 cout << "Initial ideal is: " << endl;963 /* cout << "Initial ideal is: " << endl; 894 964 idShow(initialForm); 895 //f->printFlipGB(); 896 cout << "===" << endl;965 //f->printFlipGB();*/ 966 // cout << "===" << endl; 897 967 #endif 898 968 //delete check; … … 1078 1148 */ 1079 1149 //ring dstRing=rCopyAndAddWeight(tmpRing,iv_weight); 1080 1081 // int i,j;1082 // ring dstRing=rCopy0(srcRing);1083 // i=rBlocks(srcRing);1084 //1085 // dstRing->order=(int *)omAlloc((i+1)*sizeof(int));1086 // for(j=i;j>0;j--)1087 // {1088 // dstRing->order[j]=srcRing->order[j-1];1089 // dstRing->block0[j]=srcRing->block0[j-1];1090 // dstRing->block1[j]=srcRing->block1[j-1];1091 // if (srcRing->wvhdl[j-1] != NULL)1092 // {1093 // dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);1094 // }1095 // }1096 // dstRing->order[0]=ringorder_a;1097 // dstRing->order[1]=ringorder_dp;1098 // dstRing->order[2]=ringorder_C;1099 // dstRing->wvhdl[0] =( int *)omAlloc((iv_weight->length())*sizeof(int));1100 //1101 // for (int ii=0;ii<this->numVars;ii++)1102 // {1103 // dstRing->wvhdl[0][ii]=(*iv_weight)[ii];1104 // }1105 // rComplete(dstRing);1106 1150 1107 1151 // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a … … 1283 1327 //(res)[ii] = (*res[ii])*(int)(-1); 1284 1328 //} 1285 res->show( );1329 res->show(1,0); 1286 1330 return res; 1287 1331 } … … 1693 1737 for(int ii=0;ii<this->numFacets;ii++) 1694 1738 { 1695 fNormal = fAct->getFacetNormal(); 1696 if(ii==0) 1697 { 1698 //facet *SearchListRoot = new facet(); 1699 SearchListAct = SearchListRoot; 1700 //memcpy(SearchListAct,fAct,sizeof(facet)); 1701 } 1702 else 1703 { 1704 SearchListAct->next = new facet(); 1705 SearchListAct = SearchListAct->next; 1706 //memcpy(SearchListAct,fAct,sizeof(facet)); 1707 } 1708 SearchListAct->setFacetNormal(fNormal); 1709 SearchListAct->setUCN(this->getUCN()); 1710 SearchListAct->numCodim2Facets=fAct->numCodim2Facets; 1711 1712 //Copy codim2-facets 1713 codim2Act=fAct->codim2Ptr; 1714 SearchListAct->codim2Ptr = new facet(2); 1715 sl2Root = SearchListAct->codim2Ptr; 1716 sl2Act = sl2Root; 1717 //while(codim2Act!=NULL) 1718 for(int jj=0;jj<fAct->numCodim2Facets;jj++) 1719 { 1720 f2Normal = codim2Act->getFacetNormal(); 1721 if(jj==0) 1739 //only copy flippable facets! 1740 if(fAct->isFlippable==TRUE) 1741 { 1742 fNormal = fAct->getFacetNormal(); 1743 if(ii==0) 1722 1744 { 1723 sl2Act = sl2Root;1724 sl2Act->setFacetNormal(f2Normal);1745 SearchListAct = SearchListRoot; 1746 //memcpy(SearchListAct,fAct,sizeof(facet)); 1725 1747 } 1726 1748 else 1727 1749 { 1728 sl2Act->next = new facet(2); 1729 sl2Act = sl2Act->next; 1730 sl2Act->setFacetNormal(f2Normal); 1731 } 1732 codim2Act = codim2Act->next; 1733 } 1734 fAct = fAct->next; 1750 SearchListAct->next = new facet(); 1751 SearchListAct = SearchListAct->next; 1752 //memcpy(SearchListAct,fAct,sizeof(facet)); 1753 } 1754 SearchListAct->setFacetNormal(fNormal); 1755 SearchListAct->setUCN(this->getUCN()); 1756 SearchListAct->numCodim2Facets=fAct->numCodim2Facets; 1757 SearchListAct->isFlippable=TRUE; 1758 //Copy codim2-facets 1759 codim2Act=fAct->codim2Ptr; 1760 SearchListAct->codim2Ptr = new facet(2); 1761 sl2Root = SearchListAct->codim2Ptr; 1762 sl2Act = sl2Root; 1763 //while(codim2Act!=NULL) 1764 for(int jj=0;jj<fAct->numCodim2Facets;jj++) 1765 { 1766 f2Normal = codim2Act->getFacetNormal(); 1767 if(jj==0) 1768 { 1769 sl2Act = sl2Root; 1770 sl2Act->setFacetNormal(f2Normal); 1771 } 1772 else 1773 { 1774 sl2Act->next = new facet(2); 1775 sl2Act = sl2Act->next; 1776 sl2Act->setFacetNormal(f2Normal); 1777 } 1778 codim2Act = codim2Act->next; 1779 } 1780 fAct = fAct->next; 1781 }//if(fAct->isFlippable==TRUE) 1782 else {fAct = fAct->next;} 1735 1783 } 1736 1784 … … 1763 1811 //do 1764 1812 while(fAct!=NULL) 1765 { 1813 { //Since SLA should only contain flippables there should be no need to check for that 1766 1814 gcAct->flip(gcAct->gcBasis,fAct); 1767 1815 ring rTmp=rCopy(fAct->flipRing); … … 1801 1849 UCNcounter++; 1802 1850 //SearchListAct = SearchListAct->next; 1803 SearchListAct = fAct->next; 1851 //SearchListAct = fAct->next; 1852 SearchListAct = SearchListRoot; 1804 1853 } 1805 1854 … … 1958 2007 facet *sl2Act; 1959 2008 sl2Act = f.codim2Ptr; 1960 facet *tBARoot; 1961 tBARoot = NULL; 1962 facet *tBAAct; 1963 tBAAct = NULL; 2009 1964 2010 bool doNotAdd=FALSE; 1965 2011 int ctr=0; //encountered qualities in SLA … … 1973 2019 slEndStatic = slEnd; 1974 2020 /*1st step: compare facetNormals*/ 1975 intvec *fNormal ; //= new intvec(this->numVars);1976 intvec *f2Normal ; //= new intvec(this->numVars);1977 intvec *slNormal ; //= new intvec(this->numVars);1978 intvec *sl2Normal ; //= new intvec(this->numVars);2021 intvec *fNormal=NULL; //= new intvec(this->numVars); 2022 intvec *f2Normal=NULL; //= new intvec(this->numVars); 2023 intvec *slNormal=NULL; //= new intvec(this->numVars); 2024 intvec *sl2Normal=NULL; //= new intvec(this->numVars); 1979 2025 1980 2026 while(fAct!=NULL) 1981 2027 { 1982 doNotAdd=TRUE; 1983 fNormal=fAct->getFacetNormal(); 1984 slAct = slHead; 1985 while(slAct!=NULL) 1986 { 1987 slNormal = slAct->getFacetNormal(); 1988 if(!isParallel(fNormal, slNormal)) 1989 { 1990 notParallelCtr++; 1991 } 1992 else 1993 { 1994 codim2Act = fAct->codim2Ptr; 1995 while(codim2Act!=NULL) 2028 if(fAct->isFlippable==TRUE) 2029 { 2030 doNotAdd=TRUE; 2031 fNormal=fAct->getFacetNormal(); 2032 slAct = slHead; 2033 notParallelCtr=0; 2034 while(slAct!=NULL) 2035 { 2036 slNormal = slAct->getFacetNormal(); 2037 #ifdef gfan_DEBUG 2038 cout << "Checking facet ("; 2039 fNormal->show(1,1); 2040 cout << ") against ("; 2041 slNormal->show(1,1); 2042 cout << ")" << endl; 2043 #endif 2044 if(!isParallel(fNormal, slNormal)) 1996 2045 { 1997 f2Normal = codim2Act->getFacetNormal(); 1998 sl2Act = f.codim2Ptr; 1999 while(sl2Act!=NULL) 2046 notParallelCtr++; 2047 } 2048 else 2049 { 2050 codim2Act = fAct->codim2Ptr; 2051 ctr=0; 2052 while(codim2Act!=NULL) 2000 2053 { 2001 sl2Normal = sl2Act->getFacetNormal(); 2002 if(areEqual(f2Normal, sl2Normal)) 2003 ctr++; 2004 sl2Act = sl2Act->next; 2005 }//while(sl2Act!=NULL) 2006 codim2Act = codim2Act->next; 2007 }//while(codim2Act!=NULL) 2008 if(ctr==fAct->numCodim2Facets) //facets ARE equal 2009 { 2010 if(slAct==slHead) //We want to delete the first element of SearchList 2054 f2Normal = codim2Act->getFacetNormal(); 2055 //sl2Act = f.codim2Ptr; 2056 sl2Act = slAct->codim2Ptr; 2057 while(sl2Act!=NULL) 2058 { 2059 sl2Normal = sl2Act->getFacetNormal(); 2060 if(areEqual(f2Normal, sl2Normal)) 2061 ctr++; 2062 sl2Act = sl2Act->next; 2063 }//while(sl2Act!=NULL) 2064 codim2Act = codim2Act->next; 2065 }//while(codim2Act!=NULL) 2066 if(ctr==fAct->numCodim2Facets) //facets ARE equal 2011 2067 { 2012 slHead = slAct->next; 2013 slHead->prev = NULL; 2014 //set a bool flag to mark slAct as to be deleted 2068 #ifdef gfan_DEBUG 2069 cout << "Removing facet ("; 2070 slNormal->show(1,0); 2071 cout << ") from SLA:" << endl; 2072 fAct->fDebugPrint(); 2073 slAct->fDebugPrint(); 2074 #endif 2075 if(slAct==slHead) //We want to delete the first element of SearchList 2076 { 2077 slHead = slAct->next; 2078 slHead->prev = NULL; 2079 //set a bool flag to mark slAct as to be deleted 2080 } 2081 else if(slAct==slEndStatic) 2082 { 2083 if(slEndStatic->next==NULL) 2084 { 2085 slEndStatic = slEndStatic->prev; 2086 slEndStatic->next = NULL; 2087 } 2088 else //we already added a facet after slEndStatic 2089 { 2090 slEndStatic->prev->next = slEndStatic->next; 2091 slEndStatic = slEndStatic->prev; 2092 slEnd = slEndStatic; 2093 } 2094 } 2095 else 2096 { 2097 slAct->prev->next = slAct->next; 2098 } 2099 //update lengthOfSearchList 2100 lengthOfSearchList--; 2101 break; 2102 }//if(ctr==fAct->numCodim2Facets) 2103 else //facets are NOT equal 2104 { 2105 doNotAdd=FALSE; 2106 break; 2015 2107 } 2016 else if(slAct==slEndStatic) 2017 { 2018 if(slEndStatic->next==NULL) 2019 { 2020 slEndStatic = slEndStatic->prev; 2021 slEndStatic->next = NULL; 2022 } 2023 else //we already added a facet after slEndStatic 2024 { 2025 slEndStatic->prev->next = slEndStatic->next; 2026 slEndStatic = slEndStatic->prev; 2027 slEnd = slEndStatic; 2028 } 2029 } 2030 else 2031 { 2032 slAct->prev->next = slAct->next; 2033 } 2034 //update lengthOfSearchList 2035 lengthOfSearchList--; 2036 break; 2037 }//if(ctr==fAct->numCodim2Facets) 2038 else //facets are NOT equal 2039 { 2040 doNotAdd=FALSE; 2041 break; 2042 } 2043 }//if(!isParallel(fNormal, slNormal)) 2044 slAct = slAct->next; 2045 //if slAct was marked as to be deleted, delete it here! 2046 }//while(slAct!=NULL) 2108 }//if(!isParallel(fNormal, slNormal)) 2109 slAct = slAct->next; 2110 //if slAct was marked as to be deleted, delete it here! 2111 }//while(slAct!=NULL) 2047 2112 if( (notParallelCtr==lengthOfSearchList) || (doNotAdd==FALSE) ) 2048 2113 { 2114 #ifdef gfan_DEBUG 2115 cout << "Adding facet ("; 2116 fNormal->show(1,0); 2117 cout << ") to SLA " << endl; 2118 #endif 2049 2119 //Add fAct to SLA 2050 2120 facet *marker; … … 2061 2131 2062 2132 slEnd->setUCN(this->getUCN()); 2133 slEnd->isFlippable = TRUE; 2063 2134 slEnd->setFacetNormal(fNormal); 2064 2135 slEnd->prev = marker; … … 2088 2159 } 2089 2160 fAct = fAct->next; 2161 } 2162 else 2163 { 2164 fAct = fAct->next; 2165 } 2090 2166 }//while(fAct!=NULL) 2091 2167 return slHead;
Note: See TracChangeset
for help on using the changeset viewer.