source: stable/1.16/Clp/src/ClpSolver.cpp @ 2220

Last change on this file since 2220 was 2220, checked in by forrest, 5 years ago

minor tweak for decomposition

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