source: trunk/Clp/src/ClpMain.cpp @ 1665

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 179.0 KB
Line 
1/* $Id: ClpMain.cpp 1665 2011-01-04 17:55:54Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include <cassert>
9#include <cstdio>
10#include <cstdlib>
11#include <cmath>
12#include <cfloat>
13#include <string>
14#include <iostream>
15
16int boundary_sort = 1000;
17int boundary_sort2 = 1000;
18int boundary_sort3 = 10000;
19
20#include "CoinPragma.hpp"
21#include "CoinHelperFunctions.hpp"
22#include "CoinSort.hpp"
23// History since 1.0 at end
24#include "ClpConfig.h"
25#include "CoinMpsIO.hpp"
26#include "CoinFileIO.hpp"
27
28#include "ClpFactorization.hpp"
29#include "CoinTime.hpp"
30#include "ClpSimplex.hpp"
31#include "ClpSimplexOther.hpp"
32#include "ClpSolve.hpp"
33#include "ClpPackedMatrix.hpp"
34#include "ClpPlusMinusOneMatrix.hpp"
35#include "ClpNetworkMatrix.hpp"
36#include "ClpDualRowSteepest.hpp"
37#include "ClpDualRowDantzig.hpp"
38#include "ClpLinearObjective.hpp"
39#include "ClpPrimalColumnSteepest.hpp"
40#include "ClpPrimalColumnDantzig.hpp"
41#include "ClpPresolve.hpp"
42#include "CbcOrClpParam.hpp"
43#include "CoinSignal.hpp"
44#ifdef DMALLOC
45#include "dmalloc.h"
46#endif
47#ifdef WSSMP_BARRIER
48#define FOREIGN_BARRIER
49#endif
50#ifdef UFL_BARRIER
51#define FOREIGN_BARRIER
52#endif
53#ifdef TAUCS_BARRIER
54#define FOREIGN_BARRIER
55#endif
56#ifdef MUMPS_BARRIER
57#define FOREIGN_BARRIER
58#endif
59
60static double totalTime = 0.0;
61static bool maskMatches(const int * starts, char ** masks,
62                        std::string & check);
63static ClpSimplex * currentModel = NULL;
64
65extern "C" {
66     static void
67#if defined(_MSC_VER)
68     __cdecl
69#endif // _MSC_VER
70     signal_handler(int /*whichSignal*/)
71     {
72          if (currentModel != NULL)
73               currentModel->setMaximumIterations(0); // stop at next iterations
74          return;
75     }
76}
77
78//#############################################################################
79
80#ifdef NDEBUG
81#undef NDEBUG
82#endif
83
84int mainTest (int argc, const char *argv[], int algorithm,
85              ClpSimplex empty, ClpSolve solveOptions, int switchOff, bool doVector);
86static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
87static void generateCode(const char * fileName, int type);
88// Returns next valid field
89int CbcOrClpRead_mode = 1;
90FILE * CbcOrClpReadCommand = stdin;
91extern int CbcOrClpEnvironmentIndex;
92#ifdef CLP_USER_DRIVEN1
93/* Returns true if variable sequenceOut can leave basis when
94   model->sequenceIn() enters.
95   This function may be entered several times for each sequenceOut. 
96   The first time realAlpha will be positive if going to lower bound
97   and negative if going to upper bound (scaled bounds in lower,upper) - then will be zero.
98   currentValue is distance to bound.
99   currentTheta is current theta.
100   alpha is fabs(pivot element).
101   Variable will change theta if currentValue - currentTheta*alpha < 0.0
102*/
103bool userChoiceValid1(const ClpSimplex * model,
104                      int sequenceOut,
105                      double currentValue,
106                      double currentTheta,
107                      double alpha,
108                      double realAlpha)
109{
110  return true;
111}
112/* This returns true if chosen in/out pair valid.
113   The main thing to check would be variable flipping bounds may be
114   OK.  This would be signaled by reasonable theta_ and valueOut_.
115   If you return false sequenceIn_ will be flagged as ineligible.
116*/
117bool userChoiceValid2(const ClpSimplex * model)
118{
119  return true;
120}
121/* If a good pivot then you may wish to unflag some variables.
122 */
123void userChoiceWasGood(ClpSimplex * model)
124{
125}
126#endif
127
128int
129#if defined(_MSC_VER)
130__cdecl
131#endif // _MSC_VER
132main (int argc, const char *argv[])
133{
134     // next {} is just to make sure all memory should be freed - for debug
135     {
136          double time1 = CoinCpuTime(), time2;
137          // Set up all non-standard stuff
138          //int numberModels=1;
139          ClpSimplex * models = new ClpSimplex[1];
140
141          // default action on import
142          int allowImportErrors = 0;
143          int keepImportNames = 1;
144          int doIdiot = -1;
145          int outputFormat = 2;
146          int slpValue = -1;
147          int cppValue = -1;
148          int printOptions = 0;
149          int printMode = 0;
150          int presolveOptions = 0;
151          int doCrash = 0;
152          int doVector = 0;
153          int doSprint = -1;
154          // set reasonable defaults
155          int preSolve = 5;
156          bool preSolveFile = false;
157          models->setPerturbation(50);
158          models->messageHandler()->setPrefix(false);
159          const char dirsep =  CoinFindDirSeparator();
160          std::string directory;
161          std::string dirSample;
162          std::string dirNetlib;
163          std::string dirMiplib;
164          if (dirsep == '/') {
165               directory = "./";
166               dirSample = "../../Data/Sample/";
167               dirNetlib = "../../Data/Netlib/";
168               dirMiplib = "../../Data/miplib3/";
169          } else {
170               directory = ".\\";
171#              ifdef COIN_MSVS
172               // Visual Studio builds are deeper
173               dirSample = "..\\..\\..\\..\\Data\\Sample\\";
174               dirNetlib = "..\\..\\..\\..\\Data\\Netlib\\";
175               dirMiplib = "..\\..\\..\\..\\Data\\miplib3\\";
176#              else
177               dirSample = "..\\..\\Data\\Sample\\";
178               dirNetlib = "..\\..\\Data\\Netlib\\";
179               dirMiplib = "..\\..\\Data\\miplib3\\";
180#              endif
181          }
182          std::string defaultDirectory = directory;
183          std::string importFile = "";
184          std::string exportFile = "default.mps";
185          std::string importBasisFile = "";
186          int basisHasValues = 0;
187          int substitution = 3;
188          int dualize = 3;  // dualize if looks promising
189          std::string exportBasisFile = "default.bas";
190          std::string saveFile = "default.prob";
191          std::string restoreFile = "default.prob";
192          std::string solutionFile = "stdout";
193          std::string solutionSaveFile = "solution.file";
194          std::string printMask = "";
195          CbcOrClpParam parameters[CBCMAXPARAMETERS];
196          int numberParameters ;
197          establishParams(numberParameters, parameters) ;
198          parameters[whichParam(CLP_PARAM_ACTION_BASISIN, numberParameters, parameters)].setStringValue(importBasisFile);
199          parameters[whichParam(CLP_PARAM_ACTION_BASISOUT, numberParameters, parameters)].setStringValue(exportBasisFile);
200          parameters[whichParam(CLP_PARAM_ACTION_PRINTMASK, numberParameters, parameters)].setStringValue(printMask);
201          parameters[whichParam(CLP_PARAM_ACTION_DIRECTORY, numberParameters, parameters)].setStringValue(directory);
202          parameters[whichParam(CLP_PARAM_ACTION_DIRSAMPLE, numberParameters, parameters)].setStringValue(dirSample);
203          parameters[whichParam(CLP_PARAM_ACTION_DIRNETLIB, numberParameters, parameters)].setStringValue(dirNetlib);
204          parameters[whichParam(CBC_PARAM_ACTION_DIRMIPLIB, numberParameters, parameters)].setStringValue(dirMiplib);
205          parameters[whichParam(CLP_PARAM_DBL_DUALBOUND, numberParameters, parameters)].setDoubleValue(models->dualBound());
206          parameters[whichParam(CLP_PARAM_DBL_DUALTOLERANCE, numberParameters, parameters)].setDoubleValue(models->dualTolerance());
207          parameters[whichParam(CLP_PARAM_ACTION_EXPORT, numberParameters, parameters)].setStringValue(exportFile);
208          parameters[whichParam(CLP_PARAM_INT_IDIOT, numberParameters, parameters)].setIntValue(doIdiot);
209          parameters[whichParam(CLP_PARAM_ACTION_IMPORT, numberParameters, parameters)].setStringValue(importFile);
210          parameters[whichParam(CLP_PARAM_INT_SOLVERLOGLEVEL, numberParameters, parameters)].setIntValue(models->logLevel());
211          parameters[whichParam(CLP_PARAM_INT_MAXFACTOR, numberParameters, parameters)].setIntValue(models->factorizationFrequency());
212          parameters[whichParam(CLP_PARAM_INT_MAXITERATION, numberParameters, parameters)].setIntValue(models->maximumIterations());
213          parameters[whichParam(CLP_PARAM_INT_OUTPUTFORMAT, numberParameters, parameters)].setIntValue(outputFormat);
214          parameters[whichParam(CLP_PARAM_INT_PRESOLVEPASS, numberParameters, parameters)].setIntValue(preSolve);
215          parameters[whichParam(CLP_PARAM_INT_PERTVALUE, numberParameters, parameters)].setIntValue(models->perturbation());
216          parameters[whichParam(CLP_PARAM_DBL_PRIMALTOLERANCE, numberParameters, parameters)].setDoubleValue(models->primalTolerance());
217          parameters[whichParam(CLP_PARAM_DBL_PRIMALWEIGHT, numberParameters, parameters)].setDoubleValue(models->infeasibilityCost());
218          parameters[whichParam(CLP_PARAM_ACTION_RESTORE, numberParameters, parameters)].setStringValue(restoreFile);
219          parameters[whichParam(CLP_PARAM_ACTION_SAVE, numberParameters, parameters)].setStringValue(saveFile);
220          parameters[whichParam(CLP_PARAM_DBL_TIMELIMIT, numberParameters, parameters)].setDoubleValue(models->maximumSeconds());
221          parameters[whichParam(CLP_PARAM_ACTION_SOLUTION, numberParameters, parameters)].setStringValue(solutionFile);
222          parameters[whichParam(CLP_PARAM_ACTION_SAVESOL, numberParameters, parameters)].setStringValue(solutionSaveFile);
223          parameters[whichParam(CLP_PARAM_INT_SPRINT, numberParameters, parameters)].setIntValue(doSprint);
224          parameters[whichParam(CLP_PARAM_INT_SUBSTITUTION, numberParameters, parameters)].setIntValue(substitution);
225          parameters[whichParam(CLP_PARAM_INT_DUALIZE, numberParameters, parameters)].setIntValue(dualize);
226          parameters[whichParam(CLP_PARAM_DBL_PRESOLVETOLERANCE, numberParameters, parameters)].setDoubleValue(1.0e-8);
227          int verbose = 0;
228
229          // total number of commands read
230          int numberGoodCommands = 0;
231          bool * goodModels = new bool[1];
232
233          // Hidden stuff for barrier
234          int choleskyType = 0;
235          int gamma = 0;
236          parameters[whichParam(CLP_PARAM_STR_BARRIERSCALE, numberParameters, parameters)].setCurrentOption(2);
237          int scaleBarrier = 2;
238          int doKKT = 0;
239          int crossover = 2; // do crossover unless quadratic
240
241          int iModel = 0;
242          goodModels[0] = false;
243          //models[0].scaling(1);
244          //models[0].setDualBound(1.0e6);
245          //models[0].setDualTolerance(1.0e-7);
246          //ClpDualRowSteepest steep;
247          //models[0].setDualRowPivotAlgorithm(steep);
248          //models[0].setPrimalTolerance(1.0e-7);
249          //ClpPrimalColumnSteepest steepP;
250          //models[0].setPrimalColumnPivotAlgorithm(steepP);
251          std::string field;
252          std::cout << "Coin LP version " << CLP_VERSION
253                    << ", build " << __DATE__ << std::endl;
254          // Print command line
255          if (argc > 1) {
256               printf("command line - ");
257               for (int i = 0; i < argc; i++)
258                    printf("%s ", argv[i]);
259               printf("\n");
260          }
261
262          while (1) {
263               // next command
264               field = CoinReadGetCommand(argc, argv);
265
266               // exit if null or similar
267               if (!field.length()) {
268                    if (numberGoodCommands == 1 && goodModels[0]) {
269                         // we just had file name - do dual or primal
270                         field = "either";
271                    } else if (!numberGoodCommands) {
272                         // let's give the sucker a hint
273                         std::cout
274                                   << "Clp takes input from arguments ( - switches to stdin)"
275                                   << std::endl
276                                   << "Enter ? for list of commands or help" << std::endl;
277                         field = "-";
278                    } else {
279                         break;
280                    }
281               }
282
283               // see if ? at end
284               int numberQuery = 0;
285               if (field != "?" && field != "???") {
286                    int length = field.length();
287                    int i;
288                    for (i = length - 1; i > 0; i--) {
289                         if (field[i] == '?')
290                              numberQuery++;
291                         else
292                              break;
293                    }
294                    field = field.substr(0, length - numberQuery);
295               }
296               // find out if valid command
297               int iParam;
298               int numberMatches = 0;
299               int firstMatch = -1;
300               for ( iParam = 0; iParam < numberParameters; iParam++ ) {
301                    int match = parameters[iParam].matches(field);
302                    if (match == 1) {
303                         numberMatches = 1;
304                         firstMatch = iParam;
305                         break;
306                    } else {
307                         if (match && firstMatch < 0)
308                              firstMatch = iParam;
309                         numberMatches += match >> 1;
310                    }
311               }
312               if (iParam < numberParameters && !numberQuery) {
313                    // found
314                    CbcOrClpParam found = parameters[iParam];
315                    CbcOrClpParameterType type = found.type();
316                    int valid;
317                    numberGoodCommands++;
318                    if (type == CBC_PARAM_GENERALQUERY) {
319                         std::cout << "In argument list keywords have leading - "
320                                   ", -stdin or just - switches to stdin" << std::endl;
321                         std::cout << "One command per line (and no -)" << std::endl;
322                         std::cout << "abcd? gives list of possibilities, if only one + explanation" << std::endl;
323                         std::cout << "abcd?? adds explanation, if only one fuller help" << std::endl;
324                         std::cout << "abcd without value (where expected) gives current value" << std::endl;
325                         std::cout << "abcd value sets value" << std::endl;
326                         std::cout << "Commands are:" << std::endl;
327                         int maxAcross = 10;
328                         bool evenHidden = false;
329                         int printLevel =
330                              parameters[whichParam(CLP_PARAM_STR_ALLCOMMANDS,
331                                                    numberParameters, parameters)].currentOptionAsInteger();
332                         int convertP[] = {2, 1, 0};
333                         printLevel = convertP[printLevel];
334                         if ((verbose & 8) != 0) {
335                              // even hidden
336                              evenHidden = true;
337                              verbose &= ~8;
338                         }
339                         if (verbose)
340                              maxAcross = 1;
341                         int limits[] = {1, 101, 201, 301, 401};
342                         std::vector<std::string> types;
343                         types.push_back("Double parameters:");
344                         types.push_back("Int parameters:");
345                         types.push_back("Keyword parameters:");
346                         types.push_back("Actions or string parameters:");
347                         int iType;
348                         for (iType = 0; iType < 4; iType++) {
349                              int across = 0;
350                              int lengthLine = 0;
351                              if ((verbose % 4) != 0)
352                                   std::cout << std::endl;
353                              std::cout << types[iType] << std::endl;
354                              if ((verbose & 2) != 0)
355                                   std::cout << std::endl;
356                              for ( iParam = 0; iParam < numberParameters; iParam++ ) {
357                                   int type = parameters[iParam].type();
358                                   //printf("%d type %d limits %d %d display %d\n",iParam,
359                                   //   type,limits[iType],limits[iType+1],parameters[iParam].displayThis());
360                                   if ((parameters[iParam].displayThis() >= printLevel || evenHidden) &&
361                                             type >= limits[iType]
362                                             && type < limits[iType+1]) {
363                                        if (!across) {
364                                             if ((verbose & 2) != 0)
365                                                  std::cout << "Command ";
366                                        }
367                                        int length = parameters[iParam].lengthMatchName() + 1;
368                                        if (lengthLine + length > 80) {
369                                             std::cout << std::endl;
370                                             across = 0;
371                                             lengthLine = 0;
372                                        }
373                                        std::cout << " " << parameters[iParam].matchName();
374                                        lengthLine += length ;
375                                        across++;
376                                        if (across == maxAcross) {
377                                             across = 0;
378                                             if (verbose) {
379                                                  // put out description as well
380                                                  if ((verbose & 1) != 0)
381                                                       std::cout << parameters[iParam].shortHelp();
382                                                  std::cout << std::endl;
383                                                  if ((verbose & 2) != 0) {
384                                                       std::cout << "---- description" << std::endl;
385                                                       parameters[iParam].printLongHelp();
386                                                       std::cout << "----" << std::endl << std::endl;
387                                                  }
388                                             } else {
389                                                  std::cout << std::endl;
390                                             }
391                                        }
392                                   }
393                              }
394                              if (across)
395                                   std::cout << std::endl;
396                         }
397                    } else if (type == CBC_PARAM_FULLGENERALQUERY) {
398                         std::cout << "Full list of commands is:" << std::endl;
399                         int maxAcross = 5;
400                         int limits[] = {1, 101, 201, 301, 401};
401                         std::vector<std::string> types;
402                         types.push_back("Double parameters:");
403                         types.push_back("Int parameters:");
404                         types.push_back("Keyword parameters and others:");
405                         types.push_back("Actions:");
406                         int iType;
407                         for (iType = 0; iType < 4; iType++) {
408                              int across = 0;
409                              std::cout << types[iType] << std::endl;
410                              for ( iParam = 0; iParam < numberParameters; iParam++ ) {
411                                   int type = parameters[iParam].type();
412                                   if (type >= limits[iType]
413                                             && type < limits[iType+1]) {
414                                        if (!across)
415                                             std::cout << "  ";
416                                        std::cout << parameters[iParam].matchName() << "  ";
417                                        across++;
418                                        if (across == maxAcross) {
419                                             std::cout << std::endl;
420                                             across = 0;
421                                        }
422                                   }
423                              }
424                              if (across)
425                                   std::cout << std::endl;
426                         }
427                    } else if (type < 101) {
428                         // get next field as double
429                         double value = CoinReadGetDoubleField(argc, argv, &valid);
430                         if (!valid) {
431                              parameters[iParam].setDoubleParameter(models + iModel, value);
432                         } else if (valid == 1) {
433                              std::cout << " is illegal for double parameter " << parameters[iParam].name() << " value remains " <<
434                                        parameters[iParam].doubleValue() << std::endl;
435                         } else {
436                              std::cout << parameters[iParam].name() << " has value " <<
437                                        parameters[iParam].doubleValue() << std::endl;
438                         }
439                    } else if (type < 201) {
440                         // get next field as int
441                         int value = CoinReadGetIntField(argc, argv, &valid);
442                         if (!valid) {
443                              if (parameters[iParam].type() == CLP_PARAM_INT_PRESOLVEPASS)
444                                   preSolve = value;
445                              else if (parameters[iParam].type() == CLP_PARAM_INT_IDIOT)
446                                   doIdiot = value;
447                              else if (parameters[iParam].type() == CLP_PARAM_INT_SPRINT)
448                                   doSprint = value;
449                              else if (parameters[iParam].type() == CLP_PARAM_INT_OUTPUTFORMAT)
450                                   outputFormat = value;
451                              else if (parameters[iParam].type() == CLP_PARAM_INT_SLPVALUE)
452                                   slpValue = value;
453                              else if (parameters[iParam].type() == CLP_PARAM_INT_CPP)
454                                   cppValue = value;
455                              else if (parameters[iParam].type() == CLP_PARAM_INT_PRESOLVEOPTIONS)
456                                   presolveOptions = value;
457                              else if (parameters[iParam].type() == CLP_PARAM_INT_PRINTOPTIONS)
458                                   printOptions = value;
459                              else if (parameters[iParam].type() == CLP_PARAM_INT_SUBSTITUTION)
460                                   substitution = value;
461                              else if (parameters[iParam].type() == CLP_PARAM_INT_DUALIZE)
462                                   dualize = value;
463                              else if (parameters[iParam].type() == CLP_PARAM_INT_VERBOSE)
464                                   verbose = value;
465                              parameters[iParam].setIntParameter(models + iModel, value);
466                         } else if (valid == 1) {
467                              std::cout << " is illegal for integer parameter " << parameters[iParam].name() << " value remains " <<
468                                        parameters[iParam].intValue() << std::endl;
469                         } else {
470                              std::cout << parameters[iParam].name() << " has value " <<
471                                        parameters[iParam].intValue() << std::endl;
472                         }
473                    } else if (type < 301) {
474                         // one of several strings
475                         std::string value = CoinReadGetString(argc, argv);
476                         int action = parameters[iParam].parameterOption(value);
477                         if (action < 0) {
478                              if (value != "EOL") {
479                                   // no match
480                                   parameters[iParam].printOptions();
481                              } else {
482                                   // print current value
483                                   std::cout << parameters[iParam].name() << " has value " <<
484                                             parameters[iParam].currentOption() << std::endl;
485                              }
486                         } else {
487                              parameters[iParam].setCurrentOption(action);
488                              // for now hard wired
489                              switch (type) {
490                              case CLP_PARAM_STR_DIRECTION:
491                                   if (action == 0)
492                                        models[iModel].setOptimizationDirection(1);
493                                   else if (action == 1)
494                                        models[iModel].setOptimizationDirection(-1);
495                                   else
496                                        models[iModel].setOptimizationDirection(0);
497                                   break;
498                              case CLP_PARAM_STR_DUALPIVOT:
499                                   if (action == 0) {
500                                        ClpDualRowSteepest steep(3);
501                                        models[iModel].setDualRowPivotAlgorithm(steep);
502                                   } else if (action == 1) {
503                                        //ClpDualRowDantzig dantzig;
504                                        ClpDualRowSteepest dantzig(5);
505                                        models[iModel].setDualRowPivotAlgorithm(dantzig);
506                                   } else if (action == 2) {
507                                        // partial steep
508                                        ClpDualRowSteepest steep(2);
509                                        models[iModel].setDualRowPivotAlgorithm(steep);
510                                   } else {
511                                        ClpDualRowSteepest steep;
512                                        models[iModel].setDualRowPivotAlgorithm(steep);
513                                   }
514                                   break;
515                              case CLP_PARAM_STR_PRIMALPIVOT:
516                                   if (action == 0) {
517                                        ClpPrimalColumnSteepest steep(3);
518                                        models[iModel].setPrimalColumnPivotAlgorithm(steep);
519                                   } else if (action == 1) {
520                                        ClpPrimalColumnSteepest steep(0);
521                                        models[iModel].setPrimalColumnPivotAlgorithm(steep);
522                                   } else if (action == 2) {
523                                        ClpPrimalColumnDantzig dantzig;
524                                        models[iModel].setPrimalColumnPivotAlgorithm(dantzig);
525                                   } else if (action == 3) {
526                                        ClpPrimalColumnSteepest steep(4);
527                                        models[iModel].setPrimalColumnPivotAlgorithm(steep);
528                                   } else if (action == 4) {
529                                        ClpPrimalColumnSteepest steep(1);
530                                        models[iModel].setPrimalColumnPivotAlgorithm(steep);
531                                   } else if (action == 5) {
532                                        ClpPrimalColumnSteepest steep(2);
533                                        models[iModel].setPrimalColumnPivotAlgorithm(steep);
534                                   } else if (action == 6) {
535                                        ClpPrimalColumnSteepest steep(10);
536                                        models[iModel].setPrimalColumnPivotAlgorithm(steep);
537                                   }
538                                   break;
539                              case CLP_PARAM_STR_SCALING:
540                                   models[iModel].scaling(action);
541                                   break;
542                              case CLP_PARAM_STR_AUTOSCALE:
543                                   models[iModel].setAutomaticScaling(action != 0);
544                                   break;
545                              case CLP_PARAM_STR_SPARSEFACTOR:
546                                   models[iModel].setSparseFactorization((1 - action) != 0);
547                                   break;
548                              case CLP_PARAM_STR_BIASLU:
549                                   models[iModel].factorization()->setBiasLU(action);
550                                   break;
551                              case CLP_PARAM_STR_PERTURBATION:
552                                   if (action == 0)
553                                        models[iModel].setPerturbation(50);
554                                   else
555                                        models[iModel].setPerturbation(100);
556                                   break;
557                              case CLP_PARAM_STR_ERRORSALLOWED:
558                                   allowImportErrors = action;
559                                   break;
560                              case CLP_PARAM_STR_INTPRINT:
561                                   printMode = action;
562                                   break;
563                              case CLP_PARAM_STR_KEEPNAMES:
564                                   keepImportNames = 1 - action;
565                                   break;
566                              case CLP_PARAM_STR_PRESOLVE:
567                                   if (action == 0)
568                                        preSolve = 5;
569                                   else if (action == 1)
570                                        preSolve = 0;
571                                   else if (action == 2)
572                                        preSolve = 10;
573                                   else
574                                        preSolveFile = true;
575                                   break;
576                              case CLP_PARAM_STR_PFI:
577                                   models[iModel].factorization()->setForrestTomlin(action == 0);
578                                   break;
579                              case CLP_PARAM_STR_FACTORIZATION:
580                                   models[iModel].factorization()->forceOtherFactorization(action);
581                                   break;
582                              case CLP_PARAM_STR_CRASH:
583                                   doCrash = action;
584                                   break;
585                              case CLP_PARAM_STR_VECTOR:
586                                   doVector = action;
587                                   break;
588                              case CLP_PARAM_STR_MESSAGES:
589                                   models[iModel].messageHandler()->setPrefix(action != 0);
590                                   break;
591                              case CLP_PARAM_STR_CHOLESKY:
592                                   choleskyType = action;
593                                   break;
594                              case CLP_PARAM_STR_GAMMA:
595                                   gamma = action;
596                                   break;
597                              case CLP_PARAM_STR_BARRIERSCALE:
598                                   scaleBarrier = action;
599                                   break;
600                              case CLP_PARAM_STR_KKT:
601                                   doKKT = action;
602                                   break;
603                              case CLP_PARAM_STR_CROSSOVER:
604                                   crossover = action;
605                                   break;
606                              default:
607                                   //abort();
608                                   break;
609                              }
610                         }
611                    } else {
612                         // action
613                         if (type == CLP_PARAM_ACTION_EXIT)
614                              break; // stop all
615                         switch (type) {
616                         case CLP_PARAM_ACTION_DUALSIMPLEX:
617                         case CLP_PARAM_ACTION_PRIMALSIMPLEX:
618                         case CLP_PARAM_ACTION_EITHERSIMPLEX:
619                         case CLP_PARAM_ACTION_BARRIER:
620                              // synonym for dual
621                         case CBC_PARAM_ACTION_BAB:
622                              if (goodModels[iModel]) {
623                                   double objScale =
624                                        parameters[whichParam(CLP_PARAM_DBL_OBJSCALE2, numberParameters, parameters)].doubleValue();
625                                   if (objScale != 1.0) {
626                                        int iColumn;
627                                        int numberColumns = models[iModel].numberColumns();
628                                        double * dualColumnSolution =
629                                             models[iModel].dualColumnSolution();
630                                        ClpObjective * obj = models[iModel].objectiveAsObject();
631                                        assert(dynamic_cast<ClpLinearObjective *> (obj));
632                                        double offset;
633                                        double * objective = obj->gradient(NULL, NULL, offset, true);
634                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
635                                             dualColumnSolution[iColumn] *= objScale;
636                                             objective[iColumn] *= objScale;;
637                                        }
638                                        int iRow;
639                                        int numberRows = models[iModel].numberRows();
640                                        double * dualRowSolution =
641                                             models[iModel].dualRowSolution();
642                                        for (iRow = 0; iRow < numberRows; iRow++)
643                                             dualRowSolution[iRow] *= objScale;
644                                        models[iModel].setObjectiveOffset(objScale * models[iModel].objectiveOffset());
645                                   }
646                                   ClpSolve::SolveType method;
647                                   ClpSolve::PresolveType presolveType;
648                                   ClpSimplex * model2 = models + iModel;
649                                   if (dualize) {
650                                        bool tryIt = true;
651                                        double fractionColumn = 1.0;
652                                        double fractionRow = 1.0;
653                                        if (dualize == 3) {
654                                             dualize = 1;
655                                             int numberColumns = model2->numberColumns();
656                                             int numberRows = model2->numberRows();
657                                             if (numberRows < 50000 || 5 * numberColumns > numberRows) {
658                                                  tryIt = false;
659                                             } else {
660                                                  fractionColumn = 0.1;
661                                                  fractionRow = 0.1;
662                                             }
663                                        }
664                                        if (tryIt) {
665                                             model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(fractionRow, fractionColumn);
666                                             if (model2) {
667                                                  printf("Dual of model has %d rows and %d columns\n",
668                                                         model2->numberRows(), model2->numberColumns());
669                                                  model2->setOptimizationDirection(1.0);
670                                             } else {
671                                                  model2 = models + iModel;
672                                                  dualize = 0;
673                                             }
674                                        } else {
675                                             dualize = 0;
676                                        }
677                                   }
678                                   ClpSolve solveOptions;
679                                   solveOptions.setPresolveActions(presolveOptions);
680                                   solveOptions.setSubstitution(substitution);
681                                   if (preSolve != 5 && preSolve) {
682                                        presolveType = ClpSolve::presolveNumber;
683                                        if (preSolve < 0) {
684                                             preSolve = - preSolve;
685                                             if (preSolve <= 100) {
686                                                  presolveType = ClpSolve::presolveNumber;
687                                                  printf("Doing %d presolve passes - picking up non-costed slacks\n",
688                                                         preSolve);
689                                                  solveOptions.setDoSingletonColumn(true);
690                                             } else {
691                                                  preSolve -= 100;
692                                                  presolveType = ClpSolve::presolveNumberCost;
693                                                  printf("Doing %d presolve passes - picking up costed slacks\n",
694                                                         preSolve);
695                                             }
696                                        }
697                                   } else if (preSolve) {
698                                        presolveType = ClpSolve::presolveOn;
699                                   } else {
700                                        presolveType = ClpSolve::presolveOff;
701                                   }
702                                   solveOptions.setPresolveType(presolveType, preSolve);
703                                   if (type == CLP_PARAM_ACTION_DUALSIMPLEX ||
704                                             type == CBC_PARAM_ACTION_BAB) {
705                                        method = ClpSolve::useDual;
706                                   } else if (type == CLP_PARAM_ACTION_PRIMALSIMPLEX) {
707                                        method = ClpSolve::usePrimalorSprint;
708                                   } else if (type == CLP_PARAM_ACTION_EITHERSIMPLEX) {
709                                        method = ClpSolve::automatic;
710                                   } else {
711                                        method = ClpSolve::useBarrier;
712                                        if (crossover == 1) {
713                                             method = ClpSolve::useBarrierNoCross;
714                                        } else if (crossover == 2) {
715                                             ClpObjective * obj = models[iModel].objectiveAsObject();
716                                             if (obj->type() > 1) {
717                                                  method = ClpSolve::useBarrierNoCross;
718                                                  presolveType = ClpSolve::presolveOff;
719                                                  solveOptions.setPresolveType(presolveType, preSolve);
720                                             }
721                                        }
722                                   }
723                                   solveOptions.setSolveType(method);
724                                   if (preSolveFile)
725                                        presolveOptions |= 0x40000000;
726                                   solveOptions.setSpecialOption(4, presolveOptions);
727                                   solveOptions.setSpecialOption(5, printOptions & 1);
728                                   if (doVector) {
729                                        ClpMatrixBase * matrix = models[iModel].clpMatrix();
730                                        if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
731                                             ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
732                                             clpMatrix->makeSpecialColumnCopy();
733                                        }
734                                   }
735                                   if (method == ClpSolve::useDual) {
736                                        // dual
737                                        if (doCrash)
738                                             solveOptions.setSpecialOption(0, 1, doCrash); // crash
739                                        else if (doIdiot)
740                                             solveOptions.setSpecialOption(0, 2, doIdiot); // possible idiot
741                                   } else if (method == ClpSolve::usePrimalorSprint) {
742                                        // primal
743                                        // if slp turn everything off
744                                        if (slpValue > 0) {
745                                             doCrash = false;
746                                             doSprint = 0;
747                                             doIdiot = -1;
748                                             solveOptions.setSpecialOption(1, 10, slpValue); // slp
749                                             method = ClpSolve::usePrimal;
750                                        }
751                                        if (doCrash) {
752                                             solveOptions.setSpecialOption(1, 1, doCrash); // crash
753                                        } else if (doSprint > 0) {
754                                             // sprint overrides idiot
755                                             solveOptions.setSpecialOption(1, 3, doSprint); // sprint
756                                        } else if (doIdiot > 0) {
757                                             solveOptions.setSpecialOption(1, 2, doIdiot); // idiot
758                                        } else if (slpValue <= 0) {
759                                             if (doIdiot == 0) {
760                                                  if (doSprint == 0)
761                                                       solveOptions.setSpecialOption(1, 4); // all slack
762                                                  else
763                                                       solveOptions.setSpecialOption(1, 9); // all slack or sprint
764                                             } else {
765                                                  if (doSprint == 0)
766                                                       solveOptions.setSpecialOption(1, 8); // all slack or idiot
767                                                  else
768                                                       solveOptions.setSpecialOption(1, 7); // initiative
769                                             }
770                                        }
771                                        if (basisHasValues == -1)
772                                             solveOptions.setSpecialOption(1, 11); // switch off values
773                                   } else if (method == ClpSolve::useBarrier || method == ClpSolve::useBarrierNoCross) {
774                                        int barrierOptions = choleskyType;
775                                        if (scaleBarrier) {
776                                             if ((scaleBarrier & 1) != 0)
777                                                  barrierOptions |= 8;
778                                             barrierOptions |= 2048 * (scaleBarrier >> 1);
779                                        }
780                                        if (doKKT)
781                                             barrierOptions |= 16;
782                                        if (gamma)
783                                             barrierOptions |= 32 * gamma;
784                                        if (crossover == 3)
785                                             barrierOptions |= 256; // try presolve in crossover
786                                        solveOptions.setSpecialOption(4, barrierOptions);
787                                   }
788                                   int status;
789                                   if (cppValue >= 0) {
790                                        // generate code
791                                        FILE * fp = fopen("user_driver.cpp", "w");
792                                        if (fp) {
793                                             // generate enough to do solveOptions
794                                             model2->generateCpp(fp);
795                                             solveOptions.generateCpp(fp);
796                                             fclose(fp);
797                                             // now call generate code
798                                             generateCode("user_driver.cpp", cppValue);
799                                        } else {
800                                             std::cout << "Unable to open file user_driver.cpp" << std::endl;
801                                        }
802                                   }
803#ifdef CLP_MULTIPLE_FACTORIZATIONS
804                                   int denseCode = parameters[whichParam(CBC_PARAM_INT_DENSE, numberParameters, parameters)].intValue();
805                                   model2->factorization()->setGoDenseThreshold(denseCode);
806                                   int smallCode = parameters[whichParam(CBC_PARAM_INT_SMALLFACT, numberParameters, parameters)].intValue();
807                                   model2->factorization()->setGoSmallThreshold(smallCode);
808                                   model2->factorization()->goDenseOrSmall(model2->numberRows());
809#endif
810                                   try {
811                                        status = model2->initialSolve(solveOptions);
812                                   } catch (CoinError e) {
813                                        e.print();
814                                        status = -1;
815                                   }
816                                   if (dualize) {
817                                        int returnCode = static_cast<ClpSimplexOther *> (models + iModel)->restoreFromDual(model2);
818                                        if (model2->status() == 3)
819                                             returnCode = 0;
820                                        delete model2;
821                                        if (returnCode && dualize != 2) {
822                                             currentModel = models + iModel;
823                                             // register signal handler
824                                             signal(SIGINT, signal_handler);
825                                             models[iModel].primal(1);
826                                             currentModel = NULL;
827                                        }
828                                   }
829                                   if (status >= 0)
830                                        basisHasValues = 1;
831                              } else {
832                                   std::cout << "** Current model not valid" << std::endl;
833                              }
834                              break;
835                         case CLP_PARAM_ACTION_STATISTICS:
836                              if (goodModels[iModel]) {
837                                   // If presolve on look at presolved
838                                   bool deleteModel2 = false;
839                                   ClpSimplex * model2 = models + iModel;
840                                   if (preSolve) {
841                                        ClpPresolve pinfo;
842                                        int presolveOptions2 = presolveOptions&~0x40000000;
843                                        if ((presolveOptions2 & 0xffff) != 0)
844                                             pinfo.setPresolveActions(presolveOptions2);
845                                        pinfo.setSubstitution(substitution);
846                                        if ((printOptions & 1) != 0)
847                                             pinfo.statistics();
848                                        double presolveTolerance =
849                                             parameters[whichParam(CLP_PARAM_DBL_PRESOLVETOLERANCE, numberParameters, parameters)].doubleValue();
850                                        model2 =
851                                             pinfo.presolvedModel(models[iModel], presolveTolerance,
852                                                                  true, preSolve);
853                                        if (model2) {
854                                             printf("Statistics for presolved model\n");
855                                             deleteModel2 = true;
856                                        } else {
857                                             printf("Presolved model looks infeasible - will use unpresolved\n");
858                                             model2 = models + iModel;
859                                        }
860                                   } else {
861                                        printf("Statistics for unpresolved model\n");
862                                        model2 =  models + iModel;
863                                   }
864                                   statistics(models + iModel, model2);
865                                   if (deleteModel2)
866                                        delete model2;
867                              } else {
868                                   std::cout << "** Current model not valid" << std::endl;
869                              }
870                              break;
871                         case CLP_PARAM_ACTION_TIGHTEN:
872                              if (goodModels[iModel]) {
873                                   int numberInfeasibilities = models[iModel].tightenPrimalBounds();
874                                   if (numberInfeasibilities)
875                                        std::cout << "** Analysis indicates model infeasible" << std::endl;
876                              } else {
877                                   std::cout << "** Current model not valid" << std::endl;
878                              }
879                              break;
880                         case CLP_PARAM_ACTION_PLUSMINUS:
881                              if (goodModels[iModel]) {
882                                   ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
883                                   ClpPackedMatrix* clpMatrix =
884                                        dynamic_cast< ClpPackedMatrix*>(saveMatrix);
885                                   if (clpMatrix) {
886                                        ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
887                                        if (newMatrix->getIndices()) {
888                                             models[iModel].replaceMatrix(newMatrix);
889                                             delete saveMatrix;
890                                             std::cout << "Matrix converted to +- one matrix" << std::endl;
891                                        } else {
892                                             std::cout << "Matrix can not be converted to +- 1 matrix" << std::endl;
893                                        }
894                                   } else {
895                                        std::cout << "Matrix not a ClpPackedMatrix" << std::endl;
896                                   }
897                              } else {
898                                   std::cout << "** Current model not valid" << std::endl;
899                              }
900                              break;
901                         case CLP_PARAM_ACTION_NETWORK:
902                              if (goodModels[iModel]) {
903                                   ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
904                                   ClpPackedMatrix* clpMatrix =
905                                        dynamic_cast< ClpPackedMatrix*>(saveMatrix);
906                                   if (clpMatrix) {
907                                        ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
908                                        if (newMatrix->getIndices()) {
909                                             models[iModel].replaceMatrix(newMatrix);
910                                             delete saveMatrix;
911                                             std::cout << "Matrix converted to network matrix" << std::endl;
912                                        } else {
913                                             std::cout << "Matrix can not be converted to network matrix" << std::endl;
914                                        }
915                                   } else {
916                                        std::cout << "Matrix not a ClpPackedMatrix" << std::endl;
917                                   }
918                              } else {
919                                   std::cout << "** Current model not valid" << std::endl;
920                              }
921                              break;
922                         case CLP_PARAM_ACTION_IMPORT: {
923                              // get next field
924                              field = CoinReadGetString(argc, argv);
925                              if (field == "$") {
926                                   field = parameters[iParam].stringValue();
927                              } else if (field == "EOL") {
928                                   parameters[iParam].printString();
929                                   break;
930                              } else {
931                                   parameters[iParam].setStringValue(field);
932                              }
933                              std::string fileName;
934                              bool canOpen = false;
935                              // See if gmpl file
936                              int gmpl = 0;
937                              std::string gmplData;
938                              if (field == "-") {
939                                   // stdin
940                                   canOpen = true;
941                                   fileName = "-";
942                              } else {
943                                   // See if .lp
944                                   {
945                                        const char * c_name = field.c_str();
946                                        int length = strlen(c_name);
947                                        if (length > 3 && !strncmp(c_name + length - 3, ".lp", 3))
948                                             gmpl = -1; // .lp
949                                   }
950                                   bool absolutePath;
951                                   if (dirsep == '/') {
952                                        // non Windows (or cygwin)
953                                        absolutePath = (field[0] == '/');
954                                   } else {
955                                        //Windows (non cycgwin)
956                                        absolutePath = (field[0] == '\\');
957                                        // but allow for :
958                                        if (strchr(field.c_str(), ':'))
959                                             absolutePath = true;
960                                   }
961                                   if (absolutePath) {
962                                        fileName = field;
963                                   } else if (field[0] == '~') {
964                                        char * environVar = getenv("HOME");
965                                        if (environVar) {
966                                             std::string home(environVar);
967                                             field = field.erase(0, 1);
968                                             fileName = home + field;
969                                        } else {
970                                             fileName = field;
971                                        }
972                                   } else {
973                                        fileName = directory + field;
974                                        // See if gmpl (model & data) - or even lp file
975                                        int length = field.size();
976                                        int percent = field.find('%');
977                                        if (percent < length && percent > 0) {
978                                             gmpl = 1;
979                                             fileName = directory + field.substr(0, percent);
980                                             gmplData = directory + field.substr(percent + 1);
981                                             if (percent < length - 1)
982                                                  gmpl = 2; // two files
983                                             printf("GMPL model file %s and data file %s\n",
984                                                    fileName.c_str(), gmplData.c_str());
985                                        }
986                                   }
987                                   std::string name = fileName;
988                                   if (fileCoinReadable(name)) {
989                                        // can open - lets go for it
990                                        canOpen = true;
991                                        if (gmpl == 2) {
992                                             FILE *fp;
993                                             fp = fopen(gmplData.c_str(), "r");
994                                             if (fp) {
995                                                  fclose(fp);
996                                             } else {
997                                                  canOpen = false;
998                                                  std::cout << "Unable to open file " << gmplData << std::endl;
999                                             }
1000                                        }
1001                                   } else {
1002                                        std::cout << "Unable to open file " << fileName << std::endl;
1003                                   }
1004                              }
1005                              if (canOpen) {
1006                                   int status;
1007                                   if (!gmpl)
1008                                        status = models[iModel].readMps(fileName.c_str(),
1009                                                                        keepImportNames != 0,
1010                                                                        allowImportErrors != 0);
1011                                   else if (gmpl > 0)
1012                                        status = models[iModel].readGMPL(fileName.c_str(),
1013                                                                         (gmpl == 2) ? gmplData.c_str() : NULL,
1014                                                                         keepImportNames != 0);
1015                                   else
1016                                        status = models[iModel].readLp(fileName.c_str(), 1.0e-12);
1017                                   if (!status || (status > 0 && allowImportErrors)) {
1018                                        goodModels[iModel] = true;
1019                                        // sets to all slack (not necessary?)
1020                                        models[iModel].createStatus();
1021                                        time2 = CoinCpuTime();
1022                                        totalTime += time2 - time1;
1023                                        time1 = time2;
1024                                        // Go to canned file if just input file
1025                                        if (CbcOrClpRead_mode == 2 && argc == 2) {
1026                                             // only if ends .mps
1027                                             char * find = const_cast<char *>(strstr(fileName.c_str(), ".mps"));
1028                                             if (find && find[4] == '\0') {
1029                                                  find[1] = 'p';
1030                                                  find[2] = 'a';
1031                                                  find[3] = 'r';
1032                                                  FILE *fp = fopen(fileName.c_str(), "r");
1033                                                  if (fp) {
1034                                                       CbcOrClpReadCommand = fp; // Read from that file
1035                                                       CbcOrClpRead_mode = -1;
1036                                                  }
1037                                             }
1038                                        }
1039                                   } else {
1040                                        // errors
1041                                        std::cout << "There were " << status <<
1042                                                  " errors on input" << std::endl;
1043                                   }
1044                              }
1045                         }
1046                         break;
1047                         case CLP_PARAM_ACTION_EXPORT:
1048                              if (goodModels[iModel]) {
1049                                   double objScale =
1050                                        parameters[whichParam(CLP_PARAM_DBL_OBJSCALE2, numberParameters, parameters)].doubleValue();
1051                                   if (objScale != 1.0) {
1052                                        int iColumn;
1053                                        int numberColumns = models[iModel].numberColumns();
1054                                        double * dualColumnSolution =
1055                                             models[iModel].dualColumnSolution();
1056                                        ClpObjective * obj = models[iModel].objectiveAsObject();
1057                                        assert(dynamic_cast<ClpLinearObjective *> (obj));
1058                                        double offset;
1059                                        double * objective = obj->gradient(NULL, NULL, offset, true);
1060                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1061                                             dualColumnSolution[iColumn] *= objScale;
1062                                             objective[iColumn] *= objScale;;
1063                                        }
1064                                        int iRow;
1065                                        int numberRows = models[iModel].numberRows();
1066                                        double * dualRowSolution =
1067                                             models[iModel].dualRowSolution();
1068                                        for (iRow = 0; iRow < numberRows; iRow++)
1069                                             dualRowSolution[iRow] *= objScale;
1070                                        models[iModel].setObjectiveOffset(objScale * models[iModel].objectiveOffset());
1071                                   }
1072                                   // get next field
1073                                   field = CoinReadGetString(argc, argv);
1074                                   if (field == "$") {
1075                                        field = parameters[iParam].stringValue();
1076                                   } else if (field == "EOL") {
1077                                        parameters[iParam].printString();
1078                                        break;
1079                                   } else {
1080                                        parameters[iParam].setStringValue(field);
1081                                   }
1082                                   std::string fileName;
1083                                   bool canOpen = false;
1084                                   if (field[0] == '/' || field[0] == '\\') {
1085                                        fileName = field;
1086                                   } else if (field[0] == '~') {
1087                                        char * environVar = getenv("HOME");
1088                                        if (environVar) {
1089                                             std::string home(environVar);
1090                                             field = field.erase(0, 1);
1091                                             fileName = home + field;
1092                                        } else {
1093                                             fileName = field;
1094                                        }
1095                                   } else {
1096                                        fileName = directory + field;
1097                                   }
1098                                   FILE *fp = fopen(fileName.c_str(), "w");
1099                                   if (fp) {
1100                                        // can open - lets go for it
1101                                        fclose(fp);
1102                                        canOpen = true;
1103                                   } else {
1104                                        std::cout << "Unable to open file " << fileName << std::endl;
1105                                   }
1106                                   if (canOpen) {
1107                                        // If presolve on then save presolved
1108                                        bool deleteModel2 = false;
1109                                        ClpSimplex * model2 = models + iModel;
1110                                        if (dualize && dualize < 3) {
1111                                             model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel();
1112                                             printf("Dual of model has %d rows and %d columns\n",
1113                                                    model2->numberRows(), model2->numberColumns());
1114                                             model2->setOptimizationDirection(1.0);
1115                                             preSolve = 0; // as picks up from model
1116                                        }
1117                                        if (preSolve) {
1118                                             ClpPresolve pinfo;
1119                                             int presolveOptions2 = presolveOptions&~0x40000000;
1120                                             if ((presolveOptions2 & 0xffff) != 0)
1121                                                  pinfo.setPresolveActions(presolveOptions2);
1122                                             pinfo.setSubstitution(substitution);
1123                                             if ((printOptions & 1) != 0)
1124                                                  pinfo.statistics();
1125                                             double presolveTolerance =
1126                                                  parameters[whichParam(CLP_PARAM_DBL_PRESOLVETOLERANCE, numberParameters, parameters)].doubleValue();
1127                                             model2 =
1128                                                  pinfo.presolvedModel(models[iModel], presolveTolerance,
1129                                                                       true, preSolve, false, false);
1130                                             if (model2) {
1131                                                  printf("Saving presolved model on %s\n",
1132                                                         fileName.c_str());
1133                                                  deleteModel2 = true;
1134                                             } else {
1135                                                  printf("Presolved model looks infeasible - saving original on %s\n",
1136                                                         fileName.c_str());
1137                                                  deleteModel2 = false;
1138                                                  model2 = models + iModel;
1139
1140                                             }
1141                                        } else {
1142                                             printf("Saving model on %s\n",
1143                                                    fileName.c_str());
1144                                        }
1145#if 0
1146                                        // Convert names
1147                                        int iRow;
1148                                        int numberRows = model2->numberRows();
1149                                        int iColumn;
1150                                        int numberColumns = model2->numberColumns();
1151
1152                                        char ** rowNames = NULL;
1153                                        char ** columnNames = NULL;
1154                                        if (model2->lengthNames()) {
1155                                             rowNames = new char * [numberRows];
1156                                             for (iRow = 0; iRow < numberRows; iRow++) {
1157                                                  rowNames[iRow] =
1158                                                       CoinStrdup(model2->rowName(iRow).c_str());
1159#ifdef STRIPBLANKS
1160                                                  char * xx = rowNames[iRow];
1161                                                  int i;
1162                                                  int length = strlen(xx);
1163                                                  int n = 0;
1164                                                  for (i = 0; i < length; i++) {
1165                                                       if (xx[i] != ' ')
1166                                                            xx[n++] = xx[i];
1167                                                  }
1168                                                  xx[n] = '\0';
1169#endif
1170                                             }
1171
1172                                             columnNames = new char * [numberColumns];
1173                                             for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1174                                                  columnNames[iColumn] =
1175                                                       CoinStrdup(model2->columnName(iColumn).c_str());
1176#ifdef STRIPBLANKS
1177                                                  char * xx = columnNames[iColumn];
1178                                                  int i;
1179                                                  int length = strlen(xx);
1180                                                  int n = 0;
1181                                                  for (i = 0; i < length; i++) {
1182                                                       if (xx[i] != ' ')
1183                                                            xx[n++] = xx[i];
1184                                                  }
1185                                                  xx[n] = '\0';
1186#endif
1187                                             }
1188                                        }
1189                                        CoinMpsIO writer;
1190                                        writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
1191                                                          model2->getColLower(), model2->getColUpper(),
1192                                                          model2->getObjCoefficients(),
1193                                                          (const char*) 0 /*integrality*/,
1194                                                          model2->getRowLower(), model2->getRowUpper(),
1195                                                          columnNames, rowNames);
1196                                        // Pass in array saying if each variable integer
1197                                        writer.copyInIntegerInformation(model2->integerInformation());
1198                                        writer.setObjectiveOffset(model2->objectiveOffset());
1199                                        writer.writeMps(fileName.c_str(), 0, 1, 1);
1200                                        if (rowNames) {
1201                                             for (iRow = 0; iRow < numberRows; iRow++) {
1202                                                  free(rowNames[iRow]);
1203                                             }
1204                                             delete [] rowNames;
1205                                             for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1206                                                  free(columnNames[iColumn]);
1207                                             }
1208                                             delete [] columnNames;
1209                                        }
1210#else
1211                                        model2->writeMps(fileName.c_str(), (outputFormat - 1) / 2, 1 + ((outputFormat - 1) & 1));
1212#endif
1213                                        if (deleteModel2)
1214                                             delete model2;
1215                                        time2 = CoinCpuTime();
1216                                        totalTime += time2 - time1;
1217                                        time1 = time2;
1218                                   }
1219                              } else {
1220                                   std::cout << "** Current model not valid" << std::endl;
1221                              }
1222                              break;
1223                         case CLP_PARAM_ACTION_BASISIN:
1224                              if (goodModels[iModel]) {
1225                                   // get next field
1226                                   field = CoinReadGetString(argc, argv);
1227                                   if (field == "$") {
1228                                        field = parameters[iParam].stringValue();
1229                                   } else if (field == "EOL") {
1230                                        parameters[iParam].printString();
1231                                        break;
1232                                   } else {
1233                                        parameters[iParam].setStringValue(field);
1234                                   }
1235                                   std::string fileName;
1236                                   bool canOpen = false;
1237                                   if (field == "-") {
1238                                        // stdin
1239                                        canOpen = true;
1240                                        fileName = "-";
1241                                   } else {
1242                                        if (field[0] == '/' || field[0] == '\\') {
1243                                             fileName = field;
1244                                        } else if (field[0] == '~') {
1245                                             char * environVar = getenv("HOME");
1246                                             if (environVar) {
1247                                                  std::string home(environVar);
1248                                                  field = field.erase(0, 1);
1249                                                  fileName = home + field;
1250                                             } else {
1251                                                  fileName = field;
1252                                             }
1253                                        } else {
1254                                             fileName = directory + field;
1255                                        }
1256                                        FILE *fp = fopen(fileName.c_str(), "r");
1257                                        if (fp) {
1258                                             // can open - lets go for it
1259                                             fclose(fp);
1260                                             canOpen = true;
1261                                        } else {
1262                                             std::cout << "Unable to open file " << fileName << std::endl;
1263                                        }
1264                                   }
1265                                   if (canOpen) {
1266                                        int values = models[iModel].readBasis(fileName.c_str());
1267                                        if (values == 0)
1268                                             basisHasValues = -1;
1269                                        else
1270                                             basisHasValues = 1;
1271                                   }
1272                              } else {
1273                                   std::cout << "** Current model not valid" << std::endl;
1274                              }
1275                              break;
1276                         case CLP_PARAM_ACTION_PRINTMASK:
1277                              // get next field
1278                         {
1279                              std::string name = CoinReadGetString(argc, argv);
1280                              if (name != "EOL") {
1281                                   parameters[iParam].setStringValue(name);
1282                                   printMask = name;
1283                              } else {
1284                                   parameters[iParam].printString();
1285                              }
1286                         }
1287                         break;
1288                         case CLP_PARAM_ACTION_BASISOUT:
1289                              if (goodModels[iModel]) {
1290                                   // get next field
1291                                   field = CoinReadGetString(argc, argv);
1292                                   if (field == "$") {
1293                                        field = parameters[iParam].stringValue();
1294                                   } else if (field == "EOL") {
1295                                        parameters[iParam].printString();
1296                                        break;
1297                                   } else {
1298                                        parameters[iParam].setStringValue(field);
1299                                   }
1300                                   std::string fileName;
1301                                   bool canOpen = false;
1302                                   if (field[0] == '/' || field[0] == '\\') {
1303                                        fileName = field;
1304                                   } else if (field[0] == '~') {
1305                                        char * environVar = getenv("HOME");
1306                                        if (environVar) {
1307                                             std::string home(environVar);
1308                                             field = field.erase(0, 1);
1309                                             fileName = home + field;
1310                                        } else {
1311                                             fileName = field;
1312                                        }
1313                                   } else {
1314                                        fileName = directory + field;
1315                                   }
1316                                   FILE *fp = fopen(fileName.c_str(), "w");
1317                                   if (fp) {
1318                                        // can open - lets go for it
1319                                        fclose(fp);
1320                                        canOpen = true;
1321                                   } else {
1322                                        std::cout << "Unable to open file " << fileName << std::endl;
1323                                   }
1324                                   if (canOpen) {
1325                                        ClpSimplex * model2 = models + iModel;
1326                                        model2->writeBasis(fileName.c_str(), outputFormat > 1, outputFormat - 2);
1327                                        time2 = CoinCpuTime();
1328                                        totalTime += time2 - time1;
1329                                        time1 = time2;
1330                                   }
1331                              } else {
1332                                   std::cout << "** Current model not valid" << std::endl;
1333                              }
1334                              break;
1335                         case CLP_PARAM_ACTION_PARAMETRICS:
1336                              if (goodModels[iModel]) {
1337                                   // get next field
1338                                   field = CoinReadGetString(argc, argv);
1339                                   if (field == "$") {
1340                                        field = parameters[iParam].stringValue();
1341                                   } else if (field == "EOL") {
1342                                        parameters[iParam].printString();
1343                                        break;
1344                                   } else {
1345                                        parameters[iParam].setStringValue(field);
1346                                   }
1347                                   std::string fileName;
1348                                   //bool canOpen = false;
1349                                   if (field[0] == '/' || field[0] == '\\') {
1350                                        fileName = field;
1351                                   } else if (field[0] == '~') {
1352                                        char * environVar = getenv("HOME");
1353                                        if (environVar) {
1354                                             std::string home(environVar);
1355                                             field = field.erase(0, 1);
1356                                             fileName = home + field;
1357                                        } else {
1358                                             fileName = field;
1359                                        }
1360                                   } else {
1361                                        fileName = directory + field;
1362                                   }
1363                                   ClpSimplex * model2 = models + iModel;
1364                                   static_cast<ClpSimplexOther *> (model2)->parametrics(fileName.c_str());
1365                                   time2 = CoinCpuTime();
1366                                   totalTime += time2 - time1;
1367                                   time1 = time2;
1368                              } else {
1369                                   std::cout << "** Current model not valid" << std::endl;
1370                              }
1371                              break;
1372                         case CLP_PARAM_ACTION_SAVE: {
1373                              // get next field
1374                              field = CoinReadGetString(argc, argv);
1375                              if (field == "$") {
1376                                   field = parameters[iParam].stringValue();
1377                              } else if (field == "EOL") {
1378                                   parameters[iParam].printString();
1379                                   break;
1380                              } else {
1381                                   parameters[iParam].setStringValue(field);
1382                              }
1383                              std::string fileName;
1384                              bool canOpen = false;
1385                              if (field[0] == '/' || field[0] == '\\') {
1386                                   fileName = field;
1387                              } else if (field[0] == '~') {
1388                                   char * environVar = getenv("HOME");
1389                                   if (environVar) {
1390                                        std::string home(environVar);
1391                                        field = field.erase(0, 1);
1392                                        fileName = home + field;
1393                                   } else {
1394                                        fileName = field;
1395                                   }
1396                              } else {
1397                                   fileName = directory + field;
1398                              }
1399                              FILE *fp = fopen(fileName.c_str(), "wb");
1400                              if (fp) {
1401                                   // can open - lets go for it
1402                                   fclose(fp);
1403                                   canOpen = true;
1404                              } else {
1405                                   std::cout << "Unable to open file " << fileName << std::endl;
1406                              }
1407                              if (canOpen) {
1408                                   int status;
1409                                   // If presolve on then save presolved
1410                                   bool deleteModel2 = false;
1411                                   ClpSimplex * model2 = models + iModel;
1412                                   if (preSolve) {
1413                                        ClpPresolve pinfo;
1414                                        double presolveTolerance =
1415                                             parameters[whichParam(CLP_PARAM_DBL_PRESOLVETOLERANCE, numberParameters, parameters)].doubleValue();
1416                                        model2 =
1417                                             pinfo.presolvedModel(models[iModel], presolveTolerance,
1418                                                                  false, preSolve);
1419                                        if (model2) {
1420                                             printf("Saving presolved model on %s\n",
1421                                                    fileName.c_str());
1422                                             deleteModel2 = true;
1423                                        } else {
1424                                             printf("Presolved model looks infeasible - saving original on %s\n",
1425                                                    fileName.c_str());
1426                                             deleteModel2 = false;
1427                                             model2 = models + iModel;
1428
1429                                        }
1430                                   } else {
1431                                        printf("Saving model on %s\n",
1432                                               fileName.c_str());
1433                                   }
1434                                   status = model2->saveModel(fileName.c_str());
1435                                   if (deleteModel2)
1436                                        delete model2;
1437                                   if (!status) {
1438                                        goodModels[iModel] = true;
1439                                        time2 = CoinCpuTime();
1440                                        totalTime += time2 - time1;
1441                                        time1 = time2;
1442                                   } else {
1443                                        // errors
1444                                        std::cout << "There were errors on output" << std::endl;
1445                                   }
1446                              }
1447                         }
1448                         break;
1449                         case CLP_PARAM_ACTION_RESTORE: {
1450                              // get next field
1451                              field = CoinReadGetString(argc, argv);
1452                              if (field == "$") {
1453                                   field = parameters[iParam].stringValue();
1454                              } else if (field == "EOL") {
1455                                   parameters[iParam].printString();
1456                                   break;
1457                              } else {
1458                                   parameters[iParam].setStringValue(field);
1459                              }
1460                              std::string fileName;
1461                              bool canOpen = false;
1462                              if (field[0] == '/' || field[0] == '\\') {
1463                                   fileName = field;
1464                              } else if (field[0] == '~') {
1465                                   char * environVar = getenv("HOME");
1466                                   if (environVar) {
1467                                        std::string home(environVar);
1468                                        field = field.erase(0, 1);
1469                                        fileName = home + field;
1470                                   } else {
1471                                        fileName = field;
1472                                   }
1473                              } else {
1474                                   fileName = directory + field;
1475                              }
1476                              FILE *fp = fopen(fileName.c_str(), "rb");
1477                              if (fp) {
1478                                   // can open - lets go for it
1479                                   fclose(fp);
1480                                   canOpen = true;
1481                              } else {
1482                                   std::cout << "Unable to open file " << fileName << std::endl;
1483                              }
1484                              if (canOpen) {
1485                                   int status = models[iModel].restoreModel(fileName.c_str());
1486                                   if (!status) {
1487                                        goodModels[iModel] = true;
1488                                        time2 = CoinCpuTime();
1489                                        totalTime += time2 - time1;
1490                                        time1 = time2;
1491                                   } else {
1492                                        // errors
1493                                        std::cout << "There were errors on input" << std::endl;
1494                                   }
1495                              }
1496                         }
1497                         break;
1498                         case CLP_PARAM_ACTION_MAXIMIZE:
1499                              models[iModel].setOptimizationDirection(-1);
1500                              break;
1501                         case CLP_PARAM_ACTION_MINIMIZE:
1502                              models[iModel].setOptimizationDirection(1);
1503                              break;
1504                         case CLP_PARAM_ACTION_ALLSLACK:
1505                              models[iModel].allSlackBasis(true);
1506                              break;
1507                         case CLP_PARAM_ACTION_REVERSE:
1508                              if (goodModels[iModel]) {
1509                                   int iColumn;
1510                                   int numberColumns = models[iModel].numberColumns();
1511                                   double * dualColumnSolution =
1512                                        models[iModel].dualColumnSolution();
1513                                   ClpObjective * obj = models[iModel].objectiveAsObject();
1514                                   assert(dynamic_cast<ClpLinearObjective *> (obj));
1515                                   double offset;
1516                                   double * objective = obj->gradient(NULL, NULL, offset, true);
1517                                   for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1518                                        dualColumnSolution[iColumn] = -dualColumnSolution[iColumn];
1519                                        objective[iColumn] = -objective[iColumn];
1520                                   }
1521                                   int iRow;
1522                                   int numberRows = models[iModel].numberRows();
1523                                   double * dualRowSolution =
1524                                        models[iModel].dualRowSolution();
1525                                   for (iRow = 0; iRow < numberRows; iRow++) {
1526                                        dualRowSolution[iRow] = -dualRowSolution[iRow];
1527                                   }
1528                                   models[iModel].setObjectiveOffset(-models[iModel].objectiveOffset());
1529                              }
1530                              break;
1531                         case CLP_PARAM_ACTION_DIRECTORY: {
1532                              std::string name = CoinReadGetString(argc, argv);
1533                              if (name != "EOL") {
1534                                   int length = name.length();
1535                                   if (name[length-1] == dirsep) {
1536                                        directory = name;
1537                                   } else {
1538                                        directory = name + dirsep;
1539                                   }
1540                                   parameters[iParam].setStringValue(directory);
1541                              } else {
1542                                   parameters[iParam].printString();
1543                              }
1544                         }
1545                         break;
1546                         case CLP_PARAM_ACTION_DIRSAMPLE: {
1547                              std::string name = CoinReadGetString(argc, argv);
1548                              if (name != "EOL") {
1549                                   int length = name.length();
1550                                   if (name[length-1] == dirsep) {
1551                                        dirSample = name;
1552                                   } else {
1553                                        dirSample = name + dirsep;
1554                                   }
1555                                   parameters[iParam].setStringValue(dirSample);
1556                              } else {
1557                                   parameters[iParam].printString();
1558                              }
1559                         }
1560                         break;
1561                         case CLP_PARAM_ACTION_DIRNETLIB: {
1562                              std::string name = CoinReadGetString(argc, argv);
1563                              if (name != "EOL") {
1564                                   int length = name.length();
1565                                   if (name[length-1] == dirsep) {
1566                                        dirNetlib = name;
1567                                   } else {
1568                                        dirNetlib = name + dirsep;
1569                                   }
1570                                   parameters[iParam].setStringValue(dirNetlib);
1571                              } else {
1572                                   parameters[iParam].printString();
1573                              }
1574                         }
1575                         break;
1576                         case CBC_PARAM_ACTION_DIRMIPLIB: {
1577                              std::string name = CoinReadGetString(argc, argv);
1578                              if (name != "EOL") {
1579                                   int length = name.length();
1580                                   if (name[length-1] == dirsep) {
1581                                        dirMiplib = name;
1582                                   } else {
1583                                        dirMiplib = name + dirsep;
1584                                   }
1585                                   parameters[iParam].setStringValue(dirMiplib);
1586                              } else {
1587                                   parameters[iParam].printString();
1588                              }
1589                         }
1590                         break;
1591                         case CLP_PARAM_ACTION_STDIN:
1592                              CbcOrClpRead_mode = -1;
1593                              break;
1594                         case CLP_PARAM_ACTION_NETLIB_DUAL:
1595                         case CLP_PARAM_ACTION_NETLIB_EITHER:
1596                         case CLP_PARAM_ACTION_NETLIB_BARRIER:
1597                         case CLP_PARAM_ACTION_NETLIB_PRIMAL:
1598                         case CLP_PARAM_ACTION_NETLIB_TUNE: {
1599                              // create fields for unitTest
1600                              const char * fields[4];
1601                              int nFields = 4;
1602                              fields[0] = "fake main from unitTest";
1603                              std::string mpsfield = "-dirSample=";
1604                              mpsfield += dirSample.c_str();
1605                              fields[1] = mpsfield.c_str();
1606                              std::string netfield = "-dirNetlib=";
1607                              netfield += dirNetlib.c_str();
1608                              fields[2] = netfield.c_str();
1609                              fields[3] = "-netlib";
1610                              int algorithm;
1611                              if (type == CLP_PARAM_ACTION_NETLIB_DUAL) {
1612                                   std::cerr << "Doing netlib with dual algorithm" << std::endl;
1613                                   algorithm = 0;
1614                              } else if (type == CLP_PARAM_ACTION_NETLIB_BARRIER) {
1615                                   std::cerr << "Doing netlib with barrier algorithm" << std::endl;
1616                                   algorithm = 2;
1617                              } else if (type == CLP_PARAM_ACTION_NETLIB_EITHER) {
1618                                   std::cerr << "Doing netlib with dual or primal algorithm" << std::endl;
1619                                   algorithm = 3;
1620                              } else if (type == CLP_PARAM_ACTION_NETLIB_TUNE) {
1621                                   std::cerr << "Doing netlib with best algorithm!" << std::endl;
1622                                   algorithm = 5;
1623                                   // uncomment next to get active tuning
1624                                   // algorithm=6;
1625                              } else {
1626                                   std::cerr << "Doing netlib with primal agorithm" << std::endl;
1627                                   algorithm = 1;
1628                              }
1629                              int specialOptions = models[iModel].specialOptions();
1630                              models[iModel].setSpecialOptions(0);
1631                              ClpSolve solveOptions;
1632                              ClpSolve::PresolveType presolveType;
1633                              if (preSolve)
1634                                   presolveType = ClpSolve::presolveOn;
1635                              else
1636                                   presolveType = ClpSolve::presolveOff;
1637                              solveOptions.setPresolveType(presolveType, 5);
1638                              if (doSprint >= 0 || doIdiot >= 0) {
1639                                   if (doSprint > 0) {
1640                                        // sprint overrides idiot
1641                                        solveOptions.setSpecialOption(1, 3, doSprint); // sprint
1642                                   } else if (doIdiot > 0) {
1643                                        solveOptions.setSpecialOption(1, 2, doIdiot); // idiot
1644                                   } else {
1645                                        if (doIdiot == 0) {
1646                                             if (doSprint == 0)
1647                                                  solveOptions.setSpecialOption(1, 4); // all slack
1648                                             else
1649                                                  solveOptions.setSpecialOption(1, 9); // all slack or sprint
1650                                        } else {
1651                                             if (doSprint == 0)
1652                                                  solveOptions.setSpecialOption(1, 8); // all slack or idiot
1653                                             else
1654                                                  solveOptions.setSpecialOption(1, 7); // initiative
1655                                        }
1656                                   }
1657                              }
1658                              mainTest(nFields, fields, algorithm, models[iModel],
1659                                       solveOptions, specialOptions, doVector != 0);
1660                         }
1661                         break;
1662                         case CLP_PARAM_ACTION_UNITTEST: {
1663                              // create fields for unitTest
1664                              const char * fields[2];
1665                              int nFields = 2;
1666                              fields[0] = "fake main from unitTest";
1667                              std::string dirfield = "-dirSample=";
1668                              dirfield += dirSample.c_str();
1669                              fields[1] = dirfield.c_str();
1670                              int specialOptions = models[iModel].specialOptions();
1671                              models[iModel].setSpecialOptions(0);
1672                              int algorithm = -1;
1673                              if (models[iModel].numberRows())
1674                                   algorithm = 7;
1675                              ClpSolve solveOptions;
1676                              ClpSolve::PresolveType presolveType;
1677                              if (preSolve)
1678                                   presolveType = ClpSolve::presolveOn;
1679                              else
1680                                   presolveType = ClpSolve::presolveOff;
1681                              solveOptions.setPresolveType(presolveType, 5);
1682                              mainTest(nFields, fields, algorithm, models[iModel],
1683                                       solveOptions, specialOptions, doVector != 0);
1684                         }
1685                         break;
1686                         case CLP_PARAM_ACTION_FAKEBOUND:
1687                              if (goodModels[iModel]) {
1688                                   // get bound
1689                                   double value = CoinReadGetDoubleField(argc, argv, &valid);
1690                                   if (!valid) {
1691                                        std::cout << "Setting " << parameters[iParam].name() <<
1692                                                  " to DEBUG " << value << std::endl;
1693                                        int iRow;
1694                                        int numberRows = models[iModel].numberRows();
1695                                        double * rowLower = models[iModel].rowLower();
1696                                        double * rowUpper = models[iModel].rowUpper();
1697                                        for (iRow = 0; iRow < numberRows; iRow++) {
1698                                             // leave free ones for now
1699                                             if (rowLower[iRow] > -1.0e20 || rowUpper[iRow] < 1.0e20) {
1700                                                  rowLower[iRow] = CoinMax(rowLower[iRow], -value);
1701                                                  rowUpper[iRow] = CoinMin(rowUpper[iRow], value);
1702                                             }
1703                                        }
1704                                        int iColumn;
1705                                        int numberColumns = models[iModel].numberColumns();
1706                                        double * columnLower = models[iModel].columnLower();
1707                                        double * columnUpper = models[iModel].columnUpper();
1708                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1709                                             // leave free ones for now
1710                                             if (columnLower[iColumn] > -1.0e20 ||
1711                                                       columnUpper[iColumn] < 1.0e20) {
1712                                                  columnLower[iColumn] = CoinMax(columnLower[iColumn], -value);
1713                                                  columnUpper[iColumn] = CoinMin(columnUpper[iColumn], value);
1714                                             }
1715                                        }
1716                                   } else if (valid == 1) {
1717                                        abort();
1718                                   } else {
1719                                        std::cout << "enter value for " << parameters[iParam].name() <<
1720                                                  std::endl;
1721                                   }
1722                              }
1723                              break;
1724                         case CLP_PARAM_ACTION_REALLY_SCALE:
1725                              if (goodModels[iModel]) {
1726                                   ClpSimplex newModel(models[iModel],
1727                                                       models[iModel].scalingFlag());
1728                                   printf("model really really scaled\n");
1729                                   models[iModel] = newModel;
1730                              }
1731                              break;
1732                         case CLP_PARAM_ACTION_USERCLP:
1733                              // Replace the sample code by whatever you want
1734                              if (goodModels[iModel]) {
1735                                   ClpSimplex * thisModel = &models[iModel];
1736                                   printf("Dummy user code - model has %d rows and %d columns\n",
1737                                          thisModel->numberRows(), thisModel->numberColumns());
1738                              }
1739                              break;
1740                         case CLP_PARAM_ACTION_HELP:
1741                              std::cout << "Coin LP version " << CLP_VERSION
1742                                        << ", build " << __DATE__ << std::endl;
1743                              std::cout << "Non default values:-" << std::endl;
1744                              std::cout << "Perturbation " << models[0].perturbation() << " (default 100)"
1745                                        << std::endl;
1746                              CoinReadPrintit(
1747                                   "Presolve being done with 5 passes\n\
1748Dual steepest edge steep/partial on matrix shape and factorization density\n\
1749Clpnnnn taken out of messages\n\
1750If Factorization frequency default then done on size of matrix\n\n\
1751(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
1752You can switch to interactive mode at any time so\n\
1753clp watson.mps -scaling off -primalsimplex\nis the same as\n\
1754clp watson.mps -\nscaling off\nprimalsimplex"
1755                              );
1756                              break;
1757                         case CLP_PARAM_ACTION_SOLUTION:
1758                              if (goodModels[iModel]) {
1759                                   // get next field
1760                                   field = CoinReadGetString(argc, argv);
1761                                   bool append = false;
1762                                   if (field == "append$") {
1763                                     field = "$";
1764                                     append = true;
1765                                   }
1766                                   if (field == "$") {
1767                                        field = parameters[iParam].stringValue();
1768                                   } else if (field == "EOL") {
1769                                        parameters[iParam].printString();
1770                                        break;
1771                                   } else {
1772                                        parameters[iParam].setStringValue(field);
1773                                   }
1774                                   std::string fileName;
1775                                   FILE *fp = NULL;
1776                                   if (field == "-" || field == "EOL" || field == "stdout") {
1777                                        // stdout
1778                                        fp = stdout;
1779                                        fprintf(fp, "\n");
1780                                   } else if (field == "stderr") {
1781                                        // stderr
1782                                        fp = stderr;
1783                                        fprintf(fp, "\n");
1784                                   } else {
1785                                        if (field[0] == '/' || field[0] == '\\') {
1786                                             fileName = field;
1787                                        } else if (field[0] == '~') {
1788                                             char * environVar = getenv("HOME");
1789                                             if (environVar) {
1790                                                  std::string home(environVar);
1791                                                  field = field.erase(0, 1);
1792                                                  fileName = home + field;
1793                                             } else {
1794                                                  fileName = field;
1795                                             }
1796                                        } else {
1797                                             fileName = directory + field;
1798                                        }
1799                                        if (!append)
1800                                          fp = fopen(fileName.c_str(), "w");
1801                                        else
1802                                          fp = fopen(fileName.c_str(), "a");
1803                                   }
1804                                   if (fp) {
1805                                        // Write solution header (suggested by Luigi Poderico)
1806                                        double objValue = models[iModel].getObjValue() * models[iModel].getObjSense();
1807                                        int iStat = models[iModel].status();
1808                                        if (iStat == 0) {
1809                                             fprintf(fp, "optimal\n" );
1810                                        } else if (iStat == 1) {
1811                                             // infeasible
1812                                             fprintf(fp, "infeasible\n" );
1813                                        } else if (iStat == 2) {
1814                                             // unbounded
1815                                             fprintf(fp, "unbounded\n" );
1816                                        } else if (iStat == 3) {
1817                                             fprintf(fp, "stopped on iterations or time\n" );
1818                                        } else if (iStat == 4) {
1819                                             fprintf(fp, "stopped on difficulties\n" );
1820                                        } else if (iStat == 5) {
1821                                             fprintf(fp, "stopped on ctrl-c\n" );
1822                                        } else {
1823                                             fprintf(fp, "status unknown\n" );
1824                                        }
1825                                        fprintf(fp, "Objective value %15.8g\n", objValue);
1826                                        // make fancy later on
1827                                        int iRow;
1828                                        int numberRows = models[iModel].numberRows();
1829                                        int lengthName = models[iModel].lengthNames(); // 0 if no names
1830                                        int lengthPrint = CoinMax(lengthName, 8);
1831                                        // in general I don't want to pass around massive
1832                                        // amounts of data but seems simpler here
1833                                        std::vector<std::string> rowNames =
1834                                             *(models[iModel].rowNames());
1835                                        std::vector<std::string> columnNames =
1836                                             *(models[iModel].columnNames());
1837
1838                                        double * dualRowSolution = models[iModel].dualRowSolution();
1839                                        double * primalRowSolution =
1840                                             models[iModel].primalRowSolution();
1841                                        double * rowLower = models[iModel].rowLower();
1842                                        double * rowUpper = models[iModel].rowUpper();
1843                                        double primalTolerance = models[iModel].primalTolerance();
1844                                        bool doMask = (printMask != "" && lengthName);
1845                                        int * maskStarts = NULL;
1846                                        int maxMasks = 0;
1847                                        char ** masks = NULL;
1848                                        if (doMask) {
1849                                             int nAst = 0;
1850                                             const char * pMask2 = printMask.c_str();
1851                                             char pMask[100];
1852                                             int iChar;
1853                                             int lengthMask = strlen(pMask2);
1854                                             assert (lengthMask < 100);
1855                                             if (*pMask2 == '"') {
1856                                                  if (pMask2[lengthMask-1] != '"') {
1857                                                       printf("mismatched \" in mask %s\n", pMask2);
1858                                                       break;
1859                                                  } else {
1860                                                       strcpy(pMask, pMask2 + 1);
1861                                                       *strchr(pMask, '"') = '\0';
1862                                                  }
1863                                             } else if (*pMask2 == '\'') {
1864                                                  if (pMask2[lengthMask-1] != '\'') {
1865                                                       printf("mismatched ' in mask %s\n", pMask2);
1866                                                       break;
1867                                                  } else {
1868                                                       strcpy(pMask, pMask2 + 1);
1869                                                       *strchr(pMask, '\'') = '\0';
1870                                                  }
1871                                             } else {
1872                                                  strcpy(pMask, pMask2);
1873                                             }
1874                                             if (lengthMask > lengthName) {
1875                                                  printf("mask %s too long - skipping\n", pMask);
1876                                                  break;
1877                                             }
1878                                             maxMasks = 1;
1879                                             for (iChar = 0; iChar < lengthMask; iChar++) {
1880                                                  if (pMask[iChar] == '*') {
1881                                                       nAst++;
1882                                                       maxMasks *= (lengthName + 1);
1883                                                  }
1884                                             }
1885                                             int nEntries = 1;
1886                                             maskStarts = new int[lengthName+2];
1887                                             masks = new char * [maxMasks];
1888                                             char ** newMasks = new char * [maxMasks];
1889                                             int i;
1890                                             for (i = 0; i < maxMasks; i++) {
1891                                                  masks[i] = new char[lengthName+1];
1892                                                  newMasks[i] = new char[lengthName+1];
1893                                             }
1894                                             strcpy(masks[0], pMask);
1895                                             for (int iAst = 0; iAst < nAst; iAst++) {
1896                                                  int nOldEntries = nEntries;
1897                                                  nEntries = 0;
1898                                                  for (int iEntry = 0; iEntry < nOldEntries; iEntry++) {
1899                                                       char * oldMask = masks[iEntry];
1900                                                       char * ast = strchr(oldMask, '*');
1901                                                       assert (ast);
1902                                                       int length = strlen(oldMask) - 1;
1903                                                       int nBefore = ast - oldMask;
1904                                                       int nAfter = length - nBefore;
1905                                                       // and add null
1906                                                       nAfter++;
1907                                                       for (int i = 0; i <= lengthName - length; i++) {
1908                                                            char * maskOut = newMasks[nEntries];
1909                                                            CoinMemcpyN(oldMask, nBefore, maskOut);
1910                                                            for (int k = 0; k < i; k++)
1911                                                                 maskOut[k+nBefore] = '?';
1912                                                            CoinMemcpyN(ast + 1, nAfter, maskOut + nBefore + i);
1913                                                            nEntries++;
1914                                                            assert (nEntries <= maxMasks);
1915                                                       }
1916                                                  }
1917                                                  char ** temp = masks;
1918                                                  masks = newMasks;
1919                                                  newMasks = temp;
1920                                             }
1921                                             // Now extend and sort
1922                                             int * sort = new int[nEntries];
1923                                             for (i = 0; i < nEntries; i++) {
1924                                                  char * maskThis = masks[i];
1925                                                  int length = strlen(maskThis);
1926                                                  while (maskThis[length-1] == ' ')
1927                                                       length--;
1928                                                  maskThis[length] = '\0';
1929                                                  sort[i] = length;
1930                                             }
1931                                             CoinSort_2(sort, sort + nEntries, masks);
1932                                             int lastLength = -1;
1933                                             for (i = 0; i < nEntries; i++) {
1934                                                  int length = sort[i];
1935                                                  while (length > lastLength)
1936                                                       maskStarts[++lastLength] = i;
1937                                             }
1938                                             maskStarts[++lastLength] = nEntries;
1939                                             delete [] sort;
1940                                             for (i = 0; i < maxMasks; i++)
1941                                                  delete [] newMasks[i];
1942                                             delete [] newMasks;
1943                                        }
1944                                        if (printMode > 5) {
1945                                          int numberColumns = models[iModel].numberColumns();
1946                                          // column length unless rhs ranging
1947                                          int number = numberColumns;
1948                                          switch (printMode) {
1949                                            // bound ranging
1950                                            case 6:
1951                                              fprintf(fp,"Bound ranging");
1952                                              break;
1953                                            // rhs ranging
1954                                            case 7:
1955                                              fprintf(fp,"Rhs ranging");
1956                                              number = numberRows;
1957                                              break;
1958                                            // objective ranging
1959                                            case 8:
1960                                              fprintf(fp,"Objective ranging");
1961                                              break;
1962                                          }
1963                                          if (lengthName)
1964                                            fprintf(fp,",name");
1965                                          fprintf(fp,",increase,variable,decrease,variable\n");
1966                                          int * which = new int [ number];
1967                                          if (printMode != 7) {
1968                                            if (!doMask) {
1969                                              for (int i = 0; i < number;i ++)
1970                                                which[i]=i;
1971                                            } else {
1972                                              int n = 0;
1973                                              for (int i = 0; i < number;i ++) {
1974                                                if (maskMatches(maskStarts,masks,columnNames[i]))
1975                                                  which[n++]=i;
1976                                              }
1977                                              if (n) {
1978                                                number=n;
1979                                              } else {
1980                                                printf("No names match - doing all\n");
1981                                                for (int i = 0; i < number;i ++)
1982                                                  which[i]=i;
1983                                              }
1984                                            }
1985                                          } else {
1986                                            if (!doMask) {
1987                                              for (int i = 0; i < number;i ++)
1988                                                which[i]=i+numberColumns;
1989                                            } else {
1990                                              int n = 0;
1991                                              for (int i = 0; i < number;i ++) {
1992                                                if (maskMatches(maskStarts,masks,rowNames[i]))
1993                                                  which[n++]=i+numberColumns;
1994                                              }
1995                                              if (n) {
1996                                                number=n;
1997                                              } else {
1998                                                printf("No names match - doing all\n");
1999                                                for (int i = 0; i < number;i ++)
2000                                                  which[i]=i+numberColumns;
2001                                              }
2002                                            }
2003                                          }
2004                                          double * valueIncrease = new double [ number];
2005                                          int * sequenceIncrease = new int [ number];
2006                                          double * valueDecrease = new double [ number];
2007                                          int * sequenceDecrease = new int [ number];
2008                                          switch (printMode) {
2009                                            // bound or rhs ranging
2010                                            case 6:
2011                                            case 7:
2012                                              models[iModel].primalRanging(numberRows,
2013                                                                           which, valueIncrease, sequenceIncrease,
2014                                                                           valueDecrease, sequenceDecrease);
2015                                              break;
2016                                            // objective ranging
2017                                            case 8:
2018                                              models[iModel].dualRanging(number,
2019                                                                           which, valueIncrease, sequenceIncrease,
2020                                                                           valueDecrease, sequenceDecrease);
2021                                              break;
2022                                          }
2023                                          for (int i = 0; i < number; i++) {
2024                                            int iWhich = which[i];
2025                                            fprintf(fp, "%d,", (iWhich<numberColumns) ? iWhich : iWhich-numberColumns);
2026                                            if (lengthName) {
2027                                              const char * name = (printMode==7) ? rowNames[iWhich-numberColumns].c_str() : columnNames[iWhich].c_str();
2028                                              fprintf(fp,"%s,",name);
2029                                            }
2030                                            if (valueIncrease[i]<1.0e30) {
2031                                              fprintf(fp, "%.10g,", valueIncrease[i]);
2032                                              int outSequence = sequenceIncrease[i];
2033                                              if (outSequence<numberColumns) {
2034                                                if (lengthName)
2035                                                  fprintf(fp,"%s,",columnNames[outSequence].c_str());
2036                                                else
2037                                                  fprintf(fp,"C%7.7d,",outSequence);
2038                                              } else {
2039                                                outSequence -= numberColumns;
2040                                                if (lengthName)
2041                                                  fprintf(fp,"%s,",rowNames[outSequence].c_str());
2042                                                else
2043                                                  fprintf(fp,"R%7.7d,",outSequence);
2044                                              }
2045                                            } else {
2046                                              fprintf(fp,"1.0e100,,");
2047                                            }
2048                                            if (valueDecrease[i]<1.0e30) {
2049                                              fprintf(fp, "%.10g,", valueDecrease[i]);
2050                                              int outSequence = sequenceDecrease[i];
2051                                              if (outSequence<numberColumns) {
2052                                                if (lengthName)
2053                                                  fprintf(fp,"%s",columnNames[outSequence].c_str());
2054                                                else
2055                                                  fprintf(fp,"C%7.7d",outSequence);
2056                                              } else {
2057                                                outSequence -= numberColumns;
2058                                                if (lengthName)
2059                                                  fprintf(fp,"%s",rowNames[outSequence].c_str());
2060                                                else
2061                                                  fprintf(fp,"R%7.7d",outSequence);
2062                                              }
2063                                            } else {
2064                                              fprintf(fp,"1.0e100,");
2065                                            }
2066                                            fprintf(fp,"\n");
2067                                          }
2068                                          if (fp != stdout)
2069                                            fclose(fp);
2070                                          delete [] which;
2071                                          delete [] valueIncrease;
2072                                          delete [] sequenceIncrease;
2073                                          delete [] valueDecrease;
2074                                          delete [] sequenceDecrease;
2075                                          if (masks) {
2076                                            delete [] maskStarts;
2077                                            for (int i = 0; i < maxMasks; i++)
2078                                              delete [] masks[i];
2079                                            delete [] masks;
2080                                          }
2081                                          break;
2082                                        }
2083                                        if (printMode > 2) {
2084                                             for (iRow = 0; iRow < numberRows; iRow++) {
2085                                                  int type = printMode - 3;
2086                                                  if (primalRowSolution[iRow] > rowUpper[iRow] + primalTolerance ||
2087                                                            primalRowSolution[iRow] < rowLower[iRow] - primalTolerance) {
2088                                                       fprintf(fp, "** ");
2089                                                       type = 2;
2090                                                  } else if (fabs(primalRowSolution[iRow]) > 1.0e-8) {
2091                                                       type = 1;
2092                                                  } else if (numberRows < 50) {
2093                                                       type = 3;
2094                                                  }
2095                                                  if (doMask && !maskMatches(maskStarts, masks, rowNames[iRow]))
2096                                                       type = 0;
2097                                                  if (type) {
2098                                                       fprintf(fp, "%7d ", iRow);
2099                                                       if (lengthName) {
2100                                                            const char * name = rowNames[iRow].c_str();
2101                                                            int n = strlen(name);
2102                                                            int i;
2103                                                            for (i = 0; i < n; i++)
2104                                                                 fprintf(fp, "%c", name[i]);
2105                                                            for (; i < lengthPrint; i++)
2106                                                                 fprintf(fp, " ");
2107                                                       }
2108                                                       fprintf(fp, " %15.8g        %15.8g\n", primalRowSolution[iRow],
2109                                                               dualRowSolution[iRow]);
2110                                                  }
2111                                             }
2112                                        }
2113                                        int iColumn;
2114                                        int numberColumns = models[iModel].numberColumns();
2115                                        double * dualColumnSolution =
2116                                             models[iModel].dualColumnSolution();
2117                                        double * primalColumnSolution =
2118                                             models[iModel].primalColumnSolution();
2119                                        double * columnLower = models[iModel].columnLower();
2120                                        double * columnUpper = models[iModel].columnUpper();
2121                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2122                                             int type = (printMode > 3) ? 1 : 0;
2123                                             if (primalColumnSolution[iColumn] > columnUpper[iColumn] + primalTolerance ||
2124                                                       primalColumnSolution[iColumn] < columnLower[iColumn] - primalTolerance) {
2125                                                  fprintf(fp, "** ");
2126                                                  type = 2;
2127                                             } else if (fabs(primalColumnSolution[iColumn]) > 1.0e-8) {
2128                                                  type = 1;
2129                                             } else if (numberColumns < 50) {
2130                                                  type = 3;
2131                                             }
2132                                             if (doMask && !maskMatches(maskStarts, masks,
2133                                                                        columnNames[iColumn]))
2134                                                  type = 0;
2135                                             if (type) {
2136                                                  fprintf(fp, "%7d ", iColumn);
2137                                                  if (lengthName) {
2138                                                       const char * name = columnNames[iColumn].c_str();
2139                                                       int n = strlen(name);
2140                                                       int i;
2141                                                       for (i = 0; i < n; i++)
2142                                                            fprintf(fp, "%c", name[i]);
2143                                                       for (; i < lengthPrint; i++)
2144                                                            fprintf(fp, " ");
2145                                                  }
2146                                                  fprintf(fp, " %15.8g        %15.8g\n",
2147                                                          primalColumnSolution[iColumn],
2148                                                          dualColumnSolution[iColumn]);
2149                                             }
2150                                        }
2151                                        if (fp != stdout)
2152                                             fclose(fp);
2153                                        if (masks) {
2154                                             delete [] maskStarts;
2155                                             for (int i = 0; i < maxMasks; i++)
2156                                                  delete [] masks[i];
2157                                             delete [] masks;
2158                                        }
2159                                   } else {
2160                                        std::cout << "Unable to open file " << fileName << std::endl;
2161                                   }
2162                              } else {
2163                                   std::cout << "** Current model not valid" << std::endl;
2164
2165                              }
2166
2167                              break;
2168                         case CLP_PARAM_ACTION_SAVESOL:
2169                              if (goodModels[iModel]) {
2170                                   // get next field
2171                                   field = CoinReadGetString(argc, argv);
2172                                   if (field == "$") {
2173                                        field = parameters[iParam].stringValue();
2174                                   } else if (field == "EOL") {
2175                                        parameters[iParam].printString();
2176                                        break;
2177                                   } else {
2178                                        parameters[iParam].setStringValue(field);
2179                                   }
2180                                   std::string fileName;
2181                                   if (field[0] == '/' || field[0] == '\\') {
2182                                        fileName = field;
2183                                   } else if (field[0] == '~') {
2184                                        char * environVar = getenv("HOME");
2185                                        if (environVar) {
2186                                             std::string home(environVar);
2187                                             field = field.erase(0, 1);
2188                                             fileName = home + field;
2189                                        } else {
2190                                             fileName = field;
2191                                        }
2192                                   } else {
2193                                        fileName = directory + field;
2194                                   }
2195                                   saveSolution(models + iModel, fileName);
2196                              } else {
2197                                   std::cout << "** Current model not valid" << std::endl;
2198
2199                              }
2200                              break;
2201                         case CLP_PARAM_ACTION_ENVIRONMENT:
2202                              CbcOrClpEnvironmentIndex = 0;
2203                              break;
2204                         default:
2205                              abort();
2206                         }
2207                    }
2208               } else if (!numberMatches) {
2209                    std::cout << "No match for " << field << " - ? for list of commands"
2210                              << std::endl;
2211               } else if (numberMatches == 1) {
2212                    if (!numberQuery) {
2213                         std::cout << "Short match for " << field << " - completion: ";
2214                         std::cout << parameters[firstMatch].matchName() << std::endl;
2215                    } else if (numberQuery) {
2216                         std::cout << parameters[firstMatch].matchName() << " : ";
2217                         std::cout << parameters[firstMatch].shortHelp() << std::endl;
2218                         if (numberQuery >= 2)
2219                              parameters[firstMatch].printLongHelp();
2220                    }
2221               } else {
2222                    if (!numberQuery)
2223                         std::cout << "Multiple matches for " << field << " - possible completions:"
2224                                   << std::endl;
2225                    else
2226                         std::cout << "Completions of " << field << ":" << std::endl;
2227                    for ( iParam = 0; iParam < numberParameters; iParam++ ) {
2228                         int match = parameters[iParam].matches(field);
2229                         if (match && parameters[iParam].displayThis()) {
2230                              std::cout << parameters[iParam].matchName();
2231                              if (numberQuery >= 2)
2232                                   std::cout << " : " << parameters[iParam].shortHelp();
2233                              std::cout << std::endl;
2234                         }
2235                    }
2236               }
2237          }
2238          delete [] models;
2239          delete [] goodModels;
2240     }
2241     // By now all memory should be freed
2242#ifdef DMALLOC
2243     dmalloc_log_unfreed();
2244     dmalloc_shutdown();
2245#endif
2246     return 0;
2247}
2248static void breakdown(const char * name, int numberLook, const double * region)
2249{
2250     double range[] = {
2251          -COIN_DBL_MAX,
2252          -1.0e15, -1.0e11, -1.0e8, -1.0e5, -1.0e4, -1.0e3, -1.0e2, -1.0e1,
2253          -1.0,
2254          -1.0e-1, -1.0e-2, -1.0e-3, -1.0e-4, -1.0e-5, -1.0e-8, -1.0e-11, -1.0e-15,
2255          0.0,
2256          1.0e-15, 1.0e-11, 1.0e-8, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1,
2257          1.0,
2258          1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15,
2259          COIN_DBL_MAX
2260     };
2261     int nRanges = static_cast<int> (sizeof(range) / sizeof(double));
2262     int * number = new int[nRanges];
2263     memset(number, 0, nRanges * sizeof(int));
2264     int * numberExact = new int[nRanges];
2265     memset(numberExact, 0, nRanges * sizeof(int));
2266     int i;
2267     for ( i = 0; i < numberLook; i++) {
2268          double value = region[i];
2269          for (int j = 0; j < nRanges; j++) {
2270               if (value == range[j]) {
2271                    numberExact[j]++;
2272                    break;
2273               } else if (value < range[j]) {
2274                    number[j]++;
2275                    break;
2276               }
2277          }
2278     }
2279     printf("\n%s has %d entries\n", name, numberLook);
2280     for (i = 0; i < nRanges; i++) {
2281          if (number[i])
2282               printf("%d between %g and %g", number[i], range[i-1], range[i]);
2283          if (numberExact[i]) {
2284               if (number[i])
2285                    printf(", ");
2286               printf("%d exactly at %g", numberExact[i], range[i]);
2287          }
2288          if (number[i] + numberExact[i])
2289               printf("\n");
2290     }
2291     delete [] number;
2292     delete [] numberExact;
2293}
2294void sortOnOther(int * column,
2295                 const CoinBigIndex * rowStart,
2296                 int * order,
2297                 int * other,
2298                 int nRow,
2299                 int nInRow,
2300                 int where)
2301{
2302     if (nRow < 2 || where >= nInRow)
2303          return;
2304     // do initial sort
2305     int kRow;
2306     int iRow;
2307     for ( kRow = 0; kRow < nRow; kRow++) {
2308          iRow = order[kRow];
2309          other[kRow] = column[rowStart[iRow] + where];
2310     }
2311     CoinSort_2(other, other + nRow, order);
2312     int first = 0;
2313     iRow = order[0];
2314     int firstC = column[rowStart[iRow] + where];
2315     kRow = 1;
2316     while (kRow < nRow) {
2317          int lastC = 9999999;;
2318          for (; kRow < nRow + 1; kRow++) {
2319               if (kRow < nRow) {
2320                    iRow = order[kRow];
2321                    lastC = column[rowStart[iRow] + where];
2322               } else {
2323                    lastC = 9999999;
2324               }
2325               if (lastC > firstC)
2326                    break;
2327          }
2328          // sort
2329          sortOnOther(column, rowStart, order + first, other, kRow - first,
2330                      nInRow, where + 1);
2331          firstC = lastC;
2332          first = kRow;
2333     }
2334}
2335static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
2336{
2337     int numberColumns = originalModel->numberColumns();
2338     const char * integerInformation  = originalModel->integerInformation();
2339     const double * columnLower = originalModel->columnLower();
2340     const double * columnUpper = originalModel->columnUpper();
2341     int numberIntegers = 0;
2342     int numberBinary = 0;
2343     int iRow, iColumn;
2344     if (integerInformation) {
2345          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2346               if (integerInformation[iColumn]) {
2347                    if (columnUpper[iColumn] > columnLower[iColumn]) {
2348                         numberIntegers++;
2349                         if (columnUpper[iColumn] == 0.0 && columnLower[iColumn] == 1)
2350                              numberBinary++;
2351                    }
2352               }
2353          }
2354     }
2355     numberColumns = model->numberColumns();
2356     int numberRows = model->numberRows();
2357     columnLower = model->columnLower();
2358     columnUpper = model->columnUpper();
2359     const double * rowLower = model->rowLower();
2360     const double * rowUpper = model->rowUpper();
2361     const double * objective = model->objective();
2362     CoinPackedMatrix * matrix = model->matrix();
2363     CoinBigIndex numberElements = matrix->getNumElements();
2364     const int * columnLength = matrix->getVectorLengths();
2365     //const CoinBigIndex * columnStart = matrix->getVectorStarts();
2366     const double * elementByColumn = matrix->getElements();
2367     int * number = new int[numberRows+1];
2368     memset(number, 0, (numberRows + 1)*sizeof(int));
2369     int numberObjSingletons = 0;
2370     /* cType
2371        0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
2372        8 0/1
2373     */
2374     int cType[9];
2375     std::string cName[] = {"0.0->inf,", "0.0->up,", "lo->inf,", "lo->up,", "free,", "fixed,", "-inf->0.0,",
2376                            "-inf->up,", "0.0->1.0"
2377                           };
2378     int nObjective = 0;
2379     memset(cType, 0, sizeof(cType));
2380     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2381          int length = columnLength[iColumn];
2382          if (length == 1 && objective[iColumn])
2383               numberObjSingletons++;
2384          number[length]++;
2385          if (objective[iColumn])
2386               nObjective++;
2387          if (columnLower[iColumn] > -1.0e20) {
2388               if (columnLower[iColumn] == 0.0) {
2389                    if (columnUpper[iColumn] > 1.0e20)
2390                         cType[0]++;
2391                    else if (columnUpper[iColumn] == 1.0)
2392                         cType[8]++;
2393                    else if (columnUpper[iColumn] == 0.0)
2394                         cType[5]++;
2395                    else
2396                         cType[1]++;
2397               } else {
2398                    if (columnUpper[iColumn] > 1.0e20)
2399                         cType[2]++;
2400                    else if (columnUpper[iColumn] == columnLower[iColumn])
2401                         cType[5]++;
2402                    else
2403                         cType[3]++;
2404               }
2405          } else {
2406               if (columnUpper[iColumn] > 1.0e20)
2407                    cType[4]++;
2408               else if (columnUpper[iColumn] == 0.0)
2409                    cType[6]++;
2410               else
2411                    cType[7]++;
2412          }
2413     }
2414     /* rType
2415        0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
2416        7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
2417     */
2418     int rType[13];
2419     std::string rName[] = {"E 0.0,", "E 1.0,", "E -1.0,", "E other,", "G 0.0,", "G 1.0,", "G other,",
2420                            "L 0.0,", "L 1.0,", "L other,", "Range 0.0->1.0,", "Range other,", "Free"
2421                           };
2422     memset(rType, 0, sizeof(rType));
2423     for (iRow = 0; iRow < numberRows; iRow++) {
2424          if (rowLower[iRow] > -1.0e20) {
2425               if (rowLower[iRow] == 0.0) {
2426                    if (rowUpper[iRow] > 1.0e20)
2427                         rType[4]++;
2428                    else if (rowUpper[iRow] == 1.0)
2429                         rType[10]++;
2430                    else if (rowUpper[iRow] == 0.0)
2431                         rType[0]++;
2432                    else
2433                         rType[11]++;
2434               } else if (rowLower[iRow] == 1.0) {
2435                    if (rowUpper[iRow] > 1.0e20)
2436                         rType[5]++;
2437                    else if (rowUpper[iRow] == rowLower[iRow])
2438                         rType[1]++;
2439                    else
2440                         rType[11]++;
2441               } else if (rowLower[iRow] == -1.0) {
2442                    if (rowUpper[iRow] > 1.0e20)
2443                         rType[6]++;
2444                    else if (rowUpper[iRow] == rowLower[iRow])
2445                         rType[2]++;
2446                    else
2447                         rType[11]++;
2448               } else {
2449                    if (rowUpper[iRow] > 1.0e20)
2450                         rType[6]++;
2451                    else if (rowUpper[iRow] == rowLower[iRow])
2452                         rType[3]++;
2453                    else
2454                         rType[11]++;
2455               }
2456          } else {
2457               if (rowUpper[iRow] > 1.0e20)
2458                    rType[12]++;
2459               else if (rowUpper[iRow] == 0.0)
2460                    rType[7]++;
2461               else if (rowUpper[iRow] == 1.0)
2462                    rType[8]++;
2463               else
2464                    rType[9]++;
2465          }
2466     }
2467     // Basic statistics
2468     printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
2469            numberRows, numberColumns, nObjective, numberElements);
2470     if (number[0] + number[1]) {
2471          printf("There are ");
2472          if (numberObjSingletons)
2473               printf("%d singletons with objective ", numberObjSingletons);
2474          int numberNoObj = number[1] - numberObjSingletons;
2475          if (numberNoObj)
2476               printf("%d singletons with no objective ", numberNoObj);
2477          if (number[0])
2478               printf("** %d columns have no entries", number[0]);
2479          printf("\n");
2480     }
2481     printf("Column breakdown:\n");
2482     int k;
2483     for (k = 0; k < static_cast<int> (sizeof(cType) / sizeof(int)); k++) {
2484          printf("%d of type %s ", cType[k], cName[k].c_str());
2485          if (((k + 1) % 3) == 0)
2486               printf("\n");
2487     }
2488     if ((k % 3) != 0)
2489          printf("\n");
2490     printf("Row breakdown:\n");
2491     for (k = 0; k < static_cast<int> (sizeof(rType) / sizeof(int)); k++) {
2492          printf("%d of type %s ", rType[k], rName[k].c_str());
2493          if (((k + 1) % 3) == 0)
2494               printf("\n");
2495     }
2496     if ((k % 3) != 0)
2497          printf("\n");
2498#define SYM
2499#ifndef SYM
2500     if (model->logLevel() < 2)
2501          return ;
2502#endif
2503     int kMax = model->logLevel() > 3 ? 1000000 : 10;
2504     k = 0;
2505     for (iRow = 1; iRow <= numberRows; iRow++) {
2506          if (number[iRow]) {
2507               k++;
2508               printf("%d columns have %d entries\n", number[iRow], iRow);
2509               if (k == kMax)
2510                    break;
2511          }
2512     }
2513     if (k < numberRows) {
2514          int kk = k;
2515          k = 0;
2516          for (iRow = numberRows; iRow >= 1; iRow--) {
2517               if (number[iRow]) {
2518                    k++;
2519                    if (k == kMax)
2520                         break;
2521               }
2522          }
2523          if (k > kk) {
2524               printf("\n    .........\n\n");
2525               iRow = k;
2526               k = 0;
2527               for (; iRow < numberRows; iRow++) {
2528                    if (number[iRow]) {
2529                         k++;
2530                         printf("%d columns have %d entries\n", number[iRow], iRow);
2531                         if (k == kMax)
2532                              break;
2533                    }
2534               }
2535          }
2536     }
2537     delete [] number;
2538     printf("\n\n");
2539     if (model->logLevel() == 63
2540#ifdef SYM
2541               || true
2542#endif
2543        ) {
2544          // get column copy
2545          CoinPackedMatrix columnCopy = *matrix;
2546          const int * columnLength = columnCopy.getVectorLengths();
2547          number = new int[numberRows+1];
2548          memset(number, 0, (numberRows + 1)*sizeof(int));
2549          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2550               int length = columnLength[iColumn];
2551               number[length]++;
2552          }
2553          k = 0;
2554          for (iRow = 1; iRow <= numberRows; iRow++) {
2555               if (number[iRow]) {
2556                    k++;
2557               }
2558          }
2559          int * row = columnCopy.getMutableIndices();
2560          const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2561          double * element = columnCopy.getMutableElements();
2562          int * order = new int[numberColumns];
2563          int * other = new int[numberColumns];
2564          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2565               int length = columnLength[iColumn];
2566               order[iColumn] = iColumn;
2567               other[iColumn] = length;
2568               CoinBigIndex start = columnStart[iColumn];
2569               CoinSort_2(row + start, row + start + length, element + start);
2570          }
2571          CoinSort_2(other, other + numberColumns, order);
2572          int jColumn = number[0] + number[1];
2573          for (iRow = 2; iRow <= numberRows; iRow++) {
2574               if (number[iRow]) {
2575                    printf("XX %d columns have %d entries\n", number[iRow], iRow);
2576                    int kColumn = jColumn + number[iRow];
2577                    sortOnOther(row, columnStart,
2578                                order + jColumn, other, number[iRow], iRow, 0);
2579                    // Now print etc
2580                    if (iRow < 500000) {
2581                         for (int lColumn = jColumn; lColumn < kColumn; lColumn++) {
2582                              iColumn = order[lColumn];
2583                              CoinBigIndex start = columnStart[iColumn];
2584                              if (model->logLevel() == 63) {
2585                                   printf("column %d %g <= ", iColumn, columnLower[iColumn]);
2586                                   for (CoinBigIndex i = start; i < start + iRow; i++)
2587                                        printf("( %d, %g) ", row[i], element[i]);
2588                                   printf("<= %g\n", columnUpper[iColumn]);
2589                              }
2590                         }
2591                    }
2592                    jColumn = kColumn;
2593               }
2594          }
2595          delete [] order;
2596          delete [] other;
2597          delete [] number;
2598     }
2599     // get row copy
2600     CoinPackedMatrix rowCopy = *matrix;
2601     rowCopy.reverseOrdering();
2602     const int * rowLength = rowCopy.getVectorLengths();
2603     number = new int[numberColumns+1];
2604     memset(number, 0, (numberColumns + 1)*sizeof(int));
2605     for (iRow = 0; iRow < numberRows; iRow++) {
2606          int length = rowLength[iRow];
2607          number[length]++;
2608     }
2609     if (number[0])
2610          printf("** %d rows have no entries\n", number[0]);
2611     k = 0;
2612     for (iColumn = 1; iColumn <= numberColumns; iColumn++) {
2613          if (number[iColumn]) {
2614               k++;
2615               printf("%d rows have %d entries\n", number[iColumn], iColumn);
2616               if (k == kMax)
2617                    break;
2618          }
2619     }
2620     if (k < numberColumns) {
2621          int kk = k;
2622          k = 0;
2623          for (iColumn = numberColumns; iColumn >= 1; iColumn--) {
2624               if (number[iColumn]) {
2625                    k++;
2626                    if (k == kMax)
2627                         break;
2628               }
2629          }
2630          if (k > kk) {
2631               printf("\n    .........\n\n");
2632               iColumn = k;
2633               k = 0;
2634               for (; iColumn < numberColumns; iColumn++) {
2635                    if (number[iColumn]) {
2636                         k++;
2637                         printf("%d rows have %d entries\n", number[iColumn], iColumn);
2638                         if (k == kMax)
2639                              break;
2640                    }
2641               }
2642          }
2643     }
2644     if (model->logLevel() == 63
2645#ifdef SYM
2646               || true
2647#endif
2648        ) {
2649          int * column = rowCopy.getMutableIndices();
2650          const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
2651          double * element = rowCopy.getMutableElements();
2652          int * order = new int[numberRows];
2653          int * other = new int[numberRows];
2654          for (iRow = 0; iRow < numberRows; iRow++) {
2655               int length = rowLength[iRow];
2656               order[iRow] = iRow;
2657               other[iRow] = length;
2658               CoinBigIndex start = rowStart[iRow];
2659               CoinSort_2(column + start, column + start + length, element + start);
2660          }
2661          CoinSort_2(other, other + numberRows, order);
2662          int jRow = number[0] + number[1];
2663          double * weight = new double[numberRows];
2664          double * randomColumn = new double[numberColumns+1];
2665          double * randomRow = new double [numberRows+1];
2666          int * sortRow = new int [numberRows];
2667          int * possibleRow = new int [numberRows];
2668          int * backRow = new int [numberRows];
2669          int * stackRow = new int [numberRows];
2670          int * sortColumn = new int [numberColumns];
2671          int * possibleColumn = new int [numberColumns];
2672          int * backColumn = new int [numberColumns];
2673          int * backColumn2 = new int [numberColumns];
2674          int * mapRow = new int [numberRows];
2675          int * mapColumn = new int [numberColumns];
2676          int * stackColumn = new int [numberColumns];
2677          double randomLower = CoinDrand48();
2678          double randomUpper = CoinDrand48();
2679          double randomInteger = CoinDrand48();
2680          int * startAdd = new int[numberRows+1];
2681          int * columnAdd = new int [2*numberElements];
2682          double * elementAdd = new double[2*numberElements];
2683          int nAddRows = 0;
2684          startAdd[0] = 0;
2685          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2686               randomColumn[iColumn] = CoinDrand48();
2687               backColumn2[iColumn] = -1;
2688          }
2689          for (iColumn = 2; iColumn <= numberColumns; iColumn++) {
2690               if (number[iColumn]) {
2691                    printf("XX %d rows have %d entries\n", number[iColumn], iColumn);
2692                    int kRow = jRow + number[iColumn];
2693                    sortOnOther(column, rowStart,
2694                                order + jRow, other, number[iColumn], iColumn, 0);
2695                    // Now print etc
2696                    if (iColumn < 500000) {
2697                         int nLook = 0;
2698                         for (int lRow = jRow; lRow < kRow; lRow++) {
2699                              iRow = order[lRow];
2700                              CoinBigIndex start = rowStart[iRow];
2701                              if (model->logLevel() == 63) {
2702                                   printf("row %d %g <= ", iRow, rowLower[iRow]);
2703                                   for (CoinBigIndex i = start; i < start + iColumn; i++)
2704                                        printf("( %d, %g) ", column[i], element[i]);
2705                                   printf("<= %g\n", rowUpper[iRow]);
2706                              }
2707                              int first = column[start];
2708                              double sum = 0.0;
2709                              for (CoinBigIndex i = start; i < start + iColumn; i++) {
2710                                   int jColumn = column[i];
2711                                   double value = element[i];
2712                                   jColumn -= first;
2713                                   assert (jColumn >= 0);
2714                                   sum += value * randomColumn[jColumn];
2715                              }
2716                              if (rowLower[iRow] > -1.0e30 && rowLower[iRow])
2717                                   sum += rowLower[iRow] * randomLower;
2718                              else if (!rowLower[iRow])
2719                                   sum += 1.234567e-7 * randomLower;
2720                              if (rowUpper[iRow] < 1.0e30 && rowUpper[iRow])
2721                                   sum += rowUpper[iRow] * randomUpper;
2722                              else if (!rowUpper[iRow])
2723                                   sum += 1.234567e-7 * randomUpper;
2724                              sortRow[nLook] = iRow;
2725                              randomRow[nLook++] = sum;
2726                              // best way is to number unique elements and bounds and use
2727                              if (fabs(sum) > 1.0e4)
2728                                   sum *= 1.0e-6;
2729                              weight[iRow] = sum;
2730                         }
2731                         assert (nLook <= numberRows);
2732                         CoinSort_2(randomRow, randomRow + nLook, sortRow);
2733                         randomRow[nLook] = COIN_DBL_MAX;
2734                         double last = -COIN_DBL_MAX;
2735                         int iLast = -1;
2736                         for (int iLook = 0; iLook < nLook + 1; iLook++) {
2737                              if (randomRow[iLook] > last) {
2738                                   if (iLast >= 0) {
2739                                        int n = iLook - iLast;
2740                                        if (n > 1) {
2741                                             //printf("%d rows possible?\n",n);
2742                                        }
2743                                   }
2744                                   iLast = iLook;
2745                                   last = randomRow[iLook];
2746                              }
2747                         }
2748                    }
2749                    jRow = kRow;
2750               }
2751          }
2752          CoinPackedMatrix columnCopy = *matrix;
2753          const int * columnLength = columnCopy.getVectorLengths();
2754          const int * row = columnCopy.getIndices();
2755          const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2756          const double * elementByColumn = columnCopy.getElements();
2757          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2758               int length = columnLength[iColumn];
2759               CoinBigIndex start = columnStart[iColumn];
2760               double sum = objective[iColumn];
2761               if (columnLower[iColumn] > -1.0e30 && columnLower[iColumn])
2762                    sum += columnLower[iColumn] * randomLower;
2763               else if (!columnLower[iColumn])
2764                    sum += 1.234567e-7 * randomLower;
2765               if (columnUpper[iColumn] < 1.0e30 && columnUpper[iColumn])
2766                    sum += columnUpper[iColumn] * randomUpper;
2767               else if (!columnUpper[iColumn])
2768                    sum += 1.234567e-7 * randomUpper;
2769               if (model->isInteger(iColumn))
2770                    sum += 9.87654321e-6 * randomInteger;
2771               for (CoinBigIndex i = start; i < start + length; i++) {
2772                    int iRow = row[i];
2773                    sum += elementByColumn[i] * weight[iRow];
2774               }
2775               sortColumn[iColumn] = iColumn;
2776               randomColumn[iColumn] = sum;
2777          }
2778          {
2779               CoinSort_2(randomColumn, randomColumn + numberColumns, sortColumn);
2780               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2781                    int i = sortColumn[iColumn];
2782                    backColumn[i] = iColumn;
2783               }
2784               randomColumn[numberColumns] = COIN_DBL_MAX;
2785               double last = -COIN_DBL_MAX;
2786               int iLast = -1;
2787               for (int iLook = 0; iLook < numberColumns + 1; iLook++) {
2788                    if (randomColumn[iLook] > last) {
2789                         if (iLast >= 0) {
2790                              int n = iLook - iLast;
2791                              if (n > 1) {
2792                                   //printf("%d columns possible?\n",n);
2793                              }
2794                              for (int i = iLast; i < iLook; i++) {
2795                                   possibleColumn[sortColumn[i]] = n;
2796                              }
2797                         }
2798                         iLast = iLook;
2799                         last = randomColumn[iLook];
2800                    }
2801               }
2802               for (iRow = 0; iRow < numberRows; iRow++) {
2803                    CoinBigIndex start = rowStart[iRow];
2804                    double sum = 0.0;
2805                    int length = rowLength[iRow];
2806                    for (CoinBigIndex i = start; i < start + length; i++) {
2807                         int jColumn = column[i];
2808                         double value = element[i];
2809                         jColumn = backColumn[jColumn];
2810                         sum += value * randomColumn[jColumn];
2811                         //if (iColumn==23089||iRow==23729)
2812                         //printf("row %d cola %d colb %d value %g rand %g sum %g\n",
2813                         //   iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
2814                    }
2815                    sortRow[iRow] = iRow;
2816                    randomRow[iRow] = weight[iRow];
2817                    randomRow[iRow] = sum;
2818               }
2819               CoinSort_2(randomRow, randomRow + numberRows, sortRow);
2820               for (iRow = 0; iRow < numberRows; iRow++) {
2821                    int i = sortRow[iRow];
2822                    backRow[i] = iRow;
2823               }
2824               randomRow[numberRows] = COIN_DBL_MAX;
2825               last = -COIN_DBL_MAX;
2826               iLast = -1;
2827               // Do backward indices from order
2828               for (iRow = 0; iRow < numberRows; iRow++) {
2829                    other[order[iRow]] = iRow;
2830               }
2831               for (int iLook = 0; iLook < numberRows + 1; iLook++) {
2832                    if (randomRow[iLook] > last) {
2833                         if (iLast >= 0) {
2834                              int n = iLook - iLast;
2835                              if (n > 1) {
2836                                   //printf("%d rows possible?\n",n);
2837                                   // Within group sort as for original "order"
2838                                   for (int i = iLast; i < iLook; i++) {
2839                                        int jRow = sortRow[i];
2840                                        order[i] = other[jRow];
2841                                   }
2842                                   CoinSort_2(order + iLast, order + iLook, sortRow + iLast);
2843                              }
2844                              for (int i = iLast; i < iLook; i++) {
2845                                   possibleRow[sortRow[i]] = n;
2846                              }
2847                         }
2848                         iLast = iLook;
2849                         last = randomRow[iLook];
2850                    }
2851               }
2852               // Temp out
2853               for (int iLook = 0; iLook < numberRows - 1000000; iLook++) {
2854                    iRow = sortRow[iLook];
2855                    CoinBigIndex start = rowStart[iRow];
2856                    int length = rowLength[iRow];
2857                    int numberPossible = possibleRow[iRow];
2858                    for (CoinBigIndex i = start; i < start + length; i++) {
2859                         int jColumn = column[i];
2860                         if (possibleColumn[jColumn] != numberPossible)
2861                              numberPossible = -1;
2862                    }
2863                    int n = numberPossible;
2864                    if (numberPossible > 1) {
2865                         //printf("pppppossible %d\n",numberPossible);
2866                         for (int jLook = iLook + 1; jLook < iLook + numberPossible; jLook++) {
2867                              int jRow = sortRow[jLook];
2868                              CoinBigIndex start2 = rowStart[jRow];
2869                              assert (numberPossible == possibleRow[jRow]);
2870                              assert(length == rowLength[jRow]);
2871                              for (CoinBigIndex i = start2; i < start2 + length; i++) {
2872                                   int jColumn = column[i];
2873                                   if (possibleColumn[jColumn] != numberPossible)
2874                                        numberPossible = -1;
2875                              }
2876                         }
2877                         if (numberPossible < 2) {
2878                              // switch off
2879                              for (int jLook = iLook; jLook < iLook + n; jLook++)
2880                                   possibleRow[sortRow[jLook]] = -1;
2881                         }
2882                         // skip rest
2883                         iLook += n - 1;
2884                    } else {
2885                         possibleRow[iRow] = -1;
2886                    }
2887               }
2888               for (int iLook = 0; iLook < numberRows; iLook++) {
2889                    iRow = sortRow[iLook];
2890                    int numberPossible = possibleRow[iRow];
2891                    // Only if any integers
2892                    int numberIntegers = 0;
2893                    CoinBigIndex start = rowStart[iRow];
2894                    int length = rowLength[iRow];
2895                    for (CoinBigIndex i = start; i < start + length; i++) {
2896                         int jColumn = column[i];
2897                         if (model->isInteger(jColumn))
2898                              numberIntegers++;
2899                    }
2900                    if (numberPossible > 1 && !numberIntegers) {
2901                         //printf("possible %d - but no integers\n",numberPossible);
2902                    }
2903                    if (numberPossible > 1 && (numberIntegers || false)) {
2904                         //
2905                         printf("possible %d - %d integers\n", numberPossible, numberIntegers);
2906                         int lastLook = iLook;
2907                         int nMapRow = -1;
2908                         for (int jLook = iLook + 1; jLook < iLook + numberPossible; jLook++) {
2909                              // stop if too many failures
2910                              if (jLook > iLook + 10 && nMapRow < 0)
2911                                   break;
2912                              // Create identity mapping
2913                              int i;
2914                              for (i = 0; i < numberRows; i++)
2915                                   mapRow[i] = i;
2916                              for (i = 0; i < numberColumns; i++)
2917                                   mapColumn[i] = i;
2918                              int offset = jLook - iLook;
2919                              int nStackC = 0;
2920                              // build up row and column mapping
2921                              int nStackR = 1;
2922                              stackRow[0] = iLook;
2923                              bool good = true;
2924                              while (nStackR) {
2925                                   nStackR--;
2926                                   int look1 = stackRow[nStackR];
2927                                   int look2 = look1 + offset;
2928                                   assert (randomRow[look1] == randomRow[look2]);
2929                                   int row1 = sortRow[look1];
2930                                   int row2 = sortRow[look2];
2931                                   assert (mapRow[row1] == row1);
2932                                   assert (mapRow[row2] == row2);
2933                                   mapRow[row1] = row2;
2934                                   mapRow[row2] = row1;
2935                                   CoinBigIndex start1 = rowStart[row1];
2936                                   CoinBigIndex offset2 = rowStart[row2] - start1;
2937                                   int length = rowLength[row1];
2938                                   assert( length == rowLength[row2]);
2939                                   for (CoinBigIndex i = start1; i < start1 + length; i++) {
2940                                        int jColumn1 = column[i];
2941                                        int jColumn2 = column[i+offset2];
2942                                        if (randomColumn[backColumn[jColumn1]] !=
2943                                                  randomColumn[backColumn[jColumn2]]) {
2944                                             good = false;
2945                                             break;
2946                                        }
2947                                        if (mapColumn[jColumn1] == jColumn1) {
2948                                             // not touched
2949                                             assert (mapColumn[jColumn2] == jColumn2);
2950                                             if (jColumn1 != jColumn2) {
2951                                                  // Put on stack
2952                                                  mapColumn[jColumn1] = jColumn2;
2953                                                  mapColumn[jColumn2] = jColumn1;
2954                                                  stackColumn[nStackC++] = jColumn1;
2955                                             }
2956                                        } else {
2957                                             if (mapColumn[jColumn1] != jColumn2 ||
2958                                                       mapColumn[jColumn2] != jColumn1) {
2959                                                  // bad
2960                                                  good = false;
2961                                                  printf("bad col\n");
2962                                                  break;
2963                                             }
2964                                        }
2965                                   }
2966                                   if (!good)
2967                                        break;
2968                                   while (nStackC) {
2969                                        nStackC--;
2970                                        int iColumn = stackColumn[nStackC];
2971                                        int iColumn2 = mapColumn[iColumn];
2972                                        assert (iColumn != iColumn2);
2973                                        int length = columnLength[iColumn];
2974                                        assert (length == columnLength[iColumn2]);
2975                                        CoinBigIndex start = columnStart[iColumn];
2976                                        CoinBigIndex offset2 = columnStart[iColumn2] - start;
2977                                        for (CoinBigIndex i = start; i < start + length; i++) {
2978                                             int iRow = row[i];
2979                                             int iRow2 = row[i+offset2];
2980                                             if (mapRow[iRow] == iRow) {
2981                                                  // First (but be careful)
2982                                                  if (iRow != iRow2) {
2983                                                       //mapRow[iRow]=iRow2;
2984                                                       //mapRow[iRow2]=iRow;
2985                                                       int iBack = backRow[iRow];
2986                                                       int iBack2 = backRow[iRow2];
2987                                                       if (randomRow[iBack] == randomRow[iBack2] &&
2988                                                                 iBack2 - iBack == offset) {
2989                                                            stackRow[nStackR++] = iBack;
2990                                                       } else {
2991                                                            //printf("randomRow diff - weights %g %g\n",
2992                                                            //     weight[iRow],weight[iRow2]);
2993                                                            // bad
2994                                                            good = false;
2995                                                            break;
2996                                                       }
2997                                                  }
2998                                             } else {
2999                                                  if (mapRow[iRow] != iRow2 ||
3000                                                            mapRow[iRow2] != iRow) {
3001                                                       // bad
3002                                                       good = false;
3003                                                       printf("bad row\n");
3004                                                       break;
3005                                                  }
3006                                             }
3007                                        }
3008                                        if (!good)
3009                                             break;
3010                                   }
3011                              }
3012                              // then check OK
3013                              if (good) {
3014                                   for (iRow = 0; iRow < numberRows; iRow++) {
3015                                        CoinBigIndex start = rowStart[iRow];
3016                                        int length = rowLength[iRow];
3017                                        if (mapRow[iRow] == iRow) {
3018                                             for (CoinBigIndex i = start; i < start + length; i++) {
3019                                                  int jColumn = column[i];
3020                                                  backColumn2[jColumn] = i - start;
3021                                             }
3022                                             for (CoinBigIndex i = start; i < start + length; i++) {
3023                                                  int jColumn = column[i];
3024                                                  if (mapColumn[jColumn] != jColumn) {
3025                                                       int jColumn2 = mapColumn[jColumn];
3026                                                       CoinBigIndex i2 = backColumn2[jColumn2];
3027                                                       if (i2 < 0) {
3028                                                            good = false;
3029                                                       } else if (element[i] != element[i2+start]) {
3030                                                            good = false;
3031                                                       }
3032                                                  }
3033                                             }
3034                                             for (CoinBigIndex i = start; i < start + length; i++) {
3035                                                  int jColumn = column[i];
3036                                                  backColumn2[jColumn] = -1;
3037                                             }
3038                                        } else {
3039                                             int row2 = mapRow[iRow];
3040                                             assert (iRow = mapRow[row2]);
3041                                             if (rowLower[iRow] != rowLower[row2] ||
3042                                                       rowLower[row2] != rowLower[iRow])
3043                                                  good = false;
3044                                             CoinBigIndex offset2 = rowStart[row2] - start;
3045                                             for (CoinBigIndex i = start; i < start + length; i++) {
3046                                                  int jColumn = column[i];
3047                                                  double value = element[i];
3048                                                  int jColumn2 = column[i+offset2];
3049                                                  double value2 = element[i+offset2];
3050                                                  if (value != value2 || mapColumn[jColumn] != jColumn2 ||
3051                                                            mapColumn[jColumn2] != jColumn)
3052                                                       good = false;
3053                                             }
3054                                        }
3055                                   }
3056                                   if (good) {
3057                                        // check rim
3058                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3059                                             if (mapColumn[iColumn] != iColumn) {
3060                                                  int iColumn2 = mapColumn[iColumn];
3061                                                  if (objective[iColumn] != objective[iColumn2])
3062                                                       good = false;
3063                                                  if (columnLower[iColumn] != columnLower[iColumn2])
3064                                                       good = false;
3065                                                  if (columnUpper[iColumn] != columnUpper[iColumn2])
3066                                                       good = false;
3067                                                  if (model->isInteger(iColumn) != model->isInteger(iColumn2))
3068                                                       good = false;
3069                                             }
3070                                        }
3071                                   }
3072                                   if (good) {
3073                                        // temp
3074                                        if (nMapRow < 0) {
3075                                             //const double * solution = model->primalColumnSolution();
3076                                             // find mapped
3077                                             int nMapColumn = 0;
3078                                             for (int i = 0; i < numberColumns; i++) {
3079                                                  if (mapColumn[i] > i)
3080                                                       nMapColumn++;
3081                                             }
3082                                             nMapRow = 0;
3083                                             int kRow = -1;
3084                                             for (int i = 0; i < numberRows; i++) {
3085                                                  if (mapRow[i] > i) {
3086                                                       nMapRow++;
3087                                                       kRow = i;
3088                                                  }
3089                                             }
3090                                             printf("%d columns, %d rows\n", nMapColumn, nMapRow);
3091                                             if (nMapRow == 1) {
3092                                                  CoinBigIndex start = rowStart[kRow];
3093                                                  int length = rowLength[kRow];
3094                                                  printf("%g <= ", rowLower[kRow]);
3095                                                  for (CoinBigIndex i = start; i < start + length; i++) {
3096                                                       int jColumn = column[i];
3097                                                       if (mapColumn[jColumn] != jColumn)
3098                                                            printf("* ");
3099                                                       printf("%d,%g ", jColumn, element[i]);
3100                                                  }
3101                                                  printf("<= %g\n", rowUpper[kRow]);
3102                                             }
3103                                        }
3104                                        // temp
3105                                        int row1 = sortRow[lastLook];
3106                                        int row2 = sortRow[jLook];
3107                                        lastLook = jLook;
3108                                        CoinBigIndex start1 = rowStart[row1];
3109                                        CoinBigIndex offset2 = rowStart[row2] - start1;
3110                                        int length = rowLength[row1];
3111                                        assert( length == rowLength[row2]);
3112                                        CoinBigIndex put = startAdd[nAddRows];
3113                                        double multiplier = length < 11 ? 2.0 : 1.125;
3114                                        double value = 1.0;
3115                                        for (CoinBigIndex i = start1; i < start1 + length; i++) {
3116                                             int jColumn1 = column[i];
3117                                             int jColumn2 = column[i+offset2];
3118                                             columnAdd[put] = jColumn1;
3119                                             elementAdd[put++] = value;
3120                                             columnAdd[put] = jColumn2;
3121                                             elementAdd[put++] = -value;
3122                                             value *= multiplier;
3123                                        }
3124                                        nAddRows++;
3125                                        startAdd[nAddRows] = put;
3126                                   } else {
3127                                        printf("ouch - did not check out as good\n");
3128                                   }
3129                              }
3130                         }
3131                         // skip rest
3132                         iLook += numberPossible - 1;
3133                    }
3134               }
3135          }
3136          if (nAddRows) {
3137               double * lower = new double [nAddRows];
3138               double * upper = new double[nAddRows];
3139               int i;
3140               //const double * solution = model->primalColumnSolution();
3141               for (i = 0; i < nAddRows; i++) {
3142                    lower[i] = 0.0;
3143                    upper[i] = COIN_DBL_MAX;
3144               }
3145               printf("Adding %d rows with %d elements\n", nAddRows,
3146                      startAdd[nAddRows]);
3147               //ClpSimplex newModel(*model);
3148               //newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
3149               //newModel.writeMps("modified.mps");
3150               delete [] lower;
3151               delete [] upper;
3152          }
3153          delete [] startAdd;
3154          delete [] columnAdd;
3155          delete [] elementAdd;
3156          delete [] order;
3157          delete [] other;
3158          delete [] randomColumn;
3159          delete [] weight;
3160          delete [] randomRow;
3161          delete [] sortRow;
3162          delete [] backRow;
3163          delete [] possibleRow;
3164          delete [] sortColumn;
3165          delete [] backColumn;
3166          delete [] backColumn2;
3167          delete [] possibleColumn;
3168          delete [] mapRow;
3169          delete [] mapColumn;
3170          delete [] stackRow;
3171          delete [] stackColumn;
3172     }
3173     delete [] number;
3174     // Now do breakdown of ranges
3175     breakdown("Elements", numberElements, elementByColumn);
3176     breakdown("RowLower", numberRows, rowLower);
3177     breakdown("RowUpper", numberRows, rowUpper);
3178     breakdown("ColumnLower", numberColumns, columnLower);
3179     breakdown("ColumnUpper", numberColumns, columnUpper);
3180     breakdown("Objective", numberColumns, objective);
3181}
3182static bool maskMatches(const int * starts, char ** masks,
3183                        std::string & check)
3184{
3185     // back to char as I am old fashioned
3186     const char * checkC = check.c_str();
3187     int length = strlen(checkC);
3188     while (checkC[length-1] == ' ')
3189          length--;
3190     for (int i = starts[length]; i < starts[length+1]; i++) {
3191          char * thisMask = masks[i];
3192          int k;
3193          for ( k = 0; k < length; k++) {
3194               if (thisMask[k] != '?' && thisMask[k] != checkC[k])
3195                    break;
3196          }
3197          if (k == length)
3198               return true;
3199     }
3200     return false;
3201}
3202static void clean(char * temp)
3203{
3204     char * put = temp;
3205     while (*put >= ' ')
3206          put++;
3207     *put = '\0';
3208}
3209static void generateCode(const char * fileName, int type)
3210{
3211     FILE * fp = fopen(fileName, "r");
3212     assert (fp);
3213     int numberLines = 0;
3214#define MAXLINES 500
3215#define MAXONELINE 200
3216     char line[MAXLINES][MAXONELINE];
3217     while (fgets(line[numberLines], MAXONELINE, fp)) {
3218          assert (numberLines < MAXLINES);
3219          clean(line[numberLines]);
3220          numberLines++;
3221     }
3222     fclose(fp);
3223     // add in actual solve
3224     strcpy(line[numberLines], "5  clpModel->initialSolve(clpSolve);");
3225     numberLines++;
3226     fp = fopen(fileName, "w");
3227     assert (fp);
3228     char apo = '"';
3229     char backslash = '\\';
3230
3231     fprintf(fp, "#include %cClpSimplex.hpp%c\n", apo, apo);
3232     fprintf(fp, "#include %cClpSolve.hpp%c\n", apo, apo);
3233
3234     fprintf(fp, "\nint main (int argc, const char *argv[])\n{\n");
3235     fprintf(fp, "  ClpSimplex  model;\n");
3236     fprintf(fp, "  int status=1;\n");
3237     fprintf(fp, "  if (argc<2)\n");
3238     fprintf(fp, "    fprintf(stderr,%cPlease give file name%cn%c);\n",
3239             apo, backslash, apo);
3240     fprintf(fp, "  else\n");
3241     fprintf(fp, "    status=model.readMps(argv[1],true);\n");
3242     fprintf(fp, "  if (status) {\n");
3243     fprintf(fp, "    fprintf(stderr,%cBad readMps %%s%cn%c,argv[1]);\n",
3244             apo, backslash, apo);
3245     fprintf(fp, "    exit(1);\n");
3246     fprintf(fp, "  }\n\n");
3247     fprintf(fp, "  // Now do requested saves and modifications\n");
3248     fprintf(fp, "  ClpSimplex * clpModel = & model;\n");
3249     int wanted[9];
3250     memset(wanted, 0, sizeof(wanted));
3251     wanted[0] = wanted[3] = wanted[5] = wanted[8] = 1;
3252     if (type > 0)
3253          wanted[1] = wanted[6] = 1;
3254     if (type > 1)
3255          wanted[2] = wanted[4] = wanted[7] = 1;
3256     std::string header[9] = { "", "Save values", "Redundant save of default values", "Set changed values",
3257                               "Redundant set default values", "Solve", "Restore values", "Redundant restore values", "Add to model"
3258                             };
3259     for (int iType = 0; iType < 9; iType++) {
3260          if (!wanted[iType])
3261               continue;
3262          int n = 0;
3263          int iLine;
3264          for (iLine = 0; iLine < numberLines; iLine++) {
3265               if (line[iLine][0] == '0' + iType) {
3266                    if (!n)
3267                         fprintf(fp, "\n  // %s\n\n", header[iType].c_str());
3268                    n++;
3269                    fprintf(fp, "%s\n", line[iLine] + 1);
3270               }
3271          }
3272     }
3273     fprintf(fp, "\n  // Now you would use solution etc etc\n\n");
3274     fprintf(fp, "  return 0;\n}\n");
3275     fclose(fp);
3276     printf("C++ file written to %s\n", fileName);
3277}
3278/*
3279  Version 1.00.00 October 13 2004.
3280  1.00.01 October 18.  Added basis handling helped/prodded by Thorsten Koch.
3281  Also modifications to make faster with sbb (I hope I haven't broken anything).
3282  1.00.02 March 21 2005.  Redid ClpNonLinearCost to save memory also redid
3283  createRim to try and improve cache characteristics.
3284  1.00.03 April 8 2005.  Added Volume algorithm as crash and made code more
3285  robust on testing.  Also added "either" and "tune" algorithm.
3286  1.01.01 April 12 2005.  Decided to go to different numbering.  Backups will
3287  be last 2 digits while middle 2 are for improvements.  Still take a long
3288  time to get to 2.00.01
3289  1.01.02 May 4 2005.  Will be putting in many changes - so saving stable version
3290  1.02.01 May 6 2005.  Lots of changes to try and make faster and more stable in
3291  branch and cut.
3292  1.02.02 May 19 2005.  Stuff for strong branching and some improvements to simplex
3293  1.03.01 May 24 2006.  Lots done but I can't remember what!
3294  1.03.03 June 13 2006.  For clean up after dual perturbation
3295  1.04.01 June 26 2007.  Lots of changes but I got lazy
3296  1.05.00 June 27 2007.  This is trunk so when gets to stable will be 1.5
3297  1.11.00 November 5 2009 (Guy Fawkes) - OSL factorization and better ordering
3298 */
Note: See TracBrowser for help on using the repository browser.