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

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

allow AMPL in Clp

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