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