source: stable/1.8/Clp/src/ClpMain.cpp @ 1236

Last change on this file since 1236 was 1236, checked in by ladanyi, 12 years ago

force cdecl on signal_handler() and main() if MS compiler is used (ticket #20)

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