// Copyright (c) 1997-2007 ETH Zurich (Switzerland). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // You can redistribute it and/or modify it under the terms of the GNU // General Public License as published by the Free Software Foundation, // either version 3 of the License, or (at your option) any later version. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // // Author(s) : Bernd Gaertner #ifndef CGAL_QP_FUNCTIONS_IMPL_H #define CGAL_QP_FUNCTIONS_IMPL_H #include #include #include #include #include #include namespace CGAL { namespace QP_functions_detail { // test whether the system is of the form A x == b (equations only) template bool is_in_equational_form (const R& r) { typename R::R_iterator it = r.get_r(); typename R::R_iterator end = it + r.get_m(); for (; it < end; ++it) if (*it != CGAL::EQUAL) return false; return true; } // test whether the row vectors of A that correpsond to equations // are linearly independent; this is done using type ET. The value // type of LinearInequalitySystem must be convertible to ET template bool has_linearly_independent_equations (const Ar& ar, const ET& /*dummy*/) { // we solve the following auxiliary LP, using exact type ET: // -------- // min 0 // A x r 0 // x >= 0 // -------- // Then A has linearly independent equations if and only if all // artificials have left the basis after phase I; the QP_solver // diagnostics tells us this // // auxiliary LP type typedef typename std::iterator_traits::value_type C_value; typedef typename std::iterator_traits::value_type B_value; typedef Const_oneset_iterator C_iterator; typedef Const_oneset_iterator B_iterator; typedef Nonnegative_linear_program_from_iterators LP; // auxiliary LP LP lp (ar.get_n(), ar.get_m(), ar.get_a(), B_iterator(0), ar.get_r(), C_iterator(0)); // solver Tags typedef QP_solver_impl::QP_tags< Tag_true, // Is_linear Tag_true> // Is_nonnegative Tags; // solver type typedef QP_solver Solver; // now solve auxiliary LP and compute predicate value Solver solver (lp); return !solver.diagnostics.redundant_equations; } // helper for MPS output: BOUNDS template void print_bounds (std::ostream& , const P& , CGAL::Tag_true /*is_nonnegative*/) { // nop (default bounds are nonnegative) } // helper for MPS output: BOUNDS template void print_bounds (std::ostream& out, const P& p, CGAL::Tag_false /*is_nonnegative*/) { typename P::FL_iterator fl = p.get_fl(); typename P::FU_iterator fu = p.get_fu(); typename P::L_iterator l = p.get_l(); typename P::U_iterator u = p.get_u(); int n = p.get_n(); out << "BOUNDS\n"; for (int j=0; j void print_qmatrix (std::ostream& , const P& , CGAL::Tag_true /*is_linear*/) { // nop } // helper for MPS output: DMATRIX/QMATRIX template void print_qmatrix (std::ostream& out, const P& p, CGAL::Tag_false /*is_linear*/) { typename P::D_iterator it = p.get_d(); int n = p.get_n(); bool empty_D = true; for (int i=0; i bool are_equal_qp (const Quadratic_program1 &qp1, const Quadratic_program2 &qp2) { bool return_val = true; // check n if (qp1.get_n() != qp2.get_n()) { std::cerr << "Equality test fails with n: " << qp1.get_n() << " vs. " << qp2.get_n() << std::endl; return false; // wildly wrong, abort now } // check m if (qp1.get_m() != qp2.get_m()) { std::cerr << "Equality test fails with m: " << qp1.get_m() << " vs. " << qp2.get_m() << std::endl; return false; // wildly wrong, abort now } int n = qp1.get_n(); int m = qp1.get_m(); // check A typename Quadratic_program1::A_iterator a1 = qp1.get_a(); typename Quadratic_program2::A_iterator a2 = qp2.get_a(); for (int j=0; j void print_program (std::ostream& out, const P& p, const std::string& problem_name, Is_linear is_linear, Is_nonnegative is_nonnegative) { // NAME: out << "NAME " << problem_name << "\n"; int n = p.get_n(); int m = p.get_m(); // ROWS section: typename P::R_iterator r = p.get_r(); out << "ROWS\n" << " N obj\n"; // for the objective function for (int i=0; i::value_type IT; out << "COLUMNS\n"; for (int j=0; j Quadratic_program_solution solve_program (const Program &p, const ET&, Is_linear, Is_nonnegative, const Quadratic_program_options& options) { typedef QP_solver< Program, ET, QP_solver_impl::QP_tags > Solver; const Solver* s = new Solver(p, options); Quadratic_program_solution solution(s); if (options.get_auto_validation()) { // validate solution if (options.get_verbosity() > 0) std::cout << "Validating solution...\n"; if (!solution.solves_program(p, Is_linear(), Is_nonnegative())) { // write log file std::ofstream out("QP_solver.log"); out << "Error: Program solution is invalid\n" << " (error message: " << solution.get_error() << ")\n" << "------------------\n" << "Solution function:\n" << "------------------\n"; print_solution_function (out, Is_linear(), Is_nonnegative()); out << "\n" << "--------\n" << "Program:\n" << "--------\n"; print_program (out, p, "unsolved", Is_linear(), Is_nonnegative()); out << "--------\n" << "Options:\n" << "--------\n" << options << std::endl; // print warning std::cerr << "Error: Program solution is invalid " << "(see QP_solver.log for details)" << std::endl; } } return solution; } } template Quadratic_program_solution solve_quadratic_program (const QuadraticProgram &qp, const ET& dummy, const Quadratic_program_options& options) { return QP_functions_detail:: solve_program(qp, dummy, Tag_false(), Tag_false(), options); } template Quadratic_program_solution solve_nonnegative_quadratic_program (const NonnegativeQuadraticProgram &qp, const ET& dummy, const Quadratic_program_options& options) { return QP_functions_detail:: solve_program(qp, dummy, Tag_false(), Tag_true(), options); } template Quadratic_program_solution solve_linear_program (const LinearProgram &lp, const ET& dummy, const Quadratic_program_options& options) { return QP_functions_detail:: solve_program(lp, dummy, Tag_true(), Tag_false(), options); } template Quadratic_program_solution solve_nonnegative_linear_program (const NonnegativeLinearProgram &lp, const ET& dummy, const Quadratic_program_options& options) { return QP_functions_detail:: solve_program(lp, dummy, Tag_true(), Tag_true(), options); } } //namespace CGAL #endif // CGAL_QP_FUNCTIONS_IMPL_H