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

Last change on this file since 1664 was 1664, checked in by lou, 9 years ago

Add EPL license notice in OsiClp?.

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