Annotation of OpenXM_contrib/TiGERS_0.9/gset.c, Revision 1.1
1.1 ! maekawa 1: /*
! 2: ** gset.c Birk Huber, 4/99
! 3: ** -- definitions and basic operations on gsets
! 4: **
! 5: **
! 6: ** TiGERS, Toric Groebner Basis Enumeration by Reverse Search
! 7: ** copyright (c) 1999 Birk Huber
! 8: **
! 9: */
! 10: #include<stdio.h>
! 11: #include<stdlib.h>
! 12: #include "utils.h"
! 13: #include "gset.h"
! 14: #include "matrices.h"
! 15: #include "Rsimp.h"
! 16:
! 17: /* uncomment the following to cause use sainity checks (not fully tested)*/
! 18: /*
! 19: #define GTEST 0
! 20: #include "gtesting.h"
! 21: */
! 22:
! 23: /*
! 24: ** create an empty gset, initialize all its pointers to null.
! 25: */
! 26: gset gset_new(){
! 27: gset g=0;
! 28: if ((g=(gset)malloc(sizeof(struct gset_tag)))==0){
! 29: fprintf(stderr,"gset_new(): malloc failed\n");
! 30: return 0;
! 31: }
! 32: gset_first(g)=0;
! 33: gset_cache_vtx(g)=0;
! 34: gset_cache_edge(g)=0;
! 35: return g;
! 36: }
! 37:
! 38: /*
! 39: ** delete a gset and all binomials in it.
! 40: */
! 41: void gset_free(gset g){
! 42: binomial b;
! 43: while((b=gset_first(g))!=0){
! 44: gset_first(g)=binomial_next(b);
! 45: binomial_free(b);
! 46: }
! 47: free((void *)g);
! 48: }
! 49:
! 50:
! 51: /*
! 52: ** read a gset from an input stream
! 53: ** (return TRUE if gset read in )
! 54: */
! 55: int gset_read(FILE *is,gset g){
! 56: char c;
! 57: binomial b;
! 58:
! 59: eatwhite(is);
! 60: c=getc(is);
! 61: if (c!='{') return FALSE;
! 62: c=',';
! 63: while(c==','){
! 64: b=binomial_new();
! 65: binomial_read(is,b);
! 66: gset_insert(g,b);
! 67: eatwhite(is);
! 68: c=getc(is);
! 69: }
! 70: if (c!='}'){
! 71: fprintf(stderr,"gset_read(): gset should end with a '}'\n");
! 72: return FALSE;
! 73: }
! 74: ungetc(c,is);
! 75: return TRUE;
! 76: }
! 77:
! 78: /*
! 79: ** Display a gset
! 80: */
! 81: void gset_print(FILE *of,gset g){
! 82: binomial b;
! 83:
! 84: fprintf(of,"{");
! 85: for(b=g->bottom;b!=0;b=b->next){
! 86: binomial_print(of,b);
! 87: if (b->next!=0) fprintf(of,", ");
! 88: }
! 89: fprintf(of,"}");
! 90: }
! 91:
! 92: /*
! 93: ** Display a gsets initial ideal
! 94: */
! 95: void gset_init_print(FILE *of,gset g){
! 96: binomial b;
! 97: fprintf(of,"{");
! 98: for(b=g->bottom;b!=0;b=b->next){
! 99: print_monomial(of,binomial_lead(b));
! 100: if (b->next!=0) fprintf(of,", ");
! 101: }
! 102: fprintf(of,"}");
! 103: }
! 104:
! 105: /*
! 106: ** Display a facet binomials
! 107: */
! 108: void gset_facet_print(FILE *of,gset g){
! 109: binomial b;
! 110: int tog=0;
! 111: fprintf(of,"{");
! 112: for(b=g->bottom;b!=0;b=b->next){
! 113: if (b->ff==FACET){
! 114: if (tog==1) fprintf(of,",");
! 115: else tog=1;
! 116: binomial_print(of,b);
! 117: }
! 118: }
! 119: fprintf(of,"}");
! 120: }
! 121:
! 122: /*
! 123: ** place binomial in gset in decreasing lexicographic order.
! 124: ** (using binomial_compair(,) to determine reletive order or
! 125: ** binomials)
! 126: */
! 127: void gset_insert(gset g, binomial b){
! 128: binomial t1,t2;
! 129: if (gset_first(g)==0 || binomial_compair(gset_first(g),b)==FALSE){
! 130: binomial_next(b)=gset_first(g);
! 131: gset_first(g)=b;
! 132: }
! 133: else{
! 134: t1=gset_first(g);
! 135: while(t1!=0){
! 136: t2=binomial_next(t1);
! 137: if (t2==0 || binomial_compair(t2,b)==FALSE){
! 138: binomial_next(t1)=b;
! 139: binomial_next(b)=t2;
! 140: break;
! 141: }
! 142: t1=t2;
! 143: }
! 144: }
! 145: }
! 146:
! 147: /*
! 148: ** return first facet binomial which is mis-marked
! 149: */
! 150: binomial gset_downedge(gset g){
! 151: binomial ptr=0;
! 152: for(ptr=gset_first(g);ptr!=0;ptr=binomial_next(ptr)){
! 153: if (binomial_ordered(ptr)==FALSE && gset_isfacet(g,ptr)==TRUE)break;
! 154: }
! 155: return ptr;
! 156: }
! 157:
! 158:
! 159: /*
! 160: ++void bmprint(FILE *of,binomial B, binomial M):
! 161: ++ simple output function for a generating set consisting of a binomial and
! 162: ++ a list of monomials --- usefull mainly for debugging. -- gset output
! 163: ++ routing may be modified in such a way as to use code which could double
! 164: ++ for this purpose.
! 165: ++
! 166: */
! 167: void bmprint(FILE *of,binomial B, binomial M){
! 168: fprintf(of,"{ ");
! 169: if (B!=0) binomial_print(of,B);
! 170: while(M!=0){
! 171: fprintf(of,",");
! 172: binomial_print(of,M);
! 173: M=binomial_next(M);
! 174: }
! 175: fprintf(of,"}\n");
! 176: }
! 177:
! 178: /*
! 179: **binomial bmleadreduce(binomial m1, binomial B, binomial mlist):
! 180: ** "lift" monomial m1 to a binomial using bmlist B,mlist.
! 181: */
! 182: binomial bmleadreduce(binomial m1, binomial B, binomial mlist){
! 183: binomial res,ptr;
! 184: /* case 1: the binomial is used in the decomp:*/
! 185: if (monomial_divides(binomial_lead(B),binomial_lead(m1))==TRUE){
! 186: res=binomial_new();
! 187: binomial_copy(B,res);
! 188: binomial_bumpto(res,m1);
! 189:
! 190: /* now reduce by lt(b) as many times as possible */
! 191: while(monomial_divides(binomial_lead(B),binomial_trail(res))==TRUE){
! 192: reducetrail(res,B);
! 193: }
! 194:
! 195: /* now find a monomial deviding whats left and reduce by it*/
! 196: for(ptr=mlist;ptr!=0;ptr=ptr->next){
! 197: if (monomial_divides(binomial_lead(ptr),binomial_trail(res))==TRUE){
! 198: reducetrail(res,ptr);
! 199: return res;
! 200: }
! 201: }
! 202: }
! 203: else { /* case 2 -- our monomial is a multiple of some monomial*/
! 204: for(ptr=mlist;ptr!=0;ptr=ptr->next){
! 205: if (monomial_divides(binomial_lead(ptr),binomial_lead(m1))==TRUE){
! 206: res=binomial_new();
! 207: binomial_copy(ptr,res);
! 208: binomial_set(res);
! 209: binomial_bumpto(res,m1);
! 210: return res;
! 211: }
! 212: }
! 213: }
! 214: fprintf(stderr, "Warning BMreduce does not get to zero\n");
! 215: }
! 216:
! 217:
! 218:
! 219: /*
! 220: ++ binomial find_divisor(binomial S,binomial L):
! 221: ++ search list starting at L for a monomial dividing S
! 222: ++ return this binomial or zero if no such binomial exists on list.
! 223: ++
! 224: */
! 225: binomial find_divisor(binomial S,binomial L){
! 226: binomial ptr;
! 227: for(ptr=L;ptr!=0;ptr=ptr->next){
! 228: if(ptr!=S&&monomial_divides(binomial_lead(ptr),binomial_lead(S))==TRUE){
! 229: break;
! 230: }
! 231: }
! 232: return ptr;
! 233: }
! 234:
! 235:
! 236: /*
! 237: ++ void remove_multiples(binomial S, binomial *L):
! 238: ++ Pass through list starting at *L, removing all multiples of *S
! 239: ++ encountered (except S if it is on the list).
! 240: ++
! 241: */
! 242: void remove_multiples(binomial S, binomial *L){
! 243: binomial tlist=*L,ptr;
! 244: *L=0;
! 245: while(tlist!=0){
! 246: ptr=tlist; tlist=tlist->next; ptr->next=0; /* pop tlist*/
! 247: if (S==ptr || monomial_divides(binomial_lead(S),binomial_lead(ptr))==FALSE) {
! 248: ptr->next=*L; /* push L */
! 249: *L=ptr;
! 250: }
! 251: else binomial_free(ptr);
! 252: }
! 253: }
! 254:
! 255: /*
! 256: ** void bmrgb(binomial B, binomial *ms):
! 257: ** transform a marked binomial and a list of monomials into the
! 258: ** reduced grobner basis for the ideal the generate.
! 259: **
! 260: */
! 261: void bmrgb(binomial B, binomial *ms){
! 262: int i,jk;
! 263: binomial DL,SL,tlist,S,tmp,ptr;
! 264: int szero;
! 265:
! 266: DL=*ms;
! 267: SL=0;
! 268:
! 269: while(DL!=0){
! 270: tmp=DL; DL=DL->next; /* pop DL */
! 271: tmp->next=SL; SL=tmp; /* push SL */
! 272:
! 273: /* calculate normal form of S pair with respect to {B}*/
! 274: S=monomial_spair(B,binomial_lead(tmp));
! 275:
! 276: /* check if anything on DL or on SL divides S*/
! 277: if (find_divisor(S,SL)!=0||find_divisor(S,DL)!=0)binomial_free(S);
! 278: else {
! 279: remove_multiples(S,&SL); /* remove any multiples of S from SL*/
! 280: remove_multiples(S,&DL); /* remove any multiples of S from DL*/
! 281: S->next=DL; DL=S; /* finally push S onto DL */
! 282: }
! 283: }
! 284:
! 285: remove_multiples(B,&SL); /* now that we have a grobner basis
! 286: we can remove multiples of lead(B)
! 287: from the monomial list */
! 288:
! 289: // replace list bottom pointer
! 290: *ms=SL;
! 291:
! 292: }
! 293:
! 294: /*
! 295: ** gset gset_flip(gset g1, binomial b):
! 296: ** Take a gset and one of its facet binomials and produce the new gset
! 297: ** corresponding to the grobner cone which shares the given facet binomial.
! 298: **
! 299: */
! 300: gset gset_flip(gset g1, binomial b){
! 301: binomial old=0,new=0,tmp,ptr,B,Bold;
! 302: gset gres=0, gtest=0;
! 303:
! 304: /* sainity check -- remove eventually*/
! 305: #ifdef GTEST
! 306: if (lp_isfacet(g1,b)!=FACET){
! 307: fprintf(stdout,"WARNING flip asked to flip non-binomial\n");
! 308: }
! 309: #endif
! 310:
! 311: /* copy all binomials of g's list ecept b to old and new lists*/
! 312: for(ptr=gset_first(g1);ptr!=0;ptr=ptr->next){
! 313: if (ptr!=b){
! 314: tmp=binomial_new();
! 315: binomial_copy(ptr,tmp);
! 316: monomial_set(tmp);
! 317: binomial_next(tmp)=new;
! 318: new=tmp;
! 319: tmp=binomial_new();
! 320: binomial_copy(ptr,tmp);
! 321: monomial_set(tmp);
! 322: binomial_next(tmp)=old;
! 323: old=tmp;
! 324: }
! 325: }
! 326: Bold=binomial_new();
! 327: binomial_copy(b,Bold);
! 328:
! 329: /* make fliped copy of b*/
! 330: B=binomial_new();
! 331: binomial_copy(b,B);
! 332: binomial_flip(B);
! 333:
! 334: /* calculate the reduced GB of the binomial b and the monomials in M */
! 335: bmrgb(B,&new);
! 336:
! 337: /* create new gset and insert flipped copy of b */
! 338: gres=gset_new();
! 339: gset_insert(gres,B);
! 340:
! 341: /* do lifting of elements in mlist and add them to new gset*/
! 342: while((tmp=new)!=0){
! 343: new=new->next;
! 344: ptr=bmleadreduce(tmp,Bold,old);
! 345: gset_insert(gres,ptr);
! 346: binomial_free(tmp);
! 347: }
! 348:
! 349: /* sainity check -- remove eventually*/
! 350: #ifdef GTEST
! 351: if (gset_is_gb(gres)==FALSE){
! 352: fprintf(stderr,"flipped failed to find grobner basis \n");
! 353: fprintf(stderr,"g1=");gset_print(stderr,g1);fprintf(stderr,"\n");
! 354: fprintf(stderr,"b=");binomial_print(stderr,b);fprintf(stderr,"\n");
! 355: exit(1);
! 356: }
! 357: #endif
! 358:
! 359: binomial_free(Bold);
! 360: while(old!=0){
! 361: ptr=old;
! 362: old=old->next;
! 363: binomial_free(ptr);
! 364: }
! 365:
! 366: gset_autoreduce(gres);
! 367:
! 368: /* sainity check -- remove eventually */
! 369: #ifdef GTEST
! 370: if ( gset_is_reduced(gres)==FALSE){
! 371: fprintf(stderr,"flipped failed to find reduced grobner basis \n");
! 372: fprintf(stderr,"g1=");gset_print(stderr,g1);fprintf(stderr,"\n");
! 373: fprintf(stderr,"b=");binomial_print(stderr,b);fprintf(stderr,"\n");
! 374: exit(1);
! 375: }
! 376: #endif
! 377:
! 378: return gres;
! 379: }
! 380:
! 381:
! 382: /*
! 383: ** void gset_autoreduce(gset g):
! 384: ** autoreduce g according to its marking
! 385: */
! 386: void gset_autoreduce(gset g){
! 387: binomial b1,b2;
! 388: int changed;
! 389:
! 390: /* First remove redundant elements from the grobner basis*/
! 391: while(find_divisor(gset_first(g),gset_first(g))!=0){
! 392: gset_first(g)=(b1=gset_first(g))->next;
! 393: binomial_free(b1);
! 394: }
! 395: b1=gset_first(g);
! 396: while(b1!=0 && binomial_next(b1)!=0){
! 397: if (find_divisor(binomial_next(b1),gset_first(g))!=0) {
! 398: b2=binomial_next(b1);
! 399: binomial_next(b1)=binomial_next(b2);
! 400: binomial_free(b2);
! 401: }
! 402: else b1=binomial_next(b1);
! 403: }
! 404:
! 405: /* now reduce trailing terms */
! 406: b1=gset_first(g);
! 407: while(b1!=0){
! 408: changed=1;
! 409: while(changed==1){
! 410: changed=0;
! 411: b2=gset_first(g);
! 412: while(b2!=0){
! 413: if (b2!=b1 &&
! 414: monomial_divides(binomial_lead(b2),binomial_trail(b1))==TRUE){
! 415: reducetrail(b1,b2);
! 416: changed=1;
! 417: }
! 418: b2=binomial_next(b2);
! 419: }
! 420: }
! 421: b1=binomial_next(b1);
! 422: }
! 423:
! 424: /* now reinsert terms to ensure list is in order */
! 425: b1=gset_first(g);
! 426: gset_first(g)=0;
! 427: while((b2=b1)!=0) {
! 428: b1=binomial_next(b1);
! 429: gset_insert(g,b2);
! 430: }
! 431:
! 432: }
! 433:
! 434: /*
! 435: ** void gset_rgb(gset g, int (*comp)(monomial,monomial))
! 436: ** Use buchberger's algorithm to transform g into groebner basis
! 437: ** wrt the term order implemented by comp
! 438: */
! 439: #define ordered(b) ((comp(binomial_lead(b),binomial_trail(b))>=0)?TRUE:FALSE)
! 440: int changed;
! 441: void gset_rgb(gset g, int (*comp)(monomial,monomial)){
! 442: void process_pair(binomial,binomial,binomial,binomial,binomial *,int (*comp)());
! 443: binomial ptr,b,DL=0,SL=0,NL=0;
! 444: int passes=0;
! 445:
! 446: do{
! 447: DL=0; SL=0; NL=0;
! 448: changed=FALSE;
! 449:
! 450: /* reorient g and put binomial list on DL*/
! 451: ptr=gset_first(g);
! 452: while((b=ptr)!=0){
! 453: ptr=binomial_next(ptr);
! 454: if (ordered(b)==FALSE) binomial_flip(b);
! 455: binomial_next(b)=DL;
! 456: DL=b;
! 457: }
! 458:
! 459: /* */
! 460: while(DL!=0){
! 461: NL=0;
! 462: b=DL; DL=DL->next; /* pop DL*/
! 463: b->next=SL; SL=b; /* push ptr onto SL*/
! 464: for(ptr=SL;ptr!=0;ptr=binomial_next(ptr)) process_pair(b,ptr,DL,SL,&NL,comp);
! 465: for(ptr=DL;ptr!=0;ptr=binomial_next(ptr)){
! 466: process_pair(b,ptr,DL,SL,&NL,comp);
! 467: if (binomial_next(ptr)==0){ binomial_next(ptr)=NL; NL=0;}
! 468: }
! 469: }
! 470:
! 471: /* reinsert elements in SL into gset */
! 472: ptr=SL;
! 473: gset_first(g)=0;
! 474: while(ptr!=0){
! 475: ptr=(b=ptr)->next;
! 476: gset_insert(g,b);
! 477: }
! 478:
! 479:
! 480: gset_autoreduce(g);
! 481:
! 482: if (passes>1){
! 483: if (changed==TRUE ){fprintf(stderr,"looping pass (%d)\n",passes);}
! 484: else {fprintf(stderr,"stopping after %d passes\n",passes);}
! 485: }
! 486: passes++;
! 487: }while(changed==TRUE);
! 488:
! 489: }
! 490:
! 491: /*
! 492: -- process_pair(b1,b2,DL,SL,NL,comp):
! 493: -- process spair defined by binomials b1 and b2.
! 494: -- taking spair wrt union of DL SL and *NL lists
! 495: -- if non-zero add result to NL list.
! 496: -- comp is a monomial compairison function used throuout.
! 497: */
! 498: void process_pair(binomial b1,binomial b2,binomial DL,binomial SL,binomial *NL,
! 499: int (*comp)(monomial,monomial)){
! 500: binomial S,ptr,find_divisor(binomial,binomial);
! 501: int i,tmp;
! 502:
! 503: /* test is lead terms are rel prime [buchberger's first criterion]*/
! 504: if (monomial_rel_prime(binomial_lead(b1),binomial_lead(b2))==TRUE)return;
! 505:
! 506: /* buchberger's second criterion ?*/
! 507:
! 508:
! 509: S=binomial_new();
! 510: binomial_copy(b1,S);
! 511: if (binomial_spair(S,b2)==TRUE){
! 512: binomial_free(S);
! 513: return;
! 514: }
! 515: if (ordered(S)==FALSE) binomial_flip(S);
! 516: while(((ptr=find_divisor(S,DL))!=0)||
! 517: ((ptr=find_divisor(S,SL))!=0)||
! 518: ((ptr=find_divisor(S,*NL))!=0)){
! 519: if (binomial_spair(S,ptr)==TRUE){
! 520: binomial_free(S);
! 521: return;
! 522: }
! 523: if (ordered(S)==FALSE) binomial_flip(S);
! 524: }
! 525: changed=TRUE;
! 526: binomial_next(S)=*NL;
! 527: *NL=S;
! 528: }
! 529:
! 530: /*
! 531: ** void gset_setfacets(gset g):
! 532: ** Use repeated calls to gset_isfacet to ensure that all facet flags
! 533: ** are set.
! 534: ** Also sets nfacets, nelts and deg counters (these contain garbage except
! 535: ** right after a call to set facets).
! 536: **
! 537: ** NOTE 4/20 the deg counter was just added and the code to load it is
! 538: ** cludged in here -- logically it should be part of the
! 539: ** binomial.c file -- for now it is here.
! 540: */
! 541: #define max(a,b) (((a)>(b))?(a):(b))
! 542: void gset_setfacets(gset g){
! 543: binomial b;
! 544: monomial ml,mt;
! 545: int i,d=0,dl,dt;
! 546: gset_nfacets(g)=0;
! 547: gset_nelts(g)=0;
! 548: for(b=gset_first(g);b!=0;b=binomial_next(b)){
! 549: if (gset_isfacet(g,b)==FACET) gset_nfacets(g)++;
! 550: gset_nelts(g)++;
! 551: ml=binomial_lead(b); dl=0;
! 552: mt=binomial_trail(b); dt=0;
! 553: for(i=0;i<ring_N;i++){ dl+=ml[i]; dt+=mt[i];}
! 554: d=max(d,max(dl,dt));
! 555: }
! 556: gset_deg(g)=d;
! 557: }
! 558:
! 559: /*
! 560: ** int gset_isfacet(gset g,binomial b):
! 561: ** check whether b is known to be a facet or not. if not call lp_isfacet()
! 562: ** to use a linear program to set b's facet flag.
! 563: **
! 564: */
! 565: int lptest=1; /* 1 test only LP, 2 test only flipability, 3 test both */
! 566: int gset_isflippable(gset g,binomial b);
! 567: int gset_isfacet(gset g,binomial b){
! 568: int rval=UNKNOWN;
! 569: rval=b->ff;
! 570: if (rval!=FACET && rval!=NONFACET){
! 571: switch(lptest){
! 572: case 1: rval=(b->ff=lp_isfacet(g,b));
! 573: break;
! 574: case 2: if (gset_isflippable(g,b)==FALSE) rval=(b->ff=NONFACET);
! 575: else rval=(b->ff=FACET);
! 576: break;
! 577: case 3: if (gset_isflippable(g,b)==FALSE) rval=(b->ff=NONFACET);
! 578: else rval=(b->ff=lp_isfacet(g,b));
! 579: break;
! 580: default: fprintf(stderr,"Warning lptest has invalid value in gset_isflippable()\n");
! 581: exit(4);
! 582: }
! 583: }
! 584: return rval;
! 585: }
! 586:
! 587: /*
! 588: **
! 589: */
! 590: int gset_isflippable(gset g,binomial b){
! 591: int ctg=0,ctn=1,val=TRUE;
! 592: binomial new=0,tmp,ptrg,ptrn,btmp;
! 593: /* copy all binomials of g's list ecept b to old and new lists*/
! 594: for(ptrg=gset_first(g);ptrg!=0;ptrg=ptrg->next){
! 595: ctg++;
! 596: if (ptrg!=b){
! 597: tmp=binomial_new();
! 598: binomial_copy(ptrg,tmp);
! 599: monomial_set(tmp);
! 600: binomial_next(tmp)=new;
! 601: new=tmp;
! 602: }
! 603: }
! 604: btmp=binomial_new();
! 605: binomial_copy(b,btmp);
! 606:
! 607:
! 608: /* calculate the reduced GB of the binomial b and the monomials in M */
! 609: bmrgb(btmp,&new);
! 610:
! 611: /* test if new and old have same leading terms */
! 612: for(ptrn=new ;ptrn!=0;ptrn=ptrn->next) ctn++;
! 613: if (ptrn!=ptrg) val=FALSE; /* initial ideals have diferent numbers of generators */
! 614: else {
! 615: for(ptrn=new;ptrn!=0;ptrn=ptrn->next){
! 616: for(ptrg=gset_first(g);ptrg!=0;ptrg=ptrg->next){
! 617: if (monomial_equal(binomial_lead(ptrn),binomial_lead(ptrg))==TRUE)break;
! 618: }
! 619: if (ptrg==0)break; /* no match found */
! 620: }
! 621: if (ptrn!=0)val=FALSE; /* not all terms had matches found */
! 622: }
! 623:
! 624: /* clean up */
! 625: binomial_free(btmp);
! 626: while(new!=0){
! 627: ptrn=new;
! 628: new=new->next;
! 629: binomial_free(ptrn);
! 630: }
! 631:
! 632: return val;
! 633: }
! 634:
! 635: /*
! 636: ** int lp_isfacet(gset g,binomial b):
! 637: ** set up and run linear program to set b's facet flag.
! 638: **
! 639: */
! 640: int lp_isfacet(gset g,binomial b){
! 641: binomial ptr;
! 642: int D=ring_N,M=0;
! 643: int idx=0,i;
! 644: double dtmp0,dtmp;
! 645:
! 646: /* need cheeper way to determine the size of g */
! 647: for(ptr=gset_first(g);ptr!=0;ptr=binomial_next(ptr)){
! 648: if (ptr->ff!=NONFACET)M++;
! 649: }
! 650:
! 651: /* make sure the lp structures have enough space for the program*/
! 652: LP_get_space(M,D+2*M);
! 653:
! 654: /* Set up linear program. */
! 655: ptr=gset_first(g);
! 656: idx=0;
! 657: while(ptr!=0){
! 658: if (ptr->ff!=NONFACET){
! 659: for(i=0;i<D;i++){ /* make sure supports of trail and lead disjoint? */
! 660: LP_A(idx,i)=binomial_lead(ptr)[i]-binomial_trail(ptr)[i];
! 661: }
! 662: for(i=0;i<M;i++){
! 663: LP_A(idx,D+i)=0;
! 664: LP_A(idx,D+M+i)=0;
! 665: }
! 666: if (ptr==b) {for(i=0;i<D;i++) LP_A(idx,i)*=-1;}
! 667: LP_A(idx,D+M+idx)=1;
! 668: LP_A(idx,D+idx)=-1;
! 669: LP_B[idx]=1;
! 670: idx++;
! 671: }
! 672: ptr=ptr->next;
! 673: }
! 674: for(i=0;i<D+M;i++){
! 675: LP_C[i]=0;
! 676: LP_X[i]=0;
! 677: LP_NonBasis[i]=i;
! 678: }
! 679: for(i=0;i<M;i++){
! 680: LP_Basis[i]=D+M+i;
! 681: LP_X[D+M+i]=1;
! 682: LP_C[D+M+i]=1;
! 683: }
! 684:
! 685: /* Run simplex algorithm */
! 686: Rsimp(LP_M,LP_N,LP_A,LP_B,LP_C,LP_X,LP_Basis,LP_NonBasis,
! 687: LP_R,LP_Q,LP_t1,LP_t2);
! 688:
! 689: /* use results to determine if the b is a facet*/
! 690: b->ff=FACET;
! 691: for(i=0;i<M;i++) if (LP_Basis[i]>=M+D) {b->ff=NONFACET;break;}
! 692:
! 693: return b->ff;
! 694: }
! 695:
! 696:
! 697:
! 698:
! 699:
! 700:
! 701:
! 702:
! 703:
! 704:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gset.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / TiGERS_0.9