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

Last change on this file since 2046 was 2046, checked in by tkr, 5 years ago

Merging r2045 from stable/1.15. Whoops, meant to commit it to trunk

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