restore the previous state

Signed-off-by: Lev Nachmanson <[email protected]>

src/math/lp/CMakeLists.txt

@@ -20,8 +20,11 @@ z3_add_component(lp
     lar_core_solver.cpp
     lp_core_solver_base.cpp
     lp_dual_core_solver.cpp
+    lp_dual_simplex.cpp
     lp_primal_core_solver.cpp
+    lp_primal_simplex.cpp
     lp_settings.cpp
+    lp_solver.cpp
     lu.cpp
     lp_utils.cpp
     matrix.cpp

src/math/lp/indexed_value.h

@@ -43,4 +43,15 @@ class indexed_value {
         m_value = val;
     }
 };
+#ifdef Z3DEBUG
+template <typename X>
+bool check_vector_for_small_values(indexed_vector<X> & w, lp_settings & settings) {
+    for (unsigned i : w.m_index) {
+        const X & v = w[i];
+        if ((!is_zero(v)) && settings.abs_val_is_smaller_than_drop_tolerance(v))
+            return false;
+    }
+    return true;
+}
+#endif
 }

src/math/lp/lar_core_solver.h

@@ -651,7 +651,35 @@ class lar_core_solver  {
         }
     }
 
-
+    void scale_problem_for_doubles(
+                        static_matrix<double, double>& A,        
+                        vector<double> & lower_bounds,
+                        vector<double> & upper_bounds) {
+        vector<double> column_scale_vector;
+        vector<double> right_side_vector(A.column_count());
+        settings().reps_in_scaler = 5;
+        scaler<double, double > scaler(right_side_vector,
+                                       A,
+                                       settings().scaling_minimum,
+                                       settings().scaling_maximum,
+                                       column_scale_vector,
+                                       settings());
+        if (! scaler.scale()) {
+            // the scale did not succeed, unscaling
+            A.clear();
+            create_double_matrix(A);
+        } else {
+            for (unsigned j = 0; j < A.column_count(); j++) {
+                if (m_r_solver.column_has_upper_bound(j)) {
+                    upper_bounds[j] /= column_scale_vector[j];
+                }
+                if (m_r_solver.column_has_lower_bound(j)) {
+                    lower_bounds[j] /= column_scale_vector[j];
+                }
+            }
+        }
+
+    }
     // returns the trace of basis changes
     vector<unsigned> find_solution_signature_with_doubles(lar_solution_signature & signature) {
         if (m_d_solver.m_factorization == nullptr || m_d_solver.m_factorization->get_status() != LU_status::OK) {

src/math/lp/lp_dual_simplex.cpp

@@ -0,0 +1,24 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#include "math/lp/lp_dual_simplex_def.h"
+template lp::mpq lp::lp_dual_simplex<lp::mpq, lp::mpq>::get_current_cost() const;
+template void lp::lp_dual_simplex<lp::mpq, lp::mpq>::find_maximal_solution();
+template double lp::lp_dual_simplex<double, double>::get_current_cost() const;
+template void lp::lp_dual_simplex<double, double>::find_maximal_solution();

src/math/lp/lp_dual_simplex.h

@@ -0,0 +1,93 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+#include "util/vector.h"
+#include "math/lp/lp_utils.h"
+#include "math/lp/lp_solver.h"
+#include "math/lp/lp_dual_core_solver.h"
+namespace lp {
+
+template <typename T, typename X>
+class lp_dual_simplex: public lp_solver<T, X> {
+    lp_dual_core_solver<T, X> * m_core_solver;
+    vector<T> m_b_copy;
+    vector<T> m_lower_bounds; // We don't have a convention here that all low bounds are zeros. At least it does not hold for the first stage solver
+    vector<column_type> m_column_types_of_core_solver;
+    vector<column_type> m_column_types_of_logicals;
+    vector<bool>  m_can_enter_basis;
+public:
+    ~lp_dual_simplex() override {
+        delete m_core_solver;
+    }
+
+    lp_dual_simplex() : m_core_solver(nullptr) {}
+
+
+    void decide_on_status_after_stage1();
+
+    void fix_logical_for_stage2(unsigned j);
+
+    void fix_structural_for_stage2(unsigned j);
+
+    void unmark_boxed_and_fixed_columns_and_fix_structural_costs();
+
+    void restore_right_sides();
+
+    void solve_for_stage2();
+
+    void fill_x_with_zeros();
+
+    void stage1();
+
+    void stage2();
+
+    void fill_first_stage_solver_fields();
+
+    column_type get_column_type(unsigned j);
+
+    void fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j);
+
+    void fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j);
+
+    void fill_costs_and_bounds_and_column_types_for_the_first_stage_solver();
+
+    void set_type_for_logical(unsigned j, column_type col_type) {
+        this->m_column_types_of_logicals[j - this->number_of_core_structurals()] = col_type;
+    }
+
+    void fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,
+                                                              unsigned & slack_var,
+                                                                     unsigned & artificial);
+
+    void augment_matrix_A_and_fill_x_and_allocate_some_fields();
+
+
+
+    void copy_m_b_aside_and_set_it_to_zeros();
+
+    void find_maximal_solution() override;
+
+    T get_column_value(unsigned column) const override {
+        return this->get_column_value_with_core_solver(column, m_core_solver);
+    }
+
+    T get_current_cost() const override;
+};
+}