source: stable/1.17/Clp/src/ClpSimplexNonlinear.cpp @ 2438

Last change on this file since 2438 was 2438, checked in by forrest, 4 months ago

try and fix event handler

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 152.2 KB
RevLine 
[1370]1/* $Id: ClpSimplexNonlinear.cpp 2438 2019-03-26 12:39:54Z forrest $ */
[393]2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1665]4// This code is licensed under the terms of the Eclipse Public License (EPL).
[393]5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9#include "CoinHelperFunctions.hpp"
[401]10#include "ClpHelperFunctions.hpp"
[393]11#include "ClpSimplexNonlinear.hpp"
[2025]12#include "ClpSimplexOther.hpp"
13#include "ClpSimplexDual.hpp"
[393]14#include "ClpFactorization.hpp"
15#include "ClpNonLinearCost.hpp"
16#include "ClpLinearObjective.hpp"
[1034]17#include "ClpConstraint.hpp"
[2025]18#include "ClpPresolve.hpp"
[451]19#include "ClpQuadraticObjective.hpp"
[393]20#include "CoinPackedMatrix.hpp"
21#include "CoinIndexedVector.hpp"
22#include "ClpPrimalColumnPivot.hpp"
23#include "ClpMessage.hpp"
24#include "ClpEventHandler.hpp"
25#include <cfloat>
26#include <cassert>
27#include <string>
28#include <stdio.h>
29#include <iostream>
30#ifndef NDEBUG
31#define CLP_DEBUG 1
32#endif
[1525]33// primal
[2385]34int ClpSimplexNonlinear::primal()
[393]35{
36
[2385]37  int ifValuesPass = 1;
38  algorithm_ = +3;
[393]39
[2385]40  // save data
41  ClpDataSave data = saveData();
42  matrix_->refresh(this); // make sure matrix okay
[393]43
[2385]44  // Save objective
45  ClpObjective *saveObjective = NULL;
46  if (objective_->type() > 1) {
47    // expand to full if quadratic
[451]48#ifndef NO_RTTI
[2385]49    ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
[451]50#else
[2385]51    ClpQuadraticObjective *quadraticObj = NULL;
52    if (objective_->type() == 2)
53      quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
[451]54#endif
[2385]55    // for moment only if no scaling
56    // May be faster if switched off - but can't see why
57    if (!quadraticObj->fullMatrix() && (!rowScale_ && !scalingFlag_) && objectiveScale_ == 1.0) {
58      saveObjective = objective_;
59      objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
60    }
61  }
62  double bestObjectiveWhenFlagged = COIN_DBL_MAX;
63  int pivotMode = 15;
64  //pivotMode=20;
[393]65
[2385]66  // initialize - maybe values pass and algorithm_ is +1
67  // true does something (? perturbs)
68  if (!startup(true)) {
[1525]69
[2385]70    // Set average theta
71    nonLinearCost_->setAverageTheta(1.0e3);
72    int lastCleaned = 0; // last time objective or bounds cleaned up
[1525]73
[2385]74    // Say no pivot has occurred (for steepest edge and updates)
75    pivotRow_ = -2;
[1525]76
[2385]77    // This says whether to restore things etc
78    int factorType = 0;
79    // Start check for cycles
80    progress_.startCheck();
81    /*
[1525]82            Status of problem:
83            0 - optimal
84            1 - infeasible
85            2 - unbounded
86            -1 - iterating
87            -2 - factorization wanted
88            -3 - redo checking without factorization
89            -4 - looks infeasible
90            -5 - looks unbounded
91          */
[2385]92    while (problemStatus_ < 0) {
93      int iRow, iColumn;
94      // clear
95      for (iRow = 0; iRow < 4; iRow++) {
96        rowArray_[iRow]->clear();
97      }
[1525]98
[2385]99      for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
100        columnArray_[iColumn]->clear();
101      }
[1525]102
[2385]103      // give matrix (and model costs and bounds a chance to be
104      // refreshed (normally null)
105      matrix_->refresh(this);
106      // If getting nowhere - why not give it a kick
107      // If we have done no iterations - special
108      if (lastGoodIteration_ == numberIterations_ && factorType)
109        factorType = 3;
[1525]110
[2385]111      // may factorize, checks if problem finished
112      if (objective_->type() > 1 && lastFlaggedIteration_ >= 0 && numberIterations_ > lastFlaggedIteration_ + 507) {
113        unflag();
114        lastFlaggedIteration_ = numberIterations_;
115        if (pivotMode >= 10) {
116          pivotMode--;
[393]117#ifdef CLP_DEBUG
[2385]118          if (handler_->logLevel() & 32)
119            printf("pivot mode now %d\n", pivotMode);
[393]120#endif
[2385]121          if (pivotMode == 9)
122            pivotMode = 0; // switch off fast attempt
123        }
124      }
125      statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true,
126        bestObjectiveWhenFlagged);
[393]127
[2385]128      // Say good factorization
129      factorType = 1;
[393]130
[2385]131      // Say no pivot has occurred (for steepest edge and updates)
132      pivotRow_ = -2;
[1525]133
[2385]134      // exit if victory declared
135      if (problemStatus_ >= 0)
136        break;
[1525]137
[2385]138      // test for maximum iterations
139      if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
140        problemStatus_ = 3;
141        break;
142      }
[1525]143
[2385]144      if (firstFree_ < 0) {
145        if (ifValuesPass) {
146          // end of values pass
147          ifValuesPass = 0;
148          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
149          if (status >= 0) {
150            problemStatus_ = 5;
151            secondaryStatus_ = ClpEventHandler::endOfValuesPass;
152            break;
[1525]153          }
[2385]154        }
155      }
156      // Check event
157      {
158        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
159        if (status >= 0) {
160          problemStatus_ = 5;
161          secondaryStatus_ = ClpEventHandler::endOfFactorization;
162          break;
163        }
164      }
165      // Iterate
166      whileIterating(pivotMode);
167    }
168  }
169  // if infeasible get real values
170  if (problemStatus_ == 1) {
171    infeasibilityCost_ = 0.0;
172    createRim(1 + 4);
173    delete nonLinearCost_;
174    nonLinearCost_ = new ClpNonLinearCost(this);
175    nonLinearCost_->checkInfeasibilities(0.0);
176    sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
177    numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
178    // and get good feasible duals
179    computeDuals(NULL);
180  }
181  // correct objective value
182  if (numberColumns_)
183    objectiveValue_ = nonLinearCost_->feasibleCost() + objective_->nonlinearOffset();
184  objectiveValue_ /= (objectiveScale_ * rhsScale_);
185  // clean up
186  unflag();
187  finish();
188  restoreData(data);
189  // restore objective if full
190  if (saveObjective) {
191    delete objective_;
192    objective_ = saveObjective;
193  }
194  return problemStatus_;
[393]195}
[1525]196/*  Refactorizes if necessary
[393]197    Checks if finished.  Updates status.
198    lastCleaned refers to iteration at which some objective/feasibility
199    cleaning too place.
[1525]200
[393]201    type - 0 initial so set up save arrays etc
202    - 1 normal -if good update save
[1525]203    - 2 restoring from saved
[393]204*/
[2385]205void ClpSimplexNonlinear::statusOfProblemInPrimal(int &lastCleaned, int type,
206  ClpSimplexProgress *progress,
207  bool doFactorization,
208  double &bestObjectiveWhenFlagged)
[393]209{
[2385]210  int dummy; // for use in generalExpanded
211  if (type == 2) {
212    // trouble - restore solution
213    CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_);
214    CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_);
215    CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_);
216    // restore extra stuff
217    matrix_->generalExpanded(this, 6, dummy);
218    forceFactorization_ = 1; // a bit drastic but ..
219    pivotRow_ = -1; // say no weights update
220    changeMade_++; // say change made
221  }
222  int saveThreshold = factorization_->sparseThreshold();
223  int tentativeStatus = problemStatus_;
224  int numberThrownOut = 1; // to loop round on bad factorization in values pass
225  while (numberThrownOut) {
226    if (problemStatus_ > -3 || problemStatus_ == -4) {
227      // factorize
228      // later on we will need to recover from singularities
229      // also we could skip if first time
230      // do weights
231      // This may save pivotRow_ for use
232      if (doFactorization)
233        primalColumnPivot_->saveWeights(this, 1);
[393]234
[2385]235      if (type && doFactorization) {
236        // is factorization okay?
237        int factorStatus = internalFactorize(1);
238        if (factorStatus) {
239          if (type != 1 || largestPrimalError_ > 1.0e3
240            || largestDualError_ > 1.0e3) {
241            // was ||largestDualError_>1.0e3||objective_->type()>1) {
242            // switch off dense
243            int saveDense = factorization_->denseThreshold();
244            factorization_->setDenseThreshold(0);
245            // make sure will do safe factorization
246            pivotVariable_[0] = -1;
247            internalFactorize(2);
248            factorization_->setDenseThreshold(saveDense);
249            // Go to safe
250            factorization_->pivotTolerance(0.99);
251            // restore extra stuff
252            matrix_->generalExpanded(this, 6, dummy);
[1525]253          } else {
[2385]254            // no - restore previous basis
255            CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_);
256            CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_);
257            CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_);
258            // restore extra stuff
259            matrix_->generalExpanded(this, 6, dummy);
260            matrix_->generalExpanded(this, 5, dummy);
261            forceFactorization_ = 1; // a bit drastic but ..
262            type = 2;
263            // Go to safe
264            factorization_->pivotTolerance(0.99);
265            if (internalFactorize(1) != 0)
266              largestPrimalError_ = 1.0e4; // force other type
[1525]267          }
[2385]268          // flag incoming
269          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
270            setFlagged(sequenceIn_);
271            saveStatus_[sequenceIn_] = status_[sequenceIn_];
[1525]272          }
[2385]273          changeMade_++; // say change made
274        }
275      }
276      if (problemStatus_ != -4)
277        problemStatus_ = -3;
278    }
279    // at this stage status is -3 or -5 if looks unbounded
280    // get primal and dual solutions
281    // put back original costs and then check
282    createRim(4);
283    // May need to do more if column generation
284    dummy = 4;
285    matrix_->generalExpanded(this, 9, dummy);
286    numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
287    if (numberThrownOut) {
288      problemStatus_ = tentativeStatus;
289      doFactorization = true;
290    }
291  }
292  // Double check reduced costs if no action
293  if (progress->lastIterationNumber(0) == numberIterations_) {
294    if (primalColumnPivot_->looksOptimal()) {
295      numberDualInfeasibilities_ = 0;
296      sumDualInfeasibilities_ = 0.0;
297    }
298  }
299  // Check if looping
300  int loop;
301  if (type != 2)
302    loop = progress->looping();
303  else
304    loop = -1;
305  if (loop >= 0) {
306    if (!problemStatus_) {
307      // declaring victory
308      numberPrimalInfeasibilities_ = 0;
309      sumPrimalInfeasibilities_ = 0.0;
310    } else {
311      problemStatus_ = loop; //exit if in loop
312      problemStatus_ = 10; // instead - try other algorithm
313    }
314    problemStatus_ = 10; // instead - try other algorithm
315    return;
316  } else if (loop < -1) {
317    // Is it time for drastic measures
318    if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 && progress->oddState() < 10 && progress->oddState() >= 0) {
319      progress->newOddState();
320      nonLinearCost_->zapCosts();
321    }
322    // something may have changed
323    gutsOfSolution(NULL, NULL, true);
324  }
325  // If progress then reset costs
326  if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
327    createRim(4, false); // costs back
328    delete nonLinearCost_;
329    nonLinearCost_ = new ClpNonLinearCost(this);
330    progress->endOddState();
331    gutsOfSolution(NULL, NULL, true);
332  }
333  // Flag to say whether to go to dual to clean up
334  //bool goToDual = false;
335  // really for free variables in
336  //if((progressFlag_&2)!=0)
337  //problemStatus_=-1;;
338  progressFlag_ = 0; //reset progress flag
[393]339
[2385]340  handler_->message(CLP_SIMPLEX_STATUS, messages_)
341    << numberIterations_ << nonLinearCost_->feasibleReportCost();
342  handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
343    << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
344  handler_->printing(sumDualInfeasibilities_ > 0.0)
345    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
346  handler_->printing(numberDualInfeasibilitiesWithoutFree_
347    < numberDualInfeasibilities_)
348    << numberDualInfeasibilitiesWithoutFree_;
349  handler_->message() << CoinMessageEol;
350  if (!primalFeasible()) {
351    nonLinearCost_->checkInfeasibilities(primalTolerance_);
352    gutsOfSolution(NULL, NULL, true);
353    nonLinearCost_->checkInfeasibilities(primalTolerance_);
354  }
355  double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
356  if (trueInfeasibility > 1.0) {
357    // If infeasibility going up may change weights
358    double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
359    if (progress->lastInfeasibility() < testValue) {
360      if (infeasibilityCost_ < 1.0e14) {
361        infeasibilityCost_ *= 1.5;
362        if (handler_->logLevel() == 63)
363          printf("increasing weight to %g\n", infeasibilityCost_);
364        gutsOfSolution(NULL, NULL, true);
365      }
366    }
367  }
368  // we may wish to say it is optimal even if infeasible
369  bool alwaysOptimal = (specialOptions_ & 1) != 0;
370  // give code benefit of doubt
371  if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
372    // say optimal (with these bounds etc)
373    numberDualInfeasibilities_ = 0;
374    sumDualInfeasibilities_ = 0.0;
375    numberPrimalInfeasibilities_ = 0;
376    sumPrimalInfeasibilities_ = 0.0;
377  }
378  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
379  if (dualFeasible() || problemStatus_ == -4) {
380    // see if extra helps
381    if (nonLinearCost_->numberInfeasibilities() && (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
382      && !alwaysOptimal) {
383      //may need infeasiblity cost changed
384      // we can see if we can construct a ray
385      // make up a new objective
386      double saveWeight = infeasibilityCost_;
387      // save nonlinear cost as we are going to switch off costs
388      ClpNonLinearCost *nonLinear = nonLinearCost_;
389      // do twice to make sure Primal solution has settled
390      // put non-basics to bounds in case tolerance moved
391      // put back original costs
392      createRim(4);
393      nonLinearCost_->checkInfeasibilities(0.0);
394      gutsOfSolution(NULL, NULL, true);
[1525]395
[2385]396      infeasibilityCost_ = 1.0e100;
397      // put back original costs
398      createRim(4);
399      nonLinearCost_->checkInfeasibilities(primalTolerance_);
400      // may have fixed infeasibilities - double check
401      if (nonLinearCost_->numberInfeasibilities() == 0) {
402        // carry on
403        problemStatus_ = -1;
404        infeasibilityCost_ = saveWeight;
405        nonLinearCost_->checkInfeasibilities(primalTolerance_);
406      } else {
407        nonLinearCost_ = NULL;
408        // scale
409        int i;
410        for (i = 0; i < numberRows_ + numberColumns_; i++)
411          cost_[i] *= 1.0e-95;
412        gutsOfSolution(NULL, NULL, false);
413        nonLinearCost_ = nonLinear;
414        infeasibilityCost_ = saveWeight;
415        if ((infeasibilityCost_ >= 1.0e18 || numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
416          /* goToDual = unPerturb(); // stop any further perturbation
[1525]417                         if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
418                              goToDual = false;
[1929]419                         */
[2385]420          unPerturb();
421          nonLinearCost_->checkInfeasibilities(primalTolerance_);
422          numberDualInfeasibilities_ = 1; // carry on
423          problemStatus_ = -1;
424        }
425        if (infeasibilityCost_ >= 1.0e20 || numberDualInfeasibilities_ == 0) {
426          // we are infeasible - use as ray
427          delete[] ray_;
428          ray_ = new double[numberRows_];
429          CoinMemcpyN(dual_, numberRows_, ray_);
430          // and get feasible duals
431          infeasibilityCost_ = 0.0;
432          createRim(4);
433          nonLinearCost_->checkInfeasibilities(primalTolerance_);
434          gutsOfSolution(NULL, NULL, true);
435          // so will exit
436          infeasibilityCost_ = 1.0e30;
437          // reset infeasibilities
438          sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
439          ;
440          numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
441        }
442        if (infeasibilityCost_ < 1.0e20) {
443          infeasibilityCost_ *= 5.0;
444          changeMade_++; // say change made
445          unflag();
446          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
447            << infeasibilityCost_
448            << CoinMessageEol;
449          // put back original costs and then check
450          createRim(4);
451          nonLinearCost_->checkInfeasibilities(0.0);
452          gutsOfSolution(NULL, NULL, true);
453          problemStatus_ = -1; //continue
454          //goToDual = false;
455        } else {
456          // say infeasible
457          problemStatus_ = 1;
458        }
459      }
460    } else {
461      // may be optimal
462      if (perturbation_ == 101) {
463        /* goToDual =*/unPerturb(); // stop any further perturbation
464        lastCleaned = -1; // carry on
465      }
466      bool unflagged = unflag() != 0;
467      if (lastCleaned != numberIterations_ || unflagged) {
468        handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
469          << primalTolerance_
470          << CoinMessageEol;
471        double current = nonLinearCost_->feasibleReportCost();
472        if (numberTimesOptimal_ < 4) {
473          if (bestObjectiveWhenFlagged <= current) {
474            numberTimesOptimal_++;
[393]475#ifdef CLP_DEBUG
[2385]476            if (handler_->logLevel() & 32)
477              printf("%d times optimal, current %g best %g\n", numberTimesOptimal_,
478                current, bestObjectiveWhenFlagged);
[393]479#endif
[1525]480          } else {
[2385]481            bestObjectiveWhenFlagged = current;
[1525]482          }
[2385]483          changeMade_++; // say change made
484          if (numberTimesOptimal_ == 1) {
485            // better to have small tolerance even if slower
486            factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
[1525]487          }
[2385]488          lastCleaned = numberIterations_;
489          if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
490            handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
491              << CoinMessageEol;
492          double oldTolerance = primalTolerance_;
493          primalTolerance_ = dblParam_[ClpPrimalTolerance];
494          // put back original costs and then check
495          createRim(4);
496          nonLinearCost_->checkInfeasibilities(oldTolerance);
497          gutsOfSolution(NULL, NULL, true);
498          if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
499            // say optimal (with these bounds etc)
500            numberDualInfeasibilities_ = 0;
501            sumDualInfeasibilities_ = 0.0;
502            numberPrimalInfeasibilities_ = 0;
503            sumPrimalInfeasibilities_ = 0.0;
504          }
505          if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
506            problemStatus_ = 0;
[1525]507          else
[2385]508            problemStatus_ = -1;
509        } else {
510          problemStatus_ = 0; // optimal
511          if (lastCleaned < numberIterations_) {
512            handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
513              << CoinMessageEol;
[1525]514          }
[2385]515        }
516      } else {
517        problemStatus_ = 0; // optimal
518      }
519    }
520  } else {
521    // see if looks unbounded
522    if (problemStatus_ == -5) {
523      if (nonLinearCost_->numberInfeasibilities()) {
524        if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) {
525          // back off weight
526          infeasibilityCost_ = 1.0e13;
527          unPerturb(); // stop any further perturbation
528        }
529        //we need infeasiblity cost changed
530        if (infeasibilityCost_ < 1.0e20) {
531          infeasibilityCost_ *= 5.0;
532          changeMade_++; // say change made
533          unflag();
534          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
535            << infeasibilityCost_
536            << CoinMessageEol;
537          // put back original costs and then check
538          createRim(4);
539          gutsOfSolution(NULL, NULL, true);
540          problemStatus_ = -1; //continue
541        } else {
542          // say unbounded
543          problemStatus_ = 2;
544        }
545      } else {
546        // say unbounded
547        problemStatus_ = 2;
548      }
549    } else {
550      if (type == 3 && problemStatus_ != -5)
551        unflag(); // odd
552      // carry on
553      problemStatus_ = -1;
554    }
555  }
556  // save extra stuff
557  matrix_->generalExpanded(this, 5, dummy);
558  if (type == 0 || type == 1) {
559    if (type != 1 || !saveStatus_) {
560      // create save arrays
561      delete[] saveStatus_;
562      delete[] savedSolution_;
563      saveStatus_ = new unsigned char[numberRows_ + numberColumns_];
564      savedSolution_ = new double[numberRows_ + numberColumns_];
565    }
566    // save arrays
567    CoinMemcpyN(status_, (numberColumns_ + numberRows_), saveStatus_);
568    CoinMemcpyN(rowActivityWork_, numberRows_, savedSolution_ + numberColumns_);
569    CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
570  }
571  if (doFactorization) {
572    // restore weights (if saved) - also recompute infeasibility list
573    if (tentativeStatus > -3)
574      primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
575    else
576      primalColumnPivot_->saveWeights(this, 3);
577    if (saveThreshold) {
578      // use default at present
579      factorization_->sparseThreshold(0);
580      factorization_->goSparse();
581    }
582  }
583  if (problemStatus_ < 0 && !changeMade_) {
584    problemStatus_ = 4; // unknown
585  }
586  lastGoodIteration_ = numberIterations_;
587  //if (goToDual)
588  //problemStatus_=10; // try dual
589  // Allow matrices to be sorted etc
590  int fake = -999; // signal sort
591  matrix_->correctSequence(this, fake, fake);
[393]592}
593/*
594  Reasons to come out:
595  -1 iterations etc
[1525]596  -2 inaccuracy
[393]597  -3 slight inaccuracy (and done iterations)
598  -4 end of values pass and done iterations
599  +0 looks optimal (might be infeasible - but we will investigate)
600  +2 looks unbounded
[1525]601  +3 max iterations
[393]602*/
[2385]603int ClpSimplexNonlinear::whileIterating(int &pivotMode)
[393]604{
[2385]605  // Say if values pass
606  //int ifValuesPass=(firstFree_>=0) ? 1 : 0;
607  int ifValuesPass = 1;
608  int returnCode = -1;
609  // status stays at -1 while iterating, >=0 finished, -2 to invert
610  // status -3 to go to top without an invert
611  int numberInterior = 0;
612  int nextUnflag = 10;
613  int nextUnflagIteration = numberIterations_ + 10;
614  // get two arrays
615  double *array1 = new double[2 * (numberRows_ + numberColumns_)];
616  double solutionError = -1.0;
617  while (problemStatus_ == -1) {
618    int result;
619    rowArray_[1]->clear();
620    //#define CLP_DEBUG
[393]621#if CLP_DEBUG > 1
[2385]622    rowArray_[0]->checkClear();
623    rowArray_[1]->checkClear();
624    rowArray_[2]->checkClear();
625    rowArray_[3]->checkClear();
626    columnArray_[0]->checkClear();
[393]627#endif
[2385]628    if (numberInterior >= 5) {
629      // this can go when we have better minimization
630      if (pivotMode < 10)
631        pivotMode = 1;
632      unflag();
[393]633#ifdef CLP_DEBUG
[2385]634      if (handler_->logLevel() & 32)
635        printf("interior unflag\n");
[393]636#endif
[2385]637      numberInterior = 0;
638      nextUnflag = 10;
639      nextUnflagIteration = numberIterations_ + 10;
640    } else {
641      if (numberInterior > nextUnflag && numberIterations_ > nextUnflagIteration) {
642        nextUnflagIteration = numberIterations_ + 10;
643        nextUnflag += 10;
644        unflag();
[393]645#ifdef CLP_DEBUG
[2385]646        if (handler_->logLevel() & 32)
647          printf("unflagging as interior\n");
[393]648#endif
[2385]649      }
650    }
651    pivotRow_ = -1;
652    result = pivotColumn(rowArray_[3], rowArray_[0],
653      columnArray_[0], rowArray_[1], pivotMode, solutionError,
654      array1);
655    if (result) {
656      if (result == 2 && sequenceIn_ < 0) {
657        // does not look good
658        double currentObj;
659        double thetaObj;
660        double predictedObj;
661        objective_->stepLength(this, solution_, solution_, 0.0,
662          currentObj, thetaObj, predictedObj);
663        if (currentObj == predictedObj) {
[1366]664#ifdef CLP_INVESTIGATE
[2385]665          printf("looks bad - no change in obj %g\n", currentObj);
[1366]666#endif
[2385]667          if (factorization_->pivots())
668            result = 3;
669          else
670            problemStatus_ = 0;
671        }
672      }
673      if (result == 3)
674        break; // null vector not accurate
[393]675#ifdef CLP_DEBUG
[2385]676      if (handler_->logLevel() & 32) {
677        double currentObj;
678        double thetaObj;
679        double predictedObj;
680        objective_->stepLength(this, solution_, solution_, 0.0,
681          currentObj, thetaObj, predictedObj);
682        printf("obj %g after interior move\n", currentObj);
683      }
[393]684#endif
[2385]685      // just move and try again
686      if (pivotMode < 10) {
687        pivotMode = result - 1;
688        numberInterior++;
689      }
690      continue;
691    } else {
692      if (pivotMode < 10) {
693        if (theta_ > 0.001)
694          pivotMode = 0;
695        else if (pivotMode == 2)
696          pivotMode = 1;
697      }
698      numberInterior = 0;
699      nextUnflag = 10;
700      nextUnflagIteration = numberIterations_ + 10;
701    }
702    sequenceOut_ = -1;
703    rowArray_[1]->clear();
704    if (sequenceIn_ >= 0) {
705      // we found a pivot column
706      assert(!flagged(sequenceIn_));
[393]707#ifdef CLP_DEBUG
[2385]708      if ((handler_->logLevel() & 32)) {
709        char x = isColumn(sequenceIn_) ? 'C' : 'R';
710        std::cout << "pivot column " << x << sequenceWithin(sequenceIn_) << std::endl;
711      }
[393]712#endif
[2385]713      // do second half of iteration
714      if (pivotRow_ < 0 && theta_ < 1.0e-8) {
715        assert(sequenceIn_ >= 0);
716        returnCode = pivotResult(ifValuesPass);
717      } else {
718        // specialized code
719        returnCode = pivotNonlinearResult();
720        //printf("odd pivrow %d\n",sequenceOut_);
721        if (sequenceOut_ >= 0 && theta_ < 1.0e-5) {
722          if (getStatus(sequenceOut_) != isFixed) {
723            if (getStatus(sequenceOut_) == atUpperBound)
724              solution_[sequenceOut_] = upper_[sequenceOut_];
725            else if (getStatus(sequenceOut_) == atLowerBound)
726              solution_[sequenceOut_] = lower_[sequenceOut_];
727            setFlagged(sequenceOut_);
728          }
729        }
730      }
731      if (returnCode < -1 && returnCode > -5) {
732        problemStatus_ = -2; //
733      } else if (returnCode == -5) {
734        // something flagged - continue;
735      } else if (returnCode == 2) {
736        problemStatus_ = -5; // looks unbounded
737      } else if (returnCode == 4) {
738        problemStatus_ = -2; // looks unbounded but has iterated
739      } else if (returnCode != -1) {
740        assert(returnCode == 3);
741        problemStatus_ = 3;
742      }
743    } else {
744      // no pivot column
[393]745#ifdef CLP_DEBUG
[2385]746      if (handler_->logLevel() & 32)
747        printf("** no column pivot\n");
[393]748#endif
[2385]749      if (pivotMode < 10) {
750        // looks optimal
751        primalColumnPivot_->setLooksOptimal(true);
752      } else {
753        pivotMode--;
[393]754#ifdef CLP_DEBUG
[2385]755        if (handler_->logLevel() & 32)
756          printf("pivot mode now %d\n", pivotMode);
[393]757#endif
[2385]758        if (pivotMode == 9)
759          pivotMode = 0; // switch off fast attempt
760        unflag();
761      }
762      if (nonLinearCost_->numberInfeasibilities())
763        problemStatus_ = -4; // might be infeasible
764      returnCode = 0;
765      break;
766    }
767  }
768  delete[] array1;
769  return returnCode;
[393]770}
771// Creates direction vector
[2385]772void ClpSimplexNonlinear::directionVector(CoinIndexedVector *vectorArray,
773  CoinIndexedVector *spare1, CoinIndexedVector *spare2,
774  int pivotMode2,
775  double &normFlagged, double &normUnflagged,
776  int &numberNonBasic)
[393]777{
778#if CLP_DEBUG > 1
[2385]779  vectorArray->checkClear();
780  spare1->checkClear();
781  spare2->checkClear();
[393]782#endif
[2385]783  double *array = vectorArray->denseVector();
784  int *index = vectorArray->getIndices();
785  int number = 0;
786  sequenceIn_ = -1;
787  normFlagged = 0.0;
788  normUnflagged = 1.0;
789  double dualTolerance2 = CoinMin(1.0e-8, 1.0e-2 * dualTolerance_);
790  double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
791  if (!numberNonBasic) {
792    //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
793    //printf("infeasible\n");
794    if (!pivotMode2 || pivotMode2 >= 10) {
795      normUnflagged = 0.0;
796      double bestDj = 0.0;
797      double bestSuper = 0.0;
798      double sumSuper = 0.0;
799      sequenceIn_ = -1;
800      int nSuper = 0;
801      for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
802        array[iSequence] = 0.0;
803        if (flagged(iSequence)) {
804          // accumulate  norm
805          switch (getStatus(iSequence)) {
[1525]806
[2385]807          case basic:
808          case ClpSimplex::isFixed:
809            break;
810          case atUpperBound:
811            if (dj_[iSequence] > dualTolerance3)
812              normFlagged += dj_[iSequence] * dj_[iSequence];
813            break;
814          case atLowerBound:
815            if (dj_[iSequence] < -dualTolerance3)
816              normFlagged += dj_[iSequence] * dj_[iSequence];
817            break;
818          case isFree:
819          case superBasic:
820            if (fabs(dj_[iSequence]) > dualTolerance3)
821              normFlagged += dj_[iSequence] * dj_[iSequence];
822            break;
823          }
824          continue;
825        }
826        switch (getStatus(iSequence)) {
[1525]827
[2385]828        case basic:
829        case ClpSimplex::isFixed:
830          break;
831        case atUpperBound:
832          if (dj_[iSequence] > dualTolerance_) {
833            if (dj_[iSequence] > dualTolerance3)
834              normUnflagged += dj_[iSequence] * dj_[iSequence];
835            if (pivotMode2 < 10) {
836              array[iSequence] = -dj_[iSequence];
837              index[number++] = iSequence;
838            } else {
839              if (dj_[iSequence] > bestDj) {
840                bestDj = dj_[iSequence];
841                sequenceIn_ = iSequence;
842              }
843            }
844          }
845          break;
846        case atLowerBound:
847          if (dj_[iSequence] < -dualTolerance_) {
848            if (dj_[iSequence] < -dualTolerance3)
849              normUnflagged += dj_[iSequence] * dj_[iSequence];
850            if (pivotMode2 < 10) {
851              array[iSequence] = -dj_[iSequence];
852              index[number++] = iSequence;
853            } else {
854              if (-dj_[iSequence] > bestDj) {
855                bestDj = -dj_[iSequence];
856                sequenceIn_ = iSequence;
857              }
858            }
859          }
860          break;
861        case isFree:
862        case superBasic:
863          //#define ALLSUPER
[393]864#define NOSUPER
865#ifndef ALLSUPER
[2385]866          if (fabs(dj_[iSequence]) > dualTolerance_) {
867            if (fabs(dj_[iSequence]) > dualTolerance3)
868              normUnflagged += dj_[iSequence] * dj_[iSequence];
869            nSuper++;
870            bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
871            sumSuper += fabs(dj_[iSequence]);
872          }
873          if (fabs(dj_[iSequence]) > dualTolerance2) {
874            array[iSequence] = -dj_[iSequence];
875            index[number++] = iSequence;
876          }
[393]877#else
[2385]878          array[iSequence] = -dj_[iSequence];
879          index[number++] = iSequence;
880          if (pivotMode2 >= 10)
881            bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
[393]882#endif
[2385]883          break;
884        }
885      }
[393]886#ifdef NOSUPER
[2385]887      // redo
888      bestSuper = sumSuper;
889      if (sequenceIn_ >= 0 && bestDj > bestSuper) {
890        int j;
891        // get rid of superbasics
892        for (j = 0; j < number; j++) {
893          int iSequence = index[j];
894          array[iSequence] = 0.0;
895        }
896        number = 0;
897        array[sequenceIn_] = -dj_[sequenceIn_];
898        index[number++] = sequenceIn_;
899      } else {
900        sequenceIn_ = -1;
901      }
[393]902#else
[2385]903      if (pivotMode2 >= 10 || !nSuper) {
904        bool takeBest = true;
905        if (pivotMode2 == 100 && nSuper > 1)
906          takeBest = false;
907        if (sequenceIn_ >= 0 && takeBest) {
908          if (fabs(dj_[sequenceIn_]) > bestSuper) {
909            array[sequenceIn_] = -dj_[sequenceIn_];
910            index[number++] = sequenceIn_;
911          } else {
912            sequenceIn_ = -1;
913          }
914        } else {
915          sequenceIn_ = -1;
916        }
917      }
[393]918#endif
919#ifdef CLP_DEBUG
[2385]920      if (handler_->logLevel() & 32) {
921        if (sequenceIn_ >= 0)
922          printf("%d superBasic, chosen %d - dj %g\n", nSuper, sequenceIn_,
923            dj_[sequenceIn_]);
924        else
925          printf("%d superBasic - none chosen\n", nSuper);
926      }
[393]927#endif
[2385]928    } else {
929      double bestDj = 0.0;
930      double saveDj = 0.0;
931      if (sequenceOut_ >= 0) {
932        saveDj = dj_[sequenceOut_];
933        dj_[sequenceOut_] = 0.0;
934        switch (getStatus(sequenceOut_)) {
[1525]935
[2385]936        case basic:
937          sequenceOut_ = -1;
938        case ClpSimplex::isFixed:
939          break;
940        case atUpperBound:
941          if (dj_[sequenceOut_] > dualTolerance_) {
[393]942#ifdef CLP_DEBUG
[2385]943            if (handler_->logLevel() & 32)
944              printf("after pivot out %d values %g %g %g, dj %g\n",
945                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
946                upper_[sequenceOut_], dj_[sequenceOut_]);
[393]947#endif
[2385]948          }
949          break;
950        case atLowerBound:
951          if (dj_[sequenceOut_] < -dualTolerance_) {
[393]952#ifdef CLP_DEBUG
[2385]953            if (handler_->logLevel() & 32)
954              printf("after pivot out %d values %g %g %g, dj %g\n",
955                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
956                upper_[sequenceOut_], dj_[sequenceOut_]);
[393]957#endif
[2385]958          }
959          break;
960        case isFree:
961        case superBasic:
962          if (dj_[sequenceOut_] > dualTolerance_) {
[393]963#ifdef CLP_DEBUG
[2385]964            if (handler_->logLevel() & 32)
965              printf("after pivot out %d values %g %g %g, dj %g\n",
966                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
967                upper_[sequenceOut_], dj_[sequenceOut_]);
[393]968#endif
[2385]969          } else if (dj_[sequenceOut_] < -dualTolerance_) {
[393]970#ifdef CLP_DEBUG
[2385]971            if (handler_->logLevel() & 32)
972              printf("after pivot out %d values %g %g %g, dj %g\n",
973                sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
974                upper_[sequenceOut_], dj_[sequenceOut_]);
[393]975#endif
[1525]976          }
[2385]977          break;
978        }
979      }
980      // Go for dj
981      pivotMode2 = 3;
982      for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
983        array[iSequence] = 0.0;
984        if (flagged(iSequence))
985          continue;
986        switch (getStatus(iSequence)) {
[1525]987
[2385]988        case basic:
989        case ClpSimplex::isFixed:
990          break;
991        case atUpperBound:
992          if (dj_[iSequence] > dualTolerance_) {
993            double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]);
994            double merit = distance * dj_[iSequence];
995            if (pivotMode2 == 1)
996              merit *= 1.0e-20; // discourage
997            if (pivotMode2 == 3)
998              merit = fabs(dj_[iSequence]);
999            if (merit > bestDj) {
1000              sequenceIn_ = iSequence;
1001              bestDj = merit;
1002            }
[1525]1003          }
[2385]1004          break;
1005        case atLowerBound:
1006          if (dj_[iSequence] < -dualTolerance_) {
1007            double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]);
1008            double merit = -distance * dj_[iSequence];
1009            if (pivotMode2 == 1)
1010              merit *= 1.0e-20; // discourage
1011            if (pivotMode2 == 3)
1012              merit = fabs(dj_[iSequence]);
1013            if (merit > bestDj) {
1014              sequenceIn_ = iSequence;
1015              bestDj = merit;
1016            }
[1525]1017          }
[2385]1018          break;
1019        case isFree:
1020        case superBasic:
1021          if (dj_[iSequence] > dualTolerance_) {
1022            double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]);
1023            distance = CoinMin(solution_[iSequence] - lower_[iSequence],
1024              upper_[iSequence] - solution_[iSequence]);
1025            double merit = distance * dj_[iSequence];
1026            if (pivotMode2 == 1)
1027              merit = distance;
1028            if (pivotMode2 == 3)
1029              merit = fabs(dj_[iSequence]);
1030            if (merit > bestDj) {
1031              sequenceIn_ = iSequence;
1032              bestDj = merit;
1033            }
1034          } else if (dj_[iSequence] < -dualTolerance_) {
1035            double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]);
1036            distance = CoinMin(solution_[iSequence] - lower_[iSequence],
1037              upper_[iSequence] - solution_[iSequence]);
1038            double merit = -distance * dj_[iSequence];
1039            if (pivotMode2 == 1)
1040              merit = distance;
1041            if (pivotMode2 == 3)
1042              merit = fabs(dj_[iSequence]);
1043            if (merit > bestDj) {
1044              sequenceIn_ = iSequence;
1045              bestDj = merit;
1046            }
[1525]1047          }
[2385]1048          break;
1049        }
1050      }
1051      if (sequenceOut_ >= 0) {
1052        dj_[sequenceOut_] = saveDj;
1053        sequenceOut_ = -1;
1054      }
1055      if (sequenceIn_ >= 0) {
1056        array[sequenceIn_] = -dj_[sequenceIn_];
1057        index[number++] = sequenceIn_;
1058      }
1059    }
1060    numberNonBasic = number;
1061  } else {
1062    // compute norms
1063    normUnflagged = 0.0;
1064    for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
1065      if (flagged(iSequence)) {
1066        // accumulate  norm
1067        switch (getStatus(iSequence)) {
1068
1069        case basic:
1070        case ClpSimplex::isFixed:
1071          break;
1072        case atUpperBound:
1073          if (dj_[iSequence] > dualTolerance_)
1074            normFlagged += dj_[iSequence] * dj_[iSequence];
1075          break;
1076        case atLowerBound:
1077          if (dj_[iSequence] < -dualTolerance_)
1078            normFlagged += dj_[iSequence] * dj_[iSequence];
1079          break;
1080        case isFree:
1081        case superBasic:
1082          if (fabs(dj_[iSequence]) > dualTolerance_)
1083            normFlagged += dj_[iSequence] * dj_[iSequence];
1084          break;
1085        }
1086      }
1087    }
1088    // re-use list
1089    number = 0;
1090    int j;
1091    for (j = 0; j < numberNonBasic; j++) {
1092      int iSequence = index[j];
1093      if (flagged(iSequence))
1094        continue;
1095      switch (getStatus(iSequence)) {
1096
1097      case basic:
1098      case ClpSimplex::isFixed:
1099        continue; //abort();
1100        break;
1101      case atUpperBound:
1102        if (dj_[iSequence] > dualTolerance_) {
1103          number++;
1104          normUnflagged += dj_[iSequence] * dj_[iSequence];
1105        }
1106        break;
1107      case atLowerBound:
1108        if (dj_[iSequence] < -dualTolerance_) {
1109          number++;
1110          normUnflagged += dj_[iSequence] * dj_[iSequence];
1111        }
1112        break;
1113      case isFree:
1114      case superBasic:
1115        if (fabs(dj_[iSequence]) > dualTolerance_) {
1116          number++;
1117          normUnflagged += dj_[iSequence] * dj_[iSequence];
1118        }
1119        break;
1120      }
1121      array[iSequence] = -dj_[iSequence];
1122    }
1123    // switch to large
1124    normUnflagged = 1.0;
1125    if (!number) {
1126      for (j = 0; j < numberNonBasic; j++) {
1127        int iSequence = index[j];
1128        array[iSequence] = 0.0;
1129      }
1130      numberNonBasic = 0;
1131    }
1132    number = numberNonBasic;
1133  }
1134  if (number) {
1135    int j;
1136    // Now do basic ones - how do I compensate for small basic infeasibilities?
1137    int iRow;
1138    for (iRow = 0; iRow < numberRows_; iRow++) {
1139      int iPivot = pivotVariable_[iRow];
1140      double value = 0.0;
1141      if (solution_[iPivot] > upper_[iPivot]) {
1142        value = upper_[iPivot] - solution_[iPivot];
1143      } else if (solution_[iPivot] < lower_[iPivot]) {
1144        value = lower_[iPivot] - solution_[iPivot];
1145      }
1146      //if (value)
1147      //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
1148      //     upper_[iPivot]);
1149      //value=0.0;
1150      value *= -1.0e0;
1151      if (value) {
1152        array[iPivot] = value;
1153        index[number++] = iPivot;
1154      }
1155    }
1156    double *array2 = spare1->denseVector();
1157    int *index2 = spare1->getIndices();
1158    int number2 = 0;
1159    times(-1.0, array, array2);
1160    array = array + numberColumns_;
1161    for (iRow = 0; iRow < numberRows_; iRow++) {
1162      double value = array2[iRow] + array[iRow];
1163      if (value) {
1164        array2[iRow] = value;
1165        index2[number2++] = iRow;
1166      } else {
1167        array2[iRow] = 0.0;
1168      }
1169    }
1170    array -= numberColumns_;
1171    spare1->setNumElements(number2);
1172    // Ftran
1173    factorization_->updateColumn(spare2, spare1);
1174    number2 = spare1->getNumElements();
1175    for (j = 0; j < number2; j++) {
1176      int iSequence = index2[j];
1177      double value = array2[iSequence];
1178      array2[iSequence] = 0.0;
1179      if (value) {
1180        int iPivot = pivotVariable_[iSequence];
1181        double oldValue = array[iPivot];
1182        if (!oldValue) {
1183          array[iPivot] = value;
1184          index[number++] = iPivot;
1185        } else {
1186          // something already there
1187          array[iPivot] = value + oldValue;
1188        }
1189      }
1190    }
1191    spare1->setNumElements(0);
1192  }
1193  vectorArray->setNumElements(number);
[393]1194}
[451]1195#define MINTYPE 1
[2385]1196#if MINTYPE == 2
[1525]1197static double
[2385]1198innerProductIndexed(const double *region1, int size, const double *region2, const int *which)
[451]1199{
[2385]1200  int i;
1201  double value = 0.0;
1202  for (i = 0; i < size; i++) {
1203    int j = which[i];
1204    value += region1[j] * region2[j];
1205  }
1206  return value;
[451]1207}
1208#endif
[1525]1209/*
[393]1210   Row array and column array have direction
1211   Returns 0 - can do normal iteration (basis change)
1212   1 - no basis change
1213*/
[2385]1214int ClpSimplexNonlinear::pivotColumn(CoinIndexedVector *longArray,
1215  CoinIndexedVector *rowArray,
1216  CoinIndexedVector *columnArray,
1217  CoinIndexedVector *spare,
1218  int &pivotMode,
1219  double &solutionError,
1220  double *dArray)
[393]1221{
[2385]1222  // say not optimal
1223  primalColumnPivot_->setLooksOptimal(false);
1224  double acceptablePivot = 1.0e-10;
1225  int lastSequenceIn = -1;
1226  if (pivotMode && pivotMode < 10) {
1227    acceptablePivot = 1.0e-6;
1228    if (factorization_->pivots())
1229      acceptablePivot = 1.0e-5; // if we have iterated be more strict
1230  }
1231  double acceptableBasic = 1.0e-7;
[1525]1232
[2385]1233  int number = longArray->getNumElements();
1234  int numberTotal = numberRows_ + numberColumns_;
1235  int bestSequence = -1;
1236  int bestBasicSequence = -1;
1237  double eps = 1.0e-1;
1238  eps = 1.0e-6;
1239  double basicTheta = 1.0e30;
1240  double objTheta = 0.0;
1241  bool finished = false;
1242  sequenceIn_ = -1;
1243  int nPasses = 0;
1244  int nTotalPasses = 0;
1245  int nBigPasses = 0;
1246  double djNorm0 = 0.0;
1247  double djNorm = 0.0;
1248  double normFlagged = 0.0;
1249  double normUnflagged = 0.0;
1250  int localPivotMode = pivotMode;
1251  bool allFinished = false;
1252  bool justOne = false;
1253  int returnCode = 1;
1254  double currentObj;
1255  double predictedObj;
1256  double thetaObj;
1257  objective_->stepLength(this, solution_, solution_, 0.0,
1258    currentObj, predictedObj, thetaObj);
1259  double saveObj = currentObj;
1260#if MINTYPE == 2
1261  // try Shanno's method
1262  //would be memory leak
1263  //double * saveY=new double[numberTotal];
1264  //double * saveS=new double[numberTotal];
1265  //double * saveY2=new double[numberTotal];
1266  //double * saveS2=new double[numberTotal];
1267  double saveY[100];
1268  double saveS[100];
1269  double saveY2[100];
1270  double saveS2[100];
1271  double zz[10000];
[414]1272#endif
[2385]1273  double *dArray2 = dArray + numberTotal;
1274  // big big loop
1275  while (!allFinished) {
1276    double *work = longArray->denseVector();
1277    int *which = longArray->getIndices();
1278    allFinished = true;
1279    // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
1280    //#define SMALLTHETA1 1.0e-25
1281    //#define SMALLTHETA2 1.0e-25
[393]1282#define SMALLTHETA1 1.0e-10
1283#define SMALLTHETA2 1.0e-10
1284#define CONJUGATE 2
1285#if CONJUGATE == 0
[2385]1286    int conjugate = 0;
[393]1287#elif CONJUGATE == 1
[2385]1288    int conjugate = (pivotMode < 10) ? MINTYPE : 0;
[393]1289#elif CONJUGATE == 2
[2385]1290    int conjugate = MINTYPE;
[393]1291#else
[2385]1292    int conjugate = MINTYPE;
[393]1293#endif
[1525]1294
[2385]1295    //if (pivotMode==1)
1296    //localPivotMode=11;;
[393]1297#if CLP_DEBUG > 1
[2385]1298    longArray->checkClear();
[393]1299#endif
[2385]1300    int numberNonBasic = 0;
1301    int startLocalMode = -1;
1302    while (!finished) {
1303      double simpleObjective = COIN_DBL_MAX;
1304      returnCode = 1;
1305      int iSequence;
1306      objective_->reducedGradient(this, dj_, false);
1307      // get direction vector in longArray
1308      longArray->clear();
1309      // take out comment to try slightly different idea
1310      if (nPasses + nTotalPasses > 3000 || nBigPasses > 100) {
1311        if (factorization_->pivots())
1312          returnCode = 3;
1313        break;
1314      }
1315      if (!nPasses) {
1316        numberNonBasic = 0;
1317        nBigPasses++;
1318      }
1319      // just do superbasic if in cleanup mode
1320      int local = localPivotMode;
1321      if (!local && pivotMode >= 10 && nBigPasses < 10) {
1322        local = 100;
1323      } else if (local == -1 || nBigPasses >= 10) {
1324        local = 0;
1325        localPivotMode = 0;
1326      }
1327      if (justOne) {
1328        local = 2;
1329        //local=100;
1330        justOne = false;
1331      }
1332      if (!nPasses)
1333        startLocalMode = local;
1334      directionVector(longArray, spare, rowArray, local,
1335        normFlagged, normUnflagged, numberNonBasic);
1336      {
1337        // check null vector
1338        double *rhs = spare->denseVector();
[451]1339#if CLP_DEBUG > 1
[2385]1340        spare->checkClear();
[451]1341#endif
[2385]1342        int iRow;
1343        multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, rhs, 0.0);
1344        matrix_->times(1.0, solution_, rhs, rowScale_, columnScale_);
1345        double largest = 0.0;
[1929]1346#if CLP_DEBUG > 0
[2385]1347        int iLargest = -1;
[1929]1348#endif
[2385]1349        for (iRow = 0; iRow < numberRows_; iRow++) {
1350          double value = fabs(rhs[iRow]);
1351          rhs[iRow] = 0.0;
1352          if (value > largest) {
1353            largest = value;
[1929]1354#if CLP_DEBUG > 0
[2385]1355            iLargest = iRow;
[1929]1356#endif
[2385]1357          }
1358        }
[401]1359#if CLP_DEBUG > 0
[2385]1360        if ((handler_->logLevel() & 32) && largest > 1.0e-8)
1361          printf("largest rhs error %g on row %d\n", largest, iLargest);
[401]1362#endif
[2385]1363        if (solutionError < 0.0) {
1364          solutionError = largest;
1365        } else if (largest > CoinMax(1.0e-8, 1.0e2 * solutionError) && factorization_->pivots()) {
1366          longArray->clear();
1367          pivotRow_ = -1;
1368          theta_ = 0.0;
1369          return 3;
1370        }
1371      }
1372      if (sequenceIn_ >= 0)
1373        lastSequenceIn = sequenceIn_;
1374#if MINTYPE != 2
1375      double djNormSave = djNorm;
[451]1376#endif
[2385]1377      djNorm = 0.0;
1378      int iIndex;
1379      for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
1380        int iSequence = which[iIndex];
1381        double alpha = work[iSequence];
1382        djNorm += alpha * alpha;
1383      }
1384      // go to conjugate gradient if necessary
1385      if (numberNonBasic && localPivotMode >= 10 && (nPasses > 4 || sequenceIn_ < 0)) {
1386        localPivotMode = 0;
1387        nTotalPasses += nPasses;
1388        nPasses = 0;
1389      }
[393]1390#if CONJUGATE == 2
[2385]1391      conjugate = (localPivotMode < 10) ? MINTYPE : 0;
[393]1392#endif
[2385]1393      //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
1394      //     nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
1395      if (!nPasses) {
1396        //printf("numberNon %d\n",numberNonBasic);
1397#if MINTYPE == 2
1398        assert(numberNonBasic < 100);
1399        memset(zz, 0, numberNonBasic * numberNonBasic * sizeof(double));
1400        int put = 0;
1401        for (int iVariable = 0; iVariable < numberNonBasic; iVariable++) {
1402          zz[put] = 1.0;
1403          put += numberNonBasic + 1;
1404        }
[451]1405#endif
[2385]1406        djNorm0 = CoinMax(djNorm, 1.0e-20);
1407        CoinMemcpyN(work, numberTotal, dArray);
1408        CoinMemcpyN(work, numberTotal, dArray2);
1409        if (sequenceIn_ >= 0 && numberNonBasic == 1) {
1410          // see if simple move
1411          double objTheta2 = objective_->stepLength(this, solution_, work, 1.0e30,
1412            currentObj, predictedObj, thetaObj);
1413          rowArray->clear();
[1525]1414
[2385]1415          // update the incoming column
1416          unpackPacked(rowArray);
1417          factorization_->updateColumnFT(spare, rowArray);
1418          theta_ = 0.0;
1419          double *work2 = rowArray->denseVector();
1420          int number = rowArray->getNumElements();
1421          int *which2 = rowArray->getIndices();
1422          int iIndex;
1423          bool easyMove = false;
1424          double way;
1425          if (dj_[sequenceIn_] > 0.0)
1426            way = -1.0;
1427          else
1428            way = 1.0;
1429          double largest = COIN_DBL_MAX;
[1929]1430#ifdef CLP_DEBUG
[2385]1431          int kPivot = -1;
[1929]1432#endif
[2385]1433          for (iIndex = 0; iIndex < number; iIndex++) {
1434            int iRow = which2[iIndex];
1435            double alpha = way * work2[iIndex];
1436            int iPivot = pivotVariable_[iRow];
1437            if (alpha < -1.0e-5) {
1438              double distance = upper_[iPivot] - solution_[iPivot];
1439              if (distance < -largest * alpha) {
[1929]1440#ifdef CLP_DEBUG
[2385]1441                kPivot = iPivot;
[1929]1442#endif
[2385]1443                largest = CoinMax(0.0, -distance / alpha);
1444              }
1445              if (distance < -1.0e-12 * alpha) {
1446                easyMove = true;
1447                break;
1448              }
1449            } else if (alpha > 1.0e-5) {
1450              double distance = solution_[iPivot] - lower_[iPivot];
1451              if (distance < largest * alpha) {
[1929]1452#ifdef CLP_DEBUG
[2385]1453                kPivot = iPivot;
[1929]1454#endif
[2385]1455                largest = CoinMax(0.0, distance / alpha);
1456              }
1457              if (distance < 1.0e-12 * alpha) {
1458                easyMove = true;
1459                break;
1460              }
1461            }
1462          }
1463          // But largest has to match up with change
1464          assert(work[sequenceIn_]);
1465          largest /= fabs(work[sequenceIn_]);
1466          if (largest < objTheta2) {
1467            easyMove = true;
1468          } else if (!easyMove) {
1469            double objDrop = currentObj - predictedObj;
1470            double th = objective_->stepLength(this, solution_, work, largest,
1471              currentObj, predictedObj, simpleObjective);
1472            simpleObjective = CoinMax(simpleObjective, predictedObj);
1473            double easyDrop = currentObj - simpleObjective;
1474            if (easyDrop > 1.0e-8 && easyDrop > 0.5 * objDrop) {
1475              easyMove = true;
[393]1476#ifdef CLP_DEBUG
[2385]1477              if (handler_->logLevel() & 32)
1478                printf("easy - obj drop %g, easy drop %g\n", objDrop, easyDrop);
[393]1479#endif
[2385]1480              if (easyDrop > objDrop) {
1481                // debug
1482                printf("****** th %g simple %g\n", th, simpleObjective);
1483                objective_->stepLength(this, solution_, work, 1.0e30,
1484                  currentObj, predictedObj, simpleObjective);
1485                objective_->stepLength(this, solution_, work, largest,
1486                  currentObj, predictedObj, simpleObjective);
1487              }
1488            }
1489          }
1490          rowArray->clear();
[393]1491#ifdef CLP_DEBUG
[2385]1492          if (handler_->logLevel() & 32)
1493            printf("largest %g piv %d\n", largest, kPivot);
[393]1494#endif
[2385]1495          if (easyMove) {
1496            valueIn_ = solution_[sequenceIn_];
1497            dualIn_ = dj_[sequenceIn_];
1498            lowerIn_ = lower_[sequenceIn_];
1499            upperIn_ = upper_[sequenceIn_];
1500            if (dualIn_ > 0.0)
1501              directionIn_ = -1;
1502            else
1503              directionIn_ = 1;
1504            longArray->clear();
1505            pivotRow_ = -1;
1506            theta_ = 0.0;
1507            return 0;
1508          }
1509        }
1510      } else {
1511#if MINTYPE == 1
1512        if (conjugate) {
1513          double djNorm2 = djNorm;
1514          if (numberNonBasic && false) {
1515            int iIndex;
1516            djNorm2 = 0.0;
1517            for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
1518              int iSequence = which[iIndex];
1519              double alpha = work[iSequence];
1520              //djNorm2 += alpha*alpha;
1521              double alpha2 = work[iSequence] - dArray2[iSequence];
1522              djNorm2 += alpha * alpha2;
1523            }
1524            //printf("a %.18g b %.18g\n",djNorm,djNorm2);
1525          }
1526          djNorm = djNorm2;
1527          double beta = djNorm2 / djNormSave;
1528          // reset beta every so often
1529          //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1)
1530          //beta=0.0;
1531          //printf("current norm %g, old %g - beta %g\n",
1532          //     djNorm,djNormSave,beta);
1533          for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1534            dArray[iSequence] = work[iSequence] + beta * dArray[iSequence];
1535            dArray2[iSequence] = work[iSequence];
1536          }
1537        } else {
1538          for (iSequence = 0; iSequence < numberTotal; iSequence++)
1539            dArray[iSequence] = work[iSequence];
1540        }
[451]1541#else
[2385]1542        int number2 = numberNonBasic;
1543        if (number2) {
1544          // pack down into dArray
1545          int jLast = -1;
1546          for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
1547            int j = which[iSequence];
1548            assert(j > jLast);
1549            jLast = j;
1550            double value = work[j];
1551            dArray[iSequence] = -value;
1552          }
1553          // see whether to restart
1554          // check signs - gradient
1555          double g1 = innerProduct(dArray, number2, dArray);
1556          double g2 = innerProduct(dArray, number2, saveY2);
1557          // Get differences
1558          for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
1559            saveY2[iSequence] = dArray[iSequence] - saveY2[iSequence];
1560            saveS2[iSequence] = solution_[iSequence] - saveS2[iSequence];
1561          }
1562          double g3 = innerProduct(saveS2, number2, saveY2);
1563          printf("inner %g\n", g3);
1564          //assert(g3>0);
1565          double zzz[10000];
1566          int iVariable;
1567          g2 = 1.0e50; // temp
1568          if (fabs(g2) >= 0.2 * fabs(g1)) {
1569            // restart
1570            double delta = innerProduct(saveY2, number2, saveS2) / innerProduct(saveY2, number2, saveY2);
1571            delta = 1.0; //temp
1572            memset(zz, 0, number2 * sizeof(double));
1573            int put = 0;
1574            for (iVariable = 0; iVariable < number2; iVariable++) {
1575              zz[put] = delta;
1576              put += number2 + 1;
1577            }
1578          } else {
1579          }
1580          CoinMemcpyN(zz, number2 * number2, zzz);
1581          double ww[100];
1582          // get sk -Hkyk
1583          for (iVariable = 0; iVariable < number2; iVariable++) {
1584            double value = 0.0;
1585            for (int jVariable = 0; jVariable < number2; jVariable++) {
1586              value += saveY2[jVariable] * zzz[iVariable + number2 * jVariable];
1587            }
1588            ww[iVariable] = saveS2[iVariable] - value;
1589          }
1590          double ys = innerProduct(saveY2, number2, saveS2);
1591          double multiplier1 = 1.0 / ys;
1592          double multiplier2 = innerProduct(saveY2, number2, ww) / (ys * ys);
[451]1593#if 1
[2385]1594          // and second way
1595          // Hy
1596          double h[100];
1597          for (iVariable = 0; iVariable < number2; iVariable++) {
1598            double value = 0.0;
1599            for (int jVariable = 0; jVariable < number2; jVariable++) {
1600              value += saveY2[jVariable] * zzz[iVariable + number2 * jVariable];
1601            }
1602            h[iVariable] = value;
1603          }
1604          double hh[10000];
1605          double yhy1 = innerProduct(h, number2, saveY2) * multiplier1 + 1.0;
1606          yhy1 *= multiplier1;
1607          for (iVariable = 0; iVariable < number2; iVariable++) {
1608            for (int jVariable = 0; jVariable < number2; jVariable++) {
1609              int put = iVariable + number2 * jVariable;
1610              double value = zzz[put];
1611              value += yhy1 * saveS2[iVariable] * saveS2[jVariable];
1612              hh[put] = value;
1613            }
1614          }
1615          for (iVariable = 0; iVariable < number2; iVariable++) {
1616            for (int jVariable = 0; jVariable < number2; jVariable++) {
1617              int put = iVariable + number2 * jVariable;
1618              double value = hh[put];
1619              value -= multiplier1 * (saveS2[iVariable] * h[jVariable]);
1620              value -= multiplier1 * (saveS2[jVariable] * h[iVariable]);
1621              hh[put] = value;
1622            }
1623          }
[451]1624#endif
[2385]1625          // now update H
1626          for (iVariable = 0; iVariable < number2; iVariable++) {
1627            for (int jVariable = 0; jVariable < number2; jVariable++) {
1628              int put = iVariable + number2 * jVariable;
1629              double value = zzz[put];
1630              value += multiplier1 * (ww[iVariable] * saveS2[jVariable] + ww[jVariable] * saveS2[iVariable]);
1631              value -= multiplier2 * saveS2[iVariable] * saveS2[jVariable];
1632              zzz[put] = value;
1633            }
1634          }
1635          //memcpy(zzz,hh,size*sizeof(double));
1636          // do search direction
1637          memset(dArray, 0, numberTotal * sizeof(double));
1638          for (iVariable = 0; iVariable < numberNonBasic; iVariable++) {
1639            double value = 0.0;
1640            for (int jVariable = 0; jVariable < number2; jVariable++) {
1641              int k = which[jVariable];
1642              value += work[k] * zzz[iVariable + number2 * jVariable];
1643            }
1644            int i = which[iVariable];
1645            dArray[i] = value;
1646          }
1647          // Now fill out dArray
1648          {
1649            int j;
1650            // Now do basic ones
1651            int iRow;
1652            CoinIndexedVector *spare1 = spare;
1653            CoinIndexedVector *spare2 = rowArray;
[451]1654#if CLP_DEBUG > 1
[2385]1655            spare1->checkClear();
1656            spare2->checkClear();
[451]1657#endif
[2385]1658            double *array2 = spare1->denseVector();
1659            int *index2 = spare1->getIndices();
1660            int number2 = 0;
1661            times(-1.0, dArray, array2);
1662            dArray = dArray + numberColumns_;
1663            for (iRow = 0; iRow < numberRows_; iRow++) {
1664              double value = array2[iRow] + dArray[iRow];
1665              if (value) {
1666                array2[iRow] = value;
1667                index2[number2++] = iRow;
1668              } else {
1669                array2[iRow] = 0.0;
1670              }
1671            }
1672            dArray -= numberColumns_;
1673            spare1->setNumElements(number2);
1674            // Ftran
1675            factorization_->updateColumn(spare2, spare1);
1676            number2 = spare1->getNumElements();
1677            for (j = 0; j < number2; j++) {
1678              int iSequence = index2[j];
1679              double value = array2[iSequence];
1680              array2[iSequence] = 0.0;
1681              if (value) {
1682                int iPivot = pivotVariable_[iSequence];
1683                double oldValue = dArray[iPivot];
1684                dArray[iPivot] = value + oldValue;
1685              }
1686            }
1687            spare1->setNumElements(0);
1688          }
1689          //assert (innerProductIndexed(dArray,number2,work,which)>0);
1690          //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
1691          printf("innerP1 %g innerP2 %g\n", innerProduct(dArray, numberTotal, work),
1692            innerProductIndexed(dArray, numberNonBasic, work, which));
1693          assert(innerProduct(dArray, numberTotal, work) > 0);
[451]1694#if 1
[2385]1695          {
1696            // check null vector
1697            int iRow;
1698            double qq[10];
1699            memset(qq, 0, numberRows_ * sizeof(double));
1700            multiplyAdd(dArray + numberColumns_, numberRows_, -1.0, qq, 0.0);
1701            matrix_->times(1.0, dArray, qq, rowScale_, columnScale_);
1702            double largest = 0.0;
1703            int iLargest = -1;
1704            for (iRow = 0; iRow < numberRows_; iRow++) {
1705              double value = fabs(qq[iRow]);
1706              if (value > largest) {
1707                largest = value;
1708                iLargest = iRow;
1709              }
1710            }
1711            printf("largest null error %g on row %d\n", largest, iLargest);
1712            for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1713              if (getStatus(iSequence) == basic)
1714                assert(fabs(dj_[iSequence]) < 1.0e-3);
1715            }
1716          }
[451]1717#endif
[2385]1718          CoinMemcpyN(saveY2, numberNonBasic, saveY);
1719          CoinMemcpyN(saveS2, numberNonBasic, saveS);
1720        }
[451]1721#endif
[2385]1722      }
1723#if MINTYPE == 2
1724      for (iSequence = 0; iSequence < numberNonBasic; iSequence++) {
1725        int j = which[iSequence];
1726        saveY2[iSequence] = -work[j];
1727        saveS2[iSequence] = solution_[j];
1728      }
[451]1729#endif
[2385]1730      if (djNorm < eps * djNorm0 || (nPasses > 100 && djNorm < CoinMin(1.0e-1 * djNorm0, 1.0e-12))) {
[393]1731#ifdef CLP_DEBUG
[2385]1732        if (handler_->logLevel() & 32)
1733          printf("dj norm reduced from %g to %g\n", djNorm0, djNorm);
[393]1734#endif
[2385]1735        if (pivotMode < 10 || !numberNonBasic) {
1736          finished = true;
1737        } else {
1738          localPivotMode = pivotMode;
1739          nTotalPasses += nPasses;
1740          nPasses = 0;
1741          continue;
1742        }
1743      }
1744      //if (nPasses==100||nPasses==50)
1745      //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
1746      //         normFlagged);
1747      if (nPasses > 100 && djNorm < 1.0e-2 * normFlagged && !startLocalMode) {
[393]1748#ifdef CLP_DEBUG
[2385]1749        if (handler_->logLevel() & 32)
1750          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n", djNorm0, djNorm,
1751            normFlagged);
[393]1752#endif
[2385]1753        unflag();
1754        localPivotMode = 0;
1755        nTotalPasses += nPasses;
1756        nPasses = 0;
1757        continue;
1758      }
1759      if (djNorm > 0.99 * djNorm0 && nPasses > 1500) {
1760        finished = true;
[393]1761#ifdef CLP_DEBUG
[2385]1762        if (handler_->logLevel() & 32)
1763          printf("dj norm NOT reduced from %g to %g\n", djNorm0, djNorm);
[393]1764#endif
[2385]1765        djNorm = 1.2345e-20;
1766      }
1767      number = longArray->getNumElements();
1768      if (!numberNonBasic) {
1769        // looks optimal
1770        // check if any dj look plausible
1771        int nSuper = 0;
1772        int nFlagged = 0;
1773        int nNormal = 0;
1774        for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
1775          if (flagged(iSequence)) {
1776            switch (getStatus(iSequence)) {
[1525]1777
[2385]1778            case basic:
1779            case ClpSimplex::isFixed:
1780              break;
1781            case atUpperBound:
1782              if (dj_[iSequence] > dualTolerance_)
1783                nFlagged++;
1784              break;
1785            case atLowerBound:
1786              if (dj_[iSequence] < -dualTolerance_)
1787                nFlagged++;
1788              break;
1789            case isFree:
1790            case superBasic:
1791              if (fabs(dj_[iSequence]) > dualTolerance_)
1792                nFlagged++;
1793              break;
1794            }
1795            continue;
1796          }
1797          switch (getStatus(iSequence)) {
[1525]1798
[2385]1799          case basic:
1800          case ClpSimplex::isFixed:
1801            break;
1802          case atUpperBound:
1803            if (dj_[iSequence] > dualTolerance_)
1804              nNormal++;
1805            break;
1806          case atLowerBound:
1807            if (dj_[iSequence] < -dualTolerance_)
1808              nNormal++;
1809            break;
1810          case isFree:
1811          case superBasic:
1812            if (fabs(dj_[iSequence]) > dualTolerance_)
1813              nSuper++;
1814            break;
1815          }
1816        }
[393]1817#ifdef CLP_DEBUG
[2385]1818        if (handler_->logLevel() & 32)
1819          printf("%d super, %d normal, %d flagged\n",
1820            nSuper, nNormal, nFlagged);
[393]1821#endif
[1197]1822
[2385]1823        int nFlagged2 = 1;
1824        if (lastSequenceIn < 0 && !nNormal && !nSuper) {
1825          nFlagged2 = unflag();
1826          if (pivotMode >= 10) {
1827            pivotMode--;
[393]1828#ifdef CLP_DEBUG
[2385]1829            if (handler_->logLevel() & 32)
1830              printf("pivot mode now %d\n", pivotMode);
[393]1831#endif
[2385]1832            if (pivotMode == 9)
1833              pivotMode = 0; // switch off fast attempt
1834          }
1835        } else {
1836          lastSequenceIn = -1;
1837        }
1838        if (!nFlagged2 || (normFlagged + normUnflagged < 1.0e-8)) {
1839          primalColumnPivot_->setLooksOptimal(true);
1840          return 0;
1841        } else {
1842          localPivotMode = -1;
1843          nTotalPasses += nPasses;
1844          nPasses = 0;
1845          finished = false;
1846          continue;
1847        }
1848      }
1849      bestSequence = -1;
1850      bestBasicSequence = -1;
[1525]1851
[2385]1852      // temp
1853      nPasses++;
1854      if (nPasses > 2000)
1855        finished = true;
1856      double theta = 1.0e30;
1857      basicTheta = 1.0e30;
1858      theta_ = 1.0e30;
1859      double basicTolerance = 1.0e-4 * primalTolerance_;
1860      for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1861        //if (flagged(iSequence)
1862        //  continue;
1863        double alpha = dArray[iSequence];
1864        Status thisStatus = getStatus(iSequence);
1865        double oldValue = solution_[iSequence];
1866        if (thisStatus != basic) {
1867          if (fabs(alpha) >= acceptablePivot) {
1868            if (alpha < 0.0) {
1869              // variable going towards lower bound
1870              double bound = lower_[iSequence];
1871              oldValue -= bound;
1872              if (oldValue + theta * alpha < 0.0) {
1873                bestSequence = iSequence;
1874                theta = CoinMax(0.0, oldValue / (-alpha));
1875              }
1876            } else {
1877              // variable going towards upper bound
1878              double bound = upper_[iSequence];
1879              oldValue = bound - oldValue;
1880              if (oldValue - theta * alpha < 0.0) {
1881                bestSequence = iSequence;
1882                theta = CoinMax(0.0, oldValue / alpha);
1883              }
1884            }
1885          }
1886        } else {
1887          if (fabs(alpha) >= acceptableBasic) {
1888            if (alpha < 0.0) {
1889              // variable going towards lower bound
1890              double bound = lower_[iSequence];
1891              oldValue -= bound;
1892              if (oldValue + basicTheta * alpha < -basicTolerance) {
1893                bestBasicSequence = iSequence;
1894                basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / (-alpha));
1895              }
1896            } else {
1897              // variable going towards upper bound
1898              double bound = upper_[iSequence];
1899              oldValue = bound - oldValue;
1900              if (oldValue - basicTheta * alpha < -basicTolerance) {
1901                bestBasicSequence = iSequence;
1902                basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / alpha);
1903              }
1904            }
1905          }
1906        }
1907      }
1908      theta_ = CoinMin(theta, basicTheta);
1909      // Now find minimum of function
1910      double objTheta2 = objective_->stepLength(this, solution_, dArray, CoinMin(theta, basicTheta),
1911        currentObj, predictedObj, thetaObj);
[393]1912#ifdef CLP_DEBUG
[2385]1913      if (handler_->logLevel() & 32)
1914        printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
[393]1915#endif
[2385]1916      objTheta2 = CoinMin(objTheta2, 1.0e29);
1917#if MINTYPE == 1
1918      if (conjugate) {
1919        double offset;
1920        const double *gradient = objective_->gradient(this,
1921          dArray, offset,
1922          true, 0);
1923        double product = 0.0;
1924        for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
1925          double alpha = dArray[iSequence];
1926          double value = alpha * gradient[iSequence];
1927          product += value;
1928        }
1929        //#define INCLUDESLACK
[393]1930#ifdef INCLUDESLACK
[2385]1931        for (; iSequence < numberColumns_ + numberRows_; iSequence++) {
1932          double alpha = dArray[iSequence];
1933          double value = alpha * cost_[iSequence];
1934          product += value;
1935        }
[393]1936#endif
[2385]1937        if (product > 0.0)
1938          objTheta = djNorm / product;
1939        else
1940          objTheta = COIN_DBL_MAX;
[521]1941#ifndef NDEBUG
[2385]1942        if (product < -1.0e-8 && handler_->logLevel() > 1)
1943          printf("bad product %g\n", product);
[521]1944#endif
[2385]1945        product = CoinMax(product, 0.0);
1946      } else {
1947        objTheta = objTheta2;
1948      }
[451]1949#else
[2385]1950      objTheta = objTheta2;
[451]1951#endif
[2385]1952      // if very small difference then take pivot (depends on djNorm?)
1953      // distinguish between basic and non-basic
1954      bool chooseObjTheta = objTheta < theta_;
1955      if (chooseObjTheta) {
1956        if (thetaObj < currentObj - 1.0e-12 && objTheta + 1.0e-10 > theta_)
1957          chooseObjTheta = false;
1958        //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
1959        //chooseObjTheta=false;
1960      }
1961      //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
1962      if (chooseObjTheta) {
1963        theta_ = objTheta;
1964      } else {
1965        objTheta = CoinMax(objTheta, 1.00000001 * theta_ + 1.0e-12);
1966        //if (theta+1.0e-13>basicTheta) {
1967        //theta = CoinMax(theta,1.00000001*basicTheta);
1968        //theta_ = basicTheta;
1969        //}
1970      }
1971      // always out if one variable in and zero move
1972      if (theta_ == basicTheta || (sequenceIn_ >= 0 && theta_ < 1.0e-10))
1973        finished = true; // come out
[393]1974#ifdef CLP_DEBUG
[2385]1975      if (handler_->logLevel() & 32) {
1976        printf("%d pass,", nPasses);
1977        if (sequenceIn_ >= 0)
1978          printf(" Sin %d,", sequenceIn_);
1979        if (basicTheta == theta_)
1980          printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
1981        else
1982          printf(" basicTheta %g", basicTheta);
1983        if (theta == theta_)
1984          printf(" X(%d) non-basicTheta %g", bestSequence, theta);
1985        else
1986          printf(" non-basicTheta %g", theta);
1987        printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
1988        printf(" djNorm %g\n", djNorm);
1989      }
[393]1990#endif
[2385]1991      if (handler_->logLevel() > 3 && objTheta != theta_) {
1992        printf("%d pass obj %g,", nPasses, currentObj);
1993        if (sequenceIn_ >= 0)
1994          printf(" Sin %d,", sequenceIn_);
1995        if (basicTheta == theta_)
1996          printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
1997        else
1998          printf(" basicTheta %g", basicTheta);
1999        if (theta == theta_)
2000          printf(" X(%d) non-basicTheta %g", bestSequence, theta);
2001        else
2002          printf(" non-basicTheta %g", theta);
2003        printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
2004        printf(" djNorm %g\n", djNorm);
2005      }
2006      if (objTheta != theta_) {
2007        //printf("hit boundary after %d passes\n",nPasses);
2008        nTotalPasses += nPasses;
2009        nPasses = 0; // start again
2010      }
2011      if (localPivotMode < 10 || lastSequenceIn == bestSequence) {
2012        if (theta_ == theta && theta_ < basicTheta && theta_ < 1.0e-5)
2013          setFlagged(bestSequence); // out of active set
2014      }
2015      // Update solution
2016      for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2017        //for (iIndex=0;iIndex<number;iIndex++) {
[1197]2018
[2385]2019        //int iSequence = which[iIndex];
2020        double alpha = dArray[iSequence];
2021        if (alpha) {
2022          double value = solution_[iSequence] + theta_ * alpha;
2023          solution_[iSequence] = value;
2024          switch (getStatus(iSequence)) {
[1525]2025
[2385]2026          case basic:
2027          case isFixed:
2028          case isFree:
2029          case atUpperBound:
2030          case atLowerBound:
2031            nonLinearCost_->setOne(iSequence, value);
2032            break;
2033          case superBasic:
2034            // To get correct action
2035            setStatus(iSequence, isFixed);
2036            nonLinearCost_->setOne(iSequence, value);
2037            //assert (getStatus(iSequence)!=isFixed);
2038            break;
2039          }
2040        }
2041      }
2042      if (objTheta < SMALLTHETA2 && objTheta == theta_) {
2043        if (sequenceIn_ >= 0 && basicTheta < 1.0e-9) {
2044          // try just one
2045          localPivotMode = 0;
2046          sequenceIn_ = -1;
2047          nTotalPasses += nPasses;
2048          nPasses = 0;
2049          //finished=true;
2050          //objTheta=0.0;
2051          //theta_=0.0;
2052        } else if (sequenceIn_ < 0 && nTotalPasses > 10) {
2053          if (objTheta < 1.0e-10) {
2054            finished = true;
2055            //printf("zero move\n");
2056            break;
2057          }
2058        }
2059      }
[393]2060#ifdef CLP_DEBUG
[2385]2061      if (handler_->logLevel() & 32) {
2062        objective_->stepLength(this, solution_, work, 0.0,
2063          currentObj, predictedObj, thetaObj);
2064        printf("current obj %g after update - simple was %g\n", currentObj, simpleObjective);
2065      }
[393]2066#endif
[2385]2067      if (sequenceIn_ >= 0 && !finished && objTheta > 1.0e-4) {
2068        // we made some progress - back to normal
2069        if (localPivotMode < 10) {
2070          localPivotMode = 0;
2071          sequenceIn_ = -1;
2072          nTotalPasses += nPasses;
2073          nPasses = 0;
2074        }
[393]2075#ifdef CLP_DEBUG
[2385]2076        if (handler_->logLevel() & 32)
2077          printf("some progress?\n");
[393]2078#endif
[2385]2079      }
[393]2080#if CLP_DEBUG > 1
[2385]2081      longArray->checkClean();
[393]2082#endif
[2385]2083    }
[393]2084#ifdef CLP_DEBUG
[2385]2085    if (handler_->logLevel() & 32)
2086      printf("out of loop after %d (%d) passes\n", nPasses, nTotalPasses);
[393]2087#endif
[2385]2088    if (nTotalPasses >= 1000 || (nTotalPasses > 10 && sequenceIn_ < 0 && theta_ < 1.0e-10))
2089      returnCode = 2;
2090    bool ordinaryDj = false;
2091    //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
2092    //printf("could try ordinary iteration %g\n",theta_);
2093    if (sequenceIn_ >= 0 && numberNonBasic == 1 && theta_ < 1.0e-15) {
2094      //printf("trying ordinary iteration\n");
2095      ordinaryDj = true;
2096    }
2097    if (!basicTheta && !ordinaryDj) {
2098      //returnCode=2;
2099      //objTheta=-1.0; // so we fall through
2100    }
[2438]2101    if (theta_ >= 1.0e30) // odd
2102      ordinaryDj=true;
[2385]2103    // See if we need to pivot
2104    if (theta_ == basicTheta || ordinaryDj) {
2105      if (!ordinaryDj) {
2106        // find an incoming column which will force pivot
2107        int iRow;
2108        pivotRow_ = -1;
2109        for (iRow = 0; iRow < numberRows_; iRow++) {
2110          if (pivotVariable_[iRow] == bestBasicSequence) {
2111            pivotRow_ = iRow;
2112            break;
[1525]2113          }
[2385]2114        }
2115        assert(pivotRow_ >= 0);
2116        // Get good size for pivot
2117        double acceptablePivot = 1.0e-7;
2118        if (factorization_->pivots() > 10)
2119          acceptablePivot = 1.0e-5; // if we have iterated be more strict
2120        else if (factorization_->pivots() > 5)
2121          acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
2122        // should be dArray but seems better this way!
2123        double direction = work[bestBasicSequence] > 0.0 ? -1.0 : 1.0;
2124        // create as packed
2125        rowArray->createPacked(1, &pivotRow_, &direction);
2126        factorization_->updateColumnTranspose(spare, rowArray);
2127        // put row of tableau in rowArray and columnArray
2128        matrix_->transposeTimes(this, -1.0,
2129          rowArray, spare, columnArray);
2130        // choose one futhest away from bound which has reasonable pivot
2131        // If increasing we want negative alpha
[1197]2132
[2385]2133        double *work2;
2134        int iSection;
[1197]2135
[2385]2136        sequenceIn_ = -1;
2137        double bestValue = -1.0;
2138        double bestDirection = 0.0;
2139        // First pass we take correct direction and large pivots
2140        // then correct direction
2141        // then any
2142        double check[] = { 1.0e-8, -1.0e-12, -1.0e30 };
2143        double mult[] = { 100.0, 1.0, 1.0 };
2144        for (int iPass = 0; iPass < 3; iPass++) {
2145          //if (!bestValue&&iPass==2)
2146          //bestValue=-1.0;
2147          double acceptable = acceptablePivot * mult[iPass];
2148          double checkValue = check[iPass];
2149          for (iSection = 0; iSection < 2; iSection++) {
[1197]2150
[2385]2151            int addSequence;
[1525]2152
[2385]2153            if (!iSection) {
2154              work2 = rowArray->denseVector();
2155              number = rowArray->getNumElements();
2156              which = rowArray->getIndices();
2157              addSequence = numberColumns_;
2158            } else {
2159              work2 = columnArray->denseVector();
2160              number = columnArray->getNumElements();
2161              which = columnArray->getIndices();
2162              addSequence = 0;
2163            }
2164            int i;
[1525]2165
[2385]2166            for (i = 0; i < number; i++) {
2167              int iSequence = which[i] + addSequence;
2168              if (flagged(iSequence))
2169                continue;
2170              //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
2171              //                  upper_[iSequence]-solution_[iSequence]);
2172              double alpha = work2[i];
2173              // should be dArray but seems better this way!
2174              double change = work[iSequence];
2175              Status thisStatus = getStatus(iSequence);
2176              double direction = 0;
2177              ;
2178              switch (thisStatus) {
[1525]2179
[2385]2180              case basic:
2181              case ClpSimplex::isFixed:
2182                break;
2183              case isFree:
2184              case superBasic:
2185                if (alpha < -acceptable && change > checkValue)
2186                  direction = 1.0;
2187                else if (alpha > acceptable && change < -checkValue)
2188                  direction = -1.0;
2189                break;
2190              case atUpperBound:
2191                if (alpha > acceptable && change < -checkValue)
2192                  direction = -1.0;
2193                else if (iPass == 2 && alpha < -acceptable && change < -checkValue)
2194                  direction = 1.0;
2195                break;
2196              case atLowerBound:
2197                if (alpha < -acceptable && change > checkValue)
2198                  direction = 1.0;
2199                else if (iPass == 2 && alpha > acceptable && change > checkValue)
2200                  direction = -1.0;
2201                break;
2202              }
2203              if (direction) {
2204                if (sequenceIn_ != lastSequenceIn || localPivotMode < 10) {
2205                  if (CoinMin(solution_[iSequence] - lower_[iSequence],
2206                        upper_[iSequence] - solution_[iSequence])
2207                    > bestValue) {
2208                    bestValue = CoinMin(solution_[iSequence] - lower_[iSequence],
2209                      upper_[iSequence] - solution_[iSequence]);
2210                    sequenceIn_ = iSequence;
2211                    bestDirection = direction;
2212                  }
2213                } else {
2214                  // choose
2215                  bestValue = COIN_DBL_MAX;
2216                  sequenceIn_ = iSequence;
2217                  bestDirection = direction;
2218                }
2219              }
2220            }
2221          }
2222          if (sequenceIn_ >= 0 && bestValue > 0.0)
2223            break;
2224        }
2225        if (sequenceIn_ >= 0) {
2226          valueIn_ = solution_[sequenceIn_];
2227          dualIn_ = dj_[sequenceIn_];
2228          if (bestDirection < 0.0) {
2229            // we want positive dj
2230            if (dualIn_ <= 0.0) {
2231              //printf("bad dj - xx %g\n",dualIn_);
2232              dualIn_ = 1.0;
2233            }
2234          } else {
2235            // we want negative dj
2236            if (dualIn_ >= 0.0) {
2237              //printf("bad dj - xx %g\n",dualIn_);
2238              dualIn_ = -1.0;
2239            }
2240          }
2241          lowerIn_ = lower_[sequenceIn_];
2242          upperIn_ = upper_[sequenceIn_];
2243          if (dualIn_ > 0.0)
2244            directionIn_ = -1;
2245          else
2246            directionIn_ = 1;
2247        } else {
2248          //ordinaryDj=true;
[393]2249#ifdef CLP_DEBUG
[2385]2250          if (handler_->logLevel() & 32) {
2251            printf("no easy pivot - norm %g mode %d\n", djNorm, localPivotMode);
2252            if (rowArray->getNumElements() + columnArray->getNumElements() < 12) {
2253              for (iSection = 0; iSection < 2; iSection++) {
[1197]2254
[2385]2255                int addSequence;
[1197]2256
[2385]2257                if (!iSection) {
2258                  work2 = rowArray->denseVector();
2259                  number = rowArray->getNumElements();
2260                  which = rowArray->getIndices();
2261                  addSequence = numberColumns_;
2262                } else {
2263                  work2 = columnArray->denseVector();
2264                  number = columnArray->getNumElements();
2265                  which = columnArray->getIndices();
2266                  addSequence = 0;
2267                }
2268                int i;
2269                char section[] = { 'R', 'C' };
2270                for (i = 0; i < number; i++) {
2271                  int iSequence = which[i] + addSequence;
2272                  if (flagged(iSequence)) {
2273                    printf("%c%d flagged\n", section[iSection], which[i]);
2274                    continue;
2275                  } else {
2276                    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
2277                      section[iSection], which[i], status_[iSequence],
2278                      lower_[iSequence], solution_[iSequence], upper_[iSequence],
2279                      work2[i], work[iSequence]);
2280                  }
2281                }
2282              }
2283            }
2284          }
[393]2285#endif
[2385]2286          assert(sequenceIn_ < 0);
2287          justOne = true;
2288          allFinished = false; // Round again
2289          finished = false;
2290          nTotalPasses += nPasses;
2291          nPasses = 0;
2292          if (djNorm < 0.9 * djNorm0 && djNorm < 1.0e-3 && !localPivotMode) {
[393]2293#ifdef CLP_DEBUG
[2385]2294            if (handler_->logLevel() & 32)
2295              printf("no pivot - mode %d norms %g %g - unflagging\n",
2296                localPivotMode, djNorm0, djNorm);
[393]2297#endif
[2385]2298            unflag(); //unflagging
2299            returnCode = 1;
[1525]2300          } else {
[2385]2301            returnCode = 2; // do single incoming
2302            returnCode = 1;
2303          }
2304        }
2305      }
2306      rowArray->clear();
2307      columnArray->clear();
2308      longArray->clear();
2309      if (ordinaryDj) {
2310        valueIn_ = solution_[sequenceIn_];
2311        dualIn_ = dj_[sequenceIn_];
2312        lowerIn_ = lower_[sequenceIn_];
2313        upperIn_ = upper_[sequenceIn_];
2314        if (dualIn_ > 0.0)
2315          directionIn_ = -1;
2316        else
2317          directionIn_ = 1;
2318      }
2319      if (returnCode == 1)
2320        returnCode = 0;
2321    } else {
2322      // round again
2323      longArray->clear();
2324      if (djNorm < 1.0e-3 && !localPivotMode) {
2325        if (djNorm == 1.2345e-20 && djNorm0 > 1.0e-4) {
[393]2326#ifdef CLP_DEBUG
[2385]2327          if (handler_->logLevel() & 32)
2328            printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n", djNorm0, nPasses,
2329              localPivotMode, returnCode);
[393]2330#endif
[2385]2331          //if (!localPivotMode)
2332          //returnCode=2; // force singleton
2333        } else {
[393]2334#ifdef CLP_DEBUG
[2385]2335          if (handler_->logLevel() & 32)
2336            printf("unflagging as djNorm %g %g, %d passes\n", djNorm, djNorm0, nPasses);
[393]2337#endif
[2385]2338          if (pivotMode >= 10) {
2339            pivotMode--;
[393]2340#ifdef CLP_DEBUG
[2385]2341            if (handler_->logLevel() & 32)
2342              printf("pivot mode now %d\n", pivotMode);
[393]2343#endif
[2385]2344            if (pivotMode == 9)
2345              pivotMode = 0; // switch off fast attempt
[1525]2346          }
[2385]2347          bool unflagged = unflag() != 0;
2348          if (!unflagged && djNorm < 1.0e-10) {
2349            // ? declare victory
2350            sequenceIn_ = -1;
2351            returnCode = 0;
2352          }
2353        }
2354      }
2355    }
2356  }
2357  if (djNorm0 < 1.0e-12 * normFlagged) {
[393]2358#ifdef CLP_DEBUG
[2385]2359    if (handler_->logLevel() & 32)
2360      printf("unflagging as djNorm %g %g and flagged norm %g\n", djNorm, djNorm0, normFlagged);
[393]2361#endif
[2385]2362    unflag();
2363  }
2364  if (saveObj - currentObj < 1.0e-5 && nTotalPasses > 2000) {
2365    normUnflagged = 0.0;
2366    double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
2367    for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
2368      switch (getStatus(iSequence)) {
[1197]2369
[2385]2370      case basic:
2371      case ClpSimplex::isFixed:
2372        break;
2373      case atUpperBound:
2374        if (dj_[iSequence] > dualTolerance3)
2375          normFlagged += dj_[iSequence] * dj_[iSequence];
2376        break;
2377      case atLowerBound:
2378        if (dj_[iSequence] < -dualTolerance3)
2379          normFlagged += dj_[iSequence] * dj_[iSequence];
2380        break;
2381      case isFree:
2382      case superBasic:
2383        if (fabs(dj_[iSequence]) > dualTolerance3)
2384          normFlagged += dj_[iSequence] * dj_[iSequence];
2385        break;
2386      }
2387    }
2388    if (handler_->logLevel() > 2)
2389      printf("possible optimal  %d %d %g %g\n",
2390        nBigPasses, nTotalPasses, saveObj - currentObj, normFlagged);
2391    if (normFlagged < 1.0e-5) {
2392      unflag();
2393      primalColumnPivot_->setLooksOptimal(true);
2394      returnCode = 0;
2395    }
2396  }
2397  return returnCode;
[393]2398}
2399/* Do last half of an iteration.
2400   Return codes
2401   Reasons to come out normal mode
2402   -1 normal
2403   -2 factorize now - good iteration
2404   -3 slight inaccuracy - refactorize - iteration done
2405   -4 inaccuracy - refactorize - no iteration
2406   -5 something flagged - go round again
2407   +2 looks unbounded
2408   +3 max iterations (iteration done)
[1525]2409
[393]2410*/
[2385]2411int ClpSimplexNonlinear::pivotNonlinearResult()
[393]2412{
2413
[2385]2414  int returnCode = -1;
[393]2415
[2385]2416  rowArray_[1]->clear();
[1525]2417
[2385]2418  // we found a pivot column
2419  // update the incoming column
2420  unpackPacked(rowArray_[1]);
2421  factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
2422  theta_ = 0.0;
2423  double *work = rowArray_[1]->denseVector();
2424  int number = rowArray_[1]->getNumElements();
2425  int *which = rowArray_[1]->getIndices();
2426  bool keepValue = false;
2427  double saveValue = 0.0;
2428  if (pivotRow_ >= 0) {
2429    sequenceOut_ = pivotVariable_[pivotRow_];
2430    ;
2431    valueOut_ = solution(sequenceOut_);
2432    keepValue = true;
2433    saveValue = valueOut_;
2434    lowerOut_ = lower_[sequenceOut_];
2435    upperOut_ = upper_[sequenceOut_];
2436    int iIndex;
2437    for (iIndex = 0; iIndex < number; iIndex++) {
[1525]2438
[2385]2439      int iRow = which[iIndex];
2440      if (iRow == pivotRow_) {
2441        alpha_ = work[iIndex];
2442        break;
2443      }
2444    }
2445  } else {
2446    int iIndex;
2447    double smallest = COIN_DBL_MAX;
2448    for (iIndex = 0; iIndex < number; iIndex++) {
2449      int iRow = which[iIndex];
2450      double alpha = work[iIndex];
2451      if (fabs(alpha) > 1.0e-6) {
2452        int iPivot = pivotVariable_[iRow];
2453        double distance = CoinMin(upper_[iPivot] - solution_[iPivot],
2454          solution_[iPivot] - lower_[iPivot]);
2455        if (distance < smallest) {
2456          pivotRow_ = iRow;
2457          alpha_ = alpha;
2458          smallest = distance;
2459        }
2460      }
2461    }
2462    if (smallest > primalTolerance_) {
2463      smallest = COIN_DBL_MAX;
2464      for (iIndex = 0; iIndex < number; iIndex++) {
2465        int iRow = which[iIndex];
2466        double alpha = work[iIndex];
2467        if (fabs(alpha) > 1.0e-6) {
2468          double distance = randomNumberGenerator_.randomDouble();
2469          if (distance < smallest) {
2470            pivotRow_ = iRow;
2471            alpha_ = alpha;
2472            smallest = distance;
[1525]2473          }
[2385]2474        }
2475      }
2476    }
2477    assert(pivotRow_ >= 0);
2478    sequenceOut_ = pivotVariable_[pivotRow_];
2479    ;
2480    valueOut_ = solution(sequenceOut_);
2481    lowerOut_ = lower_[sequenceOut_];
2482    upperOut_ = upper_[sequenceOut_];
2483  }
2484  double newValue = valueOut_ - theta_ * alpha_;
2485  bool isSuperBasic = false;
2486  if (valueOut_ >= upperOut_ - primalTolerance_) {
2487    directionOut_ = -1; // to upper bound
2488    upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2489    upperOut_ = newValue;
2490  } else if (valueOut_ <= lowerOut_ + primalTolerance_) {
2491    directionOut_ = 1; // to lower bound
2492    lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2493  } else {
2494    lowerOut_ = valueOut_;
2495    upperOut_ = valueOut_;
2496    isSuperBasic = true;
2497    //printf("XX superbasic out\n");
2498  }
2499  dualOut_ = dj_[sequenceOut_];
2500  // if stable replace in basis
[1525]2501
[2385]2502  int updateStatus = factorization_->replaceColumn(this,
2503    rowArray_[2],
2504    rowArray_[1],
2505    pivotRow_,
2506    alpha_);
[1525]2507
[2385]2508  // if no pivots, bad update but reasonable alpha - take and invert
2509  if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
2510    updateStatus = 4;
2511  if (updateStatus == 1 || updateStatus == 4) {
2512    // slight error
2513    if (factorization_->pivots() > 5 || updateStatus == 4) {
2514      returnCode = -3;
2515    }
2516  } else if (updateStatus == 2) {
2517    // major error
2518    // better to have small tolerance even if slower
2519    factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
2520    int maxFactor = factorization_->maximumPivots();
2521    if (maxFactor > 10) {
2522      if (forceFactorization_ < 0)
2523        forceFactorization_ = maxFactor;
2524      forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
2525    }
2526    // later we may need to unwind more e.g. fake bounds
2527    if (lastGoodIteration_ != numberIterations_) {
2528      clearAll();
2529      pivotRow_ = -1;
2530      returnCode = -4;
2531    } else {
2532      // need to reject something
2533      char x = isColumn(sequenceIn_) ? 'C' : 'R';
2534      handler_->message(CLP_SIMPLEX_FLAG, messages_)
2535        << x << sequenceWithin(sequenceIn_)
2536        << CoinMessageEol;
2537      setFlagged(sequenceIn_);
2538      progress_.clearBadTimes();
2539      lastBadIteration_ = numberIterations_; // say be more cautious
2540      clearAll();
2541      pivotRow_ = -1;
2542      sequenceOut_ = -1;
2543      returnCode = -5;
2544    }
2545    return returnCode;
2546  } else if (updateStatus == 3) {
2547    // out of memory
2548    // increase space if not many iterations
2549    if (factorization_->pivots() < 0.5 * factorization_->maximumPivots() && factorization_->pivots() < 200)
2550      factorization_->areaFactor(
2551        factorization_->areaFactor() * 1.1);
2552    returnCode = -2; // factorize now
2553  } else if (updateStatus == 5) {
2554    problemStatus_ = -2; // factorize now
2555  }
[1525]2556
[2385]2557  // update primal solution
[1525]2558
[2385]2559  double objectiveChange = 0.0;
2560  // after this rowArray_[1] is not empty - used to update djs
2561  // If pivot row >= numberRows then may be gub
2562  updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, 1);
[1525]2563
[2385]2564  double oldValue = valueIn_;
2565  if (directionIn_ == -1) {
2566    // as if from upper bound
2567    if (sequenceIn_ != sequenceOut_) {
2568      // variable becoming basic
2569      valueIn_ -= fabs(theta_);
2570    } else {
2571      valueIn_ = lowerIn_;
2572    }
2573  } else {
2574    // as if from lower bound
2575    if (sequenceIn_ != sequenceOut_) {
2576      // variable becoming basic
2577      valueIn_ += fabs(theta_);
2578    } else {
2579      valueIn_ = upperIn_;
2580    }
2581  }
2582  objectiveChange += dualIn_ * (valueIn_ - oldValue);
2583  // outgoing
2584  if (sequenceIn_ != sequenceOut_) {
2585    if (directionOut_ > 0) {
2586      valueOut_ = lowerOut_;
2587    } else {
2588      valueOut_ = upperOut_;
2589    }
2590    if (valueOut_ < lower_[sequenceOut_] - primalTolerance_)
2591      valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
2592    else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
2593      valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
2594    // may not be exactly at bound and bounds may have changed
2595    // Make sure outgoing looks feasible
2596    if (!isSuperBasic)
2597      directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
2598    solution_[sequenceOut_] = valueOut_;
2599  }
2600  // change cost and bounds on incoming if primal
2601  nonLinearCost_->setOne(sequenceIn_, valueIn_);
2602  int whatNext = housekeeping(objectiveChange);
2603  if (keepValue)
2604    solution_[sequenceOut_] = saveValue;
2605  if (isSuperBasic)
2606    setStatus(sequenceOut_, superBasic);
2607    //#define CLP_DEBUG
[393]2608#if CLP_DEBUG > 1
[2385]2609  {
2610    int ninf = matrix_->checkFeasible(this);
2611    if (ninf)
2612      printf("infeas %d\n", ninf);
2613  }
[393]2614#endif
[2385]2615  if (whatNext == 1) {
2616    returnCode = -2; // refactorize
2617  } else if (whatNext == 2) {
2618    // maximum iterations or equivalent
2619    returnCode = 3;
2620  } else if (numberIterations_ == lastGoodIteration_ + 2 * factorization_->maximumPivots()) {
2621    // done a lot of flips - be safe
2622    returnCode = -2; // refactorize
2623  }
2624  // Check event
2625  {
2626    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2627    if (status >= 0) {
2628      problemStatus_ = 5;
2629      secondaryStatus_ = ClpEventHandler::endOfIteration;
2630      returnCode = 4;
2631    }
2632  }
2633  return returnCode;
[393]2634}
[2025]2635// May use a cut approach for solving any LP
[2385]2636int ClpSimplexNonlinear::primalDualCuts(char *rowsIn, int startUp, int algorithm)
[2025]2637{
2638  if (!rowsIn) {
[2385]2639    if (algorithm > 0)
[2025]2640      return ClpSimplex::primal(startUp);
2641    else
2642      //return static_cast<ClpSimplexDual *>(static_cast<ClpSimplex *>(this))->dual(startUp);
2643      return ClpSimplex::dual(startUp);
2644  } else {
[2385]2645    int numberUsed = 0;
2646    int rowsThreshold = CoinMax(100, numberRows_ / 2);
[2025]2647    //int rowsTry=CoinMax(50,numberRows_/4);
2648    // Just add this number of rows each time in small problem
2649    int smallNumberRows = 2 * numberColumns_;
2650    smallNumberRows = CoinMin(smallNumberRows, numberRows_ / 20);
2651    // We will need arrays to choose rows to add
[2385]2652    double *weight = new double[numberRows_];
2653    int *sort = new int[numberRows_ + numberColumns_];
2654    int *whichColumns = sort + numberRows_;
[2025]2655    int numberSort = 0;
[2385]2656    for (int i = 0; i < numberRows_; i++) {
[2025]2657      if (rowsIn[i])
[2385]2658        numberUsed++;
[2025]2659    }
2660    if (numberUsed) {
2661      // normal
2662    } else {
2663      // If non slack use that information
[2385]2664      int numberBinding = 0;
2665      numberPrimalInfeasibilities_ = 0;
2666      sumPrimalInfeasibilities_ = 0.0;
[2025]2667      memset(rowActivity_, 0, numberRows_ * sizeof(double));
2668      times(1.0, columnActivity_, rowActivity_);
[2385]2669      for (int i = 0; i < numberRows_; i++) {
2670        double lowerDifference = rowActivity_[i] - rowLower_[i];
2671        double upperDifference = rowActivity_[i] - rowUpper_[i];
2672        if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
2673          numberPrimalInfeasibilities_++;
2674          if (lowerDifference < 0.0)
2675            sumPrimalInfeasibilities_ -= lowerDifference;
2676          else
2677            sumPrimalInfeasibilities_ += upperDifference;
2678          rowsIn[i] = 1;
2679        } else if (getRowStatus(i) != basic) {
2680          numberBinding++;
2681          rowsIn[i] = 1;
2682        }
[2025]2683      }
[2385]2684      if (numberBinding < rowsThreshold) {
2685        // try
[2025]2686      } else {
[2385]2687        // random?? - or give up
2688        // Set up initial list
2689        numberSort = 0;
2690        for (int iRow = 0; iRow < numberRows_; iRow++) {
[2025]2691          weight[iRow] = 1.123e50;
2692          if (rowLower_[iRow] == rowUpper_[iRow]) {
[2385]2693            sort[numberSort++] = iRow;
2694            weight[iRow] = 0.0;
[2025]2695          }
[2385]2696        }
2697        numberSort /= 2;
2698        // and pad out with random rows
2699        double ratio = ((double)(smallNumberRows - numberSort)) / ((double)numberRows_);
2700        for (int iRow = 0; iRow < numberRows_; iRow++) {
[2025]2701          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
[2385]2702            sort[numberSort++] = iRow;
2703        }
2704        // sort
2705        CoinSort_2(weight, weight + numberRows_, sort);
2706        numberSort = CoinMin(numberRows_, smallNumberRows);
2707        memset(rowsIn, 0, numberRows_);
2708        for (int iRow = 0; iRow < numberSort; iRow++)
2709          rowsIn[sort[iRow]] = 1;
[2025]2710      }
2711    }
2712    // see if feasible
[2385]2713    numberPrimalInfeasibilities_ = 0;
[2025]2714    memset(rowActivity_, 0, numberRows_ * sizeof(double));
2715    times(1.0, columnActivity_, rowActivity_);
[2385]2716    for (int i = 0; i < numberRows_; i++) {
2717      double lowerDifference = rowActivity_[i] - rowLower_[i];
2718      double upperDifference = rowActivity_[i] - rowUpper_[i];
2719      if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
2720        if (lowerDifference < 0.0)
2721          sumPrimalInfeasibilities_ -= lowerDifference;
2722        else
2723          sumPrimalInfeasibilities_ += upperDifference;
2724        numberPrimalInfeasibilities_++;
[2025]2725      }
2726    }
2727    printf("Initial infeasibilities - %g (%d)\n",
[2385]2728      sumPrimalInfeasibilities_, numberPrimalInfeasibilities_);
[2025]2729    // Just do this number of passes
2730    int maxPass = 50;
2731    // And take out slack rows until this pass
2732    int takeOutPass = 30;
2733    int iPass;
[2385]2734
2735    const CoinBigIndex *start = this->clpMatrix()->getVectorStarts();
2736    const int *length = this->clpMatrix()->getVectorLengths();
2737    const int *row = this->clpMatrix()->getIndices();
2738    problemStatus_ = 1;
[2025]2739    for (iPass = 0; iPass < maxPass; iPass++) {
2740      printf("Start of pass %d\n", iPass);
2741      int numberSmallColumns = 0;
2742      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
[2385]2743        int n = 0;
2744        for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
2745          int iRow = row[j];
2746          if (rowsIn[iRow])
2747            n++;
2748        }
2749        if (n)
2750          whichColumns[numberSmallColumns++] = iColumn;
[2025]2751      }
[2385]2752      numberSort = 0;
2753      for (int i = 0; i < numberRows_; i++) {
2754        if (rowsIn[i])
2755          sort[numberSort++] = i;
[2025]2756      }
2757      // Create small problem
2758      ClpSimplex small(this, numberSort, sort, numberSmallColumns, whichColumns);
2759      printf("Small model has %d rows, %d columns and %d elements\n",
[2385]2760        small.numberRows(), small.numberColumns(), small.getNumElements());
[2025]2761      small.setFactorizationFrequency(100 + numberSort / 200);
2762      // Solve
[2385]2763      small.setLogLevel(CoinMax(0, logLevel() - 1));
2764      if (iPass > 20) {
2765        if (sumPrimalInfeasibilities_ > 1.0e-1) {
2766          small.dual();
2767        } else {
2768          small.primal(1);
2769        }
[2025]2770      } else {
[2385]2771        // presolve
[2025]2772#if 0
2773        ClpSolve::SolveType method;
2774        ClpSolve::PresolveType presolveType = ClpSolve::presolveOn;
2775        ClpSolve solveOptions;
2776        solveOptions.setPresolveType(presolveType, 5);
2777        if (sumPrimalInfeasibilities_>1.0e-1)
2778          method = ClpSolve::useDual;
2779        else
2780          method = ClpSolve::usePrimalorSprint;
2781        solveOptions.setSolveType(method);
2782        small.initialSolve(solveOptions);
2783#else
2784#if 1
[2385]2785        ClpPresolve *pinfo = new ClpPresolve();
2786        ClpSimplex *small2 = pinfo->presolvedModel(small, 1.0e-5);
2787        assert(small2);
2788        if (sumPrimalInfeasibilities_ > 1.0e-1) {
2789          small2->dual();
2790        } else {
2791          small2->primal(1);
2792        }
2793        pinfo->postsolve(true);
2794        delete pinfo;
[2025]2795#else
[2385]2796        char *types = new char[small.numberRows() + small.numberColumns()];
2797        memset(types, 0, small.numberRows() + small.numberColumns());
2798        if (sumPrimalInfeasibilities_ > 1.0e-1) {
2799          small.miniSolve(types, types + small.numberRows(),
2800            -1, 0);
2801        } else {
2802          small.miniSolve(types, types + small.numberRows(),
2803            1, 1);
2804        }
2805        delete[] types;
[2025]2806#endif
[2385]2807        if (small.sumPrimalInfeasibilities() > 1.0)
2808          small.primal(1);
[2025]2809#endif
2810      }
2811      bool dualInfeasible = (small.status() == 2);
2812      // move solution back
[2385]2813      const double *smallSolution = small.primalColumnSolution();
[2025]2814      for (int j = 0; j < numberSmallColumns; j++) {
[2385]2815        int iColumn = whichColumns[j];
2816        columnActivity_[iColumn] = smallSolution[j];
2817        this->setColumnStatus(iColumn, small.getColumnStatus(j));
[2025]2818      }
2819      for (int iRow = 0; iRow < numberSort; iRow++) {
[2385]2820        int kRow = sort[iRow];
2821        this->setRowStatus(kRow, small.getRowStatus(iRow));
[2025]2822      }
2823      // compute full solution
2824      memset(rowActivity_, 0, numberRows_ * sizeof(double));
2825      times(1.0, columnActivity_, rowActivity_);
2826      if (iPass != maxPass - 1) {
[2385]2827        // Mark row as not looked at
2828        for (int iRow = 0; iRow < numberRows_; iRow++)
2829          weight[iRow] = 1.123e50;
2830        // Look at rows already in small problem
2831        int iSort;
2832        int numberDropped = 0;
2833        int numberKept = 0;
2834        int numberBinding = 0;
2835        numberPrimalInfeasibilities_ = 0;
2836        sumPrimalInfeasibilities_ = 0.0;
2837        bool allFeasible = small.numberIterations() == 0;
2838        for (iSort = 0; iSort < numberSort; iSort++) {
2839          int iRow = sort[iSort];
2840          //printf("%d %g %g\n",iRow,rowActivity_[iRow],small.primalRowSolution()[iSort]);
2841          if (getRowStatus(iRow) == ClpSimplex::basic) {
2842            // Basic - we can get rid of if early on
2843            if (iPass < takeOutPass && !dualInfeasible) {
2844              double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
2845                rowLower_[iRow] - rowActivity_[iRow]);
2846              weight[iRow] = -infeasibility;
2847              if (infeasibility > primalTolerance_ && !allFeasible) {
2848                numberPrimalInfeasibilities_++;
2849                sumPrimalInfeasibilities_ += infeasibility;
2850              } else {
2851                weight[iRow] = 1.0;
2852                numberDropped++;
2853              }
2854            } else {
2855              // keep
2856              weight[iRow] = -1.0e40;
2857              numberKept++;
2858            }
2859          } else {
2860            // keep
2861            weight[iRow] = -1.0e50;
2862            numberKept++;
2863            numberBinding++;
2864          }
2865        }
2866        // Now rest
2867        for (int iRow = 0; iRow < numberRows_; iRow++) {
2868          sort[iRow] = iRow;
2869          if (weight[iRow] == 1.123e50) {
2870            // not looked at yet
2871            double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
2872              rowLower_[iRow] - rowActivity_[iRow]);
2873            weight[iRow] = -infeasibility;
2874            if (infeasibility > primalTolerance_) {
2875              numberPrimalInfeasibilities_++;
2876              sumPrimalInfeasibilities_ += infeasibility;
2877            }
2878          }
2879        }
2880        // sort
2881        CoinSort_2(weight, weight + numberRows_, sort);
2882        numberSort = CoinMin(numberRows_, smallNumberRows + numberKept);
2883        memset(rowsIn, 0, numberRows_);
2884        for (int iRow = 0; iRow < numberSort; iRow++)
2885          rowsIn[sort[iRow]] = 1;
2886        printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
2887          numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns);
2888        printf("%d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
2889          sumPrimalInfeasibilities_);
2890        if (!numberPrimalInfeasibilities_) {
2891          problemStatus_ = 0;
2892          printf("Exiting as looks optimal\n");
2893          break;
2894        }
2895        numberPrimalInfeasibilities_ = 0;
2896        sumPrimalInfeasibilities_ = 0.0;
2897        for (iSort = 0; iSort < numberSort; iSort++) {
2898          if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) {
2899            numberPrimalInfeasibilities_++;
2900            sumPrimalInfeasibilities_ += -weight[iSort];
2901          }
2902        }
2903        printf("in small model %d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
2904          sumPrimalInfeasibilities_);
[2025]2905      } else {
[2385]2906        // out of passes
2907        problemStatus_ = -1;
[2025]2908      }
2909    }
[2385]2910    delete[] weight;
2911    delete[] sort;
[2025]2912    return 0;
2913  }
2914}
[393]2915// A sequential LP method
[2385]2916int ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance,
2917  int otherOptions)
[393]2918{
[2385]2919  // Are we minimizing or maximizing
2920  double whichWay = optimizationDirection();
2921  if (whichWay < 0.0)
2922    whichWay = -1.0;
2923  else if (whichWay > 0.0)
2924    whichWay = 1.0;
[393]2925
[2385]2926  int numberColumns = this->numberColumns();
2927  int numberRows = this->numberRows();
2928  double *columnLower = this->columnLower();
2929  double *columnUpper = this->columnUpper();
2930  double *solution = this->primalColumnSolution();
[393]2931
[2385]2932  if (objective_->type() < 2) {
2933    // no nonlinear part
2934    return ClpSimplex::primal(0);
2935  }
2936  // Get list of non linear columns
2937  char *markNonlinear = new char[numberColumns];
2938  memset(markNonlinear, 0, numberColumns);
2939  int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
[1525]2940
[2385]2941  if (!numberNonLinearColumns) {
2942    delete[] markNonlinear;
2943    // no nonlinear part
2944    return ClpSimplex::primal(0);
2945  }
2946  int iColumn;
2947  int *listNonLinearColumn = new int[numberNonLinearColumns];
2948  numberNonLinearColumns = 0;
2949  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2950    if (markNonlinear[iColumn])
2951      listNonLinearColumn[numberNonLinearColumns++] = iColumn;
2952  }
2953  // Replace objective
2954  ClpObjective *trueObjective = objective_;
2955  objective_ = new ClpLinearObjective(NULL, numberColumns);
2956  double *objective = this->objective();
2957  // See if user wants to use cuts
2958  char *rowsIn = NULL;
2959  if ((otherOptions & 1) != 0 || numberPasses < 0) {
2960    numberPasses = abs(numberPasses);
2961    rowsIn = new char[numberRows_];
2962    memset(rowsIn, 0, numberRows_);
2963  }
2964  // get feasible
2965  if (this->status() < 0 || numberPrimalInfeasibilities())
2966    primalDualCuts(rowsIn, 1, 1);
2967  // still infeasible
2968  if (problemStatus_) {
2969    delete[] listNonLinearColumn;
2970    return 0;
2971  }
2972  algorithm_ = 1;
2973  int jNon;
2974  int *last[3];
[1525]2975
[2385]2976  double *trust = new double[numberNonLinearColumns];
2977  double *trueLower = new double[numberNonLinearColumns];
2978  double *trueUpper = new double[numberNonLinearColumns];
2979  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
2980    iColumn = listNonLinearColumn[jNon];
2981    trust[jNon] = 0.5;
2982    trueLower[jNon] = columnLower[iColumn];
2983    trueUpper[jNon] = columnUpper[iColumn];
2984    if (solution[iColumn] < trueLower[jNon])
2985      solution[iColumn] = trueLower[jNon];
2986    else if (solution[iColumn] > trueUpper[jNon])
2987      solution[iColumn] = trueUpper[jNon];
2988  }
2989  int saveLogLevel = logLevel();
2990  int iPass;
2991  double lastObjective = 1.0e31;
2992  double *saveSolution = new double[numberColumns];
2993  double *saveRowSolution = new double[numberRows];
2994  memset(saveRowSolution, 0, numberRows * sizeof(double));
2995  double *savePi = new double[numberRows];
2996  double *safeSolution = new double[numberColumns];
2997  unsigned char *saveStatus = new unsigned char[numberRows + numberColumns];
[393]2998#define MULTIPLE 0
[2385]2999#if MULTIPLE > 2
3000  // Duplication but doesn't really matter
[1525]3001     double * saveSolutionM[MULTIPLE
[2385]3002};
3003for (jNon = 0; jNon < MULTIPLE; jNon++) {
3004  saveSolutionM[jNon] = new double[numberColumns];
3005  CoinMemcpyN(solution, numberColumns, saveSolutionM);
[1525]3006}
[393]3007#endif
[1525]3008double targetDrop = 1.0e31;
3009double objectiveOffset;
3010getDblParam(ClpObjOffset, objectiveOffset);
3011// 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
[2385]3012for (iPass = 0; iPass < 3; iPass++) {
3013  last[iPass] = new int[numberNonLinearColumns];
3014  for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
3015    last[iPass][jNon] = 0;
[1525]3016}
3017// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3018int goodMove = -2;
[2385]3019char *statusCheck = new char[numberColumns];
3020double *changeRegion = new double[numberColumns];
[1525]3021double offset = 0.0;
3022double objValue = 0.0;
3023int exitPass = 2 * numberPasses + 10;
[2385]3024for (iPass = 0; iPass < numberPasses; iPass++) {
3025  exitPass--;
3026  // redo objective
3027  offset = 0.0;
3028  objValue = -objectiveOffset;
3029  // make sure x updated
3030  trueObjective->newXValues();
3031  double theta = -1.0;
3032  double maxTheta = COIN_DBL_MAX;
3033  //maxTheta=1.0;
3034  if (iPass) {
3035    int jNon = 0;
3036    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3037      changeRegion[iColumn] = solution[iColumn] - saveSolution[iColumn];
3038      double alpha = changeRegion[iColumn];
3039      double oldValue = saveSolution[iColumn];
3040      if (markNonlinear[iColumn] == 0) {
3041        // linear
3042        if (alpha < -1.0e-15) {
3043          // variable going towards lower bound
3044          double bound = columnLower[iColumn];
3045          oldValue -= bound;
3046          if (oldValue + maxTheta * alpha < 0.0) {
3047            maxTheta = CoinMax(0.0, oldValue / (-alpha));
[1525]3048          }
[2385]3049        } else if (alpha > 1.0e-15) {
3050          // variable going towards upper bound
3051          double bound = columnUpper[iColumn];
3052          oldValue = bound - oldValue;
3053          if (oldValue - maxTheta * alpha < 0.0) {
3054            maxTheta = CoinMax(0.0, oldValue / alpha);
[1525]3055          }
[2385]3056        }
3057      } else {
3058        // nonlinear
3059        if (alpha < -1.0e-15) {
3060          // variable going towards lower bound
3061          double bound = trueLower[jNon];
3062          oldValue -= bound;
3063          if (oldValue + maxTheta * alpha < 0.0) {
3064            maxTheta = CoinMax(0.0, oldValue / (-alpha));
[1525]3065          }
[2385]3066        } else if (alpha > 1.0e-15) {
3067          // variable going towards upper bound
3068          double bound = trueUpper[jNon];
3069          oldValue = bound - oldValue;
3070          if (oldValue - maxTheta * alpha < 0.0) {
3071            maxTheta = CoinMax(0.0, oldValue / alpha);
3072          }
3073        }
3074        jNon++;
3075      }
3076    }
3077    // make sure both accurate
3078    memset(rowActivity_, 0, numberRows_ * sizeof(double));
3079    times(1.0, solution, rowActivity_);
3080    memset(saveRowSolution, 0, numberRows_ * sizeof(double));
3081    times(1.0, saveSolution, saveRowSolution);
3082    for (int iRow = 0; iRow < numberRows; iRow++) {
3083      double alpha = rowActivity_[iRow] - saveRowSolution[iRow];
3084      double oldValue = saveRowSolution[iRow];
3085      if (alpha < -1.0e-15) {
3086        // variable going towards lower bound
3087        double bound = rowLower_[iRow];
3088        oldValue -= bound;
3089        if (oldValue + maxTheta * alpha < 0.0) {
3090          maxTheta = CoinMax(0.0, oldValue / (-alpha));
3091        }
3092      } else if (alpha > 1.0e-15) {
3093        // variable going towards upper bound
3094        double bound = rowUpper_[iRow];
3095        oldValue = bound - oldValue;
3096        if (oldValue - maxTheta * alpha < 0.0) {
3097          maxTheta = CoinMax(0.0, oldValue / alpha);
3098        }
3099      }
3100    }
3101  } else {
3102    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3103      changeRegion[iColumn] = 0.0;
3104      saveSolution[iColumn] = solution[iColumn];
3105    }
3106    CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
3107  }
3108  // get current value anyway
3109  double predictedObj, thetaObj;
3110  double maxTheta2 = 2.0; // to work out a b c
3111  double theta2 = trueObjective->stepLength(this, saveSolution, changeRegion, maxTheta2,
3112    objValue, predictedObj, thetaObj);
3113  int lastMoveStatus = goodMove;
3114  if (goodMove >= 0) {
3115    theta = CoinMin(theta2, maxTheta);
[393]3116#ifdef CLP_DEBUG
[2385]3117    if (handler_->logLevel() & 32)
3118      printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
3119        theta, objValue, thetaObj, predictedObj);
[393]3120#endif
[2385]3121    if (theta > 0.0 && theta <= 1.0) {
3122      // update solution
3123      double lambda = 1.0 - theta;
3124      for (iColumn = 0; iColumn < numberColumns; iColumn++)
3125        solution[iColumn] = lambda * saveSolution[iColumn]
3126          + theta * solution[iColumn];
3127      memset(rowActivity_, 0, numberRows_ * sizeof(double));
3128      times(1.0, solution, rowActivity_);
3129      if (lambda > 0.999) {
3130        CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
3131        CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
3132      }
3133      // Do local minimization
[393]3134#define LOCAL
3135#ifdef LOCAL
[2385]3136      bool absolutelyOptimal = false;
3137      int saveScaling = scalingFlag();
3138      scaling(0);
3139      int savePerturbation = perturbation_;
3140      perturbation_ = 100;
3141      if (saveLogLevel == 1)
3142        setLogLevel(0);
3143      int status = startup(1);
3144      if (!status) {
3145        int numberTotal = numberRows_ + numberColumns_;
3146        // resize arrays
3147        for (int i = 0; i < 4; i++) {
3148          rowArray_[i]->reserve(CoinMax(numberRows_ + numberColumns_, rowArray_[i]->capacity()));
3149        }
3150        CoinIndexedVector *longArray = rowArray_[3];
3151        CoinIndexedVector *rowArray = rowArray_[0];
3152        //CoinIndexedVector * columnArray = columnArray_[0];
3153        CoinIndexedVector *spare = rowArray_[1];
3154        double *work = longArray->denseVector();
3155        int *which = longArray->getIndices();
3156        int nPass = 100;
3157        //bool conjugate=false;
3158        // Put back true bounds
3159        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3160          int iColumn = listNonLinearColumn[jNon];
3161          double value;
3162          value = trueLower[jNon];
3163          trueLower[jNon] = lower_[iColumn];
3164          lower_[iColumn] = value;
3165          value = trueUpper[jNon];
3166          trueUpper[jNon] = upper_[iColumn];
3167          upper_[iColumn] = value;
3168          switch (getStatus(iColumn)) {
[1525]3169
[2385]3170          case basic:
3171          case isFree:
3172          case superBasic:
3173            break;
3174          case isFixed:
3175          case atUpperBound:
3176          case atLowerBound:
3177            if (solution_[iColumn] > lower_[iColumn] + primalTolerance_) {
3178              if (solution_[iColumn] < upper_[iColumn] - primalTolerance_) {
3179                setStatus(iColumn, superBasic);
3180              } else {
3181                setStatus(iColumn, atUpperBound);
3182              }
3183            } else {
3184              if (solution_[iColumn] < upper_[iColumn] - primalTolerance_) {
3185                setStatus(iColumn, atLowerBound);
3186              } else {
3187                setStatus(iColumn, isFixed);
3188              }
3189            }
3190            break;
3191          }
3192        }
3193        for (int jPass = 0; jPass < nPass; jPass++) {
3194          trueObjective->reducedGradient(this, dj_, true);
3195          // get direction vector in longArray
3196          longArray->clear();
3197          double normFlagged = 0.0;
3198          double normUnflagged = 0.0;
3199          int numberNonBasic = 0;
3200          directionVector(longArray, spare, rowArray, 0,
3201            normFlagged, normUnflagged, numberNonBasic);
3202          if (normFlagged + normUnflagged < 1.0e-8) {
3203            absolutelyOptimal = true;
3204            break; //looks optimal
3205          }
3206          double djNorm = 0.0;
3207          int iIndex;
3208          for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
3209            int iSequence = which[iIndex];
3210            double alpha = work[iSequence];
3211            djNorm += alpha * alpha;
3212          }
3213          //if (!jPass)
3214          //printf("dj norm %g - %d \n",djNorm,numberNonBasic);
3215          //int number=longArray->getNumElements();
3216          if (!numberNonBasic) {
3217            // looks optimal
3218            absolutelyOptimal = true;
3219            break;
3220          }
3221          double theta = 1.0e30;
3222          int iSequence;
3223          for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3224            double alpha = work[iSequence];
3225            double oldValue = solution_[iSequence];
3226            if (alpha < -1.0e-15) {
3227              // variable going towards lower bound
3228              double bound = lower_[iSequence];
3229              oldValue -= bound;
3230              if (oldValue + theta * alpha < 0.0) {
3231                theta = CoinMax(0.0, oldValue / (-alpha));
3232              }
3233            } else if (alpha > 1.0e-15) {
3234              // variable going towards upper bound
3235              double bound = upper_[iSequence];
3236              oldValue = bound - oldValue;
3237              if (oldValue - theta * alpha < 0.0) {
3238                theta = CoinMax(0.0, oldValue / alpha);
3239              }
3240            }
3241          }
3242          // Now find minimum of function
3243          double currentObj;
3244          double predictedObj;
3245          double thetaObj;
3246          // need to use true objective
3247          double *saveCost = cost_;
3248          cost_ = NULL;
3249          double objTheta = trueObjective->stepLength(this, solution_, work, theta,
3250            currentObj, predictedObj, thetaObj);
3251          cost_ = saveCost;
[393]3252#ifdef CLP_DEBUG
[2385]3253          if (handler_->logLevel() & 32)
3254            printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
[393]3255#endif
[2385]3256          //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
3257          //printf("objTheta %g theta %g\n",objTheta,theta);
3258          if (theta > objTheta) {
3259            theta = objTheta;
3260            thetaObj = predictedObj;
3261          }
3262          // update one used outside
3263          objValue = currentObj;
3264          if (theta > 1.0e-9 && (currentObj - thetaObj < -CoinMax(1.0e-8, 1.0e-15 * fabs(currentObj)) || jPass < 5)) {
3265            // Update solution
3266            for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3267              double alpha = work[iSequence];
3268              if (alpha) {
3269                double value = solution_[iSequence] + theta * alpha;
3270                solution_[iSequence] = value;
3271                switch (getStatus(iSequence)) {
[1525]3272
[2385]3273                case basic:
3274                case isFixed:
3275                case isFree:
3276                  break;
3277                case atUpperBound:
3278                case atLowerBound:
3279                case superBasic:
3280                  if (value > lower_[iSequence] + primalTolerance_) {
3281                    if (value < upper_[iSequence] - primalTolerance_) {
3282                      setStatus(iSequence, superBasic);
3283                    } else {
3284                      setStatus(iSequence, atUpperBound);
[1525]3285                    }
[2385]3286                  } else {
3287                    if (value < upper_[iSequence] - primalTolerance_) {
3288                      setStatus(iSequence, atLowerBound);
3289                    } else {
3290                      setStatus(iSequence, isFixed);
[1525]3291                    }
[2385]3292                  }
3293                  break;
3294                }
3295              }
3296            }
3297          } else {
3298            break;
3299          }
3300        }
3301        // Put back fake bounds
3302        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3303          int iColumn = listNonLinearColumn[jNon];
3304          double value;
3305          value = trueLower[jNon];
3306          trueLower[jNon] = lower_[iColumn];
3307          lower_[iColumn] = value;
3308          value = trueUpper[jNon];
3309          trueUpper[jNon] = upper_[iColumn];
3310          upper_[iColumn] = value;
3311        }
3312      }
3313      problemStatus_ = 0;
3314      finish();
3315      scaling(saveScaling);
3316      perturbation_ = savePerturbation;
3317      setLogLevel(saveLogLevel);
[393]3318#endif
[2385]3319      // redo rowActivity
3320      memset(rowActivity_, 0, numberRows_ * sizeof(double));
3321      times(1.0, solution, rowActivity_);
3322      if (theta > 0.99999 && theta2 < 1.9 && !absolutelyOptimal) {
3323        // If big changes then tighten
3324        /*  thetaObj is objvalue + a*2*2 +b*2
[1525]3325                        predictedObj is objvalue + a*theta2*theta2 +b*theta2
3326                    */
[2385]3327        double rhs1 = thetaObj - objValue;
3328        double rhs2 = predictedObj - objValue;
3329        double subtractB = theta2 * 0.5;
3330        double a = (rhs2 - subtractB * rhs1) / (theta2 * theta2 - 4.0 * subtractB);
3331        double b = 0.5 * (rhs1 - 4.0 * a);
3332        if (fabs(a + b) > 1.0e-2) {
3333          // tighten all
3334          goodMove = -1;
3335        }
3336      }
3337    }
3338  }
3339  CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns,
3340    objective);
3341  //printf("offset comp %g orig %g\n",offset,objectiveOffset);
3342  // could do faster
3343  trueObjective->stepLength(this, solution, changeRegion, 0.0,
3344    objValue, predictedObj, thetaObj);
[1271]3345#ifdef CLP_INVESTIGATE
[2385]3346  printf("offset comp %g orig %g - obj from stepLength %g\n", offset, objectiveOffset, objValue);
[1271]3347#endif
[2385]3348  setDblParam(ClpObjOffset, objectiveOffset + offset);
3349  int *temp = last[2];
3350  last[2] = last[1];
3351  last[1] = last[0];
3352  last[0] = temp;
3353  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3354    iColumn = listNonLinearColumn[jNon];
3355    double change = solution[iColumn] - saveSolution[iColumn];
3356    if (change < -1.0e-5) {
3357      if (fabs(change + trust[jNon]) < 1.0e-5)
3358        temp[jNon] = -1;
3359      else
3360        temp[jNon] = -2;
3361    } else if (change > 1.0e-5) {
3362      if (fabs(change - trust[jNon]) < 1.0e-5)
3363        temp[jNon] = 1;
3364      else
3365        temp[jNon] = 2;
3366    } else {
3367      temp[jNon] = 0;
3368    }
3369  }
3370  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3371  double maxDelta = 0.0;
3372  if (goodMove >= 0) {
3373    if (objValue - lastObjective <= 1.0e-15 * fabs(lastObjective))
3374      goodMove = 1;
3375    else
3376      goodMove = 0;
3377  } else {
3378    maxDelta = 1.0e10;
3379  }
3380  double maxGap = 0.0;
3381  int numberSmaller = 0;
3382  int numberSmaller2 = 0;
3383  int numberLarger = 0;
3384  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3385    iColumn = listNonLinearColumn[jNon];
3386    maxDelta = CoinMax(maxDelta,
3387      fabs(solution[iColumn] - saveSolution[iColumn]));
3388    if (goodMove > 0) {
3389      if (last[0][jNon] * last[1][jNon] < 0) {
3390        // halve
3391        trust[jNon] *= 0.5;
3392        numberSmaller2++;
3393      } else {
3394        if (last[0][jNon] == last[1][jNon] && last[0][jNon] == last[2][jNon])
3395          trust[jNon] = CoinMin(1.5 * trust[jNon], 1.0e6);
3396        numberLarger++;
3397      }
3398    } else if (goodMove != -2 && trust[jNon] > 10.0 * deltaTolerance) {
3399      trust[jNon] *= 0.2;
3400      numberSmaller++;
3401    }
3402    maxGap = CoinMax(maxGap, trust[jNon]);
3403  }
[393]3404#ifdef CLP_DEBUG
[2385]3405  if (handler_->logLevel() & 32)
3406    std::cout << "largest gap is " << maxGap << " "
3407              << numberSmaller + numberSmaller2 << " reduced ("
3408              << numberSmaller << " badMove ), "
3409              << numberLarger << " increased" << std::endl;
[393]3410#endif
[2385]3411  if (iPass > 10000) {
3412    for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
3413      trust[jNon] *= 0.0001;
3414  }
3415  if (lastMoveStatus == -1 && goodMove == -1)
3416    goodMove = 1; // to force solve
3417  if (goodMove > 0) {
3418    double drop = lastObjective - objValue;
3419    handler_->message(CLP_SLP_ITER, messages_)
3420      << iPass << objValue - objectiveOffset
3421      << drop << maxDelta
3422      << CoinMessageEol;
3423    if (iPass > 20 && drop < 1.0e-12 * fabs(objValue))
3424      drop = 0.999e-4; // so will exit
3425    if (maxDelta < deltaTolerance && drop < 1.0e-4 && goodMove && theta < 0.99999) {
3426      if (handler_->logLevel() > 1)
3427        std::cout << "Exiting as maxDelta < tolerance and small drop" << std::endl;
3428      break;
3429    }
3430  } else if (!numberSmaller && iPass > 1) {
3431    if (handler_->logLevel() > 1)
3432      std::cout << "Exiting as all gaps small" << std::endl;
3433    break;
3434  }
3435  if (!iPass)
3436    goodMove = 1;
3437  targetDrop = 0.0;
3438  double *r = this->dualColumnSolution();
3439  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3440    iColumn = listNonLinearColumn[jNon];
3441    columnLower[iColumn] = CoinMax(solution[iColumn]
3442        - trust[jNon],
3443      trueLower[jNon]);
3444    columnUpper[iColumn] = CoinMin(solution[iColumn]
3445        + trust[jNon],
3446      trueUpper[jNon]);
3447  }
3448  if (iPass) {
3449    // get reduced costs
3450    this->matrix