source: trunk/Clp/src/ClpSolver.cpp @ 2371

Last change on this file since 2371 was 2371, checked in by forrest, 5 months ago

more flexibility for sos - also allow more string parameters

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