# source:trunk/ClpSimplexPrimal.cpp@284

Last change on this file since 284 was 266, checked in by forrest, 16 years ago

Trying to improve cycling and more barrier stuff

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 72.6 KB
Line
3
4
5/* Notes on implementation of primal simplex algorithm.
6
7   When primal feasible(A):
8
9   If dual feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to enter basis from a bound (B).  We now need to find an
11   outgoing variable which will leave problem primal feasible so we get
12   the column of the tableau corresponding to the incoming variable
13   (with the correct sign depending if variable will go up or down).
14
15   We now perform a ratio test to determine which outgoing variable will
16   preserve primal feasibility (C).  If no variable found then problem
17   is unbounded (in primal sense).  If there is a variable, we then
18   perform pivot and repeat.  Trivial?
19
20   -------------------------------------------
21
22   A) How do we get primal feasible?  All variables have fake costs
23   outside their feasible region so it is trivial to declare problem
24   feasible.  OSL did not have a phase 1/phase 2 approach but
25   instead effectively put an extra cost on infeasible basic variables
26   I am taking the same approach here, although it is generalized
27   to allow for non-linear costs and dual information.
28
29   In OSL, this weight was changed heuristically, here at present
30   it is only increased if problem looks finished.  If problem is
31   feasible I check for unboundedness.  If not unbounded we
32   could play with going into dual.  As long as weights increase
33   any algorithm would be finite.
34
35   B) Which incoming variable to choose is a virtual base class.
36   For difficult problems steepest edge is preferred while for
37   very easy (large) problems we will need partial scan.
38
39   C) Sounds easy, but this is hardest part of algorithm.
40      1) Instead of stopping at first choice, we may be able
41      to allow that variable to go through bound and if objective
42      still improving choose again.  These mini iterations can
43      increase speed by orders of magnitude but we may need to
44      go to more of a bucket choice of variable rather than looking
45      at them one by one (for speed).
46      2) Accuracy.  Basic infeasibilities may be less than
47      tolerance.  Pivoting on these makes objective go backwards.
48      OSL modified cost so a zero move was made, Gill et al
49      modified so a strictly positive move was made.
50      The two problems are that re-factorizations can
51      change rinfeasibilities above and below tolerances and that when
52      finished we need to reset costs and try again.
53      3) Degeneracy.  Gill et al helps but may not be enough.  We
54      may need more.  Also it can improve speed a lot if we perturb
55      the rhs and bounds significantly.
56
57  References:
58     Forrest and Goldfarb, Steepest-edge simplex algorithms for
59       linear programming - Mathematical Programming 1992
60     Forrest and Tomlin, Implementing the simplex method for
61       the Optimization Subroutine Library - IBM Systems Journal 1992
62     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
63       Procedure for Linear and Nonlinear Programming SOL report 1988
64
65
66  TODO:
67
68  a) Better recovery procedures.  At present I never check on forward
70     re-factorization frequency, but this is only on singular
71     factorizations.
72  b) Fast methods for large easy problems (and also the option for
73     the code to automatically choose which method).
74  c) We need to be able to stop in various ways for OSI - this
75     is fairly easy.
76
77 */
78
79
80#include "CoinPragma.hpp"
81
82#include <math.h>
83
84#include "CoinHelperFunctions.hpp"
85#include "ClpSimplexPrimal.hpp"
86#include "ClpFactorization.hpp"
87#include "ClpNonLinearCost.hpp"
88#include "CoinPackedMatrix.hpp"
89#include "CoinIndexedVector.hpp"
90#include "CoinWarmStartBasis.hpp"
91#include "ClpPrimalColumnPivot.hpp"
92#include "ClpMessage.hpp"
93#include <cfloat>
94#include <cassert>
95#include <string>
96#include <stdio.h>
97#include <iostream>
98// primal
99int ClpSimplexPrimal::primal (int ifValuesPass )
100{
101
102  /*
103      Method
104
105     It tries to be a single phase approach with a weight of 1.0 being
106     given to getting optimal and a weight of infeasibilityCost_ being
107     given to getting primal feasible.  In this version I have tried to
108     be clever in a stupid way.  The idea of fake bounds in dual
109     seems to work so the primal analogue would be that of getting
110     bounds on reduced costs (by a presolve approach) and using
111     these for being above or below feasible region.  I decided to waste
112     memory and keep these explicitly.  This allows for non-linear
113     costs!
114
115     The code is designed to take advantage of sparsity so arrays are
116     seldom zeroed out from scratch or gone over in their entirety.
117     The only exception is a full scan to find incoming variable for
118     Dantzig row choice.  For steepest edge we keep an updated list
119     of dual infeasibilities (actually squares).
120     On easy problems we don't need full scan - just
121     pick first reasonable.
122
123     One problem is how to tackle degeneracy and accuracy.  At present
124     I am using the modification of costs which I put in OSL and which was
125     extended by Gill et al.  I am still not sure of the exact details.
126
127     The flow of primal is three while loops as follows:
128
129     while (not finished) {
130
131       while (not clean solution) {
132
133          Factorize and/or clean up solution by changing bounds so
134          primal feasible.  If looks finished check fake primal bounds.
135          Repeat until status is iterating (-1) or finished (0,1,2)
136
137       }
138
139       while (status==-1) {
140
141         Iterate until no pivot in or out or time to re-factorize.
142
143         Flow is:
144
145         choose pivot column (incoming variable).  if none then
146         we are primal feasible so looks as if done but we need to
147         break and check bounds etc.
148
149         Get pivot column in tableau
150
151         Choose outgoing row.  If we don't find one then we look
152         primal unbounded so break and check bounds etc.  (Also the
153         pivot tolerance is larger after any iterations so that may be
154         reason)
155
156         If we do find outgoing row, we may have to adjust costs to
157         keep going forwards (anti-degeneracy).  Check pivot will be stable
158         and if unstable throw away iteration and break to re-factorize.
159         If minor error re-factorize after iteration.
160
161         Update everything (this may involve changing bounds on
162         variables to stay primal feasible.
163
164       }
165
166     }
167
168     At present we never check we are going forwards.  I overdid that in
169     OSL so will try and make a last resort.
170
171     Needs partial scan pivot in option.
172
173     May need other anti-degeneracy measures, especially if we try and use
174     loose tolerances as a way to solve in fewer iterations.
175
176     I like idea of dynamic scaling.  This gives opportunity to decouple
177     different implications of scaling for accuracy, iteration count and
178     feasibility tolerance.
179
180  */
181
182  algorithm_ = +1;
183  //specialOptions_ |= 4;
184
185  // save data
186  ClpDataSave data = saveData();
187
188  // Save so can see if doing after dual
189  int initialStatus=problemStatus_;
190  // initialize - maybe values pass and algorithm_ is +1
191  if (!startup(ifValuesPass)) {
192
193    // Set average theta
194    nonLinearCost_->setAverageTheta(1.0e3);
195    int lastCleaned=0; // last time objective or bounds cleaned up
196
197    // Say no pivot has occurred (for steepest edge and updates)
198    pivotRow_=-2;
199
200    // This says whether to restore things etc
201    int factorType=0;
202    if (problemStatus_<0&&perturbation_<100) {
203      perturb(0);
204      // Can't get here if values pass
205      assert (!ifValuesPass);
206      gutsOfSolution(NULL,NULL);
207      if (handler_->logLevel()>2) {
208        handler_->message(CLP_SIMPLEX_STATUS,messages_)
209          <<numberIterations_<<objectiveValue();
210        handler_->printing(sumPrimalInfeasibilities_>0.0)
211          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
212        handler_->printing(sumDualInfeasibilities_>0.0)
213          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
214        handler_->printing(numberDualInfeasibilitiesWithoutFree_
215                           <numberDualInfeasibilities_)
216                             <<numberDualInfeasibilitiesWithoutFree_;
217        handler_->message()<<CoinMessageEol;
218      }
219    }
220    ClpSimplex * saveModel=NULL;
221    int stopSprint=-1;
222    int sprintPass=0;
223    int reasonableSprintIteration=0;
224    int lastSprintIteration=0;
225    double lastObjectiveValue=COIN_DBL_MAX;
226    // Start check for cycles
227    progress_->startCheck();
228    /*
229      Status of problem:
230      0 - optimal
231      1 - infeasible
232      2 - unbounded
233      -1 - iterating
234      -2 - factorization wanted
235      -3 - redo checking without factorization
236      -4 - looks infeasible
237      -5 - looks unbounded
238    */
239    while (problemStatus_<0) {
240      int iRow,iColumn;
241      // clear
242      for (iRow=0;iRow<4;iRow++) {
243        rowArray_[iRow]->clear();
244      }
245
246      for (iColumn=0;iColumn<2;iColumn++) {
247        columnArray_[iColumn]->clear();
248      }
249
250      // give matrix (and model costs and bounds a chance to be
251      // refreshed (normally null)
252      matrix_->refresh(this);
253      // If getting nowhere - why not give it a kick
254#if 1
255      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
256          &&initialStatus!=10)
257        perturb(1);
258#endif
259      // If we have done no iterations - special
260      if (lastGoodIteration_==numberIterations_&&factorType)
261        factorType=3;
262      if (saveModel) {
263        // Doing sprint
264        if (sequenceIn_<0||numberIterations_>=stopSprint) {
265          problemStatus_=-1;
266          originalModel(saveModel);
267          saveModel=NULL;
268          if (sequenceIn_<0&&numberIterations_<reasonableSprintIteration&&
269              sprintPass>100)
270            primalColumnPivot_->switchOffSprint();
271          //lastSprintIteration=numberIterations_;
272          printf("End small model\n");
273        }
274      }
275
276      // may factorize, checks if problem finished
277      statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
278      // See if sprint says redo beacuse of problems
279      if (numberDualInfeasibilities_==-776) {
280        // Need new set of variables
281        problemStatus_=-1;
282        originalModel(saveModel);
283        saveModel=NULL;
284        //lastSprintIteration=numberIterations_;
285        printf("End small model after\n");
286        statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
287      }
288      int numberSprintIterations=0;
289      int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
290      if (problemStatus_==777) {
291        // problems so do one pass with normal
292        problemStatus_=-1;
293        originalModel(saveModel);
294        saveModel=NULL;
295        // Skip factorization
296        //statusOfProblemInPrimal(lastCleaned,factorType,progress_,false,saveModel);
297        statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
298      } else if (problemStatus_<0&&!saveModel&&numberSprintColumns&&firstFree_<0) {
299        int numberSort=0;
300        int numberFixed=0;
301        int numberBasic=0;
302        reasonableSprintIteration = numberIterations_ + 100;
303        int * whichColumns = new int[numberColumns_];
304        double * weight = new double[numberColumns_];
305        int numberNegative=0;
306        double sumNegative = 0.0;
307        // now massage weight so all basic in plus good djs
308        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
309          double dj = dj_[iColumn];
310          switch(getColumnStatus(iColumn)) {
311
312          case basic:
313            dj = -1.0e50;
314            numberBasic++;
315            break;
316          case atUpperBound:
317            dj = -dj;
318            break;
319          case isFixed:
320            dj=1.0e50;
321            numberFixed++;
322            break;
323          case atLowerBound:
324            dj = dj;
325            break;
326          case isFree:
327            dj = -100.0*fabs(dj);
328              break;
329          case superBasic:
330            dj = -100.0*fabs(dj);
331            break;
332          }
333          if (dj<-dualTolerance_&&dj>-1.0e50) {
334            numberNegative++;
335            sumNegative -= dj;
336          }
337          weight[iColumn]=dj;
338          whichColumns[iColumn] = iColumn;
339        }
340        handler_->message(CLP_SPRINT,messages_)
341          <<sprintPass<<numberIterations_-lastSprintIteration<<objectiveValue()<<sumNegative
342          <<numberNegative
343          <<CoinMessageEol;
344        sprintPass++;
345        lastSprintIteration=numberIterations_;
346        if (objectiveValue()*optimizationDirection_>lastObjectiveValue-1.0e-7&&sprintPass>5) {
347          // switch off
348          printf("Switching off sprint\n");
349          primalColumnPivot_->switchOffSprint();
350        } else {
351          lastObjectiveValue = objectiveValue()*optimizationDirection_;
352          // sort
353          CoinSort_2(weight,weight+numberColumns_,whichColumns);
354          numberSort = min(numberColumns_-numberFixed,numberBasic+numberSprintColumns);
355          // Sort to make consistent ?
356          std::sort(whichColumns,whichColumns+numberSort);
357          saveModel = new ClpSimplex(this,numberSort,whichColumns);
358          delete [] whichColumns;
359          delete [] weight;
360          // Skip factorization
361          //statusOfProblemInPrimal(lastCleaned,factorType,progress_,false,saveModel);
362          //statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
363          stopSprint = numberIterations_+numberSprintIterations;
364          printf("Sprint with %d columns for %d iterations\n",
365                 numberSprintColumns,numberSprintIterations);
366        }
367      }
368
369      // Say good factorization
370      factorType=1;
371
372      // Say no pivot has occurred (for steepest edge and updates)
373      pivotRow_=-2;
374
375      // exit if victory declared
376      if (problemStatus_>=0)
377        break;
378
379      // Iterate
380      whileIterating(ifValuesPass);
381    }
382  }
383  // if infeasible get real values
384  if (problemStatus_==1) {
385    infeasibilityCost_=0.0;
386    createRim(7);
387    nonLinearCost_->checkInfeasibilities(0.0);
388    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
389    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
390    // and get good feasible duals
391    computeDuals(NULL);
392  }
393  // clean up
394  unflag();
395  finish();
396  restoreData(data);
397  return problemStatus_;
398}
399/*
400  Reasons to come out:
401  -1 iterations etc
402  -2 inaccuracy
403  -3 slight inaccuracy (and done iterations)
404  -4 end of values pass and done iterations
405  +0 looks optimal (might be infeasible - but we will investigate)
406  +2 looks unbounded
407  +3 max iterations
408*/
409int
410ClpSimplexPrimal::whileIterating(int valuesOption)
411{
412  // Say if values pass
413  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
414  int returnCode=-1;
415  int superBasicType=1;
416  if (valuesOption>1)
417    superBasicType=3;
418  // status stays at -1 while iterating, >=0 finished, -2 to invert
419  // status -3 to go to top without an invert
420  while (problemStatus_==-1) {
421    //#define CLP_DEBUG 1
422#ifdef CLP_DEBUG
423    {
424      int i;
425      // not [1] as has information
426      for (i=0;i<4;i++) {
427        if (i!=1)
428          rowArray_[i]->checkClear();
429      }
430      for (i=0;i<2;i++) {
431        columnArray_[i]->checkClear();
432      }
433    }
434#endif
435#if CLP_DEBUG>2
436    // very expensive
437    if (numberIterations_>0&&numberIterations_<-2534) {
438      handler_->setLogLevel(63);
439      double saveValue = objectiveValue_;
440      double * saveRow1 = new double[numberRows_];
441      double * saveRow2 = new double[numberRows_];
442      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
443      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
444      double * saveColumn1 = new double[numberColumns_];
445      double * saveColumn2 = new double[numberColumns_];
446      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
447      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
448      createRim(7);
449      gutsOfSolution(NULL,NULL);
450      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
451             numberIterations_,
452             saveValue,objectiveValue_,sumPrimalInfeasibilities_);
453      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
454      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
455      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
456      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
457      delete [] saveRow1;
458      delete [] saveRow2;
459      delete [] saveColumn1;
460      delete [] saveColumn2;
461      objectiveValue_=saveValue;
462    }
463#endif
464    if (!ifValuesPass) {
465      // choose column to come in
466      // can use pivotRow_ to update weights
467      // pass in list of cost changes so can do row updates (rowArray_[1])
468      // NOTE rowArray_[0] is used by computeDuals which is a
469      // slow way of getting duals but might be used
470      primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
471                   columnArray_[0],columnArray_[1]);
472    } else {
473      // in values pass
474      int sequenceIn=nextSuperBasic(superBasicType,columnArray_[0]);
475      if (valuesOption>1)
476        superBasicType=2;
477      if (sequenceIn<0) {
478        // end of values pass - initialize weights etc
479        handler_->message(CLP_END_VALUES_PASS,messages_)
480          <<numberIterations_;
481        primalColumnPivot_->saveWeights(this,5);
482        problemStatus_=-2; // factorize now
483        pivotRow_=-1; // say no weights update
484        returnCode=-4;
485        // Clean up
486        int i;
487        for (i=0;i<numberRows_+numberColumns_;i++) {
488          if (getColumnStatus(i)==atLowerBound||getColumnStatus(i)==isFixed)
489            solution_[i]=lower_[i];
490          else if (getColumnStatus(i)==atUpperBound)
491            solution_[i]=upper_[i];
492        }
493        break;
494      } else {
495        // normal
496        sequenceIn_ = sequenceIn;
497        valueIn_=solution_[sequenceIn_];
498        lowerIn_=lower_[sequenceIn_];
499        upperIn_=upper_[sequenceIn_];
500        dualIn_=dj_[sequenceIn_];
501      }
502    }
503    pivotRow_=-1;
504    sequenceOut_=-1;
505    rowArray_[1]->clear();
506    if (sequenceIn_>=0) {
507      // we found a pivot column
508      assert (!flagged(sequenceIn_));
509#ifdef CLP_DEBUG
510      if ((handler_->logLevel()&32)) {
511        char x = isColumn(sequenceIn_) ? 'C' :'R';
512        std::cout<<"pivot column "<<
513          x<<sequenceWithin(sequenceIn_)<<std::endl;
514      }
515#endif
516      // do second half of iteration
517      returnCode = pivotResult(ifValuesPass);
518      if (returnCode<-1&&returnCode>-5) {
519        problemStatus_=-2; //
520      } else if (returnCode==-5) {
521        // something flagged - continue;
522      } else if (returnCode==2) {
523        problemStatus_=-5; // looks unbounded
524      } else if (returnCode==4) {
525        problemStatus_=-2; // looks unbounded but has iterated
526      } else if (returnCode!=-1) {
527        assert(returnCode==3);
528        problemStatus_=3;
529      }
530    } else {
531      // no pivot column
532#ifdef CLP_DEBUG
533      if (handler_->logLevel()&32)
534        printf("** no column pivot\n");
535#endif
536      if (nonLinearCost_->numberInfeasibilities())
537        problemStatus_=-4; // might be infeasible
538      returnCode=0;
539      break;
540    }
541  }
542  if (valuesOption>1)
543    columnArray_[0]->setNumElements(0);
544  return returnCode;
545}
546/* Checks if finished.  Updates status */
547void
548ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type,
549                                          ClpSimplexProgress * progress,
550                                          bool doFactorization,
551                                          ClpSimplex * originalModel)
552{
553  if (type==2) {
554    // trouble - restore solution
555    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
556    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
557           numberRows_*sizeof(double));
558    memcpy(columnActivityWork_,savedSolution_ ,
559           numberColumns_*sizeof(double));
560    forceFactorization_=1; // a bit drastic but ..
561    pivotRow_=-1; // say no weights update
563  }
564  int saveThreshold = factorization_->sparseThreshold();
565  int tentativeStatus = problemStatus_;
566  if (problemStatus_>-3||problemStatus_==-4) {
567    // factorize
568    // later on we will need to recover from singularities
569    // also we could skip if first time
570    // do weights
571    // This may save pivotRow_ for use
572    if (doFactorization)
573    primalColumnPivot_->saveWeights(this,1);
574
575    if (type&&doFactorization) {
576      // is factorization okay?
577      if (internalFactorize(1)) {
578        if (solveType_==2) {
579          // say odd
580          problemStatus_=5;
581          return;
582        }
583#if 1
584        // switch off dense
585        int saveDense = factorization_->denseThreshold();
586        factorization_->setDenseThreshold(0);
587        internalFactorize(2);
588        factorization_->setDenseThreshold(saveDense);
589#else
590
591        // no - restore previous basis
592        assert (type==1);
593        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
594        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
595               numberRows_*sizeof(double));
596        memcpy(columnActivityWork_,savedSolution_ ,
597               numberColumns_*sizeof(double));
598        forceFactorization_=1; // a bit drastic but ..
599        type = 2;
600        assert (internalFactorize(1)==0);
601#endif
603      }
604    }
605    if (problemStatus_!=-4)
606      problemStatus_=-3;
607  }
608  // at this stage status is -3 or -5 if looks unbounded
609  // get primal and dual solutions
610  // put back original costs and then check
611  createRim(4);
612  gutsOfSolution(NULL,NULL,(firstFree_>=0));
613  // Double check reduced costs if no action
614  if (progress->lastIterationNumber(0)==numberIterations_) {
615    if (primalColumnPivot_->looksOptimal()) {
616      numberDualInfeasibilities_ = 0;
617      sumDualInfeasibilities_ = 0.0;
618    }
619  }
620  // Check if looping
621  int loop;
622  if (type!=2)
623    loop = progress->looping();
624  else
625    loop=-1;
626  if (loop>=0) {
627    if (!problemStatus_) {
628      // declaring victory
629      numberPrimalInfeasibilities_ = 0;
630      sumPrimalInfeasibilities_ = 0.0;
631    } else {
632      problemStatus_ = loop; //exit if in loop
633      problemStatus_ = 10; // instead - try other algorithm
634    }
635    problemStatus_ = 10; // instead - try other algorithm
636    return ;
637  } else if (loop<-1) {
638    // Is it time for drastic measures
640        progress->oddState()<10&&progress->oddState()>=0) {
641      progress->newOddState();
642      nonLinearCost_->zapCosts();
643    }
644    // something may have changed
645    gutsOfSolution(NULL,NULL);
646  }
647  // If progress then reset costs
648  if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
649    createRim(4,false); // costs back
650    delete nonLinearCost_;
651    nonLinearCost_ = new ClpNonLinearCost(this);
652    progress->endOddState();
653    gutsOfSolution(NULL,NULL);
654  }
655  // Flag to say whether to go to dual to clean up
656  bool goToDual=false;
658  //if((progressFlag_&2)!=0)
659  //problemStatus_=-1;;
660  progressFlag_ = 0; //reset progress flag
661
662  handler_->message(CLP_SIMPLEX_STATUS,messages_)
663    <<numberIterations_<<nonLinearCost_->feasibleReportCost();
664  handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
665    <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
666  handler_->printing(sumDualInfeasibilities_>0.0)
667    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
668  handler_->printing(numberDualInfeasibilitiesWithoutFree_
669                     <numberDualInfeasibilities_)
670                       <<numberDualInfeasibilitiesWithoutFree_;
671  handler_->message()<<CoinMessageEol;
672  if (!primalFeasible()) {
673    nonLinearCost_->checkInfeasibilities(primalTolerance_);
674    gutsOfSolution(NULL,NULL);
675  }
676  // we may wish to say it is optimal even if infeasible
677  bool alwaysOptimal = (specialOptions_&1)!=0;
678  // give code benefit of doubt
679  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
680      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
681    // say optimal (with these bounds etc)
682    numberDualInfeasibilities_ = 0;
683    sumDualInfeasibilities_ = 0.0;
684    numberPrimalInfeasibilities_ = 0;
685    sumPrimalInfeasibilities_ = 0.0;
686    // But check if in sprint
687    if (originalModel) {
688      // Carry on and re-do
689      numberDualInfeasibilities_ = -776;
690    }
691  }
692  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
693  if (dualFeasible()||problemStatus_==-4) {
694    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
695      //may need infeasiblity cost changed
696      // we can see if we can construct a ray
697      // make up a new objective
698      double saveWeight = infeasibilityCost_;
699      // save nonlinear cost as we are going to switch off costs
700      ClpNonLinearCost * nonLinear = nonLinearCost_;
701      // do twice to make sure Primal solution has settled
702      // put non-basics to bounds in case tolerance moved
703      // put back original costs
704      createRim(4);
705      nonLinearCost_->checkInfeasibilities(0.0);
706      gutsOfSolution(NULL,NULL);
707
708      infeasibilityCost_=1.0e100;
709      // put back original costs
710      createRim(4);
711      nonLinearCost_->checkInfeasibilities(primalTolerance_);
712      // may have fixed infeasibilities - double check
713      if (nonLinearCost_->numberInfeasibilities()==0) {
714        // carry on
715        problemStatus_ = -1;
716        infeasibilityCost_=saveWeight;
717        nonLinearCost_->checkInfeasibilities(primalTolerance_);
718      } else {
719        nonLinearCost_=NULL;
720        // scale
721        int i;
722        for (i=0;i<numberRows_+numberColumns_;i++)
723          cost_[i] *= 1.0e-95;
724        gutsOfSolution(NULL,NULL);
725        nonLinearCost_=nonLinear;
726        infeasibilityCost_=saveWeight;
727        if ((infeasibilityCost_>=1.0e18||
728             numberDualInfeasibilities_==0)&&perturbation_==101) {
729          goToDual=unPerturb(); // stop any further perturbation
730          if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
731            goToDual=false;
732          nonLinearCost_->checkInfeasibilities(primalTolerance_);
733          numberDualInfeasibilities_=1; // carry on
734          problemStatus_=-1;
735        }
736        if (infeasibilityCost_>=1.0e20||
737            numberDualInfeasibilities_==0) {
738          // we are infeasible - use as ray
739          delete [] ray_;
740          ray_ = new double [numberRows_];
741          memcpy(ray_,dual_,numberRows_*sizeof(double));
742          // and get feasible duals
743          infeasibilityCost_=0.0;
744          createRim(4);
745          nonLinearCost_->checkInfeasibilities(primalTolerance_);
746          gutsOfSolution(NULL,NULL);
747          // so will exit
748          infeasibilityCost_=1.0e30;
749          // reset infeasibilities
750          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
751          numberPrimalInfeasibilities_=
752            nonLinearCost_->numberInfeasibilities();
753        }
754        if (infeasibilityCost_<1.0e20) {
755          infeasibilityCost_ *= 5.0;
757          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
758            <<infeasibilityCost_
759            <<CoinMessageEol;
760          // put back original costs and then check
761          createRim(4);
762          nonLinearCost_->checkInfeasibilities(0.0);
763          gutsOfSolution(NULL,NULL);
764          problemStatus_=-1; //continue
765          goToDual=false;
766        } else {
767          // say infeasible
768          problemStatus_ = 1;
769        }
770      }
771    } else {
772      // may be optimal
773      if (perturbation_==101) {
774        goToDual=unPerturb(); // stop any further perturbation
775        lastCleaned=-1; // carry on
776      }
777      bool unflagged = unflag();
778      if ( lastCleaned!=numberIterations_||unflagged) {
779        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
780          <<primalTolerance_
781          <<CoinMessageEol;
782        if (numberTimesOptimal_<4) {
783          numberTimesOptimal_++;
785          if (numberTimesOptimal_==1) {
786            // better to have small tolerance even if slower
787            factorization_->zeroTolerance(1.0e-15);
788          }
789          lastCleaned=numberIterations_;
790          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
791            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
792              <<CoinMessageEol;
793          double oldTolerance = primalTolerance_;
794          primalTolerance_=dblParam_[ClpPrimalTolerance];
795#if 0
796          double * xcost = new double[numberRows_+numberColumns_];
797          double * xlower = new double[numberRows_+numberColumns_];
798          double * xupper = new double[numberRows_+numberColumns_];
799          double * xdj = new double[numberRows_+numberColumns_];
800          double * xsolution = new double[numberRows_+numberColumns_];
801          memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
802          memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
803          memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
804          memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
805          memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
806#endif
807          // put back original costs and then check
808          createRim(4);
809          nonLinearCost_->checkInfeasibilities(oldTolerance);
810#if 0
811          int i;
812          for (i=0;i<numberRows_+numberColumns_;i++) {
813            if (cost_[i]!=xcost[i])
814              printf("** %d old cost %g new %g sol %g\n",
815                     i,xcost[i],cost_[i],solution_[i]);
816            if (lower_[i]!=xlower[i])
817              printf("** %d old lower %g new %g sol %g\n",
818                     i,xlower[i],lower_[i],solution_[i]);
819            if (upper_[i]!=xupper[i])
820              printf("** %d old upper %g new %g sol %g\n",
821                     i,xupper[i],upper_[i],solution_[i]);
822            if (dj_[i]!=xdj[i])
823              printf("** %d old dj %g new %g sol %g\n",
824                     i,xdj[i],dj_[i],solution_[i]);
825            if (solution_[i]!=xsolution[i])
826              printf("** %d old solution %g new %g sol %g\n",
827                     i,xsolution[i],solution_[i],solution_[i]);
828          }
829          delete [] xcost;
830          delete [] xupper;
831          delete [] xlower;
832          delete [] xdj;
833          delete [] xsolution;
834#endif
835          gutsOfSolution(NULL,NULL);
836          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
837              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
838            // say optimal (with these bounds etc)
839            numberDualInfeasibilities_ = 0;
840            sumDualInfeasibilities_ = 0.0;
841            numberPrimalInfeasibilities_ = 0;
842            sumPrimalInfeasibilities_ = 0.0;
843          }
844          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
845            problemStatus_=0;
846          else
847            problemStatus_ = -1;
848        } else {
849          problemStatus_=0; // optimal
850          if (lastCleaned<numberIterations_) {
851            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
852              <<CoinMessageEol;
853          }
854        }
855      } else {
856        problemStatus_=0; // optimal
857      }
858    }
859  } else {
860    // see if looks unbounded
861    if (problemStatus_==-5) {
862      if (nonLinearCost_->numberInfeasibilities()) {
863        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
864          // back off weight
865          infeasibilityCost_ = 1.0e13;
866          unPerturb(); // stop any further perturbation
867        }
868        //we need infeasiblity cost changed
869        if (infeasibilityCost_<1.0e20) {
870          infeasibilityCost_ *= 5.0;
872          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
873            <<infeasibilityCost_
874            <<CoinMessageEol;
875          // put back original costs and then check
876          createRim(4);
877          gutsOfSolution(NULL,NULL);
878          problemStatus_=-1; //continue
879        } else {
880          // say unbounded
881          problemStatus_ = 2;
882        }
883      } else {
884        // say unbounded
885        problemStatus_ = 2;
886      }
887    } else {
888      if(type==3&&problemStatus_!=-5)
889        unflag(); // odd
890      // carry on
891      problemStatus_ = -1;
892    }
893  }
894  if (type==0||type==1) {
895    if (type!=1||!saveStatus_) {
896      // create save arrays
897      delete [] saveStatus_;
898      delete [] savedSolution_;
899      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
900      savedSolution_ = new double [numberRows_+numberColumns_];
901    }
902    // save arrays
903    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
904    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
905           numberRows_*sizeof(double));
906    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
907  }
908  if (doFactorization) {
909    // restore weights (if saved) - also recompute infeasibility list
910    if (tentativeStatus>-3)
911      primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
912    else
913      primalColumnPivot_->saveWeights(this,3);
914    if (saveThreshold) {
915      // use default at present
916      factorization_->sparseThreshold(0);
917      factorization_->goSparse();
918    }
919  }
921    problemStatus_=4; // unknown
922  }
923  lastGoodIteration_ = numberIterations_;
924  if (goToDual)
925    problemStatus_=10; // try dual
926#if 0
927  double thisObj = progress->lastObjective(0);
928  double lastObj = progress->lastObjective(1);
929  if (lastObj<thisObj-1.0e-7*max(fabs(thisObj),fabs(lastObj))-1.0e-8
930      &&firstFree_<0) {
931    int maxFactor = factorization_->maximumPivots();
932    if (maxFactor>10) {
933      if (forceFactorization_<0)
934        forceFactorization_= maxFactor;
935      forceFactorization_ = max (1,(forceFactorization_>>1));
936      printf("Reducing factorization frequency\n");
937    }
938  }
939#endif
940}
941/*
942   Row array has pivot column
943   This chooses pivot row.
944   For speed, we may need to go to a bucket approach when many
945   variables go through bounds
946   On exit rhsArray will have changes in costs of basic variables
947*/
948void
949ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
950                            CoinIndexedVector * rhsArray,
951                            CoinIndexedVector * spareArray,
952                            CoinIndexedVector * spareArray2,
953                            int valuesPass)
954{
955  //rowArray->scanAndPack();
956  if (valuesPass) {
957    dualIn_ = cost_[sequenceIn_];
958
959    double * work=rowArray->denseVector();
960    int number=rowArray->getNumElements();
961    int * which=rowArray->getIndices();
962
963    int iIndex;
964    for (iIndex=0;iIndex<number;iIndex++) {
965
966      int iRow = which[iIndex];
967      double alpha = work[iIndex];
968      int iPivot=pivotVariable_[iRow];
969      dualIn_ -= alpha*cost_[iPivot];
970    }
971    // determine direction here
972    if (dualIn_<-dualTolerance_) {
973      directionIn_=1;
974    } else if (dualIn_>dualTolerance_) {
975      directionIn_=-1;
976    } else {
977      // towards nearest bound
978      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
979        directionIn_=-1;
980        dualIn_=dualTolerance_;
981      } else {
982        directionIn_=1;
983        dualIn_=-dualTolerance_;
984      }
985    }
986  }
987
988  // sequence stays as row number until end
989  pivotRow_=-1;
990  int numberRemaining=0;
991
992  double totalThru=0.0; // for when variables flip
993  double acceptablePivot=1.0e-7;
994  if (factorization_->pivots())
995    acceptablePivot=1.0e-5; // if we have iterated be more strict
996  double bestEverPivot=acceptablePivot;
997  int lastPivotRow = -1;
998  double lastPivot=0.0;
999  double lastTheta=1.0e50;
1000
1001  // use spareArrays to put ones looked at in
1002  // First one is list of candidates
1003  // We could compress if we really know we won't need any more
1004  // Second array has current set of pivot candidates
1005  // with a backup list saved in double * part of indexed vector
1006
1007  // pivot elements
1008  double * spare;
1009  // indices
1010  int * index;
1011  spareArray->clear();
1012  spare = spareArray->denseVector();
1013  index = spareArray->getIndices();
1014
1015  // we also need somewhere for effective rhs
1016  double * rhs=rhsArray->denseVector();
1017  //int * indexRhs = rhsArray->getIndices();
1018  //int numberFlip=0; // Those which may change if flips
1019
1020  /*
1021    First we get a list of possible pivots.  We can also see if the
1022    problem looks unbounded.
1023
1024    At first we increase theta and see what happens.  We start
1025    theta at a reasonable guess.  If in right area then we do bit by bit.
1026    We save possible pivot candidates
1027
1028   */
1029
1030  // do first pass to get possibles
1031  // We can also see if unbounded
1032
1033  double * work=rowArray->denseVector();
1034  int number=rowArray->getNumElements();
1035  int * which=rowArray->getIndices();
1036
1037  // we need to swap sign if coming in from ub
1038  double way = directionIn_;
1039  double maximumMovement;
1040  if (way>0.0)
1041    maximumMovement = min(1.0e30,upperIn_-valueIn_);
1042  else
1043    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
1044
1045  double averageTheta = nonLinearCost_->averageTheta();
1046  double tentativeTheta = min(10.0*averageTheta,maximumMovement);
1047  double upperTheta = maximumMovement;
1048  if (tentativeTheta>0.5*maximumMovement)
1049    tentativeTheta=maximumMovement;
1050
1051  double dualCheck = fabs(dualIn_);
1052  // but make a bit more pessimistic
1053  dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1054
1055  int iIndex;
1056  bool gotList=false;
1057  int pivotOne=-1;
1058  while (!gotList) {
1059    pivotOne=-1;
1060    totalThru=0.0;
1061    // We also re-compute reduced cost
1062    numberRemaining=0;
1063    dualIn_ = cost_[sequenceIn_];
1064    double tolerance = primalTolerance_*1.002;
1065    for (iIndex=0;iIndex<number;iIndex++) {
1066
1067      int iRow = which[iIndex];
1068      double alpha = work[iIndex];
1069      int iPivot=pivotVariable_[iRow];
1070      dualIn_ -= alpha*cost_[iPivot];
1071      alpha *= way;
1072      double oldValue = solution_[iPivot];
1073      // get where in bound sequence
1074      if (alpha>0.0) {
1075        // basic variable going towards lower bound
1076        double bound = lower_[iPivot];
1077        oldValue -= bound;
1078      } else {
1079        // basic variable going towards upper bound
1080        double bound = upper_[iPivot];
1081        oldValue = bound-oldValue;
1082      }
1083
1084      double value = oldValue-tentativeTheta*fabs(alpha);
1085      assert (oldValue>=-tolerance);
1086      if (value<=tolerance) {
1087        value=oldValue-upperTheta*fabs(alpha);
1088        if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot) {
1089          upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1090          pivotOne=numberRemaining;
1091        }
1093        spare[numberRemaining]=alpha;
1094        rhs[numberRemaining]=oldValue;
1095        index[numberRemaining++]=iRow;
1096        totalThru += fabs(alpha);
1097        setActive(iRow);
1098        //} else if (value<primalTolerance_*1.002) {
1099        // May change if is a flip
1100        //indexRhs[numberFlip++]=iRow;
1101      }
1102    }
1103    if (upperTheta<maximumMovement&&totalThru*infeasibilityCost_>=1.0001*dualCheck) {
1104      // Can pivot here
1105      gotList=true;
1106    } else if (tentativeTheta<maximumMovement) {
1107      //printf("Going round with average theta of %g\n",averageTheta);
1108      tentativeTheta=maximumMovement;
1109    } else {
1110      gotList=true;
1111    }
1112  }
1113  totalThru=0.0;
1114  // we need to keep where rhs non-zeros are
1115  int numberInRhs=numberRemaining;
1116  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
1117  rhsArray->setNumElements(numberInRhs);
1118  rhsArray->setPacked();
1119
1120  theta_=maximumMovement;
1121
1122  bool goBackOne = false;
1123
1124  //printf("%d remain out of %d\n",numberRemaining,number);
1125  int iTry=0;
1126#define MAXTRY 1000
1127  if (numberRemaining&&upperTheta<maximumMovement) {
1128    // First check if previously chosen one will work
1129    if (pivotOne>=0&&0) {
1130      double thruCost = infeasibilityCost_*fabs(spare[pivotOne]);
1131      if (thruCost>=0.99*fabs(dualIn_))
1132        printf("Could pivot on %d as change %g dj %g\n",
1133               index[pivotOne],thruCost,dualIn_);
1134      double alpha = spare[pivotOne];
1135      double oldValue = rhs[pivotOne];
1136      theta_ = oldValue/fabs(alpha);
1137      pivotRow_=pivotOne;
1138      // Stop loop
1139      iTry=MAXTRY;
1140    }
1141
1142    // first get ratio with tolerance
1143    for ( ;iTry<MAXTRY;iTry++) {
1144
1145      upperTheta=maximumMovement;
1146      int iBest=-1;
1147      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1148
1149        double alpha = fabs(spare[iIndex]);
1150        double oldValue = rhs[iIndex];
1151        double value = oldValue-upperTheta*alpha;
1152
1153        if (value<-primalTolerance_ && alpha>=acceptablePivot) {
1154          upperTheta = (oldValue+primalTolerance_)/alpha;
1155          iBest=iIndex; // just in case weird numbers
1156        }
1157      }
1158
1159      // now look at best in this lot
1160      double bestPivot=acceptablePivot;
1161      pivotRow_=-1;
1162      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1163
1164        int iRow = index[iIndex];
1165        double alpha = spare[iIndex];
1166        double oldValue = rhs[iIndex];
1167        double value = oldValue-upperTheta*fabs(alpha);
1168
1169        if (value<=0||iBest==iIndex) {
1170          // how much would it cost to go thru and modify bound
1171          totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha,rhs[iIndex]);
1172          setActive(iRow);
1173          if (fabs(alpha)>bestPivot) {
1174            bestPivot=fabs(alpha);
1175            theta_ = oldValue/bestPivot;
1176            pivotRow_=iIndex;
1177          }
1178        }
1179      }
1180      if (bestPivot<0.1*bestEverPivot&&
1181          bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
1182        // back to previous one
1183        goBackOne = true;
1184        break;
1185      } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1186        if (lastPivot>acceptablePivot) {
1187          // back to previous one
1188          goBackOne = true;
1189        } else {
1190          // can only get here if all pivots so far too small
1191        }
1192        break;
1193      } else if (totalThru>=dualCheck) {
1194        break; // no point trying another loop
1195      } else {
1196        lastPivotRow=pivotRow_;
1197        lastTheta = theta_;
1198        if (bestPivot>bestEverPivot)
1199          bestEverPivot=bestPivot;
1200      }
1201    }
1202    // can get here without pivotRow_ set but with lastPivotRow
1203    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1204      // back to previous one
1205      pivotRow_=lastPivotRow;
1206      theta_ = lastTheta;
1207    }
1208  } else {
1209    // mark ones which may move
1210    //for (int i=0;i<numberFlip;i++) {
1211    //int iRow= indexRhs[i];
1212    //setActive(iRow);
1213    //}
1214  }
1215  //if (iTry>50)
1216  //printf("** %d tries\n",iTry);
1217  if (pivotRow_>=0) {
1218    int position=pivotRow_; // position in list
1219    pivotRow_=index[position];
1220    //alpha_ = work[pivotRow_];
1221    alpha_ = way*spare[position];
1222    // translate to sequence
1223    sequenceOut_ = pivotVariable_[pivotRow_];
1224    valueOut_ = solution(sequenceOut_);
1225    lowerOut_=lower_[sequenceOut_];
1226    upperOut_=upper_[sequenceOut_];
1227#define MINIMUMTHETA 1.0e-12
1228    // Movement should be minimum for anti-degeneracy - unless
1229    // fixed variable out
1230    double minimumTheta;
1231    if (upperOut_>lowerOut_)
1232      minimumTheta=MINIMUMTHETA;
1233    else
1234      minimumTheta=0.0;
1235    // will we need to increase tolerance
1236    //#define CLP_DEBUG
1237#ifdef CLP_DEBUG
1238    bool found=false;
1239#endif
1240    double largestInfeasibility = primalTolerance_;
1241    if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
1242      theta_=minimumTheta;
1243      for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) {
1244        largestInfeasibility = max (largestInfeasibility,
1245                                    -(rhs[iIndex]-fabs(spare[iIndex])*theta_));
1246      }
1247#define CLP_DEBUG
1248#ifdef CLP_DEBUG
1249      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
1250        printf("Primal tolerance increased from %g to %g\n",
1251               primalTolerance_,largestInfeasibility);
1252#endif
1253#undef CLP_DEBUG
1254      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1255    }
1256    // Need to look at all in some cases
1257    if (theta_>tentativeTheta) {
1258      for (iIndex=0;iIndex<number;iIndex++)
1259        setActive(which[iIndex]);
1260    }
1261    if (way<0.0)
1262      theta_ = - theta_;
1263    double newValue = valueOut_ - theta_*alpha_;
1264    // If  4 bit set - Force outgoing variables to exact bound (primal)
1265    if (alpha_*way<0.0) {
1266      directionOut_=-1;      // to upper bound
1267      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1268        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1269      } else {
1270          upperOut_ = newValue;
1271      }
1272    } else {
1273      directionOut_=1;      // to lower bound
1274      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1275        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1276      } else {
1277        lowerOut_ = newValue;
1278      }
1279    }
1280    dualOut_ = reducedCost(sequenceOut_);
1281  } else if (maximumMovement<1.0e20) {
1282    // flip
1283    pivotRow_ = -2; // so we can tell its a flip
1284    sequenceOut_ = sequenceIn_;
1285    valueOut_ = valueIn_;
1286    dualOut_ = dualIn_;
1287    lowerOut_ = lowerIn_;
1288    upperOut_ = upperIn_;
1289    alpha_ = 0.0;
1290    if (way<0.0) {
1291      directionOut_=1;      // to lower bound
1292      theta_ = lowerOut_ - valueOut_;
1293    } else {
1294      directionOut_=-1;      // to upper bound
1295      theta_ = upperOut_ - valueOut_;
1296    }
1297  }
1298
1299  double theta1 = max(theta_,1.0e-12);
1300  double theta2 = numberIterations_*nonLinearCost_->averageTheta();
1301  // Set average theta
1302  nonLinearCost_->setAverageTheta((theta1+theta2)/((double) (numberIterations_+1)));
1303  //if (numberIterations_%1000==0)
1304  //printf("average theta is %g\n",nonLinearCost_->averageTheta());
1305
1306  // clear arrays
1307
1308  memset(spare,0,numberRemaining*sizeof(double));
1309
1310  // put back original bounds etc
1311  nonLinearCost_->goBackAll(rhsArray);
1312
1313  rhsArray->clear();
1314
1315}
1316/*
1317   Chooses primal pivot column
1318   updateArray has cost updates (also use pivotRow_ from last iteration)
1319   Would be faster with separate region to scan
1320   and will have this (with square of infeasibility) when steepest
1321   For easy problems we can just choose one of the first columns we look at
1322*/
1323void
1325                               CoinIndexedVector * spareRow1,
1326                               CoinIndexedVector * spareRow2,
1327                               CoinIndexedVector * spareColumn1,
1328                               CoinIndexedVector * spareColumn2)
1329{
1331                                               spareRow2,spareColumn1,
1332                                               spareColumn2);
1333  if (sequenceIn_>=0) {
1334    valueIn_=solution_[sequenceIn_];
1335    dualIn_=dj_[sequenceIn_];
1336    if (nonLinearCost_->lookBothWays()) {
1337      // double check
1338      ClpSimplex::Status status = getStatus(sequenceIn_);
1339
1340      switch(status) {
1341      case ClpSimplex::atUpperBound:
1342        if (dualIn_<0.0) {
1343          // move to other side
1344          printf("For %d U (%g, %g, %g) dj changed from %g",
1345                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1346                 upper_[sequenceIn_],dualIn_);
1347          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
1348          printf(" to %g\n",dualIn_);
1349          nonLinearCost_->setOne(sequenceIn_,upper_[sequenceIn_]+2.0*currentPrimalTolerance());
1350          setStatus(sequenceIn_,ClpSimplex::atLowerBound);
1351        }
1352        break;
1353      case ClpSimplex::atLowerBound:
1354        if (dualIn_>0.0) {
1355          // move to other side
1356          printf("For %d L (%g, %g, %g) dj changed from %g",
1357                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1358                 upper_[sequenceIn_],dualIn_);
1359          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
1360          printf(" to %g\n",dualIn_);
1361          nonLinearCost_->setOne(sequenceIn_,lower_[sequenceIn_]-2.0*currentPrimalTolerance());
1362          setStatus(sequenceIn_,ClpSimplex::atUpperBound);
1363        }
1364        break;
1365      default:
1366        break;
1367      }
1368    }
1369    lowerIn_=lower_[sequenceIn_];
1370    upperIn_=upper_[sequenceIn_];
1371    if (dualIn_>0.0)
1372      directionIn_ = -1;
1373    else
1374      directionIn_ = 1;
1375  } else {
1376    sequenceIn_ = -1;
1377  }
1378}
1379/* The primals are updated by the given array.
1380   Returns number of infeasibilities.
1381   After rowArray will have list of cost changes
1382*/
1383int
1384ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
1385                                        double theta,
1386                                        double & objectiveChange,
1387                                        int valuesPass)
1388{
1389  // Cost on pivot row may change - may need to change dualIn
1390  double oldCost=0.0;
1391  if (pivotRow_>=0)
1392    oldCost = cost_[sequenceOut_];
1393  //rowArray->scanAndPack();
1394  double * work=rowArray->denseVector();
1395  int number=rowArray->getNumElements();
1396  int * which=rowArray->getIndices();
1397
1398  int newNumber = 0;
1399  int pivotPosition = -1;
1400  nonLinearCost_->setChangeInCost(0.0);
1401  int iIndex;
1402  if (!valuesPass) {
1403    for (iIndex=0;iIndex<number;iIndex++) {
1404
1405      int iRow = which[iIndex];
1406      double alpha = work[iIndex];
1407      work[iIndex]=0.0;
1408      int iPivot=pivotVariable_[iRow];
1409      double change = theta*alpha;
1410      double value = solution_[iPivot] - change;
1411      solution_[iPivot]=value;
1412#ifndef NDEBUG
1413      // check if not active then okay
1414      if (!active(iRow)) {
1415        // But make sure one going out is feasible
1416        if (change>0.0) {
1417          // going down
1418          if (value<lower_[iPivot]+primalTolerance_) {
1419            if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
1420              value=lower_[iPivot];
1421            double difference = nonLinearCost_->setOne(iPivot,value);
1422            assert (!difference);
1423          }
1424        } else {
1425          // going up
1426          if (value>upper_[iPivot]-primalTolerance_) {
1427            if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
1428              value=upper_[iPivot];
1429            double difference = nonLinearCost_->setOne(iPivot,value);
1430            assert (!difference);
1431          }
1432        }
1433      }
1434#endif
1435      if (active(iRow)) {
1436        clearActive(iRow);
1437        // But make sure one going out is feasible
1438        if (change>0.0) {
1439          // going down
1440          if (value<lower_[iPivot]+primalTolerance_) {
1441            if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
1442              value=lower_[iPivot];
1443            double difference = nonLinearCost_->setOne(iPivot,value);
1444            if (difference) {
1445              if (iRow==pivotRow_)
1446                pivotPosition=newNumber;
1447              work[newNumber] = difference;
1448              //change reduced cost on this
1449              dj_[iPivot] = -difference;
1450              which[newNumber++]=iRow;
1451            }
1452          }
1453        } else {
1454          // going up
1455          if (value>upper_[iPivot]-primalTolerance_) {
1456            if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
1457              value=upper_[iPivot];
1458            double difference = nonLinearCost_->setOne(iPivot,value);
1459            if (difference) {
1460              if (iRow==pivotRow_)
1461                pivotPosition=newNumber;
1462              work[newNumber] = difference;
1463              //change reduced cost on this
1464              dj_[iPivot] = -difference;
1465              which[newNumber++]=iRow;
1466            }
1467          }
1468        }
1469      }
1470    }
1471  } else {
1472    // values pass so look at all
1473    for (iIndex=0;iIndex<number;iIndex++) {
1474
1475      int iRow = which[iIndex];
1476      double alpha = work[iIndex];
1477      work[iIndex]=0.0;
1478      int iPivot=pivotVariable_[iRow];
1479      double change = theta*alpha;
1480      double value = solution_[iPivot] - change;
1481      solution_[iPivot]=value;
1482      clearActive(iRow);
1483      // But make sure one going out is feasible
1484      if (change>0.0) {
1485        // going down
1486        if (value<lower_[iPivot]+primalTolerance_) {
1487          if (iPivot==sequenceOut_&&value>lower_[iPivot]-1.001*primalTolerance_)
1488            value=lower_[iPivot];
1489          double difference = nonLinearCost_->setOne(iPivot,value);
1490          if (difference) {
1491            if (iRow==pivotRow_)
1492              pivotPosition=newNumber;
1493            work[newNumber] = difference;
1494            //change reduced cost on this
1495            dj_[iPivot] = -difference;
1496            which[newNumber++]=iRow;
1497          }
1498        }
1499      } else {
1500        // going up
1501        if (value>upper_[iPivot]-primalTolerance_) {
1502          if (iPivot==sequenceOut_&&value<upper_[iPivot]+1.001*primalTolerance_)
1503            value=upper_[iPivot];
1504          double difference = nonLinearCost_->setOne(iPivot,value);
1505          if (difference) {
1506            if (iRow==pivotRow_)
1507              pivotPosition=newNumber;
1508            work[newNumber] = difference;
1509            //change reduced cost on this
1510            dj_[iPivot] = -difference;
1511            which[newNumber++]=iRow;
1512          }
1513        }
1514      }
1515    }
1516  }
1517  objectiveChange += nonLinearCost_->changeInCost();
1518  rowArray->setPacked();
1519#if 0
1520  rowArray->setNumElements(newNumber);
1521  rowArray->expand();
1522  if (pivotRow_>=0) {
1523    dualIn_ += (oldCost-cost_[sequenceOut_]);
1524    // update change vector to include pivot
1526    // and convert to packed
1527    rowArray->scanAndPack();
1528  } else {
1529    // and convert to packed
1530    rowArray->scanAndPack();
1531  }
1532#else
1533  if (pivotRow_>=0) {
1534    dualIn_ += (oldCost-cost_[sequenceOut_]);
1535    // update change vector to include pivot
1536    if (pivotPosition>=0) {
1537      work[pivotPosition] -= dualIn_;
1538    } else {
1539      work[newNumber]=-dualIn_;
1540      which[newNumber++]=pivotRow_;
1541    }
1542  }
1543  rowArray->setNumElements(newNumber);
1544#endif
1545  return 0;
1546}
1547// Perturbs problem
1548void
1549ClpSimplexPrimal::perturb(int type)
1550{
1551  if (perturbation_>100)
1553  if (perturbation_==100)
1554    perturbation_=50; // treat as normal
1555  int savePerturbation = perturbation_;
1556  int i;
1557  if (!numberIterations_)
1558    cleanStatus(); // make sure status okay
1559  if (!numberIterations_&&perturbation_==50) {
1560    // See if we need to perturb
1561    double * sort = new double[numberRows_];
1562    for (i=0;i<numberRows_;i++) {
1563      double lo = fabs(lower_[i]);
1564      double up = fabs(upper_[i]);
1565      double value=0.0;
1566      if (lo&&lo<1.0e20) {
1567        if (up&&up<1.0e20)
1568          value = 0.5*(lo+up);
1569        else
1570          value=lo;
1571      } else {
1572        if (up&&up<1.0e20)
1573          value = up;
1574      }
1575      sort[i]=value;
1576    }
1577    std::sort(sort,sort+numberRows_);
1578    int number=1;
1579    double last = sort[0];
1580    for (i=1;i<numberRows_;i++) {
1581      if (last!=sort[i])
1582        number++;
1583      last=sort[i];
1584    }
1585    delete [] sort;
1586    if (number*4>numberRows_) {
1587      perturbation_=100;
1588      return; // good enough
1589    }
1590  }
1591  // primal perturbation
1592  double perturbation=1.0e-20;
1593  int numberNonZero=0;
1594  // maximum fraction of rhs/bounds to perturb
1595  double maximumFraction = 1.0e-5;
1596  if (perturbation_>=50) {
1597    perturbation = 1.0e-4;
1598    for (i=0;i<numberColumns_+numberRows_;i++) {
1599      if (upper_[i]>lower_[i]+primalTolerance_) {
1600        double lowerValue, upperValue;
1601        if (lower_[i]>-1.0e20)
1602          lowerValue = fabs(lower_[i]);
1603        else
1604          lowerValue=0.0;
1605        if (upper_[i]<1.0e20)
1606          upperValue = fabs(upper_[i]);
1607        else
1608          upperValue=0.0;
1609        double value = max(fabs(lowerValue),fabs(upperValue));
1610        value = min(value,upper_[i]-lower_[i]);
1611#if 1
1612        if (value) {
1613          perturbation += value;
1614          numberNonZero++;
1615        }
1616#else
1617        perturbation = max(perturbation,value);
1618#endif
1619      }
1620    }
1621    if (numberNonZero)
1622      perturbation /= (double) numberNonZero;
1623    else
1624      perturbation = 1.0e-1;
1625  } else if (perturbation_<100) {
1626    perturbation = pow(10.0,perturbation_);
1627    // user is in charge
1628    maximumFraction = 1.0;
1629  }
1630  double largestZero=0.0;
1631  double largest=0.0;
1632  double largestPerCent=0.0;
1633  bool printOut=(handler_->logLevel()==63);
1634  printOut=false; //off
1635  // Check if all slack
1636  int number=0;
1637  int iSequence;
1638  for (iSequence=0;iSequence<numberRows_;iSequence++) {
1639    if (getRowStatus(iSequence)==basic)
1640      number++;
1641  }
1642  if (number!=numberRows_)
1643    type=1;
1644  // modify bounds
1645  // Change so at least 1.0e-5 and no more than 0.1
1646  // For now just no more than 0.1
1647  // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
1648  if (type==1) {
1649    //double multiplier = perturbation*maximumFraction;
1650    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1651      if (getStatus(iSequence)==basic) {
1652        double solutionValue = solution_[iSequence];
1653        double lowerValue = lower_[iSequence];
1654        double upperValue = upper_[iSequence];
1655        double difference = upperValue-lowerValue;
1656        difference = min(difference,perturbation);
1657        difference = min(difference,fabs(solutionValue)+1.0);
1658        double value = maximumFraction*(difference+1.0);
1659        value = min(value,0.1);
1660        value *= CoinDrand48();
1661        if (solutionValue-lowerValue<=primalTolerance_) {
1662          lower_[iSequence] -= value;
1663        } else if (upperValue-solutionValue<=primalTolerance_) {
1664          upper_[iSequence] += value;
1665        } else {
1666#if 0
1667          if (iSequence>=numberColumns_) {
1668            // may not be at bound - but still perturb (unless free)
1669            if (upperValue>1.0e30&&lowerValue<-1.0e30)
1670              value=0.0;
1671            else
1672              value = - value; // as -1.0 in matrix
1673          } else {
1674            value = 0.0;
1675          }
1676#else
1677          value=0.0;
1678#endif
1679        }
1680        if (value) {
1681          if (printOut)
1682            printf("col %d lower from %g to %g, upper from %g to %g\n",
1683                   iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
1684          if (solutionValue) {
1685            largest = max(largest,value);
1686            if (value>(fabs(solutionValue)+1.0)*largestPerCent)
1687              largestPerCent=value/(fabs(solutionValue)+1.0);
1688          } else {
1689            largestZero = max(largestZero,value);
1690          }
1691        }
1692      }
1693    }
1694  } else {
1695    double tolerance = 100.0*primalTolerance_;
1696    for (i=0;i<numberColumns_;i++) {
1697      double lowerValue=lower_[i], upperValue=upper_[i];
1698      if (upperValue>lowerValue+primalTolerance_) {
1699        double value = perturbation*maximumFraction;
1700        value = min(value,0.1);
1701        value *= CoinDrand48();
1702        if (savePerturbation!=50) {
1703          if (fabs(value)<=primalTolerance_)
1704            value=0.0;
1705          if (lowerValue>-1.0e20&&lowerValue)
1706            lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
1707          if (upperValue<1.0e20&&upperValue)
1708            upperValue += value * (max(1.0e-2,1.0e-5*fabs(upperValue)));
1709        } else if (value) {
1710          double valueL =value *(max(1.0e-2,1.0e-5*fabs(lowerValue)));
1711          // get in range
1712          if (valueL<=tolerance) {
1713            valueL *= 10.0;
1714            while (valueL<=tolerance)
1715              valueL *= 10.0;
1716          } else if (valueL>1.0) {
1717            valueL *= 0.1;
1718            while (valueL>1.0)
1719              valueL *= 0.1;
1720          }
1721          if (lowerValue>-1.0e20&&lowerValue)
1722            lowerValue -= valueL;
1723          double valueU =value *(max(1.0e-2,1.0e-5*fabs(upperValue)));
1724          // get in range
1725          if (valueU<=tolerance) {
1726            valueU *= 10.0;
1727            while (valueU<=tolerance)
1728              valueU *= 10.0;
1729          } else if (valueU>1.0) {
1730            valueU *= 0.1;
1731            while (valueU>1.0)
1732              valueU *= 0.1;
1733          }
1734          if (upperValue<1.0e20&&upperValue)
1735            upperValue += valueU;
1736        }
1737        if (lowerValue!=lower_[i]) {
1738          double difference = fabs(lowerValue-lower_[i]);
1739          largest = max(largest,difference);
1740          if (difference>fabs(lower_[i])*largestPerCent)
1741            largestPerCent=fabs(difference/lower_[i]);
1742        }
1743        if (upperValue!=upper_[i]) {
1744          double difference = fabs(upperValue-upper_[i]);
1745          largest = max(largest,difference);
1746          if (difference>fabs(upper_[i])*largestPerCent)
1747            largestPerCent=fabs(difference/upper_[i]);
1748        }
1749        if (printOut)
1750          printf("col %d lower from %g to %g, upper from %g to %g\n",
1751                 i,lower_[i],lowerValue,upper_[i],upperValue);
1752      }
1753      lower_[i]=lowerValue;
1754      upper_[i]=upperValue;
1755    }
1756    for (;i<numberColumns_+numberRows_;i++) {
1757      double lowerValue=lower_[i], upperValue=upper_[i];
1758      double value = perturbation*maximumFraction;
1759      value = min(value,0.1);
1760      value *= CoinDrand48();
1761      if (upperValue>lowerValue+tolerance) {
1762        if (savePerturbation!=50) {
1763          if (fabs(value)<=primalTolerance_)
1764            value=0.0;
1765          if (lowerValue>-1.0e20&&lowerValue)
1766            lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
1767          if (upperValue<1.0e20&&upperValue)
1768            upperValue += value * (max(1.0e-2,1.0e-5*fabs(upperValue)));
1769        } else if (value) {
1770          double valueL =value *(max(1.0e-2,1.0e-5*fabs(lowerValue)));
1771          // get in range
1772          if (valueL<=tolerance) {
1773            valueL *= 10.0;
1774            while (valueL<=tolerance)
1775              valueL *= 10.0;
1776          } else if (valueL>1.0) {
1777            valueL *= 0.1;
1778            while (valueL>1.0)
1779              valueL *= 0.1;
1780          }
1781          if (lowerValue>-1.0e20&&lowerValue)
1782            lowerValue -= valueL;
1783          double valueU =value *(max(1.0e-2,1.0e-5*fabs(upperValue)));
1784          // get in range
1785          if (valueU<=tolerance) {
1786            valueU *= 10.0;
1787            while (valueU<=tolerance)
1788              valueU *= 10.0;
1789          } else if (valueU>1.0) {
1790            valueU *= 0.1;
1791            while (valueU>1.0)
1792              valueU *= 0.1;
1793          }
1794          if (upperValue<1.0e20&&upperValue)
1795            upperValue += valueU;
1796        }
1797      } else if (upperValue>0.0) {
1798        upperValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
1799        lowerValue -= value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
1800      } else if (upperValue<0.0) {
1801        upperValue += value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
1802        lowerValue += value * (max(1.0e-2,1.0e-5*fabs(lowerValue)));
1803      } else {
1804      }
1805      if (lowerValue!=lower_[i]) {
1806        double difference = fabs(lowerValue-lower_[i]);
1807        largest = max(largest,difference);
1808        if (difference>fabs(lower_[i])*largestPerCent)
1809          largestPerCent=fabs(difference/lower_[i]);
1810      }
1811      if (upperValue!=upper_[i]) {
1812        double difference = fabs(upperValue-upper_[i]);
1813        largest = max(largest,difference);
1814        if (difference>fabs(upper_[i])*largestPerCent)
1815          largestPerCent=fabs(difference/upper_[i]);
1816      }
1817      if (printOut)
1818        printf("row %d lower from %g to %g, upper from %g to %g\n",
1819               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
1820      lower_[i]=lowerValue;
1821      upper_[i]=upperValue;
1822    }
1823  }
1824  // Clean up
1825  for (i=0;i<numberColumns_+numberRows_;i++) {
1826    switch(getStatus(i)) {
1827
1828    case basic:
1829      break;
1830    case atUpperBound:
1831      solution_[i]=upper_[i];
1832      break;
1833    case isFixed:
1834    case atLowerBound:
1835      solution_[i]=lower_[i];
1836      break;
1837    case isFree:
1838      break;
1839    case superBasic:
1840      break;
1841    }
1842  }
1843  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
1844    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
1845    <<CoinMessageEol;
1846  // redo nonlinear costs
1847  // say perturbed
1848  perturbation_=101;
1849}
1850// un perturb
1851bool
1852ClpSimplexPrimal::unPerturb()
1853{
1854  if (perturbation_!=101)
1855    return false;
1856  // put back original bounds and costs
1857  createRim(7);
1858  sanityCheck();
1859  // unflag
1860  unflag();
1861  // get a valid nonlinear cost function
1862  delete nonLinearCost_;
1863  nonLinearCost_= new ClpNonLinearCost(this);
1864  perturbation_ = 102; // stop any further perturbation
1865  // move non basic variables to new bounds
1866  nonLinearCost_->checkInfeasibilities(0.0);
1867#if 1
1868  // Try using dual
1869  return true;
1870#else
1871  gutsOfSolution(NULL,NULL);
1872  return false;
1873#endif
1874
1875}
1876// Unflag all variables and return number unflagged
1877int
1878ClpSimplexPrimal::unflag()
1879{
1880  int i;
1881  int number = numberRows_+numberColumns_;
1882  int numberFlagged=0;
1883  for (i=0;i<number;i++) {
1884    if (flagged(i)) {
1885      clearFlagged(i);
1886      numberFlagged++;
1887    }
1888  }
1889  return numberFlagged;
1890}
1891// Do not change infeasibility cost and always say optimal
1892void
1893ClpSimplexPrimal::alwaysOptimal(bool onOff)
1894{
1895  if (onOff)
1896    specialOptions_ |= 1;
1897  else
1898    specialOptions_ &= ~1;
1899}
1900bool
1901ClpSimplexPrimal::alwaysOptimal() const
1902{
1903  return (specialOptions_&1)!=0;
1904}
1905// Flatten outgoing variables i.e. - always to exact bound
1906void
1907ClpSimplexPrimal::exactOutgoing(bool onOff)
1908{
1909  if (onOff)
1910    specialOptions_ |= 4;
1911  else
1912    specialOptions_ &= ~4;
1913}
1914bool
1915ClpSimplexPrimal::exactOutgoing() const
1916{
1917  return (specialOptions_&4)!=0;
1918}
1919/*
1920  Reasons to come out (normal mode/user mode):
1921  -1 normal
1922  -2 factorize now - good iteration/ NA
1923  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
1924  -4 inaccuracy - refactorize - no iteration/ NA
1925  -5 something flagged - go round again/ pivot not possible
1926  +2 looks unbounded
1927  +3 max iterations (iteration done)
1928*/
1929int
1930ClpSimplexPrimal::pivotResult(int ifValuesPass)
1931{
1932
1933  bool roundAgain=true;
1934  int returnCode=-1;
1935
1936  // loop round if user setting and doing refactorization
1937  while (roundAgain) {
1938    roundAgain=false;
1939    returnCode=-1;
1940    pivotRow_=-1;
1941    sequenceOut_=-1;
1942    rowArray_[1]->clear();
1943#if 0
1944    {
1945      int seq[]={612,643};
1946      int k;
1947      for (k=0;k<sizeof(seq)/sizeof(int);k++) {
1948        int iSeq=seq[k];
1949        if (getColumnStatus(iSeq)!=basic) {
1950          double djval;
1951          double * work;
1952          int number;
1953          int * which;
1954
1955          int iIndex;
1956          unpack(rowArray_[1],iSeq);
1957          factorization_->updateColumn(rowArray_[2],rowArray_[1]);
1958          djval = cost_[iSeq];
1959          work=rowArray_[1]->denseVector();
1960          number=rowArray_[1]->getNumElements();
1961          which=rowArray_[1]->getIndices();
1962
1963          for (iIndex=0;iIndex<number;iIndex++) {
1964
1965            int iRow = which[iIndex];
1966            double alpha = work[iRow];
1967            int iPivot=pivotVariable_[iRow];
1968            djval -= alpha*cost_[iPivot];
1969          }
1970          double comp = 1.0e-8 + 1.0e-7*(max(fabs(dj_[iSeq]),fabs(djval)));
1971          if (fabs(djval-dj_[iSeq])>comp)
1972            printf("Bad dj %g for %d - true is %g\n",
1973                   dj_[iSeq],iSeq,djval);
1974          assert (fabs(djval)<1.0e-3||djval*dj_[iSeq]>0.0);
1975          rowArray_[1]->clear();
1976        }
1977      }
1978    }
1979#endif
1980
1981    // we found a pivot column
1982    // update the incoming column
1983    unpackPacked(rowArray_[1]);
1984    // save reduced cost
1985    double saveDj = dualIn_;
1986    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1987    // do ratio test and re-compute dj
1988    primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
1989              ifValuesPass);
1990    if (ifValuesPass) {
1991      saveDj=dualIn_;
1992      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
1993        if(fabs(dualIn_)<1.0e2*dualTolerance_) {
1994          // try other way
1995          directionIn_=-directionIn_;
1996          primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
1997                    0);
1998        }
1999        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2000          if (solveType_==1) {
2001            // reject it
2002            char x = isColumn(sequenceIn_) ? 'C' :'R';
2003            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2004              <<x<<sequenceWithin(sequenceIn_)
2005              <<CoinMessageEol;
2006            setFlagged(sequenceIn_);
2008            lastBadIteration_ = numberIterations_; // say be more cautious
2009            clearAll();
2010            pivotRow_=-1;
2011          }
2012          returnCode=-5;
2013          break;
2014        }
2015      }
2016    }
2017    double checkValue=1.0e-2;
2018    if (largestDualError_>1.0e-5)
2019      checkValue=1.0e-1;
2020    if (solveType_==1&&((saveDj*dualIn_<1.0e-20&&!ifValuesPass)||
2021        fabs(saveDj-dualIn_)>checkValue*(1.0+fabs(saveDj)))) {
2022      handler_->message(CLP_PRIMAL_DJ,messages_)
2023        <<saveDj<<dualIn_
2024        <<CoinMessageEol;
2025      if(lastGoodIteration_ != numberIterations_) {
2026        clearAll();
2027        pivotRow_=-1; // say no weights update
2028        returnCode=-4;
2029        if(lastGoodIteration_+1 == numberIterations_) {
2030          // not looking wonderful - try cleaning bounds
2031          // put non-basics to bounds in case tolerance moved
2032          nonLinearCost_->checkInfeasibilities(0.0);
2033        }
2034        sequenceOut_=-1;
2035        break;
2036      } else {
2037        // take on more relaxed criterion
2038        if (saveDj*dualIn_<1.0e-20||
2039            fabs(saveDj-dualIn_)>2.0e-1*(1.0+fabs(dualIn_))) {
2040          // need to reject something
2041          char x = isColumn(sequenceIn_) ? 'C' :'R';
2042          handler_->message(CLP_SIMPLEX_FLAG,messages_)
2043            <<x<<sequenceWithin(sequenceIn_)
2044            <<CoinMessageEol;
2045          setFlagged(sequenceIn_);
2047          lastBadIteration_ = numberIterations_; // say be more cautious
2048          clearAll();
2049          pivotRow_=-1;
2050          returnCode=-5;
2051          sequenceOut_=-1;
2052          break;
2053        }
2054      }
2055    }
2056    if (pivotRow_>=0) {
2057      if (solveType_==2) {
2058        // **** Coding for user interface
2059        // do ray
2060        primalRay(rowArray_[1]);
2061        // update duals
2062        if (pivotRow_>=0) {
2063          // as packed need to find pivot row
2064          //assert (rowArray_[1]->packedMode());
2065          //int i;
2066
2067          //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
2068          assert (fabs(alpha_)>1.0e-8);
2069          double multiplier = dualIn_/alpha_;
2070          rowArray_[0]->insert(pivotRow_,multiplier);
2071          factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
2072          // put row of tableau in rowArray[0] and columnArray[0]
2073          matrix_->transposeTimes(this,-1.0,
2074                                  rowArray_[0],columnArray_[1],columnArray_[0]);
2075          // update column djs
2076          int i;
2077          int * index = columnArray_[0]->getIndices();
2078          int number = columnArray_[0]->getNumElements();
2079          double * element = columnArray_[0]->denseVector();
2080          for (i=0;i<number;i++) {
2081            int ii = index[i];
2082            dj_[ii] += element[ii];
2083            element[ii]=0.0;
2084          }
2085          columnArray_[0]->setNumElements(0);
2086          // and row djs
2087          index = rowArray_[0]->getIndices();
2088          number = rowArray_[0]->getNumElements();
2089          element = rowArray_[0]->denseVector();
2090          for (i=0;i<number;i++) {
2091            int ii = index[i];
2092            dj_[ii+numberColumns_] += element[ii];
2093            dual_[ii] = dj_[ii+numberColumns_];
2094            element[ii]=0.0;
2095          }
2096          rowArray_[0]->setNumElements(0);
2097          // check incoming
2098          assert (fabs(dj_[sequenceIn_])<1.0e-6);
2099        }
2100      }
2101      // if stable replace in basis
2103                                                       pivotRow_,
2104                                                       alpha_);
2105      // if no pivots, bad update but reasonable alpha - take and invert
2107          lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2110        // slight error
2112          returnCode=-3;
2113        }
2114      } else if (updateStatus==2) {
2115        // major error
2116        // better to have small tolerance even if slower
2117        factorization_->zeroTolerance(1.0e-15);
2118        int maxFactor = factorization_->maximumPivots();
2119        if (maxFactor>10) {
2120          if (forceFactorization_<0)
2121            forceFactorization_= maxFactor;
2122          forceFactorization_ = max (1,(forceFactorization_>>1));
2123        }
2124        // later we may need to unwind more e.g. fake bounds
2125        if(lastGoodIteration_ != numberIterations_) {
2126          clearAll();
2127          pivotRow_=-1;
2128          if (solveType_==1) {
2129            returnCode=-4;
2130            break;
2131          } else {
2132            // user in charge - re-factorize
2133            int lastCleaned;
2134            ClpSimplexProgress dummyProgress;
2135            if (saveStatus_)
2136              statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true);
2137            else
2138              statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true);
2139            roundAgain=true;
2140            continue;
2141          }
2142        } else {
2143          // need to reject something
2144          if (solveType_==1) {
2145            char x = isColumn(sequenceIn_) ? 'C' :'R';
2146            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2147              <<x<<sequenceWithin(sequenceIn_)
2148              <<CoinMessageEol;
2149            setFlagged(sequenceIn_);
2151          }
2152          lastBadIteration_ = numberIterations_; // say be more cautious
2153          clearAll();
2154          pivotRow_=-1;
2155          sequenceOut_=-1;
2156          returnCode = -5;
2157          break;
2158
2159        }
2160      } else if (updateStatus==3) {
2161        // out of memory
2162        // increase space if not many iterations
2163        if (factorization_->pivots()<
2164            0.5*factorization_->maximumPivots()&&
2165            factorization_->pivots()<200)
2166          factorization_->areaFactor(
2167                                     factorization_->areaFactor() * 1.1);
2168        returnCode =-2; // factorize now
2169      } else if (updateStatus==5) {
2170        problemStatus_=-2; // factorize now
2171      }
2172      // here do part of steepest - ready for next iteration
2173      if (!ifValuesPass)
2174        primalColumnPivot_->updateWeights(rowArray_[1]);
2175    } else {
2176      if (pivotRow_==-1) {
2177        // no outgoing row is valid
2178        rowArray_[0]->clear();
2179        if (!factorization_->pivots()) {
2180          returnCode = 2; //say looks unbounded
2181          // do ray
2182          primalRay(rowArray_[1]);
2183        } else if (solveType_==2) {
2184          // refactorize
2185          int lastCleaned;
2186          ClpSimplexProgress dummyProgress;
2187          if (saveStatus_)
2188            statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true);
2189          else
2190            statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true);
2191          roundAgain=true;
2192          continue;
2193        } else {
2194          returnCode = 4; //say looks unbounded but has iterated
2195        }
2196        break;
2197      } else {
2198        // flipping from bound to bound
2199      }
2200    }
2201
2202
2203    // update primal solution
2204
2205    double objectiveChange=0.0;
2206    // after this rowArray_[1] is not empty - used to update djs
2207    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2208
2209    double oldValue = valueIn_;
2210    if (directionIn_==-1) {
2211      // as if from upper bound
2212      if (sequenceIn_!=sequenceOut_) {
2213        // variable becoming basic
2214        valueIn_ -= fabs(theta_);
2215      } else {
2216        valueIn_=lowerIn_;
2217      }
2218    } else {
2219      // as if from lower bound
2220      if (sequenceIn_!=sequenceOut_) {
2221        // variable becoming basic
2222        valueIn_ += fabs(theta_);
2223      } else {
2224        valueIn_=upperIn_;
2225      }
2226    }
2227    objectiveChange += dualIn_*(valueIn_-oldValue);
2228    // outgoing
2229    if (sequenceIn_!=sequenceOut_) {
2230      if (directionOut_>0) {
2231        valueOut_ = lowerOut_;
2232      } else {
2233        valueOut_ = upperOut_;
2234      }
2235      if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2236        valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2237      else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2238        valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2239      // may not be exactly at bound and bounds may have changed
2240      // Make sure outgoing looks feasible
2241      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2242      solution_[sequenceOut_]=valueOut_;
2243    }
2244    // change cost and bounds on incoming if primal
2245    nonLinearCost_->setOne(sequenceIn_,valueIn_);
2246    int whatNext=housekeeping(objectiveChange);
2247
2248#if 0
2249    if (numberIterations_==1148)
2250      whatNext=1;
2251    if (numberIterations_>1200)
2252    exit(0);
2253#endif
2254    if (whatNext==1) {
2255        returnCode =-2; // refactorize
2256    } else if (whatNext==2) {
2257      // maximum iterations or equivalent
2258      returnCode=3;
2259    } else if(numberIterations_ == lastGoodIteration_
2260              + 2 * factorization_->maximumPivots()) {
2261      // done a lot of flips - be safe
2262      returnCode =-2; // refactorize
2263    }
2264  }
2265  if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
2266    // refactorize here
2267    int lastCleaned;
2268    ClpSimplexProgress dummyProgress;
2269    if (saveStatus_)
2270      statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true);
2271    else
2272      statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true);
2273    if (problemStatus_==5) {
2274      printf("Singular basis\n");
2275      problemStatus_=-1;
2276      returnCode=5;
2277    }
2278  }
2279#ifdef CLP_DEBUG
2280  {
2281    int i;
2282    // not [1] as may have information
2283    for (i=0;i<4;i++) {
2284      if (i!=1)
2285        rowArray_[i]->checkClear();
2286    }
2287    for (i=0;i<2;i++) {
2288      columnArray_[i]->checkClear();
2289    }
2290  }
2291#endif
2292  return returnCode;
2293}
2294// Create primal ray
2295void
2296ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
2297{
2298  delete [] ray_;
2299  ray_ = new double [numberColumns_];
2300  ClpFillN(ray_,numberColumns_,0.0);
2301  int number=rowArray->getNumElements();
2302  int * index = rowArray->getIndices();
2303  double * array = rowArray->denseVector();
2304  double way=-directionIn_;
2305  int i;
2306  double zeroTolerance=1.0e-12;
2307  if (sequenceIn_<numberColumns_)
2308    ray_[sequenceIn_]=directionIn_;
2309  if (!rowArray->packedMode()) {
2310    for (i=0;i<number;i++) {
2311      int iRow=index[i];
2312      int iPivot=pivotVariable_[iRow];
2313      double arrayValue = array[iRow];
2314      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2315        ray_[iPivot] = way* arrayValue;
2316    }
2317  } else {
2318    for (i=0;i<number;i++) {
2319      int iRow=index[i];
2320      int iPivot=pivotVariable_[iRow];
2321      double arrayValue = array[i];
2322      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2323        ray_[iPivot] = way* arrayValue;
2324    }
2325  }
2326}
2327/* Get next superbasic -1 if none,
2328   Normal type is 1
2329   If type is 3 then initializes sorted list
2330   if 2 uses list.
2331*/
2332int
2333ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
2334{
2335  if (firstFree_>=0&&superBasicType) {
2336    int returnValue=firstFree_;
2337    int iColumn=firstFree_+1;
2338    if (superBasicType>1) {
2339      if (superBasicType>2) {
2340        // Initialize list
2341        // Wild guess that lower bound more natural than upper
2342        int number=0;
2343        double * work=columnArray->denseVector();
2344        int * which=columnArray->getIndices();
2345        for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2346          if (getStatus(iColumn)==superBasic) {
2347            if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2348              solution_[iColumn]=lower_[iColumn];
2349              setStatus(iColumn,atLowerBound);
2350            } else if (fabs(solution_[iColumn]-upper_[iColumn])
2351                       <=primalTolerance_) {
2352              solution_[iColumn]=upper_[iColumn];
2353              setStatus(iColumn,atUpperBound);
2354            } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2355              setStatus(iColumn,isFree);
2356              break;
2357            } else if (!flagged(iColumn)) {
2358              // put ones near bounds at end after sorting
2359              work[number]= - min(0.1*(solution_[iColumn]-lower_[iColumn]),
2360                                  upper_[iColumn]-solution_[iColumn]);
2361              which[number++] = iColumn;
2362            }
2363          }
2364        }
2365        CoinSort_2(work,work+number,which);
2366        columnArray->setNumElements(number);
2367        memset(work,0,number*sizeof(double));
2368      }
2369      int * which=columnArray->getIndices();
2370      int number = columnArray->getNumElements();
2371      if (!number) {
2372        // finished
2373        iColumn = numberRows_+numberColumns_;
2374        returnValue=-1;
2375      } else {
2376        number--;
2377        returnValue=which[number];
2378        iColumn=returnValue;
2379        columnArray->setNumElements(number);
2380      }
2381    } else {
2382      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
2383        if (getStatus(iColumn)==superBasic) {
2384          if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2385            solution_[iColumn]=lower_[iColumn];
2386            setStatus(iColumn,atLowerBound);
2387          } else if (fabs(solution_[iColumn]-upper_[iColumn])
2388                     <=primalTolerance_) {
2389            solution_[iColumn]=upper_[iColumn];
2390            setStatus(iColumn,atUpperBound);
2391          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2392            setStatus(iColumn,isFree);
2393            break;
2394          } else {
2395            break;
2396          }
2397        }
2398      }
2399    }
2400    firstFree_ = iColumn;
2401    if (firstFree_==numberRows_+numberColumns_)
2402      firstFree_=-1;
2403    return returnValue;
2404  } else {
2405    return -1;
2406  }
2407}
2408void
2409ClpSimplexPrimal::clearAll()
2410{
2411  int number=rowArray_[1]->getNumElements();
2412  int * which=rowArray_[1]->getIndices();
2413
2414  int iIndex;
2415  for (iIndex=0;iIndex<number;iIndex++) {
2416
2417    int iRow = which[iIndex];
2418    clearActive(iRow);
2419  }
2420  rowArray_[1]->clear();
2421}
2422
Note: See TracBrowser for help on using the repository browser.