Renamed LinearProgramming to Simplex
[org.ibex.core.git] / src / org / ibex / util / LinearProgramming.java
diff --git a/src/org/ibex/util/LinearProgramming.java b/src/org/ibex/util/LinearProgramming.java
deleted file mode 100644 (file)
index 2dc47d4..0000000
+++ /dev/null
@@ -1,1342 +0,0 @@
-package org.ibex.util;
-import java.io.* ;
-import java.util.* ;
-
-public class LinearProgramming {
-
-    public final static short FAIL = -1;
-    
-    public final static short NULL = 0;
-    public final static short FALSE = 0;
-    public final static short TRUE = 1;
-    
-    public final static short DEFNUMINV = 50;
-    
-    /* solve status values */
-    public final static short OPTIMAL = 0;
-    public final static short MILP_FAIL = 1;
-    public final static short INFEASIBLE = 2;
-    public final static short UNBOUNDED = 3;
-    public final static short FAILURE = 4;
-    public final static short RUNNING = 5;
-    
-    /* lag_solve extra status values */
-    public final static short FEAS_FOUND = 6;
-    public final static short NO_FEAS_FOUND = 7;
-    public final static short BREAK_BB = 8;
-    
-    public final static short FIRST_NI =       0;
-    public final static short RAND_NI = 1;
-    
-    public final static short LE = 0;
-    public final static short EQ = 1;
-    public final static short GE = 2;
-    public final static short OF = 3;
-    
-    public final static short MAX_WARN_COUNT = 20;
-    
-    public final static float DEF_INFINITE = (float)1e24; /* limit for dynamic range */
-    public final static float DEF_EPSB = (float)5.01e-7; /* for rounding RHS values to 0 determine     
-                                               infeasibility basis */
-    public final static float DEF_EPSEL = (float)1e-8; /* for rounding other values (vectors) to 0 */
-    public final static float DEF_EPSD  = (float)1e-6; /* for rounding reduced costs to zero */
-    public final static float DEF_EPSILON = (float)1e-3; /* to determine if a float value is integer */
-    
-    public final static float PREJ = (float)1e-3;  /* pivot reject (try others first) */
-    
-    public final static int ETA_START_SIZE = 10000; /* start size of array Eta. Realloced if needed */
-
-    static class Ref {
-        float value;
-        public Ref(float v) { value = v; }
-    }
-
-    public static class Simplex {
-        /* Globals used by solver */
-        short JustInverted;
-        short Status;
-        short Doiter;
-        short DoInvert;
-        short Break_bb;
-
-        public short active;           /*TRUE if the globals point to this structure*/
-        public short debug;           /* ## Print B&B information */
-        public short trace;           /* ## Print information on pivot selection */
-        public int         rows;               /* Nr of constraint rows in the problem */
-        int       rows_alloc;          /* The allocated memory for Rows sized data */
-        int       columns_alloc;  
-        int       sum;                /* The size of the variables + the slacks */
-        int       sum_alloc;
-        int       non_zeros;          /* The number of elements in the sparce matrix*/
-        int       mat_alloc;           /* The allocated size for matrix sized 
-                                           structures */
-        MatrixArray  mat;                /* mat_alloc :The sparse matrix */
-        MatrixArray  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 () */
-        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) */ 
-
-        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;
-        MatrixArray  Mat;
-        int[]     Col_no;
-        int[]     Col_end;
-        int[]     Row_end;
-        float[]    Orig_rh;
-        float[]    Rh;
-        float[]    Rhs;
-        float[]    Orig_upbo;
-        float[]    Orig_lowbo;
-        float[]    Upbo;
-        float[]    Lowbo;
-        int[]     Bas;
-        short[]   Basis;
-        short[]   Lower;
-        int     Eta_alloc; 
-        int     Eta_size;           
-        float[]    Eta_value;
-        int[]     Eta_row_nr;
-        int[]     Eta_col_end;
-        int     Num_inv;
-        float[]    Solution;
-        float[]    Best_solution;
-        float    Infinite;
-        float    Epsilon;
-        float    Epsb;
-        float    Epsd;
-        float    Epsel;
-  
-        float  TREJ;
-        float  TINV;
-  
-        short   Maximise;
-        short   Floor_first;
-        float    Extrad;
-
-        int     Warn_count; /* used in CHECK version of rounding macro */
-
-        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 ncolumns) {
-            int nsum;  
-            int nrows = 0;
-            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;
-
-            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;
-
-            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(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;
-
-            MatrixArray newmat;
-            int  elmnr;
-            int  stcol;
-
-            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(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];
-            }
-
-            for(int i =  1 ; i <= rows; i++) if(bas[i] >= rows) bas[i]++;
-
-            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) orig_rh[rows]=-rh;
-            else orig_rh[rows]=rh;  
-
-            row_end_valid=FALSE;
-            bas[rows]=rows;
-            basis[rows]=TRUE;
-            lower[rows]=TRUE;   
-            if (active != FALSE) set_globals();
-            eta_valid=FALSE;
-        }
-
-        public void bound_sum(int column1, int column2, float bound, short type, float[] scratch) {
-            for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
-            scratch[column1] = (float)1.0;
-            scratch[column2] = (float)1.0;
-            add_constraint(scratch, type, bound);
-            for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
-        }
-
-        public void bound_difference(int column1, int column2, float bound, short type, float[] scratch) {
-            for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
-            scratch[column1] = (float)1.0;
-            scratch[column2] = (float)-1.0;
-            add_constraint(scratch, type, bound);
-            for(int i=0; i<scratch.length; i++) scratch[i] = (float)0.0;
-        }
-
-        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(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(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(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_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");
-            }
-        }
-
-        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;
-        }
-
-        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);
-        }
-
-        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;
-            }
-        }
-
-        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;
-                        }
-                    }
-                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);
-        } 
-
-        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);
-            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);
-            eta_row_nr = Eta_row_nr;
-            */
-        }
-
-        private void condensecol(int row_nr, float[] pcol) {
-            int elnr;
-            elnr = Eta_col_end[Eta_size];
-            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;
-        }
-
-        private void addetacol() {
-            int k;
-            float theta;
-            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(int i = j; i < k - 1; i++) Eta_value[i] *= -theta;
-            JustInverted = FALSE;
-        }
-
-        private void setpivcol(short lower,  int varin, float[]   pcol) {
-            int colnr;
-            float f;
-            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);
-        }
-
-        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;
-                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);
-                }
-            }
-            Eta_row_nr[elnr] = row_nr;
-            Eta_value[elnr] = 1 / (float) piv;
-            elnr++;
-            theta = Rhs[row_nr] / (float) piv;
-            Rhs[row_nr] = theta;
-            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(int i = wk; i < elnr - 1; i++) Eta_value[i] /= - (float) piv;
-            Eta_col_end[Eta_size] = elnr;
-        }
-
-        private void rhsmincol(float theta, int row_nr, int varin) {
-            int varout;
-            float 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;
-        }
-
-        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(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) {
-                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) {
-                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;
-                        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(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, Rows + j);
-                    addetacol();
-                }
-            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;
-        }
-
-        private short colprim(Ref colnr, short minit, float[]   drow) {
-            int  varnr;
-            float f, dpiv;
-              dpiv = -Epsd;
-            colnr.value = 0;
-            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(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 (colnr.value == 0) {
-                Doiter   = FALSE;
-                DoInvert = FALSE;
-                Status   = OPTIMAL;
-            }
-            return(colnr.value > 0 ? (short)1 : (short)0);
-        }
-
-        private short rowprim(int colnr, Ref row_nr, Ref theta, float[] pcol) {
-            float f = 0, quot; 
-            row_nr.value = 0;
-            theta.value = Infinite;
-            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 (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 (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);
-        }
-
-        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;
-                    }
-                }
-            }
-            return(row_nr.value > 0 ? (short)1 : (short)0);
-        }
-
-        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(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;
-                    }
-                }
-                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(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 (colnr.value > 0) Doiter = TRUE;
-            return(colnr.value > 0 ? (short)1 : (short)0);
-        }
-
-        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;
-            iter++;
-            minit.value = theta.value > (up + Epsb) ? 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(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    varnr;
-            float   f = 0, theta = 0;
-            short  primal;
-            short  minit;
-            int    colnr, row_nr;
-            colnr = 0;
-            row_nr = 0;
-            short flag; 
-            Ref ref1, ref2, ref3;
-            ref1 = new Ref(0);
-            ref2 = new Ref(0);
-            ref3 = new Ref(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;
-            primal = TRUE;
-            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 {
-                Extrad = 0;
-            }
-            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 {
-                    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 (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);
-        }
-
-        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];
-            }
-        }
-
-        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 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;
-            notint = 0;
-
-            if (Break_bb != FALSE) return(BREAK_BB);
-            Level++;
-            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 (eta_valid == FALSE) {
-                for(int i = 1; i <= columns; i++)
-                    if (Lowbo[Rows + i] != 0) {
-                        theta = Lowbo[ Rows + i];
-                        if (Upbo[Rows + i]<Infinite) Upbo[Rows + i] -= theta;
-                        for(int j = Col_end[i - 1]; j < Col_end[i]; j++) Rh[get_row_nr(Mat, j)] -= theta * get_value(Mat,j);
-                    }
-                invert();
-                eta_valid = TRUE;
-            }
-            failure = solvelp();
-            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(select_not_int > 0) i++;
-                            notint = i - 1;
-                        }
-                }
-                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);
-        }
-
-        public int solve() {
-            int result;
-            if (active == FALSE) set_globals();
-            total_iter  = 0;
-            max_level   = 1;
-            total_nodes = 0;
-            if (Isvalid() != FALSE) {
-                if (Maximise != FALSE && obj_bound == Infinite) Best_solution[0]=-Infinite;
-                else if (Maximise == FALSE && obj_bound==-Infinite) Best_solution[0] = Infinite;
-                else Best_solution[0] = obj_bound;
-                Level = 0;
-                if (basis_valid == FALSE) {
-                    for(int i = 0; i <= rows; i++) {
-                        basis[i] = TRUE;
-                        bas[i] = i;
-                    }
-                    for(int i = rows+1; i <= sum; i++) basis[i] = FALSE;
-                    for(int i = 0; i <= sum; i++) lower[i] = TRUE;
-                    basis_valid = TRUE;
-                }
-                eta_valid = FALSE;
-                Break_bb      = FALSE;
-                result        = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas); 
-                eta_size  = Eta_size;
-                eta_alloc = Eta_alloc;
-                num_inv   = Num_inv;
-                return(result);
-            }
-            return(FAILURE);
-        }
-    }
-
-    private final static float round( float val, float eps) { return (Math.abs(val) < eps) ? 0 : val; }
-    static int get_row_nr(MatrixArray m, int i) { return m.row_nr[i]; }
-    static void set_row_nr(MatrixArray m, int i, int val) { m.row_nr[i] = val; }
-    static float get_value(MatrixArray m, int i) { return m.value[i]; }
-    static void set_value(MatrixArray m, int i, float val) { m.value[i] = val; }
-    public static class MatrixArray {
-        public int[] row_nr;
-        public float[] value;
-        public final int length;
-        public MatrixArray(int length) { row_nr = new int[length]; value = new float[length]; this.length = length; }
-    }
-
-}
-