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

Last change on this file since 1929 was 1929, checked in by stefan, 7 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

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