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

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