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

Last change on this file since 2449 was 2440, checked in by forrest, 8 months ago

need changes in svn as well as git

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