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 #include <assert.h> 26 #include <iostream> 27 28 #include "soplex/spxdefines.h" 29 30 namespace soplex 31 { 32 /** 33 * Here comes the ratio test for selecting a variable to leave the basis. 34 * It is assumed that Vec.delta() and fVec.idx() have been setup 35 * correctly! 36 * 37 * The leaving variable is selected such that the update of fVec() (using 38 * fVec.value() * fVec.delta()) keeps the basis feasible within 39 * solver()->entertol(). Hence, fVec.value() must be chosen such that one 40 * updated value of theFvec just reaches its bound and no other one exceeds 41 * them by more than solver()->entertol(). Further, fVec.value() must have the 42 * same sign as argument \p val. 43 * 44 * The return value of selectLeave() is the number of a variable in the 45 * basis selected to leave the basis. -1 indicates that no variable could be 46 * selected. Otherwise, parameter \p val contains the chosen fVec.value(). 47 */ 48 template <class R> 49 int SPxDefaultRT<R>::selectLeave(R& val, R, bool) 50 { 51 this->solver()->fVec().delta().setup(); 52 53 const R* vec = this->solver()->fVec().get_const_ptr(); 54 const R* upd = this->solver()->fVec().delta().values(); 55 const IdxSet& idx = this->solver()->fVec().idx(); 56 const R* ub = this->solver()->ubBound().get_const_ptr(); 57 const R* lb = this->solver()->lbBound().get_const_ptr(); 58 59 R epsilon = this->solver()->epsilon(); 60 int leave = -1; 61 62 R x; 63 int i; 64 int j; 65 66 // PARALLEL the j loop could be parallelized 67 if(val > 0) 68 { 69 // Loop over NZEs of delta vector. 70 for(j = 0; j < idx.size(); ++j) 71 { 72 i = idx.index(j); 73 x = upd[i]; 74 75 if(x > epsilon) 76 { 77 if(ub[i] < R(infinity)) 78 { 79 R y = (ub[i] - vec[i] + this->delta) / x; 80 81 if(y < val) 82 { 83 leave = i; 84 val = y; 85 } 86 } 87 } 88 else if(x < -epsilon) 89 { 90 if(lb[i] > R(-infinity)) 91 { 92 R y = (lb[i] - vec[i] - this->delta) / x; 93 94 if(y < val) 95 { 96 leave = i; 97 val = y; 98 } 99 } 100 } 101 } 102 103 if(leave >= 0) 104 { 105 x = upd[leave]; 106 107 // BH 2005-11-30: It may well happen that the basis is degenerate and the 108 // selected leaving variable is (at most this->delta) beyond its bound. (This 109 // happens for instance on LP/netlib/adlittle.mps with setting -r -t0.) 110 // In that case we do a pivot step with length zero to avoid difficulties. 111 if((x > epsilon && vec[leave] >= ub[leave]) || 112 (x < -epsilon && vec[leave] <= lb[leave])) 113 { 114 val = 0.0; 115 } 116 else 117 { 118 val = (x > epsilon) ? ub[leave] : lb[leave]; 119 val = (val - vec[leave]) / x; 120 } 121 } 122 123 SOPLEX_ASSERT_WARN("WDEFRT01", val > -epsilon); 124 } 125 else 126 { 127 for(j = 0; j < idx.size(); ++j) 128 { 129 i = idx.index(j); 130 x = upd[i]; 131 132 if(x < -epsilon) 133 { 134 if(ub[i] < R(infinity)) 135 { 136 R y = (ub[i] - vec[i] + this->delta) / x; 137 138 if(y > val) 139 { 140 leave = i; 141 val = y; 142 } 143 } 144 } 145 else if(x > epsilon) 146 { 147 if(lb[i] > R(-infinity)) 148 { 149 R y = (lb[i] - vec[i] - this->delta) / x; 150 151 if(y > val) 152 { 153 leave = i; 154 val = y; 155 } 156 } 157 } 158 } 159 160 if(leave >= 0) 161 { 162 x = upd[leave]; 163 164 // See comment above. 165 if((x < -epsilon && vec[leave] >= ub[leave]) || 166 (x > epsilon && vec[leave] <= lb[leave])) 167 { 168 val = 0.0; 169 } 170 else 171 { 172 val = (x < epsilon) ? ub[leave] : lb[leave]; 173 val = (val - vec[leave]) / x; 174 } 175 } 176 177 SOPLEX_ASSERT_WARN("WDEFRT02", val < epsilon); 178 } 179 180 return leave; 181 } 182 183 /** 184 Here comes the ratio test. It is assumed that theCoPvec.this->delta() and 185 theCoPvec.idx() have been setup correctly! 186 */ 187 template <class R> 188 SPxId SPxDefaultRT<R>::selectEnter(R& max, int, bool) 189 { 190 this->solver()->coPvec().delta().setup(); 191 this->solver()->pVec().delta().setup(); 192 193 const R* pvec = this->solver()->pVec().get_const_ptr(); 194 const R* pupd = this->solver()->pVec().delta().values(); 195 const IdxSet& pidx = this->solver()->pVec().idx(); 196 const R* lpb = this->solver()->lpBound().get_const_ptr(); 197 const R* upb = this->solver()->upBound().get_const_ptr(); 198 199 const R* cvec = this->solver()->coPvec().get_const_ptr(); 200 const R* cupd = this->solver()->coPvec().delta().values(); 201 const IdxSet& cidx = this->solver()->coPvec().idx(); 202 const R* lcb = this->solver()->lcBound().get_const_ptr(); 203 const R* ucb = this->solver()->ucBound().get_const_ptr(); 204 205 R epsilon = this->solver()->epsilon(); 206 R val = max; 207 int pnum = -1; 208 int cnum = -1; 209 210 SPxId enterId; 211 int i; 212 int j; 213 R x; 214 215 // PARALLEL the j loops could be parallelized 216 if(val > 0) 217 { 218 for(j = 0; j < pidx.size(); ++j) 219 { 220 i = pidx.index(j); 221 x = pupd[i]; 222 223 if(x > epsilon) 224 { 225 if(upb[i] < R(infinity)) 226 { 227 R y = (upb[i] - pvec[i] + this->delta) / x; 228 229 if(y < val) 230 { 231 enterId = this->solver()->id(i); 232 val = y; 233 pnum = j; 234 } 235 } 236 } 237 else if(x < -epsilon) 238 { 239 if(lpb[i] > R(-infinity)) 240 { 241 R y = (lpb[i] - pvec[i] - this->delta) / x; 242 243 if(y < val) 244 { 245 enterId = this->solver()->id(i); 246 val = y; 247 pnum = j; 248 } 249 } 250 } 251 } 252 253 for(j = 0; j < cidx.size(); ++j) 254 { 255 i = cidx.index(j); 256 x = cupd[i]; 257 258 if(x > epsilon) 259 { 260 if(ucb[i] < R(infinity)) 261 { 262 R y = (ucb[i] - cvec[i] + this->delta) / x; 263 264 if(y < val) 265 { 266 enterId = this->solver()->coId(i); 267 val = y; 268 cnum = j; 269 } 270 } 271 } 272 else if(x < -epsilon) 273 { 274 if(lcb[i] > R(-infinity)) 275 { 276 R y = (lcb[i] - cvec[i] - this->delta) / x; 277 278 if(y < val) 279 { 280 enterId = this->solver()->coId(i); 281 val = y; 282 cnum = j; 283 } 284 } 285 } 286 } 287 288 if(cnum >= 0) 289 { 290 i = cidx.index(cnum); 291 x = cupd[i]; 292 val = (x > epsilon) ? ucb[i] : lcb[i]; 293 val = (val - cvec[i]) / x; 294 } 295 else if(pnum >= 0) 296 { 297 i = pidx.index(pnum); 298 x = pupd[i]; 299 val = (x > epsilon) ? upb[i] : lpb[i]; 300 val = (val - pvec[i]) / x; 301 } 302 } 303 else 304 { 305 for(j = 0; j < pidx.size(); ++j) 306 { 307 i = pidx.index(j); 308 x = pupd[i]; 309 310 if(x > epsilon) 311 { 312 if(lpb[i] > R(-infinity)) 313 { 314 R y = (lpb[i] - pvec[i] - this->delta) / x; 315 316 if(y > val) 317 { 318 enterId = this->solver()->id(i); 319 val = y; 320 pnum = j; 321 } 322 } 323 } 324 else if(x < -epsilon) 325 { 326 if(upb[i] < R(infinity)) 327 { 328 R y = (upb[i] - pvec[i] + this->delta) / x; 329 330 if(y > val) 331 { 332 enterId = this->solver()->id(i); 333 val = y; 334 pnum = j; 335 } 336 } 337 } 338 } 339 340 for(j = 0; j < cidx.size(); ++j) 341 { 342 i = cidx.index(j); 343 x = cupd[i]; 344 345 if(x > epsilon) 346 { 347 if(lcb[i] > R(-infinity)) 348 { 349 R y = (lcb[i] - cvec[i] - this->delta) / x; 350 351 if(y > val) 352 { 353 enterId = this->solver()->coId(i); 354 val = y; 355 cnum = j; 356 } 357 } 358 } 359 else if(x < -epsilon) 360 { 361 if(ucb[i] < R(infinity)) 362 { 363 R y = (ucb[i] - cvec[i] + this->delta) / x; 364 365 if(y > val) 366 { 367 enterId = this->solver()->coId(i); 368 val = y; 369 cnum = j; 370 } 371 } 372 } 373 } 374 375 if(cnum >= 0) 376 { 377 i = cidx.index(cnum); 378 x = cupd[i]; 379 val = (x < epsilon) ? ucb[i] : lcb[i]; 380 val = (val - cvec[i]) / x; 381 } 382 else if(pnum >= 0) 383 { 384 i = pidx.index(pnum); 385 x = pupd[i]; 386 val = (x < epsilon) ? upb[i] : lpb[i]; 387 val = (val - pvec[i]) / x; 388 } 389 } 390 391 if(enterId.isValid() && this->solver()->isBasic(enterId)) 392 { 393 SPxOut::debug(this, "DDEFRT01 isValid() && isBasic(): max={}\n", max); 394 395 if(cnum >= 0) 396 this->solver()->coPvec().delta().clearNum(cnum); 397 else if(pnum >= 0) 398 this->solver()->pVec().delta().clearNum(pnum); 399 400 return SPxDefaultRT<R>::selectEnter(max, 0, false); 401 } 402 403 SPX_DEBUG( 404 405 if(!enterId.isValid()) 406 std::cout << "DDEFRT02 !isValid(): max=" << max << ", x=" << x << std::endl; 407 ) 408 max = val; 409 410 return enterId; 411 } 412 413 } // namespace soplex 414