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

Last change on this file since 2428 was 2428, checked in by stefan, 8 months ago

change reinterpret_cast of NULL to C-style case, fixes #93

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