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

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

minor changes to implement Aboca

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