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