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

Last change on this file since 1287 was 1287, checked in by forrest, 13 years ago

changes for cbc event handler and multiple factorizations

  • 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 = ((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=((ClpSimplexOther *) models+iModel)->restoreFromDual(model2);
750                if (model2->status()==3)
751                  returnCode=0;
752                delete model2;
753                if (returnCode&&dualize!=2) {
754                  CoinSighandler_t saveSignal=SIG_DFL;
755                  currentModel = models+iModel;
756                  // register signal handler
757                  saveSignal = signal(SIGINT,signal_handler);
758                  models[iModel].primal(1);
759                  currentModel=NULL;
760                }
761              }
762              if (status>=0)
763                basisHasValues=1;
764            } else {
765              std::cout<<"** Current model not valid"<<std::endl;
766            }
767            break;
768          case STATISTICS:
769            if (goodModels[iModel]) {
770              // If presolve on look at presolved
771              bool deleteModel2=false;
772              ClpSimplex * model2 = models+iModel;
773              if (preSolve) {
774                ClpPresolve pinfo;
775                int presolveOptions2 = presolveOptions&~0x40000000;
776                if ((presolveOptions2&0xffff)!=0)
777                  pinfo.setPresolveActions(presolveOptions2);
778                pinfo.setSubstitution(substitution);
779                if ((printOptions&1)!=0)
780                  pinfo.statistics();
781                double presolveTolerance = 
782                  parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
783                model2 = 
784                  pinfo.presolvedModel(models[iModel],presolveTolerance,
785                                       true,preSolve);
786                if (model2) {
787                  printf("Statistics for presolved model\n");
788                  deleteModel2=true;
789                } else {
790                  printf("Presolved model looks infeasible - will use unpresolved\n");
791                  model2 = models+iModel;
792                }
793              } else {
794                printf("Statistics for unpresolved model\n");
795                model2 =  models+iModel;
796              }
797              statistics(models+iModel,model2);
798              if (deleteModel2)
799                delete model2;
800            } else {
801              std::cout<<"** Current model not valid"<<std::endl;
802            }
803            break;
804          case TIGHTEN:
805            if (goodModels[iModel]) {
806              int numberInfeasibilities = models[iModel].tightenPrimalBounds();
807              if (numberInfeasibilities)
808                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
809            } else {
810              std::cout<<"** Current model not valid"<<std::endl;
811            }
812            break;
813          case PLUSMINUS:
814            if (goodModels[iModel]) {
815              ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
816              ClpPackedMatrix* clpMatrix =
817                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
818              if (clpMatrix) {
819                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
820                if (newMatrix->getIndices()) {
821                  models[iModel].replaceMatrix(newMatrix);
822                  delete saveMatrix;
823                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
824                } else {
825                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
826                }
827              } else {
828                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
829              }
830            } else {
831              std::cout<<"** Current model not valid"<<std::endl;
832            }
833            break;
834          case NETWORK:
835            if (goodModels[iModel]) {
836              ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
837              ClpPackedMatrix* clpMatrix =
838                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
839              if (clpMatrix) {
840                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
841                if (newMatrix->getIndices()) {
842                  models[iModel].replaceMatrix(newMatrix);
843                  delete saveMatrix;
844                  std::cout<<"Matrix converted to network matrix"<<std::endl;
845                } else {
846                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
847                }
848              } else {
849                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
850              }
851            } else {
852              std::cout<<"** Current model not valid"<<std::endl;
853            }
854            break;
855          case IMPORT:
856            {
857              // get next field
858              field = CoinReadGetString(argc,argv);
859              if (field=="$") {
860                field = parameters[iParam].stringValue();
861              } else if (field=="EOL") {
862                parameters[iParam].printString();
863                break;
864              } else {
865                parameters[iParam].setStringValue(field);
866              }
867              std::string fileName;
868              bool canOpen=false;
869              // See if gmpl file
870              int gmpl=0;
871              std::string gmplData;
872              if (field=="-") {
873                // stdin
874                canOpen=true;
875                fileName = "-";
876              } else {
877                // See if .lp
878                {
879                  const char * c_name = field.c_str();
880                  int length = strlen(c_name);
881                  if (length>3&&!strncmp(c_name+length-3,".lp",3))
882                    gmpl=-1; // .lp
883                }
884                bool absolutePath;
885                if (dirsep=='/') {
886                  // non Windows (or cygwin)
887                  absolutePath=(field[0]=='/');
888                } else {
889                  //Windows (non cycgwin)
890                  absolutePath=(field[0]=='\\');
891                  // but allow for :
892                  if (strchr(field.c_str(),':'))
893                    absolutePath=true;
894                }
895                if (absolutePath) {
896                  fileName = field;
897                } else if (field[0]=='~') {
898                  char * environVar = getenv("HOME");
899                  if (environVar) {
900                    std::string home(environVar);
901                    field=field.erase(0,1);
902                    fileName = home+field;
903                  } else {
904                    fileName=field;
905                  }
906                } else {
907                  fileName = directory+field;
908                  // See if gmpl (model & data) - or even lp file
909                  int length = field.size();
910                  int percent = field.find('%');
911                  if (percent<length&&percent>0) {
912                    gmpl=1;
913                    fileName = directory+field.substr(0,percent);
914                    gmplData = directory+field.substr(percent+1);
915                    if (percent<length-1)
916                      gmpl=2; // two files
917                    printf("GMPL model file %s and data file %s\n",
918                           fileName.c_str(),gmplData.c_str());
919                  }
920                }
921                std::string name=fileName;
922                if (fileCoinReadable(name)) {
923                  // can open - lets go for it
924                  canOpen=true;
925                  if (gmpl==2) {
926                    FILE *fp;
927                    fp=fopen(gmplData.c_str(),"r");
928                    if (fp) {
929                      fclose(fp);
930                    } else {
931                      canOpen=false;
932                      std::cout<<"Unable to open file "<<gmplData<<std::endl;
933                    }
934                  }
935                } else {
936                  std::cout<<"Unable to open file "<<fileName<<std::endl;
937                }
938              }
939              if (canOpen) {
940                int status;
941                if (!gmpl)
942                  status =models[iModel].readMps(fileName.c_str(),
943                                                 keepImportNames!=0,
944                                                 allowImportErrors!=0);
945                else if (gmpl>0)
946                  status= models[iModel].readGMPL(fileName.c_str(),
947                                                  (gmpl==2) ? gmplData.c_str() : NULL,
948                                                  keepImportNames!=0);
949                else
950                  status= models[iModel].readLp(fileName.c_str(),1.0e-12);
951                if (!status||(status>0&&allowImportErrors)) {
952                  goodModels[iModel]=true;
953                  // sets to all slack (not necessary?)
954                  models[iModel].createStatus();
955                  time2 = CoinCpuTime();
956                  totalTime += time2-time1;
957                  time1=time2;
958                  // Go to canned file if just input file
959                  if (CbcOrClpRead_mode==2&&argc==2) {
960                    // only if ends .mps
961                    char * find = (char *)strstr(fileName.c_str(),".mps");
962                    if (find&&find[4]=='\0') {
963                      find[1]='p'; find[2]='a';find[3]='r';
964                      FILE *fp=fopen(fileName.c_str(),"r");
965                      if (fp) {
966                        CbcOrClpReadCommand=fp; // Read from that file
967                        CbcOrClpRead_mode=-1;
968                      }
969                    }
970                  }
971                } else {
972                  // errors
973                  std::cout<<"There were "<<status<<
974                    " errors on input"<<std::endl;
975                }
976              }
977            }
978            break;
979          case EXPORT:
980            if (goodModels[iModel]) {
981              double objScale = 
982                parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
983              if (objScale!=1.0) {
984                int iColumn;
985                int numberColumns=models[iModel].numberColumns();
986                double * dualColumnSolution = 
987                  models[iModel].dualColumnSolution();
988                ClpObjective * obj = models[iModel].objectiveAsObject();
989                assert(dynamic_cast<ClpLinearObjective *> (obj));
990                double offset;
991                double * objective = obj->gradient(NULL,NULL,offset,true);
992                for (iColumn=0;iColumn<numberColumns;iColumn++) {
993                  dualColumnSolution[iColumn] *= objScale;
994                  objective[iColumn] *= objScale;;
995                }
996                int iRow;
997                int numberRows=models[iModel].numberRows();
998                double * dualRowSolution = 
999                  models[iModel].dualRowSolution();
1000                for (iRow=0;iRow<numberRows;iRow++) 
1001                  dualRowSolution[iRow] *= objScale;
1002                models[iModel].setObjectiveOffset(objScale*models[iModel].objectiveOffset());
1003              }
1004              // get next field
1005              field = CoinReadGetString(argc,argv);
1006              if (field=="$") {
1007                field = parameters[iParam].stringValue();
1008              } else if (field=="EOL") {
1009                parameters[iParam].printString();
1010                break;
1011              } else {
1012                parameters[iParam].setStringValue(field);
1013              }
1014              std::string fileName;
1015              bool canOpen=false;
1016              if (field[0]=='/'||field[0]=='\\') {
1017                fileName = field;
1018              } else if (field[0]=='~') {
1019                char * environVar = getenv("HOME");
1020                if (environVar) {
1021                  std::string home(environVar);
1022                  field=field.erase(0,1);
1023                  fileName = home+field;
1024                } else {
1025                  fileName=field;
1026                }
1027              } else {
1028                fileName = directory+field;
1029              }
1030              FILE *fp=fopen(fileName.c_str(),"w");
1031              if (fp) {
1032                // can open - lets go for it
1033                fclose(fp);
1034                canOpen=true;
1035              } else {
1036                std::cout<<"Unable to open file "<<fileName<<std::endl;
1037              }
1038              if (canOpen) {
1039                // If presolve on then save presolved
1040                bool deleteModel2=false;
1041                ClpSimplex * model2 = models+iModel;
1042                if (dualize&&dualize<3) {
1043                  model2 = ((ClpSimplexOther *) model2)->dualOfModel();
1044                  printf("Dual of model has %d rows and %d columns\n",
1045                         model2->numberRows(),model2->numberColumns());
1046                  model2->setOptimizationDirection(1.0);
1047                  preSolve=0; // as picks up from model
1048                }
1049                if (preSolve) {
1050                  ClpPresolve pinfo;
1051                  int presolveOptions2 = presolveOptions&~0x40000000;
1052                  if ((presolveOptions2&0xffff)!=0)
1053                    pinfo.setPresolveActions(presolveOptions2);
1054                  pinfo.setSubstitution(substitution);
1055                  if ((printOptions&1)!=0)
1056                    pinfo.statistics();
1057                  double presolveTolerance = 
1058                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1059                  model2 = 
1060                    pinfo.presolvedModel(models[iModel],presolveTolerance,
1061                                         true,preSolve,false,false);
1062                  if (model2) {
1063                    printf("Saving presolved model on %s\n",
1064                           fileName.c_str());
1065                    deleteModel2=true;
1066                  } else {
1067                    printf("Presolved model looks infeasible - saving original on %s\n",
1068                           fileName.c_str());
1069                    deleteModel2=false;
1070                    model2 = models+iModel;
1071
1072                  }
1073                } else {
1074                  printf("Saving model on %s\n",
1075                           fileName.c_str());
1076                }
1077#if 0
1078                // Convert names
1079                int iRow;
1080                int numberRows=model2->numberRows();
1081                int iColumn;
1082                int numberColumns=model2->numberColumns();
1083
1084                char ** rowNames = NULL;
1085                char ** columnNames = NULL;
1086                if (model2->lengthNames()) {
1087                  rowNames = new char * [numberRows];
1088                  for (iRow=0;iRow<numberRows;iRow++) {
1089                    rowNames[iRow] =
1090                      CoinStrdup(model2->rowName(iRow).c_str());
1091#ifdef STRIPBLANKS
1092                    char * xx = rowNames[iRow];
1093                    int i;
1094                    int length = strlen(xx);
1095                    int n=0;
1096                    for (i=0;i<length;i++) {
1097                      if (xx[i]!=' ')
1098                        xx[n++]=xx[i];
1099                    }
1100                    xx[n]='\0';
1101#endif
1102                  }
1103                 
1104                  columnNames = new char * [numberColumns];
1105                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1106                    columnNames[iColumn] =
1107                      CoinStrdup(model2->columnName(iColumn).c_str());
1108#ifdef STRIPBLANKS
1109                    char * xx = columnNames[iColumn];
1110                    int i;
1111                    int length = strlen(xx);
1112                    int n=0;
1113                    for (i=0;i<length;i++) {
1114                      if (xx[i]!=' ')
1115                        xx[n++]=xx[i];
1116                    }
1117                    xx[n]='\0';
1118#endif
1119                  }
1120                }
1121                CoinMpsIO writer;
1122                writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
1123                                  model2->getColLower(), model2->getColUpper(),
1124                                  model2->getObjCoefficients(),
1125                                  (const char*) 0 /*integrality*/,
1126                                  model2->getRowLower(), model2->getRowUpper(),
1127                                  columnNames, rowNames);
1128                // Pass in array saying if each variable integer
1129                writer.copyInIntegerInformation(model2->integerInformation());
1130                writer.setObjectiveOffset(model2->objectiveOffset());
1131                writer.writeMps(fileName.c_str(),0,1,1);
1132                if (rowNames) {
1133                  for (iRow=0;iRow<numberRows;iRow++) {
1134                    free(rowNames[iRow]);
1135                  }
1136                  delete [] rowNames;
1137                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1138                    free(columnNames[iColumn]);
1139                  }
1140                  delete [] columnNames;
1141                }
1142#else
1143                model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
1144#endif
1145                if (deleteModel2)
1146                  delete model2;
1147                time2 = CoinCpuTime();
1148                totalTime += time2-time1;
1149                time1=time2;
1150              }
1151            } else {
1152              std::cout<<"** Current model not valid"<<std::endl;
1153            }
1154            break;
1155          case BASISIN:
1156            if (goodModels[iModel]) {
1157              // get next field
1158              field = CoinReadGetString(argc,argv);
1159              if (field=="$") {
1160                field = parameters[iParam].stringValue();
1161              } else if (field=="EOL") {
1162                parameters[iParam].printString();
1163                break;
1164              } else {
1165                parameters[iParam].setStringValue(field);
1166              }
1167              std::string fileName;
1168              bool canOpen=false;
1169              if (field=="-") {
1170                // stdin
1171                canOpen=true;
1172                fileName = "-";
1173              } else {
1174                if (field[0]=='/'||field[0]=='\\') {
1175                  fileName = field;
1176                } else if (field[0]=='~') {
1177                  char * environVar = getenv("HOME");
1178                  if (environVar) {
1179                    std::string home(environVar);
1180                    field=field.erase(0,1);
1181                    fileName = home+field;
1182                  } else {
1183                    fileName=field;
1184                  }
1185                } else {
1186                  fileName = directory+field;
1187                }
1188                FILE *fp=fopen(fileName.c_str(),"r");
1189                if (fp) {
1190                  // can open - lets go for it
1191                  fclose(fp);
1192                  canOpen=true;
1193                } else {
1194                  std::cout<<"Unable to open file "<<fileName<<std::endl;
1195                }
1196              }
1197              if (canOpen) {
1198                int values = models[iModel].readBasis(fileName.c_str());
1199                if (values==0)
1200                  basisHasValues=-1;
1201                else
1202                  basisHasValues=1;
1203              }
1204            } else {
1205              std::cout<<"** Current model not valid"<<std::endl;
1206            }
1207            break;
1208          case PRINTMASK:
1209            // get next field
1210            {
1211              std::string name = CoinReadGetString(argc,argv);
1212              if (name!="EOL") {
1213                parameters[iParam].setStringValue(name);
1214                printMask = name;
1215              } else {
1216                parameters[iParam].printString();
1217              }
1218            }
1219            break;
1220          case BASISOUT:
1221            if (goodModels[iModel]) {
1222              // get next field
1223              field = CoinReadGetString(argc,argv);
1224              if (field=="$") {
1225                field = parameters[iParam].stringValue();
1226              } else if (field=="EOL") {
1227                parameters[iParam].printString();
1228                break;
1229              } else {
1230                parameters[iParam].setStringValue(field);
1231              }
1232              std::string fileName;
1233              bool canOpen=false;
1234              if (field[0]=='/'||field[0]=='\\') {
1235                fileName = field;
1236              } else if (field[0]=='~') {
1237                char * environVar = getenv("HOME");
1238                if (environVar) {
1239                  std::string home(environVar);
1240                  field=field.erase(0,1);
1241                  fileName = home+field;
1242                } else {
1243                  fileName=field;
1244                }
1245              } else {
1246                fileName = directory+field;
1247              }
1248              FILE *fp=fopen(fileName.c_str(),"w");
1249              if (fp) {
1250                // can open - lets go for it
1251                fclose(fp);
1252                canOpen=true;
1253              } else {
1254                std::cout<<"Unable to open file "<<fileName<<std::endl;
1255              }
1256              if (canOpen) {
1257                ClpSimplex * model2 = models+iModel;
1258                model2->writeBasis(fileName.c_str(),outputFormat>1,outputFormat-2);
1259                time2 = CoinCpuTime();
1260                totalTime += time2-time1;
1261                time1=time2;
1262              }
1263            } else {
1264              std::cout<<"** Current model not valid"<<std::endl;
1265            }
1266            break;
1267          case SAVE:
1268            {
1269              // get next field
1270              field = CoinReadGetString(argc,argv);
1271              if (field=="$") {
1272                field = parameters[iParam].stringValue();
1273              } else if (field=="EOL") {
1274                parameters[iParam].printString();
1275                break;
1276              } else {
1277                parameters[iParam].setStringValue(field);
1278              }
1279              std::string fileName;
1280              bool canOpen=false;
1281              if (field[0]=='/'||field[0]=='\\') {
1282                fileName = field;
1283              } else if (field[0]=='~') {
1284                char * environVar = getenv("HOME");
1285                if (environVar) {
1286                  std::string home(environVar);
1287                  field=field.erase(0,1);
1288                  fileName = home+field;
1289                } else {
1290                  fileName=field;
1291                }
1292              } else {
1293                fileName = directory+field;
1294              }
1295              FILE *fp=fopen(fileName.c_str(),"wb");
1296              if (fp) {
1297                // can open - lets go for it
1298                fclose(fp);
1299                canOpen=true;
1300              } else {
1301                std::cout<<"Unable to open file "<<fileName<<std::endl;
1302              }
1303              if (canOpen) {
1304                int status;
1305                // If presolve on then save presolved
1306                bool deleteModel2=false;
1307                ClpSimplex * model2 = models+iModel;
1308                if (preSolve) {
1309                  ClpPresolve pinfo;
1310                  double presolveTolerance = 
1311                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1312                  model2 = 
1313                    pinfo.presolvedModel(models[iModel],presolveTolerance,
1314                                         false,preSolve);
1315                  if (model2) {
1316                    printf("Saving presolved model on %s\n",
1317                           fileName.c_str());
1318                    deleteModel2=true;
1319                  } else {
1320                    printf("Presolved model looks infeasible - saving original on %s\n",
1321                           fileName.c_str());
1322                    deleteModel2=false;
1323                    model2 = models+iModel;
1324
1325                  }
1326                } else {
1327                  printf("Saving model on %s\n",
1328                           fileName.c_str());
1329                }
1330                status =model2->saveModel(fileName.c_str());
1331                if (deleteModel2)
1332                  delete model2;
1333                if (!status) {
1334                  goodModels[iModel]=true;
1335                  time2 = CoinCpuTime();
1336                  totalTime += time2-time1;
1337                  time1=time2;
1338                } else {
1339                  // errors
1340                  std::cout<<"There were errors on output"<<std::endl;
1341                }
1342              }
1343            }
1344            break;
1345          case RESTORE:
1346            {
1347              // get next field
1348              field = CoinReadGetString(argc,argv);
1349              if (field=="$") {
1350                field = parameters[iParam].stringValue();
1351              } else if (field=="EOL") {
1352                parameters[iParam].printString();
1353                break;
1354              } else {
1355                parameters[iParam].setStringValue(field);
1356              }
1357              std::string fileName;
1358              bool canOpen=false;
1359              if (field[0]=='/'||field[0]=='\\') {
1360                fileName = field;
1361              } else if (field[0]=='~') {
1362                char * environVar = getenv("HOME");
1363                if (environVar) {
1364                  std::string home(environVar);
1365                  field=field.erase(0,1);
1366                  fileName = home+field;
1367                } else {
1368                  fileName=field;
1369                }
1370              } else {
1371                fileName = directory+field;
1372              }
1373              FILE *fp=fopen(fileName.c_str(),"rb");
1374              if (fp) {
1375                // can open - lets go for it
1376                fclose(fp);
1377                canOpen=true;
1378              } else {
1379                std::cout<<"Unable to open file "<<fileName<<std::endl;
1380              }
1381              if (canOpen) {
1382                int status =models[iModel].restoreModel(fileName.c_str());
1383                if (!status) {
1384                  goodModels[iModel]=true;
1385                  time2 = CoinCpuTime();
1386                  totalTime += time2-time1;
1387                  time1=time2;
1388                } else {
1389                  // errors
1390                  std::cout<<"There were errors on input"<<std::endl;
1391                }
1392              }
1393            }
1394            break;
1395          case MAXIMIZE:
1396            models[iModel].setOptimizationDirection(-1);
1397            break;
1398          case MINIMIZE:
1399            models[iModel].setOptimizationDirection(1);
1400            break;
1401          case ALLSLACK:
1402            models[iModel].allSlackBasis(true);
1403            break;
1404          case REVERSE:
1405            if (goodModels[iModel]) {
1406              int iColumn;
1407              int numberColumns=models[iModel].numberColumns();
1408              double * dualColumnSolution = 
1409                models[iModel].dualColumnSolution();
1410              ClpObjective * obj = models[iModel].objectiveAsObject();
1411              assert(dynamic_cast<ClpLinearObjective *> (obj));
1412              double offset;
1413              double * objective = obj->gradient(NULL,NULL,offset,true);
1414              for (iColumn=0;iColumn<numberColumns;iColumn++) {
1415                dualColumnSolution[iColumn] = -dualColumnSolution[iColumn];
1416                objective[iColumn] = -objective[iColumn];
1417              }
1418              int iRow;
1419              int numberRows=models[iModel].numberRows();
1420              double * dualRowSolution = 
1421                models[iModel].dualRowSolution();
1422              for (iRow=0;iRow<numberRows;iRow++) {
1423                dualRowSolution[iRow] = -dualRowSolution[iRow];
1424              }
1425              models[iModel].setObjectiveOffset(-models[iModel].objectiveOffset());
1426            }
1427            break;
1428          case DIRECTORY:
1429            {
1430              std::string name = CoinReadGetString(argc,argv);
1431              if (name!="EOL") {
1432                int length=name.length();
1433                if (name[length-1]==dirsep) {
1434                  directory = name;
1435                } else {
1436                  directory = name+dirsep;
1437                }
1438                parameters[iParam].setStringValue(directory);
1439              } else {
1440                parameters[iParam].printString();
1441              }
1442            }
1443            break;
1444          case DIRSAMPLE:
1445            {
1446              std::string name = CoinReadGetString(argc,argv);
1447              if (name!="EOL") {
1448                int length=name.length();
1449                if (name[length-1]==dirsep) {
1450                  dirSample = name;
1451                } else {
1452                  dirSample = name+dirsep;
1453                }
1454                parameters[iParam].setStringValue(dirSample);
1455              } else {
1456                parameters[iParam].printString();
1457              }
1458            }
1459            break;
1460          case DIRNETLIB:
1461            {
1462              std::string name = CoinReadGetString(argc,argv);
1463              if (name!="EOL") {
1464                int length=name.length();
1465                if (name[length-1]==dirsep) {
1466                  dirNetlib = name;
1467                } else {
1468                  dirNetlib = name+dirsep;
1469                }
1470                parameters[iParam].setStringValue(dirNetlib);
1471              } else {
1472                parameters[iParam].printString();
1473              }
1474            }
1475            break;
1476          case DIRMIPLIB:
1477            {
1478              std::string name = CoinReadGetString(argc,argv);
1479              if (name!="EOL") {
1480                int length=name.length();
1481                if (name[length-1]==dirsep) {
1482                  dirMiplib = name;
1483                } else {
1484                  dirMiplib = name+dirsep;
1485                }
1486                parameters[iParam].setStringValue(dirMiplib);
1487              } else {
1488                parameters[iParam].printString();
1489              }
1490            }
1491            break;
1492          case STDIN:
1493            CbcOrClpRead_mode=-1;
1494            break;
1495          case NETLIB_DUAL:
1496          case NETLIB_EITHER:
1497          case NETLIB_BARRIER:
1498          case NETLIB_PRIMAL:
1499          case NETLIB_TUNE:
1500            {
1501              // create fields for unitTest
1502              const char * fields[4];
1503              int nFields=4;
1504              fields[0]="fake main from unitTest";
1505              std::string mpsfield = "-dirSample=";
1506              mpsfield += dirSample.c_str();
1507              fields[1]=mpsfield.c_str();
1508              std::string netfield = "-dirNetlib=";
1509              netfield += dirNetlib.c_str();
1510              fields[2]=netfield.c_str();
1511              fields[3]="-netlib";
1512              int algorithm;
1513              if (type==NETLIB_DUAL) {
1514                std::cerr<<"Doing netlib with dual algorithm"<<std::endl;
1515                algorithm =0;
1516              } else if (type==NETLIB_BARRIER) {
1517                std::cerr<<"Doing netlib with barrier algorithm"<<std::endl;
1518                algorithm =2;
1519              } else if (type==NETLIB_EITHER) {
1520                std::cerr<<"Doing netlib with dual or primal algorithm"<<std::endl;
1521                algorithm =3;
1522              } else if (type==NETLIB_TUNE) {
1523                std::cerr<<"Doing netlib with best algorithm!"<<std::endl;
1524                algorithm =5;
1525                // uncomment next to get active tuning
1526                // algorithm=6;
1527              } else {
1528                std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
1529                algorithm=1;
1530              }
1531              int specialOptions = models[iModel].specialOptions();
1532              models[iModel].setSpecialOptions(0);
1533              mainTest(nFields,fields,algorithm,models[iModel],
1534                       (preSolve!=0),specialOptions,doVector!=0);
1535            }
1536            break;
1537          case UNITTEST:
1538            {
1539              // create fields for unitTest
1540              const char * fields[2];
1541              int nFields=2;
1542              fields[0]="fake main from unitTest";
1543              std::string dirfield = "-dirSample=";
1544              dirfield += dirSample.c_str();
1545              fields[1]=dirfield.c_str();
1546              int specialOptions = models[iModel].specialOptions();
1547              models[iModel].setSpecialOptions(0);
1548              int algorithm=-1;
1549              if (models[iModel].numberRows())
1550                algorithm=7;
1551              mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions,doVector!=0);
1552            }
1553            break;
1554          case FAKEBOUND:
1555            if (goodModels[iModel]) {
1556              // get bound
1557              double value = CoinReadGetDoubleField(argc,argv,&valid);
1558              if (!valid) {
1559                std::cout<<"Setting "<<parameters[iParam].name()<<
1560                  " to DEBUG "<<value<<std::endl;
1561                int iRow;
1562                int numberRows=models[iModel].numberRows();
1563                double * rowLower = models[iModel].rowLower();
1564                double * rowUpper = models[iModel].rowUpper();
1565                for (iRow=0;iRow<numberRows;iRow++) {
1566                  // leave free ones for now
1567                  if (rowLower[iRow]>-1.0e20||rowUpper[iRow]<1.0e20) {
1568                    rowLower[iRow]=CoinMax(rowLower[iRow],-value);
1569                    rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
1570                  }
1571                }
1572                int iColumn;
1573                int numberColumns=models[iModel].numberColumns();
1574                double * columnLower = models[iModel].columnLower();
1575                double * columnUpper = models[iModel].columnUpper();
1576                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1577                  // leave free ones for now
1578                  if (columnLower[iColumn]>-1.0e20||
1579                      columnUpper[iColumn]<1.0e20) {
1580                    columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
1581                    columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
1582                  }
1583                }
1584              } else if (valid==1) {
1585                abort();
1586              } else {
1587                std::cout<<"enter value for "<<parameters[iParam].name()<<
1588                  std::endl;
1589              }
1590            }
1591            break;
1592          case REALLY_SCALE:
1593            if (goodModels[iModel]) {
1594              ClpSimplex newModel(models[iModel],
1595                                  models[iModel].scalingFlag());
1596              printf("model really really scaled\n");
1597              models[iModel]=newModel;
1598            }
1599            break;
1600          case USERCLP:
1601            // Replace the sample code by whatever you want
1602            if (goodModels[iModel]) {
1603              ClpSimplex * thisModel = &models[iModel];
1604              printf("Dummy user code - model has %d rows and %d columns\n",
1605                     thisModel->numberRows(),thisModel->numberColumns());
1606            }
1607            break;
1608          case HELP:
1609            std::cout<<"Coin LP version "<<CLPVERSION
1610                     <<", build "<<__DATE__<<std::endl;
1611            std::cout<<"Non default values:-"<<std::endl;
1612            std::cout<<"Perturbation "<<models[0].perturbation()<<" (default 100)"
1613                     <<std::endl;
1614            CoinReadPrintit(
1615                    "Presolve being done with 5 passes\n\
1616Dual steepest edge steep/partial on matrix shape and factorization density\n\
1617Clpnnnn taken out of messages\n\
1618If Factorization frequency default then done on size of matrix\n\n\
1619(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
1620You can switch to interactive mode at any time so\n\
1621clp watson.mps -scaling off -primalsimplex\nis the same as\n\
1622clp watson.mps -\nscaling off\nprimalsimplex"
1623                    );
1624            break;
1625          case SOLUTION:
1626            if (goodModels[iModel]) {
1627              // get next field
1628              field = CoinReadGetString(argc,argv);
1629              if (field=="$") {
1630                field = parameters[iParam].stringValue();
1631              } else if (field=="EOL") {
1632                parameters[iParam].printString();
1633                break;
1634              } else {
1635                parameters[iParam].setStringValue(field);
1636              }
1637              std::string fileName;
1638              FILE *fp=NULL;
1639              if (field=="-"||field=="EOL"||field=="stdout") {
1640                // stdout
1641                fp=stdout;
1642                fprintf(fp,"\n");
1643              } else if (field=="stderr") {
1644                // stderr
1645                fp=stderr;
1646                fprintf(fp,"\n");
1647              } else {
1648                if (field[0]=='/'||field[0]=='\\') {
1649                  fileName = field;
1650                } else if (field[0]=='~') {
1651                  char * environVar = getenv("HOME");
1652                  if (environVar) {
1653                    std::string home(environVar);
1654                    field=field.erase(0,1);
1655                    fileName = home+field;
1656                  } else {
1657                    fileName=field;
1658                  }
1659                } else {
1660                  fileName = directory+field;
1661                }
1662                fp=fopen(fileName.c_str(),"w");
1663              }
1664              if (fp) {
1665                // Write solution header (suggested by Luigi Poderico)
1666                double objValue = models[iModel].getObjValue()*models[iModel].getObjSense();
1667                int iStat = models[iModel].status();
1668                if (iStat==0) {
1669                  fprintf(fp, "optimal\n" );
1670                } else if (iStat==1) {
1671                  // infeasible
1672                  fprintf(fp, "infeasible\n" );
1673                } else if (iStat==2) {
1674                  // unbounded
1675                  fprintf(fp, "unbounded\n" );
1676                } else if (iStat==3) {
1677                  fprintf(fp, "stopped on iterations or time\n" );
1678                } else if (iStat==4) {
1679                  fprintf(fp, "stopped on difficulties\n" );
1680                } else if (iStat==5) {
1681                  fprintf(fp, "stopped on ctrl-c\n" );
1682                } else {
1683                  fprintf(fp, "status unknown\n" );
1684                }
1685                fprintf(fp, "Objective value %15.8g\n", objValue);
1686                // make fancy later on
1687                int iRow;
1688                int numberRows=models[iModel].numberRows();
1689                int lengthName = models[iModel].lengthNames(); // 0 if no names
1690                // in general I don't want to pass around massive
1691                // amounts of data but seems simpler here
1692                std::vector<std::string> rowNames =
1693                  *(models[iModel].rowNames());
1694                std::vector<std::string> columnNames =
1695                  *(models[iModel].columnNames());
1696
1697                double * dualRowSolution = models[iModel].dualRowSolution();
1698                double * primalRowSolution = 
1699                  models[iModel].primalRowSolution();
1700                double * rowLower = models[iModel].rowLower();
1701                double * rowUpper = models[iModel].rowUpper();
1702                double primalTolerance = models[iModel].primalTolerance();
1703                char format[6];
1704                sprintf(format,"%%-%ds",CoinMax(lengthName,8));
1705                bool doMask = (printMask!=""&&lengthName);
1706                int * maskStarts=NULL;
1707                int maxMasks=0;
1708                char ** masks =NULL;
1709                if (doMask) {
1710                  int nAst =0;
1711                  const char * pMask2 = printMask.c_str();
1712                  char pMask[100];
1713                  int iChar;
1714                  int lengthMask = strlen(pMask2);
1715                  assert (lengthMask<100);
1716                  if (*pMask2=='"') {
1717                    if (pMask2[lengthMask-1]!='"') {
1718                      printf("mismatched \" in mask %s\n",pMask2);
1719                      break;
1720                    } else {
1721                      strcpy(pMask,pMask2+1);
1722                      *strchr(pMask,'"')='\0';
1723                    }
1724                  } else if (*pMask2=='\'') {
1725                    if (pMask2[lengthMask-1]!='\'') {
1726                      printf("mismatched ' in mask %s\n",pMask2);
1727                      break;
1728                    } else {
1729                      strcpy(pMask,pMask2+1);
1730                      *strchr(pMask,'\'')='\0';
1731                    }
1732                  } else {
1733                    strcpy(pMask,pMask2);
1734                  }
1735                  if (lengthMask>lengthName) {
1736                    printf("mask %s too long - skipping\n",pMask);
1737                    break;
1738                  }
1739                  maxMasks = 1;
1740                  for (iChar=0;iChar<lengthMask;iChar++) {
1741                    if (pMask[iChar]=='*') {
1742                      nAst++;
1743                      maxMasks *= (lengthName+1);
1744                    }
1745                  }
1746                  int nEntries = 1;
1747                  maskStarts = new int[lengthName+2];
1748                  masks = new char * [maxMasks];
1749                  char ** newMasks = new char * [maxMasks];
1750                  int i;
1751                  for (i=0;i<maxMasks;i++) {
1752                    masks[i] = new char[lengthName+1];
1753                    newMasks[i] = new char[lengthName+1];
1754                  }
1755                  strcpy(masks[0],pMask);
1756                  for (int iAst=0;iAst<nAst;iAst++) {
1757                    int nOldEntries = nEntries;
1758                    nEntries=0;
1759                    for (int iEntry = 0;iEntry<nOldEntries;iEntry++) {
1760                      char * oldMask = masks[iEntry];
1761                      char * ast = strchr(oldMask,'*');
1762                      assert (ast);
1763                      int length = strlen(oldMask)-1;
1764                      int nBefore = ast-oldMask;
1765                      int nAfter = length-nBefore;
1766                      // and add null
1767                      nAfter++;
1768                      for (int i=0;i<=lengthName-length;i++) {
1769                        char * maskOut = newMasks[nEntries];
1770   CoinMemcpyN(oldMask,nBefore,maskOut);
1771                        for (int k=0;k<i;k++) 
1772                          maskOut[k+nBefore]='?';
1773   CoinMemcpyN(ast+1,nAfter,maskOut+nBefore+i);
1774                        nEntries++;
1775                        assert (nEntries<=maxMasks);
1776                      }
1777                    }
1778                    char ** temp = masks;
1779                    masks = newMasks;
1780                    newMasks = temp;
1781                  }
1782                  // Now extend and sort
1783                  int * sort = new int[nEntries];
1784                  for (i=0;i<nEntries;i++) {
1785                    char * maskThis = masks[i];
1786                    int length = strlen(maskThis);
1787                    while (maskThis[length-1]==' ')
1788                      length--;
1789                    maskThis[length]='\0';
1790                    sort[i]=length;
1791                  }
1792                  CoinSort_2(sort,sort+nEntries,masks);
1793                  int lastLength=-1;
1794                  for (i=0;i<nEntries;i++) {
1795                    int length = sort[i];
1796                    while (length>lastLength) 
1797                      maskStarts[++lastLength] = i;
1798                  }
1799                  maskStarts[++lastLength]=nEntries;
1800                  delete [] sort;
1801                  for (i=0;i<maxMasks;i++)
1802                    delete [] newMasks[i];
1803                  delete [] newMasks;
1804                }
1805                if (printMode>2) {
1806                  for (iRow=0;iRow<numberRows;iRow++) {
1807                    int type=printMode-3;
1808                    if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
1809                        primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
1810                      fprintf(fp,"** ");
1811                      type=2;
1812                    } else if (fabs(primalRowSolution[iRow])>1.0e-8) {
1813                      type=1;
1814                    } else if (numberRows<50) {
1815                      type=3;
1816                    } 
1817                    if (doMask&&!maskMatches(maskStarts,masks,rowNames[iRow]))
1818                      type=0;
1819                    if (type) {
1820                      fprintf(fp,"%7d ",iRow);
1821                      if (lengthName)
1822                        fprintf(fp,format,rowNames[iRow].c_str());
1823                      fprintf(fp,"%15.8g        %15.8g\n",primalRowSolution[iRow],
1824                              dualRowSolution[iRow]);
1825                    }
1826                  }
1827                }
1828                int iColumn;
1829                int numberColumns=models[iModel].numberColumns();
1830                double * dualColumnSolution = 
1831  models[iModel].dualColumnSolution();
1832                double * primalColumnSolution = 
1833  models[iModel].primalColumnSolution();
1834                double * columnLower = models[iModel].columnLower();
1835                double * columnUpper = models[iModel].columnUpper();
1836                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1837                  int type=(printMode>3) ? 1 : 0;
1838                  if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
1839                      primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
1840                    fprintf(fp,"** ");
1841                    type=2;
1842                  } else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
1843                    type=1;
1844                  } else if (numberColumns<50) {
1845                    type=3;
1846                  }
1847                  if (doMask&&!maskMatches(maskStarts,masks,
1848                                           columnNames[iColumn]))
1849                    type =0;
1850                  if (type) {
1851                    fprintf(fp,"%7d ",iColumn);
1852                    if (lengthName)
1853                      fprintf(fp,format,columnNames[iColumn].c_str());
1854                    fprintf(fp,"%15.8g        %15.8g\n",
1855                            primalColumnSolution[iColumn],
1856                            dualColumnSolution[iColumn]);
1857                  }
1858                }
1859                if (fp!=stdout)
1860                  fclose(fp);
1861                if (masks) {
1862                  delete [] maskStarts;
1863                  for (int i=0;i<maxMasks;i++)
1864                    delete [] masks[i];
1865                  delete [] masks;
1866                }
1867              } else {
1868                std::cout<<"Unable to open file "<<fileName<<std::endl;
1869              }
1870            } else {
1871              std::cout<<"** Current model not valid"<<std::endl;
1872             
1873            }
1874         
1875            break;
1876          case SAVESOL:
1877            if (goodModels[iModel]) {
1878              // get next field
1879              field = CoinReadGetString(argc,argv);
1880              if (field=="$") {
1881                field = parameters[iParam].stringValue();
1882              } else if (field=="EOL") {
1883                parameters[iParam].printString();
1884                break;
1885              } else {
1886                parameters[iParam].setStringValue(field);
1887              }
1888              std::string fileName;
1889              if (field[0]=='/'||field[0]=='\\') {
1890                fileName = field;
1891              } else if (field[0]=='~') {
1892                char * environVar = getenv("HOME");
1893                if (environVar) {
1894                  std::string home(environVar);
1895                  field=field.erase(0,1);
1896                  fileName = home+field;
1897                } else {
1898                  fileName=field;
1899                }
1900              } else {
1901                fileName = directory+field;
1902              }
1903              saveSolution(models+iModel,fileName);
1904            } else {
1905              std::cout<<"** Current model not valid"<<std::endl;
1906             
1907            }
1908            break;
1909          default:
1910            abort();
1911          }
1912        } 
1913      } else if (!numberMatches) {
1914        std::cout<<"No match for "<<field<<" - ? for list of commands"
1915                 <<std::endl;
1916      } else if (numberMatches==1) {
1917        if (!numberQuery) {
1918          std::cout<<"Short match for "<<field<<" - completion: ";
1919          std::cout<<parameters[firstMatch].matchName()<<std::endl;
1920        } else if (numberQuery) {
1921          std::cout<<parameters[firstMatch].matchName()<<" : ";
1922          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
1923          if (numberQuery>=2) 
1924            parameters[firstMatch].printLongHelp();
1925        }
1926      } else {
1927        if (!numberQuery) 
1928          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
1929                   <<std::endl;
1930        else
1931          std::cout<<"Completions of "<<field<<":"<<std::endl;
1932        for ( iParam=0; iParam<numberParameters; iParam++ ) {
1933          int match = parameters[iParam].matches(field);
1934          if (match&&parameters[iParam].displayThis()) {
1935            std::cout<<parameters[iParam].matchName();
1936            if (numberQuery>=2) 
1937              std::cout<<" : "<<parameters[iParam].shortHelp();
1938            std::cout<<std::endl;
1939          }
1940        }
1941      }
1942    }
1943    delete [] models;
1944    delete [] goodModels;
1945  }
1946  // By now all memory should be freed
1947#ifdef DMALLOC
1948  dmalloc_log_unfreed();
1949  dmalloc_shutdown();
1950#endif
1951  return 0;
1952}   
1953static void breakdown(const char * name, int numberLook, const double * region)
1954{
1955  double range[] = {
1956    -COIN_DBL_MAX,
1957    -1.0e15,-1.0e11,-1.0e8,-1.0e5,-1.0e4,-1.0e3,-1.0e2,-1.0e1,
1958    -1.0,
1959    -1.0e-1,-1.0e-2,-1.0e-3,-1.0e-4,-1.0e-5,-1.0e-8,-1.0e-11,-1.0e-15,
1960    0.0,
1961    1.0e-15,1.0e-11,1.0e-8,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,
1962    1.0,
1963    1.0e1,1.0e2,1.0e3,1.0e4,1.0e5,1.0e8,1.0e11,1.0e15,
1964    COIN_DBL_MAX};
1965  int nRanges = (int) (sizeof(range)/sizeof(double));
1966  int * number = new int[nRanges];
1967  memset(number,0,nRanges*sizeof(int));
1968  int * numberExact = new int[nRanges];
1969  memset(numberExact,0,nRanges*sizeof(int));
1970  int i;
1971  for ( i=0;i<numberLook;i++) {
1972    double value = region[i];
1973    for (int j=0;j<nRanges;j++) {
1974      if (value==range[j]) {
1975        numberExact[j]++;
1976        break;
1977      } else if (value<range[j]) {
1978        number[j]++;
1979        break;
1980      }
1981    }
1982  }
1983  printf("\n%s has %d entries\n",name,numberLook);
1984  for (i=0;i<nRanges;i++) {
1985    if (number[i]) 
1986      printf("%d between %g and %g",number[i],range[i-1],range[i]);
1987    if (numberExact[i]) {
1988      if (number[i])
1989        printf(", ");
1990      printf("%d exactly at %g",numberExact[i],range[i]);
1991    }
1992    if (number[i]+numberExact[i])
1993      printf("\n");
1994  }
1995  delete [] number;
1996  delete [] numberExact;
1997}
1998void sortOnOther(int * column,
1999                  const CoinBigIndex * rowStart,
2000                  int * order,
2001                  int * other,
2002                  int nRow,
2003                  int nInRow,
2004                  int where)
2005{
2006  if (nRow<2||where>=nInRow)
2007    return;
2008  // do initial sort
2009  int kRow;
2010  int iRow;
2011  for ( kRow=0;kRow<nRow;kRow++) {
2012    iRow = order[kRow];
2013    other[kRow]=column[rowStart[iRow]+where];
2014  }
2015  CoinSort_2(other,other+nRow,order);
2016  int first=0;
2017  iRow=order[0];
2018  int firstC=column[rowStart[iRow]+where];
2019  kRow=1;
2020  while (kRow<nRow) {
2021    int lastC=9999999;;
2022    for (;kRow<nRow+1;kRow++) {
2023      if (kRow<nRow) {
2024        iRow=order[kRow];
2025        lastC=column[rowStart[iRow]+where];
2026      } else {
2027        lastC=9999999;
2028      }
2029      if (lastC>firstC) 
2030        break;
2031    }
2032    // sort
2033    sortOnOther(column,rowStart,order+first,other,kRow-first,
2034                 nInRow,where+1);
2035    firstC=lastC;
2036    first=kRow;
2037  }
2038}
2039static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
2040{
2041  int numberColumns = originalModel->numberColumns();
2042  const char * integerInformation  = originalModel->integerInformation(); 
2043  const double * columnLower = originalModel->columnLower();
2044  const double * columnUpper = originalModel->columnUpper();
2045  int numberIntegers=0;
2046  int numberBinary=0;
2047  int iRow,iColumn;
2048  if (integerInformation) {
2049    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2050      if (integerInformation[iColumn]) {
2051        if (columnUpper[iColumn]>columnLower[iColumn]) {
2052          numberIntegers++;
2053          if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==1) 
2054            numberBinary++;
2055        }
2056      }
2057    }
2058  }
2059  numberColumns = model->numberColumns();
2060  int numberRows = model->numberRows();
2061  columnLower = model->columnLower();
2062  columnUpper = model->columnUpper();
2063  const double * rowLower = model->rowLower();
2064  const double * rowUpper = model->rowUpper();
2065  const double * objective = model->objective();
2066  CoinPackedMatrix * matrix = model->matrix();
2067  CoinBigIndex numberElements = matrix->getNumElements();
2068  const int * columnLength = matrix->getVectorLengths();
2069  //const CoinBigIndex * columnStart = matrix->getVectorStarts();
2070  const double * elementByColumn = matrix->getElements();
2071  int * number = new int[numberRows+1];
2072  memset(number,0,(numberRows+1)*sizeof(int));
2073  int numberObjSingletons=0;
2074  /* cType
2075     0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
2076     8 0/1
2077  */ 
2078  int cType[9];
2079  std::string cName[]={"0.0->inf,","0.0->up,","lo->inf,","lo->up,","free,","fixed,","-inf->0.0,",
2080                       "-inf->up,","0.0->1.0"};
2081  int nObjective=0;
2082  memset(cType,0,sizeof(cType));
2083  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2084    int length=columnLength[iColumn];
2085    if (length==1&&objective[iColumn])
2086      numberObjSingletons++;
2087    number[length]++;
2088    if (objective[iColumn])
2089      nObjective++;
2090    if (columnLower[iColumn]>-1.0e20) {
2091      if (columnLower[iColumn]==0.0) {
2092        if (columnUpper[iColumn]>1.0e20)
2093          cType[0]++;
2094        else if (columnUpper[iColumn]==1.0)
2095          cType[8]++;
2096        else if (columnUpper[iColumn]==0.0)
2097          cType[5]++;
2098        else
2099          cType[1]++;
2100      } else {
2101        if (columnUpper[iColumn]>1.0e20) 
2102          cType[2]++;
2103        else if (columnUpper[iColumn]==columnLower[iColumn])
2104          cType[5]++;
2105        else
2106          cType[3]++;
2107      }
2108    } else {
2109      if (columnUpper[iColumn]>1.0e20) 
2110        cType[4]++;
2111      else if (columnUpper[iColumn]==0.0) 
2112        cType[6]++;
2113      else
2114        cType[7]++;
2115    }
2116  }
2117  /* rType
2118     0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
2119     7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
2120  */ 
2121  int rType[13];
2122  std::string rName[]={"E 0.0,","E 1.0,","E -1.0,","E other,","G 0.0,","G 1.0,","G other,",
2123                       "L 0.0,","L 1.0,","L other,","Range 0.0->1.0,","Range other,","Free"};
2124  memset(rType,0,sizeof(rType));
2125  for (iRow=0;iRow<numberRows;iRow++) {
2126    if (rowLower[iRow]>-1.0e20) {
2127      if (rowLower[iRow]==0.0) {
2128        if (rowUpper[iRow]>1.0e20)
2129          rType[4]++;
2130        else if (rowUpper[iRow]==1.0)
2131          rType[10]++;
2132        else if (rowUpper[iRow]==0.0)
2133          rType[0]++;
2134        else
2135          rType[11]++;
2136      } else if (rowLower[iRow]==1.0) {
2137        if (rowUpper[iRow]>1.0e20) 
2138          rType[5]++;
2139        else if (rowUpper[iRow]==rowLower[iRow])
2140          rType[1]++;
2141        else
2142          rType[11]++;
2143      } else if (rowLower[iRow]==-1.0) {
2144        if (rowUpper[iRow]>1.0e20) 
2145          rType[6]++;
2146        else if (rowUpper[iRow]==rowLower[iRow])
2147          rType[2]++;
2148        else
2149          rType[11]++;
2150      } else {
2151        if (rowUpper[iRow]>1.0e20) 
2152          rType[6]++;
2153        else if (rowUpper[iRow]==rowLower[iRow])
2154          rType[3]++;
2155        else
2156          rType[11]++;
2157      }
2158    } else {
2159      if (rowUpper[iRow]>1.0e20) 
2160        rType[12]++;
2161      else if (rowUpper[iRow]==0.0) 
2162        rType[7]++;
2163      else if (rowUpper[iRow]==1.0) 
2164        rType[8]++;
2165      else
2166        rType[9]++;
2167    }
2168  }
2169  // Basic statistics
2170  printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
2171         numberRows,numberColumns,nObjective,numberElements);
2172  if (number[0]+number[1]) {
2173    printf("There are ");
2174    if (numberObjSingletons)
2175      printf("%d singletons with objective ",numberObjSingletons);
2176    int numberNoObj = number[1]-numberObjSingletons;
2177    if (numberNoObj)
2178      printf("%d singletons with no objective ",numberNoObj);
2179    if (number[0])
2180      printf("** %d columns have no entries",number[0]);
2181    printf("\n");
2182  }
2183  printf("Column breakdown:\n");
2184  int k;
2185  for (k=0;k<(int) (sizeof(cType)/sizeof(int));k++) {
2186    printf("%d of type %s ",cType[k],cName[k].c_str());
2187    if (((k+1)%3)==0)
2188      printf("\n");
2189  }
2190  if ((k%3)!=0)
2191    printf("\n");
2192  printf("Row breakdown:\n");
2193  for (k=0;k<(int) (sizeof(rType)/sizeof(int));k++) {
2194    printf("%d of type %s ",rType[k],rName[k].c_str());
2195    if (((k+1)%3)==0)
2196      printf("\n");
2197  }
2198  if ((k%3)!=0)
2199    printf("\n");
2200#define SYM
2201#ifndef SYM
2202  if (model->logLevel()<2)
2203    return ;
2204#endif
2205  int kMax = model->logLevel()>3 ? 1000000 : 10;
2206  k=0;
2207  for (iRow=1;iRow<=numberRows;iRow++) {
2208    if (number[iRow]) {
2209      k++;
2210      printf("%d columns have %d entries\n",number[iRow],iRow);
2211      if (k==kMax)
2212        break;
2213    }
2214  }
2215  if (k<numberRows) {
2216    int kk=k;
2217    k=0;
2218    for (iRow=numberRows;iRow>=1;iRow--) {
2219      if (number[iRow]) {
2220        k++;
2221        if (k==kMax)
2222          break;
2223      }
2224    }
2225    if (k>kk) {
2226      printf("\n    .........\n\n");
2227      iRow=k;
2228      k=0;
2229      for (;iRow<numberRows;iRow++) {
2230        if (number[iRow]) {
2231          k++;
2232          printf("%d columns have %d entries\n",number[iRow],iRow);
2233          if (k==kMax)
2234            break;
2235        }
2236      }
2237    }
2238  }
2239  delete [] number;
2240  printf("\n\n");
2241  if (model->logLevel()==63
2242#ifdef SYM
2243      ||true
2244#endif
2245      ) {
2246    // get column copy
2247    CoinPackedMatrix columnCopy = *matrix;
2248    const int * columnLength = columnCopy.getVectorLengths();
2249    number = new int[numberRows+1];
2250    memset(number,0,(numberRows+1)*sizeof(int));
2251    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2252      int length=columnLength[iColumn];
2253      number[length]++;
2254    }
2255    k=0;
2256    for (iRow=1;iRow<=numberRows;iRow++) {
2257      if (number[iRow]) {
2258        k++;
2259      }
2260    }
2261    int * row = columnCopy.getMutableIndices();
2262    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2263    double * element = columnCopy.getMutableElements();
2264    int * order = new int[numberColumns];
2265    int * other = new int[numberColumns];
2266    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2267      int length=columnLength[iColumn];
2268      order[iColumn]=iColumn;
2269      other[iColumn]=length;
2270      CoinBigIndex start = columnStart[iColumn];
2271      CoinSort_2(row+start,row+start+length,element+start);
2272    }
2273    CoinSort_2(other,other+numberColumns,order);
2274    int jColumn=number[0]+number[1];
2275    for (iRow=2;iRow<=numberRows;iRow++) {
2276      if (number[iRow]) {
2277        printf("XX %d columns have %d entries\n",number[iRow],iRow);
2278        int kColumn = jColumn+number[iRow];
2279        sortOnOther(row,columnStart,
2280                     order+jColumn,other,number[iRow],iRow,0);
2281        // Now print etc
2282        if (iRow<500000) {
2283          for (int lColumn =jColumn;lColumn<kColumn;lColumn++) {
2284            iColumn = order[lColumn];
2285            CoinBigIndex start = columnStart[iColumn];
2286            if (model->logLevel()==63) {
2287              printf("column %d %g <= ",iColumn,columnLower[iColumn]);
2288              for (CoinBigIndex i=start;i<start+iRow;i++)
2289                printf("( %d, %g) ",row[i],element[i]);
2290              printf("<= %g\n",columnUpper[iColumn]);
2291            }
2292          }
2293        }
2294        jColumn =kColumn;
2295      }
2296    }
2297    delete [] order;
2298    delete [] other;
2299    delete [] number;
2300  }
2301  // get row copy
2302  CoinPackedMatrix rowCopy = *matrix;
2303  rowCopy.reverseOrdering();
2304  const int * rowLength = rowCopy.getVectorLengths();
2305  number = new int[numberColumns+1];
2306  memset(number,0,(numberColumns+1)*sizeof(int));
2307  for (iRow=0;iRow<numberRows;iRow++) {
2308    int length=rowLength[iRow];
2309    number[length]++;
2310  }
2311  if (number[0])
2312    printf("** %d rows have no entries\n",number[0]);
2313  k=0;
2314  for (iColumn=1;iColumn<=numberColumns;iColumn++) {
2315    if (number[iColumn]) {
2316      k++;
2317      printf("%d rows have %d entries\n",number[iColumn],iColumn);
2318      if (k==kMax)
2319        break;
2320    }
2321  }
2322  if (k<numberColumns) {
2323    int kk=k;
2324    k=0;
2325    for (iColumn=numberColumns;iColumn>=1;iColumn--) {
2326      if (number[iColumn]) {
2327        k++;
2328        if (k==kMax)
2329          break;
2330      }
2331    }
2332    if (k>kk) {
2333      printf("\n    .........\n\n");
2334      iColumn=k;
2335      k=0;
2336      for (;iColumn<numberColumns;iColumn++) {
2337        if (number[iColumn]) {
2338          k++;
2339          printf("%d rows have %d entries\n",number[iColumn],iColumn);
2340          if (k==kMax)
2341            break;
2342        }
2343      }
2344    }
2345  }
2346  if (model->logLevel()==63
2347#ifdef SYM
2348      ||true
2349#endif
2350      ) {
2351    int * column = rowCopy.getMutableIndices();
2352    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
2353    double * element = rowCopy.getMutableElements();
2354    int * order = new int[numberRows];
2355    int * other = new int[numberRows];
2356    for (iRow=0;iRow<numberRows;iRow++) {
2357      int length=rowLength[iRow];
2358      order[iRow]=iRow;
2359      other[iRow]=length;
2360      CoinBigIndex start = rowStart[iRow];
2361      CoinSort_2(column+start,column+start+length,element+start);
2362    }
2363    CoinSort_2(other,other+numberRows,order);
2364    int jRow=number[0]+number[1];
2365    double * weight = new double[numberRows];
2366    double * randomColumn = new double[numberColumns+1];
2367    double * randomRow = new double [numberRows+1];
2368    int * sortRow = new int [numberRows];
2369    int * possibleRow = new int [numberRows];
2370    int * backRow = new int [numberRows];
2371    int * stackRow = new int [numberRows];
2372    int * sortColumn = new int [numberColumns];
2373    int * possibleColumn = new int [numberColumns];
2374    int * backColumn = new int [numberColumns];
2375    int * backColumn2 = new int [numberColumns];
2376    int * mapRow = new int [numberRows];
2377    int * mapColumn = new int [numberColumns];
2378    int * stackColumn = new int [numberColumns];
2379    double randomLower = CoinDrand48();
2380    double randomUpper = CoinDrand48();
2381    double randomInteger = CoinDrand48();
2382    int * startAdd = new int[numberRows+1];
2383    int * columnAdd = new int [2*numberElements];
2384    double * elementAdd = new double[2*numberElements];
2385    int nAddRows=0;
2386    startAdd[0]=0;
2387    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2388      randomColumn[iColumn] = CoinDrand48();
2389      backColumn2[iColumn]=-1;
2390    }
2391    for (iColumn=2;iColumn<=numberColumns;iColumn++) {
2392      if (number[iColumn]) {
2393        printf("XX %d rows have %d entries\n",number[iColumn],iColumn);
2394        int kRow = jRow+number[iColumn];
2395        sortOnOther(column,rowStart,
2396                     order+jRow,other,number[iColumn],iColumn,0);
2397        // Now print etc
2398        if (iColumn<500000) {
2399          int nLook=0;
2400          for (int lRow =jRow;lRow<kRow;lRow++) {
2401            iRow = order[lRow];
2402            CoinBigIndex start = rowStart[iRow];
2403            if (model->logLevel()==63) {
2404              printf("row %d %g <= ",iRow,rowLower[iRow]);
2405              for (CoinBigIndex i=start;i<start+iColumn;i++) 
2406                printf("( %d, %g) ",column[i],element[i]);
2407              printf("<= %g\n",rowUpper[iRow]);
2408            }
2409            int first = column[start];
2410            double sum=0.0;
2411            for (CoinBigIndex i=start;i<start+iColumn;i++) {
2412              int jColumn = column[i];
2413              double value = element[i];
2414              jColumn -= first;
2415              assert (jColumn>=0);
2416              sum += value*randomColumn[jColumn];
2417            }
2418            if (rowLower[iRow]>-1.0e30&&rowLower[iRow])
2419              sum += rowLower[iRow]*randomLower;
2420            else if (!rowLower[iRow])
2421              sum += 1.234567e-7*randomLower;
2422            if (rowUpper[iRow]<1.0e30&&rowUpper[iRow])
2423              sum += rowUpper[iRow]*randomUpper;
2424            else if (!rowUpper[iRow])
2425              sum += 1.234567e-7*randomUpper;
2426            sortRow[nLook]=iRow;
2427            randomRow[nLook++]=sum;
2428            // best way is to number unique elements and bounds and use
2429            if (fabs(sum)>1.0e4)
2430              sum *= 1.0e-6;
2431            weight[iRow]=sum;
2432          }
2433          assert (nLook<=numberRows);
2434          CoinSort_2(randomRow,randomRow+nLook,sortRow);
2435          randomRow[nLook]=COIN_DBL_MAX;
2436          double last=-COIN_DBL_MAX;
2437          int iLast=-1;
2438          for (int iLook=0;iLook<nLook+1;iLook++) {
2439            if (randomRow[iLook]>last) {
2440              if (iLast>=0) {
2441                int n=iLook-iLast;
2442                if (n>1) {
2443                  //printf("%d rows possible?\n",n);
2444                }
2445              }
2446              iLast=iLook;
2447              last=randomRow[iLook];
2448            }
2449          }
2450        }
2451        jRow =kRow;
2452      }
2453    }
2454    CoinPackedMatrix columnCopy = *matrix;
2455    const int * columnLength = columnCopy.getVectorLengths();
2456    const int * row = columnCopy.getIndices();
2457    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2458    const double * elementByColumn = columnCopy.getElements();
2459    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2460      int length=columnLength[iColumn];
2461      CoinBigIndex start = columnStart[iColumn];
2462      double sum = objective[iColumn];
2463      if (columnLower[iColumn]>-1.0e30&&columnLower[iColumn])
2464        sum += columnLower[iColumn]*randomLower;
2465      else if (!columnLower[iColumn])
2466        sum += 1.234567e-7*randomLower;
2467      if (columnUpper[iColumn]<1.0e30&&columnUpper[iColumn])
2468        sum += columnUpper[iColumn]*randomUpper;
2469      else if (!columnUpper[iColumn])
2470        sum += 1.234567e-7*randomUpper;
2471      if (model->isInteger(iColumn))
2472        sum += 9.87654321e-6*randomInteger;
2473      for (CoinBigIndex i=start;i<start+length;i++) {
2474        int iRow = row[i];
2475        sum += elementByColumn[i]*weight[iRow];
2476      }
2477      sortColumn[iColumn]=iColumn;
2478      randomColumn[iColumn]=sum;
2479    }
2480    {
2481      CoinSort_2(randomColumn,randomColumn+numberColumns,sortColumn);
2482      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2483        int i=sortColumn[iColumn];
2484        backColumn[i]=iColumn;
2485      }
2486      randomColumn[numberColumns]=COIN_DBL_MAX;
2487      double last=-COIN_DBL_MAX;
2488      int iLast=-1;
2489      for (int iLook=0;iLook<numberColumns+1;iLook++) {
2490        if (randomColumn[iLook]>last) {
2491          if (iLast>=0) {
2492            int n=iLook-iLast;
2493            if (n>1) {
2494              //printf("%d columns possible?\n",n);
2495            }
2496            for (int i=iLast;i<iLook;i++) {
2497              possibleColumn[sortColumn[i]]=n;
2498            }
2499          }
2500          iLast=iLook;
2501          last=randomColumn[iLook];
2502        }
2503      }
2504      for (iRow =0;iRow<numberRows;iRow++) {
2505        CoinBigIndex start = rowStart[iRow];
2506        double sum=0.0;
2507        int length=rowLength[iRow];
2508        for (CoinBigIndex i=start;i<start+length;i++) {
2509          int jColumn = column[i];
2510          double value = element[i];
2511          jColumn=backColumn[jColumn];
2512          sum += value*randomColumn[jColumn];
2513          //if (iColumn==23089||iRow==23729)
2514          //printf("row %d cola %d colb %d value %g rand %g sum %g\n",
2515          //   iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
2516        }
2517        sortRow[iRow]=iRow;
2518        randomRow[iRow]=weight[iRow];
2519        randomRow[iRow]=sum;
2520      }
2521      CoinSort_2(randomRow,randomRow+numberRows,sortRow);
2522      for (iRow=0;iRow<numberRows;iRow++) {
2523        int i=sortRow[iRow];
2524        backRow[i]=iRow;
2525      }
2526      randomRow[numberRows]=COIN_DBL_MAX;
2527      last=-COIN_DBL_MAX;
2528      iLast=-1;
2529      // Do backward indices from order
2530      for (iRow=0;iRow<numberRows;iRow++) {
2531        other[order[iRow]]=iRow;
2532      }
2533      for (int iLook=0;iLook<numberRows+1;iLook++) {
2534        if (randomRow[iLook]>last) {
2535          if (iLast>=0) {
2536            int n=iLook-iLast;
2537            if (n>1) {
2538              //printf("%d rows possible?\n",n);
2539              // Within group sort as for original "order"
2540              for (int i=iLast;i<iLook;i++) {
2541                int jRow=sortRow[i];
2542                order[i]=other[jRow];
2543              }
2544              CoinSort_2(order+iLast,order+iLook,sortRow+iLast);
2545            }
2546            for (int i=iLast;i<iLook;i++) {
2547              possibleRow[sortRow[i]]=n;
2548            }
2549          }
2550          iLast=iLook;
2551          last=randomRow[iLook];
2552        }
2553      }
2554      // Temp out
2555      for (int iLook=0;iLook<numberRows-1000000;iLook++) {
2556        iRow=sortRow[iLook];
2557        CoinBigIndex start = rowStart[iRow];
2558        int length=rowLength[iRow];
2559        int numberPossible = possibleRow[iRow];
2560        for (CoinBigIndex i=start;i<start+length;i++) {
2561          int jColumn = column[i];
2562          if (possibleColumn[jColumn]!=numberPossible)
2563            numberPossible=-1;
2564        }
2565        int n=numberPossible;
2566        if (numberPossible>1) {
2567          //printf("pppppossible %d\n",numberPossible);
2568          for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2569            int jRow=sortRow[jLook];
2570            CoinBigIndex start2 = rowStart[jRow];
2571            assert (numberPossible==possibleRow[jRow]);
2572            assert(length==rowLength[jRow]);
2573            for (CoinBigIndex i=start2;i<start2+length;i++) {
2574              int jColumn = column[i];
2575              if (possibleColumn[jColumn]!=numberPossible)
2576                numberPossible=-1;
2577            }
2578          }
2579          if (numberPossible<2) {
2580            // switch off
2581            for (int jLook=iLook;jLook<iLook+n;jLook++) 
2582              possibleRow[sortRow[jLook]]=-1;
2583          }
2584          // skip rest
2585          iLook += n-1;
2586        } else {
2587          possibleRow[iRow]=-1;
2588        }
2589      }
2590      for (int iLook=0;iLook<numberRows;iLook++) {
2591        iRow=sortRow[iLook];
2592        int numberPossible = possibleRow[iRow];
2593        // Only if any integers
2594        int numberIntegers=0;
2595        CoinBigIndex start = rowStart[iRow];
2596        int length=rowLength[iRow];
2597        for (CoinBigIndex i=start;i<start+length;i++) {
2598          int jColumn = column[i];
2599          if (model->isInteger(jColumn))
2600            numberIntegers++;
2601        }
2602        if (numberPossible>1&&!numberIntegers) {
2603          //printf("possible %d - but no integers\n",numberPossible);
2604        }
2605        if (numberPossible>1&&(numberIntegers||true)) {
2606          //
2607          printf("possible %d - %d integers\n",numberPossible,numberIntegers);
2608          int lastLook=iLook;
2609          int nMapRow=-1;
2610          for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2611            // stop if too many failures
2612            if (jLook>iLook+10&&nMapRow<0)
2613              break;
2614            // Create identity mapping
2615            int i;
2616            for (i=0;i<numberRows;i++)
2617              mapRow[i]=i;
2618            for (i=0;i<numberColumns;i++)
2619              mapColumn[i]=i;
2620            int offset=jLook-iLook;
2621            int nStackC=0;
2622            // build up row and column mapping
2623            int nStackR=1;
2624            stackRow[0]=iLook;
2625            bool good=true;
2626            while (nStackR) {
2627              nStackR--;
2628              int look1 = stackRow[nStackR];
2629              int look2 = look1 + offset;
2630              assert (randomRow[look1]==randomRow[look2]);
2631              int row1=sortRow[look1];
2632              int row2=sortRow[look2];
2633              assert (mapRow[row1]==row1);
2634              assert (mapRow[row2]==row2);
2635              mapRow[row1]=row2;
2636              mapRow[row2]=row1;
2637              CoinBigIndex start1 = rowStart[row1];
2638              CoinBigIndex offset2 = rowStart[row2]-start1;
2639              int length=rowLength[row1];
2640              assert( length==rowLength[row2]);
2641              for (CoinBigIndex i=start1;i<start1+length;i++) {
2642                int jColumn1 = column[i];
2643                int jColumn2 = column[i+offset2];
2644                if (randomColumn[backColumn[jColumn1]]!=
2645                    randomColumn[backColumn[jColumn2]]) {
2646                  good=false;
2647                  break;
2648                }
2649                if (mapColumn[jColumn1]==jColumn1) {
2650                  // not touched
2651                  assert (mapColumn[jColumn2]==jColumn2);
2652                  if (jColumn1!=jColumn2) {
2653                    // Put on stack
2654                    mapColumn[jColumn1]=jColumn2;
2655                    mapColumn[jColumn2]=jColumn1;
2656                    stackColumn[nStackC++]=jColumn1;
2657                  }
2658                } else {
2659                  if(mapColumn[jColumn1]!=jColumn2||
2660                     mapColumn[jColumn2]!=jColumn1) {
2661                    // bad
2662                    good=false;
2663                    printf("bad col\n");
2664                    break;
2665                  }
2666                }
2667              }
2668              if (!good)
2669                break;
2670              while (nStackC) {
2671                nStackC--;
2672                int iColumn = stackColumn[nStackC];
2673                int iColumn2 = mapColumn[iColumn];
2674                assert (iColumn!=iColumn2);
2675                int length=columnLength[iColumn];
2676                assert (length==columnLength[iColumn2]);
2677                CoinBigIndex start = columnStart[iColumn];
2678                CoinBigIndex offset2 = columnStart[iColumn2]-start;
2679                for (CoinBigIndex i=start;i<start+length;i++) {
2680                  int iRow = row[i];
2681                  int iRow2 = row[i+offset2];
2682                  if (mapRow[iRow]==iRow) {
2683                    // First (but be careful)
2684                    if (iRow!=iRow2) {
2685                      //mapRow[iRow]=iRow2;
2686                      //mapRow[iRow2]=iRow;
2687                      int iBack=backRow[iRow];
2688                      int iBack2=backRow[iRow2];
2689                      if (randomRow[iBack]==randomRow[iBack2]&&
2690                          iBack2-iBack==offset) {
2691                        stackRow[nStackR++]=iBack;
2692                      } else {
2693                        //printf("randomRow diff - weights %g %g\n",
2694                        //     weight[iRow],weight[iRow2]);
2695                        // bad
2696                        good=false;
2697                        break;
2698                      }
2699                    }
2700                  } else {
2701                    if(mapRow[iRow]!=iRow2||
2702                       mapRow[iRow2]!=iRow) {
2703                      // bad
2704                      good=false;
2705                      printf("bad row\n");
2706                      break;
2707                    }
2708                  }
2709                }
2710                if (!good)
2711                  break;
2712              }
2713            }
2714            // then check OK
2715            if (good) {
2716              for (iRow =0;iRow<numberRows;iRow++) {
2717                CoinBigIndex start = rowStart[iRow];
2718                int length=rowLength[iRow];
2719                if (mapRow[iRow]==iRow) {
2720                  for (CoinBigIndex i=start;i<start+length;i++) {
2721                    int jColumn = column[i];
2722                    backColumn2[jColumn]=i-start;
2723                  }
2724                  for (CoinBigIndex i=start;i<start+length;i++) {
2725                    int jColumn = column[i];
2726                    if (mapColumn[jColumn]!=jColumn) {
2727                      int jColumn2 =mapColumn[jColumn];
2728                      CoinBigIndex i2 = backColumn2[jColumn2];
2729                      if (i2<0) {
2730                        good=false;
2731                      } else if (element[i]!=element[i2+start]) {
2732                        good=false;
2733                      }
2734                    }
2735                  }
2736                  for (CoinBigIndex i=start;i<start+length;i++) {
2737                    int jColumn = column[i];
2738                    backColumn2[jColumn]=-1;
2739                  }
2740                } else {
2741                  int row2 = mapRow[iRow];
2742                  assert (iRow = mapRow[row2]);
2743                  if (rowLower[iRow]!=rowLower[row2]||
2744                      rowLower[row2]!=rowLower[iRow])
2745                    good=false;
2746                  CoinBigIndex offset2 = rowStart[row2]-start;
2747                  for (CoinBigIndex i=start;i<start+length;i++) {
2748                    int jColumn = column[i];
2749                    double value = element[i];
2750                    int jColumn2 = column[i+offset2];
2751                    double value2 = element[i+offset2];
2752                    if (value!=value2||mapColumn[jColumn]!=jColumn2||
2753                        mapColumn[jColumn2]!=jColumn)
2754                      good=false;
2755                  }
2756                }
2757              }
2758              if (good) {
2759                // check rim
2760                for (iColumn=0;iColumn<numberColumns;iColumn++) {
2761                  if (mapColumn[iColumn]!=iColumn) {
2762                    int iColumn2 = mapColumn[iColumn];
2763                    if(objective[iColumn]!=objective[iColumn2])
2764                      good=false;
2765                    if (columnLower[iColumn]!=columnLower[iColumn2])
2766                      good=false;
2767                    if (columnUpper[iColumn]!=columnUpper[iColumn2])
2768                      good=false;
2769                    if (model->isInteger(iColumn)!=model->isInteger(iColumn2))
2770                      good=false;
2771                  }
2772                }
2773              }
2774              if (good) {
2775                // temp
2776                if (nMapRow<0) {
2777                  //const double * solution = model->primalColumnSolution();
2778                  // find mapped
2779                  int nMapColumn=0;
2780                  for (int i=0;i<numberColumns;i++) {
2781                    if (mapColumn[i]>i) 
2782                      nMapColumn++;
2783                  }
2784                  nMapRow=0;
2785                  int kRow=-1;
2786                  for (int i=0;i<numberRows;i++) {
2787                    if (mapRow[i]>i) { 
2788                      nMapRow++;
2789                      kRow=i;
2790                    }
2791                  }
2792                  printf("%d columns, %d rows\n",nMapColumn,nMapRow);
2793                  if (nMapRow==1) {
2794                    CoinBigIndex start = rowStart[kRow];
2795                    int length=rowLength[kRow];
2796                    printf("%g <= ",rowLower[kRow]);
2797                    for (CoinBigIndex i=start;i<start+length;i++) {
2798                      int jColumn = column[i];
2799                      if (mapColumn[jColumn]!=jColumn) 
2800                        printf("* ");
2801                      printf("%d,%g ",jColumn,element[i]);
2802                    }
2803                    printf("<= %g\n",rowUpper[kRow]);
2804                  }
2805                }
2806                // temp
2807                int row1=sortRow[lastLook];
2808                int row2=sortRow[jLook];
2809                lastLook=jLook;
2810                CoinBigIndex start1 = rowStart[row1];
2811                CoinBigIndex offset2 = rowStart[row2]-start1;
2812                int length=rowLength[row1];
2813                assert( length==rowLength[row2]);
2814                CoinBigIndex put=startAdd[nAddRows];
2815                double multiplier=length<11 ? 2.0 : 1.125;
2816                double value=1.0;
2817                for (CoinBigIndex i=start1;i<start1+length;i++) {
2818                  int jColumn1 = column[i];
2819                  int jColumn2 = column[i+offset2];
2820                  columnAdd[put]=jColumn1;
2821                  elementAdd[put++]=value;
2822                  columnAdd[put]=jColumn2;
2823                  elementAdd[put++]=-value;
2824                  value *= multiplier;
2825                }
2826                nAddRows++;
2827                startAdd[nAddRows]=put;
2828              } else {
2829                printf("ouch - did not check out as good\n");
2830              }
2831            }
2832          }
2833          // skip rest
2834          iLook += numberPossible-1;
2835        }
2836      }
2837    }
2838    if (nAddRows) {
2839      double * lower = new double [nAddRows];
2840      double * upper = new double[nAddRows];
2841      int i;
2842      //const double * solution = model->primalColumnSolution();
2843      for (i=0;i<nAddRows;i++) {
2844        lower[i]=0.0;
2845        upper[i]=COIN_DBL_MAX;
2846      }
2847      printf("Adding %d rows with %d elements\n",nAddRows,
2848             startAdd[nAddRows]);
2849      //ClpSimplex newModel(*model);
2850      //newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
2851      //newModel.writeMps("modified.mps");
2852      delete [] lower;
2853      delete [] upper;
2854    }
2855    delete [] startAdd;
2856    delete [] columnAdd;
2857    delete [] elementAdd;
2858    delete [] order;
2859    delete [] other;
2860    delete [] randomColumn;
2861    delete [] weight;
2862    delete [] randomRow;
2863    delete [] sortRow;
2864    delete [] backRow;
2865    delete [] possibleRow;
2866    delete [] sortColumn;
2867    delete [] backColumn;
2868    delete [] backColumn2;
2869    delete [] possibleColumn;
2870    delete [] mapRow;
2871    delete [] mapColumn;
2872    delete [] stackRow;
2873    delete [] stackColumn;
2874  }
2875  delete [] number;
2876  // Now do breakdown of ranges
2877  breakdown("Elements",numberElements,elementByColumn);
2878  breakdown("RowLower",numberRows,rowLower);
2879  breakdown("RowUpper",numberRows,rowUpper);
2880  breakdown("ColumnLower",numberColumns,columnLower);
2881  breakdown("ColumnUpper",numberColumns,columnUpper);
2882  breakdown("Objective",numberColumns,objective);
2883}
2884static bool maskMatches(const int * starts, char ** masks,
2885                        std::string & check)
2886{
2887  // back to char as I am old fashioned
2888  const char * checkC = check.c_str();
2889  int length = strlen(checkC);
2890  while (checkC[length-1]==' ')
2891    length--;
2892  for (int i=starts[length];i<starts[length+1];i++) {
2893    char * thisMask = masks[i];
2894    int k;
2895    for ( k=0;k<length;k++) {
2896      if (thisMask[k]!='?'&&thisMask[k]!=checkC[k]) 
2897        break;
2898    }
2899    if (k==length)
2900      return true;
2901  }
2902  return false;
2903}
2904static void clean(char * temp)
2905{
2906  char * put = temp;
2907  while (*put>=' ')
2908    put++;
2909  *put='\0';
2910}
2911static void generateCode(const char * fileName,int type)
2912{
2913  FILE * fp = fopen(fileName,"r");
2914  assert (fp);
2915  int numberLines=0;
2916#define MAXLINES 500
2917#define MAXONELINE 200
2918  char line[MAXLINES][MAXONELINE];
2919  while (fgets(line[numberLines],MAXONELINE,fp)) {
2920    assert (numberLines<MAXLINES);
2921    clean(line[numberLines]);
2922    numberLines++;
2923  }
2924  fclose(fp);
2925  // add in actual solve
2926  strcpy(line[numberLines],"5  clpModel->initialSolve(clpSolve);");
2927  numberLines++;
2928  fp = fopen(fileName,"w");
2929  assert (fp);
2930  char apo='"';   
2931  char backslash = '\\';
2932
2933  fprintf(fp,"#include %cClpSimplex.hpp%c\n",apo,apo);
2934  fprintf(fp,"#include %cClpSolve.hpp%c\n",apo,apo);
2935
2936  fprintf(fp,"\nint main (int argc, const char *argv[])\n{\n");
2937  fprintf(fp,"  ClpSimplex  model;\n");
2938  fprintf(fp,"  int status=1;\n");
2939  fprintf(fp,"  if (argc<2)\n");
2940  fprintf(fp,"    fprintf(stderr,%cPlease give file name%cn%c);\n",
2941          apo,backslash,apo);
2942  fprintf(fp,"  else\n");
2943  fprintf(fp,"    status=model.readMps(argv[1],true);\n");
2944  fprintf(fp,"  if (status) {\n");
2945  fprintf(fp,"    fprintf(stderr,%cBad readMps %%s%cn%c,argv[1]);\n",
2946                apo,backslash,apo);
2947  fprintf(fp,"    exit(1);\n");
2948  fprintf(fp,"  }\n\n");
2949  fprintf(fp,"  // Now do requested saves and modifications\n");
2950  fprintf(fp,"  ClpSimplex * clpModel = & model;\n");
2951  int wanted[9];
2952  memset(wanted,0,sizeof(wanted));
2953  wanted[0]=wanted[3]=wanted[5]=wanted[8]=1;
2954  if (type>0) 
2955    wanted[1]=wanted[6]=1;
2956  if (type>1) 
2957    wanted[2]=wanted[4]=wanted[7]=1;
2958  std::string header[9]=
2959  { "","Save values","Redundant save of default values","Set changed values",
2960    "Redundant set default values","Solve","Restore values","Redundant restore values","Add to model"};
2961  for (int iType=0;iType<9;iType++) {
2962    if (!wanted[iType])
2963      continue;
2964    int n=0;
2965    int iLine;
2966    for (iLine=0;iLine<numberLines;iLine++) {
2967      if (line[iLine][0]=='0'+iType) {
2968        if (!n)
2969          fprintf(fp,"\n  // %s\n\n",header[iType].c_str());
2970        n++;
2971        fprintf(fp,"%s\n",line[iLine]+1);
2972      }
2973    }
2974  }
2975  fprintf(fp,"\n  // Now you would use solution etc etc\n\n");
2976  fprintf(fp,"  return 0;\n}\n");
2977  fclose(fp);
2978  printf("C++ file written to %s\n",fileName);
2979}
2980/*
2981  Version 1.00.00 October 13 2004.
2982  1.00.01 October 18.  Added basis handling helped/prodded by Thorsten Koch.
2983  Also modifications to make faster with sbb (I hope I haven't broken anything).
2984  1.00.02 March 21 2005.  Redid ClpNonLinearCost to save memory also redid
2985  createRim to try and improve cache characteristics.
2986  1.00.03 April 8 2005.  Added Volume algorithm as crash and made code more
2987  robust on testing.  Also added "either" and "tune" algorithm.
2988  1.01.01 April 12 2005.  Decided to go to different numbering.  Backups will
2989  be last 2 digits while middle 2 are for improvements.  Still take a long
2990  time to get to 2.00.01
2991  1.01.02 May 4 2005.  Will be putting in many changes - so saving stable version
2992  1.02.01 May 6 2005.  Lots of changes to try and make faster and more stable in
2993  branch and cut.
2994  1.02.02 May 19 2005.  Stuff for strong branching and some improvements to simplex
2995  1.03.01 May 24 2006.  Lots done but I can't remember what!
2996  1.03.03 June 13 2006.  For clean up after dual perturbation
2997  1.04.01 June 26 2007.  Lots of changes but I got lazy
2998  1.05.00 June 27 2007.  This is trunk so when gets to stable will be 1.5
2999 */
Note: See TracBrowser for help on using the repository browser.