source: trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp @ 1584

Last change on this file since 1584 was 1584, checked in by lou, 10 years ago

Change signature of getDualRays from (int) to (int,bool) to allow choice of
partial (row components) or full (row and column components) dual ray. Matches
Osi trunk revision 1551.

File size: 289.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#include <cassert>
4#ifdef CBC_STATISTICS
5extern int osi_crunch;
6extern int osi_primal;
7extern int osi_dual;
8extern int osi_hot;
9#endif
10#include "CoinTime.hpp"
11#include "CoinHelperFunctions.hpp"
12#include "CoinIndexedVector.hpp"
13#include "CoinModel.hpp"
14#include "CoinMpsIO.hpp"
15#include "CoinSort.hpp"
16#include "ClpDualRowSteepest.hpp"
17#include "ClpPrimalColumnSteepest.hpp"
18#include "ClpPackedMatrix.hpp"
19#include "ClpDualRowDantzig.hpp"
20#include "ClpPrimalColumnDantzig.hpp"
21#include "ClpFactorization.hpp"
22#include "ClpObjective.hpp"
23#include "ClpSimplex.hpp"
24#include "ClpSimplexOther.hpp"
25#include "ClpSimplexPrimal.hpp"
26#include "ClpSimplexDual.hpp"
27#include "ClpNonLinearCost.hpp"
28#include "OsiClpSolverInterface.hpp"
29#include "OsiBranchingObject.hpp"
30#include "OsiCuts.hpp"
31#include "OsiRowCut.hpp"
32#include "OsiColCut.hpp"
33#include "ClpPresolve.hpp"
34#include "CoinLpIO.hpp"
35static double totalTime=0.0;
36//#define SAVE_MODEL 1
37#ifdef SAVE_MODEL
38static int resolveTry=0;
39static int loResolveTry=0;
40static int hiResolveTry=9999999;
41#endif
42//#############################################################################
43// Solve methods
44//#############################################################################
45void OsiClpSolverInterface::initialSolve()
46{
47#define KEEP_SMALL
48#ifdef KEEP_SMALL
49  if (smallModel_) {
50    delete [] spareArrays_;
51    spareArrays_ = NULL;
52    delete smallModel_;
53    smallModel_=NULL;
54  }
55#endif
56  bool deleteSolver;
57  ClpSimplex * solver;
58  double time1 = CoinCpuTime();
59  int userFactorizationFrequency = modelPtr_->factorization()->maximumPivots();
60  int totalIterations=0;
61  bool abortSearch=false;
62  ClpObjective * savedObjective=NULL;
63  double savedDualLimit=modelPtr_->dblParam_[ClpDualObjectiveLimit];
64  if (fakeObjective_) {
65    modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~128));
66    // See if all with costs fixed
67    int numberColumns = modelPtr_->numberColumns_;
68    const double * obj = modelPtr_->objective();
69    const double * lower = modelPtr_->columnLower();
70    const double * upper = modelPtr_->columnUpper();
71    int i;
72    for (i=0;i<numberColumns;i++) {
73      double objValue = obj[i];
74      if (objValue) {
75        if (lower[i]!=upper[i])
76          break;
77      }
78    }
79    if (i==numberColumns) {
80      if ((specialOptions_&524288)==0) {
81        // Set fake
82        savedObjective=modelPtr_->objective_;
83        modelPtr_->objective_=fakeObjective_;
84        modelPtr_->dblParam_[ClpDualObjectiveLimit]=COIN_DBL_MAX;
85      } else {
86        modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()|128);
87      }
88    }
89  }
90  if ((specialOptions_&1024)==0) {
91    solver = new ClpSimplex(true);
92    deleteSolver=true;
93    solver->borrowModel(*modelPtr_);
94    // See if user set factorization frequency
95    // borrowModel does not move
96    solver->factorization()->maximumPivots(userFactorizationFrequency);
97  } else {
98    solver = modelPtr_;
99    deleteSolver=false;
100  }
101  // Treat as if user simplex not enabled
102  int saveSolveType=solver->solveType();
103  bool doingPrimal = solver->algorithm()>0;
104  if (saveSolveType==2) {
105    disableSimplexInterface();
106    solver->setSolveType(1);
107  }
108  int saveOptions = solver->specialOptions();
109  solver->setSpecialOptions(saveOptions|64|32768); // go as far as possible
110  // get original log levels
111  int saveMessageLevel=modelPtr_->logLevel();
112  int messageLevel=messageHandler()->logLevel();
113  int saveMessageLevel2 = messageLevel;
114  // Set message handler
115  solver->passInMessageHandler(handler_);
116  // But keep log level
117  solver->messageHandler()->setLogLevel(saveMessageLevel);
118  // set reasonable defaults
119  bool takeHint;
120  OsiHintStrength strength;
121  // Switch off printing if asked to
122  bool gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength));
123  assert (gotHint);
124  if (strength!=OsiHintIgnore&&takeHint) {
125    if (messageLevel>0)
126      messageLevel--;
127  }
128  if (messageLevel<saveMessageLevel)
129    solver->messageHandler()->setLogLevel(messageLevel);
130  // Allow for specialOptions_==1+8 forcing saving factorization
131  int startFinishOptions=0;
132  if ((specialOptions_&9)==(1+8)) {
133    startFinishOptions =1+2+4; // allow re-use of factorization
134  }
135  bool defaultHints=true;
136  {
137    int hint;
138    for (hint=OsiDoPresolveInInitial;hint<OsiLastHintParam;hint++) {
139      if (hint!=OsiDoReducePrint&&
140          hint!=OsiDoInBranchAndCut) {
141        bool yesNo;
142        OsiHintStrength strength;
143        getHintParam(static_cast<OsiHintParam> (hint),yesNo,strength);
144        if (yesNo) {
145          defaultHints=false;
146          break;
147        }
148        if (strength != OsiHintIgnore) {
149          defaultHints=false;
150          break;
151        }
152      }
153    }
154  }
155  /*
156    If basis then do primal (as user could do dual with resolve)
157    If not then see if dual feasible (and allow for gubs etc?)
158  */
159  bool doPrimal = (basis_.numberBasicStructurals()>0);
160  setBasis(basis_,solver);
161  bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
162  if ((!defaultHints||doPrimal)&&!solveOptions_.getSpecialOption(6)) {
163    // scaling
164    // save initial state
165    const double * rowScale1 = solver->rowScale();
166    if (modelPtr_->solveType()==1) {
167      gotHint = (getHintParam(OsiDoScale,takeHint,strength));
168      assert (gotHint);
169      if (strength==OsiHintIgnore||takeHint) {
170        if (!solver->scalingFlag())
171          solver->scaling(3);
172      } else {
173        solver->scaling(0);
174      }
175    } else {
176      solver->scaling(0);
177    }
178    //solver->setDualBound(1.0e6);
179    //solver->setDualTolerance(1.0e-7);
180   
181    //ClpDualRowSteepest steep;
182    //solver->setDualRowPivotAlgorithm(steep);
183    //solver->setPrimalTolerance(1.0e-8);
184    //ClpPrimalColumnSteepest steepP;
185    //solver->setPrimalColumnPivotAlgorithm(steepP);
186   
187    // sort out hints;
188    // algorithm 0 whatever, -1 force dual, +1 force primal
189    int algorithm = 0;
190    gotHint = (getHintParam(OsiDoDualInInitial,takeHint,strength));
191    assert (gotHint);
192    if (strength!=OsiHintIgnore)
193      algorithm = takeHint ? -1 : 1;
194    // crash 0 do lightweight if all slack, 1 do, -1 don't
195    int doCrash=0;
196    gotHint = (getHintParam(OsiDoCrash,takeHint,strength));
197    assert (gotHint);
198    if (strength!=OsiHintIgnore)
199      doCrash = takeHint ? 1 : -1;
200    // doPrimal set true if any structurals in basis so switch off crash
201    if (doPrimal)
202      doCrash = -1;
203   
204    // presolve
205    gotHint = (getHintParam(OsiDoPresolveInInitial,takeHint,strength));
206    assert (gotHint);
207    if (strength!=OsiHintIgnore&&takeHint) {
208      ClpPresolve pinfo;
209      ClpSimplex * model2 = pinfo.presolvedModel(*solver,1.0e-8);
210      if (!model2) {
211        // problem found to be infeasible - whats best?
212        model2 = solver;
213      } else {
214        model2->setSpecialOptions(solver->specialOptions());
215      }
216     
217      // change from 200 (unless changed)
218      if (modelPtr_->factorization()->maximumPivots()==200)
219        model2->factorization()->maximumPivots(100+model2->numberRows()/50);
220      else
221        model2->factorization()->maximumPivots(userFactorizationFrequency);
222      int savePerturbation = model2->perturbation();
223      if (savePerturbation==100)
224        model2->setPerturbation(50);
225      if (!doPrimal) {
226        // faster if bounds tightened
227        //int numberInfeasibilities = model2->tightenPrimalBounds();
228        model2->tightenPrimalBounds();
229        // look further
230        bool crashResult=false;
231        if (doCrash>0)
232          crashResult =  (solver->crash(1000.0,1)>0);
233        else if (doCrash==0&&algorithm>0)
234          crashResult =  (solver->crash(1000.0,1)>0);
235        doPrimal=crashResult;
236      }
237      if (algorithm<0)
238        doPrimal=false;
239      else if (algorithm>0)
240        doPrimal=true;
241      if (!doPrimal) {
242        //if (numberInfeasibilities)
243        //std::cout<<"** Analysis indicates model infeasible"
244        //       <<std::endl;
245        // up dual bound for safety
246        //model2->setDualBound(1.0e11);
247        disasterHandler_->setOsiModel(this);
248        if (inCbcOrOther) {
249          disasterHandler_->setSimplex(model2);
250          disasterHandler_->setWhereFrom(4);
251          model2->setDisasterHandler(disasterHandler_);
252        }
253        model2->dual(0);
254        totalIterations += model2->numberIterations();
255        if (inCbcOrOther) {
256          if(disasterHandler_->inTrouble()) {
257#ifdef COIN_DEVELOP
258            printf("dual trouble a\n");
259#endif
260            if (disasterHandler_->typeOfDisaster()) {
261              // We want to abort
262              abortSearch=true;
263              goto disaster;
264            }
265            // try just going back in
266            disasterHandler_->setPhase(1);
267            model2->dual();
268            totalIterations += model2->numberIterations();
269            if (disasterHandler_->inTrouble()) {
270#ifdef COIN_DEVELOP
271              printf("dual trouble b\n");
272#endif
273              if (disasterHandler_->typeOfDisaster()) {
274                // We want to abort
275                abortSearch=true;
276                goto disaster;
277              }
278              // try primal with original basis
279              disasterHandler_->setPhase(2);
280              setBasis(basis_,model2);
281              model2->primal();
282              totalIterations += model2->numberIterations();
283            }
284            if(disasterHandler_->inTrouble()) {
285#ifdef COIN_DEVELOP
286              printf("disaster - treat as infeasible\n");
287#endif
288              if (disasterHandler_->typeOfDisaster()) {
289                // We want to abort
290                abortSearch=true;
291                goto disaster;
292              }
293              model2->setProblemStatus(1);
294            }
295          }
296          // reset
297          model2->setDisasterHandler(NULL);
298        }
299        // check if clp thought it was in a loop
300        if (model2->status()==3&&!model2->hitMaximumIterations()) {
301          // switch algorithm
302          disasterHandler_->setOsiModel(this);
303          if (inCbcOrOther) {
304            disasterHandler_->setSimplex(model2);
305            disasterHandler_->setWhereFrom(6);
306            model2->setDisasterHandler(disasterHandler_);
307          }
308          model2->primal();
309          totalIterations += model2->numberIterations();
310          if (inCbcOrOther) {
311            if(disasterHandler_->inTrouble()) {
312#ifdef COIN_DEVELOP
313              printf("primal trouble a\n");
314#endif
315              if (disasterHandler_->typeOfDisaster()) {
316                // We want to abort
317                abortSearch=true;
318                goto disaster;
319              }
320              // try just going back in (but with dual)
321              disasterHandler_->setPhase(1);
322              model2->dual();
323              totalIterations += model2->numberIterations();
324              if (disasterHandler_->inTrouble()) {
325#ifdef COIN_DEVELOP
326                printf("primal trouble b\n");
327#endif
328                if (disasterHandler_->typeOfDisaster()) {
329                  // We want to abort
330                  abortSearch=true;
331                  goto disaster;
332                }
333                // try primal with original basis
334                disasterHandler_->setPhase(2);
335                setBasis(basis_,model2);
336                model2->dual();
337                totalIterations += model2->numberIterations();
338              }
339              if(disasterHandler_->inTrouble()) {
340#ifdef COIN_DEVELOP
341                printf("disaster - treat as infeasible\n");
342#endif
343                if (disasterHandler_->typeOfDisaster()) {
344                  // We want to abort
345                  abortSearch=true;
346                  goto disaster;
347                }
348                model2->setProblemStatus(1);
349              }
350            }
351            // reset
352            model2->setDisasterHandler(NULL);
353          }
354        }
355      } else {
356        // up infeasibility cost for safety
357        //model2->setInfeasibilityCost(1.0e10);
358        disasterHandler_->setOsiModel(this);
359        if (inCbcOrOther) {
360          disasterHandler_->setSimplex(model2);
361          disasterHandler_->setWhereFrom(6);
362          model2->setDisasterHandler(disasterHandler_);
363        }
364        model2->primal(1);
365        totalIterations += model2->numberIterations();
366        if (inCbcOrOther) {
367          if(disasterHandler_->inTrouble()) {
368#ifdef COIN_DEVELOP
369            printf("primal trouble a\n");
370#endif
371            if (disasterHandler_->typeOfDisaster()) {
372              // We want to abort
373              abortSearch=true;
374              goto disaster;
375            }
376            // try just going back in (but with dual)
377            disasterHandler_->setPhase(1);
378            model2->dual();
379            totalIterations += model2->numberIterations();
380            if (disasterHandler_->inTrouble()) {
381#ifdef COIN_DEVELOP
382              printf("primal trouble b\n");
383#endif
384              if (disasterHandler_->typeOfDisaster()) {
385                // We want to abort
386                abortSearch=true;
387                goto disaster;
388              }
389              // try primal with original basis
390              disasterHandler_->setPhase(2);
391              setBasis(basis_,model2);
392              model2->dual();
393              totalIterations += model2->numberIterations();
394            }
395            if(disasterHandler_->inTrouble()) {
396#ifdef COIN_DEVELOP
397              printf("disaster - treat as infeasible\n");
398#endif
399              if (disasterHandler_->typeOfDisaster()) {
400                // We want to abort
401                abortSearch=true;
402                goto disaster;
403              }
404              model2->setProblemStatus(1);
405            }
406          }
407          // reset
408          model2->setDisasterHandler(NULL);
409        }
410        // check if clp thought it was in a loop
411        if (model2->status()==3&&!model2->hitMaximumIterations()) {
412          // switch algorithm
413          disasterHandler_->setOsiModel(this);
414          if (inCbcOrOther) {
415            disasterHandler_->setSimplex(model2);
416            disasterHandler_->setWhereFrom(4);
417            model2->setDisasterHandler(disasterHandler_);
418          }
419          model2->dual(0);
420          totalIterations += model2->numberIterations();
421          if (inCbcOrOther) {
422            if(disasterHandler_->inTrouble()) {
423#ifdef COIN_DEVELOP
424              printf("dual trouble a\n");
425#endif
426              if (disasterHandler_->typeOfDisaster()) {
427                // We want to abort
428                abortSearch=true;
429                goto disaster;
430              }
431              // try just going back in
432              disasterHandler_->setPhase(1);
433              model2->dual();
434              totalIterations += model2->numberIterations();
435              if (disasterHandler_->inTrouble()) {
436#ifdef COIN_DEVELOP
437                printf("dual trouble b\n");
438#endif
439                if (disasterHandler_->typeOfDisaster()) {
440                  // We want to abort
441                  abortSearch=true;
442                goto disaster;
443                }
444                // try primal with original basis
445                disasterHandler_->setPhase(2);
446                setBasis(basis_,model2);
447                model2->primal();
448                totalIterations += model2->numberIterations();
449              }
450              if(disasterHandler_->inTrouble()) {
451#ifdef COIN_DEVELOP
452                printf("disaster - treat as infeasible\n");
453#endif
454                if (disasterHandler_->typeOfDisaster()) {
455                  // We want to abort
456                  abortSearch=true;
457                  goto disaster;
458                }
459                model2->setProblemStatus(1);
460              }
461            }
462            // reset
463            model2->setDisasterHandler(NULL);
464          }
465        }
466      }
467      model2->setPerturbation(savePerturbation);
468      if (model2!=solver) {
469        int presolvedStatus = model2->status();
470        pinfo.postsolve(true);
471       
472        delete model2;
473        //printf("Resolving from postsolved model\n");
474        // later try without (1) and check duals before solve
475        if (presolvedStatus!=3) {
476          if (!inCbcOrOther||presolvedStatus!=1) {
477            disasterHandler_->setOsiModel(this);
478            if (inCbcOrOther) {
479              disasterHandler_->setSimplex(solver); // as "borrowed"
480              disasterHandler_->setWhereFrom(6);
481              solver->setDisasterHandler(disasterHandler_);
482            }
483            solver->primal(1);
484            totalIterations += solver->numberIterations();
485            if (inCbcOrOther) {
486              if(disasterHandler_->inTrouble()) {
487#ifdef COIN_DEVELOP
488                printf("primal trouble a\n");
489#endif
490                if (disasterHandler_->typeOfDisaster()) {
491                  // We want to abort
492                  abortSearch=true;
493                  goto disaster;
494                }
495                // try just going back in (but with dual)
496                disasterHandler_->setPhase(1);
497                solver->dual();
498                totalIterations += solver->numberIterations();
499                if (disasterHandler_->inTrouble()) {
500#ifdef COIN_DEVELOP
501                  printf("primal trouble b\n");
502#endif
503                  if (disasterHandler_->typeOfDisaster()) {
504                    // We want to abort
505                    abortSearch=true;
506                    goto disaster;
507                  }
508                  // try primal with original basis
509                  disasterHandler_->setPhase(2);
510                  setBasis(basis_,solver);
511                  solver->dual();
512                  totalIterations += solver->numberIterations();
513                }
514                if(disasterHandler_->inTrouble()) {
515#ifdef COIN_DEVELOP
516                  printf("disaster - treat as infeasible\n");
517#endif
518                  if (disasterHandler_->typeOfDisaster()) {
519                    // We want to abort
520                    abortSearch=true;
521                    goto disaster;
522                  }
523                  solver->setProblemStatus(1);
524                }
525              }
526              // reset
527              solver->setDisasterHandler(NULL);
528            }
529          }
530        }
531      }
532      lastAlgorithm_=1; // primal
533      //if (solver->numberIterations())
534      //printf("****** iterated %d\n",solver->numberIterations());
535    } else {
536      // do we want crash
537      if (doCrash>0)
538        solver->crash(1000.0,2);
539      else if (doCrash==0)
540        solver->crash(1000.0,0);
541      if (algorithm<0)
542        doPrimal=false;
543      else if (algorithm>0)
544        doPrimal=true;
545      disasterHandler_->setOsiModel(this);
546      disasterHandler_->setSimplex(solver); // as "borrowed"
547      bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
548      if (!doPrimal) 
549        disasterHandler_->setWhereFrom(4);
550      else
551        disasterHandler_->setWhereFrom(6);
552      if (inCbcOrOther)
553        solver->setDisasterHandler(disasterHandler_);
554      if (!doPrimal) {
555        //printf("doing dual\n");
556        solver->dual(0);
557        totalIterations += solver->numberIterations();
558        if (inCbcOrOther) {
559          if(disasterHandler_->inTrouble()) {
560#ifdef COIN_DEVELOP
561            printf("dual trouble a\n");
562#endif
563            if (disasterHandler_->typeOfDisaster()) {
564              // We want to abort
565              abortSearch=true;
566              goto disaster;
567            }
568            // try just going back in
569            disasterHandler_->setPhase(1);
570            solver->dual();
571            totalIterations += solver->numberIterations();
572            if (disasterHandler_->inTrouble()) {
573#ifdef COIN_DEVELOP
574              printf("dual trouble b\n");
575#endif
576              if (disasterHandler_->typeOfDisaster()) {
577                // We want to abort
578                abortSearch=true;
579                goto disaster;
580              }
581              // try primal with original basis
582              disasterHandler_->setPhase(2);
583              setBasis(basis_,solver);
584              solver->primal();
585              totalIterations += solver->numberIterations();
586            }
587            if(disasterHandler_->inTrouble()) {
588#ifdef COIN_DEVELOP
589              printf("disaster - treat as infeasible\n");
590#endif
591              if (disasterHandler_->typeOfDisaster()) {
592                // We want to abort
593                abortSearch=true;
594                goto disaster;
595              }
596              solver->setProblemStatus(1);
597            }
598          }
599          // reset
600          solver->setDisasterHandler(NULL);
601        }
602        lastAlgorithm_=2; // dual
603        // check if clp thought it was in a loop
604        if (solver->status()==3&&!solver->hitMaximumIterations()) {
605          // switch algorithm
606          solver->primal(0);
607          totalIterations += solver->numberIterations();
608          lastAlgorithm_=1; // primal
609        }
610      } else {
611        //printf("doing primal\n");
612        solver->primal(1);
613        totalIterations += solver->numberIterations();
614        if (inCbcOrOther) {
615          if(disasterHandler_->inTrouble()) {
616#ifdef COIN_DEVELOP
617            printf("primal trouble a\n");
618#endif
619            if (disasterHandler_->typeOfDisaster()) {
620              // We want to abort
621              abortSearch=true;
622              goto disaster;
623            }
624            // try just going back in (but with dual)
625            disasterHandler_->setPhase(1);
626            solver->dual();
627            totalIterations += solver->numberIterations();
628            if (disasterHandler_->inTrouble()) {
629#ifdef COIN_DEVELOP
630              printf("primal trouble b\n");
631#endif
632              if (disasterHandler_->typeOfDisaster()) {
633                // We want to abort
634                abortSearch=true;
635                goto disaster;
636              }
637              // try primal with original basis
638              disasterHandler_->setPhase(2);
639              setBasis(basis_,solver);
640              solver->dual();
641              totalIterations += solver->numberIterations();
642            }
643            if(disasterHandler_->inTrouble()) {
644#ifdef COIN_DEVELOP
645              printf("disaster - treat as infeasible\n");
646#endif
647              if (disasterHandler_->typeOfDisaster()) {
648                // We want to abort
649                abortSearch=true;
650                goto disaster;
651              }
652              solver->setProblemStatus(1);
653            }
654          }
655          // reset
656          solver->setDisasterHandler(NULL);
657        }
658        lastAlgorithm_=1; // primal
659        // check if clp thought it was in a loop
660        if (solver->status()==3&&!solver->hitMaximumIterations()) {
661          // switch algorithm
662          solver->dual(0);
663          totalIterations += solver->numberIterations();
664          lastAlgorithm_=2; // dual
665        }
666      }
667    }
668    // If scaled feasible but unscaled infeasible take action
669    if (!solver->status()&&cleanupScaling_) {
670      solver->cleanup(cleanupScaling_);
671    }
672    basis_ = getBasis(solver);
673    //basis_.print();
674    const double * rowScale2 = solver->rowScale();
675    solver->setSpecialOptions(saveOptions);
676    if (!rowScale1&&rowScale2) {
677      // need to release memory
678      if (!solver->savedRowScale_) {
679        solver->setRowScale(NULL);
680        solver->setColumnScale(NULL);
681      } else {
682        solver->rowScale_=NULL;
683        solver->columnScale_=NULL;
684      }
685    }
686  } else {
687    // User doing nothing and all slack basis
688    ClpSolve options=solveOptions_;
689    bool yesNo;
690    OsiHintStrength strength;
691    getHintParam(OsiDoInBranchAndCut,yesNo,strength);
692    if (yesNo) {
693      solver->setSpecialOptions(solver->specialOptions()|1024);
694    }
695    solver->initialSolve(options);
696    totalIterations += solver->numberIterations();
697    lastAlgorithm_ = 2; // say dual
698    // If scaled feasible but unscaled infeasible take action
699    if (!solver->status()&&cleanupScaling_) {
700      solver->cleanup(cleanupScaling_);
701    }
702    basis_ = getBasis(solver);
703    //basis_.print();
704  }
705  solver->messageHandler()->setLogLevel(saveMessageLevel);
706 disaster:
707  if (deleteSolver) {
708    solver->returnModel(*modelPtr_);
709    delete solver;
710  }
711  if (startFinishOptions) {
712    int save = modelPtr_->logLevel();
713    if (save<2) modelPtr_->setLogLevel(0);
714    modelPtr_->dual(0,startFinishOptions);
715    totalIterations += modelPtr_->numberIterations();
716    modelPtr_->setLogLevel(save);
717  }
718  if (saveSolveType==2) {
719    enableSimplexInterface(doingPrimal);
720  }
721  if (savedObjective) {
722    // fix up
723    modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
724    //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
725    modelPtr_->objective_=savedObjective;
726    if (!modelPtr_->problemStatus_) {
727      CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
728      CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
729      if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
730        CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
731      modelPtr_->computeObjectiveValue();
732    }
733  }
734  modelPtr_->setNumberIterations(totalIterations);
735  handler_->setLogLevel(saveMessageLevel2);
736  if (modelPtr_->problemStatus_==3&&lastAlgorithm_==2)
737    modelPtr_->computeObjectiveValue();
738  // mark so we can pick up objective value quickly
739  modelPtr_->upperIn_=0.0;
740  time1 = CoinCpuTime()-time1;
741  totalTime += time1;
742  assert (!modelPtr_->disasterHandler());
743  if (lastAlgorithm_<1||lastAlgorithm_>2)
744    lastAlgorithm_=1;
745  if (abortSearch) {
746    lastAlgorithm_=-911;
747    modelPtr_->setProblemStatus(4);
748  }
749  modelPtr_->whatsChanged_ |= 0x30000;
750  //std::cout<<time1<<" seconds - total "<<totalTime<<std::endl;
751}
752//-----------------------------------------------------------------------------
753void OsiClpSolverInterface::resolve()
754{
755#ifdef COIN_DEVELOP
756  {
757    int i;
758    int n = getNumCols();
759    const double *lower = getColLower() ;
760    const double *upper = getColUpper() ;
761    for (i=0;i<n;i++) {
762      assert (lower[i]<1.0e12);
763      assert (upper[i]>-1.0e12);
764    }
765    n = getNumRows();
766    lower = getRowLower() ;
767    upper = getRowUpper() ;
768    for (i=0;i<n;i++) {
769      assert (lower[i]<1.0e12);
770      assert (upper[i]>-1.0e12);
771    }
772  }
773#endif
774  if ((stuff_.solverOptions_&65536)!=0) {
775    modelPtr_->fastDual2(&stuff_);
776    return;
777  }
778  //void pclp(char *);
779  //pclp("res");
780  bool takeHint;
781  OsiHintStrength strength;
782  bool gotHint = (getHintParam(OsiDoInBranchAndCut,takeHint,strength));
783  assert (gotHint);
784  // mark so we can pick up objective value quickly
785  modelPtr_->upperIn_=0.0;
786  if ((specialOptions_&4096)!=0) {
787    // Quick check to see if optimal
788    modelPtr_->checkSolutionInternal();
789    if (modelPtr_->problemStatus()==0) {
790      modelPtr_->setNumberIterations(0);
791      return;
792    }
793  }
794  int totalIterations=0;
795  bool abortSearch=false;
796  ClpObjective * savedObjective=NULL;
797  double savedDualLimit=modelPtr_->dblParam_[ClpDualObjectiveLimit];
798  if (fakeObjective_) {
799    modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~128));
800    // See if all with costs fixed
801    int numberColumns = modelPtr_->numberColumns_;
802    const double * obj = modelPtr_->objective();
803    const double * lower = modelPtr_->columnLower();
804    const double * upper = modelPtr_->columnUpper();
805    int i;
806    for (i=0;i<numberColumns;i++) {
807      double objValue = obj[i];
808      if (objValue) {
809        if (lower[i]!=upper[i])
810          break;
811      }
812    }
813    if (i==numberColumns) {
814      if ((specialOptions_&524288)==0) {
815        // Set fake
816        savedObjective=modelPtr_->objective_;
817        modelPtr_->objective_=fakeObjective_;
818        modelPtr_->dblParam_[ClpDualObjectiveLimit]=COIN_DBL_MAX;
819      } else {
820        modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()|128);
821      }
822    }
823  }
824  // If using Clp initialSolve and primal - just do here
825  gotHint = (getHintParam(OsiDoDualInResolve,takeHint,strength));
826  assert (gotHint);
827  if (strength!=OsiHintIgnore&&!takeHint&&solveOptions_.getSpecialOption(6)) {
828    ClpSolve options=solveOptions_;
829    // presolve
830    getHintParam(OsiDoPresolveInResolve,takeHint,strength);
831    if (strength!=OsiHintIgnore&&!takeHint)
832      options.setPresolveType(ClpSolve::presolveOff);
833    int saveOptions = modelPtr_->specialOptions();
834    getHintParam(OsiDoInBranchAndCut,takeHint,strength);
835    if (takeHint) {
836      modelPtr_->setSpecialOptions(modelPtr_->specialOptions()|1024);
837    }
838    setBasis(basis_,modelPtr_);
839    modelPtr_->initialSolve(options);
840    lastAlgorithm_ = 1; // say primal
841    // If scaled feasible but unscaled infeasible take action
842    if (!modelPtr_->status()&&cleanupScaling_) {
843      modelPtr_->cleanup(cleanupScaling_);
844    }
845    modelPtr_->setSpecialOptions(saveOptions); // restore
846    basis_ = getBasis(modelPtr_);
847  }
848  int saveSolveType=modelPtr_->solveType();
849  bool doingPrimal = modelPtr_->algorithm()>0;
850  if (saveSolveType==2) {
851    disableSimplexInterface();
852  }
853  int saveOptions = modelPtr_->specialOptions();
854  int startFinishOptions=0;
855  if (specialOptions_!=0x80000000) {
856    if((specialOptions_&1)==0) {
857      startFinishOptions=0;
858      modelPtr_->setSpecialOptions(saveOptions|(64|1024|32768));
859    } else {
860      startFinishOptions=1+4;
861      if ((specialOptions_&8)!=0)
862        startFinishOptions +=2; // allow re-use of factorization
863      if((specialOptions_&4)==0||!takeHint) 
864        modelPtr_->setSpecialOptions(saveOptions|(64|128|512|1024|4096|32768));
865      else
866        modelPtr_->setSpecialOptions(saveOptions|(64|128|512|1024|2048|4096|32768));
867    }
868  } else {
869    modelPtr_->setSpecialOptions(saveOptions|64|32768);
870  }
871  //printf("options %d size %d\n",modelPtr_->specialOptions(),modelPtr_->numberColumns());
872  //modelPtr_->setSolveType(1);
873  // Set message handler to have same levels etc
874  int saveMessageLevel=modelPtr_->logLevel();
875  int messageLevel=messageHandler()->logLevel();
876  bool oldDefault;
877  CoinMessageHandler * saveHandler = modelPtr_->pushMessageHandler(handler_,oldDefault);
878  //printf("basis before dual\n");
879  //basis_.print();
880  setBasis(basis_,modelPtr_);
881#ifdef SAVE_MODEL
882  resolveTry++;
883#if SAVE_MODEL > 1
884  if (resolveTry>=loResolveTry&&
885      resolveTry<=hiResolveTry) {
886    char fileName[20];
887    sprintf(fileName,"save%d.mod",resolveTry);
888    modelPtr_->saveModel(fileName);
889  }
890#endif
891#endif
892  // set reasonable defaults
893  // Switch off printing if asked to
894  gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength));
895  assert (gotHint);
896  if (strength!=OsiHintIgnore&&takeHint) {
897    if (messageLevel>0)
898      messageLevel--;
899  }
900  if (messageLevel<modelPtr_->messageHandler()->logLevel())
901    modelPtr_->messageHandler()->setLogLevel(messageLevel);
902  // See if user set factorization frequency
903  int userFactorizationFrequency = modelPtr_->factorization()->maximumPivots();
904  // scaling
905  if (modelPtr_->solveType()==1) {
906    gotHint = (getHintParam(OsiDoScale,takeHint,strength));
907    assert (gotHint);
908    if (strength==OsiHintIgnore||takeHint) {
909      if (!modelPtr_->scalingFlag())
910        modelPtr_->scaling(3);
911    } else {
912      modelPtr_->scaling(0);
913    }
914  } else {
915    modelPtr_->scaling(0);
916  }
917  // sort out hints;
918  // algorithm -1 force dual, +1 force primal
919  int algorithm = -1;
920  gotHint = (getHintParam(OsiDoDualInResolve,takeHint,strength));
921  assert (gotHint);
922  if (strength!=OsiHintIgnore)
923    algorithm = takeHint ? -1 : 1;
924  //modelPtr_->saveModel("save.bad");
925  // presolve
926  gotHint = (getHintParam(OsiDoPresolveInResolve,takeHint,strength));
927  assert (gotHint);
928  if (strength!=OsiHintIgnore&&takeHint) {
929#ifdef KEEP_SMALL
930    if (smallModel_) {
931      delete [] spareArrays_;
932      spareArrays_ = NULL;
933      delete smallModel_;
934      smallModel_=NULL;
935    }
936#endif
937    ClpPresolve pinfo;
938    if ((specialOptions_&128)!=0) {
939      specialOptions_ &= ~128;
940    }
941    if ((modelPtr_->specialOptions()&1024)!=0) {
942      pinfo.setDoDual(false);
943      pinfo.setDoTripleton(false);
944      pinfo.setDoDupcol(false);
945      pinfo.setDoDuprow(false);
946      pinfo.setDoSingletonColumn(false);
947    }
948    ClpSimplex * model2 = pinfo.presolvedModel(*modelPtr_,1.0e-8);
949    if (!model2) {
950      // problem found to be infeasible - whats best?
951      model2 = modelPtr_;
952    }
953    // return number of rows
954    int * stats = reinterpret_cast<int *> (getApplicationData());
955    if (stats) {
956      stats[0]=model2->numberRows();
957      stats[1]=model2->numberColumns();
958    }
959    //printf("rows %d -> %d, columns %d -> %d\n",
960    //     modelPtr_->numberRows(),model2->numberRows(),
961    //     modelPtr_->numberColumns(),model2->numberColumns());
962    // change from 200
963    if (modelPtr_->factorization()->maximumPivots()==200)
964      model2->factorization()->maximumPivots(100+model2->numberRows()/50);
965    else
966      model2->factorization()->maximumPivots(userFactorizationFrequency);
967    if (algorithm<0) {
968      model2->dual();
969      totalIterations += model2->numberIterations();
970      // check if clp thought it was in a loop
971      if (model2->status()==3&&!model2->hitMaximumIterations()) {
972        // switch algorithm
973        model2->primal();
974        totalIterations += model2->numberIterations();
975      }
976    } else {
977      model2->primal(1);
978      totalIterations += model2->numberIterations();
979      // check if clp thought it was in a loop
980      if (model2->status()==3&&!model2->hitMaximumIterations()) {
981        // switch algorithm
982        model2->dual();
983        totalIterations += model2->numberIterations();
984      }
985    }
986    if (model2!=modelPtr_) {
987      int finalStatus=model2->status();
988      pinfo.postsolve(true);
989   
990      delete model2;
991      // later try without (1) and check duals before solve
992      if (finalStatus!=3&&(finalStatus||modelPtr_->status()==-1)) {
993        modelPtr_->primal(1);
994        totalIterations += modelPtr_->numberIterations();
995        lastAlgorithm_=1; // primal
996        //if (modelPtr_->numberIterations())
997        //printf("****** iterated %d\n",modelPtr_->numberIterations());
998      }
999    }
1000  } else {
1001    //modelPtr_->setLogLevel(63);
1002    //modelPtr_->setDualTolerance(1.0e-7);
1003    if (false&&modelPtr_->scalingFlag_>0&&!modelPtr_->rowScale_&&
1004        !modelPtr_->rowCopy_&&matrixByRow_) {
1005      assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
1006      modelPtr_->setNewRowCopy(new ClpPackedMatrix(*matrixByRow_));
1007    }
1008    if (algorithm<0) {
1009      //writeMps("try1");
1010      int savePerturbation = modelPtr_->perturbation();
1011      if ((specialOptions_&2)!=0)
1012        modelPtr_->setPerturbation(100);
1013      //modelPtr_->messageHandler()->setLogLevel(1);
1014      //writeMpsNative("bad",NULL,NULL,2,1,1.0);
1015      disasterHandler_->setOsiModel(this);
1016      bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
1017#if 0
1018      // See how many integers fixed
1019      bool skipCrunch=true;
1020      const char * integerInformation = modelPtr_->integerType_;
1021      if (integerInformation) {
1022        int numberColumns = modelPtr_->numberColumns_;
1023        const double * lower = modelPtr_->columnLower();
1024        const double * upper = modelPtr_->columnUpper();
1025        int target=CoinMax(1,numberColumns/10000);
1026        for (int i=0;i<numberColumns;i++) {
1027          if (integerInformation[i]) {
1028            if (lower[i]==upper[i]) {
1029              target--;
1030              if (!target) {
1031                skipCrunch=false;
1032                break;
1033              }
1034            }
1035          }
1036        }
1037      }
1038#endif
1039      if((specialOptions_&1)==0||(specialOptions_&2048)!=0/*||skipCrunch*/) {
1040        disasterHandler_->setWhereFrom(0); // dual
1041        if (inCbcOrOther)
1042          modelPtr_->setDisasterHandler(disasterHandler_);
1043        bool specialScale;
1044        if ((specialOptions_&131072)!=0&&!modelPtr_->rowScale_) {
1045          modelPtr_->rowScale_ = rowScale_.array();
1046          modelPtr_->columnScale_ = columnScale_.array();
1047          specialScale=true;
1048        } else {
1049          specialScale=false;
1050        }
1051#ifdef KEEP_SMALL
1052        if (smallModel_) {
1053          delete [] spareArrays_;
1054          spareArrays_ = NULL;
1055          delete smallModel_;
1056          smallModel_=NULL;
1057        }
1058#endif
1059#ifdef CBC_STATISTICS
1060        osi_dual++;
1061#endif
1062        modelPtr_->dual(0,startFinishOptions);
1063        totalIterations += modelPtr_->numberIterations();
1064        if (specialScale) {
1065          modelPtr_->rowScale_ = NULL;
1066          modelPtr_->columnScale_ = NULL;
1067        }
1068      } else {
1069#ifdef CBC_STATISTICS
1070        osi_crunch++;
1071#endif
1072        crunch();
1073        totalIterations += modelPtr_->numberIterations();
1074        if (modelPtr_->problemStatus()==4)
1075          goto disaster;
1076        // should have already been fixed if problems
1077        inCbcOrOther=false;
1078      }
1079      if (inCbcOrOther) {
1080        if(disasterHandler_->inTrouble()) {
1081          if (disasterHandler_->typeOfDisaster()) {
1082            // We want to abort
1083            abortSearch=true;
1084            goto disaster;
1085          }
1086          // try just going back in
1087          disasterHandler_->setPhase(1);
1088          modelPtr_->dual();
1089          totalIterations += modelPtr_->numberIterations();
1090          if (disasterHandler_->inTrouble()) {
1091            if (disasterHandler_->typeOfDisaster()) {
1092              // We want to abort
1093              abortSearch=true;
1094              goto disaster;
1095            }
1096            // try primal with original basis
1097            disasterHandler_->setPhase(2);
1098            setBasis(basis_,modelPtr_);
1099            modelPtr_->primal();
1100            totalIterations += modelPtr_->numberIterations();
1101          }
1102          if(disasterHandler_->inTrouble()) {
1103#ifdef COIN_DEVELOP
1104            printf("disaster - treat as infeasible\n");
1105#endif
1106            if (disasterHandler_->typeOfDisaster()) {
1107              // We want to abort
1108              abortSearch=true;
1109              goto disaster;
1110            }
1111            modelPtr_->setProblemStatus(1);
1112          }
1113        }
1114        // reset
1115        modelPtr_->setDisasterHandler(NULL);
1116      }
1117      if (modelPtr_->problemStatus()==4) {
1118        // bad bounds?
1119        modelPtr_->setProblemStatus(1);
1120      }
1121      if (!modelPtr_->problemStatus()&&0) {
1122        int numberColumns = modelPtr_->numberColumns();
1123        const double * columnLower = modelPtr_->columnLower();
1124        const double * columnUpper = modelPtr_->columnUpper();
1125        int nBad=0;
1126        for (int i=0;i<numberColumns;i++) {
1127          if (columnLower[i]==columnUpper[i]&&modelPtr_->getColumnStatus(i)==ClpSimplex::basic) {
1128            nBad++;
1129            modelPtr_->setColumnStatus(i,ClpSimplex::isFixed);
1130          }
1131        }
1132        if (nBad) {
1133          modelPtr_->primal(1);
1134    totalIterations += modelPtr_->numberIterations();
1135          printf("%d fixed basic - %d iterations\n",nBad,modelPtr_->numberIterations());
1136        }
1137      }
1138      assert (modelPtr_->objectiveValue()<1.0e100);
1139      modelPtr_->setPerturbation(savePerturbation);
1140      lastAlgorithm_=2; // dual
1141      // check if clp thought it was in a loop
1142      if (modelPtr_->status()==3&&!modelPtr_->hitMaximumIterations()) {
1143        modelPtr_->setSpecialOptions(saveOptions);
1144        // switch algorithm
1145        //modelPtr_->messageHandler()->setLogLevel(63);
1146        // Allow for catastrophe
1147        int saveMax = modelPtr_->maximumIterations();
1148        int numberIterations = modelPtr_->numberIterations();
1149        int numberRows = modelPtr_->numberRows();
1150        int numberColumns = modelPtr_->numberColumns();
1151        if (modelPtr_->maximumIterations()>100000+numberIterations)
1152          modelPtr_->setMaximumIterations(numberIterations + 1000 + 2*numberRows+numberColumns);
1153        modelPtr_->primal(0,startFinishOptions);
1154        totalIterations += modelPtr_->numberIterations();
1155        modelPtr_->setMaximumIterations(saveMax);
1156        lastAlgorithm_=1; // primal
1157        if (modelPtr_->status()==3&&!modelPtr_->hitMaximumIterations()) {
1158#ifdef COIN_DEVELOP
1159          printf("in trouble - try all slack\n");
1160#endif
1161          CoinWarmStartBasis allSlack;
1162          setBasis(allSlack,modelPtr_);
1163          modelPtr_->dual();
1164          totalIterations += modelPtr_->numberIterations();
1165          if (modelPtr_->status()==3&&!modelPtr_->hitMaximumIterations()) {
1166            if (modelPtr_->numberPrimalInfeasibilities()) {
1167#ifdef COIN_DEVELOP
1168              printf("Real real trouble - treat as infeasible\n");
1169#endif
1170              modelPtr_->setProblemStatus(1);
1171            } else {
1172#ifdef COIN_DEVELOP
1173              printf("Real real trouble - treat as optimal\n");
1174#endif
1175              modelPtr_->setProblemStatus(0);
1176            }
1177          }
1178        }
1179      }
1180      assert (modelPtr_->objectiveValue()<1.0e100);
1181    } else {
1182#ifdef KEEP_SMALL
1183      if (smallModel_) {
1184        delete [] spareArrays_;
1185        spareArrays_ = NULL;
1186        delete smallModel_;
1187        smallModel_=NULL;
1188      }
1189#endif
1190      //printf("doing primal\n");
1191#ifdef CBC_STATISTICS
1192      osi_primal++;
1193#endif
1194      modelPtr_->primal(1,startFinishOptions);
1195      totalIterations += modelPtr_->numberIterations();
1196      lastAlgorithm_=1; // primal
1197      // check if clp thought it was in a loop
1198      if (modelPtr_->status()==3&&!modelPtr_->hitMaximumIterations()) {
1199        // switch algorithm
1200        modelPtr_->dual();
1201        totalIterations += modelPtr_->numberIterations();
1202        lastAlgorithm_=2; // dual
1203      }
1204    }
1205  }
1206  // If scaled feasible but unscaled infeasible take action
1207  //if (!modelPtr_->status()&&cleanupScaling_) {
1208  if (cleanupScaling_) {
1209    modelPtr_->cleanup(cleanupScaling_);
1210  }
1211  basis_ = getBasis(modelPtr_);
1212 disaster:
1213  //printf("basis after dual\n");
1214  //basis_.print();
1215  modelPtr_->popMessageHandler(saveHandler,oldDefault);
1216  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
1217  if (saveSolveType==2) {
1218    int saveStatus = modelPtr_->problemStatus_;
1219    enableSimplexInterface(doingPrimal);
1220    modelPtr_->problemStatus_=saveStatus;
1221  }
1222#ifdef COIN_DEVELOP_x
1223  extern bool doingDoneBranch;
1224  if (doingDoneBranch) {
1225    if (modelPtr_->numberIterations())
1226      printf("***** done %d iterations after general\n",modelPtr_->numberIterations());
1227    doingDoneBranch=false;
1228  }
1229#endif
1230  modelPtr_->setNumberIterations(totalIterations);
1231  //modelPtr_->setSolveType(saveSolveType);
1232  if (abortSearch) {
1233    lastAlgorithm_=-911;
1234    modelPtr_->setProblemStatus(4);
1235  }
1236  if (savedObjective) {
1237    // fix up
1238    modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
1239    //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
1240    modelPtr_->objective_=savedObjective;
1241    if (!modelPtr_->problemStatus_) {
1242      CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
1243      CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
1244      if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
1245        CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
1246      modelPtr_->computeObjectiveValue();
1247    }
1248  }
1249  modelPtr_->setSpecialOptions(saveOptions); // restore
1250  if (modelPtr_->problemStatus_==3&&lastAlgorithm_==2)
1251    modelPtr_->computeObjectiveValue();
1252  if (lastAlgorithm_<1||lastAlgorithm_>2)
1253    lastAlgorithm_=1;
1254#ifdef SAVE_MODEL
1255  if (resolveTry>=loResolveTry&&
1256      resolveTry<=hiResolveTry) {
1257    printf("resolve %d took %d iterations - algorithm %d\n",resolveTry,modelPtr_->numberIterations(),lastAlgorithm_);
1258  }
1259#endif
1260  // Make sure whatsChanged not out of sync
1261  if (!modelPtr_->columnUpperWork_)
1262    modelPtr_->whatsChanged_ &= ~0xffff;
1263  modelPtr_->whatsChanged_ |= 0x30000;
1264}
1265// Sort of lexicographic resolve
1266void 
1267OsiClpSolverInterface::lexSolve()
1268{
1269#if 1
1270  abort();
1271#else
1272  ((ClpSimplexPrimal *) modelPtr_)->lexSolve();
1273  printf("itA %d\n",modelPtr_->numberIterations());
1274  modelPtr_->primal();
1275  printf("itB %d\n",modelPtr_->numberIterations());
1276  basis_ = getBasis(modelPtr_);
1277#endif
1278}
1279
1280/* Sets up solver for repeated use by Osi interface.
1281   The normal usage does things like keeping factorization around so can be used.
1282   Will also do things like keep scaling and row copy of matrix if
1283   matrix does not change.
1284   adventure:
1285   0 - safe stuff as above
1286   1 - will take more risks - if it does not work then bug which will be fixed
1287   2 - don't bother doing most extreme termination checks e.g. don't bother
1288       re-factorizing if less than 20 iterations.
1289   3 - Actually safer than 1 (mainly just keeps factorization)
1290
1291   printOut - -1 always skip round common messages instead of doing some work
1292               0 skip if normal defaults
1293               1 leaves
1294  */
1295void 
1296OsiClpSolverInterface::setupForRepeatedUse(int senseOfAdventure, int printOut)
1297{
1298  // First try
1299  switch (senseOfAdventure) {
1300  case 0:
1301    specialOptions_=8;
1302    break;
1303  case 1:
1304    specialOptions_=1+2+8;
1305    break;
1306  case 2:
1307    specialOptions_=1+2+4+8;
1308    break;
1309  case 3:
1310    specialOptions_=1+8;
1311    break;
1312  }
1313  bool stopPrinting=false;
1314  if (printOut<0) {
1315    stopPrinting=true;
1316  } else if (!printOut) {
1317    bool takeHint;
1318    OsiHintStrength strength;
1319    getHintParam(OsiDoReducePrint,takeHint,strength);
1320    int messageLevel=messageHandler()->logLevel();
1321    if (strength!=OsiHintIgnore&&takeHint) 
1322      messageLevel--;
1323    stopPrinting = (messageLevel<=0);
1324  }
1325#ifndef COIN_DEVELOP
1326  if (stopPrinting) {
1327    CoinMessages * messagesPointer = modelPtr_->messagesPointer();
1328    // won't even build messages
1329    messagesPointer->setDetailMessages(100,10000,reinterpret_cast<int *> (NULL));
1330  }
1331#endif
1332}
1333#ifndef NDEBUG
1334// For errors to make sure print to screen
1335// only called in debug mode
1336static void indexError(int index,
1337                        std::string methodName)
1338{
1339  std::cerr<<"Illegal index "<<index<<" in OsiClpSolverInterface::"<<methodName<<std::endl;
1340  throw CoinError("Illegal index",methodName,"OsiClpSolverInterface");
1341}
1342#endif
1343//#############################################################################
1344// Parameter related methods
1345//#############################################################################
1346
1347bool
1348OsiClpSolverInterface::setIntParam(OsiIntParam key, int value)
1349{
1350  return modelPtr_->setIntParam(static_cast<ClpIntParam> (key), value);
1351}
1352
1353//-----------------------------------------------------------------------------
1354
1355bool
1356OsiClpSolverInterface::setDblParam(OsiDblParam key, double value)
1357{
1358  if (key != OsiLastDblParam ) {
1359    if (key==OsiDualObjectiveLimit||key==OsiPrimalObjectiveLimit)
1360      value *= modelPtr_->optimizationDirection();
1361    return modelPtr_->setDblParam(static_cast<ClpDblParam> (key), value);
1362  } else {
1363    return false;
1364  }
1365}
1366
1367//-----------------------------------------------------------------------------
1368
1369bool
1370OsiClpSolverInterface::setStrParam(OsiStrParam key, const std::string & value)
1371{
1372  assert (key!=OsiSolverName);
1373  if (key != OsiLastStrParam ) {
1374    return modelPtr_->setStrParam(static_cast<ClpStrParam> (key), value);
1375  } else {
1376    return false;
1377  }
1378}
1379
1380
1381//-----------------------------------------------------------------------------
1382
1383bool
1384OsiClpSolverInterface::getIntParam(OsiIntParam key, int& value) const 
1385{
1386  return modelPtr_->getIntParam(static_cast<ClpIntParam> (key), value);
1387}
1388
1389//-----------------------------------------------------------------------------
1390
1391bool
1392OsiClpSolverInterface::getDblParam(OsiDblParam key, double& value) const
1393{
1394  if (key != OsiLastDblParam ) {
1395    bool condition =  modelPtr_->getDblParam(static_cast<ClpDblParam> (key), value);
1396    if (key==OsiDualObjectiveLimit||key==OsiPrimalObjectiveLimit)
1397      value *= modelPtr_->optimizationDirection();
1398    return condition;
1399  } else {
1400    return false;
1401  }
1402}
1403
1404//-----------------------------------------------------------------------------
1405
1406bool
1407OsiClpSolverInterface::getStrParam(OsiStrParam key, std::string & value) const
1408{
1409  if ( key==OsiSolverName ) {
1410    value = "clp";
1411    return true;
1412  }
1413  if (key != OsiLastStrParam ) {
1414    return modelPtr_->getStrParam(static_cast<ClpStrParam> (key), value);
1415  } else {
1416    return false;
1417  }
1418}
1419
1420
1421//#############################################################################
1422// Methods returning info on how the solution process terminated
1423//#############################################################################
1424
1425bool OsiClpSolverInterface::isAbandoned() const
1426{
1427  // not sure about -1 (should not happen)
1428  return (modelPtr_->status()==4||modelPtr_->status()==-1||
1429          (modelPtr_->status()==1&&modelPtr_->secondaryStatus()==8));
1430}
1431
1432bool OsiClpSolverInterface::isProvenOptimal() const
1433{
1434
1435  const int stat = modelPtr_->status();
1436  return (stat == 0);
1437}
1438
1439bool OsiClpSolverInterface::isProvenPrimalInfeasible() const
1440{
1441
1442  const int stat = modelPtr_->status();
1443  if (stat != 1)
1444     return false;
1445  return true;
1446}
1447
1448bool OsiClpSolverInterface::isProvenDualInfeasible() const
1449{
1450  const int stat = modelPtr_->status();
1451  return stat == 2;
1452}
1453/*
1454   NOTE - Coding if limit > 1.0e30 says that 1.0e29 is loose bound
1455   so all maximization tests are changed
1456*/
1457bool OsiClpSolverInterface::isPrimalObjectiveLimitReached() const
1458{
1459  double limit = 0.0;
1460  modelPtr_->getDblParam(ClpPrimalObjectiveLimit, limit);
1461  if (fabs(limit) > 1e30) {
1462    // was not ever set
1463    return false;
1464  }
1465   
1466  const double obj = modelPtr_->objectiveValue();
1467  int maxmin = static_cast<int> (modelPtr_->optimizationDirection());
1468
1469  switch (lastAlgorithm_) {
1470   case 0: // no simplex was needed
1471     return maxmin > 0 ? (obj < limit) /*minim*/ : (-obj < limit) /*maxim*/;
1472   case 2: // dual simplex
1473     if (modelPtr_->status() == 0) // optimal
1474        return maxmin > 0 ? (obj < limit) /*minim*/ : (-obj < limit) /*maxim*/;
1475     return false;
1476   case 1: // primal simplex
1477     return maxmin > 0 ? (obj < limit) /*minim*/ : (-obj < limit) /*maxim*/;
1478  }
1479  return false; // fake return
1480}
1481
1482bool OsiClpSolverInterface::isDualObjectiveLimitReached() const
1483{
1484
1485  const int stat = modelPtr_->status();
1486  if (stat == 1)
1487  return true;
1488  double limit = 0.0;
1489  modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
1490  if (fabs(limit) > 1e30) {
1491    // was not ever set
1492    return false;
1493  }
1494   
1495  const double obj = modelPtr_->objectiveValue();
1496  int maxmin = static_cast<int> (modelPtr_->optimizationDirection());
1497
1498  switch (lastAlgorithm_) {
1499   case 0: // no simplex was needed
1500     return maxmin > 0 ? (obj > limit) /*minim*/ : (-obj > limit) /*maxim*/;
1501   case 1: // primal simplex
1502     if (stat == 0) // optimal
1503        return maxmin > 0 ? (obj > limit) /*minim*/ : (-obj > limit) /*maxim*/;
1504     return false;
1505   case 2: // dual simplex
1506     if (stat != 0 && stat != 3)
1507        // over dual limit
1508        return true;
1509     return maxmin > 0 ? (obj > limit) /*minim*/ : (-obj > limit) /*maxim*/;
1510  }
1511  return false; // fake return
1512}
1513
1514bool OsiClpSolverInterface::isIterationLimitReached() const
1515{
1516  const int stat = modelPtr_->status();
1517  return (stat == 3);
1518}
1519
1520//#############################################################################
1521// WarmStart related methods
1522//#############################################################################
1523CoinWarmStart *OsiClpSolverInterface::getEmptyWarmStart () const
1524  { return (static_cast<CoinWarmStart *>(new CoinWarmStartBasis())) ; }
1525
1526CoinWarmStart* OsiClpSolverInterface::getWarmStart() const
1527{
1528
1529  return new CoinWarmStartBasis(basis_);
1530}
1531/* Get warm start information.
1532   Return warm start information for the current state of the solver
1533   interface. If there is no valid warm start information, an empty warm
1534   start object wil be returned.  This does not necessarily create an
1535   object - may just point to one.  must Delete set true if user
1536   should delete returned object.
1537   OsiClp version always returns pointer and false.
1538*/
1539CoinWarmStart* 
1540OsiClpSolverInterface::getPointerToWarmStart(bool & mustDelete) 
1541{
1542  mustDelete = false;
1543  return &basis_;
1544}
1545
1546//-----------------------------------------------------------------------------
1547
1548bool OsiClpSolverInterface::setWarmStart(const CoinWarmStart* warmstart)
1549{
1550  modelPtr_->whatsChanged_ &= 0xffff;
1551  const CoinWarmStartBasis* ws =
1552    dynamic_cast<const CoinWarmStartBasis*>(warmstart);
1553  if (ws) {
1554    basis_ = CoinWarmStartBasis(*ws);
1555    return true;
1556  } else if (!warmstart) {
1557    // create from current basis
1558    basis_ = getBasis(modelPtr_);
1559    return true;
1560  } else {
1561    return false;
1562  }
1563}
1564
1565//#############################################################################
1566// Hotstart related methods (primarily used in strong branching)
1567//#############################################################################
1568void OsiClpSolverInterface::markHotStart()
1569{
1570#ifdef CBC_STATISTICS
1571  osi_hot++;
1572#endif
1573  //printf("HotStart options %x changed %x, small %x spare %x\n",
1574  // specialOptions_,modelPtr_->whatsChanged_,
1575  // smallModel_,spareArrays_);
1576  modelPtr_->setProblemStatus(0);
1577  saveData_.perturbation_=0;
1578  saveData_.specialOptions_ = modelPtr_->specialOptions_;
1579  modelPtr_->specialOptions_ |= 0x1000000;
1580  modelPtr_->specialOptions_ = saveData_.specialOptions_; 
1581  ClpObjective * savedObjective=NULL;
1582  double savedDualLimit=modelPtr_->dblParam_[ClpDualObjectiveLimit];
1583  if (fakeObjective_) {
1584    modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~128));
1585    // See if all with costs fixed
1586    int numberColumns = modelPtr_->numberColumns_;
1587    const double * obj = modelPtr_->objective();
1588    const double * lower = modelPtr_->columnLower();
1589    const double * upper = modelPtr_->columnUpper();
1590    int i;
1591    for (i=0;i<numberColumns;i++) {
1592      double objValue = obj[i];
1593      if (objValue) {
1594        if (lower[i]!=upper[i])
1595          break;
1596      }
1597    }
1598    if (i==numberColumns) {
1599      if ((specialOptions_&524288)==0) {
1600        // Set fake
1601        savedObjective=modelPtr_->objective_;
1602        modelPtr_->objective_=fakeObjective_;
1603        modelPtr_->dblParam_[ClpDualObjectiveLimit]=COIN_DBL_MAX;
1604        saveData_.perturbation_=1;
1605      } else {
1606        modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()|128);
1607      }
1608    }
1609  }
1610#define CLEAN_HOT_START
1611#ifdef CLEAN_HOT_START
1612  if ((specialOptions_&65536)!=0) {
1613    //specialOptions_ |= 65536;
1614    saveData_.scalingFlag_=modelPtr_->logLevel();
1615    if (modelPtr_->logLevel()<2)
1616      modelPtr_->setLogLevel(0);
1617    assert ((specialOptions_&128)==0);
1618    // space for save arrays
1619    int numberColumns = modelPtr_->numberColumns();
1620    int numberRows = modelPtr_->numberRows();
1621    // Get space for strong branching
1622    int size = static_cast<int>((1+4*(numberRows+numberColumns))*sizeof(double));
1623    // and for save of original column bounds
1624    size += static_cast<int>(2*numberColumns*sizeof(double));
1625    size += static_cast<int>((1+4*numberRows+2*numberColumns)*sizeof(int));
1626    size += numberRows+numberColumns;
1627    assert (spareArrays_==NULL);
1628    spareArrays_ = new char[size];
1629    //memset(spareArrays_,0x20,size);
1630    // Setup for strong branching
1631    assert (factorization_==NULL);
1632    if ((specialOptions_&131072)!=0) {
1633      assert (lastNumberRows_>=0);
1634      if (modelPtr_->rowScale_!=rowScale_.array()) {
1635        assert(modelPtr_->columnScale_!=columnScale_.array());
1636        delete [] modelPtr_->rowScale_;
1637        modelPtr_->rowScale_=NULL;
1638        delete [] modelPtr_->columnScale_;
1639        modelPtr_->columnScale_=NULL;
1640        if (lastNumberRows_==modelPtr_->numberRows()) {
1641          // use scaling
1642          modelPtr_->rowScale_ = rowScale_.array();
1643          modelPtr_->columnScale_ = columnScale_.array();
1644        } else {
1645          specialOptions_ &= ~131072;
1646        }
1647      }
1648      lastNumberRows_ = -1 -lastNumberRows_;
1649    }
1650    factorization_ = static_cast<ClpSimplexDual *>(modelPtr_)->setupForStrongBranching(spareArrays_,numberRows,
1651                                                                                       numberColumns,true);
1652    double * arrayD = reinterpret_cast<double *> (spareArrays_);
1653    arrayD[0]=modelPtr_->objectiveValue()* modelPtr_->optimizationDirection();
1654    double * saveSolution = arrayD+1;
1655    double * saveLower = saveSolution + (numberRows+numberColumns);
1656    double * saveUpper = saveLower + (numberRows+numberColumns);
1657    double * saveObjective = saveUpper + (numberRows+numberColumns);
1658    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
1659    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
1660    CoinMemcpyN( modelPtr_->columnLower(),numberColumns, saveLowerOriginal);
1661    CoinMemcpyN( modelPtr_->columnUpper(),numberColumns, saveUpperOriginal);
1662#if 0
1663    if (whichRange_&&whichRange_[0]) {
1664      // get ranging information
1665      int numberToDo = whichRange_[0];
1666      int * which = new int [numberToDo];
1667      // Convert column numbers
1668      int * backColumn = whichColumn+numberColumns;
1669      for (int i=0;i<numberToDo;i++) {
1670        int iColumn = whichRange_[i+1];
1671        which[i]=backColumn[iColumn];
1672      }
1673      double * downRange=new double [numberToDo];
1674      double * upRange=new double [numberToDo];
1675      int * whichDown = new int [numberToDo];
1676      int * whichUp = new int [numberToDo];
1677      modelPtr_->gutsOfSolution(NULL,NULL,false);
1678      // Tell code we can increase costs in some cases
1679      modelPtr_->setCurrentDualTolerance(0.0);
1680      ((ClpSimplexOther *) modelPtr_)->dualRanging(numberToDo,which,
1681                                                   upRange, whichUp, downRange, whichDown);
1682      delete [] whichDown;
1683      delete [] whichUp;
1684      delete [] which;
1685      rowActivity_=upRange;
1686      columnActivity_=downRange;
1687    }
1688#endif
1689    if (savedObjective) {
1690      // fix up
1691      modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
1692      //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
1693      modelPtr_->objective_=savedObjective;
1694      if (!modelPtr_->problemStatus_) {
1695        CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
1696        CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
1697        if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
1698          CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
1699        modelPtr_->computeObjectiveValue();
1700      }
1701    }
1702    return;
1703  }
1704#endif
1705  if ((specialOptions_&8192)==0&&false) { // ||(specialOptions_&65536)!=0) {
1706    delete ws_;
1707    ws_ = dynamic_cast<CoinWarmStartBasis*>(getWarmStart());
1708    int numberRows = modelPtr_->numberRows();
1709    rowActivity_= new double[numberRows];
1710    CoinMemcpyN(modelPtr_->primalRowSolution(),numberRows,rowActivity_);
1711    int numberColumns = modelPtr_->numberColumns();
1712    columnActivity_= new double[numberColumns];
1713    CoinMemcpyN(modelPtr_->primalColumnSolution(),numberColumns,columnActivity_);
1714  } else {
1715#if 0
1716    int saveLevel = modelPtr_->logLevel();
1717    modelPtr_->setLogLevel(0);
1718    //modelPtr_->dual();
1719    OsiClpSolverInterface::resolve();
1720    if (modelPtr_->numberIterations()>0)
1721      printf("**** iterated large %d\n",modelPtr_->numberIterations());
1722    //else
1723    //printf("no iterations\n");
1724    modelPtr_->setLogLevel(saveLevel);
1725#endif
1726    // called from CbcNode
1727    int numberColumns = modelPtr_->numberColumns();
1728    int numberRows = modelPtr_->numberRows();
1729    // Get space for crunch and strong branching (too much)
1730    int size = static_cast<int>((1+4*(numberRows+numberColumns))*sizeof(double));
1731    // and for save of original column bounds
1732    size += static_cast<int>(2*numberColumns*sizeof(double));
1733    size += static_cast<int>((1+4*numberRows+2*numberColumns)*sizeof(int));
1734    size += numberRows+numberColumns;
1735#ifdef KEEP_SMALL
1736    if(smallModel_&&(modelPtr_->whatsChanged_&0x30000)!=0x30000) {
1737      //printf("Bounds changed ? %x\n",modelPtr_->whatsChanged_);
1738      delete smallModel_;
1739      smallModel_=NULL;
1740    }
1741    if (!smallModel_) {
1742      delete [] spareArrays_;
1743      spareArrays_ = NULL;
1744    }
1745#endif
1746    if (spareArrays_==NULL) {
1747      //if (smallModel_)
1748      //printf("small model %x\n",smallModel_);
1749      delete smallModel_;
1750      smallModel_=NULL;
1751      spareArrays_ = new char[size];
1752      //memset(spareArrays_,0x20,size);
1753    } else {
1754      double * arrayD = reinterpret_cast<double *> (spareArrays_);
1755      double * saveSolution = arrayD+1;
1756      double * saveLower = saveSolution + (numberRows+numberColumns);
1757      double * saveUpper = saveLower + (numberRows+numberColumns);
1758      double * saveObjective = saveUpper + (numberRows+numberColumns);
1759      double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
1760      double * saveUpperOriginal = saveLowerOriginal + numberColumns;
1761      double * lowerOriginal = modelPtr_->columnLower();
1762      double * upperOriginal = modelPtr_->columnUpper();
1763      arrayD = saveUpperOriginal + numberColumns;
1764      int * savePivot = reinterpret_cast<int *> (arrayD);
1765      int * whichRow = savePivot+numberRows;
1766      int * whichColumn = whichRow + 3*numberRows;
1767      int nSame=0;
1768      int nSub=0;
1769      for (int i=0;i<numberColumns;i++) {
1770        double lo = lowerOriginal[i];
1771        //char * xx = (char *) (saveLowerOriginal+i);
1772        //assert (xx[0]!=0x20||xx[1]!=0x20);
1773        //xx = (char *) (saveUpperOriginal+i);
1774        //assert (xx[0]!=0x20||xx[1]!=0x20);
1775        double loOld = saveLowerOriginal[i];
1776        //assert (!loOld||fabs(loOld)>1.0e-30);
1777        double up = upperOriginal[i];
1778        double upOld = saveUpperOriginal[i];
1779        if (lo>=loOld&&up<=upOld) {
1780          if (lo==loOld&&up==upOld) {
1781            nSame++;
1782          } else {
1783            nSub++;
1784            //if (!isInteger(i))
1785            //nSub+=10;
1786          }
1787        }
1788      }
1789      //printf("Mark Hot %d bounds same, %d interior, %d bad\n",
1790      //     nSame,nSub,numberColumns-nSame-nSub);
1791      if (nSame<numberColumns) {
1792        if (nSame+nSub<numberColumns) {
1793          delete smallModel_;
1794          smallModel_=NULL;
1795        } else {
1796          // we can fix up (but should we if large number fixed?)
1797          assert (smallModel_);
1798          double * lowerSmall = smallModel_->columnLower();
1799          double * upperSmall = smallModel_->columnUpper();
1800          int numberColumns2 = smallModel_->numberColumns();
1801          for (int i=0;i<numberColumns2;i++) {
1802            int iColumn = whichColumn[i];
1803            lowerSmall[i]=lowerOriginal[iColumn];
1804            upperSmall[i]=upperOriginal[iColumn];
1805          }
1806        }
1807      }
1808    }
1809    double * arrayD = reinterpret_cast<double *> (spareArrays_);
1810    arrayD[0]=modelPtr_->objectiveValue()* modelPtr_->optimizationDirection();
1811    double * saveSolution = arrayD+1;
1812    double * saveLower = saveSolution + (numberRows+numberColumns);
1813    double * saveUpper = saveLower + (numberRows+numberColumns);
1814    double * saveObjective = saveUpper + (numberRows+numberColumns);
1815    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
1816    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
1817    arrayD = saveUpperOriginal + numberColumns;
1818    int * savePivot = reinterpret_cast<int *> (arrayD);
1819    int * whichRow = savePivot+numberRows;
1820    int * whichColumn = whichRow + 3*numberRows;
1821    int * arrayI = whichColumn + 2*numberColumns;
1822    //unsigned char * saveStatus = (unsigned char *) (arrayI+1);
1823    // Use dual region
1824    double * rhs = modelPtr_->dualRowSolution();
1825    int nBound=0;
1826    ClpSimplex * small;
1827#ifndef KEEP_SMALL
1828    assert (!smallModel_);
1829    small = static_cast<ClpSimplexOther *> (modelPtr_)->crunch(rhs,whichRow,whichColumn,nBound,true);
1830    bool needSolveInSetupHotStart=true;
1831#else
1832    bool needSolveInSetupHotStart=true;
1833    if (!smallModel_) {
1834#ifndef NDEBUG
1835      CoinFillN(whichRow,3*numberRows+2*numberColumns,-1);
1836#endif
1837      small = static_cast<ClpSimplexOther *> (modelPtr_)->crunch(rhs,whichRow,whichColumn,nBound,true);
1838#ifndef NDEBUG
1839      int nCopy = 3*numberRows+2*numberColumns;
1840      for (int i=0;i<nCopy;i++)
1841        assert (whichRow[i]>=-CoinMax(numberRows,numberColumns)&&whichRow[i]<CoinMax(numberRows,numberColumns));
1842#endif
1843      smallModel_=small;
1844      //int hotIts = small->intParam_[ClpMaxNumIterationHotStart];
1845      //if (5*small->factorization_->maximumPivots()>
1846      //  4*hotIts)
1847      //small->factorization_->maximumPivots(hotIts+1);
1848    } else {
1849      assert((modelPtr_->whatsChanged_&0x30000)==0x30000);
1850      //delete [] spareArrays_;
1851      //spareArrays_ = NULL;
1852      assert (spareArrays_);
1853      int nCopy = 3*numberRows+2*numberColumns;
1854      nBound = whichRow[nCopy];
1855#ifndef NDEBUG
1856      for (int i=0;i<nCopy;i++)
1857        assert (whichRow[i]>=-CoinMax(numberRows,numberColumns)&&whichRow[i]<CoinMax(numberRows,numberColumns));
1858#endif
1859      needSolveInSetupHotStart=false;
1860      small = smallModel_;
1861    }
1862#endif
1863    if (small) {
1864      small->specialOptions_ |= 262144;
1865      small->specialOptions_ &= ~65536;
1866    }
1867    if (small&&(specialOptions_&131072)!=0) {
1868      assert (lastNumberRows_>=0);
1869      int numberRows2 = small->numberRows();
1870      int numberColumns2 = small->numberColumns();
1871      double * rowScale2 = new double [2*numberRows2];
1872      const double * rowScale = rowScale_.array();
1873      double * inverseScale2 = rowScale2+numberRows2;
1874      const double * inverseScale = rowScale+modelPtr_->numberRows_;
1875      int i;
1876      for (i=0;i<numberRows2;i++) {
1877        int iRow = whichRow[i];
1878        rowScale2[i]=rowScale[iRow];
1879        inverseScale2[i]=inverseScale[iRow];
1880      }
1881      small->setRowScale(rowScale2);
1882      double * columnScale2 = new double [2*numberColumns2];
1883      const double * columnScale = columnScale_.array();
1884      inverseScale2 = columnScale2+numberColumns2;
1885      inverseScale = columnScale+modelPtr_->numberColumns_;
1886      for (i=0;i<numberColumns2;i++) {
1887        int iColumn = whichColumn[i];
1888        columnScale2[i]=columnScale[iColumn];
1889        inverseScale2[i]=inverseScale[iColumn];
1890      }
1891      small->setColumnScale(columnScale2);
1892    }
1893    if (!small) {
1894      // should never be infeasible .... but
1895      delete [] spareArrays_;
1896      spareArrays_=NULL;
1897      delete ws_;
1898      ws_ = dynamic_cast<CoinWarmStartBasis*>(getWarmStart());
1899      int numberRows = modelPtr_->numberRows();
1900      rowActivity_= new double[numberRows];
1901      CoinMemcpyN(modelPtr_->primalRowSolution(),numberRows,rowActivity_);
1902      int numberColumns = modelPtr_->numberColumns();
1903      columnActivity_= new double[numberColumns];
1904      CoinMemcpyN(modelPtr_->primalColumnSolution(),numberColumns,columnActivity_);
1905      modelPtr_->setProblemStatus(1);
1906      if (savedObjective) {
1907        // fix up
1908        modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
1909        //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
1910        modelPtr_->objective_=savedObjective;
1911      }
1912      return;
1913    }
1914    int clpOptions = modelPtr_->specialOptions();
1915    clpOptions &= ~65536;
1916    if((specialOptions_&1)==0) {
1917      small->setSpecialOptions(clpOptions|(64|1024));
1918    } else {
1919      if((specialOptions_&4)==0) 
1920        small->setSpecialOptions(clpOptions|(64|128|512|1024|4096));
1921      else
1922        small->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096));
1923    }
1924    arrayI[0]=nBound;
1925    assert (smallModel_==NULL||small==smallModel_);
1926    if ((specialOptions_&256)!=0||1) {
1927      // only need to do this on second pass in CbcNode
1928      if (modelPtr_->logLevel()<2) small->setLogLevel(0);
1929      small->specialOptions_ |= 262144;
1930      small->moreSpecialOptions_ = modelPtr_->moreSpecialOptions_;
1931#define SETUP_HOT
1932#ifndef SETUP_HOT
1933      small->dual();
1934#else
1935      assert (factorization_==NULL);
1936      //needSolveInSetupHotStart=true;
1937      ClpFactorization * factorization = static_cast<ClpSimplexDual *>(small)->setupForStrongBranching(spareArrays_,
1938                       numberRows,numberColumns,
1939                                                                                                       needSolveInSetupHotStart);
1940#endif
1941      if (small->numberIterations()>0&&small->logLevel()>2)
1942        printf("**** iterated small %d\n",small->numberIterations());
1943      //small->setLogLevel(0);
1944      // Could be infeasible if forced one way (and other way stopped on iterations)
1945      if (small->status()==1) {
1946#ifndef KEEP_SMALL
1947        if (small!=modelPtr_)
1948          delete small;
1949        //delete smallModel_;
1950        //smallModel_=NULL;
1951        assert (!smallModel_);
1952#else
1953        assert (small==smallModel_);
1954        if (smallModel_!=modelPtr_) {
1955          delete smallModel_;
1956        }
1957        smallModel_=NULL;
1958#endif
1959        delete [] spareArrays_;
1960        spareArrays_=NULL;
1961        delete ws_;
1962        ws_ = dynamic_cast<CoinWarmStartBasis*>(getWarmStart());
1963        int numberRows = modelPtr_->numberRows();
1964        rowActivity_= new double[numberRows];
1965 CoinMemcpyN(modelPtr_->primalRowSolution(),    numberRows,rowActivity_);
1966        int numberColumns = modelPtr_->numberColumns();
1967        columnActivity_= new double[numberColumns];
1968 CoinMemcpyN(modelPtr_->primalColumnSolution(), numberColumns,columnActivity_);
1969        modelPtr_->setProblemStatus(1);
1970        if (savedObjective) {
1971          // fix up
1972          modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
1973          //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
1974          modelPtr_->objective_=savedObjective;
1975        }
1976        return;
1977      } else {
1978        // update model
1979        static_cast<ClpSimplexOther *> (modelPtr_)->afterCrunch(*small,whichRow,whichColumn,nBound);
1980      } 
1981#ifndef SETUP_HOT
1982      assert (factorization_==NULL);
1983      factorization_ = static_cast<ClpSimplexDual *>(small)->setupForStrongBranching(spareArrays_,numberRows,
1984                             numberColumns,false);
1985#else
1986      assert (factorization!=NULL);
1987      factorization_ = factorization;
1988#endif
1989    } else {
1990      assert (factorization_==NULL);
1991      factorization_ = static_cast<ClpSimplexDual *>(small)->setupForStrongBranching(spareArrays_,numberRows,
1992                                                                                     numberColumns,false);
1993    }
1994    smallModel_=small;
1995    if (modelPtr_->logLevel()<2) smallModel_->setLogLevel(0);
1996    // Setup for strong branching
1997    int numberColumns2 = smallModel_->numberColumns();
1998    CoinMemcpyN( modelPtr_->columnLower(),numberColumns, saveLowerOriginal);
1999    CoinMemcpyN( modelPtr_->columnUpper(),numberColumns, saveUpperOriginal);
2000    const double * smallLower = smallModel_->columnLower();
2001    const double * smallUpper = smallModel_->columnUpper();
2002    // But modify if bounds changed in small
2003    for (int i=0;i<numberColumns2;i++) {
2004      int iColumn = whichColumn[i];
2005      saveLowerOriginal[iColumn] = CoinMax(saveLowerOriginal[iColumn],
2006                                           smallLower[i]);
2007      saveUpperOriginal[iColumn] = CoinMin(saveUpperOriginal[iColumn],
2008                                           smallUpper[i]);
2009    }
2010    if (whichRange_&&whichRange_[0]) {
2011      // get ranging information
2012      int numberToDo = whichRange_[0];
2013      int * which = new int [numberToDo];
2014      // Convert column numbers
2015      int * backColumn = whichColumn+numberColumns;
2016      for (int i=0;i<numberToDo;i++) {
2017        int iColumn = whichRange_[i+1];
2018        which[i]=backColumn[iColumn];
2019      }
2020      double * downRange=new double [numberToDo];
2021      double * upRange=new double [numberToDo];
2022      int * whichDown = new int [numberToDo];
2023      int * whichUp = new int [numberToDo];
2024      smallModel_->setFactorization(*factorization_);
2025      smallModel_->gutsOfSolution(NULL,NULL,false);
2026      // Tell code we can increase costs in some cases
2027      smallModel_->setCurrentDualTolerance(0.0);
2028      static_cast<ClpSimplexOther *> (smallModel_)->dualRanging(numberToDo,which,
2029                         upRange, whichUp, downRange, whichDown);
2030      delete [] whichDown;
2031      delete [] whichUp;
2032      delete [] which;
2033      rowActivity_=upRange;
2034      columnActivity_=downRange;
2035    }
2036  }
2037    if (savedObjective) {
2038      // fix up
2039      modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2040      //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2041      modelPtr_->objective_=savedObjective;
2042      if (!modelPtr_->problemStatus_) {
2043        CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
2044        CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
2045        if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
2046          CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
2047        modelPtr_->computeObjectiveValue();
2048      }
2049    }
2050}
2051
2052void OsiClpSolverInterface::solveFromHotStart()
2053{
2054#ifdef KEEP_SMALL
2055  if (!spareArrays_) {
2056    assert (!smallModel_);
2057    assert (modelPtr_->problemStatus_==1);
2058    return;
2059  }
2060#endif
2061  ClpObjective * savedObjective=NULL;
2062  double savedDualLimit=modelPtr_->dblParam_[ClpDualObjectiveLimit];
2063  if (saveData_.perturbation_) {
2064    // Set fake
2065    savedObjective=modelPtr_->objective_;
2066    modelPtr_->objective_=fakeObjective_;
2067    modelPtr_->dblParam_[ClpDualObjectiveLimit]=COIN_DBL_MAX;
2068  }
2069  int numberRows = modelPtr_->numberRows();
2070  int numberColumns = modelPtr_->numberColumns();
2071  modelPtr_->getIntParam(ClpMaxNumIteration,itlimOrig_);
2072  int itlim;
2073  modelPtr_->getIntParam(ClpMaxNumIterationHotStart, itlim);
2074  // Is there an extra copy of scaled bounds
2075  int extraCopy = (modelPtr_->maximumRows_>0) ? modelPtr_->maximumRows_+modelPtr_->maximumColumns_ : 0; 
2076#ifdef CLEAN_HOT_START
2077  if ((specialOptions_&65536)!=0) {
2078    double * arrayD = reinterpret_cast<double *> (spareArrays_);
2079    double saveObjectiveValue = arrayD[0];
2080    double * saveSolution = arrayD+1;
2081    int number = numberRows+numberColumns;
2082    CoinMemcpyN(saveSolution,number,modelPtr_->solutionRegion());
2083    double * saveLower = saveSolution + (numberRows+numberColumns);
2084    CoinMemcpyN(saveLower,number,modelPtr_->lowerRegion());
2085    double * saveUpper = saveLower + (numberRows+numberColumns);
2086    CoinMemcpyN(saveUpper,number,modelPtr_->upperRegion());
2087    double * saveObjective = saveUpper + (numberRows+numberColumns);
2088    CoinMemcpyN(saveObjective,number,modelPtr_->costRegion());
2089    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
2090    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
2091    arrayD = saveUpperOriginal + numberColumns;
2092    int * savePivot = reinterpret_cast<int *> (arrayD);
2093    CoinMemcpyN(savePivot,numberRows,modelPtr_->pivotVariable());
2094    int * whichRow = savePivot+numberRows;
2095    int * whichColumn = whichRow + 3*numberRows;
2096    int * arrayI = whichColumn + 2*numberColumns;
2097    unsigned char * saveStatus = reinterpret_cast<unsigned char *> (arrayI+1);
2098    CoinMemcpyN(saveStatus,number,modelPtr_->statusArray());
2099    modelPtr_->setFactorization(*factorization_);
2100    double * columnLower = modelPtr_->columnLower();
2101    double * columnUpper = modelPtr_->columnUpper();
2102    // make sure whatsChanged_ has 1 set
2103    modelPtr_->setWhatsChanged(511);
2104    double * lowerInternal = modelPtr_->lowerRegion();
2105    double * upperInternal = modelPtr_->upperRegion();
2106    double rhsScale = modelPtr_->rhsScale();
2107    const double * columnScale = NULL;
2108    if (modelPtr_->scalingFlag()>0) 
2109      columnScale = modelPtr_->columnScale() ;
2110    // and do bounds in case dual needs them
2111    int iColumn;
2112    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2113      if (columnLower[iColumn]>saveLowerOriginal[iColumn]) {
2114        double value = columnLower[iColumn];
2115        value *= rhsScale;
2116        if (columnScale)
2117          value /= columnScale[iColumn];
2118        lowerInternal[iColumn]=value;
2119        if (extraCopy)
2120          lowerInternal[iColumn+extraCopy]=value;
2121      }
2122      if (columnUpper[iColumn]<saveUpperOriginal[iColumn]) {
2123        double value = columnUpper[iColumn];
2124        value *= rhsScale;
2125        if (columnScale)
2126          value /= columnScale[iColumn];
2127        upperInternal[iColumn]=value;
2128        if (extraCopy)
2129          upperInternal[iColumn+extraCopy]=value;
2130      }
2131    }
2132    // Start of fast iterations
2133    bool alwaysFinish= ((specialOptions_&32)==0) ? true : false;
2134    //modelPtr_->setLogLevel(1);
2135    modelPtr_->setIntParam(ClpMaxNumIteration,itlim);
2136    int saveNumberFake = (static_cast<ClpSimplexDual *>(modelPtr_))->numberFake_;
2137    int status = (static_cast<ClpSimplexDual *>(modelPtr_))->fastDual(alwaysFinish);
2138    (static_cast<ClpSimplexDual *>(modelPtr_))->numberFake_ = saveNumberFake;
2139   
2140    int problemStatus = modelPtr_->problemStatus();
2141    double objectiveValue =modelPtr_->objectiveValue() * modelPtr_->optimizationDirection();
2142    CoinAssert (modelPtr_->problemStatus()||modelPtr_->objectiveValue()<1.0e50);
2143    // make sure plausible
2144    double obj = CoinMax(objectiveValue,saveObjectiveValue);
2145    if (problemStatus==10||problemStatus<0) {
2146      // was trying to clean up or something odd
2147      if (problemStatus==10)
2148        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2149      status=1;
2150    }
2151    if (status) {
2152      // not finished - might be optimal
2153      modelPtr_->checkPrimalSolution(modelPtr_->solutionRegion(0),
2154                                     modelPtr_->solutionRegion(1));
2155      //modelPtr_->gutsOfSolution(NULL,NULL,0);
2156      //if (problemStatus==3)
2157      //modelPtr_->computeObjectiveValue();
2158      objectiveValue =modelPtr_->objectiveValue() *
2159        modelPtr_->optimizationDirection();
2160      obj = CoinMax(objectiveValue,saveObjectiveValue);
2161      if (!modelPtr_->numberDualInfeasibilities()) { 
2162        double limit = 0.0;
2163        modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
2164        if (modelPtr_->secondaryStatus()==1&&!problemStatus&&obj<limit) {
2165          obj=limit;
2166          problemStatus=3;
2167        }
2168        if (!modelPtr_->numberPrimalInfeasibilities()&&obj<limit) { 
2169          problemStatus=0;
2170        } else if (problemStatus==10) {
2171          problemStatus=3;
2172        } else if (!modelPtr_->numberPrimalInfeasibilities()) {
2173          problemStatus=1; // infeasible
2174        } 
2175      } else {
2176        // can't say much
2177        //if (problemStatus==3)
2178        //modelPtr_->computeObjectiveValue();
2179        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2180        problemStatus=3;
2181      }
2182    } else if (!problemStatus) {
2183      if (modelPtr_->isDualObjectiveLimitReached()) 
2184        problemStatus=1; // infeasible
2185    }
2186    if (status&&!problemStatus) {
2187      problemStatus=3; // can't be sure
2188      lastAlgorithm_=1;
2189    }
2190    if (problemStatus<0)
2191      problemStatus=3;
2192    modelPtr_->setProblemStatus(problemStatus);
2193    modelPtr_->setObjectiveValue(obj*modelPtr_->optimizationDirection());
2194    double * solution = modelPtr_->primalColumnSolution();
2195    const double * solution2 = modelPtr_->solutionRegion();
2196    // could just do changed bounds - also try double size scale so can use * not /
2197    if (!columnScale) {
2198      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2199        solution[iColumn]= solution2[iColumn];
2200      }
2201    } else {
2202      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2203        solution[iColumn]= solution2[iColumn]*columnScale[iColumn];
2204      }
2205    }
2206    CoinMemcpyN(saveLowerOriginal,numberColumns,columnLower);
2207    CoinMemcpyN(saveUpperOriginal,numberColumns,columnUpper);
2208#if 0
2209    // could combine with loop above
2210    if (modelPtr_==modelPtr_)
2211      modelPtr_->computeObjectiveValue();
2212    if (status&&!problemStatus) {
2213      memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2214      modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2215      modelPtr_->checkSolutionInternal();
2216      //modelPtr_->setLogLevel(1);
2217      //modelPtr_->allSlackBasis();
2218      //modelPtr_->primal(1);
2219      //memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2220      //modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2221      //modelPtr_->checkSolutionInternal();
2222      assert (!modelPtr_->problemStatus());
2223    }
2224#endif
2225    // and back bounds
2226    CoinMemcpyN(saveLower,number,modelPtr_->lowerRegion());
2227    CoinMemcpyN(saveUpper,number,modelPtr_->upperRegion());
2228    if (extraCopy) {
2229      CoinMemcpyN(saveLower,number,modelPtr_->lowerRegion()+extraCopy);
2230      CoinMemcpyN(saveUpper,number,modelPtr_->upperRegion()+extraCopy);
2231    }
2232    modelPtr_->setIntParam(ClpMaxNumIteration,itlimOrig_);
2233    if (savedObjective) {
2234      // fix up
2235      modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2236      //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2237      modelPtr_->objective_=savedObjective;
2238      if (!modelPtr_->problemStatus_) {
2239        CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
2240        CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
2241        if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
2242          CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
2243        modelPtr_->computeObjectiveValue();
2244      }
2245    }
2246    return;
2247  }
2248#endif
2249  if (smallModel_==NULL) {
2250    setWarmStart(ws_);
2251    CoinMemcpyN(           rowActivity_,numberRows,modelPtr_->primalRowSolution());
2252    CoinMemcpyN(columnActivity_,numberColumns,modelPtr_->primalColumnSolution());
2253    modelPtr_->setIntParam(ClpMaxNumIteration,CoinMin(itlim,9999));
2254    resolve();
2255  } else {
2256    assert (spareArrays_);
2257    double * arrayD = reinterpret_cast<double *> (spareArrays_);
2258    double saveObjectiveValue = arrayD[0];
2259    double * saveSolution = arrayD+1;
2260    int numberRows2 = smallModel_->numberRows();
2261    int numberColumns2 = smallModel_->numberColumns();
2262    int number = numberRows2+numberColumns2;
2263    CoinMemcpyN(saveSolution,number,smallModel_->solutionRegion());
2264    double * saveLower = saveSolution + (numberRows+numberColumns);
2265    CoinMemcpyN(saveLower,number,smallModel_->lowerRegion());
2266    double * saveUpper = saveLower + (numberRows+numberColumns);
2267    CoinMemcpyN(saveUpper,number,smallModel_->upperRegion());
2268    double * saveObjective = saveUpper + (numberRows+numberColumns);
2269    CoinMemcpyN(saveObjective,number,smallModel_->costRegion());
2270    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
2271    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
2272    arrayD = saveUpperOriginal + numberColumns;
2273    int * savePivot = reinterpret_cast<int *> (arrayD);
2274    CoinMemcpyN(savePivot,numberRows2,smallModel_->pivotVariable());
2275    int * whichRow = savePivot+numberRows;
2276    int * whichColumn = whichRow + 3*numberRows;
2277    int * arrayI = whichColumn + 2*numberColumns;
2278    unsigned char * saveStatus = reinterpret_cast<unsigned char *> (arrayI+1);
2279    CoinMemcpyN(saveStatus,number,smallModel_->statusArray());
2280    smallModel_->setFactorization(*factorization_);
2281    //int * backColumn = whichColumn+numberColumns;
2282    const double * lowerBig = modelPtr_->columnLower();
2283    const double * upperBig = modelPtr_->columnUpper();
2284    // make sure whatsChanged_ has 1 set
2285    smallModel_->setWhatsChanged(511);
2286    double * lowerSmall = smallModel_->lowerRegion();
2287    double * upperSmall = smallModel_->upperRegion();
2288    double * lowerSmallReal = smallModel_->columnLower();
2289    double * upperSmallReal = smallModel_->columnUpper();
2290    int i;
2291    double rhsScale = smallModel_->rhsScale();
2292    const double * columnScale = NULL;
2293    if (smallModel_->scalingFlag()>0) {
2294      columnScale = smallModel_->columnScale();
2295    }
2296    // and do bounds in case dual needs them
2297    // may be infeasible
2298    for (i=0;i<numberColumns2;i++) {
2299      int iColumn = whichColumn[i];
2300      if (lowerBig[iColumn]>saveLowerOriginal[iColumn]) {
2301        double value = lowerBig[iColumn];
2302        lowerSmallReal[i]=value;
2303        value *= rhsScale;
2304        if (columnScale)
2305          value /= columnScale[i];
2306        lowerSmall[i]=value;
2307      }
2308      if (upperBig[iColumn]<saveUpperOriginal[iColumn]) {
2309        double value = upperBig[iColumn];
2310        upperSmallReal[i]=value;
2311        value *= rhsScale;
2312        if (columnScale)
2313          value /= columnScale[i];
2314        upperSmall[i]=value;
2315      }
2316      if (upperSmall[i]<lowerSmall[i]-1.0e-8)
2317        break;
2318    }
2319    bool infeasible = (i<numberColumns2);
2320    // Start of fast iterations
2321    bool alwaysFinish= ((specialOptions_&32)==0) ? true : false;
2322    //smallModel_->setLogLevel(1);
2323    smallModel_->setIntParam(ClpMaxNumIteration,itlim);
2324    int saveNumberFake = (static_cast<ClpSimplexDual *>(smallModel_))->numberFake_;
2325    int status;
2326    if (!infeasible) {
2327      status = static_cast<ClpSimplexDual *>(smallModel_)->fastDual(alwaysFinish);
2328    } else {
2329      status=0;
2330      smallModel_->setProblemStatus(1);
2331    }
2332    (static_cast<ClpSimplexDual *>(smallModel_))->numberFake_ = saveNumberFake;
2333    if (smallModel_->numberIterations()==-98) {
2334      printf("rrrrrrrrrrrr\n");
2335      smallModel_->checkPrimalSolution(smallModel_->solutionRegion(0),
2336                                     smallModel_->solutionRegion(1));
2337      //smallModel_->gutsOfSolution(NULL,NULL,0);
2338      //if (problemStatus==3)
2339      //smallModel_->computeObjectiveValue();
2340      printf("robj %g\n",smallModel_->objectiveValue() *
2341             modelPtr_->optimizationDirection());
2342      writeMps("rr.mps");
2343      smallModel_->writeMps("rr_small.mps");
2344      ClpSimplex temp = *smallModel_;
2345      printf("small\n");
2346      temp.setLogLevel(63);
2347      temp.dual();
2348      double limit = 0.0;
2349      modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
2350      if (temp.problemStatus()==0&&temp.objectiveValue()<limit) {
2351        printf("inf obj %g, true %g - offsets %g %g\n",smallModel_->objectiveValue(),
2352               temp.objectiveValue(),
2353               smallModel_->objectiveOffset(),temp.objectiveOffset());
2354      }
2355      printf("big\n");
2356      temp = *modelPtr_;
2357      temp.dual();
2358      if (temp.problemStatus()==0&&temp.objectiveValue()<limit) {
2359        printf("inf obj %g, true %g - offsets %g %g\n",smallModel_->objectiveValue(),
2360               temp.objectiveValue(),
2361               smallModel_->objectiveOffset(),temp.objectiveOffset());
2362      }
2363    }
2364    int problemStatus = smallModel_->problemStatus();
2365    double objectiveValue =smallModel_->objectiveValue() * modelPtr_->optimizationDirection();
2366    CoinAssert (smallModel_->problemStatus()||smallModel_->objectiveValue()<1.0e50);
2367    // make sure plausible
2368    double obj = CoinMax(objectiveValue,saveObjectiveValue);
2369    if (problemStatus==10||problemStatus<0) {
2370      // was trying to clean up or something odd
2371      if (problemStatus==10)
2372        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2373      status=1;
2374    }
2375    if (status) {
2376      // not finished - might be optimal
2377      smallModel_->checkPrimalSolution(smallModel_->solutionRegion(0),
2378                                     smallModel_->solutionRegion(1));
2379      //smallModel_->gutsOfSolution(NULL,NULL,0);
2380      //if (problemStatus==3)
2381      //smallModel_->computeObjectiveValue();
2382      objectiveValue =smallModel_->objectiveValue() *
2383        modelPtr_->optimizationDirection();
2384      if (problemStatus!=10)
2385        obj = CoinMax(objectiveValue,saveObjectiveValue);
2386      if (!smallModel_->numberDualInfeasibilities()) { 
2387        double limit = 0.0;
2388        modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
2389        if (smallModel_->secondaryStatus()==1&&!problemStatus&&obj<limit) {
2390#if 0
2391          // switch off
2392          ClpSimplex temp = *smallModel_;
2393          temp.dual();
2394          if (temp.problemStatus()==0&&temp.objectiveValue()<limit) {
2395            printf("inf obj %g, true %g - offsets %g %g\n",smallModel_->objectiveValue(),
2396                   temp.objectiveValue(),
2397                   smallModel_->objectiveOffset(),temp.objectiveOffset());
2398          }
2399          lastAlgorithm_=1;
2400          obj=limit;
2401          problemStatus=10;
2402#else
2403          obj=limit;
2404          problemStatus=3;
2405#endif
2406        }
2407        if (!smallModel_->numberPrimalInfeasibilities()&&obj<limit) { 
2408          problemStatus=0;
2409#if 0
2410          ClpSimplex temp = *smallModel_;
2411          temp.dual();
2412          if (temp.numberIterations())
2413            printf("temp iterated\n");
2414          assert (temp.problemStatus()==0&&temp.objectiveValue()<limit);
2415#endif
2416        } else if (problemStatus==10) {
2417          problemStatus=3;
2418        } else if (!smallModel_->numberPrimalInfeasibilities()) {
2419          problemStatus=1; // infeasible
2420        } 
2421      } else {
2422        // can't say much
2423        //if (problemStatus==3)
2424        //smallModel_->computeObjectiveValue();
2425        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2426        problemStatus=3;
2427      }
2428    } else if (!problemStatus) {
2429      if (smallModel_->isDualObjectiveLimitReached()) 
2430        problemStatus=1; // infeasible
2431    }
2432    if (status&&!problemStatus) {
2433      problemStatus=3; // can't be sure
2434      lastAlgorithm_=1;
2435    }
2436    if (problemStatus<0)
2437      problemStatus=3;
2438    modelPtr_->setProblemStatus(problemStatus);
2439    modelPtr_->setObjectiveValue(obj*modelPtr_->optimizationDirection());
2440    modelPtr_->setSumDualInfeasibilities(smallModel_->sumDualInfeasibilities());
2441    modelPtr_->setNumberDualInfeasibilities(smallModel_->numberDualInfeasibilities());
2442    modelPtr_->setSumPrimalInfeasibilities(smallModel_->sumPrimalInfeasibilities());
2443    modelPtr_->setNumberPrimalInfeasibilities(smallModel_->numberPrimalInfeasibilities());
2444    double * solution = modelPtr_->primalColumnSolution();
2445    const double * solution2 = smallModel_->solutionRegion();
2446    if (!columnScale) {
2447      for (i=0;i<numberColumns2;i++) {
2448        int iColumn = whichColumn[i];
2449        solution[iColumn]= solution2[i];
2450        lowerSmallReal[i]=saveLowerOriginal[iColumn];
2451        upperSmallReal[i]=saveUpperOriginal[iColumn];
2452      }
2453    } else {
2454      for (i=0;i<numberColumns2;i++) {
2455        int iColumn = whichColumn[i];
2456        solution[iColumn]= solution2[i]*columnScale[i];
2457        lowerSmallReal[i]=saveLowerOriginal[iColumn];
2458        upperSmallReal[i]=saveUpperOriginal[iColumn];
2459      }
2460    }
2461    // could combine with loop above
2462    if (modelPtr_==smallModel_)
2463      modelPtr_->computeObjectiveValue();
2464#if 1
2465    if (status&&!problemStatus) {
2466      memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2467      modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2468      modelPtr_->checkSolutionInternal();
2469      //modelPtr_->setLogLevel(1);
2470      //modelPtr_->allSlackBasis();
2471      //modelPtr_->primal(1);
2472      //memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2473      //modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2474      //modelPtr_->checkSolutionInternal();
2475      assert (!modelPtr_->problemStatus());
2476    }
2477#endif
2478    modelPtr_->setNumberIterations(smallModel_->numberIterations());
2479    // and back bounds
2480    CoinMemcpyN(saveLower,number,smallModel_->lowerRegion());
2481    CoinMemcpyN(saveUpper,number,smallModel_->upperRegion());
2482  }
2483  if (savedObjective) {
2484    // fix up
2485    modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2486    //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2487    modelPtr_->objective_=savedObjective;
2488    if (!modelPtr_->problemStatus_) {
2489      CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
2490      CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
2491      if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
2492        CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
2493      modelPtr_->computeObjectiveValue();
2494    }
2495  }
2496  modelPtr_->setIntParam(ClpMaxNumIteration,itlimOrig_);
2497}
2498
2499void OsiClpSolverInterface::unmarkHotStart()
2500{
2501#ifdef CLEAN_HOT_START
2502  if ((specialOptions_&65536)!=0) {
2503    modelPtr_->setLogLevel(saveData_.scalingFlag_);
2504    modelPtr_->deleteRim(0);
2505    if (lastNumberRows_<0) {
2506      specialOptions_ |= 131072;
2507      lastNumberRows_ = -1 -lastNumberRows_;
2508      if (modelPtr_->rowScale_) {
2509        if (modelPtr_->rowScale_!=rowScale_.array()) {
2510          delete [] modelPtr_->rowScale_;
2511          delete [] modelPtr_->columnScale_;
2512        }
2513        modelPtr_->rowScale_=NULL;
2514        modelPtr_->columnScale_=NULL;
2515      }
2516    }
2517    delete factorization_;
2518    delete [] spareArrays_;
2519    smallModel_=NULL;
2520    spareArrays_=NULL;
2521    factorization_=NULL;
2522    delete [] rowActivity_;
2523    delete [] columnActivity_;
2524    rowActivity_=NULL;
2525    columnActivity_=NULL;
2526    return;
2527  }
2528#endif
2529  if (smallModel_==NULL) {
2530    setWarmStart(ws_);
2531    int numberRows = modelPtr_->numberRows();
2532    int numberColumns = modelPtr_->numberColumns();
2533    CoinMemcpyN(           rowActivity_,numberRows,modelPtr_->primalRowSolution());
2534    CoinMemcpyN(columnActivity_,numberColumns,modelPtr_->primalColumnSolution());
2535    delete ws_;
2536    ws_ = NULL;
2537  } else {
2538#ifndef KEEP_SMALL
2539    if (smallModel_!=modelPtr_)
2540      delete smallModel_;
2541    smallModel_=NULL;
2542#else
2543    if (smallModel_==modelPtr_) {
2544      smallModel_=NULL;
2545    } else if (smallModel_) {
2546      if(!spareArrays_) {
2547        delete smallModel_;
2548        smallModel_=NULL;
2549        delete factorization_;
2550        factorization_=NULL;
2551      } else {
2552        static_cast<ClpSimplexDual *> (smallModel_)->cleanupAfterStrongBranching( factorization_);
2553        if ((smallModel_->specialOptions_&4096)==0) {
2554          delete factorization_;
2555        }
2556      }
2557    }
2558#endif
2559    //delete [] spareArrays_;
2560    //spareArrays_=NULL;
2561    factorization_=NULL;
2562  }
2563  delete [] rowActivity_;
2564  delete [] columnActivity_;
2565  rowActivity_=NULL;
2566  columnActivity_=NULL;
2567  // Make sure whatsChanged not out of sync
2568  if (!modelPtr_->columnUpperWork_)
2569    modelPtr_->whatsChanged_ &= ~0xffff;
2570  modelPtr_->specialOptions_ = saveData_.specialOptions_; 
2571}
2572
2573//#############################################################################
2574// Problem information methods (original data)
2575//#############################################################################
2576
2577//------------------------------------------------------------------
2578const char * OsiClpSolverInterface::getRowSense() const
2579{
2580  extractSenseRhsRange();
2581  return rowsense_;
2582}
2583//------------------------------------------------------------------
2584const double * OsiClpSolverInterface::getRightHandSide() const
2585{
2586  extractSenseRhsRange();
2587  return rhs_;
2588}
2589//------------------------------------------------------------------
2590const double * OsiClpSolverInterface::getRowRange() const
2591{
2592  extractSenseRhsRange();
2593  return rowrange_;
2594}
2595//------------------------------------------------------------------
2596// Return information on integrality
2597//------------------------------------------------------------------
2598bool OsiClpSolverInterface::isContinuous(int colNumber) const
2599{
2600  if ( integerInformation_==NULL ) return true;
2601#ifndef NDEBUG
2602  int n = modelPtr_->numberColumns();
2603  if (colNumber<0||colNumber>=n) {
2604    indexError(colNumber,"isContinuous");
2605  }
2606#endif
2607  if ( integerInformation_[colNumber]==0 ) return true;
2608  return false;
2609}
2610bool 
2611OsiClpSolverInterface::isBinary(int colNumber) const
2612{
2613#ifndef NDEBUG
2614  int n = modelPtr_->numberColumns();
2615  if (colNumber<0||colNumber>=n) {
2616    indexError(colNumber,"isBinary");
2617  }
2618#endif
2619  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) {
2620    return false;
2621  } else {
2622    const double * cu = getColUpper();
2623    const double * cl = getColLower();
2624    if ((cu[colNumber]== 1 || cu[colNumber]== 0) && 
2625        (cl[colNumber]== 0 || cl[colNumber]==1))
2626      return true;
2627    else 
2628      return false;
2629  }
2630}
2631//-----------------------------------------------------------------------------
2632bool 
2633OsiClpSolverInterface::isInteger(int colNumber) const
2634{
2635#ifndef NDEBUG
2636  int n = modelPtr_->numberColumns();
2637  if (colNumber<0||colNumber>=n) {
2638    indexError(colNumber,"isInteger");
2639  }
2640#endif
2641  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) 
2642    return false;
2643  else
2644    return true;
2645}
2646//-----------------------------------------------------------------------------
2647bool 
2648OsiClpSolverInterface::isIntegerNonBinary(int colNumber) const
2649{
2650#ifndef NDEBUG
2651  int n = modelPtr_->numberColumns();
2652  if (colNumber<0||colNumber>=n) {
2653    indexError(colNumber,"isIntegerNonBinary");
2654  }
2655#endif
2656  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) {
2657    return false;
2658  } else {
2659    return !isBinary(colNumber);
2660  }
2661}
2662//-----------------------------------------------------------------------------
2663bool 
2664OsiClpSolverInterface::isFreeBinary(int colNumber) const
2665{
2666#ifndef NDEBUG
2667  int n = modelPtr_->numberColumns();
2668  if (colNumber<0||colNumber>=n) {
2669    indexError(colNumber,"isFreeBinary");
2670  }
2671#endif
2672  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) {
2673    return false;
2674  } else {
2675    const double * cu = getColUpper();
2676    const double * cl = getColLower();
2677    if ((cu[colNumber]== 1) && (cl[colNumber]== 0))
2678      return true;
2679    else 
2680      return false;
2681  }
2682}
2683/*  Return array of column length
2684    0 - continuous
2685    1 - binary (may get fixed later)
2686    2 - general integer (may get fixed later)
2687*/
2688const char * 
2689OsiClpSolverInterface::getColType(bool refresh) const
2690{
2691  if (!columnType_||refresh) {
2692    const int numCols = getNumCols() ;
2693    if (!columnType_)
2694      columnType_ = new char [numCols];
2695    if ( integerInformation_==NULL ) {
2696      memset(columnType_,0,numCols);
2697    } else {
2698      const double * cu = getColUpper();
2699      const double * cl = getColLower();
2700      for (int i = 0 ; i < numCols ; ++i) {
2701        if (integerInformation_[i]) {
2702          if ((cu[i]== 1 || cu[i]== 0) && 
2703              (cl[i]== 0 || cl[i]==1))
2704            columnType_[i]=1;
2705          else
2706            columnType_[i]=2;
2707        } else {
2708          columnType_[i]=0;
2709        }
2710      }
2711    }
2712  }
2713  return columnType_;
2714}
2715
2716/* Return true if column is integer but does not have to
2717   be declared as such.
2718   Note: This function returns true if the the column
2719   is binary or a general integer.
2720*/
2721bool 
2722OsiClpSolverInterface::isOptionalInteger(int colNumber) const
2723{
2724#ifndef NDEBUG
2725  int n = modelPtr_->numberColumns();
2726  if (colNumber<0||colNumber>=n) {
2727    indexError(colNumber,"isInteger");
2728  }
2729#endif
2730  if ( integerInformation_==NULL || integerInformation_[colNumber]!=2 ) 
2731    return false;
2732  else
2733    return true;
2734}
2735/* Set the index-th variable to be an optional integer variable */
2736void 
2737OsiClpSolverInterface::setOptionalInteger(int index)
2738{
2739  if (!integerInformation_) {
2740    integerInformation_ = new char[modelPtr_->numberColumns()];
2741    CoinFillN ( integerInformation_, modelPtr_->numberColumns(),static_cast<char> (0));
2742  }
2743#ifndef NDEBUG
2744  int n = modelPtr_->numberColumns();
2745  if (index<0||index>=n) {
2746    indexError(index,"setInteger");
2747  }
2748#endif
2749  integerInformation_[index]=2;
2750  modelPtr_->setInteger(index);
2751}
2752
2753//------------------------------------------------------------------
2754// Row and column copies of the matrix ...
2755//------------------------------------------------------------------
2756const CoinPackedMatrix * OsiClpSolverInterface::getMatrixByRow() const
2757{
2758  if ( matrixByRow_ == NULL ) {
2759    matrixByRow_ = new CoinPackedMatrix();
2760    matrixByRow_->setExtraGap(0.0);
2761    matrixByRow_->setExtraMajor(0.0);
2762    matrixByRow_->reverseOrderedCopyOf(*modelPtr_->matrix());
2763    //matrixByRow_->removeGaps();
2764#if 0
2765    CoinPackedMatrix back;
2766    std::cout<<"start check"<<std::endl;
2767    back.reverseOrderedCopyOf(*matrixByRow_);
2768    modelPtr_->matrix()->isEquivalent2(back);
2769    std::cout<<"stop check"<<std::endl;
2770#endif
2771  }
2772  assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
2773  return matrixByRow_;
2774}
2775
2776const CoinPackedMatrix * OsiClpSolverInterface::getMatrixByCol() const
2777{
2778  return modelPtr_->matrix();
2779}
2780 
2781// Get pointer to mutable column-wise copy of matrix (returns NULL if not meaningful)
2782CoinPackedMatrix * 
2783OsiClpSolverInterface::getMutableMatrixByCol() const 
2784{
2785  ClpPackedMatrix * matrix = dynamic_cast<ClpPackedMatrix *>(modelPtr_->matrix_) ;
2786  if (matrix)
2787    return matrix->getPackedMatrix();
2788  else
2789    return NULL;
2790}
2791
2792//------------------------------------------------------------------
2793std::vector<double*> OsiClpSolverInterface::getDualRays(int maxNumRays,
2794                                                        bool fullRay) const
2795{
2796  if (fullRay == true) {
2797    throw CoinError("Full dual rays not yet implemented.","getDualRays",
2798                    "OsiClpSolverInterface");
2799  }
2800  return std::vector<double*>(1, modelPtr_->infeasibilityRay());
2801}
2802//------------------------------------------------------------------
2803std::vector<double*> OsiClpSolverInterface::getPrimalRays(int /*maxNumRays*/) const
2804{
2805  return std::vector<double*>(1, modelPtr_->unboundedRay());
2806}
2807//#############################################################################
2808void
2809OsiClpSolverInterface::setContinuous(int index)
2810{
2811
2812  if (integerInformation_) {
2813#ifndef NDEBUG
2814    int n = modelPtr_->numberColumns();
2815    if (index<0||index>=n) {
2816      indexError(index,"setContinuous");
2817    }
2818#endif
2819    integerInformation_[index]=0;
2820  }
2821  modelPtr_->setContinuous(index);
2822}
2823//-----------------------------------------------------------------------------
2824void
2825OsiClpSolverInterface::setInteger(int index)
2826{
2827  if (!integerInformation_) {
2828    integerInformation_ = new char[modelPtr_->numberColumns()];
2829    CoinFillN ( integerInformation_, modelPtr_->numberColumns(),static_cast<char> (0));
2830  }
2831#ifndef NDEBUG
2832  int n = modelPtr_->numberColumns();
2833  if (index<0||index>=n) {
2834    indexError(index,"setInteger");
2835  }
2836#endif
2837  integerInformation_[index]=1;
2838  modelPtr_->setInteger(index);
2839}
2840//-----------------------------------------------------------------------------
2841void
2842OsiClpSolverInterface::setContinuous(const int* indices, int len)
2843{
2844  if (integerInformation_) {
2845#ifndef NDEBUG
2846    int n = modelPtr_->numberColumns();
2847#endif
2848    int i;
2849    for (i=0; i<len;i++) {
2850      int colNumber = indices[i];
2851#ifndef NDEBUG
2852      if (colNumber<0||colNumber>=n) {
2853        indexError(colNumber,"setContinuous");
2854      }
2855#endif
2856      integerInformation_[colNumber]=0;
2857      modelPtr_->setContinuous(colNumber);
2858    }
2859  }
2860}
2861//-----------------------------------------------------------------------------
2862void
2863OsiClpSolverInterface::setInteger(const int* indices, int len)
2864{
2865  if (!integerInformation_) {
2866    integerInformation_ = new char[modelPtr_->numberColumns()];
2867    CoinFillN ( integerInformation_, modelPtr_->numberColumns(),static_cast<char> (0));
2868  }
2869#ifndef NDEBUG
2870  int n = modelPtr_->numberColumns();
2871#endif
2872  int i;
2873  for (i=0; i<len;i++) {
2874    int colNumber = indices[i];
2875#ifndef NDEBUG
2876    if (colNumber<0||colNumber>=n) {
2877      indexError(colNumber,"setInteger");
2878    }
2879#endif
2880    integerInformation_[colNumber]=1;
2881    modelPtr_->setInteger(colNumber);
2882  }
2883}
2884/* Set the objective coefficients for all columns
2885    array [getNumCols()] is an array of values for the objective.
2886    This defaults to a series of set operations and is here for speed.
2887*/
2888void 
2889OsiClpSolverInterface::setObjective(const double * array)
2890{
2891  // Say can't gurantee optimal basis etc
2892  lastAlgorithm_=999;
2893  modelPtr_->whatsChanged_ &= (0xffff&~64);
2894  CoinMemcpyN(array,modelPtr_->numberColumns(),
2895                    modelPtr_->objective());
2896}
2897/* Set the lower bounds for all columns
2898    array [getNumCols()] is an array of values for the objective.
2899    This defaults to a series of set operations and is here for speed.
2900*/
2901void 
2902OsiClpSolverInterface::setColLower(const double * array)
2903{
2904  // Say can't gurantee optimal basis etc
2905  lastAlgorithm_=999;
2906  modelPtr_->whatsChanged_ &= (0x1ffff&128);
2907  CoinMemcpyN(array,modelPtr_->numberColumns(),
2908                    modelPtr_->columnLower());
2909}
2910/* Set the upper bounds for all columns
2911    array [getNumCols()] is an array of values for the objective.
2912    This defaults to a series of set operations and is here for speed.
2913*/
2914void 
2915OsiClpSolverInterface::setColUpper(const double * array)
2916{
2917  // Say can't gurantee optimal basis etc
2918  lastAlgorithm_=999;
2919  modelPtr_->whatsChanged_ &= (0x1ffff&256);
2920  CoinMemcpyN(array,modelPtr_->numberColumns(),
2921                    modelPtr_->columnUpper());
2922}
2923//-----------------------------------------------------------------------------
2924void OsiClpSolverInterface::setColSolution(const double * cs) 
2925{
2926  // Say can't gurantee optimal basis etc
2927  lastAlgorithm_=999;
2928  CoinDisjointCopyN(cs,modelPtr_->numberColumns(),
2929                    modelPtr_->primalColumnSolution());
2930  if (modelPtr_->solveType()==2) {
2931    // directly into code as well
2932    CoinDisjointCopyN(cs,modelPtr_->numberColumns(),
2933                      modelPtr_->solutionRegion(1));
2934  }
2935}
2936//-----------------------------------------------------------------------------
2937void OsiClpSolverInterface::setRowPrice(const double * rs) 
2938{
2939  CoinDisjointCopyN(rs,modelPtr_->numberRows(),
2940                    modelPtr_->dualRowSolution());
2941  if (modelPtr_->solveType()==2) {
2942    // directly into code as well (? sign )
2943    CoinDisjointCopyN(rs,modelPtr_->numberRows(),
2944                      modelPtr_->djRegion(0));
2945  }
2946}
2947
2948//#############################################################################
2949// Problem modifying methods (matrix)
2950//#############################################################################
2951void 
2952OsiClpSolverInterface::addCol(const CoinPackedVectorBase& vec,
2953                              const double collb, const double colub,   
2954                              const double obj)
2955{
2956  int numberColumns = modelPtr_->numberColumns();
2957  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
2958  modelPtr_->resize(modelPtr_->numberRows(),numberColumns+1);
2959  linearObjective_ = modelPtr_->objective();
2960  basis_.resize(modelPtr_->numberRows(),numberColumns+1);
2961  setColBounds(numberColumns,collb,colub);
2962  setObjCoeff(numberColumns,obj);
2963  if (!modelPtr_->clpMatrix())
2964    modelPtr_->createEmptyMatrix();
2965  modelPtr_->matrix()->appendCol(vec);
2966  if (integerInformation_) {
2967    char * temp = new char[numberColumns+1];
2968    CoinMemcpyN(integerInformation_,numberColumns,temp);
2969    delete [] integerInformation_;
2970    integerInformation_ = temp;
2971    integerInformation_[numberColumns]=0;
2972  }
2973  freeCachedResults();
2974}
2975//-----------------------------------------------------------------------------
2976/* Add a column (primal variable) to the problem. */
2977void 
2978OsiClpSolverInterface::addCol(int numberElements, const int * rows, const double * elements,
2979                           const double collb, const double colub,   
2980                           const double obj) 
2981{
2982  CoinPackedVector column(numberElements, rows, elements);
2983  addCol(column,collb,colub,obj);
2984}
2985// Add a named column (primal variable) to the problem.
2986void 
2987OsiClpSolverInterface::addCol(const CoinPackedVectorBase& vec,
2988                      const double collb, const double colub,   
2989                      const double obj, std::string name) 
2990{
2991  int ndx = getNumCols() ;
2992  addCol(vec,collb,colub,obj) ;
2993  setColName(ndx,name) ;
2994}
2995//-----------------------------------------------------------------------------
2996void 
2997OsiClpSolverInterface::addCols(const int numcols,
2998                               const CoinPackedVectorBase * const * cols,
2999                               const double* collb, const double* colub,   
3000                               const double* obj)
3001{
3002  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3003  int numberColumns = modelPtr_->numberColumns();
3004  modelPtr_->resize(modelPtr_->numberRows(),numberColumns+numcols);
3005  linearObjective_ = modelPtr_->objective();
3006  basis_.resize(modelPtr_->numberRows(),numberColumns+numcols);
3007  double * lower = modelPtr_->columnLower()+numberColumns;
3008  double * upper = modelPtr_->columnUpper()+numberColumns;
3009  double * objective = modelPtr_->objective()+numberColumns;
3010  int iCol;
3011  if (collb) {
3012    for (iCol = 0; iCol < numcols; iCol++) {
3013      lower[iCol]= forceIntoRange(collb[iCol], -OsiClpInfinity, OsiClpInfinity);
3014      if (lower[iCol]<-1.0e27)
3015        lower[iCol]=-COIN_DBL_MAX;
3016    }
3017  } else {
3018    CoinFillN ( lower, numcols,0.0);
3019  }
3020  if (colub) {
3021    for (iCol = 0; iCol < numcols; iCol++) {
3022      upper[iCol]= forceIntoRange(colub[iCol], -OsiClpInfinity, OsiClpInfinity);
3023      if (upper[iCol]>1.0e27)
3024        upper[iCol]=COIN_DBL_MAX;
3025    }
3026  } else {
3027    CoinFillN ( upper, numcols,COIN_DBL_MAX);
3028  }
3029  if (obj) {
3030    for (iCol = 0; iCol < numcols; iCol++) {
3031      objective[iCol] = obj[iCol];
3032    }
3033  } else {
3034    CoinFillN ( objective, numcols,0.0);
3035  }
3036  if (!modelPtr_->clpMatrix())
3037    modelPtr_->createEmptyMatrix();
3038  modelPtr_->matrix()->appendCols(numcols,cols);
3039  if (integerInformation_) {
3040    char * temp = new char[numberColumns+numcols];
3041    CoinMemcpyN(integerInformation_,numberColumns,temp);
3042    delete [] integerInformation_;
3043    integerInformation_ = temp;
3044    for (iCol = 0; iCol < numcols; iCol++) 
3045      integerInformation_[numberColumns+iCol]=0;
3046  }
3047  freeCachedResults();
3048}
3049void 
3050OsiClpSolverInterface::addCols(const int numcols,
3051                               const int * columnStarts, const int * rows, const double * elements,
3052                               const double* collb, const double* colub,   
3053                               const double* obj)
3054{
3055  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3056  int numberColumns = modelPtr_->numberColumns();
3057  modelPtr_->resize(modelPtr_->numberRows(),numberColumns+numcols);
3058  linearObjective_ = modelPtr_->objective();
3059  basis_.resize(modelPtr_->numberRows(),numberColumns+numcols);
3060  double * lower = modelPtr_->columnLower()+numberColumns;
3061  double * upper = modelPtr_->columnUpper()+numberColumns;
3062  double * objective = modelPtr_->objective()+numberColumns;
3063  int iCol;
3064  if (collb) {
3065    for (iCol = 0; iCol < numcols; iCol++) {
3066      lower[iCol]= forceIntoRange(collb[iCol], -OsiClpInfinity, OsiClpInfinity);
3067      if (lower[iCol]<-1.0e27)
3068        lower[iCol]=-COIN_DBL_MAX;
3069    }
3070  } else {
3071    CoinFillN ( lower, numcols,0.0);
3072  }
3073  if (colub) {
3074    for (iCol = 0; iCol < numcols; iCol++) {
3075      upper[iCol]= forceIntoRange(colub[iCol], -OsiClpInfinity, OsiClpInfinity);
3076      if (upper[iCol]>1.0e27)
3077        upper[iCol]=COIN_DBL_MAX;
3078    }
3079  } else {
3080    CoinFillN ( upper, numcols,COIN_DBL_MAX);
3081  }
3082  if (obj) {
3083    for (iCol = 0; iCol < numcols; iCol++) {
3084      objective[iCol] = obj[iCol];
3085    }
3086  } else {
3087    CoinFillN ( objective, numcols,0.0);
3088  }
3089  if (!modelPtr_->clpMatrix())
3090    modelPtr_->createEmptyMatrix();
3091  modelPtr_->matrix()->appendCols(numcols,columnStarts,rows,elements);
3092  if (integerInformation_) {
3093    char * temp = new char[numberColumns+numcols];
3094    CoinMemcpyN(integerInformation_,numberColumns,temp);
3095    delete [] integerInformation_;
3096    integerInformation_ = temp;
3097    for (iCol = 0; iCol < numcols; iCol++) 
3098      integerInformation_[numberColumns+iCol]=0;
3099  }
3100  freeCachedResults();
3101}
3102//-----------------------------------------------------------------------------
3103void 
3104OsiClpSolverInterface::deleteCols(const int num, const int * columnIndices)
3105{
3106  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3107  findIntegers(false);
3108  deleteBranchingInfo(num,columnIndices);
3109  modelPtr_->deleteColumns(num,columnIndices);
3110  int nameDiscipline;
3111  getIntParam(OsiNameDiscipline,nameDiscipline) ;
3112  if (num&&nameDiscipline) {
3113    // Very clumsy (and inefficient) - need to sort and then go backwards in ? chunks
3114    int * indices = CoinCopyOfArray(columnIndices,num);
3115    std::sort(indices,indices+num);
3116    int num2=num;
3117    while(num2) {
3118      int next = indices[num2-1];
3119      int firstDelete = num2-1;
3120      int i;
3121      for (i=num2-2;i>=0;i--) {
3122        if (indices[i]+1==next) {
3123          next --;
3124          firstDelete=i;
3125        } else {
3126          break;
3127        }
3128      }
3129      OsiSolverInterface::deleteColNames(indices[firstDelete],num2-firstDelete);
3130      num2 = firstDelete;
3131      assert (num2>=0);
3132    }
3133    delete [] indices;
3134  }
3135  // synchronize integers (again)
3136  if (integerInformation_) {
3137    int numberColumns = modelPtr_->numberColumns();
3138    for (int i=0;i<numberColumns;i++) {
3139      if (modelPtr_->isInteger(i))
3140        integerInformation_[i]=1;
3141      else
3142        integerInformation_[i]=0;
3143    }
3144  }
3145  basis_.deleteColumns(num,columnIndices);
3146  linearObjective_ = modelPtr_->objective();
3147  freeCachedResults();
3148}
3149//-----------------------------------------------------------------------------
3150void 
3151OsiClpSolverInterface::addRow(const CoinPackedVectorBase& vec,
3152                              const double rowlb, const double rowub)
3153{
3154  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3155  freeCachedResults0();
3156  int numberRows = modelPtr_->numberRows();
3157  modelPtr_->resize(numberRows+1,modelPtr_->numberColumns());
3158  basis_.resize(numberRows+1,modelPtr_->numberColumns());
3159  setRowBounds(numberRows,rowlb,rowub);
3160  if (!modelPtr_->clpMatrix())
3161    modelPtr_->createEmptyMatrix();
3162  modelPtr_->matrix()->appendRow(vec);
3163  freeCachedResults1();
3164}
3165//-----------------------------------------------------------------------------
3166void OsiClpSolverInterface::addRow(const CoinPackedVectorBase& vec,
3167                                const double rowlb, const double rowub,
3168                                std::string name)
3169{
3170  int ndx = getNumRows() ;
3171  addRow(vec,rowlb,rowub) ;
3172  setRowName(ndx,name) ;
3173}
3174//-----------------------------------------------------------------------------
3175void 
3176OsiClpSolverInterface::addRow(const CoinPackedVectorBase& vec,
3177                              const char rowsen, const double rowrhs,   
3178                              const double rowrng)
3179{
3180  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3181  freeCachedResults0();
3182  int numberRows = modelPtr_->numberRows();
3183  modelPtr_->resize(numberRows+1,modelPtr_->numberColumns());
3184  basis_.resize(numberRows+1,modelPtr_->numberColumns());
3185  double rowlb = 0, rowub = 0;
3186  convertSenseToBound(rowsen, rowrhs, rowrng, rowlb, rowub);
3187  setRowBounds(numberRows,rowlb,rowub);
3188  if (!modelPtr_->clpMatrix())
3189    modelPtr_->createEmptyMatrix();
3190  modelPtr_->matrix()->appendRow(vec);
3191  freeCachedResults1();
3192}
3193//-----------------------------------------------------------------------------
3194void 
3195OsiClpSolverInterface::addRow(int numberElements, const int * columns, const double * elements,
3196                           const double rowlb, const double rowub) 
3197{
3198  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3199  freeCachedResults0();
3200  int numberRows = modelPtr_->numberRows();
3201  modelPtr_->resize(numberRows+1,modelPtr_->numberColumns());
3202  basis_.resize(numberRows+1,modelPtr_->numberColumns());
3203  setRowBounds(numberRows,rowlb,rowub);
3204  if (!modelPtr_->clpMatrix())
3205    modelPtr_->createEmptyMatrix();
3206  modelPtr_->matrix()->appendRow(numberElements, columns, elements);
3207  CoinBigIndex starts[2];
3208  starts[0]=0;
3209  starts[1]=numberElements;
3210  redoScaleFactors( 1,starts, columns, elements);
3211  freeCachedResults1();
3212}
3213//-----------------------------------------------------------------------------
3214void 
3215OsiClpSolverInterface::addRows(const int numrows,
3216                               const CoinPackedVectorBase * const * rows,
3217                               const double* rowlb, const double* rowub)
3218{
3219  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3220  freeCachedResults0();
3221  int numberRows = modelPtr_->numberRows();
3222  modelPtr_->resize(numberRows+numrows,modelPtr_->numberColumns());
3223  basis_.resize(numberRows+numrows,modelPtr_->numberColumns());
3224  double * lower = modelPtr_->rowLower()+numberRows;
3225  double * upper = modelPtr_->rowUpper()+numberRows;
3226  int iRow;
3227  for (iRow = 0; iRow < numrows; iRow++) {
3228    if (rowlb) 
3229      lower[iRow]= forceIntoRange(rowlb[iRow], -OsiClpInfinity, OsiClpInfinity);
3230    else 
3231      lower[iRow]=-OsiClpInfinity;
3232    if (rowub) 
3233      upper[iRow]= forceIntoRange(rowub[iRow], -OsiClpInfinity, OsiClpInfinity);
3234    else 
3235      upper[iRow]=OsiClpInfinity;
3236    if (lower[iRow]<-1.0e27)
3237      lower[iRow]=-COIN_DBL_MAX;
3238    if (upper[iRow]>1.0e27)
3239      upper[iRow]=COIN_DBL_MAX;
3240  }
3241  if (!modelPtr_->clpMatrix())
3242    modelPtr_->createEmptyMatrix();
3243  modelPtr_->matrix()->appendRows(numrows,rows);
3244  freeCachedResults1();
3245}
3246//-----------------------------------------------------------------------------
3247void 
3248OsiClpSolverInterface::addRows(const int numrows,
3249                               const CoinPackedVectorBase * const * rows,
3250                               const char* rowsen, const double* rowrhs,   
3251                               const double* rowrng)
3252{
3253  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3254  freeCachedResults0();
3255  int numberRows = modelPtr_->numberRows();
3256  modelPtr_->resize(numberRows+numrows,modelPtr_->numberColumns());
3257  basis_.resize(numberRows+numrows,modelPtr_->numberColumns());
3258  double * lower = modelPtr_->rowLower()+numberRows;
3259  double * upper = modelPtr_->rowUpper()+numberRows;
3260  int iRow;
3261  for (iRow = 0; iRow < numrows; iRow++) {
3262    double rowlb = 0, rowub = 0;
3263    convertSenseToBound(rowsen[iRow], rowrhs[iRow], rowrng[iRow], 
3264                        rowlb, rowub);
3265    lower[iRow]= forceIntoRange(rowlb, -OsiClpInfinity, OsiClpInfinity);
3266    upper[iRow]= forceIntoRange(rowub, -OsiClpInfinity, OsiClpInfinity);
3267    if (lower[iRow]<-1.0e27)
3268      lower[iRow]=-COIN_DBL_MAX;
3269    if (upper[iRow]>1.0e27)
3270      upper[iRow]=COIN_DBL_MAX;
3271  }
3272  if (!modelPtr_->clpMatrix())
3273    modelPtr_->createEmptyMatrix();
3274  modelPtr_->matrix()->appendRows(numrows,rows);
3275  freeCachedResults1();
3276}
3277void 
3278OsiClpSolverInterface::addRows(const int numrows,
3279                               const int * rowStarts, const int * columns, const double * element,
3280                               const double* rowlb, const double* rowub)
3281{
3282  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3283  freeCachedResults0();
3284  int numberRows = modelPtr_->numberRows();
3285  modelPtr_->resize(numberRows+numrows,modelPtr_->numberColumns());
3286  basis_.resize(numberRows+numrows,modelPtr_->numberColumns());
3287  double * lower = modelPtr_->rowLower()+numberRows;
3288  double * upper = modelPtr_->rowUpper()+numberRows;
3289  int iRow;
3290  for (iRow = 0; iRow < numrows; iRow++) {
3291    if (rowlb) 
3292      lower[iRow]= forceIntoRange(rowlb[iRow], -OsiClpInfinity, OsiClpInfinity);
3293    else 
3294      lower[iRow]=-OsiClpInfinity;
3295    if (rowub) 
3296      upper[iRow]= forceIntoRange(rowub[iRow], -OsiClpInfinity, OsiClpInfinity);
3297    else 
3298      upper[iRow]=OsiClpInfinity;
3299    if (lower[iRow]<-1.0e27)
3300      lower[iRow]=-COIN_DBL_MAX;
3301    if (upper[iRow]>1.0e27)
3302      upper[iRow]=COIN_DBL_MAX;
3303  }
3304  if (!modelPtr_->clpMatrix())
3305    modelPtr_->createEmptyMatrix();
3306  modelPtr_->matrix()->appendRows(numrows,rowStarts,columns,element);
3307  redoScaleFactors( numrows,rowStarts, columns, element);
3308  freeCachedResults1();
3309}
3310//-----------------------------------------------------------------------------
3311void 
3312OsiClpSolverInterface::deleteRows(const int num, const int * rowIndices)
3313{
3314  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3315  // will still be optimal if all rows basic
3316  bool allBasic=true;
3317  int numBasis = basis_.getNumArtificial();
3318  for (int i=0;i<num;i++) {
3319    int iRow = rowIndices[i];
3320    if (iRow<numBasis) {
3321      if (basis_.getArtifStatus(iRow)!=CoinWarmStartBasis::basic) {
3322        allBasic=false;
3323        break;
3324      }
3325    }
3326  }
3327  int saveAlgorithm = allBasic ? lastAlgorithm_ : 999;
3328  modelPtr_->deleteRows(num,rowIndices);
3329  int nameDiscipline;
3330  getIntParam(OsiNameDiscipline,nameDiscipline) ;
3331  if (num&&nameDiscipline) {
3332    // Very clumsy (and inefficient) - need to sort and then go backwards in ? chunks
3333    int * indices = CoinCopyOfArray(rowIndices,num);
3334    std::sort(indices,indices+num);
3335    int num2=num;
3336    while(num2) {
3337      int next = indices[num2-1];
3338      int firstDelete = num2-1;
3339      int i;
3340      for (i=num2-2;i>=0;i--) {
3341        if (indices[i]+1==next) {
3342          next --;
3343          firstDelete=i;
3344        } else {
3345          break;
3346        }
3347      }
3348      OsiSolverInterface::deleteRowNames(indices[firstDelete],num2-firstDelete);
3349      num2 = firstDelete;
3350      assert (num2>=0);
3351    }
3352    delete [] indices;
3353  }
3354  basis_.deleteRows(num,rowIndices);
3355  CoinPackedMatrix * saveRowCopy = matrixByRow_;
3356  matrixByRow_=NULL;
3357  freeCachedResults();
3358  modelPtr_->setNewRowCopy(NULL);
3359  delete modelPtr_->scaledMatrix_;
3360  modelPtr_->scaledMatrix_=NULL;
3361  if (saveRowCopy) {
3362#if 1
3363    matrixByRow_=saveRowCopy;
3364    matrixByRow_->deleteRows(num,rowIndices);
3365    assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
3366#else
3367    delete saveRowCopy;
3368#endif
3369  }
3370  lastAlgorithm_ = saveAlgorithm;
3371  if ((specialOptions_&131072)!=0) 
3372    lastNumberRows_=modelPtr_->numberRows();
3373}
3374
3375//#############################################################################
3376// Methods to input a problem
3377//#############################################################################
3378
3379void
3380OsiClpSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
3381                                   const double* collb, const double* colub,   
3382                                   const double* obj,
3383                                   const double* rowlb, const double* rowub)
3384{
3385  modelPtr_->whatsChanged_ = 0;
3386  // Get rid of integer information (modelPtr will get rid of its copy)
3387  delete [] integerInformation_;
3388  integerInformation_=NULL;
3389  modelPtr_->loadProblem(matrix, collb, colub, obj, rowlb, rowub);
3390  linearObjective_ = modelPtr_->objective();
3391  freeCachedResults();
3392  basis_=CoinWarmStartBasis();
3393  if (ws_) {
3394     delete ws_;
3395     ws_ = 0;
3396  }
3397}
3398
3399//-----------------------------------------------------------------------------
3400
3401void
3402OsiClpSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
3403                                     double*& collb, double*& colub,
3404                                     double*& obj,
3405                                     double*& rowlb, double*& rowub)
3406{
3407  modelPtr_->whatsChanged_ = 0;
3408  // Get rid of integer information (modelPtr will get rid of its copy)
3409  loadProblem(*matrix, collb, colub, obj, rowlb, rowub);
3410  delete matrix;   matrix = NULL;
3411  delete[] collb;  collb = NULL;
3412  delete[] colub;  colub = NULL;
3413  delete[] obj;    obj = NULL;
3414  delete[] rowlb;  rowlb = NULL;
3415  delete[] rowub;  rowub = NULL;
3416}
3417
3418//-----------------------------------------------------------------------------
3419
3420void
3421OsiClpSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
3422                                   const double* collb, const double* colub,
3423                                   const double* obj,
3424                                   const char* rowsen, const double* rowrhs,   
3425                                   const double* rowrng)
3426{
3427  modelPtr_->whatsChanged_ = 0;
3428  // Get rid of integer information (modelPtr will get rid of its copy)
3429  // assert( rowsen != NULL );
3430  // assert( rowrhs != NULL );
3431  // If any of Rhs NULLs then create arrays
3432  int numrows = matrix.getNumRows();
3433  const char * rowsenUse = rowsen;
3434  if (!rowsen) {
3435    char * rowsen = new char [numrows];
3436    for (int i=0;i<numrows;i++)
3437      rowsen[i]='G';
3438    rowsenUse = rowsen;
3439  } 
3440  const double * rowrhsUse = rowrhs;
3441  if (!rowrhs) {
3442    double * rowrhs = new double [numrows];
3443    for (int i=0;i<numrows;i++)
3444      rowrhs[i]=0.0;
3445    rowrhsUse = rowrhs;
3446  }
3447  const double * rowrngUse = rowrng;
3448  if (!rowrng) {
3449    double * rowrng = new double [numrows];
3450    for (int i=0;i<numrows;i++)
3451      rowrng[i]=0.0;
3452    rowrngUse = rowrng;
3453  }
3454  double * rowlb = new double[numrows];
3455  double * rowub = new double[numrows];
3456  for (int i = numrows-1; i >= 0; --i) {   
3457    convertSenseToBound(rowsenUse[i],rowrhsUse[i],rowrngUse[i],rowlb[i],rowub[i]);
3458  }
3459  if (rowsen!=rowsenUse)
3460    delete [] rowsenUse;
3461  if (rowrhs!=rowrhsUse)
3462    delete [] rowrhsUse;
3463  if (rowrng!=rowrngUse)
3464    delete [] rowrngUse;
3465  loadProblem(matrix, collb, colub, obj, rowlb, rowub);
3466  delete [] rowlb;
3467  delete [] rowub;
3468}
3469
3470//-----------------------------------------------------------------------------
3471
3472void
3473OsiClpSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
3474                                     double*& collb, double*& colub,
3475                                     double*& obj,
3476                                     char*& rowsen, double*& rowrhs,
3477                                     double*& rowrng)
3478{
3479  modelPtr_->whatsChanged_ = 0;
3480  // Get rid of integer information (modelPtr will get rid of its copy)
3481  loadProblem(*matrix, collb, colub, obj, rowsen, rowrhs, rowrng);
3482  delete matrix;   matrix = NULL;
3483  delete[] collb;  collb = NULL;
3484  delete[] colub;  colub = NULL;
3485  delete[] obj;    obj = NULL;
3486  delete[] rowsen; rowsen = NULL;
3487  delete[] rowrhs; rowrhs = NULL;
3488  delete[] rowrng; rowrng = NULL;
3489}
3490
3491//-----------------------------------------------------------------------------
3492
3493void
3494OsiClpSolverInterface::loadProblem(const int numcols, const int numrows,
3495                                   const CoinBigIndex * start, const int* index,
3496                                   const double* value,
3497                                   const double* collb, const double* colub,
3498                                   const double* obj,
3499                                   const double* rowlb, const double* rowub)
3500{
3501  modelPtr_->whatsChanged_ = 0;
3502  // Get rid of integer information (modelPtr will get rid of its copy)
3503  delete [] integerInformation_;
3504  integerInformation_=NULL;
3505  modelPtr_->loadProblem(numcols, numrows, start,  index,
3506            value, collb, colub, obj,
3507            rowlb,  rowub);
3508  linearObjective_ = modelPtr_->objective();
3509  freeCachedResults();
3510  basis_=CoinWarmStartBasis();
3511  if (ws_) {
3512     delete ws_;
3513     ws_ = 0;
3514  }
3515}
3516//-----------------------------------------------------------------------------
3517
3518void
3519OsiClpSolverInterface::loadProblem(const int numcols, const int numrows,
3520                                   const CoinBigIndex * start, const int* index,
3521                                   const double* value,
3522                                   const double* collb, const double* colub,
3523                                   const double* obj,
3524                                   const char* rowsen, const double* rowrhs,   
3525                                   const double* rowrng)
3526{
3527  modelPtr_->whatsChanged_ = 0;
3528  // Get rid of integer information (modelPtr will get rid of its copy)
3529  // If any of Rhs NULLs then create arrays
3530  const char * rowsenUse = rowsen;
3531  if (!rowsen) {
3532    char * rowsen = new char [numrows];
3533    for (int i=0;i<numrows;i++)
3534      rowsen[i]='G';
3535    rowsenUse = rowsen;
3536  } 
3537  const double * rowrhsUse = rowrhs;
3538  if (!rowrhs) {
3539    double * rowrhs = new double [numrows];
3540    for (int i=0;i<numrows;i++)
3541      rowrhs[i]=0.0;
3542    rowrhsUse = rowrhs;
3543  }
3544  const double * rowrngUse = rowrng;
3545  if (!rowrng) {
3546    double * rowrng = new double [numrows];
3547    for (int i=0;i<numrows;i++)
3548      rowrng[i]=0.0;
3549    rowrngUse = rowrng;
3550  }
3551  double * rowlb = new double[numrows];
3552  double * rowub = new double[numrows];
3553  for (int i = numrows-1; i >= 0; --i) {   
3554    convertSenseToBound(rowsenUse[i],rowrhsUse[i],rowrngUse[i],rowlb[i],rowub[i]);
3555  }
3556  if (rowsen!=rowsenUse)
3557    delete [] rowsenUse;
3558  if (rowrhs!=rowrhsUse)
3559    delete [] rowrhsUse;
3560  if (rowrng!=rowrngUse)
3561    delete [] rowrngUse;
3562  loadProblem(numcols, numrows, start,  index, value, collb, colub, obj,
3563              rowlb,  rowub);
3564  delete[] rowlb;
3565  delete[] rowub;
3566}
3567// This loads a model from a coinModel object - returns number of errors
3568int 
3569OsiClpSolverInterface::loadFromCoinModel (  CoinModel & modelObject, bool keepSolution)
3570{
3571  modelPtr_->whatsChanged_ = 0;
3572  int numberErrors = 0;
3573  // Set arrays for normal use
3574  double * rowLower = modelObject.rowLowerArray();
3575  double * rowUpper = modelObject.rowUpperArray();
3576  double * columnLower = modelObject.columnLowerArray();
3577  double * columnUpper = modelObject.columnUpperArray();
3578  double * objective = modelObject.objectiveArray();
3579  int * integerType = modelObject.integerTypeArray();
3580  double * associated = modelObject.associatedArray();
3581  // If strings then do copies
3582  if (modelObject.stringsExist()) {
3583    numberErrors = modelObject.createArrays(rowLower, rowUpper, columnLower, columnUpper,
3584                                            objective, integerType,associated);
3585  }
3586  CoinPackedMatrix matrix;
3587  modelObject.createPackedMatrix(matrix,associated);
3588  int numberRows = modelObject.numberRows();
3589  int numberColumns = modelObject.numberColumns();
3590  CoinWarmStart * ws = getWarmStart();
3591  bool restoreBasis = keepSolution && numberRows&&numberRows==getNumRows()&&
3592    numberColumns==getNumCols();
3593  loadProblem(matrix, 
3594              columnLower, columnUpper, objective, rowLower, rowUpper);
3595  if (restoreBasis)
3596    setWarmStart(ws);
3597  delete ws;
3598  // Do names if wanted
3599  int numberItems;
3600  numberItems = modelObject.rowNames()->numberItems();
3601  if (numberItems) {
3602    const char *const * rowNames=modelObject.rowNames()->names();
3603    modelPtr_->copyRowNames(rowNames,0,numberItems);
3604  }
3605  numberItems = modelObject.columnNames()->numberItems();
3606  if (numberItems) {
3607    const char *const * columnNames=modelObject.columnNames()->names();
3608    modelPtr_->copyColumnNames(columnNames,0,numberItems);
3609  }
3610  // Do integers if wanted
3611  assert(integerType);
3612  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3613    if (integerType[iColumn])
3614      setInteger(iColumn);
3615  }
3616  if (rowLower!=modelObject.rowLowerArray()||
3617      columnLower!=modelObject.columnLowerArray()) {
3618    delete [] rowLower;
3619    delete [] rowUpper;
3620    delete [] columnLower;
3621    delete [] columnUpper;
3622    delete [] objective;
3623    delete [] integerType;
3624    delete [] associated;
3625    //if (numberErrors)
3626    //  handler_->message(CLP_BAD_STRING_VALUES,messages_)
3627    //    <<numberErrors
3628    //    <<CoinMessageEol;
3629  }
3630  modelPtr_->optimizationDirection_ = modelObject.optimizationDirection(); 
3631  return numberErrors;
3632}
3633
3634//-----------------------------------------------------------------------------
3635// Write mps files
3636//-----------------------------------------------------------------------------
3637
3638void OsiClpSolverInterface::writeMps(const char * filename,
3639                                     const char * extension,
3640                                     double objSense) const
3641{
3642  std::string f(filename);
3643  std::string e(extension);
3644  std::string fullname;
3645  if (e!="") {
3646    fullname = f + "." + e;
3647  } else {
3648    // no extension so no trailing period
3649    fullname = f;
3650  }
3651  // get names
3652  const char * const * const rowNames = modelPtr_->rowNamesAsChar();
3653  const char * const * const columnNames = modelPtr_->columnNamesAsChar();
3654  // Fall back on Osi version - possibly with names
3655  OsiSolverInterface::writeMpsNative(fullname.c_str(), 
3656                                     const_cast<const char **>(rowNames),
3657                                     const_cast<const char **>(columnNames),0,2,objSense,
3658                                     numberSOS_,setInfo_);
3659  if (rowNames) {
3660    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_+1);
3661    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
3662  }
3663}
3664
3665int 
3666OsiClpSolverInterface::writeMpsNative(const char *filename, 
3667                  const char ** rowNames, const char ** columnNames,
3668                  int formatType,int numberAcross,double objSense) const 
3669{
3670  return OsiSolverInterface::writeMpsNative(filename, rowNames, columnNames,
3671                               formatType, numberAcross,objSense,
3672                                            numberSOS_,setInfo_);
3673}
3674
3675//#############################################################################
3676// CLP specific public interfaces
3677//#############################################################################
3678
3679ClpSimplex * OsiClpSolverInterface::getModelPtr() const
3680{
3681  int saveAlgorithm = lastAlgorithm_;
3682  //freeCachedResults();
3683  lastAlgorithm_ = saveAlgorithm;
3684  //bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
3685  return modelPtr_;
3686}
3687
3688//-------------------------------------------------------------------
3689
3690//#############################################################################
3691// Constructors, destructors clone and assignment
3692//#############################################################################
3693//-------------------------------------------------------------------
3694// Default Constructor
3695//-------------------------------------------------------------------
3696OsiClpSolverInterface::OsiClpSolverInterface ()
3697:
3698OsiSolverInterface(),
3699linearObjective_(NULL),
3700rowsense_(NULL),
3701rhs_(NULL),
3702rowrange_(NULL),
3703ws_(NULL),
3704rowActivity_(NULL),
3705columnActivity_(NULL),
3706numberSOS_(0),
3707setInfo_(NULL),
3708smallModel_(NULL),
3709factorization_(NULL),
3710smallestElementInCut_(1.0e-15),
3711smallestChangeInCut_(1.0e-10),
3712largestAway_(-1.0),
3713spareArrays_(NULL),
3714matrixByRow_(NULL),
3715matrixByRowAtContinuous_(NULL), 
3716integerInformation_(NULL),
3717whichRange_(NULL),
3718cleanupScaling_(0),
3719specialOptions_(0x80000000),
3720baseModel_(NULL),
3721lastNumberRows_(0),
3722continuousModel_(NULL),
3723fakeObjective_(NULL)
3724{
3725  //printf("in default %x\n",this);
3726  modelPtr_=NULL;
3727  notOwned_=false;
3728  disasterHandler_ = new OsiClpDisasterHandler();
3729  reset();
3730}
3731
3732//-------------------------------------------------------------------
3733// Clone
3734//-------------------------------------------------------------------
3735OsiSolverInterface * OsiClpSolverInterface::clone(bool CopyData) const
3736{
3737  //printf("in clone %x\n",this);
3738  OsiClpSolverInterface * newSolver;
3739  if (CopyData) {
3740    newSolver = new OsiClpSolverInterface(*this);
3741  } else {
3742    newSolver = new OsiClpSolverInterface();
3743  }
3744#if 0
3745  const double * obj = newSolver->getObjCoefficients();
3746  const double * oldObj = getObjCoefficients();
3747  if(newSolver->getNumCols()>3787)
3748    printf("%x - obj %x (from %x) val %g\n",newSolver,obj,oldObj,obj[3787]);
3749#endif
3750  return newSolver;
3751}
3752
3753
3754//-------------------------------------------------------------------
3755// Copy constructor
3756//-------------------------------------------------------------------
3757OsiClpSolverInterface::OsiClpSolverInterface (
3758                  const OsiClpSolverInterface & rhs)
3759: OsiSolverInterface(rhs),
3760rowsense_(NULL),
3761rhs_(NULL),
3762rowrange_(NULL),
3763ws_(NULL),
3764rowActivity_(NULL),
3765columnActivity_(NULL),
3766  stuff_(rhs.stuff_),
3767numberSOS_(rhs.numberSOS_),
3768setInfo_(NULL),
3769smallModel_(NULL),
3770factorization_(NULL),
3771smallestElementInCut_(rhs.smallestElementInCut_),
3772smallestChangeInCut_(rhs.smallestChangeInCut_),
3773largestAway_(-1.0),
3774spareArrays_(NULL),
3775basis_(),
3776itlimOrig_(9999999),
3777lastAlgorithm_(0),
3778notOwned_(false),
3779matrixByRow_(NULL),
3780matrixByRowAtContinuous_(NULL), 
3781integerInformation_(NULL),
3782whichRange_(NULL)
3783{
3784  //printf("in copy %x - > %x\n",&rhs,this);
3785  if ( rhs.modelPtr_  ) 
3786    modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
3787  else
3788    modelPtr_ = new ClpSimplex();
3789  if ( rhs.baseModel_  ) 
3790    baseModel_ = new ClpSimplex(*rhs.baseModel_);
3791  else
3792    baseModel_ = NULL;
3793  if ( rhs.continuousModel_  ) 
3794    continuousModel_ = new ClpSimplex(*rhs.continuousModel_);
3795  else
3796    continuousModel_ = NULL;
3797  if (rhs.matrixByRowAtContinuous_)
3798    matrixByRowAtContinuous_ = new CoinPackedMatrix(*rhs.matrixByRowAtContinuous_);
3799  if ( rhs.disasterHandler_  ) 
3800    disasterHandler_ = dynamic_cast<OsiClpDisasterHandler *>(rhs.disasterHandler_->clone());
3801  else
3802    disasterHandler_ = NULL;
3803  if (rhs.fakeObjective_)
3804    fakeObjective_ = new ClpLinearObjective(*rhs.fakeObjective_);
3805  else
3806    fakeObjective_ = NULL;
3807  linearObjective_ = modelPtr_->objective();
3808  if ( rhs.ws_ ) 
3809    ws_ = new CoinWarmStartBasis(*rhs.ws_);
3810  basis_ = rhs.basis_;
3811  if (rhs.integerInformation_) {
3812    int numberColumns = modelPtr_->numberColumns();
3813    integerInformation_ = new char[numberColumns];
3814    CoinMemcpyN(rhs.integerInformation_,        numberColumns,integerInformation_);
3815  }
3816  saveData_ = rhs.saveData_;
3817  solveOptions_  = rhs.solveOptions_;
3818  cleanupScaling_ = rhs.cleanupScaling_;
3819  specialOptions_ = rhs.specialOptions_;
3820  lastNumberRows_ = rhs.lastNumberRows_;
3821  rowScale_ = rhs.rowScale_;
3822  columnScale_ = rhs.columnScale_;
3823  fillParamMaps();
3824  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3825  if (numberSOS_) {
3826    setInfo_ = new CoinSet[numberSOS_];
3827    for (int i=0;i<numberSOS_;i++)
3828      setInfo_[i]=rhs.setInfo_[i];
3829  }
3830}
3831
3832// Borrow constructor - only delete one copy
3833OsiClpSolverInterface::OsiClpSolverInterface (ClpSimplex * rhs,
3834                                              bool reallyOwn)
3835:
3836OsiSolverInterface(),
3837rowsense_(NULL),
3838rhs_(NULL),
3839rowrange_(NULL),
3840ws_(NULL),
3841rowActivity_(NULL),
3842columnActivity_(NULL),
3843numberSOS_(0),
3844setInfo_(NULL),
3845smallModel_(NULL),
3846factorization_(NULL),
3847smallestElementInCut_(1.0e-15),
3848smallestChangeInCut_(1.0e-10),
3849largestAway_(-1.0),
3850spareArrays_(NULL),
3851basis_(),
3852itlimOrig_(9999999),
3853lastAlgorithm_(0),
3854notOwned_(false),
3855matrixByRow_(NULL),
3856matrixByRowAtContinuous_(NULL), 
3857integerInformation_(NULL),
3858whichRange_(NULL),
3859cleanupScaling_(0),
3860specialOptions_(0x80000000),
3861baseModel_(NULL),
3862lastNumberRows_(0),
3863continuousModel_(NULL),
3864fakeObjective_(NULL)
3865{
3866  disasterHandler_ = new OsiClpDisasterHandler();
3867  //printf("in borrow %x - > %x\n",&rhs,this);
3868  modelPtr_ = rhs;
3869  basis_.resize(modelPtr_->numberRows(),modelPtr_->numberColumns());
3870  linearObjective_ = modelPtr_->objective();
3871  if (rhs) {
3872    notOwned_=!reallyOwn;
3873
3874    if (rhs->integerInformation()) {
3875      int numberColumns = modelPtr_->numberColumns();
3876      integerInformation_ = new char[numberColumns];
3877      CoinMemcpyN(rhs->integerInformation(),    numberColumns,integerInformation_);
3878    }
3879  }
3880  fillParamMaps();
3881}
3882   
3883// Releases so won't error
3884void 
3885OsiClpSolverInterface::releaseClp()
3886{
3887  modelPtr_=NULL;
3888  notOwned_=false;
3889}
3890   
3891//-------------------------------------------------------------------
3892// Destructor
3893//-------------------------------------------------------------------
3894OsiClpSolverInterface::~OsiClpSolverInterface ()
3895{
3896  //printf("in destructor %x\n",this);
3897  freeCachedResults();
3898  if (!notOwned_)
3899    delete modelPtr_;
3900  delete baseModel_;
3901  delete continuousModel_;
3902  delete disasterHandler_;
3903  delete fakeObjective_;
3904  delete ws_;
3905  delete [] rowActivity_;
3906  delete [] columnActivity_;
3907  delete [] setInfo_;
3908#ifdef KEEP_SMALL
3909  if (smallModel_) {
3910    delete [] spareArrays_;
3911    spareArrays_ = NULL;
3912    delete smallModel_;
3913    smallModel_=NULL;
3914  }
3915#endif
3916  assert(smallModel_==NULL);
3917  assert(factorization_==NULL);
3918  assert(spareArrays_==NULL);
3919  delete [] integerInformation_;
3920  delete matrixByRowAtContinuous_;
3921  delete matrixByRow_;
3922}
3923
3924//-------------------------------------------------------------------
3925// Assignment operator
3926//-------------------------------------------------------------------
3927OsiClpSolverInterface &
3928OsiClpSolverInterface::operator=(const OsiClpSolverInterface& rhs)
3929{
3930  if (this != &rhs) {   
3931    //printf("in = %x - > %x\n",&rhs,this);
3932    OsiSolverInterface::operator=(rhs);
3933    freeCachedResults();
3934    if (!notOwned_)
3935      delete modelPtr_;
3936    delete ws_;
3937    if ( rhs.modelPtr_  ) 
3938      modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
3939    delete baseModel_;
3940    if ( rhs.baseModel_  ) 
3941      baseModel_ = new ClpSimplex(*rhs.baseModel_);
3942    else
3943      baseModel_ = NULL;
3944    delete continuousModel_;
3945    if ( rhs.continuousModel_  ) 
3946      continuousModel_ = new ClpSimplex(*rhs.continuousModel_);
3947    else
3948      continuousModel_ = NULL;
3949    delete matrixByRowAtContinuous_;
3950    delete matrixByRow_;
3951    matrixByRow_=NULL;
3952    if (rhs.matrixByRowAtContinuous_)
3953      matrixByRowAtContinuous_ = new CoinPackedMatrix(*rhs.matrixByRowAtContinuous_);
3954    else
3955      matrixByRowAtContinuous_=NULL;
3956    delete disasterHandler_;
3957    if ( rhs.disasterHandler_  ) 
3958      disasterHandler_ = dynamic_cast<OsiClpDisasterHandler *>(rhs.disasterHandler_->clone());
3959    else
3960      disasterHandler_ = NULL;
3961    delete fakeObjective_;
3962    if (rhs.fakeObjective_)
3963      fakeObjective_ = new ClpLinearObjective(*rhs.fakeObjective_);
3964    else
3965      fakeObjective_ = NULL;
3966    notOwned_=false;
3967    linearObjective_ = modelPtr_->objective();
3968    saveData_ = rhs.saveData_;
3969    solveOptions_  = rhs.solveOptions_;
3970    cleanupScaling_ = rhs.cleanupScaling_;
3971    specialOptions_ = rhs.specialOptions_;
3972    lastNumberRows_ = rhs.lastNumberRows_;
3973    rowScale_ = rhs.rowScale_;
3974    columnScale_ = rhs.columnScale_;
3975    basis_ = rhs.basis_;
3976    stuff_ = rhs.stuff_;
3977    if (rhs.integerInformation_) {
3978      int numberColumns = modelPtr_->numberColumns();
3979      integerInformation_ = new char[numberColumns];
3980      CoinMemcpyN(rhs.integerInformation_,      numberColumns,integerInformation_);
3981    }
3982    if ( rhs.ws_ ) 
3983      ws_ = new CoinWarmStartBasis(*rhs.ws_);
3984    else
3985      ws_=NULL;
3986    delete [] rowActivity_;
3987    delete [] columnActivity_;
3988    rowActivity_=NULL;
3989    columnActivity_=NULL;
3990    delete [] setInfo_;
3991    numberSOS_ = rhs.numberSOS_;
3992    setInfo_=NULL;
3993    if (numberSOS_) {
3994      setInfo_ = new CoinSet[numberSOS_];
3995      for (int i=0;i<numberSOS_;i++)
3996        setInfo_[i]=rhs.setInfo_[i];
3997    }
3998    assert(smallModel_==NULL);
3999    assert(factorization_==NULL);
4000    smallestElementInCut_ = rhs.smallestElementInCut_;
4001    smallestChangeInCut_ = rhs.smallestChangeInCut_;
4002    largestAway_ = -1.0;
4003    assert(spareArrays_==NULL);
4004    basis_ = rhs.basis_;
4005    fillParamMaps();
4006    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4007  }
4008  return *this;
4009}
4010
4011//#############################################################################
4012// Applying cuts
4013//#############################################################################
4014
4015void OsiClpSolverInterface::applyRowCut( const OsiRowCut & rowCut )
4016{
4017  applyRowCuts(1, &rowCut);
4018}
4019/* Apply a collection of row cuts which are all effective.
4020   applyCuts seems to do one at a time which seems inefficient.
4021*/
4022void 
4023OsiClpSolverInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
4024{
4025  if (numberCuts) {
4026    // Say can't gurantee optimal basis etc
4027    lastAlgorithm_=999;
4028
4029    // Thanks to js
4030    const OsiRowCut * * cutsp = new const OsiRowCut * [numberCuts];
4031    for (int i=0;i<numberCuts;i++) 
4032      cutsp[i] = &cuts[i];
4033   
4034    applyRowCuts(numberCuts, cutsp);
4035   
4036    delete [] cutsp;
4037  }
4038}
4039/* Apply a collection of row cuts which are all effective.
4040   applyCuts seems to do one at a time which seems inefficient.
4041*/
4042void 
4043OsiClpSolverInterface::applyRowCuts(int numberCuts, const OsiRowCut ** cuts)
4044{
4045  int i;
4046  if (!numberCuts)
4047    return;
4048  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
4049  CoinPackedMatrix * saveRowCopy = matrixByRow_;
4050  matrixByRow_=NULL;
4051#if 0 // was #ifndef NDEBUG
4052  int nameDiscipline;
4053  getIntParam(OsiNameDiscipline,nameDiscipline) ;
4054  assert (!nameDiscipline);
4055#endif
4056  freeCachedResults0();
4057  // Say can't gurantee optimal basis etc
4058  lastAlgorithm_=999;
4059  int numberRows = modelPtr_->numberRows();
4060  modelPtr_->resize(numberRows+numberCuts,modelPtr_->numberColumns());
4061  basis_.resize(numberRows+numberCuts,modelPtr_->numberColumns());
4062  // redo as relaxed - use modelPtr_-> addRows with starts etc
4063  int size = 0;
4064  for (i=0;i<numberCuts;i++) 
4065    size += cuts[i]->row().getNumElements();
4066  CoinBigIndex * starts = new CoinBigIndex [numberCuts+1];
4067  int * indices = new int[size];
4068  double * elements = new double[size];
4069  double * lower = modelPtr_->rowLower()+numberRows;
4070  double * upper = modelPtr_->rowUpper()+numberRows;
4071  const double * columnLower = modelPtr_->columnLower();
4072  const double * columnUpper = modelPtr_->columnUpper();
4073  size=0;
4074  for (i=0;i<numberCuts;i++) {
4075    double rowLb = cuts[i]->lb();
4076    double rowUb = cuts[i]->ub();
4077    int n=cuts[i]->row().getNumElements();
4078    const int * index = cuts[i]->row().getIndices();
4079    const double * elem = cuts[i]->row().getElements();
4080    starts[i]=size;
4081    for (int j=0;j<n;j++) {
4082      double value = elem[j];
4083      int column = index[j];
4084      if (fabs(value)>=smallestChangeInCut_) {
4085        // always take
4086        indices[size]=column;
4087        elements[size++]=value;
4088      } else if (fabs(value)>=smallestElementInCut_) {
4089        double lowerValue = columnLower[column];
4090        double upperValue = columnUpper[column];
4091        double difference = upperValue-lowerValue;
4092        if (difference<1.0e20&&difference*fabs(value)<smallestChangeInCut_&&
4093            (rowLb<-1.0e20||rowUb>1.0e20)) {
4094          // Take out and adjust to relax
4095          //printf("small el %g adjusted\n",value);
4096          if (rowLb>-1.0e20) {
4097            // just lower bound on row
4098            if (value>0.0) {
4099              // pretend at upper
4100              rowLb -= value*upperValue;
4101            } else {
4102              // pretend at lower
4103              rowLb -= value*lowerValue;
4104            }
4105          } else {
4106            // just upper bound on row
4107            if (value>0.0) {
4108              // pretend at lower
4109              rowUb -= value*lowerValue;
4110            } else {
4111              // pretend at upper
4112              rowUb -= value*upperValue;
4113            }
4114          }
4115        } else {
4116          // take (unwillingly)
4117          indices[size]=column;
4118          elements[size++]=value;
4119        }
4120      } else {
4121        //printf("small el %g ignored\n",value);
4122      }
4123    }
4124    lower[i]= forceIntoRange(rowLb, -OsiClpInfinity, OsiClpInfinity);
4125    upper[i]= forceIntoRange(rowUb, -OsiClpInfinity, OsiClpInfinity);
4126    if (lower[i]<-1.0e27)
4127      lower[i]=-COIN_DBL_MAX;
4128    if (upper[i]>1.0e27)
4129      upper[i]=COIN_DBL_MAX;
4130  }
4131  starts[numberCuts]=size;
4132  if (!modelPtr_->clpMatrix())
4133    modelPtr_->createEmptyMatrix();
4134  //modelPtr_->matrix()->appendRows(numberCuts,rows);
4135  modelPtr_->clpMatrix()->appendMatrix(numberCuts,0,starts,indices,elements);
4136  modelPtr_->setNewRowCopy(NULL);
4137  freeCachedResults1();
4138  redoScaleFactors( numberCuts,starts, indices, elements);
4139  if (saveRowCopy) {
4140#if 1
4141    matrixByRow_=saveRowCopy;
4142    matrixByRow_->appendRows(numberCuts,starts,indices,elements,0);
4143    assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
4144#else
4145    delete saveRowCopy;
4146#endif
4147  }
4148  delete [] starts;
4149  delete [] indices;
4150  delete [] elements;
4151
4152}
4153//#############################################################################
4154// Apply Cuts
4155//#############################################################################
4156
4157OsiSolverInterface::ApplyCutsReturnCode
4158OsiClpSolverInterface::applyCuts( const OsiCuts & cs, double effectivenessLb ) 
4159{
4160  OsiSolverInterface::ApplyCutsReturnCode retVal;
4161  int i;
4162
4163  // Loop once for each column cut
4164  for ( i=0; i<cs.sizeColCuts(); i ++ ) {
4165    if ( cs.colCut(i).effectiveness() < effectivenessLb ) {
4166      retVal.incrementIneffective();
4167      continue;
4168    }
4169    if ( !cs.colCut(i).consistent() ) {
4170      retVal.incrementInternallyInconsistent();
4171      continue;
4172    }
4173    if ( !cs.colCut(i).consistent(*this) ) {
4174      retVal.incrementExternallyInconsistent();
4175      continue;
4176    }
4177    if ( cs.colCut(i).infeasible(*this) ) {
4178      retVal.incrementInfeasible();
4179      continue;
4180    }
4181    applyColCut( cs.colCut(i) );
4182    retVal.incrementApplied();
4183  }
4184
4185  // Loop once for each row cut
4186  const OsiRowCut ** addCuts = new const OsiRowCut* [cs.sizeRowCuts()];
4187  int nAdd=0;
4188  for ( i=0; i<cs.sizeRowCuts(); i ++ ) {
4189    if ( cs.rowCut(i).effectiveness() < effectivenessLb ) {
4190      retVal.incrementIneffective();
4191      continue;
4192    }
4193    if ( !cs.rowCut(i).consistent() ) {
4194      retVal.incrementInternallyInconsistent();
4195      continue;
4196    }
4197    if ( !cs.rowCut(i).consistent(*this) ) {
4198      retVal.incrementExternallyInconsistent();
4199      continue;
4200    }
4201    if ( cs.rowCut(i).infeasible(*this) ) {
4202      retVal.incrementInfeasible();
4203      continue;
4204    }
4205    addCuts[nAdd++] = cs.rowCutPtr(i);
4206    retVal.incrementApplied();
4207  }
4208  // now apply
4209  applyRowCuts(nAdd,addCuts);
4210  delete [] addCuts;
4211 
4212  return retVal;
4213}
4214// Extend scale factors
4215void 
4216OsiClpSolverInterface::redoScaleFactors(int numberAdd,const CoinBigIndex * starts,
4217                                        const int * indices, const double * elements)
4218{
4219  if ((specialOptions_&131072)!=0) {
4220    int numberRows = modelPtr_->numberRows()-numberAdd;
4221    assert (lastNumberRows_==numberRows); // ???
4222    int iRow;
4223    int newNumberRows = numberRows + numberAdd;
4224    rowScale_.extend(static_cast<int>(2*newNumberRows*sizeof(double)));
4225    double * rowScale = rowScale_.array();
4226    double * oldInverseScale = rowScale + lastNumberRows_;
4227    double * inverseRowScale = rowScale + newNumberRows;
4228    for (iRow=lastNumberRows_-1;iRow>=0;iRow--)
4229      inverseRowScale[iRow] = oldInverseScale[iRow] ;
4230    //int numberColumns = baseModel_->numberColumns();
4231    const double * columnScale = columnScale_.array();
4232    //const double * inverseColumnScale = columnScale + numberColumns;
4233    // Geometric mean on row scales
4234    // adjust arrays
4235    rowScale += lastNumberRows_;
4236    inverseRowScale += lastNumberRows_;
4237    for (iRow=0;iRow<numberAdd;iRow++) {
4238      CoinBigIndex j;
4239      double largest=1.0e-20;
4240      double smallest=1.0e50;
4241      for (j=starts[iRow];j<starts[iRow+1];j++) {
4242        int iColumn = indices[j];
4243        double value = fabs(elements[j]);
4244        // Don't bother with tiny elements
4245        if (value>1.0e-20) {
4246          value *= columnScale[iColumn];
4247          largest = CoinMax(largest,value);
4248          smallest = CoinMin(smallest,value);
4249        }
4250      }
4251      double scale=sqrt(smallest*largest);
4252      scale=CoinMax(1.0e-10,CoinMin(1.0e10,scale));
4253      inverseRowScale[iRow]=scale;
4254      rowScale[iRow]=1.0/scale;
4255    }
4256    lastNumberRows_=newNumberRows;
4257  }
4258}
4259// Delete all scale factor stuff and reset option
4260void OsiClpSolverInterface::deleteScaleFactors()
4261{
4262  delete baseModel_;
4263  baseModel_=NULL;
4264  lastNumberRows_=0;
4265  specialOptions_ &= ~131072;
4266}
4267//-----------------------------------------------------------------------------
4268
4269void OsiClpSolverInterface::applyColCut( const OsiColCut & cc )
4270{
4271  modelPtr_->whatsChanged_ &= (0x1ffff&~(128|256));
4272  // Say can't gurantee optimal basis etc
4273  lastAlgorithm_=999;
4274  double * lower = modelPtr_->columnLower();
4275  double * upper = modelPtr_->columnUpper();
4276  const CoinPackedVector & lbs = cc.lbs();
4277  const CoinPackedVector & ubs = cc.ubs();
4278  int i;
4279
4280  for ( i=0; i<lbs.getNumElements(); i++ ) {
4281    int iCol = lbs.getIndices()[i];
4282    double value = lbs.getElements()[i];
4283    if ( value > lower[iCol] )
4284      lower[iCol]= value;
4285  }
4286  for ( i=0; i<ubs.getNumElements(); i++ ) {
4287    int iCol = ubs.getIndices()[i];
4288    double value = ubs.getElements()[i];
4289    if ( value < upper[iCol] )
4290      upper[iCol]= value;
4291  }
4292}
4293//#############################################################################
4294// Private methods
4295//#############################################################################
4296
4297
4298//-------------------------------------------------------------------
4299
4300void OsiClpSolverInterface::freeCachedResults() const
4301{ 
4302  // Say can't gurantee optimal basis etc
4303  lastAlgorithm_=999;
4304  delete [] rowsense_;
4305  delete [] rhs_;
4306  delete [] rowrange_;
4307  delete matrixByRow_;
4308  if (modelPtr_&&modelPtr_->scaledMatrix_) {
4309    delete modelPtr_->scaledMatrix_;
4310    modelPtr_->scaledMatrix_=NULL;
4311  }
4312  //delete ws_;
4313  rowsense_=NULL;
4314  rhs_=NULL;
4315  rowrange_=NULL;
4316  matrixByRow_=NULL;
4317  //ws_ = NULL;
4318  if (modelPtr_&&modelPtr_->clpMatrix()) {
4319    modelPtr_->clpMatrix()->refresh(modelPtr_); // make sure all clean
4320#ifndef NDEBUG
4321    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (modelPtr_->clpMatrix());
4322    if (clpMatrix) {
4323      if (clpMatrix->getNumCols())
4324        assert (clpMatrix->getNumRows()==modelPtr_->getNumRows());
4325      if (clpMatrix->getNumRows())
4326        assert (clpMatrix->getNumCols()==modelPtr_->getNumCols());
4327    }
4328#endif
4329  }
4330}
4331
4332//-------------------------------------------------------------------
4333
4334void OsiClpSolverInterface::freeCachedResults0() const
4335{ 
4336  delete [] rowsense_;
4337  delete [] rhs_;
4338  delete [] rowrange_;
4339  rowsense_=NULL;
4340  rhs_=NULL;
4341  rowrange_=NULL;
4342}
4343
4344//-------------------------------------------------------------------
4345
4346void OsiClpSolverInterface::freeCachedResults1() const
4347{ 
4348  // Say can't gurantee optimal basis etc
4349  lastAlgorithm_=999;
4350  delete matrixByRow_;
4351  matrixByRow_=NULL;
4352  //ws_ = NULL;
4353  if (modelPtr_&&modelPtr_->clpMatrix()) {
4354    delete modelPtr_->scaledMatrix_;
4355    modelPtr_->scaledMatrix_=NULL;
4356    modelPtr_->clpMatrix()->refresh(modelPtr_); // make sure all clean
4357#ifndef NDEBUG
4358    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (modelPtr_->clpMatrix());
4359    if (clpMatrix) {
4360      assert (clpMatrix->getNumRows()==modelPtr_->getNumRows());
4361      assert (clpMatrix->getNumCols()==modelPtr_->getNumCols());
4362    }
4363#endif
4364  }
4365}
4366
4367//------------------------------------------------------------------
4368void OsiClpSolverInterface::extractSenseRhsRange() const
4369{
4370  if (rowsense_ == NULL) {
4371    // all three must be NULL
4372    assert ((rhs_ == NULL) && (rowrange_ == NULL));
4373   
4374    int nr=modelPtr_->numberRows();
4375    if ( nr!=0 ) {
4376      rowsense_ = new char[nr];
4377      rhs_ = new double[nr];
4378      rowrange_ = new double[nr];
4379      std::fill(rowrange_,rowrange_+nr,0.0);
4380     
4381      const double * lb = modelPtr_->rowLower();
4382      const double * ub = modelPtr_->rowUpper();
4383     
4384      int i;
4385      for ( i=0; i<nr; i++ ) {
4386        convertBoundToSense(lb[i], ub[i], rowsense_[i], rhs_[i], rowrange_[i]);
4387      }
4388    }
4389  }
4390}
4391// Set language
4392void 
4393OsiClpSolverInterface::newLanguage(CoinMessages::Language language)
4394{
4395  modelPtr_->newLanguage(language);
4396  OsiSolverInterface::newLanguage(language);
4397}
4398//#############################################################################
4399
4400void
4401OsiClpSolverInterface::fillParamMaps()
4402{
4403  assert (static_cast<int> (OsiMaxNumIteration)==        static_cast<int>(ClpMaxNumIteration));
4404  assert (static_cast<int> (OsiMaxNumIterationHotStart)==static_cast<int>(ClpMaxNumIterationHotStart));
4405  //assert (static_cast<int> (OsiLastIntParam)==           static_cast<int>(ClpLastIntParam));
4406 
4407  assert (static_cast<int> (OsiDualObjectiveLimit)==  static_cast<int>(ClpDualObjectiveLimit));
4408  assert (static_cast<int> (OsiPrimalObjectiveLimit)==static_cast<int>(ClpPrimalObjectiveLimit));
4409  assert (static_cast<int> (OsiDualTolerance)==       static_cast<int>(ClpDualTolerance));
4410  assert (static_cast<int> (OsiPrimalTolerance)==     static_cast<int>(ClpPrimalTolerance));
4411  assert (static_cast<int> (OsiObjOffset)==           static_cast<int>(ClpObjOffset));
4412  //assert (static_cast<int> (OsiLastDblParam)==        static_cast<int>(ClpLastDblParam));
4413 
4414  assert (static_cast<int> (OsiProbName)==    static_cast<int> (ClpProbName));
4415  //strParamMap_[OsiLastStrParam] = ClpLastStrParam;
4416}
4417// Sets up basis
4418void 
4419OsiClpSolverInterface::setBasis ( const CoinWarmStartBasis & basis)
4420{
4421  setBasis(basis,modelPtr_);
4422  setWarmStart(&basis); 
4423}
4424// Warm start
4425CoinWarmStartBasis
4426OsiClpSolverInterface::getBasis(ClpSimplex * model) const
4427{
4428  int iRow,iColumn;
4429  int numberRows = model->numberRows();
4430  int numberColumns = model->numberColumns();
4431  CoinWarmStartBasis basis;
4432  basis.setSize(numberColumns,numberRows);
4433  if (model->statusExists()) {
4434    // Flip slacks
4435    int lookupA[]={0,1,3,2,0,2};
4436    for (iRow=0;iRow<numberRows;iRow++) {
4437      int iStatus = model->getRowStatus(iRow);
4438      iStatus = lookupA[iStatus];
4439      basis.setArtifStatus(iRow,static_cast<CoinWarmStartBasis::Status> (iStatus));
4440    }
4441    int lookupS[]={0,1,2,3,0,3};
4442    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4443      int iStatus = model->getColumnStatus(iColumn);
4444      iStatus = lookupS[iStatus];
4445      basis.setStructStatus(iColumn,static_cast<CoinWarmStartBasis::Status> (iStatus));
4446    }
4447  }
4448  //basis.print();
4449  return basis;
4450}
4451// Warm start from statusArray
4452CoinWarmStartBasis * 
4453OsiClpSolverInterface::getBasis(const unsigned char * statusArray) const 
4454{
4455  int iRow,iColumn;
4456  int numberRows = modelPtr_->numberRows();
4457  int numberColumns = modelPtr_->numberColumns();
4458  CoinWarmStartBasis * basis = new CoinWarmStartBasis();
4459  basis->setSize(numberColumns,numberRows);
4460  // Flip slacks
4461  int lookupA[]={0,1,3,2,0,2};
4462  for (iRow=0;iRow<numberRows;iRow++) {
4463    int iStatus = statusArray[numberColumns+iRow]&7;
4464    iStatus = lookupA[iStatus];
4465    basis->setArtifStatus(iRow,static_cast<CoinWarmStartBasis::Status> (iStatus));
4466  }
4467  int lookupS[]={0,1,2,3,0,3};
4468  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4469    int iStatus = statusArray[iColumn]&7;
4470    iStatus = lookupS[iStatus];
4471    basis->setStructStatus(iColumn,static_cast<CoinWarmStartBasis::Status> (iStatus));
4472  }
4473  //basis->print();
4474  return basis;
4475}
4476// Sets up basis
4477void 
4478OsiClpSolverInterface::setBasis ( const CoinWarmStartBasis & basis,
4479                                  ClpSimplex * model)
4480{
4481  // Say can't gurantee optimal basis etc
4482  lastAlgorithm_=999;
4483  // transform basis to status arrays
4484  int iRow,iColumn;
4485  int numberRows = model->numberRows();
4486  int numberColumns = model->numberColumns();
4487  if (!model->statusExists()) {
4488    /*
4489      get status arrays
4490      ClpBasis would seem to have overheads and we will need
4491      extra bits anyway.
4492    */
4493    model->createStatus();
4494  }
4495  if (basis.getNumArtificial()!=numberRows||
4496      basis.getNumStructural()!=numberColumns) {
4497    CoinWarmStartBasis basis2 = basis;
4498    // resize
4499    basis2.resize(numberRows,numberColumns);
4500    // move status
4501    model->createStatus();
4502    // For rows lower and upper are flipped
4503    for (iRow=0;iRow<numberRows;iRow++) {
4504      int stat = basis2.getArtifStatus(iRow);
4505      if (stat>1)
4506        stat = 5 - stat; // so 2->3 and 3->2
4507      model->setRowStatus(iRow, static_cast<ClpSimplex::Status> (stat));
4508    }
4509    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4510      model->setColumnStatus(iColumn,
4511                             static_cast<ClpSimplex::Status> (basis2.getStructStatus(iColumn)));
4512    }
4513  } else {
4514    // move status
4515    model->createStatus();
4516    // For rows lower and upper are flipped
4517    for (iRow=0;iRow<numberRows;iRow++) {
4518      int stat = basis.getArtifStatus(iRow);
4519      if (stat>1)
4520        stat = 5 - stat; // so 2->3 and 3->2
4521      model->setRowStatus(iRow, static_cast<ClpSimplex::Status> (stat));
4522    }
4523    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4524      model->setColumnStatus(iColumn,
4525                             static_cast<ClpSimplex::Status> (basis.getStructStatus(iColumn)));
4526    }
4527  }
4528}
4529// Warm start difference from basis_ to statusArray
4530CoinWarmStartDiff * 
4531OsiClpSolverInterface::getBasisDiff(const unsigned char * statusArray) const
4532{
4533  int iRow,iColumn;
4534  int numberRows = modelPtr_->numberRows();
4535  int numberColumns = modelPtr_->numberColumns();
4536  CoinWarmStartBasis basis;
4537  basis.setSize(numberColumns,numberRows);
4538  assert (modelPtr_->statusExists());
4539  int lookupS[]={0,1,2,3,0,3};
4540  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4541    int iStatus = statusArray[iColumn]&7;
4542    iStatus = lookupS[iStatus];
4543    basis.setStructStatus(iColumn,static_cast<CoinWarmStartBasis::Status> (iStatus));
4544  }
4545  statusArray += numberColumns;
4546  // Flip slacks
4547  int lookupA[]={0,1,3,2,0,2};
4548  for (iRow=0;iRow<numberRows;iRow++) {
4549    int iStatus = statusArray[iRow]&7;
4550    iStatus = lookupA[iStatus];
4551    basis.setArtifStatus(iRow,static_cast<CoinWarmStartBasis::Status> (iStatus));
4552  }
4553  // Now basis is what we want while basis_ is old
4554  CoinWarmStartDiff * difference = basis.generateDiff(&basis_);
4555  return difference;
4556}
4557/* Read an mps file from the given filename (defaults to Osi reader) - returns
4558   number of errors (see OsiMpsReader class) */
4559int 
4560OsiClpSolverInterface::readMps(const char *filename,
4561                               const char *extension ) 
4562{
4563  // Get rid of integer stuff
4564  delete [] integerInformation_;
4565  integerInformation_=NULL;
4566  freeCachedResults();
4567 
4568  CoinMpsIO m;
4569  m.setInfinity(getInfinity());
4570  m.passInMessageHandler(modelPtr_->messageHandler());
4571  *m.messagesPointer()=modelPtr_->coinMessages();
4572
4573  delete [] setInfo_;
4574  setInfo_=NULL;
4575  numberSOS_=0;
4576  CoinSet ** sets=NULL;
4577  int numberErrors = m.readMps(filename,extension,numberSOS_,sets);
4578  if (numberSOS_) {
4579    setInfo_ = new CoinSet[numberSOS_];
4580    for (int i=0;i<numberSOS_;i++) {
4581      setInfo_[i]=*sets[i];
4582      delete sets[i];
4583    }
4584    delete [] sets;
4585  }
4586  handler_->message(COIN_SOLVER_MPS,messages_)
4587    <<m.getProblemName()<< numberErrors <<CoinMessageEol;
4588  if (!numberErrors) {
4589
4590    // set objective function offest
4591    setDblParam(OsiObjOffset,m.objectiveOffset());
4592
4593    // set problem name
4594    setStrParam(OsiProbName,m.getProblemName());
4595   
4596    // no errors
4597    loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
4598                m.getObjCoefficients(),m.getRowSense(),m.getRightHandSide(),
4599                m.getRowRange());
4600    const char * integer = m.integerColumns();
4601    int nCols=m.getNumCols();
4602    int nRows=m.getNumRows();
4603    if (integer) {
4604      int i,n=0;
4605      int * index = new int [nCols];
4606      for (i=0;i<nCols;i++) {
4607        if (integer[i]) {
4608          index[n++]=i;
4609        }
4610      }
4611      setInteger(index,n);
4612      delete [] index;
4613      if (n) 
4614        modelPtr_->copyInIntegerInformation(integer);
4615    }
4616
4617    // set objective name
4618    setObjName(m.getObjectiveName());
4619
4620    // Always keep names
4621    int nameDiscipline;
4622    getIntParam(OsiNameDiscipline,nameDiscipline) ;
4623    int iRow;
4624    std::vector<std::string> rowNames = std::vector<std::string> ();
4625    std::vector<std::string> columnNames = std::vector<std::string> ();
4626    rowNames.reserve(nRows);
4627    for (iRow=0;iRow<nRows;iRow++) {
4628      const char * name = m.rowName(iRow);
4629      rowNames.push_back(name);
4630      if (nameDiscipline) 
4631        OsiSolverInterface::setRowName(iRow,name) ;
4632    }
4633   
4634    int iColumn;
4635    columnNames.reserve(nCols);
4636    for (iColumn=0;iColumn<nCols;iColumn++) {
4637      const char * name = m.columnName(iColumn);
4638      columnNames.push_back(name);
4639      if (nameDiscipline) 
4640        OsiSolverInterface::setColName(iColumn,name) ;
4641    }
4642    modelPtr_->copyNames(rowNames,columnNames);
4643  }
4644  return numberErrors;
4645}
4646int 
4647OsiClpSolverInterface::readMps(const char *filename, const char*extension,
4648                            int & numberSets, CoinSet ** & sets)
4649{
4650  int numberErrors = readMps(filename,extension);
4651  numberSets= numberSOS_;
4652  sets = &setInfo_;
4653  return numberErrors;
4654}
4655/* Read an mps file from the given filename returns
4656   number of errors (see OsiMpsReader class) */
4657int 
4658OsiClpSolverInterface::readMps(const char *filename,bool keepNames,bool allowErrors)
4659{
4660  // Get rid of integer stuff
4661  delete [] integerInformation_;
4662  integerInformation_=NULL;
4663  freeCachedResults();
4664 
4665  CoinMpsIO m;
4666  m.setInfinity(getInfinity());
4667  m.passInMessageHandler(modelPtr_->messageHandler());
4668  *m.messagesPointer()=modelPtr_->coinMessages();
4669
4670  delete [] setInfo_;
4671  setInfo_=NULL;
4672  numberSOS_=0;
4673  CoinSet ** sets=NULL;
4674  int numberErrors = m.readMps(filename,"",numberSOS_,sets);
4675  if (numberSOS_) {
4676    setInfo_ = new CoinSet[numberSOS_];
4677    for (int i=0;i<numberSOS_;i++) {
4678      setInfo_[i]=*sets[i];
4679      delete sets[i];
4680    }
4681    delete [] sets;
4682  }
4683  handler_->message(COIN_SOLVER_MPS,messages_)
4684    <<m.getProblemName()<< numberErrors <<CoinMessageEol;
4685  if (!numberErrors||((numberErrors>0&&numberErrors<100000)&&allowErrors)) {
4686
4687    // set objective function offest
4688    setDblParam(OsiObjOffset,m.objectiveOffset());
4689
4690    // set problem name
4691    setStrParam(OsiProbName,m.getProblemName());
4692
4693    // set objective name
4694    setObjName(m.getObjectiveName());
4695
4696    // no errors
4697    loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
4698                m.getObjCoefficients(),m.getRowSense(),m.getRightHandSide(),
4699                m.getRowRange());
4700    int nCols=m.getNumCols();
4701    // get quadratic part
4702    if (m.reader()->whichSection (  ) == COIN_QUAD_SECTION ) {
4703      int * start=NULL;
4704      int * column = NULL;
4705      double * element = NULL;
4706      int status=m.readQuadraticMps(NULL,start,column,element,2);
4707      if (!status) 
4708        modelPtr_->loadQuadraticObjective(nCols,start,column,element);
4709      delete [] start;
4710      delete [] column;
4711      delete [] element;
4712    }
4713    const char * integer = m.integerColumns();
4714    int nRows=m.getNumRows();
4715    if (integer) {
4716      int i,n=0;
4717      int * index = new int [nCols];
4718      for (i=0;i<nCols;i++) {
4719        if (integer[i]) {
4720          index[n++]=i;
4721        }
4722      }
4723      setInteger(index,n);
4724      delete [] index;
4725      if (n) 
4726        modelPtr_->copyInIntegerInformation(integer);
4727    }
4728    if (keepNames) {
4729      // keep names
4730      int nameDiscipline;
4731      getIntParam(OsiNameDiscipline,nameDiscipline) ;
4732      int iRow;
4733      std::vector<std::string> rowNames = std::vector<std::string> ();
4734      std::vector<std::string> columnNames = std::vector<std::string> ();
4735      rowNames.reserve(nRows);
4736      for (iRow=0;iRow<nRows;iRow++) {
4737        const char * name = m.rowName(iRow);
4738        rowNames.push_back(name);
4739        if (nameDiscipline) 
4740          OsiSolverInterface::setRowName(iRow,name) ;
4741      }
4742     
4743      int iColumn;
4744      columnNames.reserve(nCols);
4745      for (iColumn=0;iColumn<nCols;iColumn++) {
4746        const char * name = m.columnName(iColumn);
4747        columnNames.push_back(name);
4748        if (nameDiscipline) 
4749          OsiSolverInterface::setColName(iColumn,name) ;
4750      }
4751      modelPtr_->copyNames(rowNames,columnNames);
4752    }
4753  }
4754  return numberErrors;
4755}
4756// Read file in LP format (with names)
4757int 
4758OsiClpSolverInterface::readLp(const char *filename, const double epsilon )
4759{
4760  CoinLpIO m;
4761  m.readLp(filename, epsilon);
4762  freeCachedResults();
4763
4764  // set objective function offest
4765  setDblParam(OsiObjOffset, 0);
4766
4767  // set problem name
4768  setStrParam(OsiProbName, m.getProblemName());
4769
4770  // set objective name
4771  setObjName(m.getObjName());
4772
4773  // no errors
4774  loadProblem(*m.getMatrixByRow(), m.getColLower(), m.getColUpper(),
4775              m.getObjCoefficients(), m.getRowLower(), m.getRowUpper());
4776
4777  const char *integer = m.integerColumns();
4778  int nCols = m.getNumCols();
4779  int nRows = m.getNumRows();
4780  if (integer) {
4781    int i, n = 0;
4782    int *index = new int [nCols];
4783    for (i=0; i<nCols; i++) {
4784      if (integer[i]) {
4785        index[n++] = i;
4786      }
4787    }
4788    setInteger(index,n);
4789    delete [] index;
4790  }
4791  // Always keep names
4792  int nameDiscipline;
4793  getIntParam(OsiNameDiscipline,nameDiscipline) ;
4794  int iRow;
4795  std::vector<std::string> rowNames = std::vector<std::string> ();
4796  std::vector<std::string> columnNames = std::vector<std::string> ();
4797  rowNames.reserve(nRows);
4798  for (iRow=0;iRow<nRows;iRow++) {
4799    const char * name = m.rowName(iRow);
4800    rowNames.push_back(name);
4801    if (nameDiscipline) 
4802      OsiSolverInterface::setRowName(iRow,name) ;
4803  }
4804 
4805  int iColumn;
4806  columnNames.reserve(nCols);
4807  for (iColumn=0;iColumn<nCols;iColumn++) {
4808    const char * name = m.columnName(iColumn);
4809    columnNames.push_back(name);
4810    if (nameDiscipline) 
4811      OsiSolverInterface::setColName(iColumn,name) ;
4812  }
4813  modelPtr_->copyNames(rowNames,columnNames);
4814  return(0);
4815}
4816/* Write the problem into an Lp file of the given filename.
4817   If objSense is non zero then -1.0 forces the code to write a
4818   maximization objective and +1.0 to write a minimization one.
4819   If 0.0 then solver can do what it wants.
4820   This version calls writeLpNative with names */
4821void 
4822OsiClpSolverInterface::writeLp(const char *filename,
4823                               const char *extension ,
4824                               double epsilon ,
4825                               int numberAcross ,
4826                               int decimals ,
4827                               double objSense ,
4828                               bool changeNameOnRange) const
4829{
4830  std::string f(filename);
4831  std::string e(extension);
4832  std::string fullname;
4833  if (e!="") {
4834    fullname = f + "." + e;
4835  } else {
4836    // no extension so no trailing period
4837    fullname = f;
4838  }
4839  // get names
4840  const char * const * const rowNames = modelPtr_->rowNamesAsChar();
4841  const char * const * const columnNames = modelPtr_->columnNamesAsChar();
4842  // Fall back on Osi version - possibly with names
4843  OsiSolverInterface::writeLpNative(fullname.c_str(), 
4844                                    rowNames,columnNames, epsilon, numberAcross,
4845                                    decimals, objSense,changeNameOnRange);
4846  if (rowNames) {
4847    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_+1);
4848    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
4849  }
4850}
4851void 
4852OsiClpSolverInterface::writeLp(FILE * fp,
4853                               double epsilon ,
4854                               int numberAcross ,
4855                               int decimals ,
4856                               double objSense ,
4857                               bool changeNameOnRange) const
4858{
4859  // get names
4860  const char * const * const rowNames = modelPtr_->rowNamesAsChar();
4861  const char * const * const columnNames = modelPtr_->columnNamesAsChar();
4862  // Fall back on Osi version - possibly with names
4863  OsiSolverInterface::writeLpNative(fp,
4864                                    rowNames,columnNames, epsilon, numberAcross,
4865                                    decimals, objSense,changeNameOnRange);
4866  if (rowNames) {
4867    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_+1);
4868    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
4869  }
4870}
4871/*
4872  I (JJF) am getting incredibly annoyed because I can't just replace a matrix.
4873  The default behavior of this is do nothing so only use where that would not matter
4874  e.g. strengthening a matrix for MIP
4875*/
4876void 
4877OsiClpSolverInterface::replaceMatrixOptional(const CoinPackedMatrix & matrix)
4878{
4879  modelPtr_->whatsChanged_ &= (0xffff&~(2|4|8));
4880  replaceMatrix(matrix);
4881}
4882// And if it does matter (not used at present)
4883void 
4884OsiClpSolverInterface::replaceMatrix(const CoinPackedMatrix & matrix)
4885{
4886  modelPtr_->whatsChanged_ &= (0xffff&~(2|4|8));
4887  delete modelPtr_->matrix_;
4888  delete modelPtr_->rowCopy_;
4889  modelPtr_->rowCopy_=NULL;
4890  if (matrix.isColOrdered()) {
4891    modelPtr_->matrix_=new ClpPackedMatrix(matrix);
4892  } else {
4893    CoinPackedMatrix matrix2;
4894    matrix2.setExtraGap(0.0);
4895    matrix2.setExtraMajor(0.0);
4896    matrix2.reverseOrderedCopyOf(matrix);
4897    modelPtr_->matrix_=new ClpPackedMatrix(matrix2);
4898  }   
4899  modelPtr_->matrix_->setDimensions(modelPtr_->numberRows_,modelPtr_->numberColumns_);
4900  freeCachedResults();
4901}
4902// Get pointer to array[getNumCols()] of primal solution vector
4903const double * 
4904OsiClpSolverInterface::getColSolution() const 
4905{ 
4906  if (modelPtr_->solveType()!=2) {
4907    return modelPtr_->primalColumnSolution();
4908  } else {
4909    // simplex interface
4910    return modelPtr_->solutionRegion(1);
4911  }
4912}
4913 
4914// Get pointer to array[getNumRows()] of dual prices
4915const double * 
4916OsiClpSolverInterface::getRowPrice() const
4917{ 
4918  if (modelPtr_->solveType()!=2) {
4919    return modelPtr_->dualRowSolution();
4920  } else {
4921    // simplex interface
4922    //return modelPtr_->djRegion(0);
4923    return modelPtr_->dualRowSolution();
4924  }
4925}
4926 
4927// Get a pointer to array[getNumCols()] of reduced costs
4928const double * 
4929OsiClpSolverInterface::getReducedCost() const 
4930{ 
4931  if (modelPtr_->solveType()!=2) {
4932    return modelPtr_->dualColumnSolution();
4933  } else {
4934    // simplex interface
4935    return modelPtr_->djRegion(1);
4936  }
4937}
4938
4939/* Get pointer to array[getNumRows()] of row activity levels (constraint
4940   matrix times the solution vector */
4941const double * 
4942OsiClpSolverInterface::getRowActivity() const 
4943{ 
4944  if (modelPtr_->solveType()!=2) {
4945    return modelPtr_->primalRowSolution();
4946  } else {
4947    // simplex interface
4948    return modelPtr_->solutionRegion(0);
4949  }
4950}
4951double 
4952OsiClpSolverInterface::getObjValue() const 
4953{
4954  if (modelPtr_->numberIterations()||modelPtr_->upperIn_!=-COIN_DBL_MAX) {
4955    // This does not pass unitTest when getObjValue is called before solve.
4956    //printf("obj a %g %g\n",modelPtr_->objectiveValue(),
4957    //     OsiSolverInterface::getObjValue());
4958    return modelPtr_->objectiveValue();
4959  } else {
4960    return OsiSolverInterface::getObjValue();
4961  }
4962}
4963
4964/* Set an objective function coefficient */
4965void 
4966OsiClpSolverInterface::setObjCoeff( int elementIndex, double elementValue )
4967{
4968  modelPtr_->whatsChanged_ &= 0xffff;
4969  // Say can't gurantee optimal basis etc
4970  lastAlgorithm_=999;
4971#ifndef NDEBUG
4972  int n = modelPtr_->numberColumns();
4973  if (elementIndex<0||elementIndex>=n) {
4974    indexError(elementIndex,"setObjCoeff");
4975  }
4976#endif
4977  modelPtr_->setObjectiveCoefficient(elementIndex,elementValue);
4978}
4979
4980/* Set a single column lower bound<br>
4981   Use -DBL_MAX for -infinity. */
4982void 
4983OsiClpSolverInterface::setColLower( int index, double elementValue )
4984{
4985  modelPtr_->whatsChanged_ &= 0x1ffff;
4986#ifndef NDEBUG
4987  int n = modelPtr_->numberColumns();
4988  if (index<0||index>=n) {
4989    indexError(index,"setColLower");
4990  }
4991#endif
4992  double currentValue = modelPtr_->columnActivity_[index];
4993  bool changed=(currentValue<elementValue-modelPtr_->primalTolerance()||
4994                index>=basis_.getNumStructural()||
4995                basis_.getStructStatus(index)==CoinWarmStartBasis::atLowerBound);
4996  // Say can't gurantee optimal basis etc
4997  if (changed)
4998    lastAlgorithm_=999;
4999  if (!modelPtr_->lower_)
5000    modelPtr_->whatsChanged_ &= ~0xffff; // switch off
5001  modelPtr_->setColumnLower(index,elementValue);
5002}
5003     
5004/* Set a single column upper bound<br>
5005   Use DBL_MAX for infinity. */
5006void 
5007OsiClpSolverInterface::setColUpper( int index, double elementValue )
5008{
5009  modelPtr_->whatsChanged_ &= 0x1ffff;
5010#ifndef NDEBUG
5011  int n = modelPtr_->numberColumns();
5012  if (index<0||index>=n) {
5013    indexError(index,"setColUpper");
5014  }
5015#endif
5016  double currentValue = modelPtr_->columnActivity_[index];
5017  bool changed=(currentValue>elementValue+modelPtr_->primalTolerance()||
5018                index>=basis_.getNumStructural()||
5019                basis_.getStructStatus(index)==CoinWarmStartBasis::atUpperBound);
5020  // Say can't gurantee optimal basis etc
5021  if (changed)
5022    lastAlgorithm_=999;
5023  if (!modelPtr_->upper_)
5024    modelPtr_->whatsChanged_ &= ~0xffff; // switch off
5025  modelPtr_->setColumnUpper(index,elementValue);
5026}
5027
5028/* Set a single column lower and upper bound */
5029void 
5030OsiClpSolverInterface::setColBounds( int elementIndex,
5031                                     double lower, double upper )
5032{
5033  modelPtr_->whatsChanged_ &= 0x1ffff;
5034  // Say can't gurantee optimal basis etc
5035  lastAlgorithm_=999;
5036#ifndef NDEBUG
5037  int n = modelPtr_->numberColumns();
5038  if (elementIndex<0||elementIndex>=n) {
5039    indexError(elementIndex,"setColBounds");
5040  }
5041#endif
5042  if (!modelPtr_->lower_)
5043    modelPtr_->whatsChanged_ &= ~0xffff; // switch off
5044  modelPtr_->setColumnBounds(elementIndex,lower,upper);
5045}
5046void OsiClpSolverInterface::setColSetBounds(const int* indexFirst,
5047                                            const int* indexLast,
5048                                            const double* boundList)
5049{
5050  modelPtr_->whatsChanged_ &= 0x1ffff;
5051  // Say can't gurantee optimal basis etc
5052  lastAlgorithm_=999;
5053#ifndef NDEBUG
5054  int n = modelPtr_->numberColumns();
5055  const int * indexFirst2=indexFirst;
5056  while (indexFirst2 != indexLast) {
5057    const int iColumn=*indexFirst2++;
5058    if (iColumn<0||iColumn>=n) {
5059      indexError(iColumn,"setColSetBounds");
5060    }
5061  }
5062#endif
5063  modelPtr_->setColSetBounds(indexFirst,indexLast,boundList);
5064}
5065//------------------------------------------------------------------
5066/* Set a single row lower bound<br>
5067   Use -DBL_MAX for -infinity. */
5068void 
5069OsiClpSolverInterface::setRowLower( int elementIndex, double elementValue ) {
5070  // Say can't gurantee optimal basis etc
5071  lastAlgorithm_=999;
5072  modelPtr_->whatsChanged_ &= 0xffff;
5073#ifndef NDEBUG
5074  int n = modelPtr_->numberRows();
5075  if (elementIndex<0||elementIndex>=n) {
5076    indexError(elementIndex,"setRowLower");
5077  }
5078#endif
5079  modelPtr_->setRowLower(elementIndex , elementValue);
5080  if (rowsense_!=NULL) {
5081    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5082    convertBoundToSense(modelPtr_->rowLower_[elementIndex], 
5083                        modelPtr_->rowUpper_[elementIndex],
5084                        rowsense_[elementIndex], rhs_[elementIndex], rowrange_[elementIndex]);
5085  }
5086}
5087     
5088/* Set a single row upper bound<br>
5089   Use DBL_MAX for infinity. */
5090void 
5091OsiClpSolverInterface::setRowUpper( int elementIndex, double elementValue ) {
5092  modelPtr_->whatsChanged_ &= 0xffff;
5093  // Say can't gurantee optimal basis etc
5094  lastAlgorithm_=999;
5095#ifndef NDEBUG
5096  int n = modelPtr_->numberRows();
5097  if (elementIndex<0||elementIndex>=n) {
5098    indexError(elementIndex,"setRowUpper");
5099  }
5100#endif
5101  modelPtr_->setRowUpper(elementIndex , elementValue);
5102  if (rowsense_!=NULL) {
5103    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5104    convertBoundToSense(modelPtr_->rowLower_[elementIndex], 
5105                        modelPtr_->rowUpper_[elementIndex],
5106                        rowsense_[elementIndex], rhs_[elementIndex], rowrange_[elementIndex]);
5107  }
5108}
5109   
5110/* Set a single row lower and upper bound */
5111void 
5112OsiClpSolverInterface::setRowBounds( int elementIndex,
5113              double lower, double upper ) {
5114  modelPtr_->whatsChanged_ &= 0xffff;
5115  // Say can't gurantee optimal basis etc
5116  lastAlgorithm_=999;
5117#ifndef NDEBUG
5118  int n = modelPtr_->numberRows();
5119  if (elementIndex<0||elementIndex>=n) {
5120    indexError(elementIndex,"setRowBounds");
5121  }
5122#endif
5123  modelPtr_->setRowBounds(elementIndex,lower,upper);
5124  if (rowsense_!=NULL) {
5125    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5126    convertBoundToSense(modelPtr_->rowLower_[elementIndex], 
5127                        modelPtr_->rowUpper_[elementIndex],
5128                        rowsense_[elementIndex], rhs_[elementIndex], rowrange_[elementIndex]);
5129  }
5130}
5131//-----------------------------------------------------------------------------
5132void
5133OsiClpSolverInterface::setRowType(int i, char sense, double rightHandSide,
5134                                  double range)
5135{
5136  modelPtr_->whatsChanged_ &= 0xffff;
5137  // Say can't gurantee optimal basis etc
5138  lastAlgorithm_=999;
5139#ifndef NDEBUG
5140  int n = modelPtr_->numberRows();
5141  if (i<0||i>=n) {
5142    indexError(i,"setRowType");
5143  }
5144#endif
5145  double lower = 0, upper = 0;
5146  convertSenseToBound(sense, rightHandSide, range, lower, upper);
5147  setRowBounds(i, lower, upper);
5148  // If user is using sense then set
5149  if (rowsense_) {
5150    rowsense_[i] = sense;
5151    rhs_[i] = rightHandSide;
5152    rowrange_[i] = range;
5153  }
5154}
5155// Set name of row
5156void 
5157//OsiClpSolverInterface::setRowName(int rowIndex, std::string & name)
5158OsiClpSolverInterface::setRowName(int rowIndex, std::string name) 
5159{
5160  if (rowIndex>=0&&rowIndex<modelPtr_->numberRows()) {
5161    int nameDiscipline;
5162    getIntParam(OsiNameDiscipline,nameDiscipline) ;
5163    if (nameDiscipline) {
5164      modelPtr_->setRowName(rowIndex,name);
5165      OsiSolverInterface::setRowName(rowIndex,name) ;
5166    }
5167  }
5168}
5169// Return name of row if one exists or Rnnnnnnn
5170// we ignore maxLen
5171std::string
5172OsiClpSolverInterface::getRowName(int rowIndex, unsigned int /*maxLen*/) const
5173{ 
5174        if (rowIndex == getNumRows())
5175                return getObjName();
5176  return modelPtr_->getRowName(rowIndex);
5177}
5178   
5179// Set name of col
5180void 
5181//OsiClpSolverInterface::setColName(int colIndex, std::string & name)
5182OsiClpSolverInterface::setColName(int colIndex, std::string name) 
5183{
5184  if (colIndex>=0&&colIndex<modelPtr_->numberColumns()) {
5185    int nameDiscipline;
5186    getIntParam(OsiNameDiscipline,nameDiscipline) ;
5187    if (nameDiscipline) {
5188      modelPtr_->setColumnName(colIndex,name);
5189      OsiSolverInterface::setColName(colIndex,name) ;
5190    }
5191  }
5192}
5193// Return name of col if one exists or Rnnnnnnn
5194std::string
5195OsiClpSolverInterface::getColName(int colIndex, unsigned int /*maxLen*/) const
5196{
5197  return modelPtr_->getColumnName(colIndex);
5198}
5199   
5200   
5201//-----------------------------------------------------------------------------
5202void OsiClpSolverInterface::setRowSetBounds(const int* indexFirst,
5203                                            const int* indexLast,
5204                                            const double* boundList)
5205{
5206  modelPtr_->whatsChanged_ &= 0xffff;
5207  // Say can't gurantee optimal basis etc
5208  lastAlgorithm_=999;
5209#ifndef NDEBUG
5210  int n = modelPtr_->numberRows();
5211  const int * indexFirst2=indexFirst;
5212  while (indexFirst2 != indexLast) {
5213    const int iColumn=*indexFirst2++;
5214    if (iColumn<0||iColumn>=n) {
5215      indexError(iColumn,"setColumnSetBounds");
5216    }
5217  }
5218#endif
5219  modelPtr_->setRowSetBounds(indexFirst,indexLast,boundList);
5220  if (rowsense_ != NULL) {
5221    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5222    double * lower = modelPtr_->rowLower();
5223    double * upper = modelPtr_->rowUpper();
5224    while (indexFirst != indexLast) {
5225      const int iRow=*indexFirst++;
5226      convertBoundToSense(lower[iRow], upper[iRow],
5227                          rowsense_[iRow], rhs_[iRow], rowrange_[iRow]);
5228    }
5229  }
5230}
5231//-----------------------------------------------------------------------------
5232void
5233OsiClpSolverInterface::setRowSetTypes(const int* indexFirst,
5234                                      const int* indexLast,
5235                                      const char* senseList,
5236                                      const double* rhsList,
5237                                      const double* rangeList)
5238{
5239  modelPtr_->whatsChanged_ &= 0xffff;
5240  // Say can't gurantee optimal basis etc
5241  lastAlgorithm_=999;
5242#ifndef NDEBUG
5243  int n = modelPtr_->numberRows();
5244#endif
5245  const int len = static_cast<int>(indexLast - indexFirst);
5246  while (indexFirst != indexLast) {
5247    const int iRow= *indexFirst++;
5248#ifndef NDEBUG
5249    if (iRow<0||iRow>=n) {
5250      indexError(iRow,"isContinuous");
5251    }
5252#endif
5253    double lowerValue = 0;
5254    double upperValue = 0;
5255    if (rangeList){
5256      convertSenseToBound(*senseList++, *rhsList++, *rangeList++,
5257                          lowerValue, upperValue);
5258    } else {
5259      convertSenseToBound(*senseList++, *rhsList++, 0,
5260                          lowerValue, upperValue);
5261    }
5262    modelPtr_->setRowBounds(iRow,lowerValue,upperValue);
5263  }
5264  if (rowsense_ != NULL) {
5265    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5266    indexFirst -= len;
5267    senseList -= len;
5268    rhsList -= len;
5269    if (rangeList)
5270       rangeList -= len;
5271    while (indexFirst != indexLast) {
5272      const int iRow=*indexFirst++;
5273      rowsense_[iRow] = *senseList++;
5274      rhs_[iRow] = *rhsList++;
5275      if (rangeList)
5276         rowrange_[iRow] = *rangeList++;
5277    }
5278  }
5279}
5280/*Enables normal operation of subsequent functions.
5281  This method is supposed to ensure that all typical things (like
5282  reduced costs, etc.) are updated when individual pivots are executed
5283  and can be queried by other methods
5284*/
5285void 
5286OsiClpSolverInterface::enableSimplexInterface(bool doingPrimal)
5287{
5288  modelPtr_->whatsChanged_ &= 0xffff;
5289  if (modelPtr_->solveType()==2)
5290    return;
5291  assert (modelPtr_->solveType()==1);
5292  int saveIts = modelPtr_->numberIterations_;
5293  modelPtr_->setSolveType(2);
5294  if (doingPrimal)
5295    modelPtr_->setAlgorithm(1);
5296  else
5297    modelPtr_->setAlgorithm(-1);
5298  // Do initialization
5299  saveData_ = modelPtr_->saveData();
5300  saveData_.scalingFlag_=modelPtr_->scalingFlag();
5301  modelPtr_->scaling(0);
5302  specialOptions_ = 0x80000000;
5303  // set infeasibility cost up
5304  modelPtr_->setInfeasibilityCost(1.0e12);
5305  ClpDualRowDantzig dantzig;
5306  modelPtr_->setDualRowPivotAlgorithm(dantzig);
5307  ClpPrimalColumnDantzig dantzigP;
5308  dantzigP.saveWeights(modelPtr_,0); // set modelPtr
5309  modelPtr_->setPrimalColumnPivotAlgorithm(dantzigP);
5310  int saveOptions = modelPtr_->specialOptions_;
5311  modelPtr_->specialOptions_ &= ~262144;
5312  delete modelPtr_->scaledMatrix_;
5313  modelPtr_->scaledMatrix_=NULL;
5314#ifdef NDEBUG
5315  modelPtr_->startup(0);
5316#else
5317  int returnCode=modelPtr_->startup(0);
5318  assert (!returnCode||returnCode==2);
5319#endif
5320  modelPtr_->specialOptions_=saveOptions;
5321  modelPtr_->numberIterations_=saveIts;
5322}
5323
5324//Undo whatever setting changes the above method had to make
5325void 
5326OsiClpSolverInterface::disableSimplexInterface()
5327{
5328  modelPtr_->whatsChanged_ &= 0xffff;
5329  assert (modelPtr_->solveType()==2);
5330  // declare optimality anyway  (for message handler)
5331  modelPtr_->setProblemStatus(0);
5332  modelPtr_->setSolveType(1);
5333  // message will not appear anyway
5334  int saveMessageLevel=modelPtr_->messageHandler()->logLevel();
5335  modelPtr_->messageHandler()->setLogLevel(0);
5336  modelPtr_->finish();
5337  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
5338  modelPtr_->restoreData(saveData_);
5339  modelPtr_->scaling(saveData_.scalingFlag_);
5340  ClpDualRowSteepest steepest;
5341  modelPtr_->setDualRowPivotAlgorithm(steepest);
5342  ClpPrimalColumnSteepest steepestP;
5343  modelPtr_->setPrimalColumnPivotAlgorithm(steepestP);
5344  basis_ = getBasis(modelPtr_);
5345  modelPtr_->setSolveType(1);
5346}
5347void 
5348OsiClpSolverInterface::enableFactorization() const
5349{
5350  specialOptions_ &= ~0x80000000;
5351  saveData_.specialOptions_=specialOptions_;
5352  int saveStatus = modelPtr_->problemStatus_;
5353  if ((specialOptions_&(1+8))!=1+8)
5354    setSpecialOptionsMutable((1+8)|specialOptions_);
5355#ifdef NDEBUG
5356  modelPtr_->startup(0);
5357#else
5358  int returnCode=modelPtr_->startup(0);
5359  assert (!returnCode||returnCode==2);
5360#endif
5361  modelPtr_->problemStatus_=saveStatus;
5362}
5363
5364//Undo whatever setting changes the above method had to make
5365void 
5366OsiClpSolverInterface::disableFactorization() const
5367{
5368  specialOptions_=saveData_.specialOptions_;
5369  // declare optimality anyway  (for message handler)
5370  modelPtr_->setProblemStatus(0);
5371  // message will not appear anyway
5372  int saveMessageLevel=modelPtr_->messageHandler()->logLevel();
5373  modelPtr_->messageHandler()->setLogLevel(0);
5374  // Should re-do - for moment save arrays
5375  double * sol = CoinCopyOfArray(modelPtr_->columnActivity_,modelPtr_->numberColumns_);
5376  double * dj = CoinCopyOfArray(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
5377  double * rsol = CoinCopyOfArray(modelPtr_->rowActivity_,modelPtr_->numberRows_);
5378  double * dual = CoinCopyOfArray(modelPtr_->dual_,modelPtr_->numberRows_);
5379  modelPtr_->finish();
5380  CoinMemcpyN(sol,modelPtr_->numberColumns_,modelPtr_->columnActivity_);
5381  CoinMemcpyN(dj,modelPtr_->numberColumns_,modelPtr_->reducedCost_);
5382  CoinMemcpyN(rsol,modelPtr_->numberRows_,modelPtr_->rowActivity_);
5383  CoinMemcpyN(dual,modelPtr_->numberRows_,modelPtr_->dual_);
5384  delete [] sol;
5385  delete [] dj;
5386  delete [] rsol;
5387  delete [] dual;
5388  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
5389}
5390/* The following two methods may be replaced by the
5391   methods of OsiSolverInterface using OsiWarmStartBasis if:
5392   1. OsiWarmStartBasis resize operation is implemented
5393   more efficiently and
5394   2. It is ensured that effects on the solver are the same
5395   
5396   Returns a basis status of the structural/artificial variables
5397*/
5398void 
5399OsiClpSolverInterface::getBasisStatus(int* cstat, int* rstat) const
5400{
5401  int iRow,iColumn;
5402  int numberRows = modelPtr_->numberRows();
5403  int numberColumns = modelPtr_->numberColumns();
5404  const double * pi = modelPtr_->dualRowSolution();
5405  const double * dj = modelPtr_->dualColumnSolution();
5406  double multiplier = modelPtr_->optimizationDirection();
5407  // Flip slacks
5408  int lookupA[]={0,1,3,2,0,3};
5409  for (iRow=0;iRow<numberRows;iRow++) {
5410    int iStatus = modelPtr_->getRowStatus(iRow);
5411    if (iStatus==5) {
5412      // Fixed - look at reduced cost
5413      if (pi[iRow]*multiplier>1.0e-7)
5414        iStatus = 3;
5415    }
5416    iStatus = lookupA[iStatus];
5417    rstat[iRow]=iStatus;
5418  }
5419  int lookupS[]={0,1,2,3,0,3};
5420  for (iColumn=0;iColumn<numberColumns;iColumn++) {
5421    int iStatus = modelPtr_->getColumnStatus(iColumn);
5422    if (iStatus==5) {
5423      // Fixed - look at reduced cost
5424      if (dj[iColumn]*multiplier<-1.0e-7)
5425        iStatus = 2;
5426    }
5427    iStatus = lookupS[iStatus];
5428    cstat[iColumn]=iStatus;
5429  }
5430}
5431
5432//Set the status of structural/artificial variables
5433//Returns 0 if OK, 1 if problem is bad e.g. duplicate elements, too large ...
5434int 
5435OsiClpSolverInterface::setBasisStatus(const int* cstat, const int* rstat)
5436{
5437  modelPtr_->whatsChanged_ &= 0xffff;
5438  // Say can't gurantee optimal basis etc
5439  lastAlgorithm_=999;
5440  modelPtr_->createStatus();
5441  int i, n;
5442  double * lower, * upper, * solution;
5443  n=modelPtr_->numberRows();
5444  lower = modelPtr_->rowLower();
5445  upper = modelPtr_->rowUpper();
5446  solution = modelPtr_->primalRowSolution();
5447  // For rows lower and upper are flipped
5448  int lookupA[]={0,1,3,2};
5449  for (i=0;i<n;i++) {
5450    int status = lookupA[rstat[i]];
5451    if (status<0||status>3)
5452      status = 3;
5453    if (lower[i]<-1.0e50&&upper[i]>1.0e50&&status!=1)
5454      status = 0; // set free if should be
5455    else if (lower[i]<-1.0e50&&status==3)
5456      status = 2; // can't be at lower bound
5457    else if (upper[i]>1.0e50&&status==2)
5458      status = 3; // can't be at upper bound
5459    switch (status) {
5460      // free or superbasic
5461    case 0:
5462      if (lower[i]<-1.0e50&&upper[i]>1.0e50) {
5463        modelPtr_->setRowStatus(i,ClpSimplex::isFree);
5464        if (fabs(solution[i])>1.0e20)
5465          solution[i]=0.0;
5466      } else {
5467        modelPtr_->setRowStatus(i,ClpSimplex::superBasic);
5468        if (fabs(solution[i])>1.0e20)
5469          solution[i]=0.0;
5470      }
5471      break;
5472    case 1:
5473      // basic
5474      modelPtr_->setRowStatus(i,ClpSimplex::basic);
5475      break;
5476    case 2:
5477      // at upper bound
5478      solution[i]=upper[i];
5479      if (upper[i]>lower[i])
5480        modelPtr_->setRowStatus(i,ClpSimplex::atUpperBound);
5481      else
5482        modelPtr_->setRowStatus(i,ClpSimplex::isFixed);
5483      break;
5484    case 3:
5485      // at lower bound
5486      solution[i]=lower[i];
5487      if (upper[i]>lower[i])
5488        modelPtr_->setRowStatus(i,ClpSimplex::atLowerBound);
5489      else
5490        modelPtr_->setRowStatus(i,ClpSimplex::isFixed);
5491      break;
5492    }
5493  }
5494  n=modelPtr_->numberColumns();
5495  lower = modelPtr_->columnLower();
5496  upper = modelPtr_->columnUpper();
5497  solution = modelPtr_->primalColumnSolution();
5498  for (i=0;i<n;i++) {
5499    int status = cstat[i];
5500    if (status<0||status>3)
5501      status = 3;
5502    if (lower[i]<-1.0e50&&upper[i]>1.0e50&&status!=1)
5503      status = 0; // set free if should be
5504    else if (lower[i]<-1.0e50&&status==3)
5505      status = 2; // can't be at lower bound
5506    else if (upper[i]>1.0e50&&status==2)
5507      status = 3; // can't be at upper bound
5508    switch (status) {
5509      // free or superbasic
5510    case 0:
5511      if (lower[i]<-1.0e50&&upper[i]>1.0e50) {
5512        modelPtr_->setColumnStatus(i,ClpSimplex::isFree);
5513        if (fabs(solution[i])>1.0e20)
5514          solution[i]=0.0;
5515      } else {
5516        modelPtr_->setColumnStatus(i,ClpSimplex::superBasic);
5517        if (fabs(solution[i])>1.0e20)
5518          solution[i]=0.0;
5519      }
5520      break;
5521    case 1:
5522      // basic
5523      modelPtr_->setColumnStatus(i,ClpSimplex::basic);
5524      break;
5525    case 2:
5526      // at upper bound
5527      solution[i]=upper[i];
5528      if (upper[i]>lower[i])
5529        modelPtr_->setColumnStatus(i,ClpSimplex::atUpperBound);
5530      else
5531        modelPtr_->setColumnStatus(i,ClpSimplex::isFixed);
5532      break;
5533    case 3:
5534      // at lower bound
5535      solution[i]=lower[i];
5536      if (upper[i]>lower[i])
5537        modelPtr_->setColumnStatus(i,ClpSimplex::atLowerBound);
5538      else
5539        modelPtr_->setColumnStatus(i,ClpSimplex::isFixed);
5540      break;
5541    }
5542  }
5543  // say first time
5544  modelPtr_->statusOfProblem(true);
5545  // May be bad model
5546  if (modelPtr_->status()==4)
5547    return 1;
5548  // Save
5549  basis_ = getBasis(modelPtr_);
5550  return 0;
5551}
5552
5553/* Perform a pivot by substituting a colIn for colOut in the basis.
5554   The status of the leaving variable is given in statOut. Where
5555   1 is to upper bound, -1 to lower bound
5556   Return code is 0 for okay,
5557   1 if inaccuracy forced re-factorization (should be okay) and
5558   -1 for singular factorization
5559*/
5560int 
5561OsiClpSolverInterface::pivot(int colIn, int colOut, int outStatus)
5562{
5563  assert (modelPtr_->solveType()==2);
5564  // convert to Clp style (what about flips?)
5565  if (colIn<0) 
5566    colIn = modelPtr_->numberColumns()+(-1-colIn);
5567  if (colOut<0) 
5568    colOut = modelPtr_->numberColumns()+(-1-colOut);
5569  // in clp direction of out is reversed
5570  outStatus = - outStatus;
5571  // set in clp
5572  modelPtr_->setDirectionOut(outStatus);
5573  modelPtr_->setSequenceIn(colIn);
5574  modelPtr_->setSequenceOut(colOut);
5575  // do pivot
5576  return modelPtr_->pivot();
5577}
5578
5579/* Obtain a result of the primal pivot
5580   Outputs: colOut -- leaving column, outStatus -- its status,
5581   t -- step size, and, if dx!=NULL, *dx -- primal ray direction.
5582   Inputs: colIn -- entering column, sign -- direction of its change (+/-1).
5583   Both for colIn and colOut, artificial variables are index by
5584   the negative of the row index minus 1.
5585   Return code (for now): 0 -- leaving variable found,
5586   -1 -- everything else?
5587   Clearly, more informative set of return values is required
5588   Primal and dual solutions are updated
5589*/
5590int 
5591OsiClpSolverInterface::primalPivotResult(int colIn, int sign, 
5592                                         int& colOut, int& outStatus, 
5593                                         double& t, CoinPackedVector* dx)
5594{
5595  assert (modelPtr_->solveType()==2);
5596  // convert to Clp style
5597  if (colIn<0) 
5598    colIn = modelPtr_->numberColumns()+(-1-colIn);
5599  // set in clp
5600  modelPtr_->setDirectionIn(sign);
5601  modelPtr_->setSequenceIn(colIn);
5602  modelPtr_->setSequenceOut(-1);
5603  int returnCode = modelPtr_->primalPivotResult();
5604  t = modelPtr_->theta();
5605  int numberColumns = modelPtr_->numberColumns();
5606  if (dx) {
5607    double * ray = modelPtr_->unboundedRay();
5608    if  (ray)
5609      dx->setFullNonZero(numberColumns,ray);
5610    else
5611      printf("No ray?\n");
5612    delete [] ray;
5613  }
5614  outStatus = - modelPtr_->directionOut();
5615  colOut = modelPtr_->sequenceOut();
5616  if (colOut>= numberColumns) 
5617    colOut = -1-(colOut - numberColumns);
5618  return returnCode;
5619}
5620
5621/* Obtain a result of the dual pivot (similar to the previous method)
5622   Differences: entering variable and a sign of its change are now
5623   the outputs, the leaving variable and its statuts -- the inputs
5624   If dx!=NULL, then *dx contains dual ray
5625   Return code: same
5626*/
5627int 
5628OsiClpSolverInterface::dualPivotResult(int& /*colIn*/, int& /*sign*/, 
5629                                       int /*colOut*/, int /*outStatus*/, 
5630                                       double& /*t*/, CoinPackedVector* /*dx*/)
5631{
5632  assert (modelPtr_->solveType()==2);
5633  abort();
5634  return 0;
5635}
5636
5637//Get the reduced gradient for the cost vector c
5638void 
5639OsiClpSolverInterface::getReducedGradient(
5640                                          double* columnReducedCosts, 
5641                                          double * duals,
5642                                          const double * c)
5643{
5644  assert (modelPtr_->solveType()==2);
5645  // could do this faster with coding inside Clp
5646  // save current costs
5647  int numberColumns = modelPtr_->numberColumns();
5648  double