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

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

fix for ray and warning messages

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