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

Last change on this file since 26 was 26, checked in by forrest, 18 years ago

Taking out dependence on CoinWarmStartBasis?

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