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

Last change on this file since 2030 was 2030, checked in by forrest, 6 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

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