source: stable/1.11/Clp/src/ClpMain.cpp @ 1489

Last change on this file since 1489 was 1489, checked in by forrest, 10 years ago

absolute path on windows

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