From: adam Date: Mon, 5 Apr 2004 00:13:34 +0000 (+0000) Subject: another performance improvement for the simplex solver X-Git-Tag: RC4~47 X-Git-Url: http://git.megacz.com/?p=org.ibex.core.git;a=commitdiff_plain;h=ae0cd6164facf843cc136f52e3f78e3843191d61 another performance improvement for the simplex solver darcs-hash:20040405001334-5007d-3aea508c4f0a878e775f8e59c4275f187c22beb9.gz --- diff --git a/Makefile.upstream b/Makefile.upstream index 48f5cb5..b4e8759 100644 --- a/Makefile.upstream +++ b/Makefile.upstream @@ -17,7 +17,7 @@ url_jikes-1.18 := http://dist.xwt.org/jikes-1.18.tgz url_libmspack-20030726 := http://www.kyz.uklinux.net/downloads/libmspack-20030726.tar.gz url_vera-1.10 := http://fgo-temp.acc.umu.se/pub/GNOME/sources/ttf-bitstream-vera/1.10/ttf-bitstream-vera-1.10.tar.gz url_WindowMaker-0.80.2 := http://windowmaker.org/pub/source/release/WindowMaker-0.80.2.tar.gz -url_bcel-5.1 := http://apache.webmeta.com/jakarta/bcel/binaries/bcel-5.1.tar.gz +url_bcel-5.1 := http://www.apache.org/dist/jakarta/bcel/binaries/bcel-5.1.tar.gz .install_binutils-2.13.2.1_powerpc-apple-darwin: .vendor rm -rf upstream/darwin-linker/src @@ -107,6 +107,18 @@ environment_gcc_3.3_$(target) += PATH=$(shell pwd)/upstream/install/bi .PRECIOUS: .vendor .download_% .configure_%_$(target) .install_%_$(target) +.download_nestedvm: + mkdir -p upstream/nestedvm + cd upstream/nestedvm; wget -nH -r http://nestedvm.ibex.org/ + touch $@ + +.build_nestedvm: .download_nestedvm + cd upstream/nestedvm; make + touch $@ + +.install_nestedvm: .build_nestedvm + touch $@ + # vendor-supplied binaries and headers; this is stuff that comes with various OSes vendor: .vendor; @true .vendor: diff --git a/doc/reference/reference.xml b/doc/reference/reference.xml index 844bad6..ac79dec 100644 --- a/doc/reference/reference.xml +++ b/doc/reference/reference.xml @@ -402,8 +402,6 @@
-
- Each box occupies a rectangular region on the surface. The visual appearance of a surface is created by rendering each box in its tree. Unless the [[clip]] attribute is [[false]], each box will @@ -459,9 +457,6 @@ These eight components plus the size of a box fully specify its appearance. Every single box you see in Ibex is drawn only on the basis of these components and its size. -
- -
The size and position of every box is determined by its properties, its childrens' sizes, and its parent's size and @@ -495,8 +490,6 @@ resizing the surface at all. However, not all platforms give Ibex enough control to do this. - - When talking about positioning, we will often refer to the @@ -547,9 +540,6 @@ - - -
@@ -582,7 +572,14 @@
-
+ + + Each non-packed box is transformed according to the parent's + [[transform]] property and then positioned so that its alignment + point is [[(child.x, child.y)]] pixels from the corresponding + edge/corner/center of its parent. + + Ibex formulates a set of constraints for placing a box's **packed** children as follows: @@ -615,14 +612,6 @@ does not (due to [[maxwidth]] or minimum width constraints), the box's will be placed so that its alignment point coincides with the alignment point of that rectangle of cells. -
- -
- Each non-packed box is transformed according to the parent's - [[transform]] property and then positioned so that its alignment - point is [[(child.x, child.y)]] pixels from the corresponding - edge/corner/center of its parent. -
@@ -663,8 +652,7 @@
- - +
diff --git a/src/org/ibex/Box.java b/src/org/ibex/Box.java index fd4e578..2ada11f 100644 --- a/src/org/ibex/Box.java +++ b/src/org/ibex/Box.java @@ -294,46 +294,43 @@ public final class Box extends JSScope implements Scheduler.Task { } private static float[] coeff = null; - private static LinearProgramming.Simplex lp_h = new LinearProgramming.Simplex(); - private static LinearProgramming.Problem lpr_h = new LinearProgramming.Problem(50, 50, 300); - private static LinearProgramming.Simplex lp_v = new LinearProgramming.Simplex(); - private static LinearProgramming.Problem lpr_v = new LinearProgramming.Problem(50, 50, 300); + private static LinearProgramming.Simplex lp_h = new LinearProgramming.Simplex(50, 50, 300); + private static LinearProgramming.Simplex lp_v = new LinearProgramming.Simplex(50, 50, 300); void place_children() { int numkids = 0; for(Box c = firstPackedChild(); c != null; c = c.nextPackedSibling()) numkids++; //#repeat col/row colspan/rowspan contentwidth/contentheight width/height \ - // maxwidth/maxheight cols/rows minwidth/minheight lp_h/lp_v lpr_h/lpr_v + // maxwidth/maxheight cols/rows minwidth/minheight lp_h/lp_v lp_h/lp_v do { int nc = numkids * 2 + cols * 3 + 1 + 2; if (coeff == null || nc+1>coeff.length) coeff = new float[nc+1]; - lpr_h.init(nc, nc); + lp_h.init(nc, nc); // objective function for(int i=0; i=child.col && i=child.col && i b, x>=0) is translated to + s + -a`x = -b with x>=0, s>=0) */ + + int nr_lagrange; /* Nr. of Langrangian relaxation constraints */ + float[][]lag_row; /* NumLagrange, columns+1:Pointer to pointer of + rows */ + float[] lag_rhs; /* NumLagrange :Pointer to pointer of Rhs */ + float[] lambda; /* NumLagrange :Lambda Values */ + short[] lag_con_type; /* NumLagrange :TRUE if constraint type EQ */ + float lag_bound; /* the lagrangian lower bound */ + + short valid; /* Has this lp pased the 'test' */ + float infinite; /* ## numercal stuff */ + float epsilon; /* ## */ + float epsb; /* ## */ + float epsd; /* ## */ + float epsel; /* ## */ + int Rows; int columns; int Sum; int Non_zeros; int Level; - Matrix[] Mat; + MatrixArray Mat; int[] Col_no; int[] Col_end; int[] Row_end; @@ -98,1772 +202,889 @@ public class LinearProgramming { int Warn_count; /* used in CHECK version of rounding macro */ - public void set_mat(Problem lp, int Row, int column, float Value) - { - int elmnr, lastelm, i; - // System.err.println("lp.mat.length = " + lp.mat.length); - - if(Row > lp.rows || Row < 0) - System.err.print("Row out of range"); - if(column > lp.columns || column < 1) - System.err.print("column out of range"); - if(lp.scaling_used != FALSE) - Value *= lp.scale[Row] * lp.scale[lp.rows + column]; - - if(true /*abs(Value) > lp.epsilon*/) - { - if (lp.basis[column] == TRUE && Row > 0) - lp.basis_valid = FALSE; - lp.eta_valid = FALSE; - elmnr = lp.col_end[column-1]; - while((elmnr < lp.col_end[column]) ? - (lp.mat[elmnr].row_nr != Row) : false) - elmnr++; - - if((elmnr != lp.col_end[column]) ? - (lp.mat[elmnr].row_nr == Row) : false ) - if (lp.scaling_used != FALSE) - { - if (lp.ch_sign[Row] != FALSE) - lp.mat[elmnr].value = - -Value * lp.scale[Row] * lp.scale[column]; - else - lp.mat[elmnr].value = - Value * lp.scale[Row] * lp.scale[column]; - } - else - { - if (lp.ch_sign[Row] != FALSE) - lp.mat[elmnr].value = -Value; - else - lp.mat[elmnr].value = Value; - } - else - { - /* check if more space is needed for matrix */ - if (lp.non_zeros + 1 > lp.mat_alloc) throw new Error("not enough mat space; this should not happen"); - - /* Shift the matrix */ - lastelm=lp.non_zeros; - for(i = lastelm; i > elmnr ; i--) - lp.mat[i]=lp.mat[i-1]; - for(i = column; i <= lp.columns; i++) - lp.col_end[i]++; - - /* Set new element */ - lp.mat[elmnr].row_nr=Row; - - if (lp.scaling_used != FALSE) - { - if (lp.ch_sign[Row] != FALSE) - lp.mat[elmnr].value=-Value*lp.scale[Row]*lp.scale[column]; - else - lp.mat[elmnr].value=Value*lp.scale[Row]*lp.scale[column]; - } - else - { - if (lp.ch_sign[Row] != FALSE) - lp.mat[elmnr].value=-Value; - else - lp.mat[elmnr].value=Value; - } - - lp.row_end_valid=FALSE; - - lp.non_zeros++; - if (lp.active != FALSE) - Non_zeros=lp.non_zeros; - } - } + public Simplex (int nrows, int ncolumns, int matalloc) { + int nsum; + nsum=nrows+ncolumns; + rows=nrows; + columns=ncolumns; + sum=nsum; + rows_alloc=rows; + columns_alloc=columns; + sum_alloc=sum; + mat_alloc=matalloc; + eta_alloc=10000; + max_num_inv=DEFNUMINV; + col_no = new int[mat_alloc]; + col_end = new int[columns + 1]; + row_end = new int[rows + 1]; + orig_rh = new float[rows + 1]; + rh = new float[rows + 1]; + rhs = new float[rows + 1]; + orig_upbo = new float[sum + 1]; + upbo = new float[sum + 1]; + orig_lowbo = new float[sum + 1]; + lowbo = new float[sum + 1]; + bas = new int[rows+1]; + basis = new short[sum + 1]; + lower = new short[sum + 1]; + eta_value = new float[eta_alloc]; + eta_row_nr = new int[eta_alloc]; + eta_col_end = new int[rows_alloc + max_num_inv]; + solution = new float[sum + 1]; + best_solution = new float[sum + 1]; + duals = new float[rows + 1]; + ch_sign = new short[rows + 1]; + mat = new MatrixArray(mat_alloc); + alternate_mat = new MatrixArray(mat_alloc); } + + public void init(int nrows, int ncolumns) { + int nsum; + nsum=nrows+ncolumns; + active=FALSE; + debug=FALSE; + trace=FALSE; + rows=nrows; + columns=ncolumns; + sum=nsum; + obj_bound=DEF_INFINITE; + infinite=DEF_INFINITE; + epsilon=DEF_EPSILON; + epsb=DEF_EPSB; + epsd=DEF_EPSD; + epsel=DEF_EPSEL; + non_zeros=0; - public void set_obj_fn(Problem lp, float[] row) - { - for(int i=row.length-1; i>0; i--) row[i] = row[i-1]; - row[0] = (float)0.0; + for(int i = 0; i < mat_alloc; i++) { set_row_nr(mat,i, 0); set_value(mat, i, 0); } + for(int i = 0; i < mat_alloc; i++) col_no[i] = 0; + for(int i = 0; i < columns + 1; i++) col_end[i] = 0; + for(int i = 0; i < rows + 1; i++) row_end[i] = 0; + for(int i = 0; i < rows + 1; i++) orig_rh[i] = 0; + for(int i = 0; i < rows + 1; i++) rh[i] = 0; + for(int i = 0; i < rows + 1; i++) rhs[i] = 0; + for(int i = 0; i <= sum; i++) orig_upbo[i]=infinite; + for(int i = 0; i < sum + 1; i++) upbo[i] = 0; + for(int i = 0; i < sum + 1; i++) orig_lowbo[i] = 0; + for(int i = 0; i < sum + 1; i++) lowbo[i] = 0; + for(int i = 0; i <= rows; i++) bas[i] = 0; + for(int i = 0; i <= sum; i++) basis[i] = 0; + for(int i = 0; i <= rows; i++) { bas[i]=i; basis[i]=TRUE; } + for(int i = rows + 1; i <= sum; i++) basis[i]=FALSE; + for(int i = 0 ; i <= sum; i++) lower[i]=TRUE; + for(int i = 0; i < eta_alloc; i++) eta_value[i] = 0; + for(int i = 0; i < eta_alloc; i++) eta_row_nr[i] = 0; + for(int i = 0; i < rows_alloc + max_num_inv; i++) eta_col_end[i] = 0; + for(int i = 0; i <= sum; i++) solution[i] = 0; + for(int i = 0; i <= sum; i++) best_solution[i] = 0; + for(int i = 0; i <= rows; i++) duals[i] = 0; + for(int i = 0; i <= rows; i++) ch_sign[i] = FALSE; - int i; - for(i = 1; i <= lp.columns; i++) - set_mat(lp, 0, i, row[i]); + row_end_valid=FALSE; + bb_rule=FIRST_NI; + break_at_int=FALSE; + break_value=0; + iter=0; + total_iter=0; + basis_valid=TRUE; + eta_valid=TRUE; + eta_size=0; + nr_lagrange=0; + maximise = FALSE; + floor_first = TRUE; + valid = FALSE; } + public void setObjective(float[] row, boolean maximize) { + for(int i=row.length-1; i>0; i--) row[i] = row[i-1]; + row[0] = (float)0.0; + for(int j = 1; j <= columns; j++) { + int Row = 0; + int column = j; + float Value = row[j]; + int elmnr, lastelm; + + if(Row > rows || Row < 0) throw new Error("row out of range"); + if(column > columns || column < 1) throw new Error("column out of range"); + + if (basis[column] == TRUE && Row > 0) basis_valid = FALSE; + eta_valid = FALSE; + elmnr = col_end[column-1]; + while((elmnr < col_end[column]) ? (get_row_nr(mat, elmnr) != Row) : false) elmnr++; + if((elmnr != col_end[column]) ? (get_row_nr(mat, elmnr) == Row) : false ) { + if (ch_sign[Row] != FALSE) set_value(mat, elmnr, -Value); + else set_value(mat, elmnr, Value); + } else { + /* check if more space is needed for matrix */ + if (non_zeros + 1 > mat_alloc) throw new Error("not enough mat space; this should not happen"); + /* Shift the matrix */ + lastelm=non_zeros; + for(int i = lastelm; i > elmnr ; i--) { + set_row_nr(mat,i,get_row_nr(mat,i-1)); + set_value(mat,i,get_value(mat,i-1)); + } + for(int i = column; i <= columns; i++) col_end[i]++; + /* Set new element */ + set_row_nr(mat,elmnr, Row); + if (ch_sign[Row] != FALSE) set_value(mat, elmnr, -Value); + else set_value(mat, elmnr, Value); + row_end_valid=FALSE; + non_zeros++; + if (active != FALSE) Non_zeros=non_zeros; + } + } + if (maximize) { + if (maximise == FALSE) { + for(int i = 0; i < non_zeros; i++) + if(get_row_nr(mat, i)==0) + set_value(mat, i, get_value(mat,i)* (float)-1.0); + eta_valid=FALSE; + } + maximise=TRUE; + ch_sign[0]=TRUE; + if (active != FALSE) Maximise=TRUE; + } else { + if (maximise==TRUE) { + for(int i = 0; i < non_zeros; i++) + if(get_row_nr(mat, i)==0) + set_value(mat, i, get_value(mat,i) * (float)-1.0); + eta_valid=FALSE; + } + maximise=FALSE; + ch_sign[0]=FALSE; + if (active != FALSE) Maximise=FALSE; + } + } - public void add_constraint(Problem lp, float[] row, short constr_type, float rh) - { + public void add_constraint(float[] row, short constr_type, float rh) { for(int i=row.length-1; i>0; i--) row[i] = row[i-1]; row[0] = (float)0.0; - Matrix[] newmat; - int i, j; + MatrixArray newmat; int elmnr; int stcol; - short[] addtoo; - - addtoo = new short[lp.columns + 1]; - - for(i = 1; i <= lp.columns; i++) - if(row[i]!=0) - { - addtoo[i]=TRUE; - lp.non_zeros++; - } - else - addtoo[i]=FALSE; - - //newmat = new Matrix[lp.non_zeros]; - - // FIXME - newmat = lp.alternate_mat; - for (i = 0; i < newmat.length; i++) { newmat[i].row_nr = 0; newmat[i].value = 0; } - - if (lp.non_zeros > lp.mat_alloc) throw new Error("not enough mat space; this should not happen"); - lp.rows++; - lp.sum++; - if(lp.rows > lp.rows_alloc) - throw new Error("not enough rows; ("+lp.rows+" needed, have "+lp.rows_alloc+") this should never happen"); - /* - if (lp.scaling_used != FALSE) - { - // shift scale - for(i=lp.sum; i > lp.rows; i--) - lp.scale[i]=lp.scale[i-1]; - lp.scale[lp.rows]=1; - } - - if (lp.names_used != FALSE) - lp.row_name[lp.rows] = new String("r_" + lp.rows); - if(lp.scaling_used != FALSE && lp.columns_scaled != FALSE) - for(i = 1; i <= lp.columns; i++) - row[i] *= lp.scale[lp.rows+i]; - */ - if(constr_type==GE) - lp.ch_sign[lp.rows] = TRUE; - else - lp.ch_sign[lp.rows] = FALSE; + newmat = alternate_mat; + for(int i = 0; i < non_zeros; i++) { set_row_nr(newmat,i, 0); set_value(newmat, i, 0); } + for(int i = 1; i <= columns; i++) if (row[i]!=0) non_zeros++; + if (non_zeros > mat_alloc) throw new Error("not enough mat space; this should not happen"); + rows++; + sum++; + if(rows > rows_alloc) throw new Error("not enough rows; this should never happen"); + if(constr_type==GE) ch_sign[rows] = TRUE; + else ch_sign[rows] = FALSE; elmnr = 0; stcol = 0; - for(i = 1; i <= lp.columns; i++) - { - for(j = stcol; j < lp.col_end[i]; j++) - { - newmat[elmnr].row_nr=lp.mat[j].row_nr; - newmat[elmnr].value=lp.mat[j].value; - elmnr++; - } - if(addtoo[i] != FALSE) - { - if(lp.ch_sign[lp.rows] != FALSE) - newmat[elmnr].value = -row[i]; - else - newmat[elmnr].value = row[i]; - newmat[elmnr].row_nr = lp.rows; - elmnr++; - } - stcol=lp.col_end[i]; - lp.col_end[i]=elmnr; - } - - lp.alternate_mat = lp.mat; - lp.mat = newmat; - - for(i=lp.sum ; i > lp.rows; i--) - { - lp.orig_upbo[i]=lp.orig_upbo[i-1]; - lp.orig_lowbo[i]=lp.orig_lowbo[i-1]; - lp.basis[i]=lp.basis[i-1]; - lp.lower[i]=lp.lower[i-1]; - } - - for(i= 1 ; i <= lp.rows; i++) - if(lp.bas[i] >= lp.rows) - lp.bas[i]++; + for(int i = 1; i <= columns; i++) { + for(int j = stcol; j < col_end[i]; j++) { + set_row_nr(newmat,elmnr, get_row_nr(mat, j)); + set_value(newmat, elmnr, get_value(mat,j)); + elmnr++; + } + if(((i>=1 && i< columns && row[i]!=0)?TRUE:FALSE) != FALSE) { + if(ch_sign[rows] != FALSE) set_value(newmat, elmnr, -row[i]); + else set_value(newmat, elmnr, row[i]); + set_row_nr(newmat,elmnr, rows); + elmnr++; + } + stcol=col_end[i]; + col_end[i]=elmnr; + } + + alternate_mat = mat; + mat = newmat; + + for(int i = sum ; i > rows; i--) { + orig_upbo[i]=orig_upbo[i-1]; + orig_lowbo[i]=orig_lowbo[i-1]; + basis[i]=basis[i-1]; + lower[i]=lower[i-1]; + } - if(constr_type==LE || constr_type==GE) - { - lp.orig_upbo[lp.rows]=lp.infinite; - } - else if(constr_type==EQ) - { - lp.orig_upbo[lp.rows]=0; - } - else - { - System.err.print("Wrong constraint type\n"); - System.exit(FAIL); - } + for(int i = 1 ; i <= rows; i++) if(bas[i] >= rows) bas[i]++; - lp.orig_lowbo[lp.rows]=0; + if(constr_type==LE || constr_type==GE) orig_upbo[rows]=infinite; + else if(constr_type==EQ) orig_upbo[rows]=0; + else throw new Error("Wrong constraint type\n"); + orig_lowbo[rows]=0; - if(constr_type==GE && rh != 0) - lp.orig_rh[lp.rows]=-rh; - else - lp.orig_rh[lp.rows]=rh; + if(constr_type==GE && rh != 0) orig_rh[rows]=-rh; + else orig_rh[rows]=rh; - lp.row_end_valid=FALSE; - - lp.bas[lp.rows]=lp.rows; - lp.basis[lp.rows]=TRUE; - lp.lower[lp.rows]=TRUE; + row_end_valid=FALSE; - if (lp.active != FALSE) - set_globals(lp); - lp.eta_valid=FALSE; - } - - public void del_constraint(Problem lp, int del_row) - { - int i, j; - int elmnr; - int startcol; - - if(del_row < 1 || del_row > lp.rows) - { - System.err.println("There is no constraint nr. " + del_row); - System.exit(FAIL); - } - - elmnr=0; - startcol=0; - - for(i = 1; i <= lp.columns; i++) - { - for(j=startcol; j < lp.col_end[i]; j++) - { - if(lp.mat[j].row_nr!=del_row) - { - lp.mat[elmnr]=lp.mat[j]; - if(lp.mat[elmnr].row_nr > del_row) - lp.mat[elmnr].row_nr--; - elmnr++; - } - else - lp.non_zeros--; - } - startcol=lp.col_end[i]; - lp.col_end[i]=elmnr; - } - for(i = del_row; i < lp.rows; i++) - { - lp.orig_rh[i] = lp.orig_rh[i+1]; - lp.ch_sign[i] = lp.ch_sign[i+1]; - lp.bas[i] = lp.bas[i+1]; - if (lp.names_used != FALSE) - lp.row_name[i] = lp.row_name[i+1]; - } - for(i = 1; i < lp.rows; i++) - if(lp.bas[i] > del_row) - lp.bas[i]--; - - for(i=del_row; i < lp.sum; i++) - { - lp.lower[i]=lp.lower[i+1]; - lp.basis[i]=lp.basis[i+1]; - lp.orig_upbo[i]=lp.orig_upbo[i+1]; - lp.orig_lowbo[i]=lp.orig_lowbo[i+1]; - if (lp.scaling_used != FALSE) - lp.scale[i]=lp.scale[i+1]; - } - - lp.rows--; - lp.sum--; - - lp.row_end_valid=FALSE; - - if (lp.active != FALSE) - set_globals(lp); - lp.eta_valid=FALSE; - lp.basis_valid=FALSE; - } - - - public void add_column(Problem lp, float[] column) - { - int i, elmnr; - - lp.columns++; - lp.sum++; - if (lp.columns > lp.columns_alloc) throw new Error("not enough cols; this should never happen"); - if (lp.non_zeros + lp.rows+1 > lp.mat_alloc) throw new Error("not enough mat space; this should not happen"); - - if (lp.scaling_used != FALSE) - { - for(i = 0; i <= lp.rows; i++) - column[i]*=lp.scale[i]; - lp.scale[lp.sum]=1; - } - - elmnr=lp.col_end[lp.columns-1]; - for(i = 0 ; i <= lp.rows ; i++) - if(column[i] != 0) - { - lp.mat[elmnr].row_nr=i; - if(lp.ch_sign[i] != FALSE) - lp.mat[elmnr].value=-column[i]; - else - lp.mat[elmnr].value=column[i]; - lp.non_zeros++; - elmnr++; - } - lp.col_end[lp.columns]=elmnr; - lp.orig_lowbo[lp.sum]=0; - lp.orig_upbo[lp.sum]=lp.infinite; - lp.lower[lp.sum]=TRUE; - lp.basis[lp.sum]=FALSE; - if (lp.names_used != FALSE) - lp.col_name[lp.columns] = 0; - + bas[rows]=rows; + basis[rows]=TRUE; + lower[rows]=TRUE; - lp.row_end_valid=FALSE; - - if (lp.active != FALSE) - { - Sum=lp.sum; - columns=lp.columns; - Non_zeros=lp.non_zeros; - } - } - - public void del_column(Problem lp, int column) - { - int i, j, from_elm, to_elm, elm_in_col; - if(column > lp.columns || column < 1) - System.err.print("column out of range in del_column"); - for(i = 1; i <= lp.rows; i++) - { - if(lp.bas[i]==lp.rows+column) - lp.basis_valid=FALSE; - else if(lp.bas[i] > lp.rows+column) - lp.bas[i]--; - } - for(i = lp.rows+column; i < lp.sum; i++) - { - if (lp.names_used != FALSE) - lp.col_name[i-lp.rows] = lp.col_name[i-lp.rows+1]; - lp.orig_upbo[i]=lp.orig_upbo[i+1]; - lp.orig_lowbo[i]=lp.orig_lowbo[i+1]; - lp.upbo[i]=lp.upbo[i+1]; - lp.lowbo[i]=lp.lowbo[i+1]; - lp.basis[i]=lp.basis[i+1]; - lp.lower[i]=lp.lower[i+1]; - if (lp.scaling_used != FALSE) - lp.scale[i]=lp.scale[i+1]; - } - for(i = 0; i < lp.nr_lagrange; i++) - for(j = column; j <= lp.columns; j++) - lp.lag_row[i][j]=lp.lag_row[i][j+1]; - to_elm=lp.col_end[column-1]; - from_elm=lp.col_end[column]; - elm_in_col=from_elm-to_elm; - for(i = from_elm; i < lp.non_zeros; i++) - { - lp.mat[to_elm]=lp.mat[i]; - to_elm++; - } - for(i = column; i < lp.columns; i++) - lp.col_end[i]=lp.col_end[i+1]-elm_in_col; - lp.non_zeros -= elm_in_col; - lp.row_end_valid=FALSE; - lp.eta_valid=FALSE; - - lp.sum--; - lp.columns--; - if (lp.active != FALSE) - set_globals(lp); + if (active != FALSE) set_globals(); + eta_valid=FALSE; } - public void bound_sum(Problem lp, int column1, int column2, float bound, short type, float[] scratch) { + public void bound_sum(int column1, int column2, float bound, short type, float[] scratch) { for(int i=0; i lp.columns || column < 1) - System.err.print("column out of range"); - if (lp.scaling_used != FALSE) - value /= lp.scale[lp.rows + column]; - if(value < lp.orig_lowbo[lp.rows + column]) - System.err.print("UpperBound must be >= lowerBound"); - lp.eta_valid = FALSE; - lp.orig_upbo[lp.rows+column] = value; + public void set_upbo(int column, float value) { + if(column > columns || column < 1) throw new Error("column out of range"); + if(value < orig_lowbo[rows + column]) throw new Error("UpperBound must be >= lowerBound"); + eta_valid = FALSE; + orig_upbo[rows+column] = value; } - public void set_lowbo(Problem lp, int column, float value) - { - if(column > lp.columns || column < 1) - System.err.print("column out of range"); - if (lp.scaling_used != FALSE) - value /= lp.scale[lp.rows + column]; - if(value > lp.orig_upbo[lp.rows + column]) - System.err.print("UpperBound must be >= lowerBound"); - lp.eta_valid = FALSE; - lp.orig_lowbo[lp.rows+column] = value; + public void set_lowbo(int column, float value) { + if(column > columns || column < 1) throw new Error("column out of range"); + if(value > orig_upbo[rows + column]) throw new Error("UpperBound must be >= lowerBound"); + eta_valid = FALSE; + orig_lowbo[rows+column] = value; } - public void set_rh(Problem lp, int row, float value) - { - if(row > lp.rows || row < 0) - System.err.print("Row out of Range"); - if(row == 0) /* setting of RHS of OF not meaningful */ - { - System.err.println("Warning: attempt to set RHS of objective function, ignored"); - return; - } - if (lp.scaling_used != FALSE) - if (lp.ch_sign[row] != FALSE) - lp.orig_rh[row] = -value * lp.scale[row]; - else - lp.orig_rh[row] = value * lp.scale[row]; - else - if (lp.ch_sign[row] != FALSE) - lp.orig_rh[row] = -value; - else - lp.orig_rh[row] = value; - lp.eta_valid = FALSE; + public void set_rh(int row, float value) { + if(row > rows || row < 0) throw new Error("Row out of Range"); + if(row == 0) throw new Error("Warning: attempt to set RHS of objective function, ignored"); + if (ch_sign[row] != FALSE) orig_rh[row] = -value; + else orig_rh[row] = value; + eta_valid = FALSE; } - public void set_rh_vec(Problem lp, float[] rh) - { - int i; - if (lp.scaling_used != FALSE) - for(i = 1; i <= lp.rows; i++) - if(lp.ch_sign[i] != FALSE) - lp.orig_rh[i]=-rh[i]*lp.scale[i]; - else - lp.orig_rh[i]=rh[i]*lp.scale[i]; - else - for(i=1; i <= lp.rows; i++) - if(lp.ch_sign[i] != FALSE) - lp.orig_rh[i]=-rh[i]; - else - lp.orig_rh[i]=rh[i]; - lp.eta_valid=FALSE; + public void set_rh_vec(float[] rh) { + for(int i=1; i <= rows; i++) + if (ch_sign[i] != FALSE) orig_rh[i]=-rh[i]; + else orig_rh[i]=rh[i]; + eta_valid=FALSE; } - public void set_maxim(Problem lp) - { - int i; - if(lp.maximise==FALSE) - { - for(i = 0; i < lp.non_zeros; i++) - if(lp.mat[i].row_nr==0) - lp.mat[i].value*=-1; - lp.eta_valid=FALSE; - } - lp.maximise=TRUE; - lp.ch_sign[0]=TRUE; - if (lp.active != FALSE) - Maximise=TRUE; - } - public void set_minim(Problem lp) - { - int i; - if(lp.maximise==TRUE) - { - for(i = 0; i < lp.non_zeros; i++) - if(lp.mat[i].row_nr==0) - lp.mat[i].value = -lp.mat[i].value; - lp.eta_valid=FALSE; - } - lp.maximise=FALSE; - lp.ch_sign[0]=FALSE; - if (lp.active != FALSE) - Maximise=FALSE; + public void set_constr_type(int row, short con_type) { + if (row > rows || row < 1) throw new Error("Row out of Range"); + switch(con_type) { + case EQ: + orig_upbo[row]=0; + basis_valid=FALSE; + if (ch_sign[row] != FALSE) { + for(int i = 0; i < non_zeros; i++) + if (get_row_nr(mat, i)==row) set_value(mat, i, get_value(mat,i) * (float)-1); + eta_valid=FALSE; + ch_sign[row]=FALSE; + if (orig_rh[row]!=0) orig_rh[row]*=-1; + } + break; + case LE: + orig_upbo[row]=infinite; + basis_valid=FALSE; + if (ch_sign[row] != FALSE) { + for(int i = 0; i < non_zeros; i++) + if (get_row_nr(mat, i)==row) set_value(mat, i, get_value(mat,i) * (float)-1); + eta_valid=FALSE; + ch_sign[row]=FALSE; + if (orig_rh[row]!=0) orig_rh[row]*=-1; + } + break; + case GE: + orig_upbo[row]=infinite; + basis_valid=FALSE; + if (ch_sign[row] == FALSE) { + for(int i = 0; i < non_zeros; i++) + if (get_row_nr(mat, i)==row) set_value(mat, i, get_value(mat,i) * (float)-1); + eta_valid=FALSE; + ch_sign[row]=TRUE; + if (orig_rh[row]!=0) orig_rh[row]*=-1; + } + break; + default: throw new Error("Constraint type not (yet) implemented"); + } } - public void set_constr_type(Problem lp, int row, short con_type) - { - int i; - if(row > lp.rows || row < 1) - System.err.print("Row out of Range"); - if(con_type==EQ) - { - lp.orig_upbo[row]=0; - lp.basis_valid=FALSE; - if (lp.ch_sign[row] != FALSE) - { - for(i = 0; i < lp.non_zeros; i++) - if(lp.mat[i].row_nr==row) - lp.mat[i].value*=-1; - lp.eta_valid=FALSE; - lp.ch_sign[row]=FALSE; - if(lp.orig_rh[row]!=0) - lp.orig_rh[row]*=-1; - } - } - else if(con_type==LE) - { - lp.orig_upbo[row]=lp.infinite; - lp.basis_valid=FALSE; - if (lp.ch_sign[row] != FALSE) - { - for(i = 0; i < lp.non_zeros; i++) - if(lp.mat[i].row_nr==row) - lp.mat[i].value*=-1; - lp.eta_valid=FALSE; - lp.ch_sign[row]=FALSE; - if(lp.orig_rh[row]!=0) - lp.orig_rh[row]*=-1; - } - } - else if(con_type==GE) - { - lp.orig_upbo[row]=lp.infinite; - lp.basis_valid=FALSE; - if(lp.ch_sign[row] == FALSE) - { - for(i = 0; i < lp.non_zeros; i++) - if(lp.mat[i].row_nr==row) - lp.mat[i].value*=-1; - lp.eta_valid=FALSE; - lp.ch_sign[row]=TRUE; - if(lp.orig_rh[row]!=0) - lp.orig_rh[row]*=-1; - } - } - else - System.err.print("Constraint type not (yet) implemented"); + void set_globals() { + Rows = rows; + columns = columns; + Sum = Rows + columns; + Non_zeros = non_zeros; + Mat = mat; + Col_no = col_no; + Col_end = col_end; + Row_end = row_end; + Rh = rh; + Rhs = rhs; + Orig_rh = orig_rh; + Orig_upbo = orig_upbo; + Orig_lowbo = orig_lowbo; + Upbo = upbo; + Lowbo = lowbo; + Bas = bas; + Basis = basis; + Lower = lower; + Eta_alloc = eta_alloc; + Eta_size = eta_size; + Num_inv = num_inv; + Eta_value = eta_value; + Eta_row_nr = eta_row_nr; + Eta_col_end = eta_col_end; + Solution = solution; + Best_solution = best_solution; + Infinite = infinite; + Epsilon = epsilon; + Epsb = epsb; + Epsd = epsd; + Epsel = epsel; + TREJ = TREJ; + TINV = TINV; + Maximise = maximise; + Floor_first = floor_first; + active = TRUE; } - public float mat_elm(Problem lp, int row, int column) - { - float value; - int elmnr; - if(row < 0 || row > lp.rows) - System.err.print("Row out of range in mat_elm"); - if(column < 1 || column > lp.columns) - System.err.print("column out of range in mat_elm"); - value=0; - elmnr=lp.col_end[column-1]; - while(lp.mat[elmnr].row_nr != row && elmnr < lp.col_end[column]) - elmnr++; - if(elmnr != lp.col_end[column]) - { - value = lp.mat[elmnr].value; - if (lp.ch_sign[row] != FALSE) - value = -value; - if (lp.scaling_used != FALSE) - value /= lp.scale[row] * lp.scale[lp.rows + column]; - } - return(value); + private void ftran(int start, int end, float[] pcol) { + int k, r; + float theta; + for(int i = start; i <= end; i++) { + k = Eta_col_end[i] - 1; + r = Eta_row_nr[k]; + theta = pcol[r]; + if (theta != 0) for(int j = Eta_col_end[i - 1]; j < k; j++) + pcol[Eta_row_nr[j]] += theta * Eta_value[j]; + pcol[r] *= Eta_value[k]; + } + for(int i = 0; i <= Rows; i++) round(pcol[i], Epsel); } - - public void get_row(Problem lp, int row_nr, float[] row) - { - int i, j; - - if(row_nr < 0 || row_nr > lp.rows) - System.err.print("Row nr. out of range in get_row"); - for(i = 1; i <= lp.columns; i++) - { - row[i]=0; - for(j=lp.col_end[i-1]; j < lp.col_end[i]; j++) - if(lp.mat[j].row_nr==row_nr) - row[i]=lp.mat[j].value; - if (lp.scaling_used != FALSE) - row[i] /= lp.scale[lp.rows+i] * lp.scale[row_nr]; - } - if(lp.ch_sign[row_nr] != FALSE) - for(i=0; i <= lp.columns; i++) - if(row[i]!=0) - row[i] = -row[i]; + private void btran(float[] row) { + int k; + float f; + for(int i = Eta_size; i >= 1; i--) { + f = 0; + k = Eta_col_end[i] - 1; + for(int j = Eta_col_end[i - 1]; j <= k; j++) f += row[Eta_row_nr[j]] * Eta_value[j]; + f = round(f, Epsel); + row[Eta_row_nr[k]] = f; + } } - public void get_column(Problem lp, int col_nr, float[] column) - { - int i; - - if(col_nr < 1 || col_nr > lp.columns) - System.err.print("Col. nr. out of range in get_column"); - for(i = 0; i <= lp.rows; i++) - column[i]=0; - for(i = lp.col_end[col_nr-1]; i < lp.col_end[col_nr]; i++) - column[lp.mat[i].row_nr]=lp.mat[i].value; - for(i = 0; i <= lp.rows; i++) - if(column[i] !=0) - { - if(lp.ch_sign[i] != FALSE) - column[i]*=-1; - if (lp.scaling_used != FALSE) - column[i]/=(lp.scale[i] * lp.scale[lp.rows+col_nr]); + static int[] num = new int[65535]; + static int[] rownum = new int[65535]; + static int[] colnum = new int[65535]; + + short Isvalid() { + int row_nr; + if (row_end_valid == FALSE) { + for(int i = 0; i <= rows; i++) { num[i] = 0; rownum[i] = 0; } + for(int i = 0; i < non_zeros; i++) rownum[get_row_nr(mat, i)]++; + row_end[0] = 0; + for(int i = 1; i <= rows; i++) row_end[i] = row_end[i - 1] + rownum[i]; + for(int i = 1; i <= columns; i++) + for(int j = col_end[i - 1]; j < col_end[i]; j++) { + row_nr = get_row_nr(mat, j); + if (row_nr != 0) { + num[row_nr]++; + col_no[row_end[row_nr - 1] + num[row_nr]] = i; + } } - } - - public void get_reduced_costs(Problem lp, float[] rc) - { - int varnr, i, j; - float f; + row_end_valid = TRUE; + } + if (valid != FALSE) return(TRUE); + for(int i = 0; i <= rows; i++) rownum[i] = 0; + for(int i = 0; i <= columns; i++) colnum[i] = 0; + for(int i = 1 ; i <= columns; i++) + for(int j = col_end[i - 1]; j < col_end[i]; j++) { + colnum[i]++; + rownum[get_row_nr(mat, j)]++; + } + for(int i = 1; i <= columns; i++) + if (colnum[i] == 0) + throw new Error("Warning: Variable " + i + " not used in any constaints\n"); + valid = TRUE; + return(TRUE); + } - if(lp.basis_valid == FALSE) - System.err.print("Not a valid basis in get_reduced_costs"); - set_globals(lp); - if(lp.eta_valid == FALSE) - invert(); - for(i = 1; i <= lp.sum; i++) - rc[i] = 0; - rc[0] = 1; - btran(rc); - for(i = 1; i <= lp.columns; i++) - { - varnr = lp.rows + i; - if(Basis[varnr] == FALSE) - if(Upbo[varnr] > 0) - { - f = 0; - for(j = Col_end[i - 1]; j < Col_end[i]; j++) - f += rc[Mat[j].row_nr] * Mat[j].value; - rc[varnr] = f; - } - } - for(i = 1; i <= Sum; i++) - rc[i] = round(rc[i], Epsd); - } - - short is_feasible(Problem lp, float[] values) - { - int i, elmnr; - float[] this_rhs; - float dist; - - if (lp.scaling_used != FALSE) - { - for(i = lp.rows+1; i <= lp.sum; i++) - if( values[i - lp.rows] < lp.orig_lowbo[i]*lp.scale[i] - || values[i - lp.rows] > lp.orig_upbo[i]*lp.scale[i]) - return(FALSE); - } - else - { - for(i = lp.rows+1; i <= lp.sum; i++) - if( values[i - lp.rows] < lp.orig_lowbo[i] - || values[i - lp.rows] > lp.orig_upbo[i]) - return(FALSE); - } - this_rhs = new float[lp.rows+1]; - for (i = 0; i <= lp.rows; i++) - this_rhs[i] = 0; - - for(i = 1; i <= lp.columns; i++) - for(elmnr = lp.col_end[i - 1]; elmnr < lp.col_end[i]; elmnr++) - this_rhs[lp.mat[elmnr].row_nr] += lp.mat[elmnr].value * values[i]; - for(i = 1; i <= lp.rows; i++) - { - dist = lp.orig_rh[i] - this_rhs[i]; - dist = round(dist, (float)0.001); - if((lp.orig_upbo[i] == 0 && dist != 0) || dist < 0) - { - return(FALSE); - } - } - return(TRUE); - } - - short column_in_lp(Problem lp, float[] testcolumn) - { - int i, j; - short ident; - float value; - - if (lp.scaling_used != FALSE) - for(i = 1; i <= lp.columns; i++) - { - ident = TRUE; - j = lp.col_end[i-1]; - while(ident != FALSE && (j < lp.col_end[i])) - { - value = lp.mat[j].value; - if(lp.ch_sign[lp.mat[j].row_nr] != FALSE) - value = -value; - value /= lp.scale[lp.rows+i]; - value /= lp.scale[lp.mat[j].row_nr]; - value -= testcolumn[lp.mat[j].row_nr]; - if(Math.abs(value) > (float)0.001) /* should be some epsilon? MB */ - ident=FALSE; - j++; - } - if(ident != FALSE) - return(TRUE); - } - else - for(i = 1; i <= lp.columns; i++) - { - ident = TRUE; - j = lp.col_end[i-1]; - while(ident != FALSE && j < lp.col_end[i]) - { - value = lp.mat[j].value; - if(lp.ch_sign[lp.mat[j].row_nr] != FALSE) - value *= -1; - value -= testcolumn[lp.mat[j].row_nr]; - if( Math.abs(value) > (float)0.001 ) - ident=FALSE; - j++; - } - if(ident != FALSE) - return(TRUE); - } - return(FALSE); - } - - private float minmax_to_scale(float min, float max) - { - float scale; - - /* should do something sensible when min or max is 0, MB */ - if((min == 0) || (max == 0)) - return((float)1); - - scale = (float)(1 / Math.pow(Math.E, (Math.log(min) + Math.log(max)) / 2)); - return(scale); - } - - public void unscale_columns(Problem lp) - { - int i, j; - - /* unscale mat */ - for(j = 1; j <= lp.columns; j++) - for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) - lp.mat[i].value /= lp.scale[lp.rows + j]; - - /* unscale Bounds as well */ - for(i = lp.rows + 1; i < lp.sum; i++) - { - if(lp.orig_lowbo[i] != 0) - lp.orig_lowbo[i] *= lp.scale[i]; - if(lp.orig_upbo[i] != lp.infinite) - lp.orig_upbo[i] *= lp.scale[i]; - } - - for(i=lp.rows+1; i<= lp.sum; i++) - lp.scale[i]=1; - lp.columns_scaled=FALSE; - lp.eta_valid=FALSE; - } - - public void unscale(Problem lp) - { - int i, j; - - if (lp.scaling_used != FALSE) - { - - /* unscale mat */ - for(j = 1; j <= lp.columns; j++) - for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) - lp.mat[i].value /= lp.scale[lp.rows + j]; - - /* unscale Bounds */ - for(i = lp.rows + 1; i < lp.sum; i++) - { - if(lp.orig_lowbo[i] != 0) - lp.orig_lowbo[i] *= lp.scale[i]; - if(lp.orig_upbo[i] != lp.infinite) - lp.orig_upbo[i] *= lp.scale[i]; - } - - /* unscale the matrix */ - for(j = 1; j <= lp.columns; j++) - for(i = lp.col_end[j-1]; i < lp.col_end[j]; i++) - lp.mat[i].value /= lp.scale[lp.mat[i].row_nr]; - - /* unscale the rhs! */ - for(i = 0; i <= lp.rows; i++) - lp.orig_rh[i] /= lp.scale[i]; - - lp.scaling_used=FALSE; - lp.eta_valid=FALSE; - } - } - - - public void auto_scale(Problem lp) - { - int i, j, row_nr, IntUsed; - float row_max[], row_min[], scalechange[], absval; - - if(lp.scaling_used == 0) - { - lp.scale = new float[lp.sum_alloc + 1]; - for(i = 0; i <= lp.sum; i++) - lp.scale[i]=1; - } - - row_max = new float[lp.rows + 1]; - row_min = new float[lp.rows + 1]; - scalechange = new float[lp.sum + 1]; - - /* initialise min and max values */ - for(i = 0; i <= lp.rows; i++) - { - row_max[i]=0; - row_min[i]=lp.infinite; - } - - /* calculate min and max absolute values of rows */ - for(j = 1; j <= lp.columns; j++) - for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) - { - row_nr = lp.mat[i].row_nr; - absval = Math.abs(lp.mat[i].value); - if(absval!=0) - { - row_max[row_nr] = Math.max(row_max[row_nr], absval); - row_min[row_nr] = Math.min(row_min[row_nr], absval); - } - } - /* calculate scale factors for rows */ - for(i = 0; i <= lp.rows; i++) - { - scalechange[i] = minmax_to_scale(row_min[i], row_max[i]); - lp.scale[i] *= scalechange[i]; - } - - /* now actually scale the matrix */ - for(j = 1; j <= lp.columns; j++) - for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) - lp.mat[i].value *= scalechange[lp.mat[i].row_nr]; - - /* and scale the rhs! */ - for(i = 0; i <= lp.rows; i++) - lp.orig_rh[i] *= scalechange[i]; - - IntUsed=FALSE; - i=lp.rows+1; - while(IntUsed == FALSE && i <= lp.sum) - { - IntUsed=FALSE; - i++; - } - if(IntUsed == FALSE) - { - float col_max, col_min; - - /* calculate scales */ - for(j = 1; j <= lp.columns; j++) - { - col_max = 0; - col_min = lp.infinite; - for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) - { - if(lp.mat[i].value!=0) - { - col_max = Math.max(col_max, Math.abs(lp.mat[i].value)); - col_min = Math.min(col_min, Math.abs(lp.mat[i].value)); - } - } - scalechange[lp.rows + j] = minmax_to_scale(col_min, col_max); - lp.scale[lp.rows + j] *= scalechange[lp.rows + j]; - } - - /* scale mat */ - for(j = 1; j <= lp.columns; j++) - for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) - lp.mat[i].value *= scalechange[lp.rows + j]; - - /* scale Bounds as well */ - for(i = lp.rows + 1; i < lp.sum; i++) - { - if(lp.orig_lowbo[i] != 0) - lp.orig_lowbo[i] /= scalechange[i]; - if(lp.orig_upbo[i] != lp.infinite) - lp.orig_upbo[i] /= scalechange[i]; - } - lp.columns_scaled=TRUE; - } - lp.scaling_used=TRUE; - lp.eta_valid=FALSE; - } - - public void reset_basis(Problem lp) { lp.basis_valid=FALSE; } - - /* Globals used by solver */ - short JustInverted; - short Status; - short Doiter; - short DoInvert; - short Break_bb; - - void set_globals(Problem lp) - { - if(Lp != null) - Lp.active = FALSE; - Lp = lp; - Rows = lp.rows; - columns = lp.columns; - Sum = Rows + columns; - Non_zeros = lp.non_zeros; - Mat = lp.mat; - Col_no = lp.col_no; - Col_end = lp.col_end; - Row_end = lp.row_end; - Rh = lp.rh; - Rhs = lp.rhs; - Orig_rh = lp.orig_rh; - Orig_upbo = lp.orig_upbo; - Orig_lowbo = lp.orig_lowbo; - Upbo = lp.upbo; - Lowbo = lp.lowbo; - Bas = lp.bas; - Basis = lp.basis; - Lower = lp.lower; - Eta_alloc = lp.eta_alloc; - Eta_size = lp.eta_size; - Num_inv = lp.num_inv; - Eta_value = lp.eta_value; - Eta_row_nr = lp.eta_row_nr; - Eta_col_end = lp.eta_col_end; - Solution = lp.solution; - Best_solution = lp.best_solution; - Infinite = lp.infinite; - Epsilon = lp.epsilon; - Epsb = lp.epsb; - Epsd = lp.epsd; - Epsel = lp.epsel; - - /* ?? MB */ - TREJ = TREJ; - TINV = TINV; - - Maximise = lp.maximise; - Floor_first = lp.floor_first; - lp.active = TRUE; - - // System.out.println("Sum = " + Sum); - } // end of set_global - - private void ftran(int start, - int end, - float[] pcol) - { - int i, j, k, r; - float theta; - - for(i = start; i <= end; i++) - { - k = Eta_col_end[i] - 1; - r = Eta_row_nr[k]; - theta = pcol[r]; - if(theta != 0) - for(j = Eta_col_end[i - 1]; j < k; j++) - pcol[Eta_row_nr[j]] += theta * Eta_value[j]; /* cpu expensive line */ - pcol[r] *= Eta_value[k]; - } - /* round small values to zero */ - for(i = 0; i <= Rows; i++) - round(pcol[i], Epsel); - } /* ftran */ - - private void btran(float[] row) - { - int i, j, k; - float f; - - for(i = Eta_size; i >= 1; i--) - { - f = 0; - k = Eta_col_end[i] - 1; - for(j = Eta_col_end[i - 1]; j <= k; j++) - f += row[Eta_row_nr[j]] * Eta_value[j]; - f = round(f, Epsel); - row[Eta_row_nr[k]] = f; - } - } /* btran */ - - - static short Isvalid(Problem lp) - { - int i, j, rownum[], colnum[]; - int num[], row_nr; - - if(lp.row_end_valid == FALSE) - { - num = new int[lp.rows + 1]; - rownum = new int[lp.rows + 1]; - for(i = 0; i <= lp.rows; i++) - { - num[i] = 0; - rownum[i] = 0; - } - for(i = 0; i < lp.non_zeros; i++) - rownum[lp.mat[i].row_nr]++; - lp.row_end[0] = 0; - for(i = 1; i <= lp.rows; i++) - lp.row_end[i] = lp.row_end[i - 1] + rownum[i]; - for(i = 1; i <= lp.columns; i++) - for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++) - { - row_nr = lp.mat[j].row_nr; - if(row_nr != 0) - { - num[row_nr]++; - lp.col_no[lp.row_end[row_nr - 1] + num[row_nr]] = i; - } - } - lp.row_end_valid = TRUE; - } - if(lp.valid != FALSE) - return(TRUE); - rownum = new int[lp.rows + 1]; - for (i = 0; i <= lp.rows; i++) - rownum[i] = 0; - colnum = new int[lp.columns + 1]; - for (i = 0; i <= lp.columns; i++) - colnum[i] = 0; - for(i = 1 ; i <= lp.columns; i++) - for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++) - { - colnum[i]++; - rownum[lp.mat[j].row_nr]++; - } - for(i = 1; i <= lp.columns; i++) - if(colnum[i] == 0) - if (lp.names_used != FALSE) - System.err.print("Warning: Variable " + lp.col_name[i] + - " not used in any constraints\n"); - else - System.err.print("Warning: Variable " + i + " not used in any constaints\n"); - lp.valid = TRUE; - return(TRUE); - } - - private void resize_eta() - { + private void resize_eta() { Eta_alloc *= 2; + throw new Error("eta undersized; this should never happen"); + /* float[] db_ptr = Eta_value; Eta_value = new float[Eta_alloc]; System.arraycopy(db_ptr, 0, Eta_value, 0, db_ptr.length); - Lp.eta_value = Eta_value; - + eta_value = Eta_value; + int[] int_ptr = Eta_row_nr; Eta_row_nr = new int[Eta_alloc]; System.arraycopy(int_ptr, 0, Eta_row_nr, 0, int_ptr.length); - Lp.eta_row_nr = Eta_row_nr; - } /* resize_eta */ - + eta_row_nr = Eta_row_nr; + */ + } - private void condensecol(int row_nr, - float[] pcol) - { - int i, elnr; - + private void condensecol(int row_nr, float[] pcol) { + int elnr; elnr = Eta_col_end[Eta_size]; - - if(elnr + Rows + 2 > Eta_alloc) /* maximum local growth of Eta */ - resize_eta(); - - for(i = 0; i <= Rows; i++) - if(i != row_nr && pcol[i] != 0) - { - Eta_row_nr[elnr] = i; - Eta_value[elnr] = pcol[i]; - elnr++; - } + if (elnr + Rows + 2 > Eta_alloc) resize_eta(); + for(int i = 0; i <= Rows; i++) + if (i != row_nr && pcol[i] != 0) { + Eta_row_nr[elnr] = i; + Eta_value[elnr] = pcol[i]; + elnr++; + } Eta_row_nr[elnr] = row_nr; Eta_value[elnr] = pcol[row_nr]; elnr++; Eta_col_end[Eta_size + 1] = elnr; - } /* condensecol */ - + } - private void addetacol() - { - int i, j, k; + private void addetacol() { + int k; float theta; - - j = Eta_col_end[Eta_size]; + int j = Eta_col_end[Eta_size]; Eta_size++; k = Eta_col_end[Eta_size]; theta = 1 / (float) Eta_value[k - 1]; Eta_value[k - 1] = theta; - for(i = j; i < k - 1; i++) - Eta_value[i] *= -theta; + for(int i = j; i < k - 1; i++) Eta_value[i] *= -theta; JustInverted = FALSE; - } /* addetacol */ + } - private void setpivcol(short lower, - int varin, - float[] pcol) - { - int i, colnr; + private void setpivcol(short lower, int varin, float[] pcol) { + int colnr; float f; - - if(lower != FALSE) - f = 1; - else - f = -1; - for(i = 0; i <= Rows; i++) - pcol[i] = 0; - if(varin > Rows) - { - colnr = varin - Rows; - for(i = Col_end[colnr - 1]; i < Col_end[colnr]; i++) - pcol[Mat[i].row_nr] = Mat[i].value * f; - pcol[0] -= Extrad * f; - } - else - if(lower != FALSE) - pcol[varin] = 1; - else - pcol[varin] = -1; + if (lower != FALSE) f = 1; + else f = -1; + for(int i = 0; i <= Rows; i++) pcol[i] = 0; + if (varin > Rows) { + colnr = varin - Rows; + for(int i = Col_end[colnr - 1]; i < Col_end[colnr]; i++) pcol[get_row_nr(Mat, i)] = get_value(Mat,i) * f; + pcol[0] -= Extrad * f; + } else { + if (lower != FALSE) pcol[varin] = 1; + else pcol[varin] = -1; + } ftran(1, Eta_size, pcol); - } /* setpivcol */ - + } - private void minoriteration(int colnr, - int row_nr) - { - int i, j, k, wk, varin, varout, elnr; + private void minoriteration(int colnr, int row_nr) { + int k, wk, varin, varout, elnr; float piv = 0, theta; - varin = colnr + Rows; elnr = Eta_col_end[Eta_size]; wk = elnr; Eta_size++; - if(Extrad != 0) - { - Eta_row_nr[elnr] = 0; - Eta_value[elnr] = -Extrad; + if (Extrad != 0) { + Eta_row_nr[elnr] = 0; + Eta_value[elnr] = -Extrad; + elnr++; + } + for(int j = Col_end[colnr - 1] ; j < Col_end[colnr]; j++) { + k = get_row_nr(Mat, j); + if (k == 0 && Extrad != 0) Eta_value[Eta_col_end[Eta_size -1]] += get_value(Mat,j); + else if (k != row_nr) { + Eta_row_nr[elnr] = k; + Eta_value[elnr] = get_value(Mat,j); elnr++; + } else { + piv = get_value(Mat,j); } - for(j = Col_end[colnr - 1] ; j < Col_end[colnr]; j++) - { - k = Mat[j].row_nr; - if(k == 0 && Extrad != 0) - Eta_value[Eta_col_end[Eta_size -1]] += Mat[j].value; - else if(k != row_nr) - { - Eta_row_nr[elnr] = k; - Eta_value[elnr] = Mat[j].value; - elnr++; - } - else - piv = Mat[j].value; - } + } Eta_row_nr[elnr] = row_nr; Eta_value[elnr] = 1 / (float) piv; elnr++; theta = Rhs[row_nr] / (float) piv; Rhs[row_nr] = theta; - for(i = wk; i < elnr - 1; i++) - Rhs[Eta_row_nr[i]] -= theta * Eta_value[i]; + for(int i = wk; i < elnr - 1; i++) Rhs[Eta_row_nr[i]] -= theta * Eta_value[i]; varout = Bas[row_nr]; Bas[row_nr] = varin; Basis[varout] = FALSE; Basis[varin] = TRUE; - for(i = wk; i < elnr - 1; i++) - Eta_value[i] /= - (float) piv; + for(int i = wk; i < elnr - 1; i++) Eta_value[i] /= - (float) piv; Eta_col_end[Eta_size] = elnr; - } /* minoriteration */ + } - private void rhsmincol(float theta, - int row_nr, - int varin) - { - int i, j, k, varout; + private void rhsmincol(float theta, int row_nr, int varin) { + int varout; float f; - - if(row_nr > Rows + 1) - { - System.err.print("Error: rhsmincol called with row_nr: " + - row_nr + ", rows: " + Rows + "\n"); - System.err.print("This indicates numerical instability\n"); - System.exit(FAIL); - } - j = Eta_col_end[Eta_size]; - k = Eta_col_end[Eta_size + 1]; - for(i = j; i < k; i++) - { - f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i]; - f = round(f, Epsb); - Rhs[Eta_row_nr[i]] = f; - } + if (row_nr > Rows + 1) { + System.err.println("Error: rhsmincol called with row_nr: " + row_nr + ", rows: " + Rows + "\n"); + System.err.println("This indicates numerical instability\n"); + } + int j = Eta_col_end[Eta_size]; + int k = Eta_col_end[Eta_size + 1]; + for(int i = j; i < k; i++) { + f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i]; + f = round(f, Epsb); + Rhs[Eta_row_nr[i]] = f; + } Rhs[row_nr] = theta; varout = Bas[row_nr]; Bas[row_nr] = varin; Basis[varout] = FALSE; Basis[varin] = TRUE; - } /* rhsmincol */ - - - void invert() - { - int i, j, v, wk, numit, varnr, row_nr, colnr, varin; - float theta; - float[] pcol; - short[] frow; - short[] fcol; - int[] rownum, col, row; - int[] colnum; - - rownum = new int[Rows + 1]; - for (i = 0; i <= Rows; i++) - rownum[i] = 0; - col = new int[Rows + 1]; - for (i = 0; i <= Rows; i++) - col[i] = 0; - row = new int[Rows + 1]; - for (i = 0; i <= Rows; i++) - row[i] = 0; - pcol = new float[Rows + 1]; - for (i = 0; i <= Rows; i++) - pcol[i] = 0; - frow = new short[Rows + 1]; - for(i = 0; i <= Rows; i++) - frow[i] = TRUE; - fcol = new short[columns + 1]; - for(i = 0; i < columns; i++) - fcol[i] = FALSE; - colnum = new int[columns + 1]; - for(i = 0; i <= columns; i++) - colnum[i] = 0; - - for(i = 0; i <= Rows; i++) - if(Bas[i] > Rows) - fcol[Bas[i] - Rows - 1] = TRUE; - else - frow[Bas[i]] = FALSE; - - for(i = 1; i <= Rows; i++) - if(frow[i] != FALSE) - for(j = Row_end[i - 1] + 1; j <= Row_end[i]; j++) - { - wk = Col_no[j]; - if(fcol[wk - 1] != FALSE) - { - colnum[wk]++; - rownum[i - 1]++; - } - } - for(i = 1; i <= Rows; i++) - Bas[i] = i; - for(i = 1; i <= Rows; i++) - Basis[i] = TRUE; - for(i = 1; i <= columns; i++) - Basis[i + Rows] = FALSE; - - for(i = 0; i <= Rows; i++) - Rhs[i] = Rh[i]; - - for(i = 1; i <= columns; i++) - { - varnr = Rows + i; - if(Lower[varnr] == FALSE) - { - theta = Upbo[varnr]; - for(j = Col_end[i - 1]; j < Col_end[i]; j++) - Rhs[Mat[j].row_nr] -= theta * Mat[j].value; + } + + private static int[] rownum_ = new int[65535]; + private static int[] colnum_ = new int[65535]; + private static int[] col = new int[65535]; + private static int[] row = new int[65535]; + private static float[] pcol = new float[65535]; + private static short[] frow = new short[65535]; + private static short[] fcol = new short[65535]; + + void invert() { + int v, wk, numit, varnr, row_nr, colnr, varin; + float theta; + + for(int i = 0; i <= Rows; i++) rownum_[i] = 0; + for(int i = 0; i <= Rows; i++) col[i] = 0; + for(int i = 0; i <= Rows; i++) row[i] = 0; + for(int i = 0; i <= Rows; i++) pcol[i] = 0; + for(int i = 0; i <= Rows; i++) frow[i] = TRUE; + for(int i = 0; i < columns; i++) fcol[i] = FALSE; + for(int i = 0; i <= columns; i++) colnum_[i] = 0; + + for(int i = 0; i <= Rows; i++) + if (Bas[i] > Rows) fcol[Bas[i] - Rows - 1] = TRUE; + else frow[Bas[i]] = FALSE; + + for(int i = 1; i <= Rows; i++) + if (frow[i] != FALSE) + for(int j = Row_end[i - 1] + 1; j <= Row_end[i]; j++) { + wk = Col_no[j]; + if (fcol[wk - 1] != FALSE) { + colnum_[wk]++; + rownum_[i - 1]++; } + } + + for(int i = 1; i <= Rows; i++) Bas[i] = i; + for(int i = 1; i <= Rows; i++) Basis[i] = TRUE; + for(int i = 1; i <= columns; i++) Basis[i + Rows] = FALSE; + for(int i = 0; i <= Rows; i++) Rhs[i] = Rh[i]; + for(int i = 1; i <= columns; i++) { + varnr = Rows + i; + if (Lower[varnr] == FALSE) { + theta = Upbo[varnr]; + for(int j = Col_end[i - 1]; j < Col_end[i]; j++) + Rhs[get_row_nr(Mat, j)] -= theta * get_value(Mat,j); } - for(i = 1; i <= Rows; i++) - if(Lower[i] == FALSE) - Rhs[i] -= Upbo[i]; + } + for(int i = 1; i <= Rows; i++) if (Lower[i] == FALSE) Rhs[i] -= Upbo[i]; Eta_size = 0; v = 0; row_nr = 0; Num_inv = 0; numit = 0; - while(v < Rows) - { - row_nr++; - if(row_nr > Rows) - row_nr = 1; - v++; - if(rownum[row_nr - 1] == 1) - if(frow[row_nr] != FALSE) - { - v = 0; - j = Row_end[row_nr - 1] + 1; - while(fcol[Col_no[j] - 1] == FALSE) - j++; - colnr = Col_no[j]; - fcol[colnr - 1] = FALSE; - colnum[colnr] = 0; - for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) - if(frow[Mat[j].row_nr] != FALSE) - rownum[Mat[j].row_nr - 1]--; - frow[row_nr] = FALSE; - minoriteration(colnr, row_nr); - } - } + while(v < Rows) { + int j; + row_nr++; + if (row_nr > Rows) row_nr = 1; + v++; + if (rownum_[row_nr - 1] == 1) + if (frow[row_nr] != FALSE) { + v = 0; + j = Row_end[row_nr - 1] + 1; + while(fcol[Col_no[j] - 1] == FALSE) j++; + colnr = Col_no[j]; + fcol[colnr - 1] = FALSE; + colnum_[colnr] = 0; + for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) + if (frow[get_row_nr(Mat, j)] != FALSE) + rownum_[get_row_nr(Mat, j) - 1]--; + frow[row_nr] = FALSE; + minoriteration(colnr, row_nr); + } + } v = 0; colnr = 0; - while(v columns) - colnr = 1; - v++; - if(colnum[colnr] == 1) - if(fcol[colnr - 1] != FALSE) - { - v = 0; - j = Col_end[colnr - 1] + 1; - while(frow[Mat[j - 1].row_nr] == FALSE) - j++; - row_nr = Mat[j - 1].row_nr; - frow[row_nr] = FALSE; - rownum[row_nr - 1] = 0; - for(j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++) - if(fcol[Col_no[j] - 1] != FALSE) - colnum[Col_no[j]]--; - fcol[colnr - 1] = FALSE; - numit++; - col[numit - 1] = colnr; - row[numit - 1] = row_nr; - } - } - for(j = 1; j <= columns; j++) - if(fcol[j - 1] != FALSE) - { - fcol[j - 1] = FALSE; - setpivcol(Lower[Rows + j], j + Rows, pcol); - row_nr = 1; - while((frow[row_nr] == FALSE || pcol[row_nr] == FALSE) - && row_nr <= Rows) - row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes - rhsmincol crash. Solved in 2.0? MB */ - if(row_nr == Rows + 1) - System.err.print("Inverting failed"); + while(v < columns) { + int j; + colnr++; + if (colnr > columns) colnr = 1; + v++; + if (colnum_[colnr] == 1) + if (fcol[colnr - 1] != FALSE) { + v = 0; + j = Col_end[colnr - 1] + 1; + while(frow[get_row_nr(Mat, j - 1)] == FALSE) j++; + row_nr = get_row_nr(Mat, j - 1); frow[row_nr] = FALSE; - condensecol(row_nr, pcol); - theta = Rhs[row_nr] / (float) pcol[row_nr]; - rhsmincol(theta, row_nr, Rows + j); - addetacol(); + rownum_[row_nr - 1] = 0; + for(j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++) + if (fcol[Col_no[j] - 1] != FALSE) + colnum_[Col_no[j]]--; + fcol[colnr - 1] = FALSE; + numit++; + col[numit - 1] = colnr; + row[numit - 1] = row_nr; } - for(i = numit - 1; i >= 0; i--) - { - colnr = col[i]; - row_nr = row[i]; - varin = colnr + Rows; - for(j = 0; j <= Rows; j++) - pcol[j] = 0; - for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) - pcol[Mat[j].row_nr] = Mat[j].value; - pcol[0] -= Extrad; + } + for(int j = 1; j <= columns; j++) + if (fcol[j - 1] != FALSE) { + fcol[j - 1] = FALSE; + setpivcol(Lower[Rows + j], j + Rows, pcol); + row_nr = 1; + while((frow[row_nr] == FALSE || pcol[row_nr] == FALSE) && row_nr <= Rows) + row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes + rhsmincol crash. Solved in 2.0? MB */ + if (row_nr == Rows + 1) throw new Error("Inverting failed"); + frow[row_nr] = FALSE; condensecol(row_nr, pcol); theta = Rhs[row_nr] / (float) pcol[row_nr]; - rhsmincol(theta, row_nr, varin); + rhsmincol(theta, row_nr, Rows + j); addetacol(); } - for(i = 1; i <= Rows; i++) - Rhs[i] = round(Rhs[i], Epsb); + for(int i = numit - 1; i >= 0; i--) { + colnr = col[i]; + row_nr = row[i]; + varin = colnr + Rows; + for(int j = 0; j <= Rows; j++) pcol[j] = 0; + for(int j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) pcol[get_row_nr(Mat, j)] = get_value(Mat,j); + pcol[0] -= Extrad; + condensecol(row_nr, pcol); + theta = Rhs[row_nr] / (float) pcol[row_nr]; + rhsmincol(theta, row_nr, varin); + addetacol(); + } + for(int i = 1; i <= Rows; i++) Rhs[i] = round(Rhs[i], Epsb); JustInverted = TRUE; DoInvert = FALSE; - } /* invert */ + } - private short colprim(Ref colnr, - short minit, - float[] drow) - { - int varnr, i, j; + private short colprim(Ref colnr, short minit, float[] drow) { + int varnr; float f, dpiv; - - dpiv = -Epsd; + dpiv = -Epsd; colnr.value = 0; - if(minit == FALSE) - { - for(i = 1; i <= Sum; i++) - drow[i] = 0; - drow[0] = 1; - btran(drow); - for(i = 1; i <= columns; i++) - { - varnr = Rows + i; - if(Basis[varnr] == FALSE) - if(Upbo[varnr] > 0) - { - f = 0; - for(j = Col_end[i - 1]; j < Col_end[i]; j++) - f += drow[Mat[j].row_nr] * Mat[j].value; - drow[varnr] = f; - } + if (minit == FALSE) { + for(int i = 1; i <= Sum; i++) drow[i] = 0; + drow[0] = 1; + btran(drow); + for(int i = 1; i <= columns; i++) { + varnr = Rows + i; + if (Basis[varnr] == FALSE) + if (Upbo[varnr] > 0) { + f = 0; + for(int j = Col_end[i - 1]; j < Col_end[i]; j++) f += drow[get_row_nr(Mat, j)] * get_value(Mat,j); + drow[varnr] = f; } - for(i = 1; i <= Sum; i++) - drow[i] = round(drow[i], Epsd); } - /* - System.out.println("dpiv = " + dpiv); - System.out.println("drow[] at colprim"); - for(i = 1; i <= Sum; i++) // DEBUG - { - System.out.print(drow[i] + " "); - } - System.out.println(); - */ - for(i = 1; i <= Sum; i++) - if(Basis[i] == FALSE) - if(Upbo[i] > 0) - { - if(Lower[i] != FALSE) - f = drow[i]; - else - f = -drow[i]; - if(f < dpiv) - { - dpiv = f; - colnr.value = i; - } + for(int i = 1; i <= Sum; i++) drow[i] = round(drow[i], Epsd); + } + for(int i = 1; i <= Sum; i++) + if (Basis[i] == FALSE) + if (Upbo[i] > 0) { + if (Lower[i] != FALSE) f = drow[i]; + else f = -drow[i]; + if (f < dpiv) { + dpiv = f; + colnr.value = i; } - if(Lp.trace != FALSE) - if(colnr.value > 0) - System.err.print("col_prim:" + colnr.value + ", reduced cost: " + dpiv + "\n"); - else - System.err.print("col_prim: no negative reduced costs found, optimality!\n"); - if(colnr.value == 0) - { - Doiter = FALSE; - DoInvert = FALSE; - Status = OPTIMAL; - } + } + if (colnr.value == 0) { + Doiter = FALSE; + DoInvert = FALSE; + Status = OPTIMAL; + } return(colnr.value > 0 ? (short)1 : (short)0); - } /* colprim */ - - private short rowprim(int colnr, - Ref row_nr, - Ref theta, - float[] pcol) - { - int i; - float f = 0, quot; + } + private short rowprim(int colnr, Ref row_nr, Ref theta, float[] pcol) { + float f = 0, quot; row_nr.value = 0; theta.value = Infinite; - for(i = 1; i <= Rows; i++) - { - f = pcol[i]; - if(Math.abs(f) < TREJ) - f = 0; - if(f != 0) - { - quot = 2 * Infinite; - if(f > 0) - quot = Rhs[i] / (float) f; - else - if(Upbo[Bas[i]] < Infinite) - quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f; - round(quot, Epsel); - if(quot < theta.value) - { - theta.value = quot; - row_nr.value = i; - } - } - } - if(row_nr.value == 0) - for(i = 1; i <= Rows; i++) - { - f = pcol[i]; - if(f != 0) - { - quot = 2 * Infinite; - if(f > 0) - quot = Rhs[i] / (float) f; - else - if(Upbo[Bas[i]] < Infinite) - quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f; - quot = round(quot, Epsel); - if(quot < theta.value) - { - theta.value = quot; - row_nr.value = i; - } - } + for(int i = 1; i <= Rows; i++) { + f = pcol[i]; + if (Math.abs(f) < TREJ) f = 0; + if (f != 0) { + quot = 2 * Infinite; + if (f > 0) quot = Rhs[i] / (float) f; + else if (Upbo[Bas[i]] < Infinite) quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f; + round(quot, Epsel); + if (quot < theta.value) { + theta.value = quot; + row_nr.value = i; } - - if(theta.value < 0) - { - System.err.println("Warning: Numerical instability, qout = " + theta.value); - System.err.println("pcol[" + row_nr.value + "] = " + f + ", rhs[" + - row_nr.value + "] = " + Rhs[(int)row_nr.value] + - " , upbo = " + Upbo[Bas[(int)row_nr.value]]); } - if(row_nr.value == 0) - { - if(Upbo[colnr] == Infinite) - { - Doiter = FALSE; - DoInvert = FALSE; - Status = UNBOUNDED; - } - else - { - i = 1; - while(pcol[i] >= 0 && i <= Rows) - i++; - if(i > Rows) /* empty column with upperBound! */ - { - Lower[colnr] = FALSE; - Rhs[0] += Upbo[colnr]*pcol[0]; - Doiter = FALSE; - DoInvert = FALSE; - } - else if(pcol[i]<0) - { - row_nr.value = i; - } + } + if (row_nr.value == 0) + for(int i = 1; i <= Rows; i++) { + f = pcol[i]; + if (f != 0) { + quot = 2 * Infinite; + if (f > 0) quot = Rhs[i] / (float) f; + else if (Upbo[Bas[i]] < Infinite) quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f; + quot = round(quot, Epsel); + if (quot < theta.value) { + theta.value = quot; + row_nr.value = i; } + } } - if(row_nr.value > 0) - Doiter = TRUE; - if(Lp.trace != FALSE) - System.out.println("row_prim:" + row_nr.value + ", pivot element:" + - pcol[(int)row_nr.value]); + if (theta.value < 0) throw new Error("Warning: Numerical instability, qout = " + theta.value); + if (row_nr.value == 0) { + if (Upbo[colnr] == Infinite) { + Doiter = FALSE; + DoInvert = FALSE; + Status = UNBOUNDED; + } else { + int i = 1; + while(pcol[i] >= 0 && i <= Rows) i++; + if (i > Rows) { + Lower[colnr] = FALSE; + Rhs[0] += Upbo[colnr]*pcol[0]; + Doiter = FALSE; + DoInvert = FALSE; + } else if (pcol[i]<0) { + row_nr.value = i; + } + } + } + if (row_nr.value > 0) Doiter = TRUE; return((row_nr.value > 0) ? (short)1 : (short)0); - } /* rowprim */ + } - private short rowdual(Ref row_nr) - { + private short rowdual(Ref row_nr) { int i; float f, g, minrhs; short artifs; - row_nr.value = 0; minrhs = -Epsb; i = 0; artifs = FALSE; - while(i < Rows && artifs == FALSE) - { - i++; - f = Upbo[Bas[i]]; - if(f == 0 && (Rhs[i] != 0)) - { - artifs = TRUE; - row_nr.value = i; - } - else - { - if(Rhs[i] < f - Rhs[i]) - g = Rhs[i]; - else - g = f - Rhs[i]; - if(g < minrhs) - { - minrhs = g; - row_nr.value = i; - } - } - } - - if(Lp.trace != FALSE) - { - if(row_nr.value > 0) - { - System.out.println("row_dual:" + row_nr.value + - ", rhs of selected row: " - + Rhs[(int)row_nr.value]); - if(Upbo[Bas[(int)row_nr.value]] < Infinite) - System.out.println(" upper Bound of basis variable: " + - Upbo[Bas[(int)row_nr.value]]); - } - else - System.out.print("row_dual: no infeasibilities found\n"); + while(i < Rows && artifs == FALSE) { + i++; + f = Upbo[Bas[i]]; + if (f == 0 && (Rhs[i] != 0)) { + artifs = TRUE; + row_nr.value = i; + } else { + if (Rhs[i] < f - Rhs[i]) g = Rhs[i]; + else g = f - Rhs[i]; + if (g < minrhs) { + minrhs = g; + row_nr.value = i; + } } - + } return(row_nr.value > 0 ? (short)1 : (short)0); - } /* rowdual */ - - private short coldual(int row_nr, - Ref colnr, - short minit, - float[] prow, - float[] drow) - { - int i, j, r, varnr; + } + + private short coldual(int row_nr, Ref colnr, short minit, float[] prow, float[] drow) { + int r, varnr; float theta, quot, pivot, d, f, g; - Doiter = FALSE; - if(minit == FALSE) - { - for(i = 0; i <= Rows; i++) - { - prow[i] = 0; - drow[i] = 0; - } - drow[0] = 1; - prow[row_nr] = 1; - for(i = Eta_size; i >= 1; i--) - { - d = 0; - f = 0; - r = Eta_row_nr[Eta_col_end[i] - 1]; - for(j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++) - { - /* this is where the program consumes most cpu time */ - f += prow[Eta_row_nr[j]] * Eta_value[j]; - d += drow[Eta_row_nr[j]] * Eta_value[j]; - } - f = round(f, Epsel); - prow[r] = f; - d = round(d, Epsel); - drow[r] = d; - } - for(i = 1; i <= columns; i++) - { - varnr = Rows + i; - if(Basis[varnr] == FALSE) - { - d = - Extrad * drow[0]; - f = 0; - for(j = Col_end[i - 1]; j < Col_end[i]; j++) - { - d = d + drow[Mat[j].row_nr] * Mat[j].value; - f = f + prow[Mat[j].row_nr] * Mat[j].value; - } - drow[varnr] = d; - prow[varnr] = f; - } - } - for(i = 0; i <= Sum; i++) - { - prow[i] = round(prow[i], Epsel); - drow[i] = round(drow[i], Epsd); + if (minit == FALSE) { + for(int i = 0; i <= Rows; i++) { + prow[i] = 0; + drow[i] = 0; + } + drow[0] = 1; + prow[row_nr] = 1; + for(int i = Eta_size; i >= 1; i--) { + d = 0; + f = 0; + r = Eta_row_nr[Eta_col_end[i] - 1]; + for(int j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++) { + /* this is where the program consumes most cpu time */ + f += prow[Eta_row_nr[j]] * Eta_value[j]; + d += drow[Eta_row_nr[j]] * Eta_value[j]; + } + f = round(f, Epsel); + prow[r] = f; + d = round(d, Epsel); + drow[r] = d; + } + for(int i = 1; i <= columns; i++) { + varnr = Rows + i; + if (Basis[varnr] == FALSE) { + d = - Extrad * drow[0]; + f = 0; + for(int j = Col_end[i - 1]; j < Col_end[i]; j++) { + d = d + drow[get_row_nr(Mat, j)] * get_value(Mat,j); + f = f + prow[get_row_nr(Mat, j)] * get_value(Mat,j); } + drow[varnr] = d; + prow[varnr] = f; + } } - if(Rhs[row_nr] > Upbo[Bas[row_nr]]) - g = -1; - else - g = 1; + for(int i = 0; i <= Sum; i++) { + prow[i] = round(prow[i], Epsel); + drow[i] = round(drow[i], Epsd); + } + } + if (Rhs[row_nr] > Upbo[Bas[row_nr]]) g = -1; + else g = 1; pivot = 0; colnr.value = 0; theta = Infinite; - for(i = 1; i <= Sum; i++) - { - if(Lower[i] != FALSE) - d = prow[i] * g; - else - d = -prow[i] * g; - if((d < 0) && (Basis[i] == FALSE) && (Upbo[i] > 0)) - { - if(Lower[i] == FALSE) - quot = -drow[i] / (float) d; - else - quot = drow[i] / (float) d; - if(quot < theta) - { - theta = quot; - pivot = d; - colnr.value = i; - } - else if((quot == theta) && (Math.abs(d) > Math.abs(pivot))) - { - pivot = d; - colnr.value = i; - } - } + for(int i = 1; i <= Sum; i++) { + if (Lower[i] != FALSE) d = prow[i] * g; + else d = -prow[i] * g; + if ((d < 0) && (Basis[i] == FALSE) && (Upbo[i] > 0)) { + if (Lower[i] == FALSE) quot = -drow[i] / (float) d; + else quot = drow[i] / (float) d; + if (quot < theta) { + theta = quot; + pivot = d; + colnr.value = i; + } else if ((quot == theta) && (Math.abs(d) > Math.abs(pivot))) { + pivot = d; + colnr.value = i; + } } - if(Lp.trace != FALSE) - System.out.println("col_dual:" + colnr.value + ", pivot element: " + - prow[(int)colnr.value]); - - if(colnr.value > 0) - Doiter = TRUE; - + } + if (colnr.value > 0) Doiter = TRUE; return(colnr.value > 0 ? (short)1 : (short)0); - } /* coldual */ - - private void iteration(int row_nr, - int varin, - Ref theta, - float up, - Ref minit, - Ref low, - short primal, - float[] pcol) - { - int i, k, varout; + } + + private void iteration(int row_nr, int varin, Ref theta, float up, Ref minit, Ref low, short primal,float[] pcol) { + int k, varout; float f; float pivot; - - Lp.iter++; + iter++; minit.value = theta.value > (up + Epsb) ? 1 : 0; - if(minit.value != 0) - { - theta.value = up; - low.value = low.value == 0 ? 1 : 0; - } + if (minit.value != 0) { + theta.value = up; + low.value = low.value == 0 ? 1 : 0; + } k = Eta_col_end[Eta_size + 1]; pivot = Eta_value[k - 1]; - for(i = Eta_col_end[Eta_size]; i < k; i++) - { - f = Rhs[Eta_row_nr[i]] - theta.value * Eta_value[i]; - f = round(f, Epsb); - Rhs[Eta_row_nr[i]] = f; - } - if(minit.value == 0) - { - Rhs[row_nr] = theta.value; - varout = Bas[row_nr]; - Bas[row_nr] = varin; - Basis[varout] = FALSE; - Basis[varin] = TRUE; - if(primal != FALSE && pivot < 0) - Lower[varout] = FALSE; - if(low.value == 0 && up < Infinite) - { - low.value = TRUE; - Rhs[row_nr] = up - Rhs[row_nr]; - for(i = Eta_col_end[Eta_size]; i < k; i++) - Eta_value[i] = -Eta_value[i]; - } - addetacol(); - Num_inv++; - } - if(Lp.trace != FALSE) - { - System.out.print("Theta = " + theta.value + " "); - if(minit.value != 0) - { - if(Lower[varin] == FALSE) - System.out.print("Iteration:" + Lp.iter + - ", variable" + varin + " changed from 0 to its upper Bound of " - + Upbo[varin] + "\n"); - else - System.out.print("Iteration:" + Lp.iter + ", variable" + varin + - " changed its upper Bound of " + Upbo[varin] + " to 0\n"); - } - else - System.out.print("Iteration:" + Lp.iter + ", variable" + varin + - " entered basis at:" + Rhs[row_nr] + "\n"); - if(primal == 0) - { - f = 0; - for(i = 1; i <= Rows; i++) - if(Rhs[i] < 0) - f -= Rhs[i]; - else - if(Rhs[i] > Upbo[Bas[i]]) - f += Rhs[i] - Upbo[Bas[i]]; - System.out.println("feasibility gap of this basis:" + (float)f); - } - else - System.out.println("objective function value of this feasible basis: " + Rhs[0]); - } - } /* iteration */ + for(int i = Eta_col_end[Eta_size]; i < k; i++) { + f = Rhs[Eta_row_nr[i]] - theta.value * Eta_value[i]; + f = round(f, Epsb); + Rhs[Eta_row_nr[i]] = f; + } + if (minit.value == 0) { + Rhs[row_nr] = theta.value; + varout = Bas[row_nr]; + Bas[row_nr] = varin; + Basis[varout] = FALSE; + Basis[varin] = TRUE; + if (primal != FALSE && pivot < 0) Lower[varout] = FALSE; + if (low.value == 0 && up < Infinite) { + low.value = TRUE; + Rhs[row_nr] = up - Rhs[row_nr]; + for(int i = Eta_col_end[Eta_size]; i < k; i++) Eta_value[i] = -Eta_value[i]; + } + addetacol(); + Num_inv++; + } + } + static float[] drow = new float[65535]; + static float[] prow = new float[65535]; + static float[] Pcol = new float[65535]; - private int solvelp() - { - int i, j, varnr; + private int solvelp() { + int varnr; float f = 0, theta = 0; short primal; - float[] drow, prow, Pcol; short minit; int colnr, row_nr; colnr = 0; @@ -1873,1015 +1094,248 @@ public class LinearProgramming { ref1 = new Ref(0); ref2 = new Ref(0); ref3 = new Ref(0); - - drow = new float[Sum + 1]; - prow = new float[Sum + 1]; - Pcol = new float[Rows + 1]; - for (i = 0; i <= Sum; i++) { - drow[i] = 0; - prow[i] = 0; - } - for (i = 0; i <= Rows; i++) - Pcol[i] = 0; - - Lp.iter = 0; + for(int i = 0; i <= Sum; i++) { drow[i] = 0; prow[i] = 0; } + for(int i = 0; i <= Rows; i++) Pcol[i] = 0; + iter = 0; minit = FALSE; Status = RUNNING; DoInvert = FALSE; Doiter = FALSE; - i = 0; primal = TRUE; - while(i != Rows && primal != FALSE) - { - i++; - primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ? - (short)1: (short)0; - } - if(Lp.trace != FALSE) - { - if(primal != FALSE) - System.out.print("Start at feasible basis\n"); - else - System.out.print("Start at infeasible basis\n"); - } - if(primal == FALSE) - { - drow[0] = 1; - for(i = 1; i <= Rows; i++) - drow[i] = 0; - Extrad = 0; - for(i = 1; i <= columns; i++) - { - varnr = Rows + i; - drow[varnr] = 0; - for(j = Col_end[i - 1]; j < Col_end[i]; j++) - if(drow[Mat[j].row_nr] != 0) - drow[varnr] += drow[Mat[j].row_nr] * Mat[j].value; - if(drow[varnr] < Extrad) - Extrad = drow[varnr]; - } + for(int i = 0; i != Rows && primal != FALSE;) { + i++; + primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ? (short)1: (short)0; + } + if (primal == FALSE) { + drow[0] = 1; + for(int i = 1; i <= Rows; i++) drow[i] = 0; + Extrad = 0; + for(int i = 1; i <= columns; i++) { + varnr = Rows + i; + drow[varnr] = 0; + for(int j = Col_end[i - 1]; j < Col_end[i]; j++) + if (drow[get_row_nr(Mat, j)] != 0) + drow[varnr] += drow[get_row_nr(Mat, j)] * get_value(Mat,j); + if (drow[varnr] < Extrad) Extrad = drow[varnr]; } - else + } else { Extrad = 0; - if(Lp.trace != FALSE) - System.out.println("Extrad = " + Extrad); + } minit = FALSE; - - while(Status == RUNNING) - { - Doiter = FALSE; - DoInvert = FALSE; - - construct_solution(Solution); - if(primal != FALSE) - { - ref1.value = colnr; - flag = colprim(ref1, minit, drow); - colnr = (int)ref1.value; - if (flag != FALSE) - { - setpivcol(Lower[colnr], colnr, Pcol); - ref1.value = row_nr; - ref2.value = theta; - flag = rowprim(colnr, ref1, ref2, Pcol); - row_nr = (int)ref1.value; - theta = ref2.value; - if (flag != FALSE) - condensecol(row_nr, Pcol); - - } - } - else /* not primal */ - { - if(minit == FALSE) { - ref1.value = row_nr; - flag = rowdual(ref1); - row_nr = (int)ref1.value; - } - if(row_nr > 0 ) - { - ref1.value = colnr; - flag = coldual(row_nr, ref1, minit, prow, drow); - colnr = (int)ref1.value; - if (flag != FALSE) - { - setpivcol(Lower[colnr], colnr, Pcol); - /* getting div by zero here ... MB */ - if(Pcol[row_nr] == 0) - { - System.err.println("An attempt was made to divide by zero (Pcol[" + - row_nr + "])"); - System.err.println("This indicates numerical instability"); - Doiter = FALSE; - if(JustInverted == FALSE) - { - System.out.println("Reinverting Eta"); - DoInvert = TRUE; - } - else - { - System.out.println("Can't reinvert, failure"); - Status = FAILURE; - } - } - else - { - condensecol(row_nr, Pcol); - f = Rhs[row_nr] - Upbo[Bas[row_nr]]; - if(f > 0) - { - theta = f / (float) Pcol[row_nr]; - if(theta <= Upbo[colnr]) - Lower[Bas[row_nr]] = - (Lower[Bas[row_nr]] == FALSE)? - (short)1:(short)0; - } - else /* f <= 0 */ - theta = Rhs[row_nr] / (float) Pcol[row_nr]; - } - } - else - Status = INFEASIBLE; - } - else - { - primal = TRUE; - Doiter = FALSE; - Extrad = 0; - DoInvert = TRUE; - } - } - if(Doiter != FALSE) { - ref1.value = theta; - ref2.value = minit; - ref3.value = Lower[colnr]; - iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol); - theta = ref1.value; - minit = (short)ref2.value; - Lower[colnr] = (short)ref3.value; + while(Status == RUNNING) { + Doiter = FALSE; + DoInvert = FALSE; + construct_solution(Solution); + if (primal != FALSE) { + ref1.value = colnr; + flag = colprim(ref1, minit, drow); + colnr = (int)ref1.value; + if (flag != FALSE) { + setpivcol(Lower[colnr], colnr, Pcol); + ref1.value = row_nr; + ref2.value = theta; + flag = rowprim(colnr, ref1, ref2, Pcol); + row_nr = (int)ref1.value; + theta = ref2.value; + if (flag != FALSE) condensecol(row_nr, Pcol); } - if(Num_inv >= Lp.max_num_inv) + } else { + if (minit == FALSE) { + ref1.value = row_nr; + flag = rowdual(ref1); + row_nr = (int)ref1.value; + } + if (row_nr > 0) { + ref1.value = colnr; + flag = coldual(row_nr, ref1, minit, prow, drow); + colnr = (int)ref1.value; + if (flag != FALSE) { + setpivcol(Lower[colnr], colnr, Pcol); + /* getting div by zero here ... MB */ + if (Pcol[row_nr] == 0) { + throw new Error("An attempt was made to divide by zero (Pcol[" + row_nr + "])"); + } else { + condensecol(row_nr, Pcol); + f = Rhs[row_nr] - Upbo[Bas[row_nr]]; + if (f > 0) { + theta = f / (float) Pcol[row_nr]; + if (theta <= Upbo[colnr]) + Lower[Bas[row_nr]] = (Lower[Bas[row_nr]] == FALSE)? (short)1:(short)0; + } else theta = Rhs[row_nr] / (float) Pcol[row_nr]; + } + } else Status = INFEASIBLE; + } else { + primal = TRUE; + Doiter = FALSE; + Extrad = 0; DoInvert = TRUE; - if(DoInvert != FALSE) - { - invert(); - } - } - - Lp.total_iter += Lp.iter; - + } + } + if (Doiter != FALSE) { + ref1.value = theta; + ref2.value = minit; + ref3.value = Lower[colnr]; + iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol); + theta = ref1.value; + minit = (short)ref2.value; + Lower[colnr] = (short)ref3.value; + } + if (Num_inv >= max_num_inv) DoInvert = TRUE; + if (DoInvert != FALSE) invert(); + } + total_iter += iter; return(Status); - } /* solvelp */ - - - private short is_int(float value) - { - float tmp; - - tmp = (float)(value - Math.floor(value)); - if(tmp < Epsilon) - return(TRUE); - if(tmp > (1 - Epsilon)) - return(TRUE); - return(FALSE); - } /* is_int */ + } - private void construct_solution(float[] sol) - { - int i, j, basi; + private void construct_solution(float[] sol) { float f; + int basi; + for(int i = 0; i <= Rows; i++) sol[i] = 0; + for(int i = Rows + 1; i <= Sum; i++) sol[i] = Lowbo[i]; + for(int i = 1; i <= Rows; i++) { + basi = Bas[i]; + if (basi > Rows) sol[basi] += Rhs[i]; + } + for(int i = Rows + 1; i <= Sum; i++) + if (Basis[i] == FALSE && Lower[i] == FALSE) + sol[i] += Upbo[i]; + for(int j = 1; j <= columns; j++) { + f = sol[Rows + j]; + if (f != 0) + for(int i = Col_end[j - 1]; i < Col_end[j]; i++) + sol[get_row_nr(Mat, i)] += f * get_value(Mat,i); + } + for(int i = 0; i <= Rows; i++) { + if (Math.abs(sol[i]) < Epsb) sol[i] = 0; + else if (ch_sign[i] != FALSE) sol[i] = -sol[i]; + } + } - /* zero all results of rows */ - for (i = 0; i <= Rows; i++) - sol[i] = 0; - - if (Lp.scaling_used != FALSE) - { - for(i = Rows + 1; i <= Sum; i++) - sol[i] = Lowbo[i] * Lp.scale[i]; - for(i = 1; i <= Rows; i++) - { - basi = Bas[i]; - if(basi > Rows) - sol[basi] += Rhs[i] * Lp.scale[basi]; - } - for(i = Rows + 1; i <= Sum; i++) - if(Basis[i] == FALSE && Lower[i] == FALSE) - sol[i] += Upbo[i] * Lp.scale[i]; - - for(j = 1; j <= columns; j++) - { - f = sol[Rows + j]; - if(f != 0) - for(i = Col_end[j - 1]; i < Col_end[j]; i++) - sol[Mat[i].row_nr] += (f / Lp.scale[Rows+j]) - * (Mat[i].value / Lp.scale[Mat[i].row_nr]); - } - - for(i = 0; i <= Rows; i++) - { - if(Math.abs(sol[i]) < Epsb) - sol[i] = 0; - else - if(Lp.ch_sign[i] != FALSE) - sol[i] = -sol[i]; - } - } - else - { - for(i = Rows + 1; i <= Sum; i++) - sol[i] = Lowbo[i]; - for(i = 1; i <= Rows; i++) - { - basi = Bas[i]; - if(basi > Rows) - sol[basi] += Rhs[i]; - } - for(i = Rows + 1; i <= Sum; i++) - if(Basis[i] == FALSE && Lower[i] == FALSE) - sol[i] += Upbo[i]; - for(j = 1; j <= columns; j++) - { - f = sol[Rows + j]; - if(f != 0) - for(i = Col_end[j - 1]; i < Col_end[j]; i++) - sol[Mat[i].row_nr] += f * Mat[i].value; - } - - for(i = 0; i <= Rows; i++) - { - if(Math.abs(sol[i]) < Epsb) - sol[i] = 0; - else - if(Lp.ch_sign[i] != FALSE) - sol[i] = -sol[i]; - } - } - } /* construct_solution */ - - private void calculate_duals() - { - int i; - - /* initialise */ - for(i = 1; i <= Rows; i++) - Lp.duals[i] = 0; - Lp.duals[0] = 1; - btran(Lp.duals); - if (Lp.scaling_used != FALSE) - for(i = 1; i <= Rows; i++) - Lp.duals[i] *= Lp.scale[i]/Lp.scale[0]; - - /* the dual values are the reduced costs of the slacks */ - /* When the slack is at its upper Bound, change the sign. Can this happen? */ - for(i = 1; i <= Rows; i++) - { - if(Lp.basis[i] != FALSE) - Lp.duals[i] = 0; - else if( Lp.ch_sign[0] == Lp.ch_sign[i]) - Lp.duals[i] = -Lp.duals[i]; - } + private void calculate_duals() { + for(int i = 1; i <= Rows; i++) duals[i] = 0; + duals[0] = 1; + btran(duals); + for(int i = 1; i <= Rows; i++) { + if (basis[i] != FALSE) duals[i] = 0; + else if ( ch_sign[0] == ch_sign[i]) duals[i] = -duals[i]; + } } - private int milpsolve(float[] upbo, - float[] lowbo, - short[] sbasis, - short[] slower, - int[] sbas) - { - int i, j, failure, notint, is_worse; + private static Random rdm = new Random(); + + private int milpsolve(float[] upbo, float[] lowbo, short[] sbasis, short[] slower, int[] sbas) { + int failure, notint, is_worse; float theta, tmpfloat; - Random rdm = new Random(); notint = 0; - if(Break_bb != FALSE) - return(BREAK_BB); + if (Break_bb != FALSE) return(BREAK_BB); Level++; - Lp.total_nodes++; - if(Level > Lp.max_level) - Lp.max_level = Level; - /* make fresh copies of upbo, lowbo, rh as solving changes them */ + total_nodes++; + if (Level > max_level) max_level = Level; System.arraycopy(upbo, 0, Upbo, 0, Sum + 1); System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1); System.arraycopy(sbasis, 0, Basis, 0, Sum + 1); System.arraycopy(slower, 0, Lower, 0, Sum + 1); System.arraycopy(sbas, 0, Bas, 0, Rows + 1); System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1); - - if(Lp.anti_degen != FALSE) - { - for(i = 1; i <= columns; i++) - { - tmpfloat = rdm.nextFloat()*(float)0.001; - if(tmpfloat > Epsb) - Lowbo[i + Rows] -= tmpfloat; - tmpfloat = rdm.nextFloat()*(float)0.001; - if(tmpfloat > Epsb) - Upbo[i + Rows] += tmpfloat; - } - Lp.eta_valid = FALSE; - } - - if(Lp.eta_valid == FALSE) - { - for(i = 1; i <= columns; i++) - if(Lowbo[Rows + i] != 0) - { - theta = Lowbo[ Rows + i]; - if(Upbo[Rows + i]= Best_solution[0]) ? 1:0; - - if(is_worse != FALSE) - { - if(Lp.verbose != FALSE) - System.out.println("level" + Level + " OPT NOB value " + Solution[0] + - " Bound " + Best_solution[0]); - Level--; - return(MILP_FAIL); - } - - /* check if solution contains enough ints */ - if(Lp.bb_rule == FIRST_NI) - { - notint = 0; + if (failure == OPTIMAL) { + construct_solution(Solution); + /* if this solution is worse than the best sofar, this branch must die */ + if (Maximise != FALSE) is_worse = (Solution[0] <= Best_solution[0]) ? 1:0; + else is_worse = (Solution[0] >= Best_solution[0]) ? 1:0; + if (is_worse != FALSE) { + Level--; + return(MILP_FAIL); + } + /* check if solution contains enough ints */ + if (bb_rule == FIRST_NI) { + notint = 0; + int i = Rows + 1; + while(i <= Sum && notint == 0) i++; + } + if (bb_rule == RAND_NI) { + int nr_not_int, select_not_int; + nr_not_int = 0; + for(int i = Rows + 1; i <= Sum; i++) + if (nr_not_int == 0) notint = 0; + else { + select_not_int=(rdm.nextInt() % nr_not_int) + 1; i = Rows + 1; - while(i <= Sum && notint == 0) - { - i++; - } - } - if(Lp.bb_rule == RAND_NI) - { - int nr_not_int, select_not_int; - nr_not_int = 0; - for(i = Rows + 1; i <= Sum; i++) - if(nr_not_int == 0) - notint = 0; - else - { - select_not_int=(rdm.nextInt() % nr_not_int) + 1; - i = Rows + 1; - while(select_not_int > 0) - { - i++; - } - notint = i - 1; - } - } - - if(Lp.verbose == TRUE) - if(notint != FALSE) - System.out.println("level " + Level + " OPT value " + Solution[0]); - else - System.out.println("level " + Level + " OPT INT value " + Solution[0]); - - if(notint != FALSE) /* there is at least one value not yet int */ - { - /* set up two new problems */ - float[] new_upbo, new_lowbo; - float new_bound; - short[] new_lower, new_basis; - int[] new_bas; - int resone, restwo; - - /* allocate room for them */ - new_upbo = new float[Sum + 1]; - new_lowbo = new float[Sum + 1]; - new_lower = new short[Sum + 1]; - new_basis = new short[Sum + 1]; - new_bas = new int[Rows + 1]; - System.arraycopy(upbo, 0, new_upbo, 0, Sum + 1); - System.arraycopy(lowbo, 0, new_lowbo, 0, Sum + 1); - System.arraycopy(Lower, 0, new_lower, 0, Sum + 1); - System.arraycopy(Basis, 0, new_basis, 0, Sum + 1); - System.arraycopy(Bas, 0, new_bas, 0, Rows +1); - - if(Floor_first != FALSE) - { - new_bound = (float)(Math.ceil(Solution[notint]) - 1); - /* this Bound might conflict */ - if(new_bound < lowbo[notint]) - { - resone = MILP_FAIL; - } - else /* Bound feasible */ - { - new_upbo[notint] = new_bound; - Lp.eta_valid = FALSE; - resone = milpsolve(new_upbo, lowbo, new_basis, new_lower, - new_bas); - Lp.eta_valid = FALSE; - } - new_bound += 1; - if(new_bound > upbo[notint]) - { - restwo = MILP_FAIL; - } - else /* Bound feasible */ - { - new_lowbo[notint] = new_bound; - Lp.eta_valid = FALSE; - restwo = milpsolve(upbo, new_lowbo, new_basis, new_lower, - new_bas); - Lp.eta_valid = FALSE; - } - } - else /* take ceiling first */ - { - new_bound = (float)Math.ceil(Solution[notint]); - /* this Bound might conflict */ - if(new_bound > upbo[notint]) - { - resone = MILP_FAIL; - } - else /* Bound feasible */ - { - new_lowbo[notint] = new_bound; - Lp.eta_valid = FALSE; - resone = milpsolve(upbo, new_lowbo, new_basis, new_lower, - new_bas); - Lp.eta_valid = FALSE; - } - new_bound -= 1; - if(new_bound < lowbo[notint]) - { - restwo = MILP_FAIL; - } - else /* Bound feasible */ - { - new_upbo[notint] = new_bound; - Lp.eta_valid = FALSE; - restwo = milpsolve(new_upbo, lowbo, new_basis, new_lower, - new_bas); - Lp.eta_valid = FALSE; - } - } - if(resone != FALSE && restwo != FALSE) - /* both failed and must have been infeasible */ - failure = INFEASIBLE; - else - failure = OPTIMAL; - - } - else /* all required values are int */ - { - if(Maximise != FALSE) - is_worse = (Solution[0] < Best_solution[0]) ? 1:0; - else - is_worse = (Solution[0] > Best_solution[0]) ? 1:0; - - if(is_worse == FALSE) /* Current solution better */ - { - System.arraycopy(Solution, 0, Best_solution, 0, Sum + 1); - calculate_duals(); - if(Lp.break_at_int != FALSE) - { - if(Maximise != FALSE && (Best_solution[0] > Lp.break_value)) - Break_bb = TRUE; - if(Maximise == FALSE && (Best_solution[0] < Lp.break_value)) - Break_bb = TRUE; - } - } + while(select_not_int > 0) i++; + notint = i - 1; } } - - /* failure can have the values OPTIMAL, UNBOUNDED and INFEASIBLE. */ + if (notint != FALSE) throw new Error("integer linear programming not supported"); + if (Maximise != FALSE) is_worse = (Solution[0] < Best_solution[0]) ? 1:0; + else is_worse = (Solution[0] > Best_solution[0]) ? 1:0; + if (is_worse == FALSE) { + System.arraycopy(Solution, 0, Best_solution, 0, Sum + 1); + calculate_duals(); + if (break_at_int != FALSE) { + if (Maximise != FALSE && (Best_solution[0] > break_value)) Break_bb = TRUE; + if (Maximise == FALSE && (Best_solution[0] < break_value)) Break_bb = TRUE; + } + } + } Level--; return(failure); - } /* milpsolve */ - - public int solve(Problem lp) - { - int result, i; - - if(lp.active == FALSE) - set_globals(lp); - - lp.total_iter = 0; - lp.max_level = 1; - lp.total_nodes = 0; - - if(Isvalid(lp) != FALSE) - { - if(Maximise != FALSE && lp.obj_bound == Infinite) - Best_solution[0]=-Infinite; - else if(Maximise == FALSE && lp.obj_bound==-Infinite) - Best_solution[0] = Infinite; - else - Best_solution[0] = lp.obj_bound; - - Level = 0; - - if(lp.basis_valid == FALSE) - { - for(i = 0; i <= lp.rows; i++) - { - lp.basis[i] = TRUE; - lp.bas[i] = i; - } - for(i = lp.rows+1; i <= lp.sum; i++) - lp.basis[i] = FALSE; - for(i = 0; i <= lp.sum; i++) - lp.lower[i] = TRUE; - lp.basis_valid = TRUE; - } - - lp.eta_valid = FALSE; - Break_bb = FALSE; - result = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas); - lp.eta_size = Eta_size; - lp.eta_alloc = Eta_alloc; - lp.num_inv = Num_inv; - - return(result); - } - - /* if we get here, Isvalid(lp) failed. I suggest we return FAILURE */ - return(FAILURE); - } - - public int lag_solve(Problem lp, float start_bound, int num_iter, short verbose) - { - int i, j, result, citer; - short status, OrigFeas, AnyFeas, same_basis; - float[] OrigObj, ModObj, SubGrad, BestFeasSol; - float Zub, Zlb, Ztmp, pie; - float rhsmod, Step, SqrsumSubGrad; - int[] old_bas; - short[] old_lower; - - /* allocate mem */ - OrigObj = new float[lp.columns + 1]; - ModObj = new float[lp.columns + 1]; - for (i = 0; i <= lp.columns; i++) - ModObj[i] = 0; - SubGrad = new float[lp.nr_lagrange]; - for (i = 0; i < lp.nr_lagrange; i++) - SubGrad[i] = 0; - BestFeasSol = new float[lp.sum + 1]; - for (i = 0; i <= lp.sum; i++) - BestFeasSol[i] = 0; - old_bas = new int[lp.rows + 1]; - System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows + 1); - old_lower = new short[lp.sum + 1]; - System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum + 1); - - get_row(lp, 0, OrigObj); - - pie = 2; - - if(lp.maximise != FALSE) - { - Zub = DEF_INFINITE; - Zlb = start_bound; - } - else - { - Zlb = -DEF_INFINITE; - Zub = start_bound; - } - status = RUNNING; - Step = 1; - OrigFeas = FALSE; - AnyFeas = FALSE; - citer = 0; - - for(i = 0 ; i < lp.nr_lagrange; i++) - lp.lambda[i] = 0; - - while(status == RUNNING) - { - citer++; - - for(i = 1; i <= lp.columns; i++) - { - ModObj[i] = OrigObj[i]; - for(j = 0; j < lp.nr_lagrange; j++) - { - if(lp.maximise != FALSE) - ModObj[i] -= lp.lambda[j] * lp.lag_row[j][i]; - else - ModObj[i] += lp.lambda[j] * lp.lag_row[j][i]; - } - } - for(i = 1; i <= lp.columns; i++) - { - set_mat(lp, 0, i, ModObj[i]); - /* fSystem.out.print(stderr, "%f ", ModObj[i]); */ - } - rhsmod = 0; - for(i = 0; i < lp.nr_lagrange; i++) - if(lp.maximise != FALSE) - rhsmod += lp.lambda[i] * lp.lag_rhs[i]; - else - rhsmod -= lp.lambda[i] * lp.lag_rhs[i]; - - if(verbose != FALSE) - { - System.out.println("Zub: " + Zub + " Zlb: " + Zlb + " Step: " + Step + - " pie: " + pie + " Feas " + OrigFeas); - for(i = 0; i < lp.nr_lagrange; i++) - System.out.println(i + " SubGrad " + SubGrad[i] + " lambda " + lp.lambda[i]); - } - - - result = solve(lp); - - same_basis = TRUE; - i = 1; - while(same_basis != FALSE && i < lp.rows) - { - same_basis = (old_bas[i] == lp.bas[i])? (short)1: (short)0; - i++; - } - i = 1; - while(same_basis != FALSE && i < lp.sum) - { - same_basis=(old_lower[i] == lp.lower[i])? (short)1:(short)0; - i++; - } - if(same_basis == FALSE) - { - System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum+1); - System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows+1); - pie *= 0.95; - } - - if(verbose != FALSE) - System.out.println("result: " + result + " same basis: " + same_basis); - - if(result == UNBOUNDED) - { - for(i = 1; i <= lp.columns; i++) - System.out.print(ModObj[i] + " "); - System.exit(FAIL); - } - - if(result == FAILURE) - status = FAILURE; - - if(result == INFEASIBLE) - status = INFEASIBLE; - - SqrsumSubGrad = 0; - for(i = 0; i < lp.nr_lagrange; i++) - { - SubGrad[i]= -lp.lag_rhs[i]; - for(j = 1; j <= lp.columns; j++) - SubGrad[i] += lp.best_solution[lp.rows + j] * lp.lag_row[i][j]; - SqrsumSubGrad += SubGrad[i] * SubGrad[i]; - } - - OrigFeas = TRUE; - for(i = 0; i < lp.nr_lagrange; i++) - if(lp.lag_con_type[i] != FALSE) - { - if(Math.abs(SubGrad[i]) > lp.epsb) - OrigFeas = FALSE; - } - else if(SubGrad[i] > lp.epsb) - OrigFeas = FALSE; - - if(OrigFeas != FALSE) - { - AnyFeas = TRUE; - Ztmp = 0; - for(i = 1; i <= lp.columns; i++) - Ztmp += lp.best_solution[lp.rows + i] * OrigObj[i]; - if((lp.maximise != FALSE) && (Ztmp > Zlb)) - { - Zlb = Ztmp; - for(i = 1; i <= lp.sum; i++) - BestFeasSol[i] = lp.best_solution[i]; - BestFeasSol[0] = Zlb; - if(verbose != FALSE) - System.out.println("Best feasible solution: " + Zlb); - } - else if(Ztmp < Zub) - { - Zub = Ztmp; - for(i = 1; i <= lp.sum; i++) - BestFeasSol[i] = lp.best_solution[i]; - BestFeasSol[0] = Zub; - if(verbose != FALSE) - System.out.println("Best feasible solution: " + Zub); - } - } - - if(lp.maximise != FALSE) - Zub = Math.min(Zub, rhsmod + lp.best_solution[0]); - else - Zlb = Math.max(Zlb, rhsmod + lp.best_solution[0]); - - if(Math.abs(Zub-Zlb) < (float)0.001) - { - status = OPTIMAL; - } - Step = (float)(pie * ((1.05*Zub) - Zlb) / SqrsumSubGrad); - - for(i = 0; i < lp.nr_lagrange; i++) - { - lp.lambda[i] += Step * SubGrad[i]; - if(lp.lag_con_type[i] == FALSE && lp.lambda[i] < 0) - lp.lambda[i] = 0; - } - - if(citer == num_iter && status==RUNNING) - if(AnyFeas != FALSE) - status = FEAS_FOUND; - else - status = NO_FEAS_FOUND; - } - - for(i = 0; i <= lp.sum; i++) - lp.best_solution[i] = BestFeasSol[i]; - - for(i = 1; i <= lp.columns; i++) - set_mat(lp, 0, i, OrigObj[i]); - - if(lp.maximise != FALSE) - lp.lag_bound = Zub; - else - lp.lag_bound = Zlb; - return(status); - } - - } // end of class solve - - - public static class Matrix { - public int row_nr; - public float value; - public Matrix(int r, float v) { row_nr = r; value = v; } - } - - - public static class Problem { - String lp_name; /* the name of the lp */ - - public short active; /*TRUE if the globals point to this structure*/ - public short verbose; /* ## Verbose flag */ - public short debug; /* ## Print B&B information */ - public short trace; /* ## Print information on pivot selection */ - public short anti_degen; /* ## Do perturbations */ - - public int rows; /* Nr of constraint rows in the problem */ - int rows_alloc; /* The allocated memory for Rows sized data */ - int columns; /* The number of columns (= variables) */ - int columns_alloc; - int sum; /* The size of the variables + the slacks */ - int sum_alloc; - - short names_used; /* Flag to indecate if names for rows and - columns are used */ - String[] row_name; /* rows_alloc+1 */ - int[] col_name; /* columns_alloc+1 */ - - /* Row[0] of the sparce matrix is the objective function */ - - int non_zeros; /* The number of elements in the sparce matrix*/ - int mat_alloc; /* The allocated size for matrix sized - structures */ - Matrix[] mat; /* mat_alloc :The sparse matrix */ - Matrix[] alternate_mat; /* mat_alloc :The sparse matrix */ - int[] col_end; /* columns_alloc+1 :Cend[i] is the index of the - first element after column i. - column[i] is stored in elements - col_end[i-1] to col_end[i]-1 */ - int[] col_no; /* mat_alloc :From Row 1 on, col_no contains the - column nr. of the - nonzero elements, row by row */ - short row_end_valid; /* true if row_end & col_no are valid */ - int[] row_end; /* rows_alloc+1 :row_end[i] is the index of the - first element in Colno after row i */ - float[] orig_rh; /* rows_alloc+1 :The RHS after scaling & sign - changing, but before `Bound transformation' */ - float[] rh; /* rows_alloc+1 :As orig_rh, but after Bound - transformation */ - float[] rhs; /* rows_alloc+1 :The RHS of the curent simplex - tableau */ - float[] orig_upbo; /* sum_alloc+1 :Bound before transformations */ - float[] orig_lowbo; /* " " */ - float[] upbo; /* " " :Upper bound after transformation - & B&B work*/ - float[] lowbo; /* " " :Lower bound after transformation - & B&B work */ - - short basis_valid; /* TRUE is the basis is still valid */ - int[] bas; /* rows_alloc+1 :The basis column list */ - short[] basis; /* sum_alloc+1 : basis[i] is TRUE if the column - is in the basis */ - short[] lower; /* " " :TRUE is the variable is at its - lower bound (or in the basis), it is FALSE - if the variable is at its upper bound */ - - short eta_valid; /* TRUE if current Eta structures are valid */ - int eta_alloc; /* The allocated memory for Eta */ - int eta_size; /* The number of Eta columns */ - int num_inv; /* The number of float pivots */ - int max_num_inv; /* ## The number of float pivots between - reinvertions */ - float[] eta_value; /* eta_alloc :The Structure containing the - values of Eta */ - int[] eta_row_nr; /* " " :The Structure containing the Row - indexes of Eta */ - int[] eta_col_end; /* rows_alloc + MaxNumInv : eta_col_end[i] is - the start index of the next Eta column */ - - short bb_rule; /* what rule for selecting B&B variables */ - - short break_at_int; /* TRUE if stop at first integer better than - break_value */ - float break_value; - - float obj_bound; /* ## Objective function bound for speedup of - B&B */ - int iter; /* The number of iterations in the simplex - solver (LP) */ - int total_iter; /* The total number of iterations (B&B) (ILP)*/ - int max_level; /* The Deepest B&B level of the last solution */ - int total_nodes; /* total number of nodes processed in b&b */ - public float[] solution; /* sum_alloc+1 :The Solution of the last LP, - 0 = The Optimal Value, - 1..rows The Slacks, - rows+1..sum The Variables */ - public float[] best_solution; /* " " :The Best 'Integer' Solution */ - float[] duals; /* rows_alloc+1 :The dual variables of the - last LP */ - - short maximise; /* TRUE if the goal is to maximise the - objective function */ - short floor_first; /* TRUE if B&B does floor bound first */ - short[] ch_sign; /* rows_alloc+1 :TRUE if the Row in the matrix - has changed sign - (a`x > b, x>=0) is translated to - s + -a`x = -b with x>=0, s>=0) */ - - short scaling_used; /* TRUE if scaling is used */ - short columns_scaled; /* TRUE is the columns are scaled too, Only use - if all variables are non-integer */ - float[] scale; /* sum_alloc+1 :0..Rows the scaling of the Rows, - Rows+1..Sum the scaling of the columns */ - - int nr_lagrange; /* Nr. of Langrangian relaxation constraints */ - float[][]lag_row; /* NumLagrange, columns+1:Pointer to pointer of - rows */ - float[] lag_rhs; /* NumLagrange :Pointer to pointer of Rhs */ - float[] lambda; /* NumLagrange :Lambda Values */ - short[] lag_con_type; /* NumLagrange :TRUE if constraint type EQ */ - float lag_bound; /* the lagrangian lower bound */ - - short valid; /* Has this lp pased the 'test' */ - float infinite; /* ## numercal stuff */ - float epsilon; /* ## */ - float epsb; /* ## */ - float epsd; /* ## */ - float epsel; /* ## */ - - - public Problem (int nrows, int ncolumns, int matalloc) { - int i, nsum; - nsum=nrows+ncolumns; - rows=nrows; - columns=ncolumns; - sum=nsum; - rows_alloc=rows; - columns_alloc=columns; - sum_alloc=sum; - mat_alloc=matalloc; - eta_alloc=10000; - max_num_inv=DEFNUMINV; - col_no = new int[mat_alloc]; - col_end = new int[columns + 1]; - row_end = new int[rows + 1]; - orig_rh = new float[rows + 1]; - rh = new float[rows + 1]; - rhs = new float[rows + 1]; - orig_upbo = new float[sum + 1]; - upbo = new float[sum + 1]; - orig_lowbo = new float[sum + 1]; - lowbo = new float[sum + 1]; - bas = new int[rows+1]; - basis = new short[sum + 1]; - lower = new short[sum + 1]; - eta_value = new float[eta_alloc]; - eta_row_nr = new int[eta_alloc]; - eta_col_end = new int[rows_alloc + max_num_inv]; - solution = new float[sum + 1]; - best_solution = new float[sum + 1]; - duals = new float[rows + 1]; - ch_sign = new short[rows + 1]; - mat = new Matrix[mat_alloc]; - for(int j=0; j