Changeset 0ad81f in git
 Timestamp:
 Apr 7, 2009, 11:44:20 AM (14 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
 Children:
 2ae96e40fc5453bcb155aec76d376d79dd549cbe
 Parents:
 f32177e9dea6c9833a7f2f8986a2cf6273f30fab
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

kernel/gfan.cc
rf32177 r0ad81f 2 2 Compute the Groebner fan of an ideal 3 3 $Author: monerjan $ 4 $Date: 2009040 6 14:57:18$5 $Header: /exports/cvsroot2/cvsroot/kernel/gfan.cc,v 1. 29 20090406 14:57:18monerjan Exp $6 $Id: gfan.cc,v 1. 29 20090406 14:57:18monerjan Exp $4 $Date: 20090407 09:44:20 $ 5 $Header: /exports/cvsroot2/cvsroot/kernel/gfan.cc,v 1.30 20090407 09:44:20 monerjan Exp $ 6 $Id: gfan.cc,v 1.30 20090407 09:44:20 monerjan Exp $ 7 7 */ 8 8 … … 267 267 cout << "Cols = " << ddcols << endl; 268 268 #endif 269 269 270 /*Write the normals into class facet*/ 270 271 #ifdef gfan_DEBUG … … 275 276 this>facetPtr = fRoot; //set variable facetPtr of class gcone to first facet 276 277 facet *fAct; //instantiate pointer to active facet 277 fAct = fRoot; //This does not seemto do the trick. fRoot and fAct have to point to the same adress!278 fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress! 278 279 std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl; 279 280 for (int kk = 0; kk<ddrows; kk++) … … 291 292 if (jj<ddcols) //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :) 292 293 { 293 fAct>next = new facet(); //If so: instantiate new facet. Otherwise this>next=NULL due to the constructor294 //fAct>next = new facet(); //If so: instantiate new facet. Otherwise this>next=NULL due to the constructor 294 295 } 295 }//for jj<ddcols 296 /*Now load should be full and we can call setFacetNormal*/ 297 fAct>setFacetNormal(load); 298 //fAct>printNormal(); 299 fAct=fAct>next; //this should definitely not be called in the above whileloop :D 296 }//for (int jj = 1; jj <ddcols; jj++) 297 /*Quick'n'dirty hack for flippability*/ 298 bool isFlippable; 299 for (int jj = 0; jj<this>numVars; jj++) 300 { 301 intvec *ivCanonical = new intvec(this>numVars); 302 (*ivCanonical)[jj]=1; 303 if (dotProduct(load,ivCanonical)>=0) 304 { 305 isFlippable=FALSE; 306 } 307 else 308 { 309 isFlippable=TRUE; 310 } 311 delete ivCanonical; 312 }//for (int jj = 0; jj<this>numVars; jj++) 313 if (isFlippable==FALSE) 314 { 315 cout << "Ignoring facet"; 316 load>show(); 317 //fAct>next=NULL; 318 } 319 else 320 { /*Now load should be full and we can call setFacetNormal*/ 321 fAct>setFacetNormal(load); 322 fAct>next = new facet(); 323 //fAct>printNormal(); 324 fAct=fAct>next; //this should definitely not be called in the above whileloop :D 325 }//if (isFlippable==FALSE) 300 326 delete load; 301 } 327 }//for (int kk = 0; kk<ddrows; kk++) 302 328 /* 303 329 Now we should have a linked list containing the facet normals of those facets that are … … 316 342 }//method getConeNormals(ideal I) 317 343 318 /*bool isParallel(int v[], intvec iv)319 {320 }*/321 344 322 345 /** \brief Compute the Groebner Basis on the other side of a shared facet … … 335 358 void flip(ideal gb, facet *f) //Compute "the other side" 336 359 { 337 338 360 intvec *fNormal = new intvec(this>numVars); //facet normal, check for parallelity 339 361 fNormal = f>getFacetNormal(); //read this>fNormal; 340 362 #ifdef gfan_DEBUG 363 cout << "===" << endl; 364 cout << "running gcone::flip" << endl; 341 365 cout << "fNormal="; 342 366 fNormal>show(); … … 384 408 cout << "" << endl; 385 409 /*Now initialFormElement must be added to (ideal)initialForm */ 386 //f>flibGB>m[ii]=(poly)initialFormElement[ii];387 //(poly)initialForm>m[ii]=pAdd(initialForm[ii],initialFormElement[ii]);388 410 initialForm>m[ii]=initialFormElement[ii]; 389 411 }//for … … 391 413 #ifdef gfan_DEBUG 392 414 cout << "Initial ideal is: " << endl; 393 //idShow(initialForm);415 idShow(initialForm); 394 416 f>printFlipGB(); 395 417 cout << "===" << endl; … … 397 419 delete check; 398 420 399 /*2nd step: lift initial ideal to a GB of the neighbouring cone*/ 400 ring tmpring=rCopy0(currRing); 401 tmpring>order[0]=ringorder_a; 402 tmpring>order[1]=ringorder_dp; 403 tmpring>order[2]=ringorder_C; 404 //rWrite(tmpring); 405 tmpring>wvhdl[0] =( int *)omAlloc((fNormal>length())*sizeof(int)); //found in Singular/ipshell.cc 421 /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/ 422 /*Substep 2.1 423 compute $G_{\alpha}(in_v(I)) 424 see journal p. 66 425 */ 426 ring srcRing=currRing; 427 ring dstRing=rCopy0(srcRing); 428 dstRing>order[0]=ringorder_a; 429 //tmpring>order[1]=ringorder_dp; 430 //tmpring>order[2]=ringorder_C; 431 dstRing>wvhdl[0] =( int *)omAlloc((fNormal>length())*sizeof(int)); //found in Singular/ipshell.cc 406 432 407 433 for (int ii=0;ii<this>numVars;ii++) 408 { 409 cout << "ping" << endl; 410 tmpring>wvhdl[0][ii]=(*fNormal)[ii]; //What exactly am I doing here? 411 cout << "pong" << endl; 412 cout << *tmpring>wvhdl[ii] << endl; 434 { 435 dstRing>wvhdl[0][ii]=(*fNormal)[ii]; //What exactly am I doing here? 436 //cout << tmpring>wvhdl[0][ii] << endl; 413 437 } 414 //tmpring>wvhdl=(int**)(fNormal); 415 rComplete(tmpring); 416 rWrite(tmpring); 417 /*setring(tmpring); 418 ideal ina=initialForm; 438 rComplete(dstRing); 439 rChangeCurrRing(dstRing); 440 map theMap=(map)idMaxIdeal(1); 441 rWrite(currRing); cout << endl; 442 ideal ina; 443 ina=fast_map(initialForm,srcRing,(ideal)theMap,dstRing); 444 cout << "ina="; 445 idShow(ina); cout << endl; 419 446 ideal H; 420 H=kstd(ina,NULL,testHomog,NULL); 421 idShow(H);*/ 447 H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous 448 idSkipZeroes(H); 449 cout << "H="; idShow(H); cout << endl; 450 /*Substep 2.2 451 do the lifting 452 */ 453 rChangeCurrRing(srcRing); 454 ideal srcRing_H; 455 ideal srcRing_HH; 456 //map theMap = (map)idMaxIdeal(1); 457 srcRing_H=fast_map(H,dstRing,(ideal)theMap,srcRing); 458 srcRing_HH=srcRing_Hstdred(srcRing_H,this>gcBasis); 459 /*Substep 2.3 460 turn the minimal basis into a reduced one 461 */ 422 462 423 463 }//void flip(ideal gb, facet *f) … … 429 469 *\return void 430 470 */ 431 void getGB( const ideal&inputIdeal)471 void getGB(ideal const &inputIdeal) 432 472 { 433 473 ideal gb; … … 435 475 idSkipZeroes(gb); 436 476 this>gcBasis=gb; //write the GB into gcBasis 437 } 477 }//void getGB 438 478 439 479 ideal GenGrbWlk(ideal, ideal); //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets … … 442 482 /** \brief Compute the dot product of two intvecs 443 483 * 444 * THIS IS WEIRD  BUT WORKS445 484 */ 446 int dotProduct(intvec **a, intvec **b)447 { 448 intvec iva=*a;449 intvec ivb=*b;485 int dotProduct(intvec const &iva, intvec const &ivb) 486 { 487 //intvec iva=a; 488 //intvec ivb=b; 450 489 int res=0; 451 490 for (int i=0;i<this>numVars;i++) … … 454 493 } 455 494 return res; 456 } 495 }//int dotProduct 457 496 458 497 /** \brief Check whether two intvecs are parallel … … 460 499 * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$ 461 500 */ 462 bool isParallel(intvec *a, intvec *b)501 bool isParallel(intvec const &a, intvec const &b) 463 502 { 464 503 int lhs,rhs; 465 lhs=dotProduct( &a,&b)*dotProduct(&a,&b);466 rhs=dotProduct( &a,&a)*dotProduct(&b,&b);504 lhs=dotProduct(a,b)*dotProduct(a,b); 505 rhs=dotProduct(a,a)*dotProduct(b,b); 467 506 cout << "LHS="<<lhs<<", RHS="<<rhs<<endl; 468 507 if (lhs==rhs) … … 475 514 } 476 515 }//bool isParallel 477 478 /** \brief Check two intvecs for equality  PROBABLY NOT NEEDED479 *480 * \f$ \alpha=\beta\Leftrightarrow\langle\alpha\beta,\alpha\beta\rangle=0 \f$481 */482 bool isEqual(intvec a, intvec b)483 {484 intvec *ivdiff;485 int res;486 487 ivdiff=ivSub(&a,&b);488 cout << "ivdiff=";489 ivdiff>show();490 cout << endl;491 //res=dotProduct(ivdiff,ivdiff);492 if (res==0)493 {494 return TRUE;495 }496 else497 {498 return FALSE;499 }500 }//bool isEqual501 516 };//class gcone 502 517 … … 566 581 => Count the cones! 567 582 */ 568 583 rChangeCurrRing(inputRing); 569 584 res=gcAct>gcBasis; 570 585 //cout << fRoot << endl;
Note: See TracChangeset
for help on using the changeset viewer.