1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2021 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 /**@file spxbasis.h
17 * @brief Simplex basis.
18 */
19 #ifndef _SPXBASIS_H_
20 #define _SPXBASIS_H_
21
22 /* undefine SOPLEX_DEBUG flag from including files; if SOPLEX_DEBUG should be defined in this file, do so below */
23 #ifdef SOPLEX_DEBUG
24 #define SOPLEX_DEBUG_SPXBASIS
25 #undef SOPLEX_DEBUG
26 #endif
27
28 #include <assert.h>
29 #include <iostream>
30 #include <iomanip>
31 #include <string.h>
32 #include <sstream>
33
34 #include "soplex/spxdefines.h"
35 #include "soplex/spxlp.h"
36 #include "soplex/svector.h"
37 #include "soplex/ssvector.h"
38 #include "soplex/dataarray.h"
39 #include "soplex/slinsolver.h"
40 #include "soplex/nameset.h"
41 #include "soplex/spxout.h"
42 #include "soplex/timerfactory.h"
43
44 //#define MEASUREUPDATETIME
45
46 namespace soplex
47 {
48 template <class R>
49 class SPxSolverBase;
50
51 /**@class SPxBasisBase
52 @brief Simplex basis.
53 @ingroup Algo
54
55 Consider the linear program as provided from class SPxLP:
56 \f[
57 \begin{array}{rl}
58 \hbox{max} & c^T x \\
59 \hbox{s.t.} & l_r \le Ax \le u_r \\
60 & l_c \le x \le u_c
61 \end{array}
62 \f]
63 where \f$c, l_c, u_c, x \in {\bf R}^n\f$, \f$l_r, u_r \in {\bf R}^m\f$ and
64 \f$A \in {\bf R}^{m \times n}\f$. Solving this LP with the simplex algorithm
65 requires the definition of a \em basis. Such can be defined as a set of
66 column vectors or a set of row vectors building a non-singular matrix. We
67 will refer to the first case as the \em columnwise \em representation and
68 the latter case will be called the \em rowwise \em representation. In both
69 cases, a \em basis is a set of vectors forming a non-singular matrix. The
70 dimension of the vectors is referred to as the basis' \em dimension,
71 whereas the number of vectors belonging to the LP is called the basis'
72 \em codimension.
73
74 Class SPxBasisBase is designed to represent a generic simplex basis, suitable
75 for both representations. At any time the representation can be changed by
76 calling method setRep().
77
78 Class SPxBasisBase provides methods for solving linear systems with the basis
79 matrix. However, SPxBasisBase does not provide a linear solver by its own.
80 Instead, a SLinSolver object must be #load%ed to a SPxBasisBase which will
81 be called for solving linear systems.
82 */
83 template <class R> // theLP gets templated
84 class SPxBasisBase
85 {
86 public:
87
88 /// basis status.
89 /** Each SPxBasisBase is assigned a status flag, which can take on of the
90 above values.
91 */
92 enum SPxStatus
93 {
94 NO_PROBLEM = -2, ///< No Problem has been loaded to the basis.
95 SINGULAR = -1, ///< Basis is singular.
96 REGULAR = 0, ///< Basis is not known to be dual nor primal feasible.
97 DUAL = 1, ///< Basis is dual feasible.
98 PRIMAL = 2, ///< Basis is primal feasible.
99 OPTIMAL = 3, ///< Basis is optimal, i.e. dual and primal feasible.
100 UNBOUNDED = 4, ///< LP has been proven to be primal unbounded.
101 INFEASIBLE = 5 ///< LP has been proven to be primal infeasible.
102 };
103
104
105 /// Basis descriptor.
106 class Desc
107 {
108 public:
109
110 //------------------------------------
111 ///@name Status
112 ///@{
113 /// Status of a variable.
114 /** A basis is described by assigning a Status to all of the LP
115 variables and covariables. This assignment is maintained by the
116 basis #Desc%riptor.
117
118 Variables and covariables (slackvariables) may have a primal or dual Status. The
119 first type specifies that a variable is set on a primal bound, while
120 the latter type indicates a dual variable to be set on a bound.
121 If a row variable has a primal status, say #P_ON_UPPER, this means
122 that the upper bound of the inequality is set to be tight. Hence,
123 in this case the upper bound must not be infinity.
124
125 Equivalently, if the status of a variable is dual, say #D_ON_UPPER,
126 it means that the dual variable corresponding to the upper bound
127 inequality of this variable is set to 0.
128
129 For a column basis, primal #Status%es correspond to nonbasic
130 variables, while dual ones are basic. This is reversed for a row
131 basis. We will now reveal in more detail the significance of
132 variable #Status%es.
133
134 <b>Primal Variables</b>
135
136 Consider a range inequality \f$l_r \le a^T x \le u_r\f$ or bounds on
137 a variable \f$l_c \le x_c \le u_c\f$. The following table reveals
138 what is implied if the corresponding variable or covariable is
139 assigned to a primal #Status:
140
141 \f[
142 \begin{array}{lcl}
143 l_c \le x_c \le u_c & \mbox{Status}(x_i) & l_r \le a^T x \le u_r \\
144 \hline
145 x_c = u_c < \infty & \mbox{P\_ON\_UPPER} & a^T x = u_r < \infty \\
146 x_c = l_c > -\infty & \mbox{P\_ON\_LOWER} & a^T x = l_r > -\infty \\
147 -\infty < l_c = x_c = u_c < \infty
148 & \mbox{P\_FIXED} &
149 -\infty < l_r = a^T x = u_r < \infty \\
150 -\infty = l_i < x_i=0 < u_i = \infty
151 & \mbox{P\_FREE} &
152 -\infty = l_r < a^T x = 0 < u_r = \infty \\
153 \end{array}
154 \f]
155
156 Note that to determine whether a variable with #Status stat is set to
157 its upper bound, one can compute the test (-stat | -#P_ON_UPPER).
158 This will yield true even if the variable is fixed, i.e., sitting on
159 both bounds at the same time.
160
161 <b>Dual Variables</b>
162
163 In principle for implementing the Simplex algorithm it would suffice
164 to use only one dual #Status. However, for performance reasons it
165 is advisable to introduce various dual status types, reflecting
166 the structure of the bounds. Given an upper bound \f$u\f$ and a lower
167 bound \f$l\f$ of a constraint or variable, the following table
168 indicates the setting of the dual Status of this variable.
169
170 \f[
171 \begin{array}{cl}
172 l \le ... \le u & \mbox{Status} \\
173 \hline
174 -\infty < l \ne u < \infty & \mbox{D\_ON\_BOTH} \\
175 -\infty < l \ne u = \infty & \mbox{D\_ON\_UPPER} \\
176 -\infty = l \ne u < \infty & \mbox{D\_ON\_LOWER} \\
177 -\infty < l = u < \infty & \mbox{D\_FREE} \\
178 -\infty = l \ne u = \infty & \mbox{D\_UNDEFINED} \\
179 \end{array}
180 \f]
181
182 Note that unbounded primal variables are reflected by an #D_UNDEFINED
183 dual variable, since no reduced costs exist for them. To facilitate
184 the assignment of dual #Status%es, class SPxBasisBase provides methods
185 #dualStatus(), #dualColStatus() and #dualRowStatus)().
186 */
187 enum Status
188 {
189 P_ON_LOWER = -4, ///< primal variable is set to its lower bound
190 P_ON_UPPER = -2, ///< primal variable is set to its upper bound
191 P_FREE = -1, ///< primal variable is left free, but unset
192 P_FIXED = P_ON_UPPER + P_ON_LOWER, ///< primal variable is fixed to both bounds
193 D_FREE = 1, ///< dual variable is left free, but unset
194 D_ON_UPPER = 2, ///< dual variable is set to its upper bound
195 D_ON_LOWER = 4, ///< dual variable is set to its lower bound
196 D_ON_BOTH = D_ON_LOWER + D_ON_UPPER, ///< dual variable has two bounds
197 D_UNDEFINED = 8 ///< primal or dual variable is undefined
198 };
199 ///@}
200
201 friend SPxBasisBase<R>;
202 template <class T> friend std::ostream& operator<< (std::ostream& os,
203 const Status& stat); //@todo is the <> required here?
204
205 private:
206
207 //------------------------------------
208 ///@name Data
209 ///@{
210 DataArray < Status > rowstat; ///< status of rows.
211 DataArray < Status > colstat; ///< status of columns.
212 DataArray < Status >* stat; ///< basis' status.
213 DataArray < Status >* costat; ///< cobasis' status.
214 ///@}
215
216 public:
217
218 //------------------------------------
219 ///@name Access / modification
220 ///@{
221 /// returns number of columns.
222 int nCols() const
223 {
224 return colstat.size();
225 }
226 /// returns number of rows.
227 int nRows() const
228 {
229 return rowstat.size();
230 }
231 /// returns dimension.
232 int dim() const
233 {
234 return stat->size();
235 }
236 /// returns codimension.
237 int coDim() const
238 {
239 return costat->size();
240 }
241 ///
242 Status& rowStatus(int i)
243 {
244 return rowstat[i];
245 }
246 /// returns status of row \p i.
247 Status rowStatus(int i) const
248 {
249 return rowstat[i];
250 }
251 /// returns the array of row \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es.
252 const Status* rowStatus(void) const
253 {
254 return rowstat.get_const_ptr();
255 }
256 ///
257 Status& colStatus(int i)
258 {
259 return colstat[i];
260 }
261 /// returns status of column \p i.
262 Status colStatus(int i) const
263 {
264 return colstat[i];
265 }
266 /// returns the array of column \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es.
267 const Status* colStatus(void) const
268 {
269 return colstat.get_const_ptr();
270 }
271 ///
272 Status& status(int i)
273 {
274 return (*stat)[i];
275 }
276 /// returns status of variable \p i.
277 Status status(int i) const
278 {
279 return (*stat)[i];
280 }
281 /// returns the array of variable \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es.
282 const Status* status(void) const
283 {
284 return stat->get_const_ptr();
285 }
286 ///
287 Status& coStatus(int i)
288 {
289 return (*costat)[i];
290 }
291 /// returns status of covariable \p i.
292 Status coStatus(int i) const
293 {
294 return (*costat)[i];
295 }
296 /// returns the array of covariable \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es.
297 const Status* coStatus(void) const
298 {
299 return costat->get_const_ptr();
300 }
301 /// resets dimensions.
302 void reSize(int rowDim, int colDim);
303 ///@}
304
305 //------------------------------------
306 ///@name Debugging
307 ///@{
308 /// Prints out status.
309 void dump() const;
310
311 /// consistency check.
312 bool isConsistent() const;
313 ///@}
314
315 //------------------------------------
316 ///@name Construction / destruction
317 ///@{
318 /// default constructor
319 Desc()
320 : stat(0)
321 , costat(0)
322 {}
323 explicit Desc(const SPxSolverBase<R>& base);
324
325 /// copy constructor
326 Desc(const Desc& old);
327 /// assignment operator
328 Desc& operator=(const Desc& rhs);
329 ///@}
330 };
331
332 protected:
333
334 //------------------------------------
335 //**@name Protected data
336 /**
337 For storing the basis matrix we keep two arrays: Array #theBaseId
338 contains the SPxId%s of the basis vectors, and #matrix the pointers to
339 the vectors themselfes. Method #loadMatrixVecs() serves for loading
340 #matrix according to the SPxId%s stored in #theBaseId. This method must
341 be called whenever the VectorBase<R> pointers may have
342 changed due to manipulations of the LP.
343 */
344 ///@{
345 /// the LP
346 SPxSolverBase<R>* theLP;
347 /// SPxId%s of basic vectors.
348 DataArray < SPxId > theBaseId;
349 /// pointers to the vectors of the basis matrix.
350 DataArray < const SVectorBase<R>* > matrix;
351 /// \c true iff the pointers in \ref soplex::SPxBasisBase<R>::matrix "matrix" are set up correctly.
352 bool matrixIsSetup;
353
354 /* @brief LU factorization of basis matrix
355 The factorization of the matrix is stored in #factor if #factorized != 0.
356 Otherwise #factor is undefined.
357 */
358 SLinSolver<R>* factor;
359 /// \c true iff \ref soplex::SPxBasisBase<R>::factor "factor" = \ref soplex::SPxBasisBase<R>::matrix "matrix" \f$^{-1}\f$.
360 bool factorized;
361
362 /// number of updates before refactorization.
363 /** When a vector of the basis matrix is exchanged by a call to method
364 #change(), the LU factorization of the matrix is updated
365 accordingly. However, after atmost #maxUpdates updates of the
366 factorization, it is recomputed in order to regain numerical
367 stability and reduce fill in.
368 */
369 int maxUpdates;
370
371 /// allowed increase of nonzeros before refactorization.
372 /** When the number of nonzeros in LU factorization exceeds
373 #nonzeroFactor times the number of nonzeros in B, the
374 basis matrix is refactorized.
375 */
376 R nonzeroFactor;
377
378 /// allowed increase in relative fill before refactorization
379 /** When the real relative fill is bigger than fillFactor times lastFill
380 * the Basis will be refactorized.
381 */
382 R fillFactor;
383
384 /// allowed total increase in memory consumption before refactorization
385 R memFactor;
386
387 /* Rank-1-updates to the basis may be performed via method #change(). In
388 this case, the factorization is updated, and the following members are
389 reset.
390 */
391 int iterCount; ///< number of calls to change() since last manipulation
392 int lastIterCount; ///< number of calls to change() before halting the simplex
393 int iterDegenCheck;///< number of calls to change() since last degeneracy check
394 int updateCount; ///< number of calls to change() since last factorize()
395 int totalUpdateCount; ///< number of updates
396 int nzCount; ///< number of nonzeros in basis matrix
397 int lastMem; ///< memory needed after last fresh factorization
398 R lastFill; ///< fill ratio that occured during last factorization
399 int lastNzCount; ///< number of nonzeros in basis matrix after last fresh factorization
400
401 Timer* theTime; ///< time spent in updates
402 Timer::TYPE timerType; ///< type of timer (user or wallclock)
403
404 SPxId lastin; ///< lastEntered(): variable entered the base last
405 SPxId lastout; ///< lastLeft(): variable left the base last
406 int lastidx; ///< lastIndex(): basis index where last update was done
407 R minStab; ///< minimum stability
408 ///@}
409
410 private:
411
412 //------------------------------------
413 //**@name Private data */
414 ///@{
415 SPxStatus thestatus; ///< current status of the basis.
416 Desc thedesc; ///< the basis' Descriptor
417 bool freeSlinSolver; ///< true iff factor should be freed inside of this object
418 SPxOut* spxout; ///< message handler
419
420 ///@}
421
422 public:
423
424 //------------------------------------------------
425 /**@name Status and Descriptor related Methods */
426 ///@{
427 /// returns current SPxStatus.
428 SPxStatus status() const
429 {
430 return thestatus;
431 }
432
433 /// sets basis SPxStatus to \p stat.
434 void setStatus(SPxStatus stat)
435 {
436
437 if(thestatus != stat)
438 {
439 #ifdef SOPLEX_DEBUG
440 MSG_DEBUG(std::cout << "DBSTAT01 SPxBasisBase<R>::setStatus(): status: "
441 << int(thestatus) << " (" << thestatus << ") -> "
442 << int(stat) << " (" << stat << ")" << std::endl;)
443 #endif
444
445 thestatus = stat;
446
447 if(stat == NO_PROBLEM)
448 invalidate();
449 }
450 }
451
452 // TODO control factorization frequency dynamically
453 /// change maximum number of iterations until a refactorization is performed
454 void setMaxUpdates(int maxUp)
455 {
456 assert(maxUp >= 0);
457 maxUpdates = maxUp;
458 }
459
460 /// returns maximum number of updates before a refactorization is performed
461 int getMaxUpdates() const
462 {
463 return maxUpdates;
464 }
465
466 ///
467 const Desc& desc() const
468 {
469 return thedesc;
470 }
471 /// returns current basis Descriptor.
472 Desc& desc()
473 {
474 return thedesc;
475 }
476
477 /// dual Status for the \p i'th column variable of the loaded LP.
478 typename Desc::Status dualColStatus(int i) const;
479
480 /// dual Status for the column variable with ID \p id of the loaded LP.
481 typename Desc::Status dualStatus(const SPxColId& id) const;
482
483 /// dual Status for the \p i'th row variable of the loaded LP.
484 typename Desc::Status dualRowStatus(int i) const;
485
486 /// dual Status for the row variable with ID \p id of the loaded LP.
487 typename Desc::Status dualStatus(const SPxRowId& id) const;
488
489 /// dual Status for the variable with ID \p id of the loaded LP.
490 /** It is automatically detected, whether the \p id is one of a
491 row or a column variable, and the correct row or column status
492 is returned.
493 */
494 typename Desc::Status dualStatus(const SPxId& id) const
495 {
496 return id.isSPxRowId()
497 ? dualStatus(SPxRowId(id))
498 : dualStatus(SPxColId(id));
499 }
500 ///@}
501
502
503 //-----------------------------------
504 /**@name Inquiry Methods */
505 ///@{
506 ///
507 inline SPxId& baseId(int i)
508 {
509 return theBaseId[i];
510 }
511 /// returns the Id of the \p i'th basis vector.
512 inline SPxId baseId(int i) const
513 {
514 return theBaseId[i];
515 }
516
517 /// returns the \p i'th basic vector.
518 const SVectorBase<R>& baseVec(int i) const
519 {
520 assert(matrixIsSetup);
521 return *matrix[i];
522 }
523
524 /// returns SPxId of last VectorBase<R> included to the basis.
525 inline SPxId lastEntered() const
526 {
527 return lastin;
528 }
529
530 /// returns SPxId of last vector that left the basis.
531 inline SPxId lastLeft() const
532 {
533 return lastout;
534 }
535
536 /// returns index in basis where last update was done.
537 inline int lastIndex() const
538 {
539 return lastidx;
540 }
541
542 /// returns number of basis changes since last refactorization.
543 inline int lastUpdate() const
544 {
545 return updateCount;
546 }
547
548 /// returns number of basis changes since last \ref soplex::SPxBasisBase<R>::load() "load()".
549 inline int iteration() const
550 {
551 return iterCount;
552 }
553
554 /// returns the number of iterations prior to the last break in execution
555 inline int prevIteration() const
556 {
557 return lastIterCount;
558 }
559
560 /// returns the number of iterations since the last degeneracy check
561 inline int lastDegenCheck() const
562 {
563 return iterDegenCheck;
564 }
565
566 /// returns loaded solver.
567 inline SPxSolverBase<R>* solver() const
568 {
569 return theLP;
570 }
571 ///@}
572
573 //-----------------------------------
574 /**@name Linear Algebra */
575 ///@{
576 /// Basis-vector product.
577 /** Depending on the representation, for an SPxBasisBase B,
578 B.multBaseWith(x) computes
579 - \f$x \leftarrow Bx\f$ in the columnwise case, and
580 - \f$x \leftarrow x^TB\f$ in the rowwise case.
581
582 Both can be seen uniformly as multiplying the basis matrix \p B with
583 a vector \p x aligned the same way as the \em vectors of \p B.
584 */
585 VectorBase<R>& multBaseWith(VectorBase<R>& x) const;
586
587 /// Basis-vector product
588 void multBaseWith(SSVectorBase<R>& x, SSVectorBase<R>& result) const;
589
590 /// Vector-basis product.
591 /** Depending on the representation, for a #SPxBasisBase B,
592 B.multWithBase(x) computes
593 - \f$x \leftarrow x^TB\f$ in the columnwise case and
594 - \f$x \leftarrow Bx\f$ in the rowwise case.
595
596 Both can be seen uniformly as multiplying the basis matrix \p B with
597 a vector \p x aligned the same way as the \em covectors of \p B.
598 */
599 VectorBase<R>& multWithBase(VectorBase<R>& x) const;
600
601 /// VectorBase<R>-basis product
602 void multWithBase(SSVectorBase<R>& x, SSVectorBase<R>& result) const;
603
604 /* compute an estimated condition number for the current basis matrix
605 * by computing estimates of the norms of B and B^-1 using the power method.
606 * maxiters and tolerance control the accuracy of the estimate.
607 */
608 R condition(int maxiters = 10, R tolerance = 1e-6);
609
610 /* wrapper to compute an estimate of the condition number of the current basis matrix */
611 R getEstimatedCondition()
612 {
613 return condition(20, 1e-6);
614 }
615
616 /* wrapper to compute the exact condition number of the current basis matrix */
617 R getExactCondition()
618 {
619 return condition(1000, 1e-9);
620 }
621
622 /** compute one of several matrix metrics based on the diagonal of the LU factorization
623 * type = 0: max/min ratio
624 * type = 1: trace of U (sum of diagonal elements)
625 * type = 2: determinant (product of diagonal elements)
626 */
627 R getMatrixMetric(int type = 0);
628
629 /// returns the stability of the basis matrix.
630 R stability() const
631 {
632 return factor->stability();
633 }
634 ///
635 void solve(VectorBase<R>& x, const VectorBase<R>& rhs)
636 {
637 if(rhs.dim() == 0)
638 {
639 x.clear();
640 return;
641 }
642
643 if(!factorized)
644 SPxBasisBase<R>::factorize();
645
646 factor->solveRight(x, rhs);
647 }
648 ///
649 void solve(SSVectorBase<R>& x, const SVectorBase<R>& rhs)
650 {
651 if(rhs.size() == 0)
652 {
653 x.clear();
654 return;
655 }
656
657 if(!factorized)
658 SPxBasisBase<R>::factorize();
659
660 factor->solveRight(x, rhs);
661 }
662 /// solves linear system with basis matrix.
663 /** Depending on the representation, for a SPxBasisBase B,
664 B.solve(x) computes
665 - \f$x \leftarrow B^{-1}rhs\f$ in the columnwise case and
666 - \f$x \leftarrow rhs^TB^{-1}\f$ in the rowwise case.
667
668 Both can be seen uniformly as solving a linear system with the basis
669 matrix \p B and a right handside vector \p x aligned the same way as
670 the \em vectors of \p B.
671 */
672 void solve4update(SSVectorBase<R>& x, const SVectorBase<R>& rhs)
673 {
674 if(rhs.size() == 0)
675 {
676 x.clear();
677 return;
678 }
679
680 if(!factorized)
681 SPxBasisBase<R>::factorize();
682
683 factor->solveRight4update(x, rhs);
684 }
685 /// solves two systems in one call.
686 void solve4update(SSVectorBase<R>& x, VectorBase<R>& y, const SVectorBase<R>& rhsx,
687 SSVectorBase<R>& rhsy)
688 {
689 if(!factorized)
690 SPxBasisBase<R>::factorize();
691
692 factor->solve2right4update(x, y, rhsx, rhsy);
693 }
694 /// solves two systems in one call using only sparse data structures
695 void solve4update(SSVectorBase<R>& x, SSVectorBase<R>& y, const SVectorBase<R>& rhsx,
696 SSVectorBase<R>& rhsy)
697 {
698 if(!factorized)
699 SPxBasisBase<R>::factorize();
700
701 factor->solve2right4update(x, y, rhsx, rhsy);
702 }
703 /// solves three systems in one call.
704 void solve4update(SSVectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& y2,
705 const SVectorBase<R>& rhsx, SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsy2)
706 {
707 if(!factorized)
708 SPxBasisBase<R>::factorize();
709
710 assert(rhsy.isSetup());
711 assert(rhsy2.isSetup());
712 factor->solve3right4update(x, y, y2, rhsx, rhsy, rhsy2);
713 }
714 /// solves three systems in one call using only sparse data structures
715 void solve4update(SSVectorBase<R>& x, SSVectorBase<R>& y, SSVectorBase<R>& y2,
716 const SVectorBase<R>& rhsx, SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsy2)
717 {
718 if(!factorized)
719 SPxBasisBase<R>::factorize();
720
721 assert(rhsy.isSetup());
722 assert(rhsy2.isSetup());
723 factor->solve3right4update(x, y, y2, rhsx, rhsy, rhsy2);
724 }
725 /// Cosolves linear system with basis matrix.
726 /** Depending on the representation, for a SPxBasisBase B,
727 B.coSolve(x) computes
728 - \f$x \leftarrow rhs^TB^{-1}\f$ in the columnwise case and
729 - \f$x \leftarrow B^{-1}rhs\f$ in the rowwise case.
730
731 Both can be seen uniformly as solving a linear system with the basis
732 matrix \p B and a right handside vector \p x aligned the same way as
733 the \em covectors of \p B.
734 */
735 void coSolve(VectorBase<R>& x, const VectorBase<R>& rhs)
736 {
737 if(rhs.dim() == 0)
738 {
739 x.clear();
740 return;
741 }
742
743 if(!factorized)
744 SPxBasisBase<R>::factorize();
745
746 factor->solveLeft(x, rhs);
747 }
748 /// Sparse version of coSolve
749 void coSolve(SSVectorBase<R>& x, const SVectorBase<R>& rhs)
750 {
751 if(rhs.size() == 0)
752 {
753 x.clear();
754 return;
755 }
756
757 if(!factorized)
758 SPxBasisBase<R>::factorize();
759
760 factor->solveLeft(x, rhs);
761 }
762 /// solves two systems in one call.
763 void coSolve(SSVectorBase<R>& x, VectorBase<R>& y, const SVectorBase<R>& rhsx,
764 SSVectorBase<R>& rhsy)
765 {
766 if(!factorized)
767 SPxBasisBase<R>::factorize();
768
769 factor->solveLeft(x, y, rhsx, rhsy);
770 }
771 /// Sparse version of solving two systems in one call.
772 void coSolve(SSVectorBase<R>& x, SSVectorBase<R>& y, const SVectorBase<R>& rhsx,
773 SSVectorBase<R>& rhsy)
774 {
775 if(!factorized)
776 SPxBasisBase<R>::factorize();
777
778 factor->solveLeft(x, y, rhsx, rhsy);
779 }
780 /// solves three systems in one call. May be improved by using just one pass through the basis.
781 void coSolve(SSVectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& z, const SVectorBase<R>& rhsx,
782 SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsz)
783 {
784 if(!factorized)
785 SPxBasisBase<R>::factorize();
786
787 factor->solveLeft(x, y, z, rhsx, rhsy, rhsz);
788 }
789 /// Sparse version of solving three systems in one call.
790 void coSolve(SSVectorBase<R>& x, SSVectorBase<R>& y, SSVectorBase<R>& z, const SVectorBase<R>& rhsx,
791 SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsz)
792 {
793 if(!factorized)
794 SPxBasisBase<R>::factorize();
795
796 factor->solveLeft(x, y, z, rhsx, rhsy, rhsz);
797 }
798 ///@}
799
800
801 //------------------------------------
802 /**@name Modification notification.
803 These methods must be called after the loaded LP has been modified.
804 */
805 ///@{
806 /// inform SPxBasisBase, that \p n new rows had been added.
807 void addedRows(int n);
808 /// inform SPxBasisBase that row \p i had been removed.
809 void removedRow(int i);
810 /// inform SPxBasisBase that rows in \p perm with negative entry were removed.
811 void removedRows(const int perm[]);
812 /// inform SPxBasisBase that \p n new columns had been added.
813 void addedCols(int n);
814 /// inform SPxBasisBase that column \p i had been removed.
815 void removedCol(int i);
816 /// inform SPxBasisBase that columns in \p perm with negative entry were removed.
817 void removedCols(const int perm[]);
818 /// inform SPxBasisBase that a row had been changed.
819 void changedRow(int);
820 /// inform SPxBasisBase that a column had been changed.
821 void changedCol(int);
822 /// inform SPxBasisBase that a matrix entry had been changed.
823 void changedElement(int, int);
824 ///@}
825
826
827 //--------------------------------
828 /**@name Miscellaneous */
829 ///@{
830 /// performs basis update.
831 /** Changes the \p i 'th vector of the basis with the vector associated to
832 \p id. This includes:
833 - updating the factorization, or recomputing it from scratch by
834 calling \ref soplex::SPxSolverBase<R>::factorize() "factorize()",
835 - resetting \ref soplex::SPxSolverBase<R>::lastEntered() "lastEntered()",
836 - resetting \ref soplex::SPxSolverBase<R>::lastIndex() "lastIndex()",
837 - resetting \ref soplex::SPxSolverBase<R>::lastLeft() "lastLeft()",
838 - resetting \ref soplex::SPxSolverBase<R>::lastUpdate() "lastUpdate()",
839 - resetting \ref soplex::SPxSolverBase<R>::iterations() "iterations()".
840
841 The basis descriptor is \em not \em modified, since #factor()
842 cannot know about how to set up the status of the involved variables
843 correctly.
844
845 A vector \p enterVec may be passed for a fast ETA update of the LU
846 factorization associated to the basis. It must be initialized with
847 the solution vector \f$x\f$ of the right linear system \f$Bx = b\f$
848 with the entering vector as right-hand side vector \f$b\f$, where \f$B\f$
849 denotes the basis matrix. This can be computed using method #solve().
850 When using FAST updates, a vector \p eta may be passed for
851 improved performance. It must be initialized by a call to
852 factor->solveRightUpdate() as described in SLinSolver. The
853 implementation hidden behind FAST updates depends on the
854 SLinSolver implementation class.
855 */
856 virtual void change(int i, SPxId& id,
857 const SVectorBase<R>* enterVec, const SSVectorBase<R>* eta = 0);
858
859 /** Load basis from \p in in MPS format. If \p rowNames and \p colNames
860 * are \c NULL, default names are used for the constraints and variables.
861 */
862 virtual bool readBasis(std::istream& in,
863 const NameSet* rowNames, const NameSet* colNames);
864
865 /** Write basis to \p os in MPS format. If \p rowNames and \p colNames are
866 * \c NULL, default names are used for the constraints and variables.
867 */
868 virtual void writeBasis(std::ostream& os,
869 const NameSet* rownames, const NameSet* colnames, const bool cpxFormat = false) const;
870
871 virtual void printMatrix() const;
872
873 /** Prints current basis matrix to a file using the MatrixMarket format:
874 * row col value
875 * The filename is basis/basis[number].mtx where number is a parameter.
876 */
877 void printMatrixMTX(int number);
878
879 /// checks if a Descriptor is valid for the current LP w.r.t. its bounds
880 virtual bool isDescValid(const Desc& ds);
881
882 /// sets up basis.
883 /** Loads a Descriptor to the basis and sets up the basis matrix and
884 all vectors accordingly. The Descriptor must have the same number of
885 rows and columns as the currently loaded LP.
886 */
887 virtual void loadDesc(const Desc&);
888
889 /// sets up linear solver to use.
890 /** If destroy is true, solver will be freed inside this object, e.g. in the destructor.
891 */
892 virtual void loadBasisSolver(SLinSolver<R>* solver, const bool destroy = false);
893
894 /// loads the LP \p lp to the basis.
895 /** This involves resetting all counters to 0 and setting up a regular
896 default basis consisting of slacks, artificial variables or bounds.
897 */
898 virtual void load(SPxSolverBase<R>* lp, bool initSlackBasis = true);
899
900 /// unloads the LP from the basis.
901 virtual void unLoad()
902 {
903 theLP = 0;
904 setStatus(NO_PROBLEM);
905 }
906
907 /// invalidates actual basis.
908 /** This method makes the basis matrix and vectors invalid. The basis will
909 be reinitialized if needed.
910 */
911 void invalidate();
912
913 /// Restores initial basis.
914 /** This method changes the basis to that present just after loading the LP
915 (see addedRows() and addedCols()). This may be necessary if a row or a
916 column is changed, since then the current basis may become singular.
917 */
918 void restoreInitialBasis();
919
920 /// output basis entries.
921 void dump();
922
923 /// consistency check.
924 bool isConsistent() const;
925
926 /// time spent in updates
927 Real getTotalUpdateTime() const
928 {
929 return theTime->time();
930 }
931 /// number of updates performed
932 int getTotalUpdateCount() const
933 {
934 return totalUpdateCount;
935 }
936
937 /// returns statistical information in form of a string.
938 std::string statistics() const
939 {
940 std::stringstream s;
941 s << factor->statistics()
942 #ifdef MEASUREUPDATETIME
943 << "Updates : " << std::setw(10) << getTotalUpdateCount() << std::endl
944 << " Time spent : " << std::setw(10) << getTotalUpdateTime() << std::endl
945 #endif
946 ;
947
948 return s.str();
949 }
950
951 void setOutstream(SPxOut& newOutstream)
952 {
953 spxout = &newOutstream;
954 }
955 ///@}
956
957 //--------------------------------------
958 /**@name Constructors / Destructors */
959 ///@{
960 /// default constructor.
961 SPxBasisBase<R>(Timer::TYPE ttype = Timer::USER_TIME);
962 /// copy constructor
963 SPxBasisBase<R>(const SPxBasisBase<R>& old);
964 /// assignment operator
965 SPxBasisBase<R>& operator=(const SPxBasisBase<R>& rhs);
966 /// destructor.
967 virtual ~SPxBasisBase<R>();
968 ///@}
969
970
971 protected:
972
973 //--------------------------------------
974 /**@name Protected helpers */
975 ///@{
976 /// loads \ref soplex::SPxBasisBase<R>::matrix "matrix" according to the SPxId%s stored in \ref soplex::SPxBasisBase<R>::theBaseId "theBaseId".
977 /** This method must be called whenever there is a chance, that the vector
978 pointers may have changed due to manipulations of the LP.
979 */
980 void loadMatrixVecs();
981
982 /// resizes internal arrays.
983 /** When a new LP is loaded, the basis matrix and vectors become invalid
984 and possibly also of the wrong dimension. Hence, after loading an
985 LP, #reDim() is called to reset all arrays etc. accoriding to the
986 dimensions of the loaded LP.
987 */
988 void reDim();
989
990 /// factorizes the basis matrix.
991 virtual void factorize();
992
993 /// sets descriptor representation according to loaded LP.
994 void setRep();
995 ///@}
996
997 };
998
999
1000 //
1001 // Auxiliary functions.
1002 //
1003
1004 /// Pretty-printing of basis status.
1005 template <class R>
1006 std::ostream& operator<<(std::ostream& os,
1007 const typename SPxBasisBase<R>::SPxStatus& status);
1008
1009
1010 /* For backwards compatibility */
1011 typedef SPxBasisBase<Real> SPxBasis;
1012
1013 } // namespace soplex
1014
1015 // General templated definitions
1016 #include "spxbasis.hpp"
1017 #include "spxdesc.hpp"
1018
1019 /* reset the SOPLEX_DEBUG flag to its original value */
1020 #undef SOPLEX_DEBUG
1021 #ifdef SOPLEX_DEBUG_SPXBASIS
1022 #define SOPLEX_DEBUG
1023 #undef SOPLEX_DEBUG_SPXBASIS
1024 #endif
1025
1026 #endif // _SPXBASIS_H_
1027