Scippy

SoPlex

Sequential object-oriented simPlex

spxfastrt.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2017 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 #include <assert.h>
17 #include <stdio.h>
18 
19 #include "spxdefines.h"
20 #include "spxfastrt.h"
21 
22 
23 /*
24  Here comes our implementation of the Harris procedure improved by shifting
25  bounds. The basic idea is to allow a slight infeasibility within |delta| to
26  allow for more freedom when selecting the leaving variable. This freedom
27  may then be used for selecting numerically stable variables yielding great
28  improvements.
29 
30  The algorithms operates in two phases. In a first phase, the maximum value
31  |val| is determined, when infeasibility within |delta| is allowed. In the second
32  phase, between all variables with value < |val| the one is selected which
33  has the largest update vector component. However, this may not
34  always yield an improvement. In that case, we shift the variable towards
35  infeasibility and retry. This avoids cycling in the shifted LP.
36  */
37 
38 namespace soplex
39 {
40 
41 #define MINSTAB 1e-5
42 #define LOWSTAB 1e-10
43 #define TRIES 2
44 #define SHORT 1e-5
45 #define DELTA_SHIFT 1e-5
46 #define EPSILON 1e-10
47 
49 {
50  // epsilon = thesolver->epsilon();
51  epsilon = EPSILON;
52  /*
53  if(thesolver->basis().stability() < 1e-4)
54  epsilon *= 1e-4 / thesolver->basis().stability();
55  */
56 }
57 
59 {
60  /*
61  if((delta > 1.99 * DELTA_SHIFT && thesolver->theShift < 1e-4) ||
62  (delta > 1e-4 && thesolver->theShift > 1e-4))
63  */
64  // if(delta > 1.99 * DELTA_SHIFT)
65  if (fastDelta >= delta + DELTA_SHIFT)
66  {
68  if (fastDelta > 1e-4)
69  fastDelta -= 2 * DELTA_SHIFT;
70  }
71 
72  if (minStab < MINSTAB)
73  {
74  minStab /= 0.90;
75  if (minStab < 1e-6)
76  minStab /= 0.90;
77  }
78 }
79 
81 {
82  minStab *= 0.95;
83  fastDelta += 3 * DELTA_SHIFT;
84  // delta += 2 * (thesolver->theShift > delta) * DELTA_SHIFT;
85 }
86 
88 {
89  if (maxabs < 1000.0)
90  return minStab;
91  return maxabs*minStab / 1000.0;
92 }
93 
94 /* The code below implements phase 1 of the ratio test. It differs from the description in the
95  * Ph.D. thesis page 57 as follows: It uses \f$\delta_i = d_i - s_i - \delta\f$ if \f$d_i > s_i\f$.
96  *
97  * This change leads to the following behavior of the source code. Consider the first case (x >
98  * epsilon, u < infinity): If u - vec[i] <= 0, vec[i] violates the upper bound. In the Harris ratio
99  * test, we would compute (u - vec[i] + delta)/upd[i]. The code computes instead delta/upd[i].
100  */
102  Real& val, /* on return: maximum step length */
103  Real& maxabs, /* on return: maximum absolute value in upd vector */
104  UpdateVector& update,
105  const Vector& lowBound,
106  const Vector& upBound,
107  int start,
108  int incr) const
109 {
110  int i, sel;
111  Real x, y, max;
112  Real u, l;
113  bool leaving = m_type == SPxSolver::LEAVE;
114 
115  Real mabs = maxabs;
116 
117  const Real* up = upBound.get_const_ptr();
118  const Real* low = lowBound.get_const_ptr();
119  const Real* vec = update.get_const_ptr();
120  const Real* upd = update.delta().values();
121  const int* idx = update.delta().indexMem();
122 
123  sel = -1;
124  max = val;
125 
126  if (update.delta().isSetup())
127  {
128  const int* last = idx + update.delta().size();
129  for (idx += start; idx < last; idx += incr)
130  {
131  i = *idx;
132 
133  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
134  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
135  continue;
136 
137  x = upd[i];
138 
139  if (x > epsilon)
140  {
141  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
142  mabs = (x > mabs) ? x : mabs;
143  u = up[i];
144  if (u < infinity)
145  {
146  y = u - vec[i];
147  if (y <= 0)
148  x = fastDelta / x;
149  else
150  x = (y + fastDelta) / x;
151 
152  if (x < max)
153  {
154  max = x;
155  sel = i;
156  }
157  }
158  }
159  else if (x < -epsilon)
160  {
161  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
162  mabs = (-x > mabs) ? -x : mabs;
163  l = low[i];
164  if (l > -infinity)
165  {
166  y = l - vec[i];
167  if ( y >= 0 )
168  x = - fastDelta / x;
169  else
170  x = ( y - fastDelta ) / x;
171 
172  if (x < max)
173  {
174  max = x;
175  sel = i;
176  }
177  }
178  }
179  }
180  }
181  else
182  {
183  /* In this case, the indices of the semi-sparse vector update.delta() are not set up and are filled below. */
184  int* l_idx = update.delta().altIndexMem();
185  Real* uval = update.delta().altValues();
186  const Real* uend = uval + update.dim();
187  assert( uval == upd );
188 
189  for (i = 0; uval < uend; ++uval, ++i)
190  {
191  x = *uval;
192  if (x != 0.0)
193  {
194  if (x >= -epsilon && x <= epsilon)
195  {
196  *uval = 0.0;
197  continue;
198  }
199  else
200  *l_idx++ = i;
201 
202  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
203  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
204  continue;
205 
206  if (x > epsilon)
207  {
208  mabs = (x > mabs) ? x : mabs;
209  u = up[i];
210  if (u < infinity)
211  {
212  y = u - vec[i];
213  if (y <= 0)
214  x = fastDelta / x;
215  else
216  x = (y + fastDelta) / x;
217 
218  if (x < max)
219  {
220  max = x;
221  sel = i;
222  }
223  }
224  }
225  else if (x < -epsilon)
226  {
227  mabs = (-x > mabs) ? -x : mabs;
228  l = low[i];
229  if (l > -infinity)
230  {
231  y = l - vec[i];
232  if ( y >= 0 )
233  x = - fastDelta / x;
234  else
235  x = ( y - fastDelta ) / x;
236 
237  if (x < max)
238  {
239  max = x;
240  sel = i;
241  }
242  }
243  }
244  }
245  }
246  update.delta().setSize(int(l_idx - update.delta().indexMem()));
247  update.delta().forceSetup();
248  }
249 
250  val = max;
251  maxabs = mabs;
252  return sel;
253 }
254 
255 /* See maxDelta() */
257  Real& val,
258  Real& maxabs,
259  UpdateVector& update,
260  const Vector& lowBound,
261  const Vector& upBound,
262  int start,
263  int incr) const
264 {
265  int i, sel;
266  Real x, y, max;
267  Real u, l;
268  bool leaving = m_type == SPxSolver::LEAVE;
269 
270  Real mabs = maxabs;
271 
272  const Real* up = upBound.get_const_ptr();
273  const Real* low = lowBound.get_const_ptr();
274  const Real* vec = update.get_const_ptr();
275  const Real* upd = update.delta().values();
276  const int* idx = update.delta().indexMem();
277 
278  sel = -1;
279  max = val;
280 
281  if (update.delta().isSetup())
282  {
283  const int* last = idx + update.delta().size();
284  for (idx += start; idx < last; idx += incr)
285  {
286  i = *idx;
287  x = upd[i];
288 
289  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
290  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
291  continue;
292 
293  if (x > epsilon)
294  {
295  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
296  mabs = (x > mabs) ? x : mabs;
297  l = low[i];
298  if (l > -infinity)
299  {
300  y = l - vec[i];
301  if ( y >= 0 )
302  x = - fastDelta / x;
303  else
304  x = ( y - fastDelta ) / x;
305 
306  if (x > max)
307  {
308  max = x;
309  sel = i;
310  }
311  }
312  }
313  else if (x < -epsilon)
314  {
315  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
316  mabs = (-x > mabs) ? -x : mabs;
317  u = up[i];
318  if (u < infinity)
319  {
320  y = u - vec[i];
321  if (y <= 0)
322  x = fastDelta / x;
323  else
324  x = (y + fastDelta) / x;
325 
326  if (x > max)
327  {
328  max = x;
329  sel = i;
330  }
331  }
332  }
333  }
334  }
335  else
336  {
337  /* In this case, the indices of the semi-sparse vector update.delta() are not set up and are filled below. */
338  int* l_idx = update.delta().altIndexMem();
339  Real* uval = update.delta().altValues();
340  const Real* uend = uval + update.dim();
341  assert( uval == upd );
342 
343  for (i = 0; uval < uend; ++uval, ++i)
344  {
345  x = *uval;
346 
347  if (x != 0.0)
348  {
349  if (x >= -epsilon && x <= epsilon)
350  {
351  *uval = 0.0;
352  continue;
353  }
354  else
355  *l_idx++ = i;
356 
357  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
358  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
359  continue;
360 
361  if (x > epsilon)
362  {
363  mabs = (x > mabs) ? x : mabs;
364  l = low[i];
365  if (l > -infinity)
366  {
367  y = l - vec[i];
368  if ( y >= 0 )
369  x = - fastDelta / x;
370  else
371  x = ( y - fastDelta ) / x;
372 
373  if (x > max)
374  {
375  max = x;
376  sel = i;
377  }
378  }
379  }
380  else if (x < -epsilon)
381  {
382  mabs = (-x > mabs) ? -x : mabs;
383  u = up[i];
384  if (u < infinity)
385  {
386  y = u - vec[i];
387  if (y <= 0)
388  x = fastDelta / x;
389  else
390  x = (y + fastDelta) / x;
391 
392  if (x > max)
393  {
394  max = x;
395  sel = i;
396  }
397  }
398  }
399  }
400  }
401  update.delta().setSize(int(l_idx - update.delta().indexMem()));
402  update.delta().forceSetup();
403  }
404 
405  val = max;
406  maxabs = mabs;
407 
408  return sel;
409 }
410 
412  Real& val,
413  Real& maxabs)
414 {
415  assert(m_type == SPxSolver::ENTER);
416  return maxDelta(val, maxabs,
417  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
418 }
419 
421  Real& val,
422  Real& maxabs)
423 {
424  assert(m_type == SPxSolver::ENTER);
425  return minDelta(val, maxabs,
426  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
427 }
428 
430  int& nr,
431  Real& max, /* on return: maximum step length */
432  Real& maxabs) /* on return: maximum absolute value in delta vector */
433 {
434  /* The following cause side effects on coPvec and pVec - both changes may be needed later in
435  maxSelect(). We can therefore not move the first function after the (indp >= 0) check. */
436  iscoid = true;
437  int indc = maxDelta(max, maxabs,
439  iscoid = false;
440  int indp = maxDelta(max, maxabs,
441  thesolver->pVec(), thesolver->lpBound(), thesolver->upBound(), 0, 1);
442 
443  if (indp >= 0)
444  {
445  nr = indp;
446  return thesolver->id(indp);
447  }
448  if (indc >= 0)
449  {
450  nr = indc;
451  return thesolver->coId(indc);
452  }
453  nr = -1;
454  return SPxId();
455 }
456 
458  int& nr,
459  Real& max,
460  Real& maxabs)
461 {
462  /* The following cause side effects on coPvec and pVec - both changes may be needed later in
463  minSelect(). We can therefore not move the first function after the (indp >= 0) check. */
464  iscoid = true;
465  const int indc = minDelta(max, maxabs,
467  iscoid = false;
468  const int indp = minDelta(max, maxabs,
469  thesolver->pVec(), thesolver->lpBound(), thesolver->upBound(), 0, 1);
470 
471  if (indp >= 0)
472  {
473  nr = indp;
474  return thesolver->id(indp);
475  }
476  if (indc >= 0)
477  {
478  nr = indc;
479  return thesolver->coId(indc);
480  }
481  nr = -1;
482  return SPxId();
483 }
484 
485 /* \p best returns the minimum update value such that the corresponding value of \p upd.delta() is
486  * at least \p stab and the update value is smaller than \p max. If no valid update value has been
487  * found \p bestDelta returns the slack to the bound corresponding to the index used for \p best. */
489  Real& val,
490  Real& stab,
491  Real& best,
492  Real& bestDelta,
493  Real max,
494  const UpdateVector& update,
495  const Vector& lowBound,
496  const Vector& upBound,
497  int start,
498  int incr) const
499 {
500  int i;
501  Real x, y;
502  bool leaving = m_type == SPxSolver::LEAVE;
503 
504  const Real* up = upBound.get_const_ptr();
505  const Real* low = lowBound.get_const_ptr();
506  const Real* vec = update.get_const_ptr();
507  const Real* upd = update.delta().values();
508  const int* idx = update.delta().indexMem();
509  const int* last = idx + update.delta().size();
510 
511  int nr = -1;
512  int bestNr = -1;
513 
514  for (idx += start; idx < last; idx += incr)
515  {
516  i = *idx;
517  x = upd[i];
518 
519  // in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables
520  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
521  continue;
522 
523  if (x > stab)
524  {
525  y = (low[i] - vec[i]) / x;
526 
527  if (y >= max)
528  {
529  val = y;
530  nr = i;
531  stab = x;
532  }
533  else if (y < best)
534  {
535  best = y;
536  bestNr = i;
537  }
538  }
539  else if (x < -stab)
540  {
541  y = (up[i] - vec[i]) / x;
542  if (y >= max)
543  {
544  val = y;
545  nr = i;
546  stab = -x;
547  }
548  else if (y < best)
549  {
550  best = y;
551  bestNr = i;
552  }
553  }
554  }
555 
556  if (nr < 0 && bestNr > 0)
557  {
558  if (upd[bestNr] < 0)
559  bestDelta = up[bestNr] - vec[bestNr];
560  else
561  bestDelta = vec[bestNr] - low[bestNr];
562  }
563  return nr;
564 }
565 
566 /* See minSelect() */
568  Real& val,
569  Real& stab,
570  Real& best,
571  Real& bestDelta,
572  Real max,
573  const UpdateVector& update,
574  const Vector& lowBound,
575  const Vector& upBound,
576  int start,
577  int incr) const
578 {
579  int i;
580  Real x, y;
581  bool leaving = m_type == SPxSolver::LEAVE;
582 
583  const Real* up = upBound.get_const_ptr();
584  const Real* low = lowBound.get_const_ptr();
585  const Real* vec = update.get_const_ptr();
586  const Real* upd = update.delta().values();
587  const int* idx = update.delta().indexMem();
588  const int* last = idx + update.delta().size();
589 
590  int nr = -1;
591  int bestNr = -1;
592 
593  for (idx += start; idx < last; idx += incr)
594  {
595  i = *idx;
596  x = upd[i];
597 
598  // in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables
599  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
600  continue;
601 
602  if (x > stab)
603  {
604  y = (up[i] - vec[i]) / x;
605  if (y <= max)
606  {
607  val = y;
608  nr = i;
609  stab = x;
610  }
611  else if (y > best)
612  {
613  best = y;
614  bestNr = i;
615  }
616  }
617  else if (x < -stab)
618  {
619  y = (low[i] - vec[i]) / x;
620  if (y <= max)
621  {
622  val = y;
623  nr = i;
624  stab = -x;
625  }
626  else if (y > best)
627  {
628  best = y;
629  bestNr = i;
630  }
631  }
632  }
633 
634  if (nr < 0 && bestNr > 0)
635  {
636  if (upd[bestNr] > 0)
637  bestDelta = up[bestNr] - vec[bestNr];
638  else
639  bestDelta = vec[bestNr] - low[bestNr];
640  }
641 
642  return nr;
643 }
644 
645 
647  Real& val,
648  Real& stab,
649  Real& bestDelta,
650  Real max)
651 {
652  Real best = -infinity;
653  bestDelta = 0.0;
654  assert(m_type == SPxSolver::ENTER);
655  return maxSelect(val, stab, best, bestDelta, max,
656  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
657 }
658 
660  int& nr,
661  Real& val,
662  Real& stab,
663  Real& bestDelta,
664  Real max
665 )
666 {
667  int indp, indc;
668  Real best = -infinity;
669  bestDelta = 0.0;
670  iscoid = true;
671  indc = maxSelect(val, stab, best, bestDelta, max,
673  iscoid = false;
674  indp = maxSelect(val, stab, best, bestDelta, max,
675  thesolver->pVec(), thesolver->lpBound(), thesolver->upBound(), 0, 1);
676 
677  if (indp >= 0)
678  {
679  nr = indp;
680  return thesolver->id(indp);
681  }
682  if (indc >= 0)
683  {
684  nr = indc;
685  return thesolver->coId(indc);
686  }
687  nr = -1;
688  return SPxId();
689 }
690 
692  Real& val,
693  Real& stab,
694  Real& bestDelta,
695  Real max)
696 {
697  Real best = infinity;
698  bestDelta = 0.0;
699  assert(m_type == SPxSolver::ENTER);
700  return minSelect(val, stab, best, bestDelta, max,
701  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
702 }
703 
705  int& nr,
706  Real& val,
707  Real& stab,
708  Real& bestDelta,
709  Real max)
710 {
711  Real best = infinity;
712  bestDelta = 0.0;
713  iscoid = true;
714  int indc = minSelect(val, stab, best, bestDelta, max,
716  iscoid = false;
717  int indp = minSelect(val, stab, best, bestDelta, max,
718  thesolver->pVec(), thesolver->lpBound(), thesolver->upBound(), 0, 1);
719 
720  if (indp >= 0)
721  {
722  nr = indp;
723  return thesolver->id(indp);
724  }
725  if (indc >= 0)
726  {
727  nr = indc;
728  return thesolver->coId(indc);
729  }
730  nr = -1;
731  return SPxId();
732 }
733 
734 
735 bool SPxFastRT::maxShortLeave(Real& sel, int leave, Real maxabs)
736 {
737  assert(leave >= 0);
738  assert(maxabs >= 0);
739 
740  sel = thesolver->fVec().delta()[leave];
741 
742  if (sel > maxabs*SHORT)
743  {
744  sel = (thesolver->ubBound()[leave] - thesolver->fVec()[leave]) / sel;
745  return true;
746  }
747 
748  if (sel < -maxabs*SHORT)
749  {
750  sel = (thesolver->lbBound()[leave] - thesolver->fVec()[leave]) / sel;
751  return true;
752  }
753 
754  return false;
755 }
756 
757 bool SPxFastRT::minShortLeave(Real& sel, int leave, Real maxabs)
758 {
759  assert(leave >= 0);
760  assert(maxabs >= 0);
761 
762  sel = thesolver->fVec().delta()[leave];
763 
764  if (sel > maxabs*SHORT)
765  {
766  sel = (thesolver->lbBound()[leave] - thesolver->fVec()[leave]) / sel;
767  return true;
768  }
769 
770  if ( sel < -maxabs*SHORT)
771  {
772  sel = (thesolver->ubBound()[leave] - thesolver->fVec()[leave]) / sel;
773  return true;
774  }
775 
776  return false;
777 }
778 
779 bool SPxFastRT::maxReLeave(Real& sel, int leave, Real maxabs, bool polish)
780 {
781  UpdateVector& vec = thesolver->fVec();
782  Vector& low = thesolver->lbBound();
783  Vector& up = thesolver->ubBound();
784 
785  if (leave < 0)
786  return true;
787 
788  if (up[leave] > low[leave])
789  {
790  Real x = vec.delta()[leave];
791 
792  if (sel < -fastDelta / maxabs)
793  {
794  sel = 0.0;
795  if( !polish && thesolver->dualStatus(thesolver->baseId(leave)) != SPxBasis::Desc::D_ON_BOTH )
796  {
797  if (x < 0.0)
798  thesolver->shiftLBbound(leave, vec[leave]);
799  else
800  thesolver->shiftUBbound(leave, vec[leave]);
801  }
802  }
803  }
804  else
805  {
806  sel = 0.0;
807  if( !polish )
808  {
809  thesolver->shiftLBbound(leave, vec[leave]);
810  thesolver->shiftUBbound(leave, vec[leave]);
811  }
812  }
813 
814  return false;
815 }
816 
817 bool SPxFastRT::minReLeave(Real& sel, int leave, Real maxabs, bool polish)
818 {
819  UpdateVector& vec = thesolver->fVec();
820  Vector& low = thesolver->lbBound();
821  Vector& up = thesolver->ubBound();
822 
823  if (leave < 0)
824  return true;
825 
826  if (up[leave] > low[leave])
827  {
828  Real x = vec.delta()[leave];
829 
830  if (sel > fastDelta / maxabs)
831  {
832  sel = 0.0;
833  if( !polish && thesolver->dualStatus(thesolver->baseId(leave)) != SPxBasis::Desc::D_ON_BOTH )
834  {
835  if (x > 0.0)
836  thesolver->shiftLBbound(leave, vec[leave]);
837  else
838  thesolver->shiftUBbound(leave, vec[leave]);
839  }
840  }
841  }
842  else
843  {
844  sel = 0.0;
845  if( !polish )
846  {
847  thesolver->shiftLBbound(leave, vec[leave]);
848  thesolver->shiftUBbound(leave, vec[leave]);
849  }
850  }
851 
852  return false;
853 }
854 
855 int SPxFastRT::selectLeave(Real& val, Real, bool polish)
856 {
857  Real maxabs, max, sel;
858  int leave = -1;
859  int cnt = 0;
860 
861  assert( m_type == SPxSolver::ENTER );
862 
863  // force instable pivot iff true (see explanation in enter.cpp and spxsolve.cpp)
864  bool instable = solver()->instableEnter;
865  Real lowstab = LOWSTAB;
866  assert(!instable || solver()->instableEnterId.isValid());
867 
868  resetTols();
869 
870  if (val > epsilon)
871  {
872  do
873  {
874  // phase 1:
875  max = val;
876  maxabs = 0.0;
877  leave = maxDelta(max, maxabs);
878 
879  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) ||
881 
882  if (max == val)
883  return -1;
884 
885  if (!maxShortLeave(sel, leave, maxabs))
886  {
887  // phase 2:
888  Real stab, bestDelta;
889 
890  stab = 100.0 * minStability(maxabs);
891 
892  // force instable pivot iff instable is true (see explanation in enter.cpp and spxsolve.cpp)
893  if (instable)
894  leave = maxSelect(sel, lowstab, bestDelta, max);
895  else
896  leave = maxSelect(sel, stab, bestDelta, max);
897 
898  if (bestDelta < DELTA_SHIFT*TRIES)
899  cnt++;
900  else
901  cnt += TRIES;
902  }
903  if (!maxReLeave(sel, leave, maxabs, polish))
904  break;
905  relax();
906  }
907  while (cnt < TRIES);
908  }
909  else if (val < -epsilon)
910  {
911  do
912  {
913  max = val;
914  maxabs = 0;
915  leave = minDelta(max, maxabs);
916 
917  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) ||
919 
920  if (max == val)
921  return -1;
922 
923  if (!minShortLeave(sel, leave, maxabs))
924  {
925  // phase 2:
926  Real stab, bestDelta;
927 
928  stab = 100.0 * minStability(maxabs);
929 
930  // force instable pivot iff instable is true (see explanation in enter.cpp and spxsolve.cpp)
931  if (instable)
932  leave = minSelect(sel, lowstab, bestDelta, max);
933  else
934  leave = minSelect(sel, stab, bestDelta, max);
935 
936  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) || thesolver->desc().colStatus(thesolver->number(SPxColId(thesolver->baseId(leave)))) != SPxBasis::Desc::P_FIXED);
937 
938  if (bestDelta < DELTA_SHIFT*TRIES)
939  cnt++;
940  else
941  cnt += TRIES;
942  }
943  if (!minReLeave(sel, leave, maxabs, polish))
944  break;
945  relax();
946  }
947  while (cnt < TRIES);
948  }
949  else
950  return -1;
951 
952  MSG_DEBUG(
953  if (leave >= 0)
954  std::cout
955  << "DFSTRT01 "
956  << thesolver->basis().iteration() << "("
957  << std::setprecision(6) << thesolver->value() << ","
958  << std::setprecision(2) << thesolver->basis().stability() << "):"
959  << leave << "\t"
960  << std::setprecision(4) << sel << " "
961  << std::setprecision(4) << thesolver->fVec().delta()[leave] << " "
962  << std::setprecision(6) << maxabs
963  << std::endl;
964  else
965  std::cout << "DFSTRT02 " << thesolver->basis().iteration()
966  << ": skipping instable pivot" << std::endl;
967  )
968 
969  if( polish && leave >= 0 )
970  {
971  assert( thesolver->rep() == SPxSolver::COLUMN );
972  SPxId leaveId = thesolver->baseId(leave);
973  // decide whether the chosen leave index contributes to the polishing objective
975  {
976  // only allow (integer) variables to leave the basis
977  if( leaveId.isSPxRowId() )
978  return -1;
979  else if( thesolver->integerVariables.size() == thesolver->nCols() )
980  {
981  if( leaveId.isSPxColId() && thesolver->integerVariables[thesolver->number(leaveId)] == 0 )
982  return -1;
983  }
984  }
986  {
987  // only allow slacks and continuous variables to leave the basis
989  {
990  if( thesolver->baseId(leave).isSPxColId() && thesolver->integerVariables[thesolver->number(leaveId)] == 1 )
991  return -1;
992  }
993  else if( thesolver->baseId(leave).isSPxColId() )
994  return -1;
995  }
996  }
997 
998  if (leave >= 0 || minStab > 2*solver()->epsilon())
999  {
1000  val = sel;
1001  if (leave >= 0)
1002  tighten();
1003  }
1004 
1005  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) || thesolver->desc().colStatus(thesolver->number(SPxColId(thesolver->baseId(leave)))) != SPxBasis::Desc::P_FIXED);
1006 
1007  return leave;
1008 }
1009 
1010 
1012  Real& sel,
1013  Real maxabs,
1014  const SPxId& id,
1015  int nr)
1016 {
1017  Real x, d;
1018  Vector* up;
1019  Vector* low;
1020 
1021  UpdateVector& pvec = thesolver->pVec();
1022  SSVector& pupd = thesolver->pVec().delta();
1023  Vector& upb = thesolver->upBound();
1024  Vector& lpb = thesolver->lpBound();
1025  UpdateVector& cvec = thesolver->coPvec();
1026  SSVector& cupd = thesolver->coPvec().delta();
1027  Vector& ucb = thesolver->ucBound();
1028  Vector& lcb = thesolver->lcBound();
1029 
1030  if (thesolver->isCoId(id))
1031  {
1032  if (thesolver->isCoBasic(nr))
1033  {
1034  cupd.clearIdx(nr);
1035  return true;
1036  }
1037 
1038  x = cvec[nr];
1039  d = cupd[nr];
1040  up = &ucb;
1041  low = &lcb;
1042 
1043  if (d < 0.0)
1044  sel = (lcb[nr] - cvec[nr]) / d;
1045  else
1046  sel = (ucb[nr] - cvec[nr]) / d;
1047  }
1048  else if (thesolver->isId(id))
1049  {
1050  pvec[nr] = thesolver->vector(nr) * cvec;
1051  if (thesolver->isBasic(nr))
1052  {
1053  pupd.clearIdx(nr);
1054  return true;
1055  }
1056 
1057  x = pvec[nr];
1058  d = pupd[nr];
1059  up = &upb;
1060  low = &lpb;
1061 
1062  if (d < 0.0)
1063  sel = (lpb[nr] - pvec[nr]) / d;
1064  else
1065  sel = (upb[nr] - pvec[nr]) / d;
1066  }
1067  else
1068  return true;
1069 
1070  if ((*up)[nr] != (*low)[nr])
1071  {
1072  if (sel < -fastDelta / maxabs)
1073  {
1074  if (d > 0.0)
1075  {
1076  thesolver->theShift -= (*up)[nr];
1077  sel = 0.0;
1078  (*up)[nr] = x + sel * d;
1079  thesolver->theShift += (*up)[nr];
1080  }
1081  else
1082  {
1083  thesolver->theShift += (*low)[nr];
1084  sel = 0.0;
1085  (*low)[nr] = x + sel * d;
1086  thesolver->theShift -= (*low)[nr];
1087  }
1088  }
1089  }
1090  else
1091  {
1092  sel = 0.0;
1093  if (x > (*up)[nr])
1094  thesolver->theShift += x - (*up)[nr];
1095  else
1096  thesolver->theShift += (*low)[nr] - x;
1097  (*up)[nr] = (*low)[nr] = x;
1098  }
1099 
1100  return false;
1101 }
1102 
1104  Real& sel,
1105  Real maxabs,
1106  const SPxId& id,
1107  int nr)
1108 {
1109  Real x, d;
1110  Vector* up;
1111  Vector* low;
1112 
1113  UpdateVector& pvec = thesolver->pVec();
1114  SSVector& pupd = thesolver->pVec().delta();
1115  Vector& upb = thesolver->upBound();
1116  Vector& lpb = thesolver->lpBound();
1117  UpdateVector& cvec = thesolver->coPvec();
1118  SSVector& cupd = thesolver->coPvec().delta();
1119  Vector& ucb = thesolver->ucBound();
1120  Vector& lcb = thesolver->lcBound();
1121 
1122  if (thesolver->isCoId(id))
1123  {
1124  if (thesolver->isCoBasic(nr))
1125  {
1126  cupd.clearIdx(nr);
1127  return true;
1128  }
1129  x = cvec[nr];
1130  d = cupd[nr];
1131  up = &ucb;
1132  low = &lcb;
1133  if (d > 0.0)
1134  sel = (thesolver->lcBound()[nr] - cvec[nr]) / d;
1135  else
1136  sel = (thesolver->ucBound()[nr] - cvec[nr]) / d;
1137  }
1138 
1139  else if (thesolver->isId(id))
1140  {
1141  pvec[nr] = thesolver->vector(nr) * cvec;
1142  if (thesolver->isBasic(nr))
1143  {
1144  pupd.clearIdx(nr);
1145  return true;
1146  }
1147  x = pvec[nr];
1148  d = pupd[nr];
1149  up = &upb;
1150  low = &lpb;
1151  if (d > 0.0)
1152  sel = (thesolver->lpBound()[nr] - pvec[nr]) / d;
1153  else
1154  sel = (thesolver->upBound()[nr] - pvec[nr]) / d;
1155  }
1156 
1157  else
1158  return true;
1159 
1160  if ((*up)[nr] != (*low)[nr])
1161  {
1162  if (sel > fastDelta / maxabs)
1163  {
1164  if (d < 0.0)
1165  {
1166  thesolver->theShift -= (*up)[nr];
1167  sel = 0.0;
1168  (*up)[nr] = x + sel * d;
1169  thesolver->theShift += (*up)[nr];
1170  }
1171  else
1172  {
1173  thesolver->theShift += (*low)[nr];
1174  sel = 0.0;
1175  (*low)[nr] = x + sel * d;
1176  thesolver->theShift -= (*low)[nr];
1177  }
1178  }
1179  }
1180  else
1181  {
1182  sel = 0.0;
1183  if (x > (*up)[nr])
1184  thesolver->theShift += x - (*up)[nr];
1185  else
1186  thesolver->theShift += (*low)[nr] - x;
1187  (*up)[nr] = (*low)[nr] = x;
1188  }
1189 
1190  return false;
1191 }
1192 
1194  const SPxId& enterId,
1195  int nr,
1196  Real max,
1197  Real maxabs) const
1198 {
1199  if (thesolver->isCoId(enterId))
1200  {
1201  if (max != 0.0)
1202  {
1203  Real x = thesolver->coPvec().delta()[nr];
1204  if (x < maxabs * SHORT && -x < maxabs * SHORT)
1205  return false;
1206  }
1207  return true;
1208  }
1209  else if (thesolver->isId(enterId))
1210  {
1211  if (max != 0.0)
1212  {
1213  Real x = thesolver->pVec().delta()[nr];
1214  if (x < maxabs * SHORT && -x < maxabs * SHORT)
1215  return false;
1216  }
1217  return true;
1218  }
1219 
1220  return false;
1221 }
1222 
1223 SPxId SPxFastRT::selectEnter(Real& val, int, bool polish)
1224 {
1225  SPxId enterId;
1226  Real max, sel;
1227  Real maxabs = 0.0;
1228  int nr;
1229  int cnt = 0;
1230 
1231  assert( m_type == SPxSolver::LEAVE );
1232 
1233  // force instable pivot iff true (see explanation in leave.cpp and spxsolve.cpp)
1234  bool instable = solver()->instableLeave;
1235  Real lowstab = LOWSTAB;
1236  assert(!instable || solver()->instableLeaveNum >= 0);
1237 
1238  resetTols();
1239  sel = 0.0;
1240 
1241  if (val > epsilon)
1242  {
1243  do
1244  {
1245  maxabs = 0.0;
1246  max = val;
1247 
1248  enterId = maxDelta(nr, max, maxabs);
1249  if (!enterId.isValid())
1250  return enterId;
1251 
1252  assert(max >= 0.0);
1253  assert(!enterId.isValid() || !solver()->isBasic(enterId));
1254 
1255  if (!shortEnter(enterId, nr, max, maxabs))
1256  {
1257  Real bestDelta, stab;
1258 
1259  stab = minStability(maxabs);
1260 
1261  // force instable pivot iff instable is true (see explanation in leave.cpp and spxsolve.cpp)
1262  if (instable)
1263  {
1264  enterId = maxSelect(nr, sel, lowstab, bestDelta, max);
1265  }
1266  else
1267  {
1268  enterId = maxSelect(nr, sel, stab, bestDelta, max);
1269  }
1270  if (bestDelta < DELTA_SHIFT*TRIES)
1271  cnt++;
1272  else
1273  cnt += TRIES;
1274  }
1275  if (!maxReEnter(sel, maxabs, enterId, nr))
1276  break;
1277  relax();
1278  }
1279  while (cnt < TRIES);
1280  }
1281  else if (val < -epsilon)
1282  {
1283  do
1284  {
1285  maxabs = 0.0;
1286  max = val;
1287  enterId = minDelta(nr, max, maxabs);
1288  if (!enterId.isValid())
1289  return enterId;
1290 
1291  assert(max <= 0.0);
1292  assert(!enterId.isValid() || !solver()->isBasic(enterId));
1293 
1294  if (!shortEnter(enterId, nr, max, maxabs))
1295  {
1296  Real bestDelta, stab;
1297 
1298  stab = minStability(maxabs);
1299 
1300  // force instable pivot iff instable is true (see explanation in leave.cpp and spxsolve.cpp)
1301  if (instable)
1302  {
1303  enterId = minSelect(nr, sel, lowstab, bestDelta, max);
1304  }
1305  else
1306  {
1307  enterId = minSelect(nr, sel, stab, bestDelta, max);
1308  }
1309  if (bestDelta < DELTA_SHIFT*TRIES)
1310  cnt++;
1311  else
1312  cnt += TRIES;
1313  }
1314  if (!minReEnter(sel, maxabs, enterId, nr))
1315  break;
1316  relax();
1317  }
1318  while (cnt < TRIES);
1319  }
1320 
1321  MSG_DEBUG(
1322  if (enterId.isValid())
1323  {
1324  assert(!enterId.isValid() || !solver()->isBasic(enterId));
1325 
1326  Real x;
1327  if (thesolver->isCoId(enterId))
1328  x = thesolver->coPvec().delta()[ thesolver->number(enterId) ];
1329  else
1330  x = thesolver->pVec().delta()[ thesolver->number(enterId) ];
1331  std::cout << "DFSTRT03 " << thesolver->basis().iteration() << ": "
1332  << sel << '\t' << x << " (" << maxabs << ")" << std::endl;
1333  }
1334  else
1335  std::cout << "DFSTRT04 " << thesolver->basis().iteration()
1336  << ": skipping instable pivot" << std::endl;
1337  )
1338 
1339  if( polish && enterId.isValid() )
1340  {
1341  assert( thesolver->rep() == SPxSolver::ROW );
1342  // decide whether the chosen entering index contributes to the polishing objective
1344  {
1345  // only allow (integer) variables to enter the basis
1346  if( enterId.isSPxRowId() )
1347  return SPxId();
1349  return SPxId();
1350  }
1352  {
1353  // only allow slacks and continuous variables to enter the basis
1355  {
1356  if( enterId.isSPxColId() && thesolver->integerVariables[thesolver->number(enterId)] == 1 )
1357  return SPxId();
1358  }
1359  else if( enterId.isSPxColId() )
1360  return SPxId();
1361  }
1362  }
1363 
1364 
1365  if (enterId.isValid() || minStab > 2*epsilon)
1366  {
1367  val = sel;
1368  if (enterId.isValid())
1369  tighten();
1370  }
1371 
1372  assert(!enterId.isValid() || !solver()->isBasic(enterId));
1373 
1374  return enterId;
1375 }
1376 
1378 {
1379  thesolver = spx;
1380  setType(spx->type());
1381 }
1382 
1384 {
1385  m_type = type;
1386 
1387  minStab = MINSTAB;
1388  fastDelta = delta;
1389 }
1390 } // namespace soplex
Fast shifting ratio test.
SPxId coId(int i) const
id of i &#39;th covector.
Definition: spxsolver.h:1089
bool maxReEnter(Real &sel, Real maxabs, const SPxId &id, int nr)
Definition: spxfastrt.cpp:1011
const Vector & ucBound() const
Definition: spxsolver.h:1389
int iteration() const
returns number of basis changes since last load().
Definition: spxbasis.h:545
bool isSetup() const
Returns setup status.
Definition: ssvectorbase.h:120
primal variable is fixed to both bounds
Definition: spxbasis.h:190
THREADLOCAL const Real infinity
Definition: spxdefines.cpp:26
Type
Algorithmic type.
Definition: spxsolver.h:124
bool maxShortLeave(Real &sel, int leave, Real maxabs)
Definition: spxfastrt.cpp:735
int minSelect(Real &val, Real &stab, Real &best, Real &bestDelta, Real max, const UpdateVector &upd, const Vector &low, const Vector &up, int start=0, int incr=1) const
selects stable index for minimizing ratio test.
Definition: spxfastrt.cpp:488
#define TRIES
Definition: spxfastrt.cpp:43
Real delta
allowed bound violation
minimize number of basic slack variables, i.e. more variables in between bounds
Definition: spxsolver.h:231
Real fastDelta
currently allowed infeasibility.
Definition: spxfastrt.h:52
#define SHORT
Definition: spxfastrt.cpp:44
int number(const SPxRowId &id) const
Returns the row number of the row with identifier id.
Definition: spxlpbase.h:522
virtual int selectLeave(Real &val, Real, bool polish=false)
Definition: spxfastrt.cpp:855
bool isId(const SPxId &p_id) const
Is p_id an SPxId ?
Definition: spxsolver.h:1107
bool minReLeave(Real &sel, int leave, Real maxabs, bool polish=false)
numerical stability tests.
Definition: spxfastrt.cpp:817
const Vector & lcBound() const
Definition: spxsolver.h:1410
void shiftLBbound(int i, Real to)
shift i &#39;th lbBound to to.
Definition: spxsolver.h:1555
Real epsilon
|value| < epsilon is considered 0.
Definition: spxfastrt.h:50
virtual void setType(SPxSolver::Type type)
Definition: spxfastrt.cpp:1383
Ids for LP columns.Class SPxColId provides DataKeys for the column indices of an SPxLP.
Definition: spxid.h:36
int maxDelta(Real &val, Real &maxabs, UpdateVector &update, const Vector &lowBound, const Vector &upBound, int start, int incr) const
Max phase 1 value.
Definition: spxfastrt.cpp:101
void relax()
relaxes stability requirements.
Definition: spxfastrt.cpp:80
Desc::Status dualStatus(const SPxColId &id) const
dual Status for the column variable with ID id of the loaded LP.
Definition: spxbasis.cpp:34
#define LOWSTAB
Definition: spxfastrt.cpp:42
virtual Real value()
current objective value.
Definition: spxsolver.cpp:887
rowwise representation.
Definition: spxsolver.h:107
int * altIndexMem()
Returns array indices.
Definition: ssvectorbase.h:311
maximize number of basic slack variables, i.e. more variables on bounds
Definition: spxsolver.h:230
const R * get_const_ptr() const
Conversion to C-style pointer.
Definition: vectorbase.h:453
const Vector & ubBound() const
upper bound for fVec.
Definition: spxsolver.h:1313
Entering Simplex.
Definition: spxsolver.h:134
Real minStability(Real maxabs)
Compute stability requirement.
Definition: spxfastrt.cpp:87
R * altValues()
Returns array values.
Definition: ssvectorbase.h:318
Generic Ids for LP rows or columns.Both SPxColIds and SPxRowIds may be treated uniformly as SPxIds: ...
Definition: spxid.h:85
UpdateVector & coPvec() const
copricing vector.
Definition: spxsolver.h:1370
void resetTols()
resets tolerances (epsilon).
Definition: spxfastrt.cpp:48
Leaving Simplex.
Definition: spxsolver.h:143
double Real
Definition: spxdefines.h:215
const R * values() const
Returns array values.
Definition: ssvectorbase.h:299
#define MSG_DEBUG(x)
Definition: spxdefines.h:129
const Vector & lbBound() const
lower bound for fVec.
Definition: spxsolver.h:1331
SPxId id(int i) const
id of i &#39;th vector.
Definition: spxsolver.h:1070
UpdateVector & fVec() const
feasibility vector.
Definition: spxsolver.h:1295
bool isSPxColId() const
is id a column id?
Definition: spxid.h:165
bool shortEnter(const SPxId &enterId, int nr, Real max, Real maxabs) const
Tests and returns whether a shortcut after phase 1 is feasible for the selected enter pivot...
Definition: spxfastrt.cpp:1193
bool isBasic(SPxBasis::Desc::Status stat) const
does stat describe a basic index ?
Definition: spxsolver.h:1235
Dense vector with semi-sparse vector for updatesIn many algorithms vectors are updated in every itera...
Definition: updatevector.h:53
void shiftUBbound(int i, Real to)
shift i &#39;th ubBound to to.
Definition: spxsolver.h:1548
Status & colStatus(int i)
Definition: spxbasis.h:254
SPxId & baseId(int i)
Definition: spxbasis.h:503
const int * indexMem() const
Returns array indices.
Definition: ssvectorbase.h:293
Real epsilon() const
values are considered to be 0.
Definition: spxsolver.h:758
Debugging, floating point type and parameter definitions.
void setSize(int n)
Sets number of nonzeros (thereby unSetup SSVectorBase).
Definition: ssvectorbase.h:580
void clearIdx(int i)
Clears element i.
Definition: ssvectorbase.h:253
bool isCoId(const SPxId &p_id) const
Is p_id a CoId.
Definition: spxsolver.h:1116
Sequential object-oriented SimPlex.SPxSolver is an LP solver class using the revised Simplex algorith...
Definition: spxsolver.h:84
SolutionPolish polishObj
objective of solution polishing
Definition: spxsolver.h:245
UpdateVector & pVec() const
pricing vector.
Definition: spxsolver.h:1450
SPxSolver * thesolver
the solver
const Vector & lpBound() const
Definition: spxsolver.h:1476
int dim() const
Dimension of vector.
Definition: vectorbase.h:215
int size() const
Returns the number of nonzeros.
Definition: ssvectorbase.h:199
Everything should be within this namespace.
void tighten()
tightens stability requirements.
Definition: spxfastrt.cpp:58
bool iscoid
flag used in methods minSelect/maxSelect to retrieve correct basis status
Definition: spxfastrt.h:54
virtual SPxSolver * solver() const
returns loaded LP solver.
int maxSelect(Real &val, Real &stab, Real &best, Real &bestDelta, Real max, const UpdateVector &upd, const Vector &low, const Vector &up, int start=0, int incr=1) const
selects stable index for maximizing ratio test.
Definition: spxfastrt.cpp:567
bool minShortLeave(Real &sel, int leave, Real maxabs)
tests for stop after phase 1.
Definition: spxfastrt.cpp:757
int minDelta(Real &val, Real &maxabs, UpdateVector &update, const Vector &lowBound, const Vector &upBound, int start, int incr) const
Min phase 1 value.
Definition: spxfastrt.cpp:256
bool maxReLeave(Real &sel, int leave, Real maxabs, bool polish=false)
Definition: spxfastrt.cpp:779
virtual void load(SPxSolver *solver)
Definition: spxfastrt.cpp:1377
const Vector & upBound() const
Definition: spxsolver.h:1455
virtual SPxId selectEnter(Real &val, int, bool polish=false)
Definition: spxfastrt.cpp:1223
void forceSetup()
Forces setup status.
Definition: ssvectorbase.h:162
bool isSPxRowId() const
is id a row id?
Definition: spxid.h:160
bool isValid() const
returns TRUE iff the id is a valid column or row identifier.
Definition: spxid.h:150
int size() const
return nr. of elements.
Definition: dataarray.h:211
Type type() const
return current Type.
Definition: spxsolver.h:475
Real theShift
sum of all shifts applied to any bound.
Definition: spxsolver.h:266
int nCols() const
Returns number of columns in LP.
Definition: spxlpbase.h:157
bool isCoBasic(int i) const
is the i &#39;th covector basic ?
Definition: spxsolver.h:1280
DataArray< int > integerVariables
supplementary variable information, 0: continous variable, 1: integer variable
Definition: spxsolver.h:429
Real minStab
parameter for computing minimum stability requirement
Definition: spxfastrt.h:48
dual variable has two bounds
Definition: spxbasis.h:194
const SVector & vector(int i) const
i &#39;th vector.
Definition: spxsolver.h:1129
SSVector & delta()
update vector , writeable
Definition: updatevector.h:122
#define MINSTAB
Definition: spxfastrt.cpp:41
const SPxBasis & basis() const
Return current basis.
Definition: spxsolver.h:1727
bool minReEnter(Real &sel, Real maxabs, const SPxId &id, int nr)
numerical stability check.
Definition: spxfastrt.cpp:1103
#define EPSILON
Definition: spxfastrt.cpp:46
Real stability() const
returns the stability of the basis matrix.
Definition: spxbasis.h:626
SPxSolver::Type m_type
internal storage of type
const Desc & desc() const
Definition: spxbasis.h:463
Representation rep() const
return the current basis representation.
Definition: spxsolver.h:469
columnwise representation.
Definition: spxsolver.h:108
#define DELTA_SHIFT
Definition: spxfastrt.cpp:45