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