source: branches/pre/ClpSimplexPrimal.cpp @ 222

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

Start of mini sprint

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