source: branches/pre/ClpSimplexPrimal.cpp @ 210

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

Trying

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