src/org/ibex/util/Simplex.java

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