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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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