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

Last change on this file since 1379 was 1379, checked in by forrest, 11 years ago

add space

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