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

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

Remove setObjectiveAndRefresh. Make getReducedGradient const. A few other
minor tweaks.

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