SoPlex Doxygen Documentation
spxsolver.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 #include <sstream>
21 
22 #include "spxdefines.h"
23 #include "soplex.h"
24 #include "spxpricer.h"
25 #include "spxratiotester.h"
26 #include "spxstarter.h"
27 #include "spxout.h"
28 #include "exceptions.h"
29 
30 namespace soplex
31 {
32 #define MAXIMUM(x,y) ((x)>(y) ? (x) : (y))
33 
34 bool SPxSolver::read(std::istream& in, NameSet* rowNames,
35  NameSet* colNames, DIdxSet* intVars)
36 {
37  METHOD( "SPxSolver::read()" );
38  clear();
39  unInit();
40  unLoad();
41 
42  if (thepricer)
43  thepricer->clear();
44 
45  if (theratiotester)
46  theratiotester->clear();
47 
48  if (!SPxLP::read(in, rowNames, colNames, intVars))
49  return false;
50 
51  SPxBasis::load(this);
52 
53  return true;
54 }
55 
56 void SPxSolver::reLoad()
57 {
58  METHOD( "SPxSolver::reLoad()" );
59  unInit();
60  unLoad();
61  theLP = this;
62  if (thepricer)
63  thepricer->clear();
64  if (theratiotester)
65  theratiotester->clear();
66 }
67 
68 void SPxSolver::loadLP(const SPxLP& lp)
69 {
70  METHOD( "SPxSolver::loadLP()" );
71  clear();
72  unInit();
73  unLoad();
74  if (thepricer)
75  thepricer->clear();
76  if (theratiotester)
77  theratiotester->clear();
78  SPxLP::operator=(lp);
79  reDim();
80  SPxBasis::load(this);
81 }
82 
83 void SPxSolver::setSolver(SLinSolver* slu, const bool destroy)
84 {
85  METHOD( "SPxSolver::setSolver()" );
86  SPxBasis::loadSolver(slu, destroy);
87 }
88 
89 void SPxSolver::loadBasis(const SPxBasis::Desc& p_desc)
90 {
91  METHOD( "SPxSolver::loadBasis()" );
92  unInit();
93  if (SPxBasis::status() == SPxBasis::NO_PROBLEM)
94  SPxBasis::load(this);
95  SPxBasis::loadDesc(p_desc);
96  setBasisStatus(SPxBasis::REGULAR);
97 }
98 
99 void SPxSolver::setPricer(SPxPricer* x, const bool destroy)
100 {
101  METHOD( "SPxSolver::setPricer()" );
102 
103  assert(!freePricer || thepricer != 0);
104 
105  if(freePricer)
106  {
107  delete thepricer;
108  thepricer = 0;
109  }
110 
111  if (x != 0 && x != thepricer)
112  {
113  setPricing(FULL);
114  if (isInitialized())
115  x->load(this);
116  else
117  x->clear();
118  }
119 
120  if (thepricer && thepricer != x)
121  thepricer->clear();
122  thepricer = x;
123 
124  freePricer = destroy;
125 }
126 
127 void SPxSolver::setTester(SPxRatioTester* x, const bool destroy)
128 {
129  METHOD( "SPxSolver::setTester()" );
130 
131  assert(!freeRatioTester || theratiotester != 0);
132 
133  if(freeRatioTester)
134  {
135  delete theratiotester;
136  theratiotester = 0;
137  }
138 
139  if (x)
140  {
141  if (isInitialized() && x != theratiotester)
142  x->load(this);
143  else
144  x->clear();
145  }
146 
147  if (theratiotester !=0 && theratiotester != x)
148  theratiotester->clear();
149  theratiotester = x;
150 
151  freeRatioTester = destroy;
152 }
153 
154 void SPxSolver::setStarter(SPxStarter* x, const bool destroy)
155 {
156  METHOD( "SPxSolver::setStarter()" );
157 
158  assert(!freeStarter || thestarter != 0);
159 
160  if(freeStarter)
161  {
162  delete thestarter;
163  thestarter = 0;
164  }
165  thestarter = x;
166 
167  freeStarter = destroy;
168 }
169 
170 void SPxSolver::setType(Type tp)
171 {
172  METHOD( "SPxSolver::setType()" );
173 
174  if (theType != tp)
175  {
176  theType = tp;
177 
178  unInit();
179 #if 0
180  else
181  {
182  if (!matrixIsSetup)
183  {
184  SPxBasis::load(this);
185  // SPxBasis::load(desc());
186  // not needed, because load(this) allready loads descriptor
187  }
188  factorized = false;
189  m_numCycle = 0;
190 #endif
191  MSG_DEBUG( spxout << "DSOLVE20 switching to "
192  << static_cast<const char*>((tp == LEAVE)
193  ? "leaving" : "entering")
194  << " algorithm" << std::endl; )
195  }
196 }
197 
198 void SPxSolver::initRep(Representation p_rep)
199 {
200  METHOD( "SPxSolver::initRep()" );
201 
202  theRep = p_rep;
203 
204  Real tmpfeastol = feastol();
205  Real tmpopttol = opttol();
206 
207  if (theRep == COLUMN)
208  {
209  thevectors = colSet();
210  thecovectors = rowSet();
211  theFrhs = &primRhs;
212  theFvec = &primVec;
213  theCoPrhs = &dualRhs;
214  theCoPvec = &dualVec;
215  thePvec = &addVec;
216  theRPvec = theCoPvec;
217  theCPvec = thePvec;
218  theUbound = &theUCbound;
219  theLbound = &theLCbound;
220  theCoUbound = &theURbound;
221  theCoLbound = &theLRbound;
222  }
223  else
224  {
225  assert(theRep == ROW);
226 
227  thevectors = rowSet();
228  thecovectors = colSet();
229  theFrhs = &dualRhs;
230  theFvec = &dualVec;
231  theCoPrhs = &primRhs;
232  theCoPvec = &primVec;
233  thePvec = &addVec;
234  theRPvec = thePvec;
235  theCPvec = theCoPvec;
236  theUbound = &theURbound;
237  theLbound = &theLRbound;
238  theCoUbound = &theUCbound;
239  theCoLbound = &theLCbound;
240  }
241  unInit();
242  reDim();
243 
244  setFeastol(tmpfeastol);
245  setOpttol(tmpopttol);
246 
247  SPxBasis::setRep();
248  if (SPxBasis::status() > SPxBasis::NO_PROBLEM)
249  SPxBasis::loadDesc(desc());
250 
251  if (thepricer && thepricer->solver() == this)
252  thepricer->setRep(p_rep);
253 }
254 
255 void SPxSolver::setRep(Representation p_rep)
256 {
257  METHOD( "SPxSolver::setRep()" );
258 
259  if (p_rep != theRep)
260  initRep(p_rep);
261 }
262 
263 // needed for strongbranching. use carefully
264 void SPxSolver::reinitializeVecs()
265 {
266  METHOD( "SPxSolver::reinitializeVecs" );
267 
268  initialized = true;
269 
270  if (type() == ENTER)
271  {
272  if (rep() == COLUMN)
273  setPrimalBounds();
274  else
275  setDualRowBounds();
276 
277  setEnterBounds();
278  computeEnterCoPrhs();
279  }
280  else
281  {
282  if (rep() == ROW)
283  setPrimalBounds();
284  else
285  setDualColBounds();
286 
287  setLeaveBounds();
288  computeLeaveCoPrhs();
289  }
290 
291  SPxBasis::coSolve(*theCoPvec, *theCoPrhs);
292  computePvec();
293  computeFrhs();
294  SPxBasis::solve(*theFvec, *theFrhs);
295 
296  theShift = 0.0;
297  lastShift = 0.0;
298 
299  if (type() == ENTER)
300  {
301  computeCoTest();
302  computeTest();
303  }
304  else
305  {
306  computeFtest();
307  }
308  assert((testBounds(), 1));
309 }
310 
311 void SPxSolver::init()
312 {
313  METHOD( "SPxSolver::init()" );
314 
315  assert(thepricer != 0);
316  assert(theratiotester != 0);
317 
318  if (!initialized)
319  {
320  initialized = true;
321  reDim();
322  if (SPxBasis::status() <= SPxBasis::NO_PROBLEM || solver() != this)
323  SPxBasis::load(this);
324  initialized = false;
325  }
326  if (!matrixIsSetup)
327  SPxBasis::loadDesc(desc());
328 
329  // Inna/Tobi: don't "upgrade" a singular basis to a regular one
330  if( SPxBasis::status() == SPxBasis::SINGULAR )
331  return;
332 
333  //factorized = false;
334  m_numCycle = 0;
335 
336  if (type() == ENTER)
337  {
338  if (rep() == COLUMN)
339  {
340  setPrimalBounds();
341  setBasisStatus(SPxBasis::PRIMAL);
342  }
343  else
344  {
345  setDualRowBounds();
346  setBasisStatus(SPxBasis::DUAL);
347  }
348  setEnterBounds();
349  computeEnterCoPrhs();
350  // prepare support vectors for sparse pricing
351  infeasibilities.setMax(dim());
352  infeasibilitiesCo.setMax(coDim());
353  isInfeasible.reSize(dim());
354  isInfeasibleCo.reSize(coDim());
355  sparsityThresholdEnter = (int) (dim() * SPARSITYTHRESHOLD);
356  sparsityThresholdEnterCo = (int) (coDim() * SPARSITYTHRESHOLD);
357  theratiotester->setDelta(entertol());
358  }
359  else
360  {
361  if (rep() == ROW)
362  {
363  setPrimalBounds();
364  setBasisStatus(SPxBasis::PRIMAL);
365  }
366  else
367  {
368  setDualColBounds();
369  setBasisStatus(SPxBasis::DUAL);
370  }
371  setLeaveBounds();
372  computeLeaveCoPrhs();
373  // prepare support vectors for sparse pricing
374  infeasibilities.setMax(dim());
375  isInfeasible.reSize(dim());
376  sparsityThresholdLeave = (int) (dim() * SPARSITYTHRESHOLD);
377  theratiotester->setDelta(leavetol());
378  }
379 
380  SPxBasis::coSolve(*theCoPvec, *theCoPrhs);
381  computePvec();
382  computeFrhs();
383  SPxBasis::solve(*theFvec, *theFrhs);
384 
385  theShift = 0.0;
386 
387  if (type() == ENTER)
388  {
389  shiftFvec();
390  lastShift = theShift + entertol();
391 
392  computeCoTest();
393  computeTest();
394  }
395  else
396  {
397  shiftPvec();
398  lastShift = theShift + leavetol();
399 
400  computeFtest();
401  }
402 
403  if (!initialized)
404  {
405  // if(thepricer->solver() != this)
406  thepricer->load(this);
407  // if(theratiotester->solver() != this)
408  theratiotester->load(this);
409  initialized = true;
410  }
411 }
412 
413 void SPxSolver::setPricing(Pricing pr)
414 {
415  METHOD( "SPxSolver::setPricing()" );
416  thePricing = pr;
417  if (initialized && type() == ENTER)
418  {
419  computePvec();
420  computeCoTest();
421  computeTest();
422  }
423 }
424 
425 /*
426  The following method resizes all vectors and arrays of |SoPlex|
427  (excluding inherited vectors).
428  */
429 void SPxSolver::reDim()
430 {
431  METHOD( "SPxSolver::reDim()" );
432 
433  int newsize = SPxLP::nCols() > SPxLP::nRows() ? SPxLP::nCols() : SPxLP::nRows();
434 
435  if (newsize > unitVecs.size())
436  {
437  unitVecs.reSize(newsize);
438 
439  while (newsize-- > 0)
440  unitVecs[newsize] = UnitVector(newsize);
441  }
442 
443  if (isInitialized())
444  {
445  theFrhs->reDim(dim());
446  theFvec->reDim(dim());
447  thePvec->reDim(coDim());
448 
449  theCoPrhs->reDim(dim());
450  theCoPvec->reDim(dim());
451 
452  theTest.reDim(coDim());
453  theCoTest.reDim(dim());
454 
455  theURbound.reDim(SPxLP::nRows());
456  theLRbound.reDim(SPxLP::nRows());
457  theUCbound.reDim(SPxLP::nCols());
458  theLCbound.reDim(SPxLP::nCols());
459  theUBbound.reDim(dim());
460  theLBbound.reDim(dim());
461  }
462 }
463 
464 void SPxSolver::clear()
465 {
466  METHOD( "SPxSolver::clear()" );
467  unitVecs.reSize(0);
468 
469  dualRhs.clear();
470  dualVec.clear();
471  primRhs.clear();
472  primVec.clear();
473  addVec.clear();
474  theURbound.clear();
475  theLRbound.clear();
476  theUCbound.clear();
477  theLCbound.clear();
478  theTest.clear();
479  theCoTest.clear();
480 
481  unInit();
482  SPxLP::clear();
483  setBasisStatus(SPxBasis::NO_PROBLEM);
484  SPxBasis::reDim();
485 
486  infeasibilities.clear();
487  infeasibilitiesCo.clear();
488  isInfeasible.clear();
489  isInfeasibleCo.clear();
490 }
491 
492 void SPxSolver::clearUpdateVecs(void)
493 {
494  METHOD( "SPxSolver::clearUpdateVecs()" );
495  theFvec->clearUpdate();
496  thePvec->clearUpdate();
497  theCoPvec->clearUpdate();
498  solveVector2 = 0;
499  solveVector3 = 0;
500  coSolveVector2 = 0;
501 }
502 
503 /*
504  When the basis matrix factorization is recomputed from scratch,
505  we also recompute the vectors.
506  */
507 void SPxSolver::factorize()
508 {
509  METHOD( "SPxSolver::factorize()" );
510 
511  MSG_INFO1( spxout << "ISOLVE01 "
512  << "iteration = " << std::setw(8) << basis().iteration()
513  << "\tlastUpdate = " << std::setw(4) << basis().lastUpdate()
514  << "\tvalue = " << (isInitialized() ? value() : 0.0)
515  << std::endl; )
516 
517  SPxBasis::factorize();
518 
519  if (SPxBasis::status() >= SPxBasis::REGULAR)
520  {
521 #ifndef NDEBUG
522  DVector ftmp(fVec());
523  DVector ptmp(pVec());
524  DVector ctmp(coPvec());
525 #endif // NDEBUG
526 
527  if (type() == LEAVE)
528  {
529  /* we have to recompute theFrhs, because roundoff errors can occur during updating, especially when
530  * columns/rows with large bounds are present
531  */
532  computeFrhs();
533  SPxBasis::solve (*theFvec, *theFrhs);
534  SPxBasis::coSolve(*theCoPvec, *theCoPrhs);
535 
536 #ifndef NDEBUG
537  ftmp -= fVec();
538  ptmp -= pVec();
539  ctmp -= coPvec();
540  if (ftmp.length() > DEFAULT_BND_VIOL)
541  {
542  MSG_DEBUG( spxout << "DSOLVE21 fVec: " << ftmp.length() << std::endl; )
543  ftmp = fVec();
544  multBaseWith(ftmp);
545  ftmp -= fRhs();
546  if (ftmp.length() > DEFAULT_BND_VIOL)
547  MSG_ERROR( spxout << "ESOLVE29 " << iteration() << ": fVec error = "
548  << ftmp.length() << " exceeding DEFAULT_BND_VIOL = " << DEFAULT_BND_VIOL << std::endl; )
549  }
550  if (ctmp.length() > DEFAULT_BND_VIOL)
551  {
552  MSG_DEBUG( spxout << "DSOLVE23 coPvec: " << ctmp.length() << std::endl; )
553  ctmp = coPvec();
554  multWithBase(ctmp);
555  ctmp -= coPrhs();
556  if (ctmp.length() > DEFAULT_BND_VIOL)
557  MSG_ERROR( spxout << "ESOLVE30 " << iteration() << ": coPvec error = "
558  << ctmp.length() << " exceeding DEFAULT_BND_VIOL = " << DEFAULT_BND_VIOL << std::endl; )
559  }
560  if (ptmp.length() > DEFAULT_BND_VIOL)
561  {
562  MSG_DEBUG( spxout << "DSOLVE24 pVec: " << ptmp.length() << std::endl; )
563  }
564 #endif // NDEBUG
565 
566  computeFtest();
567 #if 0 /* was deactivated */
568  computePvec();
569 #endif
570  }
571  else
572  {
573  assert(type() == ENTER);
574 
575  computeCoTest();
576  if (pricing() == FULL)
577  {
578 #if 0 /* was deactivated */
579  computePvec();
580 #endif
581  /* was deactivated, but this leads to warnings in testVecs() */
582  computeTest();
583  }
584  }
585  }
586 
587  if (SPxBasis::status() == SPxBasis::SINGULAR)
588  {
589  m_status = SINGULAR;
590  std::stringstream s;
591  s << "XSOLVE21 Basis is singular (numerical troubles, feastol = " << feastol() << ", opttol = " << opttol() << ")";
592  throw SPxStatusException(s.str());
593  }
594 
595 #ifdef ENABLE_ADDITIONAL_CHECKS
596  /* moved this test after the computation of fTest and coTest below, since these vectors might not be set up at top, e.g. for an initial basis */
597  if (SPxBasis::status() > SPxBasis::SINGULAR)
598  testVecs();
599 #endif
600 }
601 
602 /* We compute how much the current solution violates (primal or dual) feasibility. In the
603  row/enter or column/leave algorithm the maximum violation of dual feasibility is
604  computed. In the row/leave or column/enter algorithm the primal feasibility is checked. */
605 Real SPxSolver::maxInfeas() const
606 {
607  METHOD( "SPxSolver::maxInfeas()" );
608  Real inf = 0.0;
609 
610  if (type() == ENTER)
611  {
612  for (int i = 0; i < dim(); i++)
613  {
614  if ((*theFvec)[i] > theUBbound[i])
615  inf = MAXIMUM(inf, (*theFvec)[i] - theUBbound[i]);
616  if (theLBbound[i] > (*theFvec)[i])
617  inf = MAXIMUM(inf, theLBbound[i] - (*theFvec)[i]);
618  }
619  }
620  else
621  {
622  assert(type() == LEAVE);
623 
624  for (int i = 0; i < dim(); i++)
625  {
626  if ((*theCoPvec)[i] > (*theCoUbound)[i])
627  inf = MAXIMUM(inf, (*theCoPvec)[i] - (*theCoUbound)[i]);
628  if ((*theCoLbound)[i] > (*theCoPvec)[i])
629  inf = MAXIMUM(inf, (*theCoLbound)[i] - (*theCoPvec)[i]);
630  }
631  for (int i = 0; i < coDim(); i++)
632  {
633  if ((*thePvec)[i] > (*theUbound)[i])
634  inf = MAXIMUM(inf, (*thePvec)[i] - (*theUbound)[i]);
635  else if ((*thePvec)[i] < (*theLbound)[i])
636  inf = MAXIMUM(inf, (*theLbound)[i] - (*thePvec)[i]);
637  }
638  }
639 
640  return inf;
641 }
642 
643 Real SPxSolver::nonbasicValue() const
644 {
645  METHOD( "SPxSolver::nonbasicValue()" );
646 
647  int i;
648  Real val = 0;
649  const SPxBasis::Desc& ds = desc();
650 
651  if (rep() == COLUMN)
652  {
653  if (type() == LEAVE)
654  {
655  for (i = nCols() - 1; i >= 0; --i)
656  {
657  switch (ds.colStatus(i))
658  {
659  case SPxBasis::Desc::P_ON_UPPER :
660  val += theUCbound[i] * SPxLP::upper(i);
661  //@ val += maxObj(i) * SPxLP::upper(i);
662  break;
663  case SPxBasis::Desc::P_ON_LOWER :
664  val += theLCbound[i] * SPxLP::lower(i);
665  //@ val += maxObj(i) * SPxLP::lower(i);
666  break;
667  case SPxBasis::Desc::P_ON_UPPER + SPxBasis::Desc::P_ON_LOWER :
668  val += maxObj(i) * SPxLP::lower(i);
669  break;
670  default:
671  break;
672  }
673  }
674  for (i = nRows() - 1; i >= 0; --i)
675  {
676  switch (ds.rowStatus(i))
677  {
678  case SPxBasis::Desc::P_ON_UPPER :
679  val += theLRbound[i] * SPxLP::rhs(i);
680  break;
681  case SPxBasis::Desc::P_ON_LOWER :
682  val += theURbound[i] * SPxLP::lhs(i);
683  break;
684  default:
685  break;
686  }
687  }
688  }
689  else
690  {
691  assert(type() == ENTER);
692  for (i = nCols() - 1; i >= 0; --i)
693  {
694  switch (ds.colStatus(i))
695  {
696  case SPxBasis::Desc::P_ON_UPPER :
697  val += maxObj(i) * theUCbound[i];
698  break;
699  case SPxBasis::Desc::P_ON_LOWER :
700  val += maxObj(i) * theLCbound[i];
701  break;
702  case SPxBasis::Desc::P_ON_UPPER + SPxBasis::Desc::P_ON_LOWER :
703  assert(theLCbound[i] == theUCbound[i]);
704  val += maxObj(i) * theLCbound[i];
705  break;
706  default:
707  break;
708  }
709  }
710  }
711  }
712  else
713  {
714  assert(rep() == ROW);
715  assert(type() == ENTER);
716  for (i = nCols() - 1; i >= 0; --i)
717  {
718  switch (ds.colStatus(i))
719  {
720  case SPxBasis::Desc::D_ON_UPPER :
721  val += theUCbound[i] * lower(i);
722  break;
723  case SPxBasis::Desc::D_ON_LOWER :
724  val += theLCbound[i] * upper(i);
725  break;
726  case SPxBasis::Desc::D_ON_BOTH :
727  val += theLCbound[i] * upper(i);
728  val += theUCbound[i] * lower(i);
729  break;
730  default:
731  break;
732  }
733  }
734  for (i = nRows() - 1; i >= 0; --i)
735  {
736  switch (ds.rowStatus(i))
737  {
738  case SPxBasis::Desc::D_ON_UPPER :
739  val += theURbound[i] * lhs(i);
740  break;
741  case SPxBasis::Desc::D_ON_LOWER :
742  val += theLRbound[i] * rhs(i);
743  break;
744  case SPxBasis::Desc::D_ON_BOTH :
745  val += theLRbound[i] * rhs(i);
746  val += theURbound[i] * lhs(i);
747  break;
748  default:
749  break;
750  }
751  }
752  }
753 
754  return val;
755 }
756 
757 Real SPxSolver::value() const
758 {
759  METHOD( "SPxSolver::value()" );
760  Real x;
761 
762  assert(isInitialized());
763 
764  // calling value() without having a suitable status is an error.
765  if (!isInitialized())
766  return infinity;
767 
768  if (rep() == ROW)
769  {
770  if (type() == LEAVE)
771  x = SPxLP::spxSense() * (coPvec() * fRhs());
772  else
773  x = SPxLP::spxSense() * (nonbasicValue() + (coPvec() * fRhs()));
774  }
775  else
776  x = SPxLP::spxSense() * (nonbasicValue() + fVec() * coPrhs());
777 
778  return x;
779 }
780 
781 void SPxSolver::setFeastol(Real d)
782 {
783  METHOD( "SPxSolver::setFeastol()" );
784 
785  if( d < 0.0 )
786  throw SPxInterfaceException("XSOLVE30 Cannot set negative feastol.");
787 
788  if( !MpqRealIsExact() && d < DEFAULT_BND_VIOL * 1e-6 )
789  {
790  MSG_WARNING( spxout << "WSOLVE32 Warning: Cannot set primal feasibility tolerance smaller than " << DEFAULT_BND_VIOL * 1e-6 << " because of missing GMP support (compile with GMP=true).\n" );
791  d = DEFAULT_BND_VIOL * 1e-6;
792  }
793 
794  if( theRep == COLUMN )
795  m_entertol = d;
796  else
797  m_leavetol = d;
798 }
799 
800 void SPxSolver::setOpttol(Real d)
801 {
802  METHOD( "SPxSolver::setOpttol()" );
803 
804  if( d < 0.0 )
805  throw SPxInterfaceException("XSOLVE31 Cannot set negative opttol.");
806 
807  if( !MpqRealIsExact() && d < DEFAULT_BND_VIOL * 1e-6 )
808  {
809  MSG_WARNING( spxout << "WSOLVE33 Warning: Cannot set dual feasibility tolerance smaller than " << DEFAULT_BND_VIOL * 1e-6 << " because of missing GMP support (compile with GMP=true).\n" );
810  d = DEFAULT_BND_VIOL * 1e-6;
811  }
812 
813  if( theRep == COLUMN )
814  m_leavetol = d;
815  else
816  m_entertol = d;
817 }
818 
819 void SPxSolver::setDelta(Real d)
820 {
821  METHOD( "SPxSolver::setDelta()" );
822 
823  if( d < 0.0 )
824  throw SPxInterfaceException("XSOLVE32 Cannot set negative delta.");
825 
826  if( !MpqRealIsExact() && d < DEFAULT_BND_VIOL * 1e-6 )
827  {
828  MSG_WARNING( spxout << "WSOLVE34 Warning: Cannot set feasibility tolerance smaller than " << DEFAULT_BND_VIOL * 1e-6 << " because of missing GMP support (compile with GMP=true).\n" );
829  d = DEFAULT_BND_VIOL * 1e-6;
830  }
831 
832  m_entertol = d;
833  m_leavetol = d;
834 }
835 
836 void SPxSolver::setIrthreshold(Real d)
837 {
838  METHOD( "SPxSolver::setIrthreshold()" );
839 
840  if( d <= 0.0 )
841  throw SPxInterfaceException("XSOLVE33 Cannot set negative or zero irthreshold.");
842 
843  m_irthreshold = d;
844 }
845 
846 SPxSolver::SPxSolver(
847  Type p_type,
848  Representation p_rep )
849  : theType (p_type)
850  , thePricing(FULL)
851  , theCumulativeTime(0.0)
852  , maxIters (-1)
853  , maxRefines (100)
854  , maxTime (infinity)
855  , objLimit(infinity)
856  , m_status(UNKNOWN)
857  , theShift (0)
858  , m_maxCycle(100)
859  , m_numCycle(0)
860  , initialized (false)
861  , solveVector2 (0)
862  , solveVector3 (0)
863  , coSolveVector2(0)
864  , freePricer (false)
865  , freeRatioTester (false)
866  , freeStarter (false)
867  , unitVecs (0)
868  , primVec (0, Param::epsilon())
869  , dualVec (0, Param::epsilon())
870  , addVec (0, Param::epsilon())
871  , thepricer (0)
872  , theratiotester (0)
873  , thestarter (0)
874  , infeasibilities(0)
875  , infeasibilitiesCo(0)
876  , isInfeasible(0)
877  , isInfeasibleCo(0)
878  , sparsePricingLeave(false)
879  , sparsePricingEnter(false)
880  , sparsePricingEnterCo(false)
881  , remainingRoundsLeave(0)
882  , remainingRoundsEnter(0)
883  , remainingRoundsEnterCo(0)
884 {
885  METHOD( "SPxSolver::SPxSolver()" );
886 
887  setDelta(DEFAULT_BND_VIOL);
888  setIrthreshold(DEFAULT_BND_VIOL * 1e-6);
889 
890  theLP = this;
891  initRep(p_rep);
892 
893  // info: SPxBasis is not consistent in this moment.
894  //assert(SPxSolver::isConsistent());
895 }
896 
897 SPxSolver::~SPxSolver()
898 {
899  assert(!freePricer || thepricer != 0);
900  assert(!freeRatioTester || theratiotester != 0);
901  assert(!freeStarter || thestarter != 0);
902 
903  if(freePricer)
904  {
905  delete thepricer;
906  thepricer = 0;
907  }
908 
909  if(freeRatioTester)
910  {
911  delete theratiotester;
912  theratiotester = 0;
913  }
914 
915  if(freeStarter)
916  {
917  delete thestarter;
918  thestarter = 0;
919  }
920 }
921 
922 
923 SPxSolver& SPxSolver::operator=(const SPxSolver& base)
924 {
925  METHOD( "SPxSolver::operator=(const SPxSolver&base)" );
926  if(this != &base)
927  {
928  SPxLP::operator=(base);
929  SPxBasis::operator=(base);
930  theType = base.theType;
931  thePricing = base.thePricing;
932  theRep = base.theRep;
933  theTime = base.theTime;
934  maxIters = base.maxIters;
935  maxRefines = base.maxRefines;
936  maxTime = base.maxTime;
937  objLimit = base.objLimit;
938  m_status = base.m_status;
939  m_entertol = base.m_entertol;
940  m_leavetol = base.m_leavetol;
941  m_irthreshold = base.m_irthreshold;
942  theShift = base.theShift;
943  lastShift = base.lastShift;
944  m_maxCycle = base.m_maxCycle;
945  m_numCycle = base.m_numCycle;
946  initialized = base.initialized;
947  instableLeaveNum = base.instableLeaveNum;
948  instableLeave = base.instableLeave;
949  unitVecs = base.unitVecs;
950  primRhs = base.primRhs;
951  primVec = base.primVec;
952  dualRhs = base.dualRhs;
953  dualVec = base.dualVec;
954  addVec = base.addVec;
955  theURbound = base.theURbound;
956  theLRbound = base.theLRbound;
957  theUCbound = base.theUCbound;
958  theLCbound = base.theLCbound;
959  theUBbound = base.theUBbound;
960  theLBbound = base.theLBbound;
961  theCoTest = base.theCoTest;
962  theTest = base.theTest;
963  primalRay = base.primalRay;
964  dualFarkas = base.dualFarkas;
965  leaveCount = base.leaveCount;
966  enterCount = base.enterCount;
967  theCumulativeTime = base.theCumulativeTime;
968  infeasibilities = base.infeasibilities;
969  infeasibilitiesCo = base.infeasibilitiesCo;
970  isInfeasible = base.isInfeasible;
971  isInfeasibleCo = base.isInfeasibleCo;
972  sparsePricingLeave = base.sparsePricingLeave;
973  sparsePricingEnter = base.sparsePricingEnter;
974  sparsePricingEnterCo = base.sparsePricingEnterCo;
975  remainingRoundsLeave = base.remainingRoundsLeave;
976  remainingRoundsEnter = base.remainingRoundsEnter;
977  remainingRoundsEnterCo = base.remainingRoundsEnterCo;
978  sparsityThresholdLeave = base.sparsityThresholdLeave;
979  sparsityThresholdEnter = base.sparsityThresholdEnter;
980  sparsityThresholdEnterCo = base.sparsityThresholdEnterCo;
981 
982  if (base.theRep == COLUMN)
983  {
984  thevectors = colSet();
985  thecovectors = rowSet();
986  theFrhs = &primRhs;
987  theFvec = &primVec;
988  theCoPrhs = &dualRhs;
989  theCoPvec = &dualVec;
990  thePvec = &addVec;
991  theRPvec = theCoPvec;
992  theCPvec = thePvec;
993  theUbound = &theUCbound;
994  theLbound = &theLCbound;
995  theCoUbound = &theURbound;
996  theCoLbound = &theLRbound;
997  }
998  else
999  {
1000  assert(base.theRep == ROW);
1001 
1002  thevectors = rowSet();
1003  thecovectors = colSet();
1004  theFrhs = &dualRhs;
1005  theFvec = &dualVec;
1006  theCoPrhs = &primRhs;
1007  theCoPvec = &primVec;
1008  thePvec = &addVec;
1009  theRPvec = thePvec;
1010  theCPvec = theCoPvec;
1011  theUbound = &theURbound;
1012  theLbound = &theLRbound;
1013  theCoUbound = &theUCbound;
1014  theCoLbound = &theLCbound;
1015  }
1016 
1017  SPxBasis::theLP = this;
1018 
1019  assert(!freePricer || thepricer != 0);
1020  assert(!freeRatioTester || theratiotester != 0);
1021  assert(!freeStarter || thestarter != 0);
1022 
1023  // thepricer
1024  if(freePricer)
1025  {
1026  delete thepricer;
1027  thepricer = 0;
1028  }
1029  if(base.thepricer == 0)
1030  {
1031  thepricer = 0;
1032  freePricer = false;
1033  }
1034  else
1035  {
1036  thepricer = base.thepricer->clone();
1037  freePricer = true;
1038  thepricer->load(this);
1039  }
1040 
1041  // theratiotester
1042  if(freeRatioTester)
1043  {
1044  delete theratiotester;
1045  theratiotester = 0;
1046  }
1047  if(base.theratiotester == 0)
1048  {
1049  theratiotester = 0;
1050  freeRatioTester = false;
1051  }
1052  else
1053  {
1054  theratiotester = base.theratiotester->clone();
1055  freeRatioTester = true;
1056  theratiotester->load(this);
1057  }
1058 
1059  // thestarter
1060  if(freeStarter)
1061  {
1062  delete thestarter;
1063  thestarter = 0;
1064  }
1065  if(base.thestarter == 0)
1066  {
1067  thestarter = 0;
1068  freeStarter = false;
1069  }
1070  else
1071  {
1072  thestarter = base.thestarter->clone();
1073  freeStarter = true;
1074  }
1075 
1076  assert(SPxSolver::isConsistent());
1077  }
1078 
1079  return *this;
1080 }
1081 
1082 
1083 SPxSolver::SPxSolver(const SPxSolver& base)
1084  : SPxLP (base)
1085  , SPxBasis(base)
1086  , theType(base.theType)
1087  , thePricing(base.thePricing)
1088  , theRep(base.theRep)
1089  , theTime(base.theTime)
1090  , theCumulativeTime(base.theCumulativeTime)
1091  , maxIters(base.maxIters)
1092  , maxRefines(base.maxRefines)
1093  , maxTime(base.maxTime)
1094  , objLimit(base.objLimit)
1095  , m_status(base.m_status)
1096  , m_entertol(base.m_entertol)
1097  , m_leavetol(base.m_leavetol)
1098  , m_irthreshold(base.m_irthreshold)
1099  , theShift(base.theShift)
1100  , lastShift(base.lastShift)
1101  , m_maxCycle(base.m_maxCycle)
1102  , m_numCycle(base.m_numCycle)
1103  , initialized(base.initialized)
1104  , solveVector2 (0)
1105  , solveVector3 (0)
1106  , coSolveVector2(0)
1107  , instableLeaveNum(base.instableLeaveNum)
1108  , instableLeave(base.instableLeave)
1109  , unitVecs(base.unitVecs)
1110  , primRhs(base.primRhs)
1111  , primVec(base.primVec)
1112  , dualRhs(base.dualRhs)
1113  , dualVec(base.dualVec)
1114  , addVec(base.addVec)
1115  , theURbound(base.theURbound)
1116  , theLRbound(base.theLRbound)
1117  , theUCbound(base.theUCbound)
1118  , theLCbound(base.theLCbound)
1119  , theUBbound(base.theUBbound)
1120  , theLBbound(base.theLBbound)
1121  , theCoTest(base.theCoTest)
1122  , theTest(base.theTest)
1123  , primalRay(base.primalRay)
1124  , dualFarkas(base.dualFarkas)
1125  , leaveCount(base.leaveCount)
1126  , enterCount(base.enterCount)
1127  , infeasibilities(base.infeasibilities)
1128  , infeasibilitiesCo(base.infeasibilitiesCo)
1129  , isInfeasible(base.isInfeasible)
1130  , isInfeasibleCo(base.isInfeasibleCo)
1131  , sparsePricingLeave(base.sparsePricingLeave)
1132  , sparsePricingEnter(base.sparsePricingEnter)
1133  , sparsePricingEnterCo(base.sparsePricingEnterCo)
1134  , remainingRoundsLeave(base.remainingRoundsLeave)
1135  , remainingRoundsEnter(base.remainingRoundsEnter)
1136  , remainingRoundsEnterCo(base.remainingRoundsEnterCo)
1137  , sparsityThresholdLeave(base.sparsityThresholdLeave)
1138  , sparsityThresholdEnter(base.sparsityThresholdEnter)
1139  , sparsityThresholdEnterCo(base.sparsityThresholdEnterCo)
1140 {
1141  METHOD( "SPxSolver::SPxSolver(const SPxSolver&base)" );
1142 
1143  if (base.theRep == COLUMN)
1144  {
1145  thevectors = colSet();
1146  thecovectors = rowSet();
1147  theFrhs = &primRhs;
1148  theFvec = &primVec;
1149  theCoPrhs = &dualRhs;
1150  theCoPvec = &dualVec;
1151  thePvec = &addVec;
1152  theRPvec = theCoPvec;
1153  theCPvec = thePvec;
1154  theUbound = &theUCbound;
1155  theLbound = &theLCbound;
1156  theCoUbound = &theURbound;
1157  theCoLbound = &theLRbound;
1158  }
1159  else
1160  {
1161  assert(base.theRep == ROW);
1162 
1163  thevectors = rowSet();
1164  thecovectors = colSet();
1165  theFrhs = &dualRhs;
1166  theFvec = &dualVec;
1167  theCoPrhs = &primRhs;
1168  theCoPvec = &primVec;
1169  thePvec = &addVec;
1170  theRPvec = thePvec;
1171  theCPvec = theCoPvec;
1172  theUbound = &theURbound;
1173  theLbound = &theLRbound;
1174  theCoUbound = &theUCbound;
1175  theCoLbound = &theLCbound;
1176  }
1177 
1178  SPxBasis::theLP = this;
1179 
1180  if(base.thepricer == 0)
1181  {
1182  thepricer = 0;
1183  freePricer = false;
1184  }
1185  else
1186  {
1187  thepricer = base.thepricer->clone();
1188  freePricer = true;
1189  thepricer->clear();
1190  thepricer->load(this);
1191  }
1192 
1193  if(base.theratiotester == 0)
1194  {
1195  theratiotester = 0;
1196  freeRatioTester = false;
1197  }
1198  else
1199  {
1200  theratiotester = base.theratiotester->clone();
1201  freeRatioTester = true;
1202  theratiotester->clear();
1203  theratiotester->load(this);
1204  }
1205 
1206  if(base.thestarter == 0)
1207  {
1208  thestarter = 0;
1209  freeStarter = false;
1210  }
1211  else
1212  {
1213  thestarter = base.thestarter->clone();
1214  freeStarter = true;
1215  }
1216 
1217  assert(SPxSolver::isConsistent());
1218 }
1219 
1220 bool SPxSolver::isConsistent() const
1221 {
1222 #ifdef ENABLE_CONSISTENCY_CHECKS
1223  METHOD( "SPxSolver::isConsistent()" );
1224  if (epsilon() < 0)
1225  return MSGinconsistent("SPxSolver");
1226 
1227  if (primVec.delta().getEpsilon() != dualVec.delta().getEpsilon())
1228  return MSGinconsistent("SPxSolver");
1229 
1230  if (dualVec.delta().getEpsilon() != addVec.delta().getEpsilon())
1231  return MSGinconsistent("SPxSolver");
1232 
1233  if (unitVecs.size() < SPxLP::nCols() || unitVecs.size() < SPxLP::nRows())
1234  return MSGinconsistent("SPxSolver");
1235 
1236  if (initialized)
1237  {
1238  if (theFrhs->dim() != dim())
1239  return MSGinconsistent("SPxSolver");
1240  if (theFvec->dim() != dim())
1241  return MSGinconsistent("SPxSolver");
1242 
1243  if (theCoPrhs->dim() != dim())
1244  return MSGinconsistent("SPxSolver");
1245  if (thePvec->dim() != coDim())
1246  return MSGinconsistent("SPxSolver");
1247  if (theCoPvec->dim() != dim())
1248  return MSGinconsistent("SPxSolver");
1249 
1250  if (theTest.dim() != coDim())
1251  return MSGinconsistent("SPxSolver");
1252  if (theCoTest.dim() != dim())
1253  return MSGinconsistent("SPxSolver");
1254 
1255  if (theURbound.dim() != SPxLP::nRows())
1256  return MSGinconsistent("SPxSolver");
1257  if (theLRbound.dim() != SPxLP::nRows())
1258  return MSGinconsistent("SPxSolver");
1259  if (theUCbound.dim() != SPxLP::nCols())
1260  return MSGinconsistent("SPxSolver");
1261  if (theLCbound.dim() != SPxLP::nCols())
1262  return MSGinconsistent("SPxSolver");
1263  if (theUBbound.dim() != dim())
1264  return MSGinconsistent("SPxSolver");
1265  if (theLBbound.dim() != dim())
1266  return MSGinconsistent("SPxSolver");
1267  }
1268 
1269  if (rep() == COLUMN)
1270  {
1271  if(thecovectors !=
1272  reinterpret_cast<const SVSet*>(static_cast<const LPRowSet*>(this))
1273  || thevectors !=
1274  reinterpret_cast<const SVSet*>(static_cast<const LPColSet*>(this))
1275  || theFrhs != &primRhs ||
1276  theFvec != &primVec ||
1277  theCoPrhs != &dualRhs ||
1278  theCoPvec != &dualVec ||
1279  thePvec != &addVec ||
1280  theRPvec != theCoPvec ||
1281  theCPvec != thePvec ||
1282  theUbound != &theUCbound ||
1283  theLbound != &theLCbound ||
1284  theCoUbound != &theURbound ||
1285  theCoLbound != &theLRbound)
1286  return MSGinconsistent("SPxSolver");
1287  }
1288  else
1289  {
1290  if (thecovectors
1291  != reinterpret_cast<const SVSet*>(static_cast<const LPColSet*>(this))
1292  || thevectors
1293  != reinterpret_cast<const SVSet*>(static_cast<const LPRowSet*>(this))
1294  || theFrhs != &dualRhs ||
1295  theFvec != &dualVec ||
1296  theCoPrhs != &primRhs ||
1297  theCoPvec != &primVec ||
1298  thePvec != &addVec ||
1299  theRPvec != thePvec ||
1300  theCPvec != theCoPvec ||
1301  theUbound != &theURbound ||
1302  theLbound != &theLRbound ||
1303  theCoUbound != &theUCbound ||
1304  theCoLbound != &theLCbound)
1305  return MSGinconsistent("SPxSolver");
1306  }
1307 
1308  return SPxLP::isConsistent()
1309  && primRhs.isConsistent()
1310  && primVec.isConsistent()
1311  && dualRhs.isConsistent()
1312  && dualVec.isConsistent()
1313  && addVec.isConsistent()
1314  && theTest.isConsistent()
1315  && theCoTest.isConsistent()
1316  && theURbound.isConsistent()
1317  && theLRbound.isConsistent()
1318  && theUCbound.isConsistent()
1319  && theLCbound.isConsistent()
1320  && SPxBasis::isConsistent()
1321  ;
1322 #else
1323  return true;
1324 #endif
1325 }
1326 
1327 
1328 void SPxSolver::setTerminationTime(Real p_time)
1329 {
1330  METHOD( "SPxSolver::setTerminationTime()" );
1331  if( p_time < 0.0 )
1332  p_time = infinity;
1333  maxTime = p_time;
1334 }
1335 
1336 Real SPxSolver::terminationTime() const
1337 {
1338  METHOD( "SPxSolver::terminationTime()" );
1339  return maxTime;
1340 }
1341 
1342 void SPxSolver::setTerminationIter(int p_iteration)
1343 {
1344  METHOD( "SPxSolver::setTerminationIter()" );
1345  if( p_iteration < 0 )
1346  p_iteration = -1;
1347  maxIters = p_iteration;
1348 }
1349 
1350 int SPxSolver::terminationIter() const
1351 {
1352  METHOD( "SPxSolver::terminationIter()" );
1353  return maxIters;
1354 }
1355 
1356 void SPxSolver::setMaxRefinements(int p_maxrefinements)
1357 {
1358  METHOD( "SPxSolver::setMaxRefinements()" );
1359  if( p_maxrefinements < 0 )
1360  p_maxrefinements = -1;
1361  maxRefines = p_maxrefinements;
1362 }
1363 
1364 int SPxSolver::maxRefinements() const
1365 {
1366  METHOD( "SPxSolver::terminationIter()" );
1367  return maxRefines;
1368 }
1369 
1370 /**@todo A first version for the termination value is
1371  * implemented. Currently we check if no bound violations (shifting)
1372  * is present. It might be even possible to use this termination
1373  * value in case of bound violations (shifting) but in this case it
1374  * is quite difficult to determine if we already reached the limit.
1375  */
1376 void SPxSolver::setTerminationValue(Real p_value)
1377 {
1378  METHOD( "SPxSolver::setTerminationValue()" );
1379  objLimit = p_value;
1380 }
1381 
1382 Real SPxSolver::terminationValue() const
1383 {
1384  METHOD( "SPxSolver::terminationValue()" );
1385  return objLimit;
1386 }
1387 
1388 SPxSolver::VarStatus
1389 SPxSolver::basisStatusToVarStatus( SPxBasis::Desc::Status stat ) const
1390 {
1391  METHOD( "SPxSolver::VarStatus()" );
1392  VarStatus vstat;
1393 
1394  switch( stat )
1395  {
1396  case SPxBasis::Desc::P_ON_LOWER:
1397  vstat = ON_LOWER;
1398  break;
1399  case SPxBasis::Desc::P_ON_UPPER:
1400  vstat = ON_UPPER;
1401  break;
1402  case SPxBasis::Desc::P_FIXED:
1403  vstat = FIXED;
1404  break;
1405  case SPxBasis::Desc::P_FREE:
1406  vstat = ZERO;
1407  break;
1408  case SPxBasis::Desc::D_ON_UPPER:
1409  case SPxBasis::Desc::D_ON_LOWER:
1410  case SPxBasis::Desc::D_ON_BOTH:
1411  case SPxBasis::Desc::D_UNDEFINED:
1412  case SPxBasis::Desc::D_FREE:
1413  vstat = BASIC;
1414  break;
1415  default:
1416  MSG_ERROR( spxout << "ESOLVE26 ERROR: unknown basis status (" << stat << ")"
1417  << std::endl; )
1418  throw SPxInternalCodeException("XSOLVE22 This should never happen.");
1419  }
1420  return vstat;
1421 }
1422 
1423 SPxBasis::Desc::Status
1424 SPxSolver::varStatusToBasisStatusRow( int row, SPxSolver::VarStatus stat ) const
1425 {
1426  METHOD( "SPxSolver::varStatusToBasisStatusRow()" );
1427  SPxBasis::Desc::Status rstat;
1428 
1429  switch( stat )
1430  {
1431  case FIXED :
1432  assert(rhs(row) == lhs(row));
1433  rstat = SPxBasis::Desc::P_FIXED;
1434  break;
1435  case ON_UPPER :
1436  assert(rhs(row) < infinity);
1437  rstat = lhs(row) < rhs(row)
1438  ? SPxBasis::Desc::P_ON_UPPER
1439  : SPxBasis::Desc::P_FIXED;
1440  break;
1441  case ON_LOWER :
1442  assert(lhs(row) > -infinity);
1443  rstat = lhs(row) < rhs(row)
1444  ? SPxBasis::Desc::P_ON_LOWER
1445  : SPxBasis::Desc::P_FIXED;
1446  break;
1447  case ZERO :
1448  /* A 'free' row (i.e., infinite lower & upper bounds) does not really make sense. The user
1449  * might (think to) know better, e.g., when temporarily turning off a row. We therefore apply
1450  * the same adjustment as in the column case in varStatusToBasisStatusCol(). */
1451  if (lhs(row) <= -infinity && rhs(row) >= infinity)
1452  rstat = SPxBasis::Desc::P_FREE;
1453  else
1454  {
1455  if ( lhs(row) == rhs(row) )
1456  {
1457  assert( rhs(row) < infinity );
1458  rstat = SPxBasis::Desc::P_FIXED;
1459  }
1460  else
1461  {
1462  if ( lhs(row) > -infinity )
1463  rstat = SPxBasis::Desc::P_ON_LOWER;
1464  else
1465  {
1466  assert( rhs(row) < infinity );
1467  rstat = SPxBasis::Desc::P_ON_UPPER;
1468  }
1469  }
1470  }
1471  break;
1472  case BASIC :
1473  rstat = dualRowStatus(row);
1474  break;
1475  default:
1476  MSG_ERROR( spxout << "ESOLVE27 ERROR: unknown VarStatus (" << int(stat) << ")"
1477  << std::endl; )
1478  throw SPxInternalCodeException("XSOLVE23 This should never happen.");
1479  }
1480  return rstat;
1481 }
1482 
1483 SPxBasis::Desc::Status
1484 SPxSolver::varStatusToBasisStatusCol( int col, SPxSolver::VarStatus stat ) const
1485 {
1486  METHOD( "SPxSolver::varStatusToBasisStatusCol()" );
1487  SPxBasis::Desc::Status cstat;
1488 
1489  switch( stat )
1490  {
1491  case FIXED :
1492  if (upper(col) == lower(col))
1493  cstat = SPxBasis::Desc::P_FIXED;
1494  else if (maxObj(col) > 0.0)
1495  cstat = SPxBasis::Desc::P_ON_UPPER;
1496  else
1497  cstat = SPxBasis::Desc::P_ON_LOWER;
1498  break;
1499  case ON_UPPER :
1500  assert(upper(col) < infinity);
1501  cstat = lower(col) < upper(col)
1502  ? SPxBasis::Desc::P_ON_UPPER
1503  : SPxBasis::Desc::P_FIXED;
1504  break;
1505  case ON_LOWER :
1506  assert(lower(col) > -infinity);
1507  cstat = lower(col) < upper(col)
1508  ? SPxBasis::Desc::P_ON_LOWER
1509  : SPxBasis::Desc::P_FIXED;
1510  break;
1511  case ZERO :
1512  /* In this case the upper and lower bounds on the variable should be infinite. The bounds
1513  * might, however, have changed and we try to recover from this by changing the status to
1514  * 'resonable' settings. A solve has to find the correct values afterwards. Note that the
1515  * approach below is consistent with changesoplex.cpp (e.g., changeUpperStatus() and
1516  * changeLowerStatus() ). */
1517  if (lower(col) <= -infinity && upper(col) >= infinity)
1518  cstat = SPxBasis::Desc::P_FREE;
1519  else
1520  {
1521  if ( lower(col) == upper(col) )
1522  {
1523  assert( upper(col) < infinity );
1524  cstat = SPxBasis::Desc::P_FIXED;
1525  }
1526  else
1527  {
1528  if ( lower(col) > -infinity )
1529  cstat = SPxBasis::Desc::P_ON_LOWER;
1530  else
1531  {
1532  assert( upper(col) < infinity );
1533  cstat = SPxBasis::Desc::P_ON_UPPER;
1534  }
1535  }
1536  }
1537  break;
1538  case BASIC :
1539  cstat = dualColStatus(col);
1540  break;
1541  default:
1542  MSG_ERROR( spxout << "ESOLVE28 ERROR: unknown VarStatus (" << int(stat) << ")"
1543  << std::endl; )
1544  throw SPxInternalCodeException("XSOLVE24 This should never happen.");
1545  }
1546  return cstat;
1547 }
1548 
1549 SPxSolver::VarStatus SPxSolver::getBasisRowStatus( int row ) const
1550 {
1551  METHOD( "SPxSolver::VarStatus()" );
1552  assert( 0 <= row && row < nRows() );
1553  return basisStatusToVarStatus( desc().rowStatus( row ) );
1554 }
1555 
1556 SPxSolver::VarStatus SPxSolver::getBasisColStatus( int col ) const
1557 {
1558  METHOD( "SPxSolver::VarStatus()" );
1559  assert( 0 <= col && col < nCols() );
1560  return basisStatusToVarStatus( desc().colStatus( col ) );
1561 }
1562 
1563 SPxSolver::Status SPxSolver::getBasis(VarStatus row[], VarStatus col[]) const
1564 {
1565  METHOD( "SPxSolver::Status()" );
1566  const SPxBasis::Desc& d = desc();
1567  int i;
1568 
1569  if (col)
1570  for (i = nCols() - 1; i >= 0; --i)
1571  col[i] = basisStatusToVarStatus( d.colStatus(i) );
1572 
1573  if (row)
1574  for (i = nRows() - 1; i >= 0; --i)
1575  row[i] = basisStatusToVarStatus( d.rowStatus(i) );
1576 
1577  return status();
1578 }
1579 
1580 bool SPxSolver::isBasisValid(DataArray<VarStatus> p_rows, DataArray<VarStatus> p_cols)
1581 {
1582  METHOD( "SPxSolver::isBasisValid()" );
1583 
1584  int basisdim;
1585 
1586  if ( p_rows.size() != nRows() || p_cols.size() != nCols() )
1587  return false;
1588 
1589  basisdim = 0;
1590  for ( int row = nRows()-1; row >= 0; --row )
1591  {
1592  if ( p_rows[row] == UNDEFINED )
1593  return false;
1594  // row is basic
1595  else if ( p_rows[row] == BASIC )
1596  {
1597  basisdim++;
1598  }
1599  // row is nonbasic
1600  else
1601  {
1602  if ( (p_rows[row] == FIXED && lhs(row) != rhs(row))
1603  || (p_rows[row] == ON_UPPER && rhs(row) >= infinity)
1604  || (p_rows[row] == ON_LOWER && lhs(row) <= -infinity) )
1605  return false;
1606  }
1607  }
1608 
1609  for ( int col = nCols()-1; col >= 0; --col )
1610  {
1611  if ( p_cols[col] == UNDEFINED )
1612  return false;
1613  // col is basic
1614  else if ( p_cols[col] == BASIC )
1615  {
1616  basisdim++;
1617  }
1618  // col is nonbasic
1619  else
1620  {
1621  if ( (p_cols[col] == FIXED && lower(col) != upper(col))
1622  || (p_cols[col] == ON_UPPER && upper(col) >= infinity)
1623  || (p_cols[col] == ON_LOWER && lower(col) <= -infinity) )
1624  return false;
1625  }
1626  }
1627 
1628  if ( basisdim != dim() )
1629  return false;
1630 
1631  // basis valid
1632  return true;
1633 }
1634 
1635 void SPxSolver::setBasis(const VarStatus p_rows[], const VarStatus p_cols[])
1636 {
1637  METHOD( "SPxSolver::setBasis()" );
1638 
1639  if (SPxBasis::status() == SPxBasis::NO_PROBLEM)
1640  SPxBasis::load(this);
1641 
1642  SPxBasis::Desc ds = desc();
1643  int i;
1644 
1645  for(i = 0; i < nRows(); i++)
1646  ds.rowStatus(i) = varStatusToBasisStatusRow( i, p_rows[i] );
1647 
1648  for(i = 0; i < nCols(); i++)
1649  ds.colStatus(i) = varStatusToBasisStatusCol( i, p_cols[i] );
1650 
1651  loadBasis(ds);
1652 }
1653 
1654 //
1655 // Auxiliary functions.
1656 //
1657 
1658 // Pretty-printing of variable status.
1659 std::ostream& operator<<( std::ostream& os,
1660  const SPxSolver::VarStatus& status )
1661 {
1662  switch( status )
1663  {
1664  case SPxSolver::BASIC:
1665  os << "BASIC";
1666  break;
1667  case SPxSolver::FIXED:
1668  os << "FIXED";
1669  break;
1670  case SPxSolver::ON_LOWER:
1671  os << "ON_LOWER";
1672  break;
1673  case SPxSolver::ON_UPPER:
1674  os << "ON_UPPER";
1675  break;
1676  case SPxSolver::ZERO:
1677  os << "ZERO";
1678  break;
1679  case SPxSolver::UNDEFINED:
1680  os << "UNDEFINED";
1681  break;
1682  default:
1683  os << "?invalid?";
1684  break;
1685  }
1686  return os;
1687 }
1688 
1689 // Pretty-printing of solver status.
1690 std::ostream& operator<<( std::ostream& os,
1691  const SPxSolver::Status& status )
1692 {
1693  switch ( status )
1694  {
1695  case SPxSolver::ERROR:
1696  os << "ERROR";
1697  break;
1698  case SPxSolver::NO_RATIOTESTER:
1699  os << "NO_RATIOTESTER";
1700  break;
1701  case SPxSolver::NO_PRICER:
1702  os << "NO_PRICER";
1703  break;
1704  case SPxSolver::NO_SOLVER:
1705  os << "NO_SOLVER";
1706  break;
1707  case SPxSolver::NOT_INIT:
1708  os << "NOT_INIT";
1709  break;
1710  case SPxSolver::ABORT_CYCLING:
1711  os << "ABORT_CYCLING";
1712  break;
1713  case SPxSolver::ABORT_TIME:
1714  os << "ABORT_TIME";
1715  break;
1716  case SPxSolver::ABORT_ITER:
1717  os << "ABORT_ITER";
1718  break;
1719  case SPxSolver::ABORT_VALUE:
1720  os << "ABORT_VALUE";
1721  break;
1722  case SPxSolver::SINGULAR:
1723  os << "SINGULAR";
1724  break;
1725  case SPxSolver::NO_PROBLEM:
1726  os << "NO_PROBLEM";
1727  break;
1728  case SPxSolver::REGULAR:
1729  os << "REGULAR";
1730  break;
1731  case SPxSolver::RUNNING:
1732  os << "RUNNING";
1733  break;
1734  case SPxSolver::UNKNOWN:
1735  os << "UNKNOWN";
1736  break;
1737  case SPxSolver::OPTIMAL:
1738  os << "OPTIMAL";
1739  break;
1740  case SPxSolver::UNBOUNDED:
1741  os << "UNBOUNDED";
1742  break;
1743  case SPxSolver::INFEASIBLE:
1744  os << "INFEASIBLE";
1745  break;
1746  default:
1747  os << "?other?";
1748  break;
1749  }
1750  return os;
1751 }
1752 
1753 // Pretty-printing of algorithm.
1754 std::ostream& operator<<( std::ostream& os,
1755  const SPxSolver::Type& status )
1756 {
1757  switch ( status )
1758  {
1759  case SPxSolver::ENTER:
1760  os << "ENTER";
1761  break;
1762  case SPxSolver::LEAVE:
1763  os << "LEAVE";
1764  break;
1765  default:
1766  os << "?other?";
1767  break;
1768  }
1769  return os;
1770 }
1771 
1772 // Pretty-printing of representation.
1773 std::ostream& operator<<( std::ostream& os,
1774  const SPxSolver::Representation& status )
1775 {
1776  switch ( status )
1777  {
1778  case SPxSolver::ROW:
1779  os << "ROW";
1780  break;
1781  case SPxSolver::COLUMN:
1782  os << "COLUMN";
1783  break;
1784  default:
1785  os << "?other?";
1786  break;
1787  }
1788  return os;
1789 }
1790 
1791 
1792 } // namespace soplex
1793 
1794 //-----------------------------------------------------------------------------
1795 //Emacs Local Variables:
1796 //Emacs mode:c++
1797 //Emacs c-basic-offset:3
1798 //Emacs tab-width:8
1799 //Emacs indent-tabs-mode:nil
1800 //Emacs End:
1801 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen