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

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

allow use of dense factorization

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