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

Last change on this file since 1924 was 1924, checked in by forrest, 7 years ago

changes for miqp, parameters and presolve

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