Scippy

SoPlex

Sequential object-oriented simPlex

spxgeometsc.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-2017 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 spxgeometsc.cpp
17  * @brief Geometric mean row/column scaling.
18  */
19 #include <assert.h>
20 
21 #include "spxgeometsc.h"
22 #include "spxout.h"
23 #include "spxlpbase.h"
24 #include "spxequilisc.h"
25 
26 namespace soplex
27 {
28 
30  const SVSet* vecset,
31  const std::vector<Real>& coScaleval,
32  std::vector<Real>& scaleval)
33 {
34  Real pmax = 0.0;
35 
36  assert(scaleval.size() == unsigned(vecset->num()));
37 
38  for( int i = 0; i < vecset->num(); ++i )
39  {
40  const SVector& vec = (*vecset)[i];
41 
42  Real maxi = 0.0;
43  Real mini = infinity;
44 
45  for( int j = 0; j < vec.size(); ++j )
46  {
47  const Real x = spxAbs(vec.value(j) * coScaleval[unsigned(vec.index(j))]);
48 
49  if (!isZero(x))
50  {
51  if (x > maxi)
52  maxi = x;
53  if (x < mini)
54  mini = x;
55  }
56  }
57  // empty rows/cols are possible
58  if (mini == infinity || maxi == 0.0)
59  {
60  mini = 1.0;
61  maxi = 1.0;
62  }
63  assert(mini < infinity);
64  assert(maxi > 0.0);
65 
66  scaleval[unsigned(i)] = 1.0 / spxSqrt(mini * maxi);
67 
68  const Real p = maxi / mini;
69 
70  if (p > pmax)
71  pmax = p;
72  }
73  return pmax;
74 }
75 
76 
77 SPxGeometSC::SPxGeometSC(bool equilibrate, int maxIters, Real minImpr, Real goodEnough)
78  : SPxScaler("Geometric")
79  , postequilibration(equilibrate)
80  , m_maxIterations(maxIters)
81  , m_minImprovement(minImpr)
82  , m_goodEnoughRatio(goodEnough)
83 {
84  assert(maxIters > 0);
85  assert(minImpr > 0.0 && minImpr <= 1.0);
86  assert(goodEnough >= 0.0);
87 }
88 
90  : SPxScaler(old)
95 {
96  assert(m_maxIterations > 0);
97  assert(m_minImprovement > 0.0 && m_minImprovement <= 1.0);
98  assert(m_goodEnoughRatio >= 0.0);
99 }
100 
102 {
103  if (this != &rhs)
104  {
106  }
107 
108  return *this;
109 }
110 
111 void SPxGeometSC::scale(SPxLPBase<Real>& lp, bool persistent)
112 {
113 
114  MSG_INFO1( (*spxout), (*spxout) << "Geometric scaling LP" << (persistent ? " (persistent)" : "") << (postequilibration ? " with post-equilibration" : "") << std::endl; )
115 
116  setup(lp);
117 
118  /* We want to do that direction first, with the lower ratio.
119  * See SPxEquiliSC::scale() for a reasoning.
120  */
121  const Real colratio = maxColRatio(lp);
122  const Real rowratio = maxRowRatio(lp);
123 
124  const bool colFirst = colratio < rowratio;
125 
126  Real p0start;
127  Real p1start;
128 
129  if( colFirst )
130  {
131  p0start = colratio;
132  p1start = rowratio;
133  }
134  else
135  {
136  p0start = rowratio;
137  p1start = colratio;
138  }
139 
140  MSG_INFO2( (*spxout), (*spxout) << "before scaling:"
141  << " min= " << lp.minAbsNzo()
142  << " max= " << lp.maxAbsNzo()
143  << " col-ratio= " << colratio
144  << " row-ratio= " << rowratio
145  << std::endl; )
146 
147  // perform geometric scaling only if maximum ratio is above threshold
148  bool geoscale = p1start > m_goodEnoughRatio;
149 
150  if( !geoscale )
151  {
152  MSG_INFO2( (*spxout), (*spxout) << "No geometric scaling done, ratio good enough" << std::endl; )
153 
154  if( !postequilibration )
155  {
156  lp.setScalingInfo(true);
157  return;
158  }
159 
160  MSG_INFO2( (*spxout), (*spxout) << " ... but will still perform equilibrium scaling" << std::endl; )
161  }
162 
163  std::vector<Real> rowscale(unsigned(lp.nRows()), 1.0);
164  std::vector<Real> colscale(unsigned(lp.nCols()), 1.0);
165 
166  Real p0 = 0.0;
167  Real p1 = 0.0;
168 
169  if( geoscale )
170  {
171  Real p0prev = p0start;
172  Real p1prev = p1start;
173 
174  // we make at most maxIterations.
175  for( int count = 0; count < m_maxIterations; count++ )
176  {
177  if( colFirst )
178  {
179  p0 = computeScalingVec(lp.colSet(), rowscale, colscale);
180  p1 = computeScalingVec(lp.rowSet(), colscale, rowscale);
181  }
182  else
183  {
184  p0 = computeScalingVec(lp.rowSet(), colscale, rowscale);
185  p1 = computeScalingVec(lp.colSet(), rowscale, colscale);
186  }
187 
188  MSG_INFO3( (*spxout), (*spxout) << "Geometric scaling round " << count
189  << " col-ratio= " << (colFirst ? p0 : p1)
190  << " row-ratio= " << (colFirst ? p1 : p0)
191  << std::endl; )
192 
193  if( p0 > m_minImprovement * p0prev && p1 > m_minImprovement * p1prev )
194  break;
195 
196  p0prev = p0;
197  p1prev = p1;
198  }
199 
200  // perform geometric scaling only if there is enough (default 15%) improvement.
201  geoscale = (p0 <= m_minImprovement * p0start || p1 <= m_minImprovement * p1start);
202  }
203 
204  if( !geoscale && !postequilibration )
205  {
206  MSG_INFO2( (*spxout), (*spxout) << "No geometric scaling done." << std::endl; )
207  lp.setScalingInfo(true);
208  }
209  else
210  {
211  DataArray<int>& colscaleExp = *m_activeColscaleExp;
212  DataArray<int>& rowscaleExp = *m_activeRowscaleExp;
213 
214  if( postequilibration )
215  {
216  if( !geoscale )
217  {
218  std::fill(rowscale.begin(), rowscale.end(), 1.0);
219  std::fill(colscale.begin(), colscale.end(), 1.0);
220  }
221  SPxEquiliSC::computePostequiExpVecs(lp, rowscale, colscale, rowscaleExp, colscaleExp);
222  }
223  else
224  {
225  computeExpVec(colscale, colscaleExp);
226  computeExpVec(rowscale, rowscaleExp);
227  }
228 
229  applyScaling(lp);
230 
231  MSG_INFO3( (*spxout), (*spxout) << "Row scaling min= " << minAbsRowscale()
232  << " max= " << maxAbsRowscale()
233  << std::endl
234  << "IGEOSC06 Col scaling min= " << minAbsColscale()
235  << " max= " << maxAbsColscale()
236  << std::endl; )
237 
238  MSG_INFO2( (*spxout), (*spxout) << "after scaling: "
239  << " min= " << lp.minAbsNzo(false)
240  << " max= " << lp.maxAbsNzo(false)
241  << " col-ratio= " << maxColRatio(lp)
242  << " row-ratio= " << maxRowRatio(lp)
243  << std::endl; )
244  }
245 }
246 
247 
248 } // namespace soplex
Rational spxAbs(const Rational &r)
Absolute.
Definition: rational.cpp:3909
virtual void setup(SPxLPBase< Real > &lp)
clear and setup scaling arrays in the LP
Definition: spxscaler.cpp:124
const Real m_goodEnoughRatio
no scaling needed if ratio is less than this.
Definition: spxgeometsc.h:45
const int m_maxIterations
maximum number of scaling iterations.
Definition: spxgeometsc.h:43
virtual R minAbsNzo(bool unscaled=true) const
Absolute smallest non-zero element in (possibly scaled) LP.
THREADLOCAL const Real infinity
Definition: spxdefines.cpp:26
virtual void scale(SPxLPBase< Real > &lp, bool persistent=true) override
Scale the loaded SPxLP.
Geometric mean row/column scaling.This SPxScaler implementation performs geometric mean scaling of th...
Definition: spxgeometsc.h:35
int size() const
Number of used indices.
Definition: svectorbase.h:152
LP geometric mean scaling.
virtual Real maxRowRatio(const SPxLPBase< Real > &lp) const
maximum ratio between absolute biggest and smallest element in any row.
Definition: spxscaler.cpp:855
R & value(int n)
Reference to value of n &#39;th nonzero.
Definition: svectorbase.h:252
Saving LPs in a form suitable for SoPlex.
Wrapper for different output streams and verbosity levels.
int nRows() const
Returns number of rows in LP.
Definition: spxlpbase.h:151
DataArray< int > * m_activeColscaleExp
pointer to currently active column scaling factors
Definition: spxscaler.h:84
virtual void applyScaling(SPxLPBase< Real > &lp)
applies m_colscale and m_rowscale to the lp.
Definition: spxscaler.cpp:171
double Real
Definition: spxdefines.h:218
SPxOut * spxout
message handler
Definition: spxscaler.h:88
int & index(int n)
Reference to index of n &#39;th nonzero.
Definition: svectorbase.h:234
#define MSG_INFO2(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO2.
Definition: spxdefines.h:120
SPxGeometSC & operator=(const SPxGeometSC &)
assignment operator
Real spxSqrt(Real a)
returns square root
Definition: spxdefines.h:334
virtual Real maxAbsColscale() const
absolute biggest column scaling factor
Definition: spxscaler.cpp:773
virtual Real minAbsColscale() const
absolute smallest column scaling factor
Definition: spxscaler.cpp:760
#define MSG_INFO3(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO3.
Definition: spxdefines.h:122
void computeExpVec(const std::vector< Real > &vec, DataArray< int > &vecExp)
round vector entries to power of 2
Definition: spxscaler.cpp:889
DataArray< int > * m_activeRowscaleExp
pointer to currently active row scaling factors
Definition: spxscaler.h:85
SPxGeometSC(bool equilibrate=false, int maxIters=8, Real minImpr=0.85, Real goodEnough=1e3)
default constructor (this scaler makes no use of inherited members m_colFirst and m_doBoth) ...
Definition: spxgeometsc.cpp:77
void setScalingInfo(bool scaled)
set whether the LP is scaled or not
Definition: spxlpbase.h:145
virtual Real minAbsRowscale() const
absolute smallest row scaling factor
Definition: spxscaler.cpp:787
static void computePostequiExpVecs(const SPxLPBase< Real > &lp, const std::vector< Real > &preRowscale, const std::vector< Real > &preColscale, DataArray< int > &rowscaleExp, DataArray< int > &colscaleExp)
compute equilibrium scaling rounded to power of 2 for existing Real scaling factors (preRowscale...
Everything should be within this namespace.
LP scaler abstract base class.Instances of classes derived from SPxScaler may be loaded to SoPlex in ...
Definition: spxscaler.h:76
virtual Real maxAbsRowscale() const
absolute biggest row scaling factor
Definition: spxscaler.cpp:800
#define MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
Definition: spxdefines.h:118
static Real computeScalingVec(const SVSet *vecset, const std::vector< Real > &coScaleval, std::vector< Real > &scaleval)
Definition: spxgeometsc.cpp:29
int nCols() const
Returns number of columns in LP.
Definition: spxlpbase.h:157
const SVSetBase< R > * colSet() const
Returns the complete SVSetBase.
Definition: lpcolsetbase.h:68
bool isZero(Real a, Real eps=Param::epsilon())
returns true iff |a| <= eps
Definition: spxdefines.h:416
virtual Real maxColRatio(const SPxLPBase< Real > &lp) const
maximum ratio between absolute biggest and smallest element in any column.
Definition: spxscaler.cpp:817
LP equilibrium scaling.
const SVSetBase< R > * rowSet() const
Returns the complete SVSet.
Definition: lprowsetbase.h:69
int num() const
Current number of SVectorBases.
Definition: svsetbase.h:746
virtual R maxAbsNzo(bool unscaled=true) const
Absolute biggest non-zero element in (in rational case possibly scaled) LP.
const bool postequilibration
equilibrate after geometric scaling?
Definition: spxgeometsc.h:42
const Real m_minImprovement
improvement necessary to carry on. (Bixby said Fourer said in MP 23, 274 ff. that 0...
Definition: spxgeometsc.h:44
SPxScaler & operator=(const SPxScaler &)
assignment operator
Definition: spxscaler.cpp:84