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

Last change on this file since 1421 was 1421, checked in by stefan, 10 years ago

let configure setup CLPVERSION

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