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