source: branches/pre/ClpSimplexPrimal.cpp @ 212

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

lots of stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 63.2 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{
436  if (type==2) {
437    // trouble - restore solution
438    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
439    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
440           numberRows_*sizeof(double));
441    memcpy(columnActivityWork_,savedSolution_ ,
442           numberColumns_*sizeof(double));
443    forceFactorization_=1; // a bit drastic but ..
444    pivotRow_=-1; // say no weights update
445    changeMade_++; // say change made
446  }
447  int saveThreshold = factorization_->sparseThreshold();
448  int tentativeStatus = problemStatus_;
449  if (problemStatus_>-3||problemStatus_==-4) {
450    // factorize
451    // later on we will need to recover from singularities
452    // also we could skip if first time
453    // do weights
454    // This may save pivotRow_ for use
455    primalColumnPivot_->saveWeights(this,1);
456
457    if (type) {
458      // is factorization okay?
459      if (internalFactorize(1)) {
460        if (solveType_==2) {
461          // say odd
462          problemStatus_=5;
463          return;
464        }
465#if 1
466        // switch off dense
467        int saveDense = factorization_->denseThreshold();
468        factorization_->setDenseThreshold(0);
469        internalFactorize(2);
470        factorization_->setDenseThreshold(saveDense);
471#else
472
473        // no - restore previous basis
474        assert (type==1);
475        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
476        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
477               numberRows_*sizeof(double));
478        memcpy(columnActivityWork_,savedSolution_ ,
479               numberColumns_*sizeof(double));
480        forceFactorization_=1; // a bit drastic but ..
481        type = 2;
482        assert (internalFactorize(1)==0);
483#endif
484        changeMade_++; // say change made
485      }
486    }
487    if (problemStatus_!=-4)
488      problemStatus_=-3;
489  }
490  // at this stage status is -3 or -5 if looks unbounded
491  // get primal and dual solutions
492  // put back original costs and then check
493  createRim(4);
494  gutsOfSolution(NULL,NULL,(firstFree_>=0));
495  // Double check reduced costs if no action
496  if (progress->lastIterationNumber(0)==numberIterations_) {
497    if (primalColumnPivot_->looksOptimal()) {
498      numberDualInfeasibilities_ = 0;
499      sumDualInfeasibilities_ = 0.0;
500    }
501  }
502  // Check if looping
503  int loop;
504  if (type!=2) 
505    loop = progress->looping();
506  else
507    loop=-1;
508  if (loop>=0) {
509    problemStatus_ = loop; //exit if in loop
510    return ;
511  } else if (loop<-1) {
512    // something may have changed
513    gutsOfSolution(NULL,NULL);
514  }
515  // Flag to say whether to go to dual to clean up
516  bool goToDual=false;
517  // really for free variables in
518  //if((progressFlag_&2)!=0)
519  //problemStatus_=-1;;
520  progressFlag_ = 0; //reset progress flag
521
522  handler_->message(CLP_SIMPLEX_STATUS,messages_)
523    <<numberIterations_<<objectiveValue();
524  handler_->printing(sumPrimalInfeasibilities_>0.0)
525    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
526  handler_->printing(sumDualInfeasibilities_>0.0)
527    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
528  handler_->printing(numberDualInfeasibilitiesWithoutFree_
529                     <numberDualInfeasibilities_)
530                       <<numberDualInfeasibilitiesWithoutFree_;
531  handler_->message()<<CoinMessageEol;
532  if (!primalFeasible()) {
533    nonLinearCost_->checkInfeasibilities(primalTolerance_);
534    gutsOfSolution(NULL,NULL);
535  }
536  // we may wish to say it is optimal even if infeasible
537  bool alwaysOptimal = (specialOptions_&1)!=0;
538  // give code benefit of doubt
539  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
540      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
541    // say optimal (with these bounds etc)
542    numberDualInfeasibilities_ = 0;
543    sumDualInfeasibilities_ = 0.0;
544    numberPrimalInfeasibilities_ = 0;
545    sumPrimalInfeasibilities_ = 0.0;
546  }
547  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
548  if (dualFeasible()||problemStatus_==-4) {
549    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
550      //may need infeasiblity cost changed
551      // we can see if we can construct a ray
552      // make up a new objective
553      double saveWeight = infeasibilityCost_;
554      // save nonlinear cost as we are going to switch off costs
555      ClpNonLinearCost * nonLinear = nonLinearCost_;
556      // do twice to make sure Primal solution has settled
557      // put non-basics to bounds in case tolerance moved
558      // put back original costs
559      createRim(4);
560      nonLinearCost_->checkInfeasibilities(primalTolerance_);
561      gutsOfSolution(NULL,NULL);
562
563      infeasibilityCost_=1.0e100;
564      // put back original costs
565      createRim(4);
566      nonLinearCost_->checkInfeasibilities(primalTolerance_);
567      // may have fixed infeasibilities - double check
568      if (nonLinearCost_->numberInfeasibilities()==0) {
569        // carry on
570        problemStatus_ = -1;
571        infeasibilityCost_=saveWeight;
572      } else {
573        nonLinearCost_=NULL;
574        // scale
575        int i;
576        for (i=0;i<numberRows_+numberColumns_;i++) 
577          cost_[i] *= 1.0e-100;
578        gutsOfSolution(NULL,NULL);
579        nonLinearCost_=nonLinear;
580        infeasibilityCost_=saveWeight;
581        if ((infeasibilityCost_>=1.0e18||
582             numberDualInfeasibilities_==0)&&perturbation_==101) {
583          goToDual=unPerturb(); // stop any further perturbation
584          numberDualInfeasibilities_=1; // carry on
585        }
586        if (infeasibilityCost_>=1.0e20||
587            numberDualInfeasibilities_==0) {
588          // we are infeasible - use as ray
589          delete [] ray_;
590          ray_ = new double [numberRows_];
591          memcpy(ray_,dual_,numberRows_*sizeof(double));
592          // and get feasible duals
593          infeasibilityCost_=0.0;
594          createRim(4);
595          nonLinearCost_->checkInfeasibilities(primalTolerance_);
596          gutsOfSolution(NULL,NULL);
597          // so will exit
598          infeasibilityCost_=1.0e30;
599          // reset infeasibilities
600          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
601          numberPrimalInfeasibilities_=
602            nonLinearCost_->numberInfeasibilities();
603        }
604        if (infeasibilityCost_<1.0e20) {
605          infeasibilityCost_ *= 5.0;
606          changeMade_++; // say change made
607          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
608            <<infeasibilityCost_
609            <<CoinMessageEol;
610          // put back original costs and then check
611          createRim(4);
612          nonLinearCost_->checkInfeasibilities(0.0);
613          gutsOfSolution(NULL,NULL);
614          problemStatus_=-1; //continue
615          goToDual=false;
616        } else {
617          // say infeasible
618          problemStatus_ = 1;
619        }
620      }
621    } else {
622      // may be optimal
623      if (perturbation_==101) {
624        goToDual=unPerturb(); // stop any further perturbation
625        lastCleaned=-1; // carry on
626      }
627      bool unflagged = unflag();
628      if ( lastCleaned!=numberIterations_||unflagged) {
629        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
630          <<primalTolerance_
631          <<CoinMessageEol;
632        if (numberTimesOptimal_<4) {
633          numberTimesOptimal_++;
634          changeMade_++; // say change made
635          if (numberTimesOptimal_==1) {
636            // better to have small tolerance even if slower
637            factorization_->zeroTolerance(1.0e-15);
638          }
639          lastCleaned=numberIterations_;
640          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
641            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
642              <<CoinMessageEol;
643          double oldTolerance = primalTolerance_;
644          primalTolerance_=dblParam_[ClpPrimalTolerance];
645#if 0
646          double * xcost = new double[numberRows_+numberColumns_];
647          double * xlower = new double[numberRows_+numberColumns_];
648          double * xupper = new double[numberRows_+numberColumns_];
649          double * xdj = new double[numberRows_+numberColumns_];
650          double * xsolution = new double[numberRows_+numberColumns_];
651          memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
652          memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
653          memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
654          memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
655          memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
656#endif
657          // put back original costs and then check
658          createRim(4);
659          nonLinearCost_->checkInfeasibilities(oldTolerance);
660#if 0
661          int i;
662          for (i=0;i<numberRows_+numberColumns_;i++) {
663            if (cost_[i]!=xcost[i])
664              printf("** %d old cost %g new %g sol %g\n",
665                     i,xcost[i],cost_[i],solution_[i]);
666            if (lower_[i]!=xlower[i])
667              printf("** %d old lower %g new %g sol %g\n",
668                     i,xlower[i],lower_[i],solution_[i]);
669            if (upper_[i]!=xupper[i])
670              printf("** %d old upper %g new %g sol %g\n",
671                     i,xupper[i],upper_[i],solution_[i]);
672            if (dj_[i]!=xdj[i])
673              printf("** %d old dj %g new %g sol %g\n",
674                     i,xdj[i],dj_[i],solution_[i]);
675            if (solution_[i]!=xsolution[i])
676              printf("** %d old solution %g new %g sol %g\n",
677                     i,xsolution[i],solution_[i],solution_[i]);
678          }
679          delete [] xcost;
680          delete [] xupper;
681          delete [] xlower;
682          delete [] xdj;
683          delete [] xsolution;
684#endif
685          gutsOfSolution(NULL,NULL);
686          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
687              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
688            // say optimal (with these bounds etc)
689            numberDualInfeasibilities_ = 0;
690            sumDualInfeasibilities_ = 0.0;
691            numberPrimalInfeasibilities_ = 0;
692            sumPrimalInfeasibilities_ = 0.0;
693          }
694          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
695            problemStatus_=0;
696          else
697            problemStatus_ = -1;
698        } else {
699          problemStatus_=0; // optimal
700          if (lastCleaned<numberIterations_) {
701            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
702              <<CoinMessageEol;
703          }
704        }
705      } else {
706        problemStatus_=0; // optimal
707      }
708    }
709  } else {
710    // see if looks unbounded
711    if (problemStatus_==-5) {
712      if (nonLinearCost_->numberInfeasibilities()) {
713        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
714          // back off weight
715          infeasibilityCost_ = 1.0e13;
716          goToDual=unPerturb(); // stop any further perturbation
717        }
718        //we need infeasiblity cost changed
719        if (infeasibilityCost_<1.0e20) {
720          infeasibilityCost_ *= 5.0;
721          changeMade_++; // say change made
722          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
723            <<infeasibilityCost_
724            <<CoinMessageEol;
725          // put back original costs and then check
726          createRim(4);
727          gutsOfSolution(NULL,NULL);
728          problemStatus_=-1; //continue
729        } else {
730          // say unbounded
731          problemStatus_ = 2;
732        }
733      } else {
734        // say unbounded
735        problemStatus_ = 2;
736      } 
737    } else {
738      if(type==3&&problemStatus_!=-5)
739        unflag(); // odd
740      // carry on
741      problemStatus_ = -1;
742    }
743  }
744  if (type==0||type==1) {
745    if (type!=1||!saveStatus_) {
746      // create save arrays
747      delete [] saveStatus_;
748      delete [] savedSolution_;
749      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
750      savedSolution_ = new double [numberRows_+numberColumns_];
751    }
752    // save arrays
753    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
754    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
755           numberRows_*sizeof(double));
756    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
757  }
758  // restore weights (if saved) - also recompute infeasibility list
759  if (tentativeStatus>-3) 
760    primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
761  else
762    primalColumnPivot_->saveWeights(this,3);
763  if (problemStatus_<0&&!changeMade_) {
764    problemStatus_=4; // unknown
765  }
766  if (saveThreshold) {
767    // use default at present
768    factorization_->sparseThreshold(0);
769    factorization_->goSparse();
770  }
771  lastGoodIteration_ = numberIterations_;
772  if (goToDual) 
773    problemStatus_=10; // try dual
774#if 0
775  double thisObj = progress->lastObjective(0);
776  double lastObj = progress->lastObjective(1);
777  if (lastObj<thisObj-1.0e-7*max(fabs(thisObj),fabs(lastObj))-1.0e-8
778      &&firstFree_<0) {
779    int maxFactor = factorization_->maximumPivots();
780    if (maxFactor>10) {
781      if (forceFactorization_<0)
782        forceFactorization_= maxFactor;
783      forceFactorization_ = max (1,(forceFactorization_>>1));
784      printf("Reducing factorization frequency\n");
785    }
786  }
787#endif
788}
789/*
790   Row array has pivot column
791   This chooses pivot row.
792   For speed, we may need to go to a bucket approach when many
793   variables go through bounds
794   On exit rhsArray will have changes in costs of basic variables
795*/
796void 
797ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
798                            CoinIndexedVector * rhsArray,
799                            CoinIndexedVector * spareArray,
800                            CoinIndexedVector * spareArray2,
801                            int valuesPass)
802{
803  if (valuesPass) {
804    dualIn_ = cost_[sequenceIn_];
805
806    double * work=rowArray->denseVector();
807    int number=rowArray->getNumElements();
808    int * which=rowArray->getIndices();
809
810    int iIndex;
811    for (iIndex=0;iIndex<number;iIndex++) {
812     
813      int iRow = which[iIndex];
814      double alpha = work[iRow];
815      int iPivot=pivotVariable_[iRow];
816      dualIn_ -= alpha*cost(iPivot);
817    }
818    // determine direction here
819    if (dualIn_<-dualTolerance_) {
820      directionIn_=1;
821    } else if (dualIn_>dualTolerance_) {
822      directionIn_=-1;
823    } else {
824      // towards nearest bound
825      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
826        directionIn_=-1;
827        dualIn_=dualTolerance_;
828      } else {
829        directionIn_=1;
830        dualIn_=-dualTolerance_;
831      }
832    }
833  }
834
835  // sequence stays as row number until end
836  pivotRow_=-1;
837  int numberSwapped=0;
838  int numberRemaining=0;
839
840  int numberThru =0; // number gone thru a barrier
841  int lastThru =0; // number gone thru a barrier on last time
842 
843  double totalThru=0.0; // for when variables flip
844  double acceptablePivot=1.0e-7;
845  if (factorization_->pivots())
846    acceptablePivot=1.0e-5; // if we have iterated be more strict
847  double bestEverPivot=acceptablePivot;
848  int lastPivotRow = -1;
849  double lastPivot=0.0;
850  double lastTheta=1.0e50;
851  int lastNumberSwapped=0;
852
853  // use spareArrays to put ones looked at in
854  // First one is list of candidates
855  // We could compress if we really know we won't need any more
856  // Second array has current set of pivot candidates
857  // with a backup list saved in double * part of indexed vector
858
859  // for zeroing out arrays after
860  int maximumSwapped=0;
861  // pivot elements
862  double * spare;
863  // indices
864  int * index, * indexSwapped;
865  int * saveSwapped;
866  spareArray->clear();
867  spareArray2->clear();
868  spare = spareArray->denseVector();
869  index = spareArray->getIndices();
870  saveSwapped = (int *) spareArray2->denseVector();
871  indexSwapped = spareArray2->getIndices();
872
873  // we also need somewhere for effective rhs
874  double * rhs=rhsArray->denseVector();
875
876  /*
877    First we get a list of possible pivots.  We can also see if the
878    problem looks unbounded.
879
880    At first we increase theta and see what happens.  We start
881    theta at a reasonable guess.  If in right area then we do bit by bit.
882    We save possible pivot candidates
883
884   */
885
886  // do first pass to get possibles
887  // We can also see if unbounded
888  // We also re-compute reduced cost
889
890  dualIn_ = cost_[sequenceIn_];
891
892  double * work=rowArray->denseVector();
893  int number=rowArray->getNumElements();
894  int * which=rowArray->getIndices();
895
896  // we need to swap sign if coming in from ub
897  double way = directionIn_;
898  double maximumMovement;
899  if (way>0.0) 
900    maximumMovement = min(1.0e30,upperIn_-valueIn_);
901  else
902    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
903
904  double tentativeTheta = maximumMovement;
905  double upperTheta = maximumMovement;
906
907  int iIndex;
908#if 0
909  if (numberIterations_<=39)
910    handler_->setLogLevel(63);
911  else
912    handler_->setLogLevel(2);
913  if (numberIterations_==38)
914    printf("trouble\n");
915  assert (solution_[29176]>-1.0e20);
916#endif
917  for (iIndex=0;iIndex<number;iIndex++) {
918
919    int iRow = which[iIndex];
920    double alpha = work[iRow];
921    int iPivot=pivotVariable_[iRow];
922    dualIn_ -= alpha*cost(iPivot);
923    alpha *= way;
924    double oldValue = solution(iPivot);
925    // get where in bound sequence
926    if (alpha>0.0) {
927      // basic variable going towards lower bound
928      double bound = lower(iPivot);
929      oldValue -= bound;
930    } else {
931      // basic variable going towards upper bound
932      double bound = upper(iPivot);
933      oldValue = bound-oldValue;
934    }
935     
936    double value = oldValue-tentativeTheta*fabs(alpha);
937    assert (oldValue>=-primalTolerance_*1.002);
938    if (value<-primalTolerance_) {
939      // add to list
940      spare[numberRemaining]=alpha;
941      rhs[iRow]=oldValue;
942      index[numberRemaining++]=iRow;
943      double value=oldValue-upperTheta*fabs(alpha);
944      if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
945        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
946    }
947  }
948#if 0
949  if (numberIterations_>17701)
950    handler_->setLogLevel(63);
951  if (!valuesPass&&fabs(dualIn_-saveDj)>1.0e-1*(1.0+fabs(saveDj))) {
952    double d=0.0;
953    for (iIndex=0;iIndex<number;iIndex++) {
954     
955      int iRow = which[iIndex];
956      double alpha = work[iRow];
957      int iPivot=pivotVariable_[iRow];
958      double value = alpha*cost(iPivot);
959      d -= value;
960      if (value>1.0e7)
961        printf("%d %g\n",iRow,value);
962    }
963  }
964#endif
965  // we need to keep where rhs non-zeros are
966  int numberInRhs=numberRemaining;
967  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
968  rhsArray->setNumElements(numberInRhs);
969
970  theta_=maximumMovement;
971
972  double dualCheck = fabs(dualIn_);
973  // but make a bit more pessimistic
974  dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
975
976  bool goBackOne = false;
977
978  if (numberRemaining) {
979
980    // looks like pivoting
981    // now try until reasonable theta
982    tentativeTheta = max(10.0*upperTheta,1.0e-7);
983    tentativeTheta = min(tentativeTheta,maximumMovement);
984   
985    // loops increasing tentative theta until can't go through
986   
987    while (tentativeTheta <= maximumMovement) {
988      double thruThis = 0.0;
989     
990      double bestPivot=acceptablePivot;
991      pivotRow_ = -1;
992     
993      numberSwapped = 0;
994     
995      upperTheta = maximumMovement;
996     
997      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
998
999        int iRow = index[iIndex];
1000        double alpha = spare[iIndex];
1001        double oldValue = rhs[iRow];
1002        double value = oldValue-tentativeTheta*fabs(alpha);
1003
1004        if (value<-primalTolerance_) {
1005          // how much would it cost to go thru
1006          thruThis += alpha*
1007            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1008          // goes on swapped list (also means candidates if too many)
1009          indexSwapped[numberSwapped++]=iRow;
1010          if (fabs(alpha)>bestPivot) {
1011            bestPivot=fabs(alpha);
1012            pivotRow_ = iRow;
1013            theta_ = oldValue/bestPivot;
1014          }
1015        } else {
1016          value = oldValue-upperTheta*fabs(alpha);
1017          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1018            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1019        }
1020      }
1021     
1022      maximumSwapped = max(maximumSwapped,numberSwapped);
1023
1024      if (totalThru+thruThis>=dualCheck) {
1025        // We should be pivoting in this batch
1026        // so compress down to this lot
1027
1028        int saveNumber = numberRemaining;
1029        numberRemaining=0;
1030        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1031          int iRow = indexSwapped[iIndex];
1032          spare[numberRemaining]=way*work[iRow];
1033          index[numberRemaining++]=iRow;
1034        }
1035        memset(spare+numberRemaining,0,
1036               (saveNumber-numberRemaining)*sizeof(double));
1037        int iTry;
1038#define MAXTRY 100
1039        // first get ratio with tolerance
1040        for (iTry=0;iTry<MAXTRY;iTry++) {
1041         
1042          upperTheta=maximumMovement;
1043          numberSwapped = 0;
1044          int iBest=-1;
1045          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1046           
1047            int iRow = index[iIndex];
1048            double alpha = fabs(spare[iIndex]);
1049            double oldValue = rhs[iRow];
1050            double value = oldValue-upperTheta*alpha;
1051           
1052            if (value<-primalTolerance_ && alpha>=acceptablePivot) {
1053              upperTheta = (oldValue+primalTolerance_)/alpha;
1054              iBest=iRow; // just in case weird numbers
1055            }
1056          }
1057         
1058          // now look at best in this lot
1059          bestPivot=acceptablePivot;
1060          pivotRow_=-1;
1061          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1062           
1063            int iRow = index[iIndex];
1064            double alpha = spare[iIndex];
1065            double oldValue = rhs[iRow];
1066            double value = oldValue-upperTheta*fabs(alpha);
1067           
1068            if (value<=0||iRow==iBest) {
1069              // how much would it cost to go thru
1070              totalThru += alpha*
1071                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1072              // goes on swapped list (also means candidates if too many)
1073              indexSwapped[numberSwapped++]=iRow;
1074              if (fabs(alpha)>bestPivot) {
1075                bestPivot=fabs(alpha);
1076                theta_ = oldValue/bestPivot;
1077                pivotRow_=iRow;
1078              }
1079            } else {
1080              value = oldValue-upperTheta*fabs(alpha);
1081              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1082                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1083            }
1084          }
1085         
1086          maximumSwapped = max(maximumSwapped,numberSwapped);
1087          // had (&&bestPivot<1.0e-3||totalThru>0.1*dualCheck) as well
1088          if (bestPivot<0.1*bestEverPivot&&
1089              bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
1090            // back to previous one
1091            goBackOne = true;
1092            break;
1093          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1094            if (lastPivot>acceptablePivot) {
1095              // back to previous one
1096              goBackOne = true;
1097              //break;
1098            } else {
1099              // can only get here if all pivots so far too small
1100            }
1101            break;
1102          } else if (totalThru>=dualCheck) {
1103            break; // no point trying another loop
1104          } else {
1105            // skip this lot
1106            nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1107            lastPivotRow=pivotRow_;
1108            lastTheta = theta_;
1109            lastThru = numberThru;
1110            numberThru += numberSwapped;
1111            lastNumberSwapped = numberSwapped;
1112            memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1113            if (bestPivot>bestEverPivot)
1114              bestEverPivot=bestPivot;
1115          }
1116        }
1117        break;
1118      } else {
1119        // skip this lot
1120        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1121        lastPivotRow=pivotRow_;
1122        lastTheta = theta_;
1123        lastThru = numberThru;
1124        numberThru += numberSwapped;
1125        lastNumberSwapped = numberSwapped;
1126        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1127        if (bestPivot>bestEverPivot)
1128          bestEverPivot=bestPivot;
1129        totalThru += thruThis;
1130        tentativeTheta = 2.0*upperTheta;
1131      }
1132    }
1133    // can get here without pivotRow_ set but with lastPivotRow
1134    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1135      // back to previous one
1136      pivotRow_=lastPivotRow;
1137      theta_ = lastTheta;
1138            // undo this lot
1139      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1140      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1141      numberSwapped = lastNumberSwapped;
1142    }
1143  }
1144
1145  if (pivotRow_>=0) {
1146   
1147    alpha_ = work[pivotRow_];
1148    // translate to sequence
1149    sequenceOut_ = pivotVariable_[pivotRow_];
1150    valueOut_ = solution(sequenceOut_);
1151    lowerOut_=lower_[sequenceOut_];
1152    upperOut_=upper_[sequenceOut_];
1153#define MINIMUMTHETA 1.0e-12
1154    // Movement should be minimum for anti-degeneracy - unless
1155    // fixed variable out
1156    double minimumTheta;
1157    if (upperOut_>lowerOut_)
1158      minimumTheta=MINIMUMTHETA;
1159    else
1160      minimumTheta=0.0;
1161    // will we need to increase tolerance
1162    //#define CLP_DEBUG
1163#ifdef CLP_DEBUG
1164    bool found=false;
1165#endif
1166    double largestInfeasibility = primalTolerance_;
1167    if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
1168      theta_=minimumTheta;
1169      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1170        int iRow = indexSwapped[iIndex];
1171#ifdef CLP_DEBUG
1172        if (iRow==pivotRow_)
1173          found=true;
1174#endif
1175        largestInfeasibility = max (largestInfeasibility,
1176                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
1177      }
1178#ifdef CLP_DEBUG
1179      assert(found);
1180      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
1181        printf("Primal tolerance increased from %g to %g\n",
1182               primalTolerance_,largestInfeasibility);
1183#endif
1184#undef CLP_DEBUG
1185      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1186    }
1187
1188    if (way<0.0) 
1189      theta_ = - theta_;
1190    double newValue = valueOut_ - theta_*alpha_;
1191    // If  4 bit set - Force outgoing variables to exact bound (primal)
1192    if (alpha_*way<0.0) {
1193      directionOut_=-1;      // to upper bound
1194      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1195        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1196      } else {
1197          upperOut_ = newValue;
1198      }
1199    } else {
1200      directionOut_=1;      // to lower bound
1201      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1202        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1203      } else {
1204        lowerOut_ = newValue;
1205      }
1206    }
1207    dualOut_ = reducedCost(sequenceOut_);
1208  } else if (maximumMovement<1.0e20) {
1209    // flip
1210    pivotRow_ = -2; // so we can tell its a flip
1211    sequenceOut_ = sequenceIn_;
1212    valueOut_ = valueIn_;
1213    dualOut_ = dualIn_;
1214    lowerOut_ = lowerIn_;
1215    upperOut_ = upperIn_;
1216    alpha_ = 0.0;
1217    if (way<0.0) {
1218      directionOut_=1;      // to lower bound
1219      theta_ = lowerOut_ - valueOut_;
1220    } else {
1221      directionOut_=-1;      // to upper bound
1222      theta_ = upperOut_ - valueOut_;
1223    }
1224  }
1225
1226  // clear arrays
1227
1228  memset(spare,0,numberRemaining*sizeof(double));
1229  memset(saveSwapped,0,maximumSwapped*sizeof(int));
1230
1231  // put back original bounds etc
1232  nonLinearCost_->goBackAll(rhsArray);
1233
1234  rhsArray->clear();
1235
1236}
1237/*
1238   Chooses primal pivot column
1239   updateArray has cost updates (also use pivotRow_ from last iteration)
1240   Would be faster with separate region to scan
1241   and will have this (with square of infeasibility) when steepest
1242   For easy problems we can just choose one of the first columns we look at
1243*/
1244void 
1245ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
1246                               CoinIndexedVector * spareRow1,
1247                               CoinIndexedVector * spareRow2,
1248                               CoinIndexedVector * spareColumn1,
1249                               CoinIndexedVector * spareColumn2)
1250{
1251  sequenceIn_ = primalColumnPivot_->pivotColumn(updates,spareRow1,
1252                                               spareRow2,spareColumn1,
1253                                               spareColumn2);
1254  if (sequenceIn_>=0) {
1255    valueIn_=solution_[sequenceIn_];
1256    dualIn_=dj_[sequenceIn_];
1257    if (nonLinearCost_->lookBothWays()) {
1258      // double check
1259      ClpSimplex::Status status = getStatus(sequenceIn_);
1260     
1261      switch(status) {
1262      case ClpSimplex::atUpperBound:
1263        if (dualIn_<0.0) {
1264          // move to other side
1265          printf("For %d U (%g, %g, %g) dj changed from %g",
1266                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1267                 upper_[sequenceIn_],dualIn_);
1268          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
1269          printf(" to %g\n",dualIn_);
1270          nonLinearCost_->setOne(sequenceIn_,upper_[sequenceIn_]+2.0*currentPrimalTolerance());
1271          setStatus(sequenceIn_,ClpSimplex::atLowerBound);
1272        }
1273        break;
1274      case ClpSimplex::atLowerBound:
1275        if (dualIn_>0.0) {
1276          // move to other side
1277          printf("For %d L (%g, %g, %g) dj changed from %g",
1278                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1279                 upper_[sequenceIn_],dualIn_);
1280          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
1281          printf(" to %g\n",dualIn_);
1282          nonLinearCost_->setOne(sequenceIn_,lower_[sequenceIn_]-2.0*currentPrimalTolerance());
1283          setStatus(sequenceIn_,ClpSimplex::atUpperBound);
1284        }
1285        break;
1286      default:
1287        break;
1288      }
1289    }
1290    lowerIn_=lower_[sequenceIn_];
1291    upperIn_=upper_[sequenceIn_];
1292    if (dualIn_>0.0)
1293      directionIn_ = -1;
1294    else 
1295      directionIn_ = 1;
1296  } else {
1297    sequenceIn_ = -1;
1298  }
1299}
1300/* The primals are updated by the given array.
1301   Returns number of infeasibilities.
1302   After rowArray will have list of cost changes
1303*/
1304int 
1305ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
1306                  double theta,
1307                  double & objectiveChange)
1308{
1309  double * work=rowArray->denseVector();
1310  int number=rowArray->getNumElements();
1311  int * which=rowArray->getIndices();
1312
1313  int newNumber = 0;
1314
1315  nonLinearCost_->setChangeInCost(0.0);
1316  int iIndex;
1317
1318  for (iIndex=0;iIndex<number;iIndex++) {
1319
1320    int iRow = which[iIndex];
1321    double alpha = work[iRow];
1322    int iPivot=pivotVariable_[iRow];
1323    double & value = solutionAddress(iPivot);
1324    double change = theta*alpha;
1325    value -= change;
1326    // But make sure one going out is feasible
1327    if (change>0.0) {
1328      // going down
1329      if (value<lower(iPivot)+primalTolerance_) {
1330        if (iPivot==sequenceOut_&&value>lower(iPivot)-1.001*primalTolerance_)
1331          value=lower(iPivot);
1332        double difference = nonLinearCost_->setOne(iPivot,value);
1333        work[iRow] = difference;
1334        if (difference) {
1335          //change reduced cost on this
1336          reducedCostAddress(iPivot) = -difference;
1337          which[newNumber++]=iRow;
1338        }
1339      } else {
1340        work[iRow]=0.0;
1341      }
1342    } else {
1343      // going up
1344      if (value>upper(iPivot)-primalTolerance_) {
1345        if (iPivot==sequenceOut_&&value<upper(iPivot)+1.001*primalTolerance_)
1346          value=upper(iPivot);
1347        double difference = nonLinearCost_->setOne(iPivot,value);
1348        work[iRow] = difference;
1349        if (difference) {
1350          //change reduced cost on this
1351          reducedCostAddress(iPivot) = -difference;
1352          which[newNumber++]=iRow;
1353        }
1354      } else {
1355        work[iRow]=0.0;
1356      }
1357    }
1358  }
1359  objectiveChange += nonLinearCost_->changeInCost();
1360  rowArray->setNumElements(newNumber);
1361#if 0
1362  if (newNumber>10) {
1363    printf("in %d out %d, row %d alpha %g\n",
1364           sequenceIn_,sequenceOut_,pivotRow_,alpha_);
1365    for (iIndex=0;iIndex<newNumber;iIndex++) {
1366      int iRow = which[iIndex];
1367      double alpha = work[iRow];
1368      printf("%d %g\n",iRow,alpha);
1369    }
1370  }
1371#endif
1372
1373  return 0;
1374}
1375// Perturbs problem
1376void 
1377ClpSimplexPrimal::perturb(int type)
1378{
1379  if (perturbation_>100)
1380    return; //perturbed already
1381  if (perturbation_==100)
1382    perturbation_=50; // treat as normal
1383  int i;
1384  if (!numberIterations_)
1385    cleanStatus(); // make sure status okay
1386  if (!numberIterations_&&perturbation_==50) {
1387    // See if we need to perturb
1388    double * sort = new double[numberRows_];
1389    for (i=0;i<numberRows_;i++) {
1390      double lo = fabs(lower_[i]);
1391      double up = fabs(upper_[i]);
1392      double value=0.0;
1393      if (lo&&lo<1.0e20) {
1394        if (up&&up<1.0e20)
1395          value = 0.5*(lo+up);
1396        else
1397          value=lo;
1398      } else {
1399        if (up&&up<1.0e20)
1400          value = up;
1401      }
1402      sort[i]=value;
1403    }
1404    std::sort(sort,sort+numberRows_);
1405    int number=1;
1406    double last = sort[0];
1407    for (i=1;i<numberRows_;i++) {
1408      if (last!=sort[i])
1409        number++;
1410      last=sort[i];
1411    }
1412    delete [] sort;
1413    if (number*4>numberRows_)
1414      return; // good enough
1415  }
1416  // primal perturbation
1417  double perturbation=1.0e-20;
1418  int numberNonZero=0;
1419  // maximum fraction of rhs/bounds to perturb
1420  double maximumFraction = 1.0e-7;
1421  if (perturbation_>=50) {
1422    perturbation = 1.0e-4;
1423    for (i=0;i<numberColumns_+numberRows_;i++) {
1424      if (upper_[i]>lower_[i]+primalTolerance_) {
1425        double lowerValue, upperValue;
1426        if (lower_[i]>-1.0e20)
1427          lowerValue = fabs(lower_[i]);
1428        else
1429          lowerValue=0.0;
1430        if (upper_[i]<1.0e20)
1431          upperValue = fabs(upper_[i]);
1432        else
1433          upperValue=0.0;
1434        double value = max(fabs(lowerValue),fabs(upperValue));
1435        value = min(value,upper_[i]-lower_[i]);
1436#if 1
1437        if (value) {
1438          perturbation += value;
1439          numberNonZero++;
1440        }
1441#else
1442        perturbation = max(perturbation,value);
1443#endif
1444      }
1445    }
1446    if (numberNonZero) 
1447      perturbation /= (double) numberNonZero;
1448    else
1449      perturbation = 1.0e-1;
1450  } else if (perturbation_<100) {
1451    perturbation = pow(10.0,perturbation_);
1452    // user is in charge
1453    maximumFraction = 1.0;
1454  }
1455  double largestZero=0.0;
1456  double largest=0.0;
1457  double largestPerCent=0.0;
1458  bool printOut=(handler_->logLevel()==63);
1459  printOut=false; //off
1460  // Check if all slack
1461  int number=0;
1462  int iSequence;
1463  for (iSequence=0;iSequence<numberRows_;iSequence++) {
1464    if (getRowStatus(iSequence)==basic) 
1465      number++;
1466  }
1467  if (number!=numberRows_)
1468    type=1;
1469  // modify bounds
1470  if (type==1) {
1471    double multiplier = perturbation*maximumFraction;
1472    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1473      if (getStatus(iSequence)==basic) {
1474        double solutionValue = solution_[iSequence];
1475        double lowerValue = lower_[iSequence];
1476        double upperValue = upper_[iSequence];
1477        double difference = upperValue-lowerValue;
1478        difference = min(difference,perturbation);
1479        difference = min(difference,fabs(solutionValue)+1.0);
1480        double value = CoinDrand48()*multiplier*(difference+1.0);
1481        if (solutionValue-lowerValue<=primalTolerance_) {
1482          lower_[iSequence] -= value;
1483        } else if (upperValue-solutionValue<=primalTolerance_) {
1484          upper_[iSequence] += value;
1485        } else {
1486#if 0
1487          if (iSequence>=numberColumns_) {
1488            // may not be at bound - but still perturb (unless free)
1489            if (upperValue>1.0e30&&lowerValue<-1.0e30)
1490              value=0.0;
1491            else
1492              value = - value; // as -1.0 in matrix
1493          } else {
1494            value = 0.0;
1495          }
1496#else
1497          value=0.0;
1498#endif
1499        }
1500        if (value) {
1501          if (printOut)
1502            printf("col %d lower from %g to %g, upper from %g to %g\n",
1503                   iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
1504          if (solutionValue) {
1505            largest = max(largest,value);
1506            if (value>(fabs(solutionValue)+1.0)*largestPerCent)
1507              largestPerCent=value/(fabs(solutionValue)+1.0);
1508          } else {
1509            largestZero = max(largestZero,value);
1510          } 
1511        }
1512      }
1513    }
1514  } else {
1515    for (i=0;i<numberColumns_;i++) {
1516      double lowerValue=lower_[i], upperValue=upper_[i];
1517      if (upperValue>lowerValue+primalTolerance_) {
1518        double value = CoinDrand48()*perturbation*maximumFraction;
1519        if (lowerValue>-1.0e20&&lowerValue)
1520          lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue))); 
1521        if (upperValue<1.0e20&&upperValue)
1522          upperValue += value * (max(1.0,1.0e-5*fabs(upperValue))); 
1523        if (lowerValue!=lower_[i]) {
1524          double difference = fabs(lowerValue-lower_[i]);
1525          largest = max(largest,difference);
1526          if (difference>fabs(lower_[i])*largestPerCent)
1527            largestPerCent=fabs(difference/lower_[i]);
1528        } 
1529        if (upperValue!=upper_[i]) {
1530          double difference = fabs(upperValue-upper_[i]);
1531          largest = max(largest,difference);
1532          if (difference>fabs(upper_[i])*largestPerCent)
1533            largestPerCent=fabs(difference/upper_[i]);
1534        } 
1535        if (printOut)
1536          printf("col %d lower from %g to %g, upper from %g to %g\n",
1537                 i,lower_[i],lowerValue,upper_[i],upperValue);
1538      }
1539      lower_[i]=lowerValue;
1540      upper_[i]=upperValue;
1541    }
1542    for (;i<numberColumns_+numberRows_;i++) {
1543      double lowerValue=lower_[i], upperValue=upper_[i];
1544      double value = CoinDrand48()*perturbation*maximumFraction;
1545      if (upperValue>lowerValue+primalTolerance_) {
1546        if (lowerValue>-1.0e20&&lowerValue)
1547          lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue))); 
1548        if (upperValue<1.0e20&&upperValue)
1549          upperValue += value * (max(1.0,1.0e-5*fabs(upperValue))); 
1550      } else if (upperValue>0.0) {
1551        upperValue -= value * (max(1.0,1.0e-5*fabs(lowerValue))); 
1552        lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue))); 
1553      } else if (upperValue<0.0) {
1554        upperValue += value * (max(1.0,1.0e-5*fabs(lowerValue))); 
1555        lowerValue += value * (max(1.0,1.0e-5*fabs(lowerValue))); 
1556      } else {
1557      }
1558      if (lowerValue!=lower_[i]) {
1559        double difference = fabs(lowerValue-lower_[i]);
1560        largest = max(largest,difference);
1561        if (difference>fabs(lower_[i])*largestPerCent)
1562          largestPerCent=fabs(difference/lower_[i]);
1563      } 
1564      if (upperValue!=upper_[i]) {
1565        double difference = fabs(upperValue-upper_[i]);
1566        largest = max(largest,difference);
1567        if (difference>fabs(upper_[i])*largestPerCent)
1568          largestPerCent=fabs(difference/upper_[i]);
1569      } 
1570      if (printOut)
1571        printf("row %d lower from %g to %g, upper from %g to %g\n",
1572               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
1573      lower_[i]=lowerValue;
1574      upper_[i]=upperValue;
1575    }
1576  }
1577  // Clean up
1578  for (i=0;i<numberColumns_+numberRows_;i++) {
1579    switch(getStatus(i)) {
1580     
1581    case basic:
1582      break;
1583    case atUpperBound:
1584      solution_[i]=upper_[i];
1585      break;
1586    case isFixed:
1587    case atLowerBound:
1588      solution_[i]=lower_[i];
1589      break;
1590    case isFree:
1591      break;
1592    case superBasic:
1593      break;
1594    }
1595  }
1596  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
1597    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
1598    <<CoinMessageEol;
1599  // redo nonlinear costs
1600  // say perturbed
1601  perturbation_=101;
1602}
1603// un perturb
1604bool
1605ClpSimplexPrimal::unPerturb()
1606{
1607  if (perturbation_!=101)
1608    return false;
1609  // put back original bounds and costs
1610  createRim(7);
1611  // unflag
1612  unflag();
1613  // get a valid nonlinear cost function
1614  delete nonLinearCost_;
1615  nonLinearCost_= new ClpNonLinearCost(this);
1616  perturbation_ = 102; // stop any further perturbation
1617  // move non basic variables to new bounds
1618  nonLinearCost_->checkInfeasibilities(0.0);
1619#if 1
1620  // Try using dual
1621  return true;
1622#else
1623  gutsOfSolution(NULL,NULL);
1624  return false;
1625#endif
1626 
1627}
1628// Unflag all variables and return number unflagged
1629int 
1630ClpSimplexPrimal::unflag()
1631{
1632  int i;
1633  int number = numberRows_+numberColumns_;
1634  int numberFlagged=0;
1635  for (i=0;i<number;i++) {
1636    if (flagged(i)) {
1637      clearFlagged(i);
1638      numberFlagged++;
1639    }
1640  }
1641  return numberFlagged;
1642}
1643// Do not change infeasibility cost and always say optimal
1644void 
1645ClpSimplexPrimal::alwaysOptimal(bool onOff)
1646{
1647  if (onOff)
1648    specialOptions_ |= 1;
1649  else
1650    specialOptions_ &= ~1;
1651}
1652bool 
1653ClpSimplexPrimal::alwaysOptimal() const
1654{
1655  return (specialOptions_&1)!=0;
1656}
1657// Flatten outgoing variables i.e. - always to exact bound
1658void 
1659ClpSimplexPrimal::exactOutgoing(bool onOff)
1660{
1661  if (onOff)
1662    specialOptions_ |= 4;
1663  else
1664    specialOptions_ &= ~4;
1665}
1666bool 
1667ClpSimplexPrimal::exactOutgoing() const
1668{
1669  return (specialOptions_&4)!=0;
1670}
1671/*
1672  Reasons to come out (normal mode/user mode):
1673  -1 normal
1674  -2 factorize now - good iteration/ NA
1675  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
1676  -4 inaccuracy - refactorize - no iteration/ NA
1677  -5 something flagged - go round again/ pivot not possible
1678  +2 looks unbounded
1679  +3 max iterations (iteration done)
1680*/
1681int
1682ClpSimplexPrimal::pivotResult(int ifValuesPass)
1683{
1684
1685  bool roundAgain=true;
1686  int returnCode=-1;
1687
1688  // loop round if user setting and doing refactorization
1689  while (roundAgain) {
1690    roundAgain=false;
1691    returnCode=-1;
1692    pivotRow_=-1;
1693    sequenceOut_=-1;
1694    rowArray_[1]->clear();
1695    // we found a pivot column
1696    // update the incoming column
1697    unpack(rowArray_[1]);
1698    // save reduced cost
1699    double saveDj = dualIn_;
1700    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1701    // do ratio test and re-compute dj
1702    primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
1703              ifValuesPass);
1704    if (ifValuesPass) {
1705      saveDj=dualIn_;
1706      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
1707        if(fabs(dualIn_)<1.0e2*dualTolerance_) {
1708          // try other way
1709          directionIn_=-directionIn_;
1710          primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
1711                    0);
1712        }
1713        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
1714          if (solveType_==1) {
1715            // reject it
1716            char x = isColumn(sequenceIn_) ? 'C' :'R';
1717            handler_->message(CLP_SIMPLEX_FLAG,messages_)
1718              <<x<<sequenceWithin(sequenceIn_)
1719              <<CoinMessageEol;
1720            setFlagged(sequenceIn_);
1721            lastBadIteration_ = numberIterations_; // say be more cautious
1722            rowArray_[1]->clear();
1723            pivotRow_=-1;
1724          }
1725          returnCode=-5;
1726          break;
1727        }
1728      }
1729    }
1730    double checkValue=1.0e-2;
1731    if (largestDualError_>1.0e-5)
1732      checkValue=1.0e-1;
1733    if (solveType_==1&&((saveDj*dualIn_<1.0e-20&&!ifValuesPass)||
1734        fabs(saveDj-dualIn_)>checkValue*(1.0+fabs(saveDj)))) {
1735      handler_->message(CLP_PRIMAL_DJ,messages_)
1736        <<saveDj<<dualIn_
1737        <<CoinMessageEol;
1738      if(lastGoodIteration_ != numberIterations_) {
1739        rowArray_[1]->clear();
1740        pivotRow_=-1; // say no weights update
1741        returnCode=-4;
1742        if(lastGoodIteration_+1 == numberIterations_) {
1743          // not looking wonderful - try cleaning bounds
1744          // put non-basics to bounds in case tolerance moved
1745          nonLinearCost_->checkInfeasibilities(0.0);
1746        }
1747        sequenceOut_=-1;
1748        break;
1749      } else {
1750        // take on more relaxed criterion
1751        if (saveDj*dualIn_<1.0e-20||
1752            fabs(saveDj-dualIn_)>2.0e-1*(1.0+fabs(dualIn_))) {
1753          // need to reject something
1754          char x = isColumn(sequenceIn_) ? 'C' :'R';
1755          handler_->message(CLP_SIMPLEX_FLAG,messages_)
1756            <<x<<sequenceWithin(sequenceIn_)
1757            <<CoinMessageEol;
1758          setFlagged(sequenceIn_);
1759          lastBadIteration_ = numberIterations_; // say be more cautious
1760          rowArray_[1]->clear();
1761          pivotRow_=-1;
1762          returnCode=-5;
1763          sequenceOut_=-1;
1764          break;
1765        }
1766      }
1767    }
1768    if (pivotRow_>=0) {
1769      if (solveType_==2) {
1770        // **** Coding for user interface
1771        // do ray
1772        primalRay(rowArray_[1]);
1773        // update duals
1774        if (pivotRow_>=0) {
1775          alpha_ = rowArray_[1]->denseVector()[pivotRow_];
1776          assert (fabs(alpha_)>1.0e-8);
1777          double multiplier = dualIn_/alpha_;
1778          rowArray_[0]->insert(pivotRow_,multiplier);
1779          factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
1780          // put row of tableau in rowArray[0] and columnArray[0]
1781          matrix_->transposeTimes(this,-1.0,
1782                                  rowArray_[0],columnArray_[1],columnArray_[0]);
1783          // update column djs
1784          int i;
1785          int * index = columnArray_[0]->getIndices();
1786          int number = columnArray_[0]->getNumElements();
1787          double * element = columnArray_[0]->denseVector();
1788          for (i=0;i<number;i++) {
1789            int ii = index[i];
1790            dj_[ii] += element[ii];
1791            element[ii]=0.0;
1792          }
1793          columnArray_[0]->setNumElements(0);
1794          // and row djs
1795          index = rowArray_[0]->getIndices();
1796          number = rowArray_[0]->getNumElements();
1797          element = rowArray_[0]->denseVector();
1798          for (i=0;i<number;i++) {
1799            int ii = index[i];
1800            dj_[ii+numberColumns_] += element[ii];
1801            dual_[ii] = dj_[ii+numberColumns_];
1802            element[ii]=0.0;
1803          }
1804          rowArray_[0]->setNumElements(0);
1805          // check incoming
1806          assert (fabs(dj_[sequenceIn_])<1.0e-6);
1807        }
1808      }
1809      // if stable replace in basis
1810      int updateStatus = factorization_->replaceColumn(rowArray_[2],
1811                                                       pivotRow_,
1812                                                       alpha_);
1813      // if no pivots, bad update but reasonable alpha - take and invert
1814      if (updateStatus==2&&
1815          lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
1816        updateStatus=4;
1817      if (updateStatus==1||updateStatus==4) {
1818        // slight error
1819        if (factorization_->pivots()>5||updateStatus==4) {
1820          returnCode=-3;
1821        }
1822      } else if (updateStatus==2) {
1823        // major error
1824        // better to have small tolerance even if slower
1825        factorization_->zeroTolerance(1.0e-15);
1826        int maxFactor = factorization_->maximumPivots();
1827        if (maxFactor>10) {
1828          if (forceFactorization_<0)
1829            forceFactorization_= maxFactor;
1830          forceFactorization_ = max (1,(forceFactorization_>>1));
1831        } 
1832        // later we may need to unwind more e.g. fake bounds
1833        if(lastGoodIteration_ != numberIterations_) {
1834          rowArray_[1]->clear();
1835          pivotRow_=-1;
1836          if (solveType_==1) {
1837            returnCode=-4;
1838            break;
1839          } else {
1840            // user in charge - re-factorize
1841            int lastCleaned;
1842            ClpSimplexProgress dummyProgress;
1843            if (saveStatus_)
1844              statusOfProblemInPrimal(lastCleaned,1,&dummyProgress);
1845            else
1846              statusOfProblemInPrimal(lastCleaned,0,&dummyProgress);
1847            roundAgain=true;
1848            continue;
1849          }
1850        } else {
1851          // need to reject something
1852          if (solveType_==1) {
1853            char x = isColumn(sequenceIn_) ? 'C' :'R';
1854            handler_->message(CLP_SIMPLEX_FLAG,messages_)
1855              <<x<<sequenceWithin(sequenceIn_)
1856              <<CoinMessageEol;
1857            setFlagged(sequenceIn_);
1858          }
1859          lastBadIteration_ = numberIterations_; // say be more cautious
1860          rowArray_[1]->clear();
1861          pivotRow_=-1;
1862          sequenceOut_=-1;
1863          returnCode = -5;
1864          break;
1865
1866        }
1867      } else if (updateStatus==3) {
1868        // out of memory
1869        // increase space if not many iterations
1870        if (factorization_->pivots()<
1871            0.5*factorization_->maximumPivots()&&
1872            factorization_->pivots()<200)
1873          factorization_->areaFactor(
1874                                     factorization_->areaFactor() * 1.1);
1875        returnCode =-2; // factorize now
1876      } else if (updateStatus==5) {
1877        problemStatus_=-2; // factorize now
1878      }
1879      // here do part of steepest - ready for next iteration
1880      if (!ifValuesPass)
1881        primalColumnPivot_->updateWeights(rowArray_[1]);
1882    } else {
1883      if (pivotRow_==-1) {
1884        // no outgoing row is valid
1885        rowArray_[0]->clear();
1886        if (!factorization_->pivots()) {
1887          returnCode = 2; //say looks unbounded
1888          // do ray
1889          primalRay(rowArray_[1]);
1890        } else if (solveType_==2) {
1891          // refactorize
1892          int lastCleaned;
1893          ClpSimplexProgress dummyProgress;
1894          if (saveStatus_)
1895            statusOfProblemInPrimal(lastCleaned,1,&dummyProgress);
1896          else
1897            statusOfProblemInPrimal(lastCleaned,0,&dummyProgress);
1898          roundAgain=true;
1899          continue;
1900        } else {
1901          returnCode = 4; //say looks unbounded but has iterated
1902        }
1903        break;
1904      } else {
1905        // flipping from bound to bound
1906      }
1907    }
1908
1909   
1910    // update primal solution
1911   
1912    double objectiveChange=0.0;
1913    // Cost on pivot row may change - may need to change dualIn
1914    double oldCost=0.0;
1915    if (pivotRow_>=0)
1916      oldCost = cost(pivotVariable_[pivotRow_]);
1917    // rowArray_[1] is not empty - used to update djs
1918    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
1919    if (pivotRow_>=0)
1920      dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
1921    double oldValue = valueIn_;
1922    if (directionIn_==-1) {
1923      // as if from upper bound
1924      if (sequenceIn_!=sequenceOut_) {
1925        // variable becoming basic
1926        valueIn_ -= fabs(theta_);
1927      } else {
1928        valueIn_=lowerIn_;
1929      }
1930    } else {
1931      // as if from lower bound
1932      if (sequenceIn_!=sequenceOut_) {
1933        // variable becoming basic
1934        valueIn_ += fabs(theta_);
1935      } else {
1936        valueIn_=upperIn_;
1937      }
1938    }
1939    objectiveChange += dualIn_*(valueIn_-oldValue);
1940    // outgoing
1941    if (sequenceIn_!=sequenceOut_) {
1942      if (directionOut_>0) {
1943        valueOut_ = lowerOut_;
1944      } else {
1945        valueOut_ = upperOut_;
1946      }
1947      assert(valueOut_>=lower_[sequenceOut_]-primalTolerance_&&
1948             valueOut_<=upper_[sequenceOut_]+primalTolerance_);
1949      // may not be exactly at bound and bounds may have changed
1950      // Make sure outgoing looks feasible
1951      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
1952      solution_[sequenceOut_]=valueOut_;
1953    }
1954    // change cost and bounds on incoming if primal
1955    nonLinearCost_->setOne(sequenceIn_,valueIn_); 
1956    int whatNext=housekeeping(objectiveChange);
1957#if 0
1958    if (numberIterations_==1148)
1959      whatNext=1;
1960    if (numberIterations_>1200)
1961    exit(0);
1962#endif
1963    if (whatNext==1) {
1964        returnCode =-2; // refactorize
1965    } else if (whatNext==2) {
1966      // maximum iterations or equivalent
1967      returnCode=3;
1968    } else if(numberIterations_ == lastGoodIteration_
1969              + 2 * factorization_->maximumPivots()) {
1970      // done a lot of flips - be safe
1971      returnCode =-2; // refactorize
1972    }
1973  }
1974  if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
1975    // refactorize here
1976    int lastCleaned;
1977    ClpSimplexProgress dummyProgress;
1978    if (saveStatus_)
1979      statusOfProblemInPrimal(lastCleaned,1,&dummyProgress);
1980    else
1981      statusOfProblemInPrimal(lastCleaned,0,&dummyProgress);
1982    if (problemStatus_==5) {
1983      printf("Singular basis\n");
1984      problemStatus_=-1;
1985      returnCode=5;
1986    }
1987  }
1988#ifdef CLP_DEBUG
1989  {
1990    int i;
1991    // not [1] as may have information
1992    for (i=0;i<4;i++) {
1993      if (i!=1)
1994        rowArray_[i]->checkClear();
1995    }   
1996    for (i=0;i<2;i++) {
1997      columnArray_[i]->checkClear();
1998    }   
1999  }     
2000#endif
2001  return returnCode;
2002}
2003// Create primal ray
2004void 
2005ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
2006{
2007  delete [] ray_;
2008  ray_ = new double [numberColumns_];
2009  ClpFillN(ray_,numberColumns_,0.0);
2010  int number=rowArray->getNumElements();
2011  int * index = rowArray->getIndices();
2012  double * array = rowArray->denseVector();
2013  double way=-directionIn_;
2014  int i;
2015  double zeroTolerance=1.0e-12;
2016  if (sequenceIn_<numberColumns_)
2017    ray_[sequenceIn_]=directionIn_;
2018  for (i=0;i<number;i++) {
2019    int iRow=index[i];
2020    int iPivot=pivotVariable_[iRow];
2021    double arrayValue = array[iRow];
2022    if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2023      ray_[iPivot] = way* array[iRow];
2024  }
2025}
2026/* Get next superbasic -1 if none,
2027   Normal type is 1
2028   If type is 3 then initializes sorted list
2029   if 2 uses list.
2030*/
2031int 
2032ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
2033{
2034  if (firstFree_>=0&&superBasicType) {
2035    int returnValue=firstFree_;
2036    int iColumn=firstFree_+1;
2037    if (superBasicType>1) {
2038      if (superBasicType>2) {
2039        // Initialize list
2040        // Wild guess that lower bound more natural than upper
2041        int number=0;
2042        double * work=columnArray->denseVector();
2043        int * which=columnArray->getIndices();
2044        for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2045          if (getStatus(iColumn)==superBasic) {
2046            if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2047              solution_[iColumn]=lower_[iColumn];
2048              setStatus(iColumn,atLowerBound);
2049            } else if (fabs(solution_[iColumn]-upper_[iColumn])
2050                       <=primalTolerance_) {
2051              solution_[iColumn]=upper_[iColumn];
2052              setStatus(iColumn,atUpperBound);
2053            } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2054              setStatus(iColumn,isFree);
2055              break;
2056            } else if (!flagged(iColumn)) {
2057              // put ones near bounds at end after sorting
2058              work[number]= - min(0.1*(solution_[iColumn]-lower_[iColumn]),
2059                                  upper_[iColumn]-solution_[iColumn]);
2060              which[number++] = iColumn;
2061            }
2062          }
2063        }
2064        CoinSort_2(work,work+number,which);
2065        columnArray->setNumElements(number);
2066        memset(work,0,number*sizeof(double));
2067      }
2068      int * which=columnArray->getIndices();
2069      int number = columnArray->getNumElements();
2070      if (!number) {
2071        // finished
2072        iColumn = numberRows_+numberColumns_;
2073        returnValue=-1;
2074      } else {
2075        number--;
2076        returnValue=which[number];
2077        iColumn=returnValue;
2078        columnArray->setNumElements(number);
2079      }     
2080    } else {
2081      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
2082        if (getStatus(iColumn)==superBasic) {
2083          if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2084            solution_[iColumn]=lower_[iColumn];
2085            setStatus(iColumn,atLowerBound);
2086          } else if (fabs(solution_[iColumn]-upper_[iColumn])
2087                     <=primalTolerance_) {
2088            solution_[iColumn]=upper_[iColumn];
2089            setStatus(iColumn,atUpperBound);
2090          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2091            setStatus(iColumn,isFree);
2092            break;
2093          } else {
2094            break;
2095          }
2096        }
2097      }
2098    }
2099    firstFree_ = iColumn;
2100    if (firstFree_==numberRows_+numberColumns_)
2101      firstFree_=-1;
2102    return returnValue;
2103  } else {
2104    return -1;
2105  }
2106}
2107 
2108 
2109
Note: See TracBrowser for help on using the repository browser.