source: trunk/Clp/src/ClpMain.cpp @ 1975

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

allow custom print format

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