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

Last change on this file since 2511 was 2511, checked in by stefan, 16 months ago

compile smallModelCut and modelCut also if CONFLICT_CUTS not defined

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