SoPlex Doxygen Documentation
spxsolve.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
17 
18 #include <assert.h>
19 #include <iostream>
20 
21 #include "spxdefines.h"
22 #include "mpqreal.h"
23 #include "spxsolver.h"
24 #include "spxpricer.h"
25 #include "spxratiotester.h"
26 #include "spxdefaultrt.h"
27 #include "spxstarter.h"
28 #include "spxout.h"
29 #include "exceptions.h"
30 
31 #define MAXCYCLES 400
32 #define MAXSTALLS 10000
33 #define MAXSTALLRECOVERS 10
34 #define MAXREFACPIVOTS 10
35 
36 namespace soplex
37 {
38 
39 /// Interval for displaying iteration information.
40 long iterationInterval = 100;
41 
42 /**@todo check separately for ENTER and LEAVE algorithm */
43 bool SPxSolver::precisionReached(Real& newpricertol) const
44 {
45  Real maxViolRedCost;
46  Real sumViolRedCost;
47  Real maxViolBounds;
48  Real sumViolBounds;
49  Real maxViolConst;
50  Real sumViolConst;
51 
52  qualRedCostViolation(maxViolRedCost, sumViolRedCost);
53  qualBoundViolation(maxViolBounds, sumViolBounds);
54  qualConstraintViolation(maxViolConst, sumViolConst);
55 
56  // is the solution good enough ?
57  bool reached = maxViolRedCost < opttol() && maxViolBounds < feastol() && maxViolConst < feastol();
58 
59  if (!reached)
60  {
61  newpricertol = thepricer->epsilon() / 10.0;
62 
63  MSG_INFO3( spxout << "ISOLVE71 "
64  << "Precision not reached: Pricer tolerance = "
65  << thepricer->epsilon()
66  << " new tolerance = " << newpricertol
67  << std::endl
68  << " maxViolRedCost= " << maxViolRedCost
69  << " maxViolBounds= " << maxViolBounds
70  << " maxViolConst= " << maxViolConst
71  << std::endl
72  << " sumViolRedCost= " << sumViolRedCost
73  << " sumViolBounds= " << sumViolBounds
74  << " sumViolConst= " << sumViolConst
75  << std::endl; );
76  }
77  return reached;
78 }
79 
80 SPxSolver::Status SPxSolver::fpsolve()
81 {
82  METHOD( "SPxSolver::fpsolve()" );
83 
84  SPxId enterId;
85  int leaveNum;
86  int loopCount = 0;
87  Real minShift = infinity;
88  int cycleCount = 0;
89  bool priced = false;
90  Real lastDelta = 1;
91 
92  /* store the last (primal or dual) feasible objective value to recover/abort in case of stalling */
93  Real stallRefValue;
94  Real stallRefShift;
95  int stallRefIter;
96  int stallNumRecovers;
97 
98  if (dim() <= 0 && coDim() <= 0) // no problem loaded
99  {
100  m_status = NO_PROBLEM;
101  throw SPxStatusException("XSOLVE01 No Problem loaded");
102  }
103 
104  if (slinSolver() == 0) // linear system solver is required.
105  {
106  m_status = NO_SOLVER;
107  throw SPxStatusException("XSOLVE02 No Solver loaded");
108  }
109  if (thepricer == 0) // pricer is required.
110  {
111  m_status = NO_PRICER;
112  throw SPxStatusException("XSOLVE03 No Pricer loaded");
113  }
114  if (theratiotester == 0) // ratiotester is required.
115  {
116  m_status = NO_RATIOTESTER;
117  throw SPxStatusException("XSOLVE04 No RatioTester loaded");
118  }
119 
120  m_numCycle = 0;
121  iterCount = 0;
122  if (!isInitialized())
123  {
124  /*
125  if(SPxBasis::status() <= NO_PROBLEM)
126  SPxBasis::load(this);
127  */
128  /**@todo != REGULAR is not enough. Also OPTIMAL/DUAL/PRIMAL should
129  * be tested and acted accordingly.
130  */
131  if (thestarter != 0 && status() != REGULAR) // no basis and no starter.
132  thestarter->generate(*this); // generate start basis.
133 
134  init();
135 
136  // Inna/Tobi: init might fail, if the basis is singular
137  if( !isInitialized() )
138  {
139  assert(SPxBasis::status() == SPxBasis::SINGULAR);
140  m_status = UNKNOWN;
141  return status();
142  }
143  }
144 
145  //setType(type());
146 
147  if (!matrixIsSetup)
148  SPxBasis::load(this);
149 
150  //factorized = false;
151 
152  assert(thepricer->solver() == this);
153  assert(theratiotester->solver() == this);
154 
155  // maybe this should be done in init() ?
156  thepricer->setType(type());
157  theratiotester->setType(type());
158 
159  MSG_INFO3(
160  spxout << "ISOLVE72 starting value = " << value() << std::endl;
161  spxout << "ISOLVE73 starting shift = " << shift() << std::endl;
162  )
163 
164  if (SPxBasis::status() == SPxBasis::OPTIMAL)
165  setBasisStatus(SPxBasis::REGULAR);
166 
167  m_status = RUNNING;
168  bool stop = terminate();
169  leaveCount = 0;
170  enterCount = 0;
171  boundflips = 0;
172  totalboundflips = 0;
173 
174  stallNumRecovers = 0;
175 
176  /* if we run into a singular basis, we will retry from regulardesc with tighter tolerance in the ratio test */
177  SPxSolver::Type tightenedtype = type();
178  bool tightened = false;
179 
180  while (!stop)
181  {
182  const SPxBasis::Desc regulardesc = desc();
183 
184  // we need to reset these pointers to avoid unnecessary/wrong solves in leave() or enter()
185  solveVector2 = 0;
186  solveVector3 = 0;
187  coSolveVector2 = 0;
188 
189  try
190  {
191 
192  if (type() == ENTER)
193  {
194  int enterCycleCount = 0;
195  int enterFacPivotCount = 0;
196 
197  stallRefIter = iteration()-1;
198  stallRefShift = shift();
199  stallRefValue = value();
200 
201  /* in the entering algorithm, entertol() should be maintained by the ratio test and leavetol() should be
202  * reached by the pricer
203  */
204  Real maxpricertol = leavetol();
205  Real minpricertol = 0.01 * maxpricertol;
206 
207  thepricer->setEpsilon(maxpricertol);
208  priced = false;
209 
210  // to avoid shifts we restrict tolerances in the ratio test
211  if( loopCount > 0 )
212  {
213  lastDelta = (lastDelta < entertol()) ? lastDelta : entertol();
214  lastDelta *= 0.01;
215  theratiotester->setDelta(lastDelta);
216  assert(theratiotester->getDelta() > 0);
217  MSG_DEBUG( spxout << "decreased delta for ratiotest to: " << theratiotester->getDelta() << std::endl; )
218  }
219  else
220  {
221  lastDelta = 1;
222  theratiotester->setDelta(entertol());
223  }
224 
225  do
226  {
227  MSG_INFO3(
228  if( iteration() % iterationInterval == 0 )
229  spxout << "ISOLVE74 Enter iteration: " << iteration()
230  << ", Value = " << value()
231  << ", Shift = " << shift() << std::endl;
232  )
233  enterId = thepricer->selectEnter();
234 
235  if (!enterId.isValid())
236  {
237  // we are not infeasible and have no shift
238  if ( shift() <= epsilon()
239  && ( SPxBasis::status() == SPxBasis::REGULAR
240  || SPxBasis::status() == SPxBasis::DUAL
241  || SPxBasis::status() == SPxBasis::PRIMAL))
242  {
243  Real newpricertol = minpricertol;
244 
245  // is the solution good enough ?
246  // max three times reduced
247  if ((thepricer->epsilon() > minpricertol) && !precisionReached(newpricertol))
248  { // no!
249  // we reduce the pricer tolerance. Note that if the pricer does not find a candiate
250  // with the reduced tolerance, we quit, regardless of the violations.
251  if (newpricertol < minpricertol)
252  newpricertol = minpricertol;
253 
254  thepricer->setEpsilon(newpricertol);
255 
256  MSG_INFO2( spxout << "ISOLVE75 Setting pricer tolerance = "
257  << thepricer->epsilon()
258  << std::endl; )
259  }
260  // solution seems good, no check whether we are precise enough
261  else if (lastUpdate() == 0)
262  {
263  priced = true;
264  break;
265  }
266  // We have an iterationlimit and everything looks good? Then stop!
267  // 6 is just a number picked.
268  else if (maxIters > 0 && lastUpdate() < 6)
269  {
270  priced = true;
271  break;
272  }
273  }
274  MSG_INFO3( spxout << "ISOLVE76 solve(enter) triggers refactorization" << std::endl; )
275 
276  // if the factorization is not fresh, we better refactorize and call the pricer again; however, this can
277  // create cycling, so it is performed only a limited number of times per ENTER round
278  if( lastUpdate() > 0 && enterFacPivotCount < MAXREFACPIVOTS )
279  {
280  factorize();
281 
282  // if the factorization was found out to be singular, we have to quit
283  if( SPxBasis::status() < SPxBasis::REGULAR )
284  {
285  MSG_ERROR( spxout << "ESOLVE09 something wrong with factorization, Basis status: " << SPxBasis::status() << std::endl; )
286  stop = true;
287  break;
288  }
289 
290  // call pricer again
291  enterId = thepricer->selectEnter();
292 
293  // count how often the pricer has found something only after refactorizing
294  if( enterId.isValid() )
295  enterFacPivotCount++;
296  }
297 
298  if( !enterId.isValid() )
299  {
300  priced = true;
301  break;
302  }
303  }
304 
305  /* check if we have iterations left */
306  if (maxIters >= 0 && iterations() >= maxIters)
307  {
308  MSG_INFO2( spxout << "ISOLVE53e Maximum number of iterations (" << maxIters
309  << ") reached" << std::endl; )
310  m_status = ABORT_ITER;
311  stop = true;
312  break;
313  }
314 
315  enter(enterId);
316  assert((testBounds(), 1));
317  thepricer->entered4(lastEntered(), lastIndex());
318  stop = terminate();
319  clearUpdateVecs();
320  if( lastEntered().isValid() ) /* either a successful pivot was performed or a nonbasic variable flipped to the other bound */
321  {
322  enterCount++;
323  enterCycleCount = 0;
324  }
325  else
326  {
327  enterCycleCount++;
328  if( enterCycleCount > MAXCYCLES )
329  {
330  MSG_INFO2( spxout << "ISOLVE77 Abort solving due to cycling in "
331  << "entering algorithm" << std::endl; );
332  m_status = ABORT_CYCLING;
333  stop = true;
334  }
335  }
336 
337  /* check every MAXSTALLS iterations whether shift and objective value have not changed */
338  if( (iteration() - stallRefIter) % MAXSTALLS == 0 )
339  {
340  if( fabs(value() - stallRefValue) <= epsilon() && fabs(shift() - stallRefShift) <= epsilon() )
341  {
342  if( stallNumRecovers < MAXSTALLRECOVERS )
343  {
344  /* try to recover by unshifting/switching algorithm up to MAXSTALLRECOVERS times (just a number picked) */
345  MSG_INFO3( spxout << "ISOLVE21 Stalling detected - trying to recover by switching to LEAVING algorithm." << std::endl; )
346 
347  ++stallNumRecovers;
348  break;
349  }
350  else
351  {
352  /* giving up */
353  MSG_INFO2( spxout << "ISOLVE22 Abort solving due to stalling in entering algorithm." << std::endl; );
354 
355  m_status = ABORT_CYCLING;
356  stop = true;
357  }
358  }
359  else
360  {
361  /* merely update reference values */
362  stallRefIter = iteration()-1;
363  stallRefShift = shift();
364  stallRefValue = value();
365  }
366  }
367 
368  //@ assert(isConsistent());
369  }
370  while (!stop);
371 
372  MSG_INFO3(
373  spxout << "ISOLVE78 Enter finished. iteration: " << iteration()
374  << ", value: " << value()
375  << ", shift: " << shift()
376  << ", epsilon: " << epsilon()
377  << ", feastol: " << feastol()
378  << ", opttol: " << opttol()
379  << std::endl
380  << "ISOLVE56 stop: " << stop
381  << ", basis status: " << SPxBasis::status() << " (" << int(SPxBasis::status()) << ")"
382  << ", solver status: " << m_status << " (" << int(m_status) << ")" << std::endl;
383  )
384 
385  if (!stop)
386  {
387  /**@todo technically it would be ok to finish already when (priced && maxinfeas + shift() <= entertol()) is
388  * satisfied; maybe at least in the case when SoPlex keeps jumping back between ENTER and LEAVE always
389  * shifting (looping), we may relax this condition here;
390  * note also that unShift may increase shift() slightly due to roundoff errors
391  */
392  if (shift() <= epsilon())
393  {
394  // factorize();
395  unShift();
396 
397  Real maxinfeas = maxInfeas();
398 
399  MSG_INFO3(
400  spxout << "ISOLVE79 maxInfeas: " << maxinfeas
401  << ", shift: " << shift()
402  << ", entertol: " << entertol() << std::endl;
403  )
404 
405  if (priced && maxinfeas + shift() <= entertol())
406  {
407  setBasisStatus(SPxBasis::OPTIMAL);
408  m_status = OPTIMAL;
409  break;
410  }
411  }
412  setType(LEAVE);
413  init();
414  thepricer->setType(type());
415  theratiotester->setType(type());
416  }
417  }
418  else
419  {
420  assert(type() == LEAVE);
421 
422  int leaveCycleCount = 0;
423  int leaveFacPivotCount = 0;
424 
425  instableLeaveNum = -1;
426  instableLeave = false;
427 
428  stallRefIter = iteration()-1;
429  stallRefShift = shift();
430  stallRefValue = value();
431 
432  /* in the leaving algorithm, leavetol() should be maintained by the ratio test and entertol() should be reached
433  * by the pricer
434  */
435  Real maxpricertol = entertol();
436  Real minpricertol = 0.01 * maxpricertol;
437 
438  thepricer->setEpsilon(maxpricertol);
439  priced = false;
440 
441  // to avoid shifts we restrict tolerances in the ratio test
442  if( loopCount > 0 )
443  {
444  lastDelta = (lastDelta < leavetol()) ? lastDelta : leavetol();
445  lastDelta *= 0.01;
446  theratiotester->setDelta(lastDelta);
447  assert(theratiotester->getDelta() > 0);
448  MSG_DEBUG( spxout << "decreased delta for ratiotest to: " << theratiotester->getDelta() << std::endl; )
449  }
450  else
451  {
452  lastDelta = 1;
453  theratiotester->setDelta(leavetol());
454  }
455 
456  do
457  {
458  MSG_INFO3(
459  if( iteration() % iterationInterval == 0 )
460  spxout << "ISOLVE80 Leave Iteration: " << iteration()
461  << ", Value = " << value()
462  << ", Shift = " << shift() << std::endl;
463  )
464 
465  leaveNum = thepricer->selectLeave();
466 
467  if (leaveNum < 0 && instableLeaveNum >= 0 && lastUpdate() == 0)
468  {
469  /* no leaving variable was found, but because of instableLeaveNum >= 0 we know
470  that this is due to the scaling of theCoTest[...]. Thus, we use
471  instableLeaveNum and SPxFastRT::selectEnter shall accept even an instable
472  entering variable. */
473  MSG_INFO3(
474  spxout << "ISOLVE98 Trying instable leave iteration" << std::endl;
475  )
476 
477  leaveNum = instableLeaveNum;
478  instableLeave = true;
479  }
480  else
481  {
482  instableLeave = false;
483  }
484 
485  if (leaveNum < 0)
486  {
487  // we are not infeasible and have no shift
488  if ( shift() <= epsilon()
489  && ( SPxBasis::status() == SPxBasis::REGULAR
490  || SPxBasis::status() == SPxBasis::DUAL
491  || SPxBasis::status() == SPxBasis::PRIMAL))
492  {
493  Real newpricertol = minpricertol;
494 
495  // is the solution good enough ?
496  // max three times reduced
497  if ((thepricer->epsilon() > minpricertol) && !precisionReached(newpricertol))
498  { // no
499  // we reduce the pricer tolerance. Note that if the pricer does not find a candiate
500  // with the reduced pricer tolerance, we quit, regardless of the violations.
501  if (newpricertol < minpricertol)
502  newpricertol = minpricertol;
503 
504  thepricer->setEpsilon(newpricertol);
505 
506  MSG_INFO2( spxout << "ISOLVE81 Setting pricer tolerance = "
507  << thepricer->epsilon()
508  << std::endl; );
509  }
510  // solution seems good, no check whether we are precise enough
511  else if (lastUpdate() == 0)
512  {
513  priced = true;
514  break;
515  }
516  // We have an iteration limit and everything looks good? Then stop!
517  // 6 is just a number picked.
518  else if (maxIters > 0 && lastUpdate() < 6)
519  {
520  priced = true;
521  break;
522  }
523  }
524  MSG_INFO3( spxout << "ISOLVE82 solve(leave) triggers refactorization" << std::endl; )
525 
526  // if the factorization is not fresh, we better refactorize and call the pricer again; however, this can
527  // create cycling, so it is performed only a limited number of times per LEAVE round
528  if( lastUpdate() > 0 && leaveFacPivotCount < MAXREFACPIVOTS )
529  {
530  factorize();
531 
532  // Inna/Tobi: if the factorization was found out to be singular, we have to quit
533  if (SPxBasis::status() < SPxBasis::REGULAR)
534  {
535  MSG_ERROR( spxout << "ESOLVE10 something wrong with factorization, Basis status: " << SPxBasis::status() << std::endl; )
536  stop = true;
537  break;
538  }
539 
540  // call pricer again
541  leaveNum = thepricer->selectLeave();
542 
543  // count how often the pricer has found something only after refactorizing
544  if( leaveNum >= 0 )
545  leaveFacPivotCount++;
546  }
547 
548  if (leaveNum < 0)
549  {
550  priced = true;
551  break;
552  }
553  }
554 
555  /* check if we have iterations left */
556  if (maxIters >= 0 && iterations() >= maxIters)
557  {
558  MSG_INFO2( spxout << "ISOLVE53l Maximum number of iterations (" << maxIters
559  << ") reached" << std::endl; )
560  m_status = ABORT_ITER;
561  stop = true;
562  break;
563  }
564 
565  leave(leaveNum);
566  assert((testBounds(), 1));
567  thepricer->left4(lastIndex(), lastLeft());
568  stop = terminate();
569  clearUpdateVecs();
570  if( lastIndex() >= 0 ) /* either a successful pivot was performed or a nonbasic variable flipped to the other bound */
571  {
572  leaveCount++;
573  leaveCycleCount = 0;
574  }
575  else
576  {
577  leaveCycleCount++;
578  if( leaveCycleCount > MAXCYCLES )
579  {
580  MSG_INFO2( spxout << "ISOLVE83 Abort solving due to cycling in leaving algorithm" << std::endl; );
581  m_status = ABORT_CYCLING;
582  stop = true;
583  }
584  }
585 
586  /* check every MAXSTALLS iterations whether shift and objective value have not changed */
587  if( (iteration() - stallRefIter) % MAXSTALLS == 0 )
588  {
589  if( fabs(value() - stallRefValue) <= epsilon() && fabs(shift() - stallRefShift) <= epsilon() )
590  {
591  if( stallNumRecovers < MAXSTALLRECOVERS )
592  {
593  /* try to recover by switching algorithm up to MAXSTALLRECOVERS times */
594  MSG_INFO3( spxout << "ISOLVE24 Stalling detected - trying to recover by switching to ENTERING algorithm." << std::endl; )
595 
596  ++stallNumRecovers;
597  break;
598  }
599  else
600  {
601  /* giving up */
602  MSG_INFO2( spxout << "ISOLVE25 Abort solving due to stalling in leaving algorithm" << std::endl; );
603 
604  m_status = ABORT_CYCLING;
605  stop = true;
606  }
607  }
608  else
609  {
610  /* merely update reference values */
611  stallRefIter = iteration()-1;
612  stallRefShift = shift();
613  stallRefValue = value();
614  }
615  }
616 
617  //@ assert(isConsistent());
618  }
619  while (!stop);
620 
621  MSG_INFO3(
622  spxout << "ISOLVE84 Leave finished. iteration: " << iteration()
623  << ", value: " << value()
624  << ", shift: " << shift()
625  << ", epsilon: " << epsilon()
626  << ", feastol: " << feastol()
627  << ", opttol: " << opttol()
628  << std::endl
629  << "ISOLVE57 stop: " << stop
630  << ", basis status: " << SPxBasis::status() << " (" << int(SPxBasis::status()) << ")"
631  << ", solver status: " << m_status << " (" << int(m_status) << ")" << std::endl;
632  )
633 
634  if (!stop)
635  {
636  if( shift() < minShift )
637  {
638  minShift = shift();
639  cycleCount = 0;
640  }
641  else
642  {
643  cycleCount++;
644  if( cycleCount > MAXCYCLES )
645  {
646  m_status = ABORT_CYCLING;
647  throw SPxStatusException("XSOLVE13 Abort solving due to cycling");
648  }
649  MSG_INFO3(
650  spxout << "ISOLVE86 maxInfeas: " << maxInfeas()
651  << ", shift: " << shift()
652  << ", leavetol: " << leavetol()
653  << ", cycle count: " << cycleCount << std::endl;
654  )
655  }
656 
657  /**@todo technically it would be ok to finish already when (priced && maxinfeas + shift() <= entertol()) is
658  * satisfied; maybe at least in the case when SoPlex keeps jumping back between ENTER and LEAVE always
659  * shifting (looping), we may relax this condition here;
660  * note also that unShift may increase shift() slightly due to roundoff errors
661  */
662  if (shift() <= epsilon())
663  {
664  cycleCount = 0;
665  // factorize();
666  unShift();
667 
668  Real maxinfeas = maxInfeas();
669 
670  MSG_INFO3(
671  spxout << "ISOLVE87 maxInfeas: " << maxinfeas
672  << ", shift: " << shift()
673  << ", leavetol: " << leavetol() << std::endl;
674  )
675 
676  // We stop if we are indeed optimal, or if we have already been
677  // two times at this place. In this case it seems futile to
678  // continue.
679  if (loopCount > 2)
680  {
681  m_status = ABORT_CYCLING;
682  throw SPxStatusException("XSOLVE14 Abort solving due to looping");
683  }
684  else if (priced && maxinfeas + shift() <= leavetol())
685  {
686  setBasisStatus(SPxBasis::OPTIMAL);
687  m_status = OPTIMAL;
688  break;
689  }
690  loopCount++;
691  }
692  setType(ENTER);
693  init();
694  thepricer->setType(type());
695  theratiotester->setType(type());
696  }
697  }
698  assert(m_status != SINGULAR);
699 
700  }
701  catch( SPxException E )
702  {
703  // if we stopped due to a singular basis, we reload the original basis and try again with tighter
704  // tolerance (only once)
705  if (m_status == SINGULAR && !tightened)
706  {
707  tightenedtype = type();
708 
709  if( tightenedtype == ENTER )
710  {
711  m_entertol = 0.01 * m_entertol;
712 
713  MSG_INFO2( spxout << "ISOLVE26e basis singular: reloading basis and solving with tighter ratio test tolerance " << m_entertol << std::endl; )
714  }
715  else
716  {
717  m_leavetol = 0.01 * m_leavetol;
718 
719  MSG_INFO2( spxout << "ISOLVE26l basis singular: reloading basis and solving with tighter ratio test tolerance " << m_leavetol << std::endl; )
720  }
721 
722  // load original basis
723  int niters = iterations();
724  loadBasis(regulardesc);
725 
726  // remember iteration count
727  iterCount = niters;
728 
729  // try initializing basis (might fail if starting basis was already singular)
730  try
731  {
732  init();
733  theratiotester->setType(type());
734  }
735  catch( SPxException Ex )
736  {
737  MSG_INFO2( spxout << "ISOLVE27 reloaded basis singular, resetting original tolerances" << std::endl; )
738 
739  if( tightenedtype == ENTER )
740  m_entertol = 100.0 * m_entertol;
741  else
742  m_leavetol = 100.0 * m_leavetol;
743 
744  theratiotester->setType(type());
745 
746  throw Ex;
747  }
748 
749  // reset status and counters
750  m_status = RUNNING;
751  m_numCycle = 0;
752  leaveCount = 0;
753  enterCount = 0;
754  stallNumRecovers = 0;
755 
756  // continue
757  stop = false;
758  tightened = true;
759  }
760  // reset tolerance to its original value and pass on the exception
761  else if (tightened)
762  {
763  if( tightenedtype == ENTER )
764  m_entertol = 100.0 * m_entertol;
765  else
766  m_leavetol = 100.0 * m_leavetol;
767 
768  theratiotester->setType(type());
769 
770  throw E;
771  }
772  // pass on the exception
773  else
774  throw E;
775  }
776  }
777 
778  // reset tolerance to its original value
779  if (tightened)
780  {
781  if( tightenedtype == ENTER )
782  m_entertol = 100.0 * m_entertol;
783  else
784  m_leavetol = 100.0 * m_leavetol;
785 
786  theratiotester->setType(type());
787  }
788 
789  if (m_status == RUNNING)
790  {
791  m_status = ERROR;
792  throw SPxStatusException("XSOLVE05 Status is still RUNNING when it shouldn't be");
793  }
794 
795  MSG_INFO1(
796  spxout << "ISOLVE02 Finished solving (status=" << status()
797  << ", iters=" << iterCount
798  << ", leave=" << leaveCount
799  << ", enter=" << enterCount
800  << ", flips=" << totalboundflips;
801  if( status() == OPTIMAL )
802  spxout << ", objValue=" << value();
803  spxout << ")" << std::endl;
804  )
805 
806 #ifdef ENABLE_ADDITIONAL_CHECKS
807  /* check if solution is really feasible */
808  if( status() == OPTIMAL )
809  {
810  int c;
811  Real val;
812  DVector sol( nCols() );
813 
814  getPrimal( sol );
815 
816  for(int row = 0; row < nRows(); ++row )
817  {
818  const SVector& rowvec = rowVector( row );
819  val = 0.0;
820  for( c = 0; c < rowvec.size(); ++c )
821  val += rowvec.value( c ) * sol[rowvec.index( c )];
822 
823  if( LT( val, lhs( row ), feastol() ) ||
824  GT( val, rhs( row ), feastol() ) )
825  {
826  // Minor rhs violations happen frequently, so print these
827  // warnings only with verbose level INFO2 and higher.
828  MSG_INFO2( spxout << "WSOLVE88 Warning! Constraint " << row
829  << " is violated by solution" << std::endl
830  << " lhs:" << lhs( row )
831  << " <= val:" << val
832  << " <= rhs:" << rhs( row ) << std::endl; )
833 
834  if( type() == LEAVE && isRowBasic( row ) )
835  {
836  // find basis variable
837  for( c = 0; c < nRows(); ++c )
838  if (basis().baseId(c).isSPxRowId()
839  && (number(basis().baseId(c)) == row))
840  break;
841 
842  assert( c < nRows() );
843 
844  MSG_WARNING( spxout << "WSOLVE90 basis idx:" << c
845  << " fVec:" << fVec()[c]
846  << " fRhs:" << fRhs()[c]
847  << " fTest:" << fTest()[c] << std::endl; )
848  }
849  }
850  }
851  for(int col = 0; col < nCols(); ++col )
852  {
853  if( LT( sol[col], lower( col ), feastol() ) ||
854  GT( sol[col], upper( col ), feastol() ) )
855  {
856  // Minor bound violations happen frequently, so print these
857  // warnings only with verbose level INFO2 and higher.
858  MSG_INFO2( spxout << "WSOLVE91 Warning! Bound for column " << col
859  << " is violated by solution" << std::endl
860  << " lower:" << lower( col )
861  << " <= val:" << sol[col]
862  << " <= upper:" << upper( col ) << std::endl; )
863 
864  if( type() == LEAVE && isColBasic( col ) )
865  {
866  for( c = 0; c < nRows() ; ++c)
867  if ( basis().baseId( c ).isSPxColId()
868  && ( number( basis().baseId( c ) ) == col ))
869  break;
870 
871  assert( c < nRows() );
872  MSG_WARNING( spxout << "WSOLVE92 basis idx:" << c
873  << " fVec:" << fVec()[c]
874  << " fRhs:" << fRhs()[c]
875  << " fTest:" << fTest()[c] << std::endl; )
876  }
877  }
878  }
879  }
880 #endif // ENABLE_ADDITIONAL_CHECKS
881 
882  return status();
883 }
884 
885 SPxSolver::Status SPxSolver::solve()
886 {
887  METHOD( "SPxSolver::solve()" );
888 
889  theTime.reset();
890  theTime.start();
891 
892  /* if tolerances are not below irthreshold(), perform a standard floating point solve */
893  if( feastol() >= irthreshold() && opttol() >= irthreshold() )
894  {
895  fpsolve();
896 
897  theTime.stop();
898  theCumulativeTime += time();
899 
900  return status();
901  }
902 
903  /* perform iterative refinement only if exact arithmetic is available */
904  if( !MpqRealIsExact() )
905  {
906  MSG_WARNING( spxout << "WSOLVE35 Warning: Iterative refinement disabled because of missing GMP support (compile with GMP=true).\n" );
907 
908  fpsolve();
909 
910  theTime.stop();
911  theCumulativeTime += time();
912 
913  return status();
914  }
915 
916  /* otherwise apply iterative refinement */
917  SPxBasis::SPxStatus basisstat;
918  SPxBasis::Desc basisdesc;
919  LPColSet slackcols;
920  Real irobjlimit;
921  Real irfeastol;
922  Real iropttol;
923  int iteroffset;
924  bool refined;
925  bool precisionreached;
926 
927  /* refined solution vectors */
928  DVector_exact primal_ex;
929  DVector_exact slack_ex;
930  DVector_exact dual_ex;
931  DVector_exact redcost_ex;
932 
933  /* store target tolerances */
934  irfeastol = feastol();
935  iropttol = opttol();
936 
937  /* limit floating point tolerance */
938  if( feastol() < DEFAULT_BND_VIOL )
939  {
940  MSG_DEBUG( spxout << "relaxing feastol() to " << DEFAULT_BND_VIOL << "\n" );
941  setFeastol(DEFAULT_BND_VIOL);
942  }
943 
944  if( opttol() < DEFAULT_BND_VIOL )
945  {
946  MSG_DEBUG( spxout << "relaxing opttol() to " << DEFAULT_BND_VIOL << "\n" );
947  setOpttol(DEFAULT_BND_VIOL);
948  }
949 
950  /* deactivate objective limit; because fpsolve() uses a relaxed tolerance for checking dual feasibility, it might
951  * return ABORT_VALUE although dual feasibility is not guaranteed within opttol()
952  */
953  /**@todo make objLimit consistent in SPxSolver */
954  irobjlimit = objLimit;
955  if( objLimit < infinity )
956  {
957  MSG_DEBUG( spxout << "deactivating objective limit\n" );
958  objLimit = infinity;
959  }
960 
961  /* store basis */
962  MSG_DEBUG( spxout << "storing basis\n" );
963 
964  basisstat = basis().status();
965  basisdesc = basis().desc();
966 
967  assert(basisstat == SPxBasis::NO_PROBLEM || basis().isDescValid(basisdesc));
968 
969  /* add artificial slack variables to convert inequality to equality constraints */
970  for( int i = 0; i < nRows(); i++ )
971  {
972  if( lhs(i) != rhs(i) )
973  {
974  slackcols.add(0.0, -rhs(i), unitVecs[i], -lhs(i));
975  changeLhs(i, 0.0);
976  changeRhs(i, 0.0);
977  }
978  }
979 
980  MSG_DEBUG( spxout << "adding slack columns\n" );
981  addCols(slackcols);
982 
983  /* adjust basis */
984  if( basisstat != SPxBasis::NO_PROBLEM )
985  {
986  basisdesc.reSize(nRows(), nCols());
987  for( int i = 0; i < slackcols.num(); i++ )
988  {
989  int col;
990  int row;
991 
992  col = nCols() - slackcols.num() + i;
993  row = slackcols.colVector(i).index(0);
994  assert(row >= 0);
995  assert(row < nRows());
996 
997  switch( basisdesc.rowStatus(row) )
998  {
999  case SPxBasis::Desc::P_ON_LOWER:
1000  basisdesc.colStatus(col) = SPxBasis::Desc::P_ON_UPPER;
1001  break;
1002  case SPxBasis::Desc::P_ON_UPPER:
1003  basisdesc.colStatus(col) = SPxBasis::Desc::P_ON_LOWER;
1004  break;
1005  case SPxBasis::Desc::D_ON_UPPER:
1006  basisdesc.colStatus(col) = SPxBasis::Desc::D_ON_LOWER;
1007  break;
1008  case SPxBasis::Desc::D_ON_LOWER:
1009  basisdesc.colStatus(col) = SPxBasis::Desc::D_ON_UPPER;
1010  break;
1011  case SPxBasis::Desc::P_FREE:
1012  case SPxBasis::Desc::P_FIXED:
1013  case SPxBasis::Desc::D_FREE:
1014  case SPxBasis::Desc::D_ON_BOTH:
1015  case SPxBasis::Desc::D_UNDEFINED:
1016  default:
1017  basisdesc.colStatus(col) = basisdesc.rowStatus(row);
1018  break;
1019  }
1020 
1021  basisdesc.rowStatus(row) = SPxBasis::Desc::P_FIXED;
1022  }
1023 
1024  /* load adjusted basis */
1025  MSG_DEBUG( spxout << "loading adjusted basis\n" );
1026 
1027  loadBasis(basisdesc);
1028  }
1029 
1030  refined = false;
1031  precisionreached = false;
1032 
1033  try
1034  {
1035  /* perform initial floating point solve */
1036  MSG_INFO1( spxout << "\nstarting floating-point solve . . .\n\n" );
1037  fpsolve();
1038 
1039  /* count simplex iterations and update iteration limit */
1040  iteroffset = iterations();
1041 
1042  if( terminationIter() >= 0 )
1043  {
1044  assert(iterations() <= terminationIter());
1045  setTerminationIter(terminationIter() - iterations());
1046  }
1047 
1048  /* apply iterative refinement only if initial solve returned optimal */
1049  refined = (status() == OPTIMAL);
1050  if( refined )
1051  {
1052  /* create refined solution vectors */
1053  primal_ex = DVector_exact(nCols());
1054  slack_ex = DVector_exact(nRows());
1055  dual_ex = DVector_exact(nRows());
1056  redcost_ex = DVector_exact(nCols());
1057 
1058  precisionreached = refine(irfeastol, iropttol, primal_ex, slack_ex, dual_ex, redcost_ex, 5 * iterations());
1059 
1060  /* count simplex iterations and update iteration limit */
1061  iteroffset += iterations();
1062 
1063  if( terminationIter() >= 0 )
1064  {
1065  assert(iterations() <= terminationIter());
1066  setTerminationIter(terminationIter() - iterations());
1067  }
1068  }
1069  else
1070  {
1071  MSG_INFO1( spxout << "\nskipping refinement because of status " << status() << "\n\n" );
1072  }
1073  }
1074  catch( SPxException E )
1075  {
1076  MSG_INFO1( spxout << "iterative refinement threw exception: " << E.what() << "\n" );
1077  }
1078 
1079  /* adjust basis; restore lhs and rhs */
1080  assert(basis().status() != SPxBasis::NO_PROBLEM);
1081  basisdesc = basis().desc();
1082  for( int i = 0; i < slackcols.num(); i++ )
1083  {
1084  int col;
1085  int row;
1086 
1087  col = nCols() - slackcols.num() + i;
1088  row = slackcols.colVector(i).index(0);
1089  assert(row >= 0);
1090  assert(row < nRows());
1091  assert(basisdesc.rowStatus(row) == SPxBasis::Desc::P_FIXED || basisdesc.rowStatus(row) == SPxBasis::Desc::D_FREE);
1092 
1093  MSG_DEBUG( spxout << "removing slack column C" << col << ": status=" << basisdesc.colStatus(col) << ", lower=" << lower(col) << ", upper=" << upper(col) << "\n" );
1094 
1095  if( basisdesc.rowStatus(row) == SPxBasis::Desc::P_FIXED )
1096  {
1097  switch( basisdesc.colStatus(col) )
1098  {
1099  case SPxBasis::Desc::P_ON_LOWER:
1100  basisdesc.rowStatus(row) = SPxBasis::Desc::P_ON_UPPER;
1101  break;
1102  case SPxBasis::Desc::P_ON_UPPER:
1103  basisdesc.rowStatus(row) = SPxBasis::Desc::P_ON_LOWER;
1104  break;
1105  case SPxBasis::Desc::D_ON_UPPER:
1106  basisdesc.rowStatus(row) = SPxBasis::Desc::D_ON_LOWER;
1107  break;
1108  case SPxBasis::Desc::D_ON_LOWER:
1109  basisdesc.rowStatus(row) = SPxBasis::Desc::D_ON_UPPER;
1110  break;
1111  case SPxBasis::Desc::P_FREE:
1112  case SPxBasis::Desc::P_FIXED:
1113  case SPxBasis::Desc::D_FREE:
1114  case SPxBasis::Desc::D_ON_BOTH:
1115  case SPxBasis::Desc::D_UNDEFINED:
1116  default:
1117  basisdesc.rowStatus(row) = basisdesc.colStatus(col);
1118  break;
1119  }
1120  }
1121 
1122  changeLhs(row, -upper(col));
1123  changeRhs(row, -lower(col));
1124  }
1125 
1126  if( refined )
1127  {
1128  /* adjust slack values */
1129  for( int i = 0; i < slackcols.num(); i++ )
1130  {
1131  int col;
1132  int row;
1133 
1134  col = nCols() - slackcols.num() + i;
1135  row = slackcols.colVector(i).index(0);
1136  assert(row >= 0);
1137  assert(row < nRows());
1138 
1139  slack_ex[row] -= primal_ex[col];
1140  }
1141 
1142  /* resize primal and redcost vector */
1143  primal_ex.reDim(nCols() - slackcols.num());
1144  redcost_ex.reDim(nCols() - slackcols.num());
1145  }
1146 
1147  /* remove slack columns */
1148  MSG_DEBUG( spxout << "removing slack columns\n" );
1149  removeColRange(nCols() - slackcols.num(), nCols() - 1);
1150 
1151  try
1152  {
1153  /* load adjusted basis */
1154  MSG_DEBUG( spxout << "loading adjusted basis\n" );
1155 
1156  basisdesc.reSize(nRows(), nCols());
1157  loadBasis(basisdesc);
1158 
1159  /* initialize data structures after removing the slack columns; if refinement has been applied, we are at an
1160  * optimal basis; otherwise we have to perform a last floating point solve after removing slack columns
1161  */
1162  if( refined )
1163  init();
1164  else
1165  fpsolve();
1166  }
1167  catch( SPxException E )
1168  {
1169  MSG_INFO1( spxout << "iterative refinement threw exception: " << E.what() << "\n" );
1170  refined = false;
1171  }
1172 
1173  /* install refined solution */
1174  if( refined )
1175  {
1176  /**@todo use only one DVector as buffer */
1177  DVector primal_fp(nCols());
1178  DVector slack_fp(nRows());
1179  DVector dual_fp(nRows());
1180  DVector redcost_fp(nCols());
1181 
1182 #ifndef NDEBUG
1183  /**@todo make debug messages */
1184  DVector tmp;
1185 
1186  tmp = DVector(nCols());
1187  getPrimal(tmp);
1188  primal_fp = DVector(primal_ex);
1189  setPrimal(primal_fp);
1190  getPrimal(primal_fp);
1191  tmp -= primal_fp;
1192  if( tmp.length() > DEFAULT_BND_VIOL )
1193  {
1194  MSG_INFO3( spxout << "distance between floating point and refined primal solution is " << tmp.length() << std::endl );
1195  }
1196 
1197  tmp = DVector(nRows());
1198  getSlacks(tmp);
1199  slack_fp = DVector(slack_ex);
1200  setSlacks(slack_fp);
1201  getSlacks(slack_fp);
1202  tmp -= slack_fp;
1203  if( tmp.length() > DEFAULT_BND_VIOL )
1204  {
1205  MSG_INFO3( spxout << "distance between floating point and refined slack vector is " << tmp.length() << std::endl );
1206  }
1207 
1208  tmp = DVector(nRows());
1209  getDual(tmp);
1210  dual_fp = DVector(dual_ex);
1211  setDual(dual_fp);
1212  getDual(dual_fp);
1213  tmp -= dual_fp;
1214  if( tmp.length() > DEFAULT_BND_VIOL )
1215  {
1216  MSG_INFO3( spxout << "distance between floating point and refined dual solution is " << tmp.length() << std::endl );
1217  }
1218 
1219  tmp = DVector(nCols());
1220  getRedCost(tmp);
1221  redcost_fp = DVector(redcost_ex);
1222  setRedCost(redcost_fp);
1223  getRedCost(redcost_fp);
1224  tmp -= redcost_fp;
1225  if( tmp.length() > DEFAULT_BND_VIOL )
1226  {
1227  MSG_INFO3( spxout << "distance between floating point and refined redcost vector is " << tmp.length() << std::endl );
1228  }
1229 #else
1230  primal_fp = DVector(primal_ex);
1231  setPrimal(primal_fp);
1232 
1233  slack_fp = DVector(slack_ex);
1234  setSlacks(slack_fp);
1235 
1236  dual_fp = DVector(dual_ex);
1237  setDual(dual_fp);
1238 
1239  redcost_fp = DVector(redcost_ex);
1240  setRedCost(redcost_fp);
1241 #endif
1242  }
1243 
1244  /* restore tolerances; this must be done after init(), because the shift procedures can't handle tolerances smaller
1245  * than epsilon
1246  */
1247  MSG_DEBUG( spxout << "restoring tolerances\n" );
1248  setFeastol(irfeastol);
1249  setOpttol(iropttol);
1250 
1251  /* restore objective limit */
1252  if( objLimit != irobjlimit )
1253  {
1254  MSG_DEBUG( spxout << "restoring objective limit\n" );
1255  objLimit = irobjlimit;
1256  }
1257 
1258  /* restore number of simplex iterations and iteration limit */
1259  SPxBasis::iterCount = iteroffset;
1260  if( terminationIter() >= 0 )
1261  setTerminationIter(terminationIter() + iteroffset);
1262 
1263  /* set refinement has been applied, we are at an optimal basis */
1264  if( refined )
1265  {
1266  m_status = SPxSolver::OPTIMAL;
1267  setBasisStatus(SPxBasis::OPTIMAL);
1268  }
1269 
1270  /* print warning if iterative refinement terminated before reaching tolerances */
1271  if( !precisionreached )
1272  {
1273  MSG_INFO1( spxout << "\nwarning: iterative refinement could not reach final precision\n" );
1274  }
1275 
1276  MSG_DEBUG( spxout << "returning with status " << status() << "\n" );
1277 
1278  theTime.stop();
1279  theCumulativeTime += time();
1280 
1281  return status();
1282 }
1283 
1284 bool SPxSolver::refine(
1285  Real irfeastol, /**< primal feasibility tolerance */
1286  Real iropttol, /**< dual feasibility tolerance */
1287  Vector_exact& primal_ex, /**< buffer to return refined primal solution values */
1288  Vector_exact& slack_ex, /**< buffer to return refined slack values */
1289  Vector_exact& dual_ex, /**< buffer to return refined dual solution values */
1290  Vector_exact& redcost_ex, /**< buffer to return refined reduced cost values */
1291  int maxitersround /**< iteration limit per refinement round */
1292  )
1293 {
1294  METHOD( "SPxSolver::refine()" );
1295 
1296  assert(maxRefines == -1 || maxRefines >= 0);
1297  assert(irfeastol >= 0.0);
1298  assert(iropttol >= 0.0);
1299  assert(status() == OPTIMAL);
1300 
1301  /* original problem in floating point precision */
1302  DVector orilower_fp(nCols());
1303  DVector oriupper_fp(nCols());
1304  DVector orirhs_fp(nRows());
1305  DVector oriobj_fp(nCols());
1306 
1307  /* modified problem in floating point precision */
1308  DVector modlower_fp(nCols());
1309  DVector modupper_fp(nCols());
1310  DVector modrhs_fp(nRows());
1311  DVector modobj_fp(nCols());
1312 
1313  /* modified problem in exact precision */
1314  DVector_exact modlower_ex(nCols());
1315  DVector_exact modupper_ex(nCols());
1316  DVector_exact modrhs_ex(nRows());
1317  DVector_exact modobj_ex(nCols());
1318 
1319  /* solution to modified problem in exact precision */
1320  DVector_exact modprimal_ex(nCols());
1321  DVector_exact moddual_ex(nRows());
1322 
1323  /* solution in floating point precision needed for cast to exact precision */
1324  DVector primal_fp(nCols());
1325  DVector dual_fp(nRows());
1326 
1327  /* maximum violations */
1328  MpqReal boundsviol_ex;
1329  MpqReal sidesviol_ex;
1330  MpqReal redcostviol_ex;
1331 
1332  /* scaling factors */
1333  MpqReal primalscale_ex;
1334  MpqReal dualscale_ex;
1335 
1336  /* number of simplex iterations, refinements in total, and refinements with actual simplex iterations being performed */
1337  int iteroffset;
1338  int nrefines;
1339 #if defined(DEBUGGING)
1340  int npivotrefines;
1341 #endif
1342 
1343  /* return status */
1344  bool precisionreached = false;
1345 
1346  /* flag to mark whether problem has been modified */
1347  bool changed = false;
1348 
1349  /* store current (regular) basis */
1350  SPxBasis::Desc basisdesc = basis().desc();
1351  assert(basis().status() == SPxBasis::OPTIMAL);
1352  assert(basis().isDescValid(basisdesc));
1353 
1354  /* store original problem */
1355  orilower_fp = lower();
1356  oriupper_fp = upper();
1357  orirhs_fp = rhs();
1358  getObj(oriobj_fp);
1359 
1360  /* get floating point solution of original problem */
1361  getPrimal(primal_fp);
1362  getDual(dual_fp);
1363 
1364  /* store floating point solution in exact precision */
1365  primal_ex = primal_fp;
1366  dual_ex = dual_fp;
1367 
1368  /* initialize values */
1369  nrefines = 0;
1370  primalscale_ex = 1;
1371  dualscale_ex = 1;
1372 #if defined(DEBUGGING)
1373  npivotrefines = 0;
1374 #endif
1375 
1376  /* count iterations performed during refinement */
1377  iteroffset = 0;
1378 
1379  /* refinement loop */
1380  do
1381  {
1382  MpqReal maxscale_ex;
1383 
1384  assert(status() == OPTIMAL);
1385 
1386  MSG_DEBUG( spxout << "computing violations exactly . . .\n" );
1387 
1388  /* compute bound violations */
1389  boundsviol_ex = 0;
1390  for( int c = 0; c < nCols(); c++ )
1391  {
1392  modlower_ex[c] = orilower_fp[c];
1393  modlower_ex[c] -= primal_ex[c];
1394 
1395  modupper_ex[c] = oriupper_fp[c];
1396  modupper_ex[c] -= primal_ex[c];
1397 
1398  if( modlower_ex[c] > boundsviol_ex )
1399  boundsviol_ex = modlower_ex[c];
1400 
1401  if( modupper_ex[c] < -boundsviol_ex )
1402  boundsviol_ex = -modupper_ex[c];
1403  }
1404 
1405  /* compute sides violation */
1406  slack_ex = computePrimalActivity(primal_ex);
1407 
1408  sidesviol_ex = 0;
1409  for( int r = 0; r < nRows(); r++ )
1410  {
1411  assert(lhs(r) == rhs(r));
1412 
1413  modrhs_ex[r] = orirhs_fp[r];
1414  modrhs_ex[r] -= slack_ex[r];
1415 
1416  if( modrhs_ex[r] > sidesviol_ex )
1417  sidesviol_ex = modrhs_ex[r];
1418  else if( modrhs_ex[r] < -sidesviol_ex )
1419  sidesviol_ex = -modrhs_ex[r];
1420  }
1421  MSG_DEBUG( spxout << "\n" );
1422 
1423  /* compute reduced costs and reduced cost violation */
1424  redcost_ex = computeDualActivity(dual_ex);
1425  redcostviol_ex = 0;
1426 
1427  for( int c = 0; c < nCols(); c++ )
1428  {
1429  SPxBasis::Desc::Status basisstat = basisdesc.colStatus(c);
1430 
1431  redcost_ex[c] *= -1;
1432  redcost_ex[c] += oriobj_fp[c];
1433 
1434  modobj_ex[c] = redcost_ex[c];
1435 
1436  if( (spxSense() == MINIMIZE && basisstat != SPxBasis::Desc::P_ON_UPPER && basisstat != SPxBasis::Desc::P_FIXED)
1437  || (spxSense() == MAXIMIZE && basisstat != SPxBasis::Desc::P_ON_LOWER && basisstat != SPxBasis::Desc::P_FIXED) )
1438  {
1439  if( redcost_ex[c] < -redcostviol_ex )
1440  redcostviol_ex = -redcost_ex[c];
1441  }
1442 
1443  if( (spxSense() == MINIMIZE && basisstat != SPxBasis::Desc::P_ON_LOWER && basisstat != SPxBasis::Desc::P_FIXED)
1444  || (spxSense() == MAXIMIZE && basisstat != SPxBasis::Desc::P_ON_UPPER && basisstat != SPxBasis::Desc::P_FIXED) )
1445  {
1446  if( redcost_ex[c] > redcostviol_ex )
1447  redcostviol_ex = redcost_ex[c];
1448  }
1449  }
1450 
1451  /* at this point the vectors primal_ex, slack_ex, dual_ex, redcost_ex are computed and we may exit */
1452 
1453  /* terminate if tolerances are satisfied */
1454  if( boundsviol_ex <= irfeastol && sidesviol_ex <= irfeastol && redcostviol_ex <= iropttol )
1455  {
1456  MSG_INFO1( spxout << "\nrefinement finished: tolerances reached\n\n" );
1457  assert(status() == OPTIMAL);
1458  precisionreached = true;
1459  break;
1460  }
1461 
1462  /* terminate if maximum number of refinements is reached */
1463  if( maxRefines >= 0 && nrefines >= maxRefines )
1464  {
1465  MSG_INFO1( spxout << "\nrefinement finished: maximum number of refinement rounds reached\n\n" );
1466  assert(status() == OPTIMAL);
1467  break;
1468  }
1469 
1470  /* otherwise continue */
1471  MSG_INFO1( spxout << "\nstarting refinement round " << nrefines+1 << ": " );
1472 
1473  /* compute primal scaling factor; limit increase in scaling by tolerance used in floating point solve */
1474  maxscale_ex = primalscale_ex / feastol();
1475 
1476  primalscale_ex = boundsviol_ex > sidesviol_ex ? boundsviol_ex : sidesviol_ex;
1477  assert(primalscale_ex >= 0);
1478 
1479  if( primalscale_ex > 0 )
1480  {
1481  primalscale_ex = 1 / primalscale_ex;
1482  if( primalscale_ex > maxscale_ex )
1483  primalscale_ex = maxscale_ex;
1484  }
1485  else
1486  primalscale_ex = maxscale_ex;
1487 
1488  if( primalscale_ex < 1 )
1489  primalscale_ex = 1;
1490 
1491  MSG_INFO1( spxout << "scaling primal by " << primalscale_ex );
1492 
1493  /* compute dual scaling factor; limit increase in scaling by tolerance used in floating point solve */
1494  maxscale_ex = dualscale_ex / opttol();
1495 
1496  dualscale_ex = redcostviol_ex;
1497  assert(dualscale_ex >= 0);
1498 
1499  if( dualscale_ex > 0 )
1500  {
1501  dualscale_ex = 1 / dualscale_ex;
1502  if( dualscale_ex > maxscale_ex )
1503  dualscale_ex = maxscale_ex;
1504  }
1505  else
1506  dualscale_ex = maxscale_ex;
1507 
1508  if( dualscale_ex < 1 )
1509  dualscale_ex = 1;
1510 
1511  MSG_INFO1( spxout << ", scaling dual by " << dualscale_ex << " . . .\n\n" );
1512 
1513  /* perform primal scaling */
1514  modlower_ex *= primalscale_ex;
1515  modupper_ex *= primalscale_ex;
1516  modrhs_ex *= primalscale_ex;
1517 
1518  /* perform dual scaling */
1519  modobj_ex *= dualscale_ex;
1520 
1521  /* change bounds */
1522  modlower_fp = DVector(modlower_ex);
1523  modupper_fp = DVector(modupper_ex);
1524  changeBounds(modlower_fp, modupper_fp);
1525 
1526  /* change sides */
1527  modrhs_fp = DVector(modrhs_ex);
1528  changeRange(modrhs_fp, modrhs_fp);
1529 
1530  /* change objective function */
1531  modobj_fp = DVector(modobj_ex);
1532  changeObj(modobj_fp);
1533 
1534  /* mark problem as changed */
1535  changed = true;
1536 
1537  MSG_DEBUG( spxout << "solving modified problem . . .\n" );
1538 
1539  /* load basis */
1540  loadBasis(basisdesc);
1541 
1542  /* count refinement rounds */
1543  nrefines++;
1544 
1545  /* solve modified problem */
1546  try
1547  {
1548  int maxiters = terminationIter();
1549 
1550  /* limit number of iterations per refinement round */
1551  if( maxitersround >= 0 && (terminationIter() < 0 || terminationIter() > maxitersround) )
1552  setTerminationIter(maxitersround);
1553 
1554  fpsolve();
1555 
1556  setTerminationIter(maxiters);
1557  }
1558  catch( SPxException E )
1559  {
1560  MSG_DEBUG( spxout << "fpsolve() threw exception: " << E.what() << "\n" );
1561  }
1562 
1563 #if defined(DEBUGGING)
1564  /* remember whether we moved to a new basis */
1565  if( iterations() > 0 )
1566  npivotrefines = nrefines;
1567 #endif
1568 
1569  /* count simplex iterations and update global iteration limit */
1570  iteroffset += iterations();
1571 
1572  if( terminationIter() >= 0 )
1573  {
1574  assert(iterations() <= terminationIter());
1575  setTerminationIter(terminationIter() - iterations());
1576  }
1577 
1578  /* if modified problem was not solved to optimality, stop refinement */
1579  if( status() != OPTIMAL )
1580  {
1581  MSG_INFO1( spxout << "\nrefinement finished: modified problem terminated with status " << status() << "\n\n" );
1582 
1583  /* load last optimal basis */
1584  loadBasis(basisdesc);
1585 
1586  break;
1587  }
1588  /* otherwise continue and correct primal and dual solution */
1589  else
1590  {
1591  int nadjusted = 0;
1592 
1593  /* get basis */
1594  basisdesc = basis().desc();
1595  assert(basis().isDescValid(basisdesc));
1596 
1597  /* get floating point solution of modified problem */
1598  getPrimal(primal_fp);
1599  getDual(dual_fp);
1600 
1601  /* store floating point solution in exact precision */
1602  modprimal_ex = primal_fp;
1603  moddual_ex = dual_fp;
1604 
1605  /* correct primal solution */
1606  MSG_DEBUG( spxout << "correcting primal solution . . ." );
1607 
1608  modprimal_ex *= MpqReal(1 / primalscale_ex);
1609  primal_ex += modprimal_ex;
1610 
1611  /* force values of nonbasic variables to bounds */
1612  for( int c = 0; c < nCols(); c++ )
1613  {
1614  SPxBasis::Desc::Status basisstat = basisdesc.colStatus(c);
1615 
1616  if( basisstat == SPxBasis::Desc::P_ON_LOWER && primal_ex[c] != orilower_fp[c] )
1617  {
1618  primal_ex[c] = orilower_fp[c];
1619  nadjusted++;
1620  }
1621  else if( basisstat == SPxBasis::Desc::P_ON_UPPER && primal_ex[c] != oriupper_fp[c] )
1622  {
1623  primal_ex[c] = oriupper_fp[c];
1624  nadjusted++;
1625  }
1626  else if( basisstat == SPxBasis::Desc::P_FIXED )
1627  {
1628  assert(orilower_fp[c] == oriupper_fp[c]);
1629 
1630  if( primal_ex[c] != orilower_fp[c] )
1631  {
1632  primal_ex[c] = orilower_fp[c];
1633  nadjusted++;
1634  }
1635  }
1636  else if( basisstat == SPxBasis::Desc::P_FREE && primal_ex[c] != 0 )
1637  {
1638  primal_ex[c] = 0;
1639  nadjusted++;
1640  }
1641  }
1642 
1643  MSG_DEBUG( spxout << " adjusted " << nadjusted << " nonbasic variables to their bounds\n" );
1644 
1645  /* correct dual solution */
1646  MSG_DEBUG( spxout << "correcting dual solution . . .\n" );
1647 
1648  moddual_ex *= MpqReal(1 / dualscale_ex);
1649  dual_ex += moddual_ex;
1650  }
1651  }
1652  while( true );
1653 
1654  /* output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers */
1655  MSG_INFO3( spxout << "maximum violations: bounds=" << boundsviol_ex << ", sides=" << sidesviol_ex << ", duals/redcosts=" << redcostviol_ex << "\n" );
1656 
1657  MSG_DEBUG(
1658  MpqReal viol_ex = boundsviol_ex;
1659  if( sidesviol_ex > viol_ex )
1660  viol_ex = sidesviol_ex;
1661  if( redcostviol_ex > viol_ex )
1662  viol_ex = redcostviol_ex;
1663 
1664  spxout << "maxboundsviol: " << boundsviol_ex
1665  << " maxsidesviol: " << sidesviol_ex
1666  << " maxredcostviol: " << redcostviol_ex
1667  << " maxviol: " << viol_ex
1668  << " after ncorrections: " << nrefines
1669  << " --"
1670  << " maxcorrections: " << maxRefines
1671  << " primaltol: " << irfeastol
1672  << " dualtol: " << iropttol
1673  << " pivotiters: " << npivotrefines
1674  << "\n"
1675  );
1676 
1677  /* restore original problem and load last regular basis */
1678  if( changed )
1679  {
1680  MSG_DEBUG( spxout << "restoring original problem . . .\n" );
1681 
1682  changeBounds(orilower_fp, oriupper_fp);
1683  changeRange(orirhs_fp, orirhs_fp);
1684  changeObj(oriobj_fp);
1685 
1686  loadBasis(basisdesc);
1687  }
1688 
1689  /* remember total number of simplex iterations during refinement */
1690  SPxBasis::iterCount = iteroffset;
1691 
1692  /* restore iteration limit */
1693  if( terminationIter() >= 0 )
1694  setTerminationIter(terminationIter() + iteroffset);
1695 
1696  return precisionreached;
1697 }
1698 
1699 void SPxSolver::testVecs()
1700 {
1701  METHOD( "SPxSolver::testVecs()" );
1702 
1703  assert(SPxBasis::status() > SPxBasis::SINGULAR);
1704 
1705  DVector tmp(dim());
1706 
1707  tmp = *theCoPvec;
1708  multWithBase(tmp);
1709  tmp -= *theCoPrhs;
1710  if (tmp.length() > leavetol())
1711  {
1712  MSG_INFO3( spxout << "ISOLVE93 " << iteration() << ":\tcoP error = \t"
1713  << tmp.length() << std::endl; )
1714 
1715  tmp.clear();
1716  SPxBasis::coSolve(tmp, *theCoPrhs);
1717  multWithBase(tmp);
1718  tmp -= *theCoPrhs;
1719  MSG_INFO3( spxout << "ISOLVE94\t\t" << tmp.length() << std::endl; )
1720 
1721  tmp.clear();
1722  SPxBasis::coSolve(tmp, *theCoPrhs);
1723  tmp -= *theCoPvec;
1724  MSG_INFO3( spxout << "ISOLVE95\t\t" << tmp.length() << std::endl; )
1725  }
1726 
1727  tmp = *theFvec;
1728  multBaseWith(tmp);
1729  tmp -= *theFrhs;
1730  if (tmp.length() > entertol())
1731  {
1732  MSG_INFO3( spxout << "ISOLVE96 " << iteration() << ":\t F error = \t"
1733  << tmp.length() << std::endl; )
1734 
1735  tmp.clear();
1736  SPxBasis::solve(tmp, *theFrhs);
1737  tmp -= *theFvec;
1738  MSG_INFO3( spxout << "ISOLVE97\t\t" << tmp.length() << std::endl; )
1739  }
1740 
1741  if (type() == ENTER)
1742  {
1743  for (int i = 0; i < dim(); ++i)
1744  {
1745  if (theCoTest[i] < -leavetol() && isCoBasic(i))
1746  {
1747  /// @todo Error message "this shalt not be": shalt this be an assert (also below)?
1748  MSG_ERROR( spxout << "ESOLVE98 testVecs: theCoTest: this shalt not be!"
1749  << std::endl
1750  << " i=" << i
1751  << ", theCoTest[i]=" << theCoTest[i]
1752  << ", leavetol()=" << leavetol() << std::endl; )
1753  }
1754  }
1755 
1756  for (int i = 0; i < coDim(); ++i)
1757  {
1758  if (theTest[i] < -leavetol() && isBasic(i))
1759  {
1760  MSG_ERROR( spxout << "ESOLVE99 testVecs: theTest: this shalt not be!"
1761  << std::endl
1762  << " i=" << i
1763  << ", theTest[i]=" << theTest[i]
1764  << ", leavetol()=" << leavetol() << std::endl; )
1765  }
1766  }
1767  }
1768 }
1769 
1770 bool SPxSolver::terminate()
1771 {
1772  METHOD( "SPxSolver::terminate()" );
1773 #ifdef ENABLE_ADDITIONAL_CHECKS
1774  if (SPxBasis::status() > SPxBasis::SINGULAR)
1775  testVecs();
1776 #endif
1777 
1778  int redo = dim();
1779 
1780  if (redo < 1000)
1781  redo = 1000;
1782 
1783  if (iteration() > 10 && iteration() % redo == 0)
1784  {
1785 #ifdef ENABLE_ADDITIONAL_CHECKS
1786  DVector cr(*theCoPrhs);
1787  DVector fr(*theFrhs);
1788 #endif
1789 
1790  if (type() == ENTER)
1791  computeEnterCoPrhs();
1792  else
1793  computeLeaveCoPrhs();
1794 
1795  computeFrhs();
1796 
1797 #ifdef ENABLE_ADDITIONAL_CHECKS
1798  cr -= *theCoPrhs;
1799  fr -= *theFrhs;
1800  if (cr.length() > leavetol())
1801  MSG_WARNING( spxout << "WSOLVE50 unexpected change of coPrhs "
1802  << cr.length() << std::endl; )
1803  if (fr.length() > entertol())
1804  MSG_WARNING( spxout << "WSOLVE51 unexpected change of Frhs "
1805  << fr.length() << std::endl; )
1806 #endif
1807 
1808  if (updateCount > 1)
1809  {
1810  MSG_INFO3( spxout << "ISOLVE52 terminate triggers refactorization"
1811  << std::endl; )
1812  factorize();
1813  }
1814  SPxBasis::coSolve(*theCoPvec, *theCoPrhs);
1815  SPxBasis::solve (*theFvec, *theFrhs);
1816 
1817  if (pricing() == FULL)
1818  {
1819  computePvec();
1820  if (type() == ENTER)
1821  computeTest();
1822  }
1823 
1824  if (shift() > 0.0)
1825  unShift();
1826  }
1827 
1828  if ( maxTime >= 0 && maxTime < infinity && time() >= maxTime )
1829  {
1830  MSG_INFO2( spxout << "ISOLVE54 Timelimit (" << maxTime
1831  << ") reached" << std::endl; )
1832  m_status = ABORT_TIME;
1833  return true;
1834  }
1835 
1836  // objLimit is set and we are running DUAL:
1837  // - objLimit is set if objLimit < infinity
1838  // - DUAL is running if rep() * type() > 0 == DUAL (-1 == PRIMAL)
1839  //
1840  // In this case we have given a objective value limit, e.g, through a
1841  // MIP solver, and we want stop solving the LP if we figure out that the
1842  // optimal value of the current LP can not be better then this objective
1843  // limit. More precisely:
1844  // - MINIMIZATION Problem
1845  // We want stop the solving process if
1846  // objLimit <= current objective value of the DUAL LP
1847  // - MAXIMIZATION Problem
1848  // We want stop the solving process if
1849  // objLimit >= current objective value of the DUAL LP
1850  if (objLimit < infinity && type() * rep() > 0)
1851  {
1852  // We have no bound shifts; therefore, we can trust the current
1853  // objective value.
1854  // It might be even possible to use this termination value in case of
1855  // bound violations (shifting) but in this case it is quite difficult
1856  // to determine if we already reached the limit.
1857  if( shift() < epsilon() && maxInfeas() + shift() <= opttol() )
1858  {
1859  // SPxSense::MINIMIZE == -1, so we have sign = 1 on minimizing
1860  if( spxSense() * value() <= spxSense() * objLimit )
1861  {
1862  MSG_INFO2( spxout << "ISOLVE55 Objective value limit (" << objLimit
1863  << ") reached" << std::endl; )
1864  MSG_DEBUG(
1865  spxout << "DSOLVE56 Objective value limit reached" << std::endl
1866  << " (value: " << value()
1867  << ", limit: " << objLimit << ")" << std::endl
1868  << " (spxSense: " << int(spxSense())
1869  << ", rep: " << int(rep())
1870  << ", type: " << int(type()) << ")" << std::endl;
1871  )
1872 
1873  m_status = ABORT_VALUE;
1874  return true;
1875  }
1876  }
1877  }
1878 
1879  if( SPxBasis::status() >= SPxBasis::OPTIMAL ||
1880  SPxBasis::status() <= SPxBasis::SINGULAR )
1881  {
1882  m_status = UNKNOWN;
1883  return true;
1884  }
1885  return false;
1886 }
1887 
1888 SPxSolver::Status SPxSolver::getPrimal (Vector& p_vector) const
1889 {
1890  METHOD( "SPxSolver::getPrimal()" );
1891 
1892  if (!isInitialized())
1893  {
1894  /* exit if presolving/simplifier cleared the problem */
1895  if (status() == NO_PROBLEM)
1896  return status();
1897  throw SPxStatusException("XSOLVE06 Not Initialized");
1898  }
1899  if (rep() == ROW)
1900  p_vector = coPvec();
1901  else
1902  {
1903  const SPxBasis::Desc& ds = desc();
1904 
1905  for (int i = 0; i < nCols(); ++i)
1906  {
1907  switch (ds.colStatus(i))
1908  {
1909  case SPxBasis::Desc::P_ON_LOWER :
1910  p_vector[i] = SPxLP::lower(i);
1911  break;
1912  case SPxBasis::Desc::P_ON_UPPER :
1913  case SPxBasis::Desc::P_FIXED :
1914  p_vector[i] = SPxLP::upper(i);
1915  break;
1916  case SPxBasis::Desc::P_FREE :
1917  p_vector[i] = 0;
1918  break;
1919  case SPxBasis::Desc::D_FREE :
1920  case SPxBasis::Desc::D_ON_UPPER :
1921  case SPxBasis::Desc::D_ON_LOWER :
1922  case SPxBasis::Desc::D_ON_BOTH :
1923  case SPxBasis::Desc::D_UNDEFINED :
1924  break;
1925  default:
1926  throw SPxInternalCodeException("XSOLVE07 This should never happen.");
1927  }
1928  }
1929  for (int j = 0; j < dim(); ++j)
1930  {
1931  if (baseId(j).isSPxColId())
1932  p_vector[ number(SPxColId(baseId(j))) ] = fVec()[j];
1933  }
1934  }
1935  return status();
1936 }
1937 
1938 SPxSolver::Status SPxSolver::getDual (Vector& p_vector) const
1939 {
1940  METHOD( "SPxSolver::getDual()" );
1941 
1942  assert(isInitialized());
1943 
1944  if (!isInitialized())
1945  {
1946  /* exit if presolving/simplifier cleared the problem */
1947  if (status() == NO_PROBLEM)
1948  return status();
1949  throw SPxStatusException("XSOLVE08 No Problem loaded");
1950  }
1951 
1952  if (rep() == ROW)
1953  {
1954  int i;
1955  p_vector.clear ();
1956  for (i = nCols() - 1; i >= 0; --i)
1957  {
1958  if (baseId(i).isSPxRowId())
1959  p_vector[ number(SPxRowId(baseId(i))) ] = fVec()[i];
1960  }
1961  }
1962  else
1963  p_vector = coPvec();
1964 
1965  p_vector *= Real(spxSense());
1966 
1967  return status();
1968 }
1969 
1970 SPxSolver::Status SPxSolver::getRedCost (Vector& p_vector) const
1971 {
1972  METHOD( "SPxSolver::getRedCost()" );
1973 
1974  assert(isInitialized());
1975 
1976  if (!isInitialized())
1977  {
1978  throw SPxStatusException("XSOLVE09 No Problem loaded");
1979  // return NOT_INIT;
1980  }
1981 
1982  if (rep() == ROW)
1983  {
1984  int i;
1985  p_vector.clear();
1986  if (spxSense() == SPxLP::MINIMIZE)
1987  {
1988  for (i = dim() - 1; i >= 0; --i)
1989  {
1990  if (baseId(i).isSPxColId())
1991  p_vector[ number(SPxColId(baseId(i))) ] = -fVec()[i];
1992  }
1993  }
1994  else
1995  {
1996  for (i = dim() - 1; i >= 0; --i)
1997  {
1998  if (baseId(i).isSPxColId())
1999  p_vector[ number(SPxColId(baseId(i))) ] = fVec()[i];
2000  }
2001  }
2002  }
2003  else
2004  {
2005  p_vector = maxObj();
2006  p_vector -= pVec();
2007  if (spxSense() == SPxLP::MINIMIZE)
2008  p_vector *= -1.0;
2009  }
2010 
2011  return status();
2012 }
2013 
2014 SPxSolver::Status SPxSolver::getPrimalray (Vector& p_vector) const
2015 {
2016  METHOD( "SPxSolver::getPrimalray()" );
2017 
2018  assert(isInitialized());
2019 
2020  if (!isInitialized())
2021  {
2022  throw SPxStatusException("XSOLVE10 No Problem loaded");
2023  // return NOT_INIT;
2024  }
2025 
2026  assert(SPxBasis::status() == SPxBasis::UNBOUNDED);
2027  p_vector.clear();
2028  p_vector = primalRay;
2029 
2030  return status();
2031 }
2032 
2033 SPxSolver::Status SPxSolver::getDualfarkas (Vector& p_vector) const
2034 {
2035  METHOD( "SPxSolver::getDualfarkas()" );
2036 
2037  assert(isInitialized());
2038 
2039  if (!isInitialized())
2040  {
2041  throw SPxStatusException("XSOLVE10 No Problem loaded");
2042  // return NOT_INIT;
2043  }
2044 
2045  assert(SPxBasis::status() == SPxBasis::INFEASIBLE);
2046  p_vector.clear();
2047  p_vector = dualFarkas;
2048 
2049  return status();
2050 }
2051 
2052 SPxSolver::Status SPxSolver::getSlacks (Vector& p_vector) const
2053 {
2054  METHOD( "SPxSolver::getSlacks()" );
2055 
2056  assert(isInitialized());
2057 
2058  if (!isInitialized())
2059  {
2060  throw SPxStatusException("XSOLVE11 No Problem loaded");
2061  // return NOT_INIT;
2062  }
2063 
2064  if (rep() == COLUMN)
2065  {
2066  int i;
2067  const SPxBasis::Desc& ds = desc();
2068  for (i = nRows() - 1; i >= 0; --i)
2069  {
2070  switch (ds.rowStatus(i))
2071  {
2072  case SPxBasis::Desc::P_ON_LOWER :
2073  p_vector[i] = lhs(i);
2074  break;
2075  case SPxBasis::Desc::P_ON_UPPER :
2076  case SPxBasis::Desc::P_FIXED :
2077  p_vector[i] = rhs(i);
2078  break;
2079  case SPxBasis::Desc::P_FREE :
2080  p_vector[i] = 0;
2081  break;
2082  case SPxBasis::Desc::D_FREE :
2083  case SPxBasis::Desc::D_ON_UPPER :
2084  case SPxBasis::Desc::D_ON_LOWER :
2085  case SPxBasis::Desc::D_ON_BOTH :
2086  case SPxBasis::Desc::D_UNDEFINED :
2087  break;
2088  default:
2089  throw SPxInternalCodeException("XSOLVE12 This should never happen.");
2090  }
2091  }
2092  for (i = dim() - 1; i >= 0; --i)
2093  {
2094  if (baseId(i).isSPxRowId())
2095  p_vector[ number(SPxRowId(baseId(i))) ] = -(*theFvec)[i];
2096  }
2097  }
2098  else
2099  p_vector = pVec();
2100 
2101  return status();
2102 }
2103 
2104 void SPxSolver::setPrimal(Vector& p_vector)
2105 {
2106  METHOD( "SPxSolver::setPrimal()" );
2107 
2108  if (!isInitialized())
2109  {
2110  throw SPxStatusException("XSOLVE20 Not Initialized");
2111  }
2112 
2113  if (rep() == ROW)
2114  coPvec() = p_vector;
2115  else
2116  {
2117  for (int j = 0; j < dim(); ++j)
2118  {
2119  if (baseId(j).isSPxColId())
2120  fVec()[j] = p_vector[ number(SPxColId(baseId(j))) ];
2121  }
2122  }
2123 }
2124 
2125 void SPxSolver::setDual(Vector& p_vector)
2126 {
2127  METHOD( "SPxSolver::setDual()" );
2128 
2129  assert(isInitialized());
2130 
2131  if (!isInitialized())
2132  {
2133  throw SPxStatusException("XSOLVE21 Not Initialized");
2134  }
2135 
2136  if (rep() == ROW)
2137  {
2138  for (int i = nCols() - 1; i >= 0; --i)
2139  {
2140  if (baseId(i).isSPxRowId())
2141  {
2142  if (spxSense() == SPxLP::MAXIMIZE)
2143  fVec()[i] = p_vector[ number(SPxRowId(baseId(i))) ];
2144  else
2145  fVec()[i] = -p_vector[ number(SPxRowId(baseId(i))) ];
2146  }
2147  }
2148  }
2149  else
2150  {
2151  coPvec() = p_vector;
2152  if (spxSense() == SPxLP::MINIMIZE)
2153  coPvec() *= -1.0;
2154  }
2155 }
2156 
2157 void SPxSolver::setRedCost(Vector& p_vector)
2158 {
2159  METHOD( "SPxSolver::setRedCost()" );
2160 
2161  assert(isInitialized());
2162 
2163  if (!isInitialized())
2164  {
2165  throw SPxStatusException("XSOLVE22 Not Initialized");
2166  }
2167 
2168  if (rep() == ROW)
2169  {
2170  for (int i = dim() - 1; i >= 0; --i)
2171  {
2172  if (baseId(i).isSPxColId())
2173  {
2174  if (spxSense() == SPxLP::MINIMIZE)
2175  fVec()[i] = -p_vector[ number(SPxColId(baseId(i))) ];
2176  else
2177  fVec()[i] = p_vector[ number(SPxColId(baseId(i))) ];
2178  }
2179  }
2180  }
2181  else
2182  {
2183  pVec() = maxObj();
2184 
2185  if (spxSense() == SPxLP::MINIMIZE)
2186  pVec() += p_vector;
2187  else
2188  pVec() -= p_vector;
2189  }
2190 }
2191 
2192 void SPxSolver::setSlacks(Vector& p_vector)
2193 {
2194  METHOD( "SPxSolver::getSlacks()" );
2195 
2196  assert(isInitialized());
2197 
2198  if (!isInitialized())
2199  {
2200  throw SPxStatusException("XSOLVE23 Not Initialized");
2201  }
2202 
2203  if (rep() == COLUMN)
2204  {
2205  for (int i = dim() - 1; i >= 0; --i)
2206  {
2207  if (baseId(i).isSPxRowId())
2208  (*theFvec)[i] = -p_vector[ number(SPxRowId(baseId(i))) ];
2209  }
2210  }
2211  else
2212  pVec() = p_vector;
2213 }
2214 
2215 SPxSolver::Status SPxSolver::status() const
2216 {
2217  METHOD( "SPxSolver::status()" );
2218  switch( m_status )
2219  {
2220  case UNKNOWN :
2221  switch (SPxBasis::status())
2222  {
2223  case SPxBasis::NO_PROBLEM :
2224  return NO_PROBLEM;
2225  case SPxBasis::SINGULAR :
2226  return SINGULAR;
2227  case SPxBasis::REGULAR :
2228  case SPxBasis::DUAL :
2229  case SPxBasis::PRIMAL :
2230  return UNKNOWN;
2231  case SPxBasis::OPTIMAL :
2232  return OPTIMAL;
2233  case SPxBasis::UNBOUNDED :
2234  return UNBOUNDED;
2235  case SPxBasis::INFEASIBLE :
2236  return INFEASIBLE;
2237  default:
2238  return ERROR;
2239  }
2240  case SINGULAR :
2241  return m_status;
2242  case OPTIMAL :
2243  assert( SPxBasis::status() == SPxBasis::OPTIMAL );
2244  /*lint -fallthrough*/
2245  case ABORT_CYCLING :
2246  case ABORT_TIME :
2247  case ABORT_ITER :
2248  case ABORT_VALUE :
2249  case RUNNING :
2250  case REGULAR :
2251  case NOT_INIT :
2252  case NO_SOLVER :
2253  case NO_PRICER :
2254  case NO_RATIOTESTER :
2255  case ERROR:
2256  return m_status;
2257  default:
2258  return ERROR;
2259  }
2260 }
2261 
2262 SPxSolver::Status SPxSolver::getResult(
2263  Real* p_value,
2264  Vector* p_primal,
2265  Vector* p_slacks,
2266  Vector* p_dual,
2267  Vector* reduCosts) const
2268 {
2269  METHOD( "SPxSolver::getResult()" );
2270  if (p_value)
2271  *p_value = this->value();
2272  if (p_primal)
2273  getPrimal(*p_primal);
2274  if (p_slacks)
2275  getSlacks(*p_slacks);
2276  if (p_dual)
2277  getDual(*p_dual);
2278  if (reduCosts)
2279  getRedCost(*reduCosts);
2280  return status();
2281 }
2282 } // namespace soplex
2283 
2284 //-----------------------------------------------------------------------------
2285 //Emacs Local Variables:
2286 //Emacs mode:c++
2287 //Emacs c-basic-offset:3
2288 //Emacs tab-width:8
2289 //Emacs indent-tabs-mode:nil
2290 //Emacs End:
2291 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen