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

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

Creating new stable branch 1.11 from trunk (rev 1457)

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