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

Last change on this file since 1924 was 1924, checked in by forrest, 7 years ago

changes for miqp, parameters and presolve

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