SoPlex Doxygen Documentation
soplexmain.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 #include <assert.h>
17 #include <math.h>
18 #include <string.h>
19 #include <iostream>
20 #include <iomanip>
21 #include <fstream>
22 
23 #include "spxdefines.h"
24 #include "soplex.h"
25 #include "spxsolver.h"
26 
27 #include "timer.h"
28 #include "spxgithash.h"
29 #include "spxpricer.h"
30 #include "spxdantzigpr.h"
31 #include "spxparmultpr.h"
32 #include "spxdevexpr.h"
33 #include "spxhybridpr.h"
34 #include "spxsteeppr.h"
35 #include "spxsteepexpr.h"
36 #include "spxweightpr.h"
37 #include "spxratiotester.h"
38 #include "spxharrisrt.h"
39 #include "spxdefaultrt.h"
40 #include "spxfastrt.h"
41 #include "spxboundflippingrt.h"
42 #include "spxsimplifier.h"
43 #include "spxmainsm.h"
44 #include "spxscaler.h"
45 #include "spxequilisc.h"
46 #include "spxgeometsc.h"
47 #include "spxsumst.h"
48 #include "spxweightst.h"
49 #include "spxvectorst.h"
50 #include "slufactor.h"
51 #include "spxout.h"
52 
53 using namespace soplex;
54 
55 
56 //------------------------------------------------------------------------
57 // for simplicity: store whether we are in check mode:
58 static bool checkMode = false;
59 //------------------------------------------------------------------------
60 
61 
62 //------------------------------------------------------------------------
63 // class MySoPlex
64 //------------------------------------------------------------------------
65 
66 /** LP solver class for the command line. */
67 class MySoPlex : public SoPlex
68 {
69 public:
70  /// default constructor
71  MySoPlex( SPxSolver::Type p_type = SPxSolver::LEAVE,
72  SPxSolver::Representation p_rep = SPxSolver::COLUMN )
73  : SoPlex(p_type, p_rep)
74  {}
75  //------------------------------------------------------------------------
76  /// virtual destructor
77  virtual ~MySoPlex()
78  {}
79  //------------------------------------------------------------------------
80  void displayQuality() const
81  {
82  Real maxviol;
83  Real sumviol;
84 
85  if ( checkMode )
86  {
87  MSG_INFO1( spxout << "IEXAMP05 Violations (max/sum)" << std::endl; )
88 
89  m_solver.qualConstraintViolation(maxviol, sumviol);
90 
91  MSG_INFO1( spxout << "IEXAMP06 Constraints :"
92  << std::setw(16) << maxviol << " "
93  << std::setw(16) << sumviol << std::endl; )
94 
95  qualConstraintViolation(maxviol, sumviol);
96 
97  MSG_INFO1( spxout << "IEXAMP07 (unscaled) :"
98  << std::setw(16) << maxviol << " "
99  << std::setw(16) << sumviol << std::endl; )
100 
101  m_solver.qualBoundViolation(maxviol, sumviol);
102 
103  MSG_INFO1( spxout << "IEXAMP08 Bounds :"
104  << std::setw(16) << maxviol << " "
105  << std::setw(16) << sumviol << std::endl; )
106 
107  qualBoundViolation(maxviol, sumviol);
108 
109  MSG_INFO1( spxout << "IEXAMP09 (unscaled) :"
110  << std::setw(16) << maxviol << " "
111  << std::setw(16) << sumviol << std::endl; )
112 
113  if (!m_vanished)
114  {
115  m_solver.qualSlackViolation(maxviol, sumviol);
116 
117  MSG_INFO1( spxout << "IEXAMP10 Slacks :"
118  << std::setw(16) << maxviol << " "
119  << std::setw(16) << sumviol << std::endl; )
120 
121  m_solver.qualRedCostViolation(maxviol, sumviol);
122 
123  MSG_INFO1( spxout << "IEXAMP11 Reduced costs :"
124  << std::setw(16) << maxviol << " "
125  << std::setw(16) << sumviol << std::endl; )
126 #if 0
127  MSG_INFO1( spxout << "IEXAMP12 Proven dual bound:"
128  << std::setw(20)
129  << std::setprecision(20)
130  << m_solver.provedDualbound() << std::endl; )
131 #endif
132  }
133  }
134  else
135  {
136  MSG_INFO1( spxout << "Violations (max/sum)" << std::endl; )
137 
138  m_solver.qualConstraintViolation(maxviol, sumviol);
139 
140  MSG_INFO1( spxout << "Constraints :"
141  << std::setw(16) << maxviol << " "
142  << std::setw(16) << sumviol << std::endl; )
143 
144  qualConstraintViolation(maxviol, sumviol);
145 
146  MSG_INFO1( spxout << " (unscaled) :"
147  << std::setw(16) << maxviol << " "
148  << std::setw(16) << sumviol << std::endl; )
149 
150  m_solver.qualBoundViolation(maxviol, sumviol);
151 
152  MSG_INFO1( spxout << "Bounds :"
153  << std::setw(16) << maxviol << " "
154  << std::setw(16) << sumviol << std::endl; )
155 
156  qualBoundViolation(maxviol, sumviol);
157 
158  MSG_INFO1( spxout << " (unscaled) :"
159  << std::setw(16) << maxviol << " "
160  << std::setw(16) << sumviol << std::endl; )
161 
162  if (!m_vanished)
163  {
164  m_solver.qualSlackViolation(maxviol, sumviol);
165 
166  MSG_INFO1( spxout << "Slacks :"
167  << std::setw(16) << maxviol << " "
168  << std::setw(16) << sumviol << std::endl; )
169 
170  m_solver.qualRedCostViolation(maxviol, sumviol);
171 
172  MSG_INFO1( spxout << "Reduced costs :"
173  << std::setw(16) << maxviol << " "
174  << std::setw(16) << sumviol << std::endl; )
175 #if 0
176  MSG_INFO1( spxout << "Proven dual bound:"
177  << std::setw(20)
178  << std::setprecision(20)
179  << m_solver.provedDualbound() << std::endl; )
180 #endif
181  }
182  }
183  }
184  //------------------------------------------------------------------------
185  void displayInfeasibility() const
186  {
187  assert(m_solver.status() == SPxSolver::INFEASIBLE);
188 
189 #if 0
190  if ( checkMode )
191  {
192  if( m_solver.isProvenInfeasible() )
193  MSG_INFO1( spxout << "IEXAMP13 Infeasibility is proven." << std::endl; )
194  else
195  MSG_INFO1( spxout << "IEXAMP13 Infeasibility could not be proven!" << std::endl; )
196  }
197  else
198  {
199  if ( m_solver.isProvenInfeasible() )
200  {
201  MSG_INFO1( spxout << "Infeasibility is proven." << std::endl; )
202  }
203  else
204  {
205  MSG_INFO1( spxout << "Infeasibility could not be proven!" << std::endl; )
206  }
207  }
208 #endif
209  }
210 };
211 
212 
213 //------------------------------------------------------------------------
214 // Helpers
215 //------------------------------------------------------------------------
216 
217 static
218 void print_version_info()
219 {
220  const char* banner1 =
221  "************************************************************************\n"
222  "* *\n"
223  "* SoPlex --- the Sequential object-oriented simPlex. *\n"
224  ;
225 
226  const char* banner2 =
227  "* *\n"
228  "* Copyright (C) 1996-2012 Konrad-Zuse-Zentrum *\n"
229  "* fuer Informationstechnik Berlin *\n"
230  "* *\n"
231  "* SoPlex is distributed under the terms of the ZIB Academic Licence. *\n"
232  "* You should have received a copy of the ZIB Academic License *\n"
233  "* along with SoPlex; If not email to soplex@zib.de. *\n"
234  "* *\n"
235  "************************************************************************\n"
236  ;
237 
238  if( !checkMode )
239  std::cout << banner1;
240 
241 #if (SOPLEX_SUBVERSION > 0)
242  if( !checkMode )
243  std::cout << "* Version ";
244  else
245  std::cout << "SoPlex version ";
246  std::cout << SOPLEX_VERSION/100 << "."
247  << (SOPLEX_VERSION % 100)/10 << "."
248  << SOPLEX_VERSION % 10 << "."
249  << SOPLEX_SUBVERSION
250  << " - Githash "
251  << std::setw(13) << std::setiosflags(std::ios::left) << getGitHash();
252  if( !checkMode )
253  std::cout << " *\n" << banner2 << std::endl;
254  else
255  std::cout << "\n";
256 #else
257  if( !checkMode )
258  std::cout << "* Release ";
259  else
260  std::cout << "SoPlex release ";
261  std::cout << SOPLEX_VERSION/100 << "."
262  << (SOPLEX_VERSION % 100)/10 << "."
263  << SOPLEX_VERSION % 10
264  << " - Githash "
265  << std::setw(13) << std::setiosflags(std::ios::left) << getGitHash();
266  if( !checkMode )
267  std::cout << " *\n" << banner2 << std::endl;
268  else
269  std::cout << "\n";
270 #endif
271 
272  /// The following code block is tests and shows compilation parameters.
273  std::cout << "[NDEBUG:"
274 #ifdef NDEBUG
275  << "YES"
276 #else
277  << "NO"
278 #endif
279  << "]";
280 
281  std::cout << "[WITH_WARNINGS:"
282 #ifdef WITH_WARNINGS
283  << "YES"
284 #else
285  << "NO"
286 #endif
287  << "]";
288 
289  std::cout << "[ENABLE_ADDITIONAL_CHECKS:"
290 #ifdef ENABLE_ADDITIONAL_CHECKS
291  << "YES"
292 #else
293  << "NO"
294 #endif
295  << "]";
296 
297  std::cout << "[ENABLE_CONSISTENCY_CHECKS:"
298 #ifdef ENABLE_CONSISTENCY_CHECKS
299  << "YES"
300 #else
301  << "NO"
302 #endif
303  << "]";
304 
305  std::cout << "[SOPLEX_WITH_GMP:"
306 #ifdef SOPLEX_WITH_GMP
307  << "YES"
308 #else
309  << "NO"
310 #endif
311  << "]" << std::endl;
312 
313  std::cout << std::endl;
314 }
315 
316 #if 0
317 static
318 void print_short_version_info()
319 {
320  const char* banner1 =
321  "************************************************************************\n"
322  "* SoPlex --- the Sequential object-oriented simPlex. ";
323  const char* banner2 =
324  "* Copyright (C) 1996-2012 Zuse Institute Berlin *\n"
325  "************************************************************************\n";
326 
327  std::cout << banner1;
328 #if (SOPLEX_SUBVERSION > 0)
329  std::cout << "Version "
330  << SOPLEX_VERSION/100 << "."
331  << (SOPLEX_VERSION % 100)/10 << "."
332  << SOPLEX_VERSION % 10 << "."
333  << SOPLEX_SUBVERSION
334  << " *\n";
335 #else
336  std::cout << "Release "
337  << SOPLEX_VERSION/100 << "."
338  << (SOPLEX_VERSION % 100)/10 << "."
339  << SOPLEX_VERSION % 10
340  << " *\n";
341 #endif
342  std::cout << banner2 << std::endl;
343 }
344 #endif
345 
346 //------------------------------------------------------------------------
347 static
348 void print_usage_and_exit( const char* const argv[] )
349 {
350  const char* usage =
351  "[options] LPfile [Basfile]\n\n"
352  " LPfile can be either in MPS or LPF format\n\n"
353  "options: (*) indicates default\n"
354  " (!) indicates experimental features which may give wrong results\n"
355  " -e select entering algorithm (default is leaving)\n"
356  " -r select row wise representation (default is column)\n"
357  " -i select Eta-update (default is Forest-Tomlin)\n"
358  " -x output solution vector\n"
359  " -y output dual multipliers\n"
360  " -q display solution quality\n"
361  " -br read file with starting basis from Basfile\n"
362  " -bw write file with optimal basis to Basfile\n"
363  " -l set time limit in seconds\n"
364  " -L set iteration limit\n"
365  " -f set primal feasibility tolerance\n"
366  " -o set optimality, i.e., dual feasibility tolerance\n"
367  " -d set primal and dual feasibility tolerance to same value\n"
368  " -R set threshold for tolerances below which iterative refinement is applied\n"
369  " -zz set general zero tolerance\n"
370  " -zf set factorization zero tolerance\n"
371  " -zu set update zero tolerance\n"
372  " -v set verbosity Level: from 0 (ERROR) to 5 (INFO3), default 3 (INFO1)\n"
373  " -V show program version\n"
374  " -C check mode (for check scripts)\n"
375  " -h show this help\n\n"
376  "Simplifier: Scaler: Starter: Pricer: Ratiotester:\n"
377  " -s0 none -g0 none -c0 none* -p0 Textbook -t0 Textbook\n"
378  " -s1 Main* -g1 uni-Equi -c1 Weight -p1 ParMult -t1 Harris\n"
379  " -g2 bi-Equi* -c2 Sum -p2 Devex -t2 Fast*\n"
380  " -g3 bi-Equi+Geo1 -c3 Vector -p3 Hybrid! -t3 Bound Flipping\n"
381  " -g4 bi-Equi+Geo8 -p4 Steep*\n"
382  " -p5 Weight\n"
383  " -p6 SteepExactSetup\n"
384  ;
385 
386  std::cerr << "usage: " << argv[0] << " " << usage << std::endl;
387  exit(0);
388 }
389 
390 //------------------------------------------------------------------------
391 static
392 void check_parameter(const char param, const char* const argv[])
393 {
394  if (param == '\0')
395  print_usage_and_exit( argv );
396 }
397 
398 //------------------------------------------------------------------------
399 static
400 void print_algorithm_parameters(
401  MySoPlex& work,
402  const SPxSolver::Representation representation,
403  const SLUFactor::UpdateType update
404  )
405 {
406  if ( checkMode )
407  {
408  MSG_INFO1( spxout
409  << "IEXAMP12 Feastol = "
410  << std::setw(16) << work.feastol() << std::endl
411  << "IEXAMP52 Opttol = "
412  << std::setw(16) << work.opttol() << std::endl
413  << "IEXAMP53 Irthreshold = "
414  << std::setw(16) << work.irthreshold() << std::endl
415  << "IEXAMP13 Epsilon Zero = "
416  << std::setw(16) << Param::epsilon() << std::endl
417  << "IEXAMP37 Epsilon Factor = "
418  << std::setw(16) << Param::epsilonFactorization() << std::endl
419  << "IEXAMP38 Epsilon Update = "
420  << std::setw(16) << Param::epsilonUpdate() << std::endl
421  << "IEXAMP14 "
422  << (work.type() == SPxSolver::ENTER ? "Entering" : "Leaving")
423  << " algorithm" << std::endl
424  << "IEXAMP15 "
425  << (representation == SPxSolver::ROW ? "Row" : "Column")
426  << " representation" << std::endl
427  << "IEXAMP16 "
428  << (update == SLUFactor::ETA ? "Eta" : "Forest-Tomlin")
429  << " update" << std::endl; )
430  }
431  else
432  {
433  MSG_INFO1( spxout
434  << "SoPlex parameters: " << std::endl
435  << "Feastol = "
436  << std::setw(16) << work.feastol() << std::endl
437  << "Opttol = "
438  << std::setw(16) << work.opttol() << std::endl
439  << "Irthreshold = "
440  << std::setw(16) << work.irthreshold() << std::endl
441  << "Epsilon Zero = "
442  << std::setw(16) << Param::epsilon() << std::endl
443  << "Epsilon Factor = "
444  << std::setw(16) << Param::epsilonFactorization() << std::endl
445  << "Epsilon Update = "
446  << std::setw(16) << Param::epsilonUpdate() << std::endl
447  << std::endl
448  << "algorithm = " << (work.type() == SPxSolver::ENTER ? "Entering" : "Leaving")
449  << std::endl
450  << "representation = " << (representation == SPxSolver::ROW ? "Row" : "Column")
451  << std::endl
452  << "update = " << (update == SLUFactor::ETA ? "Eta" : "Forest-Tomlin")
453  << std::endl; )
454  }
455 }
456 
457 //------------------------------------------------------------------------
458 static
459 SPxPricer* get_pricer(const int pricing)
460 {
461  SPxPricer* pricer = 0;
462  switch(pricing)
463  {
464  case 6 :
465  pricer = new SPxSteepExPR;
466  break;
467  case 5 :
468  pricer = new SPxWeightPR;
469  break;
470  case 4 :
471  pricer = new SPxSteepPR;
472  break;
473  case 3 :
474  pricer = new SPxHybridPR;
475  break;
476  case 2 :
477  pricer = new SPxDevexPR;
478  break;
479  case 1 :
480  pricer = new SPxParMultPR;
481  break;
482  case 0 :
483  /*FALLTHROUGH*/
484  default :
485  pricer = new SPxDantzigPR;
486  break;
487  }
488 
489  assert(pricer != 0);
490  if ( checkMode )
491 #ifdef PARTIAL_PRICING
492  MSG_INFO1( spxout << "IEXAMP17 " << pricer->getName() << " pricing"
493  << " (partial, size = " << MAX_PRICING_CANDIDATES << ")"
494  << std::endl; )
495 #else
496  MSG_INFO1( spxout << "IEXAMP17 " << pricer->getName() << " pricing"
497  << std::endl; )
498 #endif
499  else
500 #ifdef PARTIAL_PRICING
501  MSG_INFO1( spxout << "pricing = " << pricer->getName()
502  << " (partial, size = " << MAX_PRICING_CANDIDATES << ")"
503  << std::endl; )
504 #else
505  MSG_INFO1( spxout << "pricing = " << pricer->getName()
506  << std::endl; )
507 #endif
508  return pricer;
509 }
510 
511 //------------------------------------------------------------------------
512 static
513 SPxRatioTester* get_ratio_tester(const int ratiotest)
514 {
515  SPxRatioTester* ratiotester = 0;
516  switch(ratiotest)
517  {
518  case 3 :
519  ratiotester = new SPxBoundFlippingRT;
520  break;
521  case 2 :
522  ratiotester = new SPxFastRT;
523  break;
524  case 1 :
525  ratiotester = new SPxHarrisRT;
526  break;
527  case 0 :
528  /*FALLTHROUGH*/
529  default:
530  ratiotester = new SPxDefaultRT;
531  break;
532  }
533 
534  assert(ratiotester != 0);
535  if ( checkMode )
536  MSG_INFO1( spxout << "IEXAMP18 " << ratiotester->getName() << " ratiotest" << std::endl; )
537  else
538  MSG_INFO1( spxout << "ratiotest = " << ratiotester->getName() << std::endl; )
539  return ratiotester;
540 }
541 
542 //------------------------------------------------------------------------
543 static
544 void get_scalers(
545  SPxScaler*& prescaler,
546  SPxScaler*& postscaler,
547  const int scaling
548  )
549 {
550  switch(scaling)
551  {
552  case 4:
553  prescaler = new SPxEquiliSC(true);
554  postscaler = new SPxGeometSC(8);
555  break;
556  case 3:
557  prescaler = new SPxEquiliSC(true);
558  postscaler = new SPxGeometSC(1);
559  break;
560  case 2 :
561  prescaler = new SPxEquiliSC(true);
562  postscaler = 0;
563  break;
564  case 1 :
565  prescaler = new SPxEquiliSC(false);
566  postscaler = 0;
567  break;
568  case 0 :
569  /*FALLTHROUGH*/
570  default :
571  prescaler = 0;
572  postscaler = 0;
573  break;
574  }
575 
576  if ( checkMode )
577  {
578  MSG_INFO1( spxout << "IEXAMP19 "
579  << ((prescaler != 0) ? prescaler->getName() : "no")
580  << " / "
581  << ((postscaler != 0) ? postscaler->getName() : "no")
582  << " scaling" << std::endl; )
583  }
584  else
585  {
586  MSG_INFO1( spxout << "scaling = "
587  << ((prescaler != 0) ? prescaler->getName() : "no")
588  << " / "
589  << ((postscaler != 0) ? postscaler->getName() : "no")
590  << std::endl; )
591  }
592 }
593 
594 //------------------------------------------------------------------------
595 static
596 SPxSimplifier* get_simplifier(const int simplifying)
597 {
598  SPxSimplifier* simplifier = 0;
599  switch(simplifying)
600  {
601  case 1 :
602  simplifier = new SPxMainSM;
603  break;
604  case 0 :
605  /*FALLTHROUGH*/
606  default :
607  assert(simplifier == 0);
608  break;
609  }
610 
611  if ( checkMode )
612  MSG_INFO1( spxout << "IEXAMP20 " << ((simplifier == 0) ? "no" : simplifier->getName()) << " simplifier" << std::endl; )
613  else
614  MSG_INFO1( spxout << "simplifier = " << ((simplifier == 0) ? "no" : simplifier->getName()) << std::endl; )
615  return simplifier;
616 }
617 
618 //------------------------------------------------------------------------
619 static
620 SPxStarter* get_starter(const int starting)
621 {
622  SPxStarter* starter = 0;
623  switch(starting)
624  {
625  case 3 :
626  starter = new SPxVectorST;
627  break;
628  case 2 :
629  starter = new SPxSumST;
630  break;
631  case 1 :
632  starter = new SPxWeightST;
633  break;
634  case 0 :
635  /*FALLTHROUGH*/
636  default :
637  break;
638  }
639 
640  if ( checkMode )
641  MSG_INFO1( spxout << "IEXAMP21 " << ((starter == 0) ? "no" : starter->getName()) << " starter" << std::endl; )
642  else
643  MSG_INFO1( spxout << "starter = " << ((starter == 0) ? "no" : starter->getName()) << std::endl; )
644 
645  return starter;
646 }
647 
648 //------------------------------------------------------------------------
649 #ifdef SEND_ALL_OUTPUT_TO_FILES
650 static
651 void redirect_output(
652  std::ostream& myerrstream,
653  std::ostream& myinfostream
654  )
655 {
656  myerrstream .setf( std::ios::scientific | std::ios::showpoint );
657  myinfostream.setf( std::ios::scientific | std::ios::showpoint );
658  spxout.setStream( SPxOut::ERROR, myerrstream );
659  spxout.setStream( SPxOut::WARNING, myerrstream );
660  spxout.setStream( SPxOut::INFO1, myinfostream );
661  spxout.setStream( SPxOut::INFO2, myinfostream );
662  spxout.setStream( SPxOut::INFO3, myinfostream );
663  spxout.setStream( SPxOut::DEBUG, myinfostream );
664 }
665 #endif
666 //------------------------------------------------------------------------
667 static
668 void read_input_file(
669  MySoPlex& work,
670  const char* filename,
671  NameSet& rownames,
672  NameSet& colnames)
673 {
674  if ( checkMode )
675  MSG_INFO1( spxout << "IEXAMP22 loading LP file " << filename << std::endl; )
676  else
677  MSG_INFO1( spxout << "\nLoading LP file " << filename << std::endl; )
678 
679  Timer timer;
680  timer.start();
681 
682  if ( ! work.readFile(filename, &rownames, &colnames, 0) )
683  {
684  if ( checkMode )
685  MSG_INFO1( spxout << "EEXAMP23 error while reading file \"" << filename << "\"" << std::endl; )
686  else
687  MSG_INFO1( spxout << "error while reading file \"" << filename << "\"" << std::endl; )
688  exit(1);
689  }
690  assert(work.isConsistent());
691 
692  timer.stop();
693 
694  if ( checkMode )
695  {
696  MSG_INFO1( spxout << "IEXAMP24 LP has "
697  << work.nRows() << " rows "
698  << work.nCols() << " columns "
699  << work.nNzos() << " nonzeros"
700  << std::endl; )
701 
702  MSG_INFO1( spxout << "IEXAMP41 LP reading time: " << timer.userTime() << std::endl; )
703  }
704  else
705  {
706  MSG_INFO1( spxout << "LP has "
707  << work.nRows() << " rows "
708  << work.nCols() << " columns "
709  << work.nNzos() << " nonzeros"
710  << std::endl; )
711 
712  MSG_INFO1(
713  std::streamsize prec = spxout.precision();
714  spxout << "LP reading time: " << std::fixed << std::setprecision(2) << timer.userTime();
715  spxout << std::scientific << std::setprecision(int(prec)) << std::endl; )
716  }
717 }
718 
719 //------------------------------------------------------------------------
720 static
721 void read_basis_file(
722  MySoPlex& work ,
723  const char* filename,
724  const NameSet* rownames,
725  const NameSet* colnames)
726 {
727  MSG_INFO1( spxout << "Reading basis from file (disables simplifier)" << std::endl; )
728  if (!work.readBasisFile(filename, rownames, colnames))
729  {
730  if ( checkMode )
731  MSG_INFO1( spxout << "EEXAMP25 error while reading file \"" << filename << "\"" << std::endl; )
732  else
733  MSG_INFO1( spxout << "Error while reading file \"" << filename << "\"" << std::endl; )
734  exit(1);
735  }
736 }
737 
738 //------------------------------------------------------------------------
739 static
740 void solve_LP(MySoPlex& work)
741 {
742  Timer timer;
743  timer.start();
744 
745  if ( checkMode )
746  MSG_INFO1( spxout << "IEXAMP26 solving LP" << std::endl; )
747  else
748  MSG_INFO1( spxout << "\nSolving LP ..." << std::endl; )
749 
750  work.solve();
751  timer.stop();
752 
753  MSG_INFO1( spxout << "\nSoPlex statistics:\n" << work.statistics(); )
754 }
755 
756 //------------------------------------------------------------------------
757 static
758 void print_solution_and_status(
759  MySoPlex& work,
760  const NameSet& rownames,
761  const NameSet& colnames,
762  const int precision,
763  const bool print_quality,
764  const bool print_solution,
765  const bool print_dual,
766  const bool write_basis,
767  const char* basisname
768  )
769 {
770  // get the solution status
771  SPxSolver::Status stat = work.status();
772 
773  if ( ! checkMode )
774  MSG_INFO1( spxout << std::endl; )
775  switch (stat)
776  {
777  case SPxSolver::OPTIMAL:
778  if ( checkMode )
779  MSG_INFO1( spxout << "IEXAMP29 solution value is: " << std::setprecision( precision ) << work.objValue() << std::endl; )
780  else
781  MSG_INFO1( spxout << "Solution value is: " << std::setprecision( precision ) << work.objValue() << std::endl; )
782 
783  if ( print_quality )
784  work.displayQuality();
785 
786  if ( print_solution )
787  {
788  DVector objx(work.nCols());
789 
790  if( work.getPrimal(objx) != SPxSolver::ERROR )
791  {
792  MSG_INFO1( spxout << std::endl << "Primal solution (name, id, value):" << std::endl; )
793  for( int i = 0; i < work.nCols(); ++i )
794  {
795  if ( isNotZero( objx[i], 0.001 * work.feastol() ) )
796  MSG_INFO1( spxout << colnames[ work.cId(i) ] << "\t"
797  << i << "\t"
798  << std::setw(17)
799  << std::setprecision( precision )
800  << objx[i] << std::endl; )
801  }
802  MSG_INFO1( spxout << "All other variables are zero (within " << std::setprecision(1) << 0.001*work.feastol() << ")." << std::endl; )
803  }
804  }
805  if ( print_dual )
806  {
807  DVector objy(work.nRows());
808  bool allzero = true;
809 
810  if( work.getDual(objy) != SPxSolver::ERROR )
811  {
812  MSG_INFO1( spxout << std::endl << "Dual multipliers (name, id, value):" << std::endl; )
813  for( int i = 0; i < work.nRows(); ++i )
814  {
815  if ( isNotZero( objy[i] , 0.001 * work.opttol() ) )
816  {
817  MSG_INFO1( spxout << rownames[ work.rId(i) ] << "\t"
818  << i << "\t"
819  << std::setw(17)
820  << std::setprecision( precision )
821  << objy[i] << std::endl; )
822  allzero = false;
823  }
824  }
825 
826  MSG_INFO1( spxout << "All " << (allzero ? "" : "other ") << "dual values are zero (within "
827  << std::setprecision(1) << 0.001*work.opttol() << ")." << std::endl; )
828 
829  if( !allzero )
830  {
831  if( work.spxSense() == SPxLP::MINIMIZE )
832  {
833  MSG_INFO1( spxout << "Minimizing: a positive/negative value corresponds to left-hand (>=) resp. right-hand (<=) side."
834  << std::endl; )
835  }
836  else
837  {
838  MSG_INFO1( spxout << "Maximizing: a positive/negative value corresponds to right-hand (<=) resp. left-hand (>=) side."
839  << std::endl; )
840  }
841  }
842  }
843  }
844  if ( write_basis )
845  {
846  MSG_INFO1( spxout << "Writing basis of original problem to file " << basisname << std::endl; )
847  if ( ! work.writeBasisFile( basisname, &rownames, &colnames ) )
848  {
849  if ( checkMode )
850  MSG_INFO1( spxout << "EEXAMP30 error while writing file \"" << basisname << "\"" << std::endl; )
851  else
852  MSG_INFO1( spxout << "Error while writing file \"" << basisname << "\"" << std::endl; )
853  }
854  }
855  break;
856  case SPxSolver::UNBOUNDED:
857  if ( checkMode )
858  MSG_INFO1( spxout << "IEXAMP31 LP is unbounded" << std::endl; )
859  else
860  MSG_INFO1( spxout << "LP is unbounded" << std::endl; )
861 
862  if ( print_solution )
863  {
864  DVector objx(work.nCols());
865  if( work.getPrimal(objx) != SPxSolver::ERROR )
866  {
867  MSG_INFO1( spxout << std::endl << "Primal solution (name, id, value):" << std::endl; )
868  for( int i = 0; i < work.nCols(); ++i )
869  {
870  if ( isNotZero( objx[i], 0.001 * work.feastol() ) )
871  MSG_INFO1( spxout << colnames[ work.cId(i) ] << "\t"
872  << i << "\t"
873  << std::setw(17)
874  << std::setprecision( precision )
875  << objx[i] << std::endl; )
876  }
877  MSG_INFO1( spxout << "All other variables are zero (within " << std::setprecision(1) << 0.001*work.feastol() << ")." << std::endl; )
878  }
879 
880  DVector objcoef(work.nCols());
881  DVector ray(work.nCols());
882  if( work.getPrimalray(ray) != SPxSolver::ERROR )
883  {
884  Real rayobjval = 0.0;
885 
886  work.getObj(objcoef);
887 
888  MSG_INFO1( spxout << std::endl << "Primal ray (name, id, value):" << std::endl; )
889  for( int i = 0; i < work.nCols(); ++i )
890  {
891  if ( isNotZero( ray[i], 0.001 * work.feastol() ) )
892  {
893  rayobjval += ray[i] * objcoef[i];
894 
895  MSG_INFO1( spxout << colnames[ work.cId(i) ] << "\t"
896  << i << "\t"
897  << std::setw(17)
898  << std::setprecision( precision )
899  << ray[i] << std::endl; )
900  }
901  }
902  MSG_INFO1( spxout << "All other variables have zero value (within " << std::setprecision(1) << 0.001*work.feastol() << ")." << std::endl; )
903  MSG_INFO1( spxout << "Objective change per unit along primal ray is " << rayobjval << "." << std::endl; )
904  }
905  }
906  break;
907  case SPxSolver::INFEASIBLE:
908  if ( checkMode )
909  MSG_INFO1( spxout << "IEXAMP32 LP is infeasible" << std::endl; )
910  else
911  MSG_INFO1( spxout << "LP is infeasible" << std::endl; )
912  if ( print_solution )
913  {
914  DVector farkasx(work.nRows());
915 
916  if( work.getDualfarkas(farkasx) != SPxSolver::ERROR )
917  {
918  DVector proofvec(work.nCols());
919  double lhs;
920  double rhs;
921 
922  lhs = 0.0;
923  rhs = 0.0;
924  proofvec.clear();
925  for( int i = 0; i < work.nRows(); ++i )
926  {
927  if ( isNotZero( farkasx[i], 0.001 * work.opttol() ) )
928  {
929  MSG_INFO1( spxout << rownames[ work.rId(i) ] << "\t"
930  << i << "\t"
931  << std::setw(16)
932  << std::setprecision( precision )
933  << farkasx[i] << "\t"; )
934  LPRow row;
935  work.getRow(i, row);
936  if( row.lhs() > -soplex::infinity )
937  {
938  MSG_INFO1( spxout << row.lhs() << " <= "; );
939  }
940  for( int j = 0; j < row.rowVector().size(); ++j )
941  {
942  if( row.rowVector().value(j) > 0 )
943  {
944  MSG_INFO1( spxout << "+"; )
945  }
946  MSG_INFO1( spxout
947  << row.rowVector().value(j) << " "
948  << colnames[ work.cId(row.rowVector().index(j)) ]
949  << " "; );
950  }
951  if( row.rhs() < soplex::infinity )
952  {
953  MSG_INFO1( spxout << "<= " << row.rhs(); );
954  }
955  MSG_INFO1( spxout << std::endl; )
956  if( farkasx[i] > 0.0 )
957  {
958  lhs += farkasx[i] * row.lhs();
959  rhs += farkasx[i] * row.rhs();
960  }
961  else
962  {
963  lhs += farkasx[i] * row.rhs();
964  rhs += farkasx[i] * row.lhs();
965  }
966  SVector vec(row.rowVector());
967  vec *= farkasx[i];
968  proofvec += vec;
969  }
970  }
971 
972  MSG_INFO1( spxout << "All other row multipliers are zero (within " << std::setprecision(1) << 0.001*work.opttol() << ")." << std::endl; )
973  MSG_INFO1( spxout << "Farkas infeasibility proof: \t"; )
974  MSG_INFO1( spxout << lhs << " <= "; )
975 
976  bool nonzerofound = false;
977  for( int i = 0; i < work.nCols(); ++i )
978  {
979  if ( isNotZero( proofvec[i], 0.001 * work.opttol() ) )
980  {
981  if( proofvec[i] > 0 )
982  {
983  MSG_INFO1( spxout << "+"; )
984  }
985  MSG_INFO1( spxout << proofvec[i] << " " << colnames[ work.cId(i) ] << " "; )
986  nonzerofound = true;
987  }
988  }
989  if( !nonzerofound )
990  {
991  MSG_INFO1( spxout << "0 "; );
992  }
993  MSG_INFO1( spxout << "<= " << rhs << std::endl; );
994  }
995  }
996  if ( print_quality )
997  work.displayInfeasibility();
998  if ( write_basis ) // write basis even if we are infeasible
999  if ( ! work.writeBasisFile( basisname, &rownames, &colnames ) )
1000  {
1001  if ( checkMode )
1002  MSG_INFO1( spxout << "EEXAMP30 error while writing file \"" << basisname << "\"" << std::endl; )
1003  else
1004  MSG_INFO1( spxout << "Error while writing file \"" << basisname << "\"" << std::endl; )
1005  }
1006  break;
1007  case SPxSolver::ABORT_CYCLING:
1008  if ( checkMode )
1009  MSG_INFO1( spxout << "EEXAMP40 aborted due to cycling" << std::endl; )
1010  else
1011  MSG_INFO1( spxout << "Aborted due to cycling" << std::endl; )
1012  break;
1013  case SPxSolver::ABORT_TIME:
1014  if ( checkMode )
1015  MSG_INFO1( spxout << "IEXAMP33 aborted due to time limit" << std::endl; )
1016  else
1017  MSG_INFO1( spxout << "Aborted due to time limit" << std::endl; )
1018  break;
1019  case SPxSolver::ABORT_ITER:
1020  if ( checkMode )
1021  MSG_INFO1( spxout << "IEXAMP34 aborted due to iteration limit" << std::endl; )
1022  else
1023  MSG_INFO1( spxout << "Aborted due to iteration limit" << std::endl; )
1024  break;
1025  case SPxSolver::ABORT_VALUE:
1026  if ( checkMode )
1027  MSG_INFO1( spxout << "IEXAMP35 aborted due to objective value limit" << std::endl; )
1028  else
1029  MSG_INFO1( spxout << "Aborted due to objective value limit" << std::endl; )
1030  break;
1031  case SPxSolver::SINGULAR:
1032  if ( checkMode )
1033  MSG_INFO1( spxout << "EEXAMP39 basis is singular" << std::endl; )
1034  else
1035  MSG_INFO1( spxout << "Basis is singular" << std::endl; )
1036  break;
1037  default:
1038  if ( checkMode )
1039  MSG_INFO1( spxout << "EEXAMP36 An error occurred during " << "the solution process" << std::endl; )
1040  else
1041  MSG_INFO1( spxout << "An error occurred during " << "the solution process" << std::endl; )
1042  break;
1043  }
1044  MSG_INFO1( spxout << std::endl; )
1045 }
1046 
1047 //------------------------------------------------------------------------
1048 static
1049 void clean_up(
1050  SPxScaler*& prescaler,
1051  SPxScaler*& postscaler,
1052  SPxSimplifier*& simplifier,
1053  SPxStarter*& starter,
1054  SPxPricer*& pricer,
1055  SPxRatioTester*& ratiotester,
1056  char*& basisname
1057  )
1058 {
1059  if ( prescaler != 0 )
1060  {
1061  delete prescaler;
1062  prescaler = 0;
1063  }
1064  if ( postscaler != 0 )
1065  {
1066  delete postscaler;
1067  postscaler = 0;
1068  }
1069  if ( simplifier != 0 )
1070  {
1071  delete simplifier;
1072  simplifier = 0;
1073  }
1074  if ( starter != 0 )
1075  {
1076  delete starter;
1077  starter = 0;
1078  }
1079 
1080  assert( pricer != 0 );
1081  delete pricer;
1082  pricer = 0;
1083 
1084  assert( ratiotester != 0 );
1085  delete ratiotester;
1086  ratiotester = 0;
1087 
1088  if ( basisname != 0 )
1089  delete [] basisname;
1090  basisname = 0;
1091 }
1092 
1093 //------------------------------------------------------------------------
1094 // main program
1095 //------------------------------------------------------------------------
1096 
1097 int main(int argc, char* argv[])
1098 {
1099  const char* filename;
1100  char* basisname = 0;
1101  SPxSolver::Type type = SPxSolver::LEAVE;
1102  SPxSolver::Representation representation = SPxSolver::COLUMN;
1103  SLUFactor::UpdateType update = SLUFactor::FOREST_TOMLIN;
1104  SPxSimplifier* simplifier = 0;
1105  SPxStarter* starter = 0;
1106  SPxPricer* pricer = 0;
1107  SPxRatioTester* ratiotester = 0;
1108  SPxScaler* prescaler = 0;
1109  SPxScaler* postscaler = 0;
1110 
1111  try {
1112  NameSet rownames;
1113  NameSet colnames;
1114  int starting = 0;
1115  int pricing = 4;
1116  int ratiotest = 2;
1117  int scaling = 2;
1118  int simplifying = 1;
1119  int iterlimit = -1;
1120  Real timelimit = -1.0;
1121  Real delta = DEFAULT_BND_VIOL;
1122  Real feastol = DEFAULT_BND_VIOL;
1123  Real opttol = DEFAULT_BND_VIOL;
1124  Real irthreshold = DEFAULT_BND_VIOL * 1e-6;
1125  Real epsilon = DEFAULT_EPS_ZERO;
1126  Real epsilon_factor = DEFAULT_EPS_FACTOR;
1127  Real epsilon_update = DEFAULT_EPS_UPDATE;
1128  int verbose = SPxOut::INFO1;
1129  bool print_solution = false;
1130  bool print_dual = false;
1131  bool print_quality = false;
1132  bool read_basis = false;
1133  bool write_basis = false;
1134  int precision;
1135  int optidx;
1136 
1137  for(optidx = 1; optidx < argc; optidx++)
1138  {
1139  if (*argv[optidx] != '-')
1140  break;
1141 
1142  switch(argv[optidx][1])
1143  {
1144  case 'b' :
1145  check_parameter(argv[optidx][2], argv); // use -b{r,w}, not -b
1146  if (argv[optidx][2] == 'r')
1147  read_basis = true;
1148  if (argv[optidx][2] == 'w')
1149  write_basis = true;
1150  break;
1151  case 'c' :
1152  check_parameter(argv[optidx][2], argv); // use -c[0-3], not -c
1153  starting = atoi(&argv[optidx][2]);
1154  break;
1155  case 'd' :
1156  check_parameter(argv[optidx][2], argv); // use -dx, not -d
1157  delta = atof(&argv[optidx][2]);
1158  break;
1159  case 'f' :
1160  check_parameter(argv[optidx][2], argv); // use -fx, not -f
1161  feastol = atof(&argv[optidx][2]);
1162  break;
1163  case 'o' :
1164  check_parameter(argv[optidx][2], argv); // use -ox, not -o
1165  opttol = atof(&argv[optidx][2]);
1166  break;
1167  case 'R' :
1168  check_parameter(argv[optidx][2], argv); // use -Rx, not -R
1169  irthreshold = atof(&argv[optidx][2]);
1170  break;
1171  case 'e':
1172  type = SPxSolver::ENTER;
1173  break;
1174  case 'g' :
1175  check_parameter(argv[optidx][2], argv); // use -g[0-5], not -g
1176  scaling = atoi(&argv[optidx][2]);
1177  break;
1178  case 'i' :
1179  update = SLUFactor::ETA;
1180  break;
1181  case 'l' :
1182  if (argv[optidx][2] == '\0' ) // use -lx, not -l
1183  print_usage_and_exit( argv );
1184  timelimit = atof(&argv[optidx][2]);
1185  break;
1186  case 'L' :
1187  if (argv[optidx][2] == '\0' ) // use -Lx, not -L
1188  print_usage_and_exit( argv );
1189  iterlimit = atoi(&argv[optidx][2]);
1190  break;
1191  case 'p' :
1192  check_parameter(argv[optidx][2], argv); // use -p[0-5], not -p
1193  pricing = atoi(&argv[optidx][2]);
1194  break;
1195  case 'q' :
1196  print_quality = true;
1197  break;
1198  case 'r' :
1199  representation = SPxSolver::ROW;
1200  break;
1201  case 's' :
1202  check_parameter(argv[optidx][2], argv); // use -s[0-4], not -s
1203  simplifying = atoi(&argv[optidx][2]);
1204  break;
1205  case 't' :
1206  check_parameter(argv[optidx][2], argv); // use -r[0-2], not -r
1207  ratiotest = atoi(&argv[optidx][2]);
1208  break;
1209  case 'v' :
1210  check_parameter(argv[optidx][2], argv); // use -v[0-5], not -v
1211  if (argv[optidx][2] >= '0' && argv[optidx][2] <= '9')
1212  verbose = argv[optidx][2] - '0';
1213  break;
1214  case 'V' :
1215  print_version_info();
1216  exit(0);
1217  case 'x' :
1218  print_solution = true;
1219  break;
1220  case 'y' :
1221  print_dual = true;
1222  break;
1223  case 'z' :
1224  check_parameter(argv[optidx][2], argv); // must not be empty
1225  check_parameter(argv[optidx][3], argv); // must not be empty
1226  switch(argv[optidx][2])
1227  {
1228  case 'z' :
1229  epsilon = atof(&argv[optidx][3]);
1230  break;
1231  case 'f' :
1232  epsilon_factor = atof(&argv[optidx][3]);
1233  break;
1234  case 'u' :
1235  epsilon_update = atof(&argv[optidx][3]);
1236  break;
1237  default :
1238  print_usage_and_exit( argv );
1239  }
1240  break;
1241  case 'C' :
1242  checkMode = true;
1243  break;
1244  case 'h' :
1245  case '?' :
1246  print_version_info();
1247  //lint -fallthrough
1248  default :
1249  print_usage_and_exit( argv );
1250  }
1251  }
1252 
1253  // print version
1254  print_version_info();
1255 
1256  // enough arguments?
1257  if ((argc - optidx) < 1 + (read_basis ? 1 : 0) + (write_basis ? 1 : 0))
1258  print_usage_and_exit( argv );
1259  filename = argv[optidx];
1260 
1261  ++optidx;
1262 
1263  // switch off simplifier when using a starting basis
1264  if ( read_basis )
1265  simplifying = 0;
1266 
1267  if ( read_basis || write_basis )
1268  basisname = strcpy( new char[strlen(argv[optidx]) + 1], argv[optidx] );
1269 
1270  // set some algorithm parameters
1271  Param::setEpsilon ( epsilon );
1272  Param::setEpsilonFactorization( epsilon_factor );
1273  Param::setEpsilonUpdate ( epsilon_update );
1274  Param::setVerbose ( verbose );
1275 
1276  // Set the output precision.
1277  precision = int(-log10(std::min(feastol, opttol))) + 1;
1278 
1279  std::cout.setf( std::ios::scientific | std::ios::showpoint );
1280  std::cerr.setf( std::ios::scientific | std::ios::showpoint );
1281 
1282 #ifdef SEND_ALL_OUTPUT_TO_FILES
1283  // Example of redirecting output to different files.
1284  // Default is cerr for errors and warnings, cout for everything else.
1285  std::ofstream myerrstream ( "errwarn.txt" );
1286  std::ofstream myinfostream( "infos.txt" );
1287  redirect_output(myerrstream, myinfostream);
1288 #endif
1289 
1290  // create an instance of MySoPlex
1291  MySoPlex work( type, representation );
1292  work.setUtype ( update );
1293  work.setFeastol ( std::min(feastol, delta) );
1294  work.setOpttol ( std::min(opttol, delta) );
1295  work.setIrthreshold ( irthreshold );
1296  work.setTerminationTime ( timelimit );
1297  work.setTerminationIter ( iterlimit );
1298  print_algorithm_parameters( work, representation, update );
1299  assert( work.isConsistent() );
1300 
1301  // set pricer, starter, simplifier, and ratio tester
1302  work.setPricer ( pricer = get_pricer (pricing) );
1303  work.setStarter ( starter = get_starter (starting) );
1304  work.setSimplifier( simplifier = get_simplifier (simplifying) );
1305  work.setTester ( ratiotester = get_ratio_tester(ratiotest) );
1306  assert(work.isConsistent());
1307 
1308  // set pre- and postscaler
1309  get_scalers(prescaler, postscaler, scaling);
1310  work.setPreScaler (prescaler);
1311  work.setPostScaler(postscaler);
1312  assert(work.isConsistent());
1313 
1314  // read the LP from an input file (.lp or .mps)
1315  read_input_file(work, filename, rownames, colnames);
1316 
1317  // read a basis file if specified
1318  if (read_basis)
1319  read_basis_file(work, basisname, &rownames, &colnames);
1320 
1321  // solve the LP
1322  solve_LP(work);
1323 
1324  // print solution, status, infeasibility system,...
1325  print_solution_and_status(work, rownames, colnames, precision, print_quality,
1326  print_solution, print_dual, write_basis, basisname);
1327 
1328  // clean up
1329  clean_up(prescaler, postscaler, simplifier, starter, pricer, ratiotester, basisname);
1330 
1331  return 0;
1332  }
1333  catch(SPxException& x) {
1334  std::cout << "exception caught : " << x.what() << std::endl;
1335  delete [] basisname;
1336  if (simplifier)
1337  delete simplifier;
1338  delete starter;
1339  delete pricer;
1340  delete ratiotester;
1341  delete prescaler;
1342  delete postscaler;
1343  }
1344 }
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen