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

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

minor changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 93.6 KB
Line 
1/* $Id: ClpMain.cpp 1434 2009-09-24 20:48:01Z 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, ClpSolve solveOptions,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              ClpSolve solveOptions;
1557              ClpSolve::PresolveType presolveType;
1558              if (preSolve)
1559                presolveType=ClpSolve::presolveOn;
1560              else
1561                presolveType=ClpSolve::presolveOff;
1562              solveOptions.setPresolveType(presolveType,5);
1563              if (doSprint>=0||doIdiot>=0) {
1564                if (doSprint>0) {
1565                  // sprint overrides idiot
1566                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
1567                } else if (doIdiot>0) {
1568                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
1569                } else {
1570                  if (doIdiot==0) {
1571                    if (doSprint==0)
1572                      solveOptions.setSpecialOption(1,4); // all slack
1573                    else
1574                      solveOptions.setSpecialOption(1,9); // all slack or sprint
1575                  } else {
1576                    if (doSprint==0)
1577                      solveOptions.setSpecialOption(1,8); // all slack or idiot
1578                    else
1579                      solveOptions.setSpecialOption(1,7); // initiative
1580                  }
1581                }
1582              }
1583              mainTest(nFields,fields,algorithm,models[iModel],
1584                       solveOptions,specialOptions,doVector!=0);
1585            }
1586            break;
1587          case UNITTEST:
1588            {
1589              // create fields for unitTest
1590              const char * fields[2];
1591              int nFields=2;
1592              fields[0]="fake main from unitTest";
1593              std::string dirfield = "-dirSample=";
1594              dirfield += dirSample.c_str();
1595              fields[1]=dirfield.c_str();
1596              int specialOptions = models[iModel].specialOptions();
1597              models[iModel].setSpecialOptions(0);
1598              int algorithm=-1;
1599              if (models[iModel].numberRows())
1600                algorithm=7;
1601              ClpSolve solveOptions;
1602              ClpSolve::PresolveType presolveType;
1603              if (preSolve)
1604                presolveType=ClpSolve::presolveOn;
1605              else
1606                presolveType=ClpSolve::presolveOff;
1607              solveOptions.setPresolveType(presolveType,5);
1608              mainTest(nFields,fields,algorithm,models[iModel],
1609                       solveOptions,specialOptions,doVector!=0);
1610            }
1611            break;
1612          case FAKEBOUND:
1613            if (goodModels[iModel]) {
1614              // get bound
1615              double value = CoinReadGetDoubleField(argc,argv,&valid);
1616              if (!valid) {
1617                std::cout<<"Setting "<<parameters[iParam].name()<<
1618                  " to DEBUG "<<value<<std::endl;
1619                int iRow;
1620                int numberRows=models[iModel].numberRows();
1621                double * rowLower = models[iModel].rowLower();
1622                double * rowUpper = models[iModel].rowUpper();
1623                for (iRow=0;iRow<numberRows;iRow++) {
1624                  // leave free ones for now
1625                  if (rowLower[iRow]>-1.0e20||rowUpper[iRow]<1.0e20) {
1626                    rowLower[iRow]=CoinMax(rowLower[iRow],-value);
1627                    rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
1628                  }
1629                }
1630                int iColumn;
1631                int numberColumns=models[iModel].numberColumns();
1632                double * columnLower = models[iModel].columnLower();
1633                double * columnUpper = models[iModel].columnUpper();
1634                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1635                  // leave free ones for now
1636                  if (columnLower[iColumn]>-1.0e20||
1637                      columnUpper[iColumn]<1.0e20) {
1638                    columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
1639                    columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
1640                  }
1641                }
1642              } else if (valid==1) {
1643                abort();
1644              } else {
1645                std::cout<<"enter value for "<<parameters[iParam].name()<<
1646                  std::endl;
1647              }
1648            }
1649            break;
1650          case REALLY_SCALE:
1651            if (goodModels[iModel]) {
1652              ClpSimplex newModel(models[iModel],
1653                                  models[iModel].scalingFlag());
1654              printf("model really really scaled\n");
1655              models[iModel]=newModel;
1656            }
1657            break;
1658          case USERCLP:
1659            // Replace the sample code by whatever you want
1660            if (goodModels[iModel]) {
1661              ClpSimplex * thisModel = &models[iModel];
1662              printf("Dummy user code - model has %d rows and %d columns\n",
1663                     thisModel->numberRows(),thisModel->numberColumns());
1664            }
1665            break;
1666          case HELP:
1667            std::cout<<"Coin LP version "<<CLPVERSION
1668                     <<", build "<<__DATE__<<std::endl;
1669            std::cout<<"Non default values:-"<<std::endl;
1670            std::cout<<"Perturbation "<<models[0].perturbation()<<" (default 100)"
1671                     <<std::endl;
1672            CoinReadPrintit(
1673                    "Presolve being done with 5 passes\n\
1674Dual steepest edge steep/partial on matrix shape and factorization density\n\
1675Clpnnnn taken out of messages\n\
1676If Factorization frequency default then done on size of matrix\n\n\
1677(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
1678You can switch to interactive mode at any time so\n\
1679clp watson.mps -scaling off -primalsimplex\nis the same as\n\
1680clp watson.mps -\nscaling off\nprimalsimplex"
1681                    );
1682            break;
1683          case SOLUTION:
1684            if (goodModels[iModel]) {
1685              // get next field
1686              field = CoinReadGetString(argc,argv);
1687              if (field=="$") {
1688                field = parameters[iParam].stringValue();
1689              } else if (field=="EOL") {
1690                parameters[iParam].printString();
1691                break;
1692              } else {
1693                parameters[iParam].setStringValue(field);
1694              }
1695              std::string fileName;
1696              FILE *fp=NULL;
1697              if (field=="-"||field=="EOL"||field=="stdout") {
1698                // stdout
1699                fp=stdout;
1700                fprintf(fp,"\n");
1701              } else if (field=="stderr") {
1702                // stderr
1703                fp=stderr;
1704                fprintf(fp,"\n");
1705              } else {
1706                if (field[0]=='/'||field[0]=='\\') {
1707                  fileName = field;
1708                } else if (field[0]=='~') {
1709                  char * environVar = getenv("HOME");
1710                  if (environVar) {
1711                    std::string home(environVar);
1712                    field=field.erase(0,1);
1713                    fileName = home+field;
1714                  } else {
1715                    fileName=field;
1716                  }
1717                } else {
1718                  fileName = directory+field;
1719                }
1720                fp=fopen(fileName.c_str(),"w");
1721              }
1722              if (fp) {
1723                // Write solution header (suggested by Luigi Poderico)
1724                double objValue = models[iModel].getObjValue()*models[iModel].getObjSense();
1725                int iStat = models[iModel].status();
1726                if (iStat==0) {
1727                  fprintf(fp, "optimal\n" );
1728                } else if (iStat==1) {
1729                  // infeasible
1730                  fprintf(fp, "infeasible\n" );
1731                } else if (iStat==2) {
1732                  // unbounded
1733                  fprintf(fp, "unbounded\n" );
1734                } else if (iStat==3) {
1735                  fprintf(fp, "stopped on iterations or time\n" );
1736                } else if (iStat==4) {
1737                  fprintf(fp, "stopped on difficulties\n" );
1738                } else if (iStat==5) {
1739                  fprintf(fp, "stopped on ctrl-c\n" );
1740                } else {
1741                  fprintf(fp, "status unknown\n" );
1742                }
1743                fprintf(fp, "Objective value %15.8g\n", objValue);
1744                // make fancy later on
1745                int iRow;
1746                int numberRows=models[iModel].numberRows();
1747                int lengthName = models[iModel].lengthNames(); // 0 if no names
1748                int lengthPrint = CoinMax(lengthName,8);
1749                // in general I don't want to pass around massive
1750                // amounts of data but seems simpler here
1751                std::vector<std::string> rowNames =
1752                  *(models[iModel].rowNames());
1753                std::vector<std::string> columnNames =
1754                  *(models[iModel].columnNames());
1755
1756                double * dualRowSolution = models[iModel].dualRowSolution();
1757                double * primalRowSolution = 
1758                  models[iModel].primalRowSolution();
1759                double * rowLower = models[iModel].rowLower();
1760                double * rowUpper = models[iModel].rowUpper();
1761                double primalTolerance = models[iModel].primalTolerance();
1762                bool doMask = (printMask!=""&&lengthName);
1763                int * maskStarts=NULL;
1764                int maxMasks=0;
1765                char ** masks =NULL;
1766                if (doMask) {
1767                  int nAst =0;
1768                  const char * pMask2 = printMask.c_str();
1769                  char pMask[100];
1770                  int iChar;
1771                  int lengthMask = strlen(pMask2);
1772                  assert (lengthMask<100);
1773                  if (*pMask2=='"') {
1774                    if (pMask2[lengthMask-1]!='"') {
1775                      printf("mismatched \" in mask %s\n",pMask2);
1776                      break;
1777                    } else {
1778                      strcpy(pMask,pMask2+1);
1779                      *strchr(pMask,'"')='\0';
1780                    }
1781                  } else if (*pMask2=='\'') {
1782                    if (pMask2[lengthMask-1]!='\'') {
1783                      printf("mismatched ' in mask %s\n",pMask2);
1784                      break;
1785                    } else {
1786                      strcpy(pMask,pMask2+1);
1787                      *strchr(pMask,'\'')='\0';
1788                    }
1789                  } else {
1790                    strcpy(pMask,pMask2);
1791                  }
1792                  if (lengthMask>lengthName) {
1793                    printf("mask %s too long - skipping\n",pMask);
1794                    break;
1795                  }
1796                  maxMasks = 1;
1797                  for (iChar=0;iChar<lengthMask;iChar++) {
1798                    if (pMask[iChar]=='*') {
1799                      nAst++;
1800                      maxMasks *= (lengthName+1);
1801                    }
1802                  }
1803                  int nEntries = 1;
1804                  maskStarts = new int[lengthName+2];
1805                  masks = new char * [maxMasks];
1806                  char ** newMasks = new char * [maxMasks];
1807                  int i;
1808                  for (i=0;i<maxMasks;i++) {
1809                    masks[i] = new char[lengthName+1];
1810                    newMasks[i] = new char[lengthName+1];
1811                  }
1812                  strcpy(masks[0],pMask);
1813                  for (int iAst=0;iAst<nAst;iAst++) {
1814                    int nOldEntries = nEntries;
1815                    nEntries=0;
1816                    for (int iEntry = 0;iEntry<nOldEntries;iEntry++) {
1817                      char * oldMask = masks[iEntry];
1818                      char * ast = strchr(oldMask,'*');
1819                      assert (ast);
1820                      int length = strlen(oldMask)-1;
1821                      int nBefore = ast-oldMask;
1822                      int nAfter = length-nBefore;
1823                      // and add null
1824                      nAfter++;
1825                      for (int i=0;i<=lengthName-length;i++) {
1826                        char * maskOut = newMasks[nEntries];
1827   CoinMemcpyN(oldMask,nBefore,maskOut);
1828                        for (int k=0;k<i;k++) 
1829                          maskOut[k+nBefore]='?';
1830   CoinMemcpyN(ast+1,nAfter,maskOut+nBefore+i);
1831                        nEntries++;
1832                        assert (nEntries<=maxMasks);
1833                      }
1834                    }
1835                    char ** temp = masks;
1836                    masks = newMasks;
1837                    newMasks = temp;
1838                  }
1839                  // Now extend and sort
1840                  int * sort = new int[nEntries];
1841                  for (i=0;i<nEntries;i++) {
1842                    char * maskThis = masks[i];
1843                    int length = strlen(maskThis);
1844                    while (maskThis[length-1]==' ')
1845                      length--;
1846                    maskThis[length]='\0';
1847                    sort[i]=length;
1848                  }
1849                  CoinSort_2(sort,sort+nEntries,masks);
1850                  int lastLength=-1;
1851                  for (i=0;i<nEntries;i++) {
1852                    int length = sort[i];
1853                    while (length>lastLength) 
1854                      maskStarts[++lastLength] = i;
1855                  }
1856                  maskStarts[++lastLength]=nEntries;
1857                  delete [] sort;
1858                  for (i=0;i<maxMasks;i++)
1859                    delete [] newMasks[i];
1860                  delete [] newMasks;
1861                }
1862                if (printMode>2) {
1863                  for (iRow=0;iRow<numberRows;iRow++) {
1864                    int type=printMode-3;
1865                    if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
1866                        primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
1867                      fprintf(fp,"** ");
1868                      type=2;
1869                    } else if (fabs(primalRowSolution[iRow])>1.0e-8) {
1870                      type=1;
1871                    } else if (numberRows<50) {
1872                      type=3;
1873                    } 
1874                    if (doMask&&!maskMatches(maskStarts,masks,rowNames[iRow]))
1875                      type=0;
1876                    if (type) {
1877                      fprintf(fp,"%7d ",iRow);
1878                      if (lengthName) {
1879                        const char * name = rowNames[iRow].c_str();
1880                        int n=strlen(name);
1881                        int i;
1882                        for (i=0;i<n;i++)
1883                          fprintf(fp,"%c",name[i]);
1884                        for (;i<lengthPrint;i++)
1885                          fprintf(fp," ");
1886                      }
1887                      fprintf(fp," %15.8g        %15.8g\n",primalRowSolution[iRow],
1888                              dualRowSolution[iRow]);
1889                    }
1890                  }
1891                }
1892                int iColumn;
1893                int numberColumns=models[iModel].numberColumns();
1894                double * dualColumnSolution = 
1895  models[iModel].dualColumnSolution();
1896                double * primalColumnSolution = 
1897  models[iModel].primalColumnSolution();
1898                double * columnLower = models[iModel].columnLower();
1899                double * columnUpper = models[iModel].columnUpper();
1900                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1901                  int type=(printMode>3) ? 1 : 0;
1902                  if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
1903                      primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
1904                    fprintf(fp,"** ");
1905                    type=2;
1906                  } else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
1907                    type=1;
1908                  } else if (numberColumns<50) {
1909                    type=3;
1910                  }
1911                  if (doMask&&!maskMatches(maskStarts,masks,
1912                                           columnNames[iColumn]))
1913                    type =0;
1914                  if (type) {
1915                    fprintf(fp,"%7d ",iColumn);
1916                    if (lengthName) {
1917                      const char * name = columnNames[iColumn].c_str();
1918                      int n=strlen(name);
1919                      int i;
1920                      for (i=0;i<n;i++)
1921                        fprintf(fp,"%c",name[i]);
1922                      for (;i<lengthPrint;i++)
1923                        fprintf(fp," ");
1924                    }
1925                    fprintf(fp," %15.8g        %15.8g\n",
1926                            primalColumnSolution[iColumn],
1927                            dualColumnSolution[iColumn]);
1928                  }
1929                }
1930                if (fp!=stdout)
1931                  fclose(fp);
1932                if (masks) {
1933                  delete [] maskStarts;
1934                  for (int i=0;i<maxMasks;i++)
1935                    delete [] masks[i];
1936                  delete [] masks;
1937                }
1938              } else {
1939                std::cout<<"Unable to open file "<<fileName<<std::endl;
1940              }
1941            } else {
1942              std::cout<<"** Current model not valid"<<std::endl;
1943             
1944            }
1945         
1946            break;
1947          case SAVESOL:
1948            if (goodModels[iModel]) {
1949              // get next field
1950              field = CoinReadGetString(argc,argv);
1951              if (field=="$") {
1952                field = parameters[iParam].stringValue();
1953              } else if (field=="EOL") {
1954                parameters[iParam].printString();
1955                break;
1956              } else {
1957                parameters[iParam].setStringValue(field);
1958              }
1959              std::string fileName;
1960              if (field[0]=='/'||field[0]=='\\') {
1961                fileName = field;
1962              } else if (field[0]=='~') {
1963                char * environVar = getenv("HOME");
1964                if (environVar) {
1965                  std::string home(environVar);
1966                  field=field.erase(0,1);
1967                  fileName = home+field;
1968                } else {
1969                  fileName=field;
1970                }
1971              } else {
1972                fileName = directory+field;
1973              }
1974              saveSolution(models+iModel,fileName);
1975            } else {
1976              std::cout<<"** Current model not valid"<<std::endl;
1977             
1978            }
1979            break;
1980          case ENVIRONMENT:
1981            CbcOrClpEnvironmentIndex=0;
1982            break;
1983          default:
1984            abort();
1985          }
1986        } 
1987      } else if (!numberMatches) {
1988        std::cout<<"No match for "<<field<<" - ? for list of commands"
1989                 <<std::endl;
1990      } else if (numberMatches==1) {
1991        if (!numberQuery) {
1992          std::cout<<"Short match for "<<field<<" - completion: ";
1993          std::cout<<parameters[firstMatch].matchName()<<std::endl;
1994        } else if (numberQuery) {
1995          std::cout<<parameters[firstMatch].matchName()<<" : ";
1996          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
1997          if (numberQuery>=2) 
1998            parameters[firstMatch].printLongHelp();
1999        }
2000      } else {
2001        if (!numberQuery) 
2002          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
2003                   <<std::endl;
2004        else
2005          std::cout<<"Completions of "<<field<<":"<<std::endl;
2006        for ( iParam=0; iParam<numberParameters; iParam++ ) {
2007          int match = parameters[iParam].matches(field);
2008          if (match&&parameters[iParam].displayThis()) {
2009            std::cout<<parameters[iParam].matchName();
2010            if (numberQuery>=2) 
2011              std::cout<<" : "<<parameters[iParam].shortHelp();
2012            std::cout<<std::endl;
2013          }
2014        }
2015      }
2016    }
2017    delete [] models;
2018    delete [] goodModels;
2019  }
2020  // By now all memory should be freed
2021#ifdef DMALLOC
2022  dmalloc_log_unfreed();
2023  dmalloc_shutdown();
2024#endif
2025  return 0;
2026}   
2027static void breakdown(const char * name, int numberLook, const double * region)
2028{
2029  double range[] = {
2030    -COIN_DBL_MAX,
2031    -1.0e15,-1.0e11,-1.0e8,-1.0e5,-1.0e4,-1.0e3,-1.0e2,-1.0e1,
2032    -1.0,
2033    -1.0e-1,-1.0e-2,-1.0e-3,-1.0e-4,-1.0e-5,-1.0e-8,-1.0e-11,-1.0e-15,
2034    0.0,
2035    1.0e-15,1.0e-11,1.0e-8,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,
2036    1.0,
2037    1.0e1,1.0e2,1.0e3,1.0e4,1.0e5,1.0e8,1.0e11,1.0e15,
2038    COIN_DBL_MAX};
2039  int nRanges = static_cast<int> (sizeof(range)/sizeof(double));
2040  int * number = new int[nRanges];
2041  memset(number,0,nRanges*sizeof(int));
2042  int * numberExact = new int[nRanges];
2043  memset(numberExact,0,nRanges*sizeof(int));
2044  int i;
2045  for ( i=0;i<numberLook;i++) {
2046    double value = region[i];
2047    for (int j=0;j<nRanges;j++) {
2048      if (value==range[j]) {
2049        numberExact[j]++;
2050        break;
2051      } else if (value<range[j]) {
2052        number[j]++;
2053        break;
2054      }
2055    }
2056  }
2057  printf("\n%s has %d entries\n",name,numberLook);
2058  for (i=0;i<nRanges;i++) {
2059    if (number[i]) 
2060      printf("%d between %g and %g",number[i],range[i-1],range[i]);
2061    if (numberExact[i]) {
2062      if (number[i])
2063        printf(", ");
2064      printf("%d exactly at %g",numberExact[i],range[i]);
2065    }
2066    if (number[i]+numberExact[i])
2067      printf("\n");
2068  }
2069  delete [] number;
2070  delete [] numberExact;
2071}
2072void sortOnOther(int * column,
2073                  const CoinBigIndex * rowStart,
2074                  int * order,
2075                  int * other,
2076                  int nRow,
2077                  int nInRow,
2078                  int where)
2079{
2080  if (nRow<2||where>=nInRow)
2081    return;
2082  // do initial sort
2083  int kRow;
2084  int iRow;
2085  for ( kRow=0;kRow<nRow;kRow++) {
2086    iRow = order[kRow];
2087    other[kRow]=column[rowStart[iRow]+where];
2088  }
2089  CoinSort_2(other,other+nRow,order);
2090  int first=0;
2091  iRow=order[0];
2092  int firstC=column[rowStart[iRow]+where];
2093  kRow=1;
2094  while (kRow<nRow) {
2095    int lastC=9999999;;
2096    for (;kRow<nRow+1;kRow++) {
2097      if (kRow<nRow) {
2098        iRow=order[kRow];
2099        lastC=column[rowStart[iRow]+where];
2100      } else {
2101        lastC=9999999;
2102      }
2103      if (lastC>firstC) 
2104        break;
2105    }
2106    // sort
2107    sortOnOther(column,rowStart,order+first,other,kRow-first,
2108                 nInRow,where+1);
2109    firstC=lastC;
2110    first=kRow;
2111  }
2112}
2113static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
2114{
2115  int numberColumns = originalModel->numberColumns();
2116  const char * integerInformation  = originalModel->integerInformation(); 
2117  const double * columnLower = originalModel->columnLower();
2118  const double * columnUpper = originalModel->columnUpper();
2119  int numberIntegers=0;
2120  int numberBinary=0;
2121  int iRow,iColumn;
2122  if (integerInformation) {
2123    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2124      if (integerInformation[iColumn]) {
2125        if (columnUpper[iColumn]>columnLower[iColumn]) {
2126          numberIntegers++;
2127          if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==1) 
2128            numberBinary++;
2129        }
2130      }
2131    }
2132  }
2133  numberColumns = model->numberColumns();
2134  int numberRows = model->numberRows();
2135  columnLower = model->columnLower();
2136  columnUpper = model->columnUpper();
2137  const double * rowLower = model->rowLower();
2138  const double * rowUpper = model->rowUpper();
2139  const double * objective = model->objective();
2140  CoinPackedMatrix * matrix = model->matrix();
2141  CoinBigIndex numberElements = matrix->getNumElements();
2142  const int * columnLength = matrix->getVectorLengths();
2143  //const CoinBigIndex * columnStart = matrix->getVectorStarts();
2144  const double * elementByColumn = matrix->getElements();
2145  int * number = new int[numberRows+1];
2146  memset(number,0,(numberRows+1)*sizeof(int));
2147  int numberObjSingletons=0;
2148  /* cType
2149     0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
2150     8 0/1
2151  */ 
2152  int cType[9];
2153  std::string cName[]={"0.0->inf,","0.0->up,","lo->inf,","lo->up,","free,","fixed,","-inf->0.0,",
2154                       "-inf->up,","0.0->1.0"};
2155  int nObjective=0;
2156  memset(cType,0,sizeof(cType));
2157  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2158    int length=columnLength[iColumn];
2159    if (length==1&&objective[iColumn])
2160      numberObjSingletons++;
2161    number[length]++;
2162    if (objective[iColumn])
2163      nObjective++;
2164    if (columnLower[iColumn]>-1.0e20) {
2165      if (columnLower[iColumn]==0.0) {
2166        if (columnUpper[iColumn]>1.0e20)
2167          cType[0]++;
2168        else if (columnUpper[iColumn]==1.0)
2169          cType[8]++;
2170        else if (columnUpper[iColumn]==0.0)
2171          cType[5]++;
2172        else
2173          cType[1]++;
2174      } else {
2175        if (columnUpper[iColumn]>1.0e20) 
2176          cType[2]++;
2177        else if (columnUpper[iColumn]==columnLower[iColumn])
2178          cType[5]++;
2179        else
2180          cType[3]++;
2181      }
2182    } else {
2183      if (columnUpper[iColumn]>1.0e20) 
2184        cType[4]++;
2185      else if (columnUpper[iColumn]==0.0) 
2186        cType[6]++;
2187      else
2188        cType[7]++;
2189    }
2190  }
2191  /* rType
2192     0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
2193     7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
2194  */ 
2195  int rType[13];
2196  std::string rName[]={"E 0.0,","E 1.0,","E -1.0,","E other,","G 0.0,","G 1.0,","G other,",
2197                       "L 0.0,","L 1.0,","L other,","Range 0.0->1.0,","Range other,","Free"};
2198  memset(rType,0,sizeof(rType));
2199  for (iRow=0;iRow<numberRows;iRow++) {
2200    if (rowLower[iRow]>-1.0e20) {
2201      if (rowLower[iRow]==0.0) {
2202        if (rowUpper[iRow]>1.0e20)
2203          rType[4]++;
2204        else if (rowUpper[iRow]==1.0)
2205          rType[10]++;
2206        else if (rowUpper[iRow]==0.0)
2207          rType[0]++;
2208        else
2209          rType[11]++;
2210      } else if (rowLower[iRow]==1.0) {
2211        if (rowUpper[iRow]>1.0e20) 
2212          rType[5]++;
2213        else if (rowUpper[iRow]==rowLower[iRow])
2214          rType[1]++;
2215        else
2216          rType[11]++;
2217      } else if (rowLower[iRow]==-1.0) {
2218        if (rowUpper[iRow]>1.0e20) 
2219          rType[6]++;
2220        else if (rowUpper[iRow]==rowLower[iRow])
2221          rType[2]++;
2222        else
2223          rType[11]++;
2224      } else {
2225        if (rowUpper[iRow]>1.0e20) 
2226          rType[6]++;
2227        else if (rowUpper[iRow]==rowLower[iRow])
2228          rType[3]++;
2229        else
2230          rType[11]++;
2231      }
2232    } else {
2233      if (rowUpper[iRow]>1.0e20) 
2234        rType[12]++;
2235      else if (rowUpper[iRow]==0.0) 
2236        rType[7]++;
2237      else if (rowUpper[iRow]==1.0) 
2238        rType[8]++;
2239      else
2240        rType[9]++;
2241    }
2242  }
2243  // Basic statistics
2244  printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
2245         numberRows,numberColumns,nObjective,numberElements);
2246  if (number[0]+number[1]) {
2247    printf("There are ");
2248    if (numberObjSingletons)
2249      printf("%d singletons with objective ",numberObjSingletons);
2250    int numberNoObj = number[1]-numberObjSingletons;
2251    if (numberNoObj)
2252      printf("%d singletons with no objective ",numberNoObj);
2253    if (number[0])
2254      printf("** %d columns have no entries",number[0]);
2255    printf("\n");
2256  }
2257  printf("Column breakdown:\n");
2258  int k;
2259  for (k=0;k<static_cast<int> (sizeof(cType)/sizeof(int));k++) {
2260    printf("%d of type %s ",cType[k],cName[k].c_str());
2261    if (((k+1)%3)==0)
2262      printf("\n");
2263  }
2264  if ((k%3)!=0)
2265    printf("\n");
2266  printf("Row breakdown:\n");
2267  for (k=0;k<static_cast<int> (sizeof(rType)/sizeof(int));k++) {
2268    printf("%d of type %s ",rType[k],rName[k].c_str());
2269    if (((k+1)%3)==0)
2270      printf("\n");
2271  }
2272  if ((k%3)!=0)
2273    printf("\n");
2274#define SYM
2275#ifndef SYM
2276  if (model->logLevel()<2)
2277    return ;
2278#endif
2279  int kMax = model->logLevel()>3 ? 1000000 : 10;
2280  k=0;
2281  for (iRow=1;iRow<=numberRows;iRow++) {
2282    if (number[iRow]) {
2283      k++;
2284      printf("%d columns have %d entries\n",number[iRow],iRow);
2285      if (k==kMax)
2286        break;
2287    }
2288  }
2289  if (k<numberRows) {
2290    int kk=k;
2291    k=0;
2292    for (iRow=numberRows;iRow>=1;iRow--) {
2293      if (number[iRow]) {
2294        k++;
2295        if (k==kMax)
2296          break;
2297      }
2298    }
2299    if (k>kk) {
2300      printf("\n    .........\n\n");
2301      iRow=k;
2302      k=0;
2303      for (;iRow<numberRows;iRow++) {
2304        if (number[iRow]) {
2305          k++;
2306          printf("%d columns have %d entries\n",number[iRow],iRow);
2307          if (k==kMax)
2308            break;
2309        }
2310      }
2311    }
2312  }
2313  delete [] number;
2314  printf("\n\n");
2315  if (model->logLevel()==63
2316#ifdef SYM
2317      ||true
2318#endif
2319      ) {
2320    // get column copy
2321    CoinPackedMatrix columnCopy = *matrix;
2322    const int * columnLength = columnCopy.getVectorLengths();
2323    number = new int[numberRows+1];
2324    memset(number,0,(numberRows+1)*sizeof(int));
2325    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2326      int length=columnLength[iColumn];
2327      number[length]++;
2328    }
2329    k=0;
2330    for (iRow=1;iRow<=numberRows;iRow++) {
2331      if (number[iRow]) {
2332        k++;
2333      }
2334    }
2335    int * row = columnCopy.getMutableIndices();
2336    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2337    double * element = columnCopy.getMutableElements();
2338    int * order = new int[numberColumns];
2339    int * other = new int[numberColumns];
2340    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2341      int length=columnLength[iColumn];
2342      order[iColumn]=iColumn;
2343      other[iColumn]=length;
2344      CoinBigIndex start = columnStart[iColumn];
2345      CoinSort_2(row+start,row+start+length,element+start);
2346    }
2347    CoinSort_2(other,other+numberColumns,order);
2348    int jColumn=number[0]+number[1];
2349    for (iRow=2;iRow<=numberRows;iRow++) {
2350      if (number[iRow]) {
2351        printf("XX %d columns have %d entries\n",number[iRow],iRow);
2352        int kColumn = jColumn+number[iRow];
2353        sortOnOther(row,columnStart,
2354                     order+jColumn,other,number[iRow],iRow,0);
2355        // Now print etc
2356        if (iRow<500000) {
2357          for (int lColumn =jColumn;lColumn<kColumn;lColumn++) {
2358            iColumn = order[lColumn];
2359            CoinBigIndex start = columnStart[iColumn];
2360            if (model->logLevel()==63) {
2361              printf("column %d %g <= ",iColumn,columnLower[iColumn]);
2362              for (CoinBigIndex i=start;i<start+iRow;i++)
2363                printf("( %d, %g) ",row[i],element[i]);
2364              printf("<= %g\n",columnUpper[iColumn]);
2365            }
2366          }
2367        }
2368        jColumn =kColumn;
2369      }
2370    }
2371    delete [] order;
2372    delete [] other;
2373    delete [] number;
2374  }
2375  // get row copy
2376  CoinPackedMatrix rowCopy = *matrix;
2377  rowCopy.reverseOrdering();
2378  const int * rowLength = rowCopy.getVectorLengths();
2379  number = new int[numberColumns+1];
2380  memset(number,0,(numberColumns+1)*sizeof(int));
2381  for (iRow=0;iRow<numberRows;iRow++) {
2382    int length=rowLength[iRow];
2383    number[length]++;
2384  }
2385  if (number[0])
2386    printf("** %d rows have no entries\n",number[0]);
2387  k=0;
2388  for (iColumn=1;iColumn<=numberColumns;iColumn++) {
2389    if (number[iColumn]) {
2390      k++;
2391      printf("%d rows have %d entries\n",number[iColumn],iColumn);
2392      if (k==kMax)
2393        break;
2394    }
2395  }
2396  if (k<numberColumns) {
2397    int kk=k;
2398    k=0;
2399    for (iColumn=numberColumns;iColumn>=1;iColumn--) {
2400      if (number[iColumn]) {
2401        k++;
2402        if (k==kMax)
2403          break;
2404      }
2405    }
2406    if (k>kk) {
2407      printf("\n    .........\n\n");
2408      iColumn=k;
2409      k=0;
2410      for (;iColumn<numberColumns;iColumn++) {
2411        if (number[iColumn]) {
2412          k++;
2413          printf("%d rows have %d entries\n",number[iColumn],iColumn);
2414          if (k==kMax)
2415            break;
2416        }
2417      }
2418    }
2419  }
2420  if (model->logLevel()==63
2421#ifdef SYM
2422      ||true
2423#endif
2424      ) {
2425    int * column = rowCopy.getMutableIndices();
2426    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
2427    double * element = rowCopy.getMutableElements();
2428    int * order = new int[numberRows];
2429    int * other = new int[numberRows];
2430    for (iRow=0;iRow<numberRows;iRow++) {
2431      int length=rowLength[iRow];
2432      order[iRow]=iRow;
2433      other[iRow]=length;
2434      CoinBigIndex start = rowStart[iRow];
2435      CoinSort_2(column+start,column+start+length,element+start);
2436    }
2437    CoinSort_2(other,other+numberRows,order);
2438    int jRow=number[0]+number[1];
2439    double * weight = new double[numberRows];
2440    double * randomColumn = new double[numberColumns+1];
2441    double * randomRow = new double [numberRows+1];
2442    int * sortRow = new int [numberRows];
2443    int * possibleRow = new int [numberRows];
2444    int * backRow = new int [numberRows];
2445    int * stackRow = new int [numberRows];
2446    int * sortColumn = new int [numberColumns];
2447    int * possibleColumn = new int [numberColumns];
2448    int * backColumn = new int [numberColumns];
2449    int * backColumn2 = new int [numberColumns];
2450    int * mapRow = new int [numberRows];
2451    int * mapColumn = new int [numberColumns];
2452    int * stackColumn = new int [numberColumns];
2453    double randomLower = CoinDrand48();
2454    double randomUpper = CoinDrand48();
2455    double randomInteger = CoinDrand48();
2456    int * startAdd = new int[numberRows+1];
2457    int * columnAdd = new int [2*numberElements];
2458    double * elementAdd = new double[2*numberElements];
2459    int nAddRows=0;
2460    startAdd[0]=0;
2461    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2462      randomColumn[iColumn] = CoinDrand48();
2463      backColumn2[iColumn]=-1;
2464    }
2465    for (iColumn=2;iColumn<=numberColumns;iColumn++) {
2466      if (number[iColumn]) {
2467        printf("XX %d rows have %d entries\n",number[iColumn],iColumn);
2468        int kRow = jRow+number[iColumn];
2469        sortOnOther(column,rowStart,
2470                     order+jRow,other,number[iColumn],iColumn,0);
2471        // Now print etc
2472        if (iColumn<500000) {
2473          int nLook=0;
2474          for (int lRow =jRow;lRow<kRow;lRow++) {
2475            iRow = order[lRow];
2476            CoinBigIndex start = rowStart[iRow];
2477            if (model->logLevel()==63) {
2478              printf("row %d %g <= ",iRow,rowLower[iRow]);
2479              for (CoinBigIndex i=start;i<start+iColumn;i++) 
2480                printf("( %d, %g) ",column[i],element[i]);
2481              printf("<= %g\n",rowUpper[iRow]);
2482            }
2483            int first = column[start];
2484            double sum=0.0;
2485            for (CoinBigIndex i=start;i<start+iColumn;i++) {
2486              int jColumn = column[i];
2487              double value = element[i];
2488              jColumn -= first;
2489              assert (jColumn>=0);
2490              sum += value*randomColumn[jColumn];
2491            }
2492            if (rowLower[iRow]>-1.0e30&&rowLower[iRow])
2493              sum += rowLower[iRow]*randomLower;
2494            else if (!rowLower[iRow])
2495              sum += 1.234567e-7*randomLower;
2496            if (rowUpper[iRow]<1.0e30&&rowUpper[iRow])
2497              sum += rowUpper[iRow]*randomUpper;
2498            else if (!rowUpper[iRow])
2499              sum += 1.234567e-7*randomUpper;
2500            sortRow[nLook]=iRow;
2501            randomRow[nLook++]=sum;
2502            // best way is to number unique elements and bounds and use
2503            if (fabs(sum)>1.0e4)
2504              sum *= 1.0e-6;
2505            weight[iRow]=sum;
2506          }
2507          assert (nLook<=numberRows);
2508          CoinSort_2(randomRow,randomRow+nLook,sortRow);
2509          randomRow[nLook]=COIN_DBL_MAX;
2510          double last=-COIN_DBL_MAX;
2511          int iLast=-1;
2512          for (int iLook=0;iLook<nLook+1;iLook++) {
2513            if (randomRow[iLook]>last) {
2514              if (iLast>=0) {
2515                int n=iLook-iLast;
2516                if (n>1) {
2517                  //printf("%d rows possible?\n",n);
2518                }
2519              }
2520              iLast=iLook;
2521              last=randomRow[iLook];
2522            }
2523          }
2524        }
2525        jRow =kRow;
2526      }
2527    }
2528    CoinPackedMatrix columnCopy = *matrix;
2529    const int * columnLength = columnCopy.getVectorLengths();
2530    const int * row = columnCopy.getIndices();
2531    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2532    const double * elementByColumn = columnCopy.getElements();
2533    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2534      int length=columnLength[iColumn];
2535      CoinBigIndex start = columnStart[iColumn];
2536      double sum = objective[iColumn];
2537      if (columnLower[iColumn]>-1.0e30&&columnLower[iColumn])
2538        sum += columnLower[iColumn]*randomLower;
2539      else if (!columnLower[iColumn])
2540        sum += 1.234567e-7*randomLower;
2541      if (columnUpper[iColumn]<1.0e30&&columnUpper[iColumn])
2542        sum += columnUpper[iColumn]*randomUpper;
2543      else if (!columnUpper[iColumn])
2544        sum += 1.234567e-7*randomUpper;
2545      if (model->isInteger(iColumn))
2546        sum += 9.87654321e-6*randomInteger;
2547      for (CoinBigIndex i=start;i<start+length;i++) {
2548        int iRow = row[i];
2549        sum += elementByColumn[i]*weight[iRow];
2550      }
2551      sortColumn[iColumn]=iColumn;
2552      randomColumn[iColumn]=sum;
2553    }
2554    {
2555      CoinSort_2(randomColumn,randomColumn+numberColumns,sortColumn);
2556      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2557        int i=sortColumn[iColumn];
2558        backColumn[i]=iColumn;
2559      }
2560      randomColumn[numberColumns]=COIN_DBL_MAX;
2561      double last=-COIN_DBL_MAX;
2562      int iLast=-1;
2563      for (int iLook=0;iLook<numberColumns+1;iLook++) {
2564        if (randomColumn[iLook]>last) {
2565          if (iLast>=0) {
2566            int n=iLook-iLast;
2567            if (n>1) {
2568              //printf("%d columns possible?\n",n);
2569            }
2570            for (int i=iLast;i<iLook;i++) {
2571              possibleColumn[sortColumn[i]]=n;
2572            }
2573          }
2574          iLast=iLook;
2575          last=randomColumn[iLook];
2576        }
2577      }
2578      for (iRow =0;iRow<numberRows;iRow++) {
2579        CoinBigIndex start = rowStart[iRow];
2580        double sum=0.0;
2581        int length=rowLength[iRow];
2582        for (CoinBigIndex i=start;i<start+length;i++) {
2583          int jColumn = column[i];
2584          double value = element[i];
2585          jColumn=backColumn[jColumn];
2586          sum += value*randomColumn[jColumn];
2587          //if (iColumn==23089||iRow==23729)
2588          //printf("row %d cola %d colb %d value %g rand %g sum %g\n",
2589          //   iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
2590        }
2591        sortRow[iRow]=iRow;
2592        randomRow[iRow]=weight[iRow];
2593        randomRow[iRow]=sum;
2594      }
2595      CoinSort_2(randomRow,randomRow+numberRows,sortRow);
2596      for (iRow=0;iRow<numberRows;iRow++) {
2597        int i=sortRow[iRow];
2598        backRow[i]=iRow;
2599      }
2600      randomRow[numberRows]=COIN_DBL_MAX;
2601      last=-COIN_DBL_MAX;
2602      iLast=-1;
2603      // Do backward indices from order
2604      for (iRow=0;iRow<numberRows;iRow++) {
2605        other[order[iRow]]=iRow;
2606      }
2607      for (int iLook=0;iLook<numberRows+1;iLook++) {
2608        if (randomRow[iLook]>last) {
2609          if (iLast>=0) {
2610            int n=iLook-iLast;
2611            if (n>1) {
2612              //printf("%d rows possible?\n",n);
2613              // Within group sort as for original "order"
2614              for (int i=iLast;i<iLook;i++) {
2615                int jRow=sortRow[i];
2616                order[i]=other[jRow];
2617              }
2618              CoinSort_2(order+iLast,order+iLook,sortRow+iLast);
2619            }
2620            for (int i=iLast;i<iLook;i++) {
2621              possibleRow[sortRow[i]]=n;
2622            }
2623          }
2624          iLast=iLook;
2625          last=randomRow[iLook];
2626        }
2627      }
2628      // Temp out
2629      for (int iLook=0;iLook<numberRows-1000000;iLook++) {
2630        iRow=sortRow[iLook];
2631        CoinBigIndex start = rowStart[iRow];
2632        int length=rowLength[iRow];
2633        int numberPossible = possibleRow[iRow];
2634        for (CoinBigIndex i=start;i<start+length;i++) {
2635          int jColumn = column[i];
2636          if (possibleColumn[jColumn]!=numberPossible)
2637            numberPossible=-1;
2638        }
2639        int n=numberPossible;
2640        if (numberPossible>1) {
2641          //printf("pppppossible %d\n",numberPossible);
2642          for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2643            int jRow=sortRow[jLook];
2644            CoinBigIndex start2 = rowStart[jRow];
2645            assert (numberPossible==possibleRow[jRow]);
2646            assert(length==rowLength[jRow]);
2647            for (CoinBigIndex i=start2;i<start2+length;i++) {
2648              int jColumn = column[i];
2649              if (possibleColumn[jColumn]!=numberPossible)
2650                numberPossible=-1;
2651            }
2652          }
2653          if (numberPossible<2) {
2654            // switch off
2655            for (int jLook=iLook;jLook<iLook+n;jLook++) 
2656              possibleRow[sortRow[jLook]]=-1;
2657          }
2658          // skip rest
2659          iLook += n-1;
2660        } else {
2661          possibleRow[iRow]=-1;
2662        }
2663      }
2664      for (int iLook=0;iLook<numberRows;iLook++) {
2665        iRow=sortRow[iLook];
2666        int numberPossible = possibleRow[iRow];
2667        // Only if any integers
2668        int numberIntegers=0;
2669        CoinBigIndex start = rowStart[iRow];
2670        int length=rowLength[iRow];
2671        for (CoinBigIndex i=start;i<start+length;i++) {
2672          int jColumn = column[i];
2673          if (model->isInteger(jColumn))
2674            numberIntegers++;
2675        }
2676        if (numberPossible>1&&!numberIntegers) {
2677          //printf("possible %d - but no integers\n",numberPossible);
2678        }
2679        if (numberPossible>1&&(numberIntegers||false)) {
2680          //
2681          printf("possible %d - %d integers\n",numberPossible,numberIntegers);
2682          int lastLook=iLook;
2683          int nMapRow=-1;
2684          for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2685            // stop if too many failures
2686            if (jLook>iLook+10&&nMapRow<0)
2687              break;
2688            // Create identity mapping
2689            int i;
2690            for (i=0;i<numberRows;i++)
2691              mapRow[i]=i;
2692            for (i=0;i<numberColumns;i++)
2693              mapColumn[i]=i;
2694            int offset=jLook-iLook;
2695            int nStackC=0;
2696            // build up row and column mapping
2697            int nStackR=1;
2698            stackRow[0]=iLook;
2699            bool good=true;
2700            while (nStackR) {
2701              nStackR--;
2702              int look1 = stackRow[nStackR];
2703              int look2 = look1 + offset;
2704              assert (randomRow[look1]==randomRow[look2]);
2705              int row1=sortRow[look1];
2706              int row2=sortRow[look2];
2707              assert (mapRow[row1]==row1);
2708              assert (mapRow[row2]==row2);
2709              mapRow[row1]=row2;
2710              mapRow[row2]=row1;
2711              CoinBigIndex start1 = rowStart[row1];
2712              CoinBigIndex offset2 = rowStart[row2]-start1;
2713              int length=rowLength[row1];
2714              assert( length==rowLength[row2]);
2715              for (CoinBigIndex i=start1;i<start1+length;i++) {
2716                int jColumn1 = column[i];
2717                int jColumn2 = column[i+offset2];
2718                if (randomColumn[backColumn[jColumn1]]!=
2719                    randomColumn[backColumn[jColumn2]]) {
2720                  good=false;
2721                  break;
2722                }
2723                if (mapColumn[jColumn1]==jColumn1) {
2724                  // not touched
2725                  assert (mapColumn[jColumn2]==jColumn2);
2726                  if (jColumn1!=jColumn2) {
2727                    // Put on stack
2728                    mapColumn[jColumn1]=jColumn2;
2729                    mapColumn[jColumn2]=jColumn1;
2730                    stackColumn[nStackC++]=jColumn1;
2731                  }
2732                } else {
2733                  if(mapColumn[jColumn1]!=jColumn2||
2734                     mapColumn[jColumn2]!=jColumn1) {
2735                    // bad
2736                    good=false;
2737                    printf("bad col\n");
2738                    break;
2739                  }
2740                }
2741              }
2742              if (!good)
2743                break;
2744              while (nStackC) {
2745                nStackC--;
2746                int iColumn = stackColumn[nStackC];
2747                int iColumn2 = mapColumn[iColumn];
2748                assert (iColumn!=iColumn2);
2749                int length=columnLength[iColumn];
2750                assert (length==columnLength[iColumn2]);
2751                CoinBigIndex start = columnStart[iColumn];
2752                CoinBigIndex offset2 = columnStart[iColumn2]-start;
2753                for (CoinBigIndex i=start;i<start+length;i++) {
2754                  int iRow = row[i];
2755                  int iRow2 = row[i+offset2];
2756                  if (mapRow[iRow]==iRow) {
2757                    // First (but be careful)
2758                    if (iRow!=iRow2) {
2759                      //mapRow[iRow]=iRow2;
2760                      //mapRow[iRow2]=iRow;
2761                      int iBack=backRow[iRow];
2762                      int iBack2=backRow[iRow2];
2763                      if (randomRow[iBack]==randomRow[iBack2]&&
2764                          iBack2-iBack==offset) {
2765                        stackRow[nStackR++]=iBack;
2766                      } else {
2767                        //printf("randomRow diff - weights %g %g\n",
2768                        //     weight[iRow],weight[iRow2]);
2769                        // bad
2770                        good=false;
2771                        break;
2772                      }
2773                    }
2774                  } else {
2775                    if(mapRow[iRow]!=iRow2||
2776                       mapRow[iRow2]!=iRow) {
2777                      // bad
2778                      good=false;
2779                      printf("bad row\n");
2780                      break;
2781                    }
2782                  }
2783                }
2784                if (!good)
2785                  break;
2786              }
2787            }
2788            // then check OK
2789            if (good) {
2790              for (iRow =0;iRow<numberRows;iRow++) {
2791                CoinBigIndex start = rowStart[iRow];
2792                int length=rowLength[iRow];
2793                if (mapRow[iRow]==iRow) {
2794                  for (CoinBigIndex i=start;i<start+length;i++) {
2795                    int jColumn = column[i];
2796                    backColumn2[jColumn]=i-start;
2797                  }
2798                  for (CoinBigIndex i=start;i<start+length;i++) {
2799                    int jColumn = column[i];
2800                    if (mapColumn[jColumn]!=jColumn) {
2801                      int jColumn2 =mapColumn[jColumn];
2802                      CoinBigIndex i2 = backColumn2[jColumn2];
2803                      if (i2<0) {
2804                        good=false;
2805                      } else if (element[i]!=element[i2+start]) {
2806                        good=false;
2807                      }
2808                    }
2809                  }
2810                  for (CoinBigIndex i=start;i<start+length;i++) {
2811                    int jColumn = column[i];
2812                    backColumn2[jColumn]=-1;
2813                  }
2814                } else {
2815                  int row2 = mapRow[iRow];
2816                  assert (iRow = mapRow[row2]);
2817                  if (rowLower[iRow]!=rowLower[row2]||
2818                      rowLower[row2]!=rowLower[iRow])
2819                    good=false;
2820                  CoinBigIndex offset2 = rowStart[row2]-start;
2821                  for (CoinBigIndex i=start;i<start+length;i++) {
2822                    int jColumn = column[i];
2823                    double value = element[i];
2824                    int jColumn2 = column[i+offset2];
2825                    double value2 = element[i+offset2];
2826                    if (value!=value2||mapColumn[jColumn]!=jColumn2||
2827                        mapColumn[jColumn2]!=jColumn)
2828                      good=false;
2829                  }
2830                }
2831              }
2832              if (good) {
2833                // check rim
2834                for (iColumn=0;iColumn<numberColumns;iColumn++) {
2835                  if (mapColumn[iColumn]!=iColumn) {
2836                    int iColumn2 = mapColumn[iColumn];
2837                    if(objective[iColumn]!=objective[iColumn2])
2838                      good=false;
2839                    if (columnLower[iColumn]!=columnLower[iColumn2])
2840                      good=false;
2841                    if (columnUpper[iColumn]!=columnUpper[iColumn2])
2842                      good=false;
2843                    if (model->isInteger(iColumn)!=model->isInteger(iColumn2))
2844                      good=false;
2845                  }
2846                }
2847              }
2848              if (good) {
2849                // temp
2850                if (nMapRow<0) {
2851                  //const double * solution = model->primalColumnSolution();
2852                  // find mapped
2853                  int nMapColumn=0;
2854                  for (int i=0;i<numberColumns;i++) {
2855                    if (mapColumn[i]>i) 
2856                      nMapColumn++;
2857                  }
2858                  nMapRow=0;
2859                  int kRow=-1;
2860                  for (int i=0;i<numberRows;i++) {
2861                    if (mapRow[i]>i) { 
2862                      nMapRow++;
2863                      kRow=i;
2864                    }
2865                  }
2866                  printf("%d columns, %d rows\n",nMapColumn,nMapRow);
2867                  if (nMapRow==1) {
2868                    CoinBigIndex start = rowStart[kRow];
2869                    int length=rowLength[kRow];
2870                    printf("%g <= ",rowLower[kRow]);
2871                    for (CoinBigIndex i=start;i<start+length;i++) {
2872                      int jColumn = column[i];
2873                      if (mapColumn[jColumn]!=jColumn) 
2874                        printf("* ");
2875                      printf("%d,%g ",jColumn,element[i]);
2876                    }
2877                    printf("<= %g\n",rowUpper[kRow]);
2878                  }
2879                }
2880                // temp
2881                int row1=sortRow[lastLook];
2882                int row2=sortRow[jLook];
2883                lastLook=jLook;
2884                CoinBigIndex start1 = rowStart[row1];
2885                CoinBigIndex offset2 = rowStart[row2]-start1;
2886                int length=rowLength[row1];
2887                assert( length==rowLength[row2]);
2888                CoinBigIndex put=startAdd[nAddRows];
2889                double multiplier=length<11 ? 2.0 : 1.125;
2890                double value=1.0;
2891                for (CoinBigIndex i=start1;i<start1+length;i++) {
2892                  int jColumn1 = column[i];
2893                  int jColumn2 = column[i+offset2];
2894                  columnAdd[put]=jColumn1;
2895                  elementAdd[put++]=value;
2896                  columnAdd[put]=jColumn2;
2897                  elementAdd[put++]=-value;
2898                  value *= multiplier;
2899                }
2900                nAddRows++;
2901                startAdd[nAddRows]=put;
2902              } else {
2903                printf("ouch - did not check out as good\n");
2904              }
2905            }
2906          }
2907          // skip rest
2908          iLook += numberPossible-1;
2909        }
2910      }
2911    }
2912    if (nAddRows) {
2913      double * lower = new double [nAddRows];
2914      double * upper = new double[nAddRows];
2915      int i;
2916      //const double * solution = model->primalColumnSolution();
2917      for (i=0;i<nAddRows;i++) {
2918        lower[i]=0.0;
2919        upper[i]=COIN_DBL_MAX;
2920      }
2921      printf("Adding %d rows with %d elements\n",nAddRows,
2922             startAdd[nAddRows]);
2923      //ClpSimplex newModel(*model);
2924      //newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
2925      //newModel.writeMps("modified.mps");
2926      delete [] lower;
2927      delete [] upper;
2928    }
2929    delete [] startAdd;
2930    delete [] columnAdd;
2931    delete [] elementAdd;
2932    delete [] order;
2933    delete [] other;
2934    delete [] randomColumn;
2935    delete [] weight;
2936    delete [] randomRow;
2937    delete [] sortRow;
2938    delete [] backRow;
2939    delete [] possibleRow;
2940    delete [] sortColumn;
2941    delete [] backColumn;
2942    delete [] backColumn2;
2943    delete [] possibleColumn;
2944    delete [] mapRow;
2945    delete [] mapColumn;
2946    delete [] stackRow;
2947    delete [] stackColumn;
2948  }
2949  delete [] number;
2950  // Now do breakdown of ranges
2951  breakdown("Elements",numberElements,elementByColumn);
2952  breakdown("RowLower",numberRows,rowLower);
2953  breakdown("RowUpper",numberRows,rowUpper);
2954  breakdown("ColumnLower",numberColumns,columnLower);
2955  breakdown("ColumnUpper",numberColumns,columnUpper);
2956  breakdown("Objective",numberColumns,objective);
2957}
2958static bool maskMatches(const int * starts, char ** masks,
2959                        std::string & check)
2960{
2961  // back to char as I am old fashioned
2962  const char * checkC = check.c_str();
2963  int length = strlen(checkC);
2964  while (checkC[length-1]==' ')
2965    length--;
2966  for (int i=starts[length];i<starts[length+1];i++) {
2967    char * thisMask = masks[i];
2968    int k;
2969    for ( k=0;k<length;k++) {
2970      if (thisMask[k]!='?'&&thisMask[k]!=checkC[k]) 
2971        break;
2972    }
2973    if (k==length)
2974      return true;
2975  }
2976  return false;
2977}
2978static void clean(char * temp)
2979{
2980  char * put = temp;
2981  while (*put>=' ')
2982    put++;
2983  *put='\0';
2984}
2985static void generateCode(const char * fileName,int type)
2986{
2987  FILE * fp = fopen(fileName,"r");
2988  assert (fp);
2989  int numberLines=0;
2990#define MAXLINES 500
2991#define MAXONELINE 200
2992  char line[MAXLINES][MAXONELINE];
2993  while (fgets(line[numberLines],MAXONELINE,fp)) {
2994    assert (numberLines<MAXLINES);
2995    clean(line[numberLines]);
2996    numberLines++;
2997  }
2998  fclose(fp);
2999  // add in actual solve
3000  strcpy(line[numberLines],"5  clpModel->initialSolve(clpSolve);");
3001  numberLines++;
3002  fp = fopen(fileName,"w");
3003  assert (fp);
3004  char apo='"';   
3005  char backslash = '\\';
3006
3007  fprintf(fp,"#include %cClpSimplex.hpp%c\n",apo,apo);
3008  fprintf(fp,"#include %cClpSolve.hpp%c\n",apo,apo);
3009
3010  fprintf(fp,"\nint main (int argc, const char *argv[])\n{\n");
3011  fprintf(fp,"  ClpSimplex  model;\n");
3012  fprintf(fp,"  int status=1;\n");
3013  fprintf(fp,"  if (argc<2)\n");
3014  fprintf(fp,"    fprintf(stderr,%cPlease give file name%cn%c);\n",
3015          apo,backslash,apo);
3016  fprintf(fp,"  else\n");
3017  fprintf(fp,"    status=model.readMps(argv[1],true);\n");
3018  fprintf(fp,"  if (status) {\n");
3019  fprintf(fp,"    fprintf(stderr,%cBad readMps %%s%cn%c,argv[1]);\n",
3020                apo,backslash,apo);
3021  fprintf(fp,"    exit(1);\n");
3022  fprintf(fp,"  }\n\n");
3023  fprintf(fp,"  // Now do requested saves and modifications\n");
3024  fprintf(fp,"  ClpSimplex * clpModel = & model;\n");
3025  int wanted[9];
3026  memset(wanted,0,sizeof(wanted));
3027  wanted[0]=wanted[3]=wanted[5]=wanted[8]=1;
3028  if (type>0) 
3029    wanted[1]=wanted[6]=1;
3030  if (type>1) 
3031    wanted[2]=wanted[4]=wanted[7]=1;
3032  std::string header[9]=
3033  { "","Save values","Redundant save of default values","Set changed values",
3034    "Redundant set default values","Solve","Restore values","Redundant restore values","Add to model"};
3035  for (int iType=0;iType<9;iType++) {
3036    if (!wanted[iType])
3037      continue;
3038    int n=0;
3039    int iLine;
3040    for (iLine=0;iLine<numberLines;iLine++) {
3041      if (line[iLine][0]=='0'+iType) {
3042        if (!n)
3043          fprintf(fp,"\n  // %s\n\n",header[iType].c_str());
3044        n++;
3045        fprintf(fp,"%s\n",line[iLine]+1);
3046      }
3047    }
3048  }
3049  fprintf(fp,"\n  // Now you would use solution etc etc\n\n");
3050  fprintf(fp,"  return 0;\n}\n");
3051  fclose(fp);
3052  printf("C++ file written to %s\n",fileName);
3053}
3054/*
3055  Version 1.00.00 October 13 2004.
3056  1.00.01 October 18.  Added basis handling helped/prodded by Thorsten Koch.
3057  Also modifications to make faster with sbb (I hope I haven't broken anything).
3058  1.00.02 March 21 2005.  Redid ClpNonLinearCost to save memory also redid
3059  createRim to try and improve cache characteristics.
3060  1.00.03 April 8 2005.  Added Volume algorithm as crash and made code more
3061  robust on testing.  Also added "either" and "tune" algorithm.
3062  1.01.01 April 12 2005.  Decided to go to different numbering.  Backups will
3063  be last 2 digits while middle 2 are for improvements.  Still take a long
3064  time to get to 2.00.01
3065  1.01.02 May 4 2005.  Will be putting in many changes - so saving stable version
3066  1.02.01 May 6 2005.  Lots of changes to try and make faster and more stable in
3067  branch and cut.
3068  1.02.02 May 19 2005.  Stuff for strong branching and some improvements to simplex
3069  1.03.01 May 24 2006.  Lots done but I can't remember what!
3070  1.03.03 June 13 2006.  For clean up after dual perturbation
3071  1.04.01 June 26 2007.  Lots of changes but I got lazy
3072  1.05.00 June 27 2007.  This is trunk so when gets to stable will be 1.5
3073 */
Note: See TracBrowser for help on using the repository browser.