Changeset 1f3450 in git
- Timestamp:
- Feb 5, 2010, 6:02:12 PM (13 years ago)
- Branches:
- (u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
- Children:
- b086f1d24cfeb668d2a93ab34efd9d448c31be12
- Parents:
- 55587862408d1e92ee43d342befb2d834ea91635
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
kernel/gfan.cc
r555878 r1f3450 174 174 /* Cleanup the codim2-structure */ 175 175 if(this->codim==2) 176 // if(this->codim2Ptr!=NULL) 176 177 { 177 178 facet *codim2Ptr; … … 185 186 } 186 187 } 188 // delete this->codim2Ptr; 187 189 } 188 190 if(this->codim2Ptr!=NULL) … … 191 193 idDelete((ideal *)&this->flipGB); 192 194 if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb) 193 //rDelete(this->flipRing);195 rDelete(this->flipRing); 194 196 // this->flipRing=NULL; 195 197 this->prev=NULL; … … 317 319 inline int facet::getUCN() 318 320 { 319 if(this!=NULL)// && this!=(facet *)0xfbfbfbfbfbfbfbfb)// || this!=(facet *)0xfbfbfbfb) ) 320 // if(this!=NULL && ( this->fNormal!=(int64vec *)0xfbfbfbfb || this->fNormal!=(int64vec *)0xfbfbfbfbfbfbfbfb) ) 321 #ifndef NDEBUG 322 #if SIZEOF_LONG==8 323 if((this!=NULL && this!=(facet * const)0xfbfbfbfbfbfbfbfb)) 324 #elif 325 if((this!=NULL && this!=(facet * const)0xfbfbfbfb)) 326 #endif 327 #endif 328 #ifdef NDEBUG 329 if(this!=NULL) 330 #endif 321 331 return this->UCN; 322 332 else … … 435 445 this->pred=gc.UCN; 436 446 this->facetPtr=NULL; 437 // this->gcBasis=idCopy(f.flipGB);438 // ring saveRing=currRing;439 // ring tmpRing=f.flipRing;440 // rComplete(f.flipRing);441 // rChangeCurrRing(f.flipRing);442 447 this->gcBasis=idrCopyR(f.flipGB, f.flipRing); 443 448 // this->inputIdeal=idCopy(this->gcBasis); 444 449 this->baseRing=rCopy(f.flipRing); 445 // rComplete(this->baseRing);446 // rChangeCurrRing(saveRing);447 450 this->numFacets=0; 448 451 this->ivIntPt=NULL; 449 // this->rootRing=NULL;450 //rComplete(this->baseRing);451 //rChangeCurrRing(this->baseRing);452 452 } 453 453 … … 462 462 // if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe) 463 463 // rDelete(this->rootRing); 464 // if(this->UCN!=1)465 //rDelete(this->baseRing);464 if(this->UCN!=1 && this->baseRing!=NULL) 465 rDelete(this->baseRing); 466 466 facet *fAct; 467 467 facet *fDel; … … 555 555 facet *fAct; 556 556 fAct = &f; 557 facet *codim2Act; 558 codim2Act = fAct->codim2Ptr; 559 560 cout << endl; 561 while(fAct!=NULL) 562 { 563 int64vec *fNormal; 564 fNormal=fAct->getFacetNormal(); 565 cout << "("; 566 fNormal->show(1,0); 567 if(fAct->isFlippable==TRUE) 568 cout << ") "; 569 else 570 cout << ")* "; 571 delete fNormal; 557 if(fAct!=NULL) 558 { 559 facet *codim2Act; 572 560 codim2Act = fAct->codim2Ptr; 573 cout << " Codim2: "; 574 while(codim2Act!=NULL) 575 { 576 int64vec *f2Normal; 577 f2Normal = codim2Act->getFacetNormal(); 561 562 cout << endl; 563 while(fAct!=NULL) 564 { 565 int64vec *fNormal; 566 fNormal=fAct->getFacetNormal(); 578 567 cout << "("; 579 f2Normal->show(1,0); 580 cout << ") "; 581 delete f2Normal; 582 codim2Act = codim2Act->next; 568 fNormal->show(1,0); 569 if(fAct->isFlippable==TRUE) 570 cout << ") "; 571 else 572 cout << ")* "; 573 delete fNormal; 574 codim2Act = fAct->codim2Ptr; 575 cout << " Codim2: "; 576 while(codim2Act!=NULL) 577 { 578 int64vec *f2Normal; 579 f2Normal = codim2Act->getFacetNormal(); 580 cout << "("; 581 f2Normal->show(1,0); 582 cout << ") "; 583 delete f2Normal; 584 codim2Act = codim2Act->next; 585 } 586 cout << "UCN = " << fAct->getUCN() << endl; 587 fAct = fAct->next; 583 588 } 584 cout << "UCN = " << fAct->getUCN() << endl;585 fAct = fAct->next;586 589 } 587 590 } … … 1640 1643 test|=Sy_bit(OPT_REDTAIL); 1641 1644 #ifdef gfan_DEBUG 1642 test|=Sy_bit(6); //OPT_DEBUG1645 // test|=Sy_bit(6); //OPT_DEBUG 1643 1646 #endif 1644 1647 ideal tmpI; … … 1752 1755 for (int ii=0;ii<size;ii++) 1753 1756 { 1754 poly temp1=pInit();1755 poly temp2=pInit();1757 // poly temp1=pInit(); 1758 // poly temp2=pInit(); 1756 1759 poly temp3=pInit();//polys to temporarily store values for pSub 1757 1760 // res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0)); 1758 temp1=pCopy(H->m[ii]);1761 // temp1=pCopy(H->m[ii]); 1759 1762 // temp2=pCopy(res->m[ii]); 1760 1763 //NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring 1761 temp2=pCopy(kNF(G, NULL,H->m[ii],0,0)); 1762 temp3=pSub(temp1, temp2); 1764 // temp2=pCopy(kNF(G, NULL,H->m[ii],0,0)); 1765 // temp3=pSub(temp1, temp2); 1766 temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0))); 1763 1767 res->m[ii]=pCopy(temp3); 1764 1768 //res->m[ii]=pSub(temp1,temp2); //buggy 1765 1769 //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]); 1766 pDelete(&temp1);1770 // pDelete(&temp1); 1767 1771 // pDelete(&temp2); 1768 //pDelete(&temp3); //NOTE does not work, so commented out1772 pDelete(&temp3); //NOTE does not work, so commented out 1769 1773 } 1770 1774 return res; … … 1803 1807 } 1804 1808 delete iv; 1809 delete ivNull; 1805 1810 } 1806 1811 //Remove duplicates of rows … … 1873 1878 inline int64vec *gcone::ivNeg(int64vec *iv) 1874 1879 { //Hm, switching to int64vec const int64vec does no longer work 1875 //NOTE: Can't this be done without new?1876 1880 int64vec *res;// = new int64vec(iv->length()); 1877 1881 res=iv64Copy(iv); … … 2201 2205 mpq_set(ratd,ddineq->matrix[ii][k]); 2202 2206 mpq_set_str(ddineq->matrix[ii][k],"0",10); 2203 for(int ss=k+1;ss< m;ss++)2207 for(int ss=k+1;ss<n;ss++) 2204 2208 { 2205 2209 mpq_t prod; mpq_init(prod); … … 2313 2317 gcAct->getExtremalRays(*gcAct); 2314 2318 2315 //Compute unique representation of codim-2-facets2319 //Compute unique representation of Facets and rays 2316 2320 gcAct->normalize(); 2317 2321 … … 2437 2441 { 2438 2442 gcTmp->writeConeToFile(*gcTmp); 2439 //The for loop is no longer needed2440 // for(int ii=0;ii<IDELEMS(gcTmp->gcBasis);ii++)2441 // {2442 // pDelete(&gcTmp->gcBasis->m[ii]);2443 // }2444 2443 idDelete((ideal*)&gcTmp->gcBasis);//Whonder why? 2445 2444 //If you use the following make sure it is uncommented in readConeFromFile 2446 //rDelete(gcTmp->baseRing);2445 rDelete(gcTmp->baseRing); 2447 2446 } 2448 2447 #ifdef gfan_DEBUG 2449 //if(SearchListRoot!=NULL)2450 //gcTmp->showSLA(*SearchListRoot);2451 #endif 2448 if(SearchListRoot!=NULL) 2449 gcTmp->showSLA(*SearchListRoot); 2450 #endif 2452 2451 rChangeCurrRing(gcAct->baseRing); 2453 2452 rDelete(rTmp); … … 2466 2465 //Search for cone with smallest UCN 2467 2466 gcNext = gcHead; 2468 while(gcNext!=NULL && SearchListRoot!=NULL && gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && gcNext!=(gcone * const)0xfbfbfbfb) 2467 #ifndef NDEBUG 2468 #if SIZEOF_LONG==8 //64 bit 2469 while(gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && SearchListRoot!=NULL) 2470 #elif 2471 while(gcNext!=(gcone * const)0xfbfbfbfb && SearchListRoot!=NULL) 2472 #endif 2473 #endif 2474 #ifdef NDEBUG 2475 while(gcNext!=NULL && SearchListRoot!=NULL) 2476 #endif 2469 2477 { 2470 2478 if( gcNext->getUCN() == SearchListRoot->getUCN() ) 2471 { //NOTE: Implement heuristic to be used!2479 { 2472 2480 if(gfanHeuristic==0) 2473 2481 { … … 2485 2493 gcDel = gcAct; 2486 2494 gcAct = gcNext; 2487 // rDelete(gcDel->baseRing); 2488 //Read st00f from file 2495 //Read st00f from file & 2489 2496 //implant the GB into gcAct 2490 2497 readConeFromFile(gcAct->getUCN(), gcAct); 2498 //Kill the baseRing but ONLY if it is not the ring the computation started in! 2499 if(gcDel->getUCN()!=1) 2500 rDelete(gcDel->baseRing); 2491 2501 // rAct=rCopy(gcAct->baseRing); 2492 2502 /*The ring change occurs in readConeFromFile, so as to … … 2575 2585 inline void gcone::normalize() 2576 2586 { 2577 int *ggT = new int;2578 *ggT=1;2587 // int *ggT = new int; 2588 // *ggT=1; 2579 2589 facet *fAct; 2580 2590 facet *codim2Act; … … 2585 2595 int64vec *fNormal; 2586 2596 fNormal = fAct->getFacetNormal(); 2597 int *ggT = new int; 2598 *ggT=1; 2587 2599 for(int ii=0;ii<this->numVars;ii++) 2588 2600 { 2589 2601 *ggT=intgcd((*ggT),(*fNormal)[ii]); 2590 2602 } 2591 for(int ii=0;ii<this->numVars;ii++) 2592 (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT); 2593 fAct->setFacetNormal(fNormal); 2603 if(*ggT>1)//We only need to do this if the ggT is non-trivial 2604 { 2605 for(int ii=0;ii<this->numVars;ii++) 2606 (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT); 2607 fAct->setFacetNormal(fNormal); 2608 } 2594 2609 delete fNormal; 2610 delete ggT; 2595 2611 /*And now for the codim2*/ 2596 2612 while(codim2Act!=NULL) … … 2598 2614 int64vec *n; 2599 2615 n=codim2Act->getFacetNormal(); 2616 int *ggT=new int; 2617 *ggT=1; 2600 2618 for(int ii=0;ii<this->numVars;ii++) 2601 2619 { 2602 2620 *ggT = intgcd((*ggT),(*n)[ii]); 2603 2621 } 2604 for(int ii=0;ii<this->numVars;ii++) 2605 { 2606 (*n)[ii] = ((*n)[ii])/(*ggT); 2607 } 2608 codim2Act->setFacetNormal(n); 2622 if(*ggT>1) 2623 { 2624 for(int ii=0;ii<this->numVars;ii++) 2625 { 2626 (*n)[ii] = ((*n)[ii])/(*ggT); 2627 } 2628 codim2Act->setFacetNormal(n); 2629 } 2609 2630 codim2Act = codim2Act->next; 2610 2631 delete n; 2632 delete ggT; 2611 2633 } 2612 2634 fAct = fAct->next; 2613 2635 } 2614 delete ggT;2615 2636 } 2616 2637 … … 2649 2670 deleteMarker = NULL; 2650 2671 2651 /** Flag to indicate a facet that should be added to SLA*/2652 // bool doNotAdd=FALSE;2653 2672 /** \brief Flag to mark a facet that might be added 2654 2673 * The following case may occur: … … 2673 2692 // gcone::lengthOfSearchList++; 2674 2693 } 2675 /*1st step: compare facetNormals*/ 2676 // int64vec *fNormal=NULL;gcone:: 2677 // int64vec *slNormal=NULL; 2678 2694 /*1st step: compare facetNormals*/ 2679 2695 while(fAct!=NULL) 2680 2696 { … … 2682 2698 { 2683 2699 int64vec *fNormal=NULL; 2684 // int64vec *slNormal=NULL;2685 2700 fNormal=fAct->getFacetNormal(); 2686 2701 slAct = slHead; … … 2710 2725 /*End of dumping into SLA*/ 2711 2726 while(slAct!=NULL) 2712 //while(slAct!=slEndStatic->next)2713 2727 { 2714 2728 int64vec *slNormal=NULL; … … 2716 2730 slNormal = slAct->getFacetNormal(); 2717 2731 #ifdef gfan_DEBUG 2718 cout << "Checking facet ("; 2719 fNormal->show(1,1); 2720 cout << ") against ("; 2721 slNormal->show(1,1); 2722 cout << ")" << endl; 2732 // cout << "Checking facet (";fNormal->show(1,1);cout << ") against (";slNormal->show(1,1);cout << ")" << endl; 2723 2733 #endif 2724 2734 if(areEqual(fAct,slAct)) … … 2729 2739 slHead = slAct->next; 2730 2740 if(slHead!=NULL) 2731 slHead->prev = NULL; 2741 slHead->prev = NULL; 2732 2742 } 2733 2743 else if (slAct==slEnd) 2734 2744 { 2735 2745 slEnd=slEnd->prev; 2736 slEnd->next = NULL; 2746 slEnd->next = NULL; 2737 2747 } 2738 2748 else 2739 2749 { 2740 2750 slAct->prev->next = slAct->next; 2741 slAct->next->prev = slAct->prev; 2751 slAct->next->prev = slAct->prev; 2742 2752 } 2743 2753 removalOccured=TRUE; 2744 // cout << gcone::lengthOfSearchList << endl;2745 2754 gcone::lengthOfSearchList--; 2746 // cout << gcone::lengthOfSearchList << endl;2747 2755 if(deleteMarker!=NULL) 2748 2756 { … … 2751 2759 } 2752 2760 #ifdef gfan_DEBUG 2753 cout << "Removing ("; 2754 fNormal->show(1,1); 2755 cout << ") from list" << endl; 2761 cout << "Removing (";fNormal->show(1,1);cout << ") from list" << endl; 2756 2762 #endif 2757 2763 delete slNormal; 2758 2764 break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct 2759 } 2760 2765 } 2761 2766 slAct = slAct->next; 2762 2767 /* NOTE The following lines must not be here but rather called inside the if-block above. … … 2776 2781 { 2777 2782 #ifdef gfan_DEBUG 2778 cout << "Adding facet ("; 2779 fNormal->show(1,0); 2780 cout << ") to SLA " << endl; 2783 cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl; 2781 2784 #endif 2782 2785 //Add fAct to SLA … … 2787 2790 2788 2791 slEnd->next = new facet(); 2789 slEnd = slEnd->next; 2792 slEnd = slEnd->next;//Update slEnd 2790 2793 facet *slEndCodim2Root; 2791 2794 facet *slEndCodim2Act; … … 2820 2823 } 2821 2824 gcone::lengthOfSearchList++; 2822 // cout << gcone::lengthOfSearchList << endl;2823 2825 }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) || 2824 2826 fAct = fAct->next; … … 3019 3021 ss << UCN; 3020 3022 string UCNstr = ss.str(); 3021 // string line;3022 // string strGcBasisLength;3023 // string strMonom, strCoeff, strCoeffNom, strCoeffDenom;3024 3023 int gcBasisLength=0; 3025 // int intCoeff=1;3026 // int intCoeffNom=1; //better (number) but that's not compatible with stringstream;3027 // int intCoeffDenom=1;3028 3029 // bool hasCoeffInQ = FALSE; //does the polynomial have rational coeff?3030 // bool hasNegCoeff = FALSE; //or a negative one?3031 3024 size_t found; //used for find_first_(not)_of 3032 3025 … … 3042 3035 ring saveRing=currRing; 3043 3036 //Comment the following if you uncomment the if(line=="RING") part below 3044 rChangeCurrRing(gc->baseRing);3037 // rChangeCurrRing(gc->baseRing); 3045 3038 3046 3039 while( !gcInputFile.eof() ) … … 3048 3041 string line; 3049 3042 getline(gcInputFile,line); 3050 // hasCoeffInQ = FALSE;3051 // hasNegCoeff = FALSE;3052 3053 3043 if(line=="RING") 3054 3044 { 3055 //getline(gcInputFile,line);3056 //found = line.find("a(");3057 //line.erase(0,found+2);3058 //string strweight;3059 //strweight=line.substr(0,line.find_first_of(")"));3060 //int64vec *iv=new int64vec(this->numVars);3061 //for(int ii=0;ii<this->numVars;ii++)3062 //{3063 //string weight;3064 //weight=line.substr(0,line.find_first_of(",)"));3065 // (*iv)[ii]=atoi(weight.c_str()); 3066 //line.erase(0,line.find_first_of(",)")+1);3067 //}3068 //ring newRing;3069 //if(currRing->order[0]!=ringorder_a)3070 //{3071 //newRing=rCopyAndAddWeight(currRing,iv);3072 //}3073 //else3074 //{3075 //newRing=rCopy0(currRing);3076 //int length=this->numVars;3077 //int *A=(int *)omAlloc0(length*sizeof(int));3078 //for(int jj=0;jj<length;jj++)3079 //{3080 // A[jj]=-(*iv)[jj];3081 //}3082 //omFree(newRing->wvhdl[0]);3083 //newRing->wvhdl[0]=(int*)A;3084 //newRing->block1[0]=length;3085 //}3086 //delete iv;3087 //rComplete(newRing);3088 //gc->baseRing=rCopy(newRing);3089 // //rDelete(newRing);3090 //rComplete(gc->baseRing);3091 //if(currRing!=gc->baseRing)3092 //rChangeCurrRing(gc->baseRing);3045 getline(gcInputFile,line); 3046 found = line.find("a("); 3047 line.erase(0,found+2); 3048 string strweight; 3049 strweight=line.substr(0,line.find_first_of(")")); 3050 int64vec *iv=new int64vec(this->numVars); 3051 for(int ii=0;ii<this->numVars;ii++) 3052 { 3053 string weight; 3054 weight=line.substr(0,line.find_first_of(",)")); 3055 (*iv)[ii]=atol(weight.c_str());//Better to long. Weight bound in Singular:2147483647 3056 line.erase(0,line.find_first_of(",)")+1); 3057 } 3058 ring newRing; 3059 if(currRing->order[0]!=ringorder_a) 3060 { 3061 newRing=rCopyAndAddWeight(currRing,iv); 3062 } 3063 else 3064 { 3065 newRing=rCopy0(currRing); 3066 int length=this->numVars; 3067 int *A=(int *)omAlloc0(length*sizeof(int)); 3068 for(int jj=0;jj<length;jj++) 3069 { 3070 A[jj]=(*iv)[jj]; 3071 } 3072 omFree(newRing->wvhdl[0]); 3073 newRing->wvhdl[0]=(int*)A; 3074 newRing->block1[0]=length; 3075 } 3076 delete iv; 3077 rComplete(newRing); 3078 gc->baseRing=rCopy(newRing); 3079 rDelete(newRing); 3080 rComplete(gc->baseRing); 3081 if(currRing!=gc->baseRing) 3082 rChangeCurrRing(gc->baseRing); 3093 3083 } 3094 3084 … … 3115 3105 while(!line.empty()) 3116 3106 { 3117 // poly strPoly =pInit();3107 // poly strPoly;//=pInit(); 3118 3108 number nCoeff=nInit(1); 3119 3109 number nCoeffNom=nInit(1); … … 3224 3214 /** \brief Gather the output 3225 3215 * List of lists 3216 * If heuristic==1 readConeFromFile() is called once more on every cone. This may slow down the computation but it also 3217 * allows us to rDelete(gcDel->baseRing) and the such in gcone::noRevS. 3226 3218 *\param *gc Pointer to gcone, preferably gcRoot ;-) 3227 3219 *\param n the number of cones … … 3406 3398 gcAct->getConeNormals(gcAct->gcBasis); 3407 3399 gcAct->noRevS(*gcAct); //Here we go! 3400 //Switch back to the ring the computation was started in 3408 3401 rChangeCurrRing(inputRing); 3409 3402 //res=gcAct->gcBasis; … … 3436 3429 /*Return result*/ 3437 3430 #ifdef gfanp 3438 cout << "t_getConeNormals:" << gcone::time_getConeNormals << endl;3431 cout << endl << "t_getConeNormals:" << gcone::time_getConeNormals << endl; 3439 3432 // cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl; 3440 3433 // cout << " t_ddMC:" << gcone::t_ddMC << endl;
Note: See TracChangeset
for help on using the changeset viewer.