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" 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