SoPlex Doxygen Documentation
spxharrisrt.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 <iostream>
20 
21 #include "spxdefines.h"
22 #include "spxharrisrt.h"
23 
24 namespace soplex
25 {
26 /**@todo suspicious: *max is not set, but it is used
27  * (with the default setting *max=1) in selectLeave and selectEnter
28  * The question might be if max shouldn't be updated with themax?
29  *
30  * numCycle and maxCycle are integers. So degeneps will be
31  * exactly delta until numCycle >= maxCycle. Then it will be
32  * 0 until numCycle >= 2 * maxCycle, after wich it becomes
33  * negative. This does not look ok.
34  */
35 Real SPxHarrisRT::degenerateEps() const
36 {
37  return solver()->delta()
38  * (1.0 - solver()->numCycle() / solver()->maxCycle());
39 }
40 
41 int SPxHarrisRT::maxDelta(
42  Real* /*max*/, /* max abs value in upd */
43  Real* val, /* initial and chosen value */
44  int num, /* # of indices in idx */
45  const int* idx, /* nonzero indices in upd */
46  const Real* upd, /* update vector for vec */
47  const Real* vec, /* current vector */
48  const Real* low, /* lower bounds for vec */
49  const Real* up, /* upper bounds for vec */
50  Real epsilon) const /* what is 0? */
51 {
52  Real x;
53  Real theval;
54  /**@todo patch suggests using *max instead of themax */
55  Real themax;
56  int sel;
57  int i;
58 
59  assert(*val >= 0);
60 
61  theval = *val;
62  themax = 0;
63  sel = -1;
64 
65  while (num--)
66  {
67  i = idx[num];
68  x = upd[i];
69  if (x > epsilon)
70  {
71  themax = (x > themax) ? x : themax;
72  x = (up[i] - vec[i] + delta) / x;
73  if (x < theval && up[i] < infinity)
74  theval = x;
75  }
76  else if (x < -epsilon)
77  {
78  themax = (-x > themax) ? -x : themax;
79  x = (low[i] - vec[i] - delta) / x;
80  if (x < theval && low[i] > -infinity)
81  theval = x;
82  }
83  }
84  *val = theval;
85  return sel;
86 }
87 
88 /**@todo suspicious: *max is not set, but it is used
89  (with the default setting *max=1)
90  in selectLeave and selectEnter
91 */
92 int SPxHarrisRT::minDelta(
93  Real* /*max*/, /* max abs value in upd */
94  Real* val, /* initial and chosen value */
95  int num, /* # of indices in idx */
96  const int* idx, /* nonzero indices in upd */
97  const Real* upd, /* update vector for vec */
98  const Real* vec, /* current vector */
99  const Real* low, /* lower bounds for vec */
100  const Real* up, /* upper bounds for vec */
101  Real epsilon) const /* what is 0? */
102 {
103  Real x;
104  Real theval;
105  /**@todo patch suggests using *max instead of themax */
106  Real themax;
107  int sel;
108  int i;
109 
110  assert(*val < 0);
111 
112  theval = *val;
113  themax = 0;
114  sel = -1;
115 
116  while (num--)
117  {
118  i = idx[num];
119  x = upd[i];
120  if (x > epsilon)
121  {
122  themax = (x > themax) ? x : themax;
123  x = (low[i] - vec[i] - delta) / x;
124  if (x > theval && low[i] > -infinity)
125  theval = x;
126  }
127  else if (x < -epsilon)
128  {
129  themax = (-x > themax) ? -x : themax;
130  x = (up[i] - vec[i] + delta) / x;
131  if (x > theval && up[i] < infinity)
132  theval = x;
133  }
134  }
135  *val = theval;
136  return sel;
137 }
138 
139 /**
140  Here comes our implementation of the Harris procedure improved by shifting
141  bounds. The basic idea is to used the tolerated infeasibility within
142  solver()->entertol() for searching numerically stable pivots.
143 
144  The algorithms operates in two phases. In a first phase, the maximum \p val
145  is determined, when infeasibility within solver()->entertol() is allowed. In the second
146  phase, between all variables with values < \p val the one is selected which
147  gives the best step forward in the simplex iteration. However, this may not
148  allways yield an improvement. In that case, we shift the variable toward
149  infeasibility and retry. This avoids cycling in the shifted LP.
150  */
151 int SPxHarrisRT::selectLeave(Real& val, SPxId enterId)
152 {
153  int i, j;
154  Real stab, x, y;
155  Real max;
156  Real sel;
157  Real lastshift;
158  Real useeps;
159  int leave = -1;
160  Real maxabs = 1;
161 
162  Real epsilon = solver()->epsilon();
163  Real degeneps = degenerateEps();
164 
165  SSVector& upd = solver()->fVec().delta();
166  Vector& vec = solver()->fVec();
167 
168  const Vector& up = solver()->ubBound();
169  const Vector& low = solver()->lbBound();
170 
171  assert(delta > epsilon);
172  assert(epsilon > 0);
173  assert(solver()->maxCycle() > 0);
174 
175  max = val;
176  lastshift = solver()->shift();
177 
178  solver()->fVec().delta().setup();
179 
180  if (max > epsilon)
181  {
182  // phase 1:
183  maxDelta(
184  &maxabs, /* max abs value in upd */
185  &max, /* initial and chosen value */
186  upd.size(), /* # of indices in upd */
187  upd.indexMem(), /* nonzero indices in upd */
188  upd.values(), /* update vector for vec */
189  vec.get_const_ptr(), /* current vector */
190  low.get_const_ptr(), /* lower bounds for vec */
191  up.get_const_ptr(), /* upper bounds for vec */
192  epsilon); /* what is 0? */
193 
194  if (max == val)
195  return -1;
196 
197 
198  // phase 2:
199  stab = 0;
200  sel = -infinity;
201  useeps = maxabs * epsilon * 0.001;
202  if (useeps < epsilon)
203  useeps = epsilon;
204  for (j = upd.size() - 1; j >= 0; --j)
205  {
206  i = upd.index(j);
207  x = upd[i];
208  if (x > useeps)
209  {
210  y = up[i] - vec[i];
211  if (y < -degeneps)
212  solver()->shiftUBbound(i, vec[i]); // ensure simplex improvement
213  else
214  {
215  y /= x;
216  if (y <= max && y > sel - epsilon && x > stab)
217  {
218  sel = y;
219  leave = i;
220  stab = x;
221  }
222  }
223  }
224  else if (x < -useeps)
225  {
226  y = low[i] - vec[i];
227  if (y > degeneps)
228  solver()->shiftLBbound(i, vec[i]); // ensure simplex improvement
229  else
230  {
231  y /= x;
232  if (y <= max && y > sel - epsilon && -x > stab)
233  {
234  sel = y;
235  leave = i;
236  stab = -x;
237  }
238  }
239  }
240  else
241  upd.clearNum(j);
242  }
243  }
244 
245 
246  else if (max < -epsilon)
247  {
248  // phase 1:
249  minDelta(
250  &maxabs, /* max abs value in upd */
251  &max, /* initial and chosen value */
252  upd.size(), /* # of indices in upd */
253  upd.indexMem(), /* nonzero indices in upd */
254  upd.values(), /* update vector for vec */
255  vec.get_const_ptr(), /* current vector */
256  low.get_const_ptr(), /* lower bounds for vec */
257  up.get_const_ptr(), /* upper bounds for vec */
258  epsilon); /* what is 0? */
259 
260  if (max == val)
261  return -1;
262 
263  // phase 2:
264  stab = 0;
265  sel = infinity;
266  useeps = maxabs * epsilon * 0.001;
267  if (useeps < epsilon)
268  useeps = epsilon;
269  for (j = upd.size() - 1; j >= 0; --j)
270  {
271  i = upd.index(j);
272  x = upd[i];
273  if (x < -useeps)
274  {
275  y = up[i] - vec[i];
276  if (y < -degeneps)
277  solver()->shiftUBbound(i, vec[i]); // ensure simplex improvement
278  else
279  {
280  y /= x;
281  if (y >= max && y < sel + epsilon && -x > stab)
282  {
283  sel = y;
284  leave = i;
285  stab = -x;
286  }
287  }
288  }
289  else if (x > useeps)
290  {
291  y = low[i] - vec[i];
292  if (y > degeneps)
293  solver()->shiftLBbound(i, vec[i]); // ensure simplex improvement
294  else
295  {
296  y /= x;
297  if (y >= max && y < sel + epsilon && x > stab)
298  {
299  sel = y;
300  leave = i;
301  stab = x;
302  }
303  }
304  }
305  else
306  upd.clearNum(j);
307  }
308  }
309 
310  else
311  return -1;
312 
313 
314  if (lastshift != solver()->shift())
315  return selectLeave(val, enterId);
316 
317  assert(leave >= 0);
318 
319  val = sel;
320  return leave;
321 }
322 
323 
324 SPxId SPxHarrisRT::selectEnter(Real& val, int)
325 {
326  int i, j;
327  SPxId enterId;
328  Real stab, x, y;
329  Real max = 0.0;
330  Real sel = 0.0;
331  Real lastshift;
332  Real cuseeps;
333  Real ruseeps;
334  Real cmaxabs = 1;
335  Real rmaxabs = 1;
336  int pnr, cnr;
337 
338  Real minStability = 0.0001;
339  Real epsilon = solver()->epsilon();
340  Real degeneps = degenerateEps();
341 
342  Vector& pvec = solver()->pVec();
343  SSVector& pupd = solver()->pVec().delta();
344 
345  Vector& cvec = solver()->coPvec();
346  SSVector& cupd = solver()->coPvec().delta();
347 
348  const Vector& upb = solver()->upBound();
349  const Vector& lpb = solver()->lpBound();
350  const Vector& ucb = solver()->ucBound();
351  const Vector& lcb = solver()->lcBound();
352 
353  assert(delta > epsilon);
354  assert(epsilon > 0);
355  assert(solver()->maxCycle() > 0);
356 
357  solver()->coPvec().delta().setup();
358  solver()->pVec().delta().setup();
359 
360  if (val > epsilon)
361  {
362  for(;;)
363  {
364  pnr = -1;
365  cnr = -1;
366  max = val;
367  sel = val;
368  lastshift = solver()->shift();
369  assert(delta > epsilon);
370 
371  // phase 1:
372  maxDelta(
373  &rmaxabs, /* max abs value in upd */
374  &max, /* initial and chosen value */
375  pupd.size(), /* # of indices in pupd */
376  pupd.indexMem(), /* nonzero indices in pupd */
377  pupd.values(), /* update vector for vec */
378  pvec.get_const_ptr(), /* current vector */
379  lpb.get_const_ptr(), /* lower bounds for vec */
380  upb.get_const_ptr(), /* upper bounds for vec */
381  epsilon); /* what is 0? */
382 
383  maxDelta(
384  &cmaxabs, /* max abs value in upd */
385  &max, /* initial and chosen value */
386  cupd.size(), /* # of indices in cupd */
387  cupd.indexMem(), /* nonzero indices in cupd */
388  cupd.values(), /* update vector for vec */
389  cvec.get_const_ptr(), /* current vector */
390  lcb.get_const_ptr(), /* lower bounds for vec */
391  ucb.get_const_ptr(), /* upper bounds for vec */
392  epsilon); /* what is 0? */
393 
394  if (max == val)
395  return enterId;
396 
397 
398  // phase 2:
399  stab = 0;
400  sel = -infinity;
401  ruseeps = rmaxabs * 0.001 * epsilon;
402  if (ruseeps < epsilon)
403  ruseeps = epsilon;
404  cuseeps = cmaxabs * 0.001 * epsilon;
405  if (cuseeps < epsilon)
406  cuseeps = epsilon;
407  for (j = pupd.size() - 1; j >= 0; --j)
408  {
409  i = pupd.index(j);
410  x = pupd[i];
411  if (x > ruseeps)
412  {
413  y = upb[i] - pvec[i];
414  if (y < -degeneps)
415  solver()->shiftUPbound(i, pvec[i] - degeneps);
416  else
417  {
418  y /= x;
419  if (y <= max && x >= stab) // && y > sel-epsilon
420  {
421  enterId = solver()->id(i);
422  sel = y;
423  pnr = i;
424  stab = x;
425  }
426  }
427  }
428  else if (x < -ruseeps)
429  {
430  y = lpb[i] - pvec[i];
431  if (y > degeneps)
432  solver()->shiftLPbound(i, pvec[i] + degeneps);
433  else
434  {
435  y /= x;
436  if (y <= max && -x >= stab) // && y > sel-epsilon
437  {
438  enterId = solver()->id(i);
439  sel = y;
440  pnr = i;
441  stab = -x;
442  }
443  }
444  }
445  else
446  {
447  MSG_DEBUG( spxout << "DHARRI01 removing value " << pupd[i] << std::endl; )
448  pupd.clearNum(j);
449  }
450  }
451  for (j = cupd.size() - 1; j >= 0; --j)
452  {
453  i = cupd.index(j);
454  x = cupd[i];
455  if (x > cuseeps)
456  {
457  y = ucb[i] - cvec[i];
458  if (y < -degeneps)
459  solver()->shiftUCbound(i, cvec[i] - degeneps);
460  else
461  {
462  y /= x;
463  if (y <= max && x >= stab) // && y > sel-epsilon
464  {
465  enterId = solver()->coId(i);
466  sel = y;
467  cnr = j;
468  stab = x;
469  }
470  }
471  }
472  else if (x < -cuseeps)
473  {
474  y = lcb[i] - cvec[i];
475  if (y > degeneps)
476  solver()->shiftLCbound(i, cvec[i] + degeneps);
477  else
478  {
479  y /= x;
480  if (y <= max && -x >= stab) // && y > sel-epsilon
481  {
482  enterId = solver()->coId(i);
483  sel = y;
484  cnr = j;
485  stab = -x;
486  }
487  }
488  }
489  else
490  {
491  MSG_DEBUG( spxout << "DHARRI02 removing value " << cupd[i] << std::endl; )
492  cupd.clearNum(j);
493  }
494  }
495 
496  if (lastshift == solver()->shift())
497  {
498  if (cnr >= 0)
499  {
500  if (solver()->isBasic(enterId))
501  {
502  cupd.clearNum(cnr);
503  continue;
504  }
505  else
506  break;
507  }
508  else if (pnr >= 0)
509  {
510  pvec[pnr] = solver()->vector(pnr) * cvec;
511  if (solver()->isBasic(enterId))
512  {
513  pupd.setValue(pnr, 0.0);
514  continue;
515  }
516  else
517  {
518  x = pupd[pnr];
519  if (x > 0)
520  {
521  sel = upb[pnr] - pvec[pnr];
522  if (x < minStability && sel < delta)
523  {
524  minStability /= 2.0;
525  solver()->shiftUPbound(pnr, pvec[pnr]);
526  continue;
527  }
528  }
529  else
530  {
531  sel = lpb[pnr] - pvec[pnr];
532  if (-x < minStability && -sel < delta)
533  {
534  minStability /= 2.0;
535  solver()->shiftLPbound(pnr, pvec[pnr]);
536  continue;
537  }
538  }
539  sel /= x;
540  }
541  }
542  else
543  {
544  val = 0;
545  enterId.inValidate();
546  return enterId;
547  }
548 
549  if (sel > max) // instability detected => recompute
550  continue; // ratio test with corrected value
551  break;
552  }
553  }
554  }
555  else if (val < -epsilon)
556  {
557  for(;;)
558  {
559  pnr = -1;
560  cnr = -1;
561  max = val;
562  sel = val;
563  lastshift = solver()->shift();
564  assert(delta > epsilon);
565 
566 
567  // phase 1:
568  minDelta
569  (
570  &rmaxabs, /* max abs value in upd */
571  &max, /* initial and chosen value */
572  pupd.size(), /* # of indices in pupd */
573  pupd.indexMem(), /* nonzero indices in pupd */
574  pupd.values(), /* update vector for vec */
575  pvec.get_const_ptr(), /* current vector */
576  lpb.get_const_ptr(), /* lower bounds for vec */
577  upb.get_const_ptr(), /* upper bounds for vec */
578  epsilon); /* what is 0? */
579 
580  minDelta
581  (
582  &cmaxabs, /* max abs value in upd */
583  &max, /* initial and chosen value */
584  cupd.size(), /* # of indices in cupd */
585  cupd.indexMem(), /* nonzero indices in cupd */
586  cupd.values(), /* update vector for vec */
587  cvec.get_const_ptr(), /* current vector */
588  lcb.get_const_ptr(), /* lower bounds for vec */
589  ucb.get_const_ptr(), /* upper bounds for vec */
590  epsilon); /* what is 0? */
591 
592  if (max == val)
593  return enterId;
594 
595 
596  // phase 2:
597  stab = 0;
598  sel = infinity;
599  ruseeps = rmaxabs * epsilon * 0.001;
600  cuseeps = cmaxabs * epsilon * 0.001;
601  for (j = pupd.size() - 1; j >= 0; --j)
602  {
603  i = pupd.index(j);
604  x = pupd[i];
605  if (x > ruseeps)
606  {
607  y = lpb[i] - pvec[i];
608  if (y > degeneps)
609  solver()->shiftLPbound(i, pvec[i]); // ensure simplex improvement
610  else
611  {
612  y /= x;
613  if (y >= max && x > stab) // && y < sel+epsilon
614  {
615  enterId = solver()->id(i);
616  sel = y;
617  pnr = i;
618  stab = x;
619  }
620  }
621  }
622  else if (x < -ruseeps)
623  {
624  y = upb[i] - pvec[i];
625  if (y < -degeneps)
626  solver()->shiftUPbound(i, pvec[i]); // ensure simplex improvement
627  else
628  {
629  y /= x;
630  if (y >= max && -x > stab) // && y < sel+epsilon
631  {
632  enterId = solver()->id(i);
633  sel = y;
634  pnr = i;
635  stab = -x;
636  }
637  }
638  }
639  else
640  {
641  MSG_DEBUG( spxout << "DHARRI03 removing value " << pupd[i] << std::endl; )
642  pupd.clearNum(j);
643  }
644  }
645  for (j = cupd.size() - 1; j >= 0; --j)
646  {
647  i = cupd.index(j);
648  x = cupd[i];
649  if (x > cuseeps)
650  {
651  y = lcb[i] - cvec[i];
652  if (y > degeneps)
653  solver()->shiftLCbound(i, cvec[i]); // ensure simplex improvement
654  else
655  {
656  y /= x;
657  if (y >= max && x > stab) // && y < sel+epsilon
658  {
659  enterId = solver()->coId(i);
660  sel = y;
661  cnr = j;
662  stab = x;
663  }
664  }
665  }
666  else if (x < -cuseeps)
667  {
668  y = ucb[i] - cvec[i];
669  if (y < -degeneps)
670  solver()->shiftUCbound(i, cvec[i]); // ensure simplex improvement
671  else
672  {
673  y /= x;
674  if (y >= max && -x > stab) // && y < sel+epsilon
675  {
676  enterId = solver()->coId(i);
677  sel = y;
678  cnr = j;
679  stab = -x;
680  }
681  }
682  }
683  else
684  {
685  MSG_DEBUG( spxout << "DHARRI04 removing value " << x << std::endl; );
686  cupd.clearNum(j);
687  }
688  }
689 
690  if (lastshift == solver()->shift())
691  {
692  if (cnr >= 0)
693  {
694  if (solver()->isBasic(enterId))
695  {
696  cupd.clearNum(cnr);
697  continue;
698  }
699  else
700  break;
701  }
702  else if (pnr >= 0)
703  {
704  pvec[pnr] = solver()->vector(pnr) * cvec;
705  if (solver()->isBasic(enterId))
706  {
707  pupd.setValue(pnr, 0.0);
708  continue;
709  }
710  else
711  {
712  x = pupd[pnr];
713  if (x > 0)
714  {
715  sel = lpb[pnr] - pvec[pnr];
716  if (x < minStability && -sel < delta)
717  {
718  minStability /= 2;
719  solver()->shiftLPbound(pnr, pvec[pnr]);
720  continue;
721  }
722  }
723  else
724  {
725  sel = upb[pnr] - pvec[pnr];
726  if (-x < minStability && sel < delta)
727  {
728  minStability /= 2;
729  solver()->shiftUPbound(pnr, pvec[pnr]);
730  continue;
731  }
732  }
733  sel /= x;
734  }
735  }
736  else
737  {
738  val = 0;
739  enterId.inValidate();
740  return enterId;
741  }
742  if (sel < max) // instability detected => recompute
743  continue; // ratio test with corrected value
744  break;
745  }
746  }
747  }
748  assert(max * val >= 0);
749  assert(enterId.type() != SPxId::INVALID);
750 
751  val = sel;
752 
753  return enterId;
754 }
755 } // namespace soplex
756 
757 //-----------------------------------------------------------------------------
758 //Emacs Local Variables:
759 //Emacs mode:c++
760 //Emacs c-basic-offset:3
761 //Emacs tab-width:8
762 //Emacs indent-tabs-mode:nil
763 //Emacs End:
764 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen