5 public class LinearProgramming {
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 HASHSIZE = 10007; /* prime number is better, MB */
48 public final static int ETA_START_SIZE = 10000; /* start size of array Eta. Realloced if needed */
49 public final static String STD_ROW_NAME_PREFIX = "r_";
53 public Ref(float v) { value = v; }
56 public static class Simplex {
58 Problem Lp; /* pointer to active problem */
85 float[] Best_solution;
99 int Warn_count; /* used in CHECK version of rounding macro */
101 void inc_mat_space(Problem lp, int maxextra)
103 if(lp.non_zeros + maxextra > lp.mat_alloc)
106 lp.mat_alloc = lp.non_zeros + maxextra;
107 Matrix[] old_mat = lp.mat;
108 lp.mat = new Matrix[lp.mat_alloc];
109 for (i = old_mat.length; i < lp.mat.length; i++)
110 lp.mat[i] = new Matrix(0,0);
111 System.arraycopy(old_mat, 0, lp.mat, 0, old_mat.length);
112 int[] old_col_no = lp.col_no;
113 lp.col_no = new int[lp.mat_alloc];
114 System.arraycopy(old_col_no, 0, lp.col_no, 0, old_col_no.length);
115 if (lp.active != FALSE)
121 } // end of inc_mat_space
123 void inc_row_space(Problem lp)
125 if(lp.rows > lp.rows_alloc)
127 lp.rows_alloc=lp.rows+10;
128 lp.sum_alloc=lp.rows_alloc+lp.columns_alloc;
130 float[] db_ptr = lp.orig_rh;
131 lp.orig_rh = new float[lp.rows_alloc + 1];
132 System.arraycopy(db_ptr, 0, lp.orig_rh, 0, db_ptr.length);
135 lp.rh = new float[lp.rows_alloc + 1];
136 System.arraycopy(db_ptr, 0, lp.rh, 0, db_ptr.length);
139 lp.rhs = new float[lp.rows_alloc + 1];
140 System.arraycopy(db_ptr, 0, lp.rhs, 0, db_ptr.length);
142 db_ptr = lp.orig_upbo;
143 lp.orig_upbo = new float[lp.sum_alloc + 1];
144 System.arraycopy(db_ptr, 0, lp.orig_upbo, 0, db_ptr.length);
147 lp.upbo = new float[lp.sum_alloc + 1];
148 System.arraycopy(db_ptr, 0, lp.upbo, 0, db_ptr.length);
150 db_ptr = lp.orig_lowbo;
151 lp.orig_lowbo = new float[lp.sum_alloc + 1];
152 System.arraycopy(db_ptr, 0, lp.orig_lowbo, 0, db_ptr.length);
155 lp.lowbo = new float[lp.sum_alloc + 1];
156 System.arraycopy(db_ptr, 0, lp.lowbo, 0, db_ptr.length);
158 db_ptr = lp.solution;
159 lp.solution = new float[lp.sum_alloc + 1];
160 System.arraycopy(db_ptr, 0, lp.solution, 0, db_ptr.length);
162 db_ptr = lp.best_solution;
163 lp.best_solution = new float[lp.sum_alloc + 1];
164 System.arraycopy(db_ptr, 0, lp.best_solution, 0, db_ptr.length);
166 int[] int_ptr = lp.row_end;
167 lp.row_end = new int[lp.rows_alloc + 1];
168 System.arraycopy(int_ptr, 0, lp.row_end, 0, int_ptr.length);
170 short[] short_ptr = lp.basis;
171 lp.basis = new short[lp.sum_alloc + 1];
172 System.arraycopy(short_ptr, 0, lp.basis, 0, short_ptr.length);
174 short_ptr = lp.lower;
175 lp.lower = new short[lp.sum_alloc + 1];
176 System.arraycopy(short_ptr, 0, lp.lower, 0, short_ptr.length);
179 lp.bas = new int[lp.rows_alloc + 1];
180 System.arraycopy(int_ptr, 0, lp.bas, 0, int_ptr.length);
183 lp.duals = new float[lp.rows_alloc + 1];
184 System.arraycopy(db_ptr, 0, lp.duals, 0, db_ptr.length);
186 short_ptr = lp.ch_sign;
187 lp.ch_sign = new short[lp.rows_alloc + 1];
188 System.arraycopy(short_ptr, 0, lp.ch_sign, 0, short_ptr.length);
190 int_ptr = lp.eta_col_end;
191 lp.eta_col_end = new int[lp.rows_alloc + lp.max_num_inv];
192 System.arraycopy(int_ptr, 0, lp.eta_col_end, 0, int_ptr.length);
194 if(lp.names_used != FALSE) {
195 String[] str_ptr = lp.row_name;
196 lp.row_name = new String[lp.rows_alloc + 1];
197 System.arraycopy(str_ptr, 0, lp.row_name, 0, str_ptr.length);
200 if(lp.scaling_used != FALSE) {
202 lp.scale = new float[lp.sum_alloc + 1];
203 System.arraycopy(db_ptr, 0, lp.scale, 0, db_ptr.length);
206 if(lp.active != FALSE)
211 void inc_col_space(Problem lp)
213 if(lp.columns >= lp.columns_alloc)
216 lp.columns_alloc=lp.columns+10;
217 lp.sum_alloc=lp.rows_alloc+lp.columns_alloc;
219 float[] float_ptr = lp.orig_upbo;
220 lp.orig_upbo = new float[lp.sum_alloc + 1];
221 System.arraycopy(float_ptr, 0, lp.orig_upbo, 0, float_ptr.length);
224 lp.upbo = new float[lp.sum_alloc + 1];
225 System.arraycopy(float_ptr, 0, lp.upbo, 0, float_ptr.length);
227 float_ptr = lp.orig_lowbo;
228 lp.orig_lowbo = new float[lp.sum_alloc + 1];
229 System.arraycopy(float_ptr, 0, lp.orig_lowbo, 0, float_ptr.length);
231 float_ptr = lp.lowbo;
232 lp.lowbo = new float[lp.sum_alloc + 1];
233 System.arraycopy(float_ptr, 0, lp.lowbo, 0, float_ptr.length);
235 float_ptr = lp.solution;
236 lp.solution = new float[lp.sum_alloc + 1];
237 System.arraycopy(float_ptr, 0, lp.solution, 0, float_ptr.length);
239 float_ptr = lp.best_solution;
240 lp.best_solution = new float[lp.sum_alloc + 1];
241 System.arraycopy(float_ptr, 0, lp.best_solution, 0, float_ptr.length);
243 short[] short_ptr = lp.basis;
244 lp.basis = new short[lp.sum_alloc + 1];
245 System.arraycopy(short_ptr, 0, lp.basis, 0, short_ptr.length);
247 short_ptr = lp.lower;
248 lp.lower = new short[lp.sum_alloc + 1];
249 System.arraycopy(short_ptr, 0, lp.lower, 0, short_ptr.length);
251 if(lp.names_used != FALSE) {
252 int[] str_ptr = lp.col_name;
253 lp.col_name = new int[lp.columns_alloc + 1];
254 System.arraycopy(str_ptr, 0, lp.col_name, 0, str_ptr.length);
257 if(lp.scaling_used != FALSE) {
258 float_ptr = lp.scale;
259 lp.scale = new float[lp.sum_alloc + 1];
260 System.arraycopy(float_ptr, 0, lp.scale, 0, float_ptr.length);
263 int_ptr = lp.col_end;
264 lp.col_end = new int[lp.columns_alloc + 1];
265 System.arraycopy(int_ptr, 0, lp.col_end, 0, int_ptr.length);
267 if(lp.active != FALSE)
270 } // end of inc_col_space
272 public void set_mat(Problem lp, int Row, int column, float Value)
274 int elmnr, lastelm, i;
275 // System.err.println("lp.mat.length = " + lp.mat.length);
277 if(Row > lp.rows || Row < 0)
278 System.err.print("Row out of range");
279 if(column > lp.columns || column < 1)
280 System.err.print("column out of range");
281 if(lp.scaling_used != FALSE)
282 Value *= lp.scale[Row] * lp.scale[lp.rows + column];
284 if(true /*abs(Value) > lp.epsilon*/)
286 if (lp.basis[column] == TRUE && Row > 0)
287 lp.basis_valid = FALSE;
288 lp.eta_valid = FALSE;
289 elmnr = lp.col_end[column-1];
290 while((elmnr < lp.col_end[column]) ?
291 (lp.mat[elmnr].row_nr != Row) : false)
294 if((elmnr != lp.col_end[column]) ?
295 (lp.mat[elmnr].row_nr == Row) : false )
296 if (lp.scaling_used != FALSE)
298 if (lp.ch_sign[Row] != FALSE)
299 lp.mat[elmnr].value =
300 -Value * lp.scale[Row] * lp.scale[column];
302 lp.mat[elmnr].value =
303 Value * lp.scale[Row] * lp.scale[column];
307 if (lp.ch_sign[Row] != FALSE)
308 lp.mat[elmnr].value = -Value;
310 lp.mat[elmnr].value = Value;
314 /* check if more space is needed for matrix */
317 /* Shift the matrix */
318 lastelm=lp.non_zeros;
319 for(i = lastelm; i > elmnr ; i--)
320 lp.mat[i]=lp.mat[i-1];
321 for(i = column; i <= lp.columns; i++)
324 /* Set new element */
325 lp.mat[elmnr].row_nr=Row;
327 if (lp.scaling_used != FALSE)
329 if (lp.ch_sign[Row] != FALSE)
330 lp.mat[elmnr].value=-Value*lp.scale[Row]*lp.scale[column];
332 lp.mat[elmnr].value=Value*lp.scale[Row]*lp.scale[column];
336 if (lp.ch_sign[Row] != FALSE)
337 lp.mat[elmnr].value=-Value;
339 lp.mat[elmnr].value=Value;
342 lp.row_end_valid=FALSE;
345 if (lp.active != FALSE)
346 Non_zeros=lp.non_zeros;
351 public void set_obj_fn(Problem lp, float[] row)
353 for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
357 for(i = 1; i <= lp.columns; i++)
358 set_mat(lp, 0, i, row[i]);
362 public void add_constraint(Problem lp, float[] row, short constr_type, float rh)
364 for(int i=row.length-1; i>0; i--) row[i] = row[i-1];
373 addtoo = new short[lp.columns + 1];
375 for(i = 1; i <= lp.columns; i++)
384 newmat = new Matrix[lp.non_zeros];
385 for (i = 0; i < newmat.length; i++)
386 newmat[i] = new Matrix(0, (float)0);
388 inc_mat_space(lp, 0);
393 if (lp.scaling_used != FALSE)
396 for(i=lp.sum; i > lp.rows; i--)
397 lp.scale[i]=lp.scale[i-1];
401 if (lp.names_used != FALSE)
402 lp.row_name[lp.rows] = new String("r_" + lp.rows);
404 if(lp.scaling_used != FALSE && lp.columns_scaled != FALSE)
405 for(i = 1; i <= lp.columns; i++)
406 row[i] *= lp.scale[lp.rows+i];
409 lp.ch_sign[lp.rows] = TRUE;
411 lp.ch_sign[lp.rows] = FALSE;
415 for(i = 1; i <= lp.columns; i++)
417 for(j = stcol; j < lp.col_end[i]; j++)
419 newmat[elmnr].row_nr=lp.mat[j].row_nr;
420 newmat[elmnr].value=lp.mat[j].value;
423 if(addtoo[i] != FALSE)
425 if(lp.ch_sign[lp.rows] != FALSE)
426 newmat[elmnr].value = -row[i];
428 newmat[elmnr].value = row[i];
429 newmat[elmnr].row_nr = lp.rows;
438 for(i=lp.sum ; i > lp.rows; i--)
440 lp.orig_upbo[i]=lp.orig_upbo[i-1];
441 lp.orig_lowbo[i]=lp.orig_lowbo[i-1];
442 lp.basis[i]=lp.basis[i-1];
443 lp.lower[i]=lp.lower[i-1];
446 for(i= 1 ; i <= lp.rows; i++)
447 if(lp.bas[i] >= lp.rows)
450 if(constr_type==LE || constr_type==GE)
452 lp.orig_upbo[lp.rows]=lp.infinite;
454 else if(constr_type==EQ)
456 lp.orig_upbo[lp.rows]=0;
460 System.err.print("Wrong constraint type\n");
464 lp.orig_lowbo[lp.rows]=0;
466 if(constr_type==GE && rh != 0)
467 lp.orig_rh[lp.rows]=-rh;
469 lp.orig_rh[lp.rows]=rh;
471 lp.row_end_valid=FALSE;
473 lp.bas[lp.rows]=lp.rows;
474 lp.basis[lp.rows]=TRUE;
475 lp.lower[lp.rows]=TRUE;
477 if (lp.active != FALSE)
482 public void del_constraint(Problem lp, int del_row)
488 if(del_row < 1 || del_row > lp.rows)
490 System.err.println("There is no constraint nr. " + del_row);
497 for(i = 1; i <= lp.columns; i++)
499 for(j=startcol; j < lp.col_end[i]; j++)
501 if(lp.mat[j].row_nr!=del_row)
503 lp.mat[elmnr]=lp.mat[j];
504 if(lp.mat[elmnr].row_nr > del_row)
505 lp.mat[elmnr].row_nr--;
511 startcol=lp.col_end[i];
514 for(i = del_row; i < lp.rows; i++)
516 lp.orig_rh[i] = lp.orig_rh[i+1];
517 lp.ch_sign[i] = lp.ch_sign[i+1];
518 lp.bas[i] = lp.bas[i+1];
519 if (lp.names_used != FALSE)
520 lp.row_name[i] = lp.row_name[i+1];
522 for(i = 1; i < lp.rows; i++)
523 if(lp.bas[i] > del_row)
526 for(i=del_row; i < lp.sum; i++)
528 lp.lower[i]=lp.lower[i+1];
529 lp.basis[i]=lp.basis[i+1];
530 lp.orig_upbo[i]=lp.orig_upbo[i+1];
531 lp.orig_lowbo[i]=lp.orig_lowbo[i+1];
532 if (lp.scaling_used != FALSE)
533 lp.scale[i]=lp.scale[i+1];
539 lp.row_end_valid=FALSE;
541 if (lp.active != FALSE)
544 lp.basis_valid=FALSE;
548 public void add_column(Problem lp, float[] column)
555 inc_mat_space(lp, lp.rows+1);
557 if (lp.scaling_used != FALSE)
559 for(i = 0; i <= lp.rows; i++)
560 column[i]*=lp.scale[i];
564 elmnr=lp.col_end[lp.columns-1];
565 for(i = 0 ; i <= lp.rows ; i++)
568 lp.mat[elmnr].row_nr=i;
569 if(lp.ch_sign[i] != FALSE)
570 lp.mat[elmnr].value=-column[i];
572 lp.mat[elmnr].value=column[i];
576 lp.col_end[lp.columns]=elmnr;
577 lp.orig_lowbo[lp.sum]=0;
578 lp.orig_upbo[lp.sum]=lp.infinite;
579 lp.lower[lp.sum]=TRUE;
580 lp.basis[lp.sum]=FALSE;
581 if (lp.names_used != FALSE)
582 lp.col_name[lp.columns] = 0;
585 lp.row_end_valid=FALSE;
587 if (lp.active != FALSE)
591 Non_zeros=lp.non_zeros;
595 public void del_column(Problem lp, int column)
597 int i, j, from_elm, to_elm, elm_in_col;
598 if(column > lp.columns || column < 1)
599 System.err.print("column out of range in del_column");
600 for(i = 1; i <= lp.rows; i++)
602 if(lp.bas[i]==lp.rows+column)
603 lp.basis_valid=FALSE;
604 else if(lp.bas[i] > lp.rows+column)
607 for(i = lp.rows+column; i < lp.sum; i++)
609 if (lp.names_used != FALSE)
610 lp.col_name[i-lp.rows] = lp.col_name[i-lp.rows+1];
611 lp.orig_upbo[i]=lp.orig_upbo[i+1];
612 lp.orig_lowbo[i]=lp.orig_lowbo[i+1];
613 lp.upbo[i]=lp.upbo[i+1];
614 lp.lowbo[i]=lp.lowbo[i+1];
615 lp.basis[i]=lp.basis[i+1];
616 lp.lower[i]=lp.lower[i+1];
617 if (lp.scaling_used != FALSE)
618 lp.scale[i]=lp.scale[i+1];
620 for(i = 0; i < lp.nr_lagrange; i++)
621 for(j = column; j <= lp.columns; j++)
622 lp.lag_row[i][j]=lp.lag_row[i][j+1];
623 to_elm=lp.col_end[column-1];
624 from_elm=lp.col_end[column];
625 elm_in_col=from_elm-to_elm;
626 for(i = from_elm; i < lp.non_zeros; i++)
628 lp.mat[to_elm]=lp.mat[i];
631 for(i = column; i < lp.columns; i++)
632 lp.col_end[i]=lp.col_end[i+1]-elm_in_col;
633 lp.non_zeros -= elm_in_col;
634 lp.row_end_valid=FALSE;
639 if (lp.active != FALSE)
643 public void bound_sum(Problem lp, int column1, int column2, float bound, short type, float[] scratch) {
644 for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
645 scratch[column1] = (float)1.0;
646 scratch[column2] = (float)1.0;
647 add_constraint(lp, scratch, type, bound);
648 for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
651 public void bound_difference(Problem lp, int column1, int column2, float bound, short type, float[] scratch) {
652 for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
653 scratch[column1] = (float)1.0;
654 scratch[column2] = (float)-1.0;
655 add_constraint(lp, scratch, type, bound);
656 for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
659 public void set_upbo(Problem lp, int column, float value)
661 if(column > lp.columns || column < 1)
662 System.err.print("column out of range");
663 if (lp.scaling_used != FALSE)
664 value /= lp.scale[lp.rows + column];
665 if(value < lp.orig_lowbo[lp.rows + column])
666 System.err.print("UpperBound must be >= lowerBound");
667 lp.eta_valid = FALSE;
668 lp.orig_upbo[lp.rows+column] = value;
671 public void set_lowbo(Problem lp, int column, float value)
673 if(column > lp.columns || column < 1)
674 System.err.print("column out of range");
675 if (lp.scaling_used != FALSE)
676 value /= lp.scale[lp.rows + column];
677 if(value > lp.orig_upbo[lp.rows + column])
678 System.err.print("UpperBound must be >= lowerBound");
679 lp.eta_valid = FALSE;
680 lp.orig_lowbo[lp.rows+column] = value;
683 public void set_rh(Problem lp, int row, float value)
685 if(row > lp.rows || row < 0)
686 System.err.print("Row out of Range");
687 if(row == 0) /* setting of RHS of OF not meaningful */
689 System.err.println("Warning: attempt to set RHS of objective function, ignored");
692 if (lp.scaling_used != FALSE)
693 if (lp.ch_sign[row] != FALSE)
694 lp.orig_rh[row] = -value * lp.scale[row];
696 lp.orig_rh[row] = value * lp.scale[row];
698 if (lp.ch_sign[row] != FALSE)
699 lp.orig_rh[row] = -value;
701 lp.orig_rh[row] = value;
702 lp.eta_valid = FALSE;
705 public void set_rh_vec(Problem lp, float[] rh)
708 if (lp.scaling_used != FALSE)
709 for(i = 1; i <= lp.rows; i++)
710 if(lp.ch_sign[i] != FALSE)
711 lp.orig_rh[i]=-rh[i]*lp.scale[i];
713 lp.orig_rh[i]=rh[i]*lp.scale[i];
715 for(i=1; i <= lp.rows; i++)
716 if(lp.ch_sign[i] != FALSE)
717 lp.orig_rh[i]=-rh[i];
723 public void set_maxim(Problem lp)
726 if(lp.maximise==FALSE)
728 for(i = 0; i < lp.non_zeros; i++)
729 if(lp.mat[i].row_nr==0)
735 if (lp.active != FALSE)
739 public void set_minim(Problem lp)
742 if(lp.maximise==TRUE)
744 for(i = 0; i < lp.non_zeros; i++)
745 if(lp.mat[i].row_nr==0)
746 lp.mat[i].value = -lp.mat[i].value;
751 if (lp.active != FALSE)
755 public void set_constr_type(Problem lp, int row, short con_type)
758 if(row > lp.rows || row < 1)
759 System.err.print("Row out of Range");
763 lp.basis_valid=FALSE;
764 if (lp.ch_sign[row] != FALSE)
766 for(i = 0; i < lp.non_zeros; i++)
767 if(lp.mat[i].row_nr==row)
770 lp.ch_sign[row]=FALSE;
771 if(lp.orig_rh[row]!=0)
775 else if(con_type==LE)
777 lp.orig_upbo[row]=lp.infinite;
778 lp.basis_valid=FALSE;
779 if (lp.ch_sign[row] != FALSE)
781 for(i = 0; i < lp.non_zeros; i++)
782 if(lp.mat[i].row_nr==row)
785 lp.ch_sign[row]=FALSE;
786 if(lp.orig_rh[row]!=0)
790 else if(con_type==GE)
792 lp.orig_upbo[row]=lp.infinite;
793 lp.basis_valid=FALSE;
794 if(lp.ch_sign[row] == FALSE)
796 for(i = 0; i < lp.non_zeros; i++)
797 if(lp.mat[i].row_nr==row)
800 lp.ch_sign[row]=TRUE;
801 if(lp.orig_rh[row]!=0)
806 System.err.print("Constraint type not (yet) implemented");
809 public float mat_elm(Problem lp, int row, int column)
813 if(row < 0 || row > lp.rows)
814 System.err.print("Row out of range in mat_elm");
815 if(column < 1 || column > lp.columns)
816 System.err.print("column out of range in mat_elm");
818 elmnr=lp.col_end[column-1];
819 while(lp.mat[elmnr].row_nr != row && elmnr < lp.col_end[column])
821 if(elmnr != lp.col_end[column])
823 value = lp.mat[elmnr].value;
824 if (lp.ch_sign[row] != FALSE)
826 if (lp.scaling_used != FALSE)
827 value /= lp.scale[row] * lp.scale[lp.rows + column];
833 public void get_row(Problem lp, int row_nr, float[] row)
837 if(row_nr < 0 || row_nr > lp.rows)
838 System.err.print("Row nr. out of range in get_row");
839 for(i = 1; i <= lp.columns; i++)
842 for(j=lp.col_end[i-1]; j < lp.col_end[i]; j++)
843 if(lp.mat[j].row_nr==row_nr)
844 row[i]=lp.mat[j].value;
845 if (lp.scaling_used != FALSE)
846 row[i] /= lp.scale[lp.rows+i] * lp.scale[row_nr];
848 if(lp.ch_sign[row_nr] != FALSE)
849 for(i=0; i <= lp.columns; i++)
854 public void get_column(Problem lp, int col_nr, float[] column)
858 if(col_nr < 1 || col_nr > lp.columns)
859 System.err.print("Col. nr. out of range in get_column");
860 for(i = 0; i <= lp.rows; i++)
862 for(i = lp.col_end[col_nr-1]; i < lp.col_end[col_nr]; i++)
863 column[lp.mat[i].row_nr]=lp.mat[i].value;
864 for(i = 0; i <= lp.rows; i++)
867 if(lp.ch_sign[i] != FALSE)
869 if (lp.scaling_used != FALSE)
870 column[i]/=(lp.scale[i] * lp.scale[lp.rows+col_nr]);
874 public void get_reduced_costs(Problem lp, float[] rc)
879 if(lp.basis_valid == FALSE)
880 System.err.print("Not a valid basis in get_reduced_costs");
882 if(lp.eta_valid == FALSE)
884 for(i = 1; i <= lp.sum; i++)
888 for(i = 1; i <= lp.columns; i++)
891 if(Basis[varnr] == FALSE)
895 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
896 f += rc[Mat[j].row_nr] * Mat[j].value;
900 for(i = 1; i <= Sum; i++)
901 rc[i] = round(rc[i], Epsd);
904 short is_feasible(Problem lp, float[] values)
910 if (lp.scaling_used != FALSE)
912 for(i = lp.rows+1; i <= lp.sum; i++)
913 if( values[i - lp.rows] < lp.orig_lowbo[i]*lp.scale[i]
914 || values[i - lp.rows] > lp.orig_upbo[i]*lp.scale[i])
919 for(i = lp.rows+1; i <= lp.sum; i++)
920 if( values[i - lp.rows] < lp.orig_lowbo[i]
921 || values[i - lp.rows] > lp.orig_upbo[i])
924 this_rhs = new float[lp.rows+1];
925 for (i = 0; i <= lp.rows; i++)
928 for(i = 1; i <= lp.columns; i++)
929 for(elmnr = lp.col_end[i - 1]; elmnr < lp.col_end[i]; elmnr++)
930 this_rhs[lp.mat[elmnr].row_nr] += lp.mat[elmnr].value * values[i];
931 for(i = 1; i <= lp.rows; i++)
933 dist = lp.orig_rh[i] - this_rhs[i];
934 dist = round(dist, (float)0.001);
935 if((lp.orig_upbo[i] == 0 && dist != 0) || dist < 0)
943 short column_in_lp(Problem lp, float[] testcolumn)
949 if (lp.scaling_used != FALSE)
950 for(i = 1; i <= lp.columns; i++)
954 while(ident != FALSE && (j < lp.col_end[i]))
956 value = lp.mat[j].value;
957 if(lp.ch_sign[lp.mat[j].row_nr] != FALSE)
959 value /= lp.scale[lp.rows+i];
960 value /= lp.scale[lp.mat[j].row_nr];
961 value -= testcolumn[lp.mat[j].row_nr];
962 if(Math.abs(value) > (float)0.001) /* should be some epsilon? MB */
970 for(i = 1; i <= lp.columns; i++)
974 while(ident != FALSE && j < lp.col_end[i])
976 value = lp.mat[j].value;
977 if(lp.ch_sign[lp.mat[j].row_nr] != FALSE)
979 value -= testcolumn[lp.mat[j].row_nr];
980 if( Math.abs(value) > (float)0.001 )
990 private float minmax_to_scale(float min, float max)
994 /* should do something sensible when min or max is 0, MB */
995 if((min == 0) || (max == 0))
998 scale = (float)(1 / Math.pow(Math.E, (Math.log(min) + Math.log(max)) / 2));
1002 public void unscale_columns(Problem lp)
1007 for(j = 1; j <= lp.columns; j++)
1008 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
1009 lp.mat[i].value /= lp.scale[lp.rows + j];
1011 /* unscale Bounds as well */
1012 for(i = lp.rows + 1; i < lp.sum; i++)
1014 if(lp.orig_lowbo[i] != 0)
1015 lp.orig_lowbo[i] *= lp.scale[i];
1016 if(lp.orig_upbo[i] != lp.infinite)
1017 lp.orig_upbo[i] *= lp.scale[i];
1020 for(i=lp.rows+1; i<= lp.sum; i++)
1022 lp.columns_scaled=FALSE;
1026 public void unscale(Problem lp)
1030 if (lp.scaling_used != FALSE)
1034 for(j = 1; j <= lp.columns; j++)
1035 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
1036 lp.mat[i].value /= lp.scale[lp.rows + j];
1038 /* unscale Bounds */
1039 for(i = lp.rows + 1; i < lp.sum; i++)
1041 if(lp.orig_lowbo[i] != 0)
1042 lp.orig_lowbo[i] *= lp.scale[i];
1043 if(lp.orig_upbo[i] != lp.infinite)
1044 lp.orig_upbo[i] *= lp.scale[i];
1047 /* unscale the matrix */
1048 for(j = 1; j <= lp.columns; j++)
1049 for(i = lp.col_end[j-1]; i < lp.col_end[j]; i++)
1050 lp.mat[i].value /= lp.scale[lp.mat[i].row_nr];
1052 /* unscale the rhs! */
1053 for(i = 0; i <= lp.rows; i++)
1054 lp.orig_rh[i] /= lp.scale[i];
1056 lp.scaling_used=FALSE;
1062 public void auto_scale(Problem lp)
1064 int i, j, row_nr, IntUsed;
1065 float row_max[], row_min[], scalechange[], absval;
1067 if(lp.scaling_used == 0)
1069 lp.scale = new float[lp.sum_alloc + 1];
1070 for(i = 0; i <= lp.sum; i++)
1074 row_max = new float[lp.rows + 1];
1075 row_min = new float[lp.rows + 1];
1076 scalechange = new float[lp.sum + 1];
1078 /* initialise min and max values */
1079 for(i = 0; i <= lp.rows; i++)
1082 row_min[i]=lp.infinite;
1085 /* calculate min and max absolute values of rows */
1086 for(j = 1; j <= lp.columns; j++)
1087 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
1089 row_nr = lp.mat[i].row_nr;
1090 absval = Math.abs(lp.mat[i].value);
1093 row_max[row_nr] = Math.max(row_max[row_nr], absval);
1094 row_min[row_nr] = Math.min(row_min[row_nr], absval);
1097 /* calculate scale factors for rows */
1098 for(i = 0; i <= lp.rows; i++)
1100 scalechange[i] = minmax_to_scale(row_min[i], row_max[i]);
1101 lp.scale[i] *= scalechange[i];
1104 /* now actually scale the matrix */
1105 for(j = 1; j <= lp.columns; j++)
1106 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
1107 lp.mat[i].value *= scalechange[lp.mat[i].row_nr];
1109 /* and scale the rhs! */
1110 for(i = 0; i <= lp.rows; i++)
1111 lp.orig_rh[i] *= scalechange[i];
1115 while(IntUsed == FALSE && i <= lp.sum)
1120 if(IntUsed == FALSE)
1122 float col_max, col_min;
1124 /* calculate scales */
1125 for(j = 1; j <= lp.columns; j++)
1128 col_min = lp.infinite;
1129 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
1131 if(lp.mat[i].value!=0)
1133 col_max = Math.max(col_max, Math.abs(lp.mat[i].value));
1134 col_min = Math.min(col_min, Math.abs(lp.mat[i].value));
1137 scalechange[lp.rows + j] = minmax_to_scale(col_min, col_max);
1138 lp.scale[lp.rows + j] *= scalechange[lp.rows + j];
1142 for(j = 1; j <= lp.columns; j++)
1143 for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++)
1144 lp.mat[i].value *= scalechange[lp.rows + j];
1146 /* scale Bounds as well */
1147 for(i = lp.rows + 1; i < lp.sum; i++)
1149 if(lp.orig_lowbo[i] != 0)
1150 lp.orig_lowbo[i] /= scalechange[i];
1151 if(lp.orig_upbo[i] != lp.infinite)
1152 lp.orig_upbo[i] /= scalechange[i];
1154 lp.columns_scaled=TRUE;
1156 lp.scaling_used=TRUE;
1160 public void reset_basis(Problem lp) { lp.basis_valid=FALSE; }
1162 /* Globals used by solver */
1169 void set_globals(Problem lp)
1175 columns = lp.columns;
1176 Sum = Rows + columns;
1177 Non_zeros = lp.non_zeros;
1180 Col_end = lp.col_end;
1181 Row_end = lp.row_end;
1184 Orig_rh = lp.orig_rh;
1185 Orig_upbo = lp.orig_upbo;
1186 Orig_lowbo = lp.orig_lowbo;
1192 Eta_alloc = lp.eta_alloc;
1193 Eta_size = lp.eta_size;
1194 Num_inv = lp.num_inv;
1195 Eta_value = lp.eta_value;
1196 Eta_row_nr = lp.eta_row_nr;
1197 Eta_col_end = lp.eta_col_end;
1198 Solution = lp.solution;
1199 Best_solution = lp.best_solution;
1200 Infinite = lp.infinite;
1201 Epsilon = lp.epsilon;
1210 Maximise = lp.maximise;
1211 Floor_first = lp.floor_first;
1214 // System.out.println("Sum = " + Sum);
1215 } // end of set_global
1217 private void ftran(int start,
1224 for(i = start; i <= end; i++)
1226 k = Eta_col_end[i] - 1;
1230 for(j = Eta_col_end[i - 1]; j < k; j++)
1231 pcol[Eta_row_nr[j]] += theta * Eta_value[j]; /* cpu expensive line */
1232 pcol[r] *= Eta_value[k];
1234 /* round small values to zero */
1235 for(i = 0; i <= Rows; i++)
1236 round(pcol[i], Epsel);
1239 private void btran(float[] row)
1244 for(i = Eta_size; i >= 1; i--)
1247 k = Eta_col_end[i] - 1;
1248 for(j = Eta_col_end[i - 1]; j <= k; j++)
1249 f += row[Eta_row_nr[j]] * Eta_value[j];
1250 f = round(f, Epsel);
1251 row[Eta_row_nr[k]] = f;
1256 static short Isvalid(Problem lp)
1258 int i, j, rownum[], colnum[];
1261 if(lp.row_end_valid == FALSE)
1263 num = new int[lp.rows + 1];
1264 rownum = new int[lp.rows + 1];
1265 for(i = 0; i <= lp.rows; i++)
1270 for(i = 0; i < lp.non_zeros; i++)
1271 rownum[lp.mat[i].row_nr]++;
1273 for(i = 1; i <= lp.rows; i++)
1274 lp.row_end[i] = lp.row_end[i - 1] + rownum[i];
1275 for(i = 1; i <= lp.columns; i++)
1276 for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++)
1278 row_nr = lp.mat[j].row_nr;
1282 lp.col_no[lp.row_end[row_nr - 1] + num[row_nr]] = i;
1285 lp.row_end_valid = TRUE;
1287 if(lp.valid != FALSE)
1289 rownum = new int[lp.rows + 1];
1290 for (i = 0; i <= lp.rows; i++)
1292 colnum = new int[lp.columns + 1];
1293 for (i = 0; i <= lp.columns; i++)
1295 for(i = 1 ; i <= lp.columns; i++)
1296 for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++)
1299 rownum[lp.mat[j].row_nr]++;
1301 for(i = 1; i <= lp.columns; i++)
1303 if (lp.names_used != FALSE)
1304 System.err.print("Warning: Variable " + lp.col_name[i] +
1305 " not used in any constraints\n");
1307 System.err.print("Warning: Variable " + i + " not used in any constaints\n");
1312 private void resize_eta()
1315 float[] db_ptr = Eta_value;
1316 Eta_value = new float[Eta_alloc];
1317 System.arraycopy(db_ptr, 0, Eta_value, 0, db_ptr.length);
1318 Lp.eta_value = Eta_value;
1320 int[] int_ptr = Eta_row_nr;
1321 Eta_row_nr = new int[Eta_alloc];
1322 System.arraycopy(int_ptr, 0, Eta_row_nr, 0, int_ptr.length);
1323 Lp.eta_row_nr = Eta_row_nr;
1327 private void condensecol(int row_nr,
1332 elnr = Eta_col_end[Eta_size];
1334 if(elnr + Rows + 2 > Eta_alloc) /* maximum local growth of Eta */
1337 for(i = 0; i <= Rows; i++)
1338 if(i != row_nr && pcol[i] != 0)
1340 Eta_row_nr[elnr] = i;
1341 Eta_value[elnr] = pcol[i];
1344 Eta_row_nr[elnr] = row_nr;
1345 Eta_value[elnr] = pcol[row_nr];
1347 Eta_col_end[Eta_size + 1] = elnr;
1351 private void addetacol()
1356 j = Eta_col_end[Eta_size];
1358 k = Eta_col_end[Eta_size];
1359 theta = 1 / (float) Eta_value[k - 1];
1360 Eta_value[k - 1] = theta;
1361 for(i = j; i < k - 1; i++)
1362 Eta_value[i] *= -theta;
1363 JustInverted = FALSE;
1366 private void setpivcol(short lower,
1377 for(i = 0; i <= Rows; i++)
1381 colnr = varin - Rows;
1382 for(i = Col_end[colnr - 1]; i < Col_end[colnr]; i++)
1383 pcol[Mat[i].row_nr] = Mat[i].value * f;
1384 pcol[0] -= Extrad * f;
1391 ftran(1, Eta_size, pcol);
1395 private void minoriteration(int colnr,
1398 int i, j, k, wk, varin, varout, elnr;
1399 float piv = 0, theta;
1401 varin = colnr + Rows;
1402 elnr = Eta_col_end[Eta_size];
1407 Eta_row_nr[elnr] = 0;
1408 Eta_value[elnr] = -Extrad;
1411 for(j = Col_end[colnr - 1] ; j < Col_end[colnr]; j++)
1414 if(k == 0 && Extrad != 0)
1415 Eta_value[Eta_col_end[Eta_size -1]] += Mat[j].value;
1416 else if(k != row_nr)
1418 Eta_row_nr[elnr] = k;
1419 Eta_value[elnr] = Mat[j].value;
1425 Eta_row_nr[elnr] = row_nr;
1426 Eta_value[elnr] = 1 / (float) piv;
1428 theta = Rhs[row_nr] / (float) piv;
1429 Rhs[row_nr] = theta;
1430 for(i = wk; i < elnr - 1; i++)
1431 Rhs[Eta_row_nr[i]] -= theta * Eta_value[i];
1432 varout = Bas[row_nr];
1433 Bas[row_nr] = varin;
1434 Basis[varout] = FALSE;
1435 Basis[varin] = TRUE;
1436 for(i = wk; i < elnr - 1; i++)
1437 Eta_value[i] /= - (float) piv;
1438 Eta_col_end[Eta_size] = elnr;
1439 } /* minoriteration */
1441 private void rhsmincol(float theta,
1445 int i, j, k, varout;
1448 if(row_nr > Rows + 1)
1450 System.err.print("Error: rhsmincol called with row_nr: " +
1451 row_nr + ", rows: " + Rows + "\n");
1452 System.err.print("This indicates numerical instability\n");
1455 j = Eta_col_end[Eta_size];
1456 k = Eta_col_end[Eta_size + 1];
1457 for(i = j; i < k; i++)
1459 f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i];
1461 Rhs[Eta_row_nr[i]] = f;
1463 Rhs[row_nr] = theta;
1464 varout = Bas[row_nr];
1465 Bas[row_nr] = varin;
1466 Basis[varout] = FALSE;
1467 Basis[varin] = TRUE;
1473 int i, j, v, wk, numit, varnr, row_nr, colnr, varin;
1478 int[] rownum, col, row;
1481 rownum = new int[Rows + 1];
1482 for (i = 0; i <= Rows; i++)
1484 col = new int[Rows + 1];
1485 for (i = 0; i <= Rows; i++)
1487 row = new int[Rows + 1];
1488 for (i = 0; i <= Rows; i++)
1490 pcol = new float[Rows + 1];
1491 for (i = 0; i <= Rows; i++)
1493 frow = new short[Rows + 1];
1494 for(i = 0; i <= Rows; i++)
1496 fcol = new short[columns + 1];
1497 for(i = 0; i < columns; i++)
1499 colnum = new int[columns + 1];
1500 for(i = 0; i <= columns; i++)
1503 for(i = 0; i <= Rows; i++)
1505 fcol[Bas[i] - Rows - 1] = TRUE;
1507 frow[Bas[i]] = FALSE;
1509 for(i = 1; i <= Rows; i++)
1510 if(frow[i] != FALSE)
1511 for(j = Row_end[i - 1] + 1; j <= Row_end[i]; j++)
1514 if(fcol[wk - 1] != FALSE)
1520 for(i = 1; i <= Rows; i++)
1522 for(i = 1; i <= Rows; i++)
1524 for(i = 1; i <= columns; i++)
1525 Basis[i + Rows] = FALSE;
1527 for(i = 0; i <= Rows; i++)
1530 for(i = 1; i <= columns; i++)
1533 if(Lower[varnr] == FALSE)
1535 theta = Upbo[varnr];
1536 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1537 Rhs[Mat[j].row_nr] -= theta * Mat[j].value;
1540 for(i = 1; i <= Rows; i++)
1541 if(Lower[i] == FALSE)
1554 if(rownum[row_nr - 1] == 1)
1555 if(frow[row_nr] != FALSE)
1558 j = Row_end[row_nr - 1] + 1;
1559 while(fcol[Col_no[j] - 1] == FALSE)
1562 fcol[colnr - 1] = FALSE;
1564 for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++)
1565 if(frow[Mat[j].row_nr] != FALSE)
1566 rownum[Mat[j].row_nr - 1]--;
1567 frow[row_nr] = FALSE;
1568 minoriteration(colnr, row_nr);
1579 if(colnum[colnr] == 1)
1580 if(fcol[colnr - 1] != FALSE)
1583 j = Col_end[colnr - 1] + 1;
1584 while(frow[Mat[j - 1].row_nr] == FALSE)
1586 row_nr = Mat[j - 1].row_nr;
1587 frow[row_nr] = FALSE;
1588 rownum[row_nr - 1] = 0;
1589 for(j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++)
1590 if(fcol[Col_no[j] - 1] != FALSE)
1591 colnum[Col_no[j]]--;
1592 fcol[colnr - 1] = FALSE;
1594 col[numit - 1] = colnr;
1595 row[numit - 1] = row_nr;
1598 for(j = 1; j <= columns; j++)
1599 if(fcol[j - 1] != FALSE)
1601 fcol[j - 1] = FALSE;
1602 setpivcol(Lower[Rows + j], j + Rows, pcol);
1604 while((frow[row_nr] == FALSE || pcol[row_nr] == FALSE)
1606 row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes
1607 rhsmincol crash. Solved in 2.0? MB */
1608 if(row_nr == Rows + 1)
1609 System.err.print("Inverting failed");
1610 frow[row_nr] = FALSE;
1611 condensecol(row_nr, pcol);
1612 theta = Rhs[row_nr] / (float) pcol[row_nr];
1613 rhsmincol(theta, row_nr, Rows + j);
1616 for(i = numit - 1; i >= 0; i--)
1620 varin = colnr + Rows;
1621 for(j = 0; j <= Rows; j++)
1623 for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++)
1624 pcol[Mat[j].row_nr] = Mat[j].value;
1626 condensecol(row_nr, pcol);
1627 theta = Rhs[row_nr] / (float) pcol[row_nr];
1628 rhsmincol(theta, row_nr, varin);
1631 for(i = 1; i <= Rows; i++)
1632 Rhs[i] = round(Rhs[i], Epsb);
1633 JustInverted = TRUE;
1637 private short colprim(Ref colnr,
1648 for(i = 1; i <= Sum; i++)
1652 for(i = 1; i <= columns; i++)
1655 if(Basis[varnr] == FALSE)
1659 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1660 f += drow[Mat[j].row_nr] * Mat[j].value;
1664 for(i = 1; i <= Sum; i++)
1665 drow[i] = round(drow[i], Epsd);
1668 System.out.println("dpiv = " + dpiv);
1669 System.out.println("drow[] at colprim");
1670 for(i = 1; i <= Sum; i++) // DEBUG
1672 System.out.print(drow[i] + " ");
1674 System.out.println();
1676 for(i = 1; i <= Sum; i++)
1677 if(Basis[i] == FALSE)
1680 if(Lower[i] != FALSE)
1690 if(Lp.trace != FALSE)
1692 System.err.print("col_prim:" + colnr.value + ", reduced cost: " + dpiv + "\n");
1694 System.err.print("col_prim: no negative reduced costs found, optimality!\n");
1695 if(colnr.value == 0)
1701 return(colnr.value > 0 ? (short)1 : (short)0);
1704 private short rowprim(int colnr,
1713 theta.value = Infinite;
1714 for(i = 1; i <= Rows; i++)
1717 if(Math.abs(f) < TREJ)
1721 quot = 2 * Infinite;
1723 quot = Rhs[i] / (float) f;
1725 if(Upbo[Bas[i]] < Infinite)
1726 quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
1728 if(quot < theta.value)
1735 if(row_nr.value == 0)
1736 for(i = 1; i <= Rows; i++)
1741 quot = 2 * Infinite;
1743 quot = Rhs[i] / (float) f;
1745 if(Upbo[Bas[i]] < Infinite)
1746 quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f;
1747 quot = round(quot, Epsel);
1748 if(quot < theta.value)
1758 System.err.println("Warning: Numerical instability, qout = " + theta.value);
1759 System.err.println("pcol[" + row_nr.value + "] = " + f + ", rhs[" +
1760 row_nr.value + "] = " + Rhs[(int)row_nr.value] +
1761 " , upbo = " + Upbo[Bas[(int)row_nr.value]]);
1763 if(row_nr.value == 0)
1765 if(Upbo[colnr] == Infinite)
1774 while(pcol[i] >= 0 && i <= Rows)
1776 if(i > Rows) /* empty column with upperBound! */
1778 Lower[colnr] = FALSE;
1779 Rhs[0] += Upbo[colnr]*pcol[0];
1789 if(row_nr.value > 0)
1791 if(Lp.trace != FALSE)
1792 System.out.println("row_prim:" + row_nr.value + ", pivot element:" +
1793 pcol[(int)row_nr.value]);
1795 return((row_nr.value > 0) ? (short)1 : (short)0);
1798 private short rowdual(Ref row_nr)
1808 while(i < Rows && artifs == FALSE)
1812 if(f == 0 && (Rhs[i] != 0))
1819 if(Rhs[i] < f - Rhs[i])
1831 if(Lp.trace != FALSE)
1833 if(row_nr.value > 0)
1835 System.out.println("row_dual:" + row_nr.value +
1836 ", rhs of selected row: "
1837 + Rhs[(int)row_nr.value]);
1838 if(Upbo[Bas[(int)row_nr.value]] < Infinite)
1839 System.out.println(" upper Bound of basis variable: " +
1840 Upbo[Bas[(int)row_nr.value]]);
1843 System.out.print("row_dual: no infeasibilities found\n");
1846 return(row_nr.value > 0 ? (short)1 : (short)0);
1849 private short coldual(int row_nr,
1856 float theta, quot, pivot, d, f, g;
1861 for(i = 0; i <= Rows; i++)
1868 for(i = Eta_size; i >= 1; i--)
1872 r = Eta_row_nr[Eta_col_end[i] - 1];
1873 for(j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++)
1875 /* this is where the program consumes most cpu time */
1876 f += prow[Eta_row_nr[j]] * Eta_value[j];
1877 d += drow[Eta_row_nr[j]] * Eta_value[j];
1879 f = round(f, Epsel);
1881 d = round(d, Epsel);
1884 for(i = 1; i <= columns; i++)
1887 if(Basis[varnr] == FALSE)
1889 d = - Extrad * drow[0];
1891 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1893 d = d + drow[Mat[j].row_nr] * Mat[j].value;
1894 f = f + prow[Mat[j].row_nr] * Mat[j].value;
1900 for(i = 0; i <= Sum; i++)
1902 prow[i] = round(prow[i], Epsel);
1903 drow[i] = round(drow[i], Epsd);
1906 if(Rhs[row_nr] > Upbo[Bas[row_nr]])
1913 for(i = 1; i <= Sum; i++)
1915 if(Lower[i] != FALSE)
1919 if((d < 0) && (Basis[i] == FALSE) && (Upbo[i] > 0))
1921 if(Lower[i] == FALSE)
1922 quot = -drow[i] / (float) d;
1924 quot = drow[i] / (float) d;
1931 else if((quot == theta) && (Math.abs(d) > Math.abs(pivot)))
1938 if(Lp.trace != FALSE)
1939 System.out.println("col_dual:" + colnr.value + ", pivot element: " +
1940 prow[(int)colnr.value]);
1945 return(colnr.value > 0 ? (short)1 : (short)0);
1948 private void iteration(int row_nr,
1962 minit.value = theta.value > (up + Epsb) ? 1 : 0;
1963 if(minit.value != 0)
1966 low.value = low.value == 0 ? 1 : 0;
1968 k = Eta_col_end[Eta_size + 1];
1969 pivot = Eta_value[k - 1];
1970 for(i = Eta_col_end[Eta_size]; i < k; i++)
1972 f = Rhs[Eta_row_nr[i]] - theta.value * Eta_value[i];
1974 Rhs[Eta_row_nr[i]] = f;
1976 if(minit.value == 0)
1978 Rhs[row_nr] = theta.value;
1979 varout = Bas[row_nr];
1980 Bas[row_nr] = varin;
1981 Basis[varout] = FALSE;
1982 Basis[varin] = TRUE;
1983 if(primal != FALSE && pivot < 0)
1984 Lower[varout] = FALSE;
1985 if(low.value == 0 && up < Infinite)
1988 Rhs[row_nr] = up - Rhs[row_nr];
1989 for(i = Eta_col_end[Eta_size]; i < k; i++)
1990 Eta_value[i] = -Eta_value[i];
1995 if(Lp.trace != FALSE)
1997 System.out.print("Theta = " + theta.value + " ");
1998 if(minit.value != 0)
2000 if(Lower[varin] == FALSE)
2001 System.out.print("Iteration:" + Lp.iter +
2002 ", variable" + varin + " changed from 0 to its upper Bound of "
2003 + Upbo[varin] + "\n");
2005 System.out.print("Iteration:" + Lp.iter + ", variable" + varin +
2006 " changed its upper Bound of " + Upbo[varin] + " to 0\n");
2009 System.out.print("Iteration:" + Lp.iter + ", variable" + varin +
2010 " entered basis at:" + Rhs[row_nr] + "\n");
2014 for(i = 1; i <= Rows; i++)
2018 if(Rhs[i] > Upbo[Bas[i]])
2019 f += Rhs[i] - Upbo[Bas[i]];
2020 System.out.println("feasibility gap of this basis:" + (float)f);
2023 System.out.println("objective function value of this feasible basis: " + Rhs[0]);
2028 private int solvelp()
2031 float f = 0, theta = 0;
2033 float[] drow, prow, Pcol;
2039 Ref ref1, ref2, ref3;
2044 drow = new float[Sum + 1];
2045 prow = new float[Sum + 1];
2046 Pcol = new float[Rows + 1];
2047 for (i = 0; i <= Sum; i++) {
2051 for (i = 0; i <= Rows; i++)
2062 while(i != Rows && primal != FALSE)
2065 primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ?
2068 if(Lp.trace != FALSE)
2071 System.out.print("Start at feasible basis\n");
2073 System.out.print("Start at infeasible basis\n");
2078 for(i = 1; i <= Rows; i++)
2081 for(i = 1; i <= columns; i++)
2085 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
2086 if(drow[Mat[j].row_nr] != 0)
2087 drow[varnr] += drow[Mat[j].row_nr] * Mat[j].value;
2088 if(drow[varnr] < Extrad)
2089 Extrad = drow[varnr];
2094 if(Lp.trace != FALSE)
2095 System.out.println("Extrad = " + Extrad);
2098 while(Status == RUNNING)
2103 construct_solution(Solution);
2107 flag = colprim(ref1, minit, drow);
2108 colnr = (int)ref1.value;
2111 setpivcol(Lower[colnr], colnr, Pcol);
2112 ref1.value = row_nr;
2114 flag = rowprim(colnr, ref1, ref2, Pcol);
2115 row_nr = (int)ref1.value;
2118 condensecol(row_nr, Pcol);
2122 else /* not primal */
2124 if(minit == FALSE) {
2125 ref1.value = row_nr;
2126 flag = rowdual(ref1);
2127 row_nr = (int)ref1.value;
2132 flag = coldual(row_nr, ref1, minit, prow, drow);
2133 colnr = (int)ref1.value;
2136 setpivcol(Lower[colnr], colnr, Pcol);
2137 /* getting div by zero here ... MB */
2138 if(Pcol[row_nr] == 0)
2140 System.err.println("An attempt was made to divide by zero (Pcol[" +
2142 System.err.println("This indicates numerical instability");
2144 if(JustInverted == FALSE)
2146 System.out.println("Reinverting Eta");
2151 System.out.println("Can't reinvert, failure");
2157 condensecol(row_nr, Pcol);
2158 f = Rhs[row_nr] - Upbo[Bas[row_nr]];
2161 theta = f / (float) Pcol[row_nr];
2162 if(theta <= Upbo[colnr])
2163 Lower[Bas[row_nr]] =
2164 (Lower[Bas[row_nr]] == FALSE)?
2168 theta = Rhs[row_nr] / (float) Pcol[row_nr];
2172 Status = INFEASIBLE;
2182 if(Doiter != FALSE) {
2185 ref3.value = Lower[colnr];
2186 iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol);
2188 minit = (short)ref2.value;
2189 Lower[colnr] = (short)ref3.value;
2191 if(Num_inv >= Lp.max_num_inv)
2193 if(DoInvert != FALSE)
2199 Lp.total_iter += Lp.iter;
2205 private short is_int(float value)
2209 tmp = (float)(value - Math.floor(value));
2212 if(tmp > (1 - Epsilon))
2217 private void construct_solution(float[] sol)
2222 /* zero all results of rows */
2223 for (i = 0; i <= Rows; i++)
2226 if (Lp.scaling_used != FALSE)
2228 for(i = Rows + 1; i <= Sum; i++)
2229 sol[i] = Lowbo[i] * Lp.scale[i];
2230 for(i = 1; i <= Rows; i++)
2234 sol[basi] += Rhs[i] * Lp.scale[basi];
2236 for(i = Rows + 1; i <= Sum; i++)
2237 if(Basis[i] == FALSE && Lower[i] == FALSE)
2238 sol[i] += Upbo[i] * Lp.scale[i];
2240 for(j = 1; j <= columns; j++)
2244 for(i = Col_end[j - 1]; i < Col_end[j]; i++)
2245 sol[Mat[i].row_nr] += (f / Lp.scale[Rows+j])
2246 * (Mat[i].value / Lp.scale[Mat[i].row_nr]);
2249 for(i = 0; i <= Rows; i++)
2251 if(Math.abs(sol[i]) < Epsb)
2254 if(Lp.ch_sign[i] != FALSE)
2260 for(i = Rows + 1; i <= Sum; i++)
2262 for(i = 1; i <= Rows; i++)
2266 sol[basi] += Rhs[i];
2268 for(i = Rows + 1; i <= Sum; i++)
2269 if(Basis[i] == FALSE && Lower[i] == FALSE)
2271 for(j = 1; j <= columns; j++)
2275 for(i = Col_end[j - 1]; i < Col_end[j]; i++)
2276 sol[Mat[i].row_nr] += f * Mat[i].value;
2279 for(i = 0; i <= Rows; i++)
2281 if(Math.abs(sol[i]) < Epsb)
2284 if(Lp.ch_sign[i] != FALSE)
2288 } /* construct_solution */
2290 private void calculate_duals()
2295 for(i = 1; i <= Rows; i++)
2299 if (Lp.scaling_used != FALSE)
2300 for(i = 1; i <= Rows; i++)
2301 Lp.duals[i] *= Lp.scale[i]/Lp.scale[0];
2303 /* the dual values are the reduced costs of the slacks */
2304 /* When the slack is at its upper Bound, change the sign. Can this happen? */
2305 for(i = 1; i <= Rows; i++)
2307 if(Lp.basis[i] != FALSE)
2309 else if( Lp.ch_sign[0] == Lp.ch_sign[i])
2310 Lp.duals[i] = -Lp.duals[i];
2314 private int milpsolve(float[] upbo,
2320 int i, j, failure, notint, is_worse;
2321 float theta, tmpfloat;
2322 Random rdm = new Random();
2325 if(Break_bb != FALSE)
2329 if(Level > Lp.max_level)
2330 Lp.max_level = Level;
2331 /* make fresh copies of upbo, lowbo, rh as solving changes them */
2332 System.arraycopy(upbo, 0, Upbo, 0, Sum + 1);
2333 System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1);
2334 System.arraycopy(sbasis, 0, Basis, 0, Sum + 1);
2335 System.arraycopy(slower, 0, Lower, 0, Sum + 1);
2336 System.arraycopy(sbas, 0, Bas, 0, Rows + 1);
2337 System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1);
2339 if(Lp.anti_degen != FALSE)
2341 for(i = 1; i <= columns; i++)
2343 tmpfloat = rdm.nextFloat()*(float)0.001;
2345 Lowbo[i + Rows] -= tmpfloat;
2346 tmpfloat = rdm.nextFloat()*(float)0.001;
2348 Upbo[i + Rows] += tmpfloat;
2350 Lp.eta_valid = FALSE;
2353 if(Lp.eta_valid == FALSE)
2355 for(i = 1; i <= columns; i++)
2356 if(Lowbo[Rows + i] != 0)
2358 theta = Lowbo[ Rows + i];
2359 if(Upbo[Rows + i]<Infinite)
2360 Upbo[Rows + i] -= theta;
2361 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
2362 Rh[Mat[j].row_nr] -= theta * Mat[j].value;
2365 Lp.eta_valid = TRUE;
2368 failure = solvelp();
2370 if(Lp.anti_degen != FALSE)
2372 System.arraycopy(upbo, 0, Upbo, 0, Sum + 1);
2373 System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1);
2374 System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1);
2376 for(i = 1; i <= columns; i++)
2377 if(Lowbo[Rows + i] != 0)
2379 theta = Lowbo[ Rows + i];
2380 if(Upbo[Rows + i]<Infinite)
2381 Upbo[Rows + i] -= theta;
2382 for(j = Col_end[i - 1]; j < Col_end[i]; j++)
2383 Rh[Mat[j].row_nr] -= theta * Mat[j].value;
2386 Lp.eta_valid = TRUE;
2387 failure = solvelp();
2390 if(failure == INFEASIBLE && Lp.verbose != FALSE)
2391 System.out.print("level" + Level + " INF\n");
2393 if(failure == OPTIMAL) /* there is a solution */
2395 construct_solution(Solution);
2397 /* if this solution is worse than the best sofar, this branch must die */
2398 if(Maximise != FALSE)
2399 is_worse = (Solution[0] <= Best_solution[0]) ? 1:0;
2400 else /* minimising! */
2401 is_worse = (Solution[0] >= Best_solution[0]) ? 1:0;
2403 if(is_worse != FALSE)
2405 if(Lp.verbose != FALSE)
2406 System.out.println("level" + Level + " OPT NOB value " + Solution[0] +
2407 " Bound " + Best_solution[0]);
2412 /* check if solution contains enough ints */
2413 if(Lp.bb_rule == FIRST_NI)
2417 while(i <= Sum && notint == 0)
2422 if(Lp.bb_rule == RAND_NI)
2424 int nr_not_int, select_not_int;
2426 for(i = Rows + 1; i <= Sum; i++)
2431 select_not_int=(rdm.nextInt() % nr_not_int) + 1;
2433 while(select_not_int > 0)
2441 if(Lp.verbose == TRUE)
2443 System.out.println("level " + Level + " OPT value " + Solution[0]);
2445 System.out.println("level " + Level + " OPT INT value " + Solution[0]);
2447 if(notint != FALSE) /* there is at least one value not yet int */
2449 /* set up two new problems */
2450 float[] new_upbo, new_lowbo;
2452 short[] new_lower, new_basis;
2456 /* allocate room for them */
2457 new_upbo = new float[Sum + 1];
2458 new_lowbo = new float[Sum + 1];
2459 new_lower = new short[Sum + 1];
2460 new_basis = new short[Sum + 1];
2461 new_bas = new int[Rows + 1];
2462 System.arraycopy(upbo, 0, new_upbo, 0, Sum + 1);
2463 System.arraycopy(lowbo, 0, new_lowbo, 0, Sum + 1);
2464 System.arraycopy(Lower, 0, new_lower, 0, Sum + 1);
2465 System.arraycopy(Basis, 0, new_basis, 0, Sum + 1);
2466 System.arraycopy(Bas, 0, new_bas, 0, Rows +1);
2468 if(Floor_first != FALSE)
2470 new_bound = (float)(Math.ceil(Solution[notint]) - 1);
2471 /* this Bound might conflict */
2472 if(new_bound < lowbo[notint])
2476 else /* Bound feasible */
2478 new_upbo[notint] = new_bound;
2479 Lp.eta_valid = FALSE;
2480 resone = milpsolve(new_upbo, lowbo, new_basis, new_lower,
2482 Lp.eta_valid = FALSE;
2485 if(new_bound > upbo[notint])
2489 else /* Bound feasible */
2491 new_lowbo[notint] = new_bound;
2492 Lp.eta_valid = FALSE;
2493 restwo = milpsolve(upbo, new_lowbo, new_basis, new_lower,
2495 Lp.eta_valid = FALSE;
2498 else /* take ceiling first */
2500 new_bound = (float)Math.ceil(Solution[notint]);
2501 /* this Bound might conflict */
2502 if(new_bound > upbo[notint])
2506 else /* Bound feasible */
2508 new_lowbo[notint] = new_bound;
2509 Lp.eta_valid = FALSE;
2510 resone = milpsolve(upbo, new_lowbo, new_basis, new_lower,
2512 Lp.eta_valid = FALSE;
2515 if(new_bound < lowbo[notint])
2519 else /* Bound feasible */
2521 new_upbo[notint] = new_bound;
2522 Lp.eta_valid = FALSE;
2523 restwo = milpsolve(new_upbo, lowbo, new_basis, new_lower,
2525 Lp.eta_valid = FALSE;
2528 if(resone != FALSE && restwo != FALSE)
2529 /* both failed and must have been infeasible */
2530 failure = INFEASIBLE;
2535 else /* all required values are int */
2537 if(Maximise != FALSE)
2538 is_worse = (Solution[0] < Best_solution[0]) ? 1:0;
2540 is_worse = (Solution[0] > Best_solution[0]) ? 1:0;
2542 if(is_worse == FALSE) /* Current solution better */
2544 System.arraycopy(Solution, 0, Best_solution, 0, Sum + 1);
2546 if(Lp.break_at_int != FALSE)
2548 if(Maximise != FALSE && (Best_solution[0] > Lp.break_value))
2550 if(Maximise == FALSE && (Best_solution[0] < Lp.break_value))
2557 /* failure can have the values OPTIMAL, UNBOUNDED and INFEASIBLE. */
2562 public int solve(Problem lp)
2566 if(lp.active == FALSE)
2573 if(Isvalid(lp) != FALSE)
2575 if(Maximise != FALSE && lp.obj_bound == Infinite)
2576 Best_solution[0]=-Infinite;
2577 else if(Maximise == FALSE && lp.obj_bound==-Infinite)
2578 Best_solution[0] = Infinite;
2580 Best_solution[0] = lp.obj_bound;
2584 if(lp.basis_valid == FALSE)
2586 for(i = 0; i <= lp.rows; i++)
2591 for(i = lp.rows+1; i <= lp.sum; i++)
2592 lp.basis[i] = FALSE;
2593 for(i = 0; i <= lp.sum; i++)
2595 lp.basis_valid = TRUE;
2598 lp.eta_valid = FALSE;
2600 result = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas);
2601 lp.eta_size = Eta_size;
2602 lp.eta_alloc = Eta_alloc;
2603 lp.num_inv = Num_inv;
2608 /* if we get here, Isvalid(lp) failed. I suggest we return FAILURE */
2612 public int lag_solve(Problem lp, float start_bound, int num_iter, short verbose)
2614 int i, j, result, citer;
2615 short status, OrigFeas, AnyFeas, same_basis;
2616 float[] OrigObj, ModObj, SubGrad, BestFeasSol;
2617 float Zub, Zlb, Ztmp, pie;
2618 float rhsmod, Step, SqrsumSubGrad;
2623 OrigObj = new float[lp.columns + 1];
2624 ModObj = new float[lp.columns + 1];
2625 for (i = 0; i <= lp.columns; i++)
2627 SubGrad = new float[lp.nr_lagrange];
2628 for (i = 0; i < lp.nr_lagrange; i++)
2630 BestFeasSol = new float[lp.sum + 1];
2631 for (i = 0; i <= lp.sum; i++)
2633 old_bas = new int[lp.rows + 1];
2634 System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows + 1);
2635 old_lower = new short[lp.sum + 1];
2636 System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum + 1);
2638 get_row(lp, 0, OrigObj);
2642 if(lp.maximise != FALSE)
2649 Zlb = -DEF_INFINITE;
2658 for(i = 0 ; i < lp.nr_lagrange; i++)
2661 while(status == RUNNING)
2665 for(i = 1; i <= lp.columns; i++)
2667 ModObj[i] = OrigObj[i];
2668 for(j = 0; j < lp.nr_lagrange; j++)
2670 if(lp.maximise != FALSE)
2671 ModObj[i] -= lp.lambda[j] * lp.lag_row[j][i];
2673 ModObj[i] += lp.lambda[j] * lp.lag_row[j][i];
2676 for(i = 1; i <= lp.columns; i++)
2678 set_mat(lp, 0, i, ModObj[i]);
2679 /* fSystem.out.print(stderr, "%f ", ModObj[i]); */
2682 for(i = 0; i < lp.nr_lagrange; i++)
2683 if(lp.maximise != FALSE)
2684 rhsmod += lp.lambda[i] * lp.lag_rhs[i];
2686 rhsmod -= lp.lambda[i] * lp.lag_rhs[i];
2688 if(verbose != FALSE)
2690 System.out.println("Zub: " + Zub + " Zlb: " + Zlb + " Step: " + Step +
2691 " pie: " + pie + " Feas " + OrigFeas);
2692 for(i = 0; i < lp.nr_lagrange; i++)
2693 System.out.println(i + " SubGrad " + SubGrad[i] + " lambda " + lp.lambda[i]);
2701 while(same_basis != FALSE && i < lp.rows)
2703 same_basis = (old_bas[i] == lp.bas[i])? (short)1: (short)0;
2707 while(same_basis != FALSE && i < lp.sum)
2709 same_basis=(old_lower[i] == lp.lower[i])? (short)1:(short)0;
2712 if(same_basis == FALSE)
2714 System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum+1);
2715 System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows+1);
2719 if(verbose != FALSE)
2720 System.out.println("result: " + result + " same basis: " + same_basis);
2722 if(result == UNBOUNDED)
2724 for(i = 1; i <= lp.columns; i++)
2725 System.out.print(ModObj[i] + " ");
2729 if(result == FAILURE)
2732 if(result == INFEASIBLE)
2733 status = INFEASIBLE;
2736 for(i = 0; i < lp.nr_lagrange; i++)
2738 SubGrad[i]= -lp.lag_rhs[i];
2739 for(j = 1; j <= lp.columns; j++)
2740 SubGrad[i] += lp.best_solution[lp.rows + j] * lp.lag_row[i][j];
2741 SqrsumSubGrad += SubGrad[i] * SubGrad[i];
2745 for(i = 0; i < lp.nr_lagrange; i++)
2746 if(lp.lag_con_type[i] != FALSE)
2748 if(Math.abs(SubGrad[i]) > lp.epsb)
2751 else if(SubGrad[i] > lp.epsb)
2754 if(OrigFeas != FALSE)
2758 for(i = 1; i <= lp.columns; i++)
2759 Ztmp += lp.best_solution[lp.rows + i] * OrigObj[i];
2760 if((lp.maximise != FALSE) && (Ztmp > Zlb))
2763 for(i = 1; i <= lp.sum; i++)
2764 BestFeasSol[i] = lp.best_solution[i];
2765 BestFeasSol[0] = Zlb;
2766 if(verbose != FALSE)
2767 System.out.println("Best feasible solution: " + Zlb);
2772 for(i = 1; i <= lp.sum; i++)
2773 BestFeasSol[i] = lp.best_solution[i];
2774 BestFeasSol[0] = Zub;
2775 if(verbose != FALSE)
2776 System.out.println("Best feasible solution: " + Zub);
2780 if(lp.maximise != FALSE)
2781 Zub = Math.min(Zub, rhsmod + lp.best_solution[0]);
2783 Zlb = Math.max(Zlb, rhsmod + lp.best_solution[0]);
2785 if(Math.abs(Zub-Zlb) < (float)0.001)
2789 Step = (float)(pie * ((1.05*Zub) - Zlb) / SqrsumSubGrad);
2791 for(i = 0; i < lp.nr_lagrange; i++)
2793 lp.lambda[i] += Step * SubGrad[i];
2794 if(lp.lag_con_type[i] == FALSE && lp.lambda[i] < 0)
2798 if(citer == num_iter && status==RUNNING)
2799 if(AnyFeas != FALSE)
2800 status = FEAS_FOUND;
2802 status = NO_FEAS_FOUND;
2805 for(i = 0; i <= lp.sum; i++)
2806 lp.best_solution[i] = BestFeasSol[i];
2808 for(i = 1; i <= lp.columns; i++)
2809 set_mat(lp, 0, i, OrigObj[i]);
2811 if(lp.maximise != FALSE)
2818 } // end of class solve
2821 public static class Matrix
2825 public Matrix(int r, float v) {
2832 public static class Problem {
2833 String lp_name; /* the name of the lp */
2835 public short active; /*TRUE if the globals point to this structure*/
2836 public short verbose; /* ## Verbose flag */
2837 public short debug; /* ## Print B&B information */
2838 public short trace; /* ## Print information on pivot selection */
2839 public short anti_degen; /* ## Do perturbations */
2841 public int rows; /* Nr of constraint rows in the problem */
2842 int rows_alloc; /* The allocated memory for Rows sized data */
2843 int columns; /* The number of columns (= variables) */
2845 int sum; /* The size of the variables + the slacks */
2848 short names_used; /* Flag to indecate if names for rows and
2850 String[] row_name; /* rows_alloc+1 */
2851 int[] col_name; /* columns_alloc+1 */
2853 /* Row[0] of the sparce matrix is the objective function */
2855 int non_zeros; /* The number of elements in the sparce matrix*/
2856 int mat_alloc; /* The allocated size for matrix sized
2858 Matrix[] mat; /* mat_alloc :The sparse matrix */
2859 int[] col_end; /* columns_alloc+1 :Cend[i] is the index of the
2860 first element after column i.
2861 column[i] is stored in elements
2862 col_end[i-1] to col_end[i]-1 */
2863 int[] col_no; /* mat_alloc :From Row 1 on, col_no contains the
2865 nonzero elements, row by row */
2866 short row_end_valid; /* true if row_end & col_no are valid */
2867 int[] row_end; /* rows_alloc+1 :row_end[i] is the index of the
2868 first element in Colno after row i */
2869 float[] orig_rh; /* rows_alloc+1 :The RHS after scaling & sign
2870 changing, but before `Bound transformation' */
2871 float[] rh; /* rows_alloc+1 :As orig_rh, but after Bound
2873 float[] rhs; /* rows_alloc+1 :The RHS of the curent simplex
2875 float[] orig_upbo; /* sum_alloc+1 :Bound before transformations */
2876 float[] orig_lowbo; /* " " */
2877 float[] upbo; /* " " :Upper bound after transformation
2879 float[] lowbo; /* " " :Lower bound after transformation
2882 short basis_valid; /* TRUE is the basis is still valid */
2883 int[] bas; /* rows_alloc+1 :The basis column list */
2884 short[] basis; /* sum_alloc+1 : basis[i] is TRUE if the column
2886 short[] lower; /* " " :TRUE is the variable is at its
2887 lower bound (or in the basis), it is FALSE
2888 if the variable is at its upper bound */
2890 short eta_valid; /* TRUE if current Eta structures are valid */
2891 int eta_alloc; /* The allocated memory for Eta */
2892 int eta_size; /* The number of Eta columns */
2893 int num_inv; /* The number of float pivots */
2894 int max_num_inv; /* ## The number of float pivots between
2896 float[] eta_value; /* eta_alloc :The Structure containing the
2898 int[] eta_row_nr; /* " " :The Structure containing the Row
2900 int[] eta_col_end; /* rows_alloc + MaxNumInv : eta_col_end[i] is
2901 the start index of the next Eta column */
2903 short bb_rule; /* what rule for selecting B&B variables */
2905 short break_at_int; /* TRUE if stop at first integer better than
2909 float obj_bound; /* ## Objective function bound for speedup of
2911 int iter; /* The number of iterations in the simplex
2913 int total_iter; /* The total number of iterations (B&B) (ILP)*/
2914 int max_level; /* The Deepest B&B level of the last solution */
2915 int total_nodes; /* total number of nodes processed in b&b */
2916 public float[] solution; /* sum_alloc+1 :The Solution of the last LP,
2917 0 = The Optimal Value,
2919 rows+1..sum The Variables */
2920 public float[] best_solution; /* " " :The Best 'Integer' Solution */
2921 float[] duals; /* rows_alloc+1 :The dual variables of the
2924 short maximise; /* TRUE if the goal is to maximise the
2925 objective function */
2926 short floor_first; /* TRUE if B&B does floor bound first */
2927 short[] ch_sign; /* rows_alloc+1 :TRUE if the Row in the matrix
2929 (a`x > b, x>=0) is translated to
2930 s + -a`x = -b with x>=0, s>=0) */
2932 short scaling_used; /* TRUE if scaling is used */
2933 short columns_scaled; /* TRUE is the columns are scaled too, Only use
2934 if all variables are non-integer */
2935 float[] scale; /* sum_alloc+1 :0..Rows the scaling of the Rows,
2936 Rows+1..Sum the scaling of the columns */
2938 int nr_lagrange; /* Nr. of Langrangian relaxation constraints */
2939 float[][]lag_row; /* NumLagrange, columns+1:Pointer to pointer of
2941 float[] lag_rhs; /* NumLagrange :Pointer to pointer of Rhs */
2942 float[] lambda; /* NumLagrange :Lambda Values */
2943 short[] lag_con_type; /* NumLagrange :TRUE if constraint type EQ */
2944 float lag_bound; /* the lagrangian lower bound */
2946 short valid; /* Has this lp pased the 'test' */
2947 float infinite; /* ## numercal stuff */
2948 float epsilon; /* ## */
2949 float epsb; /* ## */
2950 float epsd; /* ## */
2951 float epsel; /* ## */
2954 public Problem (int nrows, int ncolumns)
2958 nsum=nrows+ncolumns;
2959 if(rows < 0 || columns < 0)
2960 System.err.print("rows < 0 or columns < 0");
2962 lp_name = new String("unnamed");
2972 columns_alloc=columns;
2976 obj_bound=DEF_INFINITE;
2977 infinite=DEF_INFINITE;
2978 epsilon=DEF_EPSILON;
2985 mat = new Matrix[mat_alloc];
2986 for (i = 0; i < mat_alloc; i++)
2987 mat[i] = new Matrix(0, 0);
2989 col_no = new int[mat_alloc];
2990 for (i = 0; i < mat_alloc; i++)
2993 col_end = new int[columns + 1];
2994 for (i = 0; i < columns + 1; i++)
2997 row_end = new int[rows + 1];
2998 for (i = 0; i < rows + 1; i++)
3001 row_end_valid=FALSE;
3003 orig_rh = new float[rows + 1];
3004 for (i = 0; i < rows + 1; i++)
3007 rh = new float[rows + 1];
3008 for (i = 0; i < rows + 1; i++)
3011 rhs = new float[rows + 1];
3012 for (i = 0; i < rows + 1; i++)
3015 orig_upbo = new float[sum + 1];
3016 for(i = 0; i <= sum; i++)
3017 orig_upbo[i]=infinite;
3019 upbo = new float[sum + 1];
3020 for (i = 0; i < sum + 1; i++)
3023 orig_lowbo = new float[sum + 1];
3024 for (i = 0; i < sum + 1; i++)
3027 lowbo = new float[sum + 1];
3028 for (i = 0; i < sum + 1; i++)
3033 bas = new int[rows+1];
3034 for (i = 0; i <= rows; i++)
3037 basis = new short[sum + 1];
3038 for (i = 0; i <= sum; i++)
3041 lower = new short[sum + 1];
3042 for(i = 0; i <= rows; i++)
3047 for(i = rows + 1; i <= sum; i++)
3049 for(i = 0 ; i <= sum; i++)
3055 max_num_inv=DEFNUMINV;
3059 eta_value = new float[eta_alloc];
3060 for (i = 0; i < eta_alloc; i++)
3063 eta_row_nr = new int[eta_alloc];
3064 for (i = 0; i < eta_alloc; i++)
3067 eta_col_end = new int[rows_alloc + max_num_inv];
3068 for (i = 0; i < rows_alloc + max_num_inv; i++)
3078 solution = new float[sum + 1];
3079 for (i = 0; i <= sum; i++)
3082 best_solution = new float[sum + 1];
3083 for (i = 0; i <= sum; i++)
3084 best_solution[i] = 0;
3086 duals = new float[rows + 1];
3087 for (i = 0; i <= rows; i++)
3093 scaling_used = FALSE;
3094 columns_scaled = FALSE;
3096 ch_sign = new short[rows + 1];
3097 for(i = 0; i <= rows; i++)
3101 } // end of constructor from row and column
3103 //***************************************
3104 // return the ith member of the best_solution[]
3106 public float getBestSolution(int i) {
3107 return best_solution[i];
3110 //***************************************
3111 // get the number of rows
3113 public int getRows() {
3117 //***************************************
3118 // get the number of columns
3120 public int getcolumns() {
3124 } // end of class Problem
3126 private final static float round( float val, float eps) { return (Math.abs(val) < eps) ? 0 : val; }