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 spxleastsqsc.hpp 26 * @brief LP least squares scaling. 27 */ 28 29 #include <assert.h> 30 #include <cmath> 31 #include "soplex/spxout.h" 32 #include "soplex/basevectors.h" 33 #include "soplex/svsetbase.h" 34 #include "soplex/svectorbase.h" 35 #include "soplex/ssvectorbase.h" 36 #include <array> 37 38 namespace soplex 39 { 40 41 /* update scaling VectorBase<R> */ 42 template <class R> 43 static void updateScale( 44 const SSVectorBase<R> vecnnzeroes, 45 const SSVectorBase<R> resnvec, 46 SSVectorBase<R>& tmpvec, 47 SSVectorBase<R>*& psccurr, 48 SSVectorBase<R>*& pscprev, 49 R qcurr, 50 R qprev, 51 R eprev1, 52 R eprev2, 53 R epsilon) 54 { 55 assert(psccurr != NULL); 56 assert(pscprev != NULL); 57 assert(qcurr * qprev != 0.0); 58 59 R fac = -(eprev1 * eprev2); 60 61 SSVectorBase<R>* pssv; 62 63 *pscprev -= *psccurr; 64 65 if(isZero(fac, epsilon)) 66 (*pscprev).clear(); 67 else 68 *pscprev *= fac; 69 70 *pscprev += tmpvec.assignPWproduct4setup(resnvec, vecnnzeroes); 71 72 *pscprev *= 1.0 / (qcurr * qprev); 73 *pscprev += *psccurr; 74 75 /* swap pointers */ 76 pssv = psccurr; 77 psccurr = pscprev; 78 pscprev = pssv; 79 } 80 81 /* update scaling VectorBase<R> after main loop */ 82 template <class R> 83 static void updateScaleFinal( 84 const SSVectorBase<R> vecnnzeroes, 85 const SSVectorBase<R> resnvec, 86 SSVectorBase<R>& tmpvec, 87 SSVectorBase<R>*& psccurr, 88 SSVectorBase<R>*& pscprev, 89 R q, 90 R eprev1, 91 R eprev2, 92 R epsilon) 93 { 94 assert(q != 0); 95 assert(psccurr != NULL); 96 assert(pscprev != NULL); 97 98 R fac = -(eprev1 * eprev2); 99 100 *pscprev -= *psccurr; 101 102 if(isZero(fac, epsilon)) 103 (*pscprev).clear(); 104 else 105 *pscprev *= fac; 106 107 *pscprev += tmpvec.assignPWproduct4setup(resnvec, vecnnzeroes); 108 *pscprev *= 1.0 / q; 109 *pscprev += *psccurr; 110 111 psccurr = pscprev; 112 } 113 114 /* update residual VectorBase<R> */ 115 template <class R> 116 static inline void updateRes( 117 const SVSetBase<R> facset, 118 const SSVectorBase<R> resvecprev, 119 SSVectorBase<R>& resvec, 120 SSVectorBase<R>& tmpvec, 121 R eprev, 122 R qcurr, 123 R epsilon) 124 { 125 assert(qcurr != 0.0); 126 127 if(isZero(eprev, epsilon)) 128 resvec.clear(); 129 else 130 resvec *= eprev; 131 132 int dummy1 = 0; 133 int dummy2 = 0; 134 tmpvec.assign2product4setup(facset, resvecprev, 0, 0, dummy1, dummy2); 135 tmpvec.setup(); 136 resvec += tmpvec; 137 138 resvec *= (-1.0 / qcurr); 139 resvec.setup(); 140 } 141 142 143 /* initialize constant vectors and matrices */ 144 template <class R> 145 static void initConstVecs( 146 const SVSetBase<R>* vecset, 147 SVSetBase<R>& facset, 148 SSVectorBase<R>& veclogs, 149 SSVectorBase<R>& vecnnzinv, 150 R epsilon) 151 { 152 assert(vecset != NULL); 153 154 const int nvec = vecset->num(); 155 156 for(int k = 0; k < nvec; ++k) 157 { 158 R logsum = 0.0; 159 int nnz = 0; 160 // get kth row or column of LP 161 const SVectorBase<R>& lpvec = (*vecset)[k]; 162 const int size = lpvec.size(); 163 164 for(int i = 0; i < size; ++i) 165 { 166 const R a = lpvec.value(i); 167 168 if(!isZero(a, epsilon)) 169 { 170 logsum += log2(double(spxAbs(a))); // todo spxLog2? 171 nnz++; 172 } 173 } 174 175 R nnzinv; 176 177 if(nnz > 0) 178 { 179 nnzinv = 1.0 / nnz; 180 } 181 else 182 { 183 /* all-0 entries */ 184 logsum = 1.0; 185 nnzinv = 1.0; 186 } 187 188 veclogs.add(k, logsum); 189 vecnnzinv.add(k, nnzinv); 190 191 /* create new VectorBase<R> for facset */ 192 SVectorBase<R>& vecnew = (*(facset.create(nnz))); 193 194 for(int i = 0; i < size; ++i) 195 { 196 if(!isZero(lpvec.value(i), epsilon)) 197 vecnew.add(lpvec.index(i), nnzinv); 198 } 199 200 vecnew.sort(); 201 } 202 203 assert(veclogs.isSetup()); 204 assert(vecnnzinv.isSetup()); 205 } 206 207 /* return name of scaler */ 208 static inline const char* makename() 209 { 210 return "Least squares"; 211 } 212 213 template <class R> 214 SPxLeastSqSC<R>::SPxLeastSqSC() 215 : SPxScaler<R>(makename(), false, false) 216 {} 217 218 template <class R> 219 SPxLeastSqSC<R>::SPxLeastSqSC(const SPxLeastSqSC<R>& old) 220 : SPxScaler<R>(old), acrcydivisor(old.acrcydivisor), maxrounds(old.maxrounds) 221 {} 222 223 template <class R> 224 SPxLeastSqSC<R>& SPxLeastSqSC<R>::operator=(const SPxLeastSqSC<R>& rhs) 225 { 226 if(this != &rhs) 227 { 228 SPxScaler<R>::operator=(rhs); 229 } 230 231 return *this; 232 } 233 234 template <class R> 235 void SPxLeastSqSC<R>::setRealParam(R param, const char* name) 236 { 237 assert(param >= 1.0); 238 acrcydivisor = param; 239 } 240 241 template <class R> 242 void SPxLeastSqSC<R>::setIntParam(int param, const char* name) 243 { 244 assert(param >= 0); 245 maxrounds = param; 246 } 247 248 // todo refactor this method. Has no abstraction and is too long 249 template <class R> 250 void SPxLeastSqSC<R>::scale(SPxLPBase<R>& lp, bool persistent) 251 { 252 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Least squares LP scaling" << 253 (persistent ? " (persistent)" : "") << std::endl;) 254 255 this->setup(lp); 256 257 const int nrows = lp.nRows(); 258 const int ncols = lp.nCols(); 259 const int lpnnz = lp.nNzos(); 260 261 // is contraints matrix empty? 262 //todo don't create the scaler in this case! 263 if(lpnnz == 0) 264 { 265 // to keep the invariants, we still need to call this method 266 this->applyScaling(lp); 267 268 return; 269 } 270 271 assert(nrows > 0 && ncols > 0 && lpnnz > 0); 272 273 /* constant factor matrices; 274 * in Curtis-Reid article 275 * facnrows equals E^T M^(-1) 276 * facncols equals E N^(-1) 277 * */ 278 SVSetBase<R> facnrows(nrows, nrows, 1.1, 1.2); 279 SVSetBase<R> facncols(ncols, ncols, 1.1, 1.2); 280 281 /* column scaling factor vectors */ 282 SSVectorBase<R> colscale1(ncols, this->_tolerances); 283 SSVectorBase<R> colscale2(ncols, this->_tolerances); 284 285 /* row scaling factor vectors */ 286 SSVectorBase<R> rowscale1(nrows, this->_tolerances); 287 SSVectorBase<R> rowscale2(nrows, this->_tolerances); 288 289 /* residual vectors */ 290 SSVectorBase<R> resnrows(nrows, this->_tolerances); 291 SSVectorBase<R> resncols(ncols, this->_tolerances); 292 293 /* vectors to store temporary values */ 294 SSVectorBase<R> tmprows(nrows, this->_tolerances); 295 SSVectorBase<R> tmpcols(ncols, this->_tolerances); 296 297 /* vectors storing the row and column sums (respectively) of logarithms of 298 *(absolute values of) non-zero elements of left hand matrix of LP 299 */ 300 SSVectorBase<R> rowlogs(nrows, this->_tolerances); 301 SSVectorBase<R> collogs(ncols, this->_tolerances); 302 303 /* vectors storing the inverted number of non-zeros in each row and column 304 *(respectively) of left hand matrix of LP 305 */ 306 SSVectorBase<R> rownnzinv(nrows, this->_tolerances); 307 SSVectorBase<R> colnnzinv(ncols, this->_tolerances); 308 309 /* VectorBase<R> pointers */ 310 SSVectorBase<R>* csccurr = &colscale1; 311 SSVectorBase<R>* cscprev = &colscale2; 312 SSVectorBase<R>* rsccurr = &rowscale1; 313 SSVectorBase<R>* rscprev = &rowscale2; 314 315 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "before scaling:" 316 << " min= " << lp.minAbsNzo() 317 << " max= " << lp.maxAbsNzo() 318 << " col-ratio= " << this->maxColRatio(lp) 319 << " row-ratio= " << this->maxRowRatio(lp) 320 << std::endl;) 321 322 /* initialize scalars, vectors and matrices */ 323 324 assert(acrcydivisor > 0.0); 325 326 const R smax = lpnnz / acrcydivisor; 327 R qcurr = 1.0; 328 R qprev = 0.0; 329 330 std::array<R, 3> eprev; 331 eprev.fill(0.0); 332 333 initConstVecs(lp.rowSet(), facnrows, rowlogs, rownnzinv, R(this->tolerances()->epsilon())); 334 initConstVecs(lp.colSet(), facncols, collogs, colnnzinv, R(this->tolerances()->epsilon())); 335 336 assert(tmprows.isSetup() && tmpcols.isSetup()); 337 assert(rowscale1.isSetup() && rowscale2.isSetup()); 338 assert(colscale1.isSetup() && colscale2.isSetup()); 339 340 // compute first residual vector r0 341 int dummy1 = 0; 342 int dummy2 = 0; 343 resncols = collogs - tmpcols.assign2product4setup(facnrows, rowlogs, 0, 0, dummy1, dummy2); 344 345 resncols.setup(); 346 resnrows.setup(); 347 348 rowscale1.assignPWproduct4setup(rownnzinv, rowlogs); 349 rowscale2 = rowscale1; 350 351 R scurr = resncols * tmpcols.assignPWproduct4setup(colnnzinv, resncols); 352 353 int k; 354 355 /* conjugate gradient loop */ 356 for(k = 0; k < maxrounds; ++k) 357 { 358 const R sprev = scurr; 359 360 // termination criterion met? 361 if(scurr < smax) 362 break; 363 364 // is k even? 365 if((k % 2) == 0) 366 { 367 // not in first iteration? 368 if(k != 0) // true, then update row scaling factor vector 369 updateScale(rownnzinv, resnrows, tmprows, rsccurr, rscprev, qcurr, qprev, eprev[1], eprev[2], 370 R(this->tolerances()->epsilon())); 371 372 updateRes(facncols, resncols, resnrows, tmprows, eprev[0], qcurr, R(this->tolerances()->epsilon())); 373 scurr = resnrows * tmprows.assignPWproduct4setup(resnrows, rownnzinv); 374 } 375 else // k is odd 376 { 377 // update column scaling factor vector 378 updateScale(colnnzinv, resncols, tmpcols, csccurr, cscprev, qcurr, qprev, eprev[1], eprev[2], 379 R(this->tolerances()->epsilon())); 380 381 updateRes(facnrows, resnrows, resncols, tmpcols, eprev[0], qcurr, R(this->tolerances()->epsilon())); 382 scurr = resncols * tmpcols.assignPWproduct4setup(resncols, colnnzinv); 383 } 384 385 // shift eprev entries one to the right 386 for(unsigned l = 2; l > 0; --l) 387 eprev[l] = eprev[l - 1]; 388 389 assert(isNotZero(sprev, R(this->tolerances()->epsilon()))); 390 391 eprev[0] = (qcurr * scurr) / sprev; 392 393 const R tmp = qcurr; 394 qcurr = 1.0 - eprev[0]; 395 qprev = tmp; 396 } 397 398 if(k > 0 && (k % 2) == 0) 399 { 400 // update column scaling factor vector 401 updateScaleFinal(colnnzinv, resncols, tmpcols, csccurr, cscprev, qprev, eprev[1], eprev[2], 402 R(this->tolerances()->epsilon())); 403 } 404 else if(k > 0) 405 { 406 // update row scaling factor vector 407 updateScaleFinal(rownnzinv, resnrows, tmprows, rsccurr, rscprev, qprev, eprev[1], eprev[2], 408 R(this->tolerances()->epsilon())); 409 } 410 411 /* compute actual scaling factors */ 412 413 const SSVectorBase<R>& rowscale = *rsccurr; 414 const SSVectorBase<R>& colscale = *csccurr; 415 416 DataArray<int>& colscaleExp = *this->m_activeColscaleExp; 417 DataArray<int>& rowscaleExp = *this->m_activeRowscaleExp; 418 419 for(k = 0; k < nrows; ++k) 420 rowscaleExp[k] = -int(rowscale[k] + ((rowscale[k] >= 0.0) ? (+0.5) : (-0.5))); 421 422 for(k = 0; k < ncols; ++k) 423 colscaleExp[k] = -int(colscale[k] + ((colscale[k] >= 0.0) ? (+0.5) : (-0.5))); 424 425 // scale 426 this->applyScaling(lp); 427 428 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Row scaling min= " << this->minAbsRowscale() 429 << " max= " << this->maxAbsRowscale() 430 << std::endl 431 << "Col scaling min= " << this->minAbsColscale() 432 << " max= " << this->maxAbsColscale() 433 << std::endl;) 434 435 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "after scaling: " 436 << " min= " << lp.minAbsNzo(false) 437 << " max= " << lp.maxAbsNzo(false) 438 << " col-ratio= " << this->maxColRatio(lp) 439 << " row-ratio= " << this->maxRowRatio(lp) 440 << std::endl;) 441 } 442 443 } // namespace soplex 444