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

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

try and fix event handler

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