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

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

changes to simplex and lots of stuff and start Mumps cholesky

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