# source:trunk/ClpSimplexPrimal.cpp@344

Last change on this file since 344 was 344, checked in by forrest, 17 years ago

event handling

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