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 <iostream>
26   	
27   	#include "soplex/spxdefines.h"
28   	#include "soplex/spxbasis.h"
29   	#include "soplex/spxsolver.h"
30   	#include "soplex/spxout.h"
31   	
32   	namespace soplex
33   	{
34   	template <class R>
35   	void SPxBasisBase<R>::reDim()
36   	{
37   	
38   	   assert(theLP != 0);
39   	
40   	   SPxOut::debug(this, "DCHBAS01 SPxBasisBase<R>::reDim(): matrixIsSetup={}, factorized={}\n",
41   	                 matrixIsSetup, factorized);
42   	
43   	   thedesc.reSize(theLP->nRows(), theLP->nCols());
44   	
45   	   if(theLP->dim() != matrix.size())
46   	   {
47   	      SPX_MSG_INFO3((*this->spxout), (*this->spxout) <<
48   	                    "ICHBAS02 basis redimensioning invalidates factorization"
49   	                    << std::endl;)
50   	
51   	      matrix.reSize(theLP->dim());
52   	      theBaseId.reSize(theLP->dim());
53   	      matrixIsSetup = false;
54   	      factorized    = false;
55   	   }
56   	
57   	   SPxOut::debug(this, "DCHBAS03 SPxBasisBase<R>::reDim(): --> matrixIsSetup={}, factorized={}\n",
58   	                 matrixIsSetup, factorized);
59   	
60   	   assert(matrix.size()    >= theLP->dim());
61   	   assert(theBaseId.size() >= theLP->dim());
62   	}
63   	
64   	/* adapt basis and basis descriptor to added rows */
65   	template <class R>
66   	void SPxBasisBase<R>::addedRows(int n)
67   	{
68   	   assert(theLP != 0);
69   	
70   	   if(n > 0)
71   	   {
72   	      reDim();
73   	
74   	      if(theLP->rep() == SPxSolverBase<R>::COLUMN)
75   	      {
76   	         /* after adding rows in column representation, reDim() should set these bools to false. */
77   	         assert(!matrixIsSetup && !factorized);
78   	
79   	         for(int i = theLP->nRows() - n; i < theLP->nRows(); ++i)
80   	         {
81   	            thedesc.rowStatus(i) = dualRowStatus(i);
82   	            baseId(i) = theLP->SPxLPBase<R>::rId(i);
83   	         }
84   	      }
85   	      else
86   	      {
87   	         assert(theLP->rep() == SPxSolverBase<R>::ROW);
88   	
89   	         for(int i = theLP->nRows() - n; i < theLP->nRows(); ++i)
90   	            thedesc.rowStatus(i) = dualRowStatus(i);
91   	      }
92   	
93   	      /* If matrix was set up, load new basis vectors to the matrix.
94   	       * In the row case, the basis is not effected by adding rows. However,
95   	       * since @c matrix stores references to the rows in the LP (SPxLPBase<R>), a realloc
96   	       * in SPxLPBase<R> (e.g. due to space requirements) might invalidate these references.
97   	       * We therefore have to "reload" the matrix if it is set up. Note that reDim()
98   	       * leaves @c matrixIsSetup untouched if only row have been added, since the basis
99   	       * matrix already has the correct size. */
100  	      if(status() > NO_PROBLEM && matrixIsSetup)
101  	         loadMatrixVecs();
102  	
103  	      /* update basis status */
104  	      switch(status())
105  	      {
106  	      case PRIMAL:
107  	      case UNBOUNDED:
108  	         setStatus(REGULAR);
109  	         break;
110  	
111  	      case OPTIMAL:
112  	      case INFEASIBLE:
113  	         setStatus(DUAL);
114  	         break;
115  	
116  	      case NO_PROBLEM:
117  	      case SINGULAR:
118  	      case REGULAR:
119  	      case DUAL:
120  	         break;
121  	
122  	      default:
123  	         SPX_MSG_ERROR(std::cerr << "ECHBAS04 Unknown basis status!" << std::endl;)
124  	         throw SPxInternalCodeException("XCHBAS01 This should never happen.");
125  	      }
126  	   }
127  	}
128  	
129  	template <class R>
130  	void SPxBasisBase<R>::removedRow(int i)
131  	{
132  	
133  	   assert(status() >  NO_PROBLEM);
134  	   assert(theLP    != 0);
135  	
136  	   if(theLP->rep() == SPxSolverBase<R>::ROW)
137  	   {
138  	      if(theLP->isBasic(thedesc.rowStatus(i)))
139  	      {
140  	         setStatus(NO_PROBLEM);
141  	         factorized = false;
142  	
143  	         SPxOut::debug(this, "DCHBAS05 Warning: deleting basic row!\n");
144  	      }
145  	   }
146  	   else
147  	   {
148  	      assert(theLP->rep() == SPxSolverBase<R>::COLUMN);
149  	      factorized = false;
150  	
151  	      if(!theLP->isBasic(thedesc.rowStatus(i)))
152  	      {
153  	         setStatus(NO_PROBLEM);
154  	         SPxOut::debug(this, "DCHBAS06 Warning: deleting nonbasic row!\n");
155  	      }
156  	      else if(status() > NO_PROBLEM && matrixIsSetup)
157  	      {
158  	         for(int j = theLP->dim(); j >= 0; --j)
159  	         {
160  	            SPxId id = baseId(j);
161  	
162  	            if(id.isSPxRowId() && !theLP->has(SPxRowId(id)))
163  	            {
164  	               baseId(j) = baseId(theLP->dim());
165  	
166  	               if(j < theLP->dim())
167  	                  matrix[j] = &theLP->vector(baseId(j));
168  	
169  	               break;
170  	            }
171  	         }
172  	      }
173  	   }
174  	
175  	   thedesc.rowStatus(i) = thedesc.rowStatus(theLP->nRows());
176  	   reDim();
177  	}
178  	
179  	template <class R>
180  	void SPxBasisBase<R>::removedRows(const int perm[])
181  	{
182  	   assert(status() > NO_PROBLEM);
183  	   assert(theLP != 0);
184  	
185  	   int i;
186  	   int n = thedesc.nRows();
187  	
188  	   if(theLP->rep() == SPxSolverBase<R>::ROW)
189  	   {
190  	      for(i = 0; i < n; ++i)
191  	      {
192  	         if(perm[i] != i)
193  	         {
194  	            if(perm[i] < 0)                // row got removed
195  	            {
196  	               if(theLP->isBasic(thedesc.rowStatus(i)))
197  	               {
198  	                  setStatus(NO_PROBLEM);
199  	                  factorized = matrixIsSetup = false;
200  	                  SPxOut::debug(this, "DCHBAS07 Warning: deleting basic row!\n");
201  	               }
202  	            }
203  	            else                            // row was moved
204  	               thedesc.rowStatus(perm[i]) = thedesc.rowStatus(i);
205  	         }
206  	      }
207  	   }
208  	   else
209  	   {
210  	      assert(theLP->rep() == SPxSolverBase<R>::COLUMN);
211  	
212  	      factorized    = false;
213  	      matrixIsSetup = false;
214  	
215  	      for(i = 0; i < n; ++i)
216  	      {
217  	         if(perm[i] != i)
218  	         {
219  	            if(perm[i] < 0)                // row got removed
220  	            {
221  	               if(!theLP->isBasic(thedesc.rowStatus(i)))
222  	                  setStatus(NO_PROBLEM);
223  	            }
224  	            else                            // row was moved
225  	               thedesc.rowStatus(perm[i]) = thedesc.rowStatus(i);
226  	         }
227  	      }
228  	   }
229  	
230  	   reDim();
231  	}
232  	
233  	template <class R>
234  	static typename SPxBasisBase<R>::Desc::Status
235  	primalColStatus(int i, const SPxLPBase<R>* theLP)
236  	{
237  	   assert(theLP != 0);
238  	
239  	   if(theLP->upper(i) < R(infinity))
240  	   {
241  	      if(theLP->lower(i) > R(-infinity))
242  	      {
243  	         if(theLP->lower(i) == theLP->SPxLPBase<R>::upper(i))
244  	            return SPxBasisBase<R>::Desc::P_FIXED;
245  	         /*
246  	           else
247  	           return (-theLP->lower(i) < theLP->upper(i))
248  	           ? SPxBasisBase<R>::Desc::P_ON_LOWER
249  	           : SPxBasisBase<R>::Desc::P_ON_UPPER;
250  	         */
251  	         else if(theLP->maxObj(i) == 0)
252  	            return (-theLP->lower(i) < theLP->upper(i))
253  	                   ? SPxBasisBase<R>::Desc::P_ON_LOWER
254  	                   : SPxBasisBase<R>::Desc::P_ON_UPPER;
255  	         else
256  	            return (theLP->maxObj(i) < 0)
257  	                   ? SPxBasisBase<R>::Desc::P_ON_LOWER
258  	                   : SPxBasisBase<R>::Desc::P_ON_UPPER;
259  	      }
260  	      else
261  	         return SPxBasisBase<R>::Desc::P_ON_UPPER;
262  	   }
263  	   else if(theLP->lower(i) > R(-infinity))
264  	      return SPxBasisBase<R>::Desc::P_ON_LOWER;
265  	   else
266  	      return SPxBasisBase<R>::Desc::P_FREE;
267  	}
268  	
269  	
270  	/* adapt basis and basis descriptor to added columns */
271  	template <class R>
272  	void SPxBasisBase<R>::addedCols(int n)
273  	{
274  	   assert(theLP != 0);
275  	
276  	   if(n > 0)
277  	   {
278  	      reDim();
279  	
280  	      if(theLP->rep() == SPxSolverBase<R>::ROW)
281  	      {
282  	         /* after adding columns in row representation, reDim() should set these bools to false. */
283  	         assert(!matrixIsSetup && !factorized);
284  	
285  	         for(int i = theLP->nCols() - n; i < theLP->nCols(); ++i)
286  	         {
287  	            thedesc.colStatus(i) = primalColStatus(i, theLP);
288  	            baseId(i) = theLP->SPxLPBase<R>::cId(i);
289  	         }
290  	      }
291  	      else
292  	      {
293  	         assert(theLP->rep() == SPxSolverBase<R>::COLUMN);
294  	
295  	         for(int i = theLP->nCols() - n; i < theLP->nCols(); ++i)
296  	            thedesc.colStatus(i) = primalColStatus(i, theLP);
297  	      }
298  	
299  	      /* If matrix was set up, load new basis vectors to the matrix
300  	       * In the column case, the basis is not effected by adding columns. However,
301  	       * since @c matrix stores references to the columns in the LP (SPxLPBase<R>), a realloc
302  	       * in SPxLPBase<R> (e.g. due to space requirements) might invalidate these references.
303  	       * We therefore have to "reload" the matrix if it is set up. Note that reDim()
304  	       * leaves @c matrixIsSetup untouched if only columns have been added, since the
305  	       * basis matrix already has the correct size. */
306  	      if(status() > NO_PROBLEM && matrixIsSetup)
307  	         loadMatrixVecs();
308  	
309  	      switch(status())
310  	      {
311  	      case DUAL:
312  	      case INFEASIBLE:
313  	         setStatus(REGULAR);
314  	         break;
315  	
316  	      case OPTIMAL:
317  	      case UNBOUNDED:
318  	         setStatus(PRIMAL);
319  	         break;
320  	
321  	      case NO_PROBLEM:
322  	      case SINGULAR:
323  	      case REGULAR:
324  	      case PRIMAL:
325  	         break;
326  	
327  	      default:
328  	         SPX_MSG_ERROR(std::cerr << "ECHBAS08 Unknown basis status!" << std::endl;)
329  	         throw SPxInternalCodeException("XCHBAS02 This should never happen.");
330  	      }
331  	   }
332  	}
333  	
334  	template <class R>
335  	void SPxBasisBase<R>::removedCol(int i)
336  	{
337  	   assert(status() > NO_PROBLEM);
338  	   assert(theLP != 0);
339  	
340  	   if(theLP->rep() == SPxSolverBase<R>::COLUMN)
341  	   {
342  	      if(theLP->isBasic(thedesc.colStatus(i)))
343  	         setStatus(NO_PROBLEM);
344  	   }
345  	   else
346  	   {
347  	      assert(theLP->rep() == SPxSolverBase<R>::ROW);
348  	      factorized = false;
349  	
350  	      if(!theLP->isBasic(thedesc.colStatus(i)))
351  	         setStatus(NO_PROBLEM);
352  	      else if(status() > NO_PROBLEM)
353  	      {
354  	         for(int j = theLP->dim(); j >= 0; --j)
355  	         {
356  	            SPxId id = baseId(j);
357  	
358  	            if(id.isSPxColId() && !theLP->has(SPxColId(id)))
359  	            {
360  	               baseId(j) = baseId(theLP->dim());
361  	
362  	               if(matrixIsSetup &&
363  	                     j < theLP->dim())
364  	                  matrix[j] = &theLP->vector(baseId(j));
365  	
366  	               break;
367  	            }
368  	         }
369  	      }
370  	   }
371  	
372  	   thedesc.colStatus(i) = thedesc.colStatus(theLP->nCols());
373  	   reDim();
374  	}
375  	
376  	template <class R>
377  	void SPxBasisBase<R>::removedCols(const int perm[])
378  	{
379  	   assert(status() > NO_PROBLEM);
380  	   assert(theLP != 0);
381  	
382  	   int i;
383  	   int n = thedesc.nCols();
384  	
385  	   if(theLP->rep() == SPxSolverBase<R>::COLUMN)
386  	   {
387  	      for(i = 0; i < n; ++i)
388  	      {
389  	         if(perm[i] < 0)            // column got removed
390  	         {
391  	            if(theLP->isBasic(thedesc.colStatus(i)))
392  	               setStatus(NO_PROBLEM);
393  	         }
394  	         else                        // column was potentially moved
395  	            thedesc.colStatus(perm[i]) = thedesc.colStatus(i);
396  	      }
397  	   }
398  	   else
399  	   {
400  	      assert(theLP->rep() == SPxSolverBase<R>::ROW);
401  	      factorized = matrixIsSetup = false;
402  	
403  	      for(i = 0; i < n; ++i)
404  	      {
405  	         if(perm[i] != i)
406  	         {
407  	            if(perm[i] < 0)                // column got removed
408  	            {
409  	               if(!theLP->isBasic(thedesc.colStatus(i)))
410  	                  setStatus(NO_PROBLEM);
411  	            }
412  	            else                            // column was moved
413  	               thedesc.colStatus(perm[i]) = thedesc.colStatus(i);
414  	         }
415  	      }
416  	   }
417  	
418  	   reDim();
419  	}
420  	
421  	
422  	/**
423  	 * mark the basis as not factorized
424  	 */
425  	template <class R>
426  	void SPxBasisBase<R>::invalidate()
427  	{
428  	   if(factorized || matrixIsSetup)
429  	   {
430  	      SPX_MSG_INFO3((*this->spxout),
431  	                    (*this->spxout) << "ICHBAS09 explicit invalidation of factorization" <<
432  	                    std::endl;)
433  	   }
434  	
435  	   factorized    = false;
436  	   matrixIsSetup = false;
437  	}
438  	
439  	/**
440  	 * Create the initial slack basis descriptor and set up the basis matrix accordingly.
441  	 * This code has been adapted from SPxBasisBase<R>::addedRows() and SPxBasisBase<R>::addedCols().
442  	 */
443  	template <class R>
444  	void SPxBasisBase<R>::restoreInitialBasis()
445  	{
446  	   assert(!factorized);
447  	
448  	   SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ICHBAS10 setup slack basis" << std::endl;)
449  	
450  	   if(theLP->rep() == SPxSolverBase<R>::COLUMN)
451  	   {
452  	      for(int i = 0; i < theLP->nRows(); ++i)
453  	      {
454  	         thedesc.rowStatus(i) = dualRowStatus(i);
455  	         baseId(i) = theLP->SPxLPBase<R>::rId(i);
456  	      }
457  	
458  	      for(int i = 0; i < theLP->nCols(); ++i)
459  	         thedesc.colStatus(i) = primalColStatus(i, theLP);
460  	   }
461  	   else
462  	   {
463  	      assert(theLP->rep() == SPxSolverBase<R>::ROW);
464  	
465  	      for(int i = 0; i < theLP->nRows(); ++i)
466  	         thedesc.rowStatus(i) = dualRowStatus(i);
467  	
468  	      for(int i = 0; i < theLP->nCols(); ++i)
469  	      {
470  	         thedesc.colStatus(i) = primalColStatus(i, theLP);
471  	         baseId(i) = theLP->SPxLPBase<R>::cId(i);
472  	      }
473  	   }
474  	
475  	   /* if matrix was set up, load new basis vectors to the matrix */
476  	   if(status() > NO_PROBLEM && matrixIsSetup)
477  	      loadMatrixVecs();
478  	
479  	   /* update basis status */
480  	   setStatus(REGULAR);
481  	}
482  	
483  	/**
484  	 * @todo Implement changedRow(), changedCol(), changedElement() in a more clever
485  	 * way. For instance, the basis won't be singular (but maybe infeasible) if the
486  	 * change doesn't affect the basis rows/columns.
487  	 *
488  	 * The following methods (changedRow(), changedCol(), changedElement()) radically
489  	 * change the current basis to the original (slack) basis also present after
490  	 * loading the LP. The reason is that through the changes, the current basis may
491  	 * become singular. Going back to the initial basis is quite inefficient, but
492  	 * correct.
493  	 */
494  	
495  	/**@todo is this correctly implemented?
496  	 */
497  	template <class R>
498  	void SPxBasisBase<R>::changedRow(int /*row*/)
499  	{
500  	   invalidate();
501  	   restoreInitialBasis();
502  	}
503  	
504  	/**@todo is this correctly implemented?
505  	 */
506  	template <class R>
507  	void SPxBasisBase<R>::changedCol(int /*col*/)
508  	{
509  	   invalidate();
510  	   restoreInitialBasis();
511  	}
512  	
513  	/**@todo is this correctly implemented?
514  	 */
515  	template <class R>
516  	void SPxBasisBase<R>::changedElement(int /*row*/, int /*col*/)
517  	{
518  	   invalidate();
519  	   restoreInitialBasis();
520  	}
521  	} // namespace soplex
522