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