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 #include "soplex/spxsolver.h"
30 #include "soplex/exceptions.h"
31
32 namespace soplex
33 {
34 /** Initialize Vectors
35
36 Computing the right hand side vector for the feasibility vector depends on
37 the chosen representation and type of the basis.
38
39 In columnwise case, |theFvec| = \f$x_B = A_B^{-1} (- A_N x_N)\f$, where \f$x_N\f$
40 are either upper or lower bounds for the nonbasic variables (depending on
41 the variables |Status|). If these values remain unchanged throughout the
42 simplex algorithm, they may be taken directly from LP. However, in the
43 entering type algorith they are changed and, hence, retreived from the
44 column or row upper or lower bound vectors.
45
46 In rowwise case, |theFvec| = \f$\pi^T_B = (c^T - 0^T A_N) A_B^{-1}\f$. However,
47 this applies only to leaving type algorithm, where no bounds on dual
48 variables are altered. In entering type algorithm they are changed and,
49 hence, retreived from the column or row upper or lower bound vectors.
50 */
51 template <class R>
52 void SPxSolverBase<R>::computeFrhs()
53 {
54
55 if(rep() == COLUMN)
56 {
57 theFrhs->clear();
58
59 if(type() == LEAVE)
60 {
61 computeFrhsXtra();
62
63 for(int i = 0; i < this->nRows(); i++)
64 {
65 R x;
66
67 typename SPxBasisBase<R>::Desc::Status stat = this->desc().rowStatus(i);
68
69 if(!isBasic(stat))
70 {
71 // coverity[switch_selector_expr_is_constant]
72 switch(stat)
73 {
74 // columnwise cases:
75 case SPxBasisBase<R>::Desc::P_FREE :
76 continue;
77
78 case(SPxBasisBase<R>::Desc::P_FIXED) :
79 assert(EQ(this->lhs(i), this->rhs(i), this->epsilon()));
80
81 //lint -fallthrough
82 case SPxBasisBase<R>::Desc::P_ON_UPPER :
83 x = this->rhs(i);
84 break;
85
86 case SPxBasisBase<R>::Desc::P_ON_LOWER :
87 x = this->lhs(i);
88 break;
89
90 default:
91 SPX_MSG_ERROR(std::cerr << "ESVECS01 ERROR: "
92 << "inconsistent basis must not happen!"
93 << std::endl;)
94 throw SPxInternalCodeException("XSVECS01 This should never happen.");
95 }
96
97 assert(x < R(infinity));
98 assert(x > R(-infinity));
99 (*theFrhs)[i] += x; // slack !
100 }
101 }
102 }
103 else
104 {
105 computeFrhs1(*theUbound, *theLbound);
106 computeFrhs2(*theCoUbound, *theCoLbound);
107 }
108 }
109 else
110 {
111 assert(rep() == ROW);
112
113 if(type() == ENTER)
114 {
115 theFrhs->clear();
116 computeFrhs1(*theUbound, *theLbound);
117 computeFrhs2(*theCoUbound, *theCoLbound);
118 *theFrhs += this->maxObj();
119 }
120 else
121 {
122 ///@todo put this into a separate method
123 *theFrhs = this->maxObj();
124 const typename SPxBasisBase<R>::Desc& ds = this->desc();
125
126 for(int i = 0; i < this->nRows(); ++i)
127 {
128 typename SPxBasisBase<R>::Desc::Status stat = ds.rowStatus(i);
129
130 if(!isBasic(stat))
131 {
132 R x;
133
134 switch(stat)
135 {
136 case SPxBasisBase<R>::Desc::D_FREE :
137 continue;
138
139 case SPxBasisBase<R>::Desc::D_ON_UPPER :
140 case SPxBasisBase<R>::Desc::D_ON_LOWER :
141 case(SPxBasisBase<R>::Desc::D_ON_BOTH) :
142 x = this->maxRowObj(i);
143 break;
144
145 default:
146 assert(this->lhs(i) <= R(-infinity) && this->rhs(i) >= R(infinity));
147 x = 0.0;
148 break;
149 }
150
151 assert(x < R(infinity));
152 assert(x > R(-infinity));
153 // assert(x == 0.0);
154
155 if(x != 0.0)
156 theFrhs->multAdd(x, vector(i));
157 }
158 }
159 }
160 }
161 }
162
163 template <class R>
164 void SPxSolverBase<R>::computeFrhsXtra()
165 {
166
167 assert(rep() == COLUMN);
168 assert(type() == LEAVE);
169
170 for(int i = 0; i < this->nCols(); ++i)
171 {
172 typename SPxBasisBase<R>::Desc::Status stat = this->desc().colStatus(i);
173
174 if(!isBasic(stat))
175 {
176 R x;
177
178 // coverity[switch_selector_expr_is_constant]
179 switch(stat)
180 {
181 // columnwise cases:
182 case SPxBasisBase<R>::Desc::P_FREE :
183 continue;
184
185 case(SPxBasisBase<R>::Desc::P_FIXED) :
186 assert(EQ(SPxLPBase<R>::lower(i), SPxLPBase<R>::upper(i), this->epsilon()));
187
188 //lint -fallthrough
189 case SPxBasisBase<R>::Desc::P_ON_UPPER :
190 x = SPxLPBase<R>::upper(i);
191 break;
192
193 case SPxBasisBase<R>::Desc::P_ON_LOWER :
194 x = SPxLPBase<R>::lower(i);
195 break;
196
197 default:
198 SPX_MSG_ERROR(std::cerr << "ESVECS02 ERROR: "
199 << "inconsistent basis must not happen!"
200 << std::endl;)
201 throw SPxInternalCodeException("XSVECS02 This should never happen.");
202 }
203
204 assert(x < R(infinity));
205 assert(x > R(-infinity));
206
207 if(x != 0.0)
208 theFrhs->multAdd(-x, vector(i));
209 }
210 }
211 }
212
213
214 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as
215 specified by the |Status| of all nonbasic variables. The values of \f$x_N\f$ or
216 \f$\pi_N\f$ are taken from the passed arrays.
217 */
218 template <class R>
219 void SPxSolverBase<R>::computeFrhs1(
220 const VectorBase<R>& ufb, ///< upper feasibility bound for variables
221 const VectorBase<R>& lfb) ///< lower feasibility bound for variables
222 {
223
224 const typename SPxBasisBase<R>::Desc& ds = this->desc();
225
226 for(int i = 0; i < coDim(); ++i)
227 {
228 typename SPxBasisBase<R>::Desc::Status stat = ds.status(i);
229
230 if(!isBasic(stat))
231 {
232 R x;
233
234 // coverity[switch_selector_expr_is_constant]
235 switch(stat)
236 {
237 case SPxBasisBase<R>::Desc::D_FREE :
238 case SPxBasisBase<R>::Desc::D_UNDEFINED :
239 case SPxBasisBase<R>::Desc::P_FREE :
240 continue;
241
242 case SPxBasisBase<R>::Desc::P_ON_UPPER :
243 case SPxBasisBase<R>::Desc::D_ON_UPPER :
244 x = ufb[i];
245 break;
246
247 case SPxBasisBase<R>::Desc::P_ON_LOWER :
248 case SPxBasisBase<R>::Desc::D_ON_LOWER :
249 x = lfb[i];
250 break;
251
252 case(SPxBasisBase<R>::Desc::P_FIXED) :
253 assert(EQ(lfb[i], ufb[i], this->epsilon()));
254
255 //lint -fallthrough
256 case(SPxBasisBase<R>::Desc::D_ON_BOTH) :
257 x = lfb[i];
258 break;
259
260 default:
261 SPX_MSG_ERROR(std::cerr << "ESVECS03 ERROR: "
262 << "inconsistent basis must not happen!"
263 << std::endl;)
264 throw SPxInternalCodeException("XSVECS04 This should never happen.");
265 }
266
267 assert(x < R(infinity));
268 assert(x > R(-infinity));
269
270 if(x != 0.0)
271 theFrhs->multAdd(-x, vector(i));
272 }
273 }
274 }
275
276 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as
277 specified by the |Status| of all nonbasic variables. The values of
278 \f$x_N\f$ or \f$\pi_N\f$ are taken from the passed arrays.
279 */
280 template <class R>
281 void SPxSolverBase<R>::computeFrhs2(
282 VectorBase<R>& coufb, ///< upper feasibility bound for covariables
283 VectorBase<R>& colfb) ///< lower feasibility bound for covariables
284 {
285 const typename SPxBasisBase<R>::Desc& ds = this->desc();
286
287 for(int i = 0; i < dim(); ++i)
288 {
289 typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(i);
290
291 if(!isBasic(stat))
292 {
293 R x;
294
295 // coverity[switch_selector_expr_is_constant]
296 switch(stat)
297 {
298 case SPxBasisBase<R>::Desc::D_FREE :
299 case SPxBasisBase<R>::Desc::D_UNDEFINED :
300 case SPxBasisBase<R>::Desc::P_FREE :
301 continue;
302
303 case SPxBasisBase<R>::Desc::P_ON_LOWER : // negative slack bounds!
304 case SPxBasisBase<R>::Desc::D_ON_UPPER :
305 x = coufb[i];
306 break;
307
308 case SPxBasisBase<R>::Desc::P_ON_UPPER : // negative slack bounds!
309 case SPxBasisBase<R>::Desc::D_ON_LOWER :
310 x = colfb[i];
311 break;
312
313 case SPxBasisBase<R>::Desc::P_FIXED :
314 case SPxBasisBase<R>::Desc::D_ON_BOTH :
315
316 if(colfb[i] != coufb[i])
317 {
318 SPX_MSG_WARNING((*this->spxout),
319 (*this->spxout) << "WSVECS04 Frhs2[" << i << "]: " << static_cast<int>
320 (stat) << " "
321 << colfb[i] << " " << coufb[i]
322 << " shouldn't be" << std::endl;)
323
324 if(isZero(colfb[i], this->epsilon()) || isZero(coufb[i], this->epsilon()))
325 colfb[i] = coufb[i] = 0.0;
326 else
327 {
328 R mid = (colfb[i] + coufb[i]) / 2.0;
329 colfb[i] = coufb[i] = mid;
330 }
331 }
332
333 assert(EQ(colfb[i], coufb[i], this->epsilon()));
334 x = colfb[i];
335 break;
336
337 default:
338 SPX_MSG_ERROR(std::cerr << "ESVECS05 ERROR: "
339 << "inconsistent basis must not happen!"
340 << std::endl;)
341 throw SPxInternalCodeException("XSVECS05 This should never happen.");
342 }
343
344 assert(x < R(infinity));
345 assert(x > R(-infinity));
346
347 (*theFrhs)[i] -= x; // This is a slack, so no need to multiply a vector.
348 }
349 }
350 }
351
352 /** Computing the right hand side vector for |theCoPvec| depends on
353 the type of the simplex algorithm. In entering algorithms, the
354 values are taken from the inequality's right handside or the
355 column's objective value.
356
357 In contrast to this leaving algorithms take the values from vectors
358 |theURbound| and so on.
359
360 We reflect this difference by providing two pairs of methods
361 |computeEnterCoPrhs(n, stat)| and |computeLeaveCoPrhs(n, stat)|. The first
362 pair operates for entering algorithms, while the second one is intended for
363 leaving algorithms. The return value of these methods is the right hand
364 side value for the \f$n\f$-th row or column id, respectively, if it had the
365 passed |Status| for both.
366
367 Both methods are again split up into two methods named |...4Row(i,n)| and
368 |...4Col(i,n)|, respectively. They do their job for the |i|-th basis
369 variable, being the |n|-th row or column.
370 */
371 template <class R>
372 void SPxSolverBase<R>::computeEnterCoPrhs4Row(int i, int n)
373 {
374 assert(this->baseId(i).isSPxRowId());
375 assert(this->number(SPxRowId(this->baseId(i))) == n);
376
377 switch(this->desc().rowStatus(n))
378 {
379 // rowwise representation:
380 case SPxBasisBase<R>::Desc::P_FIXED :
381 assert(this->lhs(n) > R(-infinity));
382 assert(EQ(this->rhs(n), this->lhs(n), this->epsilon()));
383
384 //lint -fallthrough
385 case SPxBasisBase<R>::Desc::P_ON_UPPER :
386 assert(rep() == ROW);
387 assert(this->rhs(n) < R(infinity));
388 (*theCoPrhs)[i] = this->rhs(n);
389 break;
390
391 case SPxBasisBase<R>::Desc::P_ON_LOWER :
392 assert(rep() == ROW);
393 assert(this->lhs(n) > R(-infinity));
394 (*theCoPrhs)[i] = this->lhs(n);
395 break;
396
397 // columnwise representation:
398 // slacks must be left 0!
399 default:
400 (*theCoPrhs)[i] = this->maxRowObj(n);
401 break;
402 }
403 }
404
405 template <class R>
406 void SPxSolverBase<R>::computeEnterCoPrhs4Col(int i, int n)
407 {
408 assert(this->baseId(i).isSPxColId());
409 assert(this->number(SPxColId(this->baseId(i))) == n);
410
411 switch(this->desc().colStatus(n))
412 {
413 // rowwise representation:
414 case SPxBasisBase<R>::Desc::P_FIXED :
415 assert(EQ(SPxLPBase<R>::upper(n), SPxLPBase<R>::lower(n), this->epsilon()));
416 assert(SPxLPBase<R>::lower(n) > R(-infinity));
417
418 //lint -fallthrough
419 case SPxBasisBase<R>::Desc::P_ON_UPPER :
420 assert(rep() == ROW);
421 assert(SPxLPBase<R>::upper(n) < R(infinity));
422 (*theCoPrhs)[i] = SPxLPBase<R>::upper(n);
423 break;
424
425 case SPxBasisBase<R>::Desc::P_ON_LOWER :
426 assert(rep() == ROW);
427 assert(SPxLPBase<R>::lower(n) > R(-infinity));
428 (*theCoPrhs)[i] = SPxLPBase<R>::lower(n);
429 break;
430
431 // columnwise representation:
432 case SPxBasisBase<R>::Desc::D_UNDEFINED :
433 case SPxBasisBase<R>::Desc::D_ON_BOTH :
434 case SPxBasisBase<R>::Desc::D_ON_UPPER :
435 case SPxBasisBase<R>::Desc::D_ON_LOWER :
436 case SPxBasisBase<R>::Desc::D_FREE :
437 assert(rep() == COLUMN);
438 (*theCoPrhs)[i] = this->maxObj(n);
439 break;
440
441 default: // variable left 0
442 (*theCoPrhs)[i] = 0;
443 break;
444 }
445 }
446
447 template <class R>
448 void SPxSolverBase<R>::computeEnterCoPrhs()
449 {
450 assert(type() == ENTER);
451
452 for(int i = dim() - 1; i >= 0; --i)
453 {
454 SPxId l_id = this->baseId(i);
455
456 if(l_id.isSPxRowId())
457 computeEnterCoPrhs4Row(i, this->number(SPxRowId(l_id)));
458 else
459 computeEnterCoPrhs4Col(i, this->number(SPxColId(l_id)));
460 }
461 }
462
463 template <class R>
464 void SPxSolverBase<R>::computeLeaveCoPrhs4Row(int i, int n)
465 {
466 assert(this->baseId(i).isSPxRowId());
467 assert(this->number(SPxRowId(this->baseId(i))) == n);
468
469 switch(this->desc().rowStatus(n))
470 {
471 case SPxBasisBase<R>::Desc::D_ON_BOTH :
472 case SPxBasisBase<R>::Desc::P_FIXED :
473 assert(theLRbound[n] > R(-infinity));
474 assert(EQ(theURbound[n], theLRbound[n], this->epsilon()));
475
476 //lint -fallthrough
477 case SPxBasisBase<R>::Desc::D_ON_UPPER :
478 case SPxBasisBase<R>::Desc::P_ON_UPPER :
479 assert(theURbound[n] < R(infinity));
480 (*theCoPrhs)[i] = theURbound[n];
481 break;
482
483 case SPxBasisBase<R>::Desc::D_ON_LOWER :
484 case SPxBasisBase<R>::Desc::P_ON_LOWER :
485 assert(theLRbound[n] > R(-infinity));
486 (*theCoPrhs)[i] = theLRbound[n];
487 break;
488
489 default:
490 (*theCoPrhs)[i] = this->maxRowObj(n);
491 break;
492 }
493 }
494
495 template <class R>
496 void SPxSolverBase<R>::computeLeaveCoPrhs4Col(int i, int n)
497 {
498 assert(this->baseId(i).isSPxColId());
499 assert(this->number(SPxColId(this->baseId(i))) == n);
500
501 switch(this->desc().colStatus(n))
502 {
503 case SPxBasisBase<R>::Desc::D_UNDEFINED :
504 case SPxBasisBase<R>::Desc::D_ON_BOTH :
505 case SPxBasisBase<R>::Desc::P_FIXED :
506 assert(theLCbound[n] > R(-infinity));
507 assert(EQ(theUCbound[n], theLCbound[n], this->epsilon()));
508
509 //lint -fallthrough
510 case SPxBasisBase<R>::Desc::D_ON_LOWER :
511 case SPxBasisBase<R>::Desc::P_ON_UPPER :
512 assert(theUCbound[n] < R(infinity));
513 (*theCoPrhs)[i] = theUCbound[n];
514 break;
515
516 case SPxBasisBase<R>::Desc::D_ON_UPPER :
517 case SPxBasisBase<R>::Desc::P_ON_LOWER :
518 assert(theLCbound[n] > R(-infinity));
519 (*theCoPrhs)[i] = theLCbound[n];
520 break;
521
522 default:
523 (*theCoPrhs)[i] = this->maxObj(n);
524 // (*theCoPrhs)[i] = 0;
525 break;
526 }
527 }
528
529 template <class R>
530 void SPxSolverBase<R>::computeLeaveCoPrhs()
531 {
532 assert(type() == LEAVE);
533
534 for(int i = dim() - 1; i >= 0; --i)
535 {
536 SPxId l_id = this->baseId(i);
537
538 if(l_id.isSPxRowId())
539 computeLeaveCoPrhs4Row(i, this->number(SPxRowId(l_id)));
540 else
541 computeLeaveCoPrhs4Col(i, this->number(SPxColId(l_id)));
542 }
543 }
544
545
546 /** When computing the copricing vector, we expect the pricing vector to be
547 setup correctly. Then computing the copricing vector is nothing but
548 computing all scalar products of the pricing vector with the vectors of the
549 LPs constraint matrix.
550 */
551 template <class R>
552 void SPxSolverBase<R>::computePvec()
553 {
554 int i;
555
556 for(i = coDim() - 1; i >= 0; --i)
557 (*thePvec)[i] = vector(i) * (*theCoPvec);
558 }
559
560 template <class R>
561 void SPxSolverBase<R>::setupPupdate(void)
562 {
563 SSVectorBase<R>& p = thePvec->delta();
564 SSVectorBase<R>& c = theCoPvec->delta();
565
566 if(c.isSetup())
567 {
568 if(c.size() < 0.95 * theCoPvec->dim())
569 p.assign2product4setup(*thecovectors, c,
570 multTimeSparse, multTimeFull,
571 multSparseCalls, multFullCalls);
572 else
573 {
574 multTimeColwise->start();
575 p.assign2product(c, *thevectors);
576 multTimeColwise->stop();
577 ++multColwiseCalls;
578 }
579 }
580 else
581 {
582 multTimeUnsetup->start();
583 p.assign2productAndSetup(*thecovectors, c);
584 multTimeUnsetup->stop();
585 ++multUnsetupCalls;
586 }
587
588 p.setup();
589 }
590
591 template <class R>
592 void SPxSolverBase<R>::doPupdate(void)
593 {
594 theCoPvec->update();
595
596 if(pricing() == FULL)
597 {
598 // thePvec->delta().setup();
599 thePvec->update();
600 }
601 }
602 } // namespace soplex
603