source: trunk/ClpSimplexPrimal.cpp @ 225

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

This should break everything

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