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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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