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

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

stuff for vector matrix

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