source: branches/devel-1/ClpSimplexPrimal.cpp @ 2351

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

Some changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 43.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5/* Notes on implementation of primal simplex algorithm.
6
7   When primal feasible(A):
8
9   If dual feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to enter basis from a bound (B).  We now need to find an
11   outgoing variable which will leave problem primal feasible so we get
12   the column of the tableau corresponding to the incoming variable
13   (with the correct sign depending if variable will go up or down).
14
15   We now perform a ratio test to determine which outgoing variable will
16   preserve primal feasibility (C).  If no variable found then problem
17   is unbounded (in primal sense).  If there is a variable, we then
18   perform pivot and repeat.  Trivial?
19
20   -------------------------------------------
21
22   A) How do we get primal feasible?  All variables have fake costs
23   outside their feasible region so it is trivial to declare problem
24   feasible.  OSL did not have a phase 1/phase 2 approach but
25   instead effectively put an extra cost on infeasible basic variables
26   I am taking the same approach here, although it is generalized
27   to allow for non-linear costs and dual information.
28
29   In OSL, this weight was changed heuristically, here at present
30   it is only increased if problem looks finished.  If problem is
31   feasible I check for unboundedness.  If not unbounded we
32   could play with going into dual.  As long as weights increase
33   any algorithm would be finite.
34   
35   B) Which incoming variable to choose is a virtual base class.
36   For difficult problems steepest edge is preferred while for
37   very easy (large) problems we will need partial scan.
38
39   C) Sounds easy, but this is hardest part of algorithm.
40      1) Instead of stopping at first choice, we may be able
41      to allow that variable to go through bound and if objective
42      still improving choose again.  These mini iterations can
43      increase speed by orders of magnitude but we may need to
44      go to more of a bucket choice of variable rather than looking
45      at them one by one (for speed).
46      2) Accuracy.  Basic infeasibilities may be less than
47      tolerance.  Pivoting on these makes objective go backwards.
48      OSL modified cost so a zero move was made, Gill et al
49      modified so a strictly positive move was made.
50      The two problems are that re-factorizations can
51      change rinfeasibilities above and below tolerances and that when
52      finished we need to reset costs and try again.
53      3) Degeneracy.  Gill et al helps but may not be enough.  We
54      may need more.  Also it can improve speed a lot if we perturb
55      the costs significantly. 
56
57  References:
58     Forrest and Goldfarb, Steepest-edge simplex algorithms for
59       linear programming - Mathematical Programming 1992
60     Forrest and Tomlin, Implementing the simplex method for
61       the Optimization Subroutine Library - IBM Systems Journal 1992
62     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
63       Procedure for Linear and Nonlinear Programming SOL report 1988
64
65
66  TODO:
67 
68  a) Better recovery procedures.  At present I never check on forward
69     progress.  There is checkpoint/restart with reducing
70     re-factorization frequency, but this is only on singular
71     factorizations.
72  b) Fast methods for large easy problems (and also the option for
73     the code to automatically choose which method).
74  c) We need to be able to stop in various ways for OSI - this
75     is fairly easy.
76
77 */
78
79#if defined(_MSC_VER)
80// Turn off compiler warning about long names
81#  pragma warning(disable:4786)
82#endif
83
84#include <math.h>
85
86#include "CoinHelperFunctions.hpp"
87#include "ClpSimplexPrimal.hpp"
88#include "ClpFactorization.hpp"
89#include "ClpNonLinearCost.hpp"
90#include "CoinPackedMatrix.hpp"
91#include "CoinIndexedVector.hpp"
92#include "CoinWarmStartBasis.hpp"
93#include "ClpPrimalColumnPivot.hpp"
94#include "ClpMessage.hpp"
95#include <cfloat>
96#include <cassert>
97#include <string>
98#include <stdio.h>
99#include <iostream>
100// primal
101int ClpSimplexPrimal::primal (int ifValuesPass )
102{
103
104  /*
105      Method
106
107     It tries to be a single phase approach with a weight of 1.0 being
108     given to getting optimal and a weight of infeasibilityCost_ being
109     given to getting primal feasible.  In this version I have tried to
110     be clever in a stupid way.  The idea of fake bounds in dual
111     seems to work so the primal analogue would be that of getting
112     bounds on reduced costs (by a presolve approach) and using
113     these for being above or below feasible region.  I decided to waste
114     memory and keep these explicitly.  This allows for non-linear
115     costs!
116
117     The code is designed to take advantage of sparsity so arrays are
118     seldom zeroed out from scratch or gone over in their entirety.
119     The only exception is a full scan to find incoming variable for
120     Dantzig row choice.  For steepest edge we keep an updated list
121     of dual infeasibilities (actually squares). 
122     On easy problems we don't need full scan - just
123     pick first reasonable.
124
125     One problem is how to tackle degeneracy and accuracy.  At present
126     I am using the modification of costs which I put in OSL and which was
127     extended by Gill et al.  I am still not sure of the exact details.
128
129     The flow of primal is three while loops as follows:
130
131     while (not finished) {
132
133       while (not clean solution) {
134
135          Factorize and/or clean up solution by changing bounds so
136          primal feasible.  If looks finished check fake primal bounds.
137          Repeat until status is iterating (-1) or finished (0,1,2)
138
139       }
140
141       while (status==-1) {
142
143         Iterate until no pivot in or out or time to re-factorize.
144
145         Flow is:
146
147         choose pivot column (incoming variable).  if none then
148         we are primal feasible so looks as if done but we need to
149         break and check bounds etc.
150
151         Get pivot column in tableau
152
153         Choose outgoing row.  If we don't find one then we look
154         primal unbounded so break and check bounds etc.  (Also the
155         pivot tolerance is larger after any iterations so that may be
156         reason)
157
158         If we do find outgoing row, we may have to adjust costs to
159         keep going forwards (anti-degeneracy).  Check pivot will be stable
160         and if unstable throw away iteration and break to re-factorize.
161         If minor error re-factorize after iteration.
162
163         Update everything (this may involve changing bounds on
164         variables to stay primal feasible.
165
166       }
167
168     }
169
170     At present we never check we are going forwards.  I overdid that in
171     OSL so will try and make a last resort.
172
173     Needs partial scan pivot in option.
174
175     May need other anti-degeneracy measures, especially if we try and use
176     loose tolerances as a way to solve in fewer iterations.
177
178     I like idea of dynamic scaling.  This gives opportunity to decouple
179     different implications of scaling for accuracy, iteration count and
180     feasibility tolerance.
181
182  */
183
184  // sanity check
185  assert (numberRows_==matrix_->getNumRows());
186  assert (numberColumns_==matrix_->getNumCols());
187  // for moment all arrays must exist
188  assert(columnLower_);
189  assert(columnUpper_);
190  assert(rowLower_);
191  assert(rowUpper_);
192
193
194  algorithm_ = +1;
195  primalTolerance_=dblParam_[ClpPrimalTolerance];
196  dualTolerance_=dblParam_[ClpDualTolerance];
197
198  // put in standard form (and make row copy)
199  // create modifiable copies of model rim and do optional scaling
200  bool goodMatrix=createRim(7+8+16,true);
201
202  // save infeasibility cost
203  double saveInfeasibilityCost = infeasibilityCost_;
204  // Save perturbation
205  int savePerturbation = perturbation_;
206  // save if sparse factorization wanted
207  int saveSparse = factorization_->sparseThreshold();
208
209  if (goodMatrix) {
210    // Model looks okay
211    int iRow,iColumn;
212    // Do initial factorization
213    // and set certain stuff
214    // We can either set increasing rows so ...IsBasic gives pivot row
215    // or we can just increment iBasic one by one
216    // for now let ...iBasic give pivot row
217    factorization_->increasingRows(2);
218    // row activities have negative sign
219    factorization_->slackValue(-1.0);
220    factorization_->zeroTolerance(1.0e-13);
221   
222   
223    // If user asked for perturbation - do it
224   
225    if (perturbation_<100) 
226      perturb();
227   
228    // for primal we will change bounds using infeasibilityCost_
229    if (nonLinearCost_==NULL) {
230      // get a valid nonlinear cost function
231      delete nonLinearCost_;
232      nonLinearCost_= new ClpNonLinearCost(this);
233    }
234   
235    // loop round to clean up solution if values pass
236    int numberThrownOut = -1;
237    int firstSuperBasic=numberRows_+numberColumns_;
238    int totalNumberThrownOut=0;
239    while(numberThrownOut) {
240      int status = internalFactorize(0+10*ifValuesPass);
241      if (status<0)
242        return 1; // some error
243      else
244        totalNumberThrownOut+= status;
245     
246      // for this we need clean basis so it is after factorize
247      numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
248                                     ifValuesPass);
249      totalNumberThrownOut+= numberThrownOut;
250     
251      // find first superbasic - columns, then rows
252      if (ifValuesPass) {
253        nextSuperBasic(firstSuperBasic);
254        if (firstSuperBasic==numberRows_+numberColumns_)
255          ifValuesPass=0; // signal no values pass
256      }
257    }
258   
259    if (totalNumberThrownOut)
260      handler_->message(CLP_SINGULARITIES,messages_)
261        <<totalNumberThrownOut
262        <<CoinMessageEol;
263   
264    problemStatus_ = -1;
265    numberIterations_=0;
266   
267    int lastCleaned=0; // last time objective or bounds cleaned up
268   
269    // number of times we have declared optimality
270    numberTimesOptimal_=0;
271   
272    // Progress indicator
273    ClpSimplexProgress progress(this);
274   
275    // Say no pivot has occurred (for steepest edge and updates)
276    pivotRow_=-2;
277   
278    // This says whether to restore things etc
279    int factorType=0;
280    // Save iteration number
281    int saveNumber = -1;
282    /*
283      Status of problem:
284      0 - optimal
285      1 - infeasible
286      2 - unbounded
287      -1 - iterating
288      -2 - factorization wanted
289      -3 - redo checking without factorization
290      -4 - looks infeasible
291      -5 - looks unbounded
292    */
293    while (problemStatus_<0) {
294      // clear
295      for (iRow=0;iRow<4;iRow++) {
296        rowArray_[iRow]->clear();
297      }   
298     
299      for (iColumn=0;iColumn<2;iColumn++) {
300        columnArray_[iColumn]->clear();
301      }   
302     
303      // give matrix (and model costs and bounds a chance to be
304      // refreshed (normally null)
305      matrix_->refresh(this);
306      // If getting nowhere - why not give it a kick
307#if 0
308      // primal perturbation not coded yet
309      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
310        perturb();
311#endif
312      // If we have done no iterations - special
313      if (saveNumber==numberIterations_)
314        factorType=3;
315      // may factorize, checks if problem finished
316      statusOfProblemInPrimal(lastCleaned,factorType,progress);
317     
318      // Say good factorization
319      factorType=1;
320      if (saveSparse) {
321        // use default at present
322        factorization_->sparseThreshold(0);
323        factorization_->goSparse();
324      }
325     
326      // Say no pivot has occurred (for steepest edge and updates)
327      pivotRow_=-2;
328     
329      // Save iteration number
330      saveNumber = numberIterations_;
331     
332      // Iterate
333      whileIterating(firstSuperBasic);
334    }
335  }
336  // if infeasible get real values
337  if (problemStatus_==1) {
338    infeasibilityCost_=0.0;
339    createRim(7);
340    nonLinearCost_->checkInfeasibilities(true);
341    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
342    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
343    // and get good feasible duals
344    computeDuals();
345  }
346  factorization_->sparseThreshold(saveSparse);
347  // Get rid of some arrays and empty factorization
348  deleteRim();
349  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
350    <<objectiveValue()
351    <<CoinMessageEol;
352  // Restore any saved stuff
353  perturbation_ = savePerturbation;
354  infeasibilityCost_ = saveInfeasibilityCost;
355  return problemStatus_;
356}
357/*
358  Reasons to come out:
359  -1 iterations etc
360  -2 inaccuracy
361  -3 slight inaccuracy (and done iterations)
362  -4 end of values pass and done iterations
363  +0 looks optimal (might be infeasible - but we will investigate)
364  +2 looks unbounded
365  +3 max iterations
366*/
367int
368ClpSimplexPrimal::whileIterating(int & firstSuperBasic)
369{
370
371  // Say if values pass
372  int ifValuesPass=0;
373  int returnCode=-1;
374  int startIteration = numberIterations_;
375  if (firstSuperBasic<numberRows_+numberColumns_)
376    ifValuesPass=1;
377  int saveNumber = numberIterations_;
378  // status stays at -1 while iterating, >=0 finished, -2 to invert
379  // status -3 to go to top without an invert
380  while (problemStatus_==-1) {
381#ifdef CLP_DEBUG
382    {
383      int i;
384      // not [1] as has information
385      for (i=0;i<4;i++) {
386        if (i!=1)
387          rowArray_[i]->checkClear();
388      }   
389      for (i=0;i<2;i++) {
390        columnArray_[i]->checkClear();
391      }   
392    }     
393#endif
394#if CLP_DEBUG>2
395    // very expensive
396    if (numberIterations_>0&&numberIterations_<-2534) {
397      handler_->setLogLevel(63);
398      double saveValue = objectiveValue_;
399      double * saveRow1 = new double[numberRows_];
400      double * saveRow2 = new double[numberRows_];
401      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
402      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
403      double * saveColumn1 = new double[numberColumns_];
404      double * saveColumn2 = new double[numberColumns_];
405      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
406      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
407      createRim(7);
408      gutsOfSolution(rowActivityWork_,columnActivityWork_);
409      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
410             numberIterations_,
411             saveValue,objectiveValue_,sumPrimalInfeasibilities_);
412      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
413      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
414      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
415      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
416      delete [] saveRow1;
417      delete [] saveRow2;
418      delete [] saveColumn1;
419      delete [] saveColumn2;
420      objectiveValue_=saveValue;
421    }
422#endif
423    if (!ifValuesPass) {
424      // choose column to come in
425      // can use pivotRow_ to update weights
426      // pass in list of cost changes so can do row updates (rowArray_[1])
427      // NOTE rowArray_[0] is used by computeDuals which is a
428      // slow way of getting duals but might be used
429      primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
430                   columnArray_[0],columnArray_[1]);
431    } else {
432      // in values pass
433      if (ifValuesPass>0) {
434        nextSuperBasic(firstSuperBasic);
435        if (firstSuperBasic==numberRows_+numberColumns_)
436          ifValuesPass=-1; // signal end of values pass after this
437      } else {
438        // end of values pass - initialize weights etc
439        primalColumnPivot_->saveWeights(this,5);
440        ifValuesPass=0;
441        if(saveNumber != numberIterations_) {
442          problemStatus_=-2; // factorize now
443          pivotRow_=-1; // say no weights update
444          returnCode=-4;
445          break;
446        }
447
448        // and get variable
449        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
450                     columnArray_[0],columnArray_[1]);
451      }
452    }
453    pivotRow_=-1;
454    sequenceOut_=-1;
455    rowArray_[1]->clear();
456    if (sequenceIn_>=0) {
457      // we found a pivot column
458#ifdef CLP_DEBUG
459      if ((handler_->logLevel()&32)) {
460        char x = isColumn(sequenceIn_) ? 'C' :'R';
461        std::cout<<"pivot column "<<
462          x<<sequenceWithin(sequenceIn_)<<std::endl;
463      }
464#endif
465      // update the incoming column
466      unpack(rowArray_[1]);
467      // save reduced cost
468      double saveDj = dualIn_;
469      factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
470      // do ratio test and re-compute dj
471      primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
472                ifValuesPass);
473      if (ifValuesPass) {
474        saveDj=dualIn_;
475        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
476          if(fabs(dualIn_)<1.0e2*dualTolerance_) {
477            // try other way
478            directionIn_=-directionIn_;
479            primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
480                      0);
481          }
482          if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
483            // reject it
484            char x = isColumn(sequenceIn_) ? 'C' :'R';
485            handler_->message(CLP_SIMPLEX_FLAG,messages_)
486              <<x<<sequenceWithin(sequenceIn_)
487              <<CoinMessageEol;
488            setFlagged(sequenceIn_);
489            lastBadIteration_ = numberIterations_; // say be more cautious
490            rowArray_[1]->clear();
491            pivotRow_=-1;
492            continue;
493          }
494        }
495      }
496     
497#ifdef CLP_DEBUG
498      if ((handler_->logLevel()&32))
499        printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_);
500#endif
501      if ((saveDj*dualIn_<1.0e-20&&!ifValuesPass)||
502          fabs(saveDj-dualIn_)>1.0e-3*(1.0+fabs(saveDj))) {
503        handler_->message(CLP_PRIMAL_DJ,messages_)
504          <<saveDj<<dualIn_
505          <<CoinMessageEol;
506        if(saveNumber != numberIterations_) {
507          problemStatus_=-2; // factorize now
508          rowArray_[1]->clear();
509          pivotRow_=-1; // say no weights update
510          returnCode=-2;
511          if(saveNumber+1 == numberIterations_) {
512            // not looking wonderful - try cleaning bounds
513            // put non-basics to bounds in case tolerance moved
514            nonLinearCost_->checkInfeasibilities(true);
515          }
516          break;
517        } else {
518          // take on more relaxed criterion
519          if (saveDj*dualIn_<1.0e-20||
520              fabs(saveDj-dualIn_)>1.0e-1*(1.0+fabs(dualIn_))) {
521            // need to reject something
522            char x = isColumn(sequenceIn_) ? 'C' :'R';
523            handler_->message(CLP_SIMPLEX_FLAG,messages_)
524              <<x<<sequenceWithin(sequenceIn_)
525              <<CoinMessageEol;
526            setFlagged(sequenceIn_);
527            lastBadIteration_ = numberIterations_; // say be more cautious
528            rowArray_[1]->clear();
529            pivotRow_=-1;
530            continue;
531          }
532        }
533      }
534      if (pivotRow_>=0) {
535        // if stable replace in basis
536        int updateStatus = factorization_->replaceColumn(rowArray_[2],
537                                                         pivotRow_,
538                                                         alpha_);
539        if (updateStatus==1) {
540          // slight error
541          if (factorization_->pivots()>5) {
542            problemStatus_=-2; // factorize now
543            returnCode=-3;
544          }
545        } else if (updateStatus==2) {
546          // major error
547          // later we may need to unwind more e.g. fake bounds
548          if(saveNumber != numberIterations_) {
549            problemStatus_=-2; // factorize now
550            returnCode=-2;
551            break;
552          } else {
553            // need to reject something
554            char x = isColumn(sequenceIn_) ? 'C' :'R';
555            handler_->message(CLP_SIMPLEX_FLAG,messages_)
556              <<x<<sequenceWithin(sequenceIn_)
557              <<CoinMessageEol;
558            setFlagged(sequenceIn_);
559            lastBadIteration_ = numberIterations_; // say be more cautious
560            rowArray_[1]->clear();
561            pivotRow_=-1;
562            continue;
563          }
564        } else if (updateStatus==3) {
565          // out of memory
566          // increase space if not many iterations
567          if (factorization_->pivots()<
568              0.5*factorization_->maximumPivots()&&
569              factorization_->pivots()<200)
570            factorization_->areaFactor(
571                                       factorization_->areaFactor() * 1.1);
572          problemStatus_=-2; // factorize now
573        }
574        // here do part of steepest - ready for next iteration
575        primalColumnPivot_->updateWeights(rowArray_[1]);
576      } else {
577        if (pivotRow_==-1) {
578          // no outgoing row is valid
579#ifdef CLP_DEBUG
580          if (handler_->logLevel()&32)
581            printf("** no row pivot\n");
582#endif
583          if (!factorization_->pivots()) {
584            problemStatus_=-5; //say looks unbounded
585            // do ray
586            delete [] ray_;
587            ray_ = new double [numberColumns_];
588            ClpFillN(ray_,numberColumns_,0.0);
589            int number=rowArray_[1]->getNumElements();
590            int * index = rowArray_[1]->getIndices();
591            double * array = rowArray_[1]->denseVector();
592            double way=-directionIn_;
593            int i;
594            double zeroTolerance=1.0e-12;
595            if (sequenceIn_<numberColumns_)
596              ray_[sequenceIn_]=directionIn_;
597            for (i=0;i<number;i++) {
598              int iRow=index[i];
599              int iPivot=pivotVariable_[iRow];
600              double arrayValue = array[iRow];
601              if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
602                ray_[iPivot] = way* array[iRow];
603            }
604          }
605          rowArray_[0]->clear();
606          returnCode=2;
607          break;
608        } else {
609          // flipping from bound to bound
610        }
611      }
612     
613      // update primal solution
614     
615      double objectiveChange=0.0;
616      // Cost on pivot row may change - may need to change dualIn
617      double oldCost=0.0;
618      if (pivotRow_>=0)
619        oldCost = cost(pivotVariable_[pivotRow_]);
620      // rowArray_[1] is not empty - used to update djs
621      updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
622      if (pivotRow_>=0)
623        dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
624     
625      int whatNext=housekeeping(objectiveChange);
626     
627      if (whatNext==1) {
628        problemStatus_ =-2; // refactorize
629      } else if (whatNext==2) {
630        // maximum iterations or equivalent
631        problemStatus_= 3;
632        returnCode=3;
633        break;
634      } else if(numberIterations_ == startIteration
635                + 2 * factorization_->maximumPivots()) {
636        // done a lot of flips - be safe
637        problemStatus_ =-2; // refactorize
638      }
639    } else {
640      // no pivot column
641#ifdef CLP_DEBUG
642      if (handler_->logLevel()&32)
643        printf("** no column pivot\n");
644#endif
645      if (nonLinearCost_->numberInfeasibilities())
646        problemStatus_=-4; // might be infeasible
647      returnCode=0;
648      break;
649    }
650  }
651  return returnCode;
652}
653/* Checks if finished.  Updates status */
654void 
655ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type,
656                               ClpSimplexProgress &progress)
657{
658  if (type==2) {
659    // trouble - restore solution
660    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
661    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
662           numberRows_*sizeof(double));
663    memcpy(columnActivityWork_,savedSolution_ ,
664           numberColumns_*sizeof(double));
665    forceFactorization_=1; // a bit drastic but ..
666    pivotRow_=-1; // say no weights update
667    changeMade_++; // say change made
668  }
669  int tentativeStatus = problemStatus_;
670  if (problemStatus_>-3||problemStatus_==-4) {
671    // factorize
672    // later on we will need to recover from singularities
673    // also we could skip if first time
674    // do weights
675    // This may save pivotRow_ for use
676    primalColumnPivot_->saveWeights(this,1);
677
678    if (type) {
679      // is factorization okay?
680      if (internalFactorize(1)) {
681        // no - restore previous basis
682        assert (type==1);
683        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
684        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
685               numberRows_*sizeof(double));
686        memcpy(columnActivityWork_,savedSolution_ ,
687               numberColumns_*sizeof(double));
688        forceFactorization_=1; // a bit drastic but ..
689        type = 2;
690        assert (internalFactorize(1)==0);
691        changeMade_++; // say change made
692      }
693    }
694    if (problemStatus_!=-4)
695      problemStatus_=-3;
696  }
697  // at this stage status is -3 or -5 if looks unbounded
698  // get primal and dual solutions
699  // put back original bounds and then check
700  createRim(7);
701  gutsOfSolution(rowActivityWork_, columnActivityWork_);
702  // Check if looping
703  int loop = progress.looping();
704  if (loop>=0) {
705    problemStatus_ = loop; //exit if in loop
706    return ;
707  } else if (loop<-1) {
708    // something may have changed
709    gutsOfSolution(rowActivityWork_, columnActivityWork_);
710  }
711  progressFlag_ = 0; //reset progress flag
712
713  handler_->message(CLP_SIMPLEX_STATUS,messages_)
714    <<numberIterations_<<objectiveValue();
715  handler_->printing(sumPrimalInfeasibilities_>0.0)
716    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
717  handler_->printing(sumDualInfeasibilities_>0.0)
718    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
719  handler_->printing(numberDualInfeasibilitiesWithoutFree_
720                     <numberDualInfeasibilities_)
721                       <<numberDualInfeasibilities_-
722    numberDualInfeasibilitiesWithoutFree_;
723  handler_->message()<<CoinMessageEol;
724  assert (primalFeasible());
725  // we may wish to say it is optimal even if infeasible
726  bool alwaysOptimal = (specialOptions_&1)!=0;
727  // give code benefit of doubt
728  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
729      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
730    // say optimal (with these bounds etc)
731    numberDualInfeasibilities_ = 0;
732    sumDualInfeasibilities_ = 0.0;
733    numberPrimalInfeasibilities_ = 0;
734    sumPrimalInfeasibilities_ = 0.0;
735  }
736  if (dualFeasible()||problemStatus_==-4||(type==3&&problemStatus_!=-5)) {
737    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
738      //may need infeasiblity cost changed
739      // we can see if we can construct a ray
740      // make up a new objective
741      double saveWeight = infeasibilityCost_;
742      // save nonlinear cost as we are going to switch off costs
743      ClpNonLinearCost * nonLinear = nonLinearCost_;
744      // do twice to make sure Primal solution has settled
745      // put non-basics to bounds in case tolerance moved
746      // put back original bounds
747      createRim(7);
748      nonLinearCost_->checkInfeasibilities(true);
749      gutsOfSolution(rowActivityWork_, columnActivityWork_);
750
751      infeasibilityCost_=1.0e100;
752      // put back original bounds
753      createRim(7);
754      nonLinearCost_->checkInfeasibilities(true);
755      nonLinearCost_=NULL;
756      // scale
757      int i;
758      for (i=0;i<numberRows_+numberColumns_;i++) 
759        cost_[i] *= 1.0e-100;
760      gutsOfSolution(rowActivityWork_, columnActivityWork_);
761      nonLinearCost_=nonLinear;
762      infeasibilityCost_=saveWeight;
763      if (infeasibilityCost_>=1.0e20||
764          numberDualInfeasibilities_==0) {
765        // we are infeasible - use as ray
766        delete [] ray_;
767        ray_ = new double [numberRows_];
768        memcpy(ray_,dual_,numberRows_*sizeof(double));
769        // and get feasible duals
770        infeasibilityCost_=0.0;
771        createRim(7);
772        nonLinearCost_->checkInfeasibilities(true);
773        gutsOfSolution(rowActivityWork_, columnActivityWork_);
774        // so will exit
775        infeasibilityCost_=1.0e30;
776      }
777
778      if (infeasibilityCost_<1.0e20) {
779        infeasibilityCost_ *= 5.0;
780        changeMade_++; // say change made
781        handler_->message(CLP_PRIMAL_WEIGHT,messages_)
782          <<infeasibilityCost_
783          <<CoinMessageEol;
784        // put back original bounds and then check
785        createRim(7);
786        nonLinearCost_->checkInfeasibilities(true);
787        gutsOfSolution(rowActivityWork_, columnActivityWork_);
788        problemStatus_=-1; //continue
789      } else {
790        // say infeasible
791        problemStatus_ = 1;
792      }
793    } else {
794      // may be optimal
795      if ( lastCleaned!=numberIterations_) {
796        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
797          <<primalTolerance_
798          <<CoinMessageEol;
799        if (numberTimesOptimal_<4) {
800          numberTimesOptimal_++;
801          changeMade_++; // say change made
802          if (numberTimesOptimal_==1) {
803            // better to have small tolerance even if slower
804            factorization_->zeroTolerance(1.0e-15);
805          }
806          lastCleaned=numberIterations_;
807          handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
808            <<CoinMessageEol;
809          primalTolerance_=dblParam_[ClpPrimalTolerance];
810         
811          // put back original bounds and then check
812          createRim(7);
813          nonLinearCost_->checkInfeasibilities(true);
814          gutsOfSolution(rowActivityWork_, columnActivityWork_);
815          problemStatus_ = -1;
816        } else {
817          problemStatus_=0; // optimal
818          if (lastCleaned<numberIterations_) {
819            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
820              <<CoinMessageEol;
821          }
822        }
823      } else {
824        problemStatus_=0; // optimal
825      }
826    }
827  } else {
828    // see if looks unbounded
829    if (problemStatus_==-5) {
830      if (nonLinearCost_->numberInfeasibilities()) {
831        //we need infeasiblity cost changed
832        if (infeasibilityCost_<1.0e20) {
833          infeasibilityCost_ *= 5.0;
834          changeMade_++; // say change made
835          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
836            <<infeasibilityCost_
837            <<CoinMessageEol;
838          // put back original bounds and then check
839          createRim(7);
840          gutsOfSolution(rowActivityWork_, columnActivityWork_);
841          problemStatus_=-1; //continue
842        } else {
843          // say unbounded
844          problemStatus_ = 2;
845        }
846      } else {
847        // say unbounded
848        problemStatus_ = 2;
849      } 
850    } else {
851      // carry on
852      problemStatus_ = -1;
853    }
854  }
855  if (type==0||type==1) {
856    if (!type) {
857      // create save arrays
858      delete [] saveStatus_;
859      delete [] savedSolution_;
860      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
861      savedSolution_ = new double [numberRows_+numberColumns_];
862    }
863    // save arrays
864    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
865    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
866           numberRows_*sizeof(double));
867    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
868  }
869  // restore weights (if saved) - also recompute infeasibility list
870  if (tentativeStatus>-3) 
871    primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
872  else
873    primalColumnPivot_->saveWeights(this,3);
874  if (problemStatus_<0&&!changeMade_) {
875    problemStatus_=4; // unknown
876  }
877}
878/*
879   Row array has pivot column
880   This chooses pivot row.
881   For speed, we may need to go to a bucket approach when many
882   variables go through bounds
883   On exit rhsArray will have changes in costs of basic variables
884*/
885void 
886ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
887                            CoinIndexedVector * rhsArray,
888                            CoinIndexedVector * spareArray,
889                            CoinIndexedVector * spareArray2,
890                            int valuesPass)
891{
892  if (valuesPass) {
893    dualIn_ = cost_[sequenceIn_];
894
895    double * work=rowArray->denseVector();
896    int number=rowArray->getNumElements();
897    int * which=rowArray->getIndices();
898
899    int iIndex;
900
901    for (iIndex=0;iIndex<number;iIndex++) {
902     
903      int iRow = which[iIndex];
904      double alpha = work[iRow];
905      int iPivot=pivotVariable_[iRow];
906      dualIn_ -= alpha*cost(iPivot);
907    }
908    // determine direction here
909    if (dualIn_<-dualTolerance_) {
910      directionIn_=1;
911    } else if (dualIn_>dualTolerance_) {
912      directionIn_=-1;
913    } else {
914      // towards nearest bound
915      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
916        directionIn_=-1;
917        dualIn_=dualTolerance_;
918      } else {
919        directionIn_=1;
920        dualIn_=-dualTolerance_;
921      }
922    }
923  }
924
925  // sequence stays as row number until end
926  pivotRow_=-1;
927  int numberSwapped=0;
928  int numberRemaining=0;
929
930  int numberThru =0; // number gone thru a barrier
931  int lastThru =0; // number gone thru a barrier on last time
932 
933  double totalThru=0.0; // for when variables flip
934  double acceptablePivot=1.0e-7;
935  if (factorization_->pivots())
936    acceptablePivot=1.0e-5; // if we have iterated be more strict
937  double bestEverPivot=acceptablePivot;
938  int lastPivotRow = -1;
939  double lastPivot=0.0;
940  double lastTheta=1.0e50;
941  int lastNumberSwapped=0;
942
943  // use spareArrays to put ones looked at in
944  // First one is list of candidates
945  // We could compress if we really know we won't need any more
946  // Second array has current set of pivot candidates
947  // with a backup list saved in double * part of indexed vector
948
949  // for zeroing out arrays after
950  int maximumSwapped=0;
951  // pivot elements
952  double * spare;
953  // indices
954  int * index, * indexSwapped;
955  int * saveSwapped;
956  spareArray->clear();
957  spareArray2->clear();
958  spare = spareArray->denseVector();
959  index = spareArray->getIndices();
960  saveSwapped = (int *) spareArray2->denseVector();
961  indexSwapped = spareArray2->getIndices();
962
963  // we also need somewhere for effective rhs
964  double * rhs=rhsArray->denseVector();
965
966  /*
967    First we get a list of possible pivots.  We can also see if the
968    problem looks unbounded.
969
970    At first we increase theta and see what happens.  We start
971    theta at a reasonable guess.  If in right area then we do bit by bit.
972    We save possible pivot candidates
973
974   */
975
976  // do first pass to get possibles
977  // We can also see if unbounded
978  // We also re-compute reduced cost
979
980  dualIn_ = cost_[sequenceIn_];
981
982  double * work=rowArray->denseVector();
983  int number=rowArray->getNumElements();
984  int * which=rowArray->getIndices();
985
986  // we need to swap sign if coming in from ub
987  double way = directionIn_;
988  double maximumMovement;
989  if (way>0.0) 
990    maximumMovement = min(1.0e30,upperIn_-valueIn_);
991  else
992    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
993
994  double tentativeTheta = maximumMovement;
995  double upperTheta = maximumMovement;
996
997  int iIndex;
998
999  for (iIndex=0;iIndex<number;iIndex++) {
1000
1001    int iRow = which[iIndex];
1002    double alpha = work[iRow];
1003    int iPivot=pivotVariable_[iRow];
1004    dualIn_ -= alpha*cost(iPivot);
1005    alpha *= way;
1006    double oldValue = solution(iPivot);
1007    // get where in bound sequence
1008    if (alpha>0.0) {
1009      // basic variable going towards lower bound
1010      double bound = lower(iPivot);
1011      oldValue -= bound;
1012    } else if (alpha<0.0) {
1013      // basic variable going towards upper bound
1014      double bound = upper(iPivot);
1015      oldValue = bound-oldValue;
1016    }
1017    double value = oldValue-tentativeTheta*fabs(alpha);
1018    assert (oldValue>=-primalTolerance_*1.0001);
1019    if (value<-primalTolerance_) {
1020      // add to list
1021      spare[numberRemaining]=alpha;
1022      rhs[iRow]=oldValue;
1023      index[numberRemaining++]=iRow;
1024      double value=oldValue-upperTheta*fabs(alpha);
1025      if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
1026        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1027    }
1028  }
1029
1030  // we need to keep where rhs non-zeros are
1031  int numberInRhs=numberRemaining;
1032  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
1033  rhsArray->setNumElements(numberInRhs);
1034
1035  theta_=maximumMovement;
1036
1037  double dualCheck = fabs(dualIn_);
1038  // but make a bit more pessimistic
1039  dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1040
1041  bool goBackOne = false;
1042
1043  if (numberRemaining) {
1044
1045    // looks like pivoting
1046    // now try until reasonable theta
1047    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1048    tentativeTheta = min(tentativeTheta,maximumMovement);
1049   
1050    // loops increasing tentative theta until can't go through
1051   
1052    while (tentativeTheta <= maximumMovement) {
1053      double thruThis = 0.0;
1054     
1055      double bestPivot=acceptablePivot;
1056      pivotRow_ = -1;
1057     
1058      numberSwapped = 0;
1059     
1060      upperTheta = maximumMovement;
1061     
1062      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1063
1064        int iRow = index[iIndex];
1065        double alpha = spare[iIndex];
1066        double oldValue = rhs[iRow];
1067        double value = oldValue-tentativeTheta*fabs(alpha);
1068
1069        if (value<-primalTolerance_) {
1070          // how much would it cost to go thru
1071          thruThis += alpha*
1072            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1073          // goes on swapped list (also means candidates if too many)
1074          indexSwapped[numberSwapped++]=iRow;
1075          if (fabs(alpha)>bestPivot) {
1076            bestPivot=fabs(alpha);
1077            pivotRow_ = iRow;
1078            theta_ = oldValue/bestPivot;
1079          }
1080        } else {
1081          value = oldValue-upperTheta*fabs(alpha);
1082          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1083            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1084        }
1085      }
1086     
1087      maximumSwapped = max(maximumSwapped,numberSwapped);
1088
1089      if (totalThru+thruThis>=dualCheck) {
1090        // We should be pivoting in this batch
1091        // so compress down to this lot
1092
1093        int saveNumber = numberRemaining;
1094        numberRemaining=0;
1095        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1096          int iRow = indexSwapped[iIndex];
1097          spare[numberRemaining]=way*work[iRow];
1098          index[numberRemaining++]=iRow;
1099        }
1100        memset(spare+numberRemaining,0,
1101               (saveNumber-numberRemaining)*sizeof(double));
1102        int iTry;
1103#define MAXTRY 100
1104        // first get ratio with tolerance
1105        for (iTry=0;iTry<MAXTRY;iTry++) {
1106         
1107          upperTheta=maximumMovement;
1108          numberSwapped = 0;
1109         
1110          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1111           
1112            int iRow = index[iIndex];
1113            double alpha = fabs(spare[iIndex]);
1114            double oldValue = rhs[iRow];
1115            double value = oldValue-upperTheta*alpha;
1116           
1117            if (value<-primalTolerance_ && alpha>=acceptablePivot) 
1118              upperTheta = (oldValue+primalTolerance_)/alpha;
1119           
1120          }
1121         
1122          // now look at best in this lot
1123          bestPivot=acceptablePivot;
1124          pivotRow_=-1;
1125          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1126           
1127            int iRow = index[iIndex];
1128            double alpha = spare[iIndex];
1129            double oldValue = rhs[iRow];
1130            double value = oldValue-upperTheta*fabs(alpha);
1131           
1132            if (value<=0) {
1133              // how much would it cost to go thru
1134              totalThru += alpha*
1135                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1136              // goes on swapped list (also means candidates if too many)
1137              indexSwapped[numberSwapped++]=iRow;
1138              if (fabs(alpha)>bestPivot) {
1139                bestPivot=fabs(alpha);
1140                theta_ = oldValue/bestPivot;
1141                pivotRow_=iRow;
1142              }
1143            } else {
1144              value = oldValue-upperTheta*fabs(alpha);
1145              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1146                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1147            }
1148          }
1149         
1150          maximumSwapped = max(maximumSwapped,numberSwapped);
1151          if (bestPivot<0.1*bestEverPivot&&
1152              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1153            // back to previous one
1154            goBackOne = true;
1155            break;
1156          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1157            if (lastPivot>acceptablePivot) {
1158              // back to previous one
1159              goBackOne = true;
1160              //break;
1161            } else {
1162              // can only get here if all pivots so far too small
1163            }
1164            break;
1165          } else if (totalThru>=dualCheck) {
1166            break; // no point trying another loop
1167          } else {
1168            // skip this lot
1169            nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1170            lastPivotRow=pivotRow_;
1171            lastTheta = theta_;
1172            lastThru = numberThru;
1173            numberThru += numberSwapped;
1174            lastNumberSwapped = numberSwapped;
1175            memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1176            if (bestPivot>bestEverPivot)
1177              bestEverPivot=bestPivot;
1178          }
1179        }
1180        break;
1181      } else {
1182        // skip this lot
1183        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1184        lastPivotRow=pivotRow_;
1185        lastTheta = theta_;
1186        lastThru = numberThru;
1187        numberThru += numberSwapped;
1188        lastNumberSwapped = numberSwapped;
1189        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1190        if (bestPivot>bestEverPivot)
1191          bestEverPivot=bestPivot;
1192        totalThru += thruThis;
1193        tentativeTheta = 2.0*upperTheta;
1194      }
1195    }
1196    // can get here without pivotRow_ set but with lastPivotRow
1197    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1198      // back to previous one
1199      pivotRow_=lastPivotRow;
1200      theta_ = lastTheta;
1201            // undo this lot
1202      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1203      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1204      numberSwapped = lastNumberSwapped;
1205    }
1206  }
1207
1208  if (pivotRow_>=0) {
1209   
1210#define MINIMUMTHETA 1.0e-12
1211    // will we need to increase tolerance
1212#ifdef CLP_DEBUG
1213    bool found=false;
1214#endif
1215    double largestInfeasibility = primalTolerance_;
1216    if (theta_<MINIMUMTHETA) {
1217      theta_=MINIMUMTHETA;
1218      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1219        int iRow = indexSwapped[iIndex];
1220#ifdef CLP_DEBUG
1221        if (iRow==pivotRow_)
1222          found=true;
1223#endif
1224        largestInfeasibility = max (largestInfeasibility,
1225                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
1226      }
1227#ifdef CLP_DEBUG
1228      assert(found);
1229      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
1230        printf("Primal tolerance increased from %g to %g\n",
1231               primalTolerance_,largestInfeasibility);
1232#endif
1233      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1234    }
1235    alpha_ = work[pivotRow_];
1236    // translate to sequence
1237    sequenceOut_ = pivotVariable_[pivotRow_];
1238    valueOut_ = solution(sequenceOut_);
1239    lowerOut_=lower_[sequenceOut_];
1240    upperOut_=upper_[sequenceOut_];
1241
1242    if (way<0.0) 
1243      theta_ = - theta_;
1244    double newValue = valueOut_ - theta_*alpha_;
1245    if (alpha_*way<0.0) {
1246      directionOut_=-1;      // to upper bound
1247      if (fabs(theta_)>1.0e-6)
1248        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1249      else
1250        upperOut_ = newValue;
1251    } else {
1252      directionOut_=1;      // to lower bound
1253      if (fabs(theta_)>1.0e-6)
1254        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1255      else
1256        lowerOut_ = newValue;
1257    }
1258    dualOut_ = reducedCost(sequenceOut_);
1259  } else if (maximumMovement<1.0e20) {
1260    // flip
1261    pivotRow_ = -2; // so we can tell its a flip
1262    sequenceOut_ = sequenceIn_;
1263    valueOut_ = valueIn_;
1264    dualOut_ = dualIn_;
1265    lowerOut_ = lowerIn_;
1266    upperOut_ = upperIn_;
1267    alpha_ = 0.0;
1268    if (way<0.0) {
1269      directionOut_=1;      // to lower bound
1270      theta_ = lowerOut_ - valueOut_;
1271    } else {
1272      directionOut_=-1;      // to upper bound
1273      theta_ = upperOut_ - valueOut_;
1274    }
1275  }
1276
1277  // clear arrays
1278
1279  memset(spare,0,numberRemaining*sizeof(double));
1280  memset(saveSwapped,0,maximumSwapped*sizeof(int));
1281
1282  // put back original bounds etc
1283  nonLinearCost_->goBackAll(rhsArray);
1284
1285  rhsArray->clear();
1286
1287}
1288/*
1289   Chooses primal pivot column
1290   updateArray has cost updates (also use pivotRow_ from last iteration)
1291   Would be faster with separate region to scan
1292   and will have this (with square of infeasibility) when steepest
1293   For easy problems we can just choose one of the first columns we look at
1294*/
1295void 
1296ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
1297                               CoinIndexedVector * spareRow1,
1298                               CoinIndexedVector * spareRow2,
1299                               CoinIndexedVector * spareColumn1,
1300                               CoinIndexedVector * spareColumn2)
1301{
1302  sequenceIn_ = primalColumnPivot_->pivotColumn(updates,spareRow1,
1303                                               spareRow2,spareColumn1,
1304                                               spareColumn2);
1305  if (sequenceIn_>=0) {
1306    valueIn_=solution_[sequenceIn_];
1307    lowerIn_=lower_[sequenceIn_];
1308    upperIn_=upper_[sequenceIn_];
1309    dualIn_=dj_[sequenceIn_];
1310    if (dualIn_>0.0)
1311      directionIn_ = -1;
1312    else 
1313      directionIn_ = 1;
1314  } else {
1315    sequenceIn_ = -1;
1316  }
1317}
1318/* The primals are updated by the given array.
1319   Returns number of infeasibilities.
1320   After rowArray will have list of cost changes
1321*/
1322int 
1323ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
1324                  double theta,
1325                  double & objectiveChange)
1326{
1327  double * work=rowArray->denseVector();
1328  int number=rowArray->getNumElements();
1329  int * which=rowArray->getIndices();
1330
1331  int newNumber = 0;
1332
1333  nonLinearCost_->setChangeInCost(0.0);
1334  int iIndex;
1335
1336  for (iIndex=0;iIndex<number;iIndex++) {
1337
1338    int iRow = which[iIndex];
1339    double alpha = work[iRow];
1340    int iPivot=pivotVariable_[iRow];
1341    double & value = solutionAddress(iPivot);
1342    double change = theta*alpha;
1343    value -= change;
1344
1345    if (change>0.0) {
1346      // going down
1347      if (value<=lower(iPivot)+primalTolerance_) {
1348        double difference = nonLinearCost_->setOne(iPivot,value);
1349        work[iRow] = difference;
1350        if (difference) {
1351          //change reduced cost on this
1352          reducedCostAddress(iPivot) = -difference;
1353          which[newNumber++]=iRow;
1354        }
1355      } else {
1356        work[iRow]=0.0;
1357      }
1358    } else {
1359      // going up
1360      if (value>=upper(iPivot)-primalTolerance_) {
1361        double difference = nonLinearCost_->setOne(iPivot,value);
1362        work[iRow] = difference;
1363        if (difference) {
1364          //change reduced cost on this
1365          reducedCostAddress(iPivot) = -difference;
1366          which[newNumber++]=iRow;
1367        }
1368      } else {
1369        work[iRow]=0.0;
1370      }
1371    }
1372  }
1373  objectiveChange += nonLinearCost_->changeInCost();
1374  rowArray->setNumElements(newNumber);
1375  return 0;
1376}
1377void
1378ClpSimplexPrimal::nextSuperBasic(int & firstSuperBasic)
1379{
1380  int iColumn;
1381  if (firstSuperBasic==numberRows_+numberColumns_) {
1382    // initialization
1383    iColumn=0;
1384  } else {
1385    // normal
1386    sequenceIn_=firstSuperBasic;
1387    valueIn_=solution_[sequenceIn_];
1388    lowerIn_=lower_[sequenceIn_];
1389    upperIn_=upper_[sequenceIn_];
1390    dualIn_=dj_[sequenceIn_];
1391    iColumn=firstSuperBasic+1;
1392  }
1393  for (;iColumn<numberRows_+numberColumns_;iColumn++) {
1394    if (getStatus(iColumn)==superBasic) {
1395      // is it really super basic
1396      if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
1397        solution_[iColumn]=lower_[iColumn];
1398        setStatus(iColumn,atLowerBound);
1399      } else if (fabs(solution_[iColumn]-upper_[iColumn])
1400                 <=primalTolerance_) {
1401        solution_[iColumn]=upper_[iColumn];
1402        setStatus(iColumn,atUpperBound);
1403      } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
1404        setStatus(iColumn,isFree);
1405      } else {
1406        break;
1407      }
1408    }
1409  }
1410  firstSuperBasic = iColumn;
1411}
1412// Perturbs problem
1413void 
1414ClpSimplexPrimal::perturb()
1415{
1416  if (perturbation_>100)
1417    return; //perturbed already
1418  abort();
1419}
1420// Do not change infeasibility cost and always say optimal
1421void 
1422ClpSimplexPrimal::alwaysOptimal(bool onOff)
1423{
1424  if (onOff)
1425    specialOptions_ |= 1;
1426  else
1427    specialOptions_ &= ~1;
1428}
1429
Note: See TracBrowser for help on using the repository browser.