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

Last change on this file since 1402 was 1402, checked in by forrest, 12 years ago

get rid of compiler warnings

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