1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24
25 /**@file spxequilisc.hpp
26 * @brief Equilibrium row/column scaling.
27 */
28 #include <assert.h>
29
30 #include "soplex/spxout.h"
31 #include "soplex/spxlpbase.h"
32 #include "soplex/spxlp.h"
(1) Event include_recursion: |
#include file "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/soplex/src/soplex.h" includes itself: soplex.h -> spxequilisc.h -> spxequilisc.hpp -> soplex.h |
(2) Event caretline: |
^ |
33 #include "soplex.h"
34
35 namespace soplex
36 {
37 static inline const char* makename(bool doBoth)
38 {
39 return doBoth ? "bi-Equilibrium" : "uni-Equilibrium";
40 }
41
42 /// maximum ratio between absolute biggest and smallest element in any scaled row/column.
43 template <class R>
44 static R maxPrescaledRatio(const SPxLPBase<R>& lp, const std::vector<R>& coScaleval, bool rowRatio)
45 {
46 R pmax = 0.0;
47 const int n = rowRatio ? lp.nRows() : lp.nCols();
48
49 for(int i = 0; i < n; ++i)
50 {
51 const SVectorBase<R>& vec = rowRatio ? lp.rowVector(i) : lp.colVector(i);
52 R mini = R(infinity);
53 R maxi = 0.0;
54
55 for(int j = 0; j < vec.size(); ++j)
56 {
57 assert(vec.index(j) >= 0);
58 const R x = spxAbs(vec.value(j)) * coScaleval[unsigned(vec.index(j))];
59
60 if(isZero(x, lp.tolerances()->epsilon()))
61 continue;
62
63 if(x < mini)
64 mini = x;
65
66 if(x > maxi)
67 maxi = x;
68 }
69
70 if(mini == R(infinity))
71 continue;
72
73 const R p = maxi / mini;
74
75 if(p > pmax)
76 pmax = p;
77 }
78
79 return pmax;
80 }
81
82 template <class R>
83 void SPxEquiliSC<R>::computeEquiExpVec(const SVSetBase<R>* vecset, const DataArray<int>& coScaleExp,
84 DataArray<int>& scaleExp, R epsilon)
85 {
86 assert(vecset != nullptr);
87
88 for(int i = 0; i < vecset->num(); ++i)
89 {
90 const SVectorBase<R>& vec = (*vecset)[i];
91
92 R maxi = 0.0;
93
94 for(int j = 0; j < vec.size(); ++j)
95 {
96 const R x = spxAbs(spxLdexp(vec.value(j), coScaleExp[vec.index(j)]));
97
98 if(GT(x, maxi, epsilon))
99 maxi = x;
100 }
101
102 // empty rows/cols are possible
103 if(maxi == 0.0)
104 maxi = 1.0;
105
106 assert(maxi > 0.0);
107
108 spxFrexp(Real(1.0 / maxi), &(scaleExp[i]));
109
110 scaleExp[i] -= 1;
111 }
112 }
113
114 template <class R>
115 void SPxEquiliSC<R>::computeEquiExpVec(const SVSetBase<R>* vecset, const std::vector<R>& coScaleVal,
116 DataArray<int>& scaleExp, R epsilon)
117 {
118 assert(vecset != nullptr);
119
120 for(int i = 0; i < vecset->num(); ++i)
121 {
122 const SVectorBase<R>& vec = (*vecset)[i];
123
124 R maxi = 0.0;
125
126 for(int j = 0; j < vec.size(); ++j)
127 {
128 assert(vec.index(j) >= 0);
129 const R x = spxAbs(vec.value(j) * coScaleVal[unsigned(vec.index(j))]);
130
131 if(GT(x, maxi, epsilon))
132 maxi = x;
133 }
134
135 // empty rows/cols are possible
136 if(maxi == 0.0)
137 maxi = 1.0;
138
139 assert(maxi > 0.0);
140
141 spxFrexp(Real(1.0 / maxi), &(scaleExp[i]));
142
143 scaleExp[i] -= 1;
144 }
145 }
146
147 template <class R>
148 void SPxEquiliSC<R>::computePostequiExpVecs(const SPxLPBase<R>& lp,
149 const std::vector<R>& preRowscale, const std::vector<R>& preColscale,
150 DataArray<int>& rowscaleExp, DataArray<int>& colscaleExp, R epsilon)
151 {
152 const R colratio = maxPrescaledRatio(lp, preRowscale, false);
153 const R rowratio = maxPrescaledRatio(lp, preColscale, true);
154
155 const bool colFirst = colratio < rowratio;
156
157 // see SPxEquiliSC<R>::scale for reason behind this branch
158 if(colFirst)
159 {
160 computeEquiExpVec(lp.colSet(), preRowscale, colscaleExp, epsilon);
161 computeEquiExpVec(lp.rowSet(), colscaleExp, rowscaleExp, epsilon);
162 }
163 else
164 {
165 computeEquiExpVec(lp.rowSet(), preColscale, rowscaleExp, epsilon);
166 computeEquiExpVec(lp.colSet(), rowscaleExp, colscaleExp, epsilon);
167 }
168 }
169
170 template <class R>
171 SPxEquiliSC<R>::SPxEquiliSC(bool doBoth)
172 : SPxScaler<R>(makename(doBoth), false, doBoth)
173 {}
174
175 template <class R>
176 SPxEquiliSC<R>::SPxEquiliSC(const SPxEquiliSC<R>& old)
177 : SPxScaler<R>(old)
178 {}
179
180 template <class R>
181 SPxEquiliSC<R>& SPxEquiliSC<R>::operator=(const SPxEquiliSC<R>& rhs)
182 {
183 if(this != &rhs)
184 {
185 SPxScaler<R>::operator=(rhs);
186 }
187
188 return *this;
189 }
190
191 template <class R>
192 void SPxEquiliSC<R>::scale(SPxLPBase<R>& lp, bool persistent)
193 {
194
195 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Equilibrium scaling LP" <<
196 (persistent ? " (persistent)" : "") << std::endl;)
197
198 this->setup(lp);
199
200 /* We want to do the direction first, which has a lower maximal ratio,
201 * since the lowest value in the scaled matrix is bounded from below by
202 * the inverse of the maximum ratio of the direction that is done first
203 * Example:
204 *
205 * Rowratio
206 * 0.1 1 10
207 * 10 1 10
208 *
209 * Colratio 100 1
210 *
211 * Row first => Col next =>
212 * 0.1 1 0.1 1
213 * 1 0.1 1 0.1
214 *
215 * Col first => Row next =>
216 * 0.01 1 0.01 1
217 * 1 1 1 1
218 *
219 */
220 R colratio = this->maxColRatio(lp);
221 R rowratio = this->maxRowRatio(lp);
222 R epsilon = this->tolerances()->epsilon();
223
224 bool colFirst = colratio < rowratio;
225
226 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "before scaling:"
227 << " min= " << lp.minAbsNzo()
228 << " max= " << lp.maxAbsNzo()
229 << " col-ratio= " << colratio
230 << " row-ratio= " << rowratio
231 << std::endl;)
232
233 if(colFirst)
234 {
235 computeEquiExpVec(lp.colSet(), *this->m_activeRowscaleExp, *this->m_activeColscaleExp, epsilon);
236
237 if(this->m_doBoth)
238 computeEquiExpVec(lp.rowSet(), *this->m_activeColscaleExp, *this->m_activeRowscaleExp, epsilon);
239 }
240 else
241 {
242 computeEquiExpVec(lp.rowSet(), *this->m_activeColscaleExp, *this->m_activeRowscaleExp, epsilon);
243
244 if(this->m_doBoth)
245 computeEquiExpVec(lp.colSet(), *this->m_activeRowscaleExp, *this->m_activeColscaleExp, epsilon);
246 }
247
248 /* scale */
249 this->applyScaling(lp);
250
251 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Row scaling min= " << this->minAbsRowscale()
252 << " max= " << this->maxAbsRowscale()
253 << std::endl
254 << "Col scaling min= " << this->minAbsColscale()
255 << " max= " << this->maxAbsColscale()
256 << std::endl;)
257
258 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "after scaling: "
259 << " min= " << lp.minAbsNzo(false)
260 << " max= " << lp.maxAbsNzo(false)
261 << " col-ratio= " << this->maxColRatio(lp)
262 << " row-ratio= " << this->maxRowRatio(lp)
263 << std::endl;)
264
265 }
266
267 } // namespace soplex
268