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

Last change on this file since 2332 was 2332, checked in by forrest, 8 months ago

more tweaks for marginal feasible

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