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

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

changes for miqp, parameters and presolve

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