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

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

try to fix infeasibility ray,
changes as in stable (for presolve),
stuff for Cbc parameters

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