Scippy

SoPlex

Sequential object-oriented simPlex

spxproof.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2015 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file spxproof.cpp
17  * @brief provable bounds
18  */
19 
20 #ifdef WITH_EASYVAL
21 
22 #include <assert.h>
23 #include <iostream>
24 
25 #include "Easyval.hh"
26 
27 #include "spxdefines.h"
28 #include "spxsolver.h"
29 
30 //#undef DEBUG
31 //#define DEBUG(x) {x}
32 
33 namespace soplex
34 {
35 
36 Real SPxSolver::provedBound(Vector& dualsol, const Vector& objvec) const
37 {
38  DVector lhsvec(dim());
39  DVector rhsvec(dim());
40  Real eps = opttol();
41  Easyval ytb;
42  Easyval scalprod;
43  int sen = (spxSense() == MINIMIZE ? +1 : -1);
44 
45  assert(rep() == COLUMN);
46 
47  /* for minimization problems: calculate y^Tb + min{(c^T - y^TA)x | l <= x <= u}
48  * for maximization problems: calculate y^Tb + max{(c^T - y^TA)x | l <= x <= u}
49  */
50 
51  getLhs(lhsvec);
52  getRhs(rhsvec);
53 
54  /* 1. calculate y^Tb */
55  ytb = 0;
56  for( int i = 0; i < dim(); ++i )
57  {
58  Easyval b;
59 
60  if( sen * dualsol[i] > eps )
61  b = lhsvec[i];
62  else if( sen * dualsol[i] < -eps )
63  b = rhsvec[i];
64  else
65  dualsol[i] = 0.0;
66 
67  Easyval y(dualsol[i]);
68 
69  ytb += y*b;
70  MSG_DEBUG( std::cout << "DPROOF02 y[" << i << "] = " << dualsol[i]
71  << ", lhs[" << i << "] = " << lhsvec[i]
72  << ", rhs[" << i << "] = " << rhsvec[i]
73  << " -> ytb = " << ytb << std::endl; )
74  }
75 
76  /* 2. calculate min/max{(c^T - y^TA)x} */
77  scalprod = 0;
78  for( int j = 0; j < coDim(); ++j )
79  {
80  LPCol col;
81  Easyval x(lower(j), upper(j)); /* add/sub machine epsilon? */
82 
83  getCol(j, col);
84  SVector colvec = col.colVector(); /* make this faster !!! */
85  Easyval diff = objvec[j];
86  MSG_DEBUG( std::cout << "DPROOF03 obj[" << j << "] = " << objvec[j]
87  << " -> diff = " << diff << std::endl; )
88 
89  for( int i = 0; i < colvec.size(); ++i )
90  {
91  Easyval y(dualsol[colvec.index(i)]);
92  Easyval a(colvec.value(i));
93  diff -= y*a;
94  MSG_DEBUG( std::cout << "DPROOF04 y[" << colvec.index(i) << "] = "
95  << dualsol[colvec.index(i)]
96  << ", a[" << colvec.index(i) << "] = "
97  << colvec.value(i)
98  << " -> diff = " << diff << std::endl; )
99  }
100 
101  diff *= x;
102  scalprod += diff;
103  MSG_DEBUG( std::cout << "DPROOF05 x[" << j << "] = " << x
104  << " -> diff = " << diff << std::endl
105  << " ------------> scalprod = "
106  << scalprod << std::endl; )
107  }
108 
109  /* add y^Tb */
110  scalprod += ytb;
111  MSG_INFO1( spxout, spxout << "IPROOF01 proved bound = " << scalprod << std::endl; );
112 
113  /* depending on the objective sense, choose min or max */
114  if( spxSense() == MINIMIZE )
115  return scalprod.getInf();
116  else
117  return scalprod.getSup();
118 }
119 
120 Real SPxSolver::provedDualbound() const
121 {
122  DVector dualsol(dim());
123  DVector objvec(coDim());
124 
125  getDual(dualsol);
126  getObj(objvec);
127 
128  return provedBound(dualsol, objvec);
129 }
130 
131 bool SPxSolver::isProvenInfeasible() const
132 {
133  DVector dualfarkas(dim());
134  DVector zerovec(coDim());
135 
136  getDualfarkas(dualfarkas);
137  zerovec.clear();
138 
139  return (provedBound(dualfarkas, zerovec) > 0.0);
140 }
141 
142 }
143 
144 #endif // WITH_EASYVAL