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

Last change on this file since 1679 was 1679, checked in by forrest, 9 years ago

Add code to help when bad cuts are added by user i.e. ones with tiny
or zero elements. May not fix all problems as Cbc assumes users check
the cuts they generate and so does not do all checks.

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