source: trunk/Clp/src/ClpSimplexNonlinear.cpp @ 2445

Last change on this file since 2445 was 2445, checked in by forrest, 7 months ago

hopefully fix bad read for quad.mps

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 152.3 KB
Line 
1/* $Id: ClpSimplexNonlinear.cpp 2445 2019-04-06 15:56:47Z forrest $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9#include "CoinHelperFunctions.hpp"
10#include "ClpHelperFunctions.hpp"
11#include "ClpSimplexNonlinear.hpp"
12#include "ClpSimplexOther.hpp"
13#include "ClpSimplexDual.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpNonLinearCost.hpp"
16#include "ClpLinearObjective.hpp"
17#include "ClpConstraint.hpp"
18#include "ClpPresolve.hpp"
19#include "ClpQuadraticObjective.hpp"
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
33// primal
34int ClpSimplexNonlinear::primal()
35{
36
37  int ifValuesPass = 1;
38  algorithm_ = +3;
39
40  // save data
41  ClpDataSave data = saveData();
42  matrix_->refresh(this); // make sure matrix okay
43
44  // Save objective
45  ClpObjective *saveObjective = NULL;
46  if (objective_->type() > 1) {
47    // expand to full if quadratic
48#ifndef NO_RTTI
49    ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
50#else
51    ClpQuadraticObjective *quadraticObj = NULL;
52    if (objective_->type() == 2)
53      quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
54#endif
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;
65
66  // initialize - maybe values pass and algorithm_ is +1
67  // true does something (? perturbs)
68  if (!startup(true)) {
69
70    // Set average theta
71    nonLinearCost_->setAverageTheta(1.0e3);
72    int lastCleaned = 0; // last time objective or bounds cleaned up
73
74    // Say no pivot has occurred (for steepest edge and updates)
75    pivotRow_ = -2;
76
77    // This says whether to restore things etc
78    int factorType = 0;
79    // Start check for cycles
80    progress_.startCheck();
81    /*
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          */
92    while (problemStatus_ < 0) {
93      int iRow, iColumn;
94      // clear
95      for (iRow = 0; iRow < 4; iRow++) {
96        rowArray_[iRow]->clear();
97      }
98
99      for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
100        columnArray_[iColumn]->clear();
101      }
102
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;
110
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--;
117#ifdef CLP_DEBUG
118          if (handler_->logLevel() & 32)
119            printf("pivot mode now %d\n", pivotMode);
120#endif
121          if (pivotMode == 9)
122            pivotMode = 0; // switch off fast attempt
123        }
124      }
125      statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true,
126        bestObjectiveWhenFlagged);
127
128      // Say good factorization
129      factorType = 1;
130
131      // Say no pivot has occurred (for steepest edge and updates)
132      pivotRow_ = -2;
133
134      // exit if victory declared
135      if (problemStatus_ >= 0)
136        break;
137
138      // test for maximum iterations
139      if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
140        problemStatus_ = 3;
141        break;
142      }
143
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;
153          }
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_;
195}
196/*  Refactorizes if necessary
197    Checks if finished.  Updates status.
198    lastCleaned refers to iteration at which some objective/feasibility
199    cleaning too place.
200
201    type - 0 initial so set up save arrays etc
202    - 1 normal -if good update save
203    - 2 restoring from saved
204*/
205void ClpSimplexNonlinear::statusOfProblemInPrimal(int &lastCleaned, int type,
206  ClpSimplexProgress *progress,
207  bool doFactorization,
208  double &bestObjectiveWhenFlagged)
209{
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);
234
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);
253          } else {
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
267          }
268          // flag incoming
269          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
270            setFlagged(sequenceIn_);
271            saveStatus_[sequenceIn_] = status_[sequenceIn_];
272          }
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
339
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);
395
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
417                         if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
418                              goToDual = false;
419                         */
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_++;
475#ifdef CLP_DEBUG
476            if (handler_->logLevel() & 32)
477              printf("%d times optimal, current %g best %g\n", numberTimesOptimal_,
478                current, bestObjectiveWhenFlagged);
479#endif
480          } else {
481            bestObjectiveWhenFlagged = current;
482          }
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));
487          }
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;
507          else
508            problemStatus_ = -1;
509        } else {
510          problemStatus_ = 0; // optimal
511          if (lastCleaned < numberIterations_) {
512            handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
513              << CoinMessageEol;
514          }
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);
592}
593/*
594  Reasons to come out:
595  -1 iterations etc
596  -2 inaccuracy
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
601  +3 max iterations
602*/
603int ClpSimplexNonlinear::whileIterating(int &pivotMode)
604{
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
621#if CLP_DEBUG > 1
622    rowArray_[0]->checkClear();
623    rowArray_[1]->checkClear();
624    rowArray_[2]->checkClear();
625    rowArray_[3]->checkClear();
626    columnArray_[0]->checkClear();
627#endif
628    if (numberInterior >= 5) {
629      // this can go when we have better minimization
630      if (pivotMode < 10)
631        pivotMode = 1;
632      unflag();
633#ifdef CLP_DEBUG
634      if (handler_->logLevel() & 32)
635        printf("interior unflag\n");
636#endif
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();
645#ifdef CLP_DEBUG
646        if (handler_->logLevel() & 32)
647          printf("unflagging as interior\n");
648#endif
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) {
664#ifdef CLP_INVESTIGATE
665          printf("looks bad - no change in obj %g\n", currentObj);
666#endif
667          if (factorization_->pivots())
668            result = 3;
669          else
670            problemStatus_ = 0;
671        }
672      }
673      if (result == 3)
674        break; // null vector not accurate
675#ifdef CLP_DEBUG
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      }
684#endif
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_));
707#ifdef CLP_DEBUG
708      if ((handler_->logLevel() & 32)) {
709        char x = isColumn(sequenceIn_) ? 'C' : 'R';
710        std::cout << "pivot column " << x << sequenceWithin(sequenceIn_) << std::endl;
711      }
712#endif
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
745#ifdef CLP_DEBUG
746      if (handler_->logLevel() & 32)
747        printf("** no column pivot\n");
748#endif
749      if (pivotMode < 10) {
750        // looks optimal
751        primalColumnPivot_->setLooksOptimal(true);
752      } else {
753        pivotMode--;
754#ifdef CLP_DEBUG
755        if (handler_->logLevel() & 32)
756          printf("pivot mode now %d\n", pivotMode);
757#endif
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;
770}
771// Creates direction vector
772void ClpSimplexNonlinear::directionVector(CoinIndexedVector *vectorArray,
773  CoinIndexedVector *spare1, CoinIndexedVector *spare2,
774  int pivotMode2,
775  double &normFlagged, double &normUnflagged,
776  int &numberNonBasic)
777{
778#if CLP_DEBUG > 1
779  vectorArray->checkClear();
780  spare1->checkClear();
781  spare2->checkClear();
782#endif
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)) {
806
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)) {
827
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
864#define NOSUPER
865#ifndef ALLSUPER
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          }
877#else
878          array[iSequence] = -dj_[iSequence];
879          index[number++] = iSequence;
880          if (pivotMode2 >= 10)
881            bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
882#endif
883          break;
884        }
885      }
886#ifdef NOSUPER
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      }
902#else
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      }
918#endif
919#ifdef CLP_DEBUG
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      }
927#endif
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_)) {
935
936        case basic:
937          sequenceOut_ = -1;
938        case ClpSimplex::isFixed:
939          break;
940        case atUpperBound:
941          if (dj_[sequenceOut_] > dualTolerance_) {
942#ifdef CLP_DEBUG
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_]);
947#endif
948          }
949          break;
950        case atLowerBound:
951          if (dj_[sequenceOut_] < -dualTolerance_) {
952#ifdef CLP_DEBUG
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_]);
957#endif
958          }
959          break;
960        case isFree:
961        case superBasic:
962          if (dj_[sequenceOut_] > dualTolerance_) {
963#ifdef CLP_DEBUG
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_]);
968#endif
969          } else if (dj_[sequenceOut_] < -dualTolerance_) {
970#ifdef CLP_DEBUG
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_]);
975#endif
976          }
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)) {
987
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            }
1003          }
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            }
1017          }
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            }
1047          }
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);
1194}
1195#define MINTYPE 1
1196#if MINTYPE == 2
1197static double
1198innerProductIndexed(const double *region1, int size, const double *region2, const int *which)
1199{
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;
1207}
1208#endif
1209/*
1210   Row array and column array have direction
1211   Returns 0 - can do normal iteration (basis change)
1212   1 - no basis change
1213*/
1214int ClpSimplexNonlinear::pivotColumn(CoinIndexedVector *longArray,
1215  CoinIndexedVector *rowArray,
1216  CoinIndexedVector *columnArray,
1217  CoinIndexedVector *spare,
1218  int &pivotMode,
1219  double &solutionError,
1220  double *dArray)
1221{
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;
1232
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];
1272#endif
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
1282#define SMALLTHETA1 1.0e-10
1283#define SMALLTHETA2 1.0e-10
1284#define CONJUGATE 2
1285#if CONJUGATE == 0
1286    int conjugate = 0;
1287#elif CONJUGATE == 1
1288    int conjugate = (pivotMode < 10) ? MINTYPE : 0;
1289#elif CONJUGATE == 2
1290    int conjugate = MINTYPE;
1291#else
1292    int conjugate = MINTYPE;
1293#endif
1294
1295    //if (pivotMode==1)
1296    //localPivotMode=11;;
1297#if CLP_DEBUG > 1
1298    longArray->checkClear();
1299#endif
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();
1339#if CLP_DEBUG > 1
1340        spare->checkClear();
1341#endif
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;
1346#if CLP_DEBUG > 0
1347        int iLargest = -1;
1348#endif
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;
1354#if CLP_DEBUG > 0
1355            iLargest = iRow;
1356#endif
1357          }
1358        }
1359#if CLP_DEBUG > 0
1360        if ((handler_->logLevel() & 32) && largest > 1.0e-8)
1361          printf("largest rhs error %g on row %d\n", largest, iLargest);
1362#endif
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;
1376#endif
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      }
1390#if CONJUGATE == 2
1391      conjugate = (localPivotMode < 10) ? MINTYPE : 0;
1392#endif
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        }
1405#endif
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();
1414
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;
1430#ifdef CLP_DEBUG
1431          int kPivot = -1;
1432#endif
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) {
1440#ifdef CLP_DEBUG
1441                kPivot = iPivot;
1442#endif
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) {
1452#ifdef CLP_DEBUG
1453                kPivot = iPivot;
1454#endif
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;
1476#ifdef CLP_DEBUG
1477              if (handler_->logLevel() & 32)
1478                printf("easy - obj drop %g, easy drop %g\n", objDrop, easyDrop);
1479#endif
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();
1491#ifdef CLP_DEBUG
1492          if (handler_->logLevel() & 32)
1493            printf("largest %g piv %d\n", largest, kPivot);
1494#endif
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        }
1541#else
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);
1593#if 1
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          }
1624#endif
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;
1654#if CLP_DEBUG > 1
1655            spare1->checkClear();
1656            spare2->checkClear();
1657#endif
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);
1694#if 1
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          }
1717#endif
1718          CoinMemcpyN(saveY2, numberNonBasic, saveY);
1719          CoinMemcpyN(saveS2, numberNonBasic, saveS);
1720        }
1721#endif
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      }
1729#endif
1730      if (djNorm < eps * djNorm0 || (nPasses > 100 && djNorm < CoinMin(1.0e-1 * djNorm0, 1.0e-12))) {
1731#ifdef CLP_DEBUG
1732        if (handler_->logLevel() & 32)
1733          printf("dj norm reduced from %g to %g\n", djNorm0, djNorm);
1734#endif
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) {
1748#ifdef CLP_DEBUG
1749        if (handler_->logLevel() & 32)
1750          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n", djNorm0, djNorm,
1751            normFlagged);
1752#endif
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;
1761#ifdef CLP_DEBUG
1762        if (handler_->logLevel() & 32)
1763          printf("dj norm NOT reduced from %g to %g\n", djNorm0, djNorm);
1764#endif
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)) {
1777
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)) {
1798
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        }
1817#ifdef CLP_DEBUG
1818        if (handler_->logLevel() & 32)
1819          printf("%d super, %d normal, %d flagged\n",
1820            nSuper, nNormal, nFlagged);
1821#endif
1822
1823        int nFlagged2 = 1;
1824        if (lastSequenceIn < 0 && !nNormal && !nSuper) {
1825          nFlagged2 = unflag();
1826          if (pivotMode >= 10) {
1827            pivotMode--;
1828#ifdef CLP_DEBUG
1829            if (handler_->logLevel() & 32)
1830              printf("pivot mode now %d\n", pivotMode);
1831#endif
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;
1851
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);
1912#ifdef CLP_DEBUG
1913      if (handler_->logLevel() & 32)
1914        printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
1915#endif
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
1930#ifdef INCLUDESLACK
1931        for (; iSequence < numberColumns_ + numberRows_; iSequence++) {
1932          double alpha = dArray[iSequence];
1933          double value = alpha * cost_[iSequence];
1934          product += value;
1935        }
1936#endif
1937        if (product > 0.0)
1938          objTheta = djNorm / product;
1939        else
1940          objTheta = COIN_DBL_MAX;
1941#ifndef NDEBUG
1942        if (product < -1.0e-8 && handler_->logLevel() > 1)
1943          printf("bad product %g\n", product);
1944#endif
1945        product = CoinMax(product, 0.0);
1946      } else {
1947        objTheta = objTheta2;
1948      }
1949#else
1950      objTheta = objTheta2;
1951#endif
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
1974#ifdef CLP_DEBUG
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      }
1990#endif
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++) {
2018
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)) {
2025
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      }
2060#ifdef CLP_DEBUG
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      }
2066#endif
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        }
2075#ifdef CLP_DEBUG
2076        if (handler_->logLevel() & 32)
2077          printf("some progress?\n");
2078#endif
2079      }
2080#if CLP_DEBUG > 1
2081      longArray->checkClean();
2082#endif
2083    }
2084#ifdef CLP_DEBUG
2085    if (handler_->logLevel() & 32)
2086      printf("out of loop after %d (%d) passes\n", nPasses, nTotalPasses);
2087#endif
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    }
2101    if (theta_ >= 1.0e30) { // odd
2102      returnCode=2;
2103      break;
2104    }
2105    // See if we need to pivot
2106    if (theta_ == basicTheta || ordinaryDj) {
2107      if (!ordinaryDj) {
2108        // find an incoming column which will force pivot
2109        int iRow;
2110        pivotRow_ = -1;
2111        for (iRow = 0; iRow < numberRows_; iRow++) {
2112          if (pivotVariable_[iRow] == bestBasicSequence) {
2113            pivotRow_ = iRow;
2114            break;
2115          }
2116        }
2117        assert(pivotRow_ >= 0);
2118        // Get good size for pivot
2119        double acceptablePivot = 1.0e-7;
2120        if (factorization_->pivots() > 10)
2121          acceptablePivot = 1.0e-5; // if we have iterated be more strict
2122        else if (factorization_->pivots() > 5)
2123          acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
2124        // should be dArray but seems better this way!
2125        double direction = work[bestBasicSequence] > 0.0 ? -1.0 : 1.0;
2126        // create as packed
2127        rowArray->createPacked(1, &pivotRow_, &direction);
2128        factorization_->updateColumnTranspose(spare, rowArray);
2129        // put row of tableau in rowArray and columnArray
2130        matrix_->transposeTimes(this, -1.0,
2131          rowArray, spare, columnArray);
2132        // choose one futhest away from bound which has reasonable pivot
2133        // If increasing we want negative alpha
2134
2135        double *work2;
2136        int iSection;
2137
2138        sequenceIn_ = -1;
2139        double bestValue = -1.0;
2140        double bestDirection = 0.0;
2141        // First pass we take correct direction and large pivots
2142        // then correct direction
2143        // then any
2144        double check[] = { 1.0e-8, -1.0e-12, -1.0e30 };
2145        double mult[] = { 100.0, 1.0, 1.0 };
2146        for (int iPass = 0; iPass < 3; iPass++) {
2147          //if (!bestValue&&iPass==2)
2148          //bestValue=-1.0;
2149          double acceptable = acceptablePivot * mult[iPass];
2150          double checkValue = check[iPass];
2151          for (iSection = 0; iSection < 2; iSection++) {
2152
2153            int addSequence;
2154
2155            if (!iSection) {
2156              work2 = rowArray->denseVector();
2157              number = rowArray->getNumElements();
2158              which = rowArray->getIndices();
2159              addSequence = numberColumns_;
2160            } else {
2161              work2 = columnArray->denseVector();
2162              number = columnArray->getNumElements();
2163              which = columnArray->getIndices();
2164              addSequence = 0;
2165            }
2166            int i;
2167
2168            for (i = 0; i < number; i++) {
2169              int iSequence = which[i] + addSequence;
2170              if (flagged(iSequence))
2171                continue;
2172              //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
2173              //                  upper_[iSequence]-solution_[iSequence]);
2174              double alpha = work2[i];
2175              // should be dArray but seems better this way!
2176              double change = work[iSequence];
2177              Status thisStatus = getStatus(iSequence);
2178              double direction = 0;
2179              ;
2180              switch (thisStatus) {
2181
2182              case basic:
2183              case ClpSimplex::isFixed:
2184                break;
2185              case isFree:
2186              case superBasic:
2187                if (alpha < -acceptable && change > checkValue)
2188                  direction = 1.0;
2189                else if (alpha > acceptable && change < -checkValue)
2190                  direction = -1.0;
2191                break;
2192              case atUpperBound:
2193                if (alpha > acceptable && change < -checkValue)
2194                  direction = -1.0;
2195                else if (iPass == 2 && alpha < -acceptable && change < -checkValue)
2196                  direction = 1.0;
2197                break;
2198              case atLowerBound:
2199                if (alpha < -acceptable && change > checkValue)
2200                  direction = 1.0;
2201                else if (iPass == 2 && alpha > acceptable && change > checkValue)
2202                  direction = -1.0;
2203                break;
2204              }
2205              if (direction) {
2206                if (sequenceIn_ != lastSequenceIn || localPivotMode < 10) {
2207                  if (CoinMin(solution_[iSequence] - lower_[iSequence],
2208                        upper_[iSequence] - solution_[iSequence])
2209                    > bestValue) {
2210                    bestValue = CoinMin(solution_[iSequence] - lower_[iSequence],
2211                      upper_[iSequence] - solution_[iSequence]);
2212                    sequenceIn_ = iSequence;
2213                    bestDirection = direction;
2214                  }
2215                } else {
2216                  // choose
2217                  bestValue = COIN_DBL_MAX;
2218                  sequenceIn_ = iSequence;
2219                  bestDirection = direction;
2220                }
2221              }
2222            }
2223          }
2224          if (sequenceIn_ >= 0 && bestValue > 0.0)
2225            break;
2226        }
2227        if (sequenceIn_ >= 0) {
2228          valueIn_ = solution_[sequenceIn_];
2229          dualIn_ = dj_[sequenceIn_];
2230          if (bestDirection < 0.0) {
2231            // we want positive dj
2232            if (dualIn_ <= 0.0) {
2233              //printf("bad dj - xx %g\n",dualIn_);
2234              dualIn_ = 1.0;
2235            }
2236          } else {
2237            // we want negative dj
2238            if (dualIn_ >= 0.0) {
2239              //printf("bad dj - xx %g\n",dualIn_);
2240              dualIn_ = -1.0;
2241            }
2242          }
2243          lowerIn_ = lower_[sequenceIn_];
2244          upperIn_ = upper_[sequenceIn_];
2245          if (dualIn_ > 0.0)
2246            directionIn_ = -1;
2247          else
2248            directionIn_ = 1;
2249        } else {
2250          //ordinaryDj=true;
2251#ifdef CLP_DEBUG
2252          if (handler_->logLevel() & 32) {
2253            printf("no easy pivot - norm %g mode %d\n", djNorm, localPivotMode);
2254            if (rowArray->getNumElements() + columnArray->getNumElements() < 12) {
2255              for (iSection = 0; iSection < 2; iSection++) {
2256
2257                int addSequence;
2258
2259                if (!iSection) {
2260                  work2 = rowArray->denseVector();
2261                  number = rowArray->getNumElements();
2262                  which = rowArray->getIndices();
2263                  addSequence = numberColumns_;
2264                } else {
2265                  work2 = columnArray->denseVector();
2266                  number = columnArray->getNumElements();
2267                  which = columnArray->getIndices();
2268                  addSequence = 0;
2269                }
2270                int i;
2271                char section[] = { 'R', 'C' };
2272                for (i = 0; i < number; i++) {
2273                  int iSequence = which[i] + addSequence;
2274                  if (flagged(iSequence)) {
2275                    printf("%c%d flagged\n", section[iSection], which[i]);
2276                    continue;
2277                  } else {
2278                    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
2279                      section[iSection], which[i], status_[iSequence],
2280                      lower_[iSequence], solution_[iSequence], upper_[iSequence],
2281                      work2[i], work[iSequence]);
2282                  }
2283                }
2284              }
2285            }
2286          }
2287#endif
2288          assert(sequenceIn_ < 0);
2289          justOne = true;
2290          allFinished = false; // Round again
2291          finished = false;
2292          nTotalPasses += nPasses;
2293          nPasses = 0;
2294          if (djNorm < 0.9 * djNorm0 && djNorm < 1.0e-3 && !localPivotMode) {
2295#ifdef CLP_DEBUG
2296            if (handler_->logLevel() & 32)
2297              printf("no pivot - mode %d norms %g %g - unflagging\n",
2298                localPivotMode, djNorm0, djNorm);
2299#endif
2300            unflag(); //unflagging
2301            returnCode = 1;
2302          } else {
2303            returnCode = 2; // do single incoming
2304            returnCode = 1;
2305          }
2306        }
2307      }
2308      rowArray->clear();
2309      columnArray->clear();
2310      longArray->clear();
2311      if (ordinaryDj) {
2312        assert (sequenceIn_>=0);
2313        valueIn_ = solution_[sequenceIn_];
2314        dualIn_ = dj_[sequenceIn_];
2315        lowerIn_ = lower_[sequenceIn_];
2316        upperIn_ = upper_[sequenceIn_];
2317        if (dualIn_ > 0.0)
2318          directionIn_ = -1;
2319        else
2320          directionIn_ = 1;
2321      }
2322      if (returnCode == 1)
2323        returnCode = 0;
2324    } else {
2325      // round again
2326      longArray->clear();
2327      if (djNorm < 1.0e-3 && !localPivotMode) {
2328        if (djNorm == 1.2345e-20 && djNorm0 > 1.0e-4) {
2329#ifdef CLP_DEBUG
2330          if (handler_->logLevel() & 32)
2331            printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n", djNorm0, nPasses,
2332              localPivotMode, returnCode);
2333#endif
2334          //if (!localPivotMode)
2335          //returnCode=2; // force singleton
2336        } else {
2337#ifdef CLP_DEBUG
2338          if (handler_->logLevel() & 32)
2339            printf("unflagging as djNorm %g %g, %d passes\n", djNorm, djNorm0, nPasses);
2340#endif
2341          if (pivotMode >= 10) {
2342            pivotMode--;
2343#ifdef CLP_DEBUG
2344            if (handler_->logLevel() & 32)
2345              printf("pivot mode now %d\n", pivotMode);
2346#endif
2347            if (pivotMode == 9)
2348              pivotMode = 0; // switch off fast attempt
2349          }
2350          bool unflagged = unflag() != 0;
2351          if (!unflagged && djNorm < 1.0e-10) {
2352            // ? declare victory
2353            sequenceIn_ = -1;
2354            returnCode = 0;
2355          }
2356        }
2357      }
2358    }
2359  }
2360  if (djNorm0 < 1.0e-12 * normFlagged) {
2361#ifdef CLP_DEBUG
2362    if (handler_->logLevel() & 32)
2363      printf("unflagging as djNorm %g %g and flagged norm %g\n", djNorm, djNorm0, normFlagged);
2364#endif
2365    unflag();
2366  }
2367  if (saveObj - currentObj < 1.0e-5 && nTotalPasses > 2000) {
2368    normUnflagged = 0.0;
2369    double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_);
2370    for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) {
2371      switch (getStatus(iSequence)) {
2372
2373      case basic:
2374      case ClpSimplex::isFixed:
2375        break;
2376      case atUpperBound:
2377        if (dj_[iSequence] > dualTolerance3)
2378          normFlagged += dj_[iSequence] * dj_[iSequence];
2379        break;
2380      case atLowerBound:
2381        if (dj_[iSequence] < -dualTolerance3)
2382          normFlagged += dj_[iSequence] * dj_[iSequence];
2383        break;
2384      case isFree:
2385      case superBasic:
2386        if (fabs(dj_[iSequence]) > dualTolerance3)
2387          normFlagged += dj_[iSequence] * dj_[iSequence];
2388        break;
2389      }
2390    }
2391    if (handler_->logLevel() > 2)
2392      printf("possible optimal  %d %d %g %g\n",
2393        nBigPasses, nTotalPasses, saveObj - currentObj, normFlagged);
2394    if (normFlagged < 1.0e-5) {
2395      unflag();
2396      primalColumnPivot_->setLooksOptimal(true);
2397      returnCode = 0;
2398    }
2399  }
2400  return returnCode;
2401}
2402/* Do last half of an iteration.
2403   Return codes
2404   Reasons to come out normal mode
2405   -1 normal
2406   -2 factorize now - good iteration
2407   -3 slight inaccuracy - refactorize - iteration done
2408   -4 inaccuracy - refactorize - no iteration
2409   -5 something flagged - go round again
2410   +2 looks unbounded
2411   +3 max iterations (iteration done)
2412
2413*/
2414int ClpSimplexNonlinear::pivotNonlinearResult()
2415{
2416
2417  int returnCode = -1;
2418
2419  rowArray_[1]->clear();
2420
2421  // we found a pivot column
2422  // update the incoming column
2423  unpackPacked(rowArray_[1]);
2424  factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
2425  theta_ = 0.0;
2426  double *work = rowArray_[1]->denseVector();
2427  int number = rowArray_[1]->getNumElements();
2428  int *which = rowArray_[1]->getIndices();
2429  bool keepValue = false;
2430  double saveValue = 0.0;
2431  if (pivotRow_ >= 0) {
2432    sequenceOut_ = pivotVariable_[pivotRow_];
2433    ;
2434    valueOut_ = solution(sequenceOut_);
2435    keepValue = true;
2436    saveValue = valueOut_;
2437    lowerOut_ = lower_[sequenceOut_];
2438    upperOut_ = upper_[sequenceOut_];
2439    int iIndex;
2440    for (iIndex = 0; iIndex < number; iIndex++) {
2441
2442      int iRow = which[iIndex];
2443      if (iRow == pivotRow_) {
2444        alpha_ = work[iIndex];
2445        break;
2446      }
2447    }
2448  } else {
2449    int iIndex;
2450    double smallest = COIN_DBL_MAX;
2451    for (iIndex = 0; iIndex < number; iIndex++) {
2452      int iRow = which[iIndex];
2453      double alpha = work[iIndex];
2454      if (fabs(alpha) > 1.0e-6) {
2455        int iPivot = pivotVariable_[iRow];
2456        double distance = CoinMin(upper_[iPivot] - solution_[iPivot],
2457          solution_[iPivot] - lower_[iPivot]);
2458        if (distance < smallest) {
2459          pivotRow_ = iRow;
2460          alpha_ = alpha;
2461          smallest = distance;
2462        }
2463      }
2464    }
2465    if (smallest > primalTolerance_) {
2466      smallest = COIN_DBL_MAX;
2467      for (iIndex = 0; iIndex < number; iIndex++) {
2468        int iRow = which[iIndex];
2469        double alpha = work[iIndex];
2470        if (fabs(alpha) > 1.0e-6) {
2471          double distance = randomNumberGenerator_.randomDouble();
2472          if (distance < smallest) {
2473            pivotRow_ = iRow;
2474            alpha_ = alpha;
2475            smallest = distance;
2476          }
2477        }
2478      }
2479    }
2480    assert(pivotRow_ >= 0);
2481    sequenceOut_ = pivotVariable_[pivotRow_];
2482    ;
2483    valueOut_ = solution(sequenceOut_);
2484    lowerOut_ = lower_[sequenceOut_];
2485    upperOut_ = upper_[sequenceOut_];
2486  }
2487  double newValue = valueOut_ - theta_ * alpha_;
2488  bool isSuperBasic = false;
2489  if (valueOut_ >= upperOut_ - primalTolerance_) {
2490    directionOut_ = -1; // to upper bound
2491    upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2492    upperOut_ = newValue;
2493  } else if (valueOut_ <= lowerOut_ + primalTolerance_) {
2494    directionOut_ = 1; // to lower bound
2495    lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2496  } else {
2497    lowerOut_ = valueOut_;
2498    upperOut_ = valueOut_;
2499    isSuperBasic = true;
2500    //printf("XX superbasic out\n");
2501  }
2502  dualOut_ = dj_[sequenceOut_];
2503  // if stable replace in basis
2504
2505  int updateStatus = factorization_->replaceColumn(this,
2506    rowArray_[2],
2507    rowArray_[1],
2508    pivotRow_,
2509    alpha_);
2510
2511  // if no pivots, bad update but reasonable alpha - take and invert
2512  if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
2513    updateStatus = 4;
2514  if (updateStatus == 1 || updateStatus == 4) {
2515    // slight error
2516    if (factorization_->pivots() > 5 || updateStatus == 4) {
2517      returnCode = -3;
2518    }
2519  } else if (updateStatus == 2) {
2520    // major error
2521    // better to have small tolerance even if slower
2522    factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
2523    int maxFactor = factorization_->maximumPivots();
2524    if (maxFactor > 10) {
2525      if (forceFactorization_ < 0)
2526        forceFactorization_ = maxFactor;
2527      forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
2528    }
2529    // later we may need to unwind more e.g. fake bounds
2530    if (lastGoodIteration_ != numberIterations_) {
2531      clearAll();
2532      pivotRow_ = -1;
2533      returnCode = -4;
2534    } else {
2535      // need to reject something
2536      char x = isColumn(sequenceIn_) ? 'C' : 'R';
2537      handler_->message(CLP_SIMPLEX_FLAG, messages_)
2538        << x << sequenceWithin(sequenceIn_)
2539        << CoinMessageEol;
2540      setFlagged(sequenceIn_);
2541      progress_.clearBadTimes();
2542      lastBadIteration_ = numberIterations_; // say be more cautious
2543      clearAll();
2544      pivotRow_ = -1;
2545      sequenceOut_ = -1;
2546      returnCode = -5;
2547    }
2548    return returnCode;
2549  } else if (updateStatus == 3) {
2550    // out of memory
2551    // increase space if not many iterations
2552    if (factorization_->pivots() < 0.5 * factorization_->maximumPivots() && factorization_->pivots() < 200)
2553      factorization_->areaFactor(
2554        factorization_->areaFactor() * 1.1);
2555    returnCode = -2; // factorize now
2556  } else if (updateStatus == 5) {
2557    problemStatus_ = -2; // factorize now
2558  }
2559
2560  // update primal solution
2561
2562  double objectiveChange = 0.0;
2563  // after this rowArray_[1] is not empty - used to update djs
2564  // If pivot row >= numberRows then may be gub
2565  updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, 1);
2566
2567  double oldValue = valueIn_;
2568  if (directionIn_ == -1) {
2569    // as if from upper bound
2570    if (sequenceIn_ != sequenceOut_) {
2571      // variable becoming basic
2572      valueIn_ -= fabs(theta_);
2573    } else {
2574      valueIn_ = lowerIn_;
2575    }
2576  } else {
2577    // as if from lower bound
2578    if (sequenceIn_ != sequenceOut_) {
2579      // variable becoming basic
2580      valueIn_ += fabs(theta_);
2581    } else {
2582      valueIn_ = upperIn_;
2583    }
2584  }
2585  objectiveChange += dualIn_ * (valueIn_ - oldValue);
2586  // outgoing
2587  if (sequenceIn_ != sequenceOut_) {
2588    if (directionOut_ > 0) {
2589      valueOut_ = lowerOut_;
2590    } else {
2591      valueOut_ = upperOut_;
2592    }
2593    if (valueOut_ < lower_[sequenceOut_] - primalTolerance_)
2594      valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
2595    else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
2596      valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
2597    // may not be exactly at bound and bounds may have changed
2598    // Make sure outgoing looks feasible
2599    if (!isSuperBasic)
2600      directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
2601    solution_[sequenceOut_] = valueOut_;
2602  }
2603  // change cost and bounds on incoming if primal
2604  nonLinearCost_->setOne(sequenceIn_, valueIn_);
2605  int whatNext = housekeeping(objectiveChange);
2606  if (keepValue)
2607    solution_[sequenceOut_] = saveValue;
2608  if (isSuperBasic)
2609    setStatus(sequenceOut_, superBasic);
2610    //#define CLP_DEBUG
2611#if CLP_DEBUG > 1
2612  {
2613    int ninf = matrix_->checkFeasible(this);
2614    if (ninf)
2615      printf("infeas %d\n", ninf);
2616  }
2617#endif
2618  if (whatNext == 1) {
2619    returnCode = -2; // refactorize
2620  } else if (whatNext == 2) {
2621    // maximum iterations or equivalent
2622    returnCode = 3;
2623  } else if (numberIterations_ == lastGoodIteration_ + 2 * factorization_->maximumPivots()) {
2624    // done a lot of flips - be safe
2625    returnCode = -2; // refactorize
2626  }
2627  // Check event
2628  {
2629    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2630    if (status >= 0) {
2631      problemStatus_ = 5;
2632      secondaryStatus_ = ClpEventHandler::endOfIteration;
2633      returnCode = 4;
2634    }
2635  }
2636  return returnCode;
2637}
2638// May use a cut approach for solving any LP
2639int ClpSimplexNonlinear::primalDualCuts(char *rowsIn, int startUp, int algorithm)
2640{
2641  if (!rowsIn) {
2642    if (algorithm > 0)
2643      return ClpSimplex::primal(startUp);
2644    else
2645      //return static_cast<ClpSimplexDual *>(static_cast<ClpSimplex *>(this))->dual(startUp);
2646      return ClpSimplex::dual(startUp);
2647  } else {
2648    int numberUsed = 0;
2649    int rowsThreshold = CoinMax(100, numberRows_ / 2);
2650    //int rowsTry=CoinMax(50,numberRows_/4);
2651    // Just add this number of rows each time in small problem
2652    int smallNumberRows = 2 * numberColumns_;
2653    smallNumberRows = CoinMin(smallNumberRows, numberRows_ / 20);
2654    // We will need arrays to choose rows to add
2655    double *weight = new double[numberRows_];
2656    int *sort = new int[numberRows_ + numberColumns_];
2657    int *whichColumns = sort + numberRows_;
2658    int numberSort = 0;
2659    for (int i = 0; i < numberRows_; i++) {
2660      if (rowsIn[i])
2661        numberUsed++;
2662    }
2663    if (numberUsed) {
2664      // normal
2665    } else {
2666      // If non slack use that information
2667      int numberBinding = 0;
2668      numberPrimalInfeasibilities_ = 0;
2669      sumPrimalInfeasibilities_ = 0.0;
2670      memset(rowActivity_, 0, numberRows_ * sizeof(double));
2671      times(1.0, columnActivity_, rowActivity_);
2672      for (int i = 0; i < numberRows_; i++) {
2673        double lowerDifference = rowActivity_[i] - rowLower_[i];
2674        double upperDifference = rowActivity_[i] - rowUpper_[i];
2675        if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
2676          numberPrimalInfeasibilities_++;
2677          if (lowerDifference < 0.0)
2678            sumPrimalInfeasibilities_ -= lowerDifference;
2679          else
2680            sumPrimalInfeasibilities_ += upperDifference;
2681          rowsIn[i] = 1;
2682        } else if (getRowStatus(i) != basic) {
2683          numberBinding++;
2684          rowsIn[i] = 1;
2685        }
2686      }
2687      if (numberBinding < rowsThreshold) {
2688        // try
2689      } else {
2690        // random?? - or give up
2691        // Set up initial list
2692        numberSort = 0;
2693        for (int iRow = 0; iRow < numberRows_; iRow++) {
2694          weight[iRow] = 1.123e50;
2695          if (rowLower_[iRow] == rowUpper_[iRow]) {
2696            sort[numberSort++] = iRow;
2697            weight[iRow] = 0.0;
2698          }
2699        }
2700        numberSort /= 2;
2701        // and pad out with random rows
2702        double ratio = ((double)(smallNumberRows - numberSort)) / ((double)numberRows_);
2703        for (int iRow = 0; iRow < numberRows_; iRow++) {
2704          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
2705            sort[numberSort++] = iRow;
2706        }
2707        // sort
2708        CoinSort_2(weight, weight + numberRows_, sort);
2709        numberSort = CoinMin(numberRows_, smallNumberRows);
2710        memset(rowsIn, 0, numberRows_);
2711        for (int iRow = 0; iRow < numberSort; iRow++)
2712          rowsIn[sort[iRow]] = 1;
2713      }
2714    }
2715    // see if feasible
2716    numberPrimalInfeasibilities_ = 0;
2717    memset(rowActivity_, 0, numberRows_ * sizeof(double));
2718    times(1.0, columnActivity_, rowActivity_);
2719    for (int i = 0; i < numberRows_; i++) {
2720      double lowerDifference = rowActivity_[i] - rowLower_[i];
2721      double upperDifference = rowActivity_[i] - rowUpper_[i];
2722      if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) {
2723        if (lowerDifference < 0.0)
2724          sumPrimalInfeasibilities_ -= lowerDifference;
2725        else
2726          sumPrimalInfeasibilities_ += upperDifference;
2727        numberPrimalInfeasibilities_++;
2728      }
2729    }
2730    printf("Initial infeasibilities - %g (%d)\n",
2731      sumPrimalInfeasibilities_, numberPrimalInfeasibilities_);
2732    // Just do this number of passes
2733    int maxPass = 50;
2734    // And take out slack rows until this pass
2735    int takeOutPass = 30;
2736    int iPass;
2737
2738    const CoinBigIndex *start = this->clpMatrix()->getVectorStarts();
2739    const int *length = this->clpMatrix()->getVectorLengths();
2740    const int *row = this->clpMatrix()->getIndices();
2741    problemStatus_ = 1;
2742    for (iPass = 0; iPass < maxPass; iPass++) {
2743      printf("Start of pass %d\n", iPass);
2744      int numberSmallColumns = 0;
2745      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2746        int n = 0;
2747        for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
2748          int iRow = row[j];
2749          if (rowsIn[iRow])
2750            n++;
2751        }
2752        if (n)
2753          whichColumns[numberSmallColumns++] = iColumn;
2754      }
2755      numberSort = 0;
2756      for (int i = 0; i < numberRows_; i++) {
2757        if (rowsIn[i])
2758          sort[numberSort++] = i;
2759      }
2760      // Create small problem
2761      ClpSimplex small(this, numberSort, sort, numberSmallColumns, whichColumns);
2762      printf("Small model has %d rows, %d columns and %d elements\n",
2763        small.numberRows(), small.numberColumns(), small.getNumElements());
2764      small.setFactorizationFrequency(100 + numberSort / 200);
2765      // Solve
2766      small.setLogLevel(CoinMax(0, logLevel() - 1));
2767      if (iPass > 20) {
2768        if (sumPrimalInfeasibilities_ > 1.0e-1) {
2769          small.dual();
2770        } else {
2771          small.primal(1);
2772        }
2773      } else {
2774        // presolve
2775#if 0
2776        ClpSolve::SolveType method;
2777        ClpSolve::PresolveType presolveType = ClpSolve::presolveOn;
2778        ClpSolve solveOptions;
2779        solveOptions.setPresolveType(presolveType, 5);
2780        if (sumPrimalInfeasibilities_>1.0e-1)
2781          method = ClpSolve::useDual;
2782        else
2783          method = ClpSolve::usePrimalorSprint;
2784        solveOptions.setSolveType(method);
2785        small.initialSolve(solveOptions);
2786#else
2787#if 1
2788        ClpPresolve *pinfo = new ClpPresolve();
2789        ClpSimplex *small2 = pinfo->presolvedModel(small, 1.0e-5);
2790        assert(small2);
2791        if (sumPrimalInfeasibilities_ > 1.0e-1) {
2792          small2->dual();
2793        } else {
2794          small2->primal(1);
2795        }
2796        pinfo->postsolve(true);
2797        delete pinfo;
2798#else
2799        char *types = new char[small.numberRows() + small.numberColumns()];
2800        memset(types, 0, small.numberRows() + small.numberColumns());
2801        if (sumPrimalInfeasibilities_ > 1.0e-1) {
2802          small.miniSolve(types, types + small.numberRows(),
2803            -1, 0);
2804        } else {
2805          small.miniSolve(types, types + small.numberRows(),
2806            1, 1);
2807        }
2808        delete[] types;
2809#endif
2810        if (small.sumPrimalInfeasibilities() > 1.0)
2811          small.primal(1);
2812#endif
2813      }
2814      bool dualInfeasible = (small.status() == 2);
2815      // move solution back
2816      const double *smallSolution = small.primalColumnSolution();
2817      for (int j = 0; j < numberSmallColumns; j++) {
2818        int iColumn = whichColumns[j];
2819        columnActivity_[iColumn] = smallSolution[j];
2820        this->setColumnStatus(iColumn, small.getColumnStatus(j));
2821      }
2822      for (int iRow = 0; iRow < numberSort; iRow++) {
2823        int kRow = sort[iRow];
2824        this->setRowStatus(kRow, small.getRowStatus(iRow));
2825      }
2826      // compute full solution
2827      memset(rowActivity_, 0, numberRows_ * sizeof(double));
2828      times(1.0, columnActivity_, rowActivity_);
2829      if (iPass != maxPass - 1) {
2830        // Mark row as not looked at
2831        for (int iRow = 0; iRow < numberRows_; iRow++)
2832          weight[iRow] = 1.123e50;
2833        // Look at rows already in small problem
2834        int iSort;
2835        int numberDropped = 0;
2836        int numberKept = 0;
2837        int numberBinding = 0;
2838        numberPrimalInfeasibilities_ = 0;
2839        sumPrimalInfeasibilities_ = 0.0;
2840        bool allFeasible = small.numberIterations() == 0;
2841        for (iSort = 0; iSort < numberSort; iSort++) {
2842          int iRow = sort[iSort];
2843          //printf("%d %g %g\n",iRow,rowActivity_[iRow],small.primalRowSolution()[iSort]);
2844          if (getRowStatus(iRow) == ClpSimplex::basic) {
2845            // Basic - we can get rid of if early on
2846            if (iPass < takeOutPass && !dualInfeasible) {
2847              double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
2848                rowLower_[iRow] - rowActivity_[iRow]);
2849              weight[iRow] = -infeasibility;
2850              if (infeasibility > primalTolerance_ && !allFeasible) {
2851                numberPrimalInfeasibilities_++;
2852                sumPrimalInfeasibilities_ += infeasibility;
2853              } else {
2854                weight[iRow] = 1.0;
2855                numberDropped++;
2856              }
2857            } else {
2858              // keep
2859              weight[iRow] = -1.0e40;
2860              numberKept++;
2861            }
2862          } else {
2863            // keep
2864            weight[iRow] = -1.0e50;
2865            numberKept++;
2866            numberBinding++;
2867          }
2868        }
2869        // Now rest
2870        for (int iRow = 0; iRow < numberRows_; iRow++) {
2871          sort[iRow] = iRow;
2872          if (weight[iRow] == 1.123e50) {
2873            // not looked at yet
2874            double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
2875              rowLower_[iRow] - rowActivity_[iRow]);
2876            weight[iRow] = -infeasibility;
2877            if (infeasibility > primalTolerance_) {
2878              numberPrimalInfeasibilities_++;
2879              sumPrimalInfeasibilities_ += infeasibility;
2880            }
2881          }
2882        }
2883        // sort
2884        CoinSort_2(weight, weight + numberRows_, sort);
2885        numberSort = CoinMin(numberRows_, smallNumberRows + numberKept);
2886        memset(rowsIn, 0, numberRows_);
2887        for (int iRow = 0; iRow < numberSort; iRow++)
2888          rowsIn[sort[iRow]] = 1;
2889        printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
2890          numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns);
2891        printf("%d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
2892          sumPrimalInfeasibilities_);
2893        if (!numberPrimalInfeasibilities_) {
2894          problemStatus_ = 0;
2895          printf("Exiting as looks optimal\n");
2896          break;
2897        }
2898        numberPrimalInfeasibilities_ = 0;
2899        sumPrimalInfeasibilities_ = 0.0;
2900        for (iSort = 0; iSort < numberSort; iSort++) {
2901          if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) {
2902            numberPrimalInfeasibilities_++;
2903            sumPrimalInfeasibilities_ += -weight[iSort];
2904          }
2905        }
2906        printf("in small model %d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
2907          sumPrimalInfeasibilities_);
2908      } else {
2909        // out of passes
2910        problemStatus_ = -1;
2911      }
2912    }
2913    delete[] weight;
2914    delete[] sort;
2915    return 0;
2916  }
2917}
2918// A sequential LP method
2919int ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance,
2920  int otherOptions)
2921{
2922  // Are we minimizing or maximizing
2923  double whichWay = optimizationDirection();
2924  if (whichWay < 0.0)
2925    whichWay = -1.0;
2926  else if (whichWay > 0.0)
2927    whichWay = 1.0;
2928
2929  int numberColumns = this->numberColumns();
2930  int numberRows = this->numberRows();
2931  double *columnLower = this->columnLower();
2932  double *columnUpper = this->columnUpper();
2933  double *solution = this->primalColumnSolution();
2934
2935  if (objective_->type() < 2) {
2936    // no nonlinear part
2937    return ClpSimplex::primal(0);
2938  }
2939  // Get list of non linear columns
2940  char *markNonlinear = new char[numberColumns];
2941  memset(markNonlinear, 0, numberColumns);
2942  int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
2943
2944  if (!numberNonLinearColumns) {
2945    delete[] markNonlinear;
2946    // no nonlinear part
2947    return ClpSimplex::primal(0);
2948  }
2949  int iColumn;
2950  int *listNonLinearColumn = new int[numberNonLinearColumns];
2951  numberNonLinearColumns = 0;
2952  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2953    if (markNonlinear[iColumn])
2954      listNonLinearColumn[numberNonLinearColumns++] = iColumn;
2955  }
2956  // Replace objective
2957  ClpObjective *trueObjective = objective_;
2958  objective_ = new ClpLinearObjective(NULL, numberColumns);
2959  double *objective = this->objective();
2960  // See if user wants to use cuts
2961  char *rowsIn = NULL;
2962  if ((otherOptions & 1) != 0 || numberPasses < 0) {
2963    numberPasses = abs(numberPasses);
2964    rowsIn = new char[numberRows_];
2965    memset(rowsIn, 0, numberRows_);
2966  }
2967  // get feasible
2968  if (this->status() < 0 || numberPrimalInfeasibilities())
2969    primalDualCuts(rowsIn, 1, 1);
2970  // still infeasible
2971  if (problemStatus_) {
2972    delete[] listNonLinearColumn;
2973    return 0;
2974  }
2975  algorithm_ = 1;
2976  int jNon;
2977  int *last[3];
2978
2979  double *trust = new double[numberNonLinearColumns];
2980  double *trueLower = new double[numberNonLinearColumns];
2981  double *trueUpper = new double[numberNonLinearColumns];
2982  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
2983    iColumn = listNonLinearColumn[jNon];
2984    trust[jNon] = 0.5;
2985    trueLower[jNon] = columnLower[iColumn];
2986    trueUpper[jNon] = columnUpper[iColumn];
2987    if (solution[iColumn] < trueLower[jNon])
2988      solution[iColumn] = trueLower[jNon];
2989    else if (solution[iColumn] > trueUpper[jNon])
2990      solution[iColumn] = trueUpper[jNon];
2991  }
2992  int saveLogLevel = logLevel();
2993  int iPass;
2994  double lastObjective = 1.0e31;
2995  double *saveSolution = new double[numberColumns];
2996  double *saveRowSolution = new double[numberRows];
2997  memset(saveRowSolution, 0, numberRows * sizeof(double));
2998  double *savePi = new double[numberRows];
2999  double *safeSolution = new double[numberColumns];
3000  unsigned char *saveStatus = new unsigned char[numberRows + numberColumns];
3001#define MULTIPLE 0
3002#if MULTIPLE > 2
3003  // Duplication but doesn't really matter
3004     double * saveSolutionM[MULTIPLE
3005};
3006for (jNon = 0; jNon < MULTIPLE; jNon++) {
3007  saveSolutionM[jNon] = new double[numberColumns];
3008  CoinMemcpyN(solution, numberColumns, saveSolutionM);
3009}
3010#endif
3011double targetDrop = 1.0e31;
3012double objectiveOffset;
3013getDblParam(ClpObjOffset, objectiveOffset);
3014// 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
3015for (iPass = 0; iPass < 3; iPass++) {
3016  last[iPass] = new int[numberNonLinearColumns];
3017  for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
3018    last[iPass][jNon] = 0;
3019}
3020// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3021int goodMove = -2;
3022char *statusCheck = new char[numberColumns];
3023double *changeRegion = new double[numberColumns];
3024double offset = 0.0;
3025double objValue = 0.0;
3026int exitPass = 2 * numberPasses + 10;
3027for (iPass = 0; iPass < numberPasses; iPass++) {
3028  exitPass--;
3029  // redo objective
3030  offset = 0.0;
3031  objValue = -objectiveOffset;
3032  // make sure x updated
3033  trueObjective->newXValues();
3034  double theta = -1.0;
3035  double maxTheta = COIN_DBL_MAX;
3036  //maxTheta=1.0;
3037  if (iPass) {
3038    int jNon = 0;
3039    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3040      changeRegion[iColumn] = solution[iColumn] - saveSolution[iColumn];
3041      double alpha = changeRegion[iColumn];
3042      double oldValue = saveSolution[iColumn];
3043      if (markNonlinear[iColumn] == 0) {
3044        // linear
3045        if (alpha < -1.0e-15) {
3046          // variable going towards lower bound
3047          double bound = columnLower[iColumn];
3048          oldValue -= bound;
3049          if (oldValue + maxTheta * alpha < 0.0) {
3050            maxTheta = CoinMax(0.0, oldValue / (-alpha));
3051          }
3052        } else if (alpha > 1.0e-15) {
3053          // variable going towards upper bound
3054          double bound = columnUpper[iColumn];
3055          oldValue = bound - oldValue;
3056          if (oldValue - maxTheta * alpha < 0.0) {
3057            maxTheta = CoinMax(0.0, oldValue / alpha);
3058          }
3059        }
3060      } else {
3061        // nonlinear
3062        if (alpha < -1.0e-15) {
3063          // variable going towards lower bound
3064          double bound = trueLower[jNon];
3065          oldValue -= bound;
3066          if (oldValue + maxTheta * alpha < 0.0) {
3067            maxTheta = CoinMax(0.0, oldValue / (-alpha));
3068          }
3069        } else if (alpha > 1.0e-15) {
3070          // variable going towards upper bound
3071          double bound = trueUpper[jNon];
3072          oldValue = bound - oldValue;
3073          if (oldValue - maxTheta * alpha < 0.0) {
3074            maxTheta = CoinMax(0.0, oldValue / alpha);
3075          }
3076        }
3077        jNon++;
3078      }
3079    }
3080    // make sure both accurate
3081    memset(rowActivity_, 0, numberRows_ * sizeof(double));
3082    times(1.0, solution, rowActivity_);
3083    memset(saveRowSolution, 0, numberRows_ * sizeof(double));
3084    times(1.0, saveSolution, saveRowSolution);
3085    for (int iRow = 0; iRow < numberRows; iRow++) {
3086      double alpha = rowActivity_[iRow] - saveRowSolution[iRow];
3087      double oldValue = saveRowSolution[iRow];
3088      if (alpha < -1.0e-15) {
3089        // variable going towards lower bound
3090        double bound = rowLower_[iRow];
3091        oldValue -= bound;
3092        if (oldValue + maxTheta * alpha < 0.0) {
3093          maxTheta = CoinMax(0.0, oldValue / (-alpha));
3094        }
3095      } else if (alpha > 1.0e-15) {
3096        // variable going towards upper bound
3097        double bound = rowUpper_[iRow];
3098        oldValue = bound - oldValue;
3099        if (oldValue - maxTheta * alpha < 0.0) {
3100          maxTheta = CoinMax(0.0, oldValue / alpha);
3101        }
3102      }
3103    }
3104  } else {
3105    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3106      changeRegion[iColumn] = 0.0;
3107      saveSolution[iColumn] = solution[iColumn];
3108    }
3109    CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
3110  }
3111  // get current value anyway
3112  double predictedObj, thetaObj;
3113  double maxTheta2 = 2.0; // to work out a b c
3114  double theta2 = trueObjective->stepLength(this, saveSolution, changeRegion, maxTheta2,
3115    objValue, predictedObj, thetaObj);
3116  int lastMoveStatus = goodMove;
3117  if (goodMove >= 0) {
3118    theta = CoinMin(theta2, maxTheta);
3119#ifdef CLP_DEBUG
3120    if (handler_->logLevel() & 32)
3121      printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
3122        theta, objValue, thetaObj, predictedObj);
3123#endif
3124    if (theta > 0.0 && theta <= 1.0) {
3125      // update solution
3126      double lambda = 1.0 - theta;
3127      for (iColumn = 0; iColumn < numberColumns; iColumn++)
3128        solution[iColumn] = lambda * saveSolution[iColumn]
3129          + theta * solution[iColumn];
3130      memset(rowActivity_, 0, numberRows_ * sizeof(double));
3131      times(1.0, solution, rowActivity_);
3132      if (lambda > 0.999) {
3133        CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
3134        CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
3135      }
3136      // Do local minimization
3137#define LOCAL
3138#ifdef LOCAL
3139      bool absolutelyOptimal = false;
3140      int saveScaling = scalingFlag();
3141      scaling(0);
3142      int savePerturbation = perturbation_;
3143      perturbation_ = 100;
3144      if (saveLogLevel == 1)
3145        setLogLevel(0);
3146      int status = startup(1);
3147      if (!status) {
3148        int numberTotal = numberRows_ + numberColumns_;
3149        // resize arrays
3150        for (int i = 0; i < 4; i++) {
3151          rowArray_[i]->reserve(CoinMax(numberRows_ + numberColumns_, rowArray_[i]->capacity()));
3152        }
3153        CoinIndexedVector *longArray = rowArray_[3];
3154        CoinIndexedVector *rowArray = rowArray_[0];
3155        //CoinIndexedVector * columnArray = columnArray_[0];
3156        CoinIndexedVector *spare = rowArray_[1];
3157        double *work = longArray->denseVector();
3158        int *which = longArray->getIndices();
3159        int nPass = 100;
3160        //bool conjugate=false;
3161        // Put back true bounds
3162        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3163          int iColumn = listNonLinearColumn[jNon];
3164          double value;
3165          value = trueLower[jNon];
3166          trueLower[jNon] = lower_[iColumn];
3167          lower_[iColumn] = value;
3168          value = trueUpper[jNon];
3169          trueUpper[jNon] = upper_[iColumn];
3170          upper_[iColumn] = value;
3171          switch (getStatus(iColumn)) {
3172
3173          case basic:
3174          case isFree:
3175          case superBasic:
3176            break;
3177          case isFixed:
3178          case atUpperBound:
3179          case atLowerBound:
3180            if (solution_[iColumn] > lower_[iColumn] + primalTolerance_) {
3181              if (solution_[iColumn] < upper_[iColumn] - primalTolerance_) {
3182                setStatus(iColumn, superBasic);
3183              } else {
3184                setStatus(iColumn, atUpperBound);
3185              }
3186            } else {
3187              if (solution_[iColumn] < upper_[iColumn] - primalTolerance_) {
3188                setStatus(iColumn, atLowerBound);
3189              } else {
3190                setStatus(iColumn, isFixed);
3191              }
3192            }
3193            break;
3194          }
3195        }
3196        for (int jPass = 0; jPass < nPass; jPass++) {
3197          trueObjective->reducedGradient(this, dj_, true);
3198          // get direction vector in longArray
3199          longArray->clear();
3200          double normFlagged = 0.0;
3201          double normUnflagged = 0.0;
3202          int numberNonBasic = 0;
3203          directionVector(longArray, spare, rowArray, 0,
3204            normFlagged, normUnflagged, numberNonBasic);
3205          if (normFlagged + normUnflagged < 1.0e-8) {
3206            absolutelyOptimal = true;
3207            break; //looks optimal
3208          }
3209          double djNorm = 0.0;
3210          int iIndex;
3211          for (iIndex = 0; iIndex < numberNonBasic; iIndex++) {
3212            int iSequence = which[iIndex];
3213            double alpha = work[iSequence];
3214            djNorm += alpha * alpha;
3215          }
3216          //if (!jPass)
3217          //printf("dj norm %g - %d \n",djNorm,numberNonBasic);
3218          //int number=longArray->getNumElements();
3219          if (!numberNonBasic) {
3220            // looks optimal
3221            absolutelyOptimal = true;
3222            break;
3223          }
3224          double theta = 1.0e30;
3225          int iSequence;
3226          for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3227            double alpha = work[iSequence];
3228            double oldValue = solution_[iSequence];
3229            if (alpha < -1.0e-15) {
3230              // variable going towards lower bound
3231              double bound = lower_[iSequence];
3232              oldValue -= bound;
3233              if (oldValue + theta * alpha < 0.0) {
3234                theta = CoinMax(0.0, oldValue / (-alpha));
3235              }
3236            } else if (alpha > 1.0e-15) {
3237              // variable going towards upper bound
3238              double bound = upper_[iSequence];
3239              oldValue = bound - oldValue;
3240              if (oldValue - theta * alpha < 0.0) {
3241                theta = CoinMax(0.0, oldValue / alpha);
3242              }
3243            }
3244          }
3245          // Now find minimum of function
3246          double currentObj;
3247          double predictedObj;
3248          double thetaObj;
3249          // need to use true objective
3250          double *saveCost = cost_;
3251          cost_ = NULL;
3252          double objTheta = trueObjective->stepLength(this, solution_, work, theta,
3253            currentObj, predictedObj, thetaObj);
3254          cost_ = saveCost;
3255#ifdef CLP_DEBUG
3256          if (handler_->logLevel() & 32)
3257            printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
3258#endif
3259          //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
3260          //printf("objTheta %g theta %g\n",objTheta,theta);
3261          if (theta > objTheta) {
3262            theta = objTheta;
3263            thetaObj = predictedObj;
3264          }
3265          // update one used outside
3266          objValue = currentObj;
3267          if (theta > 1.0e-9 && (currentObj - thetaObj < -CoinMax(1.0e-8, 1.0e-15 * fabs(currentObj)) || jPass < 5)) {
3268            // Update solution
3269            for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3270              double alpha = work[iSequence];
3271              if (alpha) {
3272                double value = solution_[iSequence] + theta * alpha;
3273                solution_[iSequence] = value;
3274                switch (getStatus(iSequence)) {
3275
3276                case basic:
3277                case isFixed:
3278                case isFree:
3279                  break;
3280                case atUpperBound:
3281                case atLowerBound:
3282                case superBasic:
3283                  if (value > lower_[iSequence] + primalTolerance_) {
3284                    if (value < upper_[iSequence] - primalTolerance_) {
3285                      setStatus(iSequence, superBasic);
3286                    } else {
3287                      setStatus(iSequence, atUpperBound);
3288                    }
3289                  } else {
3290                    if (value < upper_[iSequence] - primalTolerance_) {
3291                      setStatus(iSequence, atLowerBound);
3292                    } else {
3293                      setStatus(iSequence, isFixed);
3294                    }
3295                  }
3296                  break;
3297                }
3298              }
3299            }
3300          } else {
3301            break;
3302          }
3303        }
3304        // Put back fake bounds
3305        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3306          int iColumn = listNonLinearColumn[jNon];
3307          double value;
3308          value = trueLower[jNon];
3309          trueLower[jNon] = lower_[iColumn];
3310          lower_[iColumn] = value;
3311          value = trueUpper[jNon];
3312          trueUpper[jNon] = upper_[iColumn];
3313          upper_[iColumn] = value;
3314        }
3315      }
3316      problemStatus_ = 0;
3317      finish();
3318      scaling(saveScaling);
3319      perturbation_ = savePerturbation;
3320      setLogLevel(saveLogLevel);
3321#endif
3322      // redo rowActivity
3323      memset(rowActivity_, 0, numberRows_ * sizeof(double));
3324      times(1.0, solution, rowActivity_);
3325      if (theta > 0.99999 && theta2 < 1.9 && !absolutelyOptimal) {
3326        // If big changes then tighten
3327        /*  thetaObj is objvalue + a*2*2 +b*2
3328                        predictedObj is objvalue + a*theta2*theta2 +b*theta2
3329                    */
3330        double rhs1 = thetaObj - objValue;
3331        double rhs2 = predictedObj - objValue;
3332        double subtractB = theta2 * 0.5;
3333        double a = (rhs2 - subtractB * rhs1) / (theta2 * theta2 - 4.0 * subtractB);
3334        double b = 0.5 * (rhs1 - 4.0 * a);
3335        if (fabs(a + b) > 1.0e-2) {
3336          // tighten all
3337          goodMove = -1;
3338        }
3339      }
3340    }
3341  }
3342  CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns,
3343    objective);
3344  //printf("offset comp %g orig %g\n",offset,objectiveOffset);
3345  // could do faster
3346  trueObjective->stepLength(this, solution, changeRegion, 0.0,
3347    objValue, predictedObj, thetaObj);
3348#ifdef CLP_INVESTIGATE
3349  printf("offset comp %g orig %g - obj from stepLength %g\n", offset, objectiveOffset, objValue);
3350#endif
3351  setDblParam(ClpObjOffset, objectiveOffset + offset);
3352  int *temp = last[2];
3353  last[2] = last[1];
3354  last[1] = last[0];
3355  last[0] = temp;
3356  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3357    iColumn = listNonLinearColumn[jNon];
3358    double change = solution[iColumn] - saveSolution[iColumn];
3359    if (change < -1.0e-5) {
3360      if (fabs(change + trust[jNon]) < 1.0e-5)
3361        temp[jNon] = -1;
3362      else
3363        temp[jNon] = -2;
3364    } else if (change > 1.0e-5) {
3365      if (fabs(change - trust[jNon]) < 1.0e-5)
3366        temp[jNon] = 1;
3367      else
3368        temp[jNon] = 2;
3369    } else {
3370      temp[jNon] = 0;
3371    }
3372  }
3373  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3374  double maxDelta = 0.0;
3375  if (goodMove >= 0) {
3376    if (objValue - lastObjective <= 1.0e-15 * fabs(lastObjective))
3377      goodMove = 1;
3378    else
3379      goodMove = 0;
3380  } else {
3381    maxDelta = 1.0e10;
3382  }
3383  double maxGap = 0.0;
3384  int numberSmaller = 0;
3385  int numberSmaller2 = 0;
3386  int numberLarger = 0;
3387  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3388    iColumn = listNonLinearColumn[jNon];
3389    maxDelta = CoinMax(maxDelta,
3390      fabs(solution[iColumn] - saveSolution[iColumn]));
3391    if (goodMove > 0) {
3392      if (last[0][jNon] * last[1][jNon] < 0) {
3393        // halve
3394        trust[jNon] *= 0.5;
3395        numberSmaller2++;
3396      } else {
3397        if (last[0][jNon] == last[1][jNon] && last[0][jNon] == last[2][jNon])
3398          trust[jNon] = CoinMin(1.5 * trust[jNon], 1.0e6);
3399        numberLarger++;
3400      }
3401    } else if (goodMove != -2 && trust[jNon] > 10.0 * deltaTolerance) {
3402      trust[jNon] *= 0.2;
3403      numberSmaller++;
3404    }
3405    maxGap = CoinMax(maxGap, trust[jNon]);
3406  }
3407#ifdef CLP_DEBUG
3408  if (handler_->logLevel() & 32)
3409    std::cout << "largest gap is " << maxGap << " "
3410              << numberSmaller + numberSmaller2 << " reduced ("
3411              << numberSmaller << " badMove ), "
3412              << numberLarger << " increased" << std::endl;
3413#endif
3414  if (iPass > 10000) {
3415    for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
3416      trust[jNon] *= 0.0001;
3417  }
3418  if (lastMoveStatus == -1 && goodMove == -1)
3419    goodMove = 1; // to force solve
3420  if (goodMove > 0) {
3421    double drop = lastObjective - objValue;
3422    handler_->message(CLP_SLP_ITER, messages_)
3423      << iPass << objValue - objectiveOffset
3424      << drop << maxDelta
3425      << CoinMessageEol;
3426    if (iPass > 20 && drop < 1.0e-12 * fabs(objValue))
3427      drop = 0.999e-4; // so will exit
3428    if (maxDelta < deltaTolerance && drop < 1.0e-4 && goodMove && theta < 0.99999) {
3429      if (handler_->logLevel() > 1)
3430        std::cout << "Exiting as maxDelta < tolerance and small drop" << std::endl;
3431      break;
3432    }
3433  } else if (!numberSmaller && iPass > 1) {
3434    if (handler_->logLevel() > 1)
3435      std::cout << "Exiting as all gaps small" << std::endl;
3436    break;
3437  }
3438  if (!iPass)
3439    goodMove = 1;
3440  targetDrop = 0.0;
3441  double *r = this->dualColumnSolution();
3442  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3443    iColumn = listNonLinearColumn[jNon];
3444    columnLower[iColumn] = CoinMax(solution[iColumn]
3445        - trust[jNon],
3446      trueLower[jNon]);
3447    columnUpper[iColumn] = CoinMin(solution[iColumn]
3448        + trust[jNon],
3449      trueUpper[jNon]);
3450  }
3451  if (iPass) {
3452    // get reduced costs
3453    this->matrix()->transposeTimes(savePi,
3454      this->dualColumnSolution());
3455    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3456      iColumn = listNonLinearColumn[jNon];
3457      double dj = objective[iColumn] - r[iColumn];
3458      r[iColumn] = dj;
3459      if (dj < -dualTolerance_)
3460        targetDrop -= dj * (columnUpper[iColumn] - solution[iColumn]);
3461      else if (dj > dualTolerance_)
3462        targetDrop -= dj * (columnLower[iColumn] - solution[iColumn]);
3463    }
3464  } else {
3465    memset(r, 0, numberColumns * sizeof(double));
3466  }
3467#if 0
3468     for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3469          iColumn = listNonLinearColumn[jNon];
3470          if (statusCheck[iColumn] == 'L' && r[iColumn] < -1.0e-4) {
3471               columnLower[iColumn] = CoinMax(solution[iColumn],
3472                                              trueLower[jNon]);
3473               columnUpper[iColumn] = CoinMin(solution[iColumn]
3474                                              + trust[jNon],
3475                                              trueUpper[jNon]);
3476          } else if (statusCheck[iColumn] == 'U' && r[iColumn] > 1.0e-4) {
3477               columnLower[iColumn] = CoinMax(solution[iColumn]
3478                                              - trust[jNon],
3479                                              trueLower[jNon]);
3480               columnUpper[iColumn] = CoinMin(solution[iColumn],
3481                                              trueUpper[jNon]);
3482          } else {
3483               columnLower[iColumn] = CoinMax(solution[iColumn]
3484                                              - trust[jNon],
3485                                              trueLower[jNon]);
3486               columnUpper[iColumn] = CoinMin(solution[iColumn]
3487                                              + trust[jNon],
3488                                              trueUpper[jNon]);
3489          }
3490     }
3491#endif
3492  if (goodMove > 0) {
3493    CoinMemcpyN(solution, numberColumns, saveSolution);
3494    CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
3495    CoinMemcpyN(this->dualRowSolution(), numberRows, savePi);
3496    CoinMemcpyN(status_, numberRows + numberColumns, saveStatus);
3497#if MULTIPLE > 2
3498    double *tempSol = saveSolutionM[0];
3499    for (jNon = 0; jNon < MULTIPLE - 1; jNon++) {
3500      saveSolutionM[jNon] = saveSolutionM[jNon + 1];
3501    }
3502    saveSolutionM[MULTIPLE - 1] = tempSol;
3503    CoinMemcpyN(solution, numberColumns, tempSol);
3504
3505#endif
3506
3507#ifdef CLP_DEBUG
3508    if (handler_->logLevel() & 32)
3509      std::cout << "Pass - " << iPass
3510                << ", target drop is " << targetDrop
3511                << std::endl;
3512#endif
3513    lastObjective = objValue;
3514    if (targetDrop < CoinMax(1.0e-8, CoinMin(1.0e-6, 1.0e-6 * fabs(objValue))) && goodMove && iPass > 3) {
3515      if (handler_->logLevel() > 1)
3516        printf("Exiting on target drop %g\n", targetDrop);
3517      break;
3518    }
3519#ifdef CLP_DEBUG
3520    {
3521      double *r = this->dualColumnSolution();
3522      for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3523        iColumn = listNonLinearColumn[jNon];
3524        if (handler_->logLevel() & 32)
3525          printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
3526            jNon, trust[jNon], iColumn, solution[iColumn], objective[iColumn],
3527            r[iColumn], statusCheck[iColumn], columnLower[iColumn],
3528            columnUpper[iColumn]);
3529      }
3530    }
3531#endif
3532    //setLogLevel(63);
3533    this->scaling(false);
3534    if (saveLogLevel == 1)
3535      setLogLevel(0);
3536    primalDualCuts(rowsIn, 1, 1);
3537    algorithm_ = 1;
3538    setLogLevel(saveLogLevel);
3539#ifdef CLP_DEBUG
3540    if (this->status()) {
3541      writeMps("xx.mps");
3542    }
3543#endif
3544    if (this->status() == 1) {
3545      // not feasible ! - backtrack and exit
3546      // use safe solution
3547      CoinMemcpyN(safeSolution, numberColumns, solution);
3548      CoinMemcpyN(solution, numberColumns, saveSolution);
3549      memset(rowActivity_, 0, numberRows_ * sizeof(double));
3550      times(1.0, solution, rowActivity_);
3551      CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
3552      CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
3553      CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
3554      for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3555        iColumn = listNonLinearColumn[jNon];
3556        columnLower[iColumn] = CoinMax(solution[iColumn]
3557            - trust[jNon],
3558          trueLower[jNon]);
3559        columnUpper[iColumn] = CoinMin(solution[iColumn]
3560            + trust[jNon],
3561          trueUpper[jNon]);
3562      }
3563      break;
3564    } else {
3565      // save in case problems
3566      CoinMemcpyN(solution, numberColumns, safeSolution);
3567    }
3568    if (problemStatus_ == 3)
3569      break; // probably user interrupt
3570    goodMove = 1;
3571  } else {
3572    // bad pass - restore solution
3573#ifdef CLP_DEBUG
3574    if (handler_->logLevel() & 32)
3575      printf("Backtracking\n");
3576#endif
3577    CoinMemcpyN(saveSolution, numberColumns, solution);
3578    CoinMemcpyN(saveRowSolution, numberRows, rowActivity_);
3579    CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
3580    CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
3581    iPass--;
3582    assert(exitPass > 0);
3583    goodMove = -1;
3584  }
3585}
3586#if MULTIPLE > 2
3587for (jNon = 0; jNon < MULTIPLE; jNon++) {
3588  delete[] saveSolutionM[jNon];
3589}
3590#endif
3591// restore solution
3592CoinMemcpyN(saveSolution, numberColumns, solution);
3593CoinMemcpyN(saveRowSolution, numberRows, rowActivity_);
3594for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3595  iColumn = listNonLinearColumn[jNon];
3596  columnLower[iColumn] = CoinMax(solution[iColumn],
3597    trueLower[jNon]);
3598  columnUpper[iColumn] = CoinMin(solution[iColumn],
3599    trueUpper[jNon]);
3600}
3601delete[] markNonlinear;
3602delete[] statusCheck;
3603delete[] savePi;
3604delete[] saveStatus;
3605// redo objective
3606CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns,
3607  objective);
3608primalDualCuts(rowsIn, 1, 1);
3609delete objective_;
3610delete[] rowsIn;
3611objective_ = trueObjective;
3612// redo values
3613setDblParam(ClpObjOffset, objectiveOffset);
3614objectiveValue_ += whichWay * offset;
3615for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3616  iColumn = listNonLinearColumn[jNon];
3617  columnLower[iColumn] = trueLower[jNon];
3618  columnUpper[iColumn] = trueUpper[jNon];
3619}
3620delete[] saveSolution;
3621delete[] safeSolution;
3622delete[] saveRowSolution;
3623for (iPass = 0; iPass < 3; iPass++)
3624  delete[] last[iPass];
3625delete[] trust;
3626delete[] trueUpper;
3627delete[] trueLower;
3628delete[] listNonLinearColumn;
3629delete[] changeRegion;
3630// temp
3631//setLogLevel(63);
3632return 0;
3633}
3634/* Primal algorithm for nonlinear constraints
3635   Using a semi-trust region approach as for pooling problem
3636   This is in because I have it lying around
3637
3638*/
3639int ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint **constraints,
3640  int numberPasses, double deltaTolerance)
3641{
3642  if (!numberConstraints) {
3643    // no nonlinear constraints - may be nonlinear objective
3644    return primalSLP(numberPasses, deltaTolerance);
3645  }
3646  // Are we minimizing or maximizing
3647  double whichWay = optimizationDirection();
3648  if (whichWay < 0.0)
3649    whichWay = -1.0;
3650  else if (whichWay > 0.0)
3651    whichWay = 1.0;
3652  // check all matrix for odd rows is empty
3653  int iConstraint;
3654  int numberBad = 0;
3655  CoinPackedMatrix *columnCopy = matrix();
3656  // Get a row copy in standard format
3657  CoinPackedMatrix copy;
3658  copy.reverseOrderedCopyOf(*columnCopy);
3659  // get matrix data pointers
3660  //const int * column = copy.getIndices();
3661  //const CoinBigIndex * rowStart = copy.getVectorStarts();
3662  const int *rowLength = copy.getVectorLengths();
3663  //const double * elementByRow = copy.getElements();
3664  int numberArtificials = 0;
3665  // We could use nonlinearcost to do segments - maybe later
3666#define SEGMENTS 3
3667  // Penalties may be adjusted by duals
3668  // Both these should be modified depending on problem
3669  // Possibly start with big bounds
3670  //double penalties[]={1.0e-3,1.0e7,1.0e9};
3671  double penalties[] = { 1.0e7, 1.0e8, 1.0e9 };
3672  double bounds[] = { 1.0e-2, 1.0e2, COIN_DBL_MAX };
3673  // see how many extra we need
3674  CoinBigIndex numberExtra = 0;
3675  for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) {
3676    ClpConstraint *constraint = constraints[iConstraint];
3677    int iRow = constraint->rowNumber();
3678    assert(iRow >= 0);
3679    int n = constraint->numberCoefficients() - rowLength[iRow];
3680    numberExtra += n;
3681    if (iRow >= numberRows_)
3682      numberBad++;
3683    else if (rowLength[iRow] && n)
3684      numberBad++;
3685    if (rowLower_[iRow] > -1.0e20)
3686      numberArtificials++;
3687    if (rowUpper_[iRow] < 1.0e20)
3688      numberArtificials++;
3689  }
3690  if (numberBad)
3691    return numberBad;
3692  ClpObjective *trueObjective = NULL;
3693  if (objective_->type() >= 2) {
3694    // Replace objective
3695    trueObjective = objective_;
3696    objective_ = new ClpLinearObjective(NULL, numberColumns_);
3697  }
3698  ClpSimplex newModel(*this);
3699  // we can put back true objective
3700  if (trueObjective) {
3701    // Replace objective
3702    delete objective_;
3703    objective_ = trueObjective;
3704  }
3705  int numberColumns2 = numberColumns_;
3706  int numberSmallGap = numberArtificials;
3707  if (numberArtificials) {
3708    numberArtificials *= SEGMENTS;
3709    numberColumns2 += numberArtificials;
3710    CoinBigIndex *addStarts = new CoinBigIndex[numberArtificials + 1];
3711    int *addRow = new int[numberArtificials];
3712    double *addElement = new double[numberArtificials];
3713    double *addUpper = new double[numberArtificials];
3714    addStarts[0] = 0;
3715    double *addCost = new double[numberArtificials];
3716    numberArtificials = 0;
3717    for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) {
3718      ClpConstraint *constraint = constraints[iConstraint];
3719      int iRow = constraint->rowNumber();
3720      if (rowLower_[iRow] > -1.0e20) {
3721        for (int k = 0; k < SEGMENTS; k++) {
3722          addRow[numberArtificials] = iRow;
3723          addElement[numberArtificials] = 1.0;
3724          addCost[numberArtificials] = penalties[k];
3725          addUpper[numberArtificials] = bounds[k];
3726          numberArtificials++;
3727          addStarts[numberArtificials] = numberArtificials;
3728        }
3729      }
3730      if (rowUpper_[iRow] < 1.0e20) {
3731        for (int k = 0; k < SEGMENTS; k++) {
3732          addRow[numberArtificials] = iRow;
3733          addElement[numberArtificials] = -1.0;
3734          addCost[numberArtificials] = penalties[k];
3735          addUpper[numberArtificials] = bounds[k];
3736          numberArtificials++;
3737          addStarts[numberArtificials] = numberArtificials;
3738        }
3739      }
3740    }
3741    newModel.addColumns(numberArtificials, NULL, addUpper, addCost,
3742      addStarts, addRow, addElement);
3743    delete[] addStarts;
3744    delete[] addRow;
3745    delete[] addElement;
3746    delete[] addUpper;
3747    delete[] addCost;
3748    //    newModel.primal(1);
3749  }
3750  // find nonlinear columns
3751  int *listNonLinearColumn = new int[numberColumns_ + numberSmallGap];
3752  char *mark = new char[numberColumns_];
3753  memset(mark, 0, numberColumns_);
3754  for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) {
3755    ClpConstraint *constraint = constraints[iConstraint];
3756    constraint->markNonlinear(mark);
3757  }
3758  if (trueObjective)
3759    trueObjective->markNonlinear(mark);
3760  int iColumn;
3761  int numberNonLinearColumns = 0;
3762  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3763    if (mark[iColumn])
3764      listNonLinearColumn[numberNonLinearColumns++] = iColumn;
3765  }
3766  // and small gap artificials
3767  for (iColumn = numberColumns_; iColumn < numberColumns2; iColumn += SEGMENTS) {
3768    listNonLinearColumn[numberNonLinearColumns++] = iColumn;
3769  }
3770  // build row copy of matrix with space for nonzeros
3771  // Get column copy
3772  columnCopy = newModel.matrix();
3773  copy.reverseOrderedCopyOf(*columnCopy);
3774  // get matrix data pointers
3775  const int *column = copy.getIndices();
3776  const CoinBigIndex *rowStart = copy.getVectorStarts();
3777  rowLength = copy.getVectorLengths();
3778  const double *elementByRow = copy.getElements();
3779  numberExtra += copy.getNumElements();
3780  CoinBigIndex *newStarts = new CoinBigIndex[numberRows_ + 1];
3781  int *newColumn = new int[numberExtra];
3782  double *newElement = new double[numberExtra];
3783  newStarts[0] = 0;
3784  int *backRow = new int[numberRows_];
3785  int iRow;
3786  for (iRow = 0; iRow < numberRows_; iRow++)
3787    backRow[iRow] = -1;
3788  for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) {
3789    ClpConstraint *constraint = constraints[iConstraint];
3790    iRow = constraint->rowNumber();
3791    backRow[iRow] = iConstraint;
3792  }
3793  numberExtra = 0;
3794  for (iRow = 0; iRow < numberRows_; iRow++) {
3795    if (backRow[iRow] < 0) {
3796      // copy normal
3797      for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow];
3798           j++) {
3799        newColumn[numberExtra] = column[j];
3800        newElement[numberExtra++] = elementByRow[j];
3801      }
3802    } else {
3803      ClpConstraint *constraint = constraints[backRow[iRow]];
3804      assert(iRow == constraint->rowNumber());
3805      int numberArtificials = 0;
3806      if (rowLower_[iRow] > -1.0e20)
3807        numberArtificials += SEGMENTS;
3808      if (rowUpper_[iRow] < 1.0e20)
3809        numberArtificials += SEGMENTS;
3810      if (numberArtificials == rowLength[iRow]) {
3811        // all possible
3812        memset(mark, 0, numberColumns_);
3813        constraint->markNonzero(mark);
3814        for (int k = 0; k < numberColumns_; k++) {
3815          if (mark[k]) {
3816            newColumn[numberExtra] = k;
3817            newElement[numberExtra++] = 1.0;
3818          }
3819        }
3820        // copy artificials
3821        for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow];
3822             j++) {
3823          newColumn[numberExtra] = column[j];
3824          newElement[numberExtra++] = elementByRow[j];
3825        }
3826      } else {
3827        // there already
3828        // copy
3829        for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow];
3830             j++) {
3831          newColumn[numberExtra] = column[j];
3832          assert(elementByRow[j]);
3833          newElement[numberExtra++] = elementByRow[j];
3834        }
3835      }
3836    }
3837    newStarts[iRow + 1] = numberExtra;
3838  }
3839  delete[] backRow;
3840  CoinPackedMatrix saveMatrix(false, numberColumns2, numberRows_,
3841    numberExtra, newElement, newColumn, newStarts, NULL, 0.0, 0.0);
3842  delete[] newStarts;
3843  delete[] newColumn;
3844  delete[] newElement;
3845  delete[] mark;
3846  // get feasible
3847  if (whichWay < 0.0) {
3848    newModel.setOptimizationDirection(1.0);
3849    double *objective = newModel.objective();
3850    for (int iColumn = 0; iColumn < numberColumns_; iColumn++)
3851      objective[iColumn] = -objective[iColumn];
3852  }
3853  newModel.primal(1);
3854  // still infeasible
3855  if (newModel.problemStatus() == 1) {
3856    delete[] listNonLinearColumn;
3857    return 0;
3858  } else if (newModel.problemStatus() == 2) {
3859    // unbounded - add bounds
3860    double *columnLower = newModel.columnLower();
3861    double *columnUpper = newModel.columnUpper();
3862    for (int i = 0; i < numberColumns_; i++) {
3863      columnLower[i] = CoinMax(-1.0e8, columnLower[i]);
3864      columnUpper[i] = CoinMin(1.0e8, columnUpper[i]);
3865    }
3866    newModel.primal(1);
3867  }
3868  int numberRows = newModel.numberRows();
3869  double *columnLower = newModel.columnLower();
3870  double *columnUpper = newModel.columnUpper();
3871  double *objective = newModel.objective();
3872  double *rowLower = newModel.rowLower();
3873  double *rowUpper = newModel.rowUpper();
3874  double *solution = newModel.primalColumnSolution();
3875  int jNon;
3876  int *last[3];
3877
3878  double *trust = new double[numberNonLinearColumns];
3879  double *trueLower = new double[numberNonLinearColumns];
3880  double *trueUpper = new double[numberNonLinearColumns];
3881  double objectiveOffset;
3882  double objectiveOffset2;
3883  getDblParam(ClpObjOffset, objectiveOffset);
3884  objectiveOffset *= whichWay;
3885  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3886    iColumn = listNonLinearColumn[jNon];
3887    double upper = columnUpper[iColumn];
3888    double lower = columnLower[iColumn];
3889    if (solution[iColumn] < lower)
3890      solution[iColumn] = lower;
3891    else if (solution[iColumn] > upper)
3892      solution[iColumn] = upper;
3893#if 0
3894          double large = CoinMax(1000.0, 10.0 * fabs(solution[iColumn]));
3895          if (upper > 1.0e10)
3896               upper = solution[iColumn] + large;
3897          if (lower < -1.0e10)
3898               lower = solution[iColumn] - large;
3899#else
3900    upper = solution[iColumn] + 0.5;
3901    lower = solution[iColumn] - 0.5;
3902#endif
3903    //columnUpper[iColumn]=upper;
3904    trust[jNon] = 0.05 * (1.0 + upper - lower);
3905    trueLower[jNon] = columnLower[iColumn];
3906    //trueUpper[jNon]=upper;
3907    trueUpper[jNon] = columnUpper[iColumn];
3908  }
3909  bool tryFix = false;
3910  int iPass;
3911  double lastObjective = 1.0e31;
3912  double lastGoodObjective = 1.0e31;
3913  double *bestSolution = NULL;
3914  double *saveSolution = new double[numberColumns2 + numberRows];
3915  char *saveStatus = new char[numberColumns2 + numberRows];
3916  double targetDrop = 1.0e31;
3917  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
3918  for (iPass = 0; iPass < 3; iPass++) {
3919    last[iPass] = new int[numberNonLinearColumns];
3920    for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
3921      last[iPass][jNon] = 0;
3922  }
3923  int numberZeroPasses = 0;
3924  bool zeroTargetDrop = false;
3925  double *gradient = new double[numberColumns_];
3926  bool goneFeasible = false;
3927  // keep sum of artificials
3928#define KEEP_SUM 5
3929  double sumArt[KEEP_SUM];
3930  for (jNon = 0; jNon < KEEP_SUM; jNon++)
3931    sumArt[jNon] = COIN_DBL_MAX;
3932#define SMALL_FIX 0.0
3933  for (iPass = 0; iPass < numberPasses; iPass++) {
3934    objectiveOffset2 = objectiveOffset;
3935    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
3936      iColumn = listNonLinearColumn[jNon];
3937      if (solution[iColumn] < trueLower[jNon])
3938        solution[iColumn] = trueLower[jNon];
3939      else if (solution[iColumn] > trueUpper[jNon])
3940        solution[iColumn] = trueUpper[jNon];
3941      columnLower[iColumn] = CoinMax(solution[iColumn]
3942          - trust[jNon],
3943        trueLower[jNon]);
3944      if (!trueLower[jNon] && columnLower[iColumn] < SMALL_FIX)
3945        columnLower[iColumn] = SMALL_FIX;
3946      columnUpper[iColumn] = CoinMin(solution[iColumn]
3947          + trust[jNon],
3948        trueUpper[jNon]);
3949      if (!trueLower[jNon])
3950        columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
3951          columnLower[iColumn] + SMALL_FIX);
3952      if (!trueLower[jNon] && tryFix && columnLower[iColumn] == SMALL_FIX && columnUpper[iColumn] < 3.0 * SMALL_FIX) {
3953        columnLower[iColumn] = 0.0;
3954        solution[iColumn] = 0.0;
3955        columnUpper[iColumn] = 0.0;
3956        printf("fixing %d\n", iColumn);
3957      }
3958    }
3959    // redo matrix
3960    double offset;
3961    CoinPackedMatrix newMatrix(saveMatrix);
3962    // get matrix data pointers
3963    column = newMatrix.getIndices();
3964    rowStart = newMatrix.getVectorStarts();
3965    rowLength = newMatrix.getVectorLengths();
3966    // make sure x updated
3967    if (numberConstraints)
3968      constraints[0]->newXValues();
3969    else
3970      trueObjective->newXValues();
3971    double *changeableElement = newMatrix.getMutableElements();
3972    if (trueObjective) {
3973      CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns_,
3974        objective);
3975    } else {
3976      CoinMemcpyN(objective_->gradient(this, solution, offset, true, 2), numberColumns_,
3977        objective);
3978    }
3979    if (whichWay < 0.0) {
3980      for (int iColumn = 0; iColumn < numberColumns_; iColumn++)
3981        objective[iColumn] = -objective[iColumn];
3982    }
3983    for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) {
3984      ClpConstraint *constraint = constraints[iConstraint];
3985      int iRow = constraint->rowNumber();
3986      double functionValue;
3987#ifndef NDEBUG
3988      int numberErrors =
3989#endif
3990        constraint->gradient(&newModel, solution, gradient, functionValue, offset);
3991      assert(!numberErrors);
3992      // double dualValue = newModel.dualRowSolution()[iRow];
3993      int numberCoefficients = constraint->numberCoefficients();
3994      for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + numberCoefficients; j++) {
3995        int iColumn = column[j];
3996        changeableElement[j] = gradient[iColumn];
3997        //objective[iColumn] -= dualValue*gradient[iColumn];
3998        gradient[iColumn] = 0.0;
3999      }
4000      for (int k = 0; k < numberColumns_; k++)
4001        assert(!gradient[k]);
4002      if (rowLower_[iRow] > -1.0e20)
4003        rowLower[iRow] = rowLower_[iRow] - offset;
4004      if (rowUpper_[iRow] < 1.0e20)
4005        rowUpper[iRow] = rowUpper_[iRow] - offset;
4006    }
4007    // Replace matrix
4008    // Get a column copy in standard format
4009    CoinPackedMatrix *columnCopy = new CoinPackedMatrix();
4010    ;
4011    columnCopy->reverseOrderedCopyOf(newMatrix);
4012    newModel.replaceMatrix(columnCopy, true);
4013    // solve
4014    newModel.primal(1);
4015    if (newModel.status() == 1) {
4016      // Infeasible!
4017      newModel.allSlackBasis();
4018      newModel.primal();
4019      newModel.writeMps("infeas.mps");
4020      assert(!newModel.status());
4021    }
4022    double sumInfeas = 0.0;
4023    int numberInfeas = 0;
4024    for (iColumn = numberColumns_; iColumn < numberColumns2; iColumn++) {
4025      if (solution[iColumn] > 1.0e-8) {
4026        numberInfeas++;
4027        sumInfeas += solution[iColumn];
4028      }
4029    }
4030    printf("%d artificial infeasibilities - summing to %g\n",
4031      numberInfeas, sumInfeas);
4032    for (jNon = 0; jNon < KEEP_SUM - 1; jNon++)
4033      sumArt[jNon] = sumArt[jNon + 1];
4034    sumArt[KEEP_SUM - 1] = sumInfeas;
4035    if (sumInfeas > 0.01 && sumInfeas * 1.1 > sumArt[0] && penalties[1] < 1.0e7) {
4036      // not doing very well - increase - be more sophisticated later
4037      lastObjective = COIN_DBL_MAX;
4038      for (jNon = 0; jNon < KEEP_SUM; jNon++)
4039        sumArt[jNon] = COIN_DBL_MAX;
4040      //for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) {
4041      //objective[iColumn+1] *= 1.5;
4042      //}
4043      penalties[1] *= 1.5;
4044      for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
4045        if (trust[jNon] > 0.1)
4046          trust[jNon] *= 2.0;
4047        else
4048          trust[jNon] = 0.1;
4049    }
4050    if (sumInfeas < 0.001 && !goneFeasible) {
4051      goneFeasible = true;
4052      penalties[0] = 1.0e-3;
4053      penalties[1] = 1.0e6;
4054      printf("Got feasible\n");
4055    }
4056    double infValue = 0.0;
4057    double objValue = 0.0;
4058    // make sure x updated
4059    if (numberConstraints)
4060      constraints[0]->newXValues();
4061    else
4062      trueObjective->newXValues();
4063    if (trueObjective) {
4064      objValue = trueObjective->objectiveValue(this, solution);
4065      printf("objective offset %g\n", offset);
4066      objectiveOffset2 = objectiveOffset + offset; // ? sign
4067      newModel.setObjectiveOffset(objectiveOffset2);
4068    } else {
4069      objValue = objective_->objectiveValue(this, solution);
4070    }
4071    objValue *= whichWay;
4072    double infPenalty = 0.0;
4073    // This penalty is for target drop
4074    double infPenalty2 = 0.0;
4075    //const int * row = columnCopy->getIndices();
4076    //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4077    //const int * columnLength = columnCopy->getVectorLengths();
4078    //const double * element = columnCopy->getElements();
4079    double *cost = newModel.objective();
4080    column = newMatrix.getIndices();
4081    rowStart = newMatrix.getVectorStarts();
4082    rowLength = newMatrix.getVectorLengths();
4083    elementByRow = newMatrix.getElements();
4084    int jColumn = numberColumns_;
4085    double objectiveAdjustment = 0.0;
4086    for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) {
4087      ClpConstraint *constraint = constraints[iConstraint];
4088      int iRow = constraint->rowNumber();
4089      double functionValue = constraint->functionValue(this, solution);
4090      double dualValue = newModel.dualRowSolution()[iRow];
4091      if (numberConstraints < -50)
4092        printf("For row %d current value is %g (row activity %g) , dual is %g\n", iRow, functionValue,
4093          newModel.primalRowSolution()[iRow],
4094          dualValue);
4095      double movement = newModel.primalRowSolution()[iRow] + constraint->offset();
4096      movement = fabs((movement - functionValue) * dualValue);
4097      infPenalty2 += movement;
4098      double sumOfActivities = 0.0;
4099      for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
4100        int iColumn = column[j];
4101        sumOfActivities += fabs(solution[iColumn] * elementByRow[j]);
4102      }
4103      if (rowLower_[iRow] > -1.0e20) {
4104        if (functionValue < rowLower_[iRow] - 1.0e-5) {
4105          double infeasibility = rowLower_[iRow] - functionValue;
4106          double thisPenalty = 0.0;
4107          infValue += infeasibility;
4108          int k;
4109          assert(dualValue >= -1.0e-5);
4110          dualValue = CoinMax(dualValue, 0.0);
4111          for (k = 0; k < SEGMENTS; k++) {
4112            if (infeasibility <= 0)
4113              break;
4114            double thisPart = CoinMin(infeasibility, bounds[k]);
4115            thisPenalty += thisPart * cost[jColumn + k];
4116            infeasibility -= thisPart;
4117          }
4118          infeasibility = functionValue - rowUpper_[iRow];
4119          double newPenalty = 0.0;
4120          for (k = 0; k < SEGMENTS; k++) {
4121            double thisPart = CoinMin(infeasibility, bounds[k]);
4122            cost[jColumn + k] = CoinMax(penalties[k], dualValue + 1.0e-3);
4123            newPenalty += thisPart * cost[jColumn + k];
4124            infeasibility -= thisPart;
4125          }
4126          infPenalty += thisPenalty;
4127          objectiveAdjustment += CoinMax(0.0, newPenalty - thisPenalty);
4128        }
4129        jColumn += SEGMENTS;
4130      }
4131      if (rowUpper_[iRow] < 1.0e20) {
4132        if (functionValue > rowUpper_[iRow] + 1.0e-5) {
4133          double infeasibility = functionValue - rowUpper_[iRow];
4134          double thisPenalty = 0.0;
4135          infValue += infeasibility;
4136          int k;
4137          dualValue = -dualValue;
4138          assert(dualValue >= -1.0e-5);
4139          dualValue = CoinMax(dualValue, 0.0);
4140          for (k = 0; k < SEGMENTS; k++) {
4141            if (infeasibility <= 0)
4142              break;
4143            double thisPart = CoinMin(infeasibility, bounds[k]);
4144            thisPenalty += thisPart * cost[jColumn + k];
4145            infeasibility -= thisPart;
4146          }
4147          infeasibility = functionValue - rowUpper_[iRow];
4148          double newPenalty = 0.0;
4149          for (k = 0; k < SEGMENTS; k++) {
4150            double thisPart = CoinMin(infeasibility, bounds[k]);
4151            cost[jColumn + k] = CoinMax(penalties[k], dualValue + 1.0e-3);
4152            newPenalty += thisPart * cost[jColumn + k];
4153            infeasibility -= thisPart;
4154          }
4155          infPenalty += thisPenalty;
4156          objectiveAdjustment += CoinMax(0.0, newPenalty - thisPenalty);
4157        }
4158        jColumn += SEGMENTS;
4159      }
4160    }
4161    // adjust last objective value
4162    lastObjective += objectiveAdjustment;
4163    if (infValue)
4164      printf("Sum infeasibilities %g - penalty %g ", infValue, infPenalty);
4165    if (objectiveOffset2)
4166      printf("offset2 %g ", objectiveOffset2);
4167    objValue -= objectiveOffset2;
4168    printf("True objective %g or maybe %g (with penalty %g) -pen2 %g %g\n", objValue,
4169      objValue + objectiveOffset2, objValue + objectiveOffset2 + infPenalty, infPenalty2, penalties[1]);
4170    double useObjValue = objValue + objectiveOffset2 + infPenalty;
4171    objValue += infPenalty + infPenalty2;
4172    objValue = useObjValue;
4173    if (iPass) {
4174      double drop = lastObjective - objValue;
4175      std::cout << "True drop was " << drop << std::endl;
4176      if (drop < -0.05 * fabs(objValue) - 1.0e-4) {
4177        // pretty bad - go back and halve
4178        CoinMemcpyN(saveSolution, numberColumns2, solution);
4179        CoinMemcpyN(saveSolution + numberColumns2,
4180          numberRows, newModel.primalRowSolution());
4181        CoinMemcpyN(reinterpret_cast< unsigned char * >(saveStatus),
4182          numberColumns2 + numberRows, newModel.statusArray());
4183        for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
4184          if (trust[jNon] > 0.1)
4185            trust[jNon] *= 0.5;
4186          else
4187            trust[jNon] *= 0.9;
4188
4189        printf("** Halving trust\n");
4190        objValue = lastObjective;
4191        continue;
4192      } else if ((iPass % 25) == -1) {
4193        for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
4194          trust[jNon] *= 2.0;
4195        printf("** Doubling trust\n");
4196      }
4197      int *temp = last[2];
4198      last[2] = last[1];
4199      last[1] = last[0];
4200      last[0] = temp;
4201      for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
4202        iColumn = listNonLinearColumn[jNon];
4203        double change = solution[iColumn] - saveSolution[iColumn];
4204        if (change < -1.0e-5) {
4205          if (fabs(change + trust[jNon]) < 1.0e-5)
4206            temp[jNon] = -1;
4207          else
4208            temp[jNon] = -2;
4209        } else if (change > 1.0e-5) {
4210          if (fabs(change - trust[jNon]) < 1.0e-5)
4211            temp[jNon] = 1;
4212          else
4213            temp[jNon] = 2;
4214        } else {
4215          temp[jNon] = 0;
4216        }
4217      }
4218      double maxDelta = 0.0;
4219      double smallestTrust = 1.0e31;
4220      double smallestNonLinearGap = 1.0e31;
4221      double smallestGap = 1.0e31;
4222      bool increasing = false;
4223      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4224        double gap = columnUpper[iColumn] - columnLower[iColumn];
4225        assert(gap >= 0.0);
4226        if (gap)
4227          smallestGap = CoinMin(smallestGap, gap);
4228      }
4229      for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
4230        iColumn = listNonLinearColumn[jNon];
4231        double gap = columnUpper[iColumn] - columnLower[iColumn];
4232        assert(gap >= 0.0);
4233        if (gap) {
4234          smallestNonLinearGap = CoinMin(smallestNonLinearGap, gap);
4235          if (gap < 1.0e-7 && iPass == 1) {
4236            printf("Small gap %d %d %g %g %g\n",
4237              jNon, iColumn, columnLower[iColumn], columnUpper[iColumn],
4238              gap);
4239            //trueUpper[jNon]=trueLower[jNon];
4240            //columnUpper[iColumn]=columnLower[iColumn];
4241          }
4242        }
4243        maxDelta = CoinMax(maxDelta,
4244          fabs(solution[iColumn] - saveSolution[iColumn]));
4245        if (last[0][jNon] * last[1][jNon] < 0) {
4246          // halve
4247          if (trust[jNon] > 1.0)
4248            trust[jNon] *= 0.5;
4249          else
4250            trust[jNon] *= 0.7;
4251        } else {
4252          // ? only increase if +=1 ?
4253          if (last[0][jNon] == last[1][jNon] && last[0][jNon] == last[2][jNon] && last[0][jNon]) {
4254            trust[jNon] *= 1.8;
4255            increasing = true;
4256          }
4257        }
4258        smallestTrust = CoinMin(smallestTrust, trust[jNon]);
4259      }
4260      std::cout << "largest delta is " << maxDelta
4261                << ", smallest trust is " << smallestTrust
4262                << ", smallest gap is " << smallestGap
4263                << ", smallest nonlinear gap is " << smallestNonLinearGap
4264                << std::endl;
4265      if (iPass > 200) {
4266        //double useObjValue = objValue+objectiveOffset2+infPenalty;
4267        if (useObjValue + 1.0e-4 > lastGoodObjective && iPass > 250) {
4268          std::cout << "Exiting as objective not changing much" << std::endl;
4269          break;
4270        } else if (useObjValue < lastGoodObjective) {
4271          lastGoodObjective = useObjValue;
4272          if (!bestSolution)
4273            bestSolution = new double[numberColumns2];
4274          CoinMemcpyN(solution, numberColumns2, bestSolution);
4275        }
4276      }
4277      if (maxDelta < deltaTolerance && !increasing && iPass > 100) {
4278        numberZeroPasses++;
4279        if (numberZeroPasses == 3) {
4280          if (tryFix) {
4281            std::cout << "Exiting" << std::endl;
4282            break;
4283          } else {
4284            tryFix = true;
4285            for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
4286              trust[jNon] = CoinMax(trust[jNon], 1.0e-1);
4287            numberZeroPasses = 0;
4288          }
4289        }
4290      } else {
4291        numberZeroPasses = 0;
4292      }
4293    }
4294    CoinMemcpyN(solution, numberColumns2, saveSolution);
4295    CoinMemcpyN(newModel.primalRowSolution(),
4296      numberRows, saveSolution + numberColumns2);
4297    CoinMemcpyN(newModel.statusArray(),
4298      numberColumns2 + numberRows,
4299      reinterpret_cast< unsigned char * >(saveStatus));
4300
4301    targetDrop = infPenalty + infPenalty2;
4302    if (iPass) {
4303      // get reduced costs
4304      const double *pi = newModel.dualRowSolution();
4305      newModel.matrix()->transposeTimes(pi,
4306        newModel.dualColumnSolution());
4307      double *r = newModel.dualColumnSolution();
4308      for (iColumn = 0; iColumn < numberColumns_; iColumn++)
4309        r[iColumn] = objective[iColumn] - r[iColumn];
4310      for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
4311        iColumn = listNonLinearColumn[jNon];
4312        double dj = r[iColumn];
4313        if (dj < -1.0e-6) {
4314          double drop = -dj * (columnUpper[iColumn] - solution[iColumn]);
4315          //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
4316          //double drop2 = -dj*(upper-solution[iColumn]);
4317#if 0
4318                         if (drop > 1.0e8 || drop2 > 100.0 * drop || (drop > 1.0e-2 && iPass > 100))
4319                              printf("Big drop %d %g %g %g %g T %g\n",
4320                                     iColumn, columnLower[iColumn], solution[iColumn],
4321                                     columnUpper[iColumn], dj, trueUpper[jNon]);
4322#endif
4323          targetDrop += drop;
4324          if (dj < -1.0e-1 && trust[jNon] < 1.0e-3
4325            && trueUpper[jNon] - solution[iColumn] > 1.0e-2) {
4326            trust[jNon] *= 1.5;
4327            //printf("Increasing trust on %d to %g\n",
4328            //     iColumn,trust[jNon]);
4329          }
4330        } else if (dj > 1.0e-6) {
4331          double drop = -dj * (columnLower[iColumn] - solution[iColumn]);
4332          //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
4333          //double drop2 = -dj*(lower-solution[iColumn]);
4334#if 0
4335                         if (drop > 1.0e8 || drop2 > 100.0 * drop || (drop > 1.0e-2))
4336                              printf("Big drop %d %g %g %g %g T %g\n",
4337                                     iColumn, columnLower[iColumn], solution[iColumn],
4338                                     columnUpper[iColumn], dj, trueLower[jNon]);
4339#endif
4340          targetDrop += drop;
4341          if (dj > 1.0e-1 && trust[jNon] < 1.0e-3
4342            && solution[iColumn] - trueLower[jNon] > 1.0e-2) {
4343            trust[jNon] *= 1.5;
4344            printf("Increasing trust on %d to %g\n",
4345              iColumn, trust[jNon]);
4346          }
4347        }
4348      }
4349    }
4350    std::cout << "Pass - " << iPass
4351              << ", target drop is " << targetDrop
4352              << std::endl;
4353    if (iPass > 1 && targetDrop < 1.0e-5 && zeroTargetDrop)
4354      break;
4355    if (iPass > 1 && targetDrop < 1.0e-5)
4356      zeroTargetDrop = true;
4357    else
4358      zeroTargetDrop = false;
4359    //if (iPass==5)
4360    //newModel.setLogLevel(63);
4361    lastObjective = objValue;
4362    // take out when ClpPackedMatrix changed
4363    //newModel.scaling(false);
4364#if 0
4365          CoinMpsIO writer;
4366          writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
4367                            newModel.getColLower(), newModel.getColUpper(),
4368                            newModel.getObjCoefficients(),
4369                            (const char*) 0 /*integrality*/,
4370                            newModel.getRowLower(), newModel.getRowUpper(),
4371                            NULL, NULL);
4372          writer.writeMps("xx.mps");
4373#endif
4374  }
4375  delete[] saveSolution;
4376  delete[] saveStatus;
4377  for (iPass = 0; iPass < 3; iPass++)
4378    delete[] last[iPass];
4379  for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
4380    iColumn = listNonLinearColumn[jNon];
4381    columnLower[iColumn] = trueLower[jNon];
4382    columnUpper[iColumn] = trueUpper[jNon];
4383  }
4384  delete[] trust;
4385  delete[] trueUpper;
4386  delete[] trueLower;
4387  objectiveValue_ = newModel.objectiveValue();
4388  if (bestSolution) {
4389    CoinMemcpyN(bestSolution, numberColumns2, solution);
4390    delete[] bestSolution;
4391    printf("restoring objective of %g\n", lastGoodObjective);
4392    objectiveValue_ = lastGoodObjective;
4393  }
4394  // Simplest way to get true row activity ?
4395  double *rowActivity = newModel.primalRowSolution();
4396  for (iRow = 0; iRow < numberRows; iRow++) {
4397    double difference;
4398    if (fabs(rowLower_[iRow]) < fabs(rowUpper_[iRow]))
4399      difference = rowLower_[iRow] - rowLower[iRow];
4400    else
4401      difference = rowUpper_[iRow] - rowUpper[iRow];
4402    rowLower[iRow] = rowLower_[iRow];
4403    rowUpper[iRow] = rowUpper_[iRow];
4404    if (difference) {
4405      if (numberRows < 50)
4406        printf("For row %d activity changes from %g to %g\n",
4407          iRow, rowActivity[iRow], rowActivity[iRow] + difference);
4408      rowActivity[iRow] += difference;
4409    }
4410  }
4411  delete[] listNonLinearColumn;
4412  delete[] gradient;
4413  printf("solution still in newModel - do objective etc!\n");
4414  numberIterations_ = newModel.numberIterations();
4415  problemStatus_ = newModel.problemStatus();
4416  secondaryStatus_ = newModel.secondaryStatus();
4417  CoinMemcpyN(newModel.primalColumnSolution(), numberColumns_, columnActivity_);
4418  // should do status region
4419  CoinZeroN(rowActivity_, numberRows_);
4420  matrix_->times(1.0, columnActivity_, rowActivity_);
4421  return 0;
4422}
4423
4424/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
4425*/
Note: See TracBrowser for help on using the repository browser.