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

Last change on this file since 2442 was 2442, checked in by forrest, 3 months ago

maximization improvement for readLp

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#ifdef CONFLICT_CUTS
2791// Return a conflict analysis cut from small model
2792OsiRowCut *
2793OsiClpSolverInterface::smallModelCut(const double *originalLower, const double *originalUpper,
2794  int numberRowsAtContinuous, const int *whichGenerator,
2795  int typeCut)
2796{
2797  if (smallModel_ && smallModel_->ray_) {
2798    //printf("Could do small cut\n");
2799    int numberRows = modelPtr_->numberRows();
2800    int numberRows2 = smallModel_->numberRows();
2801    int numberColumns = modelPtr_->numberColumns();
2802    int numberColumns2 = smallModel_->numberColumns();
2803    double *arrayD = reinterpret_cast< double * >(spareArrays_);
2804    double *saveSolution = arrayD + 1;
2805    double *saveLower = saveSolution + (numberRows + numberColumns);
2806    double *saveUpper = saveLower + (numberRows + numberColumns);
2807    double *saveObjective = saveUpper + (numberRows + numberColumns);
2808    double *saveLowerOriginal = saveObjective + (numberRows + numberColumns);
2809    double *saveUpperOriginal = saveLowerOriginal + numberColumns;
2810    //double * lowerOriginal = modelPtr_->columnLower();
2811    //double * upperOriginal = modelPtr_->columnUpper();
2812    int *savePivot = reinterpret_cast< int * >(saveUpperOriginal + numberColumns);
2813    int *whichRow = savePivot + numberRows;
2814    int *whichColumn = whichRow + 3 * numberRows;
2815    int nCopy = 3 * numberRows + 2 * numberColumns;
2816    int nBound = whichRow[nCopy];
2817    if (smallModel_->sequenceOut_ >= 0 && smallModel_->sequenceOut_ < numberColumns2)
2818      modelPtr_->sequenceOut_ = whichColumn[smallModel_->sequenceOut_];
2819    else
2820      modelPtr_->sequenceOut_ = modelPtr_->numberColumns_ + whichRow[smallModel_->sequenceOut_];
2821    unsigned char *saveStatus = CoinCopyOfArray(modelPtr_->status_,
2822      numberRows + numberColumns);
2823    // get ray to full problem
2824    for (int i = 0; i < numberColumns2; i++) {
2825      int iColumn = whichColumn[i];
2826      modelPtr_->setStatus(iColumn, smallModel_->getStatus(i));
2827    }
2828    double *ray = new double[numberRows + numberColumns2 + numberColumns];
2829    char *mark = new char[numberRows];
2830    memset(ray, 0, (numberRows + numberColumns2 + numberColumns) * sizeof(double));
2831    double *smallFarkas = ray + numberRows;
2832    double *farkas = smallFarkas + numberColumns2;
2833    double *saveScale = smallModel_->rowScale_;
2834    smallModel_->rowScale_ = NULL;
2835    smallModel_->transposeTimes(1.0, smallModel_->ray_, smallFarkas);
2836    smallModel_->rowScale_ = saveScale;
2837    for (int i = 0; i < numberColumns2; i++)
2838      farkas[whichColumn[i]] = smallFarkas[i];
2839    memset(mark, 0, numberRows);
2840    for (int i = 0; i < numberRows2; i++) {
2841      int iRow = whichRow[i];
2842      modelPtr_->setRowStatus(iRow, smallModel_->getRowStatus(i));
2843      ray[iRow] = smallModel_->ray_[i];
2844      mark[iRow] = 1;
2845    }
2846    // Column copy of matrix
2847    const double *element = getMatrixByCol()->getElements();
2848    const int *row = getMatrixByCol()->getIndices();
2849    const CoinBigIndex *columnStart = getMatrixByCol()->getVectorStarts();
2850    const int *columnLength = getMatrixByCol()->getVectorLengths();
2851    // pick up small pivot row
2852    int pivotRow = smallModel_->spareIntArray_[3];
2853    // translate
2854    if (pivotRow >= 0)
2855      pivotRow = whichRow[pivotRow];
2856    modelPtr_->spareIntArray_[3] = pivotRow;
2857    for (int jRow = nBound; jRow < 2 * numberRows; jRow++) {
2858      int iRow = whichRow[jRow];
2859      int iColumn = whichRow[jRow + numberRows];
2860      if (modelPtr_->getColumnStatus(iColumn) == ClpSimplex::basic) {
2861        double value = 0.0;
2862        double sum = 0.0;
2863        for (CoinBigIndex j = columnStart[iColumn];
2864             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2865          if (iRow == row[j]) {
2866            value = element[j];
2867          } else if (mark[row[j]]) {
2868            sum += ray[row[j]] * element[j];
2869          }
2870        }
2871        double target = farkas[iColumn];
2872        if (iRow != pivotRow) {
2873          ray[iRow] = (target - sum) / value;
2874        } else {
2875          printf("what now - direction %d wanted %g sum %g value %g\n",
2876            smallModel_->directionOut_, ray[iRow],
2877            sum, value);
2878        }
2879        mark[iRow] = 1;
2880      }
2881    }
2882    delete[] mark;
2883    for (int i = 0; i < modelPtr_->numberColumns_; i++) {
2884      if (modelPtr_->getStatus(i) != ClpSimplex::basic && modelPtr_->columnLower_[i] == modelPtr_->columnUpper_[i])
2885        modelPtr_->setStatus(i, ClpSimplex::isFixed);
2886    }
2887#if 0
2888    {
2889      int nBad=0;
2890      double * fBig = new double [2*numberColumns+2*numberRows];
2891      double * fSmall = fBig+numberColumns;
2892      double * effectiveRhs = fSmall+numberColumns;
2893      double * effectiveRhs2 = effectiveRhs+numberRows;
2894      memset(fBig,0,2*numberColumns*sizeof(double));
2895      {
2896        memcpy(fBig,modelPtr_->columnLower_,numberColumns*sizeof(double));
2897        const double * columnLower = modelPtr_->columnLower_;
2898        const double * columnUpper = modelPtr_->columnUpper_;
2899        for (int i=0;i<numberColumns2;i++) {
2900          int iBig=whichColumn[i];
2901          if (smallModel_->getStatus(i)==ClpSimplex::atLowerBound||
2902              smallModel_->getStatus(i)==ClpSimplex::isFixed||
2903              columnLower[iBig]==columnUpper[iBig]) {
2904            fSmall[i]=columnLower[iBig];
2905          } else if (smallModel_->getStatus(i)==ClpSimplex::atUpperBound) {
2906            fSmall[i]=columnUpper[iBig];
2907          } else if (i==smallModel_->sequenceOut_) {
2908            if (smallModel_->directionOut_<0) {
2909              // above upper bound
2910              fSmall[i]=columnUpper[iBig];
2911            } else {
2912              // below lower bound
2913              fSmall[i]=columnLower[iBig];
2914            }
2915          } else if (smallModel_->columnLower_[i]==smallModel_->columnUpper_[i]) {
2916            fSmall[i]=smallModel_->columnLower_[i];
2917          }
2918          fBig[whichColumn[i]]=fSmall[i];
2919        }
2920        {
2921          // only works for one test case
2922          memcpy(effectiveRhs2,modelPtr_->rowUpper_,numberRows*sizeof(double));
2923          double * saveScale = modelPtr_->rowScale_;
2924          ClpPackedMatrix * saveMatrix = modelPtr_->scaledMatrix_;
2925          modelPtr_->rowScale_=NULL;
2926          modelPtr_->scaledMatrix_=NULL;
2927          modelPtr_->times(-1.0,fBig,effectiveRhs2);
2928          modelPtr_->scaledMatrix_=saveMatrix;
2929          modelPtr_->rowScale_=saveScale;
2930          memset(fBig,0,numberColumns*sizeof(double));
2931        }
2932        const double * rowLower = smallModel_->rowLower_;
2933        const double * rowUpper = smallModel_->rowUpper_;
2934        int pivotRow = smallModel_->spareIntArray_[3];
2935        for (int i=0;i<numberRows2;i++) {
2936          if (smallModel_->getRowStatus(i)==ClpSimplex::atLowerBound||
2937              rowUpper[i]==rowLower[i]||
2938              smallModel_->getRowStatus(i)==ClpSimplex::isFixed) {
2939            effectiveRhs[i]=rowLower[i];
2940            //if (smallModel_->getRowStatus(i)==ClpSimplex::atLowerBound)
2941            //assert (ray[i]>-1.0e-3);
2942          } else if (smallModel_->getRowStatus(i)==ClpSimplex::atUpperBound) {
2943            effectiveRhs[i]=rowUpper[i];
2944            //assert (ray[i]<1.0e-3);
2945          } else {
2946            if (smallModel_->getRowStatus(i)!=ClpSimplex::basic) {
2947              assert (rowUpper[i]>1.0e30||rowLower[i]<-1.0e30); // eventually skip
2948              if (rowUpper[i]<1.0e30)
2949                effectiveRhs[i]=rowUpper[i];
2950              else
2951                effectiveRhs[i]=rowLower[i];
2952            }
2953          }
2954          if (smallModel_->getRowStatus(i)==ClpSimplex::basic) {
2955            effectiveRhs[i]=0.0;
2956            // we should leave pivot row in and use direction for bound
2957            if (fabs(smallModel_->ray_[i])>1.0e-8) {
2958              printf("Basic small slack value %g on %d - pivotRow %d\n",smallModel_->ray_[i],i,pivotRow);
2959              if (i==pivotRow) {
2960                if (smallModel_->directionOut_<0)
2961                  effectiveRhs[i]=rowUpper[i];
2962                else
2963                  effectiveRhs[i]=rowLower[i];
2964              } else {
2965                //smallModel_->ray_[i]=0.0;
2966              }
2967            }
2968          }
2969        }
2970        //ClpPackedMatrix * saveMatrix = smallModel_->scaledMatrix_;
2971        double * saveScale = smallModel_->rowScale_;
2972        smallModel_->rowScale_=NULL;
2973        smallModel_->scaledMatrix_=NULL;
2974        smallModel_->times(-1.0,fSmall,effectiveRhs);
2975        double bSum=0.0;
2976        for (int i=0;i<numberRows2;i++) {
2977          bSum += effectiveRhs[i]*smallModel_->ray_[i];
2978        }
2979        printf("Small bsum %g\n",bSum);
2980        smallModel_->rowScale_=saveScale;
2981      }
2982      double * saveScale = smallModel_->rowScale_;
2983      smallModel_->rowScale_=NULL;
2984      memset(fSmall,0,numberColumns*sizeof(double));
2985      smallModel_->transposeTimes(1.0,smallModel_->ray_,fSmall);
2986      smallModel_->rowScale_=saveScale;
2987      ClpPackedMatrix * saveMatrix = modelPtr_->scaledMatrix_;
2988      saveScale = modelPtr_->rowScale_;
2989      modelPtr_->rowScale_=NULL;
2990      modelPtr_->scaledMatrix_=NULL;
2991      modelPtr_->transposeTimes(1.0,ray,fBig);
2992      modelPtr_->scaledMatrix_=saveMatrix;
2993      modelPtr_->rowScale_=saveScale;
2994      for (int i=0;i<numberColumns2;i++) {
2995        int iBig=whichColumn[i];
2996        if (fabs(fBig[iBig]-fSmall[i])>1.0e-4) {
2997          printf("small %d %g big %d %g\n",
2998                 i,fSmall[i],iBig,fBig[iBig]);
2999          nBad++;
3000        }
3001      }
3002      delete [] fBig;
3003      if (nBad)
3004        printf("Bad %d odd %d\n",nBad,2*numberRows-nBound);
3005    }
3006#endif
3007    modelPtr_->ray_ = ray;
3008    lastAlgorithm_ = 2;
3009    modelPtr_->directionOut_ = smallModel_->directionOut_;
3010    OsiRowCut *cut = modelCut(originalLower, originalUpper,
3011      numberRowsAtContinuous, whichGenerator, typeCut);
3012    smallModel_->deleteRay();
3013    memcpy(modelPtr_->status_, saveStatus, numberRows + numberColumns);
3014    delete[] saveStatus;
3015    if (cut) {
3016      //printf("got small cut\n");
3017      //delete cut;
3018      //cut=NULL;
3019    }
3020    return cut;
3021  } else {
3022    return NULL;
3023  }
3024}
3025static int debugMode = 0;
3026// Return a conflict analysis cut from model
3027//      If type is 0 then genuine cut, if 1 then only partially processed
3028OsiRowCut *
3029OsiClpSolverInterface::modelCut(const double *originalLower, const double *originalUpper,
3030  int numberRowsAtContinuous, const int *whichGenerator,
3031  int typeCut)
3032{
3033  OsiRowCut *cut = NULL;
3034  //return NULL;
3035  if (modelPtr_->ray_) {
3036    if (lastAlgorithm_ == 2) {
3037      //printf("Could do normal cut\n");
3038      // could use existing arrays
3039      int numberRows = modelPtr_->numberRows_;
3040      int numberColumns = modelPtr_->numberColumns_;
3041      double *farkas = new double[2 * numberColumns + numberRows];
3042      double *bound = farkas + numberColumns;
3043      double *effectiveRhs = bound + numberColumns;
3044      /*const*/ double *ray = modelPtr_->ray_;
3045      // have to get rid of local cut rows
3046      whichGenerator -= numberRowsAtContinuous;
3047      int badRows = 0;
3048      for (int iRow = numberRowsAtContinuous; iRow < numberRows; iRow++) {
3049        int iType = whichGenerator[iRow];
3050        if ((iType >= 0 && iType < 20000)) {
3051          if (fabs(ray[iRow]) > 1.0e-10) {
3052            badRows++;
3053          }
3054          ray[iRow] = 0.0;
3055        }
3056      }
3057      ClpSimplex debugModel;
3058      if ((debugMode & 4) != 0)
3059        debugModel = *modelPtr_;
3060      if (badRows && (debugMode & 1) != 0)
3061        printf("%d rows from local cuts\n", badRows);
3062      // get farkas row
3063      ClpPackedMatrix *saveMatrix = modelPtr_->scaledMatrix_;
3064      double *saveScale = modelPtr_->rowScale_;
3065      modelPtr_->rowScale_ = NULL;
3066      modelPtr_->scaledMatrix_ = NULL;
3067      memset(farkas, 0, (2 * numberColumns + numberRows) * sizeof(double));
3068      modelPtr_->transposeTimes(-1.0, ray, farkas);
3069      //const char * integerInformation = modelPtr_->integerType_;
3070      //assert (integerInformation);
3071
3072      int sequenceOut = modelPtr_->sequenceOut_;
3073      // Put nonzero bounds in bound
3074      const double *columnLower = modelPtr_->columnLower_;
3075      const double *columnUpper = modelPtr_->columnUpper_;
3076      const double *columnValue = modelPtr_->columnActivity_;
3077      int numberBad = 0;
3078      int nNonzeroBasic = 0;
3079      for (int i = 0; i < numberColumns; i++) {
3080        double value = farkas[i];
3081        double boundValue = 0.0;
3082        if (modelPtr_->getStatus(i) == ClpSimplex::basic) {
3083          // treat as zero if small
3084          if (fabs(value) < 1.0e-8) {
3085            value = 0.0;
3086            farkas[i] = 0.0;
3087          }
3088          if (value) {
3089            //printf("basic %d direction %d farkas %g\n",
3090            //i,modelPtr_->directionOut_,value);
3091            nNonzeroBasic++;
3092            if (value < 0.0)
3093              boundValue = columnLower[i];
3094            else
3095              boundValue = columnUpper[i];
3096          }
3097        } else if (fabs(value) > 1.0e-10) {
3098          if (value < 0.0) {
3099            if (columnValue[i] > columnLower[i] + 1.0e-5 && value < -1.0e-7) {
3100              //printf("bad %d alpha %g not at lower\n",i,value);
3101              numberBad++;
3102            }
3103            boundValue = columnLower[i];
3104          } else {
3105            if (columnValue[i] < columnUpper[i] - 1.0e-5 && value > 1.0e-7) {
3106              //printf("bad %d alpha %g not at upper\n",i,value);
3107              numberBad++;
3108            }
3109            boundValue = columnUpper[i];
3110          }
3111        }
3112        bound[i] = boundValue;
3113        if (fabs(boundValue) > 1.0e10)
3114          numberBad++;
3115      }
3116      const double *rowLower = modelPtr_->rowLower_;
3117      const double *rowUpper = modelPtr_->rowUpper_;
3118      //int pivotRow = modelPtr_->spareIntArray_[3];
3119      //bool badPivot=pivotRow<0;
3120      for (int i = 0; i < numberRows; i++) {
3121        double value = ray[i];
3122        double rhsValue = 0.0;
3123        if (modelPtr_->getRowStatus(i) == ClpSimplex::basic) {
3124          // treat as zero if small
3125          if (fabs(value) < 1.0e-8) {
3126            value = 0.0;
3127            ray[i] = 0.0;
3128          }
3129          if (value) {
3130            //printf("row basic %d direction %d ray %g\n",
3131            //     i,modelPtr_->directionOut_,value);
3132            nNonzeroBasic++;
3133            if (value < 0.0)
3134              rhsValue = rowLower[i];
3135            else
3136              rhsValue = rowUpper[i];
3137          }
3138        } else if (fabs(value) > 1.0e-10) {
3139          if (value < 0.0)
3140            rhsValue = rowLower[i];
3141          else
3142            rhsValue = rowUpper[i];
3143        }
3144        effectiveRhs[i] = rhsValue;
3145      }
3146      {
3147        double bSum = 0.0;
3148        for (int i = 0; i < numberRows; i++) {
3149          bSum += effectiveRhs[i] * ray[i];
3150        }
3151        //printf("before bounds - bSum %g\n",bSum);
3152      }
3153      modelPtr_->times(-1.0, bound, effectiveRhs);
3154      double bSum = 0.0;
3155      for (int i = 0; i < numberRows; i++) {
3156        bSum += effectiveRhs[i] * ray[i];
3157      }
3158#if 0
3159      double bSum2=0.0;
3160#if 1
3161      {
3162        int iSequence = modelPtr_->sequenceOut_;
3163        double lower,value,upper;
3164        if (iSequence<numberColumns) {
3165          lower=columnLower[iSequence];
3166          value=columnValue[iSequence];
3167          upper=columnUpper[iSequence];
3168        } else {
3169          iSequence -= numberColumns;
3170          lower=rowLower[iSequence];
3171          value=modelPtr_->rowActivity_[iSequence];
3172          upper=rowUpper[iSequence];
3173        }
3174        if (value<lower)
3175          bSum2=value-lower;
3176        else if (value>upper)
3177          bSum2=-(value-upper);
3178        else
3179          printf("OUCHXX %g %g %g\n",
3180                 lower,value,upper);
3181      }
3182#else
3183      if (modelPtr_->valueOut_<modelPtr_->lowerOut_)
3184        bSum2=modelPtr_->valueOut_-modelPtr_->lowerOut_;
3185      else if (modelPtr_->valueOut_>modelPtr_->upperOut_)
3186        bSum2=-(modelPtr_->valueOut_-modelPtr_->upperOut_);
3187      else
3188        printf("OUCHXX %g %g %g\n",
3189               modelPtr_->lowerOut_,modelPtr_->valueOut_,modelPtr_->upperOut_);
3190#endif
3191#if 0
3192      if (fabs(bSum-bSum2)>1.0e-5)
3193        printf("XX ");
3194      printf("after bounds - bSum %g - bSum2 %g\n",bSum,bSum2);
3195#endif
3196#endif
3197      modelPtr_->scaledMatrix_ = saveMatrix;
3198      modelPtr_->rowScale_ = saveScale;
3199      if (numberBad || bSum > -1.0e-4 /*||nNonzeroBasic>1*/ /*||bSum2>-1.0e-4*/) {
3200#if PRINT_CONFLICT > 1 //ndef NDEBUG
3201        printf("bad BOUND bSum %g (bSum2 %g) - %d bad - %d basic\n",
3202          bSum, bSum2, numberBad, nNonzeroBasic);
3203#endif
3204      } else {
3205        if (numberColumns < 0)
3206          debugMode = -numberColumns;
3207        if ((debugMode & 4) != 0) {
3208          int *tempC = new int[numberColumns];
3209          double *temp = new double[numberColumns];
3210          int n = 0;
3211          for (int i = 0; i < numberColumns; i++) {
3212            if (fabs(farkas[i]) > 1.0e-12) {
3213              temp[n] = farkas[i];
3214              tempC[n++] = i;
3215            }
3216          }
3217          debugModel.addRow(n, tempC, temp, bSum, bSum);
3218          delete[] tempC;
3219          delete[] temp;
3220        }
3221        if ((debugMode & 2) != 0) {
3222          ClpSimplex dual = *modelPtr_;
3223          dual.setLogLevel(63);
3224          dual.scaling(0);
3225          dual.dual();
3226          assert(dual.problemStatus_ == 1);
3227          if (dual.ray_) {
3228            double *farkas2 = dual.reducedCost_;
3229            // Put nonzero bounds in farkas
3230            const double *columnLower = dual.columnLower_;
3231            const double *columnUpper = dual.columnUpper_;
3232            for (int i = 0; i < numberColumns; i++) {
3233              farkas2[i] = 0.0;
3234              if (dual.getStatus(i) == ClpSimplex::atLowerBound || columnLower[i] == columnUpper[i]) {
3235                farkas2[i] = columnLower[i];
3236              } else if (dual.getStatus(i) == ClpSimplex::atUpperBound) {
3237                farkas2[i] = columnUpper[i];
3238              } else if (i == dual.sequenceOut_) {
3239                if (dual.directionOut_ < 0) {
3240                  // above upper bound
3241                  farkas2[i] = columnUpper[i];
3242                } else {
3243                  // below lower bound
3244                  farkas2[i] = columnLower[i];
3245                }
3246              } else if (columnLower[i] == columnUpper[i]) {
3247                farkas2[i] = columnLower[i];
3248              }
3249            }
3250            double *effectiveRhs2 = dual.rowActivity_;
3251            const double *rowLower = dual.rowLower_;
3252            const double *rowUpper = dual.rowUpper_;
3253            int pivotRow = dual.spareIntArray_[3];
3254            for (int i = 0; i < numberRows; i++) {
3255              if (dual.getRowStatus(i) == ClpSimplex::atLowerBound || rowUpper[i] == rowLower[i] || dual.getRowStatus(i) == ClpSimplex::isFixed) {
3256                effectiveRhs2[i] = rowLower[i];
3257              } else if (dual.getRowStatus(i) == ClpSimplex::atUpperBound) {
3258                effectiveRhs2[i] = rowUpper[i];
3259              } else {
3260                if (dual.getRowStatus(i) != ClpSimplex::basic) {
3261                  assert(rowUpper[i] > 1.0e30 || rowLower[i] < -1.0e30); // eventually skip
3262                  if (rowUpper[i] < 1.0e30)
3263                    effectiveRhs2[i] = rowUpper[i];
3264                  else
3265                    effectiveRhs2[i] = rowLower[i];
3266                }
3267              }
3268              if (dual.getRowStatus(i) == ClpSimplex::basic) {
3269                effectiveRhs2[i] = 0.0;
3270                // we should leave pivot row in and use direction for bound
3271                if (fabs(dual.ray_[i]) > 1.0e-8) {
3272                  printf("Basic slack value %g on %d - pivotRow %d\n", ray[i], i, pivotRow);
3273                  if (i == pivotRow) {
3274                    if (dual.directionOut_ < 0)
3275                      effectiveRhs[i] = rowUpper[i];
3276                    else
3277                      effectiveRhs[i] = rowLower[i];
3278                  } else {
3279                    dual.ray_[i] = 0.0;
3280                  }
3281                }
3282              }
3283            }
3284            dual.times(-1.0, farkas2, effectiveRhs2);
3285            double bSum2 = 0.0;
3286            for (int i = 0; i < numberRows; i++) {
3287              bSum2 += effectiveRhs2[i] * dual.ray_[i];
3288            }
3289            printf("Alternate bSum %g\n", bSum2);
3290            memset(farkas2, 0, numberColumns * sizeof(double));
3291            dual.transposeTimes(-1.0, dual.ray_, farkas2);
3292            int nBad = 0;
3293            double largest = -1.0;
3294            double smallest = 1.0e30;
3295            for (int i = 0; i < numberColumns; i++) {
3296              //if (fabs(farkas[i])>1.0e-12)
3297              //printf("%d ptr farkas %g dual farkas %g\n",i,farkas[i],farkas2[i]);
3298              if (fabs(farkas[i] - farkas2[i]) > 1.0e-7) {
3299                nBad++;
3300                largest = CoinMax(largest, fabs(farkas[i] - farkas2[i]));
3301                smallest = CoinMin(smallest, fabs(farkas[i] - farkas2[i]));
3302                //printf("%d ptr farkas %g dual farkas %g\n",i,farkas[i],farkas2[i]);
3303              }
3304            }
3305            if (nBad)
3306              printf("%d farkas difference %g to %g\n", nBad, smallest, largest);
3307            dual.primal();
3308            assert(dual.problemStatus_ == 1);
3309            assert(!nBad);
3310          }
3311        }
3312        const char *integerInformation = modelPtr_->integerType_;
3313        assert(integerInformation);
3314        int *conflict = new int[numberColumns];
3315        double *sort = new double[numberColumns];
3316        double relax = 0.0;
3317        int nConflict = 0;
3318        int nOriginal = 0;
3319        int nFixed = 0;
3320        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3321          double thisRelax = 0.0;
3322          if (integerInformation[iColumn]) {
3323            if ((debugMode & 1) != 0)
3324              printf("%d status %d %g <= %g <=%g (orig %g, %g) farkas %g\n",
3325                iColumn, modelPtr_->getStatus(iColumn), columnLower[iColumn],
3326                modelPtr_->columnActivity_[iColumn], columnUpper[iColumn],
3327                originalLower[iColumn], originalUpper[iColumn],
3328                farkas[iColumn]);
3329            double gap = originalUpper[iColumn] - originalLower[iColumn];
3330            if (!gap)
3331              continue;
3332            if (gap == columnUpper[iColumn] - columnLower[iColumn])
3333              nOriginal++;
3334            if (columnUpper[iColumn] == columnLower[iColumn])
3335              nFixed++;
3336            if (fabs(farkas[iColumn]) < 1.0e-15) {
3337              farkas[iColumn] = 0.0;
3338              continue;
3339            }
3340            // temp
3341            if (gap >= 20000.0 && false) {
3342              // can't use
3343              if (farkas[iColumn] < 0.0) {
3344                assert(originalLower[iColumn] - columnLower[iColumn] <= 0.0);
3345                // farkas is negative - relax lower bound all way
3346                relax += farkas[iColumn] * (originalLower[iColumn] - columnLower[iColumn]);
3347              } else {
3348                assert(originalUpper[iColumn] - columnUpper[iColumn] >= 0.0);
3349                // farkas is positive - relax upper bound all way
3350                relax += farkas[iColumn] * (originalUpper[iColumn] - columnUpper[iColumn]);
3351              }
3352              continue;
3353            }
3354            if (originalLower[iColumn] == columnLower[iColumn]) {
3355              // may need to be careful if odd basic (due to local cut)
3356              if (farkas[iColumn] > 0.0 /*&&(modelPtr_->getStatus(iColumn)==ClpSimplex::atUpperBound
3357                                        ||modelPtr_->getStatus(iColumn)==ClpSimplex::isFixed
3358                                        ||iColumn==sequenceOut)*/
3359              ) {
3360                // farkas is positive - add to list
3361                gap = originalUpper[iColumn] - columnUpper[iColumn];
3362                if (gap) {
3363                  sort[nConflict] = -farkas[iColumn] * gap;
3364                  conflict[nConflict++] = iColumn;
3365                }
3366                //assert (gap>columnUpper[iColumn]-columnLower[iColumn]);
3367              }
3368            } else if (originalUpper[iColumn] == columnUpper[iColumn]) {
3369              // may need to be careful if odd basic (due to local cut)
3370              if (farkas[iColumn] < 0.0 /*&&(modelPtr_->getStatus(iColumn)==ClpSimplex::atLowerBound
3371                                        ||modelPtr_->getStatus(iColumn)==ClpSimplex::isFixed
3372                                        ||iColumn==sequenceOut)*/
3373              ) {
3374                // farkas is negative - add to list
3375                gap = columnLower[iColumn] - originalLower[iColumn];
3376                if (gap) {
3377                  sort[nConflict] = farkas[iColumn] * gap;
3378                  conflict[nConflict++] = iColumn;
3379                }
3380                //assert (gap>columnUpper[iColumn]-columnLower[iColumn]);
3381              }
3382            } else {
3383              // can't use
3384              if (farkas[iColumn] < 0.0) {
3385                assert(originalLower[iColumn] - columnLower[iColumn] <= 0.0);
3386                // farkas is negative - relax lower bound all way
3387                thisRelax = farkas[iColumn] * (originalLower[iColumn] - columnLower[iColumn]);
3388              } else {
3389                assert(originalUpper[iColumn] - columnUpper[iColumn] >= 0.0);
3390                // farkas is positive - relax upper bound all way
3391                thisRelax = farkas[iColumn] * (originalUpper[iColumn] - columnUpper[iColumn]);
3392              }
3393            }
3394          } else {
3395            // not integer - but may have been got at
3396            double gap = originalUpper[iColumn] - originalLower[iColumn];
3397            if (gap > columnUpper[iColumn] - columnLower[iColumn]) {
3398              // can't use
3399              if (farkas[iColumn] < 0.0) {
3400                assert(originalLower[iColumn] - columnLower[iColumn] <= 0.0);
3401                // farkas is negative - relax lower bound all way
3402                thisRelax = farkas[iColumn] * (originalLower[iColumn] - columnLower[iColumn]);
3403              } else {
3404                assert(originalUpper[iColumn] - columnUpper[iColumn] >= 0.0);
3405                // farkas is positive - relax upper bound all way
3406                thisRelax = farkas[iColumn] * (originalUpper[iColumn] - columnUpper[iColumn]);
3407              }
3408            }
3409          }
3410          assert(thisRelax >= 0.0);
3411          relax += thisRelax;
3412        }
3413        if (relax + bSum > -1.0e-4 || !nConflict) {
3414          if (relax + bSum > -1.0e-4) {
3415#if PRINT_CONFLICT > 1 //ndef NDEBUG
3416            printf("General integers relax bSum to %g\n", relax + bSum);
3417#endif
3418          } else {
3419#if PRINT_CONFLICT
3420            printf("All variables relaxed and still infeasible - what does this mean?\n");
3421#endif
3422            int nR = 0;
3423            for (int i = 0; i < numberRows; i++) {
3424              if (fabs(ray[i]) > 1.0e-10)
3425                nR++;
3426              else
3427                ray[i] = 0.0;
3428            }
3429            int nC = 0;
3430            for (int i = 0; i < numberColumns; i++) {
3431              if (fabs(farkas[i]) > 1.0e-10)
3432                nC++;
3433              else
3434                farkas[i] = 0.0;
3435            }
3436            if (nR < 3 && nC < 5) {
3437              printf("BAD %d nonzero rows, %d nonzero columns\n", nR, nC);
3438            }
3439          }
3440        } else if (nConflict < 1000) {
3441#if PRINT_CONFLICT > 1 //ndef NDEBUG
3442          if (nConflict < 5)
3443            printf("BOUNDS violation bSum %g (relaxed %g) - %d at original bounds, %d fixed - %d in conflict\n", bSum,
3444              relax + bSum, nOriginal, nFixed, nConflict);
3445#endif
3446          CoinSort_2(sort, sort + nConflict, conflict);
3447          if ((debugMode & 4) != 0) {
3448            double *temp = new double[numberColumns];
3449            int *tempC = new int[numberColumns];
3450            int n = 0;
3451            for (int j = 0; j < nConflict; j++) {
3452              int i = conflict[j];
3453              if (fabs(farkas[i]) > 1.0e-12) {
3454                temp[n] = farkas[i];
3455                tempC[n++] = i;
3456              }
3457            }
3458            debugModel.addRow(n, tempC, temp, bSum, bSum);
3459            delete[] tempC;
3460            delete[] temp;
3461          }
3462          int nC = nConflict;
3463          for (int i = 0; i < nConflict; i++) {
3464            int iColumn = conflict[i];
3465            if (fabs(sort[i]) != fabs(farkas[iColumn]) && originalUpper[iColumn] == 1.0)
3466              printf("odd %d %g %d %g\n", i, sort[i], iColumn, farkas[iColumn]);
3467          }
3468          bSum += relax;
3469          double saveBsum = bSum;
3470          // we can't use same values
3471          double totalChange = 0;
3472          while (nConflict) {
3473            double last = -sort[nConflict - 1];
3474            int kConflict = nConflict;
3475            while (kConflict) {
3476              double change = -sort[kConflict - 1];
3477              if (change > last + 1.0e-5)
3478                break;
3479              totalChange += change;
3480              kConflict--;
3481            }
3482            if (bSum + totalChange > -1.0e-4)
3483              break;
3484#if 0
3485            //int iColumn=conflict[nConflict-1];
3486            double change=-sort[nConflict-1];
3487            if (bSum+change>-1.0e-4)
3488              break;
3489            nConflict--;
3490            bSum += change;
3491#else
3492            nConflict = kConflict;
3493            bSum += totalChange;
3494#endif
3495          }
3496          if (!nConflict) {
3497            int nR = 0;
3498            for (int i = 0; i < numberRows; i++) {
3499              if (fabs(ray[i]) > 1.0e-10)
3500                nR++;
3501              else
3502                ray[i] = 0.0;
3503            }
3504            int nC = 0;
3505            for (int i = 0; i < numberColumns; i++) {
3506              if (fabs(farkas[i]) > 1.0e-10)
3507                nC++;
3508              else
3509                farkas[i] = 0.0;
3510            }
3511#if 1 //PRINT_CONFLICT>1 //ndef NDEBUG
3512            {
3513              printf("BAD2 - zero nConflict %d nonzero rows, %d nonzero columns\n", nR, nC);
3514            }
3515#endif
3516          }
3517          // no point doing if no reduction (or big?) ?
3518          if (nConflict < nC + 1 && nConflict < 100 && nConflict) {
3519            cut = new OsiRowCut();
3520            cut->setUb(COIN_DBL_MAX);
3521            if (!typeCut) {
3522              double lo = 1.0;
3523              for (int i = 0; i < nConflict; i++) {
3524                int iColumn = conflict[i];
3525                if (originalLower[iColumn] == columnLower[iColumn]) {
3526                  // must be at least one higher
3527                  sort[i] = 1.0;
3528                  lo += originalLower[iColumn];
3529                } else {
3530                  // must be at least one lower
3531                  sort[i] = -1.0;
3532                  lo -= originalUpper[iColumn];
3533                }
3534              }
3535              cut->setLb(lo);
3536              cut->setRow(nConflict, conflict, sort);
3537#if PRINT_CONFLICT
3538              printf("CUT has %d (started at %d) - final bSum %g\n", nConflict, nC, bSum);
3539#endif
3540              if ((debugMode & 4) != 0) {
3541                debugModel.addRow(nConflict, conflict, sort, lo, COIN_DBL_MAX);
3542                debugModel.writeMps("bad.mps");
3543              }
3544            } else {
3545              // just save for use later
3546              // first take off small
3547              int nC2 = nC;
3548              while (nC2) {
3549                //int iColumn=conflict[nConflict-1];
3550                double change = -sort[nC2 - 1];
3551                if (saveBsum + change > -1.0e-4 || change > 1.0e-4)
3552                  break;
3553                nC2--;
3554                saveBsum += change;
3555              }
3556              cut->setLb(saveBsum);
3557              for (int i = 0; i < nC2; i++) {
3558                int iColumn = conflict[i];
3559                sort[i] = farkas[iColumn];
3560              }
3561              cut->setRow(nC2, conflict, sort);
3562              // mark as globally valid
3563              cut->setGloballyValid();
3564#if PRINT_CONFLICT > 1 //ndef NDEBUG
3565              printf("Stem CUT has %d (greedy %d - with small %d) - saved bSum %g final greedy bSum %g\n",
3566                nC2, nConflict, nC, saveBsum, bSum);
3567#endif
3568            }
3569          }
3570        }
3571        delete[] conflict;
3572        delete[] sort;
3573      }
3574      delete[] farkas;
3575    } else {
3576#if PRINT_CONFLICT > 1 //ndef NDEBUG
3577      printf("Bad dual ray\n");
3578#endif
3579    }
3580    modelPtr_->deleteRay();
3581  }
3582  return cut;
3583}
3584#endif
3585
3586//#############################################################################
3587// Problem information methods (original data)
3588//#############################################################################
3589
3590//------------------------------------------------------------------
3591const char *OsiClpSolverInterface::getRowSense() const
3592{
3593  extractSenseRhsRange();
3594  return rowsense_;
3595}
3596//------------------------------------------------------------------
3597const double *OsiClpSolverInterface::getRightHandSide() const
3598{
3599  extractSenseRhsRange();
3600  return rhs_;
3601}
3602//------------------------------------------------------------------
3603const double *OsiClpSolverInterface::getRowRange() const
3604{
3605  extractSenseRhsRange();
3606  return rowrange_;
3607}
3608//------------------------------------------------------------------
3609// Return information on integrality
3610//------------------------------------------------------------------
3611bool OsiClpSolverInterface::isContinuous(int colNumber) const
3612{
3613  if (integerInformation_ == NULL)
3614    return true;
3615#ifndef NDEBUG
3616  int n = modelPtr_->numberColumns();
3617  if (colNumber < 0 || colNumber >= n) {
3618    indexError(colNumber, "isContinuous");
3619  }
3620#endif
3621  if (integerInformation_[colNumber] == 0)
3622    return true;
3623  return false;
3624}
3625bool OsiClpSolverInterface::isBinary(int colNumber) const
3626{
3627#ifndef NDEBUG
3628  int n = modelPtr_->numberColumns();
3629  if (colNumber < 0 || colNumber >= n) {
3630    indexError(colNumber, "isBinary");
3631  }
3632#endif
3633  if (integerInformation_ == NULL || integerInformation_[colNumber] == 0) {
3634    return false;
3635  } else {
3636    const double *cu = getColUpper();
3637    const double *cl = getColLower();
3638    if ((cu[colNumber] == 1 || cu[colNumber] == 0) && (cl[colNumber] == 0 || cl[colNumber] == 1))
3639      return true;
3640    else
3641      return false;
3642  }
3643}
3644//-----------------------------------------------------------------------------
3645bool OsiClpSolverInterface::isInteger(int colNumber) const
3646{
3647#ifndef NDEBUG
3648  int n = modelPtr_->numberColumns();
3649  if (colNumber < 0 || colNumber >= n) {
3650    indexError(colNumber, "isInteger");
3651  }
3652#endif
3653  if (integerInformation_ == NULL || integerInformation_[colNumber] == 0)
3654    return false;
3655  else
3656    return true;
3657}
3658//-----------------------------------------------------------------------------
3659bool OsiClpSolverInterface::isIntegerNonBinary(int colNumber) const
3660{
3661#ifndef NDEBUG
3662  int n = modelPtr_->numberColumns();
3663  if (colNumber < 0 || colNumber >= n) {
3664    indexError(colNumber, "isIntegerNonBinary");
3665  }
3666#endif
3667  if (integerInformation_ == NULL || integerInformation_[colNumber] == 0) {
3668    return false;
3669  } else {
3670    return !isBinary(colNumber);
3671  }
3672}
3673//-----------------------------------------------------------------------------
3674bool OsiClpSolverInterface::isFreeBinary(int colNumber) const
3675{
3676#ifndef NDEBUG
3677  int n = modelPtr_->numberColumns();
3678  if (colNumber < 0 || colNumber >= n) {
3679    indexError(colNumber, "isFreeBinary");
3680  }
3681#endif
3682  if (integerInformation_ == NULL || integerInformation_[colNumber] == 0) {
3683    return false;
3684  } else {
3685    const double *cu = getColUpper();
3686    const double *cl = getColLower();
3687    if ((cu[colNumber] == 1) && (cl[colNumber] == 0))
3688      return true;
3689    else
3690      return false;
3691  }
3692}
3693/*  Return array of column length
3694    0 - continuous
3695    1 - binary (may get fixed later)
3696    2 - general integer (may get fixed later)
3697*/
3698const char *
3699OsiClpSolverInterface::getColType(bool refresh) const
3700{
3701  if (!columnType_ || refresh) {
3702    const int numCols = getNumCols();
3703    if (!columnType_)
3704      columnType_ = new char[numCols];
3705    if (integerInformation_ == NULL) {
3706      memset(columnType_, 0, numCols);
3707    } else {
3708      const double *cu = getColUpper();
3709      const double *cl = getColLower();
3710      for (int i = 0; i < numCols; ++i) {
3711        if (integerInformation_[i]) {
3712          if ((cu[i] == 1 || cu[i] == 0) && (cl[i] == 0 || cl[i] == 1))
3713            columnType_[i] = 1;
3714          else
3715            columnType_[i] = 2;
3716        } else {
3717          columnType_[i] = 0;
3718        }
3719      }
3720    }
3721  }
3722  return columnType_;
3723}
3724
3725/* Return true if column is integer but does not have to
3726   be declared as such.
3727   Note: This function returns true if the the column
3728   is binary or a general integer.
3729*/
3730bool OsiClpSolverInterface::isOptionalInteger(int colNumber) const
3731{
3732#ifndef NDEBUG
3733  int n = modelPtr_->numberColumns();
3734  if (colNumber < 0 || colNumber >= n) {
3735    indexError(colNumber, "isInteger");
3736  }
3737#endif
3738  if (integerInformation_ == NULL || integerInformation_[colNumber] != 2)
3739    return false;
3740  else
3741    return true;
3742}
3743/* Set the index-th variable to be an optional integer variable */
3744void OsiClpSolverInterface::setOptionalInteger(int index)
3745{
3746  if (!integerInformation_) {
3747    integerInformation_ = new char[modelPtr_->numberColumns()];
3748    CoinFillN(integerInformation_, modelPtr_->numberColumns(), static_cast< char >(0));
3749  }
3750#ifndef NDEBUG
3751  int n = modelPtr_->numberColumns();
3752  if (index < 0 || index >= n) {
3753    indexError(index, "setInteger");
3754  }
3755#endif
3756  integerInformation_[index] = 2;
3757  modelPtr_->setInteger(index);
3758}
3759
3760//------------------------------------------------------------------
3761// Row and column copies of the matrix ...
3762//------------------------------------------------------------------
3763const CoinPackedMatrix *OsiClpSolverInterface::getMatrixByRow() const
3764{
3765  if (matrixByRow_ == NULL || matrixByRow_->getNumElements() != modelPtr_->clpMatrix()->getNumElements()) {
3766    delete matrixByRow_;
3767    matrixByRow_ = new CoinPackedMatrix();
3768    matrixByRow_->setExtraGap(0.0);
3769    matrixByRow_->setExtraMajor(0.0);
3770    matrixByRow_->reverseOrderedCopyOf(*modelPtr_->matrix());
3771    //matrixByRow_->removeGaps();
3772#if 0
3773    CoinPackedMatrix back;
3774    std::cout<<"start check"<<std::endl;
3775    back.reverseOrderedCopyOf(*matrixByRow_);
3776    modelPtr_->matrix()->isEquivalent2(back);
3777    std::cout<<"stop check"<<std::endl;
3778#endif
3779  }
3780  assert(matrixByRow_->getNumElements() == modelPtr_->clpMatrix()->getNumElements());
3781  return matrixByRow_;
3782}
3783
3784const CoinPackedMatrix *OsiClpSolverInterface::getMatrixByCol() const
3785{
3786  return modelPtr_->matrix();
3787}
3788
3789// Get pointer to mutable column-wise copy of matrix (returns NULL if not meaningful)
3790CoinPackedMatrix *
3791OsiClpSolverInterface::getMutableMatrixByCol() const
3792{
3793  ClpPackedMatrix *matrix = dynamic_cast< ClpPackedMatrix * >(modelPtr_->matrix_);
3794  if (matrix)
3795    return matrix->getPackedMatrix();
3796  else
3797    return NULL;
3798}
3799
3800//------------------------------------------------------------------
3801std::vector< double * > OsiClpSolverInterface::getDualRays(int /*maxNumRays*/,
3802  bool fullRay) const
3803{
3804  return std::vector< double * >(1, modelPtr_->infeasibilityRay(fullRay));
3805}
3806//------------------------------------------------------------------
3807std::vector< double * > OsiClpSolverInterface::getPrimalRays(int /*maxNumRays*/) const
3808{
3809  return std::vector< double * >(1, modelPtr_->unboundedRay());
3810}
3811//#############################################################################
3812void OsiClpSolverInterface::setContinuous(int index)
3813{
3814
3815  if (integerInformation_) {
3816#ifndef NDEBUG
3817    int n = modelPtr_->numberColumns();
3818    if (index < 0 || index >= n) {
3819      indexError(index, "setContinuous");
3820    }
3821#endif
3822    integerInformation_[index] = 0;
3823  }
3824  modelPtr_->setContinuous(index);
3825}
3826//-----------------------------------------------------------------------------
3827void OsiClpSolverInterface::setInteger(int index)
3828{
3829  if (!integerInformation_) {
3830    integerInformation_ = new char[modelPtr_->numberColumns()];
3831    CoinFillN(integerInformation_, modelPtr_->numberColumns(), static_cast< char >(0));
3832  }
3833#ifndef NDEBUG
3834  int n = modelPtr_->numberColumns();
3835  if (index < 0 || index >= n) {
3836    indexError(index, "setInteger");
3837  }
3838#endif
3839  integerInformation_[index] = 1;
3840  modelPtr_->setInteger(index);
3841}
3842//-----------------------------------------------------------------------------
3843void OsiClpSolverInterface::setContinuous(const int *indices, int len)
3844{
3845  if (integerInformation_) {
3846#ifndef NDEBUG
3847    int n = modelPtr_->numberColumns();
3848#endif
3849    int i;
3850    for (i = 0; i < len; i++) {
3851      int colNumber = indices[i];
3852#ifndef NDEBUG
3853      if (colNumber < 0 || colNumber >= n) {
3854        indexError(colNumber, "setContinuous");
3855      }
3856#endif
3857      integerInformation_[colNumber] = 0;
3858      modelPtr_->setContinuous(colNumber);
3859    }
3860  }
3861}
3862//-----------------------------------------------------------------------------
3863void OsiClpSolverInterface::setInteger(const int *indices, int len)
3864{
3865  if (!integerInformation_) {
3866    integerInformation_ = new char[modelPtr_->numberColumns()];
3867    CoinFillN(integerInformation_, modelPtr_->numberColumns(), static_cast< char >(0));
3868  }
3869#ifndef NDEBUG
3870  int n = modelPtr_->numberColumns();
3871#endif
3872  int i;
3873  for (i = 0; i < len; i++) {
3874    int colNumber = indices[i];
3875#ifndef NDEBUG
3876    if (colNumber < 0 || colNumber >= n) {
3877      indexError(colNumber, "setInteger");
3878    }
3879#endif
3880    integerInformation_[colNumber] = 1;
3881    modelPtr_->setInteger(colNumber);
3882  }
3883}
3884/* Set the objective coefficients for all columns
3885    array [getNumCols()] is an array of values for the objective.
3886    This defaults to a series of set operations and is here for speed.
3887*/
3888void OsiClpSolverInterface::setObjective(const double *array)
3889{
3890  // Say can't gurantee optimal basis etc
3891  lastAlgorithm_ = 999;
3892  modelPtr_->whatsChanged_ &= (0xffff & ~64);
3893  int n = modelPtr_->numberColumns();
3894  if (fakeMinInSimplex_) {
3895    std::transform(array, array + n,
3896      modelPtr_->objective(), std::negate< double >());
3897  } else {
3898    CoinMemcpyN(array, n, modelPtr_->objective());
3899  }
3900}
3901/* Set the lower bounds for all columns
3902    array [getNumCols()] is an array of values for the objective.
3903    This defaults to a series of set operations and is here for speed.
3904*/
3905void OsiClpSolverInterface::setColLower(const double *array)
3906{
3907  // Say can't gurantee optimal basis etc
3908  lastAlgorithm_ = 999;
3909  modelPtr_->whatsChanged_ &= (0x1ffff & 128);
3910  CoinMemcpyN(array, modelPtr_->numberColumns(),
3911    modelPtr_->columnLower());
3912}
3913/* Set the upper bounds for all columns
3914    array [getNumCols()] is an array of values for the objective.
3915    This defaults to a series of set operations and is here for speed.
3916*/
3917void OsiClpSolverInterface::setColUpper(const double *array)
3918{
3919  // Say can't gurantee optimal basis etc
3920  lastAlgorithm_ = 999;
3921  modelPtr_->whatsChanged_ &= (0x1ffff & 256);
3922  CoinMemcpyN(array, modelPtr_->numberColumns(),
3923    modelPtr_->columnUpper());
3924}
3925//-----------------------------------------------------------------------------
3926void OsiClpSolverInterface::setColSolution(const double *cs)
3927{
3928  // Say can't gurantee optimal basis etc
3929  lastAlgorithm_ = 999;
3930  CoinDisjointCopyN(cs, modelPtr_->numberColumns(),
3931    modelPtr_->primalColumnSolution());
3932  if (modelPtr_->solveType() == 2) {
3933    // directly into code as well
3934    CoinDisjointCopyN(cs, modelPtr_->numberColumns(),
3935      modelPtr_->solutionRegion(1));
3936  }
3937  // compute row activity
3938  memset(modelPtr_->primalRowSolution(), 0, modelPtr_->numberRows() * sizeof(double));
3939  modelPtr_->times(1.0, modelPtr_->primalColumnSolution(), modelPtr_->primalRowSolution());
3940}
3941//-----------------------------------------------------------------------------
3942void OsiClpSolverInterface::setRowPrice(const double *rs)
3943{
3944  CoinDisjointCopyN(rs, modelPtr_->numberRows(),
3945    modelPtr_->dualRowSolution());
3946  if (modelPtr_->solveType() == 2) {
3947    // directly into code as well (? sign )
3948    CoinDisjointCopyN(rs, modelPtr_->numberRows(),
3949      modelPtr_->djRegion(0));
3950  }
3951  // compute reduced costs
3952  memcpy(modelPtr_->dualColumnSolution(), modelPtr_->objective(),
3953    modelPtr_->numberColumns() * sizeof(double));
3954  modelPtr_->transposeTimes(-1.0, modelPtr_->dualRowSolution(), modelPtr_->dualColumnSolution());
3955}
3956
3957//#############################################################################
3958// Problem modifying methods (matrix)
3959//#############################################################################
3960void OsiClpSolverInterface::addCol(const CoinPackedVectorBase &vec,
3961  const double collb, const double colub,
3962  const double obj)
3963{
3964  int numberColumns = modelPtr_->numberColumns();
3965  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 8 | 64 | 128 | 256));
3966  modelPtr_->resize(modelPtr_->numberRows(), numberColumns + 1);
3967  linearObjective_ = modelPtr_->objective();
3968  basis_.resize(modelPtr_->numberRows(), numberColumns + 1);
3969  setColBounds(numberColumns, collb, colub);
3970  setObjCoeff(numberColumns, obj);
3971  if (!modelPtr_->clpMatrix())
3972    modelPtr_->createEmptyMatrix();
3973  modelPtr_->matrix()->appendCol(vec);
3974  if (integerInformation_) {
3975    char *temp = new char[numberColumns + 1];
3976    CoinMemcpyN(integerInformation_, numberColumns, temp);
3977    delete[] integerInformation_;
3978    integerInformation_ = temp;
3979    integerInformation_[numberColumns] = 0;
3980  }
3981  freeCachedResults();
3982}
3983//-----------------------------------------------------------------------------
3984/* Add a column (primal variable) to the problem. */
3985void OsiClpSolverInterface::addCol(int numberElements, const int *rows, const double *elements,
3986  const double collb, const double colub,
3987  const double obj)
3988{
3989  CoinPackedVector column(numberElements, rows, elements);
3990  addCol(column, collb, colub, obj);
3991}
3992// Add a named column (primal variable) to the problem.
3993void OsiClpSolverInterface::addCol(const CoinPackedVectorBase &vec,
3994  const double collb, const double colub,
3995  const double obj, std::string name)
3996{
3997  int ndx = getNumCols();
3998  addCol(vec, collb, colub, obj);
3999  setColName(ndx, name);
4000}
4001//-----------------------------------------------------------------------------
4002void OsiClpSolverInterface::addCols(const int numcols,
4003  const CoinPackedVectorBase *const *cols,
4004  const double *collb, const double *colub,
4005  const double *obj)
4006{
4007  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 8 | 64 | 128 | 256));
4008  int numberColumns = modelPtr_->numberColumns();
4009  modelPtr_->resize(modelPtr_->numberRows(), numberColumns + numcols);
4010  linearObjective_ = modelPtr_->objective();
4011  basis_.resize(modelPtr_->numberRows(), numberColumns + numcols);
4012  double *lower = modelPtr_->columnLower() + numberColumns;
4013  double *upper = modelPtr_->columnUpper() + numberColumns;
4014  double *objective = modelPtr_->objective() + numberColumns;
4015  int iCol;
4016  if (collb) {
4017    for (iCol = 0; iCol < numcols; iCol++) {
4018      lower[iCol] = forceIntoRange(collb[iCol], -OsiClpInfinity, OsiClpInfinity);
4019      if (lower[iCol] < -1.0e27)
4020        lower[iCol] = -COIN_DBL_MAX;
4021    }
4022  } else {
4023    CoinFillN(lower, numcols, 0.0);
4024  }
4025  if (colub) {
4026    for (iCol = 0; iCol < numcols; iCol++) {
4027      upper[iCol] = forceIntoRange(colub[iCol], -OsiClpInfinity, OsiClpInfinity);
4028      if (upper[iCol] > 1.0e27)
4029        upper[iCol] = COIN_DBL_MAX;
4030    }
4031  } else {
4032    CoinFillN(upper, numcols, COIN_DBL_MAX);
4033  }
4034  if (obj) {
4035    for (iCol = 0; iCol < numcols; iCol++) {
4036      objective[iCol] = obj[iCol];
4037    }
4038  } else {
4039    CoinFillN(objective, numcols, 0.0);
4040  }
4041  if (!modelPtr_->clpMatrix())
4042    modelPtr_->createEmptyMatrix();
4043  modelPtr_->matrix()->appendCols(numcols, cols);
4044  if (integerInformation_) {
4045    char *temp = new char[numberColumns + numcols];
4046    CoinMemcpyN(integerInformation_, numberColumns, temp);
4047    delete[] integerInformation_;
4048    integerInformation_ = temp;
4049    for (iCol = 0; iCol < numcols; iCol++)
4050      integerInformation_[numberColumns + iCol] = 0;
4051  }
4052  freeCachedResults();
4053}
4054void OsiClpSolverInterface::addCols(const int numcols,
4055  const CoinBigIndex *columnStarts, const int *rows, const double *elements,
4056  const double *collb, const double *colub,
4057  const double *obj)
4058{
4059  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 8 | 64 | 128 | 256));
4060  int numberColumns = modelPtr_->numberColumns();
4061  modelPtr_->resize(modelPtr_->numberRows(), numberColumns + numcols);
4062  linearObjective_ = modelPtr_->objective();
4063  basis_.resize(modelPtr_->numberRows(), numberColumns + numcols);
4064  double *lower = modelPtr_->columnLower() + numberColumns;
4065  double *upper = modelPtr_->columnUpper() + numberColumns;
4066  double *objective = modelPtr_->objective() + numberColumns;
4067  int iCol;
4068  if (collb) {
4069    for (iCol = 0; iCol < numcols; iCol++) {
4070      lower[iCol] = forceIntoRange(collb[iCol], -OsiClpInfinity, OsiClpInfinity);
4071      if (lower[iCol] < -1.0e27)
4072        lower[iCol] = -COIN_DBL_MAX;
4073    }
4074  } else {
4075    CoinFillN(lower, numcols, 0.0);
4076  }
4077  if (colub) {
4078    for (iCol = 0; iCol < numcols; iCol++) {
4079      upper[iCol] = forceIntoRange(colub[iCol], -OsiClpInfinity, OsiClpInfinity);
4080      if (upper[iCol] > 1.0e27)
4081        upper[iCol] = COIN_DBL_MAX;
4082    }
4083  } else {
4084    CoinFillN(upper, numcols, COIN_DBL_MAX);
4085  }
4086  if (obj) {
4087    for (iCol = 0; iCol < numcols; iCol++) {
4088      objective[iCol] = obj[iCol];
4089    }
4090  } else {
4091    CoinFillN(objective, numcols, 0.0);
4092  }
4093  if (!modelPtr_->clpMatrix())
4094    modelPtr_->createEmptyMatrix();
4095  modelPtr_->matrix()->appendCols(numcols, columnStarts, rows, elements);
4096  if (integerInformation_) {
4097    char *temp = new char[numberColumns + numcols];
4098    CoinMemcpyN(integerInformation_, numberColumns, temp);
4099    delete[] integerInformation_;
4100    integerInformation_ = temp;
4101    for (iCol = 0; iCol < numcols; iCol++)
4102      integerInformation_[numberColumns + iCol] = 0;
4103  }
4104  freeCachedResults();
4105}
4106//-----------------------------------------------------------------------------
4107void OsiClpSolverInterface::deleteCols(const int num, const int *columnIndices)
4108{
4109  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 8 | 64 | 128 | 256));
4110  findIntegers(false);
4111  deleteBranchingInfo(num, columnIndices);
4112  modelPtr_->deleteColumns(num, columnIndices);
4113  int nameDiscipline;
4114  getIntParam(OsiNameDiscipline, nameDiscipline);
4115  if (num && nameDiscipline) {
4116    // Very clumsy (and inefficient) - need to sort and then go backwards in ? chunks
4117    int *indices = CoinCopyOfArray(columnIndices, num);
4118    std::sort(indices, indices + num);
4119    int num2 = num;
4120    while (num2) {
4121      int next = indices[num2 - 1];
4122      int firstDelete = num2 - 1;
4123      int i;
4124      for (i = num2 - 2; i >= 0; i--) {
4125        if (indices[i] + 1 == next) {
4126          next--;
4127          firstDelete = i;
4128        } else {
4129          break;
4130        }
4131      }
4132      OsiSolverInterface::deleteColNames(indices[firstDelete], num2 - firstDelete);
4133      num2 = firstDelete;
4134      assert(num2 >= 0);
4135    }
4136    delete[] indices;
4137  }
4138  // synchronize integers (again)
4139  if (integerInformation_) {
4140    int numberColumns = modelPtr_->numberColumns();
4141    for (int i = 0; i < numberColumns; i++) {
4142      if (modelPtr_->isInteger(i))
4143        integerInformation_[i] = 1;
4144      else
4145        integerInformation_[i] = 0;
4146    }
4147  }
4148  basis_.deleteColumns(num, columnIndices);
4149  linearObjective_ = modelPtr_->objective();
4150  freeCachedResults();
4151}
4152//-----------------------------------------------------------------------------
4153void OsiClpSolverInterface::addRow(const CoinPackedVectorBase &vec,
4154  const double rowlb, const double rowub)
4155{
4156  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4157  freeCachedResults0();
4158  int numberRows = modelPtr_->numberRows();
4159  modelPtr_->resize(numberRows + 1, modelPtr_->numberColumns());
4160  basis_.resize(numberRows + 1, modelPtr_->numberColumns());
4161  setRowBounds(numberRows, rowlb, rowub);
4162  if (!modelPtr_->clpMatrix())
4163    modelPtr_->createEmptyMatrix();
4164  modelPtr_->matrix()->appendRow(vec);
4165  freeCachedResults1();
4166}
4167//-----------------------------------------------------------------------------
4168void OsiClpSolverInterface::addRow(const CoinPackedVectorBase &vec,
4169  const double rowlb, const double rowub,
4170  std::string name)
4171{
4172  int ndx = getNumRows();
4173  addRow(vec, rowlb, rowub);
4174  setRowName(ndx, name);
4175}
4176//-----------------------------------------------------------------------------
4177void OsiClpSolverInterface::addRow(const CoinPackedVectorBase &vec,
4178  const char rowsen, const double rowrhs,
4179  const double rowrng)
4180{
4181  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4182  freeCachedResults0();
4183  int numberRows = modelPtr_->numberRows();
4184  modelPtr_->resize(numberRows + 1, modelPtr_->numberColumns());
4185  basis_.resize(numberRows + 1, modelPtr_->numberColumns());
4186  double rowlb = 0, rowub = 0;
4187  convertSenseToBound(rowsen, rowrhs, rowrng, rowlb, rowub);
4188  setRowBounds(numberRows, rowlb, rowub);
4189  if (!modelPtr_->clpMatrix())
4190    modelPtr_->createEmptyMatrix();
4191  modelPtr_->matrix()->appendRow(vec);
4192  freeCachedResults1();
4193}
4194//-----------------------------------------------------------------------------
4195void OsiClpSolverInterface::addRow(int numberElements, const int *columns, const double *elements,
4196  const double rowlb, const double rowub)
4197{
4198  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4199  freeCachedResults0();
4200  int numberRows = modelPtr_->numberRows();
4201  modelPtr_->resize(numberRows + 1, modelPtr_->numberColumns());
4202  basis_.resize(numberRows + 1, modelPtr_->numberColumns());
4203  setRowBounds(numberRows, rowlb, rowub);
4204  if (!modelPtr_->clpMatrix())
4205    modelPtr_->createEmptyMatrix();
4206  modelPtr_->matrix()->appendRow(numberElements, columns, elements);
4207  CoinBigIndex starts[2];
4208  starts[0] = 0;
4209  starts[1] = numberElements;
4210  redoScaleFactors(1, starts, columns, elements);
4211  freeCachedResults1();
4212}
4213//-----------------------------------------------------------------------------
4214void OsiClpSolverInterface::addRows(const int numrows,
4215  const CoinPackedVectorBase *const *rows,
4216  const double *rowlb, const double *rowub)
4217{
4218  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4219  freeCachedResults0();
4220  int numberRows = modelPtr_->numberRows();
4221  modelPtr_->resize(numberRows + numrows, modelPtr_->numberColumns());
4222  basis_.resize(numberRows + numrows, modelPtr_->numberColumns());
4223  double *lower = modelPtr_->rowLower() + numberRows;
4224  double *upper = modelPtr_->rowUpper() + numberRows;
4225  int iRow;
4226  for (iRow = 0; iRow < numrows; iRow++) {
4227    if (rowlb)
4228      lower[iRow] = forceIntoRange(rowlb[iRow], -OsiClpInfinity, OsiClpInfinity);
4229    else
4230      lower[iRow] = -OsiClpInfinity;
4231    if (rowub)
4232      upper[iRow] = forceIntoRange(rowub[iRow], -OsiClpInfinity, OsiClpInfinity);
4233    else
4234      upper[iRow] = OsiClpInfinity;
4235    if (lower[iRow] < -1.0e27)
4236      lower[iRow] = -COIN_DBL_MAX;
4237    if (upper[iRow] > 1.0e27)
4238      upper[iRow] = COIN_DBL_MAX;
4239  }
4240  if (!modelPtr_->clpMatrix())
4241    modelPtr_->createEmptyMatrix();
4242  modelPtr_->matrix()->appendRows(numrows, rows);
4243  freeCachedResults1();
4244}
4245//-----------------------------------------------------------------------------
4246void OsiClpSolverInterface::addRows(const int numrows,
4247  const CoinPackedVectorBase *const *rows,
4248  const char *rowsen, const double *rowrhs,
4249  const double *rowrng)
4250{
4251  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4252  freeCachedResults0();
4253  int numberRows = modelPtr_->numberRows();
4254  modelPtr_->resize(numberRows + numrows, modelPtr_->numberColumns());
4255  basis_.resize(numberRows + numrows, modelPtr_->numberColumns());
4256  double *lower = modelPtr_->rowLower() + numberRows;
4257  double *upper = modelPtr_->rowUpper() + numberRows;
4258  int iRow;
4259  for (iRow = 0; iRow < numrows; iRow++) {
4260    double rowlb = 0, rowub = 0;
4261    convertSenseToBound(rowsen[iRow], rowrhs[iRow], rowrng[iRow],
4262      rowlb, rowub);
4263    lower[iRow] = forceIntoRange(rowlb, -OsiClpInfinity, OsiClpInfinity);
4264    upper[iRow] = forceIntoRange(rowub, -OsiClpInfinity, OsiClpInfinity);
4265    if (lower[iRow] < -1.0e27)
4266      lower[iRow] = -COIN_DBL_MAX;
4267    if (upper[iRow] > 1.0e27)
4268      upper[iRow] = COIN_DBL_MAX;
4269  }
4270  if (!modelPtr_->clpMatrix())
4271    modelPtr_->createEmptyMatrix();
4272  modelPtr_->matrix()->appendRows(numrows, rows);
4273  freeCachedResults1();
4274}
4275void OsiClpSolverInterface::addRows(const int numrows,
4276  const CoinBigIndex *rowStarts, const int *columns, const double *element,
4277  const double *rowlb, const double *rowub)
4278{
4279  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4280  freeCachedResults0();
4281  int numberRows = modelPtr_->numberRows();
4282  modelPtr_->resize(numberRows + numrows, modelPtr_->numberColumns());
4283  basis_.resize(numberRows + numrows, modelPtr_->numberColumns());
4284  double *lower = modelPtr_->rowLower() + numberRows;
4285  double *upper = modelPtr_->rowUpper() + numberRows;
4286  int iRow;
4287  for (iRow = 0; iRow < numrows; iRow++) {
4288    if (rowlb)
4289      lower[iRow] = forceIntoRange(rowlb[iRow], -OsiClpInfinity, OsiClpInfinity);
4290    else
4291      lower[iRow] = -OsiClpInfinity;
4292    if (rowub)
4293      upper[iRow] = forceIntoRange(rowub[iRow], -OsiClpInfinity, OsiClpInfinity);
4294    else
4295      upper[iRow] = OsiClpInfinity;
4296    if (lower[iRow] < -1.0e27)
4297      lower[iRow] = -COIN_DBL_MAX;
4298    if (upper[iRow] > 1.0e27)
4299      upper[iRow] = COIN_DBL_MAX;
4300  }
4301  if (!modelPtr_->clpMatrix())
4302    modelPtr_->createEmptyMatrix();
4303  modelPtr_->matrix()->appendRows(numrows, rowStarts, columns, element);
4304  redoScaleFactors(numrows, rowStarts, columns, element);
4305  freeCachedResults1();
4306}
4307//-----------------------------------------------------------------------------
4308void OsiClpSolverInterface::deleteRows(const int num, const int *rowIndices)
4309{
4310  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
4311  // will still be optimal if all rows basic
4312  bool allBasic = true;
4313  int numBasis = basis_.getNumArtificial();
4314  for (int i = 0; i < num; i++) {
4315    int iRow = rowIndices[i];
4316    if (iRow < numBasis) {
4317      if (basis_.getArtifStatus(iRow) != CoinWarmStartBasis::basic) {
4318        allBasic = false;
4319        break;
4320      }
4321    }
4322  }
4323  int saveAlgorithm = allBasic ? lastAlgorithm_ : 999;
4324  modelPtr_->deleteRows(num, rowIndices);
4325  int nameDiscipline;
4326  getIntParam(OsiNameDiscipline, nameDiscipline);
4327  if (num && nameDiscipline) {
4328    // Very clumsy (and inefficient) - need to sort and then go backwards in ? chunks
4329    int *indices = CoinCopyOfArray(rowIndices, num);
4330    std::sort(indices, indices + num);
4331    int num2 = num;
4332    while (num2) {
4333      int next = indices[num2 - 1];
4334      int firstDelete = num2 - 1;
4335      int i;
4336      for (i = num2 - 2; i >= 0; i--) {
4337        if (indices[i] + 1 == next) {
4338          next--;
4339          firstDelete = i;
4340        } else {
4341          break;
4342        }
4343      }
4344      OsiSolverInterface::deleteRowNames(indices[firstDelete], num2 - firstDelete);
4345      num2 = firstDelete;
4346      assert(num2 >= 0);
4347    }
4348    delete[] indices;
4349  }
4350  basis_.deleteRows(num, rowIndices);
4351  CoinPackedMatrix *saveRowCopy = matrixByRow_;
4352  matrixByRow_ = NULL;
4353  freeCachedResults();
4354  modelPtr_->setNewRowCopy(NULL);
4355  delete modelPtr_->scaledMatrix_;
4356  modelPtr_->scaledMatrix_ = NULL;
4357  if (saveRowCopy) {
4358    matrixByRow_ = saveRowCopy;
4359    matrixByRow_->deleteRows(num, rowIndices);
4360    if (matrixByRow_->getNumElements() != modelPtr_->clpMatrix()->getNumElements()) {
4361      delete matrixByRow_; // odd type matrix
4362      matrixByRow_ = NULL;
4363    }
4364  }
4365  lastAlgorithm_ = saveAlgorithm;
4366  if ((specialOptions_ & 131072) != 0)
4367    lastNumberRows_ = modelPtr_->numberRows();
4368}
4369
4370//#############################################################################
4371// Methods to input a problem
4372//#############################################################################
4373
4374void OsiClpSolverInterface::loadProblem(const CoinPackedMatrix &matrix,
4375  const double *collb, const double *colub,
4376  const double *obj,
4377  const double *rowlb, const double *rowub)
4378{
4379  modelPtr_->whatsChanged_ = 0;
4380  // Get rid of integer information (modelPtr will get rid of its copy)
4381  delete[] integerInformation_;
4382  integerInformation_ = NULL;
4383  modelPtr_->loadProblem(matrix, collb, colub, obj, rowlb, rowub);
4384  linearObjective_ = modelPtr_->objective();
4385  freeCachedResults();
4386  basis_ = CoinWarmStartBasis();
4387  if (ws_) {
4388    delete ws_;
4389    ws_ = 0;
4390  }
4391}
4392
4393//-----------------------------------------------------------------------------
4394
4395/*
4396  Expose the method that takes ClpMatrixBase. User request. Can't hurt, given
4397  the number of non-OSI methods already here.
4398*/
4399void OsiClpSolverInterface::loadProblem(const ClpMatrixBase &matrix,
4400  const double *collb,
4401  const double *colub,
4402  const double *obj,
4403  const double *rowlb,
4404  const double *rowub)
4405{
4406  modelPtr_->whatsChanged_ = 0;
4407  // Get rid of integer information (modelPtr will get rid of its copy)
4408  delete[] integerInformation_;
4409  integerInformation_ = NULL;
4410  modelPtr_->loadProblem(matrix, collb, colub, obj, rowlb, rowub);
4411  linearObjective_ = modelPtr_->objective();
4412  freeCachedResults();
4413  basis_ = CoinWarmStartBasis();
4414  if (ws_) {
4415    delete ws_;
4416    ws_ = 0;
4417  }
4418}
4419
4420//-----------------------------------------------------------------------------
4421
4422void OsiClpSolverInterface::assignProblem(CoinPackedMatrix *&matrix,
4423  double *&collb, double *&colub,
4424  double *&obj,
4425  double *&rowlb, double *&rowub)
4426{
4427  modelPtr_->whatsChanged_ = 0;
4428  // Get rid of integer information (modelPtr will get rid of its copy)
4429  loadProblem(*matrix, collb, colub, obj, rowlb, rowub);
4430  delete matrix;
4431  matrix = NULL;
4432  delete[] collb;
4433  collb = NULL;
4434  delete[] colub;
4435  colub = NULL;
4436  delete[] obj;
4437  obj = NULL;
4438  delete[] rowlb;
4439  rowlb = NULL;
4440  delete[] rowub;
4441  rowub = NULL;
4442}
4443
4444//-----------------------------------------------------------------------------
4445
4446void OsiClpSolverInterface::loadProblem(const CoinPackedMatrix &matrix,
4447  const double *collb, const double *colub,
4448  const double *obj,
4449  const char *rowsen, const double *rowrhs,
4450  const double *rowrng)
4451{
4452  modelPtr_->whatsChanged_ = 0;
4453  // Get rid of integer information (modelPtr will get rid of its copy)
4454  // assert( rowsen != NULL );
4455  // assert( rowrhs != NULL );
4456  // If any of Rhs NULLs then create arrays
4457  int numrows = matrix.getNumRows();
4458  const char *rowsenUse = rowsen;
4459  if (!rowsen) {
4460    char *rowsen = new char[numrows];
4461    for (int i = 0; i < numrows; i++)
4462      rowsen[i] = 'G';
4463    rowsenUse = rowsen;
4464  }
4465  const double *rowrhsUse = rowrhs;
4466  if (!rowrhs) {
4467    double *rowrhs = new double[numrows];
4468    for (int i = 0; i < numrows; i++)
4469      rowrhs[i] = 0.0;
4470    rowrhsUse = rowrhs;
4471  }
4472  const double *rowrngUse = rowrng;
4473  if (!rowrng) {
4474    double *rowrng = new double[numrows];
4475    for (int i = 0; i < numrows; i++)
4476      rowrng[i] = 0.0;
4477    rowrngUse = rowrng;
4478  }
4479  double *rowlb = new double[numrows];
4480  double *rowub = new double[numrows];
4481  for (int i = numrows - 1; i >= 0; --i) {
4482    convertSenseToBound(rowsenUse[i], rowrhsUse[i], rowrngUse[i], rowlb[i], rowub[i]);
4483  }
4484  if (rowsen != rowsenUse)
4485    delete[] rowsenUse;
4486  if (rowrhs != rowrhsUse)
4487    delete[] rowrhsUse;
4488  if (rowrng != rowrngUse)
4489    delete[] rowrngUse;
4490  loadProblem(matrix, collb, colub, obj, rowlb, rowub);
4491  delete[] rowlb;
4492  delete[] rowub;
4493}
4494
4495//-----------------------------------------------------------------------------
4496
4497void OsiClpSolverInterface::assignProblem(CoinPackedMatrix *&matrix,
4498  double *&collb, double *&colub,
4499  double *&obj,
4500  char *&rowsen, double *&rowrhs,
4501  double *&rowrng)
4502{
4503  modelPtr_->whatsChanged_ = 0;
4504  // Get rid of integer information (modelPtr will get rid of its copy)
4505  loadProblem(*matrix, collb, colub, obj, rowsen, rowrhs, rowrng);
4506  delete matrix;
4507  matrix = NULL;
4508  delete[] collb;
4509  collb = NULL;
4510  delete[] colub;
4511  colub = NULL;
4512  delete[] obj;
4513  obj = NULL;
4514  delete[] rowsen;
4515  rowsen = NULL;
4516  delete[] rowrhs;
4517  rowrhs = NULL;
4518  delete[] rowrng;
4519  rowrng = NULL;
4520}
4521
4522//-----------------------------------------------------------------------------
4523
4524void OsiClpSolverInterface::loadProblem(const int numcols, const int numrows,
4525  const CoinBigIndex *start, const int *index,
4526  const double *value,
4527  const double *collb, const double *colub,
4528  const double *obj,
4529  const double *rowlb, const double *rowub)
4530{
4531  modelPtr_->whatsChanged_ = 0;
4532  // Get rid of integer information (modelPtr will get rid of its copy)
4533  delete[] integerInformation_;
4534  integerInformation_ = NULL;
4535  modelPtr_->loadProblem(numcols, numrows, start, index,
4536    value, collb, colub, obj,
4537    rowlb, rowub);
4538  linearObjective_ = modelPtr_->objective();
4539  freeCachedResults();
4540  basis_ = CoinWarmStartBasis();
4541  if (ws_) {
4542    delete ws_;
4543    ws_ = 0;
4544  }
4545}
4546//-----------------------------------------------------------------------------
4547
4548void OsiClpSolverInterface::loadProblem(const int numcols, const int numrows,
4549  const CoinBigIndex *start, const int *index,
4550  const double *value,
4551  const double *collb, const double *colub,
4552  const double *obj,
4553  const char *rowsen, const double *rowrhs,
4554  const double *rowrng)
4555{
4556  modelPtr_->whatsChanged_ = 0;
4557  // Get rid of integer information (modelPtr will get rid of its copy)
4558  // If any of Rhs NULLs then create arrays
4559  const char *rowsenUse = rowsen;
4560  if (!rowsen) {
4561    char *rowsen = new char[numrows];
4562    for (int i = 0; i < numrows; i++)
4563      rowsen[i] = 'G';
4564    rowsenUse = rowsen;
4565  }
4566  const double *rowrhsUse = rowrhs;
4567  if (!rowrhs) {
4568    double *rowrhs = new double[numrows];
4569    for (int i = 0; i < numrows; i++)
4570      rowrhs[i] = 0.0;
4571    rowrhsUse = rowrhs;
4572  }
4573  const double *rowrngUse = rowrng;
4574  if (!rowrng) {
4575    double *rowrng = new double[numrows];
4576    for (int i = 0; i < numrows; i++)
4577      rowrng[i] = 0.0;
4578    rowrngUse = rowrng;
4579  }
4580  double *rowlb = new double[numrows];
4581  double *rowub = new double[numrows];
4582  for (int i = numrows - 1; i >= 0; --i) {
4583    convertSenseToBound(rowsenUse[i], rowrhsUse[i], rowrngUse[i], rowlb[i], rowub[i]);
4584  }
4585  if (rowsen != rowsenUse)
4586    delete[] rowsenUse;
4587  if (rowrhs != rowrhsUse)
4588    delete[] rowrhsUse;
4589  if (rowrng != rowrngUse)
4590    delete[] rowrngUse;
4591  loadProblem(numcols, numrows, start, index, value, collb, colub, obj,
4592    rowlb, rowub);
4593  delete[] rowlb;
4594  delete[] rowub;
4595}
4596// This loads a model from a coinModel object - returns number of errors
4597int OsiClpSolverInterface::loadFromCoinModel(CoinModel &modelObject, bool keepSolution)
4598{
4599  modelPtr_->whatsChanged_ = 0;
4600  int numberErrors = 0;
4601  // Set arrays for normal use
4602  double *rowLower = modelObject.rowLowerArray();
4603  double *rowUpper = modelObject.rowUpperArray();
4604  double *columnLower = modelObject.columnLowerArray();
4605  double *columnUpper = modelObject.columnUpperArray();
4606  double *objective = modelObject.objectiveArray();
4607  int *integerType = modelObject.integerTypeArray();
4608  double *associated = modelObject.associatedArray();
4609  // If strings then do copies
4610  if (modelObject.stringsExist()) {
4611    numberErrors = modelObject.createArrays(rowLower, rowUpper, columnLower, columnUpper,
4612      objective, integerType, associated);
4613  }
4614  CoinPackedMatrix matrix;
4615  modelObject.createPackedMatrix(matrix, associated);
4616  int numberRows = modelObject.numberRows();
4617  int numberColumns = modelObject.numberColumns();
4618  CoinWarmStart *ws = getWarmStart();
4619  bool restoreBasis = keepSolution && numberRows && numberRows == getNumRows() && numberColumns == getNumCols();
4620  loadProblem(matrix,
4621    columnLower, columnUpper, objective, rowLower, rowUpper);
4622  if (restoreBasis)
4623    setWarmStart(ws);
4624  delete ws;
4625  // Do names if wanted
4626  int numberItems;
4627  numberItems = modelObject.rowNames()->numberItems();
4628  if (numberItems) {
4629    const char *const *rowNames = modelObject.rowNames()->names();
4630    modelPtr_->copyRowNames(rowNames, 0, numberItems);
4631  }
4632  numberItems = modelObject.columnNames()->numberItems();
4633  if (numberItems) {
4634    const char *const *columnNames = modelObject.columnNames()->names();
4635    modelPtr_->copyColumnNames(columnNames, 0, numberItems);
4636  }
4637  // Do integers if wanted
4638  assert(integerType);
4639  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4640    if (integerType[iColumn])
4641      setInteger(iColumn);
4642  }
4643  if (rowLower != modelObject.rowLowerArray() || columnLower != modelObject.columnLowerArray()) {
4644    delete[] rowLower;
4645    delete[] rowUpper;
4646    delete[] columnLower;
4647    delete[] columnUpper;
4648    delete[] objective;
4649    delete[] integerType;
4650    delete[] associated;
4651    //if (numberErrors)
4652    //  handler_->message(CLP_BAD_STRING_VALUES,messages_)
4653    //    <<numberErrors
4654    //    <<CoinMessageEol;
4655  }
4656  modelPtr_->optimizationDirection_ = modelObject.optimizationDirection();
4657  return numberErrors;
4658}
4659
4660//-----------------------------------------------------------------------------
4661// Write mps files
4662//-----------------------------------------------------------------------------
4663
4664void OsiClpSolverInterface::writeMps(const char *filename,
4665  const char *extension,
4666  double objSense) const
4667{
4668  std::string f(filename);
4669  std::string e(extension);
4670  std::string fullname;
4671  if (e != "") {
4672    fullname = f + "." + e;
4673  } else {
4674    // no extension so no trailing period
4675    fullname = f;
4676  }
4677  // get names
4678  const char *const *const rowNames = modelPtr_->rowNamesAsChar();
4679  const char *const *const columnNames = modelPtr_->columnNamesAsChar();
4680  // Fall back on Osi version - possibly with names
4681  OsiSolverInterface::writeMpsNative(fullname.c_str(),
4682    const_cast< const char ** >(rowNames),
4683    const_cast< const char ** >(columnNames), 0, 2, objSense,
4684    numberSOS_, setInfo_);
4685  if (rowNames) {
4686    modelPtr_->deleteNamesAsChar(rowNames, modelPtr_->numberRows_ + 1);
4687    modelPtr_->deleteNamesAsChar(columnNames, modelPtr_->numberColumns_);
4688  }
4689}
4690
4691int OsiClpSolverInterface::writeMpsNative(const char *filename,
4692  const char **rowNames, const char **columnNames,
4693  int formatType, int numberAcross, double objSense) const
4694{
4695  return OsiSolverInterface::writeMpsNative(filename, rowNames, columnNames,
4696    formatType, numberAcross, objSense,
4697    numberSOS_, setInfo_);
4698}
4699
4700//#############################################################################
4701// CLP specific public interfaces
4702//#############################################################################
4703
4704ClpSimplex *OsiClpSolverInterface::getModelPtr() const
4705{
4706  int saveAlgorithm = lastAlgorithm_;
4707  //freeCachedResults();
4708  lastAlgorithm_ = saveAlgorithm;
4709  //bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
4710  return modelPtr_;
4711}
4712//-------------------------------------------------------------------
4713
4714//#############################################################################
4715// Constructors, destructors clone and assignment
4716//#############################################################################
4717//-------------------------------------------------------------------
4718// Default Constructor
4719//-------------------------------------------------------------------
4720OsiClpSolverInterface::OsiClpSolverInterface()
4721  : OsiSolverInterface()
4722  , rowsense_(NULL)
4723  , rhs_(NULL)
4724  , rowrange_(NULL)
4725  , ws_(NULL)
4726  , rowActivity_(NULL)
4727  , columnActivity_(NULL)
4728  , numberSOS_(0)
4729  , setInfo_(NULL)
4730  , smallModel_(NULL)
4731  , factorization_(NULL)
4732  , smallestElementInCut_(1.0e-15)
4733  , smallestChangeInCut_(1.0e-10)
4734  , largestAway_(-1.0)
4735  , spareArrays_(NULL)
4736  , matrixByRow_(NULL)
4737  , matrixByRowAtContinuous_(NULL)
4738  , integerInformation_(NULL)
4739  , whichRange_(NULL)
4740  , fakeMinInSimplex_(false)
4741  , linearObjective_(NULL)
4742  , cleanupScaling_(0)
4743  , specialOptions_(0x80000000)
4744  , baseModel_(NULL)
4745  , lastNumberRows_(0)
4746  , continuousModel_(NULL)
4747  , fakeObjective_(NULL)
4748{
4749  //printf("%p %d null constructor\n",this,xxxxxx);xxxxxx++;
4750  modelPtr_ = NULL;
4751  notOwned_ = false;
4752  disasterHandler_ = new OsiClpDisasterHandler();
4753  reset();
4754}
4755
4756//-------------------------------------------------------------------
4757// Clone
4758//-------------------------------------------------------------------
4759OsiSolverInterface *OsiClpSolverInterface::clone(bool CopyData) const
4760{
4761  //printf("in clone %x\n",this);
4762  OsiClpSolverInterface *newSolver;
4763  if (CopyData) {
4764    newSolver = new OsiClpSolverInterface(*this);
4765  } else {
4766    newSolver = new OsiClpSolverInterface();
4767  }
4768#if 0
4769  const double * obj = newSolver->getObjCoefficients();
4770  const double * oldObj = getObjCoefficients();
4771  if(newSolver->getNumCols()>3787)
4772    printf("%x - obj %x (from %x) val %g\n",newSolver,obj,oldObj,obj[3787]);
4773#endif
4774  return newSolver;
4775}
4776
4777//-------------------------------------------------------------------
4778// Copy constructor
4779//-------------------------------------------------------------------
4780OsiClpSolverInterface::OsiClpSolverInterface(
4781  const OsiClpSolverInterface &rhs)
4782  : OsiSolverInterface(rhs)
4783  , rowsense_(NULL)
4784  , rhs_(NULL)
4785  , rowrange_(NULL)
4786  , ws_(NULL)
4787  , rowActivity_(NULL)
4788  , columnActivity_(NULL)
4789  , stuff_(rhs.stuff_)
4790  , numberSOS_(rhs.numberSOS_)
4791  , setInfo_(NULL)
4792  , smallModel_(NULL)
4793  , factorization_(NULL)
4794  , smallestElementInCut_(rhs.smallestElementInCut_)
4795  , smallestChangeInCut_(rhs.smallestChangeInCut_)
4796  , largestAway_(-1.0)
4797  , spareArrays_(NULL)
4798  , basis_()
4799  , itlimOrig_(9999999)
4800  , lastAlgorithm_(0)
4801  , notOwned_(false)
4802  , matrixByRow_(NULL)
4803  , matrixByRowAtContinuous_(NULL)
4804  , integerInformation_(NULL)
4805  , whichRange_(NULL)
4806  , fakeMinInSimplex_(rhs.fakeMinInSimplex_)
4807{
4808  //printf("%p %d copy constructor %p\n",this,xxxxxx,&rhs);xxxxxx++;
4809  if (rhs.modelPtr_)
4810    modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
4811  else
4812    modelPtr_ = new ClpSimplex();
4813  if (rhs.baseModel_)
4814    baseModel_ = new ClpSimplex(*rhs.baseModel_);
4815  else
4816    baseModel_ = NULL;
4817  if (rhs.continuousModel_)
4818    continuousModel_ = new ClpSimplex(*rhs.continuousModel_);
4819  else
4820    continuousModel_ = NULL;
4821  if (rhs.matrixByRowAtContinuous_)
4822    matrixByRowAtContinuous_ = new CoinPackedMatrix(*rhs.matrixByRowAtContinuous_);
4823  if (rhs.disasterHandler_)
4824    disasterHandler_ = dynamic_cast< OsiClpDisasterHandler * >(rhs.disasterHandler_->clone());
4825  else
4826    disasterHandler_ = NULL;
4827  if (rhs.fakeObjective_)
4828    fakeObjective_ = new ClpLinearObjective(*rhs.fakeObjective_);
4829  else
4830    fakeObjective_ = NULL;
4831  linearObjective_ = modelPtr_->objective();
4832  if (rhs.ws_)
4833    ws_ = new CoinWarmStartBasis(*rhs.ws_);
4834  basis_ = rhs.basis_;
4835  if (rhs.integerInformation_) {
4836    int numberColumns = modelPtr_->numberColumns();
4837    integerInformation_ = new char[numberColumns];
4838    CoinMemcpyN(rhs.integerInformation_, numberColumns, integerInformation_);
4839  }
4840  saveData_ = rhs.saveData_;
4841  solveOptions_ = rhs.solveOptions_;
4842  cleanupScaling_ = rhs.cleanupScaling_;
4843  specialOptions_ = rhs.specialOptions_;
4844  lastNumberRows_ = rhs.lastNumberRows_;
4845  rowScale_ = rhs.rowScale_;
4846  columnScale_ = rhs.columnScale_;
4847  fillParamMaps();
4848  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4849  if (numberSOS_) {
4850    setInfo_ = new CoinSet[numberSOS_];
4851    for (int i = 0; i < numberSOS_; i++)
4852      setInfo_[i] = rhs.setInfo_[i];
4853  }
4854}
4855
4856// Borrow constructor - only delete one copy
4857OsiClpSolverInterface::OsiClpSolverInterface(ClpSimplex *rhs,
4858  bool reallyOwn)
4859  : OsiSolverInterface()
4860  , rowsense_(NULL)
4861  , rhs_(NULL)
4862  , rowrange_(NULL)
4863  , ws_(NULL)
4864  , rowActivity_(NULL)
4865  , columnActivity_(NULL)
4866  , numberSOS_(0)
4867  , setInfo_(NULL)
4868  , smallModel_(NULL)
4869  , factorization_(NULL)
4870  , smallestElementInCut_(1.0e-15)
4871  , smallestChangeInCut_(1.0e-10)
4872  , largestAway_(-1.0)
4873  , spareArrays_(NULL)
4874  , basis_()
4875  , itlimOrig_(9999999)
4876  , lastAlgorithm_(0)
4877  , notOwned_(false)
4878  , matrixByRow_(NULL)
4879  , matrixByRowAtContinuous_(NULL)
4880  , integerInformation_(NULL)
4881  , whichRange_(NULL)
4882  , fakeMinInSimplex_(false)
4883  , cleanupScaling_(0)
4884  , specialOptions_(0x80000000)
4885  , baseModel_(NULL)
4886  , lastNumberRows_(0)
4887  , continuousModel_(NULL)
4888  , fakeObjective_(NULL)
4889{
4890  disasterHandler_ = new OsiClpDisasterHandler();
4891  //printf("in borrow %x - > %x\n",&rhs,this);
4892  modelPtr_ = rhs;
4893  basis_.resize(modelPtr_->numberRows(), modelPtr_->numberColumns());
4894  linearObjective_ = modelPtr_->objective();
4895  if (rhs) {
4896    notOwned_ = !reallyOwn;
4897
4898    if (rhs->integerInformation()) {
4899      int numberColumns = modelPtr_->numberColumns();
4900      integerInformation_ = new char[numberColumns];
4901      CoinMemcpyN(rhs->integerInformation(), numberColumns, integerInformation_);
4902    }
4903  }
4904  fillParamMaps();
4905}
4906
4907// Releases so won't error
4908void OsiClpSolverInterface::releaseClp()
4909{
4910  modelPtr_ = NULL;
4911  notOwned_ = false;
4912}
4913
4914//-------------------------------------------------------------------
4915// Destructor
4916//-------------------------------------------------------------------
4917OsiClpSolverInterface::~OsiClpSolverInterface()
4918{
4919  //printf("%p destructor\n",this);
4920  freeCachedResults();
4921  if (!notOwned_)
4922    delete modelPtr_;
4923  delete baseModel_;
4924  delete continuousModel_;
4925  delete disasterHandler_;
4926  delete fakeObjective_;
4927  delete ws_;
4928  delete[] rowActivity_;
4929  delete[] columnActivity_;
4930  delete[] setInfo_;
4931#ifdef KEEP_SMALL
4932  if (smallModel_) {
4933    delete[] spareArrays_;
4934    spareArrays_ = NULL;
4935    delete smallModel_;
4936    smallModel_ = NULL;
4937  }
4938#endif
4939  assert(smallModel_ == NULL);
4940  assert(factorization_ == NULL);
4941  assert(spareArrays_ == NULL);
4942  delete[] integerInformation_;
4943  delete matrixByRowAtContinuous_;
4944  delete matrixByRow_;
4945}
4946
4947//-------------------------------------------------------------------
4948// Assignment operator
4949//-------------------------------------------------------------------
4950OsiClpSolverInterface &
4951OsiClpSolverInterface::operator=(const OsiClpSolverInterface &rhs)
4952{
4953  if (this != &rhs) {
4954    //printf("in = %x - > %x\n",&rhs,this);
4955    OsiSolverInterface::operator=(rhs);
4956    freeCachedResults();
4957    if (!notOwned_)
4958      delete modelPtr_;
4959    delete ws_;
4960    if (rhs.modelPtr_)
4961      modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
4962    delete baseModel_;
4963    if (rhs.baseModel_)
4964      baseModel_ = new ClpSimplex(*rhs.baseModel_);
4965    else
4966      baseModel_ = NULL;
4967    delete continuousModel_;
4968    if (rhs.continuousModel_)
4969      continuousModel_ = new ClpSimplex(*rhs.continuousModel_);
4970    else
4971      continuousModel_ = NULL;
4972    delete matrixByRowAtContinuous_;
4973    delete matrixByRow_;
4974    matrixByRow_ = NULL;
4975    if (rhs.matrixByRowAtContinuous_)
4976      matrixByRowAtContinuous_ = new CoinPackedMatrix(*rhs.matrixByRowAtContinuous_);
4977    else
4978      matrixByRowAtContinuous_ = NULL;
4979    delete disasterHandler_;
4980    if (rhs.disasterHandler_)
4981      disasterHandler_ = dynamic_cast< OsiClpDisasterHandler * >(rhs.disasterHandler_->clone());
4982    else
4983      disasterHandler_ = NULL;
4984    delete fakeObjective_;
4985    if (rhs.fakeObjective_)
4986      fakeObjective_ = new ClpLinearObjective(*rhs.fakeObjective_);
4987    else
4988      fakeObjective_ = NULL;
4989    notOwned_ = false;
4990    linearObjective_ = modelPtr_->objective();
4991    saveData_ = rhs.saveData_;
4992    solveOptions_ = rhs.solveOptions_;
4993    cleanupScaling_ = rhs.cleanupScaling_;
4994    specialOptions_ = rhs.specialOptions_;
4995    lastNumberRows_ = rhs.lastNumberRows_;
4996    rowScale_ = rhs.rowScale_;
4997    columnScale_ = rhs.columnScale_;
4998    basis_ = rhs.basis_;
4999    stuff_ = rhs.stuff_;
5000    delete[] integerInformation_;
5001    integerInformation_ = NULL;
5002    if (rhs.integerInformation_) {
5003      int numberColumns = modelPtr_->numberColumns();
5004      integerInformation_ = new char[numberColumns];
5005      CoinMemcpyN(rhs.integerInformation_, numberColumns, integerInformation_);
5006    }
5007    if (rhs.ws_)
5008      ws_ = new CoinWarmStartBasis(*rhs.ws_);
5009    else
5010      ws_ = NULL;
5011    delete[] rowActivity_;
5012    delete[] columnActivity_;
5013    rowActivity_ = NULL;
5014    columnActivity_ = NULL;
5015    delete[] setInfo_;
5016    numberSOS_ = rhs.numberSOS_;
5017    setInfo_ = NULL;
5018    if (numberSOS_) {
5019      setInfo_ = new CoinSet[numberSOS_];
5020      for (int i = 0; i < numberSOS_; i++)
5021        setInfo_[i] = rhs.setInfo_[i];
5022    }
5023    assert(smallModel_ == NULL);
5024    assert(factorization_ == NULL);
5025    smallestElementInCut_ = rhs.smallestElementInCut_;
5026    smallestChangeInCut_ = rhs.smallestChangeInCut_;
5027    largestAway_ = -1.0;
5028    assert(spareArrays_ == NULL);
5029    basis_ = rhs.basis_;
5030    fillParamMaps();
5031    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
5032  }
5033  return *this;
5034}
5035
5036//#############################################################################
5037// Applying cuts
5038//#############################################################################
5039
5040void OsiClpSolverInterface::applyRowCut(const OsiRowCut &rowCut)
5041{
5042  applyRowCuts(1, &rowCut);
5043}
5044/* Apply a collection of row cuts which are all effective.
5045   applyCuts seems to do one at a time which seems inefficient.
5046*/
5047void OsiClpSolverInterface::applyRowCuts(int numberCuts, const OsiRowCut *cuts)
5048{
5049  if (numberCuts) {
5050    // Say can't gurantee optimal basis etc
5051    lastAlgorithm_ = 999;
5052
5053    // Thanks to js
5054    const OsiRowCut **cutsp = new const OsiRowCut *[numberCuts];
5055    for (int i = 0; i < numberCuts; i++)
5056      cutsp[i] = &cuts[i];
5057
5058    applyRowCuts(numberCuts, cutsp);
5059
5060    delete[] cutsp;
5061  }
5062}
5063/* Apply a collection of row cuts which are all effective.
5064   applyCuts seems to do one at a time which seems inefficient.
5065*/
5066void OsiClpSolverInterface::applyRowCuts(int numberCuts, const OsiRowCut **cuts)
5067{
5068  int i;
5069  if (!numberCuts)
5070    return;
5071  modelPtr_->whatsChanged_ &= (0xffff & ~(1 | 2 | 4 | 16 | 32));
5072  CoinPackedMatrix *saveRowCopy = matrixByRow_;
5073  matrixByRow_ = NULL;
5074#if 0 // was #ifndef NDEBUG
5075  int nameDiscipline;
5076  getIntParam(OsiNameDiscipline,nameDiscipline) ;
5077  assert (!nameDiscipline);
5078#endif
5079  freeCachedResults0();
5080  // Say can't gurantee optimal basis etc
5081  lastAlgorithm_ = 999;
5082  int numberRows = modelPtr_->numberRows();
5083  modelPtr_->resize(numberRows + numberCuts, modelPtr_->numberColumns());
5084  basis_.resize(numberRows + numberCuts, modelPtr_->numberColumns());
5085  // redo as relaxed - use modelPtr_-> addRows with starts etc
5086  int size = 0;
5087  for (i = 0; i < numberCuts; i++)
5088    size += cuts[i]->row().getNumElements();
5089  CoinBigIndex *starts = new CoinBigIndex[numberCuts + 1];
5090  int *indices = new int[size];
5091  double *elements = new double[size];
5092  double *lower = modelPtr_->rowLower() + numberRows;
5093  double *upper = modelPtr_->rowUpper() + numberRows;
5094  const double *columnLower = modelPtr_->columnLower();
5095  const double *columnUpper = modelPtr_->columnUpper();
5096  size = 0;
5097  for (i = 0; i < numberCuts; i++) {
5098    double rowLb = cuts[i]->lb();
5099    double rowUb = cuts[i]->ub();
5100    int n = cuts[i]->row().getNumElements();
5101    const int *index = cuts[i]->row().getIndices();
5102    const double *elem = cuts[i]->row().getElements();
5103    starts[i] = size;
5104    for (int j = 0; j < n; j++) {
5105      double value = elem[j];
5106      int column = index[j];
5107      if (fabs(value) >= smallestChangeInCut_) {
5108        // always take
5109        indices[size] = column;
5110        elements[size++] = value;
5111      } else if (fabs(value) >= smallestElementInCut_) {
5112        double lowerValue = columnLower[column];
5113        double upperValue = columnUpper[column];
5114        double difference = upperValue - lowerValue;
5115        if (difference < 1.0e20 && difference * fabs(value) < smallestChangeInCut_ && (rowLb < -1.0e20 || rowUb > 1.0e20)) {
5116          // Take out and adjust to relax
5117          //printf("small el %g adjusted\n",value);
5118          if (rowLb > -1.0e20) {
5119            // just lower bound on row
5120            if (value > 0.0) {
5121              // pretend at upper
5122              rowLb -= value * upperValue;
5123            } else {
5124              // pretend at lower
5125              rowLb -= value * lowerValue;
5126            }
5127          } else {
5128            // just upper bound on row
5129            if (value > 0.0) {
5130              // pretend at lower
5131              rowUb -= value * lowerValue;
5132            } else {
5133              // pretend at upper
5134              rowUb -= value * upperValue;
5135            }
5136          }
5137        } else {
5138          // take (unwillingly)
5139          indices[size] = column;
5140          elements[size++] = value;
5141        }
5142      } else {
5143        //printf("small el %g ignored\n",value);
5144      }
5145    }
5146    lower[i] = forceIntoRange(rowLb, -OsiClpInfinity, OsiClpInfinity);
5147    upper[i] = forceIntoRange(rowUb, -OsiClpInfinity, OsiClpInfinity);
5148    if (lower[i] < -1.0e27)
5149      lower[i] = -COIN_DBL_MAX;
5150    if (upper[i] > 1.0e27)
5151      upper[i] = COIN_DBL_MAX;
5152  }
5153  starts[numberCuts] = size;
5154  if (!modelPtr_->clpMatrix())
5155    modelPtr_->createEmptyMatrix();
5156  //modelPtr_->matrix()->appendRows(numberCuts,rows);
5157  modelPtr_->clpMatrix()->appendMatrix(numberCuts, 0, starts, indices, elements);
5158  modelPtr_->setNewRowCopy(NULL);
5159  modelPtr_->setClpScaledMatrix(NULL);
5160  freeCachedResults1();
5161  redoScaleFactors(numberCuts, starts, indices, elements);
5162  if (saveRowCopy) {
5163#if 1
5164    matrixByRow_ = saveRowCopy;
5165    matrixByRow_->appendRows(numberCuts, starts, indices, elements, 0);
5166    if (matrixByRow_->getNumElements() != modelPtr_->clpMatrix()->getNumElements()) {
5167      delete matrixByRow_; // odd type matrix
5168      matrixByRow_ = NULL;
5169    }
5170#else
5171    delete saveRowCopy;
5172#endif
5173  }
5174  delete[] starts;
5175  delete[] indices;
5176  delete[] elements;
5177}
5178//#############################################################################
5179// Apply Cuts
5180//#############################################################################
5181
5182OsiSolverInterface::ApplyCutsReturnCode
5183OsiClpSolverInterface::applyCuts(const OsiCuts &cs, double effectivenessLb)
5184{
5185  OsiSolverInterface::ApplyCutsReturnCode retVal;
5186  int i;
5187
5188  // Loop once for each column cut
5189  for (i = 0; i < cs.sizeColCuts(); i++) {
5190    if (cs.colCut(i).effectiveness() < effectivenessLb) {
5191      retVal.incrementIneffective();
5192      continue;
5193    }
5194    if (!cs.colCut(i).consistent()) {
5195      retVal.incrementInternallyInconsistent();
5196      continue;
5197    }
5198    if (!cs.colCut(i).consistent(*this)) {
5199      retVal.incrementExternallyInconsistent();
5200      continue;
5201    }
5202    if (cs.colCut(i).infeasible(*this)) {
5203      retVal.incrementInfeasible();
5204      continue;
5205    }
5206    applyColCut(cs.colCut(i));
5207    retVal.incrementApplied();
5208  }
5209
5210  // Loop once for each row cut
5211  const OsiRowCut **addCuts = new const OsiRowCut *[cs.sizeRowCuts()];
5212  int nAdd = 0;
5213  for (i = 0; i < cs.sizeRowCuts(); i++) {
5214    if (cs.rowCut(i).effectiveness() < effectivenessLb) {
5215      retVal.incrementIneffective();
5216      continue;
5217    }
5218    if (!cs.rowCut(i).consistent()) {
5219      retVal.incrementInternallyInconsistent();
5220      continue;
5221    }
5222    if (!cs.rowCut(i).consistent(*this)) {
5223      retVal.incrementExternallyInconsistent();
5224      continue;
5225    }
5226    if (cs.rowCut(i).infeasible(*this)) {
5227      retVal.incrementInfeasible();
5228      continue;
5229    }
5230    addCuts[nAdd++] = cs.rowCutPtr(i);
5231    retVal.incrementApplied();
5232  }
5233  // now apply
5234  applyRowCuts(nAdd, addCuts);
5235  delete[] addCuts;
5236
5237  return retVal;
5238}
5239// Extend scale factors
5240void OsiClpSolverInterface::redoScaleFactors(int numberAdd, const CoinBigIndex *starts,
5241  const int *indices, const double *elements)
5242{
5243  if ((specialOptions_ & 131072) != 0) {
5244    int numberRows = modelPtr_->numberRows() - numberAdd;
5245    assert(lastNumberRows_ == numberRows); // ???
5246    int iRow;
5247    int newNumberRows = numberRows + numberAdd;
5248    rowScale_.extend(static_cast< int >(2 * newNumberRows * sizeof(double)));
5249    double *rowScale = rowScale_.array();
5250    double *oldInverseScale = rowScale + lastNumberRows_;
5251    double *inverseRowScale = rowScale + newNumberRows;
5252    for (iRow = lastNumberRows_ - 1; iRow >= 0; iRow--)
5253      inverseRowScale[iRow] = oldInverseScale[iRow];
5254    //int numberColumns = baseModel_->numberColumns();
5255    const double *columnScale = columnScale_.array();
5256    //const double * inverseColumnScale = columnScale + numberColumns;
5257    // Geometric mean on row scales
5258    // adjust arrays
5259    rowScale += lastNumberRows_;
5260    inverseRowScale += lastNumberRows_;
5261    for (iRow = 0; iRow < numberAdd; iRow++) {
5262      CoinBigIndex j;
5263      double largest = 1.0e-20;
5264      double smallest = 1.0e50;
5265      for (j = starts[iRow]; j < starts[iRow + 1]; j++) {
5266        int iColumn = indices[j];
5267        double value = fabs(elements[j]);
5268        // Don't bother with tiny elements
5269        if (value > 1.0e-20) {
5270          value *= columnScale[iColumn];
5271          largest = CoinMax(largest, value);
5272          smallest = CoinMin(smallest, value);
5273        }
5274      }
5275      double scale = sqrt(smallest * largest);
5276      scale = CoinMax(1.0e-10, CoinMin(1.0e10, scale));
5277      inverseRowScale[iRow] = scale;
5278      rowScale[iRow] = 1.0 / scale;
5279    }
5280    lastNumberRows_ = newNumberRows;
5281  }
5282}
5283// Delete all scale factor stuff and reset option
5284void OsiClpSolverInterface::deleteScaleFactors()
5285{
5286  delete baseModel_;
5287  baseModel_ = NULL;
5288  lastNumberRows_ = 0;
5289  specialOptions_ &= ~131072;
5290}
5291//-----------------------------------------------------------------------------
5292
5293void OsiClpSolverInterface::applyColCut(const OsiColCut &cc)
5294{
5295  modelPtr_->whatsChanged_ &= (0x1ffff & ~(128 | 256));
5296  // Say can't gurantee optimal basis etc
5297  lastAlgorithm_ = 999;
5298  double *lower = modelPtr_->columnLower();
5299  double *upper = modelPtr_->columnUpper();
5300  const CoinPackedVector &lbs = cc.lbs();
5301  const CoinPackedVector &ubs = cc.ubs();
5302  int i;
5303
5304  for (i = 0; i < lbs.getNumElements(); i++) {
5305    int iCol = lbs.getIndices()[i];
5306    double value = lbs.getElements()[i];
5307    if (value > lower[iCol])
5308      lower[iCol] = value;
5309  }
5310  for (i = 0; i < ubs.getNumElements(); i++) {
5311    int iCol = ubs.getIndices()[i];
5312    double value = ubs.getElements()[i];
5313    if (value < upper[iCol])
5314      upper[iCol] = value;
5315  }
5316}
5317//#############################################################################
5318// Private methods
5319//#############################################################################
5320
5321//-------------------------------------------------------------------
5322
5323void OsiClpSolverInterface::freeCachedResults() const
5324{
5325  // Say can't gurantee optimal basis etc
5326  lastAlgorithm_ = 999;
5327  delete[] rowsense_;
5328  delete[] rhs_;
5329  delete[] rowrange_;
5330  delete matrixByRow_;
5331  //delete ws_;
5332  rowsense_ = NULL;
5333  rhs_ = NULL;
5334  rowrange_ = NULL;
5335  matrixByRow_ = NULL;
5336  //ws_ = NULL;
5337  if (!notOwned_ && modelPtr_) {
5338    if (modelPtr_->scaledMatrix_) {
5339      delete modelPtr_->scaledMatrix_;
5340      modelPtr_->scaledMatrix_ = NULL;
5341    }
5342    if (modelPtr_->clpMatrix()) {
5343      modelPtr_->clpMatrix()->refresh(modelPtr_); // make sure all clean
5344#ifndef NDEBUG
5345      ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(modelPtr_->clpMatrix());
5346      if (clpMatrix) {
5347        if (clpMatrix->getNumCols())
5348          assert(clpMatrix->getNumRows() == modelPtr_->getNumRows());
5349        if (clpMatrix->getNumRows())
5350          assert(clpMatrix->getNumCols() == modelPtr_->getNumCols());
5351      }
5352#endif
5353    }
5354  }
5355}
5356
5357//-------------------------------------------------------------------
5358
5359void OsiClpSolverInterface::freeCachedResults0() const
5360{
5361  delete[] rowsense_;
5362  delete[] rhs_;
5363  delete[] rowrange_;
5364  rowsense_ = NULL;
5365  rhs_ = NULL;
5366  rowrange_ = NULL;
5367}
5368
5369//-------------------------------------------------------------------
5370
5371void OsiClpSolverInterface::freeCachedResults1() const
5372{
5373  // Say can't gurantee optimal basis etc
5374  lastAlgorithm_ = 999;
5375  delete matrixByRow_;
5376  matrixByRow_ = NULL;
5377  //ws_ = NULL;
5378  if (modelPtr_ && modelPtr_->clpMatrix()) {
5379    delete modelPtr_->scaledMatrix_;
5380    modelPtr_->scaledMatrix_ = NULL;
5381    modelPtr_->clpMatrix()->refresh(modelPtr_); // make sure all clean
5382#ifndef NDEBUG
5383    ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(modelPtr_->clpMatrix());
5384    if (clpMatrix) {
5385      assert(clpMatrix->getNumRows() == modelPtr_->getNumRows());
5386      assert(clpMatrix->getNumCols() == modelPtr_->getNumCols());
5387    }
5388#endif
5389  }
5390}
5391
5392//------------------------------------------------------------------
5393void OsiClpSolverInterface::extractSenseRhsRange() const
5394{
5395  if (rowsense_ == NULL) {
5396    // all three must be NULL
5397    assert((rhs_ == NULL) && (rowrange_ == NULL));
5398
5399    int nr = modelPtr_->numberRows();
5400    if (nr != 0) {
5401      rowsense_ = new char[nr];
5402      rhs_ = new double[nr];
5403      rowrange_ = new double[nr];
5404      std::fill(rowrange_, rowrange_ + nr, 0.0);
5405
5406      const double *lb = modelPtr_->rowLower();
5407      const double *ub = modelPtr_->rowUpper();
5408
5409      int i;
5410      for (i = 0; i < nr; i++) {
5411        convertBoundToSense(lb[i], ub[i], rowsense_[i], rhs_[i], rowrange_[i]);
5412      }
5413    }
5414  }
5415}
5416// Set language
5417void OsiClpSolverInterface::newLanguage(CoinMessages::Language language)
5418{
5419  modelPtr_->newLanguage(language);
5420  OsiSolverInterface::newLanguage(language);
5421}
5422//#############################################################################
5423
5424void OsiClpSolverInterface::fillParamMaps()
5425{
5426  assert(static_cast< int >(OsiMaxNumIteration) == static_cast< int >(ClpMaxNumIteration));
5427  assert(static_cast< int >(OsiMaxNumIterationHotStart) == static_cast< int >(ClpMaxNumIterationHotStart));
5428  //assert (static_cast<int> (OsiLastIntParam)==           static_cast<int>(ClpLastIntParam));
5429
5430  assert(static_cast< int >(OsiDualObjectiveLimit) == static_cast< int >(ClpDualObjectiveLimit));
5431  assert(static_cast< int >(OsiPrimalObjectiveLimit) == static_cast< int >(ClpPrimalObjectiveLimit));
5432  assert(static_cast< int >(OsiDualTolerance) == static_cast< int >(ClpDualTolerance));
5433  assert(static_cast< int >(OsiPrimalTolerance) == static_cast< int >(ClpPrimalTolerance));
5434  assert(static_cast< int >(OsiObjOffset) == static_cast< int >(ClpObjOffset));
5435  //assert (static_cast<int> (OsiLastDblParam)==        static_cast<int>(ClpLastDblParam));
5436
5437  assert(static_cast< int >(OsiProbName) == static_cast< int >(ClpProbName));
5438  //strParamMap_[OsiLastStrParam] = ClpLastStrParam;
5439}
5440// Sets up basis
5441void OsiClpSolverInterface::setBasis(const CoinWarmStartBasis &basis)
5442{
5443  setBasis(basis, modelPtr_);
5444  setWarmStart(&basis);
5445}
5446// Warm start
5447CoinWarmStartBasis
5448OsiClpSolverInterface::getBasis(ClpSimplex *model) const
5449{
5450  int iRow, iColumn;
5451  int numberRows = model->numberRows();
5452  int numberColumns = model->numberColumns();
5453  CoinWarmStartBasis basis;
5454  basis.setSize(numberColumns, numberRows);
5455  if (model->statusExists()) {
5456    // Flip slacks
5457    int lookupA[] = { 0, 1, 3, 2, 0, 2 };
5458    for (iRow = 0; iRow < numberRows; iRow++) {
5459      int iStatus = model->getRowStatus(iRow);
5460      iStatus = lookupA[iStatus];
5461      basis.setArtifStatus(iRow, static_cast< CoinWarmStartBasis::Status >(iStatus));
5462    }
5463    int lookupS[] = { 0, 1, 2, 3, 0, 3 };
5464    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5465      int iStatus = model->getColumnStatus(iColumn);
5466      iStatus = lookupS[iStatus];
5467      basis.setStructStatus(iColumn, static_cast< CoinWarmStartBasis::Status >(iStatus));
5468    }
5469  }
5470  //basis.print();
5471  return basis;
5472}
5473// Warm start from statusArray
5474CoinWarmStartBasis *
5475OsiClpSolverInterface::getBasis(const unsigned char *statusArray) const
5476{
5477  int iRow, iColumn;
5478  int numberRows = modelPtr_->numberRows();
5479  int numberColumns = modelPtr_->numberColumns();
5480  CoinWarmStartBasis *basis = new CoinWarmStartBasis();
5481  basis->setSize(numberColumns, numberRows);
5482  // Flip slacks
5483  int lookupA[] = { 0, 1, 3, 2, 0, 2 };
5484  for (iRow = 0; iRow < numberRows; iRow++) {
5485    int iStatus = statusArray[numberColumns + iRow] & 7;
5486    iStatus = lookupA[iStatus];
5487    basis->setArtifStatus(iRow, static_cast< CoinWarmStartBasis::Status >(iStatus));
5488  }
5489  int lookupS[]