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

Last change on this file since 1929 was 1929, checked in by stefan, 7 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

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