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

Last change on this file since 2254 was 2254, checked in by forrest, 2 years ago

segfault in primal (and out some messages)

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