source: branches/sandbox/Clp/src/ClpMain.cpp @ 1474

Last change on this file since 1474 was 1474, checked in by bjarni, 11 years ago

Renamed parameter constants in CbcOrClpParam? and ClpMain? to make them more readable/search-able

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