SoPlex Doxygen Documentation
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-2012 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 //#define DEBUGGING 1
17 
18 #include <assert.h>
19 #include <stdio.h>
20 
21 #include "spxdefines.h"
22 #include "spxfastrt.h"
23 
24 
25 /*
26  Here comes our implementation of the Harris procedure improved by shifting
27  bounds. The basic idea is to allow a slight infeasibility within |delta| to
28  allow for more freedom when selecting the leaving variable. This freedom
29  may then be used for selecting numerically stable variables yielding great
30  improvements.
31 
32  The algorithms operates in two phases. In a first phase, the maximum value
33  |val| is determined, when infeasibility within |delta| is allowed. In the second
34  phase, between all variables with value < |val| the one is selected which
35  has the largest update vector component. However, this may not
36  always yield an improvement. In that case, we shift the variable towards
37  infeasibility and retry. This avoids cycling in the shifted LP.
38  */
39 
40 namespace soplex
41 {
42 
43 #define MINSTAB 1e-5
44 #define LOWSTAB 1e-10
45 #define TRIES 2
46 #define SHORT 1e-5
47 #define DELTA_SHIFT 1e-5
48 #define EPSILON 1e-10
49 
50 void SPxFastRT::resetTols()
51 {
52  // epsilon = thesolver->epsilon();
53  epsilon = EPSILON;
54  /*
55  if(thesolver->basis().stability() < 1e-4)
56  epsilon *= 1e-4 / thesolver->basis().stability();
57  */
58 }
59 
60 void SPxFastRT::tighten()
61 {
62  /*
63  if((delta > 1.99 * DELTA_SHIFT && thesolver->theShift < 1e-4) ||
64  (delta > 1e-4 && thesolver->theShift > 1e-4))
65  */
66  // if(delta > 1.99 * DELTA_SHIFT)
67  if (fastDelta >= delta + DELTA_SHIFT)
68  {
69  fastDelta -= DELTA_SHIFT;
70  if (fastDelta > 1e-4)
71  fastDelta -= 2 * DELTA_SHIFT;
72  }
73 
74  if (minStab < MINSTAB)
75  {
76  minStab /= 0.90;
77  if (minStab < 1e-6)
78  minStab /= 0.90;
79  }
80 }
81 
82 void SPxFastRT::relax()
83 {
84  minStab *= 0.95;
85  fastDelta += 3 * DELTA_SHIFT;
86  // delta += 2 * (thesolver->theShift > delta) * DELTA_SHIFT;
87 }
88 
89 Real SPxFastRT::minStability(Real maxabs)
90 {
91  if (maxabs < 1000.0)
92  return minStab;
93  return maxabs*minStab / 1000.0;
94 }
95 
96 /* The code below implements phase 1 of the ratio test. It differs from the description in the
97  * 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$.
98  *
99  * This change leads to the following behavior of the source code. Consider the first case (x >
100  * epsilon, u < infinity): If u - vec[i] <= 0, vec[i] violates the upper bound. In the Harris ratio
101  * test, we would compute (u - vec[i] + delta)/upd[i]. The code computes instead delta/upd[i].
102  */
103 int SPxFastRT::maxDelta(
104  Real& val, /* on return: maximum step length */
105  Real& maxabs, /* on return: maximum absolute value in upd vector */
106  UpdateVector& update,
107  const Vector& lowBound,
108  const Vector& upBound,
109  int start,
110  int incr) const
111 {
112  int i, sel;
113  Real x, y, max;
114  Real u, l;
115  bool leaving = m_type == SPxSolver::LEAVE;
116 
117  Real mabs = maxabs;
118 
119  const Real* up = upBound.get_const_ptr();
120  const Real* low = lowBound.get_const_ptr();
121  const Real* vec = update.get_const_ptr();
122  const Real* upd = update.delta().values();
123  const int* idx = update.delta().indexMem();
124 
125  sel = -1;
126  max = val;
127 
128  if (update.delta().isSetup())
129  {
130  const int* last = idx + update.delta().size();
131  for (idx += start; idx < last; idx += incr)
132  {
133  i = *idx;
134  x = upd[i];
135 
136  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
137  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
138  continue;
139 
140  if (x > epsilon)
141  {
142  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
143  mabs = (x > mabs) ? x : mabs;
144  u = up[i];
145  if (u < infinity)
146  {
147  y = u - vec[i];
148  if (y <= 0)
149  x = fastDelta / x;
150  else
151  x = (y + fastDelta) / x;
152 
153  if (x < max)
154  {
155  max = x;
156  sel = i;
157  }
158  }
159  }
160  else if (x < -epsilon)
161  {
162  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
163  mabs = (-x > mabs) ? -x : mabs;
164  l = low[i];
165  if (l > -infinity)
166  {
167  y = l - vec[i];
168  if ( y >= 0 )
169  x = - fastDelta / x;
170  else
171  x = ( y - fastDelta ) / x;
172 
173  if (x < max)
174  {
175  max = x;
176  sel = i;
177  }
178  }
179  }
180  }
181  }
182  else
183  {
184  /* In this case, the indices of the semi-sparse vector update.delta() are not set up and are filled below. */
185  int* l_idx = update.delta().altIndexMem();
186  Real* uval = update.delta().altValues();
187  const Real* uend = uval + update.dim();
188  assert( uval == upd );
189 
190  for (i = 0; uval < uend; ++uval, ++i)
191  {
192  x = *uval;
193  if (x)
194  {
195  if (x >= -epsilon && x <= epsilon)
196  {
197  *uval = 0.0;
198  continue;
199  }
200  else
201  *l_idx++ = i;
202 
203  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
204  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
205  continue;
206 
207  if (x > epsilon)
208  {
209  mabs = (x > mabs) ? x : mabs;
210  u = up[i];
211  if (u < infinity)
212  {
213  y = u - vec[i];
214  if (y <= 0)
215  x = fastDelta / x;
216  else
217  x = (y + fastDelta) / x;
218 
219  if (x < max)
220  {
221  max = x;
222  sel = i;
223  }
224  }
225  }
226  else if (x < -epsilon)
227  {
228  mabs = (-x > mabs) ? -x : mabs;
229  l = low[i];
230  if (l > -infinity)
231  {
232  y = l - vec[i];
233  if ( y >= 0 )
234  x = - fastDelta / x;
235  else
236  x = ( y - fastDelta ) / x;
237 
238  if (x < max)
239  {
240  max = x;
241  sel = i;
242  }
243  }
244  }
245  }
246  }
247  update.delta().setSize(int(l_idx - update.delta().indexMem()));
248  update.delta().forceSetup();
249  }
250 
251  val = max;
252  maxabs = mabs;
253 
254  return sel;
255 }
256 
257 /* See maxDelta() */
258 int SPxFastRT::minDelta(
259  Real& val,
260  Real& maxabs,
261  UpdateVector& update,
262  const Vector& lowBound,
263  const Vector& upBound,
264  int start,
265  int incr) const
266 {
267  int i, sel;
268  Real x, y, max;
269  Real u, l;
270  bool leaving = m_type == SPxSolver::LEAVE;
271 
272  Real mabs = maxabs;
273 
274  const Real* up = upBound.get_const_ptr();
275  const Real* low = lowBound.get_const_ptr();
276  const Real* vec = update.get_const_ptr();
277  const Real* upd = update.delta().values();
278  const int* idx = update.delta().indexMem();
279 
280  sel = -1;
281  max = val;
282 
283  if (update.delta().isSetup())
284  {
285  const int* last = idx + update.delta().size();
286  for (idx += start; idx < last; idx += incr)
287  {
288  i = *idx;
289  x = upd[i];
290 
291  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
292  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
293  continue;
294 
295  if (x > epsilon)
296  {
297  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
298  mabs = (x > mabs) ? x : mabs;
299  l = low[i];
300  if (l > -infinity)
301  {
302  y = l - vec[i];
303  if ( y >= 0 )
304  x = - fastDelta / x;
305  else
306  x = ( y - fastDelta ) / x;
307 
308  if (x > max)
309  {
310  max = x;
311  sel = i;
312  }
313  }
314  }
315  else if (x < -epsilon)
316  {
317  // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
318  mabs = (-x > mabs) ? -x : mabs;
319  u = up[i];
320  if (u < infinity)
321  {
322  y = u - vec[i];
323  if (y <= 0)
324  x = fastDelta / x;
325  else
326  x = (y + fastDelta) / x;
327 
328  if (x > max)
329  {
330  max = x;
331  sel = i;
332  }
333  }
334  }
335  }
336  }
337  else
338  {
339  /* In this case, the indices of the semi-sparse vector update.delta() are not set up and are filled below. */
340  int* l_idx = update.delta().altIndexMem();
341  Real* uval = update.delta().altValues();
342  const Real* uend = uval + update.dim();
343  assert( uval == upd );
344 
345  for (i = 0; uval < uend; ++uval, ++i)
346  {
347  x = *uval;
348 
349  if (x)
350  {
351  if (x >= -epsilon && x <= epsilon)
352  {
353  *uval = 0.0;
354  continue;
355  }
356  else
357  *l_idx++ = i;
358 
359  /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
360  if( leaving && ((iscoid && thesolver->isCoBasic(i)) || (!iscoid && thesolver->isBasic(i))) )
361  continue;
362 
363  if (x > epsilon)
364  {
365  mabs = (x > mabs) ? x : mabs;
366  l = low[i];
367  if (l > -infinity)
368  {
369  y = l - vec[i];
370  if ( y >= 0 )
371  x = - fastDelta / x;
372  else
373  x = ( y - fastDelta ) / x;
374 
375  if (x > max)
376  {
377  max = x;
378  sel = i;
379  }
380  }
381  }
382  else if (x < -epsilon)
383  {
384  mabs = (-x > mabs) ? -x : mabs;
385  u = up[i];
386  if (u < infinity)
387  {
388  y = u - vec[i];
389  if (y <= 0)
390  x = fastDelta / x;
391  else
392  x = (y + fastDelta) / x;
393 
394  if (x > max)
395  {
396  max = x;
397  sel = i;
398  }
399  }
400  }
401  }
402  }
403  update.delta().setSize(int(l_idx - update.delta().indexMem()));
404  update.delta().forceSetup();
405  }
406 
407  val = max;
408  maxabs = mabs;
409 
410  return sel;
411 }
412 
413 int SPxFastRT::maxDelta(
414  Real& val,
415  Real& maxabs)
416 {
417  return maxDelta(val, maxabs,
418  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
419 }
420 
421 int SPxFastRT::minDelta(
422  Real& val,
423  Real& maxabs)
424 {
425  return minDelta(val, maxabs,
426  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
427 }
428 
429 SPxId SPxFastRT::maxDelta(
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,
438  thesolver->coPvec(), thesolver->lcBound(), thesolver->ucBound(), 0, 1);
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 
457 SPxId SPxFastRT::minDelta(
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,
466  thesolver->coPvec(), thesolver->lcBound(), thesolver->ucBound(), 0, 1);
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. */
488 int SPxFastRT::minSelect(
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() */
567 int SPxFastRT::maxSelect(
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 
646 int SPxFastRT::maxSelect(
647  Real& val,
648  Real& stab,
649  Real& bestDelta,
650  Real max)
651 {
652  Real best = -infinity;
653  bestDelta = 0.0;
654  return maxSelect(val, stab, best, bestDelta, max,
655  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
656 }
657 
658 SPxId SPxFastRT::maxSelect(
659  int& nr,
660  Real& val,
661  Real& stab,
662  Real& bestDelta,
663  Real max
664 )
665 {
666  int indp, indc;
667  Real best = -infinity;
668  bestDelta = 0.0;
669  iscoid = true;
670  indc = maxSelect(val, stab, best, bestDelta, max,
671  thesolver->coPvec(), thesolver->lcBound(), thesolver->ucBound(), 0, 1);
672  iscoid = false;
673  indp = maxSelect(val, stab, best, bestDelta, max,
674  thesolver->pVec(), thesolver->lpBound(), thesolver->upBound(), 0, 1);
675 
676  if (indp >= 0)
677  {
678  nr = indp;
679  return thesolver->id(indp);
680  }
681  if (indc >= 0)
682  {
683  nr = indc;
684  return thesolver->coId(indc);
685  }
686  nr = -1;
687  return SPxId();
688 }
689 
690 int SPxFastRT::minSelect(
691  Real& val,
692  Real& stab,
693  Real& bestDelta,
694  Real max)
695 {
696  Real best = infinity;
697  bestDelta = 0.0;
698  return minSelect(val, stab, best, bestDelta, max,
699  thesolver->fVec(), thesolver->lbBound(), thesolver->ubBound(), 0, 1);
700 }
701 
702 SPxId SPxFastRT::minSelect(
703  int& nr,
704  Real& val,
705  Real& stab,
706  Real& bestDelta,
707  Real max)
708 {
709  Real best = infinity;
710  bestDelta = 0.0;
711  iscoid = true;
712  int indc = minSelect(val, stab, best, bestDelta, max,
713  thesolver->coPvec(), thesolver->lcBound(), thesolver->ucBound(), 0, 1);
714  iscoid = false;
715  int indp = minSelect(val, stab, best, bestDelta, max,
716  thesolver->pVec(), thesolver->lpBound(), thesolver->upBound(), 0, 1);
717 
718  if (indp >= 0)
719  {
720  nr = indp;
721  return thesolver->id(indp);
722  }
723  if (indc >= 0)
724  {
725  nr = indc;
726  return thesolver->coId(indc);
727  }
728  nr = -1;
729  return SPxId();
730 }
731 
732 
733 bool SPxFastRT::maxShortLeave(Real& sel, int leave, Real maxabs)
734 {
735  assert(leave >= 0);
736  assert(maxabs >= 0);
737 
738  sel = thesolver->fVec().delta()[leave];
739 
740  if (sel > maxabs*SHORT)
741  {
742  sel = (thesolver->ubBound()[leave] - thesolver->fVec()[leave]) / sel;
743  return true;
744  }
745 
746  if (sel < -maxabs*SHORT)
747  {
748  sel = (thesolver->lbBound()[leave] - thesolver->fVec()[leave]) / sel;
749  return true;
750  }
751 
752  return false;
753 }
754 
755 bool SPxFastRT::minShortLeave(Real& sel, int leave, Real maxabs)
756 {
757  assert(leave >= 0);
758  assert(maxabs >= 0);
759 
760  sel = thesolver->fVec().delta()[leave];
761 
762  if (sel > maxabs*SHORT)
763  {
764  sel = (thesolver->lbBound()[leave] - thesolver->fVec()[leave]) / sel;
765  return true;
766  }
767 
768  if ( sel < -maxabs*SHORT)
769  {
770  sel = (thesolver->ubBound()[leave] - thesolver->fVec()[leave]) / sel;
771  return true;
772  }
773 
774  return false;
775 }
776 
777 bool SPxFastRT::maxReLeave(Real& sel, int leave, Real maxabs)
778 {
779  UpdateVector& vec = thesolver->fVec();
780  Vector& low = thesolver->lbBound();
781  Vector& up = thesolver->ubBound();
782 
783  if (leave < 0)
784  return true;
785 
786  if (up[leave] > low[leave])
787  {
788  Real x = vec.delta()[leave];
789 
790  if (sel < -fastDelta / maxabs)
791  {
792  sel = 0.0;
793  if (x < 0.0)
794  thesolver->shiftLBbound(leave, vec[leave]);
795  else
796  thesolver->shiftUBbound(leave, vec[leave]);
797  }
798  }
799  else
800  {
801  sel = 0.0;
802  thesolver->shiftLBbound(leave, vec[leave]);
803  thesolver->shiftUBbound(leave, vec[leave]);
804  }
805 
806  return false;
807 }
808 
809 bool SPxFastRT::minReLeave(Real& sel, int leave, Real maxabs)
810 {
811  UpdateVector& vec = thesolver->fVec();
812  Vector& low = thesolver->lbBound();
813  Vector& up = thesolver->ubBound();
814 
815  if (leave < 0)
816  return true;
817 
818  if (up[leave] > low[leave])
819  {
820  Real x = vec.delta()[leave];
821 
822  if (sel > fastDelta / maxabs)
823  {
824  if (x > 0.0)
825  {
826  thesolver->theShift += low[leave];
827  sel = 0.0;
828  low[leave] = vec[leave] + sel * x;
829  thesolver->theShift -= low[leave];
830  }
831  else
832  {
833  thesolver->theShift -= up[leave];
834  sel = 0.0;
835  up[leave] = vec[leave] + sel * x;
836  thesolver->theShift += up[leave];
837  }
838  }
839  }
840  else
841  {
842  sel = 0.0;
843  if (vec[leave] < low[leave])
844  thesolver->theShift += low[leave] - vec[leave];
845  else
846  thesolver->theShift += vec[leave] - up[leave];
847  low[leave] = up[leave] = vec[leave];
848  }
849 
850  return false;
851 }
852 
853 int SPxFastRT::selectLeave(Real& val, SPxId)
854 {
855  Real maxabs, max, sel;
856  int leave = -1;
857  int cnt = 0;
858 
859  assert( m_type == SPxSolver::ENTER );
860 
861  resetTols();
862 
863  if (val > epsilon)
864  {
865  do
866  {
867  // phase 1:
868  max = val;
869  maxabs = 0.0;
870  leave = maxDelta(max, maxabs);
871  if (max == val)
872  return -1;
873 
874  if (!maxShortLeave(sel, leave, maxabs))
875  {
876  // phase 2:
877  Real stab, bestDelta;
878  stab = 100.0 * minStability(maxabs);
879  leave = maxSelect(sel, stab, bestDelta, max);
880  if (bestDelta < DELTA_SHIFT*TRIES)
881  cnt++;
882  else
883  cnt += TRIES;
884  }
885  if (!maxReLeave(sel, leave, maxabs))
886  break;
887  relax();
888  }
889  while (cnt < TRIES);
890  }
891  else if (val < -epsilon)
892  {
893  do
894  {
895  max = val;
896  maxabs = 0;
897  leave = minDelta(max, maxabs);
898 
899  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) || thesolver->desc().colStatus(thesolver->number(SPxColId(thesolver->baseId(leave)))) != SPxBasis::Desc::P_FIXED);
900 
901  if (max == val)
902  return -1;
903 
904  if (!minShortLeave(sel, leave, maxabs))
905  {
906  // phase 2:
907  Real stab, bestDelta;
908  stab = 100.0 * minStability(maxabs);
909  leave = minSelect(sel, stab, bestDelta, max);
910 
911  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) || thesolver->desc().colStatus(thesolver->number(SPxColId(thesolver->baseId(leave)))) != SPxBasis::Desc::P_FIXED);
912 
913  if (bestDelta < DELTA_SHIFT*TRIES)
914  cnt++;
915  else
916  cnt += TRIES;
917  }
918  if (!minReLeave(sel, leave, maxabs))
919  break;
920  relax();
921  }
922  while (cnt < TRIES);
923  }
924  else
925  return -1;
926 
927  MSG_DEBUG(
928  if (leave >= 0)
929  spxout
930  << "DFSTRT01 "
931  << thesolver->basis().iteration() << "("
932  << std::setprecision(6) << thesolver->value() << ","
933  << std::setprecision(2) << thesolver->basis().stability() << "):"
934  << leave << "\t"
935  << std::setprecision(4) << sel << " "
936  << std::setprecision(4) << thesolver->fVec().delta()[leave] << " "
937  << std::setprecision(6) << maxabs
938  << std::endl;
939  else
940  spxout << "DFSTRT02 " << thesolver->basis().iteration()
941  << ": skipping instable pivot" << std::endl;
942  )
943 
944  if (leave >= 0 || minStab > 2*solver()->epsilon())
945  {
946  val = sel;
947  if (leave >= 0)
948  tighten();
949  }
950 
951  assert(leave < 0 || !(thesolver->baseId(leave).isSPxColId()) || thesolver->desc().colStatus(thesolver->number(SPxColId(thesolver->baseId(leave)))) != SPxBasis::Desc::P_FIXED);
952 
953  return leave;
954 }
955 
956 
957 bool SPxFastRT::maxReEnter(
958  Real& sel,
959  Real maxabs,
960  const SPxId& id,
961  int nr)
962 {
963  Real x, d;
964  Vector* up;
965  Vector* low;
966 
967  UpdateVector& pvec = thesolver->pVec();
968  SSVector& pupd = thesolver->pVec().delta();
969  Vector& upb = thesolver->upBound();
970  Vector& lpb = thesolver->lpBound();
971  UpdateVector& cvec = thesolver->coPvec();
972  SSVector& cupd = thesolver->coPvec().delta();
973  Vector& ucb = thesolver->ucBound();
974  Vector& lcb = thesolver->lcBound();
975 
976  if (thesolver->isCoId(id))
977  {
978  if (thesolver->isCoBasic(nr))
979  {
980  cupd.clearIdx(nr);
981  return true;
982  }
983 
984  x = cvec[nr];
985  d = cupd[nr];
986  up = &ucb;
987  low = &lcb;
988 
989  if (d < 0.0)
990  sel = (lcb[nr] - cvec[nr]) / d;
991  else
992  sel = (ucb[nr] - cvec[nr]) / d;
993  }
994  else if (thesolver->isId(id))
995  {
996  pvec[nr] = thesolver->vector(nr) * cvec;
997  if (thesolver->isBasic(nr))
998  {
999  pupd.clearIdx(nr);
1000  return true;
1001  }
1002 
1003  x = pvec[nr];
1004  d = pupd[nr];
1005  up = &upb;
1006  low = &lpb;
1007 
1008  if (d < 0.0)
1009  sel = (lpb[nr] - pvec[nr]) / d;
1010  else
1011  sel = (upb[nr] - pvec[nr]) / d;
1012  }
1013  else
1014  return true;
1015 
1016  if ((*up)[nr] != (*low)[nr])
1017  {
1018  if (sel < -fastDelta / maxabs)
1019  {
1020  if (d > 0.0)
1021  {
1022  thesolver->theShift -= (*up)[nr];
1023  sel = 0.0;
1024  (*up)[nr] = x + sel * d;
1025  thesolver->theShift += (*up)[nr];
1026  }
1027  else
1028  {
1029  thesolver->theShift += (*low)[nr];
1030  sel = 0.0;
1031  (*low)[nr] = x + sel * d;
1032  thesolver->theShift -= (*low)[nr];
1033  }
1034  }
1035  }
1036  else
1037  {
1038  sel = 0.0;
1039  if (x > (*up)[nr])
1040  thesolver->theShift += x - (*up)[nr];
1041  else
1042  thesolver->theShift += (*low)[nr] - x;
1043  (*up)[nr] = (*low)[nr] = x;
1044  }
1045 
1046  return false;
1047 }
1048 
1049 bool SPxFastRT::minReEnter(
1050  Real& sel,
1051  Real maxabs,
1052  const SPxId& id,
1053  int nr)
1054 {
1055  Real x, d;
1056  Vector* up;
1057  Vector* low;
1058 
1059  UpdateVector& pvec = thesolver->pVec();
1060  SSVector& pupd = thesolver->pVec().delta();
1061  Vector& upb = thesolver->upBound();
1062  Vector& lpb = thesolver->lpBound();
1063  UpdateVector& cvec = thesolver->coPvec();
1064  SSVector& cupd = thesolver->coPvec().delta();
1065  Vector& ucb = thesolver->ucBound();
1066  Vector& lcb = thesolver->lcBound();
1067 
1068  if (thesolver->isCoId(id))
1069  {
1070  if (thesolver->isCoBasic(nr))
1071  {
1072  cupd.clearIdx(nr);
1073  return true;
1074  }
1075  x = cvec[nr];
1076  d = cupd[nr];
1077  up = &ucb;
1078  low = &lcb;
1079  if (d > 0.0)
1080  sel = (thesolver->lcBound()[nr] - cvec[nr]) / d;
1081  else
1082  sel = (thesolver->ucBound()[nr] - cvec[nr]) / d;
1083  }
1084 
1085  else if (thesolver->isId(id))
1086  {
1087  pvec[nr] = thesolver->vector(nr) * cvec;
1088  if (thesolver->isBasic(nr))
1089  {
1090  pupd.clearIdx(nr);
1091  return true;
1092  }
1093  x = pvec[nr];
1094  d = pupd[nr];
1095  up = &upb;
1096  low = &lpb;
1097  if (d > 0.0)
1098  sel = (thesolver->lpBound()[nr] - pvec[nr]) / d;
1099  else
1100  sel = (thesolver->upBound()[nr] - pvec[nr]) / d;
1101  }
1102 
1103  else
1104  return true;
1105 
1106  if ((*up)[nr] != (*low)[nr])
1107  {
1108  if (sel > fastDelta / maxabs)
1109  {
1110  if (d < 0.0)
1111  {
1112  thesolver->theShift -= (*up)[nr];
1113  sel = 0.0;
1114  (*up)[nr] = x + sel * d;
1115  thesolver->theShift += (*up)[nr];
1116  }
1117  else
1118  {
1119  thesolver->theShift += (*low)[nr];
1120  sel = 0.0;
1121  (*low)[nr] = x + sel * d;
1122  thesolver->theShift -= (*low)[nr];
1123  }
1124  }
1125  }
1126  else
1127  {
1128  sel = 0.0;
1129  if (x > (*up)[nr])
1130  thesolver->theShift += x - (*up)[nr];
1131  else
1132  thesolver->theShift += (*low)[nr] - x;
1133  (*up)[nr] = (*low)[nr] = x;
1134  }
1135 
1136  return false;
1137 }
1138 
1139 bool SPxFastRT::shortEnter(
1140  const SPxId& enterId,
1141  int nr,
1142  Real max,
1143  Real maxabs) const
1144 {
1145  if (thesolver->isCoId(enterId))
1146  {
1147  if (max != 0.0)
1148  {
1149  Real x = thesolver->coPvec().delta()[nr];
1150  if (x < maxabs * SHORT && -x < maxabs * SHORT)
1151  return false;
1152  }
1153  return true;
1154  }
1155  else if (thesolver->isId(enterId))
1156  {
1157  if (max != 0.0)
1158  {
1159  Real x = thesolver->pVec().delta()[nr];
1160  if (x < maxabs * SHORT && -x < maxabs * SHORT)
1161  return false;
1162  }
1163  return true;
1164  }
1165 
1166  return false;
1167 }
1168 
1169 SPxId SPxFastRT::selectEnter(Real& val, int)
1170 {
1171  SPxId enterId;
1172  Real max, sel;
1173  Real maxabs = 0.0;
1174  int nr;
1175  int cnt = 0;
1176 
1177  assert( m_type == SPxSolver::LEAVE );
1178 
1179  // force instable pivot iff true (see explanation in leave.cpp and spxsolve.cpp)
1180  bool instable = solver()->instableLeave;
1181  Real lowstab = LOWSTAB;
1182  assert(!instable || solver()->instableLeaveNum >= 0);
1183 
1184  resetTols();
1185  sel = 0.0;
1186 
1187  if (val > epsilon)
1188  {
1189  do
1190  {
1191  maxabs = 0.0;
1192  max = val;
1193 
1194  enterId = maxDelta(nr, max, maxabs);
1195  if (!enterId.isValid())
1196  return enterId;
1197  assert(max >= 0.0);
1198 
1199  if (!shortEnter(enterId, nr, max, maxabs))
1200  {
1201  Real bestDelta, stab;
1202 
1203  stab = minStability(maxabs);
1204 
1205  // force instable pivot iff instable is true (see explanation in leave.cpp and spxsolve.cpp)
1206  if (instable)
1207  {
1208  enterId = maxSelect(nr, sel, lowstab, bestDelta, max);
1209  }
1210  else
1211  {
1212  enterId = maxSelect(nr, sel, stab, bestDelta, max);
1213  }
1214  if (bestDelta < DELTA_SHIFT*TRIES)
1215  cnt++;
1216  else
1217  cnt += TRIES;
1218  }
1219  if (!maxReEnter(sel, maxabs, enterId, nr))
1220  break;
1221  relax();
1222  }
1223  while (cnt < TRIES);
1224  }
1225  else if (val < -epsilon)
1226  {
1227  do
1228  {
1229  maxabs = 0.0;
1230  max = val;
1231  enterId = minDelta(nr, max, maxabs);
1232  if (!enterId.isValid())
1233  return enterId;
1234  assert(max <= 0.0);
1235 
1236  if (!shortEnter(enterId, nr, max, maxabs))
1237  {
1238  Real bestDelta, stab;
1239 
1240  stab = minStability(maxabs);
1241 
1242  // force instable pivot iff instable is true (see explanation in leave.cpp and spxsolve.cpp)
1243  if (instable)
1244  {
1245  enterId = minSelect(nr, sel, lowstab, bestDelta, max);
1246  }
1247  else
1248  {
1249  enterId = minSelect(nr, sel, stab, bestDelta, max);
1250  }
1251  if (bestDelta < DELTA_SHIFT*TRIES)
1252  cnt++;
1253  else
1254  cnt += TRIES;
1255  }
1256  if (!minReEnter(sel, maxabs, enterId, nr))
1257  break;
1258  relax();
1259  }
1260  while (cnt < TRIES);
1261  }
1262 
1263  MSG_DEBUG(
1264  if (enterId.isValid())
1265  {
1266  Real x;
1267  if (thesolver->isCoId(enterId))
1268  x = thesolver->coPvec().delta()[ thesolver->number(enterId) ];
1269  else
1270  x = thesolver->pVec().delta()[ thesolver->number(enterId) ];
1271  spxout << "DFSTRT03 " << thesolver->basis().iteration() << ": "
1272  << sel << '\t' << x << " (" << maxabs << ")" << std::endl;
1273  }
1274  else
1275  spxout << "DFSTRT04 " << thesolver->basis().iteration()
1276  << ": skipping instable pivot" << std::endl;
1277  )
1278 
1279 
1280  if (enterId.isValid() || minStab > 2*epsilon)
1281  {
1282  val = sel;
1283  if (enterId.isValid())
1284  tighten();
1285  }
1286 
1287  return enterId;
1288 }
1289 
1290 void SPxFastRT::load(SPxSolver* spx)
1291 {
1292  thesolver = spx;
1293  setType(spx->type());
1294 }
1295 
1296 void SPxFastRT::setType(SPxSolver::Type type)
1297 {
1298  m_type = type;
1299 
1300  minStab = MINSTAB;
1301  fastDelta = delta;
1302 }
1303 } // namespace soplex
1304 
1305 //-----------------------------------------------------------------------------
1306 //Emacs Local Variables:
1307 //Emacs mode:c++
1308 //Emacs c-basic-offset:3
1309 //Emacs tab-width:8
1310 //Emacs indent-tabs-mode:nil
1311 //Emacs End:
1312 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen