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

Last change on this file since 2370 was 2370, checked in by forrest, 5 months ago

fix zero row problem and null factorization

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