lib/CodeGen/PBQP.cpp

   1 //===---------------- PBQP.cpp --------- PBQP Solver ------------*- C++ -*-===//
   2 //
   3 //                     The LLVM Compiler Infrastructure
   4 //
   5 // This file is distributed under the University of Illinois Open Source
   6 // License. See LICENSE.TXT for details.
   7 //
   8 //===----------------------------------------------------------------------===//
   9 //
  10 // Developed by:                   Bernhard Scholz
  11 //                             The University of Sydney
  12 //                         http://www.it.usyd.edu.au/~scholz
  13 //===----------------------------------------------------------------------===//
  14
  15
  16 #include <limits>
  17 #include <cassert>
  18 #include <cstring>
  19
  20 #include "PBQP.h"
  21 #include "llvm/Config/alloca.h"
  22
  23 namespace llvm {
  24
  25 /**************************************************************************
  26  * Data Structures
  27  **************************************************************************/
  28
  29 /* edge of PBQP graph */
  30 typedef struct adjnode {
  31   struct adjnode *prev,      /* doubly chained list */
  32                  *succ,
  33                  *reverse;   /* reverse edge */
  34   int adj;                   /* adj. node */
  35   PBQPMatrix *costs;         /* cost matrix of edge */
  36
  37   bool tc_valid;              /* flag whether following fields are valid */
  38   int *tc_safe_regs;          /* safe registers */
  39   int tc_impact;              /* impact */
  40 } adjnode;
  41
  42 /* bucket node */
  43 typedef struct bucketnode {
  44   struct bucketnode *prev;   /* doubly chained list */
  45   struct bucketnode *succ;
  46   int u;                     /* node */
  47 } bucketnode;
  48
  49 /* data structure of partitioned boolean quadratic problem */
  50 struct pbqp {
  51   int num_nodes;             /* number of nodes */
  52   int max_deg;               /* maximal degree of a node */
  53   bool solved;               /* flag that indicates whether PBQP has been solved yet */
  54   bool optimal;              /* flag that indicates whether PBQP is optimal */
  55   PBQPNum min;
  56   bool changed;              /* flag whether graph has changed in simplification */
  57
  58                              /* node fields */
  59   PBQPVector **node_costs;   /* cost vectors of nodes */
  60   int *node_deg;             /* node degree of nodes */
  61   int *solution;             /* solution for node */
  62   adjnode **adj_list;        /* adj. list */
  63   bucketnode **bucket_ptr;   /* bucket pointer of a node */
  64
  65                              /* node stack */
  66   int *stack;                /* stack of nodes */
  67   int stack_ptr;             /* stack pointer */
  68
  69                              /* bucket fields */
  70   bucketnode **bucket_list;  /* bucket list */
  71
  72   int num_r0;                /* counters for number statistics */
  73   int num_ri;
  74   int num_rii;
  75   int num_rn;
  76   int num_rn_special;
  77 };
  78
  79 bool isInf(PBQPNum n) { return n == std::numeric_limits<PBQPNum>::infinity(); }
  80
  81 /*****************************************************************************
  82  * allocation/de-allocation of pbqp problem
  83  ****************************************************************************/
  84
  85 /* allocate new partitioned boolean quadratic program problem */
  86 pbqp *alloc_pbqp(int num_nodes)
  87 {
  88   pbqp *this_;
  89   int u;
  90
  91   assert(num_nodes > 0);
  92
  93   /* allocate memory for pbqp data structure */
  94   this_ = (pbqp *)malloc(sizeof(pbqp));
  95
  96   /* Initialize pbqp fields */
  97   this_->num_nodes = num_nodes;
  98   this_->solved = false;
  99   this_->optimal = true;
 100   this_->min = 0.0;
 101   this_->max_deg = 0;
 102   this_->changed = false;
 103   this_->num_r0 = 0;
 104   this_->num_ri = 0;
 105   this_->num_rii = 0;
 106   this_->num_rn = 0;
 107   this_->num_rn_special = 0;
 108
 109   /* initialize/allocate stack fields of pbqp */
 110   this_->stack = (int *) malloc(sizeof(int)*num_nodes);
 111   this_->stack_ptr = 0;
 112
 113   /* initialize/allocate node fields of pbqp */
 114   this_->adj_list = (adjnode **) malloc(sizeof(adjnode *)*num_nodes);
 115   this_->node_deg = (int *) malloc(sizeof(int)*num_nodes);
 116   this_->solution = (int *) malloc(sizeof(int)*num_nodes);
 117   this_->bucket_ptr = (bucketnode **) malloc(sizeof(bucketnode **)*num_nodes);
 118   this_->node_costs = (PBQPVector**) malloc(sizeof(PBQPVector*) * num_nodes);
 119   for(u=0;u<num_nodes;u++) {
 120     this_->solution[u]=-1;
 121     this_->adj_list[u]=NULL;
 122     this_->node_deg[u]=0;
 123     this_->bucket_ptr[u]=NULL;
 124     this_->node_costs[u]=NULL;
 125   }
 126
 127   /* initialize bucket list */
 128   this_->bucket_list = NULL;
 129
 130   return this_;
 131 }
 132
 133 /* free pbqp problem */
 134 void free_pbqp(pbqp *this_)
 135 {
 136   int u;
 137   int deg;
 138   adjnode *adj_ptr,*adj_next;
 139   bucketnode *bucket,*bucket_next;
 140
 141   assert(this_ != NULL);
 142
 143   /* free node cost fields */
 144   for(u=0;u < this_->num_nodes;u++) {
 145     delete this_->node_costs[u];
 146   }
 147   free(this_->node_costs);
 148
 149   /* free bucket list */
 150   for(deg=0;deg<=this_->max_deg;deg++) {
 151     for(bucket=this_->bucket_list[deg];bucket!=NULL;bucket=bucket_next) {
 152       this_->bucket_ptr[bucket->u] = NULL;
 153       bucket_next = bucket-> succ;
 154       free(bucket);
 155     }
 156   }
 157   free(this_->bucket_list);
 158
 159   /* free adj. list */
 160   assert(this_->adj_list != NULL);
 161   for(u=0;u < this_->num_nodes; u++) {
 162     for(adj_ptr = this_->adj_list[u]; adj_ptr != NULL; adj_ptr = adj_next) {
 163       adj_next = adj_ptr -> succ;
 164       if (u < adj_ptr->adj) {
 165         assert(adj_ptr != NULL);
 166         delete adj_ptr->costs;
 167       }
 168       if (adj_ptr -> tc_safe_regs != NULL) {
 169            free(adj_ptr -> tc_safe_regs);
 170       }
 171       free(adj_ptr);
 172     }
 173   }
 174   free(this_->adj_list);
 175
 176   /* free other node fields */
 177   free(this_->node_deg);
 178   free(this_->solution);
 179   free(this_->bucket_ptr);
 180
 181   /* free stack */
 182   free(this_->stack);
 183
 184   /* free pbqp data structure itself */
 185   free(this_);
 186 }
 187
 188
 189 /****************************************************************************
 190  * adj. node routines
 191  ****************************************************************************/
 192
 193 /* find data structure of adj. node of a given node */
 194 static
 195 adjnode *find_adjnode(pbqp *this_,int u,int v)
 196 {
 197   adjnode *adj_ptr;
 198
 199   assert (this_ != NULL);
 200   assert (u >= 0 && u < this_->num_nodes);
 201   assert (v >= 0 && v < this_->num_nodes);
 202   assert(this_->adj_list != NULL);
 203
 204   for(adj_ptr = this_ -> adj_list[u];adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
 205     if (adj_ptr->adj == v) {
 206       return adj_ptr;
 207     }
 208   }
 209   return NULL;
 210 }
 211
 212 /* allocate a new data structure for adj. node */
 213 static
 214 adjnode *alloc_adjnode(pbqp *this_,int u, PBQPMatrix *costs)
 215 {
 216   adjnode *p;
 217
 218   assert(this_ != NULL);
 219   assert(costs != NULL);
 220   assert(u >= 0 && u < this_->num_nodes);
 221
 222   p = (adjnode *)malloc(sizeof(adjnode));
 223   assert(p != NULL);
 224
 225   p->adj = u;
 226   p->costs = costs;
 227
 228   p->tc_valid= false;
 229   p->tc_safe_regs = NULL;
 230   p->tc_impact = 0;
 231
 232   return p;
 233 }
 234
 235 /* insert adjacence node to adj. list */
 236 static
 237 void insert_adjnode(pbqp *this_, int u, adjnode *adj_ptr)
 238 {
 239
 240   assert(this_ != NULL);
 241   assert(adj_ptr != NULL);
 242   assert(u >= 0 && u < this_->num_nodes);
 243
 244   /* if adjacency list of node is not empty -> update
 245      first node of the list */
 246   if (this_ -> adj_list[u] != NULL) {
 247     assert(this_->adj_list[u]->prev == NULL);
 248     this_->adj_list[u] -> prev = adj_ptr;
 249   }
 250
 251   /* update doubly chained list pointers of pointers */
 252   adj_ptr -> succ = this_->adj_list[u];
 253   adj_ptr -> prev = NULL;
 254
 255   /* update adjacency list pointer of node u */
 256   this_->adj_list[u] = adj_ptr;
 257 }
 258
 259 /* remove entry in an adj. list */
 260 static
 261 void remove_adjnode(pbqp *this_, int u, adjnode *adj_ptr)
 262 {
 263   assert(this_!= NULL);
 264   assert(u >= 0 && u <= this_->num_nodes);
 265   assert(this_->adj_list != NULL);
 266   assert(adj_ptr != NULL);
 267
 268   if (adj_ptr -> prev == NULL) {
 269     this_->adj_list[u] = adj_ptr -> succ;
 270   } else {
 271     adj_ptr -> prev -> succ = adj_ptr -> succ;
 272   }
 273
 274   if (adj_ptr -> succ != NULL) {
 275     adj_ptr -> succ -> prev = adj_ptr -> prev;
 276   }
 277
 278   if(adj_ptr->reverse != NULL) {
 279     adjnode *rev = adj_ptr->reverse;
 280     rev->reverse = NULL;
 281   }
 282
 283   if (adj_ptr -> tc_safe_regs != NULL) {
 284      free(adj_ptr -> tc_safe_regs);
 285   }
 286
 287   free(adj_ptr);
 288 }
 289
 290 /*****************************************************************************
 291  * node functions
 292  ****************************************************************************/
 293
 294 /* get degree of a node */
 295 static
 296 int get_deg(pbqp *this_,int u)
 297 {
 298   adjnode *adj_ptr;
 299   int deg = 0;
 300
 301   assert(this_ != NULL);
 302   assert(u >= 0 && u < this_->num_nodes);
 303   assert(this_->adj_list != NULL);
 304
 305   for(adj_ptr = this_ -> adj_list[u];adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
 306     deg ++;
 307   }
 308   return deg;
 309 }
 310
 311 /* reinsert node */
 312 static
 313 void reinsert_node(pbqp *this_,int u)
 314 {
 315   adjnode *adj_u,
 316           *adj_v;
 317
 318   assert(this_!= NULL);
 319   assert(u >= 0 && u <= this_->num_nodes);
 320   assert(this_->adj_list != NULL);
 321
 322   for(adj_u = this_ -> adj_list[u]; adj_u != NULL; adj_u = adj_u -> succ) {
 323     int v = adj_u -> adj;
 324     adj_v = alloc_adjnode(this_,u,adj_u->costs);
 325     insert_adjnode(this_,v,adj_v);
 326   }
 327 }
 328
 329 /* remove node */
 330 static
 331 void remove_node(pbqp *this_,int u)
 332 {
 333   adjnode *adj_ptr;
 334
 335   assert(this_!= NULL);
 336   assert(u >= 0 && u <= this_->num_nodes);
 337   assert(this_->adj_list != NULL);
 338
 339   for(adj_ptr = this_ -> adj_list[u]; adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
 340     remove_adjnode(this_,adj_ptr->adj,adj_ptr -> reverse);
 341   }
 342 }
 343
 344 /*****************************************************************************
 345  * edge functions
 346  ****************************************************************************/
 347
 348 /* insert edge to graph */
 349 /* (does not check whether edge exists in graph */
 350 static
 351 void insert_edge(pbqp *this_, int u, int v, PBQPMatrix *costs)
 352 {
 353   adjnode *adj_u,
 354           *adj_v;
 355
 356   /* create adjanceny entry for u */
 357   adj_u = alloc_adjnode(this_,v,costs);
 358   insert_adjnode(this_,u,adj_u);
 359
 360
 361   /* create adjanceny entry for v */
 362   adj_v = alloc_adjnode(this_,u,costs);
 363   insert_adjnode(this_,v,adj_v);
 364
 365   /* create link for reverse edge */
 366   adj_u -> reverse = adj_v;
 367   adj_v -> reverse = adj_u;
 368 }
 369
 370 /* delete edge */
 371 static
 372 void delete_edge(pbqp *this_,int u,int v)
 373 {
 374   adjnode *adj_ptr;
 375   adjnode *rev;
 376
 377   assert(this_ != NULL);
 378   assert( u >= 0 && u < this_->num_nodes);
 379   assert( v >= 0 && v < this_->num_nodes);
 380
 381   adj_ptr=find_adjnode(this_,u,v);
 382   assert(adj_ptr != NULL);
 383   assert(adj_ptr->reverse != NULL);
 384
 385   delete adj_ptr -> costs;
 386
 387   rev = adj_ptr->reverse;
 388   remove_adjnode(this_,u,adj_ptr);
 389   remove_adjnode(this_,v,rev);
 390 }
 391
 392 /*****************************************************************************
 393  * cost functions
 394  ****************************************************************************/
 395
 396 /* Note: Since cost(u,v) = transpose(cost(v,u)), it would be necessary to store
 397    two matrices for both edges (u,v) and (v,u). However, we only store the
 398    matrix for the case u < v. For the other case we transpose the stored matrix
 399    if required.
 400 */
 401
 402 /* add costs to cost vector of a node */
 403 void add_pbqp_nodecosts(pbqp *this_,int u, PBQPVector *costs)
 404 {
 405   assert(this_ != NULL);
 406   assert(costs != NULL);
 407   assert(u >= 0 && u <= this_->num_nodes);
 408
 409   if (!this_->node_costs[u]) {
 410     this_->node_costs[u] = new PBQPVector(*costs);
 411   } else {
 412     *this_->node_costs[u] += *costs;
 413   }
 414 }
 415
 416 /* get cost matrix ptr */
 417 static
 418 PBQPMatrix *get_costmatrix_ptr(pbqp *this_, int u, int v)
 419 {
 420   adjnode *adj_ptr;
 421   PBQPMatrix *m = NULL;
 422
 423   assert (this_ != NULL);
 424   assert (u >= 0 && u < this_->num_nodes);
 425   assert (v >= 0 && v < this_->num_nodes);
 426
 427   adj_ptr = find_adjnode(this_,u,v);
 428
 429   if (adj_ptr != NULL) {
 430     m = adj_ptr -> costs;
 431   }
 432
 433   return m;
 434 }
 435
 436 /* get cost matrix ptr */
 437 /* Note: only the pointer is returned for
 438    cost(u,v), if u < v.
 439 */
 440 static
 441 PBQPMatrix *pbqp_get_costmatrix(pbqp *this_, int u, int v)
 442 {
 443   adjnode *adj_ptr = find_adjnode(this_,u,v);
 444
 445   if (adj_ptr != NULL) {
 446     if ( u < v) {
 447       return new PBQPMatrix(*adj_ptr->costs);
 448     } else {
 449       return new PBQPMatrix(adj_ptr->costs->transpose());
 450     }
 451   } else {
 452     return NULL;
 453   }
 454 }
 455
 456 /* add costs to cost matrix of an edge */
 457 void add_pbqp_edgecosts(pbqp *this_,int u,int v, PBQPMatrix *costs)
 458 {
 459   PBQPMatrix *adj_costs;
 460
 461   assert(this_!= NULL);
 462   assert(costs != NULL);
 463   assert(u >= 0 && u <= this_->num_nodes);
 464   assert(v >= 0 && v <= this_->num_nodes);
 465
 466   /* does the edge u-v exists ? */
 467   if (u == v) {
 468     PBQPVector *diag = new PBQPVector(costs->diagonalize());
 469     add_pbqp_nodecosts(this_,v,diag);
 470     delete diag;
 471   } else if ((adj_costs = get_costmatrix_ptr(this_,u,v))!=NULL) {
 472     if ( u < v) {
 473       *adj_costs += *costs;
 474     } else {
 475       *adj_costs += costs->transpose();
 476     }
 477   } else {
 478     adj_costs = new PBQPMatrix((u < v) ? *costs : costs->transpose());
 479     insert_edge(this_,u,v,adj_costs);
 480   }
 481 }
 482
 483 /* remove bucket from bucket list */
 484 static
 485 void pbqp_remove_bucket(pbqp *this_, bucketnode *bucket)
 486 {
 487   int u = bucket->u;
 488
 489   assert(this_ != NULL);
 490   assert(u >= 0 && u < this_->num_nodes);
 491   assert(this_->bucket_list != NULL);
 492   assert(this_->bucket_ptr[u] != NULL);
 493
 494   /* update predecessor node in bucket list
 495      (if no preceeding bucket exists, then
 496      the bucket_list pointer needs to be
 497      updated.)
 498   */
 499   if (bucket->prev != NULL) {
 500     bucket->prev-> succ = bucket->succ;
 501   } else {
 502     this_->bucket_list[this_->node_deg[u]] = bucket -> succ;
 503   }
 504
 505   /* update successor node in bucket list */
 506   if (bucket->succ != NULL) {
 507     bucket->succ-> prev = bucket->prev;
 508   }
 509 }
 510
 511 /**********************************************************************************
 512  * pop functions
 513  **********************************************************************************/
 514
 515 /* pop node of given degree */
 516 static
 517 int pop_node(pbqp *this_,int deg)
 518 {
 519   bucketnode *bucket;
 520   int u;
 521
 522   assert(this_ != NULL);
 523   assert(deg >= 0 && deg <= this_->max_deg);
 524   assert(this_->bucket_list != NULL);
 525
 526   /* get first bucket of bucket list */
 527   bucket = this_->bucket_list[deg];
 528   assert(bucket != NULL);
 529
 530   /* remove bucket */
 531   pbqp_remove_bucket(this_,bucket);
 532   u = bucket->u;
 533   free(bucket);
 534   return u;
 535 }
 536
 537 /**********************************************************************************
 538  * reorder functions
 539  **********************************************************************************/
 540
 541 /* add bucket to bucketlist */
 542 static
 543 void add_to_bucketlist(pbqp *this_,bucketnode *bucket, int deg)
 544 {
 545   bucketnode *old_head;
 546
 547   assert(bucket != NULL);
 548   assert(this_ != NULL);
 549   assert(deg >= 0 && deg <= this_->max_deg);
 550   assert(this_->bucket_list != NULL);
 551
 552   /* store node degree (for re-ordering purposes)*/
 553   this_->node_deg[bucket->u] = deg;
 554
 555   /* put bucket to front of doubly chained list */
 556   old_head = this_->bucket_list[deg];
 557   bucket -> prev = NULL;
 558   bucket -> succ = old_head;
 559   this_ -> bucket_list[deg] = bucket;
 560   if (bucket -> succ != NULL ) {
 561     assert ( old_head -> prev == NULL);
 562     old_head -> prev = bucket;
 563   }
 564 }
 565
 566
 567 /* reorder node in bucket list according to
 568    current node degree */
 569 static
 570 void reorder_node(pbqp *this_, int u)
 571 {
 572   int deg;
 573
 574   assert(this_ != NULL);
 575   assert(u>= 0 && u < this_->num_nodes);
 576   assert(this_->bucket_list != NULL);
 577   assert(this_->bucket_ptr[u] != NULL);
 578
 579   /* get current node degree */
 580   deg = get_deg(this_,u);
 581
 582   /* remove bucket from old bucket list only
 583      if degree of node has changed. */
 584   if (deg != this_->node_deg[u]) {
 585     pbqp_remove_bucket(this_,this_->bucket_ptr[u]);
 586     add_to_bucketlist(this_,this_->bucket_ptr[u],deg);
 587   }
 588 }
 589
 590 /* reorder adj. nodes of a node */
 591 static
 592 void reorder_adjnodes(pbqp *this_,int u)
 593 {
 594   adjnode *adj_ptr;
 595
 596   assert(this_!= NULL);
 597   assert(u >= 0 && u <= this_->num_nodes);
 598   assert(this_->adj_list != NULL);
 599
 600   for(adj_ptr = this_ -> adj_list[u]; adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
 601     reorder_node(this_,adj_ptr->adj);
 602   }
 603 }
 604
 605 /**********************************************************************************
 606  * creation functions
 607  **********************************************************************************/
 608
 609 /* create new bucket entry */
 610 /* consistency of the bucket list is not checked! */
 611 static
 612 void create_bucket(pbqp *this_,int u,int deg)
 613 {
 614   bucketnode *bucket;
 615
 616   assert(this_ != NULL);
 617   assert(u >= 0 && u < this_->num_nodes);
 618   assert(this_->bucket_list != NULL);
 619
 620   bucket = (bucketnode *)malloc(sizeof(bucketnode));
 621   assert(bucket != NULL);
 622
 623   bucket -> u = u;
 624   this_->bucket_ptr[u] = bucket;
 625
 626   add_to_bucketlist(this_,bucket,deg);
 627 }
 628
 629 /* create bucket list */
 630 static
 631 void create_bucketlist(pbqp *this_)
 632 {
 633   int u;
 634   int max_deg;
 635   int deg;
 636
 637   assert(this_ != NULL);
 638   assert(this_->bucket_list == NULL);
 639
 640   /* determine max. degree of the nodes */
 641   max_deg = 2;  /* at least of degree two! */
 642   for(u=0;u<this_->num_nodes;u++) {
 643     deg = this_->node_deg[u] = get_deg(this_,u);
 644     if (deg > max_deg) {
 645       max_deg = deg;
 646     }
 647   }
 648   this_->max_deg = max_deg;
 649
 650   /* allocate bucket list */
 651   this_ -> bucket_list = (bucketnode **)malloc(sizeof(bucketnode *)*(max_deg + 1));
 652   memset(this_->bucket_list,0,sizeof(bucketnode *)*(max_deg + 1));
 653   assert(this_->bucket_list != NULL);
 654
 655   /* insert nodes to the list */
 656   for(u=0;u<this_->num_nodes;u++) {
 657     create_bucket(this_,u,this_->node_deg[u]);
 658   }
 659 }
 660
 661 /*****************************************************************************
 662  * PBQP simplification for trivial nodes
 663  ****************************************************************************/
 664
 665 /* remove trivial node with cost vector length of one */
 666 static
 667 void disconnect_trivialnode(pbqp *this_,int u)
 668 {
 669   int v;
 670   adjnode *adj_ptr,
 671           *next;
 672   PBQPMatrix *c_uv;
 673   PBQPVector *c_v;
 674
 675   assert(this_ != NULL);
 676   assert(this_->node_costs != NULL);
 677   assert(u >= 0 && u < this_ -> num_nodes);
 678   assert(this_->node_costs[u]->getLength() == 1);
 679
 680   /* add edge costs to node costs of adj. nodes */
 681   for(adj_ptr = this_->adj_list[u]; adj_ptr != NULL; adj_ptr = next){
 682     next = adj_ptr -> succ;
 683     v = adj_ptr -> adj;
 684     assert(v >= 0 && v < this_ -> num_nodes);
 685
 686     /* convert matrix to cost vector offset for adj. node */
 687     c_uv = pbqp_get_costmatrix(this_,u,v);
 688     c_v = new PBQPVector(c_uv->getRowAsVector(0));
 689     *this_->node_costs[v] += *c_v;
 690
 691     /* delete edge & free vec/mat */
 692     delete c_v;
 693     delete c_uv;
 694     delete_edge(this_,u,v);
 695   }
 696 }
 697
 698 /* find all trivial nodes and disconnect them */
 699 static
 700 void eliminate_trivial_nodes(pbqp *this_)
 701 {
 702    int u;
 703
 704    assert(this_ != NULL);
 705    assert(this_ -> node_costs != NULL);
 706
 707    for(u=0;u < this_ -> num_nodes; u++) {
 708      if (this_->node_costs[u]->getLength() == 1) {
 709        disconnect_trivialnode(this_,u);
 710      }
 711    }
 712 }
 713
 714 /*****************************************************************************
 715  * Normal form for PBQP
 716  ****************************************************************************/
 717
 718 /* simplify a cost matrix. If the matrix
 719    is independent, then simplify_matrix
 720    returns true - otherwise false. In
 721    vectors u and v the offset values of
 722    the decomposition are stored.
 723 */
 724
 725 static
 726 bool normalize_matrix(PBQPMatrix *m, PBQPVector *u, PBQPVector *v)
 727 {
 728   assert( m != NULL);
 729   assert( u != NULL);
 730   assert( v != NULL);
 731   assert( u->getLength() > 0);
 732   assert( v->getLength() > 0);
 733
 734   assert(m->getRows() == u->getLength());
 735   assert(m->getCols() == v->getLength());
 736
 737   /* determine u vector */
 738   for(unsigned r = 0; r < m->getRows(); ++r) {
 739     PBQPNum min = m->getRowMin(r);
 740     (*u)[r] += min;
 741     if (!isInf(min)) {
 742       m->subFromRow(r, min);
 743     } else {
 744       m->setRow(r, 0);
 745     }
 746   }
 747
 748   /* determine v vector */
 749   for(unsigned c = 0; c < m->getCols(); ++c) {
 750     PBQPNum min = m->getColMin(c);
 751     (*v)[c] += min;
 752     if (!isInf(min)) {
 753       m->subFromCol(c, min);
 754     } else {
 755       m->setCol(c, 0);
 756     }
 757   }
 758
 759   /* determine whether matrix is
 760      independent or not.
 761     */
 762   return m->isZero();
 763 }
 764
 765 /* simplify single edge */
 766 static
 767 void simplify_edge(pbqp *this_,int u,int v)
 768 {
 769   PBQPMatrix *costs;
 770   bool is_zero;
 771
 772   assert (this_ != NULL);
 773   assert (u >= 0 && u <this_->num_nodes);
 774   assert (v >= 0 && v <this_->num_nodes);
 775   assert (u != v);
 776
 777   /* swap u and v  if u > v in order to avoid un-necessary
 778      tranpositions of the cost matrix */
 779
 780   if (u > v) {
 781     int swap = u;
 782     u = v;
 783     v = swap;
 784   }
 785
 786   /* get cost matrix and simplify it */
 787   costs = get_costmatrix_ptr(this_,u,v);
 788   is_zero=normalize_matrix(costs,this_->node_costs[u],this_->node_costs[v]);
 789
 790   /* delete edge */
 791   if(is_zero){
 792     delete_edge(this_,u,v);
 793     this_->changed = true;
 794   }
 795 }
 796
 797 /* normalize cost matrices and remove
 798    edges in PBQP if they ary independent,
 799    i.e. can be decomposed into two
 800    cost vectors.
 801 */
 802 static
 803 void eliminate_independent_edges(pbqp *this_)
 804 {
 805   int u,v;
 806   adjnode *adj_ptr,*next;
 807
 808   assert(this_ != NULL);
 809   assert(this_ -> adj_list != NULL);
 810
 811   this_->changed = false;
 812   for(u=0;u < this_->num_nodes;u++) {
 813     for (adj_ptr = this_ -> adj_list[u]; adj_ptr != NULL; adj_ptr = next) {
 814       next = adj_ptr -> succ;
 815       v = adj_ptr -> adj;
 816       assert(v >= 0 && v < this_->num_nodes);
 817       if (u < v) {
 818         simplify_edge(this_,u,v);
 819       }
 820     }
 821   }
 822 }
 823
 824
 825 /*****************************************************************************
 826  * PBQP reduction rules
 827  ****************************************************************************/
 828
 829 /* RI reduction
 830    This reduction rule is applied for nodes
 831    of degree one. */
 832
 833 static
 834 void apply_RI(pbqp *this_,int x)
 835 {
 836   int y;
 837   unsigned xlen,
 838            ylen;
 839   PBQPMatrix *c_yx;
 840   PBQPVector *c_x, *delta;
 841
 842   assert(this_ != NULL);
 843   assert(x >= 0 && x < this_->num_nodes);
 844   assert(this_ -> adj_list[x] != NULL);
 845   assert(this_ -> adj_list[x] -> succ == NULL);
 846
 847   /* get adjacence matrix */
 848   y = this_ -> adj_list[x] -> adj;
 849   assert(y >= 0 && y < this_->num_nodes);
 850
 851   /* determine length of cost vectors for node x and y */
 852   xlen = this_ -> node_costs[x]->getLength();
 853   ylen = this_ -> node_costs[y]->getLength();
 854
 855   /* get cost vector c_x and matrix c_yx */
 856   c_x = this_ -> node_costs[x];
 857   c_yx = pbqp_get_costmatrix(this_,y,x);
 858   assert (c_yx != NULL);
 859
 860
 861   /* allocate delta vector */
 862   delta = new PBQPVector(ylen);
 863
 864   /* compute delta vector */
 865   for(unsigned i = 0; i < ylen; ++i) {
 866     PBQPNum min =  (*c_yx)[i][0] + (*c_x)[0];
 867     for(unsigned j = 1; j < xlen; ++j) {
 868       PBQPNum c =  (*c_yx)[i][j] + (*c_x)[j];
 869       if ( c < min )
 870          min = c;
 871     }
 872     (*delta)[i] = min;
 873   }
 874
 875   /* add delta vector */
 876   *this_ -> node_costs[y] += *delta;
 877
 878   /* delete node x */
 879   remove_node(this_,x);
 880
 881   /* reorder adj. nodes of node x */
 882   reorder_adjnodes(this_,x);
 883
 884   /* push node x on stack */
 885   assert(this_ -> stack_ptr < this_ -> num_nodes);
 886   this_->stack[this_ -> stack_ptr++] = x;
 887
 888   /* free vec/mat */
 889   delete c_yx;
 890   delete delta;
 891
 892   /* increment counter for number statistic */
 893   this_->num_ri++;
 894 }
 895
 896 /* RII reduction
 897    This reduction rule is applied for nodes
 898    of degree two. */
 899
 900 static
 901 void apply_RII(pbqp *this_,int x)
 902 {
 903   int y,z;
 904   unsigned xlen,ylen,zlen;
 905   adjnode *adj_yz;
 906
 907   PBQPMatrix *c_yx, *c_zx;
 908   PBQPVector *cx;
 909   PBQPMatrix *delta;
 910
 911   assert(this_ != NULL);
 912   assert(x >= 0 && x < this_->num_nodes);
 913   assert(this_ -> adj_list[x] != NULL);
 914   assert(this_ -> adj_list[x] -> succ != NULL);
 915   assert(this_ -> adj_list[x] -> succ -> succ == NULL);
 916
 917   /* get adjacence matrix */
 918   y = this_ -> adj_list[x] -> adj;
 919   z = this_ -> adj_list[x] -> succ -> adj;
 920   assert(y >= 0 && y < this_->num_nodes);
 921   assert(z >= 0 && z < this_->num_nodes);
 922
 923   /* determine length of cost vectors for node x and y */
 924   xlen = this_ -> node_costs[x]->getLength();
 925   ylen = this_ -> node_costs[y]->getLength();
 926   zlen = this_ -> node_costs[z]->getLength();
 927
 928   /* get cost vector c_x and matrix c_yx */
 929   cx = this_ -> node_costs[x];
 930   c_yx = pbqp_get_costmatrix(this_,y,x);
 931   c_zx = pbqp_get_costmatrix(this_,z,x);
 932   assert(c_yx != NULL);
 933   assert(c_zx != NULL);
 934
 935   /* Colour Heuristic */
 936   if ( (adj_yz = find_adjnode(this_,y,z)) != NULL) {
 937     adj_yz->tc_valid = false;
 938     adj_yz->reverse->tc_valid = false;
 939   }
 940
 941   /* allocate delta matrix */
 942   delta = new PBQPMatrix(ylen, zlen);
 943
 944   /* compute delta matrix */
 945   for(unsigned i=0;i<ylen;i++) {
 946     for(unsigned j=0;j<zlen;j++) {
 947       PBQPNum min = (*c_yx)[i][0] + (*c_zx)[j][0] + (*cx)[0];
 948       for(unsigned k=1;k<xlen;k++) {
 949         PBQPNum c = (*c_yx)[i][k] + (*c_zx)[j][k] + (*cx)[k];
 950         if ( c < min ) {
 951           min = c;
 952         }
 953       }
 954       (*delta)[i][j] = min;
 955     }
 956   }
 957
 958   /* add delta matrix */
 959   add_pbqp_edgecosts(this_,y,z,delta);
 960
 961   /* delete node x */
 962   remove_node(this_,x);
 963
 964   /* simplify cost matrix c_yz */
 965   simplify_edge(this_,y,z);
 966
 967   /* reorder adj. nodes */
 968   reorder_adjnodes(this_,x);
 969
 970   /* push node x on stack */
 971   assert(this_ -> stack_ptr < this_ -> num_nodes);
 972   this_->stack[this_ -> stack_ptr++] = x;
 973
 974   /* free vec/mat */
 975   delete c_yx;
 976   delete c_zx;
 977   delete delta;
 978
 979   /* increment counter for number statistic */
 980   this_->num_rii++;
 981
 982 }
 983
 984 /* RN reduction */
 985 static
 986 void apply_RN(pbqp *this_,int x)
 987 {
 988   unsigned xlen;
 989
 990   assert(this_ != NULL);
 991   assert(x >= 0 && x < this_->num_nodes);
 992   assert(this_ -> node_costs[x] != NULL);
 993
 994   xlen = this_ -> node_costs[x] -> getLength();
 995
 996   /* after application of RN rule no optimality
 997      can be guaranteed! */
 998   this_ -> optimal = false;
 999
1000   /* push node x on stack */
1001   assert(this_ -> stack_ptr < this_ -> num_nodes);
1002   this_->stack[this_ -> stack_ptr++] = x;
1003
1004   /* delete node x */
1005   remove_node(this_,x);
1006
1007   /* reorder adj. nodes of node x */
1008   reorder_adjnodes(this_,x);
1009
1010   /* increment counter for number statistic */
1011   this_->num_rn++;
1012 }
1013
1014
1015 static
1016 void compute_tc_info(pbqp *this_, adjnode *p)
1017 {
1018    adjnode *r;
1019    PBQPMatrix *m;
1020    int x,y;
1021    PBQPVector *c_x, *c_y;
1022    int *row_inf_counts;
1023
1024    assert(p->reverse != NULL);
1025
1026    /* set flags */
1027    r = p->reverse;
1028    p->tc_valid = true;
1029    r->tc_valid = true;
1030
1031    /* get edge */
1032    x = r->adj;
1033    y = p->adj;
1034
1035    /* get cost vectors */
1036    c_x = this_ -> node_costs[x];
1037    c_y = this_ -> node_costs[y];
1038
1039    /* get cost matrix */
1040    m = pbqp_get_costmatrix(this_, x, y);
1041
1042
1043    /* allocate allowed set for edge (x,y) and (y,x) */
1044    if (p->tc_safe_regs == NULL) {
1045      p->tc_safe_regs = (int *) malloc(sizeof(int) * c_x->getLength());
1046    }
1047
1048    if (r->tc_safe_regs == NULL ) {
1049      r->tc_safe_regs = (int *) malloc(sizeof(int) * c_y->getLength());
1050    }
1051
1052    p->tc_impact = r->tc_impact = 0;
1053
1054    row_inf_counts = (int *) alloca(sizeof(int) * c_x->getLength());
1055
1056    /* init arrays */
1057    p->tc_safe_regs[0] = 0;
1058    row_inf_counts[0] = 0;
1059    for(unsigned i = 1; i < c_x->getLength(); ++i){
1060      p->tc_safe_regs[i] = 1;
1061      row_inf_counts[i] = 0;
1062    }
1063
1064    r->tc_safe_regs[0] = 0;
1065    for(unsigned j = 1; j < c_y->getLength(); ++j){
1066      r->tc_safe_regs[j] = 1;
1067    }
1068
1069    for(unsigned j = 0; j < c_y->getLength(); ++j) {
1070       int col_inf_counts = 0;
1071       for (unsigned i = 0; i < c_x->getLength(); ++i) {
1072          if (isInf((*m)[i][j])) {
1073               ++col_inf_counts;
1074               ++row_inf_counts[i];
1075
1076               p->tc_safe_regs[i] = 0;
1077               r->tc_safe_regs[j] = 0;
1078          }
1079       }
1080       if (col_inf_counts > p->tc_impact) {
1081            p->tc_impact = col_inf_counts;
1082       }
1083    }
1084
1085    for(unsigned i = 0; i < c_x->getLength(); ++i){
1086      if (row_inf_counts[i] > r->tc_impact)
1087      {
1088            r->tc_impact = row_inf_counts[i];
1089      }
1090    }
1091
1092    delete m;
1093 }
1094
1095 /*
1096  * Checks whether node x can be locally coloured.
1097  */
1098 static
1099 int is_colorable(pbqp *this_,int x)
1100 {
1101   adjnode *adj_ptr;
1102   PBQPVector *c_x;
1103   int result = 1;
1104   int *allowed;
1105   int num_allowed = 0;
1106   unsigned total_impact = 0;
1107
1108   assert(this_ != NULL);
1109   assert(x >= 0 && x < this_->num_nodes);
1110   assert(this_ -> node_costs[x] != NULL);
1111
1112   c_x = this_ -> node_costs[x];
1113
1114   /* allocate allowed set */
1115   allowed = (int *)malloc(sizeof(int) * c_x->getLength());
1116   for(unsigned i = 0; i < c_x->getLength(); ++i){
1117     if (!isInf((*c_x)[i]) && i > 0) {
1118       allowed[i] = 1;
1119       ++num_allowed;
1120     } else {
1121       allowed[i] = 0;
1122     }
1123   }
1124
1125   /* determine local minimum */
1126   for(adj_ptr=this_->adj_list[x] ;adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
1127       if (!adj_ptr -> tc_valid) {
1128           compute_tc_info(this_, adj_ptr);
1129       }
1130
1131       total_impact += adj_ptr->tc_impact;
1132
1133       if (num_allowed > 0) {
1134           for (unsigned i = 1; i < c_x->getLength(); ++i){
1135             if (allowed[i]){
1136               if (!adj_ptr->tc_safe_regs[i]){
1137                 allowed[i] = 0;
1138                 --num_allowed;
1139                 if (num_allowed == 0)
1140                     break;
1141               }
1142             }
1143           }
1144       }
1145
1146       if ( total_impact >= c_x->getLength() - 1 && num_allowed == 0 ) {
1147          result = 0;
1148          break;
1149       }
1150   }
1151   free(allowed);
1152
1153   return result;
1154 }
1155
1156 /* use briggs heuristic
1157   note: this_ is not a general heuristic. it only is useful for
1158   interference graphs.
1159  */
1160 int pop_colorablenode(pbqp *this_)
1161 {
1162   int deg;
1163   bucketnode *min_bucket=NULL;
1164   PBQPNum min = std::numeric_limits<PBQPNum>::infinity();
1165
1166   /* select node where the number of colors is less than the node degree */
1167   for(deg=this_->max_deg;deg > 2;deg--) {
1168     bucketnode *bucket;
1169     for(bucket=this_->bucket_list[deg];bucket!= NULL;bucket = bucket -> succ) {
1170       int u = bucket->u;
1171       if (is_colorable(this_,u)) {
1172         pbqp_remove_bucket(this_,bucket);
1173         this_->num_rn_special++;
1174         free(bucket);
1175         return u;
1176       }
1177     }
1178   }
1179
1180   /* select node with minimal ratio between average node costs and degree of node */
1181   for(deg=this_->max_deg;deg >2; deg--) {
1182     bucketnode *bucket;
1183     for(bucket=this_->bucket_list[deg];bucket!= NULL;bucket = bucket -> succ) {
1184       PBQPNum h;
1185       int u;
1186
1187       u = bucket->u;
1188       assert(u>=0 && u < this_->num_nodes);
1189       h = (*this_->node_costs[u])[0] / (PBQPNum) deg;
1190       if (h < min) {
1191         min_bucket = bucket;
1192         min = h;
1193       }
1194     }
1195   }
1196
1197   /* return node and free bucket */
1198   if (min_bucket != NULL) {
1199     int u;
1200
1201     pbqp_remove_bucket(this_,min_bucket);
1202     u = min_bucket->u;
1203     free(min_bucket);
1204     return u;
1205   } else {
1206     return -1;
1207   }
1208 }
1209
1210
1211 /*****************************************************************************
1212  * PBQP graph parsing
1213  ****************************************************************************/
1214
1215 /* reduce pbqp problem (first phase) */
1216 static
1217 void reduce_pbqp(pbqp *this_)
1218 {
1219   int u;
1220
1221   assert(this_ != NULL);
1222   assert(this_->bucket_list != NULL);
1223
1224   for(;;){
1225
1226     if (this_->bucket_list[1] != NULL) {
1227       u = pop_node(this_,1);
1228       apply_RI(this_,u);
1229     } else if (this_->bucket_list[2] != NULL) {
1230       u = pop_node(this_,2);
1231       apply_RII(this_,u);
1232     } else if ((u = pop_colorablenode(this_)) != -1) {
1233       apply_RN(this_,u);
1234     } else {
1235       break;
1236     }
1237   }
1238 }
1239
1240 /*****************************************************************************
1241  * PBQP back propagation
1242  ****************************************************************************/
1243
1244 /* determine solution of a reduced node. Either
1245    RI or RII was applied for this_ node. */
1246 static
1247 void determine_solution(pbqp *this_,int x)
1248 {
1249   PBQPVector *v = new PBQPVector(*this_ -> node_costs[x]);
1250   adjnode *adj_ptr;
1251
1252   assert(this_ != NULL);
1253   assert(x >= 0 && x < this_->num_nodes);
1254   assert(this_ -> adj_list != NULL);
1255   assert(this_ -> solution != NULL);
1256
1257   for(adj_ptr=this_->adj_list[x] ;adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
1258     int y = adj_ptr -> adj;
1259     int y_sol = this_ -> solution[y];
1260
1261     PBQPMatrix *c_yx = pbqp_get_costmatrix(this_,y,x);
1262     assert(y_sol >= 0 && y_sol < (int)this_->node_costs[y]->getLength());
1263     (*v) += c_yx->getRowAsVector(y_sol);
1264     delete c_yx;
1265   }
1266   this_ -> solution[x] = v->minIndex();
1267
1268   delete v;
1269 }
1270
1271 /* back popagation phase of PBQP */
1272 static
1273 void back_propagate(pbqp *this_)
1274 {
1275    int i;
1276
1277    assert(this_ != NULL);
1278    assert(this_->stack != NULL);
1279    assert(this_->stack_ptr < this_->num_nodes);
1280
1281    for(i=this_ -> stack_ptr-1;i>=0;i--) {
1282       int x = this_ -> stack[i];
1283       assert( x >= 0 && x < this_ -> num_nodes);
1284       reinsert_node(this_,x);
1285       determine_solution(this_,x);
1286    }
1287 }
1288
1289 /* solve trivial nodes of degree zero */
1290 static
1291 void determine_trivialsolution(pbqp *this_)
1292 {
1293   int u;
1294   PBQPNum delta;
1295
1296   assert( this_ != NULL);
1297   assert( this_ -> bucket_list != NULL);
1298
1299   /* determine trivial solution */
1300   while (this_->bucket_list[0] != NULL) {
1301     u = pop_node(this_,0);
1302
1303     assert( u >= 0 && u < this_ -> num_nodes);
1304
1305     this_->solution[u] = this_->node_costs[u]->minIndex();
1306     delta = (*this_->node_costs[u])[this_->solution[u]];
1307     this_->min = this_->min + delta;
1308
1309     /* increment counter for number statistic */
1310     this_->num_r0++;
1311   }
1312 }
1313
1314 /*****************************************************************************
1315  * debug facilities
1316  ****************************************************************************/
1317 static
1318 void check_pbqp(pbqp *this_)
1319 {
1320   int u,v;
1321   PBQPMatrix *costs;
1322   adjnode *adj_ptr;
1323
1324   assert( this_ != NULL);
1325
1326   for(u=0;u< this_->num_nodes; u++) {
1327     assert (this_ -> node_costs[u] != NULL);
1328     for(adj_ptr = this_ -> adj_list[u];adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
1329       v = adj_ptr -> adj;
1330       assert( v>= 0 && v < this_->num_nodes);
1331       if (u < v ) {
1332         costs = adj_ptr -> costs;
1333         assert( costs->getRows() == this_->node_costs[u]->getLength() &&
1334                 costs->getCols() == this_->node_costs[v]->getLength());
1335       }
1336     }
1337   }
1338 }
1339
1340 /*****************************************************************************
1341  * PBQP solve routines
1342  ****************************************************************************/
1343
1344 /* solve PBQP problem */
1345 void solve_pbqp(pbqp *this_)
1346 {
1347   assert(this_ != NULL);
1348   assert(!this_->solved);
1349
1350   /* check vector & matrix dimensions */
1351   check_pbqp(this_);
1352
1353   /* simplify PBQP problem */
1354
1355   /* eliminate trivial nodes, i.e.
1356      nodes with cost vectors of length one.  */
1357   eliminate_trivial_nodes(this_);
1358
1359   /* eliminate edges with independent
1360      cost matrices and normalize matrices */
1361   eliminate_independent_edges(this_);
1362
1363   /* create bucket list for graph parsing */
1364   create_bucketlist(this_);
1365
1366   /* reduce phase */
1367   reduce_pbqp(this_);
1368
1369   /* solve trivial nodes */
1370   determine_trivialsolution(this_);
1371
1372   /* back propagation phase */
1373   back_propagate(this_);
1374
1375   this_->solved = true;
1376 }
1377
1378 /* get solution of a node */
1379 int get_pbqp_solution(pbqp *this_,int x)
1380 {
1381   assert(this_ != NULL);
1382   assert(this_->solution != NULL);
1383   assert(this_ -> solved);
1384
1385   return this_->solution[x];
1386 }
1387
1388 /* is solution optimal? */
1389 bool is_pbqp_optimal(pbqp *this_)
1390 {
1391   assert(this_ -> solved);
1392   return this_->optimal;
1393 }
1394
1395 }
1396
1397 /* end of pbqp.c */