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

Last change on this file since 2383 was 2383, checked in by forrest, 4 months ago

allow compressed .lp files

File size: 224.7 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                                            (length > 6 && !strncmp(c_name + length - 6, ".lp.gz", 6)) ||
1454                                            (length > 7 && !strncmp(c_name + length - 7, ".lp.bz2", 7)))
1455                                            gmpl = -1; // .lp
1456                                   }
1457                                   bool absolutePath;
1458                                   if (dirsep == '/') {
1459                                        // non Windows (or cygwin)
1460                                        absolutePath = (field[0] == '/');
1461                                   } else {
1462                                        //Windows (non cycgwin)
1463                                        absolutePath = (field[0] == '\\');
1464                                        // but allow for :
1465                                        if (strchr(field.c_str(), ':'))
1466                                             absolutePath = true;
1467                                   }
1468                                   if (absolutePath) {
1469                                        fileName = field;
1470                                        size_t length = field.size();
1471                                        size_t percent = field.find('%');
1472                                        if (percent < length && percent > 0) {
1473                                             gmpl = 1;
1474                                             fileName = field.substr(0, percent);
1475                                             gmplData = field.substr(percent + 1);
1476                                             if (percent < length - 1)
1477                                                  gmpl = 2; // two files
1478                                             printf("GMPL model file %s and data file %s\n",
1479                                                    fileName.c_str(), gmplData.c_str());
1480                                        }
1481                                   } else if (field[0] == '~') {
1482                                        char * environVar = getenv("HOME");
1483                                        if (environVar) {
1484                                             std::string home(environVar);
1485                                             field = field.erase(0, 1);
1486                                             fileName = home + field;
1487                                        } else {
1488                                             fileName = field;
1489                                        }
1490                                   } else {
1491                                        fileName = directory + field;
1492                                        // See if gmpl (model & data) - or even lp file
1493                                        size_t length = field.size();
1494                                        size_t percent = field.find('%');
1495                                        if (percent < length && percent > 0) {
1496                                             gmpl = 1;
1497                                             fileName = directory + field.substr(0, percent);
1498                                             gmplData = directory + field.substr(percent + 1);
1499                                             if (percent < length - 1)
1500                                                  gmpl = 2; // two files
1501                                             printf("GMPL model file %s and data file %s\n",
1502                                                    fileName.c_str(), gmplData.c_str());
1503                                        }
1504                                   }
1505                                   std::string name = fileName;
1506                                   if (fileCoinReadable(name)) {
1507                                        // can open - lets go for it
1508                                        canOpen = true;
1509                                        if (gmpl == 2) {
1510                                             FILE *fp;
1511                                             fp = fopen(gmplData.c_str(), "r");
1512                                             if (fp) {
1513                                                  fclose(fp);
1514                                             } else {
1515                                                  canOpen = false;
1516                                                  std::cout << "Unable to open file " << gmplData << std::endl;
1517                                             }
1518                                        }
1519                                   } else {
1520                                        std::cout << "Unable to open file " << fileName << std::endl;
1521                                   }
1522                              }
1523                              if (canOpen) {
1524                                   int status;
1525#ifdef CLP_USEFUL_PRINTOUT
1526                                   mpsFile=fileName;
1527#endif
1528                                   if (!gmpl)
1529                                        status = models[iModel].readMps(fileName.c_str(),
1530                                                                        keepImportNames != 0,
1531                                                                        allowImportErrors != 0);
1532                                   else if (gmpl > 0)
1533                                        status = models[iModel].readGMPL(fileName.c_str(),
1534                                                                         (gmpl == 2) ? gmplData.c_str() : NULL,
1535                                                                         keepImportNames != 0);
1536                                   else
1537#ifdef KILL_ZERO_READLP
1538                                     status = models[iModel].readLp(fileName.c_str(), models[iModel].getSmallElementValue());
1539#else
1540                                     status = models[iModel].readLp(fileName.c_str(), 1.0e-12);
1541#endif
1542                                   if (!status || (status > 0 && allowImportErrors)) {
1543                                        goodModels[iModel] = true;
1544                                        // sets to all slack (not necessary?)
1545                                        thisModel->createStatus();
1546                                        time2 = CoinCpuTime();
1547                                        totalTime += time2 - time1;
1548                                        time1 = time2;
1549                                        // Go to canned file if just input file
1550                                        if (CbcOrClpRead_mode == 2 && argc == 2) {
1551                                             // only if ends .mps
1552                                             char * find = const_cast<char *>(strstr(fileName.c_str(), ".mps"));
1553                                             if (find && find[4] == '\0') {
1554                                                  find[1] = 'p';
1555                                                  find[2] = 'a';
1556                                                  find[3] = 'r';
1557                                                  FILE *fp = fopen(fileName.c_str(), "r");
1558                                                  if (fp) {
1559                                                       CbcOrClpReadCommand = fp; // Read from that file
1560                                                       CbcOrClpRead_mode = -1;
1561                                                  }
1562                                             }
1563                                        }
1564                                   } else {
1565                                        // errors
1566                                        std::cout << "There were " << status <<
1567                                                  " errors on input" << std::endl;
1568                                   }
1569                              }
1570                         }
1571                         break;
1572                         case CLP_PARAM_ACTION_EXPORT:
1573                              if (goodModels[iModel]) {
1574                                   double objScale =
1575                                        parameters[whichParam(CLP_PARAM_DBL_OBJSCALE2, numberParameters, parameters)].doubleValue();
1576                                   if (objScale != 1.0) {
1577                                        int iColumn;
1578                                        int numberColumns = models[iModel].numberColumns();
1579                                        double * dualColumnSolution =
1580                                             models[iModel].dualColumnSolution();
1581                                        ClpObjective * obj = models[iModel].objectiveAsObject();
1582                                        assert(dynamic_cast<ClpLinearObjective *> (obj));
1583                                        double offset;
1584                                        double * objective = obj->gradient(NULL, NULL, offset, true);
1585                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1586                                             dualColumnSolution[iColumn] *= objScale;
1587                                             objective[iColumn] *= objScale;;
1588                                        }
1589                                        int iRow;
1590                                        int numberRows = models[iModel].numberRows();
1591                                        double * dualRowSolution =
1592                                             models[iModel].dualRowSolution();
1593                                        for (iRow = 0; iRow < numberRows; iRow++)
1594                                             dualRowSolution[iRow] *= objScale;
1595                                        models[iModel].setObjectiveOffset(objScale * models[iModel].objectiveOffset());
1596                                   }
1597                                   // get next field
1598                                   field = CoinReadGetString(argc, argv);
1599                                   if (field == "$") {
1600                                        field = parameters[iParam].stringValue();
1601                                   } else if (field == "EOL") {
1602                                        parameters[iParam].printString();
1603                                        break;
1604                                   } else {
1605                                        parameters[iParam].setStringValue(field);
1606                                   }
1607                                   std::string fileName;
1608                                   bool canOpen = false;
1609                                   if (field[0] == '/' || field[0] == '\\') {
1610                                        fileName = field;
1611                                   } else if (field[0] == '~') {
1612                                        char * environVar = getenv("HOME");
1613                                        if (environVar) {
1614                                             std::string home(environVar);
1615                                             field = field.erase(0, 1);
1616                                             fileName = home + field;
1617                                        } else {
1618                                             fileName = field;
1619                                        }
1620                                   } else {
1621                                        fileName = directory + field;
1622                                   }
1623                                   FILE *fp = fopen(fileName.c_str(), "w");
1624                                   if (fp) {
1625                                        // can open - lets go for it
1626                                        fclose(fp);
1627                                        canOpen = true;
1628                                   } else {
1629                                        std::cout << "Unable to open file " << fileName << std::endl;
1630                                   }
1631                                   if (canOpen) {
1632                                        // If presolve on then save presolved
1633                                        bool deleteModel2 = false;
1634                                        ClpSimplex * model2 = models + iModel;
1635                                        if (dualize && dualize < 3) {
1636                                             model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel();
1637                                             printf("Dual of model has %d rows and %d columns\n",
1638                                                    model2->numberRows(), model2->numberColumns());
1639                                             model2->setOptimizationDirection(1.0);
1640                                             preSolve = 0; // as picks up from model
1641                                        }
1642                                        if (preSolve) {
1643                                             ClpPresolve pinfo;
1644                                             int presolveOptions2 = presolveOptions&~0x40000000;
1645                                             if ((presolveOptions2 & 0xffff) != 0)
1646                                                  pinfo.setPresolveActions(presolveOptions2);
1647                                             pinfo.setSubstitution(substitution);
1648                                             if ((printOptions & 1) != 0)
1649                                                  pinfo.statistics();
1650                                             double presolveTolerance =
1651                                                  parameters[whichParam(CLP_PARAM_DBL_PRESOLVETOLERANCE, numberParameters, parameters)].doubleValue();
1652                                             model2 =
1653                                                  pinfo.presolvedModel(models[iModel], presolveTolerance,
1654                                                                       true, preSolve, false, false);
1655                                             if (model2) {
1656                                                  printf("Saving presolved model on %s\n",
1657                                                         fileName.c_str());
1658                                                  deleteModel2 = true;
1659                                             } else {
1660                                                  printf("Presolved model looks infeasible - saving original on %s\n",
1661                                                         fileName.c_str());
1662                                                  deleteModel2 = false;
1663                                                  model2 = models + iModel;
1664
1665                                             }
1666                                        } else {
1667                                             printf("Saving model on %s\n",
1668                                                    fileName.c_str());
1669                                        }
1670#if 0
1671                                        // Convert names
1672                                        int iRow;
1673                                        int numberRows = model2->numberRows();
1674                                        int iColumn;
1675                                        int numberColumns = model2->numberColumns();
1676
1677                                        char ** rowNames = NULL;
1678                                        char ** columnNames = NULL;
1679                                        if (model2->lengthNames()) {
1680                                             rowNames = new char * [numberRows];
1681                                             for (iRow = 0; iRow < numberRows; iRow++) {
1682                                                  rowNames[iRow] =
1683                                                       CoinStrdup(model2->rowName(iRow).c_str());
1684#ifdef STRIPBLANKS
1685                                                  char * xx = rowNames[iRow];
1686                                                  int i;
1687                                                  int length = strlen(xx);
1688                                                  int n = 0;
1689                                                  for (i = 0; i < length; i++) {
1690                                                       if (xx[i] != ' ')
1691                                                            xx[n++] = xx[i];
1692                                                  }
1693                                                  xx[n] = '\0';
1694#endif
1695                                             }
1696
1697                                             columnNames = new char * [numberColumns];
1698                                             for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1699                                                  columnNames[iColumn] =
1700                                                       CoinStrdup(model2->columnName(iColumn).c_str());
1701#ifdef STRIPBLANKS
1702                                                  char * xx = columnNames[iColumn];
1703                                                  int i;
1704                                                  int length = strlen(xx);
1705                                                  int n = 0;
1706                                                  for (i = 0; i < length; i++) {
1707                                                       if (xx[i] != ' ')
1708                                                            xx[n++] = xx[i];
1709                                                  }
1710                                                  xx[n] = '\0';
1711#endif
1712                                             }
1713                                        }
1714                                        CoinMpsIO writer;
1715                                        writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
1716                                                          model2->getColLower(), model2->getColUpper(),
1717                                                          model2->getObjCoefficients(),
1718                                                          (const char*) 0 /*integrality*/,
1719                                                          model2->getRowLower(), model2->getRowUpper(),
1720                                                          columnNames, rowNames);
1721                                        // Pass in array saying if each variable integer
1722                                        writer.copyInIntegerInformation(model2->integerInformation());
1723                                        writer.setObjectiveOffset(model2->objectiveOffset());
1724                                        writer.writeMps(fileName.c_str(), 0, 1, 1);
1725                                        if (rowNames) {
1726                                             for (iRow = 0; iRow < numberRows; iRow++) {
1727                                                  free(rowNames[iRow]);
1728                                             }
1729                                             delete [] rowNames;
1730                                             for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1731                                                  free(columnNames[iColumn]);
1732                                             }
1733                                             delete [] columnNames;
1734                                        }
1735#else
1736                                        model2->writeMps(fileName.c_str(), (outputFormat - 1) / 2, 1 + ((outputFormat - 1) & 1));
1737#endif
1738                                        if (deleteModel2)
1739                                             delete model2;
1740                                        time2 = CoinCpuTime();
1741                                        totalTime += time2 - time1;
1742                                        time1 = time2;
1743                                   }
1744                              } else {
1745                                   std::cout << "** Current model not valid" << std::endl;
1746                              }
1747                              break;
1748                         case CLP_PARAM_ACTION_BASISIN:
1749                              if (goodModels[iModel]) {
1750                                   // get next field
1751                                   field = CoinReadGetString(argc, argv);
1752                                   if (field == "$") {
1753                                        field = parameters[iParam].stringValue();
1754                                   } else if (field == "EOL") {
1755                                        parameters[iParam].printString();
1756                                        break;
1757                                   } else {
1758                                        parameters[iParam].setStringValue(field);
1759                                   }
1760                                   std::string fileName;
1761                                   bool canOpen = false;
1762                                   if (field == "-") {
1763                                        // stdin
1764                                        canOpen = true;
1765                                        fileName = "-";
1766                                   } else {
1767                                        if (field[0] == '/' || field[0] == '\\') {
1768                                             fileName = field;
1769                                        } else if (field[0] == '~') {
1770                                             char * environVar = getenv("HOME");
1771                                             if (environVar) {
1772                                                  std::string home(environVar);
1773                                                  field = field.erase(0, 1);
1774                                                  fileName = home + field;
1775                                             } else {
1776                                                  fileName = field;
1777                                             }
1778                                        } else {
1779                                             fileName = directory + field;
1780                                        }
1781                                        FILE *fp = fopen(fileName.c_str(), "r");
1782                                        if (fp) {
1783                                             // can open - lets go for it
1784                                             fclose(fp);
1785                                             canOpen = true;
1786                                        } else {
1787                                             std::cout << "Unable to open file " << fileName << std::endl;
1788                                        }
1789                                   }
1790                                   if (canOpen) {
1791                                        int values = thisModel->readBasis(fileName.c_str());
1792                                        if (values == 0)
1793                                             basisHasValues = -1;
1794                                        else
1795                                             basisHasValues = 1;
1796                                   }
1797                              } else {
1798                                   std::cout << "** Current model not valid" << std::endl;
1799                              }
1800                              break;
1801                         case CLP_PARAM_ACTION_PRINTMASK:
1802                              // get next field
1803                         {
1804                              std::string name = CoinReadGetString(argc, argv);
1805                              if (name != "EOL") {
1806                                   parameters[iParam].setStringValue(name);
1807                                   printMask = name;
1808                              } else {
1809                                   parameters[iParam].printString();
1810                              }
1811                         }
1812                         break;
1813                         case CLP_PARAM_ACTION_BASISOUT:
1814                              if (goodModels[iModel]) {
1815                                   // get next field
1816                                   field = CoinReadGetString(argc, argv);
1817                                   if (field == "$") {
1818                                        field = parameters[iParam].stringValue();
1819                                   } else if (field == "EOL") {
1820                                        parameters[iParam].printString();
1821                                        break;
1822                                   } else {
1823                                        parameters[iParam].setStringValue(field);
1824                                   }
1825                                   std::string fileName;
1826                                   bool canOpen = false;
1827                                   if (field[0] == '/' || field[0] == '\\') {
1828                                        fileName = field;
1829                                   } else if (field[0] == '~') {
1830                                        char * environVar = getenv("HOME");
1831                                        if (environVar) {
1832                                             std::string home(environVar);
1833                                             field = field.erase(0, 1);
1834                                             fileName = home + field;
1835                                        } else {
1836                                             fileName = field;
1837                                        }
1838                                   } else {
1839                                        fileName = directory + field;
1840                                   }
1841                                   FILE *fp = fopen(fileName.c_str(), "w");
1842                                   if (fp) {
1843                                        // can open - lets go for it
1844                                        fclose(fp);
1845                                        canOpen = true;
1846                                   } else {
1847                                        std::cout << "Unable to open file " << fileName << std::endl;
1848                                   }
1849                                   if (canOpen) {
1850                                        ClpSimplex * model2 = models + iModel;
1851                                        model2->writeBasis(fileName.c_str(), outputFormat > 1, outputFormat - 2);
1852                                        time2 = CoinCpuTime();
1853                                        totalTime += time2 - time1;
1854                                        time1 = time2;
1855                                   }
1856                              } else {
1857                                   std::cout << "** Current model not valid" << std::endl;
1858                              }
1859                              break;
1860                         case CLP_PARAM_ACTION_PARAMETRICS:
1861                              if (goodModels[iModel]) {
1862                                   // get next field
1863                                   field = CoinReadGetString(argc, argv);
1864                                   if (field == "$") {
1865                                        field = parameters[iParam].stringValue();
1866                                   } else if (field == "EOL") {
1867                                        parameters[iParam].printString();
1868                                        break;
1869                                   } else {
1870                                        parameters[iParam].setStringValue(field);
1871                                   }
1872                                   std::string fileName;
1873                                   //bool canOpen = false;
1874                                   if (field[0] == '/' || field[0] == '\\') {
1875                                        fileName = field;
1876                                   } else if (field[0] == '~') {
1877                                        char * environVar = getenv("HOME");
1878                                        if (environVar) {
1879                                             std::string home(environVar);
1880                                             field = field.erase(0, 1);
1881                                             fileName = home + field;
1882                                        } else {
1883                                             fileName = field;
1884                                        }
1885                                   } else {
1886                                        fileName = directory + field;
1887                                   }
1888                                   ClpSimplex * model2 = models + iModel;
1889                                   static_cast<ClpSimplexOther *> (model2)->parametrics(fileName.c_str());
1890                                   time2 = CoinCpuTime();
1891                                   totalTime += time2 - time1;
1892                                   time1 = time2;
1893                              } else {
1894                                   std::cout << "** Current model not valid" << std::endl;
1895                              }
1896                              break;
1897                         case CLP_PARAM_ACTION_SAVE: {
1898                              // get next field
1899                              field = CoinReadGetString(argc, argv);
1900                              if (field == "$") {
1901                                   field = parameters[iParam].stringValue();
1902                              } else if (field == "EOL") {
1903                                   parameters[iParam].printString();
1904                                   break;
1905                              } else {
1906                                   parameters[iParam].setStringValue(field);
1907                              }
1908                              std::string fileName;
1909                              bool canOpen = false;
1910                              if (field[0] == '/' || field[0] == '\\') {
1911                                   fileName = field;
1912                              } else if (field[0] == '~') {
1913                                   char * environVar = getenv("HOME");
1914                                   if (environVar) {
1915                                        std::string home(environVar);
1916                                        field = field.erase(0, 1);
1917                                        fileName = home + field;
1918                                   } else {
1919                                        fileName = field;
1920                                   }
1921                              } else {
1922                                   fileName = directory + field;
1923                              }
1924                              FILE *fp = fopen(fileName.c_str(), "wb");
1925                              if (fp) {
1926                                   // can open - lets go for it
1927                                   fclose(fp);
1928                                   canOpen = true;
1929                              } else {
1930                                   std::cout << "Unable to open file " << fileName << std::endl;
1931                              }
1932                              if (canOpen) {
1933                                   int status;
1934                                   // If presolve on then save presolved
1935                                   bool deleteModel2 = false;
1936                                   ClpSimplex * model2 = models + iModel;
1937                                   if (preSolve) {
1938                                        ClpPresolve pinfo;
1939                                        double presolveTolerance =
1940                                             parameters[whichParam(CLP_PARAM_DBL_PRESOLVETOLERANCE, numberParameters, parameters)].doubleValue();
1941                                        model2 =
1942                                             pinfo.presolvedModel(models[iModel], presolveTolerance,
1943                                                                  false, preSolve);
1944                                        if (model2) {
1945                                             printf("Saving presolved model on %s\n",
1946                                                    fileName.c_str());
1947                                             deleteModel2 = true;
1948                                        } else {
1949                                             printf("Presolved model looks infeasible - saving original on %s\n",
1950                                                    fileName.c_str());
1951                                             deleteModel2 = false;
1952                                             model2 = models + iModel;
1953
1954                                        }
1955                                   } else {
1956                                        printf("Saving model on %s\n",
1957                                               fileName.c_str());
1958                                   }
1959                                   status = model2->saveModel(fileName.c_str());
1960                                   if (deleteModel2)
1961                                        delete model2;
1962                                   if (!status) {
1963                                        goodModels[iModel] = true;
1964                                        time2 = CoinCpuTime();
1965                                        totalTime += time2 - time1;
1966                                        time1 = time2;
1967                                   } else {
1968                                        // errors
1969                                        std::cout << "There were errors on output" << std::endl;
1970                                   }
1971                              }
1972                         }
1973                         break;
1974                         case CLP_PARAM_ACTION_RESTORE: {
1975                              // get next field
1976                              field = CoinReadGetString(argc, argv);
1977                              if (field == "$") {
1978                                   field = parameters[iParam].stringValue();
1979                              } else if (field == "EOL") {
1980                                   parameters[iParam].printString();
1981                                   break;
1982                              } else {
1983                                   parameters[iParam].setStringValue(field);
1984                              }
1985                              std::string fileName;
1986                              bool canOpen = false;
1987                              if (field[0] == '/' || field[0] == '\\') {
1988                                   fileName = field;
1989                              } else if (field[0] == '~') {
1990                                   char * environVar = getenv("HOME");
1991                                   if (environVar) {
1992                                        std::string home(environVar);
1993                                        field = field.erase(0, 1);
1994                                        fileName = home + field;
1995                                   } else {
1996                                        fileName = field;
1997                                   }
1998                              } else {
1999                                   fileName = directory + field;
2000                              }
2001                              FILE *fp = fopen(fileName.c_str(), "rb");
2002                              if (fp) {
2003                                   // can open - lets go for it
2004                                   fclose(fp);
2005                                   canOpen = true;
2006                              } else {
2007                                   std::cout << "Unable to open file " << fileName << std::endl;
2008                              }
2009                              if (canOpen) {
2010                                   int status = models[iModel].restoreModel(fileName.c_str());
2011                                   if (!status) {
2012                                        goodModels[iModel] = true;
2013                                        time2 = CoinCpuTime();
2014                                        totalTime += time2 - time1;
2015                                        time1 = time2;
2016                                   } else {
2017                                        // errors
2018                                        std::cout << "There were errors on input" << std::endl;
2019                                   }
2020                              }
2021                         }
2022                         break;
2023                         case CLP_PARAM_ACTION_MAXIMIZE:
2024                              models[iModel].setOptimizationDirection(-1);
2025#ifdef ABC_INHERIT
2026                              thisModel->setOptimizationDirection(-1);
2027#endif
2028                              break;
2029                         case CLP_PARAM_ACTION_MINIMIZE:
2030                              models[iModel].setOptimizationDirection(1);
2031#ifdef ABC_INHERIT
2032                              thisModel->setOptimizationDirection(1);
2033#endif
2034                              break;
2035                         case CLP_PARAM_ACTION_ALLSLACK:
2036                              thisModel->allSlackBasis(true);
2037#ifdef ABC_INHERIT
2038                              models[iModel].allSlackBasis();
2039#endif
2040                              break;
2041                         case CLP_PARAM_ACTION_REVERSE:
2042                              if (goodModels[iModel]) {
2043                                   int iColumn;
2044                                   int numberColumns = models[iModel].numberColumns();
2045                                   double * dualColumnSolution =
2046                                        models[iModel].dualColumnSolution();
2047                                   ClpObjective * obj = models[iModel].objectiveAsObject();
2048                                   assert(dynamic_cast<ClpLinearObjective *> (obj));
2049                                   double offset;
2050                                   double * objective = obj->gradient(NULL, NULL, offset, true);
2051                                   for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2052                                        dualColumnSolution[iColumn] = -dualColumnSolution[iColumn];
2053                                        objective[iColumn] = -objective[iColumn];
2054                                   }
2055                                   int iRow;
2056                                   int numberRows = models[iModel].numberRows();
2057                                   double * dualRowSolution =
2058                                        models[iModel].dualRowSolution();
2059                                   for (iRow = 0; iRow < numberRows; iRow++) {
2060                                        dualRowSolution[iRow] = -dualRowSolution[iRow];
2061                                   }
2062                                   models[iModel].setObjectiveOffset(-models[iModel].objectiveOffset());
2063                              }
2064                              break;
2065                         case CLP_PARAM_ACTION_DIRECTORY: {
2066                              std::string name = CoinReadGetString(argc, argv);
2067                              if (name != "EOL") {
2068                                   size_t length = name.length();
2069                                   if (length > 0 && name[length-1] == dirsep) {
2070                                        directory = name;
2071                                   } else {
2072                                        directory = name + dirsep;
2073                                   }
2074                                   parameters[iParam].setStringValue(directory);
2075                              } else {
2076                                   parameters[iParam].printString();
2077                              }
2078                         }
2079                         break;
2080                         case CLP_PARAM_ACTION_DIRSAMPLE: {
2081                              std::string name = CoinReadGetString(argc, argv);
2082                              if (name != "EOL") {
2083                                   size_t length = name.length();
2084                                   if (length > 0 && name[length-1] == dirsep) {
2085                                        dirSample = name;
2086                                   } else {
2087                                        dirSample = name + dirsep;
2088                                   }
2089                                   parameters[iParam].setStringValue(dirSample);
2090                              } else {
2091                                   parameters[iParam].printString();
2092                              }
2093                         }
2094                         break;
2095                         case CLP_PARAM_ACTION_DIRNETLIB: {
2096                              std::string name = CoinReadGetString(argc, argv);
2097                              if (name != "EOL") {
2098                                   size_t length = name.length();
2099                                   if (length > 0 && name[length-1] == dirsep) {
2100                                        dirNetlib = name;
2101                                   } else {
2102                                        dirNetlib = name + dirsep;
2103                                   }
2104                                   parameters[iParam].setStringValue(dirNetlib);
2105                              } else {
2106                                   parameters[iParam].printString();
2107                              }
2108                         }
2109                         break;
2110                         case CBC_PARAM_ACTION_DIRMIPLIB: {
2111                              std::string name = CoinReadGetString(argc, argv);
2112                              if (name != "EOL") {
2113                                   size_t length = name.length();
2114                                   if (length > 0 && name[length-1] == dirsep) {
2115                                        dirMiplib = name;
2116                                   } else {
2117                                        dirMiplib = name + dirsep;
2118                                   }
2119                                   parameters[iParam].setStringValue(dirMiplib);
2120                              } else {
2121                                   parameters[iParam].printString();
2122                              }
2123                         }
2124                         break;
2125                         case CLP_PARAM_ACTION_STDIN:
2126                              CbcOrClpRead_mode = -1;
2127                              break;
2128                         case CLP_PARAM_ACTION_NETLIB_DUAL:
2129                         case CLP_PARAM_ACTION_NETLIB_EITHER:
2130                         case CLP_PARAM_ACTION_NETLIB_BARRIER:
2131                         case CLP_PARAM_ACTION_NETLIB_PRIMAL:
2132                         case CLP_PARAM_ACTION_NETLIB_TUNE: {
2133                              // create fields for unitTest
2134                              const char * fields[4];
2135                              int nFields = 4;
2136                              fields[0] = "fake main from unitTest";
2137                              std::string mpsfield = "-dirSample=";
2138                              mpsfield += dirSample.c_str();
2139                              fields[1] = mpsfield.c_str();
2140                              std::string netfield = "-dirNetlib=";
2141                              netfield += dirNetlib.c_str();
2142                              fields[2] = netfield.c_str();
2143                              fields[3] = "-netlib";
2144                              int algorithm;
2145                              if (type == CLP_PARAM_ACTION_NETLIB_DUAL) {
2146                                   std::cerr << "Doing netlib with dual algorithm" << std::endl;
2147                                   algorithm = 0;
2148                                   models[iModel].setMoreSpecialOptions(models[iModel].moreSpecialOptions()|32768);
2149                              } else if (type == CLP_PARAM_ACTION_NETLIB_BARRIER) {
2150                                   std::cerr << "Doing netlib with barrier algorithm" << std::endl;
2151                                   algorithm = 2;
2152                              } else if (type == CLP_PARAM_ACTION_NETLIB_EITHER) {
2153                                   std::cerr << "Doing netlib with dual or primal algorithm" << std::endl;
2154                                   algorithm = 3;
2155                              } else if (type == CLP_PARAM_ACTION_NETLIB_TUNE) {
2156                                   std::cerr << "Doing netlib with best algorithm!" << std::endl;
2157                                   algorithm = 5;
2158                                   // uncomment next to get active tuning
2159                                   // algorithm=6;
2160                              } else {
2161                                   std::cerr << "Doing netlib with primal algorithm" << std::endl;
2162                                   algorithm = 1;
2163                              }
2164                              //int specialOptions = models[iModel].specialOptions();
2165                              models[iModel].setSpecialOptions(0);
2166                              ClpSolve solveOptions;
2167                              ClpSolve::PresolveType presolveType;
2168                              if (preSolve)
2169                                   presolveType = ClpSolve::presolveOn;
2170                              else
2171                                   presolveType = ClpSolve::presolveOff;
2172                              solveOptions.setPresolveType(presolveType, 5);
2173                              if (doSprint >= 0 || doIdiot >= 0) {
2174                                   if (doSprint > 0) {
2175                                        // sprint overrides idiot
2176                                        solveOptions.setSpecialOption(1, 3, doSprint); // sprint
2177                                   } else if (doIdiot > 0) {
2178                                        solveOptions.setSpecialOption(1, 2, doIdiot); // idiot
2179                                   } else {
2180                                        if (doIdiot == 0) {
2181                                             if (doSprint == 0)
2182                                                  solveOptions.setSpecialOption(1, 4); // all slack
2183                                             else
2184                                                  solveOptions.setSpecialOption(1, 9); // all slack or sprint
2185                                        } else {
2186                                             if (doSprint == 0)
2187                                                  solveOptions.setSpecialOption(1, 8); // all slack or idiot
2188                                             else
2189                                                  solveOptions.setSpecialOption(1, 7); // initiative
2190                                        }
2191                                   }
2192                              }
2193#if FACTORIZATION_STATISTICS
2194                              {
2195                                extern int loSizeX;
2196                                extern int hiSizeX;
2197                                for (int i=0;i<argc;i++) {
2198                                  if (!strcmp(argv[i],"-losize")) {
2199                                    int size=atoi(argv[i+1]);
2200                                    if (size>0)
2201                                      loSizeX=size;
2202                                  }
2203                                  if (!strcmp(argv[i],"-hisize")) {
2204                                    int size=atoi(argv[i+1]);
2205                                    if (size>loSizeX)
2206                                      hiSizeX=size;
2207                                  }
2208                                }
2209                                if (loSizeX!=-1||hiSizeX!=1000000)
2210                                  printf("Solving problems %d<= and <%d\n",loSizeX,hiSizeX);
2211                              }
2212#endif
2213                              // for moment then back to models[iModel]
2214#ifndef ABC_INHERIT
2215                              int specialOptions = models[iModel].specialOptions();
2216                              mainTest(nFields, fields, algorithm, *thisModel,
2217                                       solveOptions, specialOptions, doVector != 0);
2218#else
2219                              //if (!processTune) {
2220                              //specialOptions=0;
2221                              //models->setSpecialOptions(models->specialOptions()&~65536);
2222                              // }
2223                              mainTest(nFields, fields, algorithm, *models,
2224                                       solveOptions, 0, doVector != 0);
2225#endif
2226                         }
2227                         break;
2228                         case CLP_PARAM_ACTION_UNITTEST: {
2229                              // create fields for unitTest
2230                              const char * fields[2];
2231                              int nFields = 2;
2232                              fields[0] = "fake main from unitTest";
2233                              std::string dirfield = "-dirSample=";
2234                              dirfield += dirSample.c_str();
2235                              fields[1] = dirfield.c_str();
2236                              int specialOptions = models[iModel].specialOptions();
2237                              models[iModel].setSpecialOptions(0);
2238                              int algorithm = -1;
2239                              if (models[iModel].numberRows())
2240                                   algorithm = 7;
2241                              ClpSolve solveOptions;
2242                              ClpSolve::PresolveType presolveType;
2243                              if (preSolve)
2244                                   presolveType = ClpSolve::presolveOn;
2245                              else
2246                                   presolveType = ClpSolve::presolveOff;
2247                              solveOptions.setPresolveType(presolveType, 5);
2248#ifndef ABC_INHERIT
2249                              mainTest(nFields, fields, algorithm, *thisModel,
2250                                       solveOptions, specialOptions, doVector != 0);
2251#else
2252                              mainTest(nFields, fields, algorithm, *models,
2253                                       solveOptions, specialOptions, doVector != 0);
2254#endif
2255                         }
2256                         break;
2257                         case CLP_PARAM_ACTION_FAKEBOUND:
2258                              if (goodModels[iModel]) {
2259                                   // get bound
2260                                   double value = CoinReadGetDoubleField(argc, argv, &valid);
2261                                   if (!valid) {
2262                                        std::cout << "Setting " << parameters[iParam].name() <<
2263                                                  " to DEBUG " << value << std::endl;
2264                                        int iRow;
2265                                        int numberRows = models[iModel].numberRows();
2266                                        double * rowLower = models[iModel].rowLower();
2267                                        double * rowUpper = models[iModel].rowUpper();
2268                                        for (iRow = 0; iRow < numberRows; iRow++) {
2269                                             // leave free ones for now
2270                                             if (rowLower[iRow] > -1.0e20 || rowUpper[iRow] < 1.0e20) {
2271                                                  rowLower[iRow] = CoinMax(rowLower[iRow], -value);
2272                                                  rowUpper[iRow] = CoinMin(rowUpper[iRow], value);
2273                                             }
2274                                        }
2275                                        int iColumn;
2276                                        int numberColumns = models[iModel].numberColumns();
2277                                        double * columnLower = models[iModel].columnLower();
2278                                        double * columnUpper = models[iModel].columnUpper();
2279                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2280                                             // leave free ones for now
2281                                             if (columnLower[iColumn] > -1.0e20 ||
2282                                                       columnUpper[iColumn] < 1.0e20) {
2283                                                  columnLower[iColumn] = CoinMax(columnLower[iColumn], -value);
2284                                                  columnUpper[iColumn] = CoinMin(columnUpper[iColumn], value);
2285                                             }
2286                                        }
2287                                   } else if (valid == 1) {
2288                                        abort();
2289                                   } else {
2290                                        std::cout << "enter value for " << parameters[iParam].name() <<
2291                                                  std::endl;
2292                                   }
2293                              }
2294                              break;
2295                         case CLP_PARAM_ACTION_REALLY_SCALE:
2296                              if (goodModels[iModel]) {
2297                                   ClpSimplex newModel(models[iModel],
2298                                                       models[iModel].scalingFlag());
2299                                   printf("model really really scaled\n");
2300                                   models[iModel] = newModel;
2301                              }
2302                              break;
2303                         case CLP_PARAM_ACTION_USERCLP:
2304                              // Replace the sample code by whatever you want
2305                              if (goodModels[iModel]) {
2306                                   ClpSimplex * thisModel = &models[iModel];
2307                                   printf("Dummy user code - model has %d rows and %d columns\n",
2308                                          thisModel->numberRows(), thisModel->numberColumns());
2309                              }
2310                              break;
2311                         case CLP_PARAM_ACTION_HELP:
2312                              std::cout << "Coin LP version " << CLP_VERSION
2313                                        << ", build " << __DATE__ << std::endl;
2314                              std::cout << "Non default values:-" << std::endl;
2315                              std::cout << "Perturbation " << models[0].perturbation() << " (default 100)"
2316                                        << std::endl;
2317                              CoinReadPrintit(
2318                                   "Presolve being done with 5 passes\n\
2319Dual steepest edge steep/partial on matrix shape and factorization density\n\
2320Clpnnnn taken out of messages\n\
2321If Factorization frequency default then done on size of matrix\n\n\
2322(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
2323You can switch to interactive mode at any time so\n\
2324clp watson.mps -scaling off -primalsimplex\nis the same as\n\
2325clp watson.mps -\nscaling off\nprimalsimplex"
2326                              );
2327                              break;
2328                         case CLP_PARAM_ACTION_SOLUTION:
2329                         case CLP_PARAM_ACTION_GMPL_SOLUTION:
2330                              if (goodModels[iModel]) {
2331                                   // get next field
2332                                   field = CoinReadGetString(argc, argv);
2333                                   bool append = false;
2334                                   if (field == "append$") {
2335                                     field = "$";
2336                                     append = true;
2337                                   }
2338                                   if (field == "$") {
2339                                        field = parameters[iParam].stringValue();
2340                                   } else if (field == "EOL") {
2341                                        parameters[iParam].printString();
2342                                        break;
2343                                   } else {
2344                                        parameters[iParam].setStringValue(field);
2345                                   }
2346                                   std::string fileName;
2347                                   FILE *fp = NULL;
2348                                   if (field == "-" || field == "EOL" || field == "stdout") {
2349                                        // stdout
2350                                        fp = stdout;
2351                                        fprintf(fp, "\n");
2352                                   } else if (field == "stderr") {
2353                                        // stderr
2354                                        fp = stderr;
2355                                        fprintf(fp, "\n");
2356                                   } else {
2357                                        if (field[0] == '/' || field[0] == '\\') {
2358                                             fileName = field;
2359                                        } else if (field[0] == '~') {
2360                                             char * environVar = getenv("HOME");
2361                                             if (environVar) {
2362                                                  std::string home(environVar);
2363                                                  field = field.erase(0, 1);
2364                                                  fileName = home + field;
2365                                             } else {
2366                                                  fileName = field;
2367                                             }
2368                                        } else {
2369                                             fileName = directory + field;
2370                                        }
2371                                        if (!append)
2372                                          fp = fopen(fileName.c_str(), "w");
2373                                        else
2374                                          fp = fopen(fileName.c_str(), "a");
2375                                   }
2376                                   if (fp) {
2377                                     // See if Glpk
2378                                     if (type == CLP_PARAM_ACTION_GMPL_SOLUTION) {
2379                                       int numberRows = models[iModel].getNumRows();
2380                                       int numberColumns = models[iModel].getNumCols();
2381                                       int numberGlpkRows=numberRows+1;
2382#ifdef COIN_HAS_GLPK
2383                                       if (cbc_glp_prob) {
2384                                         // from gmpl
2385                                         numberGlpkRows=glp_get_num_rows(cbc_glp_prob);
2386                                         if (numberGlpkRows!=numberRows)
2387                                           printf("Mismatch - cbc %d rows, glpk %d\n",
2388                                                  numberRows,numberGlpkRows);
2389                                       }
2390#endif
2391                                       fprintf(fp,"%d %d\n",numberGlpkRows,
2392                                               numberColumns);
2393                                       int iStat = models[iModel].status();
2394                                       int iStat2 = GLP_UNDEF;
2395                                       if (iStat == 0) {
2396                                         // optimal
2397                                         iStat2 = GLP_FEAS;
2398                                       } else if (iStat == 1) {
2399                                         // infeasible
2400                                         iStat2 = GLP_NOFEAS;
2401                                       } else if (iStat == 2) {
2402                                         // unbounded
2403                                         // leave as 1
2404                                       } else if (iStat >= 3 && iStat <= 5) {
2405                                         iStat2 = GLP_FEAS;
2406                                       }
2407                                       double objValue = models[iModel].getObjValue() 
2408                                         * models[iModel].getObjSense();
2409                                       fprintf(fp,"%d 2 %g\n",iStat2,objValue);
2410                                       if (numberGlpkRows > numberRows) {
2411                                         // objective as row
2412                                         fprintf(fp,"4 %g 1.0\n",objValue);
2413                                       }
2414                                       int lookup[6]=
2415                                         {4,1,3,2,4,5};
2416                                       const double * primalRowSolution =
2417                                         models[iModel].primalRowSolution();
2418                                       const double * dualRowSolution =
2419                                         models[iModel].dualRowSolution();
2420                                       for (int i=0;i<numberRows;i++) {
2421                                         fprintf(fp,"%d %g %g\n",lookup[models[iModel].getRowStatus(i)],
2422                                                 primalRowSolution[i],dualRowSolution[i]);
2423                                       }
2424                                       const double * primalColumnSolution =
2425                                         models[iModel].primalColumnSolution();
2426                                       const double * dualColumnSolution =
2427                                         models[iModel].dualColumnSolution();
2428                                       for (int i=0;i<numberColumns;i++) {
2429                                         fprintf(fp,"%d %g %g\n",lookup[models[iModel].getColumnStatus(i)],
2430                                                 primalColumnSolution[i],dualColumnSolution[i]);
2431                                       }
2432                                       fclose(fp);
2433#ifdef COIN_HAS_GLPK
2434                                       if (cbc_glp_prob) {
2435                                         glp_read_sol(cbc_glp_prob,fileName.c_str());
2436                                         glp_mpl_postsolve(cbc_glp_tran,
2437                                                           cbc_glp_prob,
2438                                                           GLP_SOL);
2439                                         // free up as much as possible
2440                                         glp_free(cbc_glp_prob);
2441                                         glp_mpl_free_wksp(cbc_glp_tran);
2442                                         cbc_glp_prob = NULL;
2443                                         cbc_glp_tran = NULL;
2444                                         //gmp_free_mem();
2445                                         /* check that no memory blocks are still allocated */
2446                                         glp_free_env();
2447                                       }
2448#endif
2449                                       break;
2450                                     }
2451                                        // Write solution header (suggested by Luigi Poderico)
2452                                        double objValue = models[iModel].getObjValue() * models[iModel].getObjSense();
2453                                        int iStat = models[iModel].status();
2454                                        if (iStat == 0) {
2455                                             fprintf(fp, "optimal\n" );
2456                                        } else if (iStat == 1) {
2457                                             // infeasible
2458                                             fprintf(fp, "infeasible\n" );
2459                                        } else if (iStat == 2) {
2460                                             // unbounded
2461                                             fprintf(fp, "unbounded\n" );
2462                                        } else if (iStat == 3) {
2463                                             fprintf(fp, "stopped on iterations or time\n" );
2464                                        } else if (iStat == 4) {
2465                                             fprintf(fp, "stopped on difficulties\n" );
2466                                        } else if (iStat == 5) {
2467                                             fprintf(fp, "stopped on ctrl-c\n" );
2468                                        } else {
2469                                             fprintf(fp, "status unknown\n" );
2470                                        }
2471                                        char printFormat[50];
2472                                        sprintf(printFormat,"Objective value %s\n",
2473                                                CLP_QUOTE(CLP_OUTPUT_FORMAT));
2474                                        fprintf(fp, printFormat, objValue);
2475                                        if (printMode==9) {
2476                                          // just statistics
2477                                          int numberRows = models[iModel].numberRows();
2478                                          double * dualRowSolution = models[iModel].dualRowSolution();
2479                                          double * primalRowSolution =
2480                                            models[iModel].primalRowSolution();
2481                                          double * rowLower = models[iModel].rowLower();
2482                                          double * rowUpper = models[iModel].rowUpper();
2483                                          double highestPrimal;
2484                                          double lowestPrimal;
2485                                          double highestDual;
2486                                          double lowestDual;
2487                                          double largestAway;
2488                                          int numberAtLower;
2489                                          int numberAtUpper;
2490                                          int numberBetween;
2491                                          highestPrimal=-COIN_DBL_MAX;
2492                                          lowestPrimal=COIN_DBL_MAX;
2493                                          highestDual=-COIN_DBL_MAX;
2494                                          lowestDual=COIN_DBL_MAX;
2495                                          largestAway=0.0;;
2496                                          numberAtLower=0;
2497                                          numberAtUpper=0;
2498                                          numberBetween=0;
2499                                          for (int iRow=0;iRow<numberRows;iRow++) {
2500                                            double primal=primalRowSolution[iRow];
2501                                            double lower=rowLower[iRow];
2502                                            double upper=rowUpper[iRow];
2503                                            double dual=dualRowSolution[iRow];
2504                                            highestPrimal=CoinMax(highestPrimal,primal);
2505                                            lowestPrimal=CoinMin(lowestPrimal,primal);
2506                                            highestDual=CoinMax(highestDual,dual);
2507                                            lowestDual=CoinMin(lowestDual,dual);
2508                                            if (primal<lower+1.0e-6) {
2509                                              numberAtLower++;
2510                                            } else if (primal>upper-1.0e-6) {
2511                                              numberAtUpper++;
2512                                            } else {
2513                                              numberBetween++;
2514                                              largestAway=CoinMax(largestAway,
2515                                                                  CoinMin(primal-lower,upper-primal));
2516                                            }
2517                                          }
2518                                          printf("For rows %d at lower, %d between, %d at upper - lowest %g, highest %g most away %g - highest dual %g lowest %g\n",
2519                                                 numberAtLower,numberBetween,numberAtUpper,
2520                                                 lowestPrimal,highestPrimal,largestAway,
2521                                                 lowestDual,highestDual);
2522                                          int numberColumns = models[iModel].numberColumns();
2523                                          double * dualColumnSolution = models[iModel].dualColumnSolution();
2524                                          double * primalColumnSolution =
2525                                            models[iModel].primalColumnSolution();
2526                                          double * columnLower = models[iModel].columnLower();
2527                                          double * columnUpper = models[iModel].columnUpper();
2528                                          highestPrimal=-COIN_DBL_MAX;
2529                                          lowestPrimal=COIN_DBL_MAX;
2530                                          highestDual=-COIN_DBL_MAX;
2531                                          lowestDual=COIN_DBL_MAX;
2532                                          largestAway=0.0;;
2533                                          numberAtLower=0;
2534                                          numberAtUpper=0;
2535                                          numberBetween=0;
2536                                          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2537                                            double primal=primalColumnSolution[iColumn];
2538                                            double lower=columnLower[iColumn];
2539                                            double upper=columnUpper[iColumn];
2540                                            double dual=dualColumnSolution[iColumn];
2541                                            highestPrimal=CoinMax(highestPrimal,primal);
2542                                            lowestPrimal=CoinMin(lowestPrimal,primal);
2543                                            highestDual=CoinMax(highestDual,dual);
2544                                            lowestDual=CoinMin(lowestDual,dual);
2545                                            if (primal<lower+1.0e-6) {
2546                                              numberAtLower++;
2547                                            } else if (primal>upper-1.0e-6) {
2548                                              numberAtUpper++;
2549                                            } else {
2550                                              numberBetween++;
2551                                              largestAway=CoinMax(largestAway,
2552                                                                  CoinMin(primal-lower,upper-primal));
2553                                            }
2554                                          }
2555                                          printf("For columns %d at lower, %d between, %d at upper - lowest %g, highest %g most away %g - highest dual %g lowest %g\n",
2556                                                 numberAtLower,numberBetween,numberAtUpper,
2557                                                 lowestPrimal,highestPrimal,largestAway,
2558                                                 lowestDual,highestDual);
2559                                          break;
2560                                        }
2561                                        // make fancy later on
2562                                        int iRow;
2563                                        int numberRows = models[iModel].numberRows();
2564                                        int lengthName = models[iModel].lengthNames(); // 0 if no names
2565                                        int lengthPrint = CoinMax(lengthName, 8);
2566                                        // in general I don't want to pass around massive
2567                                        // amounts of data but seems simpler here
2568                                        std::vector<std::string> rowNames =
2569                                             *(models[iModel].rowNames());
2570                                        std::vector<std::string> columnNames =
2571                                             *(models[iModel].columnNames());
2572
2573                                        double * dualRowSolution = models[iModel].dualRowSolution();
2574                                        double * primalRowSolution =
2575                                             models[iModel].primalRowSolution();
2576                                        double * rowLower = models[iModel].rowLower();
2577                                        double * rowUpper = models[iModel].rowUpper();
2578                                        double primalTolerance = models[iModel].primalTolerance();
2579                                        bool doMask = (printMask != "" && lengthName);
2580                                        int * maskStarts = NULL;
2581                                        int maxMasks = 0;
2582                                        char ** masks = NULL;
2583                                        if (doMask) {
2584                                             int nAst = 0;
2585                                             const char * pMask2 = printMask.c_str();
2586                                             char pMask[100];
2587                                             size_t lengthMask = strlen(pMask2);
2588                                             assert (lengthMask < 100);
2589                                             if (*pMask2 == '"') {
2590                                                  if (pMask2[lengthMask-1] != '"') {
2591                                                       printf("mismatched \" in mask %s\n", pMask2);
2592                                                       break;
2593                                                  } else {
2594                                                       strcpy(pMask, pMask2 + 1);
2595                                                       *strchr(pMask, '"') = '\0';
2596                                                  }
2597                                             } else if (*pMask2 == '\'') {
2598                                                  if (pMask2[lengthMask-1] != '\'') {
2599                                                       printf("mismatched ' in mask %s\n", pMask2);
2600                                                       break;
2601                                                  } else {
2602                                                       strcpy(pMask, pMask2 + 1);
2603                                                       *strchr(pMask, '\'') = '\0';
2604                                                  }
2605                                             } else {
2606                                                  strcpy(pMask, pMask2);
2607                                             }
2608                                             if (lengthMask > static_cast<size_t>(lengthName)) {
2609                                                  printf("mask %s too long - skipping\n", pMask);
2610                                                  break;
2611                                             }
2612                                             maxMasks = 1;
2613                                             for (size_t iChar = 0; iChar < lengthMask; iChar++) {
2614                                                  if (pMask[iChar] == '*') {
2615                                                       nAst++;
2616                                                       maxMasks *= (lengthName + 1);
2617                                                  }
2618                                             }
2619                                             int nEntries = 1;
2620                                             maskStarts = new int[lengthName+2];
2621                                             masks = new char * [maxMasks];
2622                                             char ** newMasks = new char * [maxMasks];
2623                                             int i;
2624                                             for (i = 0; i < maxMasks; i++) {
2625                                                  masks[i] = new char[lengthName+1];
2626                                                  newMasks[i] = new char[lengthName+1];
2627                                             }
2628                                             strcpy(masks[0], pMask);
2629                                             for (int iAst = 0; iAst < nAst; iAst++) {
2630                                                  int nOldEntries = nEntries;
2631                                                  nEntries = 0;
2632                                                  for (int iEntry = 0; iEntry < nOldEntries; iEntry++) {
2633                                                       char * oldMask = masks[iEntry];
2634                                                       char * ast = strchr(oldMask, '*');
2635                                                       assert (ast);
2636                                                       size_t length = strlen(oldMask) - 1;
2637                                                       size_t nBefore = ast - oldMask;
2638                                                       size_t nAfter = length - nBefore;
2639                                                       // and add null
2640                                                       nAfter++;
2641                                                       for (size_t i = 0; i <= lengthName - length; i++) {
2642                                                            char * maskOut = newMasks[nEntries];
2643                                                            CoinMemcpyN(oldMask, static_cast<int>(nBefore), maskOut);
2644                                                            for (size_t k = 0; k < i; k++)
2645                                                                 maskOut[k+nBefore] = '?';
2646                                                            CoinMemcpyN(ast + 1, static_cast<int>(nAfter), maskOut + nBefore + i);
2647                                                            nEntries++;
2648                                                            assert (nEntries <= maxMasks);
2649                                                       }
2650                                                  }
2651                                                  char ** temp = masks;
2652                                                  masks = newMasks;
2653                                                  newMasks = temp;
2654                                             }
2655                                             // Now extend and sort
2656                                             int * sort = new int[nEntries];
2657                                             for (i = 0; i < nEntries; i++) {
2658                                                  char * maskThis = masks[i];
2659                                                  size_t length = strlen(maskThis);
2660                                                  while (length > 0 && maskThis[length-1] == ' ')
2661                                                       length--;
2662                                                  maskThis[length] = '\0';
2663                                                  sort[i] = static_cast<int>(length);
2664                                             }
2665                                             CoinSort_2(sort, sort + nEntries, masks);
2666                                             int lastLength = -1;
2667                                             for (i = 0; i < nEntries; i++) {
2668                                                  int length = sort[i];
2669                                                  while (length > lastLength)
2670                                                       maskStarts[++lastLength] = i;
2671                                             }
2672                                             maskStarts[++lastLength] = nEntries;
2673                                             delete [] sort;
2674                                             for (i = 0; i < maxMasks; i++)
2675                                                  delete [] newMasks[i];
2676                                             delete [] newMasks;
2677                                        }
2678                                        if (printMode > 5) {
2679                                          int numberColumns = models[iModel].numberColumns();
2680                                          // column length unless rhs ranging
2681                                          int number = numberColumns;
2682                                          switch (printMode) {
2683                                            // bound ranging
2684                                            case 6:
2685                                              fprintf(fp,"Bound ranging");
2686                                              break;
2687                                            // rhs ranging
2688                                            case 7:
2689                                              fprintf(fp,"Rhs ranging");
2690                                              number = numberRows;
2691                                              break;
2692                                            // objective ranging
2693                                            case 8:
2694                                              fprintf(fp,"Objective ranging");
2695                                              break;
2696                                          }
2697                                          if (lengthName)
2698                                            fprintf(fp,",name");
2699                                          fprintf(fp,",increase,variable,decrease,variable\n");
2700                                          int * which = new int [ number];
2701                                          if (printMode != 7) {
2702                                            if (!doMask) {
2703                                              for (int i = 0; i < number;i ++)
2704                                                which[i]=i;
2705                                            } else {
2706                                              int n = 0;
2707                                              for (int i = 0; i < number;i ++) {
2708                                                if (maskMatches(maskStarts,masks,columnNames[i]))
2709                                                  which[n++]=i;
2710                                              }
2711                                              if (n) {
2712                                                number=n;
2713                                              } else {
2714                                                printf("No names match - doing all\n");
2715                                                for (int i = 0; i < number;i ++)
2716                                                  which[i]=i;
2717                                              }
2718                                            }
2719                                          } else {
2720                                            if (!doMask) {
2721                                              for (int i = 0; i < number;i ++)
2722                                                which[i]=i+numberColumns;
2723                                            } else {
2724                                              int n = 0;
2725                                              for (int i = 0; i < number;i ++) {
2726                                                if (maskMatches(maskStarts,masks,rowNames[i]))
2727                                                  which[n++]=i+numberColumns;
2728                                              }
2729                                              if (n) {
2730                                                number=n;
2731                                              } else {
2732                                                printf("No names match - doing all\n");
2733                                                for (int i = 0; i < number;i ++)
2734                                                  which[i]=i+numberColumns;
2735                                              }
2736                                            }
2737                                          }
2738                                          double * valueIncrease = new double [ number];
2739                                          int * sequenceIncrease = new int [ number];
2740                                          double * valueDecrease = new double [ number];
2741                                          int * sequenceDecrease = new int [ number];
2742                                          switch (printMode) {
2743                                            // bound or rhs ranging
2744                                            case 6:
2745                                            case 7:
2746                                              models[iModel].primalRanging(numberRows,
2747                                                                           which, valueIncrease, sequenceIncrease,
2748                                                                           valueDecrease, sequenceDecrease);
2749                                              break;
2750                                            // objective ranging
2751                                            case 8:
2752                                              models[iModel].dualRanging(number,
2753                                                                           which, valueIncrease, sequenceIncrease,
2754                                                                           valueDecrease, sequenceDecrease);
2755                                              break;
2756                                          }
2757                                          for (int i = 0; i < number; i++) {
2758                                            int iWhich = which[i];
2759                                            fprintf(fp, "%d,", (iWhich<numberColumns) ? iWhich : iWhich-numberColumns);
2760                                            if (lengthName) {
2761                                              const char * name = (printMode==7) ? rowNames[iWhich-numberColumns].c_str() : columnNames[iWhich].c_str();
2762                                              fprintf(fp,"%s,",name);
2763                                            }
2764                                            if (valueIncrease[i]<1.0e30) {
2765                                              fprintf(fp, "%.10g,", valueIncrease[i]);
2766                                              int outSequence = sequenceIncrease[i];
2767                                              if (outSequence<numberColumns) {
2768                                                if (lengthName)
2769                                                  fprintf(fp,"%s,",columnNames[outSequence].c_str());
2770                                                else
2771                                                  fprintf(fp,"C%7.7d,",outSequence);
2772                                              } else {
2773                                                outSequence -= numberColumns;
2774                                                if (lengthName)
2775                                                  fprintf(fp,"%s,",rowNames[outSequence].c_str());
2776                                                else
2777                                                  fprintf(fp,"R%7.7d,",outSequence);
2778                                              }
2779                                            } else {
2780                                              fprintf(fp,"1.0e100,,");
2781                                            }
2782                                            if (valueDecrease[i]<1.0e30) {
2783                                              fprintf(fp, "%.10g,", valueDecrease[i]);
2784                                              int outSequence = sequenceDecrease[i];
2785                                              if (outSequence<numberColumns) {
2786                                                if (lengthName)
2787                                                  fprintf(fp,"%s",columnNames[outSequence].c_str());
2788                                                else
2789                                                  fprintf(fp,"C%7.7d",outSequence);
2790                                              } else {
2791                                                outSequence -= numberColumns;
2792                                                if (lengthName)
2793                                                  fprintf(fp,"%s",rowNames[outSequence].c_str());
2794                                                else
2795                                                  fprintf(fp,"R%7.7d",outSequence);
2796                                              }
2797                                            } else {
2798                                              fprintf(fp,"1.0e100,");
2799                                            }
2800                                            fprintf(fp,"\n");
2801                                          }
2802                                          if (fp != stdout)
2803                                            fclose(fp);
2804                                          delete [] which;
2805                                          delete [] valueIncrease;
2806                                          delete [] sequenceIncrease;
2807                                          delete [] valueDecrease;
2808                                          delete [] sequenceDecrease;
2809                                          if (masks) {
2810                                            delete [] maskStarts;
2811                                            for (int i = 0; i < maxMasks; i++)
2812                                              delete [] masks[i];
2813                                            delete [] masks;
2814                                          }
2815                                          break;
2816                                        }
2817                                        sprintf(printFormat," %s         %s\n",
2818                                                CLP_QUOTE(CLP_OUTPUT_FORMAT),
2819                                                CLP_QUOTE(CLP_OUTPUT_FORMAT));
2820                                        if (printMode > 2) {
2821                                             for (iRow = 0; iRow < numberRows; iRow++) {
2822                                                  int type = printMode - 3;
2823                                                  if (primalRowSolution[iRow] > rowUpper[iRow] + primalTolerance ||
2824                                                            primalRowSolution[iRow] < rowLower[iRow] - primalTolerance) {
2825                                                       fprintf(fp, "** ");
2826                                                       type = 2;
2827                                                  } else if (fabs(primalRowSolution[iRow]) > 1.0e-8) {
2828                                                       type = 1;
2829                                                  } else if (numberRows < 50) {
2830                                                       type = 3;
2831                                                  }
2832                                                  if (doMask && !maskMatches(maskStarts, masks, rowNames[iRow]))
2833                                                       type = 0;
2834                                                  if (type) {
2835                                                       fprintf(fp, "%7d ", iRow);
2836                                                       if (lengthName) {
2837                                                            const char * name = rowNames[iRow].c_str();
2838                                                            size_t n = strlen(name);
2839                                                            size_t i;
2840                                                            for (i = 0; i < n; i++)
2841                                                                 fprintf(fp, "%c", name[i]);
2842                                                            for (; i < static_cast<size_t>(lengthPrint); i++)
2843                                                                 fprintf(fp, " ");
2844                                                       }
2845                                                       fprintf(fp, printFormat, primalRowSolution[iRow],
2846                                                               dualRowSolution[iRow]);
2847                                                  }
2848                                             }
2849                                        }
2850                                        int iColumn;
2851                                        int numberColumns = models[iModel].numberColumns();
2852                                        double * dualColumnSolution =
2853                                             models[iModel].dualColumnSolution();
2854                                        double * primalColumnSolution =
2855                                             models[iModel].primalColumnSolution();
2856                                        double * columnLower = models[iModel].columnLower();
2857                                        double * columnUpper = models[iModel].columnUpper();
2858                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2859                                             int type = (printMode > 3) ? 1 : 0;
2860                                             if (primalColumnSolution[iColumn] > columnUpper[iColumn] + primalTolerance ||
2861                                                       primalColumnSolution[iColumn] < columnLower[iColumn] - primalTolerance) {
2862                                                  fprintf(fp, "** ");
2863                                                  type = 2;
2864                                             } else if (fabs(primalColumnSolution[iColumn]) > 1.0e-8) {
2865                                                  type = 1;
2866                                             } else if (numberColumns < 50) {
2867                                                  type = 3;
2868                                             }
2869                                             if (doMask && !maskMatches(maskStarts, masks,
2870                                                                        columnNames[iColumn]))
2871                                                  type = 0;
2872                                             if (type) {
2873                                                  fprintf(fp, "%7d ", iColumn);
2874                                                  if (lengthName) {
2875                                                       const char * name = columnNames[iColumn].c_str();
2876                                                       size_t n = strlen(name);
2877                                                       size_t i;
2878                                                       for (i = 0; i < n; i++)
2879                                                            fprintf(fp, "%c", name[i]);
2880                                                       for (; i < static_cast<size_t>(lengthPrint); i++)
2881                                                            fprintf(fp, " ");
2882                                                  }
2883                                                  fprintf(fp, printFormat, 
2884                                                          primalColumnSolution[iColumn],
2885                                                          dualColumnSolution[iColumn]);
2886                                             }
2887                                        }
2888                                        if (fp != stdout)
2889                                             fclose(fp);
2890                                        if (masks) {
2891                                             delete [] maskStarts;
2892                                             for (int i = 0; i < maxMasks; i++)
2893                                                  delete [] masks[i];
2894                                             delete [] masks;
2895                                        }
2896                                   } else {
2897                                        std::cout << "Unable to open file " << fileName << std::endl;
2898                                   }
2899                              } else {
2900                                   std::cout << "** Current model not valid" << std::endl;
2901
2902                              }
2903
2904                              break;
2905                         case CLP_PARAM_ACTION_SAVESOL:
2906                              if (goodModels[iModel]) {
2907                                   // get next field
2908                                   field = CoinReadGetString(argc, argv);
2909                                   if (field == "$") {
2910                                        field = parameters[iParam].stringValue();
2911                                   } else if (field == "EOL") {
2912                                        parameters[iParam].printString();
2913                                        break;
2914                                   } else {
2915                                        parameters[iParam].setStringValue(field);
2916                                   }
2917                                   std::string fileName;
2918                                   if (field[0] == '/' || field[0] == '\\') {
2919                                        fileName = field;
2920                                   } else if (field[0] == '~') {
2921                                        char * environVar = getenv("HOME");
2922                                        if (environVar) {
2923                                             std::string home(environVar);
2924                                             field = field.erase(0, 1);
2925                                             fileName = home + field;
2926                                        } else {
2927                                             fileName = field;
2928                                        }
2929                                   } else {
2930                                        fileName = directory + field;
2931                                   }
2932                                   saveSolution(models + iModel, fileName);
2933                              } else {
2934                                   std::cout << "** Current model not valid" << std::endl;
2935
2936                              }
2937                              break;
2938                         case CLP_PARAM_ACTION_ENVIRONMENT:
2939                              CbcOrClpEnvironmentIndex = 0;
2940                              break;
2941                         default:
2942                              abort();
2943                         }
2944                    }
2945               } else if (!numberMatches) {
2946                    std::cout << "No match for " << field << " - ? for list of commands"
2947                              << std::endl;
2948               } else if (numberMatches == 1) {
2949                    if (!numberQuery) {
2950                         std::cout << "Short match for " << field << " - completion: ";
2951                         std::cout << parameters[firstMatch].matchName() << std::endl;
2952                    } else if (numberQuery) {
2953                         std::cout << parameters[firstMatch].matchName() << " : ";
2954                         std::cout << parameters[firstMatch].shortHelp() << std::endl;
2955                         if (numberQuery >= 2)
2956                              parameters[firstMatch].printLongHelp();
2957                    }
2958               } else {
2959                    if (!numberQuery)
2960                         std::cout << "Multiple matches for " << field << " - possible completions:"
2961                                   << std::endl;
2962                    else
2963                         std::cout << "Completions of " << field << ":" << std::endl;
2964                    for ( iParam = 0; iParam < numberParameters; iParam++ ) {
2965                         int match = parameters[iParam].matches(field);
2966                         if (match && parameters[iParam].displayThis()) {
2967                              std::cout << parameters[iParam].matchName();
2968                              if (numberQuery >= 2)
2969                                   std::cout << " : " << parameters[iParam].shortHelp();
2970                              std::cout << std::endl;
2971                         }
2972                    }
2973               }
2974          }
2975          delete [] goodModels;
2976#ifdef COIN_HAS_GLPK
2977     if (cbc_glp_prob) {
2978       // free up as much as possible
2979       glp_free(cbc_glp_prob);
2980       glp_mpl_free_wksp(cbc_glp_tran);
2981       glp_free_env(); 
2982       cbc_glp_prob = NULL;
2983       cbc_glp_tran = NULL;
2984     }
2985#endif
2986     // By now all memory should be freed
2987#ifdef DMALLOC
2988     dmalloc_log_unfreed();
2989     dmalloc_shutdown();
2990#endif
2991#ifdef CLP_USEFUL_PRINTOUT
2992     printf("BENCHMARK_RESULT %s took %g cpu seconds (%g elapsed - %d iterations) - %d row, %d columns, %d elements - reduced to %d, %d, %d\n",
2993              mpsFile.c_str(),CoinCpuTime()-startCpu,CoinGetTimeOfDay()-startElapsed,
2994            debugInt[23],
2995              debugInt[0],debugInt[1],debugInt[2],
2996              debugInt[3],debugInt[4],debugInt[5]);
2997            const char * method[]={"","Dual","Primal","Sprint","Idiot"};
2998            printf("BENCHMARK using method %s (%d)\n",
2999              method[debugInt[6]],debugInt[7]);
3000#endif
3001     return 0;
3002}
3003static void breakdown(const char * name, CoinBigIndex numberLook, const double * region)
3004{
3005     double range[] = {
3006          -COIN_DBL_MAX,
3007          -1.0e15, -1.0e11, -1.0e8, -1.0e5, -1.0e4, -1.0e3, -1.0e2, -1.0e1,
3008          -1.0,
3009          -1.0e-1, -1.0e-2, -1.0e-3, -1.0e-4, -1.0e-5, -1.0e-8, -1.0e-11, -1.0e-15,
3010          0.0,
3011          1.0e-15, 1.0e-11, 1.0e-8, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1,
3012          1.0,
3013          1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15,
3014          COIN_DBL_MAX
3015     };
3016     int nRanges = static_cast<int> (sizeof(range) / sizeof(double));
3017     int * number = new int[nRanges];
3018     memset(number, 0, nRanges * sizeof(int));
3019     int * numberExact = new int[nRanges];
3020     memset(numberExact, 0, nRanges * sizeof(int));
3021     int i;
3022     for ( i = 0; i < numberLook; i++) {
3023          double value = region[i];
3024          for (int j = 0; j < nRanges; j++) {
3025               if (value == range[j]) {
3026                    numberExact[j]++;
3027                    break;
3028               } else if (value < range[j]) {
3029                    number[j]++;
3030                    break;
3031               }
3032          }
3033     }
3034     printf("\n%s has %d entries\n", name, numberLook);
3035     for (i = 0; i < nRanges; i++) {
3036          if (number[i])
3037               printf("%d between %g and %g", number[i], range[i-1], range[i]);
3038          if (numberExact[i]) {
3039               if (number[i])
3040                    printf(", ");
3041               printf("%d exactly at %g", numberExact[i], range[i]);
3042          }
3043          if (number[i] + numberExact[i])
3044               printf("\n");
3045     }
3046     delete [] number;
3047     delete [] numberExact;
3048}
3049void sortOnOther(int * column,
3050                 const CoinBigIndex * rowStart,
3051                 int * order,
3052                 int * other,
3053                 int nRow,
3054                 int nInRow,
3055                 int where)
3056{
3057     if (nRow < 2 || where >= nInRow)
3058          return;
3059     // do initial sort
3060     int kRow;
3061     int iRow;
3062     for ( kRow = 0; kRow < nRow; kRow++) {
3063          iRow = order[kRow];
3064          other[kRow] = column[rowStart[iRow] + where];
3065     }
3066     CoinSort_2(other, other + nRow, order);
3067     int first = 0;
3068     iRow = order[0];
3069     int firstC = column[rowStart[iRow] + where];
3070     kRow = 1;
3071     while (kRow < nRow) {
3072          int lastC = 9999999;;
3073          for (; kRow < nRow + 1; kRow++) {
3074               if (kRow < nRow) {
3075                    iRow = order[kRow];
3076                    lastC = column[rowStart[iRow] + where];
3077               } else {
3078                    lastC = 9999999;
3079               }
3080               if (lastC > firstC)
3081                    break;
3082          }
3083          // sort
3084          sortOnOther(column, rowStart, order + first, other, kRow - first,
3085                      nInRow, where + 1);
3086          firstC = lastC;
3087          first = kRow;
3088     }
3089}
3090static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
3091{
3092     int numberColumns = originalModel->numberColumns();
3093     const char * integerInformation  = originalModel->integerInformation();
3094     const double * columnLower = originalModel->columnLower();
3095     const double * columnUpper = originalModel->columnUpper();
3096     int numberIntegers = 0;
3097     int numberBinary = 0;
3098     int iRow, iColumn;
3099     if (integerInformation) {
3100          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3101               if (integerInformation[iColumn]) {
3102                    if (columnUpper[iColumn] > columnLower[iColumn]) {
3103                         numberIntegers++;
3104                         if (columnLower[iColumn] == 0.0 && columnUpper[iColumn] == 1)
3105                              numberBinary++;
3106                    }
3107               }
3108          }
3109          printf("Original problem has %d integers (%d of which binary)\n",
3110                 numberIntegers,numberBinary);
3111     }
3112     numberColumns = model->numberColumns();
3113     int numberRows = model->numberRows();
3114     columnLower = model->columnLower();
3115     columnUpper = model->columnUpper();
3116     const double * rowLower = model->rowLower();
3117     const double * rowUpper = model->rowUpper();
3118     const double * objective = model->objective();
3119     if (model->integerInformation()) {
3120       const char * integerInformation  = model->integerInformation();
3121       int numberIntegers = 0;
3122       int numberBinary = 0;
3123       double * obj = new double [numberColumns];
3124       int * which = new int [numberColumns];
3125       for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3126         if (columnUpper[iColumn] > columnLower[iColumn]) {
3127           if (integerInformation[iColumn]) {
3128             numberIntegers++;
3129             if (columnLower[iColumn] == 0.0 && columnUpper[iColumn] == 1)
3130               numberBinary++;
3131           }
3132         }
3133       }
3134       if(numberColumns != originalModel->numberColumns())
3135         printf("Presolved problem has %d integers (%d of which binary)\n",
3136                numberIntegers,numberBinary);
3137       for (int ifInt=0;ifInt<2;ifInt++) {
3138         for (int ifAbs=0;ifAbs<2;ifAbs++) {
3139           int numberSort=0;
3140           int numberZero=0;
3141           int numberDifferentObj=0;
3142           for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3143             if (columnUpper[iColumn] > columnLower[iColumn]) {
3144               if (!ifInt||integerInformation[iColumn]) {
3145                 obj[numberSort]=(ifAbs) ? fabs(objective[iColumn]) :
3146                   objective[iColumn];
3147                 which[numberSort++]=iColumn;
3148                 if (!objective[iColumn])
3149                   numberZero++;
3150               }
3151             }
3152           }
3153           CoinSort_2(obj,obj+numberSort,which);
3154           double last=obj[0];
3155           for (int jColumn = 1; jColumn < numberSort; jColumn++) {
3156             if (fabs(obj[jColumn]-last)>1.0e-12) {
3157               numberDifferentObj++;
3158               last=obj[jColumn];
3159             }
3160           }
3161           numberDifferentObj++;
3162           printf("==== ");
3163           if (ifInt)
3164             printf("for integers ");
3165           if (!ifAbs)
3166             printf("%d zero objective ",numberZero);
3167           else
3168             printf("absolute objective values ");
3169           printf("%d different\n",numberDifferentObj);
3170           bool saveModel=false;
3171           int target=model->logLevel();
3172           if (target>10000) {
3173             if (ifInt&&!ifAbs)
3174               saveModel=true;
3175             target-=10000;
3176           }
3177
3178           if (target<=100)
3179             target=12;
3180           else
3181             target-=100;
3182           if (numberDifferentObj<target) {
3183             int iLast=0;
3184             double last=obj[0];
3185             for (int jColumn = 1; jColumn < numberSort; jColumn++) {
3186               if (fabs(obj[jColumn]-last)>1.0e-12) {
3187                 printf("%d variables have objective of %g\n",
3188                        jColumn-iLast,last);
3189                 iLast=jColumn;
3190                 last=obj[jColumn];
3191               }
3192             }
3193             printf("%d variables have objective of %g\n",
3194                    numberSort-iLast,last);
3195             if (saveModel) {
3196               int spaceNeeded=numberSort+numberDifferentObj;
3197               int * columnAdd = new int[spaceNeeded];
3198               CoinBigIndex * startAdd = new CoinBigIndex[numberDifferentObj+1];
3199               double * elementAdd = new double[spaceNeeded];
3200               CoinBigIndex * rowAdd = new CoinBigIndex[2*numberDifferentObj+1];
3201               int * newIsInteger = reinterpret_cast<int *>(rowAdd+numberDifferentObj+1);
3202               double * objectiveNew = new double[3*numberDifferentObj];
3203               double * lowerNew = objectiveNew+numberDifferentObj;
3204               double * upperNew = lowerNew+numberDifferentObj;
3205               memset(startAdd,0,
3206                      (numberDifferentObj+1)*sizeof(CoinBigIndex));
3207               ClpSimplex tempModel=*model;
3208               int iLast=0;
3209               double last=obj[0];
3210               numberDifferentObj=0;
3211               int numberElements=0;
3212               rowAdd[0]=0;
3213               double * objective = tempModel.objective();
3214               for (int jColumn = 1; jColumn < numberSort+1; jColumn++) {
3215                 if (jColumn==numberSort||fabs(obj[jColumn]-last)>1.0e-12) {
3216                   // not if just one
3217                   if (jColumn-iLast>1) {
3218                     bool allInteger=integerInformation!=NULL;
3219                     int iColumn=which[iLast];
3220                     objectiveNew[numberDifferentObj]=objective[iColumn];
3221                     double lower=0.0;
3222                     double upper=0.0;
3223                     for (int kColumn=iLast;kColumn<jColumn;kColumn++) {
3224                       iColumn=which[kColumn];
3225                       objective[iColumn]=0.0;
3226                       double lowerValue=columnLower[iColumn];
3227                       double upperValue=columnUpper[iColumn];
3228                       double elementValue=-1.0;
3229                       if (objectiveNew[numberDifferentObj]*objective[iColumn]<0.0) {
3230                         lowerValue=-columnUpper[iColumn];
3231                         upperValue=-columnLower[iColumn];
3232                         elementValue=1.0;
3233                       }
3234                       columnAdd[numberElements]=iColumn;
3235                       elementAdd[numberElements++]=elementValue;
3236                       if (integerInformation&&!integerInformation[iColumn])
3237                         allInteger=false;
3238                       if (lower!=-COIN_DBL_MAX) {
3239                         if (lowerValue!=-COIN_DBL_MAX)
3240                           lower += lowerValue;
3241                         else
3242                           lower=-COIN_DBL_MAX;
3243                       }
3244                       if (upper!=COIN_DBL_MAX) {
3245                         if (upperValue!=COIN_DBL_MAX)
3246                           upper += upperValue;
3247                         else
3248                           upper=COIN_DBL_MAX;
3249                       }
3250                     }
3251                     columnAdd[numberElements]=numberColumns+numberDifferentObj;
3252                     elementAdd[numberElements++]=1.0;
3253                     newIsInteger[numberDifferentObj]= (allInteger) ? 1 : 0;
3254                     lowerNew[numberDifferentObj]=lower;
3255                     upperNew[numberDifferentObj]=upper;
3256                     numberDifferentObj++;
3257                     rowAdd[numberDifferentObj]=numberElements;
3258                   }
3259                   iLast=jColumn;
3260                   last=obj[jColumn];
3261                 }
3262               }
3263               // add columns
3264               tempModel.addColumns(numberDifferentObj, lowerNew, upperNew,
3265                                    objectiveNew,
3266                                    startAdd, NULL, NULL);
3267               // add constraints and make integer if all integer in group
3268               for (int iObj=0; iObj < numberDifferentObj; iObj++) {
3269                 lowerNew[iObj]=0.0;
3270                 upperNew[iObj]=0.0;
3271                 if (newIsInteger[iObj])
3272                   tempModel.setInteger(numberColumns+iObj);
3273               }
3274               tempModel.addRows(numberDifferentObj, lowerNew, upperNew,
3275                                 rowAdd,columnAdd,elementAdd);
3276               delete [] columnAdd;
3277               delete [] startAdd;
3278               delete [] elementAdd;
3279               delete [] rowAdd;
3280               delete [] objectiveNew;
3281               // save
3282               std::string tempName = model->problemName();
3283               if (ifInt)
3284                 tempName += "_int";
3285               if (ifAbs)
3286                 tempName += "_abs";
3287               tempName += ".mps";
3288               tempModel.writeMps(tempName.c_str());
3289             }
3290           }
3291         }
3292       }
3293       delete [] which;
3294       delete [] obj;
3295       printf("===== end objective counts\n");
3296     }
3297     CoinPackedMatrix * matrix = model->matrix();
3298     CoinBigIndex numberElements = matrix->getNumElements();
3299     const int * columnLength = matrix->getVectorLengths();
3300     //const CoinBigIndex * columnStart = matrix->getVectorStarts();
3301     const double * elementByColumn = matrix->getElements();
3302     int * number = new int[numberRows+1];
3303     memset(number, 0, (numberRows + 1)*sizeof(int));
3304     int numberObjSingletons = 0;
3305     /* cType
3306        0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
3307        8 0/1
3308     */
3309     int cType[9];
3310     std::string cName[] = {"0.0->inf,", "0.0->up,", "lo->inf,", "lo->up,", "free,", "fixed,", "-inf->0.0,",
3311                            "-inf->up,", "0.0->1.0"
3312                           };
3313     int nObjective = 0;
3314     memset(cType, 0, sizeof(cType));
3315     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3316          int length = columnLength[iColumn];
3317          if (length == 1 && objective[iColumn])
3318               numberObjSingletons++;
3319          number[length]++;
3320          if (objective[iColumn])
3321               nObjective++;
3322          if (columnLower[iColumn] > -1.0e20) {
3323               if (columnLower[iColumn] == 0.0) {
3324                    if (columnUpper[iColumn] > 1.0e20)
3325                         cType[0]++;
3326                    else if (columnUpper[iColumn] == 1.0)
3327                         cType[8]++;
3328                    else if (columnUpper[iColumn] == 0.0)
3329                         cType[5]++;
3330                    else
3331                         cType[1]++;
3332               } else {
3333                    if (columnUpper[iColumn] > 1.0e20)
3334                         cType[2]++;
3335                    else if (columnUpper[iColumn] == columnLower[iColumn])
3336                         cType[5]++;
3337                    else
3338                         cType[3]++;
3339               }
3340          } else {
3341               if (columnUpper[iColumn] > 1.0e20)
3342                    cType[4]++;
3343               else if (columnUpper[iColumn] == 0.0)
3344                    cType[6]++;
3345               else
3346                    cType[7]++;
3347          }
3348     }
3349     /* rType
3350        0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
3351        7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
3352     */
3353     int rType[13];
3354     std::string rName[] = {"E 0.0,", "E 1.0,", "E -1.0,", "E other,", "G 0.0,", "G 1.0,", "G other,",
3355                            "L 0.0,", "L 1.0,", "L other,", "Range 0.0->1.0,", "Range other,", "Free"
3356                           };
3357     memset(rType, 0, sizeof(rType));
3358     for (iRow = 0; iRow < numberRows; iRow++) {
3359          if (rowLower[iRow] > -1.0e20) {
3360               if (rowLower[iRow] == 0.0) {
3361                    if (rowUpper[iRow] > 1.0e20)
3362                         rType[4]++;
3363                    else if (rowUpper[iRow] == 1.0)
3364                         rType[10]++;
3365                    else if (rowUpper[iRow] == 0.0)
3366                         rType[0]++;
3367                    else
3368                         rType[11]++;
3369               } else if (rowLower[iRow] == 1.0) {
3370                    if (rowUpper[iRow] > 1.0e20)
3371                         rType[5]++;
3372                    else if (rowUpper[iRow] == rowLower[iRow])
3373                         rType[1]++;
3374                    else
3375                         rType[11]++;
3376               } else if (rowLower[iRow] == -1.0) {
3377                    if (rowUpper[iRow] > 1.0e20)
3378                         rType[6]++;
3379                    else if (rowUpper[iRow] == rowLower[iRow])
3380                         rType[2]++;
3381                    else
3382                         rType[11]++;
3383               } else {
3384                    if (rowUpper[iRow] > 1.0e20)
3385                         rType[6]++;
3386                    else if (rowUpper[iRow] == rowLower[iRow])
3387                         rType[3]++;
3388                    else
3389                         rType[11]++;
3390               }
3391          } else {
3392               if (rowUpper[iRow] > 1.0e20)
3393                    rType[12]++;
3394               else if (rowUpper[iRow] == 0.0)
3395                    rType[7]++;
3396               else if (rowUpper[iRow] == 1.0)
3397                    rType[8]++;
3398               else
3399                    rType[9]++;
3400          }
3401     }
3402     // Basic statistics
3403     printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
3404            numberRows, numberColumns, nObjective, numberElements);
3405     if (number[0] + number[1]) {
3406          printf("There are ");
3407          if (numberObjSingletons)
3408               printf("%d singletons with objective ", numberObjSingletons);
3409          int numberNoObj = number[1] - numberObjSingletons;
3410          if (numberNoObj)
3411               printf("%d singletons with no objective ", numberNoObj);
3412          if (number[0])
3413               printf("** %d columns have no entries", number[0]);
3414          printf("\n");
3415     }
3416     printf("Column breakdown:\n");
3417     int k;
3418     for (k = 0; k < static_cast<int> (sizeof(cType) / sizeof(int)); k++) {
3419          printf("%d of type %s ", cType[k], cName[k].c_str());
3420          if (((k + 1) % 3) == 0)
3421               printf("\n");
3422     }
3423     if ((k % 3) != 0)
3424          printf("\n");
3425     printf("Row breakdown:\n");
3426     for (k = 0; k < static_cast<int> (sizeof(rType) / sizeof(int)); k++) {
3427          printf("%d of type %s ", rType[k], rName[k].c_str());
3428          if (((k + 1) % 3) == 0)
3429               printf("\n");
3430     }
3431     if ((k % 3) != 0)
3432          printf("\n");
3433     //#define SYM
3434#ifndef SYM
3435     if (model->logLevel() < 2)
3436          return ;
3437#endif
3438     int kMax = model->logLevel() > 3 ? 10000000 : 10;
3439     k = 0;
3440     for (iRow = 1; iRow <= numberRows; iRow++) {
3441          if (number[iRow]) {
3442               k++;
3443               printf("%d columns have %d entries\n", number[iRow], iRow);
3444               if (k == kMax)
3445                    break;
3446          }
3447     }
3448     {
3449       int n = 0;
3450       int nLast=-1;
3451       int iLast=-1;
3452       int jRow=iRow;
3453       iRow++;
3454       for (; iRow <numberRows; iRow++) {
3455         if (number[iRow]) {
3456           n += number[iRow];
3457           nLast = number[iRow];
3458           iLast=iRow;
3459         }
3460       }
3461       if (n) {
3462         printf("\n    ... set logLevel >3 to see all ......\n\n");
3463         printf("%d columns > %d entries < %d\n", n-nLast, jRow,iLast);
3464         printf("%d column%s %d entries\n", nLast,
3465                nLast>1 ? "s have" : " has",iLast);
3466       }
3467     }
3468     delete [] number;
3469     printf("\n\n");
3470     if (model->logLevel() == 63
3471#ifdef SYM
3472               || true
3473#endif
3474        ) {
3475          // get column copy
3476          CoinPackedMatrix columnCopy = *matrix;
3477          const int * columnLength = columnCopy.getVectorLengths();
3478          number = new int[numberRows+1];
3479          memset(number, 0, (numberRows + 1)*sizeof(int));
3480          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3481               int length = columnLength[iColumn];
3482               number[length]++;
3483          }
3484          k = 0;
3485          for (iRow = 1; iRow <= numberRows; iRow++) {
3486               if (number[iRow]) {
3487                    k++;
3488               }
3489          }
3490          int * row = columnCopy.getMutableIndices();
3491          const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
3492          double * element = columnCopy.getMutableElements();
3493          int * order = new int[numberColumns];
3494          int * other = new int[numberColumns];
3495          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3496               int length = columnLength[iColumn];
3497               order[iColumn] = iColumn;
3498               other[iColumn] = length;
3499               CoinBigIndex start = columnStart[iColumn];
3500               CoinSort_2(row + start, row + start + length, element + start);
3501          }
3502          CoinSort_2(other, other + numberColumns, order);
3503          int jColumn = number[0] + number[1];
3504          for (iRow = 2; iRow <= numberRows; iRow++) {
3505               if (number[iRow]) {
3506                    printf("XX %d columns have %d entries\n", number[iRow], iRow);
3507                    int kColumn = jColumn + number[iRow];
3508                    sortOnOther(row, columnStart,
3509                                order + jColumn, other, number[iRow], iRow, 0);
3510                    // Now print etc
3511                    if (iRow < 500000) {
3512                         for (int lColumn = jColumn; lColumn < kColumn; lColumn++) {
3513                              iColumn = order[lColumn];
3514                              CoinBigIndex start = columnStart[iColumn];
3515                              if (model->logLevel() == 63) {
3516                                   printf("column %d %g <= ", iColumn, columnLower[iColumn]);
3517                                   for (CoinBigIndex i = start; i < start + iRow; i++)
3518                                        printf("( %d, %g) ", row[i], element[i]);
3519                                   printf("<= %g\n", columnUpper[iColumn]);
3520                              }
3521                         }
3522                    }
3523                    jColumn = kColumn;
3524               }
3525          }
3526          delete [] order;
3527          delete [] other;
3528          delete [] number;
3529     }
3530     // get row copy
3531     CoinPackedMatrix rowCopy = *matrix;
3532     rowCopy.reverseOrdering();
3533     const int * rowLength = rowCopy.getVectorLengths();
3534     number = new int[numberColumns+1];
3535     memset(number, 0, (numberColumns + 1)*sizeof(int));
3536     int logLevel=model->logLevel();
3537     bool reduceMaster=(logLevel&16)!=0;
3538     bool writeMatrices=(logLevel&8)!=0;
3539     bool morePrint=(logLevel&32)!=0;
3540     if (logLevel > 3) {
3541       // See what is needed
3542       // get column copy
3543       CoinPackedMatrix columnCopy = *matrix;
3544       const int * columnLength = columnCopy.getVectorLengths();
3545       const int * row = columnCopy.getIndices();
3546       const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
3547       const double * element = columnCopy.getElements();
3548       const double * elementByRow = rowCopy.getElements();
3549       const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
3550       const int * column = rowCopy.getIndices();
3551       int nPossibleZeroCost=0;
3552       int nPossibleNonzeroCost=0;
3553       for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3554         int length = columnLength[iColumn];
3555         if (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30) {
3556           if (length==1) {
3557             printf("Singleton free %d - cost %g\n",iColumn,objective[iColumn]);
3558           } else if (length==2) {
3559             int iRow0=row[columnStart[iColumn]];
3560             int iRow1=row[columnStart[iColumn]+1];
3561             double element0=element[columnStart[iColumn]];
3562             double element1=element[columnStart[iColumn]+1];
3563             int n0=rowLength[iRow0];
3564             int n1=rowLength[iRow1];
3565             printf("Doubleton free %d - cost %g - %g in %srow with %d entries and %g in %srow with %d entries\n",
3566                    iColumn,objective[iColumn],element0,(rowLower[iRow0]==rowUpper[iRow0]) ? "==" : "",n0,
3567                    element1,(rowLower[iRow1]==rowUpper[iRow1]) ? "==" : "",n1);
3568
3569           }
3570         }
3571         if (length==1) {
3572           int iRow=row[columnStart[iColumn]];
3573           double value=COIN_DBL_MAX;
3574           for (CoinBigIndex i=rowStart[iRow];i<rowStart[iRow]+rowLength[iRow];i++) {
3575             int jColumn=column[i];
3576             if (jColumn!=iColumn) {
3577               if (value!=elementByRow[i]) {
3578                 if (value==COIN_DBL_MAX) {
3579                   value=elementByRow[i];
3580                 } else {
3581                   value = -COIN_DBL_MAX;
3582                   break;
3583                 }
3584               }
3585             }
3586           }
3587           if (!objective[iColumn]) {
3588             if (morePrint) 
3589             printf("Singleton %d with no objective in row with %d elements - rhs %g,%g\n",iColumn,rowLength[iRow],rowLower[iRow],rowUpper[iRow]);
3590             nPossibleZeroCost++;
3591           } else if (value!=-COIN_DBL_MAX) {
3592             if (morePrint) 
3593               printf("Singleton %d (%s) with objective in row %d (%s) with %d equal elements - rhs %g,%g\n",iColumn,model->getColumnName(iColumn).c_str(),
3594                      iRow,model->getRowName(iRow).c_str(),
3595                      rowLength[iRow],rowLower[iRow],rowUpper[iRow]);
3596             nPossibleNonzeroCost++;
3597           }
3598         }
3599       }
3600       if (nPossibleZeroCost||nPossibleNonzeroCost)
3601         printf("%d singletons with zero cost, %d with valid cost\n",
3602                nPossibleZeroCost,nPossibleNonzeroCost);
3603       // look for DW
3604       int * blockStart = new int [3*(numberRows+numberColumns)+1];
3605       int * columnBlock = blockStart+numberRows;
3606       int * nextColumn = columnBlock+numberColumns;
3607       int * blockCount = nextColumn+numberColumns;
3608       int * blockEls = blockCount+numberRows+1;
3609       int * countIntegers = blockEls+numberRows;
3610       memset(countIntegers,0,numberColumns*sizeof(int));
3611       int direction[2]={-1,1};
3612       int bestBreak=-1;
3613       double bestValue=0.0;
3614       int iPass=0;
3615       int halfway=(numberRows+1)/2;
3616       int firstMaster=-1;
3617       int lastMaster=-2;
3618       while (iPass<2) {
3619         int increment=direction[iPass];
3620         int start= increment>0 ? 0 : numberRows-1;
3621         int stop=increment>0 ? numberRows : -1;
3622         int numberBlocks=0;
3623         int thisBestBreak=-1;
3624         double thisBestValue=COIN_DBL_MAX;
3625         int numberRowsDone=0;
3626         int numberMarkedColumns=0;
3627         int maximumBlockSize=0;
3628         for (int i=0;i<numberRows+2*numberColumns;i++) 
3629           blockStart[i]=-1;
3630         for (int i=0;i<numberRows+1;i++)
3631           blockCount[i]=0;
3632         for (int iRow=start;iRow!=stop;iRow+=increment) {
3633           int iBlock = -1;
3634           for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3635             int iColumn=column[j];
3636             int whichColumnBlock=columnBlock[iColumn];
3637             if (whichColumnBlock>=0) {
3638               // column marked
3639               if (iBlock<0) {
3640                 // put row in that block
3641                 iBlock=whichColumnBlock;
3642               } else if (iBlock!=whichColumnBlock) {
3643                 // merge
3644                 blockCount[iBlock]+=blockCount[whichColumnBlock];
3645                 blockCount[whichColumnBlock]=0;
3646                 int jColumn=blockStart[whichColumnBlock];
3647                 while (jColumn>=0) {
3648                   columnBlock[jColumn]=iBlock;
3649                   iColumn=jColumn;
3650                   jColumn=nextColumn[jColumn];
3651                 }
3652                 nextColumn[iColumn]=blockStart[iBlock];
3653                 blockStart[iBlock]=blockStart[whichColumnBlock];
3654                 blockStart[whichColumnBlock]=-1;
3655               }
3656             }
3657           }
3658           int n=numberMarkedColumns;
3659           if (iBlock<0) {
3660             //new block
3661             if (rowLength[iRow]) {
3662               numberBlocks++;
3663               iBlock=numberBlocks;
3664               int jColumn=column[rowStart[iRow]];
3665               columnBlock[jColumn]=iBlock;
3666               blockStart[iBlock]=jColumn;
3667               numberMarkedColumns++;
3668               for (CoinBigIndex j=rowStart[iRow]+1;j<rowStart[iRow]+rowLength[iRow];j++) {
3669                 int iColumn=column[j];
3670                 columnBlock[iColumn]=iBlock;
3671                 numberMarkedColumns++;
3672                 nextColumn[jColumn]=iColumn;
3673                 jColumn=iColumn;
3674               }
3675               blockCount[iBlock]=numberMarkedColumns-n;
3676             } else {
3677               // empty
3678               iBlock=numberRows;
3679             }
3680           } else {
3681             // put in existing block
3682             int jColumn=blockStart[iBlock];
3683             for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3684               int iColumn=column[j];
3685               assert (columnBlock[iColumn]<0||columnBlock[iColumn]==iBlock);
3686               if (columnBlock[iColumn]<0) {
3687                 columnBlock[iColumn]=iBlock;
3688                 numberMarkedColumns++;
3689                 nextColumn[iColumn]=jColumn;
3690                 jColumn=iColumn;
3691               }
3692             }
3693             blockStart[iBlock]=jColumn;
3694             blockCount[iBlock]+=numberMarkedColumns-n;
3695           }
3696           maximumBlockSize=CoinMax(maximumBlockSize,blockCount[iBlock]);
3697           numberRowsDone++;
3698           if (thisBestValue*numberRowsDone > maximumBlockSize&&numberRowsDone>halfway) { 
3699             thisBestBreak=iRow;
3700             thisBestValue=static_cast<double>(maximumBlockSize)/
3701               static_cast<double>(numberRowsDone);
3702           }
3703         }
3704         // If wanted minimize master rows
3705         if (reduceMaster)
3706           thisBestValue= (increment<0) ? thisBestBreak :
3707             numberRows-thisBestBreak;
3708         if (thisBestBreak==stop)
3709           thisBestValue=COIN_DBL_MAX;
3710         iPass++;
3711         if (iPass==1) {
3712           bestBreak=thisBestBreak;
3713           bestValue=thisBestValue;
3714         } else {
3715           if (bestValue<thisBestValue) {
3716             firstMaster=0;
3717             lastMaster=bestBreak;
3718           } else {
3719             firstMaster=thisBestBreak+1;
3720             lastMaster=numberRows;
3721           }
3722         }
3723       }
3724       if (firstMaster<=lastMaster) {
3725         if (firstMaster==lastMaster)
3726           printf("Total decomposition! - ");
3727         printf("%d master rows %d <= < %d\n",lastMaster-firstMaster,
3728                firstMaster,lastMaster);
3729         for (int i=0;i<numberRows+2*numberColumns;i++) 
3730           blockStart[i]=-1;
3731         for (int i=firstMaster;i<lastMaster;i++)
3732           blockStart[i]=-2;
3733         int firstRow=0;
3734         int numberBlocks=-1;
3735         while (true) {
3736           for (;firstRow<numberRows;firstRow++) {
3737             if (blockStart[firstRow]==-1)
3738               break;
3739           }
3740           if (firstRow==numberRows)
3741             break;
3742           int nRows=0;
3743           numberBlocks++;
3744           int numberStack=1;
3745           blockCount[0] = firstRow;
3746           while (numberStack) {
3747             int iRow=blockCount[--numberStack];
3748             for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3749               int iColumn=column[j];
3750               int iBlock=columnBlock[iColumn];
3751               if (iBlock<0) {
3752                 columnBlock[iColumn]=numberBlocks;
3753                 for (CoinBigIndex k=columnStart[iColumn];
3754                      k<columnStart[iColumn]+columnLength[iColumn];k++) {
3755                   int jRow=row[k];
3756                   int rowBlock=blockStart[jRow];
3757                   if (rowBlock==-1) {
3758                     nRows++;
3759                     blockStart[jRow]=numberBlocks;
3760                     blockCount[numberStack++]=jRow;
3761                   }
3762                 }
3763               }
3764             }
3765           }
3766           if (!nRows) {
3767             // empty!!
3768             numberBlocks--;
3769           }
3770           firstRow++;
3771         }
3772         // adjust
3773         numberBlocks++;
3774         for (int i=0;i<numberBlocks;i++) {
3775           blockCount[i]=0;
3776           nextColumn[i]=0;
3777         }
3778         int numberEmpty=0;
3779         int numberMaster=0;
3780         memset(blockEls,0,numberBlocks*sizeof(int));
3781         for (int iRow = 0; iRow < numberRows; iRow++) {
3782           int iBlock=blockStart[iRow];
3783           if (iBlock>=0) {
3784             blockCount[iBlock]++;
3785             blockEls[iBlock]+=rowLength[iRow];
3786           } else {
3787             if (iBlock==-2)
3788               numberMaster++;
3789             else
3790               numberEmpty++;
3791           }
3792         }
3793         int numberEmptyColumns=0;
3794         int numberMasterColumns=0;
3795         int numberMasterIntegers=0;
3796         for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3797           int iBlock=columnBlock[iColumn];
3798           bool isInteger = (model->isInteger(iColumn));
3799           if (iBlock>=0) {
3800             nextColumn[iBlock]++;
3801             if (isInteger)
3802               countIntegers[iBlock]++;
3803           } else {
3804             if (isInteger)
3805               numberMasterIntegers++;
3806             if (columnLength[iColumn])
3807               numberMasterColumns++;
3808             else
3809               numberEmptyColumns++;
3810           }
3811         }
3812         int largestRows=0;
3813         int largestColumns=0;
3814         for (int i=0;i<numberBlocks;i++) {
3815           if (blockCount[i]+nextColumn[i]>largestRows+largestColumns) {
3816             largestRows=blockCount[i];
3817             largestColumns=nextColumn[i];
3818           }
3819         }
3820         bool useful=true;
3821         if (numberMaster>halfway||largestRows*3>numberRows)
3822           useful=false;
3823         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",
3824                useful ? "**Useful" : "NoGood",
3825                numberBlocks,largestRows,largestColumns,numberMaster,numberEmpty,numberRows,
3826                numberMasterColumns,numberEmptyColumns,numberMasterIntegers,
3827                numberColumns);
3828         for (int i=0;i<numberBlocks;i++) 
3829           printf("Block %d has %d rows and %d columns (%d elements, %d integers)\n",
3830                  i,blockCount[i],nextColumn[i],blockEls[i],countIntegers[i]);
3831         FILE * fpBlocks = fopen("blocks.data","wb");
3832         printf("Blocks data on file blocks.data\n");
3833         int stats[3];
3834         stats[0]=numberRows;
3835         stats[1]=numberColumns;
3836         stats[2]=numberBlocks;
3837         size_t numberWritten;
3838         numberWritten=fwrite(stats,sizeof(int),3,fpBlocks);
3839         assert (numberWritten==3);
3840         numberWritten=fwrite(blockStart,sizeof(int),numberRows,fpBlocks);
3841         assert (numberWritten==numberRows);
3842         numberWritten=fwrite(columnBlock,sizeof(int),numberColumns,fpBlocks);
3843         assert (numberWritten==numberColumns);
3844         fclose(fpBlocks);
3845         if (writeMatrices) {
3846           int * whichRows=new int[numberRows+numberColumns];
3847           int * whichColumns=whichRows+numberRows;
3848           char name[20];
3849           for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
3850             sprintf(name,"block%d.mps",iBlock);
3851             int nRows=0;
3852             for (int iRow=0;iRow<numberRows;iRow++) {
3853               if (blockStart[iRow]==iBlock)
3854                 whichRows[nRows++]=iRow;
3855             }
3856             int nColumns=0;
3857             for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3858               if (columnBlock[iColumn]==iBlock)
3859                 whichColumns[nColumns++]=iColumn;
3860             }
3861             ClpSimplex subset(model,nRows,whichRows,nColumns,whichColumns);
3862             for (int jRow=0;jRow<nRows;jRow++) {
3863               int iRow = whichRows[jRow];
3864               std::string name = model->getRowName(iRow);
3865               subset.setRowName(jRow,name);
3866             }
3867             for (int jColumn=0;jColumn<nColumns;jColumn++) {
3868               int iColumn = whichColumns[jColumn];
3869               if (model->isInteger(iColumn))
3870                 subset.setInteger(jColumn);
3871               std::string name = model->getColumnName(iColumn);
3872               subset.setColumnName(jColumn,name);
3873             }
3874             subset.writeMps(name,0,1);
3875           }
3876           delete [] whichRows;
3877         } 
3878       }
3879       delete [] blockStart;
3880     }
3881     for (iRow = 0; iRow < numberRows; iRow++) {
3882          int length = rowLength[iRow];
3883          number[length]++;
3884     }
3885     if (number[0])
3886          printf("** %d rows have no entries\n", number[0]);
3887     k = 0;
3888     for (iColumn = 1; iColumn <= numberColumns; iColumn++) {
3889          if (number[iColumn]) {
3890               k++;
3891               printf("%d rows have %d entries\n", number[iColumn], iColumn);
3892               if (k == kMax)
3893                    break;
3894          }
3895     }
3896     {
3897       int n = 0;
3898       int nLast=-1;
3899       int iLast=-1;
3900       int jColumn=iColumn;
3901       iColumn++;
3902       for (; iColumn <numberColumns; iColumn++) {
3903         if (number[iColumn]) {
3904           n += number[iColumn];
3905           nLast = number[iColumn];
3906           iLast=iColumn;
3907         }
3908       }
3909       if (n) {
3910         printf("\n    ... set logLevel >3 to see all ......\n\n");
3911         printf("%d rows > %d entries < %d\n", n-nLast, jColumn,iLast);
3912         printf("%d row%s %d entries\n", nLast,
3913                nLast>1 ? "s have" : " has",iLast);
3914       }
3915     }
3916     if (morePrint
3917#ifdef SYM
3918               || true
3919#endif
3920        ) {
3921          int * column = rowCopy.getMutableIndices();
3922          const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
3923          double * element = rowCopy.getMutableElements();
3924          int * order = new int[numberRows];
3925          int * other = new int[numberRows];
3926          for (iRow = 0; iRow < numberRows; iRow++) {
3927               int length = rowLength[iRow];
3928               order[iRow] = iRow;
3929               other[iRow] = length;
3930               CoinBigIndex start = rowStart[iRow];
3931               CoinSort_2(column + start, column + start + length, element + start);
3932          }
3933          CoinSort_2(other, other + numberRows, order);
3934          int jRow = number[0] + number[1];
3935          double * weight = new double[numberRows];
3936          double * randomColumn = new double[numberColumns+1];
3937          double * randomRow = new double [numberRows+1];
3938          int * sortRow = new int [numberRows];
3939          int * possibleRow = new int [numberRows];
3940          int * backRow = new int [numberRows];
3941          int * stackRow = new int [numberRows];
3942          int * sortColumn = new int [numberColumns];
3943          int * possibleColumn = new int [numberColumns];
3944          int * backColumn = new int [numberColumns];
3945          int * backColumn2 = new int [numberColumns];
3946          int * mapRow = new int [numberRows];
3947          int * mapColumn = new int [numberColumns];
3948          int * stackColumn = new int [numberColumns];
3949          double randomLower = CoinDrand48();
3950          double randomUpper = CoinDrand48();
3951          double randomInteger = CoinDrand48();
3952          CoinBigIndex * startAdd = new CoinBigIndex[numberRows+1];
3953          int * columnAdd = new int [2*numberElements];
3954          double * elementAdd = new double[2*numberElements];
3955          int nAddRows = 0;
3956          startAdd[0] = 0;
3957          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3958               randomColumn[iColumn] = CoinDrand48();
3959               backColumn2[iColumn] = -1;
3960          }
3961          for (iColumn = 2; iColumn <= numberColumns; iColumn++) {
3962               if (number[iColumn]) {
3963                    printf("XX %d rows have %d entries\n", number[iColumn], iColumn);
3964                    int kRow = jRow + number[iColumn];
3965                    sortOnOther(column, rowStart,
3966                                order + jRow, other, number[iColumn], iColumn, 0);
3967                    // Now print etc
3968                    if (iColumn < 500000) {
3969                         int nLook = 0;
3970                         for (int lRow = jRow; lRow < kRow; lRow++) {
3971                              iRow = order[lRow];
3972                              CoinBigIndex start = rowStart[iRow];
3973                              if (morePrint) {
3974                                   printf("row %d %g <= ", iRow, rowLower[iRow]);
3975                                   for (CoinBigIndex i = start; i < start + iColumn; i++)
3976                                        printf("( %d, %g) ", column[i], element[i]);
3977                                   printf("<= %g\n", rowUpper[iRow]);
3978                              }
3979                              int first = column[start];
3980                              double sum = 0.0;
3981                              for (CoinBigIndex i = start; i < start + iColumn; i++) {
3982                                   int jColumn = column[i];
3983                                   double value = element[i];
3984                                   jColumn -= first;
3985                                   assert (jColumn >= 0);
3986                                   sum += value * randomColumn[jColumn];
3987                              }
3988                              if (rowLower[iRow] > -1.0e30 && rowLower[iRow])
3989                                   sum += rowLower[iRow] * randomLower;
3990                              else if (!rowLower[iRow])
3991                                   sum += 1.234567e-7 * randomLower;
3992                              if (rowUpper[iRow] < 1.0e30 && rowUpper[iRow])
3993                                   sum += rowUpper[iRow] * randomUpper;
3994                              else if (!rowUpper[iRow])
3995                                   sum += 1.234567e-7 * randomUpper;
3996                              sortRow[nLook] = iRow;
3997                              randomRow[nLook++] = sum;
3998                              // best way is to number unique elements and bounds and use
3999                              if (fabs(sum) > 1.0e4)
4000                                   sum *= 1.0e-6;
4001                              weight[iRow] = sum;
4002                         }
4003                         assert (nLook <= numberRows);
4004                         CoinSort_2(randomRow, randomRow + nLook, sortRow);
4005                         randomRow[nLook] = COIN_DBL_MAX;
4006                         double last = -COIN_DBL_MAX;
4007                         int iLast = -1;
4008                         for (int iLook = 0; iLook < nLook + 1; iLook++) {
4009                              if (randomRow[iLook] > last) {
4010                                   if (iLast >= 0) {
4011                                        int n = iLook - iLast;
4012                                        if (n > 1) {
4013                                             //printf("%d rows possible?\n",n);
4014                                        }
4015                                   }
4016                                   iLast = iLook;
4017                                   last = randomRow[iLook];
4018                              }
4019                         }
4020                    }
4021                    jRow = kRow;
4022               }
4023          }
4024          CoinPackedMatrix columnCopy = *matrix;
4025          const int * columnLength = columnCopy.getVectorLengths();
4026          const int * row = columnCopy.getIndices();
4027          const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
4028          const double * elementByColumn = columnCopy.getElements();
4029          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4030               int length = columnLength[iColumn];
4031               CoinBigIndex start = columnStart[iColumn];
4032               double sum = objective[iColumn];
4033               if (columnLower[iColumn] > -1.0e30 && columnLower[iColumn])
4034                    sum += columnLower[iColumn] * randomLower;
4035               else if (!columnLower[iColumn])
4036                    sum += 1.234567e-7 * randomLower;
4037               if (columnUpper[iColumn] < 1.0e30 && columnUpper[iColumn])
4038                    sum += columnUpper[iColumn] * randomUpper;
4039               else if (!columnUpper[iColumn])
4040                    sum += 1.234567e-7 * randomUpper;
4041               if (model->isInteger(iColumn))
4042                    sum += 9.87654321e-6 * randomInteger;
4043               for (CoinBigIndex i = start; i < start + length; i++) {
4044                    int iRow = row[i];
4045                    sum += elementByColumn[i] * weight[iRow];
4046               }
4047               sortColumn[iColumn] = iColumn;
4048               randomColumn[iColumn] = sum;
4049          }
4050          {
4051               CoinSort_2(randomColumn, randomColumn + numberColumns, sortColumn);
4052               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4053                    int i = sortColumn[iColumn];
4054                    backColumn[i] = iColumn;
4055               }
4056               randomColumn[numberColumns] = COIN_DBL_MAX;
4057               double last = -COIN_DBL_MAX;
4058               int iLast = -1;
4059               for (int iLook = 0; iLook < numberColumns + 1; iLook++) {
4060                    if (randomColumn[iLook] > last) {
4061                         if (iLast >= 0) {
4062                              int n = iLook - iLast;
4063                              if (n > 1) {
4064                                   //printf("%d columns possible?\n",n);
4065                              }
4066                              for (int i = iLast; i < iLook; i++) {
4067                                   possibleColumn[sortColumn[i]] = n;
4068                              }
4069                         }
4070                         iLast = iLook;
4071                         last = randomColumn[iLook];
4072                    }
4073               }
4074               for (iRow = 0; iRow < numberRows; iRow++) {
4075                    CoinBigIndex start = rowStart[iRow];
4076                    double sum = 0.0;
4077                    int length = rowLength[iRow];
4078                    for (CoinBigIndex i = start; i < start + length; i++) {
4079                         int jColumn = column[i];
4080                         double value = element[i];
4081                         jColumn = backColumn[jColumn];
4082                         sum += value * randomColumn[jColumn];
4083                         //if (iColumn==23089||iRow==23729)
4084                         //printf("row %d cola %d colb %d value %g rand %g sum %g\n",
4085                         //   iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
4086                    }
4087                    sortRow[iRow] = iRow;
4088                    randomRow[iRow] = weight[iRow];
4089                    randomRow[iRow] = sum;
4090               }
4091               CoinSort_2(randomRow, randomRow + numberRows, sortRow);
4092               for (iRow = 0; iRow < numberRows; iRow++) {
4093                    int i = sortRow[iRow];
4094                    backRow[i] = iRow;
4095               }
4096               randomRow[numberRows] = COIN_DBL_MAX;
4097               last = -COIN_DBL_MAX;
4098               iLast = -1;
4099               // Do backward indices from order
4100               for (iRow = 0; iRow < numberRows; iRow++) {
4101                    other[order[iRow]] = iRow;
4102               }
4103               for (int iLook = 0; iLook < numberRows + 1; iLook++) {
4104                    if (randomRow[iLook] > last) {
4105                         if (iLast >= 0) {
4106                              int n = iLook - iLast;
4107                              if (n > 1) {
4108                                   //printf("%d rows possible?\n",n);
4109                                   // Within group sort as for original "order"
4110                                   for (int i = iLast; i < iLook; i++) {
4111                                        int jRow = sortRow[i];
4112                                        order[i] = other[jRow];
4113                                   }
4114                                   CoinSort_2(order + iLast, order + iLook, sortRow + iLast);
4115                              }
4116                              for (int i = iLast; i < iLook; i++) {
4117                                   possibleRow[sortRow[i]] = n;
4118                              }
4119                         }
4120                         iLast = iLook;
4121                         last = randomRow[iLook];
4122                    }
4123               }
4124               // Temp out
4125               for (int iLook = 0; iLook < numberRows - 1000000; iLook++) {
4126                    iRow = sortRow[iLook];
4127                    CoinBigIndex start = rowStart[iRow];
4128                    int length = rowLength[iRow];
4129                    int numberPossible = possibleRow[iRow];
4130                    for (CoinBigIndex i = start; i < start + length; i++) {
4131                         int jColumn = column[i];
4132                         if (possibleColumn[jColumn] != numberPossible)
4133                              numberPossible = -1;
4134                    }
4135                    int n = numberPossible;
4136                    if (numberPossible > 1) {
4137                         //printf("pppppossible %d\n",numberPossible);
4138                         for (int jLook = iLook + 1; jLook < iLook + numberPossible; jLook++) {
4139                              int jRow = sortRow[jLook];
4140                              CoinBigIndex start2 = rowStart[jRow];
4141                              assert (numberPossible == possibleRow[jRow]);
4142                              assert(length == rowLength[jRow]);
4143                              for (CoinBigIndex i = start2; i < start2 + length; i++) {
4144                                   int jColumn = column[i];
4145                                   if (possibleColumn[jColumn] != numberPossible)
4146                                        numberPossible = -1;
4147                              }
4148                         }
4149                         if (numberPossible < 2) {
4150                              // switch off
4151                              for (int jLook = iLook; jLook < iLook + n; jLook++)
4152                                   possibleRow[sortRow[jLook]] = -1;
4153                         }
4154                         // skip rest
4155                         iLook += n - 1;
4156                    } else {
4157                         possibleRow[iRow] = -1;
4158                    }
4159               }
4160               for (int iLook = 0; iLook < numberRows; iLook++) {
4161                    iRow = sortRow[iLook];
4162                    int numberPossible = possibleRow[iRow];
4163                    // Only if any integers
4164                    int numberIntegers = 0;
4165                    CoinBigIndex start = rowStart[iRow];
4166                    int length = rowLength[iRow];
4167                    for (CoinBigIndex i = start; i < start + length; i++) {
4168                         int jColumn = column[i];
4169                         if (model->isInteger(jColumn))
4170                              numberIntegers++;
4171                    }
4172                    if (numberPossible > 1 && !numberIntegers) {
4173                         //printf("possible %d - but no integers\n",numberPossible);
4174                    }
4175                    if (numberPossible > 1 && (numberIntegers || false)) {
4176                         //
4177                         printf("possible %d - %d integers\n", numberPossible, numberIntegers);
4178                         int lastLook = iLook;
4179                         int nMapRow = -1;
4180                         for (int jLook = iLook + 1; jLook < iLook + numberPossible; jLook++) {
4181                              // stop if too many failures
4182                              if (jLook > iLook + 10 && nMapRow < 0)
4183                                   break;
4184                              // Create identity mapping
4185                              int i;
4186                              for (i = 0; i < numberRows; i++)
4187                                   mapRow[i] = i;
4188                              for (i = 0; i < numberColumns; i++)
4189                                   mapColumn[i] = i;
4190                              int offset = jLook - iLook;
4191                              int nStackC = 0;
4192                              // build up row and column mapping
4193                              int nStackR = 1;
4194                              stackRow[0] = iLook;
4195                              bool good = true;
4196                              while (nStackR) {
4197                                   nStackR--;
4198                                   int look1 = stackRow[nStackR];
4199                                   int look2 = look1 + offset;
4200                                   assert (randomRow[look1] == randomRow[look2]);
4201                                   int row1 = sortRow[look1];
4202                                   int row2 = sortRow[look2];
4203                                   assert (mapRow[row1] == row1);
4204                                   assert (mapRow[row2] == row2);
4205                                   mapRow[row1] = row2;
4206                                   mapRow[row2] = row1;
4207                                   CoinBigIndex start1 = rowStart[row1];
4208                                   CoinBigIndex offset2 = rowStart[row2] - start1;
4209                                   int length = rowLength[row1];
4210                                   assert( length == rowLength[row2]);
4211                                   for (CoinBigIndex i = start1; i < start1 + length; i++) {
4212                                        int jColumn1 = column[i];
4213                                        int jColumn2 = column[i+offset2];
4214                                        if (randomColumn[backColumn[jColumn1]] !=
4215                                                  randomColumn[backColumn[jColumn2]]) {
4216                                             good = false;
4217                                             break;
4218                                        }
4219                                        if (mapColumn[jColumn1] == jColumn1) {
4220                                             // not touched
4221                                             assert (mapColumn[jColumn2] == jColumn2);
4222                                             if (jColumn1 != jColumn2) {
4223                                                  // Put on stack
4224                                                  mapColumn[jColumn1] = jColumn2;
4225                                                  mapColumn[jColumn2] = jColumn1;
4226                                                  stackColumn[nStackC++] = jColumn1;
4227                                             }
4228                                        } else {
4229                                             if (mapColumn[jColumn1] != jColumn2 ||
4230                                                       mapColumn[jColumn2] != jColumn1) {
4231                                                  // bad
4232                                                  good = false;
4233                                                  printf("bad col\n");
4234                                                  break;
4235                                             }
4236                                        }
4237                                   }
4238                                   if (!good)
4239                                        break;
4240                                   while (nStackC) {
4241                                        nStackC--;
4242                                        int iColumn = stackColumn[nStackC];
4243                                        int iColumn2 = mapColumn[iColumn];
4244                                        assert (iColumn != iColumn2);
4245                                        int length = columnLength[iColumn];
4246                                        assert (length == columnLength[iColumn2]);
4247                                        CoinBigIndex start = columnStart[iColumn];
4248                                        CoinBigIndex offset2 = columnStart[iColumn2] - start;
4249                                        for (CoinBigIndex i = start; i < start + length; i++) {
4250                                             int iRow = row[i];
4251                                             int iRow2 = row[i+offset2];
4252                                             if (mapRow[iRow] == iRow) {
4253                                                  // First (but be careful)
4254                                                  if (iRow != iRow2) {
4255                                                       //mapRow[iRow]=iRow2;
4256                                                       //mapRow[iRow2]=iRow;
4257                                                       int iBack = backRow[iRow];
4258                                                       int iBack2 = backRow[iRow2];
4259                                                       if (randomRow[iBack] == randomRow[iBack2] &&
4260                                                                 iBack2 - iBack == offset) {
4261                                                            stackRow[nStackR++] = iBack;
4262                                                       } else {
4263                                                            //printf("randomRow diff - weights %g %g\n",
4264                                                            //     weight[iRow],weight[iRow2]);
4265                                                            // bad
4266                                                            good = false;
4267                                                            break;
4268                                                       }
4269                                                  }
4270                                             } else {
4271                                                  if (mapRow[iRow] != iRow2 ||
4272                                                            mapRow[iRow2] != iRow) {
4273                                                       // bad
4274                                                       good = false;
4275                                                       printf("bad row\n");
4276                                                       break;
4277                                                  }
4278                                             }
4279                                        }
4280                                        if (!good)
4281                                             break;
4282                                   }
4283                              }
4284                              // then check OK
4285                              if (good) {
4286                                   for (iRow = 0; iRow < numberRows; iRow++) {
4287                                        CoinBigIndex start = rowStart[iRow];
4288                                        int length = rowLength[iRow];
4289                                        if (mapRow[iRow] == iRow) {
4290                                             for (CoinBigIndex i = start; i < start + length; i++) {
4291                                                  int jColumn = column[i];
4292                                                  backColumn2[jColumn] = static_cast<int>(i - start);
4293                                             }
4294                                             for (CoinBigIndex i = start; i < start + length; i++) {
4295                                                  int jColumn = column[i];
4296                                                  if (mapColumn[jColumn] != jColumn) {
4297                                                       int jColumn2 = mapColumn[jColumn];
4298                                                       CoinBigIndex i2 = backColumn2[jColumn2];
4299                                                       if (i2 < 0) {
4300                                                            good = false;
4301                                                       } else if (element[i] != element[i2+start]) {
4302                                                            good = false;
4303                                                       }
4304                                                  }
4305                                             }
4306                                             for (CoinBigIndex i = start; i < start + length; i++) {
4307                                                  int jColumn = column[i];
4308                                                  backColumn2[jColumn] = -1;
4309                                             }
4310                                        } else {
4311                                             int row2 = mapRow[iRow];
4312                                             assert (iRow = mapRow[row2]);
4313                                             if (rowLower[iRow] != rowLower[row2] ||
4314                                                       rowLower[row2] != rowLower[iRow])
4315                                                  good = false;
4316                                             CoinBigIndex offset2 = rowStart[row2] - start;
4317                                             for (CoinBigIndex i = start; i < start + length; i++) {
4318                                                  int jColumn = column[i];
4319                                                  double value = element[i];
4320                                                  int jColumn2 = column[i+offset2];
4321                                                  double value2 = element[i+offset2];
4322                                                  if (value != value2 || mapColumn[jColumn] != jColumn2 ||
4323                                                            mapColumn[jColumn2] != jColumn)
4324                                                       good = false;
4325                                             }
4326                                        }
4327                                   }
4328                                   if (good) {
4329                                        // check rim
4330                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4331                                             if (mapColumn[iColumn] != iColumn) {
4332                                                  int iColumn2 = mapColumn[iColumn];
4333                                                  if (objective[iColumn] != objective[iColumn2])
4334                                                       good = false;
4335                                                  if (columnLower[iColumn] != columnLower[iColumn2])
4336                                                       good = false;
4337                                                  if (columnUpper[iColumn] != columnUpper[iColumn2])
4338                                                       good = false;
4339                                                  if (model->isInteger(iColumn) != model->isInteger(iColumn2))
4340                                                       good = false;
4341                                             }
4342                                        }
4343                                   }
4344                                   if (good) {
4345                                        // temp
4346                                        if (nMapRow < 0) {
4347                                             //const double * solution = model->primalColumnSolution();
4348                                             // find mapped
4349                                             int nMapColumn = 0;
4350                                             for (int i = 0; i < numberColumns; i++) {
4351                                                  if (mapColumn[i] > i)
4352                                                       nMapColumn++;
4353                                             }
4354                                             nMapRow = 0;
4355                                             int kRow = -1;
4356                                             for (int i = 0; i < numberRows; i++) {
4357                                                  if (mapRow[i] > i) {
4358                                                       nMapRow++;
4359                                                       kRow = i;
4360                                                  }
4361                                             }
4362                                             printf("%d columns, %d rows\n", nMapColumn, nMapRow);
4363                                             if (nMapRow == 1) {
4364                                                  CoinBigIndex start = rowStart[kRow];
4365                                                  int length = rowLength[kRow];
4366                                                  printf("%g <= ", rowLower[kRow]);
4367                                                  for (CoinBigIndex i = start; i < start + length; i++) {
4368                                                       int jColumn = column[i];
4369                                                       if (mapColumn[jColumn] != jColumn)
4370                                                            printf("* ");
4371                                                       printf("%d,%g ", jColumn, element[i]);
4372                                                  }
4373                                                  printf("<= %g\n", rowUpper[kRow]);
4374                                             }
4375                                        }
4376                                        // temp
4377                                        int row1 = sortRow[lastLook];
4378                                        int row2 = sortRow[jLook];
4379                                        lastLook = jLook;
4380                                        CoinBigIndex start1 = rowStart[row1];
4381                                        CoinBigIndex offset2 = rowStart[row2] - start1;
4382                                        int length = rowLength[row1];
4383                                        assert( length == rowLength[row2]);
4384                                        CoinBigIndex put = startAdd[nAddRows];
4385                                        double multiplier = length < 11 ? 2.0 : 1.125;
4386                                        double value = 1.0;
4387                                        for (CoinBigIndex i = start1; i < start1 + length; i++) {
4388                                             int jColumn1 = column[i];
4389                                             int jColumn2 = column[i+offset2];
4390                                             columnAdd[put] = jColumn1;
4391                                             elementAdd[put++] = value;
4392                                             columnAdd[put] = jColumn2;
4393                                             elementAdd[put++] = -value;
4394                                             value *= multiplier;
4395                                        }
4396                                        nAddRows++;
4397                                        startAdd[nAddRows] = put;
4398                                   } else {
4399                                        printf("ouch - did not check out as good\n");
4400                                   }
4401                              }
4402                         }
4403                         // skip rest
4404                         iLook += numberPossible - 1;
4405                    }
4406               }
4407          }
4408          if (nAddRows) {
4409               double * lower = new double [nAddRows];
4410               double * upper = new double[nAddRows];
4411               int i;
4412               //const double * solution = model->primalColumnSolution();
4413               for (i = 0; i < nAddRows; i++) {
4414                    lower[i] = 0.0;
4415                    upper[i] = COIN_DBL_MAX;
4416               }
4417               printf("Adding %d rows with %d elements\n", nAddRows,
4418                      startAdd[nAddRows]);
4419               //ClpSimplex newModel(*model);
4420               //newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
4421               //newModel.writeMps("modified.mps");
4422               delete [] lower;
4423               delete [] upper;
4424          }
4425          delete [] startAdd;
4426          delete [] columnAdd;
4427          delete [] elementAdd;
4428          delete [] order;
4429          delete [] other;
4430          delete [] randomColumn;
4431          delete [] weight;
4432          delete [] randomRow;
4433          delete [] sortRow;
4434          delete [] backRow;
4435          delete [] possibleRow;
4436          delete [] sortColumn;
4437          delete [] backColumn;
4438          delete [] backColumn2;
4439          delete [] possibleColumn;
4440          delete [] mapRow;
4441          delete [] mapColumn;
4442          delete [] stackRow;
4443          delete [] stackColumn;
4444     }
4445     delete [] number;
4446     // Now do breakdown of ranges
4447     breakdown("Elements", numberElements, elementByColumn);
4448     breakdown("RowLower", numberRows, rowLower);
4449     breakdown("RowUpper", numberRows, rowUpper);
4450     breakdown("ColumnLower", numberColumns, columnLower);
4451     breakdown("ColumnUpper", numberColumns, columnUpper);
4452     breakdown("Objective", numberColumns, objective);
4453     // do integer objective
4454     double * obj = CoinCopyOfArray(objective,numberColumns);
4455     int n=0;
4456     //#define FIX_COSTS 1.0
4457#ifdef FIX_COSTS
4458     double * obj2   = originalModel->objective();
4459     double * upper2   = originalModel->columnUpper();
4460#endif
4461     for (int i=0;i<numberColumns;i++) {
4462       if (integerInformation&&integerInformation[i]) {
4463         obj[n++]=obj[i];
4464#ifdef FIX_COSTS
4465         if (obj2[i]>FIX_COSTS)
4466         upper2[i]=0.0;
4467#endif
4468       }
4469     }
4470     if (n) {
4471       breakdown("Integer objective", n, obj);
4472     }
4473}
4474static bool maskMatches(const int * starts, char ** masks,
4475                        std::string & check)
4476{
4477     // back to char as I am old fashioned
4478     const char * checkC = check.c_str();
4479     size_t length = strlen(checkC);
4480     while (checkC[length-1] == ' ')
4481          length--;
4482     for (int i = starts[length]; i < starts[length+1]; i++) {
4483          char * thisMask = masks[i];
4484          size_t k;
4485          for ( k = 0; k < length; k++) {
4486               if (thisMask[k] != '?' && thisMask[k] != checkC[k])
4487                    break;
4488          }
4489          if (k == length)
4490               return true;
4491     }
4492     return false;
4493}
4494static void clean(char * temp)
4495{
4496     char * put = temp;
4497     while (*put >= ' ')
4498          put++;
4499     *put = '\0';
4500}
4501static void generateCode(const char * fileName, int type)
4502{
4503     FILE * fp = fopen(fileName, "r");
4504     assert (fp);
4505     int numberLines = 0;
4506#define MAXLINES 500
4507#define MAXONELINE 200
4508     char line[MAXLINES][MAXONELINE];
4509     while (fgets(line[numberLines], MAXONELINE, fp)) {
4510          assert (numberLines < MAXLINES);
4511          clean(line[numberLines]);
4512          numberLines++;
4513     }
4514     fclose(fp);
4515     // add in actual solve
4516     strcpy(line[numberLines], "5  clpModel->initialSolve(clpSolve);");
4517     numberLines++;
4518     fp = fopen(fileName, "w");
4519     assert (fp);
4520     char apo = '"';
4521     char backslash = '\\';
4522
4523     fprintf(fp, "#include %cClpSimplex.hpp%c\n", apo, apo);
4524     fprintf(fp, "#include %cClpSolve.hpp%c\n", apo, apo);
4525
4526     fprintf(fp, "\nint main (int argc, const char *argv[])\n{\n");
4527     fprintf(fp, "  ClpSimplex  model;\n");
4528     fprintf(fp, "  int status=1;\n");
4529     fprintf(fp, "  if (argc<2)\n");
4530     fprintf(fp, "    fprintf(stderr,%cPlease give file name%cn%c);\n",
4531             apo, backslash, apo);
4532     fprintf(fp, "  else\n");
4533     fprintf(fp, "    status=model.readMps(argv[1],true);\n");
4534     fprintf(fp, "  if (status) {\n");
4535     fprintf(fp, "    fprintf(stderr,%cBad readMps %%s%cn%c,argv[1]);\n",
4536             apo, backslash, apo);
4537     fprintf(fp, "    exit(1);\n");
4538     fprintf(fp, "  }\n\n");
4539     fprintf(fp, "  // Now do requested saves and modifications\n");
4540     fprintf(fp, "  ClpSimplex * clpModel = & model;\n");
4541     int wanted[9];
4542     memset(wanted, 0, sizeof(wanted));
4543     wanted[0] = wanted[3] = wanted[5] = wanted[8] = 1;
4544     if (type > 0)
4545          wanted[1] = wanted[6] = 1;
4546     if (type > 1)
4547          wanted[2] = wanted[4] = wanted[7] = 1;
4548     std::string header[9] = { "", "Save values", "Redundant save of default values", "Set changed values",
4549                               "Redundant set default values", "Solve", "Restore values", "Redundant restore values", "Add to model"
4550                             };
4551     for (int iType = 0; iType < 9; iType++) {
4552          if (!wanted[iType])
4553               continue;
4554          int n = 0;
4555          int iLine;
4556          for (iLine = 0; iLine < numberLines; iLine++) {
4557               if