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

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

Allow a strategy for initial solve where code analyzes problem and guesses at parameters

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