source: stable/1.15/Clp/src/ClpMain.cpp @ 1949

Last change on this file since 1949 was 1949, checked in by forrest, 7 years ago

fix for postsolve should be in stable as well

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