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

Last change on this file since 1322 was 1322, checked in by jpfasano, 11 years ago

Change to source so will compile with MS Visual Studio Express V9.
Error message indicated that reinterpret cast can not be used to remove constness.

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