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