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

Last change on this file since 1424 was 1424, checked in by forrest, 12 years ago

fixes for various things

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 92.5 KB
Line 
1/* $Id: ClpMain.cpp 1424 2009-08-21 16:15:31Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4   
5#include "CoinPragma.hpp"
6
7#include <cassert>
8#include <cstdio>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12#include <string>
13#include <iostream>
14
15int boundary_sort=1000;
16int boundary_sort2=1000;
17int boundary_sort3=10000;
18
19#include "CoinPragma.hpp"
20#include "CoinHelperFunctions.hpp"
21#include "CoinSort.hpp"
22// History since 1.0 at end
23#include "CoinMpsIO.hpp"
24#include "CoinFileIO.hpp"
25
26#include "ClpFactorization.hpp"
27#include "CoinTime.hpp"
28#include "ClpSimplex.hpp"
29#include "ClpSimplexOther.hpp"
30#include "ClpSolve.hpp"
31#include "ClpPackedMatrix.hpp"
32#include "ClpPlusMinusOneMatrix.hpp"
33#include "ClpNetworkMatrix.hpp"
34#include "ClpDualRowSteepest.hpp"
35#include "ClpDualRowDantzig.hpp"
36#include "ClpLinearObjective.hpp"
37#include "ClpPrimalColumnSteepest.hpp"
38#include "ClpPrimalColumnDantzig.hpp"
39#include "ClpPresolve.hpp"
40#include "CbcOrClpParam.hpp"
41#include "CoinSignal.hpp"
42#ifdef DMALLOC
43#include "dmalloc.h"
44#endif
45#ifdef WSSMP_BARRIER
46#define FOREIGN_BARRIER
47#endif
48#ifdef UFL_BARRIER
49#define FOREIGN_BARRIER
50#endif
51#ifdef TAUCS_BARRIER
52#define FOREIGN_BARRIER
53#endif
54#ifdef MUMPS_BARRIER
55#define FOREIGN_BARRIER
56#endif
57
58static double totalTime=0.0;
59static bool maskMatches(const int * starts, char ** masks,
60                        std::string & check);
61static ClpSimplex * currentModel = NULL;
62
63extern "C" {
64   static void 
65#if defined(_MSC_VER)
66   __cdecl
67#endif // _MSC_VER
68   signal_handler(int /*whichSignal*/)
69   {
70      if (currentModel!=NULL) 
71         currentModel->setMaximumIterations(0); // stop at next iterations
72      return;
73   }
74}
75
76//#############################################################################
77
78#ifdef NDEBUG
79#undef NDEBUG
80#endif
81
82int mainTest (int argc, const char *argv[],int algorithm,
83              ClpSimplex empty, bool doPresolve,int switchOff,bool doVector);
84static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
85static void generateCode(const char * fileName,int type);
86// Returns next valid field
87int CbcOrClpRead_mode=1;
88FILE * CbcOrClpReadCommand=stdin;
89extern int CbcOrClpEnvironmentIndex;
90
91int
92#if defined(_MSC_VER)
93__cdecl
94#endif // _MSC_VER
95main (int argc, const char *argv[])
96{
97  // next {} is just to make sure all memory should be freed - for debug
98  {
99    double time1 = CoinCpuTime(),time2;
100    // Set up all non-standard stuff
101    //int numberModels=1;
102    ClpSimplex * models = new ClpSimplex[1];
103   
104    // default action on import
105    int allowImportErrors=0;
106    int keepImportNames=1;
107    int doIdiot=-1;
108    int outputFormat=2;
109    int slpValue=-1;
110    int cppValue=-1;
111    int printOptions=0;
112    int printMode=0;
113    int presolveOptions=0;
114    int doCrash=0;
115    int doVector=0;
116    int doSprint=-1;
117    // set reasonable defaults
118    int preSolve=5;
119    bool preSolveFile=false;
120    models->setPerturbation(50);
121    models->messageHandler()->setPrefix(false);
122    const char dirsep =  CoinFindDirSeparator();
123    std::string directory;
124    std::string dirSample;
125    std::string dirNetlib;
126    std::string dirMiplib;
127    if (dirsep == '/') {
128      directory = "./";
129      dirSample = "../../Data/Sample/";
130      dirNetlib = "../../Data/Netlib/";
131      dirMiplib = "../../Data/miplib3/";
132    } else {
133      directory = ".\\";
134      dirSample = "..\\..\\Data\\Sample\\";
135      dirNetlib = "..\\..\\Data\\Netlib\\";
136      dirMiplib = "..\\..\\Data\\miplib3\\";
137    }
138    std::string defaultDirectory = directory;
139    std::string importFile ="";
140    std::string exportFile ="default.mps";
141    std::string importBasisFile ="";
142    int basisHasValues=0;
143    int substitution=3;
144    int dualize=3;  // dualize if looks promising
145    std::string exportBasisFile ="default.bas";
146    std::string saveFile ="default.prob";
147    std::string restoreFile ="default.prob";
148    std::string solutionFile ="stdout";
149    std::string solutionSaveFile ="solution.file";
150    std::string printMask="";
151    CbcOrClpParam parameters[CBCMAXPARAMETERS];
152    int numberParameters ;
153    establishParams(numberParameters,parameters) ;
154    parameters[whichParam(BASISIN,numberParameters,parameters)].setStringValue(importBasisFile);
155    parameters[whichParam(BASISOUT,numberParameters,parameters)].setStringValue(exportBasisFile);
156    parameters[whichParam(PRINTMASK,numberParameters,parameters)].setStringValue(printMask);
157    parameters[whichParam(DIRECTORY,numberParameters,parameters)].setStringValue(directory);
158    parameters[whichParam(DIRSAMPLE,numberParameters,parameters)].setStringValue(dirSample);
159    parameters[whichParam(DIRNETLIB,numberParameters,parameters)].setStringValue(dirNetlib);
160    parameters[whichParam(DIRMIPLIB,numberParameters,parameters)].setStringValue(dirMiplib);
161    parameters[whichParam(DUALBOUND,numberParameters,parameters)].setDoubleValue(models->dualBound());
162    parameters[whichParam(DUALTOLERANCE,numberParameters,parameters)].setDoubleValue(models->dualTolerance());
163    parameters[whichParam(EXPORT,numberParameters,parameters)].setStringValue(exportFile);
164    parameters[whichParam(IDIOT,numberParameters,parameters)].setIntValue(doIdiot);
165    parameters[whichParam(IMPORT,numberParameters,parameters)].setStringValue(importFile);
166    parameters[whichParam(SOLVERLOGLEVEL,numberParameters,parameters)].setIntValue(models->logLevel());
167    parameters[whichParam(MAXFACTOR,numberParameters,parameters)].setIntValue(models->factorizationFrequency());
168    parameters[whichParam(MAXITERATION,numberParameters,parameters)].setIntValue(models->maximumIterations());
169    parameters[whichParam(OUTPUTFORMAT,numberParameters,parameters)].setIntValue(outputFormat);
170    parameters[whichParam(PRESOLVEPASS,numberParameters,parameters)].setIntValue(preSolve);
171    parameters[whichParam(PERTVALUE,numberParameters,parameters)].setIntValue(models->perturbation());
172    parameters[whichParam(PRIMALTOLERANCE,numberParameters,parameters)].setDoubleValue(models->primalTolerance());
173    parameters[whichParam(PRIMALWEIGHT,numberParameters,parameters)].setDoubleValue(models->infeasibilityCost());
174    parameters[whichParam(RESTORE,numberParameters,parameters)].setStringValue(restoreFile);
175    parameters[whichParam(SAVE,numberParameters,parameters)].setStringValue(saveFile);
176    parameters[whichParam(TIMELIMIT,numberParameters,parameters)].setDoubleValue(models->maximumSeconds());
177    parameters[whichParam(SOLUTION,numberParameters,parameters)].setStringValue(solutionFile);
178    parameters[whichParam(SAVESOL,numberParameters,parameters)].setStringValue(solutionSaveFile);
179    parameters[whichParam(SPRINT,numberParameters,parameters)].setIntValue(doSprint);
180    parameters[whichParam(SUBSTITUTION,numberParameters,parameters)].setIntValue(substitution);
181    parameters[whichParam(DUALIZE,numberParameters,parameters)].setIntValue(dualize);
182    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].setDoubleValue(1.0e-8);
183    int verbose=0;
184   
185    // total number of commands read
186    int numberGoodCommands=0;
187    bool * goodModels = new bool[1];
188
189    // Hidden stuff for barrier
190    int choleskyType = 0;
191    int gamma=0;
192    parameters[whichParam(BARRIERSCALE,numberParameters,parameters)].setCurrentOption(2);
193    int scaleBarrier=2;
194    int doKKT=0;
195    int crossover=2; // do crossover unless quadratic
196   
197    int iModel=0;
198    goodModels[0]=false;
199    //models[0].scaling(1);
200    //models[0].setDualBound(1.0e6);
201    //models[0].setDualTolerance(1.0e-7);
202    //ClpDualRowSteepest steep;
203    //models[0].setDualRowPivotAlgorithm(steep);
204    //models[0].setPrimalTolerance(1.0e-7);
205    //ClpPrimalColumnSteepest steepP;
206    //models[0].setPrimalColumnPivotAlgorithm(steepP);
207    std::string field;
208    std::cout<<"Coin LP version "<<CLPVERSION
209             <<", build "<<__DATE__<<std::endl;
210    // Print command line
211    if (argc>1) {
212      printf("command line - ");
213      for (int i=0;i<argc;i++)
214        printf("%s ",argv[i]);
215      printf("\n");
216    }
217   
218    while (1) {
219      // next command
220      field=CoinReadGetCommand(argc,argv);
221     
222      // exit if null or similar
223      if (!field.length()) {
224        if (numberGoodCommands==1&&goodModels[0]) {
225          // we just had file name - do dual or primal
226          field="either";
227        } else if (!numberGoodCommands) {
228          // let's give the sucker a hint
229          std::cout
230            <<"Clp takes input from arguments ( - switches to stdin)"
231            <<std::endl
232            <<"Enter ? for list of commands or help"<<std::endl;
233          field="-";
234        } else {
235          break;
236        }
237      }
238     
239      // see if ? at end
240      int numberQuery=0;
241      if (field!="?"&&field!="???") {
242        int length = field.length();
243        int i;
244        for (i=length-1;i>0;i--) {
245          if (field[i]=='?') 
246            numberQuery++;
247          else
248            break;
249        }
250        field=field.substr(0,length-numberQuery);
251      }
252      // find out if valid command
253      int iParam;
254      int numberMatches=0;
255      int firstMatch=-1;
256      for ( iParam=0; iParam<numberParameters; iParam++ ) {
257        int match = parameters[iParam].matches(field);
258        if (match==1) {
259          numberMatches = 1;
260          firstMatch=iParam;
261          break;
262        } else {
263          if (match&&firstMatch<0)
264            firstMatch=iParam;
265          numberMatches += match>>1;
266        }
267      }
268      if (iParam<numberParameters&&!numberQuery) {
269        // found
270        CbcOrClpParam found = parameters[iParam];
271        CbcOrClpParameterType type = found.type();
272        int valid;
273        numberGoodCommands++;
274        if (type==GENERALQUERY) {
275          std::cout<<"In argument list keywords have leading - "
276            ", -stdin or just - switches to stdin"<<std::endl;
277          std::cout<<"One command per line (and no -)"<<std::endl;
278          std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
279          std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
280          std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
281          std::cout<<"abcd value sets value"<<std::endl;
282          std::cout<<"Commands are:"<<std::endl;
283          int maxAcross=10;
284          bool evenHidden=false;
285          int printLevel = 
286            parameters[whichParam(ALLCOMMANDS,
287                                   numberParameters,parameters)].currentOptionAsInteger();
288          int convertP[]={2,1,0};
289          printLevel=convertP[printLevel];
290          if ((verbose&8)!=0) {
291            // even hidden
292            evenHidden = true;
293            verbose &= ~8;
294          }
295          if (verbose)
296            maxAcross=1;
297          int limits[]={1,101,201,301,401};
298          std::vector<std::string> types;
299          types.push_back("Double parameters:");
300          types.push_back("Int parameters:");
301          types.push_back("Keyword parameters:");
302          types.push_back("Actions or string parameters:");
303          int iType;
304          for (iType=0;iType<4;iType++) {
305            int across=0;
306            int lengthLine=0;
307            if ((verbose%4)!=0)
308              std::cout<<std::endl;
309            std::cout<<types[iType]<<std::endl;
310            if ((verbose&2)!=0)
311              std::cout<<std::endl;
312            for ( iParam=0; iParam<numberParameters; iParam++ ) {
313              int type = parameters[iParam].type();
314              //printf("%d type %d limits %d %d display %d\n",iParam,
315              //   type,limits[iType],limits[iType+1],parameters[iParam].displayThis());
316              if ((parameters[iParam].displayThis()>=printLevel||evenHidden)&&
317                  type>=limits[iType]
318                  &&type<limits[iType+1]) {
319                if (!across) {
320                  if ((verbose&2)!=0) 
321                    std::cout<<"Command ";
322                }
323                int length = parameters[iParam].lengthMatchName()+1;
324                if (lengthLine+length>80) {
325                  std::cout<<std::endl;
326                  across=0;
327                  lengthLine=0;
328                }
329                std::cout<<" "<<parameters[iParam].matchName();
330                lengthLine += length ;
331                across++;
332                if (across==maxAcross) {
333                  across=0;
334                  if (verbose) {
335                    // put out description as well
336                    if ((verbose&1)!=0) 
337                      std::cout<<parameters[iParam].shortHelp();
338                    std::cout<<std::endl;
339                    if ((verbose&2)!=0) {
340                      std::cout<<"---- description"<<std::endl;
341                      parameters[iParam].printLongHelp();
342                      std::cout<<"----"<<std::endl<<std::endl;
343                    }
344                  } else {
345                    std::cout<<std::endl;
346                  }
347                }
348              }
349            }
350            if (across)
351              std::cout<<std::endl;
352          }
353        } else if (type==FULLGENERALQUERY) {
354          std::cout<<"Full list of commands is:"<<std::endl;
355          int maxAcross=5;
356          int limits[]={1,101,201,301,401};
357          std::vector<std::string> types;
358          types.push_back("Double parameters:");
359          types.push_back("Int parameters:");
360          types.push_back("Keyword parameters and others:");
361          types.push_back("Actions:");
362          int iType;
363          for (iType=0;iType<4;iType++) {
364            int across=0;
365            std::cout<<types[iType]<<std::endl;
366            for ( iParam=0; iParam<numberParameters; iParam++ ) {
367              int type = parameters[iParam].type();
368              if (type>=limits[iType]
369                  &&type<limits[iType+1]) {
370                if (!across)
371                  std::cout<<"  ";
372                std::cout<<parameters[iParam].matchName()<<"  ";
373                across++;
374                if (across==maxAcross) {
375                  std::cout<<std::endl;
376                  across=0;
377                }
378              }
379            }
380            if (across)
381              std::cout<<std::endl;
382          }
383        } else if (type<101) {
384          // get next field as double
385          double value = CoinReadGetDoubleField(argc,argv,&valid);
386          if (!valid) {
387            parameters[iParam].setDoubleParameter(models+iModel,value);
388          } else if (valid==1) {
389            std::cout<<" is illegal for double parameter "<<parameters[iParam].name()<<" value remains "<<
390              parameters[iParam].doubleValue()<<std::endl;
391          } else {
392            std::cout<<parameters[iParam].name()<<" has value "<<
393              parameters[iParam].doubleValue()<<std::endl;
394          }
395        } else if (type<201) {
396          // get next field as int
397          int value = CoinReadGetIntField(argc,argv,&valid);
398          if (!valid) {
399            if (parameters[iParam].type()==PRESOLVEPASS)
400              preSolve = value;
401            else if (parameters[iParam].type()==IDIOT)
402              doIdiot = value;
403            else if (parameters[iParam].type()==SPRINT)
404              doSprint = value;
405            else if (parameters[iParam].type()==OUTPUTFORMAT)
406              outputFormat = value;
407            else if (parameters[iParam].type()==SLPVALUE)
408              slpValue = value;
409            else if (parameters[iParam].type()==CPP)
410              cppValue = value;
411            else if (parameters[iParam].type()==PRESOLVEOPTIONS)
412              presolveOptions = value;
413            else if (parameters[iParam].type()==PRINTOPTIONS)
414              printOptions = value;
415            else if (parameters[iParam].type()==SUBSTITUTION)
416              substitution = value;
417            else if (parameters[iParam].type()==DUALIZE)
418              dualize = value;
419            else if (parameters[iParam].type()==VERBOSE)
420              verbose = value;
421            parameters[iParam].setIntParameter(models+iModel,value);
422          } else if (valid==1) {
423            std::cout<<" is illegal for integer parameter "<<parameters[iParam].name()<<" value remains "<<
424              parameters[iParam].intValue()<<std::endl;
425          } else {
426            std::cout<<parameters[iParam].name()<<" has value "<<
427              parameters[iParam].intValue()<<std::endl;
428          }
429        } else if (type<301) {
430          // one of several strings
431          std::string value = CoinReadGetString(argc,argv);
432          int action = parameters[iParam].parameterOption(value);
433          if (action<0) {
434            if (value!="EOL") {
435              // no match
436              parameters[iParam].printOptions();
437            } else {
438              // print current value
439              std::cout<<parameters[iParam].name()<<" has value "<<
440                parameters[iParam].currentOption()<<std::endl;
441            }
442          } else {
443            parameters[iParam].setCurrentOption(action);
444            // for now hard wired
445            switch (type) {
446            case DIRECTION:
447              if (action==0)
448                models[iModel].setOptimizationDirection(1);
449              else if (action==1)
450                models[iModel].setOptimizationDirection(-1);
451              else
452                models[iModel].setOptimizationDirection(0);
453              break;
454            case DUALPIVOT:
455              if (action==0) {
456                ClpDualRowSteepest steep(3);
457                models[iModel].setDualRowPivotAlgorithm(steep);
458              } else if (action==1) {
459                //ClpDualRowDantzig dantzig;
460                ClpDualRowSteepest dantzig(5);
461                models[iModel].setDualRowPivotAlgorithm(dantzig);
462              } else if (action==2) {
463                // partial steep
464                ClpDualRowSteepest steep(2);
465                models[iModel].setDualRowPivotAlgorithm(steep);
466              } else {
467                ClpDualRowSteepest steep;
468                models[iModel].setDualRowPivotAlgorithm(steep);
469              }
470              break;
471            case PRIMALPIVOT:
472              if (action==0) {
473                ClpPrimalColumnSteepest steep(3);
474                models[iModel].setPrimalColumnPivotAlgorithm(steep);
475              } else if (action==1) {
476                ClpPrimalColumnSteepest steep(0);
477                models[iModel].setPrimalColumnPivotAlgorithm(steep);
478              } else if (action==2) {
479                ClpPrimalColumnDantzig dantzig;
480                models[iModel].setPrimalColumnPivotAlgorithm(dantzig);
481              } else if (action==3) {
482                ClpPrimalColumnSteepest steep(4);
483                models[iModel].setPrimalColumnPivotAlgorithm(steep);
484              } else if (action==4) {
485                ClpPrimalColumnSteepest steep(1);
486                models[iModel].setPrimalColumnPivotAlgorithm(steep);
487              } else if (action==5) {
488                ClpPrimalColumnSteepest steep(2);
489                models[iModel].setPrimalColumnPivotAlgorithm(steep);
490              } else if (action==6) {
491                ClpPrimalColumnSteepest steep(10);
492                models[iModel].setPrimalColumnPivotAlgorithm(steep);
493              }
494              break;
495            case SCALING:
496              models[iModel].scaling(action);
497              break;
498            case AUTOSCALE:
499              models[iModel].setAutomaticScaling(action!=0);
500              break;
501            case SPARSEFACTOR:
502              models[iModel].setSparseFactorization((1-action)!=0);
503              break;
504            case BIASLU:
505              models[iModel].factorization()->setBiasLU(action);
506              break;
507            case PERTURBATION:
508              if (action==0)
509                models[iModel].setPerturbation(50);
510              else
511                models[iModel].setPerturbation(100);
512              break;
513            case ERRORSALLOWED:
514              allowImportErrors = action;
515              break;
516            case INTPRINT:
517              printMode=action;
518              break;
519            case KEEPNAMES:
520              keepImportNames = 1-action;
521              break;
522            case PRESOLVE:
523              if (action==0)
524                preSolve = 5;
525              else if (action==1)
526                preSolve=0;
527              else if (action==2)
528                preSolve=10;
529              else
530                preSolveFile=true;
531              break;
532            case PFI:
533              models[iModel].factorization()->setForrestTomlin(action==0);
534              break;
535            case FACTORIZATION:
536              models[iModel].factorization()->forceOtherFactorization(action);
537              break;
538            case CRASH:
539              doCrash=action;
540              break;
541            case VECTOR:
542              doVector=action;
543              break;
544            case MESSAGES:
545              models[iModel].messageHandler()->setPrefix(action!=0);
546              break;
547            case CHOLESKY:
548              choleskyType = action;
549              break;
550            case GAMMA:
551              gamma=action;
552              break;
553            case BARRIERSCALE:
554              scaleBarrier=action;
555              break;
556            case KKT:
557              doKKT=action;
558              break;
559            case CROSSOVER:
560              crossover=action;
561              break;
562            default:
563              //abort();
564              break;
565            }
566          }
567        } else {
568          // action
569          if (type==EXIT)
570            break; // stop all
571          switch (type) {
572          case DUALSIMPLEX:
573          case PRIMALSIMPLEX:
574          case EITHERSIMPLEX:
575          case BARRIER:
576            // synonym for dual
577          case BAB:
578            if (goodModels[iModel]) {
579              double objScale = 
580                parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
581              if (objScale!=1.0) {
582                int iColumn;
583                int numberColumns=models[iModel].numberColumns();
584                double * dualColumnSolution = 
585                  models[iModel].dualColumnSolution();
586                ClpObjective * obj = models[iModel].objectiveAsObject();
587                assert(dynamic_cast<ClpLinearObjective *> (obj));
588                double offset;
589                double * objective = obj->gradient(NULL,NULL,offset,true);
590                for (iColumn=0;iColumn<numberColumns;iColumn++) {
591                  dualColumnSolution[iColumn] *= objScale;
592                  objective[iColumn] *= objScale;;
593                }
594                int iRow;
595                int numberRows=models[iModel].numberRows();
596                double * dualRowSolution = 
597                  models[iModel].dualRowSolution();
598                for (iRow=0;iRow<numberRows;iRow++) 
599                  dualRowSolution[iRow] *= objScale;
600                models[iModel].setObjectiveOffset(objScale*models[iModel].objectiveOffset());
601              }
602              ClpSolve::SolveType method;
603              ClpSolve::PresolveType presolveType;
604              ClpSimplex * model2 = models+iModel;
605              if (dualize) {
606                bool tryIt=true;
607                double fractionColumn=1.0;
608                double fractionRow=1.0;
609                if (dualize==3) {
610                  dualize=1;
611                  int numberColumns=model2->numberColumns();
612                  int numberRows=model2->numberRows();
613                  if (numberRows<50000||5*numberColumns>numberRows) {
614                    tryIt=false;
615                  } else {
616                    fractionColumn=0.1;
617                    fractionRow=0.1;
618                  }
619                }
620                if (tryIt) {
621                  model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(fractionRow,fractionColumn);
622                  if (model2) {
623                    printf("Dual of model has %d rows and %d columns\n",
624                           model2->numberRows(),model2->numberColumns());
625                    model2->setOptimizationDirection(1.0);
626                  } else {
627                    model2 = models+iModel;
628                    dualize=0;
629                  }
630                } else {
631                  dualize=0;
632                }
633              }
634              ClpSolve solveOptions;
635              solveOptions.setPresolveActions(presolveOptions);
636              solveOptions.setSubstitution(substitution);
637              if (preSolve!=5&&preSolve) {
638                presolveType=ClpSolve::presolveNumber;
639                if (preSolve<0) {
640                  preSolve = - preSolve;
641                  if (preSolve<=100) {
642                    presolveType=ClpSolve::presolveNumber;
643                    printf("Doing %d presolve passes - picking up non-costed slacks\n",
644                           preSolve);
645                    solveOptions.setDoSingletonColumn(true);
646                  } else {
647                    preSolve -=100;
648                    presolveType=ClpSolve::presolveNumberCost;
649                    printf("Doing %d presolve passes - picking up costed slacks\n",
650                           preSolve);
651                  }
652                } 
653              } else if (preSolve) {
654                presolveType=ClpSolve::presolveOn;
655              } else {
656                presolveType=ClpSolve::presolveOff;
657              }
658              solveOptions.setPresolveType(presolveType,preSolve);
659              if (type==DUALSIMPLEX||type==BAB) {
660                method=ClpSolve::useDual;
661              } else if (type==PRIMALSIMPLEX) {
662                method=ClpSolve::usePrimalorSprint;
663              } else if (type==EITHERSIMPLEX) {
664                method=ClpSolve::automatic;
665              } else {
666                method = ClpSolve::useBarrier;
667                if (crossover==1) {
668                  method=ClpSolve::useBarrierNoCross;
669                } else if (crossover==2) {
670                  ClpObjective * obj = models[iModel].objectiveAsObject();
671                  if (obj->type()>1) {
672                    method=ClpSolve::useBarrierNoCross;
673                    presolveType=ClpSolve::presolveOff;
674                    solveOptions.setPresolveType(presolveType,preSolve);
675                  } 
676                }
677              }
678              solveOptions.setSolveType(method);
679              if(preSolveFile)
680                presolveOptions |= 0x40000000;
681              solveOptions.setSpecialOption(4,presolveOptions);
682              solveOptions.setSpecialOption(5,printOptions&1);
683              if (doVector) {
684                ClpMatrixBase * matrix = models[iModel].clpMatrix();
685                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
686                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
687                  clpMatrix->makeSpecialColumnCopy();
688                }
689              }
690              if (method==ClpSolve::useDual) {
691                // dual
692                if (doCrash)
693                  solveOptions.setSpecialOption(0,1,doCrash); // crash
694                else if (doIdiot)
695                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
696              } else if (method==ClpSolve::usePrimalorSprint) {
697                // primal
698                // if slp turn everything off
699                if (slpValue>0) {
700                  doCrash=false;
701                  doSprint=0;
702                  doIdiot=-1;
703                  solveOptions.setSpecialOption(1,10,slpValue); // slp
704                  method=ClpSolve::usePrimal;
705                }
706                if (doCrash) {
707                  solveOptions.setSpecialOption(1,1,doCrash); // crash
708                } else if (doSprint>0) {
709                  // sprint overrides idiot
710                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
711                } else if (doIdiot>0) {
712                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
713                } else if (slpValue<=0) {
714                  if (doIdiot==0) {
715                    if (doSprint==0)
716                      solveOptions.setSpecialOption(1,4); // all slack
717                    else
718                      solveOptions.setSpecialOption(1,9); // all slack or sprint
719                  } else {
720                    if (doSprint==0)
721                      solveOptions.setSpecialOption(1,8); // all slack or idiot
722                    else
723                      solveOptions.setSpecialOption(1,7); // initiative
724                  }
725                }
726                if (basisHasValues==-1)
727                  solveOptions.setSpecialOption(1,11); // switch off values
728              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
729                int barrierOptions = choleskyType;
730                if (scaleBarrier) {
731                  if ((scaleBarrier&1)!=0) 
732                    barrierOptions |= 8;
733                  barrierOptions |= 2048*(scaleBarrier>>1);
734                }
735                if (doKKT)
736                  barrierOptions |= 16;
737                if (gamma)
738                  barrierOptions |= 32*gamma;
739                if (crossover==3) 
740                  barrierOptions |= 256; // try presolve in crossover
741                solveOptions.setSpecialOption(4,barrierOptions);
742              }
743              int status;
744              if (cppValue>=0) {
745                // generate code
746                FILE * fp = fopen("user_driver.cpp","w");
747                if (fp) {
748                  // generate enough to do solveOptions
749                  model2->generateCpp(fp);
750                  solveOptions.generateCpp(fp);
751                  fclose(fp);
752                  // now call generate code
753                  generateCode("user_driver.cpp",cppValue);
754                } else {
755                  std::cout<<"Unable to open file user_driver.cpp"<<std::endl;
756                }
757              }
758#ifdef CLP_MULTIPLE_FACTORIZATIONS   
759              int denseCode = parameters[whichParam(DENSE,numberParameters,parameters)].intValue();
760              model2->factorization()->setGoDenseThreshold(denseCode);
761              int smallCode = parameters[whichParam(SMALLFACT,numberParameters,parameters)].intValue();
762              model2->factorization()->setGoSmallThreshold(smallCode);
763              model2->factorization()->goDenseOrSmall(model2->numberRows());
764#endif
765              try {
766                status=model2->initialSolve(solveOptions);
767              }
768              catch (CoinError e) {
769                e.print();
770                status=-1;
771              }
772              if (dualize) {
773                int returnCode=static_cast<ClpSimplexOther *> (models+iModel)->restoreFromDual(model2);
774                if (model2->status()==3)
775                  returnCode=0;
776                delete model2;
777                if (returnCode&&dualize!=2) {
778                  currentModel = models+iModel;
779                  // register signal handler
780                  signal(SIGINT,signal_handler);
781                  models[iModel].primal(1);
782                  currentModel=NULL;
783                }
784              }
785              if (status>=0)
786                basisHasValues=1;
787            } else {
788              std::cout<<"** Current model not valid"<<std::endl;
789            }
790            break;
791          case STATISTICS:
792            if (goodModels[iModel]) {
793              // If presolve on look at presolved
794              bool deleteModel2=false;
795              ClpSimplex * model2 = models+iModel;
796              if (preSolve) {
797                ClpPresolve pinfo;
798                int presolveOptions2 = presolveOptions&~0x40000000;
799                if ((presolveOptions2&0xffff)!=0)
800                  pinfo.setPresolveActions(presolveOptions2);
801                pinfo.setSubstitution(substitution);
802                if ((printOptions&1)!=0)
803                  pinfo.statistics();
804                double presolveTolerance = 
805                  parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
806                model2 = 
807                  pinfo.presolvedModel(models[iModel],presolveTolerance,
808                                       true,preSolve);
809                if (model2) {
810                  printf("Statistics for presolved model\n");
811                  deleteModel2=true;
812                } else {
813                  printf("Presolved model looks infeasible - will use unpresolved\n");
814                  model2 = models+iModel;
815                }
816              } else {
817                printf("Statistics for unpresolved model\n");
818                model2 =  models+iModel;
819              }
820              statistics(models+iModel,model2);
821              if (deleteModel2)
822                delete model2;
823            } else {
824              std::cout<<"** Current model not valid"<<std::endl;
825            }
826            break;
827          case TIGHTEN:
828            if (goodModels[iModel]) {
829              int numberInfeasibilities = models[iModel].tightenPrimalBounds();
830              if (numberInfeasibilities)
831                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
832            } else {
833              std::cout<<"** Current model not valid"<<std::endl;
834            }
835            break;
836          case PLUSMINUS:
837            if (goodModels[iModel]) {
838              ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
839              ClpPackedMatrix* clpMatrix =
840                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
841              if (clpMatrix) {
842                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
843                if (newMatrix->getIndices()) {
844                  models[iModel].replaceMatrix(newMatrix);
845                  delete saveMatrix;
846                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
847                } else {
848                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
849                }
850              } else {
851                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
852              }
853            } else {
854              std::cout<<"** Current model not valid"<<std::endl;
855            }
856            break;
857          case NETWORK:
858            if (goodModels[iModel]) {
859              ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
860              ClpPackedMatrix* clpMatrix =
861                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
862              if (clpMatrix) {
863                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
864                if (newMatrix->getIndices()) {
865                  models[iModel].replaceMatrix(newMatrix);
866                  delete saveMatrix;
867                  std::cout<<"Matrix converted to network matrix"<<std::endl;
868                } else {
869                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
870                }
871              } else {
872                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
873              }
874            } else {
875              std::cout<<"** Current model not valid"<<std::endl;
876            }
877            break;
878          case IMPORT:
879            {
880              // get next field
881              field = CoinReadGetString(argc,argv);
882              if (field=="$") {
883                field = parameters[iParam].stringValue();
884              } else if (field=="EOL") {
885                parameters[iParam].printString();
886                break;
887              } else {
888                parameters[iParam].setStringValue(field);
889              }
890              std::string fileName;
891              bool canOpen=false;
892              // See if gmpl file
893              int gmpl=0;
894              std::string gmplData;
895              if (field=="-") {
896                // stdin
897                canOpen=true;
898                fileName = "-";
899              } else {
900                // See if .lp
901                {
902                  const char * c_name = field.c_str();
903                  int length = strlen(c_name);
904                  if (length>3&&!strncmp(c_name+length-3,".lp",3))
905                    gmpl=-1; // .lp
906                }
907                bool absolutePath;
908                if (dirsep=='/') {
909                  // non Windows (or cygwin)
910                  absolutePath=(field[0]=='/');
911                } else {
912                  //Windows (non cycgwin)
913                  absolutePath=(field[0]=='\\');
914                  // but allow for :
915                  if (strchr(field.c_str(),':'))
916                    absolutePath=true;
917                }
918                if (absolutePath) {
919                  fileName = field;
920                } else if (field[0]=='~') {
921                  char * environVar = getenv("HOME");
922                  if (environVar) {
923                    std::string home(environVar);
924                    field=field.erase(0,1);
925                    fileName = home+field;
926                  } else {
927                    fileName=field;
928                  }
929                } else {
930                  fileName = directory+field;
931                  // See if gmpl (model & data) - or even lp file
932                  int length = field.size();
933                  int percent = field.find('%');
934                  if (percent<length&&percent>0) {
935                    gmpl=1;
936                    fileName = directory+field.substr(0,percent);
937                    gmplData = directory+field.substr(percent+1);
938                    if (percent<length-1)
939                      gmpl=2; // two files
940                    printf("GMPL model file %s and data file %s\n",
941                           fileName.c_str(),gmplData.c_str());
942                  }
943                }
944                std::string name=fileName;
945                if (fileCoinReadable(name)) {
946                  // can open - lets go for it
947                  canOpen=true;
948                  if (gmpl==2) {
949                    FILE *fp;
950                    fp=fopen(gmplData.c_str(),"r");
951                    if (fp) {
952                      fclose(fp);
953                    } else {
954                      canOpen=false;
955                      std::cout<<"Unable to open file "<<gmplData<<std::endl;
956                    }
957                  }
958                } else {
959                  std::cout<<"Unable to open file "<<fileName<<std::endl;
960                }
961              }
962              if (canOpen) {
963                int status;
964                if (!gmpl)
965                  status =models[iModel].readMps(fileName.c_str(),
966                                                 keepImportNames!=0,
967                                                 allowImportErrors!=0);
968                else if (gmpl>0)
969                  status= models[iModel].readGMPL(fileName.c_str(),
970                                                  (gmpl==2) ? gmplData.c_str() : NULL,
971                                                  keepImportNames!=0);
972                else
973                  status= models[iModel].readLp(fileName.c_str(),1.0e-12);
974                if (!status||(status>0&&allowImportErrors)) {
975                  goodModels[iModel]=true;
976                  // sets to all slack (not necessary?)
977                  models[iModel].createStatus();
978                  time2 = CoinCpuTime();
979                  totalTime += time2-time1;
980                  time1=time2;
981                  // Go to canned file if just input file
982                  if (CbcOrClpRead_mode==2&&argc==2) {
983                    // only if ends .mps
984                    char * find = const_cast<char *>(strstr(fileName.c_str(),".mps"));
985                    if (find&&find[4]=='\0') {
986                      find[1]='p'; find[2]='a';find[3]='r';
987                      FILE *fp=fopen(fileName.c_str(),"r");
988                      if (fp) {
989                        CbcOrClpReadCommand=fp; // Read from that file
990                        CbcOrClpRead_mode=-1;
991                      }
992                    }
993                  }
994                } else {
995                  // errors
996                  std::cout<<"There were "<<status<<
997                    " errors on input"<<std::endl;
998                }
999              }
1000            }
1001            break;
1002          case EXPORT:
1003            if (goodModels[iModel]) {
1004              double objScale = 
1005                parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
1006              if (objScale!=1.0) {
1007                int iColumn;
1008                int numberColumns=models[iModel].numberColumns();
1009                double * dualColumnSolution = 
1010                  models[iModel].dualColumnSolution();
1011                ClpObjective * obj = models[iModel].objectiveAsObject();
1012                assert(dynamic_cast<ClpLinearObjective *> (obj));
1013                double offset;
1014                double * objective = obj->gradient(NULL,NULL,offset,true);
1015                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1016                  dualColumnSolution[iColumn] *= objScale;
1017                  objective[iColumn] *= objScale;;
1018                }
1019                int iRow;
1020                int numberRows=models[iModel].numberRows();
1021                double * dualRowSolution = 
1022                  models[iModel].dualRowSolution();
1023                for (iRow=0;iRow<numberRows;iRow++) 
1024                  dualRowSolution[iRow] *= objScale;
1025                models[iModel].setObjectiveOffset(objScale*models[iModel].objectiveOffset());
1026              }
1027              // get next field
1028              field = CoinReadGetString(argc,argv);
1029              if (field=="$") {
1030                field = parameters[iParam].stringValue();
1031              } else if (field=="EOL") {
1032                parameters[iParam].printString();
1033                break;
1034              } else {
1035                parameters[iParam].setStringValue(field);
1036              }
1037              std::string fileName;
1038              bool canOpen=false;
1039              if (field[0]=='/'||field[0]=='\\') {
1040                fileName = field;
1041              } else if (field[0]=='~') {
1042                char * environVar = getenv("HOME");
1043                if (environVar) {
1044                  std::string home(environVar);
1045                  field=field.erase(0,1);
1046                  fileName = home+field;
1047                } else {
1048                  fileName=field;
1049                }
1050              } else {
1051                fileName = directory+field;
1052              }
1053              FILE *fp=fopen(fileName.c_str(),"w");
1054              if (fp) {
1055                // can open - lets go for it
1056                fclose(fp);
1057                canOpen=true;
1058              } else {
1059                std::cout<<"Unable to open file "<<fileName<<std::endl;
1060              }
1061              if (canOpen) {
1062                // If presolve on then save presolved
1063                bool deleteModel2=false;
1064                ClpSimplex * model2 = models+iModel;
1065                if (dualize&&dualize<3) {
1066                  model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel();
1067                  printf("Dual of model has %d rows and %d columns\n",
1068                         model2->numberRows(),model2->numberColumns());
1069                  model2->setOptimizationDirection(1.0);
1070                  preSolve=0; // as picks up from model
1071                }
1072                if (preSolve) {
1073                  ClpPresolve pinfo;
1074                  int presolveOptions2 = presolveOptions&~0x40000000;
1075                  if ((presolveOptions2&0xffff)!=0)
1076                    pinfo.setPresolveActions(presolveOptions2);
1077                  pinfo.setSubstitution(substitution);
1078                  if ((printOptions&1)!=0)
1079                    pinfo.statistics();
1080                  double presolveTolerance = 
1081                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1082                  model2 = 
1083                    pinfo.presolvedModel(models[iModel],presolveTolerance,
1084                                         true,preSolve,false,false);
1085                  if (model2) {
1086                    printf("Saving presolved model on %s\n",
1087                           fileName.c_str());
1088                    deleteModel2=true;
1089                  } else {
1090                    printf("Presolved model looks infeasible - saving original on %s\n",
1091                           fileName.c_str());
1092                    deleteModel2=false;
1093                    model2 = models+iModel;
1094
1095                  }
1096                } else {
1097                  printf("Saving model on %s\n",
1098                           fileName.c_str());
1099                }
1100#if 0
1101                // Convert names
1102                int iRow;
1103                int numberRows=model2->numberRows();
1104                int iColumn;
1105                int numberColumns=model2->numberColumns();
1106
1107                char ** rowNames = NULL;
1108                char ** columnNames = NULL;
1109                if (model2->lengthNames()) {
1110                  rowNames = new char * [numberRows];
1111                  for (iRow=0;iRow<numberRows;iRow++) {
1112                    rowNames[iRow] =
1113                      CoinStrdup(model2->rowName(iRow).c_str());
1114#ifdef STRIPBLANKS
1115                    char * xx = rowNames[iRow];
1116                    int i;
1117                    int length = strlen(xx);
1118                    int n=0;
1119                    for (i=0;i<length;i++) {
1120                      if (xx[i]!=' ')
1121                        xx[n++]=xx[i];
1122                    }
1123                    xx[n]='\0';
1124#endif
1125                  }
1126                 
1127                  columnNames = new char * [numberColumns];
1128                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1129                    columnNames[iColumn] =
1130                      CoinStrdup(model2->columnName(iColumn).c_str());
1131#ifdef STRIPBLANKS
1132                    char * xx = columnNames[iColumn];
1133                    int i;
1134                    int length = strlen(xx);
1135                    int n=0;
1136                    for (i=0;i<length;i++) {
1137                      if (xx[i]!=' ')
1138                        xx[n++]=xx[i];
1139                    }
1140                    xx[n]='\0';
1141#endif
1142                  }
1143                }
1144                CoinMpsIO writer;
1145                writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
1146                                  model2->getColLower(), model2->getColUpper(),
1147                                  model2->getObjCoefficients(),
1148                                  (const char*) 0 /*integrality*/,
1149                                  model2->getRowLower(), model2->getRowUpper(),
1150                                  columnNames, rowNames);
1151                // Pass in array saying if each variable integer
1152                writer.copyInIntegerInformation(model2->integerInformation());
1153                writer.setObjectiveOffset(model2->objectiveOffset());
1154                writer.writeMps(fileName.c_str(),0,1,1);
1155                if (rowNames) {
1156                  for (iRow=0;iRow<numberRows;iRow++) {
1157                    free(rowNames[iRow]);
1158                  }
1159                  delete [] rowNames;
1160                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1161                    free(columnNames[iColumn]);
1162                  }
1163                  delete [] columnNames;
1164                }
1165#else
1166                model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
1167#endif
1168                if (deleteModel2)
1169                  delete model2;
1170                time2 = CoinCpuTime();
1171                totalTime += time2-time1;
1172                time1=time2;
1173              }
1174            } else {
1175              std::cout<<"** Current model not valid"<<std::endl;
1176            }
1177            break;
1178          case BASISIN:
1179            if (goodModels[iModel]) {
1180              // get next field
1181              field = CoinReadGetString(argc,argv);
1182              if (field=="$") {
1183                field = parameters[iParam].stringValue();
1184              } else if (field=="EOL") {
1185                parameters[iParam].printString();
1186                break;
1187              } else {
1188                parameters[iParam].setStringValue(field);
1189              }
1190              std::string fileName;
1191              bool canOpen=false;
1192              if (field=="-") {
1193                // stdin
1194                canOpen=true;
1195                fileName = "-";
1196              } else {
1197                if (field[0]=='/'||field[0]=='\\') {
1198                  fileName = field;
1199                } else if (field[0]=='~') {
1200                  char * environVar = getenv("HOME");
1201                  if (environVar) {
1202                    std::string home(environVar);
1203                    field=field.erase(0,1);
1204                    fileName = home+field;
1205                  } else {
1206                    fileName=field;
1207                  }
1208                } else {
1209                  fileName = directory+field;
1210                }
1211                FILE *fp=fopen(fileName.c_str(),"r");
1212                if (fp) {
1213                  // can open - lets go for it
1214                  fclose(fp);
1215                  canOpen=true;
1216                } else {
1217                  std::cout<<"Unable to open file "<<fileName<<std::endl;
1218                }
1219              }
1220              if (canOpen) {
1221                int values = models[iModel].readBasis(fileName.c_str());
1222                if (values==0)
1223                  basisHasValues=-1;
1224                else
1225                  basisHasValues=1;
1226              }
1227            } else {
1228              std::cout<<"** Current model not valid"<<std::endl;
1229            }
1230            break;
1231          case PRINTMASK:
1232            // get next field
1233            {
1234              std::string name = CoinReadGetString(argc,argv);
1235              if (name!="EOL") {
1236                parameters[iParam].setStringValue(name);
1237                printMask = name;
1238              } else {
1239                parameters[iParam].printString();
1240              }
1241            }
1242            break;
1243          case BASISOUT:
1244            if (goodModels[iModel]) {
1245              // get next field
1246              field = CoinReadGetString(argc,argv);
1247              if (field=="$") {
1248                field = parameters[iParam].stringValue();
1249              } else if (field=="EOL") {
1250                parameters[iParam].printString();
1251                break;
1252              } else {
1253                parameters[iParam].setStringValue(field);
1254              }
1255              std::string fileName;
1256              bool canOpen=false;
1257              if (field[0]=='/'||field[0]=='\\') {
1258                fileName = field;
1259              } else if (field[0]=='~') {
1260                char * environVar = getenv("HOME");
1261                if (environVar) {
1262                  std::string home(environVar);
1263                  field=field.erase(0,1);
1264                  fileName = home+field;
1265                } else {
1266                  fileName=field;
1267                }
1268              } else {
1269                fileName = directory+field;
1270              }
1271              FILE *fp=fopen(fileName.c_str(),"w");
1272              if (fp) {
1273                // can open - lets go for it
1274                fclose(fp);
1275                canOpen=true;
1276              } else {
1277                std::cout<<"Unable to open file "<<fileName<<std::endl;
1278              }
1279              if (canOpen) {
1280                ClpSimplex * model2 = models+iModel;
1281                model2->writeBasis(fileName.c_str(),outputFormat>1,outputFormat-2);
1282                time2 = CoinCpuTime();
1283                totalTime += time2-time1;
1284                time1=time2;
1285              }
1286            } else {
1287              std::cout<<"** Current model not valid"<<std::endl;
1288            }
1289            break;
1290          case SAVE:
1291            {
1292              // get next field
1293              field = CoinReadGetString(argc,argv);
1294              if (field=="$") {
1295                field = parameters[iParam].stringValue();
1296              } else if (field=="EOL") {
1297                parameters[iParam].printString();
1298                break;
1299              } else {
1300                parameters[iParam].setStringValue(field);
1301              }
1302              std::string fileName;
1303              bool canOpen=false;
1304              if (field[0]=='/'||field[0]=='\\') {
1305                fileName = field;
1306              } else if (field[0]=='~') {
1307                char * environVar = getenv("HOME");
1308                if (environVar) {
1309                  std::string home(environVar);
1310                  field=field.erase(0,1);
1311                  fileName = home+field;
1312                } else {
1313                  fileName=field;
1314                }
1315              } else {
1316                fileName = directory+field;
1317              }
1318              FILE *fp=fopen(fileName.c_str(),"wb");
1319              if (fp) {
1320                // can open - lets go for it
1321                fclose(fp);
1322                canOpen=true;
1323              } else {
1324                std::cout<<"Unable to open file "<<fileName<<std::endl;
1325              }
1326              if (canOpen) {
1327                int status;
1328                // If presolve on then save presolved
1329                bool deleteModel2=false;
1330                ClpSimplex * model2 = models+iModel;
1331                if (preSolve) {
1332                  ClpPresolve pinfo;
1333                  double presolveTolerance = 
1334                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1335                  model2 = 
1336                    pinfo.presolvedModel(models[iModel],presolveTolerance,
1337                                         false,preSolve);
1338                  if (model2) {
1339                    printf("Saving presolved model on %s\n",
1340                           fileName.c_str());
1341                    deleteModel2=true;
1342                  } else {
1343                    printf("Presolved model looks infeasible - saving original on %s\n",
1344                           fileName.c_str());
1345                    deleteModel2=false;
1346                    model2 = models+iModel;
1347
1348                  }
1349                } else {
1350                  printf("Saving model on %s\n",
1351                           fileName.c_str());
1352                }
1353                status =model2->saveModel(fileName.c_str());
1354                if (deleteModel2)
1355                  delete model2;
1356                if (!status) {
1357                  goodModels[iModel]=true;
1358                  time2 = CoinCpuTime();
1359                  totalTime += time2-time1;
1360                  time1=time2;
1361                } else {
1362                  // errors
1363                  std::cout<<"There were errors on output"<<std::endl;
1364                }
1365              }
1366            }
1367            break;
1368          case RESTORE:
1369            {
1370              // get next field
1371              field = CoinReadGetString(argc,argv);
1372              if (field=="$") {
1373                field = parameters[iParam].stringValue();
1374              } else if (field=="EOL") {
1375                parameters[iParam].printString();
1376                break;
1377              } else {
1378                parameters[iParam].setStringValue(field);
1379              }
1380              std::string fileName;
1381              bool canOpen=false;
1382              if (field[0]=='/'||field[0]=='\\') {
1383                fileName = field;
1384              } else if (field[0]=='~') {
1385                char * environVar = getenv("HOME");
1386                if (environVar) {
1387                  std::string home(environVar);
1388                  field=field.erase(0,1);
1389                  fileName = home+field;
1390                } else {
1391                  fileName=field;
1392                }
1393              } else {
1394                fileName = directory+field;
1395              }
1396              FILE *fp=fopen(fileName.c_str(),"rb");
1397              if (fp) {
1398                // can open - lets go for it
1399                fclose(fp);
1400                canOpen=true;
1401              } else {
1402                std::cout<<"Unable to open file "<<fileName<<std::endl;
1403              }
1404              if (canOpen) {
1405                int status =models[iModel].restoreModel(fileName.c_str());
1406                if (!status) {
1407                  goodModels[iModel]=true;
1408                  time2 = CoinCpuTime();
1409                  totalTime += time2-time1;
1410                  time1=time2;
1411                } else {
1412                  // errors
1413                  std::cout<<"There were errors on input"<<std::endl;
1414                }
1415              }
1416            }
1417            break;
1418          case MAXIMIZE:
1419            models[iModel].setOptimizationDirection(-1);
1420            break;
1421          case MINIMIZE:
1422            models[iModel].setOptimizationDirection(1);
1423            break;
1424          case ALLSLACK:
1425            models[iModel].allSlackBasis(true);
1426            break;
1427          case REVERSE:
1428            if (goodModels[iModel]) {
1429              int iColumn;
1430              int numberColumns=models[iModel].numberColumns();
1431              double * dualColumnSolution = 
1432                models[iModel].dualColumnSolution();
1433              ClpObjective * obj = models[iModel].objectiveAsObject();
1434              assert(dynamic_cast<ClpLinearObjective *> (obj));
1435              double offset;
1436              double * objective = obj->gradient(NULL,NULL,offset,true);
1437              for (iColumn=0;iColumn<numberColumns;iColumn++) {
1438                dualColumnSolution[iColumn] = -dualColumnSolution[iColumn];
1439                objective[iColumn] = -objective[iColumn];
1440              }
1441              int iRow;
1442              int numberRows=models[iModel].numberRows();
1443              double * dualRowSolution = 
1444                models[iModel].dualRowSolution();
1445              for (iRow=0;iRow<numberRows;iRow++) {
1446                dualRowSolution[iRow] = -dualRowSolution[iRow];
1447              }
1448              models[iModel].setObjectiveOffset(-models[iModel].objectiveOffset());
1449            }
1450            break;
1451          case DIRECTORY:
1452            {
1453              std::string name = CoinReadGetString(argc,argv);
1454              if (name!="EOL") {
1455                int length=name.length();
1456                if (name[length-1]==dirsep) {
1457                  directory = name;
1458                } else {
1459                  directory = name+dirsep;
1460                }
1461                parameters[iParam].setStringValue(directory);
1462              } else {
1463                parameters[iParam].printString();
1464              }
1465            }
1466            break;
1467          case DIRSAMPLE:
1468            {
1469              std::string name = CoinReadGetString(argc,argv);
1470              if (name!="EOL") {
1471                int length=name.length();
1472                if (name[length-1]==dirsep) {
1473                  dirSample = name;
1474                } else {
1475                  dirSample = name+dirsep;
1476                }
1477                parameters[iParam].setStringValue(dirSample);
1478              } else {
1479                parameters[iParam].printString();
1480              }
1481            }
1482            break;
1483          case DIRNETLIB:
1484            {
1485              std::string name = CoinReadGetString(argc,argv);
1486              if (name!="EOL") {
1487                int length=name.length();
1488                if (name[length-1]==dirsep) {
1489                  dirNetlib = name;
1490                } else {
1491                  dirNetlib = name+dirsep;
1492                }
1493                parameters[iParam].setStringValue(dirNetlib);
1494              } else {
1495                parameters[iParam].printString();
1496              }
1497            }
1498            break;
1499          case DIRMIPLIB:
1500            {
1501              std::string name = CoinReadGetString(argc,argv);
1502              if (name!="EOL") {
1503                int length=name.length();
1504                if (name[length-1]==dirsep) {
1505                  dirMiplib = name;
1506                } else {
1507                  dirMiplib = name+dirsep;
1508                }
1509                parameters[iParam].setStringValue(dirMiplib);
1510              } else {
1511                parameters[iParam].printString();
1512              }
1513            }
1514            break;
1515          case STDIN:
1516            CbcOrClpRead_mode=-1;
1517            break;
1518          case NETLIB_DUAL:
1519          case NETLIB_EITHER:
1520          case NETLIB_BARRIER:
1521          case NETLIB_PRIMAL:
1522          case NETLIB_TUNE:
1523            {
1524              // create fields for unitTest
1525              const char * fields[4];
1526              int nFields=4;
1527              fields[0]="fake main from unitTest";
1528              std::string mpsfield = "-dirSample=";
1529              mpsfield += dirSample.c_str();
1530              fields[1]=mpsfield.c_str();
1531              std::string netfield = "-dirNetlib=";
1532              netfield += dirNetlib.c_str();
1533              fields[2]=netfield.c_str();
1534              fields[3]="-netlib";
1535              int algorithm;
1536              if (type==NETLIB_DUAL) {
1537                std::cerr<<"Doing netlib with dual algorithm"<<std::endl;
1538                algorithm =0;
1539              } else if (type==NETLIB_BARRIER) {
1540                std::cerr<<"Doing netlib with barrier algorithm"<<std::endl;
1541                algorithm =2;
1542              } else if (type==NETLIB_EITHER) {
1543                std::cerr<<"Doing netlib with dual or primal algorithm"<<std::endl;
1544                algorithm =3;
1545              } else if (type==NETLIB_TUNE) {
1546                std::cerr<<"Doing netlib with best algorithm!"<<std::endl;
1547                algorithm =5;
1548                // uncomment next to get active tuning
1549                // algorithm=6;
1550              } else {
1551                std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
1552                algorithm=1;
1553              }
1554              int specialOptions = models[iModel].specialOptions();
1555              models[iModel].setSpecialOptions(0);
1556              mainTest(nFields,fields,algorithm,models[iModel],
1557                       (preSolve!=0),specialOptions,doVector!=0);
1558            }
1559            break;
1560          case UNITTEST:
1561            {
1562              // create fields for unitTest
1563              const char * fields[2];
1564              int nFields=2;
1565              fields[0]="fake main from unitTest";
1566              std::string dirfield = "-dirSample=";
1567              dirfield += dirSample.c_str();
1568              fields[1]=dirfield.c_str();
1569              int specialOptions = models[iModel].specialOptions();
1570              models[iModel].setSpecialOptions(0);
1571              int algorithm=-1;
1572              if (models[iModel].numberRows())
1573                algorithm=7;
1574              mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions,doVector!=0);
1575            }
1576            break;
1577          case FAKEBOUND:
1578            if (goodModels[iModel]) {
1579              // get bound
1580              double value = CoinReadGetDoubleField(argc,argv,&valid);
1581              if (!valid) {
1582                std::cout<<"Setting "<<parameters[iParam].name()<<
1583                  " to DEBUG "<<value<<std::endl;
1584                int iRow;
1585                int numberRows=models[iModel].numberRows();
1586                double * rowLower = models[iModel].rowLower();
1587                double * rowUpper = models[iModel].rowUpper();
1588                for (iRow=0;iRow<numberRows;iRow++) {
1589                  // leave free ones for now
1590                  if (rowLower[iRow]>-1.0e20||rowUpper[iRow]<1.0e20) {
1591                    rowLower[iRow]=CoinMax(rowLower[iRow],-value);
1592                    rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
1593                  }
1594                }
1595                int iColumn;
1596                int numberColumns=models[iModel].numberColumns();
1597                double * columnLower = models[iModel].columnLower();
1598                double * columnUpper = models[iModel].columnUpper();
1599                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1600                  // leave free ones for now
1601                  if (columnLower[iColumn]>-1.0e20||
1602                      columnUpper[iColumn]<1.0e20) {
1603                    columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
1604                    columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
1605                  }
1606                }
1607              } else if (valid==1) {
1608                abort();
1609              } else {
1610                std::cout<<"enter value for "<<parameters[iParam].name()<<
1611                  std::endl;
1612              }
1613            }
1614            break;
1615          case REALLY_SCALE:
1616            if (goodModels[iModel]) {
1617              ClpSimplex newModel(models[iModel],
1618                                  models[iModel].scalingFlag());
1619              printf("model really really scaled\n");
1620              models[iModel]=newModel;
1621            }
1622            break;
1623          case USERCLP:
1624            // Replace the sample code by whatever you want
1625            if (goodModels[iModel]) {
1626              ClpSimplex * thisModel = &models[iModel];
1627              printf("Dummy user code - model has %d rows and %d columns\n",
1628                     thisModel->numberRows(),thisModel->numberColumns());
1629            }
1630            break;
1631          case HELP:
1632            std::cout<<"Coin LP version "<<CLPVERSION
1633                     <<", build "<<__DATE__<<std::endl;
1634            std::cout<<"Non default values:-"<<std::endl;
1635            std::cout<<"Perturbation "<<models[0].perturbation()<<" (default 100)"
1636                     <<std::endl;
1637            CoinReadPrintit(
1638                    "Presolve being done with 5 passes\n\
1639Dual steepest edge steep/partial on matrix shape and factorization density\n\
1640Clpnnnn taken out of messages\n\
1641If Factorization frequency default then done on size of matrix\n\n\
1642(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
1643You can switch to interactive mode at any time so\n\
1644clp watson.mps -scaling off -primalsimplex\nis the same as\n\
1645clp watson.mps -\nscaling off\nprimalsimplex"
1646                    );
1647            break;
1648          case SOLUTION:
1649            if (goodModels[iModel]) {
1650              // get next field
1651              field = CoinReadGetString(argc,argv);
1652              if (field=="$") {
1653                field = parameters[iParam].stringValue();
1654              } else if (field=="EOL") {
1655                parameters[iParam].printString();
1656                break;
1657              } else {
1658                parameters[iParam].setStringValue(field);
1659              }
1660              std::string fileName;
1661              FILE *fp=NULL;
1662              if (field=="-"||field=="EOL"||field=="stdout") {
1663                // stdout
1664                fp=stdout;
1665                fprintf(fp,"\n");
1666              } else if (field=="stderr") {
1667                // stderr
1668                fp=stderr;
1669                fprintf(fp,"\n");
1670              } else {
1671                if (field[0]=='/'||field[0]=='\\') {
1672                  fileName = field;
1673                } else if (field[0]=='~') {
1674                  char * environVar = getenv("HOME");
1675                  if (environVar) {
1676                    std::string home(environVar);
1677                    field=field.erase(0,1);
1678                    fileName = home+field;
1679                  } else {
1680                    fileName=field;
1681                  }
1682                } else {
1683                  fileName = directory+field;
1684                }
1685                fp=fopen(fileName.c_str(),"w");
1686              }
1687              if (fp) {
1688                // Write solution header (suggested by Luigi Poderico)
1689                double objValue = models[iModel].getObjValue()*models[iModel].getObjSense();
1690                int iStat = models[iModel].status();
1691                if (iStat==0) {
1692                  fprintf(fp, "optimal\n" );
1693                } else if (iStat==1) {
1694                  // infeasible
1695                  fprintf(fp, "infeasible\n" );
1696                } else if (iStat==2) {
1697                  // unbounded
1698                  fprintf(fp, "unbounded\n" );
1699                } else if (iStat==3) {
1700                  fprintf(fp, "stopped on iterations or time\n" );
1701                } else if (iStat==4) {
1702                  fprintf(fp, "stopped on difficulties\n" );
1703                } else if (iStat==5) {
1704                  fprintf(fp, "stopped on ctrl-c\n" );
1705                } else {
1706                  fprintf(fp, "status unknown\n" );
1707                }
1708                fprintf(fp, "Objective value %15.8g\n", objValue);
1709                // make fancy later on
1710                int iRow;
1711                int numberRows=models[iModel].numberRows();
1712                int lengthName = models[iModel].lengthNames(); // 0 if no names
1713                int lengthPrint = CoinMax(lengthName,8);
1714                // in general I don't want to pass around massive
1715                // amounts of data but seems simpler here
1716                std::vector<std::string> rowNames =
1717                  *(models[iModel].rowNames());
1718                std::vector<std::string> columnNames =
1719                  *(models[iModel].columnNames());
1720
1721                double * dualRowSolution = models[iModel].dualRowSolution();
1722                double * primalRowSolution = 
1723                  models[iModel].primalRowSolution();
1724                double * rowLower = models[iModel].rowLower();
1725                double * rowUpper = models[iModel].rowUpper();
1726                double primalTolerance = models[iModel].primalTolerance();
1727                bool doMask = (printMask!=""&&lengthName);
1728                int * maskStarts=NULL;
1729                int maxMasks=0;
1730                char ** masks =NULL;
1731                if (doMask) {
1732                  int nAst =0;
1733                  const char * pMask2 = printMask.c_str();
1734                  char pMask[100];
1735                  int iChar;
1736                  int lengthMask = strlen(pMask2);
1737                  assert (lengthMask<100);
1738                  if (*pMask2=='"') {
1739                    if (pMask2[lengthMask-1]!='"') {
1740                      printf("mismatched \" in mask %s\n",pMask2);
1741                      break;
1742                    } else {
1743                      strcpy(pMask,pMask2+1);
1744                      *strchr(pMask,'"')='\0';
1745                    }
1746                  } else if (*pMask2=='\'') {
1747                    if (pMask2[lengthMask-1]!='\'') {
1748                      printf("mismatched ' in mask %s\n",pMask2);
1749                      break;
1750                    } else {
1751                      strcpy(pMask,pMask2+1);
1752                      *strchr(pMask,'\'')='\0';
1753                    }
1754                  } else {
1755                    strcpy(pMask,pMask2);
1756                  }
1757                  if (lengthMask>lengthName) {
1758                    printf("mask %s too long - skipping\n",pMask);
1759                    break;
1760                  }
1761                  maxMasks = 1;
1762                  for (iChar=0;iChar<lengthMask;iChar++) {
1763                    if (pMask[iChar]=='*') {
1764                      nAst++;
1765                      maxMasks *= (lengthName+1);
1766                    }
1767                  }
1768                  int nEntries = 1;
1769                  maskStarts = new int[lengthName+2];
1770                  masks = new char * [maxMasks];
1771                  char ** newMasks = new char * [maxMasks];
1772                  int i;
1773                  for (i=0;i<maxMasks;i++) {
1774                    masks[i] = new char[lengthName+1];
1775                    newMasks[i] = new char[lengthName+1];
1776                  }
1777                  strcpy(masks[0],pMask);
1778                  for (int iAst=0;iAst<nAst;iAst++) {
1779                    int nOldEntries = nEntries;
1780                    nEntries=0;
1781                    for (int iEntry = 0;iEntry<nOldEntries;iEntry++) {
1782                      char * oldMask = masks[iEntry];
1783                      char * ast = strchr(oldMask,'*');
1784                      assert (ast);
1785                      int length = strlen(oldMask)-1;
1786                      int nBefore = ast-oldMask;
1787                      int nAfter = length-nBefore;
1788                      // and add null
1789                      nAfter++;
1790                      for (int i=0;i<=lengthName-length;i++) {
1791                        char * maskOut = newMasks[nEntries];
1792   CoinMemcpyN(oldMask,nBefore,maskOut);
1793                        for (int k=0;k<i;k++) 
1794                          maskOut[k+nBefore]='?';
1795   CoinMemcpyN(ast+1,nAfter,maskOut+nBefore+i);
1796                        nEntries++;
1797                        assert (nEntries<=maxMasks);
1798                      }
1799                    }
1800                    char ** temp = masks;
1801                    masks = newMasks;
1802                    newMasks = temp;
1803                  }
1804                  // Now extend and sort
1805                  int * sort = new int[nEntries];
1806                  for (i=0;i<nEntries;i++) {
1807                    char * maskThis = masks[i];
1808                    int length = strlen(maskThis);
1809                    while (maskThis[length-1]==' ')
1810                      length--;
1811                    maskThis[length]='\0';
1812                    sort[i]=length;
1813                  }
1814                  CoinSort_2(sort,sort+nEntries,masks);
1815                  int lastLength=-1;
1816                  for (i=0;i<nEntries;i++) {
1817                    int length = sort[i];
1818                    while (length>lastLength) 
1819                      maskStarts[++lastLength] = i;
1820                  }
1821                  maskStarts[++lastLength]=nEntries;
1822                  delete [] sort;
1823                  for (i=0;i<maxMasks;i++)
1824                    delete [] newMasks[i];
1825                  delete [] newMasks;
1826                }
1827                if (printMode>2) {
1828                  for (iRow=0;iRow<numberRows;iRow++) {
1829                    int type=printMode-3;
1830                    if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
1831                        primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
1832                      fprintf(fp,"** ");
1833                      type=2;
1834                    } else if (fabs(primalRowSolution[iRow])>1.0e-8) {
1835                      type=1;
1836                    } else if (numberRows<50) {
1837                      type=3;
1838                    } 
1839                    if (doMask&&!maskMatches(maskStarts,masks,rowNames[iRow]))
1840                      type=0;
1841                    if (type) {
1842                      fprintf(fp,"%7d ",iRow);
1843                      if (lengthName) {
1844                        const char * name = rowNames[iRow].c_str();
1845                        int n=strlen(name);
1846                        int i;
1847                        for (i=0;i<n;i++)
1848                          fprintf(fp,"%c",name[i]);
1849                        for (;i<lengthPrint;i++)
1850                          fprintf(fp," ");
1851                      }
1852                      fprintf(fp," %15.8g        %15.8g\n",primalRowSolution[iRow],
1853                              dualRowSolution[iRow]);
1854                    }
1855                  }
1856                }
1857                int iColumn;
1858                int numberColumns=models[iModel].numberColumns();
1859                double * dualColumnSolution = 
1860  models[iModel].dualColumnSolution();
1861                double * primalColumnSolution = 
1862  models[iModel].primalColumnSolution();
1863                double * columnLower = models[iModel].columnLower();
1864                double * columnUpper = models[iModel].columnUpper();
1865                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1866                  int type=(printMode>3) ? 1 : 0;
1867                  if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
1868                      primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
1869                    fprintf(fp,"** ");
1870                    type=2;
1871                  } else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
1872                    type=1;
1873                  } else if (numberColumns<50) {
1874                    type=3;
1875                  }
1876                  if (doMask&&!maskMatches(maskStarts,masks,
1877                                           columnNames[iColumn]))
1878                    type =0;
1879                  if (type) {
1880                    fprintf(fp,"%7d ",iColumn);
1881                    if (lengthName) {
1882                      const char * name = columnNames[iColumn].c_str();
1883                      int n=strlen(name);
1884                      int i;
1885                      for (i=0;i<n;i++)
1886                        fprintf(fp,"%c",name[i]);
1887                      for (;i<lengthPrint;i++)
1888                        fprintf(fp," ");
1889                    }
1890                    fprintf(fp," %15.8g        %15.8g\n",
1891                            primalColumnSolution[iColumn],
1892                            dualColumnSolution[iColumn]);
1893                  }
1894                }
1895                if (fp!=stdout)
1896                  fclose(fp);
1897                if (masks) {
1898                  delete [] maskStarts;
1899                  for (int i=0;i<maxMasks;i++)
1900                    delete [] masks[i];
1901                  delete [] masks;
1902                }
1903              } else {
1904                std::cout<<"Unable to open file "<<fileName<<std::endl;
1905              }
1906            } else {
1907              std::cout<<"** Current model not valid"<<std::endl;
1908             
1909            }
1910         
1911            break;
1912          case SAVESOL:
1913            if (goodModels[iModel]) {
1914              // get next field
1915              field = CoinReadGetString(argc,argv);
1916              if (field=="$") {
1917                field = parameters[iParam].stringValue();
1918              } else if (field=="EOL") {
1919                parameters[iParam].printString();
1920                break;
1921              } else {
1922                parameters[iParam].setStringValue(field);
1923              }
1924              std::string fileName;
1925              if (field[0]=='/'||field[0]=='\\') {
1926                fileName = field;
1927              } else if (field[0]=='~') {
1928                char * environVar = getenv("HOME");
1929                if (environVar) {
1930                  std::string home(environVar);
1931                  field=field.erase(0,1);
1932                  fileName = home+field;
1933                } else {
1934                  fileName=field;
1935                }
1936              } else {
1937                fileName = directory+field;
1938              }
1939              saveSolution(models+iModel,fileName);
1940            } else {
1941              std::cout<<"** Current model not valid"<<std::endl;
1942             
1943            }
1944            break;
1945          case ENVIRONMENT:
1946            CbcOrClpEnvironmentIndex=0;
1947            break;
1948          default:
1949            abort();
1950          }
1951        } 
1952      } else if (!numberMatches) {
1953        std::cout<<"No match for "<<field<<" - ? for list of commands"
1954                 <<std::endl;
1955      } else if (numberMatches==1) {
1956        if (!numberQuery) {
1957          std::cout<<"Short match for "<<field<<" - completion: ";
1958          std::cout<<parameters[firstMatch].matchName()<<std::endl;
1959        } else if (numberQuery) {
1960          std::cout<<parameters[firstMatch].matchName()<<" : ";
1961          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
1962          if (numberQuery>=2) 
1963            parameters[firstMatch].printLongHelp();
1964        }
1965      } else {
1966        if (!numberQuery) 
1967          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
1968                   <<std::endl;
1969        else
1970          std::cout<<"Completions of "<<field<<":"<<std::endl;
1971        for ( iParam=0; iParam<numberParameters; iParam++ ) {
1972          int match = parameters[iParam].matches(field);
1973          if (match&&parameters[iParam].displayThis()) {
1974            std::cout<<parameters[iParam].matchName();
1975            if (numberQuery>=2) 
1976              std::cout<<" : "<<parameters[iParam].shortHelp();
1977            std::cout<<std::endl;
1978          }
1979        }
1980      }
1981    }
1982    delete [] models;
1983    delete [] goodModels;
1984  }
1985  // By now all memory should be freed
1986#ifdef DMALLOC
1987  dmalloc_log_unfreed();
1988  dmalloc_shutdown();
1989#endif
1990  return 0;
1991}   
1992static void breakdown(const char * name, int numberLook, const double * region)
1993{
1994  double range[] = {
1995    -COIN_DBL_MAX,
1996    -1.0e15,-1.0e11,-1.0e8,-1.0e5,-1.0e4,-1.0e3,-1.0e2,-1.0e1,
1997    -1.0,
1998    -1.0e-1,-1.0e-2,-1.0e-3,-1.0e-4,-1.0e-5,-1.0e-8,-1.0e-11,-1.0e-15,
1999    0.0,
2000    1.0e-15,1.0e-11,1.0e-8,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,
2001    1.0,
2002    1.0e1,1.0e2,1.0e3,1.0e4,1.0e5,1.0e8,1.0e11,1.0e15,
2003    COIN_DBL_MAX};
2004  int nRanges = static_cast<int> (sizeof(range)/sizeof(double));
2005  int * number = new int[nRanges];
2006  memset(number,0,nRanges*sizeof(int));
2007  int * numberExact = new int[nRanges];
2008  memset(numberExact,0,nRanges*sizeof(int));
2009  int i;
2010  for ( i=0;i<numberLook;i++) {
2011    double value = region[i];
2012    for (int j=0;j<nRanges;j++) {
2013      if (value==range[j]) {
2014        numberExact[j]++;
2015        break;
2016      } else if (value<range[j]) {
2017        number[j]++;
2018        break;
2019      }
2020    }
2021  }
2022  printf("\n%s has %d entries\n",name,numberLook);
2023  for (i=0;i<nRanges;i++) {
2024    if (number[i]) 
2025      printf("%d between %g and %g",number[i],range[i-1],range[i]);
2026    if (numberExact[i]) {
2027      if (number[i])
2028        printf(", ");
2029      printf("%d exactly at %g",numberExact[i],range[i]);
2030    }
2031    if (number[i]+numberExact[i])
2032      printf("\n");
2033  }
2034  delete [] number;
2035  delete [] numberExact;
2036}
2037void sortOnOther(int * column,
2038                  const CoinBigIndex * rowStart,
2039                  int * order,
2040                  int * other,
2041                  int nRow,
2042                  int nInRow,
2043                  int where)
2044{
2045  if (nRow<2||where>=nInRow)
2046    return;
2047  // do initial sort
2048  int kRow;
2049  int iRow;
2050  for ( kRow=0;kRow<nRow;kRow++) {
2051    iRow = order[kRow];
2052    other[kRow]=column[rowStart[iRow]+where];
2053  }
2054  CoinSort_2(other,other+nRow,order);
2055  int first=0;
2056  iRow=order[0];
2057  int firstC=column[rowStart[iRow]+where];
2058  kRow=1;
2059  while (kRow<nRow) {
2060    int lastC=9999999;;
2061    for (;kRow<nRow+1;kRow++) {
2062      if (kRow<nRow) {
2063        iRow=order[kRow];
2064        lastC=column[rowStart[iRow]+where];
2065      } else {
2066        lastC=9999999;
2067      }
2068      if (lastC>firstC) 
2069        break;
2070    }
2071    // sort
2072    sortOnOther(column,rowStart,order+first,other,kRow-first,
2073                 nInRow,where+1);
2074    firstC=lastC;
2075    first=kRow;
2076  }
2077}
2078static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
2079{
2080  int numberColumns = originalModel->numberColumns();
2081  const char * integerInformation  = originalModel->integerInformation(); 
2082  const double * columnLower = originalModel->columnLower();
2083  const double * columnUpper = originalModel->columnUpper();
2084  int numberIntegers=0;
2085  int numberBinary=0;
2086  int iRow,iColumn;
2087  if (integerInformation) {
2088    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2089      if (integerInformation[iColumn]) {
2090        if (columnUpper[iColumn]>columnLower[iColumn]) {
2091          numberIntegers++;
2092          if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==1) 
2093            numberBinary++;
2094        }
2095      }
2096    }
2097  }
2098  numberColumns = model->numberColumns();
2099  int numberRows = model->numberRows();
2100  columnLower = model->columnLower();
2101  columnUpper = model->columnUpper();
2102  const double * rowLower = model->rowLower();
2103  const double * rowUpper = model->rowUpper();
2104  const double * objective = model->objective();
2105  CoinPackedMatrix * matrix = model->matrix();
2106  CoinBigIndex numberElements = matrix->getNumElements();
2107  const int * columnLength = matrix->getVectorLengths();
2108  //const CoinBigIndex * columnStart = matrix->getVectorStarts();
2109  const double * elementByColumn = matrix->getElements();
2110  int * number = new int[numberRows+1];
2111  memset(number,0,(numberRows+1)*sizeof(int));
2112  int numberObjSingletons=0;
2113  /* cType
2114     0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
2115     8 0/1
2116  */ 
2117  int cType[9];
2118  std::string cName[]={"0.0->inf,","0.0->up,","lo->inf,","lo->up,","free,","fixed,","-inf->0.0,",
2119                       "-inf->up,","0.0->1.0"};
2120  int nObjective=0;
2121  memset(cType,0,sizeof(cType));
2122  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2123    int length=columnLength[iColumn];
2124    if (length==1&&objective[iColumn])
2125      numberObjSingletons++;
2126    number[length]++;
2127    if (objective[iColumn])
2128      nObjective++;
2129    if (columnLower[iColumn]>-1.0e20) {
2130      if (columnLower[iColumn]==0.0) {
2131        if (columnUpper[iColumn]>1.0e20)
2132          cType[0]++;
2133        else if (columnUpper[iColumn]==1.0)
2134          cType[8]++;
2135        else if (columnUpper[iColumn]==0.0)
2136          cType[5]++;
2137        else
2138          cType[1]++;
2139      } else {
2140        if (columnUpper[iColumn]>1.0e20) 
2141          cType[2]++;
2142        else if (columnUpper[iColumn]==columnLower[iColumn])
2143          cType[5]++;
2144        else
2145          cType[3]++;
2146      }
2147    } else {
2148      if (columnUpper[iColumn]>1.0e20) 
2149        cType[4]++;
2150      else if (columnUpper[iColumn]==0.0) 
2151        cType[6]++;
2152      else
2153        cType[7]++;
2154    }
2155  }
2156  /* rType
2157     0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
2158     7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
2159  */ 
2160  int rType[13];
2161  std::string rName[]={"E 0.0,","E 1.0,","E -1.0,","E other,","G 0.0,","G 1.0,","G other,",
2162                       "L 0.0,","L 1.0,","L other,","Range 0.0->1.0,","Range other,","Free"};
2163  memset(rType,0,sizeof(rType));
2164  for (iRow=0;iRow<numberRows;iRow++) {
2165    if (rowLower[iRow]>-1.0e20) {
2166      if (rowLower[iRow]==0.0) {
2167        if (rowUpper[iRow]>1.0e20)
2168          rType[4]++;
2169        else if (rowUpper[iRow]==1.0)
2170          rType[10]++;
2171        else if (rowUpper[iRow]==0.0)
2172          rType[0]++;
2173        else
2174          rType[11]++;
2175      } else if (rowLower[iRow]==1.0) {
2176        if (rowUpper[iRow]>1.0e20) 
2177          rType[5]++;
2178        else if (rowUpper[iRow]==rowLower[iRow])
2179          rType[1]++;
2180        else
2181          rType[11]++;
2182      } else if (rowLower[iRow]==-1.0) {
2183        if (rowUpper[iRow]>1.0e20) 
2184          rType[6]++;
2185        else if (rowUpper[iRow]==rowLower[iRow])
2186          rType[2]++;
2187        else
2188          rType[11]++;
2189      } else {
2190        if (rowUpper[iRow]>1.0e20) 
2191          rType[6]++;
2192        else if (rowUpper[iRow]==rowLower[iRow])
2193          rType[3]++;
2194        else
2195          rType[11]++;
2196      }
2197    } else {
2198      if (rowUpper[iRow]>1.0e20) 
2199        rType[12]++;
2200      else if (rowUpper[iRow]==0.0) 
2201        rType[7]++;
2202      else if (rowUpper[iRow]==1.0) 
2203        rType[8]++;
2204      else
2205        rType[9]++;
2206    }
2207  }
2208  // Basic statistics
2209  printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
2210         numberRows,numberColumns,nObjective,numberElements);
2211  if (number[0]+number[1]) {
2212    printf("There are ");
2213    if (numberObjSingletons)
2214      printf("%d singletons with objective ",numberObjSingletons);
2215    int numberNoObj = number[1]-numberObjSingletons;
2216    if (numberNoObj)
2217      printf("%d singletons with no objective ",numberNoObj);
2218    if (number[0])
2219      printf("** %d columns have no entries",number[0]);
2220    printf("\n");
2221  }
2222  printf("Column breakdown:\n");
2223  int k;
2224  for (k=0;k<static_cast<int> (sizeof(cType)/sizeof(int));k++) {
2225    printf("%d of type %s ",cType[k],cName[k].c_str());
2226    if (((k+1)%3)==0)
2227      printf("\n");
2228  }
2229  if ((k%3)!=0)
2230    printf("\n");
2231  printf("Row breakdown:\n");
2232  for (k=0;k<static_cast<int> (sizeof(rType)/sizeof(int));k++) {
2233    printf("%d of type %s ",rType[k],rName[k].c_str());
2234    if (((k+1)%3)==0)
2235      printf("\n");
2236  }
2237  if ((k%3)!=0)
2238    printf("\n");
2239#define SYM
2240#ifndef SYM
2241  if (model->logLevel()<2)
2242    return ;
2243#endif
2244  int kMax = model->logLevel()>3 ? 1000000 : 10;
2245  k=0;
2246  for (iRow=1;iRow<=numberRows;iRow++) {
2247    if (number[iRow]) {
2248      k++;
2249      printf("%d columns have %d entries\n",number[iRow],iRow);
2250      if (k==kMax)
2251        break;
2252    }
2253  }
2254  if (k<numberRows) {
2255    int kk=k;
2256    k=0;
2257    for (iRow=numberRows;iRow>=1;iRow--) {
2258      if (number[iRow]) {
2259        k++;
2260        if (k==kMax)
2261          break;
2262      }
2263    }
2264    if (k>kk) {
2265      printf("\n    .........\n\n");
2266      iRow=k;
2267      k=0;
2268      for (;iRow<numberRows;iRow++) {
2269        if (number[iRow]) {
2270          k++;
2271          printf("%d columns have %d entries\n",number[iRow],iRow);
2272          if (k==kMax)
2273            break;
2274        }
2275      }
2276    }
2277  }
2278  delete [] number;
2279  printf("\n\n");
2280  if (model->logLevel()==63
2281#ifdef SYM
2282      ||true
2283#endif
2284      ) {
2285    // get column copy
2286    CoinPackedMatrix columnCopy = *matrix;
2287    const int * columnLength = columnCopy.getVectorLengths();
2288    number = new int[numberRows+1];
2289    memset(number,0,(numberRows+1)*sizeof(int));
2290    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2291      int length=columnLength[iColumn];
2292      number[length]++;
2293    }
2294    k=0;
2295    for (iRow=1;iRow<=numberRows;iRow++) {
2296      if (number[iRow]) {
2297        k++;
2298      }
2299    }
2300    int * row = columnCopy.getMutableIndices();
2301    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2302    double * element = columnCopy.getMutableElements();
2303    int * order = new int[numberColumns];
2304    int * other = new int[numberColumns];
2305    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2306      int length=columnLength[iColumn];
2307      order[iColumn]=iColumn;
2308      other[iColumn]=length;
2309      CoinBigIndex start = columnStart[iColumn];
2310      CoinSort_2(row+start,row+start+length,element+start);
2311    }
2312    CoinSort_2(other,other+numberColumns,order);
2313    int jColumn=number[0]+number[1];
2314    for (iRow=2;iRow<=numberRows;iRow++) {
2315      if (number[iRow]) {
2316        printf("XX %d columns have %d entries\n",number[iRow],iRow);
2317        int kColumn = jColumn+number[iRow];
2318        sortOnOther(row,columnStart,
2319                     order+jColumn,other,number[iRow],iRow,0);
2320        // Now print etc
2321        if (iRow<500000) {
2322          for (int lColumn =jColumn;lColumn<kColumn;lColumn++) {
2323            iColumn = order[lColumn];
2324            CoinBigIndex start = columnStart[iColumn];
2325            if (model->logLevel()==63) {
2326              printf("column %d %g <= ",iColumn,columnLower[iColumn]);
2327              for (CoinBigIndex i=start;i<start+iRow;i++)
2328                printf("( %d, %g) ",row[i],element[i]);
2329              printf("<= %g\n",columnUpper[iColumn]);
2330            }
2331          }
2332        }
2333        jColumn =kColumn;
2334      }
2335    }
2336    delete [] order;
2337    delete [] other;
2338    delete [] number;
2339  }
2340  // get row copy
2341  CoinPackedMatrix rowCopy = *matrix;
2342  rowCopy.reverseOrdering();
2343  const int * rowLength = rowCopy.getVectorLengths();
2344  number = new int[numberColumns+1];
2345  memset(number,0,(numberColumns+1)*sizeof(int));
2346  for (iRow=0;iRow<numberRows;iRow++) {
2347    int length=rowLength[iRow];
2348    number[length]++;
2349  }
2350  if (number[0])
2351    printf("** %d rows have no entries\n",number[0]);
2352  k=0;
2353  for (iColumn=1;iColumn<=numberColumns;iColumn++) {
2354    if (number[iColumn]) {
2355      k++;
2356      printf("%d rows have %d entries\n",number[iColumn],iColumn);
2357      if (k==kMax)
2358        break;
2359    }
2360  }
2361  if (k<numberColumns) {
2362    int kk=k;
2363    k=0;
2364    for (iColumn=numberColumns;iColumn>=1;iColumn--) {
2365      if (number[iColumn]) {
2366        k++;
2367        if (k==kMax)
2368          break;
2369      }
2370    }
2371    if (k>kk) {
2372      printf("\n    .........\n\n");
2373      iColumn=k;
2374      k=0;
2375      for (;iColumn<numberColumns;iColumn++) {
2376        if (number[iColumn]) {
2377          k++;
2378          printf("%d rows have %d entries\n",number[iColumn],iColumn);
2379          if (k==kMax)
2380            break;
2381        }
2382      }
2383    }
2384  }
2385  if (model->logLevel()==63
2386#ifdef SYM
2387      ||true
2388#endif
2389      ) {
2390    int * column = rowCopy.getMutableIndices();
2391    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
2392    double * element = rowCopy.getMutableElements();
2393    int * order = new int[numberRows];
2394    int * other = new int[numberRows];
2395    for (iRow=0;iRow<numberRows;iRow++) {
2396      int length=rowLength[iRow];
2397      order[iRow]=iRow;
2398      other[iRow]=length;
2399      CoinBigIndex start = rowStart[iRow];
2400      CoinSort_2(column+start,column+start+length,element+start);
2401    }
2402    CoinSort_2(other,other+numberRows,order);
2403    int jRow=number[0]+number[1];
2404    double * weight = new double[numberRows];
2405    double * randomColumn = new double[numberColumns+1];
2406    double * randomRow = new double [numberRows+1];
2407    int * sortRow = new int [numberRows];
2408    int * possibleRow = new int [numberRows];
2409    int * backRow = new int [numberRows];
2410    int * stackRow = new int [numberRows];
2411    int * sortColumn = new int [numberColumns];
2412    int * possibleColumn = new int [numberColumns];
2413    int * backColumn = new int [numberColumns];
2414    int * backColumn2 = new int [numberColumns];
2415    int * mapRow = new int [numberRows];
2416    int * mapColumn = new int [numberColumns];
2417    int * stackColumn = new int [numberColumns];
2418    double randomLower = CoinDrand48();
2419    double randomUpper = CoinDrand48();
2420    double randomInteger = CoinDrand48();
2421    int * startAdd = new int[numberRows+1];
2422    int * columnAdd = new int [2*numberElements];
2423    double * elementAdd = new double[2*numberElements];
2424    int nAddRows=0;
2425    startAdd[0]=0;
2426    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2427      randomColumn[iColumn] = CoinDrand48();
2428      backColumn2[iColumn]=-1;
2429    }
2430    for (iColumn=2;iColumn<=numberColumns;iColumn++) {
2431      if (number[iColumn]) {
2432        printf("XX %d rows have %d entries\n",number[iColumn],iColumn);
2433        int kRow = jRow+number[iColumn];
2434        sortOnOther(column,rowStart,
2435                     order+jRow,other,number[iColumn],iColumn,0);
2436        // Now print etc
2437        if (iColumn<500000) {
2438          int nLook=0;
2439          for (int lRow =jRow;lRow<kRow;lRow++) {
2440            iRow = order[lRow];
2441            CoinBigIndex start = rowStart[iRow];
2442            if (model->logLevel()==63) {
2443              printf("row %d %g <= ",iRow,rowLower[iRow]);
2444              for (CoinBigIndex i=start;i<start+iColumn;i++) 
2445                printf("( %d, %g) ",column[i],element[i]);
2446              printf("<= %g\n",rowUpper[iRow]);
2447            }
2448            int first = column[start];
2449            double sum=0.0;
2450            for (CoinBigIndex i=start;i<start+iColumn;i++) {
2451              int jColumn = column[i];
2452              double value = element[i];
2453              jColumn -= first;
2454              assert (jColumn>=0);
2455              sum += value*randomColumn[jColumn];
2456            }
2457            if (rowLower[iRow]>-1.0e30&&rowLower[iRow])
2458              sum += rowLower[iRow]*randomLower;
2459            else if (!rowLower[iRow])
2460              sum += 1.234567e-7*randomLower;
2461            if (rowUpper[iRow]<1.0e30&&rowUpper[iRow])
2462              sum += rowUpper[iRow]*randomUpper;
2463            else if (!rowUpper[iRow])
2464              sum += 1.234567e-7*randomUpper;
2465            sortRow[nLook]=iRow;
2466            randomRow[nLook++]=sum;
2467            // best way is to number unique elements and bounds and use
2468            if (fabs(sum)>1.0e4)
2469              sum *= 1.0e-6;
2470            weight[iRow]=sum;
2471          }
2472          assert (nLook<=numberRows);
2473          CoinSort_2(randomRow,randomRow+nLook,sortRow);
2474          randomRow[nLook]=COIN_DBL_MAX;
2475          double last=-COIN_DBL_MAX;
2476          int iLast=-1;
2477          for (int iLook=0;iLook<nLook+1;iLook++) {
2478            if (randomRow[iLook]>last) {
2479              if (iLast>=0) {
2480                int n=iLook-iLast;
2481                if (n>1) {
2482                  //printf("%d rows possible?\n",n);
2483                }
2484              }
2485              iLast=iLook;
2486              last=randomRow[iLook];
2487            }
2488          }
2489        }
2490        jRow =kRow;
2491      }
2492    }
2493    CoinPackedMatrix columnCopy = *matrix;
2494    const int * columnLength = columnCopy.getVectorLengths();
2495    const int * row = columnCopy.getIndices();
2496    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2497    const double * elementByColumn = columnCopy.getElements();
2498    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2499      int length=columnLength[iColumn];
2500      CoinBigIndex start = columnStart[iColumn];
2501      double sum = objective[iColumn];
2502      if (columnLower[iColumn]>-1.0e30&&columnLower[iColumn])
2503        sum += columnLower[iColumn]*randomLower;
2504      else if (!columnLower[iColumn])
2505        sum += 1.234567e-7*randomLower;
2506      if (columnUpper[iColumn]<1.0e30&&columnUpper[iColumn])
2507        sum += columnUpper[iColumn]*randomUpper;
2508      else if (!columnUpper[iColumn])
2509        sum += 1.234567e-7*randomUpper;
2510      if (model->isInteger(iColumn))
2511        sum += 9.87654321e-6*randomInteger;
2512      for (CoinBigIndex i=start;i<start+length;i++) {
2513        int iRow = row[i];
2514        sum += elementByColumn[i]*weight[iRow];
2515      }
2516      sortColumn[iColumn]=iColumn;
2517      randomColumn[iColumn]=sum;
2518    }
2519    {
2520      CoinSort_2(randomColumn,randomColumn+numberColumns,sortColumn);
2521      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2522        int i=sortColumn[iColumn];
2523        backColumn[i]=iColumn;
2524      }
2525      randomColumn[numberColumns]=COIN_DBL_MAX;
2526      double last=-COIN_DBL_MAX;
2527      int iLast=-1;
2528      for (int iLook=0;iLook<numberColumns+1;iLook++) {
2529        if (randomColumn[iLook]>last) {
2530          if (iLast>=0) {
2531            int n=iLook-iLast;
2532            if (n>1) {
2533              //printf("%d columns possible?\n",n);
2534            }
2535            for (int i=iLast;i<iLook;i++) {
2536              possibleColumn[sortColumn[i]]=n;
2537            }
2538          }
2539          iLast=iLook;
2540          last=randomColumn[iLook];
2541        }
2542      }
2543      for (iRow =0;iRow<numberRows;iRow++) {
2544        CoinBigIndex start = rowStart[iRow];
2545        double sum=0.0;
2546        int length=rowLength[iRow];
2547        for (CoinBigIndex i=start;i<start+length;i++) {
2548          int jColumn = column[i];
2549          double value = element[i];
2550          jColumn=backColumn[jColumn];
2551          sum += value*randomColumn[jColumn];
2552          //if (iColumn==23089||iRow==23729)
2553          //printf("row %d cola %d colb %d value %g rand %g sum %g\n",
2554          //   iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
2555        }
2556        sortRow[iRow]=iRow;
2557        randomRow[iRow]=weight[iRow];
2558        randomRow[iRow]=sum;
2559      }
2560      CoinSort_2(randomRow,randomRow+numberRows,sortRow);
2561      for (iRow=0;iRow<numberRows;iRow++) {
2562        int i=sortRow[iRow];
2563        backRow[i]=iRow;
2564      }
2565      randomRow[numberRows]=COIN_DBL_MAX;
2566      last=-COIN_DBL_MAX;
2567      iLast=-1;
2568      // Do backward indices from order
2569      for (iRow=0;iRow<numberRows;iRow++) {
2570        other[order[iRow]]=iRow;
2571      }
2572      for (int iLook=0;iLook<numberRows+1;iLook++) {
2573        if (randomRow[iLook]>last) {
2574          if (iLast>=0) {
2575            int n=iLook-iLast;
2576            if (n>1) {
2577              //printf("%d rows possible?\n",n);
2578              // Within group sort as for original "order"
2579              for (int i=iLast;i<iLook;i++) {
2580                int jRow=sortRow[i];
2581                order[i]=other[jRow];
2582              }
2583              CoinSort_2(order+iLast,order+iLook,sortRow+iLast);
2584            }
2585            for (int i=iLast;i<iLook;i++) {
2586              possibleRow[sortRow[i]]=n;
2587            }
2588          }
2589          iLast=iLook;
2590          last=randomRow[iLook];
2591        }
2592      }
2593      // Temp out
2594      for (int iLook=0;iLook<numberRows-1000000;iLook++) {
2595        iRow=sortRow[iLook];
2596        CoinBigIndex start = rowStart[iRow];
2597        int length=rowLength[iRow];
2598        int numberPossible = possibleRow[iRow];
2599        for (CoinBigIndex i=start;i<start+length;i++) {
2600          int jColumn = column[i];
2601          if (possibleColumn[jColumn]!=numberPossible)
2602            numberPossible=-1;
2603        }
2604        int n=numberPossible;
2605        if (numberPossible>1) {
2606          //printf("pppppossible %d\n",numberPossible);
2607          for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2608            int jRow=sortRow[jLook];
2609            CoinBigIndex start2 = rowStart[jRow];
2610            assert (numberPossible==possibleRow[jRow]);
2611            assert(length==rowLength[jRow]);
2612            for (CoinBigIndex i=start2;i<start2+length;i++) {
2613              int jColumn = column[i];
2614              if (possibleColumn[jColumn]!=numberPossible)
2615                numberPossible=-1;
2616            }
2617          }
2618          if (numberPossible<2) {
2619            // switch off
2620            for (int jLook=iLook;jLook<iLook+n;jLook++) 
2621              possibleRow[sortRow[jLook]]=-1;
2622          }
2623          // skip rest
2624          iLook += n-1;
2625        } else {
2626          possibleRow[iRow]=-1;
2627        }
2628      }
2629      for (int iLook=0;iLook<numberRows;iLook++) {
2630        iRow=sortRow[iLook];
2631        int numberPossible = possibleRow[iRow];
2632        // Only if any integers
2633        int numberIntegers=0;
2634        CoinBigIndex start = rowStart[iRow];
2635        int length=rowLength[iRow];
2636        for (CoinBigIndex i=start;i<start+length;i++) {
2637          int jColumn = column[i];
2638          if (model->isInteger(jColumn))
2639            numberIntegers++;
2640        }
2641        if (numberPossible>1&&!numberIntegers) {
2642          //printf("possible %d - but no integers\n",numberPossible);
2643        }
2644        if (numberPossible>1&&(numberIntegers||false)) {
2645          //
2646          printf("possible %d - %d integers\n",numberPossible,numberIntegers);
2647          int lastLook=iLook;
2648          int nMapRow=-1;
2649          for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2650            // stop if too many failures
2651            if (jLook>iLook+10&&nMapRow<0)
2652              break;
2653            // Create identity mapping
2654            int i;
2655            for (i=0;i<numberRows;i++)
2656              mapRow[i]=i;
2657            for (i=0;i<numberColumns;i++)
2658              mapColumn[i]=i;
2659            int offset=jLook-iLook;
2660            int nStackC=0;
2661            // build up row and column mapping
2662            int nStackR=1;
2663            stackRow[0]=iLook;
2664            bool good=true;
2665            while (nStackR) {
2666              nStackR--;
2667              int look1 = stackRow[nStackR];
2668              int look2 = look1 + offset;
2669              assert (randomRow[look1]==randomRow[look2]);
2670              int row1=sortRow[look1];
2671              int row2=sortRow[look2];
2672              assert (mapRow[row1]==row1);
2673              assert (mapRow[row2]==row2);
2674              mapRow[row1]=row2;
2675              mapRow[row2]=row1;
2676              CoinBigIndex start1 = rowStart[row1];
2677              CoinBigIndex offset2 = rowStart[row2]-start1;
2678              int length=rowLength[row1];
2679              assert( length==rowLength[row2]);
2680              for (CoinBigIndex i=start1;i<start1+length;i++) {
2681                int jColumn1 = column[i];
2682                int jColumn2 = column[i+offset2];
2683                if (randomColumn[backColumn[jColumn1]]!=
2684                    randomColumn[backColumn[jColumn2]]) {
2685                  good=false;
2686                  break;
2687                }
2688                if (mapColumn[jColumn1]==jColumn1) {
2689                  // not touched
2690                  assert (mapColumn[jColumn2]==jColumn2);
2691                  if (jColumn1!=jColumn2) {
2692                    // Put on stack
2693                    mapColumn[jColumn1]=jColumn2;
2694                    mapColumn[jColumn2]=jColumn1;
2695                    stackColumn[nStackC++]=jColumn1;
2696                  }
2697                } else {
2698                  if(mapColumn[jColumn1]!=jColumn2||
2699                     mapColumn[jColumn2]!=jColumn1) {
2700                    // bad
2701                    good=false;
2702                    printf("bad col\n");
2703                    break;
2704                  }
2705                }
2706              }
2707              if (!good)
2708                break;
2709              while (nStackC) {
2710                nStackC--;
2711                int iColumn = stackColumn[nStackC];
2712                int iColumn2 = mapColumn[iColumn];
2713                assert (iColumn!=iColumn2);
2714                int length=columnLength[iColumn];
2715                assert (length==columnLength[iColumn2]);
2716                CoinBigIndex start = columnStart[iColumn];
2717                CoinBigIndex offset2 = columnStart[iColumn2]-start;
2718                for (CoinBigIndex i=start;i<start+length;i++) {
2719                  int iRow = row[i];
2720                  int iRow2 = row[i+offset2];
2721                  if (mapRow[iRow]==iRow) {
2722                    // First (but be careful)
2723                    if (iRow!=iRow2) {
2724                      //mapRow[iRow]=iRow2;
2725                      //mapRow[iRow2]=iRow;
2726                      int iBack=backRow[iRow];
2727                      int iBack2=backRow[iRow2];
2728                      if (randomRow[iBack]==randomRow[iBack2]&&
2729                          iBack2-iBack==offset) {
2730                        stackRow[nStackR++]=iBack;
2731                      } else {
2732                        //printf("randomRow diff - weights %g %g\n",
2733                        //     weight[iRow],weight[iRow2]);
2734                        // bad
2735                        good=false;
2736                        break;
2737                      }
2738                    }
2739                  } else {
2740                    if(mapRow[iRow]!=iRow2||
2741                       mapRow[iRow2]!=iRow) {
2742                      // bad
2743                      good=false;
2744                      printf("bad row\n");
2745                      break;
2746                    }
2747                  }
2748                }
2749                if (!good)
2750                  break;
2751              }
2752            }
2753            // then check OK
2754            if (good) {
2755              for (iRow =0;iRow<numberRows;iRow++) {
2756                CoinBigIndex start = rowStart[iRow];
2757                int length=rowLength[iRow];
2758                if (mapRow[iRow]==iRow) {
2759                  for (CoinBigIndex i=start;i<start+length;i++) {
2760                    int jColumn = column[i];
2761                    backColumn2[jColumn]=i-start;
2762                  }
2763                  for (CoinBigIndex i=start;i<start+length;i++) {
2764                    int jColumn = column[i];
2765                    if (mapColumn[jColumn]!=jColumn) {
2766                      int jColumn2 =mapColumn[jColumn];
2767                      CoinBigIndex i2 = backColumn2[jColumn2];
2768                      if (i2<0) {
2769                        good=false;
2770                      } else if (element[i]!=element[i2+start]) {
2771                        good=false;
2772                      }
2773                    }
2774                  }
2775                  for (CoinBigIndex i=start;i<start+length;i++) {
2776                    int jColumn = column[i];
2777                    backColumn2[jColumn]=-1;
2778                  }
2779                } else {
2780                  int row2 = mapRow[iRow];
2781                  assert (iRow = mapRow[row2]);
2782                  if (rowLower[iRow]!=rowLower[row2]||
2783                      rowLower[row2]!=rowLower[iRow])
2784                    good=false;
2785                  CoinBigIndex offset2 = rowStart[row2]-start;
2786                  for (CoinBigIndex i=start;i<start+length;i++) {
2787                    int jColumn = column[i];
2788                    double value = element[i];
2789                    int jColumn2 = column[i+offset2];
2790                    double value2 = element[i+offset2];
2791                    if (value!=value2||mapColumn[jColumn]!=jColumn2||
2792                        mapColumn[jColumn2]!=jColumn)
2793                      good=false;
2794                  }
2795                }
2796              }
2797              if (good) {
2798                // check rim
2799                for (iColumn=0;iColumn<numberColumns;iColumn++) {
2800                  if (mapColumn[iColumn]!=iColumn) {
2801                    int iColumn2 = mapColumn[iColumn];
2802                    if(objective[iColumn]!=objective[iColumn2])
2803                      good=false;
2804                    if (columnLower[iColumn]!=columnLower[iColumn2])
2805                      good=false;
2806                    if (columnUpper[iColumn]!=columnUpper[iColumn2])
2807                      good=false;
2808                    if (model->isInteger(iColumn)!=model->isInteger(iColumn2))
2809                      good=false;
2810                  }
2811                }
2812              }
2813              if (good) {
2814                // temp
2815                if (nMapRow<0) {
2816                  //const double * solution = model->primalColumnSolution();
2817                  // find mapped
2818                  int nMapColumn=0;
2819                  for (int i=0;i<numberColumns;i++) {
2820                    if (mapColumn[i]>i) 
2821                      nMapColumn++;
2822                  }
2823                  nMapRow=0;
2824                  int kRow=-1;
2825                  for (int i=0;i<numberRows;i++) {
2826                    if (mapRow[i]>i) { 
2827                      nMapRow++;
2828                      kRow=i;
2829                    }
2830                  }
2831                  printf("%d columns, %d rows\n",nMapColumn,nMapRow);
2832                  if (nMapRow==1) {
2833                    CoinBigIndex start = rowStart[kRow];
2834                    int length=rowLength[kRow];
2835                    printf("%g <= ",rowLower[kRow]);
2836                    for (CoinBigIndex i=start;i<start+length;i++) {
2837                      int jColumn = column[i];
2838                      if (mapColumn[jColumn]!=jColumn) 
2839                        printf("* ");
2840                      printf("%d,%g ",jColumn,element[i]);
2841                    }
2842                    printf("<= %g\n",rowUpper[kRow]);
2843                  }
2844                }
2845                // temp
2846                int row1=sortRow[lastLook];
2847                int row2=sortRow[jLook];
2848                lastLook=jLook;
2849                CoinBigIndex start1 = rowStart[row1];
2850                CoinBigIndex offset2 = rowStart[row2]-start1;
2851                int length=rowLength[row1];
2852                assert( length==rowLength[row2]);
2853                CoinBigIndex put=startAdd[nAddRows];
2854                double multiplier=length<11 ? 2.0 : 1.125;
2855                double value=1.0;
2856                for (CoinBigIndex i=start1;i<start1+length;i++) {
2857                  int jColumn1 = column[i];
2858                  int jColumn2 = column[i+offset2];
2859                  columnAdd[put]=jColumn1;
2860                  elementAdd[put++]=value;
2861                  columnAdd[put]=jColumn2;
2862                  elementAdd[put++]=-value;
2863                  value *= multiplier;
2864                }
2865                nAddRows++;
2866                startAdd[nAddRows]=put;
2867              } else {
2868                printf("ouch - did not check out as good\n");
2869              }
2870            }
2871          }
2872          // skip rest
2873          iLook += numberPossible-1;
2874        }
2875      }
2876    }
2877    if (nAddRows) {
2878      double * lower = new double [nAddRows];
2879      double * upper = new double[nAddRows];
2880      int i;
2881      //const double * solution = model->primalColumnSolution();
2882      for (i=0;i<nAddRows;i++) {
2883        lower[i]=0.0;
2884        upper[i]=COIN_DBL_MAX;
2885      }
2886      printf("Adding %d rows with %d elements\n",nAddRows,
2887             startAdd[nAddRows]);
2888      //ClpSimplex newModel(*model);
2889      //newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
2890      //newModel.writeMps("modified.mps");
2891      delete [] lower;
2892      delete [] upper;
2893    }
2894    delete [] startAdd;
2895    delete [] columnAdd;
2896    delete [] elementAdd;
2897    delete [] order;
2898    delete [] other;
2899    delete [] randomColumn;
2900    delete [] weight;
2901    delete [] randomRow;
2902    delete [] sortRow;
2903    delete [] backRow;
2904    delete [] possibleRow;
2905    delete [] sortColumn;
2906    delete [] backColumn;
2907    delete [] backColumn2;
2908    delete [] possibleColumn;
2909    delete [] mapRow;
2910    delete [] mapColumn;
2911    delete [] stackRow;
2912    delete [] stackColumn;
2913  }
2914  delete [] number;
2915  // Now do breakdown of ranges
2916  breakdown("Elements",numberElements,elementByColumn);
2917  breakdown("RowLower",numberRows,rowLower);
2918  breakdown("RowUpper",numberRows,rowUpper);
2919  breakdown("ColumnLower",numberColumns,columnLower);
2920  breakdown("ColumnUpper",numberColumns,columnUpper);
2921  breakdown("Objective",numberColumns,objective);
2922}
2923static bool maskMatches(const int * starts, char ** masks,
2924                        std::string & check)
2925{
2926  // back to char as I am old fashioned
2927  const char * checkC = check.c_str();
2928  int length = strlen(checkC);
2929  while (checkC[length-1]==' ')
2930    length--;
2931  for (int i=starts[length];i<starts[length+1];i++) {
2932    char * thisMask = masks[i];
2933    int k;
2934    for ( k=0;k<length;k++) {
2935      if (thisMask[k]!='?'&&thisMask[k]!=checkC[k]) 
2936        break;
2937    }
2938    if (k==length)
2939      return true;
2940  }
2941  return false;
2942}
2943static void clean(char * temp)
2944{
2945  char * put = temp;
2946  while (*put>=' ')
2947    put++;
2948  *put='\0';
2949}
2950static void generateCode(const char * fileName,int type)
2951{
2952  FILE * fp = fopen(fileName,"r");
2953  assert (fp);
2954  int numberLines=0;
2955#define MAXLINES 500
2956#define MAXONELINE 200
2957  char line[MAXLINES][MAXONELINE];
2958  while (fgets(line[numberLines],MAXONELINE,fp)) {
2959    assert (numberLines<MAXLINES);
2960    clean(line[numberLines]);
2961    numberLines++;
2962  }
2963  fclose(fp);
2964  // add in actual solve
2965  strcpy(line[numberLines],"5  clpModel->initialSolve(clpSolve);");
2966  numberLines++;
2967  fp = fopen(fileName,"w");
2968  assert (fp);
2969  char apo='"';   
2970  char backslash = '\\';
2971
2972  fprintf(fp,"#include %cClpSimplex.hpp%c\n",apo,apo);
2973  fprintf(fp,"#include %cClpSolve.hpp%c\n",apo,apo);
2974
2975  fprintf(fp,"\nint main (int argc, const char *argv[])\n{\n");
2976  fprintf(fp,"  ClpSimplex  model;\n");
2977  fprintf(fp,"  int status=1;\n");
2978  fprintf(fp,"  if (argc<2)\n");
2979  fprintf(fp,"    fprintf(stderr,%cPlease give file name%cn%c);\n",
2980          apo,backslash,apo);
2981  fprintf(fp,"  else\n");
2982  fprintf(fp,"    status=model.readMps(argv[1],true);\n");
2983  fprintf(fp,"  if (status) {\n");
2984  fprintf(fp,"    fprintf(stderr,%cBad readMps %%s%cn%c,argv[1]);\n",
2985                apo,backslash,apo);
2986  fprintf(fp,"    exit(1);\n");
2987  fprintf(fp,"  }\n\n");
2988  fprintf(fp,"  // Now do requested saves and modifications\n");
2989  fprintf(fp,"  ClpSimplex * clpModel = & model;\n");
2990  int wanted[9];
2991  memset(wanted,0,sizeof(wanted));
2992  wanted[0]=wanted[3]=wanted[5]=wanted[8]=1;
2993  if (type>0) 
2994    wanted[1]=wanted[6]=1;
2995  if (type>1) 
2996    wanted[2]=wanted[4]=wanted[7]=1;
2997  std::string header[9]=
2998  { "","Save values","Redundant save of default values","Set changed values",
2999    "Redundant set default values","Solve","Restore values","Redundant restore values","Add to model"};
3000  for (int iType=0;iType<9;iType++) {
3001    if (!wanted[iType])
3002      continue;
3003    int n=0;
3004    int iLine;
3005    for (iLine=0;iLine<numberLines;iLine++) {
3006      if (line[iLine][0]=='0'+iType) {
3007        if (!n)
3008          fprintf(fp,"\n  // %s\n\n",header[iType].c_str());
3009        n++;
3010        fprintf(fp,"%s\n",line[iLine]+1);
3011      }
3012    }
3013  }
3014  fprintf(fp,"\n  // Now you would use solution etc etc\n\n");
3015  fprintf(fp,"  return 0;\n}\n");
3016  fclose(fp);
3017  printf("C++ file written to %s\n",fileName);
3018}
3019/*
3020  Version 1.00.00 October 13 2004.
3021  1.00.01 October 18.  Added basis handling helped/prodded by Thorsten Koch.
3022  Also modifications to make faster with sbb (I hope I haven't broken anything).
3023  1.00.02 March 21 2005.  Redid ClpNonLinearCost to save memory also redid
3024  createRim to try and improve cache characteristics.
3025  1.00.03 April 8 2005.  Added Volume algorithm as crash and made code more
3026  robust on testing.  Also added "either" and "tune" algorithm.
3027  1.01.01 April 12 2005.  Decided to go to different numbering.  Backups will
3028  be last 2 digits while middle 2 are for improvements.  Still take a long
3029  time to get to 2.00.01
3030  1.01.02 May 4 2005.  Will be putting in many changes - so saving stable version
3031  1.02.01 May 6 2005.  Lots of changes to try and make faster and more stable in
3032  branch and cut.
3033  1.02.02 May 19 2005.  Stuff for strong branching and some improvements to simplex
3034  1.03.01 May 24 2006.  Lots done but I can't remember what!
3035  1.03.03 June 13 2006.  For clean up after dual perturbation
3036  1.04.01 June 26 2007.  Lots of changes but I got lazy
3037  1.05.00 June 27 2007.  This is trunk so when gets to stable will be 1.5
3038 */
Note: See TracBrowser for help on using the repository browser.