src/baseparsuif/suifmath/linear_ineq.cc

   1 /* file "linear_ineq.cc" */
   2
   3 /*  Copyright (c) 1994 Stanford University
   4
   5     All rights reserved.
   6
   7     This software is provided under the terms described in
   8     the "suif_copyright.h" include file. */
   9
  10 #include <suif_copyright.h>
  11
  12 #define _MODULE_ "libsuifmath.a"
  13 #pragma implementation "linear_ineq.h"
  14
  15 #include <stdlib.h>
  16 #include <string.h>
  17 #include <math.h>
  18 #include <suif1.h>
  19 #include "suifmath.h"
  20
  21
  22 no_integer_result_hook_f NIR_hook = NULL;
  23 integer_solver_hook_f    IS_hook  = NULL;
  24
  25
  26
  27 /* ##################################################
  28    #####   constraint                           #####
  29    ################################################## */
  30
  31
  32
  33 /********************************************************************************
  34  * complement                                                                   *
  35  *                                                                              *
  36  * when the raw represent R>=0  -> R>-1 compliment is R<=-1 -> -1-R>=0          *
  37  ********************************************************************************/
  38 void constraint::complement()
  39 {
  40     (*this) *= -1;
  41     (*this)[0]--;
  42 }
  43
  44
  45
  46 /********************************************************************************
  47  * normalize                                                                    *
  48  ********************************************************************************/
  49 /* if (!norm_bounds), this /= row_gcd();
  50    if (norm_bounds),
  51       g = row_gcd(all but 0)
  52       this /= g; (take floor of this[0]/g);
  53 */
  54
  55 boolean constraint::normalize(boolean norm_bounds)
  56 {
  57     int g = -1;
  58     for(int i=(norm_bounds)?1:0; i<n(); i++)
  59         if((*this)[i] != 0) {
  60             if(g == -1)
  61                 g = ABS((*this)[i]);
  62             else
  63                 g = gcd(g, ABS((*this)[i]));
  64         }
  65
  66     if(g > 1) {
  67         if(norm_bounds) {
  68             for(int i=1; i<n(); i++)
  69                 R[i] /= g;
  70             int new0 = (int) floor(((double)R[0])/(double)g);
  71             boolean chg = ((g*new0) != R[0]);
  72             R[0] = new0;
  73             return chg;
  74         } else
  75             (*this) /= g;
  76         return FALSE;
  77     } else
  78         return FALSE;
  79 }
  80
  81
  82
  83 /********************************************************************************
  84  * highest non-zero value                                                       *
  85  ********************************************************************************/
  86 int constraint::rank() const
  87 {
  88
  89     const int * curr = &(R[n()-1]);
  90     int i = n();
  91     while(i-->0)
  92         if(*curr-- != 0) return i;
  93     return 0;
  94
  95 /* speed hack above
  96     int rnk = 0;
  97     for(int cc = 0; cc<n(); cc++)
  98         if(c(cc) != 0) rnk=cc;
  99     return rnk; */
 100 }
 101
 102 /********************************************************************************
 103  * highest order                                                                *
 104  * return val s.t. For all c R[c]!=0 val +=sort_order[c]                        *
 105  ********************************************************************************/
 106 int constraint::highest_order(const constraint & sort_order) const
 107 {
 108     assert(sort_order.n() == this->n());
 109
 110     int hi = 0;
 111     for(int cc = 0; cc<n(); cc++)
 112         if(c(cc) != 0) hi = hi + sort_order.c(cc);
 113
 114     return hi;
 115 }
 116
 117 /********************************************************************************
 118  * Only one non-zero variable (1st location is a const, thus can have any val.  *
 119  ********************************************************************************/
 120 boolean constraint::unique() const
 121 {
 122     boolean found = FALSE;
 123     for(int cc = 1; cc<n(); cc++)
 124         if(c(cc) != 0) {
 125             if (found)
 126                 return FALSE; // more than one
 127             else
 128                 found = TRUE; // found first
 129         }
 130     if(found)
 131         return TRUE;  // last one was the only var
 132     else
 133         return FALSE; // Only constant
 134 }
 135
 136
 137 /********************************************************************************
 138  * gcd and lcm                                                                  *
 139  ********************************************************************************/
 140 int constraint::row_gcd() const
 141 {
 142     int g = 0;
 143     for(int i = 0; i<n(); i++)
 144         if(c(i) != 0) {
 145             if(g==0)
 146                 g = ABS(c(i));
 147             else
 148                 g = gcd(g, ABS(c(i)));
 149         }
 150     return g;
 151 }
 152
 153 int constraint::row_lcm() const
 154 {
 155     int l = 0;
 156     for(int i = 0; i<n(); i++)
 157         if(c(i) != 0) {
 158             if(l==0)
 159                 l = ABS(c(i));
 160             else
 161                 l = lcm(l, ABS(c(i)));
 162         }
 163     return l;
 164 }
 165
 166
 167 /*
 168    if (type==0)
 169    then result[i]=1 if lower<=this[i]<=upper
 170    else result[i]=1 if lower<=leftofthis[i]<=upper
 171                     where leftofthis[i] = sum(0,i,this[i]==type)-1
 172 */
 173
 174 constraint constraint::make_filter(int type, int lower, int upper) const
 175 {
 176     constraint reply(n());
 177
 178     if(type == 0) { // array vars;
 179         for(int i = 0; i<n(); i++)
 180             reply[i] = ((c(i) >= lower)&&(c(i) <= upper))?1:0;
 181     } else {
 182         int count = -1;
 183         for(int i = 0; i<n(); i++) {
 184             if(type==c(i)) count++;
 185             reply[i] = ((count >= lower)&&(count <= upper))?1:0;
 186         }
 187     }
 188
 189     return reply;
 190 }
 191
 192 /*
 193    returns AND(0<=i<n(),kernal[i] && same_sign(this[i],sgn))
 194 */
 195
 196 int constraint::filter(const constraint & kernal, int sgn) const
 197 {
 198     for(int i = 0; i<n(); i++)
 199         if(kernal.c(i) == 1)                      // mark in kernal
 200             if(c(i) != 0)                 // entry in constraint
 201                 if(c(i)*sgn >= 0)        // sign +/- check sign of entry
 202                     return 1;                   // if sign=0 any sign is ok
 203     return 0;
 204 }
 205
 206
 207
 208 /* ##################################################
 209    #####   lin_ineq                             #####
 210    ################################################## */
 211
 212
 213 integer_row *lin_ineq::alloc_rows(int num_rows)
 214 {
 215     realmm = num_rows;
 216     return new constraint[num_rows];
 217 }
 218
 219
 220 /********************************************************************************
 221  * Column Operations                                                            *
 222  ********************************************************************************/
 223 int lin_ineq::col_lcm(int col) const
 224 {
 225     assert((col>=0)&&(col<nn));
 226     int l = 0;
 227     for(int i = 0; i<mm; i++)
 228         if(r(i).c(col) != 0) {
 229             if(l==0)
 230                 l = ABS(r(i).c(col));
 231             else
 232                 l = lcm(l, ABS(r(i).c(col)));
 233         }
 234     return l;
 235 }
 236
 237
 238 int lin_ineq::col_gcd(int col) const
 239 {
 240     assert((col>=0)&&(col<nn));
 241     int g = 0;
 242     for(int i = 0; i<mm; i++)
 243         if(r(i).c(col) != 0) {
 244             if(g==0)
 245                 g = ABS(r(i).c(col));
 246             else
 247                 g = gcd(g, ABS(r(i).c(col)));
 248         }
 249     return g;
 250 }
 251
 252
 253
 254 /********************************************************************************
 255  * Operators                                                                       *
 256  ********************************************************************************/
 257 lin_ineq lin_ineq::operator&&(const lin_ineq & mat) const
 258 {
 259     if(m() == 0) return mat;
 260     if(mat.m() == 0) return *this;
 261
 262     assert(n() == mat.n());
 263
 264     lin_ineq result(this, m()+mat.m());
 265     for(int i=0; i<mat.m(); i++)
 266         result[i+m()] = mat.r(i);
 267
 268     return result;
 269 }
 270
 271 lin_ineq &lin_ineq::operator&=(const lin_ineq & mat)
 272 {
 273     if(mat.m() == 0) return *this;
 274     if(m() == 0) { init(mat); return *this; };
 275
 276     assert(n() == mat.n());
 277
 278     int oldm = m();
 279     do_add_rows(mat.m());
 280
 281     for(int i=0; i<mat.m(); i++)
 282         (*this)[i+oldm] = mat.r(i);
 283
 284     return *this;
 285 }
 286
 287 // for debugging operator~();
 288 // static int soln_count = 0;
 289
 290 int lin_ineq::operator~() const
 291 {
 292     assert_msg(NIR_hook, ("No integer solver given, "
 293                           "(link in the dependence library)"));
 294 //    fprintf(stderr,"[%d]lin_ineq::~ applied to\n",++soln_count);
 295 //    print(stderr);
 296     boolean ans = (*NIR_hook)(*this, NULL);
 297 //    fprintf(stderr,"result == %d\n",(int)ans);
 298     return ans;
 299 }
 300
 301 boolean lin_ineq::operator>=(const lin_ineq & mat) const
 302 {
 303     assert(n() == mat.n());
 304
 305     if (~mat) return 1;
 306
 307     lin_ineq testln(&mat, mat.m()+1);
 308
 309     for(int i=0; i<m(); i++) {
 310         testln[mat.m()] = -r(i);
 311         if((~testln)==0) return 0;
 312     }
 313
 314     lin_ineq tmp = mat||(*this);
 315     if (~tmp) {
 316         fprintf(stderr,"Error in dependence test!!!\n THIS = \n");
 317         print(stderr);
 318         fprintf(stderr,"MAT = \n");
 319         mat.print(stderr);
 320         fprintf(stderr,"tmp = \n");
 321         tmp.print(stderr);
 322         assert(0);
 323     }
 324
 325     return 1;
 326 }
 327
 328 lin_ineq lin_ineq::operator%(const constraint & row) const
 329 {
 330     return lin_ineq(integer_matrix::operator%(row));
 331 }
 332
 333 /********************************************************************************
 334  * output                                                                       *
 335  ********************************************************************************/
 336 void lin_ineq::print(FILE *fp, int flag) const
 337 {
 338
 339     for (int i=0; i<m(); i++)  {
 340         if(flag == 1) fprintf(fp,"%2d:  ", i);
 341         for (int j=0; j<n(); j++)  {
 342             fprintf(fp, " %6d", r(i).c(j));
 343         }
 344         fprintf(fp,"\n");
 345     }
 346 }
 347
 348
 349 int lin_ineq::max_rank() const
 350 {
 351     int rank = 0;
 352     for(int i=0; i<m(); i++) rank = MAX(rank, r(i).rank());
 353
 354     return rank;
 355 }
 356
 357
 358 /********************************************************************************
 359  * zero                                                                         *
 360  *                                                                              *
 361  ********************************************************************************/
 362 void lin_ineq::del_zeros()
 363 {
 364     int *lst = new int[m()];
 365
 366     int i;
 367     int cnt = 0;
 368     for(i=0; i<m(); i++)
 369         if((r(i).rank() != 0) || (r(i).c(0) < 0)) {
 370             cnt++;
 371             lst[i] = 1;
 372         } else
 373             lst[i] = 0;
 374
 375     if(cnt != m()) {
 376         int src, dst;
 377         for (src=0, dst=0; src < m(); src++) {
 378             if (lst[src]) {
 379                 if (src != dst)
 380                     (*this)[dst] ^= (*this)[src];
 381                 dst++;
 382             }
 383         }
 384         mm = cnt;
 385     }
 386     delete[] lst;
 387 }
 388
 389 /********************************************************************************
 390  * zero                                                                         *
 391  *                                                                              *
 392  ********************************************************************************/
 393 void lin_ineq::del_unused_var()
 394 {
 395     lin_ineq tmpM(this);
 396
 397     integer_row lst(tmpM.n());
 398
 399     int i, j;
 400
 401     lst[0] = 1;
 402     for(i=0; i<tmpM.m(); i++)
 403         for(j=0; j<tmpM.n(); j++) if(tmpM[i][j]) lst[j] = 1;
 404
 405     do_del_columns(lst);
 406 }
 407
 408
 409 /********************************************************************************
 410  * normalization                                                                *
 411  *                                                                              *
 412  ********************************************************************************/
 413
 414 void lin_ineq::normalize(boolean norm_bounds)
 415 {
 416     for(int i=0; i < this->m(); i++)
 417         (*this)[i].normalize(norm_bounds);
 418
 419 }
 420
 421 /********************************************************************************
 422  * repatitions                                                                  *
 423  *                                                                              *
 424  ********************************************************************************/
 425 // rem identical rows;
 426 // as compared to itself and (0,a-1);
 427 void lin_ineq::del_repetition(int a, int b, integer_row *deleted)
 428 {
 429
 430     if ((m() > 0) && (a <= b)) {
 431         assert(0<=a);
 432         assert(b<m());
 433         int *rank = new int[m()];
 434         int i;
 435         for(i=0; i < m(); i++) {
 436             rank[i] = (*this)[i].rank();
 437         }
 438
 439         integer_row *del_rows = deleted ? deleted : new integer_row(m());
 440         if (deleted)
 441             assert(deleted->n() == m());
 442         for(i=a; i <= b; i++)
 443             for(int j=0; j < i; j++) {
 444                 if(rank[i] == rank[j])
 445                     if((*this)[i] == (*this)[j]) {
 446                         (*del_rows)[j] = 1; // two equals
 447                         rank[i] = 0;
 448                         break;
 449                     }
 450             }
 451         delete[] rank;
 452         do_del_rows(*del_rows);
 453         if (!deleted) delete del_rows;
 454     }
 455 }
 456
 457 // rem identical rows in (incl)range;
 458 void lin_ineq::del_repetition()
 459 {
 460     if(m() <= 1) return;
 461
 462     del_repetition(0, m()-1);
 463 }
 464
 465 // rem otherwise identical row with bigger constant;
 466 // in (incl) range given;
 467 void lin_ineq::del_loose(int r1, int r2, integer_row *deleted)
 468 {
 469     if (r1 > r2) return;
 470
 471     integer_row *del_rows = deleted ? deleted : new integer_row(m());
 472     if (deleted) {
 473         assert (deleted->n() == m());
 474     }
 475     int *this_rank = new int[m()];
 476
 477     assert ((0<=r1) && (r1<=r2) && (r2<m()));
 478
 479     for(int j = 0; j<m(); j++)
 480         this_rank[j] = (*this)[j].rank();
 481
 482     for (int i=r1; i <= r2; i++)
 483         if (((*del_rows)[i]==0) && ((*this)[i].rank() != 0)) {
 484             constraint A = (*this)[i];
 485             int c1 = A[0];
 486             int A_rank = A.rank();
 487             A[0] = 0;
 488             for (int j=0; j < i; j++) {
 489                 if (A_rank == this_rank[j]) {
 490                     constraint B = (*this)[j];
 491                     int c2 = (*this)[j][0];
 492                     B[0] = 0;
 493                     if (A == B) {
 494                         if (c2 > c1)
 495                             (*del_rows)[j] = 1;
 496                         else {
 497                             (*del_rows)[i] = 1;
 498                             break;
 499                         };
 500                     };
 501                 };
 502             };
 503         };
 504     delete[] this_rank;
 505     do_del_rows(*del_rows);
 506     if (!deleted) delete del_rows;
 507 }
 508
 509 // rem otherwise identical row with bigger constant;
 510 void lin_ineq::del_loose()
 511 {
 512     if (m() <= 1) return;
 513     del_loose(0,m()-1);
 514 }
 515
 516 void lin_ineq::row_swap(int i, int j)
 517 {
 518     assert((i>=0)&&(i<m()));
 519     assert((j>=0)&&(j<m()));
 520
 521     A[i] ^= A[j];
 522 }
 523
 524
 525
 526
 527 /********************************************************************************
 528  * sort                                                                         *
 529  *                                                                              *
 530  * If sort order is given sort in the increasing order of f(M[i])               *
 531  *   f(M[i]) = val s.t. For all c M[i][c]!=0 val=MAX(s_order[c])                *
 532  * else Sort S in the order of increasing constraint rank                       *
 533  ********************************************************************************/
 534 struct rh {
 535     int rnk;
 536     int high;
 537
 538     void rh_swap(rh & o) {
 539         int tmp;
 540         tmp = rnk;  rnk  = o.rnk;  o.rnk  = tmp;
 541         tmp = high; high = o.high; o.high = tmp;
 542     }
 543 };
 544
 545 // if s_order, sort according to
 546 //   row.highest_order(s_order)+row.rank();
 547 //   (recall row.highest_order(s_order) = Sum{i s.t. row[i]!=0} s_order[i]);
 548 // dir > 0 => sort into row[i] <= row[i+1];
 549 void lin_ineq::sort(const constraint * s_order, int dir)
 550 {
 551     if(this->m() <= 1) return;
 552
 553     if(s_order) assert(s_order->n() == this->n());
 554
 555     rh * RH = new rh[this->m()];
 556
 557     int i;
 558     for(i=0; i<this->m(); i++) {
 559         RH[i].rnk  = (*this)[i].rank();
 560         RH[i].high = (s_order)?(*this)[i].highest_order(*s_order):0;
 561     }
 562
 563 //    s_order->print(stderr);
 564 //    for(int f=0; f<this->m(); f++)
 565 //        fprintf(stderr,"(%d %d+%d) ", f, (*this)[f].highest_order(*s_order),
 566 //                                      (*this)[f].rank());
 567 //    fprintf(stderr,"\n");
 568
 569     for(i=0; i<this->m()-1; i++) {
 570         int swap = i;
 571         for(int j=i+1; j<this->m(); j++)
 572             if(s_order) {
 573                 if(dir > 0) {
 574                     if(RH[swap].high+RH[swap].rnk >  RH[j].high+RH[j].rnk)
 575                         swap = j;
 576                 } else {
 577                     if(RH[swap].high+RH[swap].rnk <  RH[j].high+RH[j].rnk)
 578                         swap = j;
 579                 }
 580             } else {
 581                 if(dir > 0) {
 582                     if(RH[swap].rnk > RH[j].rnk) swap = j;
 583                 } else {
 584                     if(RH[swap].rnk < RH[j].rnk) swap = j;
 585                 }
 586             }
 587         if(swap != i) {
 588             RH[i].rh_swap(RH[swap]);
 589             this->row_swap(i, swap);
 590         }
 591     }
 592
 593     delete[] RH;
 594 }
 595
 596 /********************************************************************************
 597  * filter                                                                       *
 598  ********************************************************************************/
 599
 600 /*
 601    deletes rows for which filter(kernal,sgn) is not true;
 602    */
 603 void lin_ineq::do_filter_thru(const constraint & kernal, int sign)
 604 {
 605     if (m() > 0) {
 606         integer_row del_rows(m());
 607
 608         for(int i=0; i<m(); i++)
 609             if(!(*this)[i].filter(kernal, sign))
 610                 del_rows[i] = 1;
 611         do_del_rows(del_rows);
 612     }
 613 }
 614
 615 lin_ineq lin_ineq::filter_thru(const constraint & kernal, int sign) const
 616 {
 617     lin_ineq result(this);
 618     result.do_filter_thru(kernal,sign);
 619     return result;
 620 }
 621
 622 /*
 623    deletes rows for which filter(kernal,sgn) is true;
 624    */
 625 void lin_ineq::do_filter_away(const constraint & kernal, int sign)
 626 {
 627     if (m() > 0) {
 628         integer_row del_rows(m());
 629
 630         for(int i=0; i<m(); i++)
 631             if((*this)[i].filter(kernal, sign))
 632                 del_rows[i] = 1;
 633         do_del_rows(del_rows);
 634     }
 635 }
 636
 637 lin_ineq lin_ineq::filter_away(const constraint & kernal, int sign) const
 638 {
 639     lin_ineq result(this);
 640     result.do_filter_away(kernal,sign);
 641     return result;
 642 }
 643
 644 /*
 645    Takes a row with a non-0 for every variable whose bounds matter.
 646    Modifies the row to take the transitive closure over this leq,
 647    adding every variable related to those vars.
 648    */
 649 void lin_ineq::do_bounds(integer_row &bounds0) const
 650 {
 651     assert(bounds0.n() == n());
 652
 653     boolean changed = TRUE;
 654     while (changed) {
 655         changed = FALSE;
 656
 657         int i;
 658         for (i=1; i<n(); i++) {
 659             if (bounds0[i]) {
 660                 for (int j=0; j<m(); j++) {
 661                     const integer_row &thisrow = r(j);
 662                     if (thisrow.c(i)) {
 663                         for (int i2=1; i2<n(); i2++) {
 664                             if ((i2 != i) && thisrow.c(i2)) {
 665                                 if (!bounds0[i2]) {
 666                                     bounds0[i2] = 1;
 667                                     changed = TRUE;
 668                                 };
 669                             };
 670                         }
 671                     }
 672                 }
 673             }
 674         }
 675     }
 676 }
 677
 678 /* ##################################################
 679    #####   lin_ineq_negate_iter                 #####
 680    ################################################## */
 681 void lin_ineq_negate_iter::init(const lin_ineq & l)
 682 {
 683     assert(l.m() < 32);
 684     L = l;
 685
 686     curr = 0x01;
 687     prev = 0x00;
 688
 689     done_mark = 0x01 << L.m();
 690 }
 691
 692
 693 lin_ineq * lin_ineq_negate_iter::step()
 694 {
 695     unsigned long chang;
 696
 697     chang = curr ^ prev;
 698
 699     for(int i=0; i<L.m(); i++) {
 700         if(chang & 0x01) L[i] = -L[i];
 701         chang >>= 1;
 702     }
 703     prev = curr;
 704     curr = curr + 0x01;
 705
 706     return &L;
 707 }
 708
 709
 710 /* ##################################################
 711    #####   lin_ineq_difference_iter             #####
 712    ################################################## */
 713
 714 void lin_ineq_difference_iter::init(const lin_ineq & A, const lin_ineq & B)
 715 {
 716     D = A & B;
 717     if(D.m()==0) {
 718         indep = TRUE;
 719         D = A;
 720     } else {
 721         indep = FALSE;
 722
 723         D.del_zeros();
 724         D.del_repetition();
 725         D.del_loose();
 726         D.sort();
 727
 728         // remove constraints of D that are also bounds of A
 729         integer_matrix m1A(1,A.n());
 730         lin_ineq A1 = Compose(2, 1,
 731                               & A,
 732                               & m1A);
 733         for(int i=0; i<D.m(); i++) {
 734             A1[A1.m()-1] = -D[i];
 735             if(~A1) D[i] = 0;
 736         }
 737         D.del_zeros();
 738     }
 739 }
 740
 741
 742 lin_ineq * lin_ineq_difference_iter::step()
 743 {
 744     curr++;
 745     if(indep) {
 746         return new lin_ineq(D);
 747     } else {
 748         lin_ineq * ret = new lin_ineq(curr, D.n());
 749
 750         ret = new lin_ineq(D.resize(0, curr, 0, D.n()));
 751         (*ret)[curr] = -(*ret)[curr];
 752
 753         return ret;
 754     }
 755 }
 756
 757 /* ##################################################
 758    #####   lin_ineq_op                          #####
 759    ################################################## */
 760
 761
 762 /********************************************************************************
 763  *  constructor                                                                 *
 764  ********************************************************************************/
 765 lin_ineq_op::lin_ineq_op(int sz, const int * k)
 766 {
 767     kind.init(sz);
 768
 769     kind[0] = NM_CONSTANT;
 770
 771     for(int i=1; i<kind.n(); i++)
 772         kind[i] = k[i-1];
 773 }
 774
 775 lin_ineq_op::lin_ineq_op(int a, int b)
 776 {
 777     kind.init(a+b);
 778
 779     kind[0] = NM_CONSTANT;
 780     int i;
 781     for(i=1; i<a; i++) kind[i] = NM_SYMBOLS;
 782     for(i=a; i<a+b; i++)   kind[i] = NM_LOCATIONS;
 783 }
 784
 785
 786 /********************************************************************************
 787  *  filtering                                                                   *
 788  ********************************************************************************/
 789
 790 lin_ineq lin_ineq_op::get(const lin_ineq & LEQ, int nloc, int sgn) const
 791 {
 792     /* flt[i] = 1 if there are nloc+1 NM_LOCATIONS left of column i in kind */
 793     constraint flt = kind.make_filter(NM_LOCATIONS, nloc);
 794
 795     /* removes rows not containing right signed value in column of flt */
 796     lin_ineq result = LEQ.filter_thru(flt, sgn);
 797
 798     return result;
 799 }
 800
 801 lin_ineq lin_ineq_op::rem(const lin_ineq & LEQ, int nloc, int sgn) const
 802 {
 803     constraint flt = kind.make_filter(NM_LOCATIONS, nloc);
 804
 805     lin_ineq result = LEQ.filter_away(flt, sgn);
 806
 807     return result;
 808 }
 809
 810
 811 // keeps only rows with a + value in a column i;
 812 // with at least nloc+1 NM_LOCATIONS to the left of it in kind;
 813 lin_ineq lin_ineq_op::get_lesser(const lin_ineq & LEQ, int nloc) const
 814 {
 815     constraint flt = kind.make_filter(NM_LOCATIONS, nloc, LEQ.n());
 816
 817     lin_ineq result = LEQ.filter_away(flt, 0);
 818
 819     return result;
 820 }
 821
 822 // opposite of get_lesser();
 823 lin_ineq lin_ineq_op::get_greater(const lin_ineq & LEQ, int nloc) const
 824 {
 825     constraint flt = kind.make_filter(NM_LOCATIONS, nloc, LEQ.n());
 826
 827     lin_ineq result = LEQ.filter_thru(flt, 0);
 828
 829     return result;
 830 }
 831
 832
 833 lin_ineq lin_ineq_op::get_exclusive(const lin_ineq & LEQ, int nloc, int sgn)
 834         const
 835 {
 836     // get rows with a value of right sign in nloc+1 NM_LOCATIONS;
 837     lin_ineq result = get(LEQ, nloc, sgn);
 838
 839     constraint flt = kind.make_filter(NM_LOCATIONS) - kind.make_filter(NM_LOCATIONS, nloc);
 840
 841     // get rid of such rows with any non-0 value in another NM_LOCATIONS col;
 842     result = result.filter_away(flt, 0);
 843
 844     return result;
 845 }
 846
 847
 848
 849 /********************************************************************************
 850  *  overlap or contigous ranges in given dim                                    *
 851  *                                                                              *
 852  *  - get both ineq's with given array dim present and no other array dims      *
 853  *  - If no subset of ineq`s exist like that, cannot get an indipendent         *
 854  *    variable which might overlap or fuse toghter -> good bye                  *
 855  *  - If there is an overlap of the ind. dimention all the time can merge       *
 856  *  - Now check if the two dimentions can be fused toghter. i.e. no gap in the  *
 857  *    middle of the two ranges.                                                 *
 858  ********************************************************************************/
 859 lin_ineq lin_ineq_op::overlap_or_contig(const lin_ineq & a, const lin_ineq & b,
 860                                         int nloc) const
 861 {
 862     lin_ineq result;
 863
 864     lin_ineq curr_a;
 865     lin_ineq curr_b;
 866
 867     // get rows from a/b with non-0 only in column of nloc+1 NM_LOCATIONS;
 868     // but in no other NM_LOCATIONS columns;
 869     curr_a = get_exclusive(a, nloc, 0);
 870     curr_b = get_exclusive(b, nloc, 0);
 871
 872     if((curr_a.is_empty())||(curr_b.is_empty())) return result;
 873
 874     // Now check if a is after b
 875
 876     // get rows from a with + values in nloc+1 columns;
 877     curr_a = get_exclusive_lower(a, nloc);
 878     // get rows from b with - values in nloc+1 columns;
 879     curr_b = get_exclusive_upper(b, nloc);
 880
 881     int gap1 = check_gap(curr_a, curr_b);
 882     if(gap1 < 0) return result;
 883
 884     // vice versa;
 885     curr_a = get_exclusive_upper(a, nloc);
 886     curr_b = get_exclusive_lower(b, nloc);
 887
 888     int gap2 = check_gap(curr_a, curr_b);
 889     if(gap2 < 0) return result;
 890
 891     if(gap1 < gap2)
 892         result = rem_lower(a, nloc) || rem_upper(b, nloc);  // remove lower_a and upper_b
 893     else
 894         result = rem_upper(a, nloc) || rem_lower(b, nloc);  // remove upper_a and lower_b
 895
 896     return result;
 897 }
 898
 899
 900 /********************************************************************************
 901  *  check_gap                                                                   *
 902  *                                                                              *
 903  * input:       for a single array dim. a set with upper bounds and the         *
 904  *              other set with lower bounds                                     *
 905  *              Order does not matter (if a is LB then b is UB vice vesa)       *
 906  *                                                                              *
 907  * Theory:                                                                      *
 908  *         ------------|---|-------                                             *
 909  *                     a   b                                                    *
 910  *      From LB  i <= a     a - i >= 0                                          *
 911  *      From UB  i >= b    -b + i >= 0                                          *
 912  *              add:        a - b                                               *
 913  *      The gap between rangers cannot be grater than 1   a >= b - 1            *
 914  *                          a - b >= -1                                         *
 915  *                                                                              *
 916  * Method:                                                                      *
 917  *    For each constraint L in Lower bound set of ineqs                         *
 918  *        For each constraint U in Upper bound set of ineqs                     *
 919  *            k = U + L                                                         *
 920  *            if k only is a constant (rank(k) == 0)                            *
 921  *                if const(k) >= -1                                             *
 922  *                    dist = min(dist, cosnt(k))                                *
 923  *                    delete U from set of upper bounds                         *
 924  *    if gap is no grater than 1 U shoud be empty.                              *
 925  *    Now do the same to L                                                      *
 926  ********************************************************************************/
 927 int lin_ineq_op::check_gap(const lin_ineq & a, const lin_ineq & b) const
 928 {
 929     lin_ineq liq_a(a);
 930     lin_ineq liq_b(b);
 931     constraint k;
 932     int dist = 10000;
 933
 934     for(int i = 0; i < a.num_constraint(); i++)
 935         for(int j = 0; j < b.num_constraint(); j++) {
 936             k = a.r(i) + b.r(j);
 937             if(k.rank() == 0)
 938                 if(k.c(0) >= -1) {
 939                     dist = MIN(dist+2, k.c(0));
 940                     liq_a[i] = 0;
 941                     liq_b[j] = 0;
 942                 }
 943         }
 944
 945     liq_a.del_zeros();
 946     liq_b.del_zeros();
 947
 948     if((!liq_a.is_empty())||(!liq_a.is_empty())) return -1;
 949
 950     return dist;
 951 }