1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2022 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.hpp
17 * @brief Geometric mean row/column scaling.
18 */
19 #include <assert.h>
20
|
(1) Event include_recursion: |
#include file "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scipoptsuite-8.0.1/soplex/src/soplex/spxgeometsc.h" includes itself: spxgeometsc.h -> spxgeometsc.hpp -> spxgeometsc.h |
|
(2) Event caretline: |
^ |
21 #include "soplex/spxgeometsc.h"
22 #include "soplex/spxout.h"
23 #include "soplex/spxlpbase.h"
24 #include "soplex/spxequilisc.h"
25
26 namespace soplex
27 {
28
29 template <class R>
30 static R computeScalingVec(
31 const SVSetBase<R>* vecset,
32 const std::vector<R>& coScaleval,
33 std::vector<R>& scaleval)
34 {
35 R pmax = 0.0;
36
37 assert(scaleval.size() == unsigned(vecset->num()));
38
39 for(int i = 0; i < vecset->num(); ++i)
40 {
41 const SVectorBase<R>& vec = (*vecset)[i];
42
43 R maxi = 0.0;
44 R mini = R(infinity);
45
46 for(int j = 0; j < vec.size(); ++j)
47 {
48 const R x = spxAbs(vec.value(j) * coScaleval[unsigned(vec.index(j))]);
49
50 if(!isZero(x))
51 {
52 if(x > maxi)
53 maxi = x;
54
55 if(x < mini)
56 mini = x;
57 }
58 }
59
60 // empty rows/cols are possible
61 if(mini == R(infinity) || maxi == 0.0)
62 {
63 mini = 1.0;
64 maxi = 1.0;
65 }
66
67 assert(mini < R(infinity));
68 assert(maxi > 0.0);
69
70 scaleval[unsigned(i)] = 1.0 / spxSqrt(mini * maxi);
71
72 const R p = maxi / mini;
73
74 if(p > pmax)
75 pmax = p;
76 }
77
78 return pmax;
79 }
80
81 template <class R>
82 SPxGeometSC<R>::SPxGeometSC(bool equilibrate, int maxIters, R minImpr, R goodEnough)
83 : SPxScaler<R>("Geometric")
84 , postequilibration(equilibrate)
85 , m_maxIterations(maxIters)
86 , m_minImprovement(minImpr)
87 , m_goodEnoughRatio(goodEnough)
88 {
89 assert(maxIters > 0);
90 assert(minImpr > 0.0 && minImpr <= 1.0);
91 assert(goodEnough >= 0.0);
92 }
93
94 template <class R>
95 SPxGeometSC<R>::SPxGeometSC(const SPxGeometSC<R>& old)
96 : SPxScaler<R>(old)
97 , postequilibration(old.postequilibration)
98 , m_maxIterations(old.m_maxIterations)
99 , m_minImprovement(old.m_minImprovement)
100 , m_goodEnoughRatio(old.m_goodEnoughRatio)
101 {
102 assert(m_maxIterations > 0);
103 assert(m_minImprovement > 0.0 && m_minImprovement <= 1.0);
104 assert(m_goodEnoughRatio >= 0.0);
105 }
106
107 template <class R>
108 SPxGeometSC<R>& SPxGeometSC<R>::operator=(const SPxGeometSC<R>& rhs)
109 {
110 if(this != &rhs)
111 {
112 SPxScaler<R>::operator=(rhs);
113 }
114
115 return *this;
116 }
117
118 template <class R>
119 void SPxGeometSC<R>::scale(SPxLPBase<R>& lp, bool persistent)
120 {
121
122 MSG_INFO1((*this->spxout), (*this->spxout) << "Geometric scaling LP" <<
123 (persistent ? " (persistent)" : "") << (postequilibration ? " with post-equilibration" : "") <<
124 std::endl;)
125
126 this->setup(lp);
127
128 /* We want to do that direction first, with the lower ratio.
129 * See SPxEquiliSC<R>::scale() for a reasoning.
130 */
131 const R colratio = this->maxColRatio(lp);
132 const R rowratio = this->maxRowRatio(lp);
133
134 const bool colFirst = colratio < rowratio;
135
136 R p0start;
137 R p1start;
138
139 if(colFirst)
140 {
141 p0start = colratio;
142 p1start = rowratio;
143 }
144 else
145 {
146 p0start = rowratio;
147 p1start = colratio;
148 }
149
150 MSG_INFO2((*this->spxout), (*this->spxout) << "before scaling:"
151 << " min= " << lp.minAbsNzo()
152 << " max= " << lp.maxAbsNzo()
153 << " col-ratio= " << colratio
154 << " row-ratio= " << rowratio
155 << std::endl;)
156
157 // perform geometric scaling only if maximum ratio is above threshold
158 bool geoscale = p1start > m_goodEnoughRatio;
159
160 if(!geoscale)
161 {
162 MSG_INFO2((*this->spxout), (*this->spxout) << "No geometric scaling done, ratio good enough" <<
163 std::endl;)
164
165 if(!postequilibration)
166 {
167 lp.setScalingInfo(true);
168 return;
169 }
170
171 MSG_INFO2((*this->spxout), (*this->spxout) << " ... but will still perform equilibrium scaling" <<
172 std::endl;)
173 }
174
175 std::vector<R> rowscale(unsigned(lp.nRows()), 1.0);
176 std::vector<R> colscale(unsigned(lp.nCols()), 1.0);
177
178 R p0 = 0.0;
179 R p1 = 0.0;
180
181 if(geoscale)
182 {
183 R p0prev = p0start;
184 R p1prev = p1start;
185
186 // we make at most maxIterations.
187 for(int count = 0; count < m_maxIterations; count++)
188 {
189 if(colFirst)
190 {
191 p0 = computeScalingVec(lp.colSet(), rowscale, colscale);
192 p1 = computeScalingVec(lp.rowSet(), colscale, rowscale);
193 }
194 else
195 {
196 p0 = computeScalingVec(lp.rowSet(), colscale, rowscale);
197 p1 = computeScalingVec(lp.colSet(), rowscale, colscale);
198 }
199
200 MSG_INFO3((*this->spxout), (*this->spxout) << "Geometric scaling round " << count
201 << " col-ratio= " << (colFirst ? p0 : p1)
202 << " row-ratio= " << (colFirst ? p1 : p0)
203 << std::endl;)
204
205 if(p0 > m_minImprovement * p0prev && p1 > m_minImprovement * p1prev)
206 break;
207
208 p0prev = p0;
209 p1prev = p1;
210 }
211
212 // perform geometric scaling only if there is enough (default 15%) improvement.
213 geoscale = (p0 <= m_minImprovement * p0start || p1 <= m_minImprovement * p1start);
214 }
215
216 if(!geoscale && !postequilibration)
217 {
218 MSG_INFO2((*this->spxout), (*this->spxout) << "No geometric scaling done." << std::endl;)
219 lp.setScalingInfo(true);
220 }
221 else
222 {
223 DataArray<int>& colscaleExp = *this->m_activeColscaleExp;
224 DataArray<int>& rowscaleExp = *this->m_activeRowscaleExp;
225
226 if(postequilibration)
227 {
228 if(!geoscale)
229 {
230 std::fill(rowscale.begin(), rowscale.end(), 1.0);
231 std::fill(colscale.begin(), colscale.end(), 1.0);
232 }
233
234 SPxEquiliSC<R>::computePostequiExpVecs(lp, rowscale, colscale, rowscaleExp, colscaleExp);
235 }
236 else
237 {
238 this->computeExpVec(colscale, colscaleExp);
239 this->computeExpVec(rowscale, rowscaleExp);
240 }
241
242 this->applyScaling(lp);
243
244 MSG_INFO3((*this->spxout), (*this->spxout) << "Row scaling min= " << this->minAbsRowscale()
245 << " max= " << this->maxAbsRowscale()
246 << std::endl
247 << "IGEOSC06 Col scaling min= " << this->minAbsColscale()
248 << " max= " << this->maxAbsColscale()
249 << std::endl;)
250
251 MSG_INFO2((*this->spxout), (*this->spxout) << "after scaling: "
252 << " min= " << lp.minAbsNzo(false)
253 << " max= " << lp.maxAbsNzo(false)
254 << " col-ratio= " << this->maxColRatio(lp)
255 << " row-ratio= " << this->maxRowRatio(lp)
256 << std::endl;)
257 }
258 }
259
260
261 } // namespace soplex
262