performance boost for simplex solver
[org.ibex.core.git] / src / org / ibex / util / LinearProgramming.java
index 869737f..d6b6176 100644 (file)
@@ -98,177 +98,6 @@ public class LinearProgramming {
 
         int     Warn_count; /* used in CHECK version of rounding macro */
 
-        void inc_mat_space(Problem lp, int maxextra)
-        {
-            if(lp.non_zeros + maxextra > lp.mat_alloc)
-                {
-                    int i;
-                    lp.mat_alloc = lp.non_zeros + maxextra;
-                    Matrix[] old_mat = lp.mat;
-                    lp.mat = new Matrix[lp.mat_alloc];
-                    for (i = old_mat.length; i < lp.mat.length; i++)
-                        lp.mat[i] = new Matrix(0,0);
-                    System.arraycopy(old_mat, 0, lp.mat, 0, old_mat.length);
-                    int[] old_col_no = lp.col_no;
-                    lp.col_no = new int[lp.mat_alloc];
-                    System.arraycopy(old_col_no, 0, lp.col_no, 0, old_col_no.length);
-                    if (lp.active != FALSE)
-                        {
-                            Mat=lp.mat;
-                            Col_no=lp.col_no;
-                        }
-                }
-        } // end of inc_mat_space
-        void inc_row_space(Problem lp)
-        {
-            if(lp.rows > lp.rows_alloc)
-                {
-                    lp.rows_alloc=lp.rows+10;
-                    lp.sum_alloc=lp.rows_alloc+lp.columns_alloc;
-
-                    float[] db_ptr = lp.orig_rh;
-                    lp.orig_rh = new float[lp.rows_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.orig_rh, 0, db_ptr.length);
-      
-                    db_ptr = lp.rh;
-                    lp.rh = new float[lp.rows_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.rh, 0, db_ptr.length);
-
-                    db_ptr = lp.rhs;
-                    lp.rhs = new float[lp.rows_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.rhs, 0, db_ptr.length);
-
-                    db_ptr = lp.orig_upbo;
-                    lp.orig_upbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.orig_upbo, 0, db_ptr.length);
-
-                    db_ptr = lp.upbo;
-                    lp.upbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.upbo, 0, db_ptr.length);
-
-                    db_ptr = lp.orig_lowbo;
-                    lp.orig_lowbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.orig_lowbo, 0, db_ptr.length);
-
-                    db_ptr = lp.lowbo;
-                    lp.lowbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.lowbo, 0, db_ptr.length);
-
-                    db_ptr = lp.solution;
-                    lp.solution = new float[lp.sum_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.solution, 0, db_ptr.length);
-
-                    db_ptr = lp.best_solution;
-                    lp.best_solution = new float[lp.sum_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.best_solution, 0, db_ptr.length);
-
-                    int[] int_ptr = lp.row_end;
-                    lp.row_end = new int[lp.rows_alloc + 1];
-                    System.arraycopy(int_ptr, 0, lp.row_end, 0, int_ptr.length);
-
-                    short[] short_ptr = lp.basis;
-                    lp.basis = new short[lp.sum_alloc + 1];
-                    System.arraycopy(short_ptr, 0, lp.basis, 0, short_ptr.length);
-
-                    short_ptr = lp.lower;
-                    lp.lower = new short[lp.sum_alloc + 1];
-                    System.arraycopy(short_ptr, 0, lp.lower, 0, short_ptr.length);
-
-                    int_ptr = lp.bas;
-                    lp.bas = new int[lp.rows_alloc + 1];
-                    System.arraycopy(int_ptr, 0, lp.bas, 0, int_ptr.length);
-
-                    db_ptr = lp.duals;
-                    lp.duals = new float[lp.rows_alloc + 1];
-                    System.arraycopy(db_ptr, 0, lp.duals, 0, db_ptr.length);
-
-                    short_ptr = lp.ch_sign;
-                    lp.ch_sign = new short[lp.rows_alloc + 1];
-                    System.arraycopy(short_ptr, 0, lp.ch_sign, 0, short_ptr.length);
-
-                    int_ptr = lp.eta_col_end;
-                    lp.eta_col_end = new int[lp.rows_alloc + lp.max_num_inv];
-                    System.arraycopy(int_ptr, 0, lp.eta_col_end, 0, int_ptr.length);
-
-                    if(lp.names_used != FALSE) {
-                        String[] str_ptr = lp.row_name;
-                        lp.row_name = new String[lp.rows_alloc + 1];
-                        System.arraycopy(str_ptr, 0, lp.row_name, 0, str_ptr.length);
-                    }
-
-                    if(lp.scaling_used != FALSE) {
-                        db_ptr = lp.scale;
-                        lp.scale = new float[lp.sum_alloc + 1];
-                        System.arraycopy(db_ptr, 0, lp.scale, 0, db_ptr.length);
-                    }
-
-                    if(lp.active != FALSE)
-                        set_globals(lp); 
-                }
-        }
-
-        void inc_col_space(Problem lp)
-        {
-            if(lp.columns >= lp.columns_alloc)
-                {
-                    int[] int_ptr;
-                    lp.columns_alloc=lp.columns+10;
-                    lp.sum_alloc=lp.rows_alloc+lp.columns_alloc;
-
-                    float[] float_ptr = lp.orig_upbo;
-                    lp.orig_upbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(float_ptr, 0, lp.orig_upbo, 0, float_ptr.length);
-
-                    float_ptr = lp.upbo;
-                    lp.upbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(float_ptr, 0, lp.upbo, 0, float_ptr.length);
-
-                    float_ptr = lp.orig_lowbo;
-                    lp.orig_lowbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(float_ptr, 0, lp.orig_lowbo, 0, float_ptr.length);
-
-                    float_ptr = lp.lowbo;
-                    lp.lowbo = new float[lp.sum_alloc + 1];
-                    System.arraycopy(float_ptr, 0, lp.lowbo, 0, float_ptr.length);
-
-                    float_ptr = lp.solution;
-                    lp.solution = new float[lp.sum_alloc + 1];
-                    System.arraycopy(float_ptr, 0, lp.solution, 0, float_ptr.length);
-
-                    float_ptr = lp.best_solution;
-                    lp.best_solution = new float[lp.sum_alloc + 1];
-                    System.arraycopy(float_ptr, 0, lp.best_solution, 0, float_ptr.length);
-
-                    short[] short_ptr = lp.basis;
-                    lp.basis = new short[lp.sum_alloc + 1];
-                    System.arraycopy(short_ptr, 0, lp.basis, 0, short_ptr.length);
-
-                    short_ptr = lp.lower;
-                    lp.lower = new short[lp.sum_alloc + 1];
-                    System.arraycopy(short_ptr, 0, lp.lower, 0, short_ptr.length);
-
-                    if(lp.names_used != FALSE) {
-                        int[] str_ptr = lp.col_name;
-                        lp.col_name = new int[lp.columns_alloc + 1];
-                        System.arraycopy(str_ptr, 0, lp.col_name, 0, str_ptr.length);
-                    }
-
-                    if(lp.scaling_used != FALSE) {
-                        float_ptr = lp.scale;
-                        lp.scale = new float[lp.sum_alloc + 1];
-                        System.arraycopy(float_ptr, 0, lp.scale, 0, float_ptr.length);
-                    }
-
-                    int_ptr = lp.col_end;
-                    lp.col_end = new int[lp.columns_alloc + 1];
-                    System.arraycopy(int_ptr, 0, lp.col_end, 0, int_ptr.length);
-
-                    if(lp.active != FALSE)
-                        set_globals(lp);
-                }
-        } // end of inc_col_space
-
         public void set_mat(Problem lp, int Row, int column, float Value)
         {
             int elmnr, lastelm, i;
@@ -312,7 +141,7 @@ public class LinearProgramming {
                     else
                         {
                             /* check if more space is needed for matrix */
-                            inc_mat_space(lp,1);
+                            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; 
@@ -380,19 +209,22 @@ public class LinearProgramming {
                     }
                 else
                     addtoo[i]=FALSE;
+            
+            //newmat = new Matrix[lp.non_zeros];
 
-            newmat = new Matrix[lp.non_zeros];
-            for (i = 0; i < newmat.length; i++)
-                newmat[i] = new Matrix(0, (float)0);
+            // FIXME
+            newmat = lp.alternate_mat;
+            for (i = 0; i < newmat.length; i++) { newmat[i].row_nr = 0; newmat[i].value = 0; }
 
-            inc_mat_space(lp, 0);
+            if (lp.non_zeros > lp.mat_alloc) throw new Error("not enough mat space; this should not happen");
             lp.rows++;
             lp.sum++;
-            inc_row_space(lp);
-
+            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 */
+                    // shift scale
                     for(i=lp.sum; i > lp.rows; i--)
                         lp.scale[i]=lp.scale[i-1];
                     lp.scale[lp.rows]=1;
@@ -404,7 +236,7 @@ public class LinearProgramming {
             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
@@ -432,7 +264,8 @@ public class LinearProgramming {
                     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--)
@@ -551,8 +384,8 @@ public class LinearProgramming {
 
             lp.columns++;
             lp.sum++;
-            inc_col_space(lp);
-            inc_mat_space(lp, lp.rows+1);
+            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)
                 {
@@ -2818,14 +2651,10 @@ public class LinearProgramming {
     } // end of class solve
 
 
-    public static class Matrix
-    {
-        int row_nr;
-        float value;
-        public Matrix(int r, float v) {
-            row_nr = r;
-            value = v;
-        }
+    public static class Matrix {
+        public int row_nr;
+        public float value;
+        public Matrix(int r, float v) { row_nr = r; value = v; }
     }
 
 
@@ -2856,6 +2685,7 @@ public class LinearProgramming {
         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 
@@ -2951,177 +2781,106 @@ public class LinearProgramming {
         float    epsel;              /* ## */
 
 
-        public Problem (int nrows, int ncolumns)
-        {
+        public Problem (int nrows, int ncolumns, int matalloc) {
             int i, nsum;  
-
             nsum=nrows+ncolumns;
-            if(rows < 0 || columns < 0)
-                System.err.print("rows < 0 or columns < 0");
-
-            lp_name = new String("unnamed");
-            active=FALSE;
-            verbose=FALSE;
-            debug=FALSE;
-            trace=FALSE;
-
             rows=nrows;
             columns=ncolumns;
             sum=nsum;
             rows_alloc=rows;
             columns_alloc=columns;
             sum_alloc=sum;
-            names_used=FALSE;
-
-            obj_bound=DEF_INFINITE;
-            infinite=DEF_INFINITE;
-            epsilon=DEF_EPSILON;
-            epsb=DEF_EPSB;
-            epsd=DEF_EPSD;
-            epsel=DEF_EPSEL;
-            non_zeros=0;
-            mat_alloc=1;
-
-            mat = new Matrix[mat_alloc];
-            for (i = 0; i < mat_alloc; i++)
-                mat[i] = new Matrix(0, 0);
-
+            mat_alloc=matalloc;
+            eta_alloc=10000;
+            max_num_inv=DEFNUMINV;
             col_no = new int[mat_alloc];
-            for (i = 0; i < mat_alloc; i++)
-                col_no[i] = 0;
-
             col_end = new int[columns + 1];
-            for (i = 0; i < columns + 1; i++)
-                col_end[i] = 0;
-
             row_end = new int[rows + 1];
-            for (i = 0; i < rows + 1; i++)
-                row_end[i] = 0;
-
-            row_end_valid=FALSE;
-
             orig_rh = new float[rows + 1];
-            for (i = 0; i < rows + 1; i++)
-                orig_rh[i] = 0;
-
             rh = new float[rows + 1];
-            for (i = 0; i < rows + 1; i++)
-                rh[i] = 0;
-
             rhs = new float[rows + 1];
-            for (i = 0; i < rows + 1; i++)
-                rhs[i] = 0;
-
             orig_upbo = new float[sum + 1];
-            for(i = 0; i <= sum; i++)
-                orig_upbo[i]=infinite;
-
             upbo = new float[sum + 1];
-            for (i = 0; i < sum + 1; i++)
-                upbo[i] = 0;
-
             orig_lowbo = new float[sum + 1];
-            for (i = 0; i < sum + 1; i++)
-                orig_lowbo[i] = 0;
-
             lowbo = new float[sum + 1];
-            for (i = 0; i < sum + 1; i++)
-                lowbo[i] = 0;
-
-            basis_valid=TRUE;
-
             bas = new int[rows+1];
-            for (i = 0; i <= rows; i++)
-                bas[i] = 0;
-
             basis = new short[sum + 1];
-            for (i = 0; i <= sum; i++)
-                basis[i] = 0;
-
             lower = new short[sum + 1];
-            for(i = 0; i <= rows; i++)
-                {
-                    bas[i]=i;
-                    basis[i]=TRUE;
-                }
-            for(i = rows + 1; i <= sum; i++)
-                basis[i]=FALSE;
-            for(i = 0 ; i <= sum; i++)
-                lower[i]=TRUE;
-            eta_valid=TRUE;
-            eta_size=0;
-            eta_alloc=10000;
-            max_num_inv=DEFNUMINV;
-
-            nr_lagrange=0;
-
             eta_value = new float[eta_alloc];
-            for (i = 0; i < eta_alloc; i++)
-                eta_value[i] = 0;
-
             eta_row_nr = new int[eta_alloc];
-            for (i = 0; i < eta_alloc; i++)
-                eta_row_nr[i] = 0;
-
             eta_col_end = new int[rows_alloc + max_num_inv];
-            for (i = 0; i < rows_alloc + max_num_inv; i++)
-                eta_col_end[i] = 0;
+            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<mat.length; j++) mat[j] = new Matrix(0, 0);
+            alternate_mat = new Matrix[mat_alloc];
+            for(int j=0; j<alternate_mat.length; j++) alternate_mat[j] = new Matrix(0, 0);
+        }
+
+        public void init(int nrows, int ncolumns) {
+            int i, nsum;  
+            nsum=nrows+ncolumns;
+            lp_name = new String("unnamed");
+            active=FALSE;
+            verbose=FALSE;
+            debug=FALSE;
+            trace=FALSE;
+            rows=nrows;
+            columns=ncolumns;
+            sum=nsum;
+            names_used=FALSE;
+            obj_bound=DEF_INFINITE;
+            infinite=DEF_INFINITE;
+            epsilon=DEF_EPSILON;
+            epsb=DEF_EPSB;
+            epsd=DEF_EPSD;
+            epsel=DEF_EPSEL;
+            non_zeros=0;
 
+            for (i = 0; i < mat_alloc; i++) { mat[i].row_nr = 0; mat[i].value = 0; }
+            for (i = 0; i < mat_alloc; i++)   col_no[i] = 0;
+            for (i = 0; i < columns + 1; i++) col_end[i] = 0;
+            for (i = 0; i < rows + 1; i++)    row_end[i] = 0;
+            for (i = 0; i < rows + 1; i++)   orig_rh[i] = 0;
+            for (i = 0; i < rows + 1; i++)   rh[i] = 0;
+            for (i = 0; i < rows + 1; i++)   rhs[i] = 0;
+            for (i = 0; i <= sum; i++)       orig_upbo[i]=infinite;
+            for (i = 0; i < sum + 1; i++)    upbo[i] = 0;
+            for (i = 0; i < sum + 1; i++)    orig_lowbo[i] = 0;
+            for (i = 0; i < sum + 1; i++)    lowbo[i] = 0;
+            for (i = 0; i <= rows; i++)      bas[i] = 0;
+            for (i = 0; i <= sum; i++)       basis[i] = 0;
+            for (i = 0; i <= rows; i++)     { bas[i]=i; basis[i]=TRUE; }
+            for (i = rows + 1; i <= sum; i++) basis[i]=FALSE;
+            for (i = 0 ; i <= sum; i++)       lower[i]=TRUE;
+            for (i = 0; i < eta_alloc; i++) eta_value[i] = 0;
+            for (i = 0; i < eta_alloc; i++) eta_row_nr[i] = 0;
+            for (i = 0; i < rows_alloc + max_num_inv; i++) eta_col_end[i] = 0;
+            for (i = 0; i <= sum; i++) solution[i] = 0;
+            for (i = 0; i <= sum; i++) best_solution[i] = 0;
+            for (i = 0; i <= rows; i++) duals[i] = 0;
+            for (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;
-
-            solution = new float[sum + 1];
-            for (i = 0; i <= sum; i++)
-                solution[i] = 0;
-
-            best_solution = new float[sum + 1];
-            for (i = 0; i <= sum; i++)
-                best_solution[i] = 0;
-
-            duals = new float[rows + 1];
-            for (i = 0; i <= rows; i++)
-                duals[i] = 0;
-
+            basis_valid=TRUE; 
+            eta_valid=TRUE;
+            eta_size=0;
+            nr_lagrange=0;
             maximise = FALSE;
             floor_first = TRUE;
-
             scaling_used = FALSE;
             columns_scaled = FALSE;
-
-            ch_sign = new short[rows + 1];
-            for(i = 0; i <= rows; i++)
-                ch_sign[i] = FALSE;
-
             valid = FALSE; 
-        } // end of constructor from row and column
-
-        //***************************************
-        // return the ith member of the best_solution[]
-        //
-        public float getBestSolution(int i) {
-            return best_solution[i];
-        }
-
-        //***************************************
-        // get the number of rows
-        //
-        public int getRows() {
-            return rows;
-        }
-
-        //***************************************
-        // get the number of columns
-        //
-        public int getcolumns() {
-            return columns;
         }
 
-    } // end of class Problem
+    }
 
     private final static float round( float val, float eps) { return (Math.abs(val) < eps) ? 0 : val; }