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 ETA_START_SIZE = 10000; /* start size of array Eta. Realloced if needed */
  48
  49     static class Ref {
  50         float value;
  51         public Ref(float v) { value = v; }
  52     }
  53
  54     public static class Simplex {
  55         /* Globals used by solver */
  56         short JustInverted;
  57         short Status;
  58         short Doiter;
  59         short DoInvert;
  60         short Break_bb;
  61
  62         public short active;            /*TRUE if the globals point to this structure*/
  63         public short debug;           /* ## Print B&B information */
  64         public short trace;           /* ## Print information on pivot selection */
  65         public int          rows;               /* Nr of constraint rows in the problem */
  66         int       rows_alloc;           /* The allocated memory for Rows sized data */
  67         int       columns_alloc;
  68         int       sum;                /* The size of the variables + the slacks */
  69         int       sum_alloc;
  70         int       non_zeros;          /* The number of elements in the sparce matrix*/
  71         int       mat_alloc;            /* The allocated size for matrix sized
  72                                            structures */
  73         MatrixArray  mat;                /* mat_alloc :The sparse matrix */
  74         MatrixArray  alternate_mat;                /* mat_alloc :The sparse matrix */
  75         int[]     col_end;            /* columns_alloc+1 :Cend[i] is the index of the
  76                                          first element after column i.
  77                                          column[i] is stored in elements
  78                                          col_end[i-1] to col_end[i]-1 */
  79         int[]     col_no;             /* mat_alloc :From Row 1 on, col_no contains the
  80                                          column nr. of the
  81                                          nonzero elements, row by row */
  82         short     row_end_valid;        /* true if row_end & col_no are valid */
  83         int[]     row_end;            /* rows_alloc+1 :row_end[i] is the index of the
  84                                          first element in Colno after row i */
  85         float[]  orig_rh;            /* rows_alloc+1 :The RHS after scaling & sign
  86                                          changing, but before `Bound transformation' */
  87         float[]  rh;                    /* rows_alloc+1 :As orig_rh, but after Bound
  88                                            transformation */
  89         float[]  rhs;           /* rows_alloc+1 :The RHS of the curent simplex
  90                                    tableau */
  91         float[]  orig_upbo;          /* sum_alloc+1 :Bound before transformations */
  92         float[]  orig_lowbo;            /*  "       "                   */
  93         float[]  upbo;               /*  "       "  :Upper bound after transformation
  94                                           & B&B work*/
  95         float[]  lowbo;              /*  "       "  :Lower bound after transformation
  96                                           & B&B work */
  97
  98         short     basis_valid;        /* TRUE is the basis is still valid */
  99         int[]     bas;                /* rows_alloc+1 :The basis column list */
 100         short[]   basis;              /* sum_alloc+1 : basis[i] is TRUE if the column
 101                                          is in the basis */
 102         short[]   lower;              /*  "       "  :TRUE is the variable is at its
 103                                           lower bound (or in the basis), it is FALSE
 104                                           if the variable is at its upper bound */
 105
 106         short     eta_valid;          /* TRUE if current Eta structures are valid */
 107         int       eta_alloc;          /* The allocated memory for Eta */
 108         int       eta_size;           /* The number of Eta columns */
 109         int       num_inv;            /* The number of float pivots */
 110         int       max_num_inv;        /* ## The number of float pivots between
 111                                          reinvertions */
 112         float[]  eta_value;          /* eta_alloc :The Structure containing the
 113                                          values of Eta */
 114         int[]     eta_row_nr;         /*  "     "  :The Structure containing the Row
 115                                           indexes of Eta */
 116         int[]     eta_col_end;        /* rows_alloc + MaxNumInv : eta_col_end[i] is
 117                                          the start index of the next Eta column */
 118
 119         short       bb_rule;            /* what rule for selecting B&B variables */
 120
 121         short     break_at_int;       /* TRUE if stop at first integer better than
 122                                          break_value */
 123         float    break_value;
 124
 125         float    obj_bound;          /* ## Objective function bound for speedup of
 126                                          B&B */
 127         int       iter;               /* The number of iterations in the simplex
 128                                          solver () */
 129         int       total_iter;         /* The total number of iterations (B&B) (ILP)*/
 130         int       max_level;          /* The Deepest B&B level of the last solution */
 131         int         total_nodes;        /* total number of nodes processed in b&b */
 132         public float[]  solution;           /* sum_alloc+1 :The Solution of the last LP,
 133                                          0 = The Optimal Value,
 134                                          1..rows The Slacks,
 135                                          rows+1..sum The Variables */
 136         public float[]  best_solution;      /*  "       "  :The Best 'Integer' Solution */
 137         float[]  duals;              /* rows_alloc+1 :The dual variables of the
 138                                          last LP */
 139
 140         short     maximise;           /* TRUE if the goal is to maximise the
 141                                          objective function */
 142         short     floor_first;        /* TRUE if B&B does floor bound first */
 143         short[]   ch_sign;            /* rows_alloc+1 :TRUE if the Row in the matrix
 144                                          has changed sign
 145                                          (a`x > b, x>=0) is translated to
 146                                          s + -a`x = -b with x>=0, s>=0) */
 147
 148         int         nr_lagrange;        /* Nr. of Langrangian relaxation constraints */
 149         float[][]lag_row;               /* NumLagrange, columns+1:Pointer to pointer of
 150                                            rows */
 151         float[]  lag_rhs;               /* NumLagrange :Pointer to pointer of Rhs */
 152         float[]  lambda;                /* NumLagrange :Lambda Values */
 153         short[]   lag_con_type;       /* NumLagrange :TRUE if constraint type EQ */
 154         float    lag_bound;             /* the lagrangian lower bound */
 155
 156         short     valid;                /* Has this lp pased the 'test' */
 157         float    infinite;           /* ## numercal stuff */
 158         float    epsilon;            /* ## */
 159         float    epsb;               /* ## */
 160         float    epsd;               /* ## */
 161         float    epsel;              /* ## */
 162
 163         int     Rows;
 164         int     columns;
 165         int     Sum;
 166         int     Non_zeros;
 167         int     Level;
 168         MatrixArray  Mat;
 169         int[]     Col_no;
 170         int[]     Col_end;
 171         int[]     Row_end;
 172         float[]    Orig_rh;
 173         float[]    Rh;
 174         float[]    Rhs;
 175         float[]    Orig_upbo;
 176         float[]    Orig_lowbo;
 177         float[]    Upbo;
 178         float[]    Lowbo;
 179         int[]     Bas;
 180         short[]   Basis;
 181         short[]   Lower;
 182         int     Eta_alloc;
 183         int     Eta_size;
 184         float[]    Eta_value;
 185         int[]     Eta_row_nr;
 186         int[]     Eta_col_end;
 187         int     Num_inv;
 188         float[]    Solution;
 189         float[]    Best_solution;
 190         float    Infinite;
 191         float    Epsilon;
 192         float    Epsb;
 193         float    Epsd;
 194         float    Epsel;
 195
 196         float   TREJ;
 197         float   TINV;
 198
 199         short   Maximise;
 200         short   Floor_first;
 201         float    Extrad;
 202
 203         int     Warn_count; /* used in CHECK version of rounding macro */
 204
 205         public Simplex (int nrows, int ncolumns, int matalloc) {
 206             int nsum;
 207             nsum=nrows+ncolumns;
 208             rows=nrows;
 209             columns=ncolumns;
 210             sum=nsum;
 211             rows_alloc=rows;
 212             columns_alloc=columns;
 213             sum_alloc=sum;
 214             mat_alloc=matalloc;
 215             eta_alloc=10000;
 216             max_num_inv=DEFNUMINV;
 217             col_no = new int[mat_alloc];
 218             col_end = new int[columns + 1];
 219             row_end = new int[rows + 1];
 220             orig_rh = new float[rows + 1];
 221             rh = new float[rows + 1];
 222             rhs = new float[rows + 1];
 223             orig_upbo = new float[sum + 1];
 224             upbo = new float[sum + 1];
 225             orig_lowbo = new float[sum + 1];
 226             lowbo = new float[sum + 1];
 227             bas = new int[rows+1];
 228             basis = new short[sum + 1];
 229             lower = new short[sum + 1];
 230             eta_value = new float[eta_alloc];
 231             eta_row_nr = new int[eta_alloc];
 232             eta_col_end = new int[rows_alloc + max_num_inv];
 233             solution = new float[sum + 1];
 234             best_solution = new float[sum + 1];
 235             duals = new float[rows + 1];
 236             ch_sign = new short[rows + 1];
 237             mat = new MatrixArray(mat_alloc);
 238             alternate_mat = new MatrixArray(mat_alloc);
 239         }
 240
 241         public void init(int nrows, int ncolumns) {
 242             int nsum;
 243             nsum=nrows+ncolumns;
 244             active=FALSE;
 245             debug=FALSE;
 246             trace=FALSE;
 247             rows=nrows;
 248             columns=ncolumns;
 249             sum=nsum;
 250             obj_bound=DEF_INFINITE;
 251             infinite=DEF_INFINITE;
 252             epsilon=DEF_EPSILON;
 253             epsb=DEF_EPSB;
 254             epsd=DEF_EPSD;
 255             epsel=DEF_EPSEL;
 256             non_zeros=0;
 257
 258             for(int i = 0; i < mat_alloc; i++) { set_row_nr(mat,i, 0); set_value(mat, i, 0); }
 259             for(int i = 0; i < mat_alloc; i++)   col_no[i] = 0;
 260             for(int i = 0; i < columns + 1; i++) col_end[i] = 0;
 261             for(int i = 0; i < rows + 1; i++)    row_end[i] = 0;
 262             for(int i = 0; i < rows + 1; i++)   orig_rh[i] = 0;
 263             for(int i = 0; i < rows + 1; i++)   rh[i] = 0;
 264             for(int i = 0; i < rows + 1; i++)   rhs[i] = 0;
 265             for(int i = 0; i <= sum; i++)       orig_upbo[i]=infinite;
 266             for(int i = 0; i < sum + 1; i++)    upbo[i] = 0;
 267             for(int i = 0; i < sum + 1; i++)    orig_lowbo[i] = 0;
 268             for(int i = 0; i < sum + 1; i++)    lowbo[i] = 0;
 269             for(int i = 0; i <= rows; i++)      bas[i] = 0;
 270             for(int i = 0; i <= sum; i++)       basis[i] = 0;
 271             for(int i = 0; i <= rows; i++)     { bas[i]=i; basis[i]=TRUE; }
 272             for(int i = rows + 1; i <= sum; i++) basis[i]=FALSE;
 273             for(int i = 0 ; i <= sum; i++)       lower[i]=TRUE;
 274             for(int i = 0; i < eta_alloc; i++) eta_value[i] = 0;
 275             for(int i = 0; i < eta_alloc; i++) eta_row_nr[i] = 0;
 276             for(int i = 0; i < rows_alloc + max_num_inv; i++) eta_col_end[i] = 0;
 277             for(int i = 0; i <= sum; i++) solution[i] = 0;
 278             for(int i = 0; i <= sum; i++) best_solution[i] = 0;
 279             for(int i = 0; i <= rows; i++) duals[i] = 0;
 280             for(int i = 0; i <= rows; i++) ch_sign[i] = FALSE;
 281
 282             row_end_valid=FALSE;
 283             bb_rule=FIRST_NI;
 284             break_at_int=FALSE;
 285             break_value=0;
 286             iter=0;
 287             total_iter=0;
 288             basis_valid=TRUE;
 289             eta_valid=TRUE;
 290             eta_size=0;
 291             nr_lagrange=0;
 292             maximise = FALSE;
 293             floor_first = TRUE;
 294             valid = FALSE;
 295         }
 296
 297         public void setObjective(float[] row, boolean maximize) {
 298             for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
 299             row[0] = (float)0.0;
 300             for(int j = 1; j <= columns; j++) {
 301                 int Row = 0;
 302                 int column = j;
 303                 float Value = row[j];
 304                 int elmnr, lastelm;
 305
 306                 if(Row > rows || Row < 0) throw new Error("row out of range");
 307                 if(column > columns || column < 1) throw new Error("column out of range");
 308
 309                 if (basis[column] == TRUE && Row > 0) basis_valid = FALSE;
 310                 eta_valid = FALSE;
 311                 elmnr = col_end[column-1];
 312                 while((elmnr < col_end[column]) ? (get_row_nr(mat, elmnr) != Row) : false) elmnr++;
 313                 if((elmnr != col_end[column]) ? (get_row_nr(mat, elmnr) == Row) : false ) {
 314                     if (ch_sign[Row] != FALSE) set_value(mat, elmnr, -Value);
 315                     else set_value(mat, elmnr, Value);
 316                 } else {
 317                     /* check if more space is needed for matrix */
 318                     if (non_zeros + 1 > mat_alloc) throw new Error("not enough mat space; this should not happen");
 319                     /* Shift the matrix */
 320                     lastelm=non_zeros;
 321                     for(int i = lastelm; i > elmnr ; i--) {
 322                         set_row_nr(mat,i,get_row_nr(mat,i-1));
 323                         set_value(mat,i,get_value(mat,i-1));
 324                     }
 325                     for(int i = column; i <= columns; i++) col_end[i]++;
 326                     /* Set new element */
 327                     set_row_nr(mat,elmnr, Row);
 328                     if (ch_sign[Row] != FALSE) set_value(mat, elmnr, -Value);
 329                     else set_value(mat, elmnr, Value);
 330                     row_end_valid=FALSE;
 331                     non_zeros++;
 332                     if (active != FALSE) Non_zeros=non_zeros;
 333                 }
 334             }
 335             if (maximize) {
 336                 if (maximise == FALSE) {
 337                     for(int i = 0; i < non_zeros; i++)
 338                         if(get_row_nr(mat, i)==0)
 339                             set_value(mat, i, get_value(mat,i)* (float)-1.0);
 340                     eta_valid=FALSE;
 341                 }
 342                 maximise=TRUE;
 343                 ch_sign[0]=TRUE;
 344                 if (active != FALSE) Maximise=TRUE;
 345             } else {
 346                 if (maximise==TRUE) {
 347                     for(int i = 0; i < non_zeros; i++)
 348                         if(get_row_nr(mat, i)==0)
 349                             set_value(mat, i, get_value(mat,i) * (float)-1.0);
 350                     eta_valid=FALSE;
 351                 }
 352                 maximise=FALSE;
 353                 ch_sign[0]=FALSE;
 354                 if (active != FALSE) Maximise=FALSE;
 355             }
 356         }
 357
 358         public void add_constraint(float[] row, short constr_type, float rh) {
 359             for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
 360             row[0] = (float)0.0;
 361
 362             MatrixArray newmat;
 363             int  elmnr;
 364             int  stcol;
 365
 366             newmat = alternate_mat;
 367             for(int i = 0; i < non_zeros; i++) { set_row_nr(newmat,i, 0); set_value(newmat, i, 0); }
 368             for(int i = 1; i <= columns; i++) if (row[i]!=0) non_zeros++;
 369             if (non_zeros > mat_alloc) throw new Error("not enough mat space; this should not happen");
 370             rows++;
 371             sum++;
 372             if(rows > rows_alloc) throw new Error("not enough rows; this should never happen");
 373             if(constr_type==GE) ch_sign[rows] = TRUE;
 374             else ch_sign[rows] = FALSE;
 375
 376             elmnr = 0;
 377             stcol = 0;
 378             for(int i = 1; i <= columns; i++) {
 379                 for(int j = stcol; j < col_end[i]; j++) {
 380                     set_row_nr(newmat,elmnr, get_row_nr(mat, j));
 381                     set_value(newmat, elmnr, get_value(mat,j));
 382                     elmnr++;
 383                 }
 384                 if(((i>=1 && i< columns && row[i]!=0)?TRUE:FALSE) != FALSE) {
 385                     if(ch_sign[rows] != FALSE) set_value(newmat, elmnr, -row[i]);
 386                     else set_value(newmat, elmnr, row[i]);
 387                     set_row_nr(newmat,elmnr, rows);
 388                     elmnr++;
 389                 }
 390                 stcol=col_end[i];
 391                 col_end[i]=elmnr;
 392             }
 393
 394             alternate_mat = mat;
 395             mat = newmat;
 396
 397             for(int i = sum ; i > rows; i--) {
 398                 orig_upbo[i]=orig_upbo[i-1];
 399                 orig_lowbo[i]=orig_lowbo[i-1];
 400                 basis[i]=basis[i-1];
 401                 lower[i]=lower[i-1];
 402             }
 403
 404             for(int i =  1 ; i <= rows; i++) if(bas[i] >= rows) bas[i]++;
 405
 406             if(constr_type==LE || constr_type==GE) orig_upbo[rows]=infinite;
 407             else if(constr_type==EQ) orig_upbo[rows]=0;
 408             else throw new Error("Wrong constraint type\n");
 409             orig_lowbo[rows]=0;
 410
 411             if(constr_type==GE && rh != 0) orig_rh[rows]=-rh;
 412             else orig_rh[rows]=rh;
 413
 414             row_end_valid=FALSE;
 415
 416             bas[rows]=rows;
 417             basis[rows]=TRUE;
 418             lower[rows]=TRUE;
 419
 420             if (active != FALSE) set_globals();
 421             eta_valid=FALSE;
 422         }
 423
 424         public void bound_sum(int column1, int column2, float bound, short type, float[] scratch) {
 425             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 426             scratch[column1] = (float)1.0;
 427             scratch[column2] = (float)1.0;
 428             add_constraint(scratch, type, bound);
 429             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 430         }
 431
 432         public void bound_difference(int column1, int column2, float bound, short type, float[] scratch) {
 433             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 434             scratch[column1] = (float)1.0;
 435             scratch[column2] = (float)-1.0;
 436             add_constraint(scratch, type, bound);
 437             for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
 438         }
 439
 440         public void set_upbo(int column, float value) {
 441             if(column > columns || column < 1) throw new Error("column out of range");
 442             if(value < orig_lowbo[rows + column]) throw new Error("UpperBound must be >= lowerBound");
 443             eta_valid = FALSE;
 444             orig_upbo[rows+column] = value;
 445         }
 446
 447         public void set_lowbo(int column, float value) {
 448             if(column > columns || column < 1) throw new Error("column out of range");
 449             if(value > orig_upbo[rows + column]) throw new Error("UpperBound must be >= lowerBound");
 450             eta_valid = FALSE;
 451             orig_lowbo[rows+column] = value;
 452         }
 453
 454         public void set_rh(int row, float value) {
 455             if(row > rows || row < 0) throw new Error("Row out of Range");
 456             if(row == 0) throw new Error("Warning: attempt to set RHS of objective function, ignored");
 457             if (ch_sign[row] != FALSE) orig_rh[row] = -value;
 458             else orig_rh[row] = value;
 459             eta_valid = FALSE;
 460         }
 461
 462         public void set_rh_vec(float[] rh) {
 463             for(int i=1; i <= rows; i++)
 464                 if (ch_sign[i] != FALSE) orig_rh[i]=-rh[i];
 465                 else orig_rh[i]=rh[i];
 466             eta_valid=FALSE;
 467         }
 468
 469
 470         public void set_constr_type(int row, short con_type) {
 471             if (row > rows || row < 1) throw new Error("Row out of Range");
 472             switch(con_type) {
 473                 case EQ:
 474                     orig_upbo[row]=0;
 475                     basis_valid=FALSE;
 476                     if (ch_sign[row] != FALSE) {
 477                         for(int i = 0; i < non_zeros; i++)
 478                             if (get_row_nr(mat, i)==row) set_value(mat, i, get_value(mat,i) * (float)-1);
 479                         eta_valid=FALSE;
 480                         ch_sign[row]=FALSE;
 481                         if (orig_rh[row]!=0) orig_rh[row]*=-1;
 482                     }
 483                     break;
 484                 case LE:
 485                     orig_upbo[row]=infinite;
 486                     basis_valid=FALSE;
 487                     if (ch_sign[row] != FALSE) {
 488                         for(int i = 0; i < non_zeros; i++)
 489                             if (get_row_nr(mat, i)==row) set_value(mat, i, get_value(mat,i) * (float)-1);
 490                         eta_valid=FALSE;
 491                         ch_sign[row]=FALSE;
 492                         if (orig_rh[row]!=0) orig_rh[row]*=-1;
 493                     }
 494                     break;
 495                 case GE:
 496                     orig_upbo[row]=infinite;
 497                     basis_valid=FALSE;
 498                     if (ch_sign[row] == FALSE) {
 499                         for(int i = 0; i < non_zeros; i++)
 500                             if (get_row_nr(mat, i)==row) set_value(mat, i, get_value(mat,i) * (float)-1);
 501                         eta_valid=FALSE;
 502                         ch_sign[row]=TRUE;
 503                         if (orig_rh[row]!=0) orig_rh[row]*=-1;
 504                     }
 505                     break;
 506                 default: throw new Error("Constraint type not (yet) implemented");
 507             }
 508         }
 509
 510         void set_globals() {
 511             Rows = rows;
 512             columns = columns;
 513             Sum = Rows + columns;
 514             Non_zeros = non_zeros;
 515             Mat = mat;
 516             Col_no = col_no;
 517             Col_end = col_end;
 518             Row_end = row_end;
 519             Rh = rh;
 520             Rhs = rhs;
 521             Orig_rh = orig_rh;
 522             Orig_upbo = orig_upbo;
 523             Orig_lowbo = orig_lowbo;
 524             Upbo = upbo;
 525             Lowbo = lowbo;
 526             Bas = bas;
 527             Basis = basis;
 528             Lower = lower;
 529             Eta_alloc = eta_alloc;
 530             Eta_size = eta_size;
 531             Num_inv = num_inv;
 532             Eta_value = eta_value;
 533             Eta_row_nr = eta_row_nr;
 534             Eta_col_end = eta_col_end;
 535             Solution = solution;
 536             Best_solution = best_solution;
 537             Infinite = infinite;
 538             Epsilon = epsilon;
 539             Epsb = epsb;
 540             Epsd = epsd;
 541             Epsel = epsel;
 542             TREJ = TREJ;
 543             TINV = TINV;
 544             Maximise = maximise;
 545             Floor_first = floor_first;
 546             active = TRUE;
 547         }
 548
 549         private void ftran(int start, int end, float[] pcol) {
 550             int k, r;
 551             float theta;
 552             for(int i = start; i <= end; i++) {
 553                 k = Eta_col_end[i] - 1;
 554                 r = Eta_row_nr[k];
 555                 theta = pcol[r];
 556                 if (theta != 0) for(int j = Eta_col_end[i - 1]; j < k; j++)
 557                     pcol[Eta_row_nr[j]] += theta * Eta_value[j];
 558                 pcol[r] *= Eta_value[k];
 559             }
 560             for(int i = 0; i <= Rows; i++) round(pcol[i], Epsel);
 561         }
 562
 563         private void btran(float[] row) {
 564             int k;
 565             float f;
 566             for(int i = Eta_size; i >= 1; i--) {
 567                 f = 0;
 568                 k = Eta_col_end[i] - 1;
 569                 for(int j = Eta_col_end[i - 1]; j <= k; j++) f += row[Eta_row_nr[j]] * Eta_value[j];
 570                 f = round(f, Epsel);
 571                 row[Eta_row_nr[k]] = f;
 572             }
 573         }
 574
 575         static int[] num = new int[65535];
 576         static int[] rownum = new int[65535];
 577         static int[] colnum = new int[65535];
 578
 579         short Isvalid() {
 580             int row_nr;
 581             if (row_end_valid == FALSE) {
 582                 for(int i = 0; i <= rows; i++) { num[i] = 0; rownum[i] = 0; }
 583                 for(int i = 0; i < non_zeros; i++) rownum[get_row_nr(mat, i)]++;
 584                 row_end[0] = 0;
 585                 for(int i = 1; i <= rows; i++) row_end[i] = row_end[i - 1] + rownum[i];
 586                 for(int i = 1; i <= columns; i++)
 587                     for(int j = col_end[i - 1]; j < col_end[i]; j++) {
 588                         row_nr = get_row_nr(mat, j);
 589                         if (row_nr != 0) {
 590                             num[row_nr]++;
 591                             col_no[row_end[row_nr - 1] + num[row_nr]] = i;
 592                         }
 593                     }
 594                 row_end_valid = TRUE;
 595             }
 596             if (valid != FALSE) return(TRUE);
 597             for(int i = 0; i <= rows; i++) rownum[i] = 0;
 598             for(int i = 0; i <= columns; i++) colnum[i] = 0;
 599             for(int i = 1 ; i <= columns; i++)
 600                 for(int j = col_end[i - 1]; j < col_end[i]; j++) {
 601                     colnum[i]++;
 602                     rownum[get_row_nr(mat, j)]++;
 603                 }
 604             for(int i = 1; i <= columns; i++)
 605                 if (colnum[i] == 0)
 606                     throw new Error("Warning: Variable " + i + " not used in any constaints\n");
 607             valid = TRUE;
 608             return(TRUE);
 609         }
 610
 611         private void resize_eta() {
 612             Eta_alloc *= 2;
 613             throw new Error("eta undersized; this should never happen");
 614             /*
 615             float[] db_ptr = Eta_value;
 616             Eta_value = new float[Eta_alloc];
 617             System.arraycopy(db_ptr, 0, Eta_value, 0, db_ptr.length);
 618             eta_value = Eta_value;
 619
 620             int[] int_ptr = Eta_row_nr;
 621             Eta_row_nr = new int[Eta_alloc];
 622             System.arraycopy(int_ptr, 0, Eta_row_nr, 0, int_ptr.length);
 623             eta_row_nr = Eta_row_nr;
 624             */
 625         }
 626
 627         private void condensecol(int row_nr, float[] pcol) {
 628             int elnr;
 629             elnr = Eta_col_end[Eta_size];
 630             if (elnr + Rows + 2 > Eta_alloc) resize_eta();
 631             for(int i = 0; i <= Rows; i++)
 632                 if (i != row_nr && pcol[i] != 0) {
 633                     Eta_row_nr[elnr] = i;
 634                     Eta_value[elnr] = pcol[i];
 635                     elnr++;
 636                 }
 637             Eta_row_nr[elnr] = row_nr;
 638             Eta_value[elnr] = pcol[row_nr];
 639             elnr++;
 640             Eta_col_end[Eta_size + 1] = elnr;
 641         }
 642
 643         private void addetacol() {
 644             int k;
 645             float theta;
 646             int j = Eta_col_end[Eta_size];
 647             Eta_size++;
 648             k = Eta_col_end[Eta_size];
 649             theta = 1 / (float) Eta_value[k - 1];
 650             Eta_value[k - 1] = theta;
 651             for(int i = j; i < k - 1; i++) Eta_value[i] *= -theta;
 652             JustInverted = FALSE;
 653         }
 654
 655         private void setpivcol(short lower,  int varin, float[]   pcol) {
 656             int colnr;
 657             float f;
 658             if (lower != FALSE) f = 1;
 659             else f = -1;
 660             for(int i = 0; i <= Rows; i++) pcol[i] = 0;
 661             if (varin > Rows) {
 662                 colnr = varin - Rows;
 663                 for(int i = Col_end[colnr - 1]; i < Col_end[colnr]; i++) pcol[get_row_nr(Mat, i)] = get_value(Mat,i) * f;
 664                 pcol[0] -= Extrad * f;
 665             } else {
 666                 if (lower != FALSE) pcol[varin] = 1;
 667                 else pcol[varin] = -1;
 668             }
 669             ftran(1, Eta_size, pcol);
 670         }
 671
 672         private void minoriteration(int colnr, int row_nr) {
 673             int k, wk, varin, varout, elnr;
 674             float piv = 0, theta;
 675             varin = colnr + Rows;
 676             elnr = Eta_col_end[Eta_size];
 677             wk = elnr;
 678             Eta_size++;
 679             if (Extrad != 0) {
 680                 Eta_row_nr[elnr] = 0;
 681                 Eta_value[elnr] = -Extrad;
 682                 elnr++;
 683             }
 684             for(int j = Col_end[colnr - 1] ; j < Col_end[colnr]; j++) {
 685                 k = get_row_nr(Mat, j);
 686                 if (k == 0 && Extrad != 0) Eta_value[Eta_col_end[Eta_size -1]] += get_value(Mat,j);
 687                 else if (k != row_nr) {
 688                     Eta_row_nr[elnr] = k;
 689                     Eta_value[elnr] = get_value(Mat,j);
 690                     elnr++;
 691                 } else {
 692                     piv = get_value(Mat,j);
 693                 }
 694             }
 695             Eta_row_nr[elnr] = row_nr;
 696             Eta_value[elnr] = 1 / (float) piv;
 697             elnr++;
 698             theta = Rhs[row_nr] / (float) piv;
 699             Rhs[row_nr] = theta;
 700             for(int i = wk; i < elnr - 1; i++) Rhs[Eta_row_nr[i]] -= theta * Eta_value[i];
 701             varout = Bas[row_nr];
 702             Bas[row_nr] = varin;
 703             Basis[varout] = FALSE;
 704             Basis[varin] = TRUE;
 705             for(int i = wk; i < elnr - 1; i++) Eta_value[i] /= - (float) piv;
 706             Eta_col_end[Eta_size] = elnr;
 707         }
 708
 709         private void rhsmincol(float theta, int row_nr, int varin) {
 710             int varout;
 711             float f;
 712             if (row_nr > Rows + 1) {
 713                 System.err.println("Error: rhsmincol called with row_nr: " + row_nr + ", rows: " + Rows + "\n");
 714                 System.err.println("This indicates numerical instability\n");
 715             }
 716             int j = Eta_col_end[Eta_size];
 717             int k = Eta_col_end[Eta_size + 1];
 718             for(int i = j; i < k; i++) {
 719                 f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i];
 720                 f = round(f, Epsb);
 721                 Rhs[Eta_row_nr[i]] = f;
 722             }
 723             Rhs[row_nr] = theta;
 724             varout = Bas[row_nr];
 725             Bas[row_nr] = varin;
 726             Basis[varout] = FALSE;
 727             Basis[varin] = TRUE;
 728         }
 729
 730         private static int[] rownum_ = new int[65535];
 731         private static int[] colnum_ = new int[65535];
 732         private static int[] col = new int[65535];
 733         private static int[] row = new int[65535];
 734         private static float[] pcol = new float[65535];
 735         private static short[] frow = new short[65535];
 736         private static short[] fcol = new short[65535];
 737
 738         void invert() {
 739             int    v, wk, numit, varnr, row_nr, colnr, varin;
 740             float    theta;
 741
 742             for(int i = 0; i <= Rows; i++) rownum_[i] = 0;
 743             for(int i = 0; i <= Rows; i++) col[i] = 0;
 744             for(int i = 0; i <= Rows; i++) row[i] = 0;
 745             for(int i = 0; i <= Rows; i++) pcol[i] = 0;
 746             for(int i = 0; i <= Rows; i++) frow[i] = TRUE;
 747             for(int i = 0; i < columns; i++) fcol[i] = FALSE;
 748             for(int i = 0; i <= columns; i++) colnum_[i] = 0;
 749
 750             for(int i = 0; i <= Rows; i++)
 751                 if (Bas[i] > Rows) fcol[Bas[i] - Rows - 1] = TRUE;
 752                 else frow[Bas[i]] = FALSE;
 753
 754             for(int i = 1; i <= Rows; i++)
 755                 if (frow[i] != FALSE)
 756                     for(int j = Row_end[i - 1] + 1; j <= Row_end[i]; j++) {
 757                         wk = Col_no[j];
 758                         if (fcol[wk - 1] != FALSE) {
 759                             colnum_[wk]++;
 760                             rownum_[i - 1]++;
 761                         }
 762                     }
 763
 764             for(int i = 1; i <= Rows; i++) Bas[i] = i;
 765             for(int i = 1; i <= Rows; i++) Basis[i] = TRUE;
 766             for(int i = 1; i <= columns; i++) Basis[i + Rows] = FALSE;
 767             for(int i = 0; i <= Rows; i++) Rhs[i] = Rh[i];
 768             for(int i = 1; i <= columns; i++) {
 769                 varnr = Rows + i;
 770                 if (Lower[varnr] == FALSE) {
 771                     theta = Upbo[varnr];
 772                     for(int j = Col_end[i - 1]; j < Col_end[i]; j++)
 773                         Rhs[get_row_nr(Mat, j)] -= theta * get_value(Mat,j);
 774                 }
 775             }
 776             for(int i = 1; i <= Rows; i++) if (Lower[i] == FALSE) Rhs[i] -= Upbo[i];
 777             Eta_size = 0;
 778             v = 0;
 779             row_nr = 0;
 780             Num_inv = 0;
 781             numit = 0;
 782             while(v < Rows) {
 783                 int j;
 784                 row_nr++;
 785                 if (row_nr > Rows) row_nr = 1;
 786                 v++;
 787                 if (rownum_[row_nr - 1] == 1)
 788                     if (frow[row_nr] != FALSE) {
 789                         v = 0;
 790                         j = Row_end[row_nr - 1] + 1;
 791                         while(fcol[Col_no[j] - 1] == FALSE) j++;
 792                         colnr = Col_no[j];
 793                         fcol[colnr - 1] = FALSE;
 794                         colnum_[colnr] = 0;
 795                         for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++)
 796                             if (frow[get_row_nr(Mat, j)] != FALSE)
 797                                 rownum_[get_row_nr(Mat, j) - 1]--;
 798                         frow[row_nr] = FALSE;
 799                         minoriteration(colnr, row_nr);
 800                     }
 801             }
 802             v = 0;
 803             colnr = 0;
 804             while(v < columns) {
 805                 int j;
 806                 colnr++;
 807                 if (colnr > columns) colnr = 1;
 808                 v++;
 809                 if (colnum_[colnr] == 1)
 810                     if (fcol[colnr - 1] != FALSE) {
 811                         v = 0;
 812                         j = Col_end[colnr - 1] + 1;
 813                         while(frow[get_row_nr(Mat, j - 1)] == FALSE) j++;
 814                         row_nr = get_row_nr(Mat, j - 1);
 815                         frow[row_nr] = FALSE;
 816                         rownum_[row_nr - 1] = 0;
 817                         for(j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++)
 818                             if (fcol[Col_no[j] - 1] != FALSE)
 819                                 colnum_[Col_no[j]]--;
 820                         fcol[colnr - 1] = FALSE;
 821                         numit++;
 822                         col[numit - 1] = colnr;
 823                         row[numit - 1] = row_nr;
 824                     }
 825             }
 826             for(int j = 1; j <= columns; j++)
 827                 if (fcol[j - 1] != FALSE) {
 828                     fcol[j - 1] = FALSE;
 829                     setpivcol(Lower[Rows + j], j + Rows, pcol);
 830                     row_nr = 1;
 831                     while((frow[row_nr] == FALSE || pcol[row_nr] == FALSE) && row_nr <= Rows)
 832                         row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes
 833                                      rhsmincol crash. Solved in 2.0? MB */
 834                     if (row_nr == Rows + 1) throw new Error("Inverting failed");
 835                     frow[row_nr] = FALSE;
 836                     condensecol(row_nr, pcol);
 837                     theta = Rhs[row_nr] / (float) pcol[row_nr];
 838                     rhsmincol(theta, row_nr, Rows + j);
 839                     addetacol();
 840                 }
 841             for(int i = numit - 1; i >= 0; i--) {
 842                 colnr = col[i];
 843                 row_nr = row[i];
 844                 varin = colnr + Rows;
 845                 for(int j = 0; j <= Rows; j++) pcol[j] = 0;
 846                 for(int j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) pcol[get_row_nr(Mat, j)] = get_value(Mat,j);
 847                 pcol[0] -= Extrad;
 848                 condensecol(row_nr, pcol);
 849                 theta = Rhs[row_nr] / (float) pcol[row_nr];
 850                 rhsmincol(theta, row_nr, varin);
 851                 addetacol();
 852             }
 853             for(int i = 1; i <= Rows; i++) Rhs[i] = round(Rhs[i], Epsb);
 854             JustInverted = TRUE;
 855             DoInvert = FALSE;
 856         }
 857
 858         private short colprim(Ref colnr, short minit, float[]   drow) {
 859             int  varnr;
 860             float f, dpiv;
 861               dpiv = -Epsd;
 862             colnr.value = 0;
 863             if (minit == FALSE) {
 864                 for(int i = 1; i <= Sum; i++) drow[i] = 0;
 865                 drow[0] = 1;
 866                 btran(drow);
 867                 for(int i = 1; i <= columns; i++) {
 868                     varnr = Rows + i;
 869                     if (Basis[varnr] == FALSE)
 870                         if (Upbo[varnr] > 0) {
 871                             f = 0;
 872                             for(int j = Col_end[i - 1]; j < Col_end[i]; j++) f += drow[get_row_nr(Mat, j)] * get_value(Mat,j);
 873                             drow[varnr] = f;
 874                         }
 875                 }
 876                 for(int i = 1; i <= Sum; i++) drow[i] = round(drow[i], Epsd);
 877             }
 878             for(int i = 1; i <= Sum; i++)
 879                 if (Basis[i] == FALSE)
 880                     if (Upbo[i] > 0) {
 881                         if (Lower[i] != FALSE) f = drow[i];
 882                         else f = -drow[i];
 883                         if (f < dpiv) {
 884                             dpiv = f;
 885                             colnr.value = i;
 886                         }
 887                     }
 888             if (colnr.value == 0) {
 889                 Doiter   = FALSE;
 890                 DoInvert = FALSE;
 891                 Status   = OPTIMAL;
 892             }
 893             return(colnr.value > 0 ? (short)1 : (short)0);
 894         }
 895
 896         private short rowprim(int colnr, Ref row_nr, Ref theta, float[] pcol) {
 897             float f = 0, quot;
 898             row_nr.value = 0;
 899             theta.value = Infinite;
 900             for(int i = 1; i <= Rows; i++) {
 901                 f = pcol[i];
 902                 if (Math.abs(f) < TREJ) f = 0;
 903                 if (f != 0) {
 904                     quot = 2 * Infinite;
 905                     if (f > 0) quot = Rhs[i] / (float) f;
 906                     else if (Upbo[Bas[i]] < Infinite) quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
 907                     round(quot, Epsel);
 908                     if (quot < theta.value) {
 909                         theta.value = quot;
 910                         row_nr.value = i;
 911                     }
 912                 }
 913             }
 914             if (row_nr.value == 0)
 915                 for(int i = 1; i <= Rows; i++) {
 916                     f = pcol[i];
 917                     if (f != 0) {
 918                         quot = 2 * Infinite;
 919                         if (f > 0) quot = Rhs[i] / (float) f;
 920                         else if (Upbo[Bas[i]] < Infinite) quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
 921                         quot = round(quot, Epsel);
 922                         if (quot < theta.value) {
 923                             theta.value = quot;
 924                             row_nr.value = i;
 925                         }
 926                     }
 927                 }
 928
 929             if (theta.value < 0) throw new Error("Warning: Numerical instability, qout = " + theta.value);
 930             if (row_nr.value == 0) {
 931                 if (Upbo[colnr] == Infinite) {
 932                     Doiter   = FALSE;
 933                     DoInvert = FALSE;
 934                     Status   = UNBOUNDED;
 935                 } else {
 936                     int i = 1;
 937                     while(pcol[i] >= 0 && i <= Rows) i++;
 938                     if (i > Rows) {
 939                         Lower[colnr] = FALSE;
 940                         Rhs[0] += Upbo[colnr]*pcol[0];
 941                         Doiter = FALSE;
 942                         DoInvert = FALSE;
 943                     } else if (pcol[i]<0) {
 944                         row_nr.value = i;
 945                     }
 946                 }
 947             }
 948             if (row_nr.value > 0) Doiter = TRUE;
 949             return((row_nr.value > 0) ? (short)1 : (short)0);
 950         }
 951
 952         private short rowdual(Ref row_nr) {
 953             int   i;
 954             float  f, g, minrhs;
 955             short artifs;
 956             row_nr.value = 0;
 957             minrhs = -Epsb;
 958             i = 0;
 959             artifs = FALSE;
 960             while(i < Rows && artifs == FALSE) {
 961                 i++;
 962                 f = Upbo[Bas[i]];
 963                 if (f == 0 && (Rhs[i] != 0)) {
 964                     artifs = TRUE;
 965                     row_nr.value = i;
 966                 } else {
 967                     if (Rhs[i] < f - Rhs[i]) g = Rhs[i];
 968                     else g = f - Rhs[i];
 969                     if (g < minrhs) {
 970                         minrhs = g;
 971                         row_nr.value = i;
 972                     }
 973                 }
 974             }
 975             return(row_nr.value > 0 ? (short)1 : (short)0);
 976         }
 977
 978         private short coldual(int row_nr, Ref colnr, short minit, float[] prow, float[] drow) {
 979             int r, varnr;
 980             float theta, quot, pivot, d, f, g;
 981             Doiter = FALSE;
 982             if (minit == FALSE) {
 983                 for(int i = 0; i <= Rows; i++) {
 984                     prow[i] = 0;
 985                     drow[i] = 0;
 986                 }
 987                 drow[0] = 1;
 988                 prow[row_nr] = 1;
 989                 for(int i = Eta_size; i >= 1; i--) {
 990                     d = 0;
 991                     f = 0;
 992                     r = Eta_row_nr[Eta_col_end[i] - 1];
 993                     for(int j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++) {
 994                         /* this is where the program consumes most cpu time */
 995                         f += prow[Eta_row_nr[j]] * Eta_value[j];
 996                         d += drow[Eta_row_nr[j]] * Eta_value[j];
 997                     }
 998                     f = round(f, Epsel);
 999                     prow[r] = f;
1000                     d = round(d, Epsel);
1001                     drow[r] = d;
1002                 }
1003                 for(int i = 1; i <= columns; i++) {
1004                     varnr = Rows + i;
1005                     if (Basis[varnr] == FALSE) {
1006                         d = - Extrad * drow[0];
1007                         f = 0;
1008                         for(int j = Col_end[i - 1]; j < Col_end[i]; j++) {
1009                             d = d + drow[get_row_nr(Mat, j)] * get_value(Mat,j);
1010                             f = f + prow[get_row_nr(Mat, j)] * get_value(Mat,j);
1011                         }
1012                         drow[varnr] = d;
1013                         prow[varnr] = f;
1014                     }
1015                 }
1016                 for(int i = 0; i <= Sum; i++) {
1017                     prow[i] = round(prow[i], Epsel);
1018                     drow[i] = round(drow[i], Epsd);
1019                 }
1020             }
1021             if (Rhs[row_nr] > Upbo[Bas[row_nr]]) g = -1;
1022             else g = 1;
1023             pivot = 0;
1024             colnr.value = 0;
1025             theta = Infinite;
1026             for(int i = 1; i <= Sum; i++) {
1027                 if (Lower[i] != FALSE) d = prow[i] * g;
1028                 else d = -prow[i] * g;
1029                 if ((d < 0) && (Basis[i] == FALSE) && (Upbo[i] > 0)) {
1030                     if (Lower[i] == FALSE) quot = -drow[i] / (float) d;
1031                     else quot = drow[i] / (float) d;
1032                     if (quot < theta) {
1033                         theta = quot;
1034                         pivot = d;
1035                         colnr.value = i;
1036                     } else if ((quot == theta) && (Math.abs(d) > Math.abs(pivot))) {
1037                         pivot = d;
1038                         colnr.value = i;
1039                     }
1040                 }
1041             }
1042             if (colnr.value > 0) Doiter = TRUE;
1043             return(colnr.value > 0 ? (short)1 : (short)0);
1044         }
1045
1046         private void iteration(int row_nr, int varin, Ref theta, float up, Ref minit, Ref low, short primal,float[] pcol) {
1047             int k, varout;
1048             float f;
1049             float pivot;
1050             iter++;
1051             minit.value = theta.value > (up + Epsb) ? 1 : 0;
1052             if (minit.value != 0) {
1053                 theta.value = up;
1054                 low.value = low.value == 0 ? 1 : 0;
1055             }
1056             k = Eta_col_end[Eta_size + 1];
1057             pivot = Eta_value[k - 1];
1058             for(int i = Eta_col_end[Eta_size]; i < k; i++) {
1059                 f = Rhs[Eta_row_nr[i]] - theta.value * Eta_value[i];
1060                 f = round(f, Epsb);
1061                 Rhs[Eta_row_nr[i]] = f;
1062             }
1063             if (minit.value == 0) {
1064                 Rhs[row_nr] = theta.value;
1065                 varout = Bas[row_nr];
1066                 Bas[row_nr] = varin;
1067                 Basis[varout] = FALSE;
1068                 Basis[varin] = TRUE;
1069                 if (primal != FALSE && pivot < 0) Lower[varout] = FALSE;
1070                 if (low.value == 0 && up < Infinite) {
1071                     low.value = TRUE;
1072                     Rhs[row_nr] = up - Rhs[row_nr];
1073                     for(int i = Eta_col_end[Eta_size]; i < k; i++) Eta_value[i] = -Eta_value[i];
1074                 }
1075                 addetacol();
1076                 Num_inv++;
1077             }
1078         }
1079
1080         static float[] drow = new float[65535];
1081         static float[] prow = new float[65535];
1082         static float[] Pcol = new float[65535];
1083
1084         private int solvelp() {
1085             int    varnr;
1086             float   f = 0, theta = 0;
1087             short  primal;
1088             short  minit;
1089             int    colnr, row_nr;
1090             colnr = 0;
1091             row_nr = 0;
1092             short flag;
1093             Ref ref1, ref2, ref3;
1094             ref1 = new Ref(0);
1095             ref2 = new Ref(0);
1096             ref3 = new Ref(0);
1097
1098             for(int i = 0; i <= Sum; i++) { drow[i] = 0; prow[i] = 0; }
1099             for(int i = 0; i <= Rows; i++) Pcol[i] = 0;
1100             iter = 0;
1101             minit = FALSE;
1102             Status = RUNNING;
1103             DoInvert = FALSE;
1104             Doiter = FALSE;
1105             primal = TRUE;
1106             for(int i = 0; i != Rows && primal != FALSE;) {
1107                 i++;
1108                 primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ? (short)1: (short)0;
1109             }
1110             if (primal == FALSE) {
1111                 drow[0] = 1;
1112                 for(int i = 1; i <= Rows; i++) drow[i] = 0;
1113                 Extrad = 0;
1114                 for(int i = 1; i <= columns; i++) {
1115                     varnr = Rows + i;
1116                     drow[varnr] = 0;
1117                     for(int j = Col_end[i - 1]; j < Col_end[i]; j++)
1118                         if (drow[get_row_nr(Mat, j)] != 0)
1119                             drow[varnr] += drow[get_row_nr(Mat, j)] * get_value(Mat,j);
1120                     if (drow[varnr] < Extrad) Extrad = drow[varnr];
1121                 }
1122             } else {
1123                 Extrad = 0;
1124             }
1125             minit = FALSE;
1126             while(Status == RUNNING) {
1127                 Doiter = FALSE;
1128                 DoInvert = FALSE;
1129                 construct_solution(Solution);
1130                 if (primal != FALSE) {
1131                     ref1.value = colnr;
1132                     flag = colprim(ref1, minit, drow);
1133                     colnr = (int)ref1.value;
1134                     if (flag != FALSE) {
1135                         setpivcol(Lower[colnr], colnr, Pcol);
1136                         ref1.value = row_nr;
1137                         ref2.value = theta;
1138                         flag = rowprim(colnr, ref1, ref2, Pcol);
1139                         row_nr = (int)ref1.value;
1140                         theta = ref2.value;
1141                         if (flag != FALSE) condensecol(row_nr, Pcol);
1142                     }
1143                 } else {
1144                     if (minit == FALSE) {
1145                         ref1.value = row_nr;
1146                         flag = rowdual(ref1);
1147                         row_nr = (int)ref1.value;
1148                     }
1149                     if (row_nr > 0) {
1150                         ref1.value = colnr;
1151                         flag = coldual(row_nr, ref1, minit, prow, drow);
1152                         colnr = (int)ref1.value;
1153                         if (flag != FALSE) {
1154                             setpivcol(Lower[colnr], colnr, Pcol);
1155                             /* getting div by zero here ... MB */
1156                             if (Pcol[row_nr] == 0) {
1157                                 throw new Error("An attempt was made to divide by zero (Pcol[" + row_nr + "])");
1158                             } else {
1159                                 condensecol(row_nr, Pcol);
1160                                 f = Rhs[row_nr] - Upbo[Bas[row_nr]];
1161                                 if (f > 0) {
1162                                     theta = f / (float) Pcol[row_nr];
1163                                     if (theta <= Upbo[colnr])
1164                                         Lower[Bas[row_nr]] = (Lower[Bas[row_nr]] == FALSE)? (short)1:(short)0;
1165                                 } else theta = Rhs[row_nr] / (float) Pcol[row_nr];
1166                             }
1167                         } else Status = INFEASIBLE;
1168                     } else {
1169                         primal   = TRUE;
1170                         Doiter   = FALSE;
1171                         Extrad   = 0;
1172                         DoInvert = TRUE;
1173                     }
1174                 }
1175                 if (Doiter != FALSE) {
1176                     ref1.value = theta;
1177                     ref2.value = minit;
1178                     ref3.value = Lower[colnr];
1179                     iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol);
1180                     theta = ref1.value;
1181                     minit = (short)ref2.value;
1182                     Lower[colnr] = (short)ref3.value;
1183                 }
1184                 if (Num_inv >= max_num_inv) DoInvert = TRUE;
1185                 if (DoInvert != FALSE) invert();
1186             }
1187             total_iter += iter;
1188             return(Status);
1189         }
1190
1191         private void construct_solution(float[]   sol) {
1192             float   f;
1193             int basi;
1194             for(int i = 0; i <= Rows; i++) sol[i] = 0;
1195             for(int i = Rows + 1; i <= Sum; i++) sol[i] = Lowbo[i];
1196             for(int i = 1; i <= Rows; i++) {
1197                 basi = Bas[i];
1198                 if (basi > Rows) sol[basi] += Rhs[i];
1199             }
1200             for(int i = Rows + 1; i <= Sum; i++)
1201                 if (Basis[i] == FALSE && Lower[i] == FALSE)
1202                     sol[i] += Upbo[i];
1203             for(int j = 1; j <= columns; j++) {
1204                 f = sol[Rows + j];
1205                 if (f != 0)
1206                     for(int i = Col_end[j - 1]; i < Col_end[j]; i++)
1207                         sol[get_row_nr(Mat, i)] += f * get_value(Mat,i);
1208             }
1209             for(int i = 0; i <= Rows; i++) {
1210                 if (Math.abs(sol[i]) < Epsb) sol[i] = 0;
1211                 else if (ch_sign[i] != FALSE) sol[i] = -sol[i];
1212             }
1213         }
1214
1215         private void calculate_duals() {
1216             for(int i = 1; i <= Rows; i++) duals[i] = 0;
1217             duals[0] = 1;
1218             btran(duals);
1219             for(int i = 1; i <= Rows; i++) {
1220                 if (basis[i] != FALSE) duals[i] = 0;
1221                 else if ( ch_sign[0] == ch_sign[i]) duals[i] = -duals[i];
1222             }
1223         }
1224
1225         private static Random rdm = new Random();
1226
1227         private int milpsolve(float[]   upbo, float[]   lowbo, short[]  sbasis, short[]  slower, int[]    sbas) {
1228             int failure, notint, is_worse;
1229             float theta, tmpfloat;
1230             notint = 0;
1231
1232             if (Break_bb != FALSE) return(BREAK_BB);
1233             Level++;
1234             total_nodes++;
1235             if (Level > max_level) max_level = Level;
1236             System.arraycopy(upbo, 0, Upbo, 0, Sum + 1);
1237             System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1);
1238             System.arraycopy(sbasis, 0, Basis, 0, Sum + 1);
1239             System.arraycopy(slower, 0, Lower, 0, Sum + 1);
1240             System.arraycopy(sbas, 0, Bas, 0, Rows + 1);
1241             System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1);
1242             if (eta_valid == FALSE) {
1243                 for(int i = 1; i <= columns; i++)
1244                     if (Lowbo[Rows + i] != 0) {
1245                         theta = Lowbo[ Rows + i];
1246                         if (Upbo[Rows + i]<Infinite) Upbo[Rows + i] -= theta;
1247                         for(int j = Col_end[i - 1]; j < Col_end[i]; j++) Rh[get_row_nr(Mat, j)] -= theta * get_value(Mat,j);
1248                     }
1249                 invert();
1250                 eta_valid = TRUE;
1251             }
1252             failure = solvelp();
1253             if (failure == OPTIMAL) {
1254                 construct_solution(Solution);
1255                 /* if this solution is worse than the best sofar, this branch must die */
1256                 if (Maximise != FALSE) is_worse = (Solution[0] <= Best_solution[0]) ? 1:0;
1257                 else is_worse = (Solution[0] >= Best_solution[0]) ? 1:0;
1258                 if (is_worse != FALSE) {
1259                     Level--;
1260                     return(MILP_FAIL);
1261                 }
1262                 /* check if solution contains enough ints */
1263                 if (bb_rule == FIRST_NI) {
1264                     notint = 0;
1265                     int i = Rows + 1;
1266                     while(i <= Sum && notint == 0) i++;
1267                 }
1268                 if (bb_rule == RAND_NI) {
1269                     int nr_not_int, select_not_int;
1270                     nr_not_int = 0;
1271                     for(int i = Rows + 1; i <= Sum; i++)
1272                         if (nr_not_int == 0) notint = 0;
1273                         else {
1274                             select_not_int=(rdm.nextInt() % nr_not_int) + 1;
1275                             i = Rows + 1;
1276                             while(select_not_int > 0) i++;
1277                             notint = i - 1;
1278                         }
1279                 }
1280                 if (notint != FALSE) throw new Error("integer linear programming not supported");
1281                 if (Maximise != FALSE) is_worse = (Solution[0] < Best_solution[0]) ? 1:0;
1282                 else is_worse = (Solution[0] > Best_solution[0]) ? 1:0;
1283                 if (is_worse == FALSE) {
1284                     System.arraycopy(Solution, 0, Best_solution, 0, Sum + 1);
1285                     calculate_duals();
1286                     if (break_at_int != FALSE) {
1287                         if (Maximise != FALSE &&  (Best_solution[0] > break_value)) Break_bb = TRUE;
1288                         if (Maximise == FALSE &&  (Best_solution[0] < break_value)) Break_bb = TRUE;
1289                     }
1290                 }
1291             }
1292             Level--;
1293             return(failure);
1294         }
1295
1296         public int solve() {
1297             int result;
1298             if (active == FALSE) set_globals();
1299             total_iter  = 0;
1300             max_level   = 1;
1301             total_nodes = 0;
1302             if (Isvalid() != FALSE) {
1303                 if (Maximise != FALSE && obj_bound == Infinite) Best_solution[0]=-Infinite;
1304                 else if (Maximise == FALSE && obj_bound==-Infinite) Best_solution[0] = Infinite;
1305                 else Best_solution[0] = obj_bound;
1306                 Level = 0;
1307                 if (basis_valid == FALSE) {
1308                     for(int i = 0; i <= rows; i++) {
1309                         basis[i] = TRUE;
1310                         bas[i] = i;
1311                     }
1312                     for(int i = rows+1; i <= sum; i++) basis[i] = FALSE;
1313                     for(int i = 0; i <= sum; i++) lower[i] = TRUE;
1314                     basis_valid = TRUE;
1315                 }
1316                 eta_valid = FALSE;
1317                 Break_bb      = FALSE;
1318                 result        = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas);
1319                 eta_size  = Eta_size;
1320                 eta_alloc = Eta_alloc;
1321                 num_inv   = Num_inv;
1322                 return(result);
1323             }
1324             return(FAILURE);
1325         }
1326     }
1327
1328     private final static float round( float val, float eps) { return (Math.abs(val) < eps) ? 0 : val; }
1329     static int get_row_nr(MatrixArray m, int i) { return m.row_nr[i]; }
1330     static void set_row_nr(MatrixArray m, int i, int val) { m.row_nr[i] = val; }
1331     static float get_value(MatrixArray m, int i) { return m.value[i]; }
1332     static void set_value(MatrixArray m, int i, float val) { m.value[i] = val; }
1333     public static class MatrixArray {
1334         public int[] row_nr;
1335         public float[] value;
1336         public final int length;
1337         public MatrixArray(int length) { row_nr = new int[length]; value = new float[length]; this.length = length; }
1338     }
1339
1340 }
1341