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

Last change on this file since 2007 was 2007, checked in by tkr, 6 years ago

Syncing with trunk at r2006

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