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

Last change on this file since 1742 was 1742, checked in by forrest, 8 years ago

setRowPrice setColSolution to pass tests

File size: 297.4 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      if (small->status()==1) {
2042#ifndef KEEP_SMALL
2043        if (small!=modelPtr_)
2044          delete small;
2045        //delete smallModel_;
2046        //smallModel_=NULL;
2047        assert (!smallModel_);
2048#else
2049        assert (small==smallModel_);
2050        if (smallModel_!=modelPtr_) {
2051          delete smallModel_;
2052        }
2053        smallModel_=NULL;
2054#endif
2055        delete [] spareArrays_;
2056        spareArrays_=NULL;
2057        delete ws_;
2058        ws_ = dynamic_cast<CoinWarmStartBasis*>(getWarmStart());
2059        int numberRows = modelPtr_->numberRows();
2060        rowActivity_= new double[numberRows];
2061 CoinMemcpyN(modelPtr_->primalRowSolution(),    numberRows,rowActivity_);
2062        int numberColumns = modelPtr_->numberColumns();
2063        columnActivity_= new double[numberColumns];
2064 CoinMemcpyN(modelPtr_->primalColumnSolution(), numberColumns,columnActivity_);
2065        modelPtr_->setProblemStatus(1);
2066        if (savedObjective) {
2067          // fix up
2068          modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2069          //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2070          modelPtr_->objective_=savedObjective;
2071        }
2072        return;
2073      } else {
2074        // update model
2075        static_cast<ClpSimplexOther *> (modelPtr_)->afterCrunch(*small,whichRow,whichColumn,nBound);
2076      } 
2077#ifndef SETUP_HOT
2078      assert (factorization_==NULL);
2079      factorization_ = static_cast<ClpSimplexDual *>(small)->setupForStrongBranching(spareArrays_,numberRows,
2080                             numberColumns,false);
2081#else
2082      assert (factorization!=NULL || small->problemStatus_ );
2083      factorization_ = factorization;
2084#endif
2085    } else {
2086      assert (factorization_==NULL);
2087      factorization_ = static_cast<ClpSimplexDual *>(small)->setupForStrongBranching(spareArrays_,numberRows,
2088                                                                                     numberColumns,false);
2089    }
2090    smallModel_=small;
2091    if (modelPtr_->logLevel()<2) smallModel_->setLogLevel(0);
2092    // Setup for strong branching
2093    int numberColumns2 = smallModel_->numberColumns();
2094    CoinMemcpyN( modelPtr_->columnLower(),numberColumns, saveLowerOriginal);
2095    CoinMemcpyN( modelPtr_->columnUpper(),numberColumns, saveUpperOriginal);
2096    const double * smallLower = smallModel_->columnLower();
2097    const double * smallUpper = smallModel_->columnUpper();
2098    // But modify if bounds changed in small
2099    for (int i=0;i<numberColumns2;i++) {
2100      int iColumn = whichColumn[i];
2101      saveLowerOriginal[iColumn] = CoinMax(saveLowerOriginal[iColumn],
2102                                           smallLower[i]);
2103      saveUpperOriginal[iColumn] = CoinMin(saveUpperOriginal[iColumn],
2104                                           smallUpper[i]);
2105    }
2106    if (whichRange_&&whichRange_[0]) {
2107      // get ranging information
2108      int numberToDo = whichRange_[0];
2109      int * which = new int [numberToDo];
2110      // Convert column numbers
2111      int * backColumn = whichColumn+numberColumns;
2112      for (int i=0;i<numberToDo;i++) {
2113        int iColumn = whichRange_[i+1];
2114        which[i]=backColumn[iColumn];
2115      }
2116      double * downRange=new double [numberToDo];
2117      double * upRange=new double [numberToDo];
2118      int * whichDown = new int [numberToDo];
2119      int * whichUp = new int [numberToDo];
2120      smallModel_->setFactorization(*factorization_);
2121      smallModel_->gutsOfSolution(NULL,NULL,false);
2122      // Tell code we can increase costs in some cases
2123      smallModel_->setCurrentDualTolerance(0.0);
2124      static_cast<ClpSimplexOther *> (smallModel_)->dualRanging(numberToDo,which,
2125                         upRange, whichUp, downRange, whichDown);
2126      delete [] whichDown;
2127      delete [] whichUp;
2128      delete [] which;
2129      rowActivity_=upRange;
2130      columnActivity_=downRange;
2131    }
2132  }
2133    if (savedObjective) {
2134      // fix up
2135      modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2136      //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2137      modelPtr_->objective_=savedObjective;
2138      if (!modelPtr_->problemStatus_) {
2139        CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
2140        CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
2141        if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
2142          CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
2143        modelPtr_->computeObjectiveValue();
2144      }
2145    }
2146}
2147
2148void OsiClpSolverInterface::solveFromHotStart()
2149{
2150#ifdef KEEP_SMALL
2151  if (!spareArrays_) {
2152    assert (!smallModel_);
2153    assert (modelPtr_->problemStatus_==1);
2154    return;
2155  }
2156#endif
2157  ClpObjective * savedObjective=NULL;
2158  double savedDualLimit=modelPtr_->dblParam_[ClpDualObjectiveLimit];
2159  if (saveData_.perturbation_) {
2160    // Set fake
2161    savedObjective=modelPtr_->objective_;
2162    modelPtr_->objective_=fakeObjective_;
2163    modelPtr_->dblParam_[ClpDualObjectiveLimit]=COIN_DBL_MAX;
2164  }
2165  int numberRows = modelPtr_->numberRows();
2166  int numberColumns = modelPtr_->numberColumns();
2167  modelPtr_->getIntParam(ClpMaxNumIteration,itlimOrig_);
2168  int itlim;
2169  modelPtr_->getIntParam(ClpMaxNumIterationHotStart, itlim);
2170  // Is there an extra copy of scaled bounds
2171  int extraCopy = (modelPtr_->maximumRows_>0) ? modelPtr_->maximumRows_+modelPtr_->maximumColumns_ : 0; 
2172#ifdef CLEAN_HOT_START
2173  if ((specialOptions_&65536)!=0) {
2174    double * arrayD = reinterpret_cast<double *> (spareArrays_);
2175    double saveObjectiveValue = arrayD[0];
2176    double * saveSolution = arrayD+1;
2177    int number = numberRows+numberColumns;
2178    CoinMemcpyN(saveSolution,number,modelPtr_->solutionRegion());
2179    double * saveLower = saveSolution + (numberRows+numberColumns);
2180    CoinMemcpyN(saveLower,number,modelPtr_->lowerRegion());
2181    double * saveUpper = saveLower + (numberRows+numberColumns);
2182    CoinMemcpyN(saveUpper,number,modelPtr_->upperRegion());
2183    double * saveObjective = saveUpper + (numberRows+numberColumns);
2184    CoinMemcpyN(saveObjective,number,modelPtr_->costRegion());
2185    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
2186    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
2187    arrayD = saveUpperOriginal + numberColumns;
2188    int * savePivot = reinterpret_cast<int *> (arrayD);
2189    CoinMemcpyN(savePivot,numberRows,modelPtr_->pivotVariable());
2190    int * whichRow = savePivot+numberRows;
2191    int * whichColumn = whichRow + 3*numberRows;
2192    int * arrayI = whichColumn + 2*numberColumns;
2193    unsigned char * saveStatus = reinterpret_cast<unsigned char *> (arrayI+1);
2194    CoinMemcpyN(saveStatus,number,modelPtr_->statusArray());
2195    modelPtr_->setFactorization(*factorization_);
2196    double * columnLower = modelPtr_->columnLower();
2197    double * columnUpper = modelPtr_->columnUpper();
2198    // make sure whatsChanged_ has 1 set
2199    modelPtr_->setWhatsChanged(511);
2200    double * lowerInternal = modelPtr_->lowerRegion();
2201    double * upperInternal = modelPtr_->upperRegion();
2202    double rhsScale = modelPtr_->rhsScale();
2203    const double * columnScale = NULL;
2204    if (modelPtr_->scalingFlag()>0) 
2205      columnScale = modelPtr_->columnScale() ;
2206    // and do bounds in case dual needs them
2207    int iColumn;
2208    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2209      if (columnLower[iColumn]>saveLowerOriginal[iColumn]) {
2210        double value = columnLower[iColumn];
2211        value *= rhsScale;
2212        if (columnScale)
2213          value /= columnScale[iColumn];
2214        lowerInternal[iColumn]=value;
2215        if (extraCopy)
2216          lowerInternal[iColumn+extraCopy]=value;
2217      }
2218      if (columnUpper[iColumn]<saveUpperOriginal[iColumn]) {
2219        double value = columnUpper[iColumn];
2220        value *= rhsScale;
2221        if (columnScale)
2222          value /= columnScale[iColumn];
2223        upperInternal[iColumn]=value;
2224        if (extraCopy)
2225          upperInternal[iColumn+extraCopy]=value;
2226      }
2227    }
2228    // Start of fast iterations
2229    bool alwaysFinish= ((specialOptions_&32)==0) ? true : false;
2230    //modelPtr_->setLogLevel(1);
2231    modelPtr_->setIntParam(ClpMaxNumIteration,itlim);
2232    int saveNumberFake = (static_cast<ClpSimplexDual *>(modelPtr_))->numberFake_;
2233    int status = (static_cast<ClpSimplexDual *>(modelPtr_))->fastDual(alwaysFinish);
2234    (static_cast<ClpSimplexDual *>(modelPtr_))->numberFake_ = saveNumberFake;
2235   
2236    int problemStatus = modelPtr_->problemStatus();
2237    double objectiveValue =modelPtr_->objectiveValue() * modelPtr_->optimizationDirection();
2238    CoinAssert (modelPtr_->problemStatus()||modelPtr_->objectiveValue()<1.0e50);
2239    // make sure plausible
2240    double obj = CoinMax(objectiveValue,saveObjectiveValue);
2241    if (problemStatus==10||problemStatus<0) {
2242      // was trying to clean up or something odd
2243      if (problemStatus==10)
2244        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2245      status=1;
2246    }
2247    if (status) {
2248      // not finished - might be optimal
2249      modelPtr_->checkPrimalSolution(modelPtr_->solutionRegion(0),
2250                                     modelPtr_->solutionRegion(1));
2251      //modelPtr_->gutsOfSolution(NULL,NULL,0);
2252      //if (problemStatus==3)
2253      //modelPtr_->computeObjectiveValue();
2254      objectiveValue =modelPtr_->objectiveValue() *
2255        modelPtr_->optimizationDirection();
2256      obj = CoinMax(objectiveValue,saveObjectiveValue);
2257      if (!modelPtr_->numberDualInfeasibilities()) { 
2258        double limit = 0.0;
2259        modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
2260        if (modelPtr_->secondaryStatus()==1&&!problemStatus&&obj<limit) {
2261          obj=limit;
2262          problemStatus=3;
2263        }
2264        if (!modelPtr_->numberPrimalInfeasibilities()&&obj<limit) { 
2265          problemStatus=0;
2266        } else if (problemStatus==10) {
2267          problemStatus=3;
2268        } else if (!modelPtr_->numberPrimalInfeasibilities()) {
2269          problemStatus=1; // infeasible
2270        } 
2271      } else {
2272        // can't say much
2273        //if (problemStatus==3)
2274        //modelPtr_->computeObjectiveValue();
2275        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2276        problemStatus=3;
2277      }
2278    } else if (!problemStatus) {
2279      if (modelPtr_->isDualObjectiveLimitReached()) 
2280        problemStatus=1; // infeasible
2281    }
2282    if (status&&!problemStatus) {
2283      problemStatus=3; // can't be sure
2284      lastAlgorithm_=1;
2285    }
2286    if (problemStatus<0)
2287      problemStatus=3;
2288    modelPtr_->setProblemStatus(problemStatus);
2289    modelPtr_->setObjectiveValue(obj*modelPtr_->optimizationDirection());
2290    double * solution = modelPtr_->primalColumnSolution();
2291    const double * solution2 = modelPtr_->solutionRegion();
2292    // could just do changed bounds - also try double size scale so can use * not /
2293    if (!columnScale) {
2294      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2295        solution[iColumn]= solution2[iColumn];
2296      }
2297    } else {
2298      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2299        solution[iColumn]= solution2[iColumn]*columnScale[iColumn];
2300      }
2301    }
2302    CoinMemcpyN(saveLowerOriginal,numberColumns,columnLower);
2303    CoinMemcpyN(saveUpperOriginal,numberColumns,columnUpper);
2304#if 0
2305    // could combine with loop above
2306    if (modelPtr_==modelPtr_)
2307      modelPtr_->computeObjectiveValue();
2308    if (status&&!problemStatus) {
2309      memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2310      modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2311      modelPtr_->checkSolutionInternal();
2312      //modelPtr_->setLogLevel(1);
2313      //modelPtr_->allSlackBasis();
2314      //modelPtr_->primal(1);
2315      //memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2316      //modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2317      //modelPtr_->checkSolutionInternal();
2318      assert (!modelPtr_->problemStatus());
2319    }
2320#endif
2321    // and back bounds
2322    CoinMemcpyN(saveLower,number,modelPtr_->lowerRegion());
2323    CoinMemcpyN(saveUpper,number,modelPtr_->upperRegion());
2324    if (extraCopy) {
2325      CoinMemcpyN(saveLower,number,modelPtr_->lowerRegion()+extraCopy);
2326      CoinMemcpyN(saveUpper,number,modelPtr_->upperRegion()+extraCopy);
2327    }
2328    modelPtr_->setIntParam(ClpMaxNumIteration,itlimOrig_);
2329    if (savedObjective) {
2330      // fix up
2331      modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2332      //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2333      modelPtr_->objective_=savedObjective;
2334      if (!modelPtr_->problemStatus_) {
2335        CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
2336        CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
2337        if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
2338          CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
2339        modelPtr_->computeObjectiveValue();
2340      }
2341    }
2342    return;
2343  }
2344#endif
2345  if (smallModel_==NULL) {
2346    setWarmStart(ws_);
2347    CoinMemcpyN(           rowActivity_,numberRows,modelPtr_->primalRowSolution());
2348    CoinMemcpyN(columnActivity_,numberColumns,modelPtr_->primalColumnSolution());
2349    modelPtr_->setIntParam(ClpMaxNumIteration,CoinMin(itlim,9999));
2350    resolve();
2351  } else {
2352    assert (spareArrays_);
2353    double * arrayD = reinterpret_cast<double *> (spareArrays_);
2354    double saveObjectiveValue = arrayD[0];
2355    double * saveSolution = arrayD+1;
2356    int numberRows2 = smallModel_->numberRows();
2357    int numberColumns2 = smallModel_->numberColumns();
2358    int number = numberRows2+numberColumns2;
2359    CoinMemcpyN(saveSolution,number,smallModel_->solutionRegion());
2360    double * saveLower = saveSolution + (numberRows+numberColumns);
2361    CoinMemcpyN(saveLower,number,smallModel_->lowerRegion());
2362    double * saveUpper = saveLower + (numberRows+numberColumns);
2363    CoinMemcpyN(saveUpper,number,smallModel_->upperRegion());
2364    double * saveObjective = saveUpper + (numberRows+numberColumns);
2365    CoinMemcpyN(saveObjective,number,smallModel_->costRegion());
2366    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
2367    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
2368    arrayD = saveUpperOriginal + numberColumns;
2369    int * savePivot = reinterpret_cast<int *> (arrayD);
2370    CoinMemcpyN(savePivot,numberRows2,smallModel_->pivotVariable());
2371    int * whichRow = savePivot+numberRows;
2372    int * whichColumn = whichRow + 3*numberRows;
2373    int * arrayI = whichColumn + 2*numberColumns;
2374    unsigned char * saveStatus = reinterpret_cast<unsigned char *> (arrayI+1);
2375    CoinMemcpyN(saveStatus,number,smallModel_->statusArray());
2376    smallModel_->setFactorization(*factorization_);
2377    //int * backColumn = whichColumn+numberColumns;
2378    const double * lowerBig = modelPtr_->columnLower();
2379    const double * upperBig = modelPtr_->columnUpper();
2380    // make sure whatsChanged_ has 1 set
2381    smallModel_->setWhatsChanged(511);
2382    double * lowerSmall = smallModel_->lowerRegion();
2383    double * upperSmall = smallModel_->upperRegion();
2384    double * lowerSmallReal = smallModel_->columnLower();
2385    double * upperSmallReal = smallModel_->columnUpper();
2386    int i;
2387    double rhsScale = smallModel_->rhsScale();
2388    const double * columnScale = NULL;
2389    if (smallModel_->scalingFlag()>0) {
2390      columnScale = smallModel_->columnScale();
2391    }
2392    // and do bounds in case dual needs them
2393    // may be infeasible
2394    for (i=0;i<numberColumns2;i++) {
2395      int iColumn = whichColumn[i];
2396      if (lowerBig[iColumn]>saveLowerOriginal[iColumn]) {
2397        double value = lowerBig[iColumn];
2398        lowerSmallReal[i]=value;
2399        value *= rhsScale;
2400        if (columnScale)
2401          value /= columnScale[i];
2402        lowerSmall[i]=value;
2403      }
2404      if (upperBig[iColumn]<saveUpperOriginal[iColumn]) {
2405        double value = upperBig[iColumn];
2406        upperSmallReal[i]=value;
2407        value *= rhsScale;
2408        if (columnScale)
2409          value /= columnScale[i];
2410        upperSmall[i]=value;
2411      }
2412      if (upperSmall[i]<lowerSmall[i]-1.0e-8)
2413        break;
2414    }
2415    bool infeasible = (i<numberColumns2);
2416    // Start of fast iterations
2417    bool alwaysFinish= ((specialOptions_&32)==0) ? true : false;
2418    //smallModel_->setLogLevel(1);
2419    smallModel_->setIntParam(ClpMaxNumIteration,itlim);
2420    int saveNumberFake = (static_cast<ClpSimplexDual *>(smallModel_))->numberFake_;
2421    int status;
2422    if (!infeasible) {
2423      status = static_cast<ClpSimplexDual *>(smallModel_)->fastDual(alwaysFinish);
2424    } else {
2425      status=0;
2426      smallModel_->setProblemStatus(1);
2427    }
2428    (static_cast<ClpSimplexDual *>(smallModel_))->numberFake_ = saveNumberFake;
2429    if (smallModel_->numberIterations()==-98) {
2430      printf("rrrrrrrrrrrr\n");
2431      smallModel_->checkPrimalSolution(smallModel_->solutionRegion(0),
2432                                     smallModel_->solutionRegion(1));
2433      //smallModel_->gutsOfSolution(NULL,NULL,0);
2434      //if (problemStatus==3)
2435      //smallModel_->computeObjectiveValue();
2436      printf("robj %g\n",smallModel_->objectiveValue() *
2437             modelPtr_->optimizationDirection());
2438      writeMps("rr.mps");
2439      smallModel_->writeMps("rr_small.mps");
2440      ClpSimplex temp = *smallModel_;
2441      printf("small\n");
2442      temp.setLogLevel(63);
2443      temp.dual();
2444      double limit = 0.0;
2445      modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
2446      if (temp.problemStatus()==0&&temp.objectiveValue()<limit) {
2447        printf("inf obj %g, true %g - offsets %g %g\n",smallModel_->objectiveValue(),
2448               temp.objectiveValue(),
2449               smallModel_->objectiveOffset(),temp.objectiveOffset());
2450      }
2451      printf("big\n");
2452      temp = *modelPtr_;
2453      temp.dual();
2454      if (temp.problemStatus()==0&&temp.objectiveValue()<limit) {
2455        printf("inf obj %g, true %g - offsets %g %g\n",smallModel_->objectiveValue(),
2456               temp.objectiveValue(),
2457               smallModel_->objectiveOffset(),temp.objectiveOffset());
2458      }
2459    }
2460    int problemStatus = smallModel_->problemStatus();
2461    double objectiveValue =smallModel_->objectiveValue() * modelPtr_->optimizationDirection();
2462    CoinAssert (smallModel_->problemStatus()||smallModel_->objectiveValue()<1.0e50);
2463    // make sure plausible
2464    double obj = CoinMax(objectiveValue,saveObjectiveValue);
2465    if (problemStatus==10||problemStatus<0) {
2466      // was trying to clean up or something odd
2467      if (problemStatus==10)
2468        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2469      status=1;
2470    }
2471    if (status) {
2472      // not finished - might be optimal
2473      smallModel_->checkPrimalSolution(smallModel_->solutionRegion(0),
2474                                     smallModel_->solutionRegion(1));
2475      //smallModel_->gutsOfSolution(NULL,NULL,0);
2476      //if (problemStatus==3)
2477      //smallModel_->computeObjectiveValue();
2478      objectiveValue =smallModel_->objectiveValue() *
2479        modelPtr_->optimizationDirection();
2480      if (problemStatus!=10)
2481        obj = CoinMax(objectiveValue,saveObjectiveValue);
2482      if (!smallModel_->numberDualInfeasibilities()) { 
2483        double limit = 0.0;
2484        modelPtr_->getDblParam(ClpDualObjectiveLimit, limit);
2485        if (smallModel_->secondaryStatus()==1&&!problemStatus&&obj<limit) {
2486#if 0
2487          // switch off
2488          ClpSimplex temp = *smallModel_;
2489          temp.dual();
2490          if (temp.problemStatus()==0&&temp.objectiveValue()<limit) {
2491            printf("inf obj %g, true %g - offsets %g %g\n",smallModel_->objectiveValue(),
2492                   temp.objectiveValue(),
2493                   smallModel_->objectiveOffset(),temp.objectiveOffset());
2494          }
2495          lastAlgorithm_=1;
2496          obj=limit;
2497          problemStatus=10;
2498#else
2499          obj=limit;
2500          problemStatus=3;
2501#endif
2502        }
2503        if (!smallModel_->numberPrimalInfeasibilities()&&obj<limit) { 
2504          problemStatus=0;
2505#if 0
2506          ClpSimplex temp = *smallModel_;
2507          temp.dual();
2508          if (temp.numberIterations())
2509            printf("temp iterated\n");
2510          assert (temp.problemStatus()==0&&temp.objectiveValue()<limit);
2511#endif
2512        } else if (problemStatus==10) {
2513          problemStatus=3;
2514        } else if (!smallModel_->numberPrimalInfeasibilities()) {
2515          problemStatus=1; // infeasible
2516        } 
2517      } else {
2518        // can't say much
2519        //if (problemStatus==3)
2520        //smallModel_->computeObjectiveValue();
2521        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
2522        problemStatus=3;
2523      }
2524    } else if (!problemStatus) {
2525      if (smallModel_->isDualObjectiveLimitReached()) 
2526        problemStatus=1; // infeasible
2527    }
2528    if (status&&!problemStatus) {
2529      problemStatus=3; // can't be sure
2530      lastAlgorithm_=1;
2531    }
2532    if (problemStatus<0)
2533      problemStatus=3;
2534    modelPtr_->setProblemStatus(problemStatus);
2535    modelPtr_->setObjectiveValue(obj*modelPtr_->optimizationDirection());
2536    modelPtr_->setSumDualInfeasibilities(smallModel_->sumDualInfeasibilities());
2537    modelPtr_->setNumberDualInfeasibilities(smallModel_->numberDualInfeasibilities());
2538    modelPtr_->setSumPrimalInfeasibilities(smallModel_->sumPrimalInfeasibilities());
2539    modelPtr_->setNumberPrimalInfeasibilities(smallModel_->numberPrimalInfeasibilities());
2540    double * solution = modelPtr_->primalColumnSolution();
2541    const double * solution2 = smallModel_->solutionRegion();
2542    if (!columnScale) {
2543      for (i=0;i<numberColumns2;i++) {
2544        int iColumn = whichColumn[i];
2545        solution[iColumn]= solution2[i];
2546        lowerSmallReal[i]=saveLowerOriginal[iColumn];
2547        upperSmallReal[i]=saveUpperOriginal[iColumn];
2548      }
2549    } else {
2550      for (i=0;i<numberColumns2;i++) {
2551        int iColumn = whichColumn[i];
2552        solution[iColumn]= solution2[i]*columnScale[i];
2553        lowerSmallReal[i]=saveLowerOriginal[iColumn];
2554        upperSmallReal[i]=saveUpperOriginal[iColumn];
2555      }
2556    }
2557    // could combine with loop above
2558    if (modelPtr_==smallModel_)
2559      modelPtr_->computeObjectiveValue();
2560#if 1
2561    if (status&&!problemStatus) {
2562      memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2563      modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2564      modelPtr_->checkSolutionInternal();
2565      //modelPtr_->setLogLevel(1);
2566      //modelPtr_->allSlackBasis();
2567      //modelPtr_->primal(1);
2568      //memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
2569      //modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
2570      //modelPtr_->checkSolutionInternal();
2571      assert (!modelPtr_->problemStatus());
2572    }
2573#endif
2574    modelPtr_->setNumberIterations(smallModel_->numberIterations());
2575    // and back bounds
2576    CoinMemcpyN(saveLower,number,smallModel_->lowerRegion());
2577    CoinMemcpyN(saveUpper,number,smallModel_->upperRegion());
2578  }
2579  if (savedObjective) {
2580    // fix up
2581    modelPtr_->dblParam_[ClpDualObjectiveLimit]=savedDualLimit;
2582    //modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~32));
2583    modelPtr_->objective_=savedObjective;
2584    if (!modelPtr_->problemStatus_) {
2585      CoinZeroN(modelPtr_->dual_,modelPtr_->numberRows_);
2586      CoinZeroN(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
2587      if (modelPtr_->dj_&&(modelPtr_->whatsChanged_&1)!=0)
2588        CoinZeroN(modelPtr_->dj_,modelPtr_->numberColumns_+modelPtr_->numberRows_);
2589      modelPtr_->computeObjectiveValue();
2590    }
2591  }
2592  modelPtr_->setIntParam(ClpMaxNumIteration,itlimOrig_);
2593}
2594
2595void OsiClpSolverInterface::unmarkHotStart()
2596{
2597#ifdef CLEAN_HOT_START
2598  if ((specialOptions_&65536)!=0) {
2599    modelPtr_->setLogLevel(saveData_.scalingFlag_);
2600    modelPtr_->deleteRim(0);
2601    if (lastNumberRows_<0) {
2602      specialOptions_ |= 131072;
2603      lastNumberRows_ = -1 -lastNumberRows_;
2604      if (modelPtr_->rowScale_) {
2605        if (modelPtr_->rowScale_!=rowScale_.array()) {
2606          delete [] modelPtr_->rowScale_;
2607          delete [] modelPtr_->columnScale_;
2608        }
2609        modelPtr_->rowScale_=NULL;
2610        modelPtr_->columnScale_=NULL;
2611      }
2612    }
2613    delete factorization_;
2614    delete [] spareArrays_;
2615    smallModel_=NULL;
2616    spareArrays_=NULL;
2617    factorization_=NULL;
2618    delete [] rowActivity_;
2619    delete [] columnActivity_;
2620    rowActivity_=NULL;
2621    columnActivity_=NULL;
2622    return;
2623  }
2624#endif
2625  if (smallModel_==NULL) {
2626    setWarmStart(ws_);
2627    int numberRows = modelPtr_->numberRows();
2628    int numberColumns = modelPtr_->numberColumns();
2629    CoinMemcpyN(           rowActivity_,numberRows,modelPtr_->primalRowSolution());
2630    CoinMemcpyN(columnActivity_,numberColumns,modelPtr_->primalColumnSolution());
2631    delete ws_;
2632    ws_ = NULL;
2633  } else {
2634#ifndef KEEP_SMALL
2635    if (smallModel_!=modelPtr_)
2636      delete smallModel_;
2637    smallModel_=NULL;
2638#else
2639    if (smallModel_==modelPtr_) {
2640      smallModel_=NULL;
2641    } else if (smallModel_) {
2642      if(!spareArrays_) {
2643        delete smallModel_;
2644        smallModel_=NULL;
2645        delete factorization_;
2646        factorization_=NULL;
2647      } else {
2648        static_cast<ClpSimplexDual *> (smallModel_)->cleanupAfterStrongBranching( factorization_);
2649        if ((smallModel_->specialOptions_&4096)==0) {
2650          delete factorization_;
2651        }
2652      }
2653    }
2654#endif
2655    //delete [] spareArrays_;
2656    //spareArrays_=NULL;
2657    factorization_=NULL;
2658  }
2659  delete [] rowActivity_;
2660  delete [] columnActivity_;
2661  rowActivity_=NULL;
2662  columnActivity_=NULL;
2663  // Make sure whatsChanged not out of sync
2664  if (!modelPtr_->columnUpperWork_)
2665    modelPtr_->whatsChanged_ &= ~0xffff;
2666  modelPtr_->specialOptions_ = saveData_.specialOptions_; 
2667}
2668
2669//#############################################################################
2670// Problem information methods (original data)
2671//#############################################################################
2672
2673//------------------------------------------------------------------
2674const char * OsiClpSolverInterface::getRowSense() const
2675{
2676  extractSenseRhsRange();
2677  return rowsense_;
2678}
2679//------------------------------------------------------------------
2680const double * OsiClpSolverInterface::getRightHandSide() const
2681{
2682  extractSenseRhsRange();
2683  return rhs_;
2684}
2685//------------------------------------------------------------------
2686const double * OsiClpSolverInterface::getRowRange() const
2687{
2688  extractSenseRhsRange();
2689  return rowrange_;
2690}
2691//------------------------------------------------------------------
2692// Return information on integrality
2693//------------------------------------------------------------------
2694bool OsiClpSolverInterface::isContinuous(int colNumber) const
2695{
2696  if ( integerInformation_==NULL ) return true;
2697#ifndef NDEBUG
2698  int n = modelPtr_->numberColumns();
2699  if (colNumber<0||colNumber>=n) {
2700    indexError(colNumber,"isContinuous");
2701  }
2702#endif
2703  if ( integerInformation_[colNumber]==0 ) return true;
2704  return false;
2705}
2706bool 
2707OsiClpSolverInterface::isBinary(int colNumber) const
2708{
2709#ifndef NDEBUG
2710  int n = modelPtr_->numberColumns();
2711  if (colNumber<0||colNumber>=n) {
2712    indexError(colNumber,"isBinary");
2713  }
2714#endif
2715  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) {
2716    return false;
2717  } else {
2718    const double * cu = getColUpper();
2719    const double * cl = getColLower();
2720    if ((cu[colNumber]== 1 || cu[colNumber]== 0) && 
2721        (cl[colNumber]== 0 || cl[colNumber]==1))
2722      return true;
2723    else 
2724      return false;
2725  }
2726}
2727//-----------------------------------------------------------------------------
2728bool 
2729OsiClpSolverInterface::isInteger(int colNumber) const
2730{
2731#ifndef NDEBUG
2732  int n = modelPtr_->numberColumns();
2733  if (colNumber<0||colNumber>=n) {
2734    indexError(colNumber,"isInteger");
2735  }
2736#endif
2737  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) 
2738    return false;
2739  else
2740    return true;
2741}
2742//-----------------------------------------------------------------------------
2743bool 
2744OsiClpSolverInterface::isIntegerNonBinary(int colNumber) const
2745{
2746#ifndef NDEBUG
2747  int n = modelPtr_->numberColumns();
2748  if (colNumber<0||colNumber>=n) {
2749    indexError(colNumber,"isIntegerNonBinary");
2750  }
2751#endif
2752  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) {
2753    return false;
2754  } else {
2755    return !isBinary(colNumber);
2756  }
2757}
2758//-----------------------------------------------------------------------------
2759bool 
2760OsiClpSolverInterface::isFreeBinary(int colNumber) const
2761{
2762#ifndef NDEBUG
2763  int n = modelPtr_->numberColumns();
2764  if (colNumber<0||colNumber>=n) {
2765    indexError(colNumber,"isFreeBinary");
2766  }
2767#endif
2768  if ( integerInformation_==NULL || integerInformation_[colNumber]==0 ) {
2769    return false;
2770  } else {
2771    const double * cu = getColUpper();
2772    const double * cl = getColLower();
2773    if ((cu[colNumber]== 1) && (cl[colNumber]== 0))
2774      return true;
2775    else 
2776      return false;
2777  }
2778}
2779/*  Return array of column length
2780    0 - continuous
2781    1 - binary (may get fixed later)
2782    2 - general integer (may get fixed later)
2783*/
2784const char * 
2785OsiClpSolverInterface::getColType(bool refresh) const
2786{
2787  if (!columnType_||refresh) {
2788    const int numCols = getNumCols() ;
2789    if (!columnType_)
2790      columnType_ = new char [numCols];
2791    if ( integerInformation_==NULL ) {
2792      memset(columnType_,0,numCols);
2793    } else {
2794      const double * cu = getColUpper();
2795      const double * cl = getColLower();
2796      for (int i = 0 ; i < numCols ; ++i) {
2797        if (integerInformation_[i]) {
2798          if ((cu[i]== 1 || cu[i]== 0) && 
2799              (cl[i]== 0 || cl[i]==1))
2800            columnType_[i]=1;
2801          else
2802            columnType_[i]=2;
2803        } else {
2804          columnType_[i]=0;
2805        }
2806      }
2807    }
2808  }
2809  return columnType_;
2810}
2811
2812/* Return true if column is integer but does not have to
2813   be declared as such.
2814   Note: This function returns true if the the column
2815   is binary or a general integer.
2816*/
2817bool 
2818OsiClpSolverInterface::isOptionalInteger(int colNumber) const
2819{
2820#ifndef NDEBUG
2821  int n = modelPtr_->numberColumns();
2822  if (colNumber<0||colNumber>=n) {
2823    indexError(colNumber,"isInteger");
2824  }
2825#endif
2826  if ( integerInformation_==NULL || integerInformation_[colNumber]!=2 ) 
2827    return false;
2828  else
2829    return true;
2830}
2831/* Set the index-th variable to be an optional integer variable */
2832void 
2833OsiClpSolverInterface::setOptionalInteger(int index)
2834{
2835  if (!integerInformation_) {
2836    integerInformation_ = new char[modelPtr_->numberColumns()];
2837    CoinFillN ( integerInformation_, modelPtr_->numberColumns(),static_cast<char> (0));
2838  }
2839#ifndef NDEBUG
2840  int n = modelPtr_->numberColumns();
2841  if (index<0||index>=n) {
2842    indexError(index,"setInteger");
2843  }
2844#endif
2845  integerInformation_[index]=2;
2846  modelPtr_->setInteger(index);
2847}
2848
2849//------------------------------------------------------------------
2850// Row and column copies of the matrix ...
2851//------------------------------------------------------------------
2852const CoinPackedMatrix * OsiClpSolverInterface::getMatrixByRow() const
2853{
2854  if ( matrixByRow_ == NULL ||
2855       matrixByRow_->getNumElements() != 
2856       modelPtr_->clpMatrix()->getNumElements() ) {
2857    delete matrixByRow_;
2858    matrixByRow_ = new CoinPackedMatrix();
2859    matrixByRow_->setExtraGap(0.0);
2860    matrixByRow_->setExtraMajor(0.0);
2861    matrixByRow_->reverseOrderedCopyOf(*modelPtr_->matrix());
2862    //matrixByRow_->removeGaps();
2863#if 0
2864    CoinPackedMatrix back;
2865    std::cout<<"start check"<<std::endl;
2866    back.reverseOrderedCopyOf(*matrixByRow_);
2867    modelPtr_->matrix()->isEquivalent2(back);
2868    std::cout<<"stop check"<<std::endl;
2869#endif
2870  }
2871  assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
2872  return matrixByRow_;
2873}
2874
2875const CoinPackedMatrix * OsiClpSolverInterface::getMatrixByCol() const
2876{
2877  return modelPtr_->matrix();
2878}
2879 
2880// Get pointer to mutable column-wise copy of matrix (returns NULL if not meaningful)
2881CoinPackedMatrix * 
2882OsiClpSolverInterface::getMutableMatrixByCol() const 
2883{
2884  ClpPackedMatrix * matrix = dynamic_cast<ClpPackedMatrix *>(modelPtr_->matrix_) ;
2885  if (matrix)
2886    return matrix->getPackedMatrix();
2887  else
2888    return NULL;
2889}
2890
2891//------------------------------------------------------------------
2892std::vector<double*> OsiClpSolverInterface::getDualRays(int /*maxNumRays*/,
2893                                                        bool fullRay) const
2894{
2895  if (fullRay == true) {
2896    throw CoinError("Full dual rays not yet implemented.","getDualRays",
2897                    "OsiClpSolverInterface");
2898  }
2899  return std::vector<double*>(1, modelPtr_->infeasibilityRay());
2900}
2901//------------------------------------------------------------------
2902std::vector<double*> OsiClpSolverInterface::getPrimalRays(int /*maxNumRays*/) const
2903{
2904  return std::vector<double*>(1, modelPtr_->unboundedRay());
2905}
2906//#############################################################################
2907void
2908OsiClpSolverInterface::setContinuous(int index)
2909{
2910
2911  if (integerInformation_) {
2912#ifndef NDEBUG
2913    int n = modelPtr_->numberColumns();
2914    if (index<0||index>=n) {
2915      indexError(index,"setContinuous");
2916    }
2917#endif
2918    integerInformation_[index]=0;
2919  }
2920  modelPtr_->setContinuous(index);
2921}
2922//-----------------------------------------------------------------------------
2923void
2924OsiClpSolverInterface::setInteger(int index)
2925{
2926  if (!integerInformation_) {
2927    integerInformation_ = new char[modelPtr_->numberColumns()];
2928    CoinFillN ( integerInformation_, modelPtr_->numberColumns(),static_cast<char> (0));
2929  }
2930#ifndef NDEBUG
2931  int n = modelPtr_->numberColumns();
2932  if (index<0||index>=n) {
2933    indexError(index,"setInteger");
2934  }
2935#endif
2936  integerInformation_[index]=1;
2937  modelPtr_->setInteger(index);
2938}
2939//-----------------------------------------------------------------------------
2940void
2941OsiClpSolverInterface::setContinuous(const int* indices, int len)
2942{
2943  if (integerInformation_) {
2944#ifndef NDEBUG
2945    int n = modelPtr_->numberColumns();
2946#endif
2947    int i;
2948    for (i=0; i<len;i++) {
2949      int colNumber = indices[i];
2950#ifndef NDEBUG
2951      if (colNumber<0||colNumber>=n) {
2952        indexError(colNumber,"setContinuous");
2953      }
2954#endif
2955      integerInformation_[colNumber]=0;
2956      modelPtr_->setContinuous(colNumber);
2957    }
2958  }
2959}
2960//-----------------------------------------------------------------------------
2961void
2962OsiClpSolverInterface::setInteger(const int* indices, int len)
2963{
2964  if (!integerInformation_) {
2965    integerInformation_ = new char[modelPtr_->numberColumns()];
2966    CoinFillN ( integerInformation_, modelPtr_->numberColumns(),static_cast<char> (0));
2967  }
2968#ifndef NDEBUG
2969  int n = modelPtr_->numberColumns();
2970#endif
2971  int i;
2972  for (i=0; i<len;i++) {
2973    int colNumber = indices[i];
2974#ifndef NDEBUG
2975    if (colNumber<0||colNumber>=n) {
2976      indexError(colNumber,"setInteger");
2977    }
2978#endif
2979    integerInformation_[colNumber]=1;
2980    modelPtr_->setInteger(colNumber);
2981  }
2982}
2983/* Set the objective coefficients for all columns
2984    array [getNumCols()] is an array of values for the objective.
2985    This defaults to a series of set operations and is here for speed.
2986*/
2987void 
2988OsiClpSolverInterface::setObjective(const double * array)
2989{
2990  // Say can't gurantee optimal basis etc
2991  lastAlgorithm_=999;
2992  modelPtr_->whatsChanged_ &= (0xffff&~64);
2993  int n = modelPtr_->numberColumns() ;
2994  if (fakeMinInSimplex_) {
2995    std::transform(array,array+n,
2996                   modelPtr_->objective(),std::negate<double>()) ;
2997  } else {
2998    CoinMemcpyN(array,n,modelPtr_->objective());
2999  }
3000}
3001/* Set the lower bounds for all columns
3002    array [getNumCols()] is an array of values for the objective.
3003    This defaults to a series of set operations and is here for speed.
3004*/
3005void 
3006OsiClpSolverInterface::setColLower(const double * array)
3007{
3008  // Say can't gurantee optimal basis etc
3009  lastAlgorithm_=999;
3010  modelPtr_->whatsChanged_ &= (0x1ffff&128);
3011  CoinMemcpyN(array,modelPtr_->numberColumns(),
3012                    modelPtr_->columnLower());
3013}
3014/* Set the upper bounds for all columns
3015    array [getNumCols()] is an array of values for the objective.
3016    This defaults to a series of set operations and is here for speed.
3017*/
3018void 
3019OsiClpSolverInterface::setColUpper(const double * array)
3020{
3021  // Say can't gurantee optimal basis etc
3022  lastAlgorithm_=999;
3023  modelPtr_->whatsChanged_ &= (0x1ffff&256);
3024  CoinMemcpyN(array,modelPtr_->numberColumns(),
3025                    modelPtr_->columnUpper());
3026}
3027//-----------------------------------------------------------------------------
3028void OsiClpSolverInterface::setColSolution(const double * cs) 
3029{
3030  // Say can't gurantee optimal basis etc
3031  lastAlgorithm_=999;
3032  CoinDisjointCopyN(cs,modelPtr_->numberColumns(),
3033                    modelPtr_->primalColumnSolution());
3034  if (modelPtr_->solveType()==2) {
3035    // directly into code as well
3036    CoinDisjointCopyN(cs,modelPtr_->numberColumns(),
3037                      modelPtr_->solutionRegion(1));
3038  }
3039  // compute row activity
3040  memset(modelPtr_->primalRowSolution(),0,modelPtr_->numberRows()*sizeof(double));
3041  modelPtr_->times(1.0,modelPtr_->primalColumnSolution(),modelPtr_->primalRowSolution());
3042}
3043//-----------------------------------------------------------------------------
3044void OsiClpSolverInterface::setRowPrice(const double * rs) 
3045{
3046  CoinDisjointCopyN(rs,modelPtr_->numberRows(),
3047                    modelPtr_->dualRowSolution());
3048  if (modelPtr_->solveType()==2) {
3049    // directly into code as well (? sign )
3050    CoinDisjointCopyN(rs,modelPtr_->numberRows(),
3051                      modelPtr_->djRegion(0));
3052  }
3053  // compute reduced costs
3054  memcpy(modelPtr_->dualColumnSolution(),modelPtr_->objective(),
3055         modelPtr_->numberColumns()*sizeof(double));
3056  modelPtr_->transposeTimes(-1.0,modelPtr_->dualRowSolution(),modelPtr_->dualColumnSolution());
3057}
3058
3059//#############################################################################
3060// Problem modifying methods (matrix)
3061//#############################################################################
3062void 
3063OsiClpSolverInterface::addCol(const CoinPackedVectorBase& vec,
3064                              const double collb, const double colub,   
3065                              const double obj)
3066{
3067  int numberColumns = modelPtr_->numberColumns();
3068  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3069  modelPtr_->resize(modelPtr_->numberRows(),numberColumns+1);
3070  linearObjective_ = modelPtr_->objective();
3071  basis_.resize(modelPtr_->numberRows(),numberColumns+1);
3072  setColBounds(numberColumns,collb,colub);
3073  setObjCoeff(numberColumns,obj);
3074  if (!modelPtr_->clpMatrix())
3075    modelPtr_->createEmptyMatrix();
3076  modelPtr_->matrix()->appendCol(vec);
3077  if (integerInformation_) {
3078    char * temp = new char[numberColumns+1];
3079    CoinMemcpyN(integerInformation_,numberColumns,temp);
3080    delete [] integerInformation_;
3081    integerInformation_ = temp;
3082    integerInformation_[numberColumns]=0;
3083  }
3084  freeCachedResults();
3085}
3086//-----------------------------------------------------------------------------
3087/* Add a column (primal variable) to the problem. */
3088void 
3089OsiClpSolverInterface::addCol(int numberElements, const int * rows, const double * elements,
3090                           const double collb, const double colub,   
3091                           const double obj) 
3092{
3093  CoinPackedVector column(numberElements, rows, elements);
3094  addCol(column,collb,colub,obj);
3095}
3096// Add a named column (primal variable) to the problem.
3097void 
3098OsiClpSolverInterface::addCol(const CoinPackedVectorBase& vec,
3099                      const double collb, const double colub,   
3100                      const double obj, std::string name) 
3101{
3102  int ndx = getNumCols() ;
3103  addCol(vec,collb,colub,obj) ;
3104  setColName(ndx,name) ;
3105}
3106//-----------------------------------------------------------------------------
3107void 
3108OsiClpSolverInterface::addCols(const int numcols,
3109                               const CoinPackedVectorBase * const * cols,
3110                               const double* collb, const double* colub,   
3111                               const double* obj)
3112{
3113  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3114  int numberColumns = modelPtr_->numberColumns();
3115  modelPtr_->resize(modelPtr_->numberRows(),numberColumns+numcols);
3116  linearObjective_ = modelPtr_->objective();
3117  basis_.resize(modelPtr_->numberRows(),numberColumns+numcols);
3118  double * lower = modelPtr_->columnLower()+numberColumns;
3119  double * upper = modelPtr_->columnUpper()+numberColumns;
3120  double * objective = modelPtr_->objective()+numberColumns;
3121  int iCol;
3122  if (collb) {
3123    for (iCol = 0; iCol < numcols; iCol++) {
3124      lower[iCol]= forceIntoRange(collb[iCol], -OsiClpInfinity, OsiClpInfinity);
3125      if (lower[iCol]<-1.0e27)
3126        lower[iCol]=-COIN_DBL_MAX;
3127    }
3128  } else {
3129    CoinFillN ( lower, numcols,0.0);
3130  }
3131  if (colub) {
3132    for (iCol = 0; iCol < numcols; iCol++) {
3133      upper[iCol]= forceIntoRange(colub[iCol], -OsiClpInfinity, OsiClpInfinity);
3134      if (upper[iCol]>1.0e27)
3135        upper[iCol]=COIN_DBL_MAX;
3136    }
3137  } else {
3138    CoinFillN ( upper, numcols,COIN_DBL_MAX);
3139  }
3140  if (obj) {
3141    for (iCol = 0; iCol < numcols; iCol++) {
3142      objective[iCol] = obj[iCol];
3143    }
3144  } else {
3145    CoinFillN ( objective, numcols,0.0);
3146  }
3147  if (!modelPtr_->clpMatrix())
3148    modelPtr_->createEmptyMatrix();
3149  modelPtr_->matrix()->appendCols(numcols,cols);
3150  if (integerInformation_) {
3151    char * temp = new char[numberColumns+numcols];
3152    CoinMemcpyN(integerInformation_,numberColumns,temp);
3153    delete [] integerInformation_;
3154    integerInformation_ = temp;
3155    for (iCol = 0; iCol < numcols; iCol++) 
3156      integerInformation_[numberColumns+iCol]=0;
3157  }
3158  freeCachedResults();
3159}
3160void 
3161OsiClpSolverInterface::addCols(const int numcols,
3162                               const int * columnStarts, const int * rows, const double * elements,
3163                               const double* collb, const double* colub,   
3164                               const double* obj)
3165{
3166  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3167  int numberColumns = modelPtr_->numberColumns();
3168  modelPtr_->resize(modelPtr_->numberRows(),numberColumns+numcols);
3169  linearObjective_ = modelPtr_->objective();
3170  basis_.resize(modelPtr_->numberRows(),numberColumns+numcols);
3171  double * lower = modelPtr_->columnLower()+numberColumns;
3172  double * upper = modelPtr_->columnUpper()+numberColumns;
3173  double * objective = modelPtr_->objective()+numberColumns;
3174  int iCol;
3175  if (collb) {
3176    for (iCol = 0; iCol < numcols; iCol++) {
3177      lower[iCol]= forceIntoRange(collb[iCol], -OsiClpInfinity, OsiClpInfinity);
3178      if (lower[iCol]<-1.0e27)
3179        lower[iCol]=-COIN_DBL_MAX;
3180    }
3181  } else {
3182    CoinFillN ( lower, numcols,0.0);
3183  }
3184  if (colub) {
3185    for (iCol = 0; iCol < numcols; iCol++) {
3186      upper[iCol]= forceIntoRange(colub[iCol], -OsiClpInfinity, OsiClpInfinity);
3187      if (upper[iCol]>1.0e27)
3188        upper[iCol]=COIN_DBL_MAX;
3189    }
3190  } else {
3191    CoinFillN ( upper, numcols,COIN_DBL_MAX);
3192  }
3193  if (obj) {
3194    for (iCol = 0; iCol < numcols; iCol++) {
3195      objective[iCol] = obj[iCol];
3196    }
3197  } else {
3198    CoinFillN ( objective, numcols,0.0);
3199  }
3200  if (!modelPtr_->clpMatrix())
3201    modelPtr_->createEmptyMatrix();
3202  modelPtr_->matrix()->appendCols(numcols,columnStarts,rows,elements);
3203  if (integerInformation_) {
3204    char * temp = new char[numberColumns+numcols];
3205    CoinMemcpyN(integerInformation_,numberColumns,temp);
3206    delete [] integerInformation_;
3207    integerInformation_ = temp;
3208    for (iCol = 0; iCol < numcols; iCol++) 
3209      integerInformation_[numberColumns+iCol]=0;
3210  }
3211  freeCachedResults();
3212}
3213//-----------------------------------------------------------------------------
3214void 
3215OsiClpSolverInterface::deleteCols(const int num, const int * columnIndices)
3216{
3217  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|8|64|128|256));
3218  findIntegers(false);
3219  deleteBranchingInfo(num,columnIndices);
3220  modelPtr_->deleteColumns(num,columnIndices);
3221  int nameDiscipline;
3222  getIntParam(OsiNameDiscipline,nameDiscipline) ;
3223  if (num&&nameDiscipline) {
3224    // Very clumsy (and inefficient) - need to sort and then go backwards in ? chunks
3225    int * indices = CoinCopyOfArray(columnIndices,num);
3226    std::sort(indices,indices+num);
3227    int num2=num;
3228    while(num2) {
3229      int next = indices[num2-1];
3230      int firstDelete = num2-1;
3231      int i;
3232      for (i=num2-2;i>=0;i--) {
3233        if (indices[i]+1==next) {
3234          next --;
3235          firstDelete=i;
3236        } else {
3237          break;
3238        }
3239      }
3240      OsiSolverInterface::deleteColNames(indices[firstDelete],num2-firstDelete);
3241      num2 = firstDelete;
3242      assert (num2>=0);
3243    }
3244    delete [] indices;
3245  }
3246  // synchronize integers (again)
3247  if (integerInformation_) {
3248    int numberColumns = modelPtr_->numberColumns();
3249    for (int i=0;i<numberColumns;i++) {
3250      if (modelPtr_->isInteger(i))
3251        integerInformation_[i]=1;
3252      else
3253        integerInformation_[i]=0;
3254    }
3255  }
3256  basis_.deleteColumns(num,columnIndices);
3257  linearObjective_ = modelPtr_->objective();
3258  freeCachedResults();
3259}
3260//-----------------------------------------------------------------------------
3261void 
3262OsiClpSolverInterface::addRow(const CoinPackedVectorBase& vec,
3263                              const double rowlb, const double rowub)
3264{
3265  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3266  freeCachedResults0();
3267  int numberRows = modelPtr_->numberRows();
3268  modelPtr_->resize(numberRows+1,modelPtr_->numberColumns());
3269  basis_.resize(numberRows+1,modelPtr_->numberColumns());
3270  setRowBounds(numberRows,rowlb,rowub);
3271  if (!modelPtr_->clpMatrix())
3272    modelPtr_->createEmptyMatrix();
3273  modelPtr_->matrix()->appendRow(vec);
3274  freeCachedResults1();
3275}
3276//-----------------------------------------------------------------------------
3277void OsiClpSolverInterface::addRow(const CoinPackedVectorBase& vec,
3278                                const double rowlb, const double rowub,
3279                                std::string name)
3280{
3281  int ndx = getNumRows() ;
3282  addRow(vec,rowlb,rowub) ;
3283  setRowName(ndx,name) ;
3284}
3285//-----------------------------------------------------------------------------
3286void 
3287OsiClpSolverInterface::addRow(const CoinPackedVectorBase& vec,
3288                              const char rowsen, const double rowrhs,   
3289                              const double rowrng)
3290{
3291  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3292  freeCachedResults0();
3293  int numberRows = modelPtr_->numberRows();
3294  modelPtr_->resize(numberRows+1,modelPtr_->numberColumns());
3295  basis_.resize(numberRows+1,modelPtr_->numberColumns());
3296  double rowlb = 0, rowub = 0;
3297  convertSenseToBound(rowsen, rowrhs, rowrng, rowlb, rowub);
3298  setRowBounds(numberRows,rowlb,rowub);
3299  if (!modelPtr_->clpMatrix())
3300    modelPtr_->createEmptyMatrix();
3301  modelPtr_->matrix()->appendRow(vec);
3302  freeCachedResults1();
3303}
3304//-----------------------------------------------------------------------------
3305void 
3306OsiClpSolverInterface::addRow(int numberElements, const int * columns, const double * elements,
3307                           const double rowlb, const double rowub) 
3308{
3309  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3310  freeCachedResults0();
3311  int numberRows = modelPtr_->numberRows();
3312  modelPtr_->resize(numberRows+1,modelPtr_->numberColumns());
3313  basis_.resize(numberRows+1,modelPtr_->numberColumns());
3314  setRowBounds(numberRows,rowlb,rowub);
3315  if (!modelPtr_->clpMatrix())
3316    modelPtr_->createEmptyMatrix();
3317  modelPtr_->matrix()->appendRow(numberElements, columns, elements);
3318  CoinBigIndex starts[2];
3319  starts[0]=0;
3320  starts[1]=numberElements;
3321  redoScaleFactors( 1,starts, columns, elements);
3322  freeCachedResults1();
3323}
3324//-----------------------------------------------------------------------------
3325void 
3326OsiClpSolverInterface::addRows(const int numrows,
3327                               const CoinPackedVectorBase * const * rows,
3328                               const double* rowlb, const double* rowub)
3329{
3330  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3331  freeCachedResults0();
3332  int numberRows = modelPtr_->numberRows();
3333  modelPtr_->resize(numberRows+numrows,modelPtr_->numberColumns());
3334  basis_.resize(numberRows+numrows,modelPtr_->numberColumns());
3335  double * lower = modelPtr_->rowLower()+numberRows;
3336  double * upper = modelPtr_->rowUpper()+numberRows;
3337  int iRow;
3338  for (iRow = 0; iRow < numrows; iRow++) {
3339    if (rowlb) 
3340      lower[iRow]= forceIntoRange(rowlb[iRow], -OsiClpInfinity, OsiClpInfinity);
3341    else 
3342      lower[iRow]=-OsiClpInfinity;
3343    if (rowub) 
3344      upper[iRow]= forceIntoRange(rowub[iRow], -OsiClpInfinity, OsiClpInfinity);
3345    else 
3346      upper[iRow]=OsiClpInfinity;
3347    if (lower[iRow]<-1.0e27)
3348      lower[iRow]=-COIN_DBL_MAX;
3349    if (upper[iRow]>1.0e27)
3350      upper[iRow]=COIN_DBL_MAX;
3351  }
3352  if (!modelPtr_->clpMatrix())
3353    modelPtr_->createEmptyMatrix();
3354  modelPtr_->matrix()->appendRows(numrows,rows);
3355  freeCachedResults1();
3356}
3357//-----------------------------------------------------------------------------
3358void 
3359OsiClpSolverInterface::addRows(const int numrows,
3360                               const CoinPackedVectorBase * const * rows,
3361                               const char* rowsen, const double* rowrhs,   
3362                               const double* rowrng)
3363{
3364  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3365  freeCachedResults0();
3366  int numberRows = modelPtr_->numberRows();
3367  modelPtr_->resize(numberRows+numrows,modelPtr_->numberColumns());
3368  basis_.resize(numberRows+numrows,modelPtr_->numberColumns());
3369  double * lower = modelPtr_->rowLower()+numberRows;
3370  double * upper = modelPtr_->rowUpper()+numberRows;
3371  int iRow;
3372  for (iRow = 0; iRow < numrows; iRow++) {
3373    double rowlb = 0, rowub = 0;
3374    convertSenseToBound(rowsen[iRow], rowrhs[iRow], rowrng[iRow], 
3375                        rowlb, rowub);
3376    lower[iRow]= forceIntoRange(rowlb, -OsiClpInfinity, OsiClpInfinity);
3377    upper[iRow]= forceIntoRange(rowub, -OsiClpInfinity, OsiClpInfinity);
3378    if (lower[iRow]<-1.0e27)
3379      lower[iRow]=-COIN_DBL_MAX;
3380    if (upper[iRow]>1.0e27)
3381      upper[iRow]=COIN_DBL_MAX;
3382  }
3383  if (!modelPtr_->clpMatrix())
3384    modelPtr_->createEmptyMatrix();
3385  modelPtr_->matrix()->appendRows(numrows,rows);
3386  freeCachedResults1();
3387}
3388void 
3389OsiClpSolverInterface::addRows(const int numrows,
3390                               const int * rowStarts, const int * columns, const double * element,
3391                               const double* rowlb, const double* rowub)
3392{
3393  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3394  freeCachedResults0();
3395  int numberRows = modelPtr_->numberRows();
3396  modelPtr_->resize(numberRows+numrows,modelPtr_->numberColumns());
3397  basis_.resize(numberRows+numrows,modelPtr_->numberColumns());
3398  double * lower = modelPtr_->rowLower()+numberRows;
3399  double * upper = modelPtr_->rowUpper()+numberRows;
3400  int iRow;
3401  for (iRow = 0; iRow < numrows; iRow++) {
3402    if (rowlb) 
3403      lower[iRow]= forceIntoRange(rowlb[iRow], -OsiClpInfinity, OsiClpInfinity);
3404    else 
3405      lower[iRow]=-OsiClpInfinity;
3406    if (rowub) 
3407      upper[iRow]= forceIntoRange(rowub[iRow], -OsiClpInfinity, OsiClpInfinity);
3408    else 
3409      upper[iRow]=OsiClpInfinity;
3410    if (lower[iRow]<-1.0e27)
3411      lower[iRow]=-COIN_DBL_MAX;
3412    if (upper[iRow]>1.0e27)
3413      upper[iRow]=COIN_DBL_MAX;
3414  }
3415  if (!modelPtr_->clpMatrix())
3416    modelPtr_->createEmptyMatrix();
3417  modelPtr_->matrix()->appendRows(numrows,rowStarts,columns,element);
3418  redoScaleFactors( numrows,rowStarts, columns, element);
3419  freeCachedResults1();
3420}
3421//-----------------------------------------------------------------------------
3422void 
3423OsiClpSolverInterface::deleteRows(const int num, const int * rowIndices)
3424{
3425  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
3426  // will still be optimal if all rows basic
3427  bool allBasic=true;
3428  int numBasis = basis_.getNumArtificial();
3429  for (int i=0;i<num;i++) {
3430    int iRow = rowIndices[i];
3431    if (iRow<numBasis) {
3432      if (basis_.getArtifStatus(iRow)!=CoinWarmStartBasis::basic) {
3433        allBasic=false;
3434        break;
3435      }
3436    }
3437  }
3438  int saveAlgorithm = allBasic ? lastAlgorithm_ : 999;
3439  modelPtr_->deleteRows(num,rowIndices);
3440  int nameDiscipline;
3441  getIntParam(OsiNameDiscipline,nameDiscipline) ;
3442  if (num&&nameDiscipline) {
3443    // Very clumsy (and inefficient) - need to sort and then go backwards in ? chunks
3444    int * indices = CoinCopyOfArray(rowIndices,num);
3445    std::sort(indices,indices+num);
3446    int num2=num;
3447    while(num2) {
3448      int next = indices[num2-1];
3449      int firstDelete = num2-1;
3450      int i;
3451      for (i=num2-2;i>=0;i--) {
3452        if (indices[i]+1==next) {
3453          next --;
3454          firstDelete=i;
3455        } else {
3456          break;
3457        }
3458      }
3459      OsiSolverInterface::deleteRowNames(indices[firstDelete],num2-firstDelete);
3460      num2 = firstDelete;
3461      assert (num2>=0);
3462    }
3463    delete [] indices;
3464  }
3465  basis_.deleteRows(num,rowIndices);
3466  CoinPackedMatrix * saveRowCopy = matrixByRow_;
3467  matrixByRow_=NULL;
3468  freeCachedResults();
3469  modelPtr_->setNewRowCopy(NULL);
3470  delete modelPtr_->scaledMatrix_;
3471  modelPtr_->scaledMatrix_=NULL;
3472  if (saveRowCopy) {
3473#if 1
3474    matrixByRow_=saveRowCopy;
3475    matrixByRow_->deleteRows(num,rowIndices);
3476    assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
3477#else
3478    delete saveRowCopy;
3479#endif
3480  }
3481  lastAlgorithm_ = saveAlgorithm;
3482  if ((specialOptions_&131072)!=0) 
3483    lastNumberRows_=modelPtr_->numberRows();
3484}
3485
3486//#############################################################################
3487// Methods to input a problem
3488//#############################################################################
3489
3490void
3491OsiClpSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
3492                                   const double* collb, const double* colub,   
3493                                   const double* obj,
3494                                   const double* rowlb, const double* rowub)
3495{
3496  modelPtr_->whatsChanged_ = 0;
3497  // Get rid of integer information (modelPtr will get rid of its copy)
3498  delete [] integerInformation_;
3499  integerInformation_=NULL;
3500  modelPtr_->loadProblem(matrix, collb, colub, obj, rowlb, rowub);
3501  linearObjective_ = modelPtr_->objective();
3502  freeCachedResults();
3503  basis_=CoinWarmStartBasis();
3504  if (ws_) {
3505     delete ws_;
3506     ws_ = 0;
3507  }
3508}
3509
3510//-----------------------------------------------------------------------------
3511
3512void
3513OsiClpSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
3514                                     double*& collb, double*& colub,
3515                                     double*& obj,
3516                                     double*& rowlb, double*& rowub)
3517{
3518  modelPtr_->whatsChanged_ = 0;
3519  // Get rid of integer information (modelPtr will get rid of its copy)
3520  loadProblem(*matrix, collb, colub, obj, rowlb, rowub);
3521  delete matrix;   matrix = NULL;
3522  delete[] collb;  collb = NULL;
3523  delete[] colub;  colub = NULL;
3524  delete[] obj;    obj = NULL;
3525  delete[] rowlb;  rowlb = NULL;
3526  delete[] rowub;  rowub = NULL;
3527}
3528
3529//-----------------------------------------------------------------------------
3530
3531void
3532OsiClpSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
3533                                   const double* collb, const double* colub,
3534                                   const double* obj,
3535                                   const char* rowsen, const double* rowrhs,   
3536                                   const double* rowrng)
3537{
3538  modelPtr_->whatsChanged_ = 0;
3539  // Get rid of integer information (modelPtr will get rid of its copy)
3540  // assert( rowsen != NULL );
3541  // assert( rowrhs != NULL );
3542  // If any of Rhs NULLs then create arrays
3543  int numrows = matrix.getNumRows();
3544  const char * rowsenUse = rowsen;
3545  if (!rowsen) {
3546    char * rowsen = new char [numrows];
3547    for (int i=0;i<numrows;i++)
3548      rowsen[i]='G';
3549    rowsenUse = rowsen;
3550  } 
3551  const double * rowrhsUse = rowrhs;
3552  if (!rowrhs) {
3553    double * rowrhs = new double [numrows];
3554    for (int i=0;i<numrows;i++)
3555      rowrhs[i]=0.0;
3556    rowrhsUse = rowrhs;
3557  }
3558  const double * rowrngUse = rowrng;
3559  if (!rowrng) {
3560    double * rowrng = new double [numrows];
3561    for (int i=0;i<numrows;i++)
3562      rowrng[i]=0.0;
3563    rowrngUse = rowrng;
3564  }
3565  double * rowlb = new double[numrows];
3566  double * rowub = new double[numrows];
3567  for (int i = numrows-1; i >= 0; --i) {   
3568    convertSenseToBound(rowsenUse[i],rowrhsUse[i],rowrngUse[i],rowlb[i],rowub[i]);
3569  }
3570  if (rowsen!=rowsenUse)
3571    delete [] rowsenUse;
3572  if (rowrhs!=rowrhsUse)
3573    delete [] rowrhsUse;
3574  if (rowrng!=rowrngUse)
3575    delete [] rowrngUse;
3576  loadProblem(matrix, collb, colub, obj, rowlb, rowub);
3577  delete [] rowlb;
3578  delete [] rowub;
3579}
3580
3581//-----------------------------------------------------------------------------
3582
3583void
3584OsiClpSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
3585                                     double*& collb, double*& colub,
3586                                     double*& obj,
3587                                     char*& rowsen, double*& rowrhs,
3588                                     double*& rowrng)
3589{
3590  modelPtr_->whatsChanged_ = 0;
3591  // Get rid of integer information (modelPtr will get rid of its copy)
3592  loadProblem(*matrix, collb, colub, obj, rowsen, rowrhs, rowrng);
3593  delete matrix;   matrix = NULL;
3594  delete[] collb;  collb = NULL;
3595  delete[] colub;  colub = NULL;
3596  delete[] obj;    obj = NULL;
3597  delete[] rowsen; rowsen = NULL;
3598  delete[] rowrhs; rowrhs = NULL;
3599  delete[] rowrng; rowrng = NULL;
3600}
3601
3602//-----------------------------------------------------------------------------
3603
3604void
3605OsiClpSolverInterface::loadProblem(const int numcols, const int numrows,
3606                                   const CoinBigIndex * start, const int* index,
3607                                   const double* value,
3608                                   const double* collb, const double* colub,
3609                                   const double* obj,
3610                                   const double* rowlb, const double* rowub)
3611{
3612  modelPtr_->whatsChanged_ = 0;
3613  // Get rid of integer information (modelPtr will get rid of its copy)
3614  delete [] integerInformation_;
3615  integerInformation_=NULL;
3616  modelPtr_->loadProblem(numcols, numrows, start,  index,
3617            value, collb, colub, obj,
3618            rowlb,  rowub);
3619  linearObjective_ = modelPtr_->objective();
3620  freeCachedResults();
3621  basis_=CoinWarmStartBasis();
3622  if (ws_) {
3623     delete ws_;
3624     ws_ = 0;
3625  }
3626}
3627//-----------------------------------------------------------------------------
3628
3629void
3630OsiClpSolverInterface::loadProblem(const int numcols, const int numrows,
3631                                   const CoinBigIndex * start, const int* index,
3632                                   const double* value,
3633                                   const double* collb, const double* colub,
3634                                   const double* obj,
3635                                   const char* rowsen, const double* rowrhs,   
3636                                   const double* rowrng)
3637{
3638  modelPtr_->whatsChanged_ = 0;
3639  // Get rid of integer information (modelPtr will get rid of its copy)
3640  // If any of Rhs NULLs then create arrays
3641  const char * rowsenUse = rowsen;
3642  if (!rowsen) {
3643    char * rowsen = new char [numrows];
3644    for (int i=0;i<numrows;i++)
3645      rowsen[i]='G';
3646    rowsenUse = rowsen;
3647  } 
3648  const double * rowrhsUse = rowrhs;
3649  if (!rowrhs) {
3650    double * rowrhs = new double [numrows];
3651    for (int i=0;i<numrows;i++)
3652      rowrhs[i]=0.0;
3653    rowrhsUse = rowrhs;
3654  }
3655  const double * rowrngUse = rowrng;
3656  if (!rowrng) {
3657    double * rowrng = new double [numrows];
3658    for (int i=0;i<numrows;i++)
3659      rowrng[i]=0.0;
3660    rowrngUse = rowrng;
3661  }
3662  double * rowlb = new double[numrows];
3663  double * rowub = new double[numrows];
3664  for (int i = numrows-1; i >= 0; --i) {   
3665    convertSenseToBound(rowsenUse[i],rowrhsUse[i],rowrngUse[i],rowlb[i],rowub[i]);
3666  }
3667  if (rowsen!=rowsenUse)
3668    delete [] rowsenUse;
3669  if (rowrhs!=rowrhsUse)
3670    delete [] rowrhsUse;
3671  if (rowrng!=rowrngUse)
3672    delete [] rowrngUse;
3673  loadProblem(numcols, numrows, start,  index, value, collb, colub, obj,
3674              rowlb,  rowub);
3675  delete[] rowlb;
3676  delete[] rowub;
3677}
3678// This loads a model from a coinModel object - returns number of errors
3679int 
3680OsiClpSolverInterface::loadFromCoinModel (  CoinModel & modelObject, bool keepSolution)
3681{
3682  modelPtr_->whatsChanged_ = 0;
3683  int numberErrors = 0;
3684  // Set arrays for normal use
3685  double * rowLower = modelObject.rowLowerArray();
3686  double * rowUpper = modelObject.rowUpperArray();
3687  double * columnLower = modelObject.columnLowerArray();
3688  double * columnUpper = modelObject.columnUpperArray();
3689  double * objective = modelObject.objectiveArray();
3690  int * integerType = modelObject.integerTypeArray();
3691  double * associated = modelObject.associatedArray();
3692  // If strings then do copies
3693  if (modelObject.stringsExist()) {
3694    numberErrors = modelObject.createArrays(rowLower, rowUpper, columnLower, columnUpper,
3695                                            objective, integerType,associated);
3696  }
3697  CoinPackedMatrix matrix;
3698  modelObject.createPackedMatrix(matrix,associated);
3699  int numberRows = modelObject.numberRows();
3700  int numberColumns = modelObject.numberColumns();
3701  CoinWarmStart * ws = getWarmStart();
3702  bool restoreBasis = keepSolution && numberRows&&numberRows==getNumRows()&&
3703    numberColumns==getNumCols();
3704  loadProblem(matrix, 
3705              columnLower, columnUpper, objective, rowLower, rowUpper);
3706  if (restoreBasis)
3707    setWarmStart(ws);
3708  delete ws;
3709  // Do names if wanted
3710  int numberItems;
3711  numberItems = modelObject.rowNames()->numberItems();
3712  if (numberItems) {
3713    const char *const * rowNames=modelObject.rowNames()->names();
3714    modelPtr_->copyRowNames(rowNames,0,numberItems);
3715  }
3716  numberItems = modelObject.columnNames()->numberItems();
3717  if (numberItems) {
3718    const char *const * columnNames=modelObject.columnNames()->names();
3719    modelPtr_->copyColumnNames(columnNames,0,numberItems);
3720  }
3721  // Do integers if wanted
3722  assert(integerType);
3723  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3724    if (integerType[iColumn])
3725      setInteger(iColumn);
3726  }
3727  if (rowLower!=modelObject.rowLowerArray()||
3728      columnLower!=modelObject.columnLowerArray()) {
3729    delete [] rowLower;
3730    delete [] rowUpper;
3731    delete [] columnLower;
3732    delete [] columnUpper;
3733    delete [] objective;
3734    delete [] integerType;
3735    delete [] associated;
3736    //if (numberErrors)
3737    //  handler_->message(CLP_BAD_STRING_VALUES,messages_)
3738    //    <<numberErrors
3739    //    <<CoinMessageEol;
3740  }
3741  modelPtr_->optimizationDirection_ = modelObject.optimizationDirection(); 
3742  return numberErrors;
3743}
3744
3745//-----------------------------------------------------------------------------
3746// Write mps files
3747//-----------------------------------------------------------------------------
3748
3749void OsiClpSolverInterface::writeMps(const char * filename,
3750                                     const char * extension,
3751                                     double objSense) const
3752{
3753  std::string f(filename);
3754  std::string e(extension);
3755  std::string fullname;
3756  if (e!="") {
3757    fullname = f + "." + e;
3758  } else {
3759    // no extension so no trailing period
3760    fullname = f;
3761  }
3762  // get names
3763  const char * const * const rowNames = modelPtr_->rowNamesAsChar();
3764  const char * const * const columnNames = modelPtr_->columnNamesAsChar();
3765  // Fall back on Osi version - possibly with names
3766  OsiSolverInterface::writeMpsNative(fullname.c_str(), 
3767                                     const_cast<const char **>(rowNames),
3768                                     const_cast<const char **>(columnNames),0,2,objSense,
3769                                     numberSOS_,setInfo_);
3770  if (rowNames) {
3771    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_+1);
3772    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
3773  }
3774}
3775
3776int 
3777OsiClpSolverInterface::writeMpsNative(const char *filename, 
3778                  const char ** rowNames, const char ** columnNames,
3779                  int formatType,int numberAcross,double objSense) const 
3780{
3781  return OsiSolverInterface::writeMpsNative(filename, rowNames, columnNames,
3782                               formatType, numberAcross,objSense,
3783                                            numberSOS_,setInfo_);
3784}
3785
3786//#############################################################################
3787// CLP specific public interfaces
3788//#############################################################################
3789
3790ClpSimplex * OsiClpSolverInterface::getModelPtr() const
3791{
3792  int saveAlgorithm = lastAlgorithm_;
3793  //freeCachedResults();
3794  lastAlgorithm_ = saveAlgorithm;
3795  //bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
3796  return modelPtr_;
3797}
3798
3799//-------------------------------------------------------------------
3800
3801//#############################################################################
3802// Constructors, destructors clone and assignment
3803//#############################################################################
3804//-------------------------------------------------------------------
3805// Default Constructor
3806//-------------------------------------------------------------------
3807OsiClpSolverInterface::OsiClpSolverInterface ()
3808:
3809OsiSolverInterface(),
3810rowsense_(NULL),
3811rhs_(NULL),
3812rowrange_(NULL),
3813ws_(NULL),
3814rowActivity_(NULL),
3815columnActivity_(NULL),
3816numberSOS_(0),
3817setInfo_(NULL),
3818smallModel_(NULL),
3819factorization_(NULL),
3820smallestElementInCut_(1.0e-15),
3821smallestChangeInCut_(1.0e-10),
3822largestAway_(-1.0),
3823spareArrays_(NULL),
3824matrixByRow_(NULL),
3825matrixByRowAtContinuous_(NULL), 
3826integerInformation_(NULL),
3827whichRange_(NULL),
3828fakeMinInSimplex_(false),
3829linearObjective_(NULL),
3830cleanupScaling_(0),
3831specialOptions_(0x80000000),
3832baseModel_(NULL),
3833lastNumberRows_(0),
3834continuousModel_(NULL),
3835fakeObjective_(NULL)
3836{
3837  //printf("in default %x\n",this);
3838  modelPtr_=NULL;
3839  notOwned_=false;
3840  disasterHandler_ = new OsiClpDisasterHandler();
3841  reset();
3842}
3843
3844//-------------------------------------------------------------------
3845// Clone
3846//-------------------------------------------------------------------
3847OsiSolverInterface * OsiClpSolverInterface::clone(bool CopyData) const
3848{
3849  //printf("in clone %x\n",this);
3850  OsiClpSolverInterface * newSolver;
3851  if (CopyData) {
3852    newSolver = new OsiClpSolverInterface(*this);
3853  } else {
3854    newSolver = new OsiClpSolverInterface();
3855  }
3856#if 0
3857  const double * obj = newSolver->getObjCoefficients();
3858  const double * oldObj = getObjCoefficients();
3859  if(newSolver->getNumCols()>3787)
3860    printf("%x - obj %x (from %x) val %g\n",newSolver,obj,oldObj,obj[3787]);
3861#endif
3862  return newSolver;
3863}
3864
3865
3866//-------------------------------------------------------------------
3867// Copy constructor
3868//-------------------------------------------------------------------
3869OsiClpSolverInterface::OsiClpSolverInterface (
3870                  const OsiClpSolverInterface & rhs)
3871: OsiSolverInterface(rhs),
3872rowsense_(NULL),
3873rhs_(NULL),
3874rowrange_(NULL),
3875ws_(NULL),
3876rowActivity_(NULL),
3877columnActivity_(NULL),
3878  stuff_(rhs.stuff_),
3879numberSOS_(rhs.numberSOS_),
3880setInfo_(NULL),
3881smallModel_(NULL),
3882factorization_(NULL),
3883smallestElementInCut_(rhs.smallestElementInCut_),
3884smallestChangeInCut_(rhs.smallestChangeInCut_),
3885largestAway_(-1.0),
3886spareArrays_(NULL),
3887basis_(),
3888itlimOrig_(9999999),
3889lastAlgorithm_(0),
3890notOwned_(false),
3891matrixByRow_(NULL),
3892matrixByRowAtContinuous_(NULL), 
3893integerInformation_(NULL),
3894whichRange_(NULL),
3895fakeMinInSimplex_(rhs.fakeMinInSimplex_)
3896{
3897  //printf("in copy %x - > %x\n",&rhs,this);
3898  if ( rhs.modelPtr_  ) 
3899    modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
3900  else
3901    modelPtr_ = new ClpSimplex();
3902  if ( rhs.baseModel_  ) 
3903    baseModel_ = new ClpSimplex(*rhs.baseModel_);
3904  else
3905    baseModel_ = NULL;
3906  if ( rhs.continuousModel_  ) 
3907    continuousModel_ = new ClpSimplex(*rhs.continuousModel_);
3908  else
3909    continuousModel_ = NULL;
3910  if (rhs.matrixByRowAtContinuous_)
3911    matrixByRowAtContinuous_ = new CoinPackedMatrix(*rhs.matrixByRowAtContinuous_);
3912  if ( rhs.disasterHandler_  ) 
3913    disasterHandler_ = dynamic_cast<OsiClpDisasterHandler *>(rhs.disasterHandler_->clone());
3914  else
3915    disasterHandler_ = NULL;
3916  if (rhs.fakeObjective_)
3917    fakeObjective_ = new ClpLinearObjective(*rhs.fakeObjective_);
3918  else
3919    fakeObjective_ = NULL;
3920  linearObjective_ = modelPtr_->objective();
3921  if ( rhs.ws_ ) 
3922    ws_ = new CoinWarmStartBasis(*rhs.ws_);
3923  basis_ = rhs.basis_;
3924  if (rhs.integerInformation_) {
3925    int numberColumns = modelPtr_->numberColumns();
3926    integerInformation_ = new char[numberColumns];
3927    CoinMemcpyN(rhs.integerInformation_,        numberColumns,integerInformation_);
3928  }
3929  saveData_ = rhs.saveData_;
3930  solveOptions_  = rhs.solveOptions_;
3931  cleanupScaling_ = rhs.cleanupScaling_;
3932  specialOptions_ = rhs.specialOptions_;
3933  lastNumberRows_ = rhs.lastNumberRows_;
3934  rowScale_ = rhs.rowScale_;
3935  columnScale_ = rhs.columnScale_;
3936  fillParamMaps();
3937  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3938  if (numberSOS_) {
3939    setInfo_ = new CoinSet[numberSOS_];
3940    for (int i=0;i<numberSOS_;i++)
3941      setInfo_[i]=rhs.setInfo_[i];
3942  }
3943}
3944
3945// Borrow constructor - only delete one copy
3946OsiClpSolverInterface::OsiClpSolverInterface (ClpSimplex * rhs,
3947                                              bool reallyOwn)
3948:
3949OsiSolverInterface(),
3950rowsense_(NULL),
3951rhs_(NULL),
3952rowrange_(NULL),
3953ws_(NULL),
3954rowActivity_(NULL),
3955columnActivity_(NULL),
3956numberSOS_(0),
3957setInfo_(NULL),
3958smallModel_(NULL),
3959factorization_(NULL),
3960smallestElementInCut_(1.0e-15),
3961smallestChangeInCut_(1.0e-10),
3962largestAway_(-1.0),
3963spareArrays_(NULL),
3964basis_(),
3965itlimOrig_(9999999),
3966lastAlgorithm_(0),
3967notOwned_(false),
3968matrixByRow_(NULL),
3969matrixByRowAtContinuous_(NULL), 
3970integerInformation_(NULL),
3971whichRange_(NULL),
3972fakeMinInSimplex_(false),
3973cleanupScaling_(0),
3974specialOptions_(0x80000000),
3975baseModel_(NULL),
3976lastNumberRows_(0),
3977continuousModel_(NULL),
3978fakeObjective_(NULL)
3979{
3980  disasterHandler_ = new OsiClpDisasterHandler();
3981  //printf("in borrow %x - > %x\n",&rhs,this);
3982  modelPtr_ = rhs;
3983  basis_.resize(modelPtr_->numberRows(),modelPtr_->numberColumns());
3984  linearObjective_ = modelPtr_->objective();
3985  if (rhs) {
3986    notOwned_=!reallyOwn;
3987
3988    if (rhs->integerInformation()) {
3989      int numberColumns = modelPtr_->numberColumns();
3990      integerInformation_ = new char[numberColumns];
3991      CoinMemcpyN(rhs->integerInformation(),    numberColumns,integerInformation_);
3992    }
3993  }
3994  fillParamMaps();
3995}
3996   
3997// Releases so won't error
3998void 
3999OsiClpSolverInterface::releaseClp()
4000{
4001  modelPtr_=NULL;
4002  notOwned_=false;
4003}
4004   
4005//-------------------------------------------------------------------
4006// Destructor
4007//-------------------------------------------------------------------
4008OsiClpSolverInterface::~OsiClpSolverInterface ()
4009{
4010  //printf("in destructor %x\n",this);
4011  freeCachedResults();
4012  if (!notOwned_)
4013    delete modelPtr_;
4014  delete baseModel_;
4015  delete continuousModel_;
4016  delete disasterHandler_;
4017  delete fakeObjective_;
4018  delete ws_;
4019  delete [] rowActivity_;
4020  delete [] columnActivity_;
4021  delete [] setInfo_;
4022#ifdef KEEP_SMALL
4023  if (smallModel_) {
4024    delete [] spareArrays_;
4025    spareArrays_ = NULL;
4026    delete smallModel_;
4027    smallModel_=NULL;
4028  }
4029#endif
4030  assert(smallModel_==NULL);
4031  assert(factorization_==NULL);
4032  assert(spareArrays_==NULL);
4033  delete [] integerInformation_;
4034  delete matrixByRowAtContinuous_;
4035  delete matrixByRow_;
4036}
4037
4038//-------------------------------------------------------------------
4039// Assignment operator
4040//-------------------------------------------------------------------
4041OsiClpSolverInterface &
4042OsiClpSolverInterface::operator=(const OsiClpSolverInterface& rhs)
4043{
4044  if (this != &rhs) {   
4045    //printf("in = %x - > %x\n",&rhs,this);
4046    OsiSolverInterface::operator=(rhs);
4047    freeCachedResults();
4048    if (!notOwned_)
4049      delete modelPtr_;
4050    delete ws_;
4051    if ( rhs.modelPtr_  ) 
4052      modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
4053    delete baseModel_;
4054    if ( rhs.baseModel_  ) 
4055      baseModel_ = new ClpSimplex(*rhs.baseModel_);
4056    else
4057      baseModel_ = NULL;
4058    delete continuousModel_;
4059    if ( rhs.continuousModel_  ) 
4060      continuousModel_ = new ClpSimplex(*rhs.continuousModel_);
4061    else
4062      continuousModel_ = NULL;
4063    delete matrixByRowAtContinuous_;
4064    delete matrixByRow_;
4065    matrixByRow_=NULL;
4066    if (rhs.matrixByRowAtContinuous_)
4067      matrixByRowAtContinuous_ = new CoinPackedMatrix(*rhs.matrixByRowAtContinuous_);
4068    else
4069      matrixByRowAtContinuous_=NULL;
4070    delete disasterHandler_;
4071    if ( rhs.disasterHandler_  ) 
4072      disasterHandler_ = dynamic_cast<OsiClpDisasterHandler *>(rhs.disasterHandler_->clone());
4073    else
4074      disasterHandler_ = NULL;
4075    delete fakeObjective_;
4076    if (rhs.fakeObjective_)
4077      fakeObjective_ = new ClpLinearObjective(*rhs.fakeObjective_);
4078    else
4079      fakeObjective_ = NULL;
4080    notOwned_=false;
4081    linearObjective_ = modelPtr_->objective();
4082    saveData_ = rhs.saveData_;
4083    solveOptions_  = rhs.solveOptions_;
4084    cleanupScaling_ = rhs.cleanupScaling_;
4085    specialOptions_ = rhs.specialOptions_;
4086    lastNumberRows_ = rhs.lastNumberRows_;
4087    rowScale_ = rhs.rowScale_;
4088    columnScale_ = rhs.columnScale_;
4089    basis_ = rhs.basis_;
4090    stuff_ = rhs.stuff_;
4091    if (rhs.integerInformation_) {
4092      int numberColumns = modelPtr_->numberColumns();
4093      integerInformation_ = new char[numberColumns];
4094      CoinMemcpyN(rhs.integerInformation_,      numberColumns,integerInformation_);
4095    }
4096    if ( rhs.ws_ ) 
4097      ws_ = new CoinWarmStartBasis(*rhs.ws_);
4098    else
4099      ws_=NULL;
4100    delete [] rowActivity_;
4101    delete [] columnActivity_;
4102    rowActivity_=NULL;
4103    columnActivity_=NULL;
4104    delete [] setInfo_;
4105    numberSOS_ = rhs.numberSOS_;
4106    setInfo_=NULL;
4107    if (numberSOS_) {
4108      setInfo_ = new CoinSet[numberSOS_];
4109      for (int i=0;i<numberSOS_;i++)
4110        setInfo_[i]=rhs.setInfo_[i];
4111    }
4112    assert(smallModel_==NULL);
4113    assert(factorization_==NULL);
4114    smallestElementInCut_ = rhs.smallestElementInCut_;
4115    smallestChangeInCut_ = rhs.smallestChangeInCut_;
4116    largestAway_ = -1.0;
4117    assert(spareArrays_==NULL);
4118    basis_ = rhs.basis_;
4119    fillParamMaps();
4120    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4121  }
4122  return *this;
4123}
4124
4125//#############################################################################
4126// Applying cuts
4127//#############################################################################
4128
4129void OsiClpSolverInterface::applyRowCut( const OsiRowCut & rowCut )
4130{
4131  applyRowCuts(1, &rowCut);
4132}
4133/* Apply a collection of row cuts which are all effective.
4134   applyCuts seems to do one at a time which seems inefficient.
4135*/
4136void 
4137OsiClpSolverInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
4138{
4139  if (numberCuts) {
4140    // Say can't gurantee optimal basis etc
4141    lastAlgorithm_=999;
4142
4143    // Thanks to js
4144    const OsiRowCut * * cutsp = new const OsiRowCut * [numberCuts];
4145    for (int i=0;i<numberCuts;i++) 
4146      cutsp[i] = &cuts[i];
4147   
4148    applyRowCuts(numberCuts, cutsp);
4149   
4150    delete [] cutsp;
4151  }
4152}
4153/* Apply a collection of row cuts which are all effective.
4154   applyCuts seems to do one at a time which seems inefficient.
4155*/
4156void 
4157OsiClpSolverInterface::applyRowCuts(int numberCuts, const OsiRowCut ** cuts)
4158{
4159  int i;
4160  if (!numberCuts)
4161    return;
4162  modelPtr_->whatsChanged_ &= (0xffff&~(1|2|4|16|32));
4163  CoinPackedMatrix * saveRowCopy = matrixByRow_;
4164  matrixByRow_=NULL;
4165#if 0 // was #ifndef NDEBUG
4166  int nameDiscipline;
4167  getIntParam(OsiNameDiscipline,nameDiscipline) ;
4168  assert (!nameDiscipline);
4169#endif
4170  freeCachedResults0();
4171  // Say can't gurantee optimal basis etc
4172  lastAlgorithm_=999;
4173  int numberRows = modelPtr_->numberRows();
4174  modelPtr_->resize(numberRows+numberCuts,modelPtr_->numberColumns());
4175  basis_.resize(numberRows+numberCuts,modelPtr_->numberColumns());
4176  // redo as relaxed - use modelPtr_-> addRows with starts etc
4177  int size = 0;
4178  for (i=0;i<numberCuts;i++) 
4179    size += cuts[i]->row().getNumElements();
4180  CoinBigIndex * starts = new CoinBigIndex [numberCuts+1];
4181  int * indices = new int[size];
4182  double * elements = new double[size];
4183  double * lower = modelPtr_->rowLower()+numberRows;
4184  double * upper = modelPtr_->rowUpper()+numberRows;
4185  const double * columnLower = modelPtr_->columnLower();
4186  const double * columnUpper = modelPtr_->columnUpper();
4187  size=0;
4188  for (i=0;i<numberCuts;i++) {
4189    double rowLb = cuts[i]->lb();
4190    double rowUb = cuts[i]->ub();
4191    int n=cuts[i]->row().getNumElements();
4192    const int * index = cuts[i]->row().getIndices();
4193    const double * elem = cuts[i]->row().getElements();
4194    starts[i]=size;
4195    for (int j=0;j<n;j++) {
4196      double value = elem[j];
4197      int column = index[j];
4198      if (fabs(value)>=smallestChangeInCut_) {
4199        // always take
4200        indices[size]=column;
4201        elements[size++]=value;
4202      } else if (fabs(value)>=smallestElementInCut_) {
4203        double lowerValue = columnLower[column];
4204        double upperValue = columnUpper[column];
4205        double difference = upperValue-lowerValue;
4206        if (difference<1.0e20&&difference*fabs(value)<smallestChangeInCut_&&
4207            (rowLb<-1.0e20||rowUb>1.0e20)) {
4208          // Take out and adjust to relax
4209          //printf("small el %g adjusted\n",value);
4210          if (rowLb>-1.0e20) {
4211            // just lower bound on row
4212            if (value>0.0) {
4213              // pretend at upper
4214              rowLb -= value*upperValue;
4215            } else {
4216              // pretend at lower
4217              rowLb -= value*lowerValue;
4218            }
4219          } else {
4220            // just upper bound on row
4221            if (value>0.0) {
4222              // pretend at lower
4223              rowUb -= value*lowerValue;
4224            } else {
4225              // pretend at upper
4226              rowUb -= value*upperValue;
4227            }
4228          }
4229        } else {
4230          // take (unwillingly)
4231          indices[size]=column;
4232          elements[size++]=value;
4233        }
4234      } else {
4235        //printf("small el %g ignored\n",value);
4236      }
4237    }
4238    lower[i]= forceIntoRange(rowLb, -OsiClpInfinity, OsiClpInfinity);
4239    upper[i]= forceIntoRange(rowUb, -OsiClpInfinity, OsiClpInfinity);
4240    if (lower[i]<-1.0e27)
4241      lower[i]=-COIN_DBL_MAX;
4242    if (upper[i]>1.0e27)
4243      upper[i]=COIN_DBL_MAX;
4244  }
4245  starts[numberCuts]=size;
4246  if (!modelPtr_->clpMatrix())
4247    modelPtr_->createEmptyMatrix();
4248  //modelPtr_->matrix()->appendRows(numberCuts,rows);
4249  modelPtr_->clpMatrix()->appendMatrix(numberCuts,0,starts,indices,elements);
4250  modelPtr_->setNewRowCopy(NULL);
4251  freeCachedResults1();
4252  redoScaleFactors( numberCuts,starts, indices, elements);
4253  if (saveRowCopy) {
4254#if 1
4255    matrixByRow_=saveRowCopy;
4256    matrixByRow_->appendRows(numberCuts,starts,indices,elements,0);
4257    assert (matrixByRow_->getNumElements()==modelPtr_->clpMatrix()->getNumElements());
4258#else
4259    delete saveRowCopy;
4260#endif
4261  }
4262  delete [] starts;
4263  delete [] indices;
4264  delete [] elements;
4265
4266}
4267//#############################################################################
4268// Apply Cuts
4269//#############################################################################
4270
4271OsiSolverInterface::ApplyCutsReturnCode
4272OsiClpSolverInterface::applyCuts( const OsiCuts & cs, double effectivenessLb ) 
4273{
4274  OsiSolverInterface::ApplyCutsReturnCode retVal;
4275  int i;
4276
4277  // Loop once for each column cut
4278  for ( i=0; i<cs.sizeColCuts(); i ++ ) {
4279    if ( cs.colCut(i).effectiveness() < effectivenessLb ) {
4280      retVal.incrementIneffective();
4281      continue;
4282    }
4283    if ( !cs.colCut(i).consistent() ) {
4284      retVal.incrementInternallyInconsistent();
4285      continue;
4286    }
4287    if ( !cs.colCut(i).consistent(*this) ) {
4288      retVal.incrementExternallyInconsistent();
4289      continue;
4290    }
4291    if ( cs.colCut(i).infeasible(*this) ) {
4292      retVal.incrementInfeasible();
4293      continue;
4294    }
4295    applyColCut( cs.colCut(i) );
4296    retVal.incrementApplied();
4297  }
4298
4299  // Loop once for each row cut
4300  const OsiRowCut ** addCuts = new const OsiRowCut* [cs.sizeRowCuts()];
4301  int nAdd=0;
4302  for ( i=0; i<cs.sizeRowCuts(); i ++ ) {
4303    if ( cs.rowCut(i).effectiveness() < effectivenessLb ) {
4304      retVal.incrementIneffective();
4305      continue;
4306    }
4307    if ( !cs.rowCut(i).consistent() ) {
4308      retVal.incrementInternallyInconsistent();
4309      continue;
4310    }
4311    if ( !cs.rowCut(i).consistent(*this) ) {
4312      retVal.incrementExternallyInconsistent();
4313      continue;
4314    }
4315    if ( cs.rowCut(i).infeasible(*this) ) {
4316      retVal.incrementInfeasible();
4317      continue;
4318    }
4319    addCuts[nAdd++] = cs.rowCutPtr(i);
4320    retVal.incrementApplied();
4321  }
4322  // now apply
4323  applyRowCuts(nAdd,addCuts);
4324  delete [] addCuts;
4325 
4326  return retVal;
4327}
4328// Extend scale factors
4329void 
4330OsiClpSolverInterface::redoScaleFactors(int numberAdd,const CoinBigIndex * starts,
4331                                        const int * indices, const double * elements)
4332{
4333  if ((specialOptions_&131072)!=0) {
4334    int numberRows = modelPtr_->numberRows()-numberAdd;
4335    assert (lastNumberRows_==numberRows); // ???
4336    int iRow;
4337    int newNumberRows = numberRows + numberAdd;
4338    rowScale_.extend(static_cast<int>(2*newNumberRows*sizeof(double)));
4339    double * rowScale = rowScale_.array();
4340    double * oldInverseScale = rowScale + lastNumberRows_;
4341    double * inverseRowScale = rowScale + newNumberRows;
4342    for (iRow=lastNumberRows_-1;iRow>=0;iRow--)
4343      inverseRowScale[iRow] = oldInverseScale[iRow] ;
4344    //int numberColumns = baseModel_->numberColumns();
4345    const double * columnScale = columnScale_.array();
4346    //const double * inverseColumnScale = columnScale + numberColumns;
4347    // Geometric mean on row scales
4348    // adjust arrays
4349    rowScale += lastNumberRows_;
4350    inverseRowScale += lastNumberRows_;
4351    for (iRow=0;iRow<numberAdd;iRow++) {
4352      CoinBigIndex j;
4353      double largest=1.0e-20;
4354      double smallest=1.0e50;
4355      for (j=starts[iRow];j<starts[iRow+1];j++) {
4356        int iColumn = indices[j];
4357        double value = fabs(elements[j]);
4358        // Don't bother with tiny elements
4359        if (value>1.0e-20) {
4360          value *= columnScale[iColumn];
4361          largest = CoinMax(largest,value);
4362          smallest = CoinMin(smallest,value);
4363        }
4364      }
4365      double scale=sqrt(smallest*largest);
4366      scale=CoinMax(1.0e-10,CoinMin(1.0e10,scale));
4367      inverseRowScale[iRow]=scale;
4368      rowScale[iRow]=1.0/scale;
4369    }
4370    lastNumberRows_=newNumberRows;
4371  }
4372}
4373// Delete all scale factor stuff and reset option
4374void OsiClpSolverInterface::deleteScaleFactors()
4375{
4376  delete baseModel_;
4377  baseModel_=NULL;
4378  lastNumberRows_=0;
4379  specialOptions_ &= ~131072;
4380}
4381//-----------------------------------------------------------------------------
4382
4383void OsiClpSolverInterface::applyColCut( const OsiColCut & cc )
4384{
4385  modelPtr_->whatsChanged_ &= (0x1ffff&~(128|256));
4386  // Say can't gurantee optimal basis etc
4387  lastAlgorithm_=999;
4388  double * lower = modelPtr_->columnLower();
4389  double * upper = modelPtr_->columnUpper();
4390  const CoinPackedVector & lbs = cc.lbs();
4391  const CoinPackedVector & ubs = cc.ubs();
4392  int i;
4393
4394  for ( i=0; i<lbs.getNumElements(); i++ ) {
4395    int iCol = lbs.getIndices()[i];
4396    double value = lbs.getElements()[i];
4397    if ( value > lower[iCol] )
4398      lower[iCol]= value;
4399  }
4400  for ( i=0; i<ubs.getNumElements(); i++ ) {
4401    int iCol = ubs.getIndices()[i];
4402    double value = ubs.getElements()[i];
4403    if ( value < upper[iCol] )
4404      upper[iCol]= value;
4405  }
4406}
4407//#############################################################################
4408// Private methods
4409//#############################################################################
4410
4411
4412//-------------------------------------------------------------------
4413
4414void OsiClpSolverInterface::freeCachedResults() const
4415{ 
4416  // Say can't gurantee optimal basis etc
4417  lastAlgorithm_=999;
4418  delete [] rowsense_;
4419  delete [] rhs_;
4420  delete [] rowrange_;
4421  delete matrixByRow_;
4422  if (modelPtr_&&modelPtr_->scaledMatrix_) {
4423    delete modelPtr_->scaledMatrix_;
4424    modelPtr_->scaledMatrix_=NULL;
4425  }
4426  //delete ws_;
4427  rowsense_=NULL;
4428  rhs_=NULL;
4429  rowrange_=NULL;
4430  matrixByRow_=NULL;
4431  //ws_ = NULL;
4432  if (modelPtr_&&modelPtr_->clpMatrix()) {
4433    modelPtr_->clpMatrix()->refresh(modelPtr_); // make sure all clean
4434#ifndef NDEBUG
4435    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (modelPtr_->clpMatrix());
4436    if (clpMatrix) {
4437      if (clpMatrix->getNumCols())
4438        assert (clpMatrix->getNumRows()==modelPtr_->getNumRows());
4439      if (clpMatrix->getNumRows())
4440        assert (clpMatrix->getNumCols()==modelPtr_->getNumCols());
4441    }
4442#endif
4443  }
4444}
4445
4446//-------------------------------------------------------------------
4447
4448void OsiClpSolverInterface::freeCachedResults0() const
4449{ 
4450  delete [] rowsense_;
4451  delete [] rhs_;
4452  delete [] rowrange_;
4453  rowsense_=NULL;
4454  rhs_=NULL;
4455  rowrange_=NULL;
4456}
4457
4458//-------------------------------------------------------------------
4459
4460void OsiClpSolverInterface::freeCachedResults1() const
4461{ 
4462  // Say can't gurantee optimal basis etc
4463  lastAlgorithm_=999;
4464  delete matrixByRow_;
4465  matrixByRow_=NULL;
4466  //ws_ = NULL;
4467  if (modelPtr_&&modelPtr_->clpMatrix()) {
4468    delete modelPtr_->scaledMatrix_;
4469    modelPtr_->scaledMatrix_=NULL;
4470    modelPtr_->clpMatrix()->refresh(modelPtr_); // make sure all clean
4471#ifndef NDEBUG
4472    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (modelPtr_->clpMatrix());
4473    if (clpMatrix) {
4474      assert (clpMatrix->getNumRows()==modelPtr_->getNumRows());
4475      assert (clpMatrix->getNumCols()==modelPtr_->getNumCols());
4476    }
4477#endif
4478  }
4479}
4480
4481//------------------------------------------------------------------
4482void OsiClpSolverInterface::extractSenseRhsRange() const
4483{
4484  if (rowsense_ == NULL) {
4485    // all three must be NULL
4486    assert ((rhs_ == NULL) && (rowrange_ == NULL));
4487   
4488    int nr=modelPtr_->numberRows();
4489    if ( nr!=0 ) {
4490      rowsense_ = new char[nr];
4491      rhs_ = new double[nr];
4492      rowrange_ = new double[nr];
4493      std::fill(rowrange_,rowrange_+nr,0.0);
4494     
4495      const double * lb = modelPtr_->rowLower();
4496      const double * ub = modelPtr_->rowUpper();
4497     
4498      int i;
4499      for ( i=0; i<nr; i++ ) {
4500        convertBoundToSense(lb[i], ub[i], rowsense_[i], rhs_[i], rowrange_[i]);
4501      }
4502    }
4503  }
4504}
4505// Set language
4506void 
4507OsiClpSolverInterface::newLanguage(CoinMessages::Language language)
4508{
4509  modelPtr_->newLanguage(language);
4510  OsiSolverInterface::newLanguage(language);
4511}
4512//#############################################################################
4513
4514void
4515OsiClpSolverInterface::fillParamMaps()
4516{
4517  assert (static_cast<int> (OsiMaxNumIteration)==        static_cast<int>(ClpMaxNumIteration));
4518  assert (static_cast<int> (OsiMaxNumIterationHotStart)==static_cast<int>(ClpMaxNumIterationHotStart));
4519  //assert (static_cast<int> (OsiLastIntParam)==           static_cast<int>(ClpLastIntParam));
4520 
4521  assert (static_cast<int> (OsiDualObjectiveLimit)==  static_cast<int>(ClpDualObjectiveLimit));
4522  assert (static_cast<int> (OsiPrimalObjectiveLimit)==static_cast<int>(ClpPrimalObjectiveLimit));
4523  assert (static_cast<int> (OsiDualTolerance)==       static_cast<int>(ClpDualTolerance));
4524  assert (static_cast<int> (OsiPrimalTolerance)==     static_cast<int>(ClpPrimalTolerance));
4525  assert (static_cast<int> (OsiObjOffset)==           static_cast<int>(ClpObjOffset));
4526  //assert (static_cast<int> (OsiLastDblParam)==        static_cast<int>(ClpLastDblParam));
4527 
4528  assert (static_cast<int> (OsiProbName)==    static_cast<int> (ClpProbName));
4529  //strParamMap_[OsiLastStrParam] = ClpLastStrParam;
4530}
4531// Sets up basis
4532void 
4533OsiClpSolverInterface::setBasis ( const CoinWarmStartBasis & basis)
4534{
4535  setBasis(basis,modelPtr_);
4536  setWarmStart(&basis); 
4537}
4538// Warm start
4539CoinWarmStartBasis
4540OsiClpSolverInterface::getBasis(ClpSimplex * model) const
4541{
4542  int iRow,iColumn;
4543  int numberRows = model->numberRows();
4544  int numberColumns = model->numberColumns();
4545  CoinWarmStartBasis basis;
4546  basis.setSize(numberColumns,numberRows);
4547  if (model->statusExists()) {
4548    // Flip slacks
4549    int lookupA[]={0,1,3,2,0,2};
4550    for (iRow=0;iRow<numberRows;iRow++) {
4551      int iStatus = model->getRowStatus(iRow);
4552      iStatus = lookupA[iStatus];
4553      basis.setArtifStatus(iRow,static_cast<CoinWarmStartBasis::Status> (iStatus));
4554    }
4555    int lookupS[]={0,1,2,3,0,3};
4556    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4557      int iStatus = model->getColumnStatus(iColumn);
4558      iStatus = lookupS[iStatus];
4559      basis.setStructStatus(iColumn,static_cast<CoinWarmStartBasis::Status> (iStatus));
4560    }
4561  }
4562  //basis.print();
4563  return basis;
4564}
4565// Warm start from statusArray
4566CoinWarmStartBasis * 
4567OsiClpSolverInterface::getBasis(const unsigned char * statusArray) const 
4568{
4569  int iRow,iColumn;
4570  int numberRows = modelPtr_->numberRows();
4571  int numberColumns = modelPtr_->numberColumns();
4572  CoinWarmStartBasis * basis = new CoinWarmStartBasis();
4573  basis->setSize(numberColumns,numberRows);
4574  // Flip slacks
4575  int lookupA[]={0,1,3,2,0,2};
4576  for (iRow=0;iRow<numberRows;iRow++) {
4577    int iStatus = statusArray[numberColumns+iRow]&7;
4578    iStatus = lookupA[iStatus];
4579    basis->setArtifStatus(iRow,static_cast<CoinWarmStartBasis::Status> (iStatus));
4580  }
4581  int lookupS[]={0,1,2,3,0,3};
4582  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4583    int iStatus = statusArray[iColumn]&7;
4584    iStatus = lookupS[iStatus];
4585    basis->setStructStatus(iColumn,static_cast<CoinWarmStartBasis::Status> (iStatus));
4586  }
4587  //basis->print();
4588  return basis;
4589}
4590// Sets up basis
4591void 
4592OsiClpSolverInterface::setBasis ( const CoinWarmStartBasis & basis,
4593                                  ClpSimplex * model)
4594{
4595  // Say can't gurantee optimal basis etc
4596  lastAlgorithm_=999;
4597  // transform basis to status arrays
4598  int iRow,iColumn;
4599  int numberRows = model->numberRows();
4600  int numberColumns = model->numberColumns();
4601  if (!model->statusExists()) {
4602    /*
4603      get status arrays
4604      ClpBasis would seem to have overheads and we will need
4605      extra bits anyway.
4606    */
4607    model->createStatus();
4608  }
4609  if (basis.getNumArtificial()!=numberRows||
4610      basis.getNumStructural()!=numberColumns) {
4611    CoinWarmStartBasis basis2 = basis;
4612    // resize
4613    basis2.resize(numberRows,numberColumns);
4614    // move status
4615    model->createStatus();
4616    // For rows lower and upper are flipped
4617    for (iRow=0;iRow<numberRows;iRow++) {
4618      int stat = basis2.getArtifStatus(iRow);
4619      if (stat>1)
4620        stat = 5 - stat; // so 2->3 and 3->2
4621      model->setRowStatus(iRow, static_cast<ClpSimplex::Status> (stat));
4622    }
4623    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4624      model->setColumnStatus(iColumn,
4625                             static_cast<ClpSimplex::Status> (basis2.getStructStatus(iColumn)));
4626    }
4627  } else {
4628    // move status
4629    model->createStatus();
4630    // For rows lower and upper are flipped
4631    for (iRow=0;iRow<numberRows;iRow++) {
4632      int stat = basis.getArtifStatus(iRow);
4633      if (stat>1)
4634        stat = 5 - stat; // so 2->3 and 3->2
4635      model->setRowStatus(iRow, static_cast<ClpSimplex::Status> (stat));
4636    }
4637    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4638      model->setColumnStatus(iColumn,
4639                             static_cast<ClpSimplex::Status> (basis.getStructStatus(iColumn)));
4640    }
4641  }
4642}
4643// Warm start difference from basis_ to statusArray
4644CoinWarmStartDiff * 
4645OsiClpSolverInterface::getBasisDiff(const unsigned char * statusArray) const
4646{
4647  int iRow,iColumn;
4648  int numberRows = modelPtr_->numberRows();
4649  int numberColumns = modelPtr_->numberColumns();
4650  CoinWarmStartBasis basis;
4651  basis.setSize(numberColumns,numberRows);
4652  assert (modelPtr_->statusExists());
4653  int lookupS[]={0,1,2,3,0,3};
4654  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4655    int iStatus = statusArray[iColumn]&7;
4656    iStatus = lookupS[iStatus];
4657    basis.setStructStatus(iColumn,static_cast<CoinWarmStartBasis::Status> (iStatus));
4658  }
4659  statusArray += numberColumns;
4660  // Flip slacks
4661  int lookupA[]={0,1,3,2,0,2};
4662  for (iRow=0;iRow<numberRows;iRow++) {
4663    int iStatus = statusArray[iRow]&7;
4664    iStatus = lookupA[iStatus];
4665    basis.setArtifStatus(iRow,static_cast<CoinWarmStartBasis::Status> (iStatus));
4666  }
4667  // Now basis is what we want while basis_ is old
4668  CoinWarmStartDiff * difference = basis.generateDiff(&basis_);
4669  return difference;
4670}
4671/*
4672  Read an mps file from the given filename - returns number of errors
4673  (see CoinMpsIO class)
4674*/
4675int 
4676OsiClpSolverInterface::readMps(const char *filename,
4677                               const char *extension ) 
4678{
4679  // Get rid of integer stuff
4680  delete [] integerInformation_;
4681  integerInformation_=NULL;
4682  freeCachedResults();
4683 
4684  CoinMpsIO m;
4685  m.setInfinity(getInfinity());
4686  m.passInMessageHandler(modelPtr_->messageHandler());
4687  *m.messagesPointer()=modelPtr_->coinMessages();
4688
4689  delete [] setInfo_;
4690  setInfo_=NULL;
4691  numberSOS_=0;
4692  CoinSet ** sets=NULL;
4693  // Temporarily reduce log level to get CoinMpsIO to shut up.
4694  int saveLogLevel = modelPtr_->messageHandler()->logLevel() ;
4695  modelPtr_->messageHandler()->setLogLevel(0) ;
4696  int numberErrors = m.readMps(filename,extension,numberSOS_,sets);
4697  modelPtr_->messageHandler()->setLogLevel(saveLogLevel) ;
4698  if (numberSOS_) {
4699    setInfo_ = new CoinSet[numberSOS_];
4700    for (int i=0;i<numberSOS_;i++) {
4701      setInfo_[i]=*sets[i];
4702      delete sets[i];
4703    }
4704    delete [] sets;
4705  }
4706  handler_->message(COIN_SOLVER_MPS,messages_)
4707    <<m.getProblemName()<< numberErrors <<CoinMessageEol;
4708  if (!numberErrors) {
4709
4710    // set objective function offest
4711    setDblParam(OsiObjOffset,m.objectiveOffset());
4712
4713    // set problem name
4714    setStrParam(OsiProbName,m.getProblemName());
4715   
4716    // no errors
4717    loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
4718                m.getObjCoefficients(),m.getRowSense(),m.getRightHandSide(),
4719                m.getRowRange());
4720    const char * integer = m.integerColumns();
4721    int nCols=m.getNumCols();
4722    int nRows=m.getNumRows();
4723    if (integer) {
4724      int i,n=0;
4725      int * index = new int [nCols];
4726      for (i=0;i<nCols;i++) {
4727        if (integer[i]) {
4728          index[n++]=i;
4729        }
4730      }
4731      setInteger(index,n);
4732      delete [] index;
4733      if (n) 
4734        modelPtr_->copyInIntegerInformation(integer);
4735    }
4736
4737    // set objective name
4738    setObjName(m.getObjectiveName());
4739
4740    // Always keep names
4741    int nameDiscipline;
4742    getIntParam(OsiNameDiscipline,nameDiscipline) ;
4743    int iRow;
4744    std::vector<std::string> rowNames = std::vector<std::string> ();
4745    std::vector<std::string> columnNames = std::vector<std::string> ();
4746    rowNames.reserve(nRows);
4747    for (iRow=0;iRow<nRows;iRow++) {
4748      const char * name = m.rowName(iRow);
4749      rowNames.push_back(name);
4750      if (nameDiscipline) 
4751        OsiSolverInterface::setRowName(iRow,name) ;
4752    }
4753   
4754    int iColumn;
4755    columnNames.reserve(nCols);
4756    for (iColumn=0;iColumn<nCols;iColumn++) {
4757      const char * name = m.columnName(iColumn);
4758      columnNames.push_back(name);
4759      if (nameDiscipline) 
4760        OsiSolverInterface::setColName(iColumn,name) ;
4761    }
4762    modelPtr_->copyNames(rowNames,columnNames);
4763  }
4764  return numberErrors;
4765}
4766int 
4767OsiClpSolverInterface::readMps(const char *filename, const char*extension,
4768                            int & numberSets, CoinSet ** & sets)
4769{
4770  int numberErrors = readMps(filename,extension);
4771  numberSets= numberSOS_;
4772  sets = &setInfo_;
4773  return numberErrors;
4774}
4775/* Read an mps file from the given filename returns
4776   number of errors (see OsiMpsReader class) */
4777int 
4778OsiClpSolverInterface::readMps(const char *filename,bool keepNames,bool allowErrors)
4779{
4780  // Get rid of integer stuff
4781  delete [] integerInformation_;
4782  integerInformation_=NULL;
4783  freeCachedResults();
4784 
4785  CoinMpsIO m;
4786  m.setInfinity(getInfinity());
4787  m.passInMessageHandler(modelPtr_->messageHandler());
4788  *m.messagesPointer()=modelPtr_->coinMessages();
4789  m.setSmallElementValue(CoinMax(modelPtr_->getSmallElementValue(), 
4790                                 m.getSmallElementValue()));
4791
4792  delete [] setInfo_;
4793  setInfo_=NULL;
4794  numberSOS_=0;
4795  CoinSet ** sets=NULL;
4796  int numberErrors = m.readMps(filename,"",numberSOS_,sets);
4797  if (numberSOS_) {
4798    setInfo_ = new CoinSet[numberSOS_];
4799    for (int i=0;i<numberSOS_;i++) {
4800      setInfo_[i]=*sets[i];
4801      delete sets[i];
4802    }
4803    delete [] sets;
4804  }
4805  handler_->message(COIN_SOLVER_MPS,messages_)
4806    <<m.getProblemName()<< numberErrors <<CoinMessageEol;
4807  if (!numberErrors||((numberErrors>0&&numberErrors<100000)&&allowErrors)) {
4808
4809    // set objective function offest
4810    setDblParam(OsiObjOffset,m.objectiveOffset());
4811
4812    // set problem name
4813    setStrParam(OsiProbName,m.getProblemName());
4814
4815    // set objective name
4816    setObjName(m.getObjectiveName());
4817
4818    // no errors
4819    loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
4820                m.getObjCoefficients(),m.getRowSense(),m.getRightHandSide(),
4821                m.getRowRange());
4822    int nCols=m.getNumCols();
4823    // get quadratic part
4824    if (m.reader()->whichSection (  ) == COIN_QUAD_SECTION ) {
4825      int * start=NULL;
4826      int * column = NULL;
4827      double * element = NULL;
4828      int status=m.readQuadraticMps(NULL,start,column,element,2);
4829      if (!status) 
4830        modelPtr_->loadQuadraticObjective(nCols,start,column,element);
4831      delete [] start;
4832      delete [] column;
4833      delete [] element;
4834    }
4835    const char * integer = m.integerColumns();
4836    int nRows=m.getNumRows();
4837    if (integer) {
4838      int i,n=0;
4839      int * index = new int [nCols];
4840      for (i=0;i<nCols;i++) {
4841        if (integer[i]) {
4842          index[n++]=i;
4843        }
4844      }
4845      setInteger(index,n);
4846      delete [] index;
4847      if (n) 
4848        modelPtr_->copyInIntegerInformation(integer);
4849    }
4850    if (keepNames) {
4851      // keep names
4852      int nameDiscipline;
4853      getIntParam(OsiNameDiscipline,nameDiscipline) ;
4854      int iRow;
4855      std::vector<std::string> rowNames = std::vector<std::string> ();
4856      std::vector<std::string> columnNames = std::vector<std::string> ();
4857      rowNames.reserve(nRows);
4858      for (iRow=0;iRow<nRows;iRow++) {
4859        const char * name = m.rowName(iRow);
4860        rowNames.push_back(name);
4861        if (nameDiscipline) 
4862          OsiSolverInterface::setRowName(iRow,name) ;
4863      }
4864     
4865      int iColumn;
4866      columnNames.reserve(nCols);
4867      for (iColumn=0;iColumn<nCols;iColumn++) {
4868        const char * name = m.columnName(iColumn);
4869        columnNames.push_back(name);
4870        if (nameDiscipline) 
4871          OsiSolverInterface::setColName(iColumn,name) ;
4872      }
4873      modelPtr_->copyNames(rowNames,columnNames);
4874    }
4875  }
4876  return numberErrors;
4877}
4878// Read file in LP format (with names)
4879int 
4880OsiClpSolverInterface::readLp(const char *filename, const double epsilon )
4881{
4882  CoinLpIO m;
4883  m.passInMessageHandler(modelPtr_->messageHandler());
4884  *m.messagesPointer()=modelPtr_->coinMessages();
4885  m.readLp(filename, epsilon);
4886  freeCachedResults();
4887
4888  // set objective function offest
4889  setDblParam(OsiObjOffset, 0);
4890
4891  // set problem name
4892  setStrParam(OsiProbName, m.getProblemName());
4893
4894  // set objective name
4895  setObjName(m.getObjName());
4896
4897  // no errors
4898  loadProblem(*m.getMatrixByRow(), m.getColLower(), m.getColUpper(),
4899              m.getObjCoefficients(), m.getRowLower(), m.getRowUpper());
4900
4901  const char *integer = m.integerColumns();
4902  int nCols = m.getNumCols();
4903  int nRows = m.getNumRows();
4904  if (integer) {
4905    int i, n = 0;
4906    int *index = new int [nCols];
4907    for (i=0; i<nCols; i++) {
4908      if (integer[i]) {
4909        index[n++] = i;
4910      }
4911    }
4912    setInteger(index,n);
4913    delete [] index;
4914  }
4915  // Always keep names
4916  int nameDiscipline;
4917  getIntParam(OsiNameDiscipline,nameDiscipline) ;
4918  int iRow;
4919  std::vector<std::string> rowNames = std::vector<std::string> ();
4920  std::vector<std::string> columnNames = std::vector<std::string> ();
4921  rowNames.reserve(nRows);
4922  for (iRow=0;iRow<nRows;iRow++) {
4923    const char * name = m.rowName(iRow);
4924    rowNames.push_back(name);
4925    if (nameDiscipline) 
4926      OsiSolverInterface::setRowName(iRow,name) ;
4927  }
4928 
4929  int iColumn;
4930  columnNames.reserve(nCols);
4931  for (iColumn=0;iColumn<nCols;iColumn++) {
4932    const char * name = m.columnName(iColumn);
4933    columnNames.push_back(name);
4934    if (nameDiscipline) 
4935      OsiSolverInterface::setColName(iColumn,name) ;
4936  }
4937  modelPtr_->copyNames(rowNames,columnNames);
4938  return(0);
4939}
4940/* Write the problem into an Lp file of the given filename.
4941   If objSense is non zero then -1.0 forces the code to write a
4942   maximization objective and +1.0 to write a minimization one.
4943   If 0.0 then solver can do what it wants.
4944   This version calls writeLpNative with names */
4945void 
4946OsiClpSolverInterface::writeLp(const char *filename,
4947                               const char *extension ,
4948                               double epsilon ,
4949                               int numberAcross ,
4950                               int decimals ,
4951                               double objSense ,
4952                               bool changeNameOnRange) const
4953{
4954  std::string f(filename);
4955  std::string e(extension);
4956  std::string fullname;
4957  if (e!="") {
4958    fullname = f + "." + e;
4959  } else {
4960    // no extension so no trailing period
4961    fullname = f;
4962  }
4963  // get names
4964  const char * const * const rowNames = modelPtr_->rowNamesAsChar();
4965  const char * const * const columnNames = modelPtr_->columnNamesAsChar();
4966  // Fall back on Osi version - possibly with names
4967  OsiSolverInterface::writeLpNative(fullname.c_str(), 
4968                                    rowNames,columnNames, epsilon, numberAcross,
4969                                    decimals, objSense,changeNameOnRange);
4970  if (rowNames) {
4971    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_+1);
4972    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
4973  }
4974}
4975void 
4976OsiClpSolverInterface::writeLp(FILE * fp,
4977                               double epsilon ,
4978                               int numberAcross ,
4979                               int decimals ,
4980                               double objSense ,
4981                               bool changeNameOnRange) const
4982{
4983  // get names
4984  const char * const * const rowNames = modelPtr_->rowNamesAsChar();
4985  const char * const * const columnNames = modelPtr_->columnNamesAsChar();
4986  // Fall back on Osi version - possibly with names
4987  OsiSolverInterface::writeLpNative(fp,
4988                                    rowNames,columnNames, epsilon, numberAcross,
4989                                    decimals, objSense,changeNameOnRange);
4990  if (rowNames) {
4991    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_+1);
4992    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
4993  }
4994}
4995/*
4996  I (JJF) am getting incredibly annoyed because I can't just replace a matrix.
4997  The default behavior of this is do nothing so only use where that would not matter
4998  e.g. strengthening a matrix for MIP
4999*/
5000void 
5001OsiClpSolverInterface::replaceMatrixOptional(const CoinPackedMatrix & matrix)
5002{
5003  modelPtr_->whatsChanged_ &= (0xffff&~(2|4|8));
5004  replaceMatrix(matrix);
5005}
5006// And if it does matter (not used at present)
5007void 
5008OsiClpSolverInterface::replaceMatrix(const CoinPackedMatrix & matrix)
5009{
5010  modelPtr_->whatsChanged_ &= (0xffff&~(2|4|8));
5011  delete modelPtr_->matrix_;
5012  delete modelPtr_->rowCopy_;
5013  modelPtr_->rowCopy_=NULL;
5014  if (matrix.isColOrdered()) {
5015    modelPtr_->matrix_=new ClpPackedMatrix(matrix);
5016  } else {
5017    CoinPackedMatrix matrix2;
5018    matrix2.setExtraGap(0.0);
5019    matrix2.setExtraMajor(0.0);
5020    matrix2.reverseOrderedCopyOf(matrix);
5021    modelPtr_->matrix_=new ClpPackedMatrix(matrix2);
5022  }   
5023  modelPtr_->matrix_->setDimensions(modelPtr_->numberRows_,modelPtr_->numberColumns_);
5024  freeCachedResults();
5025}
5026// Get pointer to array[getNumCols()] of primal solution vector
5027const double * 
5028OsiClpSolverInterface::getColSolution() const 
5029{ 
5030  if (modelPtr_->solveType()!=2) {
5031    return modelPtr_->primalColumnSolution();
5032  } else {
5033    // simplex interface
5034    return modelPtr_->solutionRegion(1);
5035  }
5036}
5037 
5038// Get pointer to array[getNumRows()] of dual prices
5039const double * 
5040OsiClpSolverInterface::getRowPrice() const
5041{ 
5042  if (modelPtr_->solveType()!=2) {
5043    return modelPtr_->dualRowSolution();
5044  } else {
5045    // simplex interface
5046    //return modelPtr_->djRegion(0);
5047    return modelPtr_->dualRowSolution();
5048  }
5049}
5050 
5051// Get a pointer to array[getNumCols()] of reduced costs
5052const double * 
5053OsiClpSolverInterface::getReducedCost() const 
5054{ 
5055  if (modelPtr_->solveType()!=2) {
5056    return modelPtr_->dualColumnSolution();
5057  } else {
5058    // simplex interface
5059    return modelPtr_->djRegion(1);
5060  }
5061}
5062
5063/* Get pointer to array[getNumRows()] of row activity levels (constraint
5064   matrix times the solution vector */
5065const double * 
5066OsiClpSolverInterface::getRowActivity() const 
5067{ 
5068  if (modelPtr_->solveType()!=2) {
5069    return modelPtr_->primalRowSolution();
5070  } else {
5071    // simplex interface
5072    return modelPtr_->solutionRegion(0);
5073  }
5074}
5075double 
5076OsiClpSolverInterface::getObjValue() const 
5077{
5078  if (modelPtr_->numberIterations()||modelPtr_->upperIn_!=-COIN_DBL_MAX) {
5079    // This does not pass unitTest when getObjValue is called before solve.
5080    //printf("obj a %g %g\n",modelPtr_->objectiveValue(),
5081    //     OsiSolverInterface::getObjValue());
5082    if (fakeMinInSimplex_)
5083      return -modelPtr_->objectiveValue() ;
5084    else
5085      return modelPtr_->objectiveValue();
5086  } else {
5087    return OsiSolverInterface::getObjValue();
5088  }
5089}
5090
5091/* Set an objective function coefficient */
5092void 
5093OsiClpSolverInterface::setObjCoeff( int elementIndex, double elementValue )
5094{
5095  modelPtr_->whatsChanged_ &= 0xffff;
5096  // Say can't gurantee optimal basis etc
5097  lastAlgorithm_=999;
5098#ifndef NDEBUG
5099  int n = modelPtr_->numberColumns();
5100  if (elementIndex<0||elementIndex>=n) {
5101    indexError(elementIndex,"setObjCoeff");
5102  }
5103#endif
5104  modelPtr_->setObjectiveCoefficient(elementIndex,
5105    ((fakeMinInSimplex_)?-elementValue:elementValue));
5106}
5107
5108/* Set a single column lower bound<br>
5109   Use -DBL_MAX for -infinity. */
5110void 
5111OsiClpSolverInterface::setColLower( int index, double elementValue )
5112{
5113  modelPtr_->whatsChanged_ &= 0x1ffff;
5114#ifndef NDEBUG
5115  int n = modelPtr_->numberColumns();
5116  if (index<0||index>=n) {
5117    indexError(index,"setColLower");
5118  }
5119#endif
5120  double currentValue = modelPtr_->columnActivity_[index];
5121  bool changed=(currentValue<elementValue-modelPtr_->primalTolerance()||
5122                index>=basis_.getNumStructural()||
5123                basis_.getStructStatus(index)==CoinWarmStartBasis::atLowerBound);
5124  // Say can't gurantee optimal basis etc
5125  if (changed)
5126    lastAlgorithm_=999;
5127  if (!modelPtr_->lower_)
5128    modelPtr_->whatsChanged_ &= ~0xffff; // switch off
5129  modelPtr_->setColumnLower(index,elementValue);
5130}
5131     
5132/* Set a single column upper bound<br>
5133   Use DBL_MAX for infinity. */
5134void 
5135OsiClpSolverInterface::setColUpper( int index, double elementValue )
5136{
5137  modelPtr_->whatsChanged_ &= 0x1ffff;
5138#ifndef NDEBUG
5139  int n = modelPtr_->numberColumns();
5140  if (index<0||index>=n) {
5141    indexError(index,"setColUpper");
5142  }
5143#endif
5144  double currentValue = modelPtr_->columnActivity_[index];
5145  bool changed=(currentValue>elementValue+modelPtr_->primalTolerance()||
5146                index>=basis_.getNumStructural()||
5147                basis_.getStructStatus(index)==CoinWarmStartBasis::atUpperBound);
5148  // Say can't gurantee optimal basis etc
5149  if (changed)
5150    lastAlgorithm_=999;
5151  if (!modelPtr_->upper_)
5152    modelPtr_->whatsChanged_ &= ~0xffff; // switch off
5153  modelPtr_->setColumnUpper(index,elementValue);
5154}
5155
5156/* Set a single column lower and upper bound */
5157void 
5158OsiClpSolverInterface::setColBounds( int elementIndex,
5159                                     double lower, double upper )
5160{
5161  modelPtr_->whatsChanged_ &= 0x1ffff;
5162  // Say can't gurantee optimal basis etc
5163  lastAlgorithm_=999;
5164#ifndef NDEBUG
5165  int n = modelPtr_->numberColumns();
5166  if (elementIndex<0||elementIndex>=n) {
5167    indexError(elementIndex,"setColBounds");
5168  }
5169#endif
5170  if (!modelPtr_->lower_)
5171    modelPtr_->whatsChanged_ &= ~0xffff; // switch off
5172  modelPtr_->setColumnBounds(elementIndex,lower,upper);
5173}
5174void OsiClpSolverInterface::setColSetBounds(const int* indexFirst,
5175                                            const int* indexLast,
5176                                            const double* boundList)
5177{
5178  modelPtr_->whatsChanged_ &= 0x1ffff;
5179  // Say can't gurantee optimal basis etc
5180  lastAlgorithm_=999;
5181#ifndef NDEBUG
5182  int n = modelPtr_->numberColumns();
5183  const int * indexFirst2=indexFirst;
5184  while (indexFirst2 != indexLast) {
5185    const int iColumn=*indexFirst2++;
5186    if (iColumn<0||iColumn>=n) {
5187      indexError(iColumn,"setColSetBounds");
5188    }
5189  }
5190#endif
5191  modelPtr_->setColSetBounds(indexFirst,indexLast,boundList);
5192}
5193//------------------------------------------------------------------
5194/* Set a single row lower bound<br>
5195   Use -DBL_MAX for -infinity. */
5196void 
5197OsiClpSolverInterface::setRowLower( int elementIndex, double elementValue ) {
5198  // Say can't gurantee optimal basis etc
5199  lastAlgorithm_=999;
5200  modelPtr_->whatsChanged_ &= 0xffff;
5201#ifndef NDEBUG
5202  int n = modelPtr_->numberRows();
5203  if (elementIndex<0||elementIndex>=n) {
5204    indexError(elementIndex,"setRowLower");
5205  }
5206#endif
5207  modelPtr_->setRowLower(elementIndex , elementValue);
5208  if (rowsense_!=NULL) {
5209    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5210    convertBoundToSense(modelPtr_->rowLower_[elementIndex], 
5211                        modelPtr_->rowUpper_[elementIndex],
5212                        rowsense_[elementIndex], rhs_[elementIndex], rowrange_[elementIndex]);
5213  }
5214}
5215     
5216/* Set a single row upper bound<br>
5217   Use DBL_MAX for infinity. */
5218void 
5219OsiClpSolverInterface::setRowUpper( int elementIndex, double elementValue ) {
5220  modelPtr_->whatsChanged_ &= 0xffff;
5221  // Say can't gurantee optimal basis etc
5222  lastAlgorithm_=999;
5223#ifndef NDEBUG
5224  int n = modelPtr_->numberRows();
5225  if (elementIndex<0||elementIndex>=n) {
5226    indexError(elementIndex,"setRowUpper");
5227  }
5228#endif
5229  modelPtr_->setRowUpper(elementIndex , elementValue);
5230  if (rowsense_!=NULL) {
5231    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5232    convertBoundToSense(modelPtr_->rowLower_[elementIndex], 
5233                        modelPtr_->rowUpper_[elementIndex],
5234                        rowsense_[elementIndex], rhs_[elementIndex], rowrange_[elementIndex]);
5235  }
5236}
5237   
5238/* Set a single row lower and upper bound */
5239void 
5240OsiClpSolverInterface::setRowBounds( int elementIndex,
5241              double lower, double upper ) {
5242  modelPtr_->whatsChanged_ &= 0xffff;
5243  // Say can't gurantee optimal basis etc
5244  lastAlgorithm_=999;
5245#ifndef NDEBUG
5246  int n = modelPtr_->numberRows();
5247  if (elementIndex<0||elementIndex>=n) {
5248    indexError(elementIndex,"setRowBounds");
5249  }
5250#endif
5251  modelPtr_->setRowBounds(elementIndex,lower,upper);
5252  if (rowsense_!=NULL) {
5253    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5254    convertBoundToSense(modelPtr_->rowLower_[elementIndex], 
5255                        modelPtr_->rowUpper_[elementIndex],
5256                        rowsense_[elementIndex], rhs_[elementIndex], rowrange_[elementIndex]);
5257  }
5258}
5259//-----------------------------------------------------------------------------
5260void
5261OsiClpSolverInterface::setRowType(int i, char sense, double rightHandSide,
5262                                  double range)
5263{
5264  modelPtr_->whatsChanged_ &= 0xffff;
5265  // Say can't gurantee optimal basis etc
5266  lastAlgorithm_=999;
5267#ifndef NDEBUG
5268  int n = modelPtr_->numberRows();
5269  if (i<0||i>=n) {
5270    indexError(i,"setRowType");
5271  }
5272#endif
5273  double lower = 0, upper = 0;
5274  convertSenseToBound(sense, rightHandSide, range, lower, upper);
5275  setRowBounds(i, lower, upper);
5276  // If user is using sense then set
5277  if (rowsense_) {
5278    rowsense_[i] = sense;
5279    rhs_[i] = rightHandSide;
5280    rowrange_[i] = range;
5281  }
5282}
5283// Set name of row
5284void 
5285//OsiClpSolverInterface::setRowName(int rowIndex, std::string & name)
5286OsiClpSolverInterface::setRowName(int rowIndex, std::string name) 
5287{
5288  if (rowIndex>=0&&rowIndex<modelPtr_->numberRows()) {
5289    int nameDiscipline;
5290    getIntParam(OsiNameDiscipline,nameDiscipline) ;
5291    if (nameDiscipline) {
5292      modelPtr_->setRowName(rowIndex,name);
5293      OsiSolverInterface::setRowName(rowIndex,name) ;
5294    }
5295  }
5296}
5297// Return name of row if one exists or Rnnnnnnn
5298// we ignore maxLen
5299std::string
5300OsiClpSolverInterface::getRowName(int rowIndex, unsigned int /*maxLen*/) const
5301{ 
5302        if (rowIndex == getNumRows())
5303                return getObjName();
5304  return modelPtr_->getRowName(rowIndex);
5305}
5306   
5307// Set name of col
5308void 
5309//OsiClpSolverInterface::setColName(int colIndex, std::string & name)
5310OsiClpSolverInterface::setColName(int colIndex, std::string name) 
5311{
5312  if (colIndex>=0&&colIndex<modelPtr_->numberColumns()) {
5313    int nameDiscipline;
5314    getIntParam(OsiNameDiscipline,nameDiscipline) ;
5315    if (nameDiscipline) {
5316      modelPtr_->setColumnName(colIndex,name);
5317      OsiSolverInterface::setColName(colIndex,name) ;
5318    }
5319  }
5320}
5321// Return name of col if one exists or Rnnnnnnn
5322std::string
5323OsiClpSolverInterface::getColName(int colIndex, unsigned int /*maxLen*/) const
5324{
5325  return modelPtr_->getColumnName(colIndex);
5326}
5327   
5328   
5329//-----------------------------------------------------------------------------
5330void OsiClpSolverInterface::setRowSetBounds(const int* indexFirst,
5331                                            const int* indexLast,
5332                                            const double* boundList)
5333{
5334  modelPtr_->whatsChanged_ &= 0xffff;
5335  // Say can't gurantee optimal basis etc
5336  lastAlgorithm_=999;
5337#ifndef NDEBUG
5338  int n = modelPtr_->numberRows();
5339  const int * indexFirst2=indexFirst;
5340  while (indexFirst2 != indexLast) {
5341    const int iColumn=*indexFirst2++;
5342    if (iColumn<0||iColumn>=n) {
5343      indexError(iColumn,"setColumnSetBounds");
5344    }
5345  }
5346#endif
5347  modelPtr_->setRowSetBounds(indexFirst,indexLast,boundList);
5348  if (rowsense_ != NULL) {
5349    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5350    double * lower = modelPtr_->rowLower();
5351    double * upper = modelPtr_->rowUpper();
5352    while (indexFirst != indexLast) {
5353      const int iRow=*indexFirst++;
5354      convertBoundToSense(lower[iRow], upper[iRow],
5355                          rowsense_[iRow], rhs_[iRow], rowrange_[iRow]);
5356    }
5357  }
5358}
5359//-----------------------------------------------------------------------------
5360void
5361OsiClpSolverInterface::setRowSetTypes(const int* indexFirst,
5362                                      const int* indexLast,
5363                                      const char* senseList,
5364                                      const double* rhsList,
5365                                      const double* rangeList)
5366{
5367  modelPtr_->whatsChanged_ &= 0xffff;
5368  // Say can't gurantee optimal basis etc
5369  lastAlgorithm_=999;
5370#ifndef NDEBUG
5371  int n = modelPtr_->numberRows();
5372#endif
5373  const int len = static_cast<int>(indexLast - indexFirst);
5374  while (indexFirst != indexLast) {
5375    const int iRow= *indexFirst++;
5376#ifndef NDEBUG
5377    if (iRow<0||iRow>=n) {
5378      indexError(iRow,"isContinuous");
5379    }
5380#endif
5381    double lowerValue = 0;
5382    double upperValue = 0;
5383    if (rangeList){
5384      convertSenseToBound(*senseList++, *rhsList++, *rangeList++,
5385                          lowerValue, upperValue);
5386    } else {
5387      convertSenseToBound(*senseList++, *rhsList++, 0,
5388                          lowerValue, upperValue);
5389    }
5390    modelPtr_->setRowBounds(iRow,lowerValue,upperValue);
5391  }
5392  if (rowsense_ != NULL) {
5393    assert ((rhs_ != NULL) && (rowrange_ != NULL));
5394    indexFirst -= len;
5395    senseList -= len;
5396    rhsList -= len;
5397    if (rangeList)
5398       rangeList -= len;
5399    while (indexFirst != indexLast) {
5400      const int iRow=*indexFirst++;
5401      rowsense_[iRow] = *senseList++;
5402      rhs_[iRow] = *rhsList++;
5403      if (rangeList)
5404         rowrange_[iRow] = *rangeList++;
5405    }
5406  }
5407}
5408
5409/*
5410  Clp's copy-in/copy-out design paradigm is a challenge for the simplex modes.
5411  Normal operation goes like this:
5412    * startup() loads clp's work arrays, performing scaling for numerical
5413      stability and compensating for max.
5414    * clp solves the problem
5415    * finish() unloads the work arrays into answer arrays, undoing scaling
5416      and max compensation.
5417  There are two solutions: undo scaling and max on demand, or make them
5418  into noops. The various getBInv* methods undo scaling on demand (but
5419  see special option 512) and do not need to worry about max. Other get
5420  solution methods are not coded to do this, so the second approach is
5421  used. For simplex modes, turn off scaling (necessary for both primal and
5422  dual solutions) and temporarily convert max to min (necessary for dual
5423  solution). This makes the unscaling in getBInv* superfluous, but don't
5424  remove it. Arguably the better solution here would be to go through and
5425  add unscaling and max compensation to the get solution methods. Look for
5426  fakeMinInSimplex to see the places this propagates to.
5427
5428  TODO: setRowPrice never has worked properly, and I didn't try to fix it in
5429        this go-round.
5430
5431  As of 100907, change applied to [enable|disable]Factorization (mode 1).
5432  Limitation of [enable|disable]SimplexInterface (mode 2) noted in
5433  documentation.  -- lh, 100907 --
5434*/
5435/*
5436  Enables normal operation of subsequent functions.  This method is supposed
5437  to ensure that all typical things (like reduced costs, etc.) are updated
5438  when individual pivots are executed and can be queried by other methods
5439*/
5440void 
5441OsiClpSolverInterface::enableSimplexInterface(bool doingPrimal)
5442{
5443  modelPtr_->whatsChanged_ &= 0xffff;
5444  if (modelPtr_->solveType()==2)
5445    return;
5446  assert (modelPtr_->solveType()==1);
5447  int saveIts = modelPtr_->numberIterations_;
5448  modelPtr_->setSolveType(2);
5449  if (doingPrimal)
5450    modelPtr_->setAlgorithm(1);
5451  else
5452    modelPtr_->setAlgorithm(-1);
5453  // Do initialization
5454  saveData_ = modelPtr_->saveData();
5455  saveData_.scalingFlag_=modelPtr_->scalingFlag();
5456  modelPtr_->scaling(0);
5457  specialOptions_ = 0x80000000;
5458  // set infeasibility cost up
5459  modelPtr_->setInfeasibilityCost(1.0e12);
5460  ClpDualRowDantzig dantzig;
5461  modelPtr_->setDualRowPivotAlgorithm(dantzig);
5462  ClpPrimalColumnDantzig dantzigP;
5463  dantzigP.saveWeights(modelPtr_,0); // set modelPtr
5464  modelPtr_->setPrimalColumnPivotAlgorithm(dantzigP);
5465  int saveOptions = modelPtr_->specialOptions_;
5466  modelPtr_->specialOptions_ &= ~262144;
5467  delete modelPtr_->scaledMatrix_;
5468  modelPtr_->scaledMatrix_=NULL;
5469#ifdef NDEBUG
5470  modelPtr_->startup(0);
5471#else
5472  int returnCode=modelPtr_->startup(0);
5473  assert (!returnCode||returnCode==2);
5474#endif
5475  modelPtr_->specialOptions_=saveOptions;
5476  modelPtr_->numberIterations_=saveIts;
5477}
5478
5479//Undo whatever setting changes the above method had to make
5480void 
5481OsiClpSolverInterface::disableSimplexInterface()
5482{
5483  modelPtr_->whatsChanged_ &= 0xffff;
5484  assert (modelPtr_->solveType()==2);
5485  // declare optimality anyway  (for message handler)
5486  modelPtr_->setProblemStatus(0);
5487  modelPtr_->setSolveType(1);
5488  // message will not appear anyway
5489  int saveMessageLevel=modelPtr_->messageHandler()->logLevel();
5490  modelPtr_->messageHandler()->setLogLevel(0);
5491  modelPtr_->finish();
5492  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
5493  modelPtr_->restoreData(saveData_);
5494  modelPtr_->scaling(saveData_.scalingFlag_);
5495  ClpDualRowSteepest steepest;
5496  modelPtr_->setDualRowPivotAlgorithm(steepest);
5497  ClpPrimalColumnSteepest steepestP;
5498  modelPtr_->setPrimalColumnPivotAlgorithm(steepestP);
5499  basis_ = getBasis(modelPtr_);
5500  modelPtr_->setSolveType(1);
5501}
5502
5503/*
5504  Force scaling off. If the client thinks we're maximising, arrange it so
5505  that clp sees minimisation while the client still sees maximisation. In
5506  keeping with the spirit of the getBInv methods, special option 512 will
5507  leave all work to the client.
5508*/
5509void 
5510OsiClpSolverInterface::enableFactorization() const
5511{
5512  saveData_.specialOptions_=specialOptions_;
5513  // Try to preserve work regions, reuse factorization
5514  if ((specialOptions_&(1+8))!=1+8)
5515    setSpecialOptionsMutable((1+8)|specialOptions_);
5516  // Are we allowed to make the output sensible to mere mortals?
5517  if ((specialOptions_&512)==0) {
5518    // Force scaling to off
5519    saveData_.scalingFlag_ = modelPtr_->scalingFlag() ;
5520    modelPtr_->scaling(0) ;
5521    // Temporarily force to min but keep a copy of original objective.
5522    if (getObjSense() < 0.0) {
5523      fakeMinInSimplex_ = true ;
5524      modelPtr_->setOptimizationDirection(1.0) ;
5525      double *c = modelPtr_->objective() ;
5526      int n = getNumCols() ;
5527      linearObjective_ = new double[n] ;
5528      CoinMemcpyN(c,n,linearObjective_) ;
5529      std::transform(c,c+n,c,std::negate<double>()) ;
5530    }
5531  }
5532  int saveStatus = modelPtr_->problemStatus_;
5533#ifdef NDEBUG
5534  modelPtr_->startup(0);
5535#else
5536  int returnCode=modelPtr_->startup(0);
5537  assert (!returnCode||returnCode==2);
5538#endif
5539  modelPtr_->problemStatus_=saveStatus;
5540}
5541
5542/*
5543  Undo enableFactorization. Retrieve the special options and scaling and
5544  remove the temporary objective used to fake minimisation in clp.
5545*/
5546void 
5547OsiClpSolverInterface::disableFactorization() const
5548{
5549  specialOptions_=saveData_.specialOptions_;
5550  // declare optimality anyway  (for message handler)
5551  modelPtr_->setProblemStatus(0);
5552  // message will not appear anyway
5553  int saveMessageLevel=modelPtr_->messageHandler()->logLevel();
5554  modelPtr_->messageHandler()->setLogLevel(0);
5555  modelPtr_->finish();
5556  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
5557  // Client asked for transforms on the way in, so back out.
5558  if ((specialOptions_&512)==0) {
5559    modelPtr_->scaling(saveData_.scalingFlag_) ;
5560    if (fakeMinInSimplex_ == true) {
5561      fakeMinInSimplex_ = false ;
5562      modelPtr_->setOptimizationDirection(-1.0) ;
5563      double *c = modelPtr_->objective() ;
5564      int n = getNumCols() ;
5565      std::transform(c,c+n,c,std::negate<double>()) ;
5566      delete[] linearObjective_ ;
5567    }
5568  }
5569}
5570
5571
5572
5573/* The following two methods may be replaced by the
5574   methods of OsiSolverInterface using OsiWarmStartBasis if:
5575   1. OsiWarmStartBasis resize operation is implemented
5576   more efficiently and
5577   2. It is ensured that effects on the solver are the same
5578   
5579   Returns a basis status of the structural/artificial variables
5580*/
5581void 
5582OsiClpSolverInterface::getBasisStatus(int* cstat, int* rstat) const
5583{
5584  int iRow,iColumn;
5585  int numberRows = modelPtr_->numberRows();
5586  int numberColumns = modelPtr_->numberColumns();
5587  const double * pi = modelPtr_->dualRowSolution();
5588  const double * dj = modelPtr_->dualColumnSolution();
5589  double multiplier = modelPtr_->optimizationDirection();
5590  // Flip slacks
5591  int lookupA[]={0,1,3,2,0,3};
5592  for (iRow=0;iRow<numberRows;iRow++) {
5593    int iStatus = modelPtr_->getRowStatus(iRow);
5594    if (iStatus==5) {
5595      // Fixed - look at reduced cost
5596      if (pi[iRow]*multiplier>1.0e-7)
5597        iStatus = 3;
5598    }
5599    iStatus = lookupA[iStatus];
5600    rstat[iRow]=iStatus;
5601  }
5602  int lookupS[]={0,1,2,3,0,3};
5603  for (iColumn=0;iColumn<numberColumns;iColumn++) {
5604    int iStatus = modelPtr_->getColumnStatus(iColumn);
5605    if (iStatus==5) {
5606      // Fixed - look at reduced cost
5607      if (dj[iColumn]*multiplier<-1.0e-7)
5608        iStatus = 2;
5609    }
5610    iStatus = lookupS[iStatus];
5611    cstat[iColumn]=iStatus;
5612  }
5613}
5614
5615//Set the status of structural/artificial variables
5616//Returns 0 if OK, 1 if problem is bad e.g. duplicate elements, too large ...
5617int 
5618OsiClpSolverInterface::setBasisStatus(const int* cstat, const int* rstat)
5619{
5620  modelPtr_->whatsChanged_ &= 0xffff;
5621  // Say can't gurantee optimal basis etc
5622  lastAlgorithm_=999;
5623  modelPtr_->createStatus();
5624  int i, n;
5625  double * lower, * upper, * solution;
5626  n=modelPtr_->numberRows();
5627  lower = modelPtr_->rowLower();
5628  upper = modelPtr_->rowUpper();
5629  solution = modelPtr_->primalRowSolution();
5630  // For rows lower and upper are flipped
5631  int lookupA[]={0,1,3,2};
5632  for (i=0;i<n;i++) {
5633    int status = lookupA[rstat[i]];
5634    if (status<0||status>3)
5635      status = 3;
5636    if (lower[i]<-1.0e50&&upper[i]>1.0e50&&status!=1)
5637      status = 0; // set free if should be
5638    else if (lower[i]<-1.0e50&&status==3)
5639      status = 2; // can't be at lower bound
5640    else if (upper[i]>1.0e50&&status==2)
5641      status = 3; // can't be at upper bound
5642    switch (status) {
5643      // free or superbasic
5644    case 0:
5645      if (lower[i]<-1.0e50&&upper[i]>1.0e50) {
5646        mod