src/org/ibex/util/LinearProgramming.java

   1 package org.ibex.util;
   2 import java.io.* ;
   3 import java.util.* ;
   4
   5 public class LinearProgramming {
   6
   7     public final static short FAIL = -1;
   8
   9     public final static short NULL = 0;
  10     public final static short FALSE = 0;
  11     public final static short TRUE = 1;
  12
  13     public final static short DEFNUMINV = 50;
  14
  15     /* solve status values */
  16     public final static short OPTIMAL = 0;
  17     public final static short MILP_FAIL = 1;
  18     public final static short INFEASIBLE = 2;
  19     public final static short UNBOUNDED = 3;
  20     public final static short FAILURE = 4;
  21     public final static short RUNNING = 5;
  22
  23     /* lag_solve extra status values */
  24     public final static short FEAS_FOUND = 6;
  25     public final static short NO_FEAS_FOUND = 7;
  26     public final static short BREAK_BB = 8;
  27
  28     public final static short FIRST_NI =        0;
  29     public final static short RAND_NI = 1;
  30
  31     public final static short LE = 0;
  32     public final static short EQ = 1;
  33     public final static short GE = 2;
  34     public final static short OF = 3;
  35
  36     public final static short MAX_WARN_COUNT = 20;
  37
  38     public final static float DEF_INFINITE = (float)1e24; /* limit for dynamic range */
  39     public final static float DEF_EPSB = (float)5.01e-7; /* for rounding RHS values to 0 determine
  40                                                infeasibility basis */
  41     public final static float DEF_EPSEL = (float)1e-8; /* for rounding other values (vectors) to 0 */
  42     public final static float DEF_EPSD  = (float)1e-6; /* for rounding reduced costs to zero */
  43     public final static float DEF_EPSILON = (float)1e-3; /* to determine if a float value is integer */
  44
  45     public final static float PREJ = (float)1e-3;  /* pivot reject (try others first) */
  46
  47     public final static int HASHSIZE = 10007; /* prime number is better, MB */
  48     public final static int ETA_START_SIZE = 10000; /* start size of array Eta. Realloced if needed */
  49     public final static String STD_ROW_NAME_PREFIX = "r_";
  50
  51     static class Ref {
  52         float value;
  53         public Ref(float v) { value = v; }
  54     }
  55
  56     public static class Simplex {
  57         /* Globals */
  58         Problem   Lp; /* pointer to active problem */
  59         int     Rows;
  60         int     columns;
  61         int     Sum;
  62         int     Non_zeros;
  63         int     Level;
  64         Matrix[]  Mat;
  65         int[]     Col_no;
  66         int[]     Col_end;
  67         int[]     Row_end;
  68         float[]    Orig_rh;
  69         float[]    Rh;
  70         float[]    Rhs;
  71         float[]    Orig_upbo;
  72         float[]    Orig_lowbo;
  73         float[]    Upbo;
  74         float[]    Lowbo;
  75         int[]     Bas;
  76         short[]   Basis;
  77         short[]   Lower;
  78         int     Eta_alloc;
  79         int     Eta_size;
  80         float[]    Eta_value;
  81         int[]     Eta_row_nr;
  82         int[]     Eta_col_end;
  83         int     Num_inv;
  84         float[]    Solution;
  85         float[]    Best_solution;
  86         float    Infinite;
  87         float    Epsilon;
  88         float    Epsb;
  89         float    Epsd;
  90         float    Epsel;
  91
  92         float   TREJ;
  93         float   TINV;
  94
  95         short   Maximise;
  96         short   Floor_first;
  97         float    Extrad;
  98
  99         int     Warn_count; /* used in CHECK version of rounding macro */
 100
 101         public void set_mat(Problem lp, int Row, int column, float Value)
 102         {
 103             int elmnr, lastelm, i;
 104             //  System.err.println("lp.mat.length = " + lp.mat.length);
 105
 106             if(Row > lp.rows || Row < 0)
 107                 System.err.print("Row out of range");
 108             if(column > lp.columns || column < 1)
 109                 System.err.print("column out of range");
 110             if(lp.scaling_used != FALSE)
 111                 Value *= lp.scale[Row] * lp.scale[lp.rows + column];
 112
 113             if(true /*abs(Value) > lp.epsilon*/)
 114                 {
 115                     if (lp.basis[column] == TRUE && Row > 0)
 116                         lp.basis_valid = FALSE;
 117                     lp.eta_valid = FALSE;
 118                     elmnr = lp.col_end[column-1];
 119                     while((elmnr < lp.col_end[column]) ?
 120                           (lp.mat[elmnr].row_nr != Row) : false)
 121                         elmnr++;
 122
 123                     if((elmnr != lp.col_end[column]) ?
 124                        (lp.mat[elmnr].row_nr == Row) : false )
 125                         if (lp.scaling_used != FALSE)
 126                             {
 127                                 if (lp.ch_sign[Row] != FALSE)
 128                                     lp.mat[elmnr].value =
 129                                         -Value * lp.scale[Row] * lp.scale[column];
 130                                 else
 131                                     lp.mat[elmnr].value =
 132                                         Value * lp.scale[Row] * lp.scale[column];
 133                             }
 134                         else
 135                             {
 136                                 if (lp.ch_sign[Row] != FALSE)
 137                                     lp.mat[elmnr].value = -Value;
 138                                 else
 139                                     lp.mat[elmnr].value = Value;
 140                             }
 141                     else
 142                         {
 143                             /* check if more space is needed for matrix */
 144                             if (lp.non_zeros + 1 > lp.mat_alloc) throw new Error("not enough mat space; this should not happen");
 145
 146                             /* Shift the matrix */
 147                             lastelm=lp.non_zeros;
 148                             for(i = lastelm; i > elmnr ; i--)
 149                                 lp.mat[i]=lp.mat[i-1];
 150                             for(i = column; i <= lp.columns; i++)
 151                                 lp.col_end[i]++;
 152
 153                             /* Set new element */
 154                             lp.mat[elmnr].row_nr=Row;
 155
 156                             if (lp.scaling_used != FALSE)
 157                                 {
 158                                     if (lp.ch_sign[Row] != FALSE)
 159                                         lp.mat[elmnr].value=-Value*lp.scale[Row]*lp.scale[column];
 160                                     else
 161                                         lp.mat[elmnr].value=Value*lp.scale[Row]*lp.scale[column];
 162                                 }
 163                             else
 164                                 {
 165                                     if (lp.ch_sign[Row] != FALSE)
 166                                         lp.mat[elmnr].value=-Value;
 167                                     else
 168                                         lp.mat[elmnr].value=Value;
 169                                 }
 170
 171                             lp.row_end_valid=FALSE;
 172
 173                             lp.non_zeros++;
 174                             if (lp.active != FALSE)
 175                                 Non_zeros=lp.non_zeros;
 176                         }
 177                 }
 178         }
 179
 180         public void set_obj_fn(Problem lp, float[] row)
 181         {
 182             for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
 183             row[0] = (float)0.0;
 184
 185             int i;
 186             for(i = 1; i <= lp.columns; i++)
 187                 set_mat(lp, 0, i, row[i]);
 188         }
 189
 190
 191         public void add_constraint(Problem lp, float[] row, short constr_type, float rh)
 192         {
 193             for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
 194             row[0] = (float)0.0;
 195
 196             Matrix[] newmat;
 197             int  i, j;
 198             int  elmnr;
 199             int  stcol;
 200             short[]  addtoo;
 201
 202             addtoo = new short[lp.columns + 1];
 203
 204             for(i = 1; i <= lp.columns; i++)
 205                 if(row[i]!=0)
 206                     {
 207                         addtoo[i]=TRUE;
 208                         lp.non_zeros++;
 209                     }
 210                 else
 211                     addtoo[i]=FALSE;
 212
 213             //newmat = new Matrix[lp.non_zeros];
 214
 215             // FIXME
 216             newmat = lp.alternate_mat;
 217             for (i = 0; i < newmat.length; i++) { newmat[i].row_nr = 0; newmat[i].value = 0; }
 218
 219             if (lp.non_zeros > lp.mat_alloc) throw new Error("not enough mat space; this should not happen");
 220             lp.rows++;
 221             lp.sum++;
 222             if(lp.rows > lp.rows_alloc)
 223                 throw new Error("not enough rows; ("+lp.rows+" needed, have "+lp.rows_alloc+") this should never happen");
 224             /*
 225             if (lp.scaling_used != FALSE)
 226                 {
 227                     // shift scale
 228                     for(i=lp.sum; i > lp.rows; i--)
 229                         lp.scale[i]=lp.scale[i-1];
 230                     lp.scale[lp.rows]=1;
 231                 }
 232
 233             if (lp.names_used != FALSE)
 234                 lp.row_name[lp.rows] = new String("r_" + lp.rows);
 235
 236             if(lp.scaling_used != FALSE && lp.columns_scaled != FALSE)
 237                 for(i = 1; i <= lp.columns; i++)
 238                     row[i] *= lp.scale[lp.rows+i];
 239      */
 240             if(constr_type==GE)
 241                 lp.ch_sign[lp.rows] = TRUE;
 242             else
 243                 lp.ch_sign[lp.rows] = FALSE;
 244
 245             elmnr = 0;
 246             stcol = 0;
 247             for(i = 1; i <= lp.columns; i++)
 248                 {
 249                     for(j = stcol; j < lp.col_end[i]; j++)
 250                         {
 251                             newmat[elmnr].row_nr=lp.mat[j].row_nr;
 252                             newmat[elmnr].value=lp.mat[j].value;
 253                             elmnr++;
 254                         }
 255                     if(addtoo[i] != FALSE)
 256                         {
 257                             if(lp.ch_sign[lp.rows] != FALSE)
 258                                 newmat[elmnr].value = -row[i];
 259                             else
 260                                 newmat[elmnr].value = row[i];
 261                             newmat[elmnr].row_nr = lp.rows;
 262                             elmnr++;
 263                         }
 264                     stcol=lp.col_end[i];
 265                     lp.col_end[i]=elmnr;
 266                 }
 267
 268             lp.alternate_mat = lp.mat;
 269             lp.mat = newmat;
 270
 271             for(i=lp.sum ; i > lp.rows; i--)
 272                 {
 273                     lp.orig_upbo[i]=lp.orig_upbo[i-1];
 274                     lp.orig_lowbo[i]=lp.orig_lowbo[i-1];
 275                     lp.basis[i]=lp.basis[i-1];
 276                     lp.lower[i]=lp.lower[i-1];
 277                 }
 278
 279             for(i= 1 ; i <= lp.rows; i++)
 280                 if(lp.bas[i] >= lp.rows)
 281                     lp.bas[i]++;
 282
 283             if(constr_type==LE || constr_type==GE)
 284                 {
 285                     lp.orig_upbo[lp.rows]=lp.infinite;
 286                 }
 287             else if(constr_type==EQ)
 288                 {
 289                     lp.orig_upbo[lp.rows]=0;
 290                 }
 291             else
 292                 {
 293                     System.err.print("Wrong constraint type\n");
 294                     System.exit(FAIL);
 295                 }
 296
 297             lp.orig_lowbo[lp.rows]=0;
 298
 299             if(constr_type==GE && rh != 0)
 300                 lp.orig_rh[lp.rows]=-rh;
 301             else
 302                 lp.orig_rh[lp.rows]=rh;
 303
 304             lp.row_end_valid=FALSE;
 305
 306             lp.bas[lp.rows]=lp.rows;
 307             lp.basis[lp.rows]=TRUE;
 308             lp.lower[lp.rows]=TRUE;
 309
 310             if (lp.active != FALSE)
 311                 set_globals(lp);
 312             lp.eta_valid=FALSE;
 313         }
 314
 315         public void del_constraint(Problem lp, int del_row)
 316         {
 317             int i, j;
 318             int elmnr;
 319             int startcol;
 320
 321             if(del_row < 1 || del_row > lp.rows)
 322                 {
 323                     System.err.println("There is no constraint nr. " + del_row);
 324                     System.exit(FAIL);
 325                 }
 326
 327             elmnr=0;
 328             startcol=0;
 329
 330             for(i = 1; i <= lp.columns; i++)
 331                 {
 332                     for(j=startcol; j < lp.col_end[i]; j++)
 333                         {
 334                             if(lp.mat[j].row_nr!=del_row)
 335                                 {
 336                                     lp.mat[elmnr]=lp.mat[j];
 337                                     if(lp.mat[elmnr].row_nr > del_row)
 338                                         lp.mat[elmnr].row_nr--;
 339                                     elmnr++;
 340                                 }
 341                             else
 342                                 lp.non_zeros--;
 343                         }
 344                     startcol=lp.col_end[i];
 345                     lp.col_end[i]=elmnr;
 346                 }
 347             for(i = del_row; i < lp.rows; i++)
 348                 {
 349                     lp.orig_rh[i] = lp.orig_rh[i+1];
 350                     lp.ch_sign[i] = lp.ch_sign[i+1];
 351                     lp.bas[i] = lp.bas[i+1];
 352                     if (lp.names_used != FALSE)
 353                         lp.row_name[i] = lp.row_name[i+1];
 354                 }
 355             for(i = 1; i < lp.rows; i++)
 356                 if(lp.bas[i] >  del_row)
 357                     lp.bas[i]--;
 358
 359             for(i=del_row; i < lp.sum; i++)
 360                 {
 361                     lp.lower[i]=lp.lower[i+1];
 362                     lp.basis[i]=lp.basis[i+1];
 363                     lp.orig_upbo[i]=lp.orig_upbo[i+1];
 364                     lp.orig_lowbo[i]=lp.orig_lowbo[i+1];
 365                     if (lp.scaling_used != FALSE)
 366                         lp.scale[i]=lp.scale[i+1];
 367                 }
 368
 369             lp.rows--;
 370             lp.sum--;
 371
 372             lp.row_end_valid=FALSE;
 373
 374             if (lp.active != FALSE)
 375                 set_globals(lp);
 376             lp.eta_valid=FALSE;
 377             lp.basis_valid=FALSE;
 378         }
 379
 380
 381         public void add_column(Problem lp, float[] column)
 382         {
 383             int i, elmnr;
 384
 385             lp.columns++;
 386             lp.sum++;
 387             if (lp.columns > lp.columns_alloc) throw new Error("not enough cols; this should never happen");
 388             if (lp.non_zeros + lp.rows+1 > lp.mat_alloc) throw new Error("not enough mat space; this should not happen");
 389
 390             if (lp.scaling_used != FALSE)
 391                 {
 392                     for(i = 0; i <= lp.rows; i++)
 393                         column[i]*=lp.scale[i];
 394                     lp.scale[lp.sum]=1;
 395                 }
 396
 397             elmnr=lp.col_end[lp.columns-1];
 398             for(i = 0 ; i <= lp.rows ; i++)
 399                 if(column[i] != 0)
 400                     {
 401                         lp.mat[elmnr].row_nr=i;
 402                         if(lp.ch_sign[i] != FALSE)
 403                             lp.mat[elmnr].value=-column[i];
 404                         else
 405                             lp.mat[elmnr].value=column[i];
 406                         lp.non_zeros++;
 407                         elmnr++;
 408                     }
 409             lp.col_end[lp.columns]=elmnr;
 410             lp.orig_lowbo[lp.sum]=0;
 411             lp.orig_upbo[lp.sum]=lp.infinite;
 412             lp.lower[lp.sum]=TRUE;
 413             lp.basis[lp.sum]=FALSE;
 414             if (lp.names_used != FALSE)
 415                 lp.col_name[lp.columns] = 0;
 416
 417
 418             lp.row_end_valid=FALSE;
 419
 420             if (lp.active != FALSE)
 421                 {
 422                     Sum=lp.sum;
 423                     columns=lp.columns;
 424                     Non_zeros=lp.non_zeros;
 425                 }
 426         }
 427
 428         public void del_column(Problem lp, int column)
 429         {
 430             int i, j, from_elm, to_elm, elm_in_col;
 431             if(column > lp.columns || column < 1)
 432                 System.err.print("column out of range in del_column");
 433             for(i = 1; i <= lp.rows; i++)
 434                 {
 435                     if(lp.bas[i]==lp.rows+column)
 436                         lp.basis_valid=FALSE;
 437                     else if(lp.bas[i] > lp.rows+column)
 438                         lp.bas[i]--;
 439                 }
 440             for(i = lp.rows+column; i < lp.sum; i++)
 441                 {
 442                     if (lp.names_used != FALSE)
 443                         lp.col_name[i-lp.rows] =  lp.col_name[i-lp.rows+1];
 444                     lp.orig_upbo[i]=lp.orig_upbo[i+1];
 445                     lp.orig_lowbo[i]=lp.orig_lowbo[i+1];
 446                     lp.upbo[i]=lp.upbo[i+1];
 447                     lp.lowbo[i]=lp.lowbo[i+1];
 448                     lp.basis[i]=lp.basis[i+1];
 449                     lp.lower[i]=lp.lower[i+1];
 450                     if (lp.scaling_used != FALSE)
 451                         lp.scale[i]=lp.scale[i+1];
 452                 }
 453             for(i = 0; i < lp.nr_lagrange; i++)
 454                 for(j = column; j <= lp.columns; j++)
 455                     lp.lag_row[i][j]=lp.lag_row[i][j+1];
 456             to_elm=lp.col_end[column-1];
 457             from_elm=lp.col_end[column];
 458             elm_in_col=from_elm-to_elm;
 459             for(i = from_elm; i < lp.non_zeros; i++)
 460                 {
 461                     lp.mat[to_elm]=lp.mat[i];
 462                     to_elm++;
 463                 }
 464             for(i = column; i < lp.columns; i++)
 465                 lp.col_end[i]=lp.col_end[i+1]-elm_in_col;
 466             lp.non_zeros -= elm_in_col;
 467             lp.row_end_valid=FALSE;
 468             lp.eta_valid=FALSE;
 469
 470             lp.sum--;
 471             lp.columns--;
 472             if (lp.active != FALSE)
 473                 set_globals(lp);
 474         }
 475
 476         public void bound_sum(Problem lp, int column1, int column2, float bound, short type, float[] scratch) {
 477             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 478             scratch[column1] = (float)1.0;
 479             scratch[column2] = (float)1.0;
 480             add_constraint(lp, scratch, type, bound);
 481             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 482         }
 483
 484         public void bound_difference(Problem lp, int column1, int column2, float bound, short type, float[] scratch) {
 485             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 486             scratch[column1] = (float)1.0;
 487             scratch[column2] = (float)-1.0;
 488             add_constraint(lp, scratch, type, bound);
 489             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 490         }
 491
 492         public void set_upbo(Problem lp, int column, float value)
 493         {
 494             if(column > lp.columns || column < 1)
 495                 System.err.print("column out of range");
 496             if (lp.scaling_used != FALSE)
 497                 value /= lp.scale[lp.rows + column];
 498             if(value < lp.orig_lowbo[lp.rows + column])
 499                 System.err.print("UpperBound must be >= lowerBound");
 500             lp.eta_valid = FALSE;
 501             lp.orig_upbo[lp.rows+column] = value;
 502         }
 503
 504         public void set_lowbo(Problem lp, int column, float value)
 505         {
 506             if(column > lp.columns || column < 1)
 507                 System.err.print("column out of range");
 508             if (lp.scaling_used != FALSE)
 509                 value /= lp.scale[lp.rows + column];
 510             if(value > lp.orig_upbo[lp.rows + column])
 511                 System.err.print("UpperBound must be >= lowerBound");
 512             lp.eta_valid = FALSE;
 513             lp.orig_lowbo[lp.rows+column] = value;
 514         }
 515
 516         public void set_rh(Problem lp, int row, float value)
 517         {
 518             if(row > lp.rows || row < 0)
 519                 System.err.print("Row out of Range");
 520             if(row == 0)                        /* setting of RHS of OF not meaningful */
 521                 {
 522                     System.err.println("Warning: attempt to set RHS of objective function, ignored");
 523                     return;
 524                 }
 525             if (lp.scaling_used != FALSE)
 526                 if (lp.ch_sign[row] != FALSE)
 527                     lp.orig_rh[row] = -value * lp.scale[row];
 528                 else
 529                     lp.orig_rh[row] = value * lp.scale[row];
 530             else
 531                 if (lp.ch_sign[row] != FALSE)
 532                     lp.orig_rh[row] = -value;
 533                 else
 534                     lp.orig_rh[row] = value;
 535             lp.eta_valid = FALSE;
 536         }
 537
 538         public void set_rh_vec(Problem lp, float[] rh)
 539         {
 540             int i;
 541             if (lp.scaling_used != FALSE)
 542                 for(i = 1; i <= lp.rows; i++)
 543                     if(lp.ch_sign[i] != FALSE)
 544                         lp.orig_rh[i]=-rh[i]*lp.scale[i];
 545                     else
 546                         lp.orig_rh[i]=rh[i]*lp.scale[i];
 547             else
 548                 for(i=1; i <= lp.rows; i++)
 549                     if(lp.ch_sign[i] != FALSE)
 550                         lp.orig_rh[i]=-rh[i];
 551                     else
 552                         lp.orig_rh[i]=rh[i];
 553             lp.eta_valid=FALSE;
 554         }
 555
 556         public void set_maxim(Problem lp)
 557         {
 558             int i;
 559             if(lp.maximise==FALSE)
 560                 {
 561                     for(i = 0; i < lp.non_zeros; i++)
 562                         if(lp.mat[i].row_nr==0)
 563                             lp.mat[i].value*=-1;
 564                     lp.eta_valid=FALSE;
 565                 }
 566             lp.maximise=TRUE;
 567             lp.ch_sign[0]=TRUE;
 568             if (lp.active != FALSE)
 569                 Maximise=TRUE;
 570         }
 571
 572         public void set_minim(Problem lp)
 573         {
 574             int i;
 575             if(lp.maximise==TRUE)
 576                 {
 577                     for(i = 0; i < lp.non_zeros; i++)
 578                         if(lp.mat[i].row_nr==0)
 579                             lp.mat[i].value = -lp.mat[i].value;
 580                     lp.eta_valid=FALSE;
 581                 }
 582             lp.maximise=FALSE;
 583             lp.ch_sign[0]=FALSE;
 584             if (lp.active != FALSE)
 585                 Maximise=FALSE;
 586         }
 587
 588         public void set_constr_type(Problem lp, int row, short con_type)
 589         {
 590             int i;
 591             if(row > lp.rows || row < 1)
 592                 System.err.print("Row out of Range");
 593             if(con_type==EQ)
 594                 {
 595                     lp.orig_upbo[row]=0;
 596                     lp.basis_valid=FALSE;
 597                     if (lp.ch_sign[row] != FALSE)
 598                         {
 599                             for(i = 0; i < lp.non_zeros; i++)
 600                                 if(lp.mat[i].row_nr==row)
 601                                     lp.mat[i].value*=-1;
 602                             lp.eta_valid=FALSE;
 603                             lp.ch_sign[row]=FALSE;
 604                             if(lp.orig_rh[row]!=0)
 605                                 lp.orig_rh[row]*=-1;
 606                         }
 607                 }
 608             else if(con_type==LE)
 609                 {
 610                     lp.orig_upbo[row]=lp.infinite;
 611                     lp.basis_valid=FALSE;
 612                     if (lp.ch_sign[row] != FALSE)
 613                         {
 614                             for(i = 0; i < lp.non_zeros; i++)
 615                                 if(lp.mat[i].row_nr==row)
 616                                     lp.mat[i].value*=-1;
 617                             lp.eta_valid=FALSE;
 618                             lp.ch_sign[row]=FALSE;
 619                             if(lp.orig_rh[row]!=0)
 620                                 lp.orig_rh[row]*=-1;
 621                         }
 622                 }
 623             else if(con_type==GE)
 624                 {
 625                     lp.orig_upbo[row]=lp.infinite;
 626                     lp.basis_valid=FALSE;
 627                     if(lp.ch_sign[row] == FALSE)
 628                         {
 629                             for(i = 0; i < lp.non_zeros; i++)
 630                                 if(lp.mat[i].row_nr==row)
 631                                     lp.mat[i].value*=-1;
 632                             lp.eta_valid=FALSE;
 633                             lp.ch_sign[row]=TRUE;
 634                             if(lp.orig_rh[row]!=0)
 635                                 lp.orig_rh[row]*=-1;
 636                         }
 637                 }
 638             else
 639                 System.err.print("Constraint type not (yet) implemented");
 640         }
 641
 642         public float mat_elm(Problem lp, int row, int column)
 643         {
 644             float value;
 645             int elmnr;
 646             if(row < 0 || row > lp.rows)
 647                 System.err.print("Row out of range in mat_elm");
 648             if(column < 1 || column > lp.columns)
 649                 System.err.print("column out of range in mat_elm");
 650             value=0;
 651             elmnr=lp.col_end[column-1];
 652             while(lp.mat[elmnr].row_nr != row && elmnr < lp.col_end[column])
 653                 elmnr++;
 654             if(elmnr != lp.col_end[column])
 655                 {
 656                     value = lp.mat[elmnr].value;
 657                     if (lp.ch_sign[row] != FALSE)
 658                         value = -value;
 659                     if (lp.scaling_used != FALSE)
 660                         value /= lp.scale[row] * lp.scale[lp.rows + column];
 661                 }
 662             return(value);
 663         }
 664
 665
 666         public void get_row(Problem lp, int row_nr, float[] row)
 667         {
 668             int i, j;
 669
 670             if(row_nr < 0 || row_nr > lp.rows)
 671                 System.err.print("Row nr. out of range in get_row");
 672             for(i = 1; i <= lp.columns; i++)
 673                 {
 674                     row[i]=0;
 675                     for(j=lp.col_end[i-1]; j < lp.col_end[i]; j++)
 676                         if(lp.mat[j].row_nr==row_nr)
 677                             row[i]=lp.mat[j].value;
 678                     if (lp.scaling_used != FALSE)
 679                         row[i] /= lp.scale[lp.rows+i] * lp.scale[row_nr];
 680                 }
 681             if(lp.ch_sign[row_nr] != FALSE)
 682                 for(i=0; i <= lp.columns; i++)
 683                     if(row[i]!=0)
 684                         row[i] = -row[i];
 685         }
 686
 687         public void get_column(Problem lp, int col_nr, float[] column)
 688         {
 689             int i;
 690
 691             if(col_nr < 1 || col_nr > lp.columns)
 692                 System.err.print("Col. nr. out of range in get_column");
 693             for(i = 0; i <= lp.rows; i++)
 694                 column[i]=0;
 695             for(i = lp.col_end[col_nr-1]; i < lp.col_end[col_nr]; i++)
 696                 column[lp.mat[i].row_nr]=lp.mat[i].value;
 697             for(i = 0; i <= lp.rows; i++)
 698                 if(column[i] !=0)
 699                     {
 700                         if(lp.ch_sign[i] != FALSE)
 701                             column[i]*=-1;
 702                         if (lp.scaling_used != FALSE)
 703                             column[i]/=(lp.scale[i] * lp.scale[lp.rows+col_nr]);
 704                     }
 705         }
 706
 707         public void get_reduced_costs(Problem lp, float[] rc)
 708         {
 709             int varnr, i, j;
 710             float f;
 711
 712             if(lp.basis_valid == FALSE)
 713                 System.err.print("Not a valid basis in get_reduced_costs");
 714             set_globals(lp);
 715             if(lp.eta_valid == FALSE)
 716                 invert();
 717             for(i = 1; i <= lp.sum; i++)
 718                 rc[i] = 0;
 719             rc[0] = 1;
 720             btran(rc);
 721             for(i = 1; i <= lp.columns; i++)
 722                 {
 723                     varnr = lp.rows + i;
 724                     if(Basis[varnr] == FALSE)
 725                         if(Upbo[varnr] > 0)
 726                             {
 727                                 f = 0;
 728                                 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
 729                                     f += rc[Mat[j].row_nr] * Mat[j].value;
 730                                 rc[varnr] = f;
 731                             }
 732                 }
 733             for(i = 1; i <= Sum; i++)
 734                 rc[i] = round(rc[i], Epsd);
 735         }
 736
 737         short is_feasible(Problem lp, float[] values)
 738         {
 739             int i, elmnr;
 740             float[] this_rhs;
 741             float dist;
 742
 743             if (lp.scaling_used != FALSE)
 744                 {
 745                     for(i = lp.rows+1; i <= lp.sum; i++)
 746                         if( values[i - lp.rows] < lp.orig_lowbo[i]*lp.scale[i]
 747                             || values[i - lp.rows] > lp.orig_upbo[i]*lp.scale[i])
 748                             return(FALSE);
 749                 }
 750             else
 751                 {
 752                     for(i = lp.rows+1; i <= lp.sum; i++)
 753                         if(   values[i - lp.rows] < lp.orig_lowbo[i]
 754                               || values[i - lp.rows] > lp.orig_upbo[i])
 755                             return(FALSE);
 756                 }
 757             this_rhs = new float[lp.rows+1];
 758             for (i = 0; i <= lp.rows; i++)
 759                 this_rhs[i] = 0;
 760
 761             for(i = 1; i <= lp.columns; i++)
 762                 for(elmnr = lp.col_end[i - 1]; elmnr < lp.col_end[i]; elmnr++)
 763                     this_rhs[lp.mat[elmnr].row_nr] += lp.mat[elmnr].value * values[i];
 764             for(i = 1; i <= lp.rows; i++)
 765                 {
 766                     dist = lp.orig_rh[i] - this_rhs[i];
 767                     dist = round(dist, (float)0.001);
 768                     if((lp.orig_upbo[i] == 0 && dist != 0) || dist < 0)
 769                         {
 770                             return(FALSE);
 771                         }
 772                 }
 773             return(TRUE);
 774         }
 775
 776         short column_in_lp(Problem lp, float[] testcolumn)
 777         {
 778             int i, j;
 779             short ident;
 780             float value;
 781
 782             if (lp.scaling_used != FALSE)
 783                 for(i = 1; i <= lp.columns; i++)
 784                     {
 785                         ident = TRUE;
 786                         j = lp.col_end[i-1];
 787                         while(ident != FALSE && (j < lp.col_end[i]))
 788                             {
 789                                 value = lp.mat[j].value;
 790                                 if(lp.ch_sign[lp.mat[j].row_nr] != FALSE)
 791                                     value = -value;
 792                                 value /= lp.scale[lp.rows+i];
 793                                 value /= lp.scale[lp.mat[j].row_nr];
 794                                 value -= testcolumn[lp.mat[j].row_nr];
 795                                 if(Math.abs(value) >  (float)0.001) /* should be some epsilon? MB */
 796                                     ident=FALSE;
 797                                 j++;
 798                             }
 799                         if(ident != FALSE)
 800                             return(TRUE);
 801                     }
 802             else
 803                 for(i = 1; i <= lp.columns; i++)
 804                     {
 805                         ident = TRUE;
 806                         j = lp.col_end[i-1];
 807                         while(ident != FALSE && j < lp.col_end[i])
 808                             {
 809                                 value = lp.mat[j].value;
 810                                 if(lp.ch_sign[lp.mat[j].row_nr] != FALSE)
 811                                     value *= -1;
 812                                 value -= testcolumn[lp.mat[j].row_nr];
 813                                 if( Math.abs(value) >  (float)0.001 )
 814                                     ident=FALSE;
 815                                 j++;
 816                             }
 817                         if(ident != FALSE)
 818                             return(TRUE);
 819                     }
 820             return(FALSE);
 821         }
 822
 823         private float minmax_to_scale(float min, float max)
 824         {
 825             float scale;
 826
 827             /* should do something sensible when min or max is 0, MB */
 828             if((min == 0) || (max == 0))
 829                 return((float)1);
 830
 831             scale = (float)(1 / Math.pow(Math.E, (Math.log(min) + Math.log(max)) / 2));
 832             return(scale);
 833         }
 834
 835         public void unscale_columns(Problem lp)
 836         {
 837             int i, j;
 838
 839             /* unscale mat */
 840             for(j = 1; j <= lp.columns; j++)
 841                 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
 842                     lp.mat[i].value /= lp.scale[lp.rows + j];
 843
 844             /* unscale Bounds as well */
 845             for(i = lp.rows + 1; i < lp.sum; i++)
 846                 {
 847                     if(lp.orig_lowbo[i] != 0)
 848                         lp.orig_lowbo[i] *= lp.scale[i];
 849                     if(lp.orig_upbo[i] != lp.infinite)
 850                         lp.orig_upbo[i] *= lp.scale[i];
 851                 }
 852
 853             for(i=lp.rows+1; i<= lp.sum; i++)
 854                 lp.scale[i]=1;
 855             lp.columns_scaled=FALSE;
 856             lp.eta_valid=FALSE;
 857         }
 858
 859         public void unscale(Problem lp)
 860         {
 861             int i, j;
 862
 863             if (lp.scaling_used != FALSE)
 864                 {
 865
 866                     /* unscale mat */
 867                     for(j = 1; j <= lp.columns; j++)
 868                         for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
 869                             lp.mat[i].value /= lp.scale[lp.rows + j];
 870
 871                     /* unscale Bounds */
 872                     for(i = lp.rows + 1; i < lp.sum; i++)
 873                         {
 874                             if(lp.orig_lowbo[i] != 0)
 875                                 lp.orig_lowbo[i] *= lp.scale[i];
 876                             if(lp.orig_upbo[i] != lp.infinite)
 877                                 lp.orig_upbo[i] *= lp.scale[i];
 878                         }
 879
 880                     /* unscale the matrix */
 881                     for(j = 1; j <= lp.columns; j++)
 882                         for(i = lp.col_end[j-1]; i < lp.col_end[j]; i++)
 883                             lp.mat[i].value /= lp.scale[lp.mat[i].row_nr];
 884
 885                     /* unscale the rhs! */
 886                     for(i = 0; i <= lp.rows; i++)
 887                         lp.orig_rh[i] /= lp.scale[i];
 888
 889                     lp.scaling_used=FALSE;
 890                     lp.eta_valid=FALSE;
 891                 }
 892         }
 893
 894
 895         public void auto_scale(Problem lp)
 896         {
 897             int i, j, row_nr, IntUsed;
 898             float row_max[], row_min[], scalechange[], absval;
 899
 900             if(lp.scaling_used == 0)
 901                 {
 902                     lp.scale = new float[lp.sum_alloc + 1];
 903                     for(i = 0; i <= lp.sum; i++)
 904                         lp.scale[i]=1;
 905                 }
 906
 907             row_max = new float[lp.rows + 1];
 908             row_min = new float[lp.rows + 1];
 909             scalechange = new float[lp.sum + 1];
 910
 911             /* initialise min and max values */
 912             for(i = 0; i <= lp.rows; i++)
 913                 {
 914                     row_max[i]=0;
 915                     row_min[i]=lp.infinite;
 916                 }
 917
 918             /* calculate min and max absolute values of rows */
 919             for(j = 1; j <= lp.columns; j++)
 920                 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
 921                     {
 922                         row_nr = lp.mat[i].row_nr;
 923                         absval = Math.abs(lp.mat[i].value);
 924                         if(absval!=0)
 925                             {
 926                                 row_max[row_nr] = Math.max(row_max[row_nr], absval);
 927                                 row_min[row_nr] = Math.min(row_min[row_nr], absval);
 928                             }
 929                     }
 930             /* calculate scale factors for rows */
 931             for(i = 0; i <= lp.rows; i++)
 932                 {
 933                     scalechange[i] = minmax_to_scale(row_min[i], row_max[i]);
 934                     lp.scale[i] *= scalechange[i];
 935                 }
 936
 937             /* now actually scale the matrix */
 938             for(j = 1; j <= lp.columns; j++)
 939                 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
 940                     lp.mat[i].value *= scalechange[lp.mat[i].row_nr];
 941
 942             /* and scale the rhs! */
 943             for(i = 0; i <= lp.rows; i++)
 944                 lp.orig_rh[i] *= scalechange[i];
 945
 946             IntUsed=FALSE;
 947             i=lp.rows+1;
 948             while(IntUsed == FALSE && i <= lp.sum)
 949                 {
 950                     IntUsed=FALSE;
 951                     i++;
 952                 }
 953             if(IntUsed == FALSE)
 954                 {
 955                     float col_max, col_min;
 956
 957                     /* calculate scales */
 958                     for(j = 1; j <= lp.columns; j++)
 959                         {
 960                             col_max = 0;
 961                             col_min = lp.infinite;
 962                             for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
 963                                 {
 964                                     if(lp.mat[i].value!=0)
 965                                         {
 966                                             col_max = Math.max(col_max, Math.abs(lp.mat[i].value));
 967                                             col_min = Math.min(col_min, Math.abs(lp.mat[i].value));
 968                                         }
 969                                 }
 970                             scalechange[lp.rows + j]  = minmax_to_scale(col_min, col_max);
 971                             lp.scale[lp.rows + j] *= scalechange[lp.rows + j];
 972                         }
 973
 974                     /* scale mat */
 975                     for(j = 1; j <= lp.columns; j++)
 976                         for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
 977                             lp.mat[i].value *= scalechange[lp.rows + j];
 978
 979                     /* scale Bounds as well */
 980                     for(i = lp.rows + 1; i < lp.sum; i++)
 981                         {
 982                             if(lp.orig_lowbo[i] != 0)
 983                                 lp.orig_lowbo[i] /= scalechange[i];
 984                             if(lp.orig_upbo[i] != lp.infinite)
 985                                 lp.orig_upbo[i] /= scalechange[i];
 986                         }
 987                     lp.columns_scaled=TRUE;
 988                 }
 989             lp.scaling_used=TRUE;
 990             lp.eta_valid=FALSE;
 991         }
 992
 993         public void reset_basis(Problem lp) { lp.basis_valid=FALSE; }
 994
 995         /* Globals used by solver */
 996         short JustInverted;
 997         short Status;
 998         short Doiter;
 999         short DoInvert;
1000         short Break_bb;
1001
1002         void set_globals(Problem lp)
1003         {
1004             if(Lp != null)
1005                 Lp.active = FALSE;
1006             Lp = lp;
1007             Rows = lp.rows;
1008             columns = lp.columns;
1009             Sum = Rows + columns;
1010             Non_zeros = lp.non_zeros;
1011             Mat = lp.mat;
1012             Col_no = lp.col_no;
1013             Col_end = lp.col_end;
1014             Row_end = lp.row_end;
1015             Rh = lp.rh;
1016             Rhs = lp.rhs;
1017             Orig_rh = lp.orig_rh;
1018             Orig_upbo = lp.orig_upbo;
1019             Orig_lowbo = lp.orig_lowbo;
1020             Upbo = lp.upbo;
1021             Lowbo = lp.lowbo;
1022             Bas = lp.bas;
1023             Basis = lp.basis;
1024             Lower = lp.lower;
1025             Eta_alloc = lp.eta_alloc;
1026             Eta_size = lp.eta_size;
1027             Num_inv = lp.num_inv;
1028             Eta_value = lp.eta_value;
1029             Eta_row_nr = lp.eta_row_nr;
1030             Eta_col_end = lp.eta_col_end;
1031             Solution = lp.solution;
1032             Best_solution = lp.best_solution;
1033             Infinite = lp.infinite;
1034             Epsilon = lp.epsilon;
1035             Epsb = lp.epsb;
1036             Epsd = lp.epsd;
1037             Epsel = lp.epsel;
1038
1039             /* ?? MB */
1040             TREJ = TREJ;
1041             TINV = TINV;
1042
1043             Maximise = lp.maximise;
1044             Floor_first = lp.floor_first;
1045             lp.active = TRUE;
1046
1047             //  System.out.println("Sum = " + Sum);
1048         } // end of set_global
1049
1050         private void ftran(int start,
1051                            int end,
1052                            float[] pcol)
1053         {
1054             int  i, j, k, r;
1055             float theta;
1056
1057             for(i = start; i <= end; i++)
1058                 {
1059                     k = Eta_col_end[i] - 1;
1060                     r = Eta_row_nr[k];
1061                     theta = pcol[r];
1062                     if(theta != 0)
1063                         for(j = Eta_col_end[i - 1]; j < k; j++)
1064                             pcol[Eta_row_nr[j]] += theta * Eta_value[j]; /* cpu expensive line */
1065                     pcol[r] *= Eta_value[k];
1066                 }
1067             /* round small values to zero */
1068             for(i = 0; i <= Rows; i++)
1069                 round(pcol[i], Epsel);
1070         } /* ftran */
1071
1072         private void btran(float[] row)
1073         {
1074             int  i, j, k;
1075             float f;
1076
1077             for(i = Eta_size; i >= 1; i--)
1078                 {
1079                     f = 0;
1080                     k = Eta_col_end[i] - 1;
1081                     for(j = Eta_col_end[i - 1]; j <= k; j++)
1082                         f += row[Eta_row_nr[j]] * Eta_value[j];
1083                     f = round(f, Epsel);
1084                     row[Eta_row_nr[k]] = f;
1085                 }
1086         } /* btran */
1087
1088
1089         static short Isvalid(Problem lp)
1090         {
1091             int i, j, rownum[], colnum[];
1092             int num[], row_nr;
1093
1094             if(lp.row_end_valid == FALSE)
1095                 {
1096                     num = new int[lp.rows + 1];
1097                     rownum = new int[lp.rows + 1];
1098                     for(i = 0; i <= lp.rows; i++)
1099                         {
1100                             num[i] = 0;
1101                             rownum[i] = 0;
1102                         }
1103                     for(i = 0; i < lp.non_zeros; i++)
1104                         rownum[lp.mat[i].row_nr]++;
1105                     lp.row_end[0] = 0;
1106                     for(i = 1; i <= lp.rows; i++)
1107                         lp.row_end[i] = lp.row_end[i - 1] + rownum[i];
1108                     for(i = 1; i <= lp.columns; i++)
1109                         for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++)
1110                             {
1111                                 row_nr = lp.mat[j].row_nr;
1112                                 if(row_nr != 0)
1113                                     {
1114                                         num[row_nr]++;
1115                                         lp.col_no[lp.row_end[row_nr - 1] + num[row_nr]] = i;
1116                                     }
1117                             }
1118                     lp.row_end_valid = TRUE;
1119                 }
1120             if(lp.valid != FALSE)
1121                 return(TRUE);
1122             rownum = new int[lp.rows + 1];
1123             for (i = 0; i <= lp.rows; i++)
1124                 rownum[i] = 0;
1125             colnum = new int[lp.columns + 1];
1126             for (i = 0; i <= lp.columns; i++)
1127                 colnum[i] = 0;
1128             for(i = 1 ; i <= lp.columns; i++)
1129                 for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++)
1130                     {
1131                         colnum[i]++;
1132                         rownum[lp.mat[j].row_nr]++;
1133                     }
1134             for(i = 1; i <= lp.columns; i++)
1135                 if(colnum[i] == 0)
1136                     if (lp.names_used != FALSE)
1137                         System.err.print("Warning: Variable " + lp.col_name[i] +
1138                                          " not used in any constraints\n");
1139                     else
1140                         System.err.print("Warning: Variable " + i + " not used in any constaints\n");
1141             lp.valid = TRUE;
1142             return(TRUE);
1143         }
1144
1145         private void resize_eta()
1146         {
1147             Eta_alloc *= 2;
1148             float[] db_ptr = Eta_value;
1149             Eta_value = new float[Eta_alloc];
1150             System.arraycopy(db_ptr, 0, Eta_value, 0, db_ptr.length);
1151             Lp.eta_value = Eta_value;
1152
1153             int[] int_ptr = Eta_row_nr;
1154             Eta_row_nr = new int[Eta_alloc];
1155             System.arraycopy(int_ptr, 0, Eta_row_nr, 0, int_ptr.length);
1156             Lp.eta_row_nr = Eta_row_nr;
1157         } /* resize_eta */
1158
1159
1160         private void condensecol(int row_nr,
1161                                  float[] pcol)
1162         {
1163             int i, elnr;
1164
1165             elnr = Eta_col_end[Eta_size];
1166
1167             if(elnr + Rows + 2 > Eta_alloc) /* maximum local growth of Eta */
1168                 resize_eta();
1169
1170             for(i = 0; i <= Rows; i++)
1171                 if(i != row_nr && pcol[i] != 0)
1172                     {
1173                         Eta_row_nr[elnr] = i;
1174                         Eta_value[elnr] = pcol[i];
1175                         elnr++;
1176                     }
1177             Eta_row_nr[elnr] = row_nr;
1178             Eta_value[elnr] = pcol[row_nr];
1179             elnr++;
1180             Eta_col_end[Eta_size + 1] = elnr;
1181         } /* condensecol */
1182
1183
1184         private void addetacol()
1185         {
1186             int  i, j, k;
1187             float theta;
1188
1189             j = Eta_col_end[Eta_size];
1190             Eta_size++;
1191             k = Eta_col_end[Eta_size];
1192             theta = 1 / (float) Eta_value[k - 1];
1193             Eta_value[k - 1] = theta;
1194             for(i = j; i < k - 1; i++)
1195                 Eta_value[i] *= -theta;
1196             JustInverted = FALSE;
1197         } /* addetacol */
1198
1199         private void setpivcol(short lower,
1200                                int varin,
1201                                float[]   pcol)
1202         {
1203             int  i, colnr;
1204             float f;
1205
1206             if(lower != FALSE)
1207                 f = 1;
1208             else
1209                 f = -1;
1210             for(i = 0; i <= Rows; i++)
1211                 pcol[i] = 0;
1212             if(varin > Rows)
1213                 {
1214                     colnr = varin - Rows;
1215                     for(i = Col_end[colnr - 1]; i < Col_end[colnr]; i++)
1216                         pcol[Mat[i].row_nr] = Mat[i].value * f;
1217                     pcol[0] -= Extrad * f;
1218                 }
1219             else
1220                 if(lower != FALSE)
1221                     pcol[varin] = 1;
1222                 else
1223                     pcol[varin] = -1;
1224             ftran(1, Eta_size, pcol);
1225         } /* setpivcol */
1226
1227
1228         private void minoriteration(int colnr,
1229                                     int row_nr)
1230         {
1231             int  i, j, k, wk, varin, varout, elnr;
1232             float piv = 0, theta;
1233
1234             varin = colnr + Rows;
1235             elnr = Eta_col_end[Eta_size];
1236             wk = elnr;
1237             Eta_size++;
1238             if(Extrad != 0)
1239                 {
1240                     Eta_row_nr[elnr] = 0;
1241                     Eta_value[elnr] = -Extrad;
1242                     elnr++;
1243                 }
1244             for(j = Col_end[colnr - 1] ; j < Col_end[colnr]; j++)
1245                 {
1246                     k = Mat[j].row_nr;
1247                     if(k == 0 && Extrad != 0)
1248                         Eta_value[Eta_col_end[Eta_size -1]] += Mat[j].value;
1249                     else if(k != row_nr)
1250                         {
1251                             Eta_row_nr[elnr] = k;
1252                             Eta_value[elnr] = Mat[j].value;
1253                             elnr++;
1254                         }
1255                     else
1256                         piv = Mat[j].value;
1257                 }
1258             Eta_row_nr[elnr] = row_nr;
1259             Eta_value[elnr] = 1 / (float) piv;
1260             elnr++;
1261             theta = Rhs[row_nr] / (float) piv;
1262             Rhs[row_nr] = theta;
1263             for(i = wk; i < elnr - 1; i++)
1264                 Rhs[Eta_row_nr[i]] -= theta * Eta_value[i];
1265             varout = Bas[row_nr];
1266             Bas[row_nr] = varin;
1267             Basis[varout] = FALSE;
1268             Basis[varin] = TRUE;
1269             for(i = wk; i < elnr - 1; i++)
1270                 Eta_value[i] /= - (float) piv;
1271             Eta_col_end[Eta_size] = elnr;
1272         } /* minoriteration */
1273
1274         private void rhsmincol(float theta,
1275                                int row_nr,
1276                                int varin)
1277         {
1278             int  i, j, k, varout;
1279             float f;
1280
1281             if(row_nr > Rows + 1)
1282                 {
1283                     System.err.print("Error: rhsmincol called with row_nr: " +
1284                                      row_nr + ", rows: " + Rows + "\n");
1285                     System.err.print("This indicates numerical instability\n");
1286                     System.exit(FAIL);
1287                 }
1288             j = Eta_col_end[Eta_size];
1289             k = Eta_col_end[Eta_size + 1];
1290             for(i = j; i < k; i++)
1291                 {
1292                     f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i];
1293                     f = round(f, Epsb);
1294                     Rhs[Eta_row_nr[i]] = f;
1295                 }
1296             Rhs[row_nr] = theta;
1297             varout = Bas[row_nr];
1298             Bas[row_nr] = varin;
1299             Basis[varout] = FALSE;
1300             Basis[varin] = TRUE;
1301         } /* rhsmincol */
1302
1303
1304         void invert()
1305         {
1306             int    i, j, v, wk, numit, varnr, row_nr, colnr, varin;
1307             float   theta;
1308             float[]   pcol;
1309             short[]  frow;
1310             short[]  fcol;
1311             int[]    rownum, col, row;
1312             int[]    colnum;
1313
1314             rownum = new int[Rows + 1];
1315             for (i = 0; i <= Rows; i++)
1316                 rownum[i] = 0;
1317             col = new int[Rows + 1];
1318             for (i = 0; i <= Rows; i++)
1319                 col[i] = 0;
1320             row = new int[Rows + 1];
1321             for (i = 0; i <= Rows; i++)
1322                 row[i] = 0;
1323             pcol = new float[Rows + 1];
1324             for (i = 0; i <= Rows; i++)
1325                 pcol[i] = 0;
1326             frow = new short[Rows + 1];
1327             for(i = 0; i <= Rows; i++)
1328                 frow[i] = TRUE;
1329             fcol = new short[columns + 1];
1330             for(i = 0; i < columns; i++)
1331                 fcol[i] = FALSE;
1332             colnum = new int[columns + 1];
1333             for(i = 0; i <= columns; i++)
1334                 colnum[i] = 0;
1335
1336             for(i = 0; i <= Rows; i++)
1337                 if(Bas[i] > Rows)
1338                     fcol[Bas[i] - Rows - 1] = TRUE;
1339                 else
1340                     frow[Bas[i]] = FALSE;
1341
1342             for(i = 1; i <= Rows; i++)
1343                 if(frow[i] != FALSE)
1344                     for(j = Row_end[i - 1] + 1; j <= Row_end[i]; j++)
1345                         {
1346                             wk = Col_no[j];
1347                             if(fcol[wk - 1] != FALSE)
1348                                 {
1349                                     colnum[wk]++;
1350                                     rownum[i - 1]++;
1351                                 }
1352                         }
1353             for(i = 1; i <= Rows; i++)
1354                 Bas[i] = i;
1355             for(i = 1; i <= Rows; i++)
1356                 Basis[i] = TRUE;
1357             for(i = 1; i <= columns; i++)
1358                 Basis[i + Rows] = FALSE;
1359
1360             for(i = 0; i <= Rows; i++)
1361                 Rhs[i] = Rh[i];
1362
1363             for(i = 1; i <= columns; i++)
1364                 {
1365                     varnr = Rows + i;
1366                     if(Lower[varnr] == FALSE)
1367                         {
1368                             theta = Upbo[varnr];
1369                             for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1370                                 Rhs[Mat[j].row_nr] -= theta * Mat[j].value;
1371                         }
1372                 }
1373             for(i = 1; i <= Rows; i++)
1374                 if(Lower[i] == FALSE)
1375                     Rhs[i] -= Upbo[i];
1376             Eta_size = 0;
1377             v = 0;
1378             row_nr = 0;
1379             Num_inv = 0;
1380             numit = 0;
1381             while(v < Rows)
1382                 {
1383                     row_nr++;
1384                     if(row_nr > Rows)
1385                         row_nr = 1;
1386                     v++;
1387                     if(rownum[row_nr - 1] == 1)
1388                         if(frow[row_nr] != FALSE)
1389                             {
1390                                 v = 0;
1391                                 j = Row_end[row_nr - 1] + 1;
1392                                 while(fcol[Col_no[j] - 1] == FALSE)
1393                                     j++;
1394                                 colnr = Col_no[j];
1395                                 fcol[colnr - 1] = FALSE;
1396                                 colnum[colnr] = 0;
1397                                 for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++)
1398                                     if(frow[Mat[j].row_nr] != FALSE)
1399                                         rownum[Mat[j].row_nr - 1]--;
1400                                 frow[row_nr] = FALSE;
1401                                 minoriteration(colnr, row_nr);
1402                             }
1403                 }
1404             v = 0;
1405             colnr = 0;
1406             while(v <columns)
1407                 {
1408                     colnr++;
1409                     if(colnr > columns)
1410                         colnr = 1;
1411                     v++;
1412                     if(colnum[colnr] == 1)
1413                         if(fcol[colnr - 1] != FALSE)
1414                             {
1415                                 v = 0;
1416                                 j = Col_end[colnr - 1] + 1;
1417                                 while(frow[Mat[j - 1].row_nr] == FALSE)
1418                                     j++;
1419                                 row_nr = Mat[j - 1].row_nr;
1420                                 frow[row_nr] = FALSE;
1421                                 rownum[row_nr - 1] = 0;
1422                                 for(j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++)
1423                                     if(fcol[Col_no[j] - 1] != FALSE)
1424                                         colnum[Col_no[j]]--;
1425                                 fcol[colnr - 1] = FALSE;
1426                                 numit++;
1427                                 col[numit - 1] = colnr;
1428                                 row[numit - 1] = row_nr;
1429                             }
1430                 }
1431             for(j = 1; j <= columns; j++)
1432                 if(fcol[j - 1] != FALSE)
1433                     {
1434                         fcol[j - 1] = FALSE;
1435                         setpivcol(Lower[Rows + j], j + Rows, pcol);
1436                         row_nr = 1;
1437                         while((frow[row_nr] == FALSE || pcol[row_nr] == FALSE)
1438                               && row_nr <= Rows)
1439                             row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes
1440                                          rhsmincol crash. Solved in 2.0? MB */
1441                         if(row_nr == Rows + 1)
1442                             System.err.print("Inverting failed");
1443                         frow[row_nr] = FALSE;
1444                         condensecol(row_nr, pcol);
1445                         theta = Rhs[row_nr] / (float) pcol[row_nr];
1446                         rhsmincol(theta, row_nr, Rows + j);
1447                         addetacol();
1448                     }
1449             for(i = numit - 1; i >= 0; i--)
1450                 {
1451                     colnr = col[i];
1452                     row_nr = row[i];
1453                     varin = colnr + Rows;
1454                     for(j = 0; j <= Rows; j++)
1455                         pcol[j] = 0;
1456                     for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++)
1457                         pcol[Mat[j].row_nr] = Mat[j].value;
1458                     pcol[0] -= Extrad;
1459                     condensecol(row_nr, pcol);
1460                     theta = Rhs[row_nr] / (float) pcol[row_nr];
1461                     rhsmincol(theta, row_nr, varin);
1462                     addetacol();
1463                 }
1464             for(i = 1; i <= Rows; i++)
1465                 Rhs[i] = round(Rhs[i], Epsb);
1466             JustInverted = TRUE;
1467             DoInvert = FALSE;
1468         } /* invert */
1469
1470         private short colprim(Ref colnr,
1471                               short minit,
1472                               float[]   drow)
1473         {
1474             int  varnr, i, j;
1475             float f, dpiv;
1476
1477             dpiv = -Epsd;
1478             colnr.value = 0;
1479             if(minit == FALSE)
1480                 {
1481                     for(i = 1; i <= Sum; i++)
1482                         drow[i] = 0;
1483                     drow[0] = 1;
1484                     btran(drow);
1485                     for(i = 1; i <= columns; i++)
1486                         {
1487                             varnr = Rows + i;
1488                             if(Basis[varnr] == FALSE)
1489                                 if(Upbo[varnr] > 0)
1490                                     {
1491                                         f = 0;
1492                                         for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1493                                             f += drow[Mat[j].row_nr] * Mat[j].value;
1494                                         drow[varnr] = f;
1495                                     }
1496                         }
1497                     for(i = 1; i <= Sum; i++)
1498                         drow[i] = round(drow[i], Epsd);
1499                 }
1500             /*
1501               System.out.println("dpiv = " + dpiv);
1502               System.out.println("drow[] at colprim");
1503               for(i = 1; i <= Sum; i++) // DEBUG
1504               {
1505               System.out.print(drow[i] + " ");
1506               }
1507               System.out.println();
1508             */
1509             for(i = 1; i <= Sum; i++)
1510                 if(Basis[i] == FALSE)
1511                     if(Upbo[i] > 0)
1512                         {
1513                             if(Lower[i] != FALSE)
1514                                 f = drow[i];
1515                             else
1516                                 f = -drow[i];
1517                             if(f < dpiv)
1518                                 {
1519                                     dpiv = f;
1520                                     colnr.value = i;
1521                                 }
1522                         }
1523             if(Lp.trace != FALSE)
1524                 if(colnr.value > 0)
1525                     System.err.print("col_prim:" + colnr.value + ", reduced cost: " + dpiv + "\n");
1526                 else
1527                     System.err.print("col_prim: no negative reduced costs found, optimality!\n");
1528             if(colnr.value == 0)
1529                 {
1530                     Doiter   = FALSE;
1531                     DoInvert = FALSE;
1532                     Status   = OPTIMAL;
1533                 }
1534             return(colnr.value > 0 ? (short)1 : (short)0);
1535         } /* colprim */
1536
1537         private short rowprim(int colnr,
1538                               Ref row_nr,
1539                               Ref theta,
1540                               float[] pcol)
1541         {
1542             int  i;
1543             float f = 0, quot;
1544
1545             row_nr.value = 0;
1546             theta.value = Infinite;
1547             for(i = 1; i <= Rows; i++)
1548                 {
1549                     f = pcol[i];
1550                     if(Math.abs(f) < TREJ)
1551                         f = 0;
1552                     if(f != 0)
1553                         {
1554                             quot = 2 * Infinite;
1555                             if(f > 0)
1556                                 quot = Rhs[i] / (float) f;
1557                             else
1558                                 if(Upbo[Bas[i]] < Infinite)
1559                                     quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
1560                             round(quot, Epsel);
1561                             if(quot < theta.value)
1562                                 {
1563                                     theta.value = quot;
1564                                     row_nr.value = i;
1565                                 }
1566                         }
1567                 }
1568             if(row_nr.value == 0)
1569                 for(i = 1; i <= Rows; i++)
1570                     {
1571                         f = pcol[i];
1572                         if(f != 0)
1573                             {
1574                                 quot = 2 * Infinite;
1575                                 if(f > 0)
1576                                     quot = Rhs[i] / (float) f;
1577                                 else
1578                                     if(Upbo[Bas[i]] < Infinite)
1579                                         quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
1580                                 quot = round(quot, Epsel);
1581                                 if(quot < theta.value)
1582                                     {
1583                                         theta.value = quot;
1584                                         row_nr.value = i;
1585                                     }
1586                             }
1587                     }
1588
1589             if(theta.value < 0)
1590                 {
1591                     System.err.println("Warning: Numerical instability, qout = " + theta.value);
1592                     System.err.println("pcol[" + row_nr.value + "] = " + f + ", rhs[" +
1593                                        row_nr.value + "] = " + Rhs[(int)row_nr.value] +
1594                                        " , upbo = " + Upbo[Bas[(int)row_nr.value]]);
1595                 }
1596             if(row_nr.value == 0)
1597                 {
1598                     if(Upbo[colnr] == Infinite)
1599                         {
1600                             Doiter   = FALSE;
1601                             DoInvert = FALSE;
1602                             Status   = UNBOUNDED;
1603                         }
1604                     else
1605                         {
1606                             i = 1;
1607                             while(pcol[i] >= 0 && i <= Rows)
1608                                 i++;
1609                             if(i > Rows) /* empty column with upperBound! */
1610                                 {
1611                                     Lower[colnr] = FALSE;
1612                                     Rhs[0] += Upbo[colnr]*pcol[0];
1613                                     Doiter = FALSE;
1614                                     DoInvert = FALSE;
1615                                 }
1616                             else if(pcol[i]<0)
1617                                 {
1618                                     row_nr.value = i;
1619                                 }
1620                         }
1621                 }
1622             if(row_nr.value > 0)
1623                 Doiter = TRUE;
1624             if(Lp.trace != FALSE)
1625                 System.out.println("row_prim:" + row_nr.value + ", pivot element:" +
1626                                    pcol[(int)row_nr.value]);
1627
1628             return((row_nr.value > 0) ? (short)1 : (short)0);
1629         } /* rowprim */
1630
1631         private short rowdual(Ref row_nr)
1632         {
1633             int   i;
1634             float  f, g, minrhs;
1635             short artifs;
1636
1637             row_nr.value = 0;
1638             minrhs = -Epsb;
1639             i = 0;
1640             artifs = FALSE;
1641             while(i < Rows && artifs == FALSE)
1642                 {
1643                     i++;
1644                     f = Upbo[Bas[i]];
1645                     if(f == 0 && (Rhs[i] != 0))
1646                         {
1647                             artifs = TRUE;
1648                             row_nr.value = i;
1649                         }
1650                     else
1651                         {
1652                             if(Rhs[i] < f - Rhs[i])
1653                                 g = Rhs[i];
1654                             else
1655                                 g = f - Rhs[i];
1656                             if(g < minrhs)
1657                                 {
1658                                     minrhs = g;
1659                                     row_nr.value = i;
1660                                 }
1661                         }
1662                 }
1663
1664             if(Lp.trace != FALSE)
1665                 {
1666                     if(row_nr.value > 0)
1667                         {
1668                             System.out.println("row_dual:" + row_nr.value +
1669                                                ", rhs of selected row:           "
1670                                                + Rhs[(int)row_nr.value]);
1671                             if(Upbo[Bas[(int)row_nr.value]] < Infinite)
1672                                 System.out.println("               upper Bound of basis variable:    " +
1673                                                    Upbo[Bas[(int)row_nr.value]]);
1674                         }
1675                     else
1676                         System.out.print("row_dual: no infeasibilities found\n");
1677                 }
1678
1679             return(row_nr.value > 0 ? (short)1 : (short)0);
1680         } /* rowdual */
1681
1682         private short coldual(int row_nr,
1683                               Ref colnr,
1684                               short minit,
1685                               float[] prow,
1686                               float[] drow)
1687         {
1688             int  i, j, r, varnr;
1689             float theta, quot, pivot, d, f, g;
1690
1691             Doiter = FALSE;
1692             if(minit == FALSE)
1693                 {
1694                     for(i = 0; i <= Rows; i++)
1695                         {
1696                             prow[i] = 0;
1697                             drow[i] = 0;
1698                         }
1699                     drow[0] = 1;
1700                     prow[row_nr] = 1;
1701                     for(i = Eta_size; i >= 1; i--)
1702                         {
1703                             d = 0;
1704                             f = 0;
1705                             r = Eta_row_nr[Eta_col_end[i] - 1];
1706                             for(j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++)
1707                                 {
1708                                     /* this is where the program consumes most cpu time */
1709                                     f += prow[Eta_row_nr[j]] * Eta_value[j];
1710                                     d += drow[Eta_row_nr[j]] * Eta_value[j];
1711                                 }
1712                             f = round(f, Epsel);
1713                             prow[r] = f;
1714                             d = round(d, Epsel);
1715                             drow[r] = d;
1716                         }
1717                     for(i = 1; i <= columns; i++)
1718                         {
1719                             varnr = Rows + i;
1720                             if(Basis[varnr] == FALSE)
1721                                 {
1722                                     d = - Extrad * drow[0];
1723                                     f = 0;
1724                                     for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1725                                         {
1726                                             d = d + drow[Mat[j].row_nr] * Mat[j].value;
1727                                             f = f + prow[Mat[j].row_nr] * Mat[j].value;
1728                                         }
1729                                     drow[varnr] = d;
1730                                     prow[varnr] = f;
1731                                 }
1732                         }
1733                     for(i = 0; i <= Sum; i++)
1734                         {
1735                             prow[i] = round(prow[i], Epsel);
1736                             drow[i] = round(drow[i], Epsd);
1737                         }
1738                 }
1739             if(Rhs[row_nr] > Upbo[Bas[row_nr]])
1740                 g = -1;
1741             else
1742                 g = 1;
1743             pivot = 0;
1744             colnr.value = 0;
1745             theta = Infinite;
1746             for(i = 1; i <= Sum; i++)
1747                 {
1748                     if(Lower[i] != FALSE)
1749                         d = prow[i] * g;
1750                     else
1751                         d = -prow[i] * g;
1752                     if((d < 0) && (Basis[i] == FALSE) && (Upbo[i] > 0))
1753                         {
1754                             if(Lower[i] == FALSE)
1755                                 quot = -drow[i] / (float) d;
1756                             else
1757                                 quot = drow[i] / (float) d;
1758                             if(quot < theta)
1759                                 {
1760                                     theta = quot;
1761                                     pivot = d;
1762                                     colnr.value = i;
1763                                 }
1764                             else if((quot == theta) && (Math.abs(d) > Math.abs(pivot)))
1765                                 {
1766                                     pivot = d;
1767                                     colnr.value = i;
1768                                 }
1769                         }
1770                 }
1771             if(Lp.trace != FALSE)
1772                 System.out.println("col_dual:" + colnr.value + ", pivot element:  " +
1773                                    prow[(int)colnr.value]);
1774
1775             if(colnr.value > 0)
1776                 Doiter = TRUE;
1777
1778             return(colnr.value > 0 ? (short)1 : (short)0);
1779         } /* coldual */
1780
1781         private void iteration(int row_nr,
1782                                int varin,
1783                                Ref theta,
1784                                float up,
1785                                Ref minit,
1786                                Ref low,
1787                                short primal,
1788                                float[] pcol)
1789         {
1790             int  i, k, varout;
1791             float f;
1792             float pivot;
1793
1794             Lp.iter++;
1795             minit.value = theta.value > (up + Epsb) ? 1 : 0;
1796             if(minit.value != 0)
1797                 {
1798                     theta.value = up;
1799                     low.value = low.value == 0 ? 1 : 0;
1800                 }
1801             k = Eta_col_end[Eta_size + 1];
1802             pivot = Eta_value[k - 1];
1803             for(i = Eta_col_end[Eta_size]; i < k; i++)
1804                 {
1805                     f = Rhs[Eta_row_nr[i]] - theta.value * Eta_value[i];
1806                     f = round(f, Epsb);
1807                     Rhs[Eta_row_nr[i]] = f;
1808                 }
1809             if(minit.value == 0)
1810                 {
1811                     Rhs[row_nr] = theta.value;
1812                     varout = Bas[row_nr];
1813                     Bas[row_nr] = varin;
1814                     Basis[varout] = FALSE;
1815                     Basis[varin] = TRUE;
1816                     if(primal != FALSE && pivot < 0)
1817                         Lower[varout] = FALSE;
1818                     if(low.value == 0 && up < Infinite)
1819                         {
1820                             low.value = TRUE;
1821                             Rhs[row_nr] = up - Rhs[row_nr];
1822                             for(i = Eta_col_end[Eta_size]; i < k; i++)
1823                                 Eta_value[i] = -Eta_value[i];
1824                         }
1825                     addetacol();
1826                     Num_inv++;
1827                 }
1828             if(Lp.trace != FALSE)
1829                 {
1830                     System.out.print("Theta = " + theta.value + " ");
1831                     if(minit.value != 0)
1832                         {
1833                             if(Lower[varin] == FALSE)
1834                                 System.out.print("Iteration:" + Lp.iter +
1835                                                  ", variable" + varin + " changed from 0 to its upper Bound of "
1836                                                  + Upbo[varin] + "\n");
1837                             else
1838                                 System.out.print("Iteration:" + Lp.iter + ", variable" + varin +
1839                                                  " changed its upper Bound of " + Upbo[varin] + " to 0\n");
1840                         }
1841                     else
1842                         System.out.print("Iteration:" + Lp.iter + ", variable" + varin +
1843                                          " entered basis at:" + Rhs[row_nr] + "\n");
1844                     if(primal == 0)
1845                         {
1846                             f = 0;
1847                             for(i = 1; i <= Rows; i++)
1848                                 if(Rhs[i] < 0)
1849                                     f -= Rhs[i];
1850                                 else
1851                                     if(Rhs[i] > Upbo[Bas[i]])
1852                                         f += Rhs[i] - Upbo[Bas[i]];
1853                             System.out.println("feasibility gap of this basis:" + (float)f);
1854                         }
1855                     else
1856                         System.out.println("objective function value of this feasible basis: " + Rhs[0]);
1857                 }
1858         } /* iteration */
1859
1860
1861         private int solvelp()
1862         {
1863             int    i, j, varnr;
1864             float   f = 0, theta = 0;
1865             short  primal;
1866             float[]   drow, prow, Pcol;
1867             short  minit;
1868             int    colnr, row_nr;
1869             colnr = 0;
1870             row_nr = 0;
1871             short flag;
1872             Ref ref1, ref2, ref3;
1873             ref1 = new Ref(0);
1874             ref2 = new Ref(0);
1875             ref3 = new Ref(0);
1876
1877             drow = new float[Sum + 1];
1878             prow = new float[Sum + 1];
1879             Pcol = new float[Rows + 1];
1880             for (i = 0; i <= Sum; i++) {
1881                 drow[i] = 0;
1882                 prow[i] = 0;
1883             }
1884             for (i = 0; i <= Rows; i++)
1885                 Pcol[i] = 0;
1886
1887
1888             Lp.iter = 0;
1889             minit = FALSE;
1890             Status = RUNNING;
1891             DoInvert = FALSE;
1892             Doiter = FALSE;
1893             i = 0;
1894             primal = TRUE;
1895             while(i != Rows && primal != FALSE)
1896                 {
1897                     i++;
1898                     primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ?
1899                         (short)1: (short)0;
1900                 }
1901             if(Lp.trace != FALSE)
1902                 {
1903                     if(primal != FALSE)
1904                         System.out.print("Start at feasible basis\n");
1905                     else
1906                         System.out.print("Start at infeasible basis\n");
1907                 }
1908             if(primal == FALSE)
1909                 {
1910                     drow[0] = 1;
1911                     for(i = 1; i <= Rows; i++)
1912                         drow[i] = 0;
1913                     Extrad = 0;
1914                     for(i = 1; i <= columns; i++)
1915                         {
1916                             varnr = Rows + i;
1917                             drow[varnr] = 0;
1918                             for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1919                                 if(drow[Mat[j].row_nr] != 0)
1920                                     drow[varnr] += drow[Mat[j].row_nr] * Mat[j].value;
1921                             if(drow[varnr] < Extrad)
1922                                 Extrad = drow[varnr];
1923                         }
1924                 }
1925             else
1926                 Extrad = 0;
1927             if(Lp.trace != FALSE)
1928                 System.out.println("Extrad = " + Extrad);
1929             minit = FALSE;
1930
1931             while(Status == RUNNING)
1932                 {
1933                     Doiter = FALSE;
1934                     DoInvert = FALSE;
1935
1936                     construct_solution(Solution);
1937                     if(primal != FALSE)
1938                         {
1939                             ref1.value = colnr;
1940                             flag = colprim(ref1, minit, drow);
1941                             colnr = (int)ref1.value;
1942                             if (flag != FALSE)
1943                                 {
1944                                     setpivcol(Lower[colnr], colnr, Pcol);
1945                                     ref1.value = row_nr;
1946                                     ref2.value = theta;
1947                                     flag = rowprim(colnr, ref1, ref2, Pcol);
1948                                     row_nr = (int)ref1.value;
1949                                     theta = ref2.value;
1950                                     if (flag != FALSE)
1951                                         condensecol(row_nr, Pcol);
1952
1953                                 }
1954                         }
1955                     else /* not primal */
1956                         {
1957                             if(minit == FALSE) {
1958                                 ref1.value = row_nr;
1959                                 flag = rowdual(ref1);
1960                                 row_nr = (int)ref1.value;
1961                             }
1962                             if(row_nr > 0 )
1963                                 {
1964                                     ref1.value = colnr;
1965                                     flag = coldual(row_nr, ref1, minit, prow, drow);
1966                                     colnr = (int)ref1.value;
1967                                     if (flag != FALSE)
1968                                         {
1969                                             setpivcol(Lower[colnr], colnr, Pcol);
1970                                             /* getting div by zero here ... MB */
1971                                             if(Pcol[row_nr] == 0)
1972                                                 {
1973                                                     System.err.println("An attempt was made to divide by zero (Pcol[" +
1974                                                                        row_nr + "])");
1975                                                     System.err.println("This indicates numerical instability");
1976                                                     Doiter = FALSE;
1977                                                     if(JustInverted == FALSE)
1978                                                         {
1979                                                             System.out.println("Reinverting Eta");
1980                                                             DoInvert = TRUE;
1981                                                         }
1982                                                     else
1983                                                         {
1984                                                             System.out.println("Can't reinvert, failure");
1985                                                             Status = FAILURE;
1986                                                         }
1987                                                 }
1988                                             else
1989                                                 {
1990                                                     condensecol(row_nr, Pcol);
1991                                                     f = Rhs[row_nr] - Upbo[Bas[row_nr]];
1992                                                     if(f > 0)
1993                                                         {
1994                                                             theta = f / (float) Pcol[row_nr];
1995                                                             if(theta <= Upbo[colnr])
1996                                                                 Lower[Bas[row_nr]] =
1997                                                                     (Lower[Bas[row_nr]] == FALSE)?
1998                                                                     (short)1:(short)0;
1999                                                         }
2000                                                     else /* f <= 0 */
2001                                                         theta = Rhs[row_nr] / (float) Pcol[row_nr];
2002                                                 }
2003                                         }
2004                                     else
2005                                         Status = INFEASIBLE;
2006                                 }
2007                             else
2008                                 {
2009                                     primal   = TRUE;
2010                                     Doiter   = FALSE;
2011                                     Extrad   = 0;
2012                                     DoInvert = TRUE;
2013                                 }
2014                         }
2015                     if(Doiter != FALSE) {
2016                         ref1.value = theta;
2017                         ref2.value = minit;
2018                         ref3.value = Lower[colnr];
2019                         iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol);
2020                         theta = ref1.value;
2021                         minit = (short)ref2.value;
2022                         Lower[colnr] = (short)ref3.value;
2023                     }
2024                     if(Num_inv >= Lp.max_num_inv)
2025                         DoInvert = TRUE;
2026                     if(DoInvert != FALSE)
2027                         {
2028                             invert();
2029                         }
2030                 }
2031
2032             Lp.total_iter += Lp.iter;
2033
2034             return(Status);
2035         } /* solvelp */
2036
2037
2038         private short is_int(float value)
2039         {
2040             float   tmp;
2041
2042             tmp = (float)(value - Math.floor(value));
2043             if(tmp < Epsilon)
2044                 return(TRUE);
2045             if(tmp > (1 - Epsilon))
2046                 return(TRUE);
2047             return(FALSE);
2048         } /* is_int */
2049
2050         private void construct_solution(float[]   sol)
2051         {
2052             int    i, j, basi;
2053             float   f;
2054
2055             /* zero all results of rows */
2056             for (i = 0; i <= Rows; i++)
2057                 sol[i] = 0;
2058
2059             if (Lp.scaling_used != FALSE)
2060                 {
2061                     for(i = Rows + 1; i <= Sum; i++)
2062                         sol[i] = Lowbo[i] * Lp.scale[i];
2063                     for(i = 1; i <= Rows; i++)
2064                         {
2065                             basi = Bas[i];
2066                             if(basi > Rows)
2067                                 sol[basi] += Rhs[i] * Lp.scale[basi];
2068                         }
2069                     for(i = Rows + 1; i <= Sum; i++)
2070                         if(Basis[i] == FALSE && Lower[i] == FALSE)
2071                             sol[i] += Upbo[i] * Lp.scale[i];
2072
2073                     for(j = 1; j <= columns; j++)
2074                         {
2075                             f = sol[Rows + j];
2076                             if(f != 0)
2077                                 for(i = Col_end[j - 1]; i < Col_end[j]; i++)
2078                                     sol[Mat[i].row_nr] += (f / Lp.scale[Rows+j])
2079                                         * (Mat[i].value / Lp.scale[Mat[i].row_nr]);
2080                         }
2081
2082                     for(i = 0; i <= Rows; i++)
2083                         {
2084                             if(Math.abs(sol[i]) < Epsb)
2085                                 sol[i] = 0;
2086                             else
2087                                 if(Lp.ch_sign[i] != FALSE)
2088                                     sol[i] = -sol[i];
2089                         }
2090                 }
2091             else
2092                 {
2093                     for(i = Rows + 1; i <= Sum; i++)
2094                         sol[i] = Lowbo[i];
2095                     for(i = 1; i <= Rows; i++)
2096                         {
2097                             basi = Bas[i];
2098                             if(basi > Rows)
2099                                 sol[basi] += Rhs[i];
2100                         }
2101                     for(i = Rows + 1; i <= Sum; i++)
2102                         if(Basis[i] == FALSE && Lower[i] == FALSE)
2103                             sol[i] += Upbo[i];
2104                     for(j = 1; j <= columns; j++)
2105                         {
2106                             f = sol[Rows + j];
2107                             if(f != 0)
2108                                 for(i = Col_end[j - 1]; i < Col_end[j]; i++)
2109                                     sol[Mat[i].row_nr] += f * Mat[i].value;
2110                         }
2111
2112                     for(i = 0; i <= Rows; i++)
2113                         {
2114                             if(Math.abs(sol[i]) < Epsb)
2115                                 sol[i] = 0;
2116                             else
2117                                 if(Lp.ch_sign[i] != FALSE)
2118                                     sol[i] = -sol[i];
2119                         }
2120                 }
2121         } /* construct_solution */
2122
2123         private void calculate_duals()
2124         {
2125             int i;
2126
2127             /* initialise */
2128             for(i = 1; i <= Rows; i++)
2129                 Lp.duals[i] = 0;
2130             Lp.duals[0] = 1;
2131             btran(Lp.duals);
2132             if (Lp.scaling_used != FALSE)
2133                 for(i = 1; i <= Rows; i++)
2134                     Lp.duals[i] *= Lp.scale[i]/Lp.scale[0];
2135
2136             /* the dual values are the reduced costs of the slacks */
2137             /* When the slack is at its upper Bound, change the sign. Can this happen? */
2138             for(i = 1; i <= Rows; i++)
2139                 {
2140                     if(Lp.basis[i] != FALSE)
2141                         Lp.duals[i] = 0;
2142                     else if( Lp.ch_sign[0] == Lp.ch_sign[i])
2143                         Lp.duals[i] = -Lp.duals[i];
2144                 }
2145         }
2146
2147         private int milpsolve(float[]   upbo,
2148                               float[]   lowbo,
2149                               short[]  sbasis,
2150                               short[]  slower,
2151                               int[]    sbas)
2152         {
2153             int i, j, failure, notint, is_worse;
2154             float theta, tmpfloat;
2155             Random rdm = new Random();
2156             notint = 0;
2157
2158             if(Break_bb != FALSE)
2159                 return(BREAK_BB);
2160             Level++;
2161             Lp.total_nodes++;
2162             if(Level > Lp.max_level)
2163                 Lp.max_level = Level;
2164             /* make fresh copies of upbo, lowbo, rh as solving changes them */
2165             System.arraycopy(upbo, 0, Upbo, 0, Sum + 1);
2166             System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1);
2167             System.arraycopy(sbasis, 0, Basis, 0, Sum + 1);
2168             System.arraycopy(slower, 0, Lower, 0, Sum + 1);
2169             System.arraycopy(sbas, 0, Bas, 0, Rows + 1);
2170             System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1);
2171
2172             if(Lp.anti_degen != FALSE)
2173                 {
2174                     for(i = 1; i <= columns; i++)
2175                         {
2176                             tmpfloat = rdm.nextFloat()*(float)0.001;
2177                             if(tmpfloat > Epsb)
2178                                 Lowbo[i + Rows] -= tmpfloat;
2179                             tmpfloat = rdm.nextFloat()*(float)0.001;
2180                             if(tmpfloat > Epsb)
2181                                 Upbo[i + Rows] += tmpfloat;
2182                         }
2183                     Lp.eta_valid = FALSE;
2184                 }
2185
2186             if(Lp.eta_valid == FALSE)
2187                 {
2188                     for(i = 1; i <= columns; i++)
2189                         if(Lowbo[Rows + i] != 0)
2190                             {
2191                                 theta = Lowbo[ Rows + i];
2192                                 if(Upbo[Rows + i]<Infinite)
2193                                     Upbo[Rows + i] -= theta;
2194                                 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
2195                                     Rh[Mat[j].row_nr] -= theta * Mat[j].value;
2196                             }
2197                     invert();
2198                     Lp.eta_valid = TRUE;
2199                 }
2200
2201             failure = solvelp();
2202
2203             if(Lp.anti_degen != FALSE)
2204                 {
2205                     System.arraycopy(upbo, 0, Upbo, 0, Sum + 1);
2206                     System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1);
2207                     System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1);
2208
2209                     for(i = 1; i <= columns; i++)
2210                         if(Lowbo[Rows + i] != 0)
2211                             {
2212                                 theta = Lowbo[ Rows + i];
2213                                 if(Upbo[Rows + i]<Infinite)
2214                                     Upbo[Rows + i] -= theta;
2215                                 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
2216                                     Rh[Mat[j].row_nr] -= theta * Mat[j].value;
2217                             }
2218                     invert();
2219                     Lp.eta_valid = TRUE;
2220                     failure = solvelp();
2221                 }
2222
2223             if(failure == INFEASIBLE && Lp.verbose != FALSE)
2224                 System.out.print("level" + Level + " INF\n");
2225
2226             if(failure == OPTIMAL)      /* there is a solution */
2227                 {
2228                     construct_solution(Solution);
2229
2230                     /* if this solution is worse than the best sofar, this branch must die */
2231                     if(Maximise != FALSE)
2232                         is_worse = (Solution[0] <= Best_solution[0]) ? 1:0;
2233                     else                        /* minimising! */
2234                         is_worse = (Solution[0] >= Best_solution[0]) ? 1:0;
2235
2236                     if(is_worse != FALSE)
2237                         {
2238                             if(Lp.verbose != FALSE)
2239                                 System.out.println("level" + Level + " OPT NOB value " + Solution[0] +
2240                                                    " Bound " + Best_solution[0]);
2241                             Level--;
2242                             return(MILP_FAIL);
2243                         }
2244
2245                     /* check if solution contains enough ints */
2246                     if(Lp.bb_rule == FIRST_NI)
2247                         {
2248                             notint = 0;
2249                             i = Rows + 1;
2250                             while(i <= Sum && notint == 0)
2251                                 {
2252                                     i++;
2253                                 }
2254                         }
2255                     if(Lp.bb_rule == RAND_NI)
2256                         {
2257                             int nr_not_int, select_not_int;
2258                             nr_not_int = 0;
2259                             for(i = Rows + 1; i <= Sum; i++)
2260                             if(nr_not_int == 0)
2261                                 notint = 0;
2262                             else
2263                                 {
2264                                     select_not_int=(rdm.nextInt() % nr_not_int) + 1;
2265                                     i = Rows + 1;
2266                                     while(select_not_int > 0)
2267                                         {
2268                                             i++;
2269                                         }
2270                                     notint = i - 1;
2271                                 }
2272                         }
2273
2274                     if(Lp.verbose == TRUE)
2275                         if(notint != FALSE)
2276                             System.out.println("level " + Level + " OPT     value " +  Solution[0]);
2277                         else
2278                             System.out.println("level " + Level + " OPT INT value " +  Solution[0]);
2279
2280                     if(notint != FALSE)         /* there is at least one value not yet int */
2281                         {
2282                             /* set up two new problems */
2283                             float[]   new_upbo, new_lowbo;
2284                             float   new_bound;
2285                             short[]  new_lower, new_basis;
2286                             int[]    new_bas;
2287                             int     resone, restwo;
2288
2289                             /* allocate room for them */
2290                             new_upbo = new float[Sum + 1];
2291                             new_lowbo = new float[Sum + 1];
2292                             new_lower = new short[Sum + 1];
2293                             new_basis = new short[Sum + 1];
2294                             new_bas = new int[Rows + 1];
2295                             System.arraycopy(upbo, 0, new_upbo, 0, Sum + 1);
2296                             System.arraycopy(lowbo, 0, new_lowbo, 0, Sum + 1);
2297                             System.arraycopy(Lower, 0, new_lower, 0, Sum + 1);
2298                             System.arraycopy(Basis, 0, new_basis, 0, Sum + 1);
2299                             System.arraycopy(Bas, 0, new_bas, 0, Rows +1);
2300
2301                             if(Floor_first != FALSE)
2302                                 {
2303                                     new_bound = (float)(Math.ceil(Solution[notint]) - 1);
2304                                     /* this Bound might conflict */
2305                                     if(new_bound < lowbo[notint])
2306                                         {
2307                                             resone = MILP_FAIL;
2308                                         }
2309                                     else                /* Bound feasible */
2310                                         {
2311                                             new_upbo[notint] = new_bound;
2312                                             Lp.eta_valid = FALSE;
2313                                             resone = milpsolve(new_upbo, lowbo, new_basis, new_lower,
2314                                                                new_bas);
2315                                             Lp.eta_valid = FALSE;
2316                                         }
2317                                     new_bound += 1;
2318                                     if(new_bound > upbo[notint])
2319                                         {
2320                                             restwo = MILP_FAIL;
2321                                         }
2322                                     else                /* Bound feasible */
2323                                         {
2324                                             new_lowbo[notint] = new_bound;
2325                                             Lp.eta_valid = FALSE;
2326                                             restwo = milpsolve(upbo, new_lowbo, new_basis, new_lower,
2327                                                                new_bas);
2328                                             Lp.eta_valid = FALSE;
2329                                         }
2330                                 }
2331                             else                        /* take ceiling first */
2332                                 {
2333                                     new_bound = (float)Math.ceil(Solution[notint]);
2334                                     /* this Bound might conflict */
2335                                     if(new_bound > upbo[notint])
2336                                         {
2337                                             resone = MILP_FAIL;
2338                                         }
2339                                     else                /* Bound feasible */
2340                                         {
2341                                             new_lowbo[notint] = new_bound;
2342                                             Lp.eta_valid = FALSE;
2343                                             resone = milpsolve(upbo, new_lowbo, new_basis, new_lower,
2344                                                                new_bas);
2345                                             Lp.eta_valid = FALSE;
2346                                         }
2347                                     new_bound -= 1;
2348                                     if(new_bound < lowbo[notint])
2349                                         {
2350                                             restwo = MILP_FAIL;
2351                                         }
2352                                     else                /* Bound feasible */
2353                                         {
2354                                             new_upbo[notint] = new_bound;
2355                                             Lp.eta_valid = FALSE;
2356                                             restwo = milpsolve(new_upbo, lowbo, new_basis, new_lower,
2357                                                                new_bas);
2358                                             Lp.eta_valid = FALSE;
2359                                         }
2360                                 }
2361                             if(resone != FALSE && restwo != FALSE)
2362                                 /* both failed and must have been infeasible */
2363                                 failure = INFEASIBLE;
2364                             else
2365                                 failure = OPTIMAL;
2366
2367                         }
2368                     else /* all required values are int */
2369                         {
2370                             if(Maximise != FALSE)
2371                                 is_worse = (Solution[0] < Best_solution[0]) ? 1:0;
2372                             else
2373                                 is_worse = (Solution[0] > Best_solution[0]) ? 1:0;
2374
2375                             if(is_worse == FALSE) /* Current solution better */
2376                                 {
2377                                     System.arraycopy(Solution, 0, Best_solution, 0, Sum + 1);
2378                                     calculate_duals();
2379                                     if(Lp.break_at_int != FALSE)
2380                                         {
2381                                             if(Maximise != FALSE &&  (Best_solution[0] > Lp.break_value))
2382                                                 Break_bb = TRUE;
2383                                             if(Maximise == FALSE &&  (Best_solution[0] < Lp.break_value))
2384                                                 Break_bb = TRUE;
2385                                         }
2386                                 }
2387                         }
2388                 }
2389
2390             /* failure can have the values OPTIMAL, UNBOUNDED and INFEASIBLE. */
2391             Level--;
2392             return(failure);
2393         } /* milpsolve */
2394
2395         public int solve(Problem lp)
2396         {
2397             int result, i;
2398
2399             if(lp.active == FALSE)
2400                 set_globals(lp);
2401
2402             lp.total_iter  = 0;
2403             lp.max_level   = 1;
2404             lp.total_nodes = 0;
2405
2406             if(Isvalid(lp) != FALSE)
2407                 {
2408                     if(Maximise != FALSE && lp.obj_bound == Infinite)
2409                         Best_solution[0]=-Infinite;
2410                     else if(Maximise == FALSE && lp.obj_bound==-Infinite)
2411                         Best_solution[0] = Infinite;
2412                     else
2413                         Best_solution[0] = lp.obj_bound;
2414
2415                     Level = 0;
2416
2417                     if(lp.basis_valid == FALSE)
2418                         {
2419                             for(i = 0; i <= lp.rows; i++)
2420                                 {
2421                                     lp.basis[i] = TRUE;
2422                                     lp.bas[i] = i;
2423                                 }
2424                             for(i = lp.rows+1; i <= lp.sum; i++)
2425                                 lp.basis[i] = FALSE;
2426                             for(i = 0; i <= lp.sum; i++)
2427                                 lp.lower[i] = TRUE;
2428                             lp.basis_valid = TRUE;
2429                         }
2430
2431                     lp.eta_valid = FALSE;
2432                     Break_bb      = FALSE;
2433                     result        = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas);
2434                     lp.eta_size  = Eta_size;
2435                     lp.eta_alloc = Eta_alloc;
2436                     lp.num_inv   = Num_inv;
2437
2438                     return(result);
2439                 }
2440
2441             /* if we get here, Isvalid(lp) failed. I suggest we return FAILURE */
2442             return(FAILURE);
2443         }
2444
2445         public int lag_solve(Problem lp, float start_bound, int num_iter, short verbose)
2446         {
2447             int i, j, result, citer;
2448             short status, OrigFeas, AnyFeas, same_basis;
2449             float[] OrigObj, ModObj, SubGrad, BestFeasSol;
2450             float Zub, Zlb, Ztmp, pie;
2451             float rhsmod, Step, SqrsumSubGrad;
2452             int[]   old_bas;
2453             short[] old_lower;
2454
2455             /* allocate mem */
2456             OrigObj = new float[lp.columns + 1];
2457             ModObj = new float[lp.columns + 1];
2458             for (i = 0; i <= lp.columns; i++)
2459                 ModObj[i] = 0;
2460             SubGrad = new float[lp.nr_lagrange];
2461             for (i = 0; i < lp.nr_lagrange; i++)
2462                 SubGrad[i] = 0;
2463             BestFeasSol = new float[lp.sum + 1];
2464             for (i = 0; i <= lp.sum; i++)
2465                 BestFeasSol[i] = 0;
2466             old_bas = new int[lp.rows + 1];
2467             System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows + 1);
2468             old_lower = new short[lp.sum + 1];
2469             System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum + 1);
2470
2471             get_row(lp, 0, OrigObj);
2472
2473             pie = 2;
2474
2475             if(lp.maximise != FALSE)
2476                 {
2477                     Zub = DEF_INFINITE;
2478                     Zlb = start_bound;
2479                 }
2480             else
2481                 {
2482                     Zlb = -DEF_INFINITE;
2483                     Zub = start_bound;
2484                 }
2485             status   = RUNNING;
2486             Step     = 1;
2487             OrigFeas = FALSE;
2488             AnyFeas  = FALSE;
2489             citer    = 0;
2490
2491             for(i = 0 ; i < lp.nr_lagrange; i++)
2492                 lp.lambda[i] = 0;
2493
2494             while(status == RUNNING)
2495                 {
2496                     citer++;
2497
2498                     for(i = 1; i <= lp.columns; i++)
2499                         {
2500                             ModObj[i] = OrigObj[i];
2501                             for(j = 0; j < lp.nr_lagrange; j++)
2502                                 {
2503                                     if(lp.maximise != FALSE)
2504                                         ModObj[i] -= lp.lambda[j] * lp.lag_row[j][i];
2505                                     else
2506                                         ModObj[i] += lp.lambda[j] * lp.lag_row[j][i];
2507                                 }
2508                         }
2509                     for(i = 1; i <= lp.columns; i++)
2510                         {
2511                             set_mat(lp, 0, i, ModObj[i]);
2512                             /* fSystem.out.print(stderr, "%f ", ModObj[i]); */
2513                         }
2514                     rhsmod = 0;
2515                     for(i = 0; i < lp.nr_lagrange; i++)
2516                         if(lp.maximise != FALSE)
2517                             rhsmod += lp.lambda[i] * lp.lag_rhs[i];
2518                         else
2519                             rhsmod -= lp.lambda[i] * lp.lag_rhs[i];
2520
2521                     if(verbose != FALSE)
2522                         {
2523                             System.out.println("Zub: " + Zub + " Zlb: " + Zlb + " Step: " + Step +
2524                                                " pie: " + pie + " Feas " + OrigFeas);
2525                             for(i = 0; i < lp.nr_lagrange; i++)
2526                                 System.out.println(i + " SubGrad " + SubGrad[i] + " lambda " + lp.lambda[i]);
2527                         }
2528
2529
2530                     result = solve(lp);
2531
2532                     same_basis = TRUE;
2533                     i = 1;
2534                     while(same_basis != FALSE && i < lp.rows)
2535                         {
2536                             same_basis = (old_bas[i] == lp.bas[i])? (short)1: (short)0;
2537                             i++;
2538                         }
2539                     i = 1;
2540                     while(same_basis != FALSE && i < lp.sum)
2541                         {
2542                             same_basis=(old_lower[i] == lp.lower[i])? (short)1:(short)0;
2543                             i++;
2544                         }
2545                     if(same_basis == FALSE)
2546                         {
2547                             System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum+1);
2548                             System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows+1);
2549                             pie *= 0.95;
2550                         }
2551
2552                     if(verbose != FALSE)
2553                         System.out.println("result: " + result + "  same basis: " + same_basis);
2554
2555                     if(result == UNBOUNDED)
2556                         {
2557                             for(i = 1; i <= lp.columns; i++)
2558                                 System.out.print(ModObj[i] + " ");
2559                             System.exit(FAIL);
2560                         }
2561
2562                     if(result == FAILURE)
2563                         status = FAILURE;
2564
2565                     if(result == INFEASIBLE)
2566                         status = INFEASIBLE;
2567
2568                     SqrsumSubGrad = 0;
2569                     for(i = 0; i < lp.nr_lagrange; i++)
2570                         {
2571                             SubGrad[i]= -lp.lag_rhs[i];
2572                             for(j = 1; j <= lp.columns; j++)
2573                                 SubGrad[i] += lp.best_solution[lp.rows + j] * lp.lag_row[i][j];
2574                             SqrsumSubGrad += SubGrad[i] * SubGrad[i];
2575                         }
2576
2577                     OrigFeas = TRUE;
2578                     for(i = 0; i < lp.nr_lagrange; i++)
2579                         if(lp.lag_con_type[i] != FALSE)
2580                             {
2581                                 if(Math.abs(SubGrad[i]) > lp.epsb)
2582                                     OrigFeas = FALSE;
2583                             }
2584                         else if(SubGrad[i] > lp.epsb)
2585                             OrigFeas = FALSE;
2586
2587                     if(OrigFeas != FALSE)
2588                         {
2589                             AnyFeas = TRUE;
2590                             Ztmp = 0;
2591                             for(i = 1; i <= lp.columns; i++)
2592                                 Ztmp += lp.best_solution[lp.rows + i] * OrigObj[i];
2593                             if((lp.maximise != FALSE) && (Ztmp > Zlb))
2594                                 {
2595                                     Zlb = Ztmp;
2596                                     for(i = 1; i <= lp.sum; i++)
2597                                         BestFeasSol[i] = lp.best_solution[i];
2598                                     BestFeasSol[0] = Zlb;
2599                                     if(verbose != FALSE)
2600                                         System.out.println("Best feasible solution: " + Zlb);
2601                                 }
2602                             else if(Ztmp < Zub)
2603                                 {
2604                                     Zub = Ztmp;
2605                                     for(i = 1; i <= lp.sum; i++)
2606                                         BestFeasSol[i] = lp.best_solution[i];
2607                                     BestFeasSol[0] = Zub;
2608                                     if(verbose != FALSE)
2609                                         System.out.println("Best feasible solution: " + Zub);
2610                                 }
2611                         }
2612
2613                     if(lp.maximise != FALSE)
2614                         Zub = Math.min(Zub, rhsmod + lp.best_solution[0]);
2615                     else
2616                         Zlb = Math.max(Zlb, rhsmod + lp.best_solution[0]);
2617
2618                     if(Math.abs(Zub-Zlb) < (float)0.001)
2619                         {
2620                             status = OPTIMAL;
2621                         }
2622                     Step = (float)(pie * ((1.05*Zub) - Zlb) / SqrsumSubGrad);
2623
2624                     for(i = 0; i < lp.nr_lagrange; i++)
2625                         {
2626                             lp.lambda[i] += Step * SubGrad[i];
2627                             if(lp.lag_con_type[i] == FALSE && lp.lambda[i] < 0)
2628                                 lp.lambda[i] = 0;
2629                         }
2630
2631                     if(citer == num_iter && status==RUNNING)
2632                         if(AnyFeas != FALSE)
2633                             status = FEAS_FOUND;
2634                         else
2635                             status = NO_FEAS_FOUND;
2636                 }
2637
2638             for(i = 0; i <= lp.sum; i++)
2639                 lp.best_solution[i] = BestFeasSol[i];
2640
2641             for(i = 1; i <= lp.columns; i++)
2642                 set_mat(lp, 0, i, OrigObj[i]);
2643
2644             if(lp.maximise != FALSE)
2645                 lp.lag_bound = Zub;
2646             else
2647                 lp.lag_bound = Zlb;
2648             return(status);
2649         }
2650
2651     } // end of class solve
2652
2653
2654     public static class Matrix {
2655         public int row_nr;
2656         public float value;
2657         public Matrix(int r, float v) { row_nr = r; value = v; }
2658     }
2659
2660
2661     public static class Problem {
2662         String   lp_name;               /* the name of the lp */
2663
2664         public short active;            /*TRUE if the globals point to this structure*/
2665         public short verbose;         /* ## Verbose flag */
2666         public short debug;           /* ## Print B&B information */
2667         public short trace;           /* ## Print information on pivot selection */
2668         public short anti_degen;        /* ## Do perturbations */
2669
2670         public int          rows;               /* Nr of constraint rows in the problem */
2671         int       rows_alloc;           /* The allocated memory for Rows sized data */
2672         int       columns;            /* The number of columns (= variables) */
2673         int       columns_alloc;
2674         int       sum;                /* The size of the variables + the slacks */
2675         int       sum_alloc;
2676
2677         short     names_used;         /* Flag to indecate if names for rows and
2678                                          columns are used */
2679         String[]  row_name;             /* rows_alloc+1 */
2680         int[]  col_name;                /* columns_alloc+1 */
2681
2682         /* Row[0] of the sparce matrix is the objective function */
2683
2684         int       non_zeros;          /* The number of elements in the sparce matrix*/
2685         int       mat_alloc;            /* The allocated size for matrix sized
2686                                            structures */
2687         Matrix[]  mat;                /* mat_alloc :The sparse matrix */
2688         Matrix[]  alternate_mat;                /* mat_alloc :The sparse matrix */
2689         int[]     col_end;            /* columns_alloc+1 :Cend[i] is the index of the
2690                                          first element after column i.
2691                                          column[i] is stored in elements
2692                                          col_end[i-1] to col_end[i]-1 */
2693         int[]     col_no;             /* mat_alloc :From Row 1 on, col_no contains the
2694                                          column nr. of the
2695                                          nonzero elements, row by row */
2696         short     row_end_valid;        /* true if row_end & col_no are valid */
2697         int[]     row_end;            /* rows_alloc+1 :row_end[i] is the index of the
2698                                          first element in Colno after row i */
2699         float[]  orig_rh;            /* rows_alloc+1 :The RHS after scaling & sign
2700                                          changing, but before `Bound transformation' */
2701         float[]  rh;                    /* rows_alloc+1 :As orig_rh, but after Bound
2702                                            transformation */
2703         float[]  rhs;           /* rows_alloc+1 :The RHS of the curent simplex
2704                                    tableau */
2705         float[]  orig_upbo;          /* sum_alloc+1 :Bound before transformations */
2706         float[]  orig_lowbo;            /*  "       "                   */
2707         float[]  upbo;               /*  "       "  :Upper bound after transformation
2708                                           & B&B work*/
2709         float[]  lowbo;              /*  "       "  :Lower bound after transformation
2710                                           & B&B work */
2711
2712         short     basis_valid;        /* TRUE is the basis is still valid */
2713         int[]     bas;                /* rows_alloc+1 :The basis column list */
2714         short[]   basis;              /* sum_alloc+1 : basis[i] is TRUE if the column
2715                                          is in the basis */
2716         short[]   lower;              /*  "       "  :TRUE is the variable is at its
2717                                           lower bound (or in the basis), it is FALSE
2718                                           if the variable is at its upper bound */
2719
2720         short     eta_valid;          /* TRUE if current Eta structures are valid */
2721         int       eta_alloc;          /* The allocated memory for Eta */
2722         int       eta_size;           /* The number of Eta columns */
2723         int       num_inv;            /* The number of float pivots */
2724         int       max_num_inv;        /* ## The number of float pivots between
2725                                          reinvertions */
2726         float[]  eta_value;          /* eta_alloc :The Structure containing the
2727                                          values of Eta */
2728         int[]     eta_row_nr;         /*  "     "  :The Structure containing the Row
2729                                           indexes of Eta */
2730         int[]     eta_col_end;        /* rows_alloc + MaxNumInv : eta_col_end[i] is
2731                                          the start index of the next Eta column */
2732
2733         short       bb_rule;            /* what rule for selecting B&B variables */
2734
2735         short     break_at_int;       /* TRUE if stop at first integer better than
2736                                          break_value */
2737         float    break_value;
2738
2739         float    obj_bound;          /* ## Objective function bound for speedup of
2740                                          B&B */
2741         int       iter;               /* The number of iterations in the simplex
2742                                          solver (LP) */
2743         int       total_iter;         /* The total number of iterations (B&B) (ILP)*/
2744         int       max_level;          /* The Deepest B&B level of the last solution */
2745         int         total_nodes;        /* total number of nodes processed in b&b */
2746         public float[]  solution;           /* sum_alloc+1 :The Solution of the last LP,
2747                                          0 = The Optimal Value,
2748                                          1..rows The Slacks,
2749                                          rows+1..sum The Variables */
2750         public float[]  best_solution;      /*  "       "  :The Best 'Integer' Solution */
2751         float[]  duals;              /* rows_alloc+1 :The dual variables of the
2752                                          last LP */
2753
2754         short     maximise;           /* TRUE if the goal is to maximise the
2755                                          objective function */
2756         short     floor_first;        /* TRUE if B&B does floor bound first */
2757         short[]   ch_sign;            /* rows_alloc+1 :TRUE if the Row in the matrix
2758                                          has changed sign
2759                                          (a`x > b, x>=0) is translated to
2760                                          s + -a`x = -b with x>=0, s>=0) */
2761
2762         short     scaling_used; /* TRUE if scaling is used */
2763         short     columns_scaled;     /* TRUE is the columns are scaled too, Only use
2764                                          if all variables are non-integer */
2765         float[]  scale;              /* sum_alloc+1 :0..Rows the scaling of the Rows,
2766                                          Rows+1..Sum the scaling of the columns */
2767
2768         int         nr_lagrange;        /* Nr. of Langrangian relaxation constraints */
2769         float[][]lag_row;               /* NumLagrange, columns+1:Pointer to pointer of
2770                                            rows */
2771         float[]  lag_rhs;               /* NumLagrange :Pointer to pointer of Rhs */
2772         float[]  lambda;                /* NumLagrange :Lambda Values */
2773         short[]   lag_con_type;       /* NumLagrange :TRUE if constraint type EQ */
2774         float    lag_bound;             /* the lagrangian lower bound */
2775
2776         short     valid;                /* Has this lp pased the 'test' */
2777         float    infinite;           /* ## numercal stuff */
2778         float    epsilon;            /* ## */
2779         float    epsb;               /* ## */
2780         float    epsd;               /* ## */
2781         float    epsel;              /* ## */
2782
2783
2784         public Problem (int nrows, int ncolumns, int matalloc) {
2785             int i, nsum;
2786             nsum=nrows+ncolumns;
2787             rows=nrows;
2788             columns=ncolumns;
2789             sum=nsum;
2790             rows_alloc=rows;
2791             columns_alloc=columns;
2792             sum_alloc=sum;
2793             mat_alloc=matalloc;
2794             eta_alloc=10000;
2795             max_num_inv=DEFNUMINV;
2796             col_no = new int[mat_alloc];
2797             col_end = new int[columns + 1];
2798             row_end = new int[rows + 1];
2799             orig_rh = new float[rows + 1];
2800             rh = new float[rows + 1];
2801             rhs = new float[rows + 1];
2802             orig_upbo = new float[sum + 1];
2803             upbo = new float[sum + 1];
2804             orig_lowbo = new float[sum + 1];
2805             lowbo = new float[sum + 1];
2806             bas = new int[rows+1];
2807             basis = new short[sum + 1];
2808             lower = new short[sum + 1];
2809             eta_value = new float[eta_alloc];
2810             eta_row_nr = new int[eta_alloc];
2811             eta_col_end = new int[rows_alloc + max_num_inv];
2812             solution = new float[sum + 1];
2813             best_solution = new float[sum + 1];
2814             duals = new float[rows + 1];
2815             ch_sign = new short[rows + 1];
2816             mat = new Matrix[mat_alloc];
2817             for(int j=0; j<mat.length; j++) mat[j] = new Matrix(0, 0);
2818             alternate_mat = new Matrix[mat_alloc];
2819             for(int j=0; j<alternate_mat.length; j++) alternate_mat[j] = new Matrix(0, 0);
2820         }
2821
2822         public void init(int nrows, int ncolumns) {
2823             int i, nsum;
2824             nsum=nrows+ncolumns;
2825             lp_name = new String("unnamed");
2826             active=FALSE;
2827             verbose=FALSE;
2828             debug=FALSE;
2829             trace=FALSE;
2830             rows=nrows;
2831             columns=ncolumns;
2832             sum=nsum;
2833             names_used=FALSE;
2834             obj_bound=DEF_INFINITE;
2835             infinite=DEF_INFINITE;
2836             epsilon=DEF_EPSILON;
2837             epsb=DEF_EPSB;
2838             epsd=DEF_EPSD;
2839             epsel=DEF_EPSEL;
2840             non_zeros=0;
2841
2842             for (i = 0; i < mat_alloc; i++) { mat[i].row_nr = 0; mat[i].value = 0; }
2843             for (i = 0; i < mat_alloc; i++)   col_no[i] = 0;
2844             for (i = 0; i < columns + 1; i++) col_end[i] = 0;
2845             for (i = 0; i < rows + 1; i++)    row_end[i] = 0;
2846             for (i = 0; i < rows + 1; i++)   orig_rh[i] = 0;
2847             for (i = 0; i < rows + 1; i++)   rh[i] = 0;
2848             for (i = 0; i < rows + 1; i++)   rhs[i] = 0;
2849             for (i = 0; i <= sum; i++)       orig_upbo[i]=infinite;
2850             for (i = 0; i < sum + 1; i++)    upbo[i] = 0;
2851             for (i = 0; i < sum + 1; i++)    orig_lowbo[i] = 0;
2852             for (i = 0; i < sum + 1; i++)    lowbo[i] = 0;
2853             for (i = 0; i <= rows; i++)      bas[i] = 0;
2854             for (i = 0; i <= sum; i++)       basis[i] = 0;
2855             for (i = 0; i <= rows; i++)     { bas[i]=i; basis[i]=TRUE; }
2856             for (i = rows + 1; i <= sum; i++) basis[i]=FALSE;
2857             for (i = 0 ; i <= sum; i++)       lower[i]=TRUE;
2858             for (i = 0; i < eta_alloc; i++) eta_value[i] = 0;
2859             for (i = 0; i < eta_alloc; i++) eta_row_nr[i] = 0;
2860             for (i = 0; i < rows_alloc + max_num_inv; i++) eta_col_end[i] = 0;
2861             for (i = 0; i <= sum; i++) solution[i] = 0;
2862             for (i = 0; i <= sum; i++) best_solution[i] = 0;
2863             for (i = 0; i <= rows; i++) duals[i] = 0;
2864             for (i = 0; i <= rows; i++) ch_sign[i] = FALSE;
2865
2866             row_end_valid=FALSE;
2867             bb_rule=FIRST_NI;
2868             break_at_int=FALSE;
2869             break_value=0;
2870             iter=0;
2871             total_iter=0;
2872             basis_valid=TRUE;
2873             eta_valid=TRUE;
2874             eta_size=0;
2875             nr_lagrange=0;
2876             maximise = FALSE;
2877             floor_first = TRUE;
2878             scaling_used = FALSE;
2879             columns_scaled = FALSE;
2880             valid = FALSE;
2881         }
2882
2883     }
2884
2885     private final static float round( float val, float eps) { return (Math.abs(val) < eps) ? 0 : val; }
2886
2887 }