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