7 public final static short FAIL = -1;
9 public final static short NULL = 0;
10 public final static short FALSE = 0;
11 public final static short TRUE = 1;
13 public final static short DEFNUMINV = 50;
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;
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;
28 public final static short FIRST_NI = 0;
29 public final static short RAND_NI = 1;
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;
36 public final static short MAX_WARN_COUNT = 20;
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 */
45 public final static float PREJ = (float)1e-3; /* pivot reject (try others first) */
47 public final static int ETA_START_SIZE = 10000; /* start size of array Eta. Realloced if needed */
51 public Ref(float v) { value = v; }
54 //public static class Simplex {
55 /* Globals used by solver */
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 */
68 int sum; /* The size of the variables + the slacks */
70 int non_zeros; /* The number of elements in the sparce matrix*/
71 int mat_alloc; /* The allocated size for matrix sized
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
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
89 float[] rhs; /* rows_alloc+1 :The RHS of the curent simplex
91 float[] orig_upbo; /* sum_alloc+1 :Bound before transformations */
92 float[] orig_lowbo; /* " " */
93 float[] upbo; /* " " :Upper bound after transformation
95 float[] lowbo; /* " " :Lower bound after transformation
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
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 */
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
112 float[] eta_value; /* eta_alloc :The Structure containing the
114 int[] eta_row_nr; /* " " :The Structure containing the Row
116 int[] eta_col_end; /* rows_alloc + MaxNumInv : eta_col_end[i] is
117 the start index of the next Eta column */
119 short bb_rule; /* what rule for selecting B&B variables */
121 short break_at_int; /* TRUE if stop at first integer better than
125 float obj_bound; /* ## Objective function bound for speedup of
127 int iter; /* The number of iterations in the simplex
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,
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
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
145 (a`x > b, x>=0) is translated to
146 s + -a`x = -b with x>=0, s>=0) */
148 int nr_lagrange; /* Nr. of Langrangian relaxation constraints */
149 float[][]lag_row; /* NumLagrange, columns+1:Pointer to pointer of
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 */
156 short valid; /* Has this lp pased the 'test' */
157 float infinite; /* ## numercal stuff */
158 float epsilon; /* ## */
161 float epsel; /* ## */
189 public float[] Best_solution;
203 int Warn_count; /* used in CHECK version of rounding macro */
205 public Simplex (int nrows, int ncolumns, int matalloc) {
212 columns_alloc=columns;
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);
241 public void init(int ncolumns) {
251 obj_bound=DEF_INFINITE;
252 infinite=DEF_INFINITE;
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;
293 public void setObjective(float[] row, boolean maximize) {
294 for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
296 for(int j = 1; j <= columns; j++) {
299 float Value = row[j];
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");
305 if (basis[column] == TRUE && Row > 0) basis_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);
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 */
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));
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);
328 if (active != FALSE) Non_zeros=non_zeros;
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);
340 if (active != FALSE) Maximise=TRUE;
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);
350 if (active != FALSE) Maximise=FALSE;
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];
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");
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;
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));
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);
393 for(int i = sum ; i > rows; i--) {
394 orig_upbo[i]=orig_upbo[i-1];
395 orig_lowbo[i]=orig_lowbo[i-1];
400 for(int i = 1 ; i <= rows; i++) if(bas[i] >= rows) bas[i]++;
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");
407 if(constr_type==GE && rh != 0) orig_rh[rows]=-rh;
408 else orig_rh[rows]=rh;
416 if (active != FALSE) set_globals();
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;
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;
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");
440 orig_upbo[rows+column] = value;
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");
447 orig_lowbo[rows+column] = value;
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;
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];
466 public void set_constr_type(int row, short con_type) {
467 if (row > rows || row < 1) throw new Error("Row out of Range");
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);
477 if (orig_rh[row]!=0) orig_rh[row]*=-1;
481 orig_upbo[row]=infinite;
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);
488 if (orig_rh[row]!=0) orig_rh[row]*=-1;
492 orig_upbo[row]=infinite;
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);
499 if (orig_rh[row]!=0) orig_rh[row]*=-1;
502 default: throw new Error("Constraint type not (yet) implemented");
509 Sum = Rows + columns;
510 Non_zeros = non_zeros;
518 Orig_upbo = orig_upbo;
519 Orig_lowbo = orig_lowbo;
525 Eta_alloc = eta_alloc;
528 Eta_value = eta_value;
529 Eta_row_nr = eta_row_nr;
530 Eta_col_end = eta_col_end;
532 Best_solution = best_solution;
541 Floor_first = floor_first;
545 private void ftran(int start, int end, float[] pcol) {
548 for(int i = start; i <= end; i++) {
549 k = Eta_col_end[i] - 1;
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];
556 for(int i = 0; i <= Rows; i++) round(pcol[i], Epsel);
559 private void btran(float[] row) {
562 for(int i = Eta_size; i >= 1; i--) {
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];
567 row[Eta_row_nr[k]] = f;
571 static int[] num = new int[65535];
572 static int[] rownum = new int[65535];
573 static int[] colnum = new int[65535];
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)]++;
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);
587 col_no[row_end[row_nr - 1] + num[row_nr]] = i;
590 row_end_valid = TRUE;
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++) {
598 rownum[get_row_nr(mat, j)]++;
600 for(int i = 1; i <= columns; i++)
602 throw new Error("Warning: Variable " + i + " not used in any constaints\n");
607 private void resize_eta() {
609 throw new Error("eta undersized; this should never happen");
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;
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;
623 private void condensecol(int row_nr, float[] pcol) {
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];
633 Eta_row_nr[elnr] = row_nr;
634 Eta_value[elnr] = pcol[row_nr];
636 Eta_col_end[Eta_size + 1] = elnr;
639 private void addetacol() {
642 int j = Eta_col_end[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;
651 private void setpivcol(short lower, int varin, float[] pcol) {
654 if (lower != FALSE) f = 1;
656 for(int i = 0; i <= Rows; i++) pcol[i] = 0;
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;
662 if (lower != FALSE) pcol[varin] = 1;
663 else pcol[varin] = -1;
665 ftran(1, Eta_size, pcol);
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];
676 Eta_row_nr[elnr] = 0;
677 Eta_value[elnr] = -Extrad;
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);
688 piv = get_value(Mat,j);
691 Eta_row_nr[elnr] = row_nr;
692 Eta_value[elnr] = 1 / (float) piv;
694 theta = Rhs[row_nr] / (float) piv;
696 for(int i = wk; i < elnr - 1; i++) Rhs[Eta_row_nr[i]] -= theta * Eta_value[i];
697 varout = Bas[row_nr];
699 Basis[varout] = FALSE;
701 for(int i = wk; i < elnr - 1; i++) Eta_value[i] /= - (float) piv;
702 Eta_col_end[Eta_size] = elnr;
705 private void rhsmincol(float theta, int row_nr, int varin) {
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");
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];
717 Rhs[Eta_row_nr[i]] = f;
720 varout = Bas[row_nr];
722 Basis[varout] = FALSE;
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];
735 int v, wk, numit, varnr, row_nr, colnr, varin;
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;
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;
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++) {
754 if (fcol[wk - 1] != FALSE) {
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++) {
766 if (Lower[varnr] == FALSE) {
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);
772 for(int i = 1; i <= Rows; i++) if (Lower[i] == FALSE) Rhs[i] -= Upbo[i];
781 if (row_nr > Rows) row_nr = 1;
783 if (rownum_[row_nr - 1] == 1)
784 if (frow[row_nr] != FALSE) {
786 j = Row_end[row_nr - 1] + 1;
787 while(fcol[Col_no[j] - 1] == FALSE) j++;
789 fcol[colnr - 1] = FALSE;
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);
803 if (colnr > columns) colnr = 1;
805 if (colnum_[colnr] == 1)
806 if (fcol[colnr - 1] != FALSE) {
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;
818 col[numit - 1] = colnr;
819 row[numit - 1] = row_nr;
822 for(int j = 1; j <= columns; j++)
823 if (fcol[j - 1] != FALSE) {
825 setpivcol(Lower[Rows + j], j + Rows, pcol);
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);
837 for(int i = numit - 1; i >= 0; 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);
844 condensecol(row_nr, pcol);
845 theta = Rhs[row_nr] / (float) pcol[row_nr];
846 rhsmincol(theta, row_nr, varin);
849 for(int i = 1; i <= Rows; i++) Rhs[i] = round(Rhs[i], Epsb);
854 private short colprim(Ref colnr, short minit, float[] drow) {
859 if (minit == FALSE) {
860 for(int i = 1; i <= Sum; i++) drow[i] = 0;
863 for(int i = 1; i <= columns; i++) {
865 if (Basis[varnr] == FALSE)
866 if (Upbo[varnr] > 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);
872 for(int i = 1; i <= Sum; i++) drow[i] = round(drow[i], Epsd);
874 for(int i = 1; i <= Sum; i++)
875 if (Basis[i] == FALSE)
877 if (Lower[i] != FALSE) f = drow[i];
884 if (colnr.value == 0) {
889 return(colnr.value > 0 ? (short)1 : (short)0);
892 private short rowprim(int colnr, Ref row_nr, Ref theta, float[] pcol) {
895 theta.value = Infinite;
896 for(int i = 1; i <= Rows; i++) {
898 if (Math.abs(f) < TREJ) f = 0;
901 if (f > 0) quot = Rhs[i] / (float) f;
902 else if (Upbo[Bas[i]] < Infinite) quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
904 if (quot < theta.value) {
910 if (row_nr.value == 0)
911 for(int i = 1; i <= Rows; i++) {
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) {
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) {
933 while(pcol[i] >= 0 && i <= Rows) i++;
935 Lower[colnr] = FALSE;
936 Rhs[0] += Upbo[colnr]*pcol[0];
939 } else if (pcol[i]<0) {
944 if (row_nr.value > 0) Doiter = TRUE;
945 return((row_nr.value > 0) ? (short)1 : (short)0);
948 private short rowdual(Ref row_nr) {
956 while(i < Rows && artifs == FALSE) {
959 if (f == 0 && (Rhs[i] != 0)) {
963 if (Rhs[i] < f - Rhs[i]) g = Rhs[i];
971 return(row_nr.value > 0 ? (short)1 : (short)0);
974 private short coldual(int row_nr, Ref colnr, short minit, float[] prow, float[] drow) {
976 float theta, quot, pivot, d, f, g;
978 if (minit == FALSE) {
979 for(int i = 0; i <= Rows; i++) {
985 for(int i = Eta_size; i >= 1; i--) {
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];
999 for(int i = 1; i <= columns; i++) {
1001 if (Basis[varnr] == FALSE) {
1002 d = - Extrad * drow[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);
1012 for(int i = 0; i <= Sum; i++) {
1013 prow[i] = round(prow[i], Epsel);
1014 drow[i] = round(drow[i], Epsd);
1017 if (Rhs[row_nr] > Upbo[Bas[row_nr]]) g = -1;
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;
1032 } else if ((quot == theta) && (Math.abs(d) > Math.abs(pivot))) {
1038 if (colnr.value > 0) Doiter = TRUE;
1039 return(colnr.value > 0 ? (short)1 : (short)0);
1042 private void iteration(int row_nr, int varin, Ref theta, float up, Ref minit, Ref low, short primal,float[] pcol) {
1047 minit.value = theta.value > (up + Epsb) ? 1 : 0;
1048 if (minit.value != 0) {
1050 low.value = low.value == 0 ? 1 : 0;
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];
1057 Rhs[Eta_row_nr[i]] = f;
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) {
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];
1076 static float[] drow = new float[65535];
1077 static float[] prow = new float[65535];
1078 static float[] Pcol = new float[65535];
1080 private int solvelp() {
1082 float f = 0, theta = 0;
1089 Ref ref1, ref2, ref3;
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;
1102 for(int i = 0; i != Rows && primal != FALSE;) {
1104 primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ? (short)1: (short)0;
1106 if (primal == FALSE) {
1108 for(int i = 1; i <= Rows; i++) drow[i] = 0;
1110 for(int i = 1; i <= columns; i++) {
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];
1122 while(Status == RUNNING) {
1125 construct_solution(Solution);
1126 if (primal != FALSE) {
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;
1134 flag = rowprim(colnr, ref1, ref2, Pcol);
1135 row_nr = (int)ref1.value;
1137 if (flag != FALSE) condensecol(row_nr, Pcol);
1140 if (minit == FALSE) {
1141 ref1.value = row_nr;
1142 flag = rowdual(ref1);
1143 row_nr = (int)ref1.value;
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 + "])");
1155 condensecol(row_nr, Pcol);
1156 f = Rhs[row_nr] - Upbo[Bas[row_nr]];
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];
1163 } else Status = INFEASIBLE;
1171 if (Doiter != FALSE) {
1174 ref3.value = Lower[colnr];
1175 iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol);
1177 minit = (short)ref2.value;
1178 Lower[colnr] = (short)ref3.value;
1180 if (Num_inv >= max_num_inv) DoInvert = TRUE;
1181 if (DoInvert != FALSE) invert();
1187 private void construct_solution(float[] sol) {
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++) {
1194 if (basi > Rows) sol[basi] += Rhs[i];
1196 for(int i = Rows + 1; i <= Sum; i++)
1197 if (Basis[i] == FALSE && Lower[i] == FALSE)
1199 for(int j = 1; j <= columns; j++) {
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);
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];
1211 private void calculate_duals() {
1212 for(int i = 1; i <= Rows; i++) duals[i] = 0;
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];
1221 private static Random rdm = new Random();
1223 private int milpsolve(float[] upbo, float[] lowbo, short[] sbasis, short[] slower, int[] sbas) {
1224 int failure, notint, is_worse;
1225 float theta, tmpfloat;
1228 if (Break_bb != FALSE) return(BREAK_BB);
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);
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) {
1258 /* check if solution contains enough ints */
1259 if (bb_rule == FIRST_NI) {
1262 while(i <= Sum && notint == 0) i++;
1264 if (bb_rule == RAND_NI) {
1265 int nr_not_int, select_not_int;
1267 for(int i = Rows + 1; i <= Sum; i++)
1268 if (nr_not_int == 0) notint = 0;
1270 select_not_int=(rdm.nextInt() % nr_not_int) + 1;
1272 while(select_not_int > 0) i++;
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);
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;
1292 public int solve() {
1293 int result = FAILURE;
1294 if (active == FALSE) set_globals();
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;
1303 if (basis_valid == FALSE) {
1304 for(int i = 0; i <= rows; i++) {
1308 for(int i = rows+1; i <= sum; i++) basis[i] = FALSE;
1309 for(int i = 0; i <= sum; i++) lower[i] = TRUE;
1314 result = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas);
1315 eta_size = Eta_size;
1316 eta_alloc = Eta_alloc;
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;
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;
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; }