source: stable/1.16/Clp/src/ClpSimplexNonlinear.cpp @ 2128

Last change on this file since 2128 was 2128, checked in by forrest, 4 years ago

out mac warnings

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