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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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