source: branches/pre/ClpPrimalColumnSteepest.cpp @ 221

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

Still trying to go faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 96.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include "ClpSimplex.hpp"
7#include "ClpPrimalColumnSteepest.hpp"
8#include "CoinIndexedVector.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpNonLinearCost.hpp"
11#include "ClpMessage.hpp"
12#include "CoinHelperFunctions.hpp"
13#include <stdio.h>
14
15
16// bias for free variables
17#define FREE_BIAS 1.0e1
18// Acceptance criteria for free variables
19#define FREE_ACCEPT 1.0e2
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode) 
28  : ClpPrimalColumnPivot(),
29    devex_(0.0),
30    weights_(NULL),
31    infeasible_(NULL),
32    alternateWeights_(NULL),
33    savedWeights_(NULL),
34    reference_(NULL),
35    state_(-1),
36    mode_(mode),
37    numberSwitched_(0),
38    pivotSequence_(-1),
39    savedPivotSequence_(-1),
40    savedSequenceOut_(-1)
41{
42  type_=2+64*mode;
43}
44
45//-------------------------------------------------------------------
46// Copy constructor
47//-------------------------------------------------------------------
48ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs) 
49: ClpPrimalColumnPivot(rhs)
50{ 
51  state_=rhs.state_;
52  mode_ = rhs.mode_;
53  numberSwitched_ = rhs.numberSwitched_;
54  model_ = rhs.model_;
55  pivotSequence_ = rhs.pivotSequence_;
56  savedPivotSequence_ = rhs.savedPivotSequence_;
57  savedSequenceOut_ = rhs.savedSequenceOut_;
58  devex_ = rhs.devex_;
59  if (rhs.infeasible_) {
60    infeasible_= new CoinIndexedVector(rhs.infeasible_);
61  } else {
62    infeasible_=NULL;
63  }
64  reference_=NULL;
65  if (rhs.weights_) {
66    assert(model_);
67    int number = model_->numberRows()+model_->numberColumns();
68    weights_= new double[number];
69    ClpDisjointCopyN(rhs.weights_,number,weights_);
70    savedWeights_= new double[number];
71    ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
72    if (mode_!=1) {
73      reference_ = new unsigned int[(number+31)>>5];
74      memcpy(reference_,rhs.reference_,((number+31)>>5)*sizeof(unsigned int));
75    }
76  } else {
77    weights_=NULL;
78    savedWeights_=NULL;
79  }
80  if (rhs.alternateWeights_) {
81    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
82  } else {
83    alternateWeights_=NULL;
84  }
85}
86
87//-------------------------------------------------------------------
88// Destructor
89//-------------------------------------------------------------------
90ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
91{
92  delete [] weights_;
93  delete infeasible_;
94  delete alternateWeights_;
95  delete [] savedWeights_;
96  delete [] reference_;
97}
98
99//----------------------------------------------------------------
100// Assignment operator
101//-------------------------------------------------------------------
102ClpPrimalColumnSteepest &
103ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
104{
105  if (this != &rhs) {
106    ClpPrimalColumnPivot::operator=(rhs);
107    state_=rhs.state_;
108    mode_ = rhs.mode_;
109    numberSwitched_ = rhs.numberSwitched_;
110    model_ = rhs.model_;
111    pivotSequence_ = rhs.pivotSequence_;
112    savedPivotSequence_ = rhs.savedPivotSequence_;
113    savedSequenceOut_ = rhs.savedSequenceOut_;
114    devex_ = rhs.devex_;
115    delete [] weights_;
116    delete [] reference_;
117    reference_=NULL;
118    delete infeasible_;
119    delete alternateWeights_;
120    delete [] savedWeights_;
121    savedWeights_ = NULL;
122    if (rhs.infeasible_!=NULL) {
123      infeasible_= new CoinIndexedVector(rhs.infeasible_);
124    } else {
125      infeasible_=NULL;
126    }
127    if (rhs.weights_!=NULL) {
128      assert(model_);
129      int number = model_->numberRows()+model_->numberColumns();
130      weights_= new double[number];
131      ClpDisjointCopyN(rhs.weights_,number,weights_);
132      savedWeights_= new double[number];
133      ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
134      if (mode_!=1) {
135        reference_ = new unsigned int[(number+31)>>5];
136        memcpy(reference_,rhs.reference_,
137               ((number+31)>>5)*sizeof(unsigned int));
138      }
139    } else {
140      weights_=NULL;
141    }
142    if (rhs.alternateWeights_!=NULL) {
143      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
144    } else {
145      alternateWeights_=NULL;
146    }
147  }
148  return *this;
149}
150
151#define TRY_NORM 1.0e-4
152#define ADD_ONE 1.0
153// Returns pivot column, -1 if none
154int 
155ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
156                                    CoinIndexedVector * spareRow1,
157                                    CoinIndexedVector * spareRow2,
158                                    CoinIndexedVector * spareColumn1,
159                                    CoinIndexedVector * spareColumn2)
160{
161  assert(model_);
162  if (model_->nonLinearCost()->lookBothWays()||model_->algorithm()==2) {
163    // Do old way
164    return pivotColumnOldMethod(updates,spareRow1,spareRow2,
165                                spareColumn1,spareColumn2);
166  }
167  int number=0;
168  int * index;
169  // dj could be very small (or even zero - take care)
170  double dj = model_->dualIn();
171  double tolerance=model_->currentDualTolerance();
172  // we can't really trust infeasibilities if there is dual error
173  // this coding has to mimic coding in checkDualSolution
174  double error = min(1.0e-3,model_->largestDualError());
175  // allow tolerance at least slightly bigger than standard
176  tolerance = tolerance +  error;
177  int pivotRow = model_->pivotRow();
178  int anyUpdates;
179  double * infeas = infeasible_->denseVector();
180
181  // Local copy of mode so can decide what to do
182  int switchType;
183  if (mode_==4) 
184    switchType = 5-numberSwitched_;
185  else
186    switchType = mode_;
187  /* switchType -
188     0 - all exact devex
189     1 - all steepest
190     2 - some exact devex
191     3 - auto some exact devex
192     4 - devex
193     5 - dantzig
194  */
195  if (updates->getNumElements()) {
196    // would have to have two goes for devex, three for steepest
197    anyUpdates=2;
198    // add in pivot contribution
199    if (pivotRow>=0) 
200      updates->add(pivotRow,-dj);
201  } else if (pivotRow>=0) {
202    if (fabs(dj)>1.0e-15) {
203      // some dj
204      updates->insert(pivotRow,-dj);
205      if (fabs(dj)>1.0e-6) {
206        // reasonable size
207        anyUpdates=1;
208      } else {
209        // too small
210        anyUpdates=2;
211      }
212    } else {
213      // zero dj
214      anyUpdates=-1;
215    }
216  } else if (pivotSequence_>=0){
217    // just after re-factorization
218    anyUpdates=-1;
219  } else {
220    // sub flip - nothing to do
221    anyUpdates=0;
222  }
223  int sequenceOut = model_->sequenceOut();
224  if (switchType==5) {
225    if (anyUpdates>0) {
226      justDjs(updates,spareRow1,spareRow2,
227        spareColumn1,spareColumn2);
228    }
229  } else if (anyUpdates==1) {
230    if (switchType<4) {
231      // exact etc when can use dj
232      djsAndSteepest(updates,spareRow1,spareRow2,
233        spareColumn1,spareColumn2);
234    } else {
235      // devex etc when can use dj
236      djsAndDevex(updates,spareRow1,spareRow2,
237        spareColumn1,spareColumn2);
238    }
239  } else if (anyUpdates==-1) {
240    if (switchType<4) {
241      // exact etc when djs okay
242      justSteepest(updates,spareRow1,spareRow2,
243        spareColumn1,spareColumn2);
244    } else {
245      // devex etc when djs okay
246      justDevex(updates,spareRow1,spareRow2,
247        spareColumn1,spareColumn2);
248    }
249  } else if (anyUpdates==2) {
250    if (switchType<4) {
251      // exact etc when have to use pivot
252      djsAndSteepest2(updates,spareRow1,spareRow2,
253        spareColumn1,spareColumn2);
254    } else {
255      // devex etc when have to use pivot
256      djsAndDevex2(updates,spareRow1,spareRow2,
257        spareColumn1,spareColumn2);
258    }
259  } 
260  // make sure outgoing from last iteration okay
261  if (sequenceOut>=0) {
262    ClpSimplex::Status status = model_->getStatus(sequenceOut);
263    double value = model_->reducedCost(sequenceOut);
264   
265    switch(status) {
266     
267    case ClpSimplex::basic:
268    case ClpSimplex::isFixed:
269      break;
270    case ClpSimplex::isFree:
271    case ClpSimplex::superBasic:
272      if (fabs(value)>FREE_ACCEPT*tolerance) { 
273        // we are going to bias towards free (but only if reasonable)
274        value *= FREE_BIAS;
275        // store square in list
276        if (infeas[sequenceOut])
277          infeas[sequenceOut] = value*value; // already there
278        else
279          infeasible_->quickAdd(sequenceOut,value*value);
280      } else {
281        infeasible_->zero(sequenceOut);
282      }
283      break;
284    case ClpSimplex::atUpperBound:
285      if (value>tolerance) {
286        // store square in list
287        if (infeas[sequenceOut])
288          infeas[sequenceOut] = value*value; // already there
289        else
290          infeasible_->quickAdd(sequenceOut,value*value);
291      } else {
292        infeasible_->zero(sequenceOut);
293      }
294      break;
295    case ClpSimplex::atLowerBound:
296      if (value<-tolerance) {
297        // store square in list
298        if (infeas[sequenceOut])
299          infeas[sequenceOut] = value*value; // already there
300        else
301          infeasible_->quickAdd(sequenceOut,value*value);
302      } else {
303        infeasible_->zero(sequenceOut);
304      }
305    }
306  }
307  // update of duals finished - now do pricing
308  // See what sort of pricing
309  int numberWanted=0;
310  number = infeasible_->getNumElements();
311  int numberColumns = model_->numberColumns();
312  if (switchType==5) {
313    pivotSequence_=-1;
314    pivotRow=-1;
315    // See if to switch
316    int numberElements = model_->factorization()->numberElements();
317    int numberRows = model_->numberRows();
318    // ratio is done on number of columns here
319    //double ratio = (double) numberElements/(double) numberColumns;
320    double ratio = (double) numberElements/(double) numberRows;
321    //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
322    if (ratio<0.1) {
323      numberWanted = max(100,number/200);
324    } else if (ratio<0.15) {
325      numberWanted = max(500,number/40);
326    } else if (ratio<0.2) {
327      numberWanted = max(2000,number/10);
328      numberWanted = max(numberWanted,numberColumns/30);
329    } else {
330      switchType=4;
331      // initialize
332      numberSwitched_++;
333      // Make sure will re-do
334      delete [] weights_;
335      weights_=NULL;
336      saveWeights(model_,4);
337      printf("switching to devex %d nel ratio %g\n",numberElements,ratio);
338    }
339    //if (model_->numberIterations()%1000==0)
340    //printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
341  }
342  if(switchType==4) {
343    // Still in devex mode
344    int numberElements = model_->factorization()->numberElements();
345    int numberRows = model_->numberRows();
346    // ratio is done on number of rows here
347    double ratio = (double) numberElements/(double) numberRows;
348    // Go to steepest if lot of iterations?
349    if (ratio<1.0) {
350      numberWanted = max(2000,number/20);
351    } else if (ratio<2.0) {
352      numberWanted = max(2000,number/10);
353      numberWanted = max(numberWanted,numberColumns/20);
354    } else {
355      // we can zero out
356      updates->clear();
357      spareColumn1->clear();
358      switchType=3;
359      // initialize
360      pivotSequence_=-1;
361      pivotRow=-1;
362      numberSwitched_++;
363      // Make sure will re-do
364      delete [] weights_;
365      weights_=NULL;
366      saveWeights(model_,4);
367      printf("switching to exact %d nel ratio %g\n",numberElements,ratio);
368      updates->clear();
369    }
370    //if (model_->numberIterations()%1000==0)
371    //printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
372  } 
373  if (switchType<4) {
374    if (switchType<2 ) {
375      numberWanted = number+1;
376    } else if (switchType==2) {
377      numberWanted = max(2000,number/8);
378    } else {
379      int numberElements = model_->factorization()->numberElements();
380      double ratio = (double) numberElements/(double) model_->numberRows();
381      if (ratio<1.0) {
382        numberWanted = max(2000,number/20);
383      } else if (ratio<5.0) {
384        numberWanted = max(2000,number/10);
385        numberWanted = max(numberWanted,numberColumns/20);
386      } else if (ratio<10.0) {
387        numberWanted = max(2000,number/8);
388        numberWanted = max(numberWanted,numberColumns/20);
389      } else {
390        ratio = number * (ratio/80.0);
391        if (ratio>number) {
392          numberWanted=number+1;
393        } else {
394          numberWanted = max(2000,(int) ratio);
395          numberWanted = max(numberWanted,numberColumns/10);
396        }
397      }
398    }
399  }
400
401
402  double bestDj = 1.0e-30;
403  int bestSequence=-1;
404
405  int i,iSequence;
406  index = infeasible_->getIndices();
407  number = infeasible_->getNumElements();
408  // Re-sort infeasible every 100 pivots
409  if (0&&model_->factorization()->pivots()>0&&
410      (model_->factorization()->pivots()%100)==0) {
411    int nLook = model_->numberRows()+numberColumns;
412    number=0;
413    for (i=0;i<nLook;i++) {
414      if (infeas[i]) { 
415        if (fabs(infeas[i])>COIN_INDEXED_TINY_ELEMENT) 
416          index[number++]=i;
417        else
418          infeas[i]=0.0;
419      }
420    }
421    infeasible_->setNumElements(number);
422  }
423  if(model_->numberIterations()<model_->lastBadIteration()+200) {
424    // we can't really trust infeasibilities if there is dual error
425    double checkTolerance = 1.0e-8;
426    if (!model_->factorization()->pivots())
427      checkTolerance = 1.0e-6;
428    if (model_->largestDualError()>checkTolerance)
429      tolerance *= model_->largestDualError()/checkTolerance;
430    // But cap
431    tolerance = min(1000.0,tolerance);
432  }
433#ifdef CLP_DEBUG
434  if (model_->numberDualInfeasibilities()==1) 
435    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
436           model_->largestDualError(),model_,model_->messageHandler(),
437           number);
438#endif
439  // stop last one coming immediately
440  double saveOutInfeasibility=0.0;
441  if (sequenceOut>=0) {
442    saveOutInfeasibility = infeas[sequenceOut];
443    infeas[sequenceOut]=0.0;
444  }
445  tolerance *= tolerance; // as we are using squares
446
447  int iPass;
448  // Setup two passes
449  int start[4];
450  start[1]=number;
451  start[2]=0;
452  double dstart = ((double) number) * CoinDrand48();
453  start[0]=(int) dstart;
454  start[3]=start[0];
455  //double largestWeight=0.0;
456  //double smallestWeight=1.0e100;
457  for (iPass=0;iPass<2;iPass++) {
458    int end = start[2*iPass+1];
459    if (switchType<5) {
460      for (i=start[2*iPass];i<end;i++) {
461        iSequence = index[i];
462        double value = infeas[iSequence];
463        if (value>tolerance) {
464          double weight = weights_[iSequence];
465          //weight=1.0;
466          if (value>bestDj*weight) {
467            // check flagged variable and correct dj
468            if (!model_->flagged(iSequence)) {
469              bestDj=value/weight;
470              bestSequence = iSequence;
471            } else {
472              // just to make sure we don't exit before got something
473              numberWanted++;
474            }
475          }
476        }
477        numberWanted--;
478        if (!numberWanted)
479          break;
480      }
481    } else {
482      // Dantzig
483      for (i=start[2*iPass];i<end;i++) {
484        iSequence = index[i];
485        double value = infeas[iSequence];
486        if (value>tolerance) {
487          if (value>bestDj) {
488            // check flagged variable and correct dj
489            if (!model_->flagged(iSequence)) {
490              bestDj=value;
491              bestSequence = iSequence;
492            } else {
493              // just to make sure we don't exit before got something
494              numberWanted++;
495            }
496          }
497        }
498        numberWanted--;
499        if (!numberWanted)
500          break;
501      }
502    }
503    if (!numberWanted)
504      break;
505  }
506  if (sequenceOut>=0) {
507    infeas[sequenceOut]=saveOutInfeasibility;
508  }
509  /*if (model_->numberIterations()%100==0)
510    printf("%d best %g\n",bestSequence,bestDj);*/
511 
512#ifdef CLP_DEBUG
513  if (bestSequence>=0) {
514    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
515      assert(model_->reducedCost(bestSequence)<0.0);
516    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
517      assert(model_->reducedCost(bestSequence)>0.0);
518  }
519#endif
520  return bestSequence;
521}
522// Just update djs
523void 
524ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
525               CoinIndexedVector * spareRow1,
526               CoinIndexedVector * spareRow2,
527               CoinIndexedVector * spareColumn1,
528               CoinIndexedVector * spareColumn2)
529{
530  int iSection,j;
531  int number=0;
532  int * index;
533  double * updateBy;
534  double * reducedCost;
535  double tolerance=model_->currentDualTolerance();
536  // we can't really trust infeasibilities if there is dual error
537  // this coding has to mimic coding in checkDualSolution
538  double error = min(1.0e-3,model_->largestDualError());
539  // allow tolerance at least slightly bigger than standard
540  tolerance = tolerance +  error;
541  int pivotRow = model_->pivotRow();
542  double * infeas = infeasible_->denseVector();
543  model_->factorization()->updateColumnTranspose(spareRow2,updates);
544 
545  // put row of tableau in rowArray and columnArray
546  model_->clpMatrix()->transposeTimes(model_,-1.0,
547                                      updates,spareColumn2,spareColumn1);
548  // normal
549  for (iSection=0;iSection<2;iSection++) {
550   
551    reducedCost=model_->djRegion(iSection);
552    int addSequence;
553   
554    if (!iSection) {
555      number = updates->getNumElements();
556      index = updates->getIndices();
557      updateBy = updates->denseVector();
558      addSequence = model_->numberColumns();
559    } else {
560      number = spareColumn1->getNumElements();
561      index = spareColumn1->getIndices();
562      updateBy = spareColumn1->denseVector();
563      addSequence = 0;
564    }
565   
566    for (j=0;j<number;j++) {
567      int iSequence = index[j];
568      double value = reducedCost[iSequence];
569      value -= updateBy[iSequence];
570      updateBy[iSequence]=0.0;
571      reducedCost[iSequence] = value;
572      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
573     
574      switch(status) {
575       
576      case ClpSimplex::basic:
577        infeasible_->zero(iSequence+addSequence);
578      case ClpSimplex::isFixed:
579        break;
580      case ClpSimplex::isFree:
581      case ClpSimplex::superBasic:
582        if (fabs(value)>FREE_ACCEPT*tolerance) {
583          // we are going to bias towards free (but only if reasonable)
584          value *= FREE_BIAS;
585          // store square in list
586          if (infeas[iSequence+addSequence])
587            infeas[iSequence+addSequence] = value*value; // already there
588          else
589            infeasible_->quickAdd(iSequence+addSequence,value*value);
590        } else {
591          infeasible_->zero(iSequence+addSequence);
592        }
593        break;
594      case ClpSimplex::atUpperBound:
595        if (value>tolerance) {
596          // store square in list
597          if (infeas[iSequence+addSequence])
598            infeas[iSequence+addSequence] = value*value; // already there
599          else
600            infeasible_->quickAdd(iSequence+addSequence,value*value);
601        } else {
602          infeasible_->zero(iSequence+addSequence);
603        }
604        break;
605      case ClpSimplex::atLowerBound:
606        if (value<-tolerance) {
607          // store square in list
608          if (infeas[iSequence+addSequence])
609            infeas[iSequence+addSequence] = value*value; // already there
610          else
611            infeasible_->quickAdd(iSequence+addSequence,value*value);
612        } else {
613          infeasible_->zero(iSequence+addSequence);
614        }
615      }
616    }
617  }
618  updates->setNumElements(0);
619  spareColumn1->setNumElements(0);
620  if (pivotRow>=0) {
621    // make sure infeasibility on incoming is 0.0
622    int sequenceIn = model_->sequenceIn();
623    infeasible_->zero(sequenceIn);
624  }
625}
626// Update djs, weights for Devex
627void 
628ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
629               CoinIndexedVector * spareRow1,
630               CoinIndexedVector * spareRow2,
631               CoinIndexedVector * spareColumn1,
632               CoinIndexedVector * spareColumn2)
633{
634  int j;
635  int number=0;
636  int * index;
637  double * updateBy;
638  double * reducedCost;
639  // dj could be very small (or even zero - take care)
640  double dj = model_->dualIn();
641  double tolerance=model_->currentDualTolerance();
642  // we can't really trust infeasibilities if there is dual error
643  // this coding has to mimic coding in checkDualSolution
644  double error = min(1.0e-3,model_->largestDualError());
645  // allow tolerance at least slightly bigger than standard
646  tolerance = tolerance +  error;
647  // for weights update we use pivotSequence
648  // unset in case sub flip
649  assert (pivotSequence_>=0);
650  pivotSequence_=-1;
651  double * infeas = infeasible_->denseVector();
652  model_->factorization()->updateColumnTranspose(spareRow2,updates);
653  // and we can see if reference
654  double referenceIn=0.0;
655  int sequenceIn = model_->sequenceIn();
656  if (mode_!=1&&reference(sequenceIn))
657    referenceIn=1.0;
658  // save outgoing weight round update
659  double outgoingWeight=0.0;
660  int sequenceOut = model_->sequenceOut();
661  if (sequenceOut>=0)
662    outgoingWeight=weights_[sequenceOut];
663   
664  // put row of tableau in rowArray and columnArray
665  model_->clpMatrix()->transposeTimes(model_,-1.0,
666                                      updates,spareColumn2,spareColumn1);
667  // update weights
668  double * weight;
669  int numberColumns = model_->numberColumns();
670  double scaleFactor = -1.0/dj; // as formula is with 1.0
671  // rows
672  reducedCost=model_->djRegion(0);
673  int addSequence=model_->numberColumns();;
674   
675  number = updates->getNumElements();
676  index = updates->getIndices();
677  updateBy = updates->denseVector();
678  weight = weights_+numberColumns;
679  // Devex
680  for (j=0;j<number;j++) {
681    double thisWeight;
682    double pivot;
683    double value3;
684    int iSequence = index[j];
685    double value = reducedCost[iSequence];
686    double value2 = updateBy[iSequence];
687    updateBy[iSequence]=0.0;
688    value -= value2;
689    reducedCost[iSequence] = value;
690    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
691   
692    switch(status) {
693     
694    case ClpSimplex::basic:
695      infeasible_->zero(iSequence+addSequence);
696    case ClpSimplex::isFixed:
697      break;
698    case ClpSimplex::isFree:
699    case ClpSimplex::superBasic:
700      thisWeight = weight[iSequence];
701      // row has -1
702      pivot = value2*scaleFactor;
703      value3 = pivot * pivot*devex_;
704      if (reference(iSequence+numberColumns))
705        value3 += 1.0;
706      weight[iSequence] = max(0.99*thisWeight,value3);
707      if (fabs(value)>FREE_ACCEPT*tolerance) {
708        // we are going to bias towards free (but only if reasonable)
709        value *= FREE_BIAS;
710        // store square in list
711        if (infeas[iSequence+addSequence])
712          infeas[iSequence+addSequence] = value*value; // already there
713        else
714          infeasible_->quickAdd(iSequence+addSequence,value*value);
715      } else {
716        infeasible_->zero(iSequence+addSequence);
717      }
718      break;
719    case ClpSimplex::atUpperBound:
720      thisWeight = weight[iSequence];
721      // row has -1
722      pivot = value2*scaleFactor;
723      value3 = pivot * pivot*devex_;
724      if (reference(iSequence+numberColumns))
725        value3 += 1.0;
726      weight[iSequence] = max(0.99*thisWeight,value3);
727      if (value>tolerance) {
728        // store square in list
729        if (infeas[iSequence+addSequence])
730          infeas[iSequence+addSequence] = value*value; // already there
731        else
732          infeasible_->quickAdd(iSequence+addSequence,value*value);
733      } else {
734        infeasible_->zero(iSequence+addSequence);
735      }
736      break;
737    case ClpSimplex::atLowerBound:
738      thisWeight = weight[iSequence];
739      // row has -1
740      pivot = value2*scaleFactor;
741      value3 = pivot * pivot*devex_;
742      if (reference(iSequence+numberColumns))
743        value3 += 1.0;
744      weight[iSequence] = max(0.99*thisWeight,value3);
745      if (value<-tolerance) {
746        // store square in list
747        if (infeas[iSequence+addSequence])
748          infeas[iSequence+addSequence] = value*value; // already there
749        else
750          infeasible_->quickAdd(iSequence+addSequence,value*value);
751      } else {
752        infeasible_->zero(iSequence+addSequence);
753      }
754    }
755  }
756 
757  // columns
758  weight = weights_;
759 
760  scaleFactor = -scaleFactor;
761  reducedCost=model_->djRegion(1);
762  number = spareColumn1->getNumElements();
763  index = spareColumn1->getIndices();
764  updateBy = spareColumn1->denseVector();
765
766  // Devex
767 
768  for (j=0;j<number;j++) {
769    double thisWeight;
770    double pivot;
771    double value3;
772    int iSequence = index[j];
773    double value = reducedCost[iSequence];
774    double value2 = updateBy[iSequence];
775    value -= value2;
776    updateBy[iSequence]=0.0;
777    reducedCost[iSequence] = value;
778    ClpSimplex::Status status = model_->getStatus(iSequence);
779   
780    switch(status) {
781     
782    case ClpSimplex::basic:
783      infeasible_->zero(iSequence);
784    case ClpSimplex::isFixed:
785      break;
786    case ClpSimplex::isFree:
787    case ClpSimplex::superBasic:
788      thisWeight = weight[iSequence];
789      // row has -1
790      pivot = value2*scaleFactor;
791      value3 = pivot * pivot*devex_;
792      if (reference(iSequence))
793        value3 += 1.0;
794      weight[iSequence] = max(0.99*thisWeight,value3);
795      if (fabs(value)>FREE_ACCEPT*tolerance) {
796        // we are going to bias towards free (but only if reasonable)
797        value *= FREE_BIAS;
798        // store square in list
799        if (infeas[iSequence])
800          infeas[iSequence] = value*value; // already there
801        else
802          infeasible_->quickAdd(iSequence,value*value);
803      } else {
804        infeasible_->zero(iSequence);
805      }
806      break;
807    case ClpSimplex::atUpperBound:
808      thisWeight = weight[iSequence];
809      // row has -1
810      pivot = value2*scaleFactor;
811      value3 = pivot * pivot*devex_;
812      if (reference(iSequence))
813        value3 += 1.0;
814      weight[iSequence] = max(0.99*thisWeight,value3);
815      if (value>tolerance) {
816        // store square in list
817        if (infeas[iSequence])
818          infeas[iSequence] = value*value; // already there
819        else
820          infeasible_->quickAdd(iSequence,value*value);
821      } else {
822        infeasible_->zero(iSequence);
823      }
824      break;
825    case ClpSimplex::atLowerBound:
826      thisWeight = weight[iSequence];
827      // row has -1
828      pivot = value2*scaleFactor;
829      value3 = pivot * pivot*devex_;
830      if (reference(iSequence))
831        value3 += 1.0;
832      weight[iSequence] = max(0.99*thisWeight,value3);
833      if (value<-tolerance) {
834        // store square in list
835        if (infeas[iSequence])
836          infeas[iSequence] = value*value; // already there
837        else
838          infeasible_->quickAdd(iSequence,value*value);
839      } else {
840        infeasible_->zero(iSequence);
841      }
842    }
843  }
844  // restore outgoing weight
845  if (sequenceOut>=0)
846    weights_[sequenceOut]=outgoingWeight;
847  // make sure infeasibility on incoming is 0.0
848  infeasible_->zero(sequenceIn);
849  spareRow2->setNumElements(0);
850  //#define SOME_DEBUG_1
851#ifdef SOME_DEBUG_1
852  // check for accuracy
853  int iCheck=229;
854  //printf("weight for iCheck is %g\n",weights_[iCheck]);
855  int numberRows = model_->numberRows();
856  //int numberColumns = model_->numberColumns();
857  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
858    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
859        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
860      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
861  }
862#endif
863  updates->setNumElements(0);
864  spareColumn1->setNumElements(0);
865}
866// Update djs, weights for Steepest
867void 
868ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates,
869               CoinIndexedVector * spareRow1,
870               CoinIndexedVector * spareRow2,
871               CoinIndexedVector * spareColumn1,
872               CoinIndexedVector * spareColumn2)
873{
874  int j;
875  int number=0;
876  int * index;
877  double * updateBy;
878  double * reducedCost;
879  // dj could be very small (or even zero - take care)
880  double dj = model_->dualIn();
881  double tolerance=model_->currentDualTolerance();
882  // we can't really trust infeasibilities if there is dual error
883  // this coding has to mimic coding in checkDualSolution
884  double error = min(1.0e-3,model_->largestDualError());
885  // allow tolerance at least slightly bigger than standard
886  tolerance = tolerance +  error;
887  // for weights update we use pivotSequence
888  // unset in case sub flip
889  assert (pivotSequence_>=0);
890  pivotSequence_=-1;
891  double * infeas = infeasible_->denseVector();
892  model_->factorization()->updateColumnTranspose(spareRow2,updates);
893
894  model_->factorization()->updateColumnTranspose(spareRow2,
895                                                 alternateWeights_);
896  // and we can see if reference
897  double referenceIn=0.0;
898  int sequenceIn = model_->sequenceIn();
899  if (mode_!=1&&reference(sequenceIn))
900    referenceIn=1.0;
901  // save outgoing weight round update
902  double outgoingWeight=0.0;
903  int sequenceOut = model_->sequenceOut();
904  if (sequenceOut>=0)
905    outgoingWeight=weights_[sequenceOut];
906   
907  // put row of tableau in rowArray and columnArray
908  model_->clpMatrix()->transposeTimes(model_,-1.0,
909                                      updates,spareColumn2,spareColumn1);
910    // get subset which have nonzero tableau elements
911    // Luckily spareRow2 is long enough (rowArray_[3])
912    model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
913                                              spareColumn1,
914                                              spareRow2);
915  // update weights
916  double * weight;
917  double * other = alternateWeights_->denseVector();
918  int numberColumns = model_->numberColumns();
919  double scaleFactor = -1.0/dj; // as formula is with 1.0
920  // rows
921  reducedCost=model_->djRegion(0);
922  int addSequence=model_->numberColumns();;
923   
924  number = updates->getNumElements();
925  index = updates->getIndices();
926  updateBy = updates->denseVector();
927  weight = weights_+numberColumns;
928 
929  for (j=0;j<number;j++) {
930    double thisWeight;
931    double pivot;
932    double modification;
933    double pivotSquared;
934    int iSequence = index[j];
935    double value = reducedCost[iSequence];
936    double value2 = updateBy[iSequence];
937    value -= value2;
938    reducedCost[iSequence] = value;
939    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
940    updateBy[iSequence]=0.0;
941   
942    switch(status) {
943     
944    case ClpSimplex::basic:
945      infeasible_->zero(iSequence+addSequence);
946    case ClpSimplex::isFixed:
947      break;
948    case ClpSimplex::isFree:
949    case ClpSimplex::superBasic:
950      thisWeight = weight[iSequence];
951      // row has -1
952      pivot = value2*scaleFactor;
953      modification = other[iSequence];
954      pivotSquared = pivot * pivot;
955     
956      thisWeight += pivotSquared * devex_ + pivot * modification;
957      if (thisWeight<TRY_NORM) {
958        if (mode_==1) {
959          // steepest
960          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
961        } else {
962          // exact
963          thisWeight = referenceIn*pivotSquared;
964          if (reference(iSequence+numberColumns))
965            thisWeight += 1.0;
966          thisWeight = max(thisWeight,TRY_NORM);
967        }
968      }
969      weight[iSequence] = thisWeight;
970      if (fabs(value)>FREE_ACCEPT*tolerance) {
971        // we are going to bias towards free (but only if reasonable)
972        value *= FREE_BIAS;
973        // store square in list
974        if (infeas[iSequence+addSequence])
975          infeas[iSequence+addSequence] = value*value; // already there
976        else
977          infeasible_->quickAdd(iSequence+addSequence,value*value);
978      } else {
979        infeasible_->zero(iSequence+addSequence);
980      }
981      break;
982    case ClpSimplex::atUpperBound:
983      thisWeight = weight[iSequence];
984      // row has -1
985      pivot = value2*scaleFactor;
986      modification = other[iSequence];
987      pivotSquared = pivot * pivot;
988     
989      thisWeight += pivotSquared * devex_ + pivot * modification;
990      if (thisWeight<TRY_NORM) {
991        if (mode_==1) {
992          // steepest
993          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
994        } else {
995          // exact
996          thisWeight = referenceIn*pivotSquared;
997          if (reference(iSequence+numberColumns))
998            thisWeight += 1.0;
999          thisWeight = max(thisWeight,TRY_NORM);
1000        }
1001      }
1002      weight[iSequence] = thisWeight;
1003      if (value>tolerance) {
1004        // store square in list
1005        if (infeas[iSequence+addSequence])
1006          infeas[iSequence+addSequence] = value*value; // already there
1007        else
1008          infeasible_->quickAdd(iSequence+addSequence,value*value);
1009      } else {
1010        infeasible_->zero(iSequence+addSequence);
1011      }
1012      break;
1013    case ClpSimplex::atLowerBound:
1014      thisWeight = weight[iSequence];
1015      // row has -1
1016      pivot = value2*scaleFactor;
1017      modification = other[iSequence];
1018      pivotSquared = pivot * pivot;
1019     
1020      thisWeight += pivotSquared * devex_ + pivot * modification;
1021      if (thisWeight<TRY_NORM) {
1022        if (mode_==1) {
1023          // steepest
1024          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1025        } else {
1026          // exact
1027          thisWeight = referenceIn*pivotSquared;
1028          if (reference(iSequence+numberColumns))
1029            thisWeight += 1.0;
1030          thisWeight = max(thisWeight,TRY_NORM);
1031        }
1032      }
1033      weight[iSequence] = thisWeight;
1034      if (value<-tolerance) {
1035        // store square in list
1036        if (infeas[iSequence+addSequence])
1037          infeas[iSequence+addSequence] = value*value; // already there
1038        else
1039          infeasible_->quickAdd(iSequence+addSequence,value*value);
1040      } else {
1041        infeasible_->zero(iSequence+addSequence);
1042      }
1043    }
1044  }
1045 
1046  // columns
1047  weight = weights_;
1048 
1049  scaleFactor = -scaleFactor;
1050  reducedCost=model_->djRegion(1);
1051  number = spareColumn1->getNumElements();
1052  index = spareColumn1->getIndices();
1053  updateBy = spareColumn1->denseVector();
1054  double * updateBy2 = spareRow2->denseVector();
1055
1056  for (j=0;j<number;j++) {
1057    double thisWeight;
1058    double pivot;
1059    double modification;
1060    double pivotSquared;
1061    int iSequence = index[j];
1062    double value = reducedCost[iSequence];
1063    double value2 = updateBy[iSequence];
1064    updateBy[iSequence]=0.0;
1065    value -= value2;
1066    reducedCost[iSequence] = value;
1067    ClpSimplex::Status status = model_->getStatus(iSequence);
1068   
1069    switch(status) {
1070     
1071    case ClpSimplex::basic:
1072      infeasible_->zero(iSequence);
1073    case ClpSimplex::isFixed:
1074      updateBy2[iSequence]=0.0;
1075      break;
1076    case ClpSimplex::isFree:
1077    case ClpSimplex::superBasic:
1078      thisWeight = weight[iSequence];
1079      pivot = value2*scaleFactor;
1080      modification = updateBy2[iSequence];
1081      updateBy2[iSequence]=0.0;
1082      pivotSquared = pivot * pivot;
1083     
1084      thisWeight += pivotSquared * devex_ + pivot * modification;
1085      if (thisWeight<TRY_NORM) {
1086        if (mode_==1) {
1087          // steepest
1088          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1089        } else {
1090          // exact
1091          thisWeight = referenceIn*pivotSquared;
1092          if (reference(iSequence))
1093            thisWeight += 1.0;
1094          thisWeight = max(thisWeight,TRY_NORM);
1095        }
1096      }
1097      weight[iSequence] = thisWeight;
1098      if (fabs(value)>FREE_ACCEPT*tolerance) {
1099        // we are going to bias towards free (but only if reasonable)
1100        value *= FREE_BIAS;
1101        // store square in list
1102        if (infeas[iSequence])
1103          infeas[iSequence] = value*value; // already there
1104        else
1105          infeasible_->quickAdd(iSequence,value*value);
1106      } else {
1107        infeasible_->zero(iSequence);
1108      }
1109      break;
1110    case ClpSimplex::atUpperBound:
1111      thisWeight = weight[iSequence];
1112      pivot = value2*scaleFactor;
1113      modification = updateBy2[iSequence];
1114      updateBy2[iSequence]=0.0;
1115      pivotSquared = pivot * pivot;
1116     
1117      thisWeight += pivotSquared * devex_ + pivot * modification;
1118      if (thisWeight<TRY_NORM) {
1119        if (mode_==1) {
1120          // steepest
1121          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1122        } else {
1123          // exact
1124          thisWeight = referenceIn*pivotSquared;
1125          if (reference(iSequence))
1126            thisWeight += 1.0;
1127          thisWeight = max(thisWeight,TRY_NORM);
1128        }
1129      }
1130      weight[iSequence] = thisWeight;
1131      if (value>tolerance) {
1132        // store square in list
1133        if (infeas[iSequence])
1134          infeas[iSequence] = value*value; // already there
1135        else
1136          infeasible_->quickAdd(iSequence,value*value);
1137      } else {
1138        infeasible_->zero(iSequence);
1139      }
1140      break;
1141    case ClpSimplex::atLowerBound:
1142      thisWeight = weight[iSequence];
1143      pivot = value2*scaleFactor;
1144      modification = updateBy2[iSequence];
1145      updateBy2[iSequence]=0.0;
1146      pivotSquared = pivot * pivot;
1147     
1148      thisWeight += pivotSquared * devex_ + pivot * modification;
1149      if (thisWeight<TRY_NORM) {
1150        if (mode_==1) {
1151          // steepest
1152          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1153        } else {
1154          // exact
1155          thisWeight = referenceIn*pivotSquared;
1156          if (reference(iSequence))
1157            thisWeight += 1.0;
1158          thisWeight = max(thisWeight,TRY_NORM);
1159        }
1160      }
1161      weight[iSequence] = thisWeight;
1162      if (value<-tolerance) {
1163        // store square in list
1164        if (infeas[iSequence])
1165          infeas[iSequence] = value*value; // already there
1166        else
1167          infeasible_->quickAdd(iSequence,value*value);
1168      } else {
1169        infeasible_->zero(iSequence);
1170      }
1171    }
1172  }
1173  // restore outgoing weight
1174  if (sequenceOut>=0)
1175    weights_[sequenceOut]=outgoingWeight;
1176  // make sure infeasibility on incoming is 0.0
1177  infeasible_->zero(sequenceIn);
1178  alternateWeights_->clear();
1179  spareRow2->setNumElements(0);
1180  //#define SOME_DEBUG_1
1181#ifdef SOME_DEBUG_1
1182  // check for accuracy
1183  int iCheck=229;
1184  //printf("weight for iCheck is %g\n",weights_[iCheck]);
1185  int numberRows = model_->numberRows();
1186  //int numberColumns = model_->numberColumns();
1187  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
1188    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
1189        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
1190      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
1191  }
1192#endif
1193  updates->setNumElements(0);
1194  spareColumn1->setNumElements(0);
1195}
1196// Update djs, weights for Devex
1197void 
1198ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
1199               CoinIndexedVector * spareRow1,
1200               CoinIndexedVector * spareRow2,
1201               CoinIndexedVector * spareColumn1,
1202               CoinIndexedVector * spareColumn2)
1203{
1204  int iSection,j;
1205  int number=0;
1206  int * index;
1207  double * updateBy;
1208  double * reducedCost;
1209  // dj could be very small (or even zero - take care)
1210  double dj = model_->dualIn();
1211  double tolerance=model_->currentDualTolerance();
1212  // we can't really trust infeasibilities if there is dual error
1213  // this coding has to mimic coding in checkDualSolution
1214  double error = min(1.0e-3,model_->largestDualError());
1215  // allow tolerance at least slightly bigger than standard
1216  tolerance = tolerance +  error;
1217  int pivotRow = model_->pivotRow();
1218  double * infeas = infeasible_->denseVector();
1219  model_->factorization()->updateColumnTranspose(spareRow2,updates);
1220 
1221  // put row of tableau in rowArray and columnArray
1222  model_->clpMatrix()->transposeTimes(model_,-1.0,
1223                                      updates,spareColumn2,spareColumn1);
1224  // normal
1225  for (iSection=0;iSection<2;iSection++) {
1226   
1227    reducedCost=model_->djRegion(iSection);
1228    int addSequence;
1229   
1230    if (!iSection) {
1231      number = updates->getNumElements();
1232      index = updates->getIndices();
1233      updateBy = updates->denseVector();
1234      addSequence = model_->numberColumns();
1235    } else {
1236      number = spareColumn1->getNumElements();
1237      index = spareColumn1->getIndices();
1238      updateBy = spareColumn1->denseVector();
1239      addSequence = 0;
1240    }
1241   
1242    for (j=0;j<number;j++) {
1243      int iSequence = index[j];
1244      double value = reducedCost[iSequence];
1245      value -= updateBy[iSequence];
1246      reducedCost[iSequence] = value;
1247      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
1248     
1249      switch(status) {
1250       
1251      case ClpSimplex::basic:
1252        infeasible_->zero(iSequence+addSequence);
1253      case ClpSimplex::isFixed:
1254        break;
1255      case ClpSimplex::isFree:
1256      case ClpSimplex::superBasic:
1257        if (fabs(value)>FREE_ACCEPT*tolerance) {
1258          // we are going to bias towards free (but only if reasonable)
1259          value *= FREE_BIAS;
1260          // store square in list
1261          if (infeas[iSequence+addSequence])
1262            infeas[iSequence+addSequence] = value*value; // already there
1263          else
1264            infeasible_->quickAdd(iSequence+addSequence,value*value);
1265        } else {
1266          infeasible_->zero(iSequence+addSequence);
1267        }
1268        break;
1269      case ClpSimplex::atUpperBound:
1270        if (value>tolerance) {
1271          // store square in list
1272          if (infeas[iSequence+addSequence])
1273            infeas[iSequence+addSequence] = value*value; // already there
1274          else
1275            infeasible_->quickAdd(iSequence+addSequence,value*value);
1276        } else {
1277          infeasible_->zero(iSequence+addSequence);
1278        }
1279        break;
1280      case ClpSimplex::atLowerBound:
1281        if (value<-tolerance) {
1282          // store square in list
1283          if (infeas[iSequence+addSequence])
1284            infeas[iSequence+addSequence] = value*value; // already there
1285          else
1286            infeasible_->quickAdd(iSequence+addSequence,value*value);
1287        } else {
1288          infeasible_->zero(iSequence+addSequence);
1289        }
1290      }
1291    }
1292  }
1293  // we can zero out as will have to get pivot row
1294  updates->clear();
1295  spareColumn1->clear();
1296  // make sure infeasibility on incoming is 0.0
1297  int sequenceIn = model_->sequenceIn();
1298  infeasible_->zero(sequenceIn);
1299  // for weights update we use pivotSequence
1300  assert (pivotSequence_>=0);
1301  pivotRow = pivotSequence_;
1302  // unset in case sub flip
1303  pivotSequence_=-1;
1304  // make sure infeasibility on incoming is 0.0
1305  const int * pivotVariable = model_->pivotVariable();
1306  sequenceIn = pivotVariable[pivotRow];
1307  infeasible_->zero(sequenceIn);
1308  // and we can see if reference
1309  double referenceIn=0.0;
1310  if (mode_!=1&&reference(sequenceIn))
1311    referenceIn=1.0;
1312  // save outgoing weight round update
1313  double outgoingWeight=0.0;
1314  int sequenceOut = model_->sequenceOut();
1315  if (sequenceOut>=0)
1316    outgoingWeight=weights_[sequenceOut];
1317  // update weights
1318  updates->setNumElements(0);
1319  spareColumn1->setNumElements(0);
1320  // might as well set dj to 1
1321  dj = 1.0;
1322  updates->insert(pivotRow,-dj);
1323  model_->factorization()->updateColumnTranspose(spareRow2,updates);
1324  // put row of tableau in rowArray and columnArray
1325  model_->clpMatrix()->transposeTimes(model_,-1.0,
1326                                      updates,spareColumn2,spareColumn1);
1327  double * weight;
1328  int numberColumns = model_->numberColumns();
1329  double scaleFactor = -1.0/dj; // as formula is with 1.0
1330  // rows
1331  number = updates->getNumElements();
1332  index = updates->getIndices();
1333  updateBy = updates->denseVector();
1334  weight = weights_+numberColumns;
1335 
1336  assert (devex_>0.0);
1337  for (j=0;j<number;j++) {
1338    int iSequence = index[j];
1339    double thisWeight = weight[iSequence];
1340    // row has -1
1341    double pivot = updateBy[iSequence]*scaleFactor;
1342    updateBy[iSequence]=0.0;
1343    double value = pivot * pivot*devex_;
1344    if (reference(iSequence+numberColumns))
1345      value += 1.0;
1346    weight[iSequence] = max(0.99*thisWeight,value);
1347  }
1348 
1349  // columns
1350  weight = weights_;
1351 
1352  scaleFactor = -scaleFactor;
1353 
1354  number = spareColumn1->getNumElements();
1355  index = spareColumn1->getIndices();
1356  updateBy = spareColumn1->denseVector();
1357  for (j=0;j<number;j++) {
1358    int iSequence = index[j];
1359    double thisWeight = weight[iSequence];
1360    // row has -1
1361    double pivot = updateBy[iSequence]*scaleFactor;
1362    updateBy[iSequence]=0.0;
1363    double value = pivot * pivot*devex_;
1364    if (reference(iSequence))
1365      value += 1.0;
1366    weight[iSequence] = max(0.99*thisWeight,value);
1367  }
1368  // restore outgoing weight
1369  if (sequenceOut>=0)
1370    weights_[sequenceOut]=outgoingWeight;
1371  spareColumn2->setNumElements(0);
1372  //#define SOME_DEBUG_1
1373#ifdef SOME_DEBUG_1
1374  // check for accuracy
1375  int iCheck=229;
1376  //printf("weight for iCheck is %g\n",weights_[iCheck]);
1377  int numberRows = model_->numberRows();
1378  //int numberColumns = model_->numberColumns();
1379  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
1380    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
1381        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
1382      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
1383  }
1384#endif
1385  updates->setNumElements(0);
1386  spareColumn1->setNumElements(0);
1387}
1388// Update djs, weights for Steepest
1389void 
1390ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates,
1391               CoinIndexedVector * spareRow1,
1392               CoinIndexedVector * spareRow2,
1393               CoinIndexedVector * spareColumn1,
1394               CoinIndexedVector * spareColumn2)
1395{
1396  int iSection,j;
1397  int number=0;
1398  int * index;
1399  double * updateBy;
1400  double * reducedCost;
1401  // dj could be very small (or even zero - take care)
1402  double dj = model_->dualIn();
1403  double tolerance=model_->currentDualTolerance();
1404  // we can't really trust infeasibilities if there is dual error
1405  // this coding has to mimic coding in checkDualSolution
1406  double error = min(1.0e-3,model_->largestDualError());
1407  // allow tolerance at least slightly bigger than standard
1408  tolerance = tolerance +  error;
1409  int pivotRow = model_->pivotRow();
1410  double * infeas = infeasible_->denseVector();
1411  model_->factorization()->updateColumnTranspose(spareRow2,updates);
1412 
1413  // put row of tableau in rowArray and columnArray
1414  model_->clpMatrix()->transposeTimes(model_,-1.0,
1415                                      updates,spareColumn2,spareColumn1);
1416  // normal
1417  for (iSection=0;iSection<2;iSection++) {
1418   
1419    reducedCost=model_->djRegion(iSection);
1420    int addSequence;
1421   
1422    if (!iSection) {
1423      number = updates->getNumElements();
1424      index = updates->getIndices();
1425      updateBy = updates->denseVector();
1426      addSequence = model_->numberColumns();
1427    } else {
1428      number = spareColumn1->getNumElements();
1429      index = spareColumn1->getIndices();
1430      updateBy = spareColumn1->denseVector();
1431      addSequence = 0;
1432    }
1433   
1434    for (j=0;j<number;j++) {
1435      int iSequence = index[j];
1436      double value = reducedCost[iSequence];
1437      value -= updateBy[iSequence];
1438      reducedCost[iSequence] = value;
1439      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
1440     
1441      switch(status) {
1442       
1443      case ClpSimplex::basic:
1444        infeasible_->zero(iSequence+addSequence);
1445      case ClpSimplex::isFixed:
1446        break;
1447      case ClpSimplex::isFree:
1448      case ClpSimplex::superBasic:
1449        if (fabs(value)>FREE_ACCEPT*tolerance) {
1450          // we are going to bias towards free (but only if reasonable)
1451          value *= FREE_BIAS;
1452          // store square in list
1453          if (infeas[iSequence+addSequence])
1454            infeas[iSequence+addSequence] = value*value; // already there
1455          else
1456            infeasible_->quickAdd(iSequence+addSequence,value*value);
1457        } else {
1458          infeasible_->zero(iSequence+addSequence);
1459        }
1460        break;
1461      case ClpSimplex::atUpperBound:
1462        if (value>tolerance) {
1463          // store square in list
1464          if (infeas[iSequence+addSequence])
1465            infeas[iSequence+addSequence] = value*value; // already there
1466          else
1467            infeasible_->quickAdd(iSequence+addSequence,value*value);
1468        } else {
1469          infeasible_->zero(iSequence+addSequence);
1470        }
1471        break;
1472      case ClpSimplex::atLowerBound:
1473        if (value<-tolerance) {
1474          // store square in list
1475          if (infeas[iSequence+addSequence])
1476            infeas[iSequence+addSequence] = value*value; // already there
1477          else
1478            infeasible_->quickAdd(iSequence+addSequence,value*value);
1479        } else {
1480          infeasible_->zero(iSequence+addSequence);
1481        }
1482      }
1483    }
1484  }
1485  // we can zero out as will have to get pivot row
1486  // ***** move
1487  updates->clear();
1488  spareColumn1->clear();
1489  if (pivotRow>=0) {
1490    // make sure infeasibility on incoming is 0.0
1491    int sequenceIn = model_->sequenceIn();
1492    infeasible_->zero(sequenceIn);
1493  }
1494  // for weights update we use pivotSequence
1495  pivotRow = pivotSequence_;
1496  // unset in case sub flip
1497  pivotSequence_=-1;
1498  if (pivotRow>=0) {
1499    // make sure infeasibility on incoming is 0.0
1500    const int * pivotVariable = model_->pivotVariable();
1501    int sequenceIn = pivotVariable[pivotRow];
1502    infeasible_->zero(sequenceIn);
1503    // and we can see if reference
1504    double referenceIn=0.0;
1505    if (mode_!=1&&reference(sequenceIn))
1506      referenceIn=1.0;
1507    // save outgoing weight round update
1508    double outgoingWeight=0.0;
1509    int sequenceOut = model_->sequenceOut();
1510    if (sequenceOut>=0)
1511      outgoingWeight=weights_[sequenceOut];
1512    // update weights
1513    updates->setNumElements(0);
1514    spareColumn1->setNumElements(0);
1515    // might as well set dj to 1
1516    dj = 1.0;
1517    updates->insert(pivotRow,-dj);
1518    model_->factorization()->updateColumnTranspose(spareRow2,updates);
1519    // put row of tableau in rowArray and columnArray
1520    model_->clpMatrix()->transposeTimes(model_,-1.0,
1521                                        updates,spareColumn2,spareColumn1);
1522    double * weight;
1523    double * other = alternateWeights_->denseVector();
1524    int numberColumns = model_->numberColumns();
1525    double scaleFactor = -1.0/dj; // as formula is with 1.0
1526    // rows
1527    number = updates->getNumElements();
1528    index = updates->getIndices();
1529    updateBy = updates->denseVector();
1530    weight = weights_+numberColumns;
1531   
1532    if (mode_<4||numberSwitched_>1) {
1533      // Exact
1534      // now update weight update array
1535      model_->factorization()->updateColumnTranspose(spareRow2,
1536                                                     alternateWeights_);
1537      for (j=0;j<number;j++) {
1538        int iSequence = index[j];
1539        double thisWeight = weight[iSequence];
1540        // row has -1
1541        double pivot = updateBy[iSequence]*scaleFactor;
1542        updateBy[iSequence]=0.0;
1543        double modification = other[iSequence];
1544        double pivotSquared = pivot * pivot;
1545       
1546        thisWeight += pivotSquared * devex_ + pivot * modification;
1547        if (thisWeight<TRY_NORM) {
1548          if (mode_==1) {
1549            // steepest
1550            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1551          } else {
1552            // exact
1553            thisWeight = referenceIn*pivotSquared;
1554            if (reference(iSequence+numberColumns))
1555              thisWeight += 1.0;
1556            thisWeight = max(thisWeight,TRY_NORM);
1557          }
1558        }
1559        weight[iSequence] = thisWeight;
1560      }
1561    } else if (mode_==4) {
1562      // Devex
1563      assert (devex_>0.0);
1564      for (j=0;j<number;j++) {
1565        int iSequence = index[j];
1566        double thisWeight = weight[iSequence];
1567        // row has -1
1568        double pivot = updateBy[iSequence]*scaleFactor;
1569        updateBy[iSequence]=0.0;
1570        double value = pivot * pivot*devex_;
1571        if (reference(iSequence+numberColumns))
1572          value += 1.0;
1573        weight[iSequence] = max(0.99*thisWeight,value);
1574      }
1575    }
1576   
1577    // columns
1578    weight = weights_;
1579   
1580    scaleFactor = -scaleFactor;
1581   
1582    number = spareColumn1->getNumElements();
1583    index = spareColumn1->getIndices();
1584    updateBy = spareColumn1->denseVector();
1585    if (mode_<4||numberSwitched_>1) {
1586      // Exact
1587      // get subset which have nonzero tableau elements
1588      model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
1589                                                spareColumn1,
1590                                                spareColumn2);
1591      double * updateBy2 = spareColumn2->denseVector();
1592      for (j=0;j<number;j++) {
1593        int iSequence = index[j];
1594        double thisWeight = weight[iSequence];
1595        double pivot = updateBy[iSequence]*scaleFactor;
1596        updateBy[iSequence]=0.0;
1597        double modification = updateBy2[iSequence];
1598        updateBy2[iSequence]=0.0;
1599        double pivotSquared = pivot * pivot;
1600       
1601        thisWeight += pivotSquared * devex_ + pivot * modification;
1602        if (thisWeight<TRY_NORM) {
1603          if (mode_==1) {
1604            // steepest
1605            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1606          } else {
1607            // exact
1608            thisWeight = referenceIn*pivotSquared;
1609            if (reference(iSequence))
1610              thisWeight += 1.0;
1611            thisWeight = max(thisWeight,TRY_NORM);
1612          }
1613        }
1614        weight[iSequence] = thisWeight;
1615      }
1616    } else if (mode_==4) {
1617      // Devex
1618      for (j=0;j<number;j++) {
1619        int iSequence = index[j];
1620        double thisWeight = weight[iSequence];
1621        // row has -1
1622        double pivot = updateBy[iSequence]*scaleFactor;
1623        updateBy[iSequence]=0.0;
1624        double value = pivot * pivot*devex_;
1625        if (reference(iSequence))
1626          value += 1.0;
1627        weight[iSequence] = max(0.99*thisWeight,value);
1628      }
1629    }
1630    // restore outgoing weight
1631    if (sequenceOut>=0)
1632      weights_[sequenceOut]=outgoingWeight;
1633    alternateWeights_->clear();
1634    spareColumn2->setNumElements(0);
1635    //#define SOME_DEBUG_1
1636#ifdef SOME_DEBUG_1
1637    // check for accuracy
1638    int iCheck=229;
1639    //printf("weight for iCheck is %g\n",weights_[iCheck]);
1640    int numberRows = model_->numberRows();
1641    //int numberColumns = model_->numberColumns();
1642    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
1643      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
1644          !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
1645        checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
1646    }
1647#endif
1648  }
1649  updates->setNumElements(0);
1650  spareColumn1->setNumElements(0);
1651}
1652// Update weights for Devex
1653void 
1654ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
1655               CoinIndexedVector * spareRow1,
1656               CoinIndexedVector * spareRow2,
1657               CoinIndexedVector * spareColumn1,
1658               CoinIndexedVector * spareColumn2)
1659{
1660  int j;
1661  int number=0;
1662  int * index;
1663  double * updateBy;
1664  // dj could be very small (or even zero - take care)
1665  double dj = model_->dualIn();
1666  double tolerance=model_->currentDualTolerance();
1667  // we can't really trust infeasibilities if there is dual error
1668  // this coding has to mimic coding in checkDualSolution
1669  double error = min(1.0e-3,model_->largestDualError());
1670  // allow tolerance at least slightly bigger than standard
1671  tolerance = tolerance +  error;
1672  int pivotRow = model_->pivotRow();
1673  // for weights update we use pivotSequence
1674  pivotRow = pivotSequence_;
1675  assert (pivotRow>=0);
1676  // make sure infeasibility on incoming is 0.0
1677  const int * pivotVariable = model_->pivotVariable();
1678  int sequenceIn = pivotVariable[pivotRow];
1679  infeasible_->zero(sequenceIn);
1680  // and we can see if reference
1681  double referenceIn=0.0;
1682  if (mode_!=1&&reference(sequenceIn))
1683    referenceIn=1.0;
1684  // save outgoing weight round update
1685  double outgoingWeight=0.0;
1686  int sequenceOut = model_->sequenceOut();
1687  if (sequenceOut>=0)
1688    outgoingWeight=weights_[sequenceOut];
1689  assert (!updates->getNumElements());
1690  assert (!spareColumn1->getNumElements());
1691  // unset in case sub flip
1692  pivotSequence_=-1;
1693  // might as well set dj to 1
1694  dj = 1.0;
1695  updates->insert(pivotRow,-dj);
1696  model_->factorization()->updateColumnTranspose(spareRow2,updates);
1697  // put row of tableau in rowArray and columnArray
1698  model_->clpMatrix()->transposeTimes(model_,-1.0,
1699                                      updates,spareColumn2,spareColumn1);
1700  double * weight;
1701  int numberColumns = model_->numberColumns();
1702  double scaleFactor = -1.0/dj; // as formula is with 1.0
1703  // rows
1704  number = updates->getNumElements();
1705  index = updates->getIndices();
1706  updateBy = updates->denseVector();
1707  weight = weights_+numberColumns;
1708 
1709  // Devex
1710  assert (devex_>0.0);
1711  for (j=0;j<number;j++) {
1712    int iSequence = index[j];
1713    double thisWeight = weight[iSequence];
1714    // row has -1
1715    double pivot = updateBy[iSequence]*scaleFactor;
1716    updateBy[iSequence]=0.0;
1717    double value = pivot * pivot*devex_;
1718    if (reference(iSequence+numberColumns))
1719      value += 1.0;
1720    weight[iSequence] = max(0.99*thisWeight,value);
1721  }
1722 
1723  // columns
1724  weight = weights_;
1725 
1726  scaleFactor = -scaleFactor;
1727 
1728  number = spareColumn1->getNumElements();
1729  index = spareColumn1->getIndices();
1730  updateBy = spareColumn1->denseVector();
1731  // Devex
1732  for (j=0;j<number;j++) {
1733    int iSequence = index[j];
1734    double thisWeight = weight[iSequence];
1735    // row has -1
1736    double pivot = updateBy[iSequence]*scaleFactor;
1737    updateBy[iSequence]=0.0;
1738    double value = pivot * pivot*devex_;
1739    if (reference(iSequence))
1740      value += 1.0;
1741    weight[iSequence] = max(0.99*thisWeight,value);
1742  }
1743  // restore outgoing weight
1744  if (sequenceOut>=0)
1745    weights_[sequenceOut]=outgoingWeight;
1746  spareColumn2->setNumElements(0);
1747  //#define SOME_DEBUG_1
1748#ifdef SOME_DEBUG_1
1749  // check for accuracy
1750  int iCheck=229;
1751  //printf("weight for iCheck is %g\n",weights_[iCheck]);
1752  int numberRows = model_->numberRows();
1753  //int numberColumns = model_->numberColumns();
1754  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
1755    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
1756        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
1757      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
1758  }
1759#endif
1760  updates->setNumElements(0);
1761  spareColumn1->setNumElements(0);
1762}
1763// Update weights for Steepest
1764void 
1765ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates,
1766               CoinIndexedVector * spareRow1,
1767               CoinIndexedVector * spareRow2,
1768               CoinIndexedVector * spareColumn1,
1769               CoinIndexedVector * spareColumn2)
1770{
1771  int j;
1772  int number=0;
1773  int * index;
1774  double * updateBy;
1775  // dj could be very small (or even zero - take care)
1776  double dj = model_->dualIn();
1777  double tolerance=model_->currentDualTolerance();
1778  // we can't really trust infeasibilities if there is dual error
1779  // this coding has to mimic coding in checkDualSolution
1780  double error = min(1.0e-3,model_->largestDualError());
1781  // allow tolerance at least slightly bigger than standard
1782  tolerance = tolerance +  error;
1783  int pivotRow = model_->pivotRow();
1784  // for weights update we use pivotSequence
1785  pivotRow = pivotSequence_;
1786  // unset in case sub flip
1787  pivotSequence_=-1;
1788  assert (pivotRow>=0);
1789  // make sure infeasibility on incoming is 0.0
1790  const int * pivotVariable = model_->pivotVariable();
1791  int sequenceIn = pivotVariable[pivotRow];
1792  infeasible_->zero(sequenceIn);
1793  // and we can see if reference
1794  double referenceIn=0.0;
1795  if (mode_!=1&&reference(sequenceIn))
1796    referenceIn=1.0;
1797  // save outgoing weight round update
1798  double outgoingWeight=0.0;
1799  int sequenceOut = model_->sequenceOut();
1800  if (sequenceOut>=0)
1801    outgoingWeight=weights_[sequenceOut];
1802  assert (!updates->getNumElements());
1803  assert (!spareColumn1->getNumElements());
1804  // update weights
1805  //updates->setNumElements(0);
1806  //spareColumn1->setNumElements(0);
1807  // might as well set dj to 1
1808  dj = 1.0;
1809  updates->insert(pivotRow,-dj);
1810  model_->factorization()->updateColumnTranspose(spareRow2,updates);
1811  // put row of tableau in rowArray and columnArray
1812  model_->clpMatrix()->transposeTimes(model_,-1.0,
1813                                      updates,spareColumn2,spareColumn1);
1814  double * weight;
1815  double * other = alternateWeights_->denseVector();
1816  int numberColumns = model_->numberColumns();
1817  double scaleFactor = -1.0/dj; // as formula is with 1.0
1818  // rows
1819  number = updates->getNumElements();
1820  index = updates->getIndices();
1821  updateBy = updates->denseVector();
1822  weight = weights_+numberColumns;
1823 
1824  // Exact
1825  // now update weight update array
1826  model_->factorization()->updateColumnTranspose(spareRow2,
1827                                                 alternateWeights_);
1828  for (j=0;j<number;j++) {
1829    int iSequence = index[j];
1830    double thisWeight = weight[iSequence];
1831    // row has -1
1832    double pivot = updateBy[iSequence]*scaleFactor;
1833    updateBy[iSequence]=0.0;
1834    double modification = other[iSequence];
1835    double pivotSquared = pivot * pivot;
1836   
1837    thisWeight += pivotSquared * devex_ + pivot * modification;
1838    if (thisWeight<TRY_NORM) {
1839      if (mode_==1) {
1840        // steepest
1841        thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1842      } else {
1843        // exact
1844        thisWeight = referenceIn*pivotSquared;
1845        if (reference(iSequence+numberColumns))
1846          thisWeight += 1.0;
1847        thisWeight = max(thisWeight,TRY_NORM);
1848      }
1849    }
1850    weight[iSequence] = thisWeight;
1851  }
1852 
1853  // columns
1854  weight = weights_;
1855 
1856  scaleFactor = -scaleFactor;
1857 
1858  number = spareColumn1->getNumElements();
1859  index = spareColumn1->getIndices();
1860  updateBy = spareColumn1->denseVector();
1861  // Exact
1862  // get subset which have nonzero tableau elements
1863  model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
1864                                            spareColumn1,
1865                                            spareColumn2);
1866  double * updateBy2 = spareColumn2->denseVector();
1867  for (j=0;j<number;j++) {
1868    int iSequence = index[j];
1869    double thisWeight = weight[iSequence];
1870    double pivot = updateBy[iSequence]*scaleFactor;
1871    updateBy[iSequence]=0.0;
1872    double modification = updateBy2[iSequence];
1873    updateBy2[iSequence]=0.0;
1874    double pivotSquared = pivot * pivot;
1875   
1876    thisWeight += pivotSquared * devex_ + pivot * modification;
1877    if (thisWeight<TRY_NORM) {
1878      if (mode_==1) {
1879        // steepest
1880        thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
1881      } else {
1882        // exact
1883        thisWeight = referenceIn*pivotSquared;
1884        if (reference(iSequence))
1885          thisWeight += 1.0;
1886        thisWeight = max(thisWeight,TRY_NORM);
1887      }
1888    }
1889    weight[iSequence] = thisWeight;
1890  }
1891  // restore outgoing weight
1892  if (sequenceOut>=0)
1893    weights_[sequenceOut]=outgoingWeight;
1894  alternateWeights_->clear();
1895  spareColumn2->setNumElements(0);
1896  //#define SOME_DEBUG_1
1897#ifdef SOME_DEBUG_1
1898  // check for accuracy
1899  int iCheck=229;
1900  //printf("weight for iCheck is %g\n",weights_[iCheck]);
1901  int numberRows = model_->numberRows();
1902  //int numberColumns = model_->numberColumns();
1903  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
1904    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
1905        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
1906      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
1907  }
1908#endif
1909  updates->setNumElements(0);
1910  spareColumn1->setNumElements(0);
1911}
1912// Returns pivot column, -1 if none
1913int 
1914ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates,
1915                                    CoinIndexedVector * spareRow1,
1916                                    CoinIndexedVector * spareRow2,
1917                                    CoinIndexedVector * spareColumn1,
1918                                    CoinIndexedVector * spareColumn2)
1919{
1920  assert(model_);
1921  int iSection,j;
1922  int number=0;
1923  int * index;
1924  double * updateBy;
1925  double * reducedCost;
1926  // dj could be very small (or even zero - take care)
1927  double dj = model_->dualIn();
1928  double tolerance=model_->currentDualTolerance();
1929  // we can't really trust infeasibilities if there is dual error
1930  // this coding has to mimic coding in checkDualSolution
1931  double error = min(1.0e-3,model_->largestDualError());
1932  // allow tolerance at least slightly bigger than standard
1933  tolerance = tolerance +  error;
1934  int pivotRow = model_->pivotRow();
1935  int anyUpdates;
1936  double * infeas = infeasible_->denseVector();
1937
1938  // Local copy of mode so can decide what to do
1939  int switchType;
1940  if (mode_==4) 
1941    switchType = 5-numberSwitched_;
1942  else
1943    switchType = mode_;
1944  /* switchType -
1945     0 - all exact devex
1946     1 - all steepest
1947     2 - some exact devex
1948     3 - auto some exact devex
1949     4 - devex
1950     5 - dantzig
1951  */
1952  if (updates->getNumElements()) {
1953    // would have to have two goes for devex, three for steepest
1954    anyUpdates=2;
1955    // add in pivot contribution
1956    if (pivotRow>=0) 
1957      updates->add(pivotRow,-dj);
1958  } else if (pivotRow>=0) {
1959    if (fabs(dj)>1.0e-15) {
1960      // some dj
1961      updates->insert(pivotRow,-dj);
1962      if (fabs(dj)>1.0e-6) {
1963        // reasonable size
1964        anyUpdates=1;
1965      } else {
1966        // too small
1967        anyUpdates=2;
1968      }
1969    } else {
1970      // zero dj
1971      anyUpdates=-1;
1972    }
1973  } else if (pivotSequence_>=0){
1974    // just after re-factorization
1975    anyUpdates=-1;
1976  } else {
1977    // sub flip - nothing to do
1978    anyUpdates=0;
1979  }
1980
1981  if (anyUpdates>0) {
1982    model_->factorization()->updateColumnTranspose(spareRow2,updates);
1983   
1984    // put row of tableau in rowArray and columnArray
1985    model_->clpMatrix()->transposeTimes(model_,-1.0,
1986                                        updates,spareColumn2,spareColumn1);
1987    // normal
1988    for (iSection=0;iSection<2;iSection++) {
1989     
1990      reducedCost=model_->djRegion(iSection);
1991      int addSequence;
1992     
1993      if (!iSection) {
1994        number = updates->getNumElements();
1995        index = updates->getIndices();
1996        updateBy = updates->denseVector();
1997        addSequence = model_->numberColumns();
1998      } else {
1999        number = spareColumn1->getNumElements();
2000        index = spareColumn1->getIndices();
2001        updateBy = spareColumn1->denseVector();
2002        addSequence = 0;
2003      }
2004      if (!model_->nonLinearCost()->lookBothWays()) {
2005
2006        for (j=0;j<number;j++) {
2007          int iSequence = index[j];
2008          double value = reducedCost[iSequence];
2009          value -= updateBy[iSequence];
2010          reducedCost[iSequence] = value;
2011          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
2012         
2013          switch(status) {
2014           
2015          case ClpSimplex::basic:
2016            infeasible_->zero(iSequence+addSequence);
2017          case ClpSimplex::isFixed:
2018            break;
2019          case ClpSimplex::isFree:
2020          case ClpSimplex::superBasic:
2021            if (fabs(value)>FREE_ACCEPT*tolerance) {
2022              // we are going to bias towards free (but only if reasonable)
2023              value *= FREE_BIAS;
2024              // store square in list
2025              if (infeas[iSequence+addSequence])
2026                infeas[iSequence+addSequence] = value*value; // already there
2027              else
2028                infeasible_->quickAdd(iSequence+addSequence,value*value);
2029            } else {
2030              infeasible_->zero(iSequence+addSequence);
2031            }
2032            break;
2033          case ClpSimplex::atUpperBound:
2034            if (value>tolerance) {
2035              // store square in list
2036              if (infeas[iSequence+addSequence])
2037                infeas[iSequence+addSequence] = value*value; // already there
2038              else
2039                infeasible_->quickAdd(iSequence+addSequence,value*value);
2040            } else {
2041              infeasible_->zero(iSequence+addSequence);
2042            }
2043            break;
2044          case ClpSimplex::atLowerBound:
2045            if (value<-tolerance) {
2046              // store square in list
2047              if (infeas[iSequence+addSequence])
2048                infeas[iSequence+addSequence] = value*value; // already there
2049              else
2050                infeasible_->quickAdd(iSequence+addSequence,value*value);
2051            } else {
2052              infeasible_->zero(iSequence+addSequence);
2053            }
2054          }
2055        }
2056      } else {
2057        ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
2058        // We can go up OR down
2059        for (j=0;j<number;j++) {
2060          int iSequence = index[j];
2061          double value = reducedCost[iSequence];
2062          value -= updateBy[iSequence];
2063          reducedCost[iSequence] = value;
2064          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
2065         
2066          switch(status) {
2067           
2068          case ClpSimplex::basic:
2069            infeasible_->zero(iSequence+addSequence);
2070          case ClpSimplex::isFixed:
2071            break;
2072          case ClpSimplex::isFree:
2073          case ClpSimplex::superBasic:
2074            if (fabs(value)>FREE_ACCEPT*tolerance) {
2075              // we are going to bias towards free (but only if reasonable)
2076              value *= FREE_BIAS;
2077              // store square in list
2078              if (infeas[iSequence+addSequence])
2079                infeas[iSequence+addSequence] = value*value; // already there
2080              else
2081                infeasible_->quickAdd(iSequence+addSequence,value*value);
2082            } else {
2083              infeasible_->zero(iSequence+addSequence);
2084            }
2085            break;
2086          case ClpSimplex::atUpperBound:
2087            if (value>tolerance) {
2088              // store square in list
2089              if (infeas[iSequence+addSequence])
2090                infeas[iSequence+addSequence] = value*value; // already there
2091              else
2092                infeasible_->quickAdd(iSequence+addSequence,value*value);
2093            } else {
2094              // look other way - change up should be negative
2095              value -= nonLinear->changeUpInCost(iSequence+addSequence);
2096              if (value>-tolerance) {
2097                infeasible_->zero(iSequence+addSequence);
2098              } else {
2099                // store square in list
2100                if (infeas[iSequence+addSequence])
2101                  infeas[iSequence+addSequence] = value*value; // already there
2102                else
2103                  infeasible_->quickAdd(iSequence+addSequence,value*value);
2104              }
2105            }
2106            break;
2107          case ClpSimplex::atLowerBound:
2108            if (value<-tolerance) {
2109              // store square in list
2110              if (infeas[iSequence+addSequence])
2111                infeas[iSequence+addSequence] = value*value; // already there
2112              else
2113                infeasible_->quickAdd(iSequence+addSequence,value*value);
2114            } else {
2115              // look other way - change down should be positive
2116              value -= nonLinear->changeDownInCost(iSequence+addSequence);
2117              if (value<tolerance) {
2118                infeasible_->zero(iSequence+addSequence);
2119              } else {
2120                // store square in list
2121                if (infeas[iSequence+addSequence])
2122                  infeas[iSequence+addSequence] = value*value; // already there
2123                else
2124                  infeasible_->quickAdd(iSequence+addSequence,value*value);
2125              }
2126            }
2127          }
2128        }
2129      }
2130    }
2131    if (anyUpdates==2) {
2132      // we can zero out as will have to get pivot row
2133      updates->clear();
2134      spareColumn1->clear();
2135    }
2136    if (pivotRow>=0) {
2137      // make sure infeasibility on incoming is 0.0
2138      int sequenceIn = model_->sequenceIn();
2139      infeasible_->zero(sequenceIn);
2140    }
2141  }
2142  // make sure outgoing from last iteration okay
2143  int sequenceOut = model_->sequenceOut();
2144  if (sequenceOut>=0) {
2145    ClpSimplex::Status status = model_->getStatus(sequenceOut);
2146    double value = model_->reducedCost(sequenceOut);
2147   
2148    switch(status) {
2149     
2150    case ClpSimplex::basic:
2151    case ClpSimplex::isFixed:
2152      break;
2153    case ClpSimplex::isFree:
2154    case ClpSimplex::superBasic:
2155      if (fabs(value)>FREE_ACCEPT*tolerance) { 
2156        // we are going to bias towards free (but only if reasonable)
2157        value *= FREE_BIAS;
2158        // store square in list
2159        if (infeas[sequenceOut])
2160          infeas[sequenceOut] = value*value; // already there
2161        else
2162          infeasible_->quickAdd(sequenceOut,value*value);
2163      } else {
2164        infeasible_->zero(sequenceOut);
2165      }
2166      break;
2167    case ClpSimplex::atUpperBound:
2168      if (value>tolerance) {
2169        // store square in list
2170        if (infeas[sequenceOut])
2171          infeas[sequenceOut] = value*value; // already there
2172        else
2173          infeasible_->quickAdd(sequenceOut,value*value);
2174      } else {
2175        infeasible_->zero(sequenceOut);
2176      }
2177      break;
2178    case ClpSimplex::atLowerBound:
2179      if (value<-tolerance) {
2180        // store square in list
2181        if (infeas[sequenceOut])
2182          infeas[sequenceOut] = value*value; // already there
2183        else
2184          infeasible_->quickAdd(sequenceOut,value*value);
2185      } else {
2186        infeasible_->zero(sequenceOut);
2187      }
2188    }
2189  }
2190
2191  // If in quadratic re-compute all
2192  if (model_->algorithm()==2) {
2193    for (iSection=0;iSection<2;iSection++) {
2194     
2195      reducedCost=model_->djRegion(iSection);
2196      int addSequence;
2197      int iSequence;
2198     
2199      if (!iSection) {
2200        number = model_->numberRows();
2201        addSequence = model_->numberColumns();
2202      } else {
2203        number = model_->numberColumns();
2204        addSequence = 0;
2205      }
2206     
2207      if (!model_->nonLinearCost()->lookBothWays()) {
2208        for (iSequence=0;iSequence<number;iSequence++) {
2209          double value = reducedCost[iSequence];
2210          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
2211         
2212          switch(status) {
2213           
2214          case ClpSimplex::basic:
2215            infeasible_->zero(iSequence+addSequence);
2216          case ClpSimplex::isFixed:
2217            break;
2218          case ClpSimplex::isFree:
2219          case ClpSimplex::superBasic:
2220            if (fabs(value)>tolerance) {
2221              // we are going to bias towards free (but only if reasonable)
2222              value *= FREE_BIAS;
2223              // store square in list
2224              if (infeas[iSequence+addSequence])
2225                infeas[iSequence+addSequence] = value*value; // already there
2226              else
2227                infeasible_->quickAdd(iSequence+addSequence,value*value);
2228            } else {
2229              infeasible_->zero(iSequence+addSequence);
2230            }
2231            break;
2232          case ClpSimplex::atUpperBound:
2233            if (value>tolerance) {
2234              // store square in list
2235              if (infeas[iSequence+addSequence])
2236                infeas[iSequence+addSequence] = value*value; // already there
2237              else
2238                infeasible_->quickAdd(iSequence+addSequence,value*value);
2239            } else {
2240              infeasible_->zero(iSequence+addSequence);
2241            }
2242            break;
2243          case ClpSimplex::atLowerBound:
2244            if (value<-tolerance) {
2245              // store square in list
2246              if (infeas[iSequence+addSequence])
2247                infeas[iSequence+addSequence] = value*value; // already there
2248              else
2249                infeasible_->quickAdd(iSequence+addSequence,value*value);
2250            } else {
2251              infeasible_->zero(iSequence+addSequence);
2252            }
2253          }
2254        }
2255      } else {
2256        // we can go both ways
2257        ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
2258        for (iSequence=0;iSequence<number;iSequence++) {
2259          double value = reducedCost[iSequence];
2260          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
2261         
2262          switch(status) {
2263           
2264          case ClpSimplex::basic:
2265            infeasible_->zero(iSequence+addSequence);
2266          case ClpSimplex::isFixed:
2267            break;
2268          case ClpSimplex::isFree:
2269          case ClpSimplex::superBasic:
2270            if (fabs(value)>tolerance) {
2271              // we are going to bias towards free (but only if reasonable)
2272              value *= FREE_BIAS;
2273              // store square in list
2274              if (infeas[iSequence+addSequence])
2275                infeas[iSequence+addSequence] = value*value; // already there
2276              else
2277                infeasible_->quickAdd(iSequence+addSequence,value*value);
2278            } else {
2279              infeasible_->zero(iSequence+addSequence);
2280            }
2281            break;
2282          case ClpSimplex::atUpperBound:
2283            if (value>tolerance) {
2284              // store square in list
2285              if (infeas[iSequence+addSequence])
2286                infeas[iSequence+addSequence] = value*value; // already there
2287              else
2288                infeasible_->quickAdd(iSequence+addSequence,value*value);
2289            } else {
2290              // look other way - change up should be negative
2291              value -= nonLinear->changeUpInCost(iSequence+addSequence);
2292              if (value>-tolerance) {
2293                infeasible_->zero(iSequence+addSequence);
2294              } else {
2295                // store square in list
2296                if (infeas[iSequence+addSequence])
2297                  infeas[iSequence+addSequence] = value*value; // already there
2298                else
2299                  infeasible_->quickAdd(iSequence+addSequence,value*value);
2300              }
2301            }
2302            break;
2303          case ClpSimplex::atLowerBound:
2304            if (value<-tolerance) {
2305              // store square in list
2306              if (infeas[iSequence+addSequence])
2307                infeas[iSequence+addSequence] = value*value; // already there
2308              else
2309                infeasible_->quickAdd(iSequence+addSequence,value*value);
2310            } else {
2311              // look other way - change down should be positive
2312              value -= nonLinear->changeDownInCost(iSequence+addSequence);
2313              if (value<tolerance) {
2314                infeasible_->zero(iSequence+addSequence);
2315              } else {
2316                // store square in list
2317                if (infeas[iSequence+addSequence])
2318                  infeas[iSequence+addSequence] = value*value; // already there
2319                else
2320                  infeasible_->quickAdd(iSequence+addSequence,value*value);
2321              }
2322            }
2323          }
2324        }
2325      }
2326    }
2327  }
2328  // See what sort of pricing
2329  int numberWanted=0;
2330  number = infeasible_->getNumElements();
2331  int numberColumns = model_->numberColumns();
2332  if (switchType==5) {
2333    // we can zero out
2334    updates->clear();
2335    spareColumn1->clear();
2336    pivotSequence_=-1;
2337    pivotRow=-1;
2338    // See if to switch
2339    int numberElements = model_->factorization()->numberElements();
2340    int numberRows = model_->numberRows();
2341    // ratio is done on number of columns here
2342    //double ratio = (double) numberElements/(double) numberColumns;
2343    double ratio = (double) numberElements/(double) numberRows;
2344    //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
2345    if (ratio<0.1) {
2346      numberWanted = max(100,number/200);
2347    } else if (ratio<0.3) {
2348      numberWanted = max(500,number/40);
2349    } else if (ratio<0.5) {
2350      numberWanted = max(2000,number/10);
2351      numberWanted = max(numberWanted,numberColumns/30);
2352    } else {
2353      switchType=4;
2354      // initialize
2355      numberSwitched_++;
2356      // Make sure will re-do
2357      delete [] weights_;
2358      weights_=NULL;
2359      saveWeights(model_,4);
2360      printf("switching to devex %d nel ratio %g\n",numberElements,ratio);
2361      updates->clear();
2362    }
2363    if (model_->numberIterations()%1000==0)
2364      printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
2365  }
2366  if(switchType==4) {
2367    // Still in devex mode
2368    int numberElements = model_->factorization()->numberElements();
2369    int numberRows = model_->numberRows();
2370    // ratio is done on number of rows here
2371    double ratio = (double) numberElements/(double) numberRows;
2372    // Go to steepest if lot of iterations?
2373    if (ratio<1.0) {
2374      numberWanted = max(2000,number/20);
2375    } else if (ratio<5.0) {
2376      numberWanted = max(2000,number/10);
2377      numberWanted = max(numberWanted,numberColumns/20);
2378    } else {
2379      // we can zero out
2380      updates->clear();
2381      spareColumn1->clear();
2382      switchType=3;
2383      // initialize
2384      pivotSequence_=-1;
2385      pivotRow=-1;
2386      numberSwitched_++;
2387      // Make sure will re-do
2388      delete [] weights_;
2389      weights_=NULL;
2390      saveWeights(model_,4);
2391      printf("switching to exact %d nel ratio %g\n",numberElements,ratio);
2392      updates->clear();
2393    }
2394    if (model_->numberIterations()%1000==0)
2395      printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
2396  } 
2397  if (switchType<4) {
2398    if (switchType<2 ) {
2399      numberWanted = number+1;
2400    } else if (switchType==2) {
2401      numberWanted = max(2000,number/8);
2402    } else {
2403      int numberElements = model_->factorization()->numberElements();
2404      double ratio = (double) numberElements/(double) model_->numberRows();
2405      if (ratio<1.0) {
2406        numberWanted = max(2000,number/20);
2407      } else if (ratio<5.0) {
2408        numberWanted = max(2000,number/10);
2409        numberWanted = max(numberWanted,numberColumns/20);
2410      } else if (ratio<10.0) {
2411        numberWanted = max(2000,number/8);
2412        numberWanted = max(numberWanted,numberColumns/20);
2413      } else {
2414        ratio = number * (ratio/80.0);
2415        if (ratio>number) {
2416          numberWanted=number+1;
2417        } else {
2418          numberWanted = max(2000,(int) ratio);
2419          numberWanted = max(numberWanted,numberColumns/10);
2420        }
2421      }
2422    }
2423  }
2424  // for weights update we use pivotSequence
2425  pivotRow = pivotSequence_;
2426  // unset in case sub flip
2427  pivotSequence_=-1;
2428  if (pivotRow>=0) {
2429    // make sure infeasibility on incoming is 0.0
2430    const int * pivotVariable = model_->pivotVariable();
2431    int sequenceIn = pivotVariable[pivotRow];
2432    infeasible_->zero(sequenceIn);
2433    // and we can see if reference
2434    double referenceIn=0.0;
2435    if (switchType!=1&&reference(sequenceIn))
2436      referenceIn=1.0;
2437    // save outgoing weight round update
2438    double outgoingWeight=0.0;
2439    if (sequenceOut>=0)
2440      outgoingWeight=weights_[sequenceOut];
2441    // update weights
2442    if (anyUpdates!=1) {
2443      updates->setNumElements(0);
2444      spareColumn1->setNumElements(0);
2445      // might as well set dj to 1
2446      dj = 1.0;
2447      updates->insert(pivotRow,-dj);
2448      model_->factorization()->updateColumnTranspose(spareRow2,updates);
2449      // put row of tableau in rowArray and columnArray
2450      model_->clpMatrix()->transposeTimes(model_,-1.0,
2451                                          updates,spareColumn2,spareColumn1);
2452    }
2453    double * weight;
2454    double * other = alternateWeights_->denseVector();
2455    int numberColumns = model_->numberColumns();
2456    double scaleFactor = -1.0/dj; // as formula is with 1.0
2457    // rows
2458    number = updates->getNumElements();
2459    index = updates->getIndices();
2460    updateBy = updates->denseVector();
2461    weight = weights_+numberColumns;
2462     
2463    if (switchType<4) {
2464      // Exact
2465      // now update weight update array
2466      model_->factorization()->updateColumnTranspose(spareRow2,
2467                                                     alternateWeights_);
2468      for (j=0;j<number;j++) {
2469        int iSequence = index[j];
2470        double thisWeight = weight[iSequence];
2471        // row has -1
2472        double pivot = updateBy[iSequence]*scaleFactor;
2473        updateBy[iSequence]=0.0;
2474        double modification = other[iSequence];
2475        double pivotSquared = pivot * pivot;
2476       
2477        thisWeight += pivotSquared * devex_ + pivot * modification;
2478        if (thisWeight<TRY_NORM) {
2479          if (switchType==1) {
2480            // steepest
2481            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
2482          } else {
2483            // exact
2484            thisWeight = referenceIn*pivotSquared;
2485            if (reference(iSequence+numberColumns))
2486              thisWeight += 1.0;
2487            thisWeight = max(thisWeight,TRY_NORM);
2488          }
2489        }
2490        weight[iSequence] = thisWeight;
2491      }
2492    } else if (switchType==4) {
2493      // Devex
2494      assert (devex_>0.0);
2495      for (j=0;j<number;j++) {
2496        int iSequence = index[j];
2497        double thisWeight = weight[iSequence];
2498        // row has -1
2499        double pivot = updateBy[iSequence]*scaleFactor;
2500        updateBy[iSequence]=0.0;
2501        double value = pivot * pivot*devex_;
2502        if (reference(iSequence+numberColumns))
2503          value += 1.0;
2504        weight[iSequence] = max(0.99*thisWeight,value);
2505      }
2506    }
2507   
2508    // columns
2509    weight = weights_;
2510   
2511    scaleFactor = -scaleFactor;
2512   
2513    number = spareColumn1->getNumElements();
2514    index = spareColumn1->getIndices();
2515    updateBy = spareColumn1->denseVector();
2516    if (switchType<4) {
2517      // Exact
2518      // get subset which have nonzero tableau elements
2519      model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
2520                                                spareColumn1,
2521                                                spareColumn2);
2522      double * updateBy2 = spareColumn2->denseVector();
2523      for (j=0;j<number;j++) {
2524        int iSequence = index[j];
2525        double thisWeight = weight[iSequence];
2526        double pivot = updateBy[iSequence]*scaleFactor;
2527        updateBy[iSequence]=0.0;
2528        double modification = updateBy2[iSequence];
2529        updateBy2[iSequence]=0.0;
2530        double pivotSquared = pivot * pivot;
2531       
2532        thisWeight += pivotSquared * devex_ + pivot * modification;
2533        if (thisWeight<TRY_NORM) {
2534          if (switchType==1) {
2535            // steepest
2536            thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
2537          } else {
2538            // exact
2539            thisWeight = referenceIn*pivotSquared;
2540            if (reference(iSequence))
2541              thisWeight += 1.0;
2542            thisWeight = max(thisWeight,TRY_NORM);
2543          }
2544        }
2545        weight[iSequence] = thisWeight;
2546      }
2547    } else if (switchType==4) {
2548      // Devex
2549      for (j=0;j<number;j++) {
2550        int iSequence = index[j];
2551        double thisWeight = weight[iSequence];
2552        // row has -1
2553        double pivot = updateBy[iSequence]*scaleFactor;
2554        updateBy[iSequence]=0.0;
2555        double value = pivot * pivot*devex_;
2556        if (reference(iSequence))
2557          value += 1.0;
2558        weight[iSequence] = max(0.99*thisWeight,value);
2559      }
2560    }
2561    // restore outgoing weight
2562    if (sequenceOut>=0)
2563      weights_[sequenceOut]=outgoingWeight;
2564    alternateWeights_->clear();
2565    spareColumn2->setNumElements(0);
2566    //#define SOME_DEBUG_1
2567#ifdef SOME_DEBUG_1
2568    // check for accuracy
2569    int iCheck=229;
2570    //printf("weight for iCheck is %g\n",weights_[iCheck]);
2571    int numberRows = model_->numberRows();
2572    //int numberColumns = model_->numberColumns();
2573    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
2574      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
2575          !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
2576        checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
2577    }
2578#endif
2579    updates->setNumElements(0);
2580    spareColumn1->setNumElements(0);
2581  }
2582
2583  // update of duals finished - now do pricing
2584
2585
2586  double bestDj = 1.0e-30;
2587  int bestSequence=-1;
2588
2589  int i,iSequence;
2590  index = infeasible_->getIndices();
2591  number = infeasible_->getNumElements();
2592  if(model_->numberIterations()<model_->lastBadIteration()+200) {
2593    // we can't really trust infeasibilities if there is dual error
2594    double checkTolerance = 1.0e-8;
2595    if (!model_->factorization()->pivots())
2596      checkTolerance = 1.0e-6;
2597    if (model_->largestDualError()>checkTolerance)
2598      tolerance *= model_->largestDualError()/checkTolerance;
2599    // But cap
2600    tolerance = min(1000.0,tolerance);
2601  }
2602#ifdef CLP_DEBUG
2603  if (model_->numberDualInfeasibilities()==1) 
2604    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
2605           model_->largestDualError(),model_,model_->messageHandler(),
2606           number);
2607#endif
2608  // stop last one coming immediately
2609  double saveOutInfeasibility=0.0;
2610  if (sequenceOut>=0) {
2611    saveOutInfeasibility = infeas[sequenceOut];
2612    infeas[sequenceOut]=0.0;
2613  }
2614  tolerance *= tolerance; // as we are using squares
2615
2616  int iPass;
2617  // Setup two passes
2618  int start[4];
2619  start[1]=number;
2620  start[2]=0;
2621  double dstart = ((double) number) * CoinDrand48();
2622  start[0]=(int) dstart;
2623  start[3]=start[0];
2624  //double largestWeight=0.0;
2625  //double smallestWeight=1.0e100;
2626  for (iPass=0;iPass<2;iPass++) {
2627    int end = start[2*iPass+1];
2628    if (switchType<5) {
2629      for (i=start[2*iPass];i<end;i++) {
2630        iSequence = index[i];
2631        double value = infeas[iSequence];
2632        if (value>tolerance) {
2633          double weight = weights_[iSequence];
2634          //weight=1.0;
2635          if (value>bestDj*weight) {
2636            // check flagged variable and correct dj
2637            if (!model_->flagged(iSequence)) {
2638              bestDj=value/weight;
2639              bestSequence = iSequence;
2640            } else {
2641              // just to make sure we don't exit before got something
2642              numberWanted++;
2643            }
2644          }
2645        }
2646        numberWanted--;
2647        if (!numberWanted)
2648          break;
2649      }
2650    } else {
2651      // Dantzig
2652      for (i=start[2*iPass];i<end;i++) {
2653        iSequence = index[i];
2654        double value = infeas[iSequence];
2655        if (value>tolerance) {
2656          if (value>bestDj) {
2657            // check flagged variable and correct dj
2658            if (!model_->flagged(iSequence)) {
2659              bestDj=value;
2660              bestSequence = iSequence;
2661            } else {
2662              // just to make sure we don't exit before got something
2663              numberWanted++;
2664            }
2665          }
2666        }
2667        numberWanted--;
2668        if (!numberWanted)
2669          break;
2670      }
2671    }
2672    if (!numberWanted)
2673      break;
2674  }
2675  if (sequenceOut>=0) {
2676    infeas[sequenceOut]=saveOutInfeasibility;
2677  }
2678  /*if (model_->numberIterations()%100==0)
2679    printf("%d best %g\n",bestSequence,bestDj);*/
2680 
2681#ifdef CLP_DEBUG
2682  if (bestSequence>=0) {
2683    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
2684      assert(model_->reducedCost(bestSequence)<0.0);
2685    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
2686      assert(model_->reducedCost(bestSequence)>0.0);
2687  }
2688#endif
2689  return bestSequence;
2690}
2691/*
2692   1) before factorization
2693   2) after factorization
2694   3) just redo infeasibilities
2695   4) restore weights
2696*/
2697void 
2698ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
2699{
2700  model_ = model;
2701  if (mode_==4) {
2702    if (mode==1&&!weights_) 
2703      numberSwitched_=0; // Reset
2704  }
2705  // alternateWeights_ is defined as indexed but is treated oddly
2706  // at times
2707  int numberRows = model_->numberRows();
2708  int numberColumns = model_->numberColumns();
2709  const int * pivotVariable = model_->pivotVariable();
2710  if (mode==1&&weights_) {
2711    if (pivotSequence_>=0) {
2712      // save pivot order
2713      memcpy(alternateWeights_->getIndices(),pivotVariable,
2714             numberRows*sizeof(int));
2715      // change from pivot row number to sequence number
2716      pivotSequence_=pivotVariable[pivotSequence_];
2717    } else {
2718      alternateWeights_->clear();
2719    }
2720    state_=1;
2721  } else if (mode==2||mode==4||mode==5) {
2722    // restore
2723    if (!weights_||state_==-1||mode==5) {
2724      // initialize weights
2725      delete [] weights_;
2726      delete alternateWeights_;
2727      weights_ = new double[numberRows+numberColumns];
2728      alternateWeights_ = new CoinIndexedVector();
2729      // enough space so can use it for factorization
2730      alternateWeights_->reserve(numberRows+
2731                                 model_->factorization()->maximumPivots());
2732      initializeWeights();
2733      // create saved weights
2734      delete [] savedWeights_;
2735      savedWeights_ = new double[numberRows+numberColumns];
2736      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
2737             sizeof(double));
2738      savedPivotSequence_=-2;
2739      savedSequenceOut_ = -2;
2740     
2741    } else {
2742      if (mode!=4) {
2743        // save
2744        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
2745               sizeof(double));
2746        savedPivotSequence_= pivotSequence_;
2747        savedSequenceOut_ = model_->sequenceOut();
2748      } else {
2749        // restore
2750        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
2751               sizeof(double));
2752        // was - but I think should not be
2753        //pivotSequence_= savedPivotSequence_;
2754        //model_->setSequenceOut(savedSequenceOut_);
2755        pivotSequence_= -1;
2756        model_->setSequenceOut(-1); 
2757      }
2758    }
2759    state_=0;
2760    // set up infeasibilities
2761    if (!infeasible_) {
2762      infeasible_ = new CoinIndexedVector();
2763      infeasible_->reserve(numberColumns+numberRows);
2764    }
2765  }
2766  if (mode>=2&&mode!=5) {
2767    if (mode!=3&&pivotSequence_>=0) {
2768      // restore pivot row
2769      int iRow;
2770      // permute alternateWeights
2771      double * temp = new double[numberRows+numberColumns];
2772      double * work = alternateWeights_->denseVector();
2773      int * oldPivotOrder = alternateWeights_->getIndices();
2774      for (iRow=0;iRow<numberRows;iRow++) {
2775        int iPivot=oldPivotOrder[iRow];
2776        temp[iPivot]=work[iRow];
2777      }
2778      int number=0;
2779      int found=-1;
2780      int * which = oldPivotOrder;
2781      // find pivot row and re-create alternateWeights
2782      for (iRow=0;iRow<numberRows;iRow++) {
2783        int iPivot=pivotVariable[iRow];
2784        if (iPivot==pivotSequence_) 
2785          found=iRow;
2786        work[iRow]=temp[iPivot];
2787        if (work[iRow])
2788          which[number++]=iRow;
2789      }
2790      alternateWeights_->setNumElements(number);
2791#ifdef CLP_DEBUG
2792      // Can happen but I should clean up
2793      assert(found>=0);
2794#endif
2795      pivotSequence_ = found;
2796      delete [] temp;
2797    }
2798    infeasible_->clear();
2799    double tolerance=model_->currentDualTolerance();
2800    int number = model_->numberRows() + model_->numberColumns();
2801    int iSequence;
2802   
2803    double * reducedCost = model_->djRegion();
2804     
2805    if (!model_->nonLinearCost()->lookBothWays()) {
2806      for (iSequence=0;iSequence<number;iSequence++) {
2807        double value = reducedCost[iSequence];
2808        ClpSimplex::Status status = model_->getStatus(iSequence);
2809
2810        switch(status) {
2811         
2812        case ClpSimplex::basic:
2813        case ClpSimplex::isFixed:
2814          break;
2815        case ClpSimplex::isFree:
2816        case ClpSimplex::superBasic:
2817          if (fabs(value)>FREE_ACCEPT*tolerance) { 
2818            // we are going to bias towards free (but only if reasonable)
2819            value *= FREE_BIAS;
2820            // store square in list
2821            infeasible_->quickAdd(iSequence,value*value);
2822          }
2823          break;
2824        case ClpSimplex::atUpperBound:
2825          if (value>tolerance) {
2826            infeasible_->quickAdd(iSequence,value*value);
2827          }
2828          break;
2829        case ClpSimplex::atLowerBound:
2830          if (value<-tolerance) {
2831            infeasible_->quickAdd(iSequence,value*value);
2832          }
2833        }
2834      }
2835    } else {
2836      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
2837      // can go both ways
2838      for (iSequence=0;iSequence<number;iSequence++) {
2839        double value = reducedCost[iSequence];
2840        ClpSimplex::Status status = model_->getStatus(iSequence);
2841
2842        switch(status) {
2843         
2844        case ClpSimplex::basic:
2845        case ClpSimplex::isFixed:
2846          break;
2847        case ClpSimplex::isFree:
2848        case ClpSimplex::superBasic:
2849          if (fabs(value)>FREE_ACCEPT*tolerance) { 
2850            // we are going to bias towards free (but only if reasonable)
2851            value *= FREE_BIAS;
2852            // store square in list
2853            infeasible_->quickAdd(iSequence,value*value);
2854          }
2855          break;
2856        case ClpSimplex::atUpperBound:
2857          if (value>tolerance) {
2858            infeasible_->quickAdd(iSequence,value*value);
2859          } else {
2860            // look other way - change up should be negative
2861            value -= nonLinear->changeUpInCost(iSequence);
2862            if (value<-tolerance) {
2863              // store square in list
2864              infeasible_->quickAdd(iSequence,value*value);
2865            }
2866          }
2867          break;
2868        case ClpSimplex::atLowerBound:
2869          if (value<-tolerance) {
2870            infeasible_->quickAdd(iSequence,value*value);
2871          } else {
2872            // look other way - change down should be positive
2873            value -= nonLinear->changeDownInCost(iSequence);
2874            if (value>tolerance) {
2875              // store square in list
2876              infeasible_->quickAdd(iSequence,value*value);
2877            }
2878          }
2879        }
2880      }
2881    }
2882  }
2883}
2884// Gets rid of last update
2885void 
2886ClpPrimalColumnSteepest::unrollWeights()
2887{
2888  if (mode_==4&&!numberSwitched_)
2889    return;
2890  double * saved = alternateWeights_->denseVector();
2891  int number = alternateWeights_->getNumElements();
2892  int * which = alternateWeights_->getIndices();
2893  int i;
2894  for (i=0;i<number;i++) {
2895    int iRow = which[i];
2896    weights_[iRow]=saved[iRow];
2897    saved[iRow]=0.0;
2898  }
2899  alternateWeights_->setNumElements(0);
2900}
2901 
2902//-------------------------------------------------------------------
2903// Clone
2904//-------------------------------------------------------------------
2905ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
2906{
2907  if (CopyData) {
2908    return new ClpPrimalColumnSteepest(*this);
2909  } else {
2910    return new ClpPrimalColumnSteepest();
2911  }
2912}
2913void
2914ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
2915{
2916  // Local copy of mode so can decide what to do
2917  int switchType = mode_;
2918  if (mode_==4&&numberSwitched_)
2919    switchType=3;
2920  else if (mode_==4)
2921    return;
2922  int number=input->getNumElements();
2923  int * which = input->getIndices();
2924  double * work = input->denseVector();
2925  int newNumber = 0;
2926  int * newWhich = alternateWeights_->getIndices();
2927  double * newWork = alternateWeights_->denseVector();
2928  int i;
2929  int sequenceIn = model_->sequenceIn();
2930  int sequenceOut = model_->sequenceOut();
2931  const int * pivotVariable = model_->pivotVariable();
2932
2933  int pivotRow = model_->pivotRow();
2934  pivotSequence_ = pivotRow;
2935
2936  devex_ =0.0;
2937
2938  if (pivotRow>=0) {
2939    if (switchType==1) {
2940      for (i=0;i<number;i++) {
2941        int iRow = which[i];
2942        devex_ += work[iRow]*work[iRow];
2943        newWork[iRow] = -2.0*work[iRow];
2944      }
2945      newWork[pivotRow] = -2.0*max(devex_,0.0);
2946      devex_+=ADD_ONE;
2947      weights_[sequenceOut]=1.0+ADD_ONE;
2948      memcpy(newWhich,which,number*sizeof(int));
2949      alternateWeights_->setNumElements(number);
2950    } else {
2951      if (mode_!=4||numberSwitched_>1) {
2952        for (i=0;i<number;i++) {
2953          int iRow = which[i];
2954          int iPivot = pivotVariable[iRow];
2955          if (reference(iPivot)) {
2956            devex_ += work[iRow]*work[iRow];
2957            newWork[iRow] = -2.0*work[iRow];
2958            newWhich[newNumber++]=iRow;
2959          }
2960        }
2961        if (!newWork[pivotRow])
2962          newWhich[newNumber++]=pivotRow; // add if not already in
2963        newWork[pivotRow] = -2.0*max(devex_,0.0);
2964      } else {
2965        for (i=0;i<number;i++) {
2966          int iRow = which[i];
2967          int iPivot = pivotVariable[iRow];
2968          if (reference(iPivot)) 
2969            devex_ += work[iRow]*work[iRow];
2970        }
2971      }
2972      if (reference(sequenceIn)) {
2973        devex_+=1.0;
2974      } else {
2975      }
2976      if (reference(sequenceOut)) {
2977        weights_[sequenceOut]=1.0+1.0;
2978      } else {
2979        weights_[sequenceOut]=1.0;
2980      }
2981      alternateWeights_->setNumElements(newNumber);
2982    }
2983  } else {
2984    if (switchType==1) {
2985      for (i=0;i<number;i++) {
2986        int iRow = which[i];
2987        devex_ += work[iRow]*work[iRow];
2988      }
2989      devex_ += ADD_ONE;
2990    } else {
2991      for (i=0;i<number;i++) {
2992        int iRow = which[i];
2993        int iPivot = pivotVariable[iRow];
2994        if (reference(iPivot)) {
2995          devex_ += work[iRow]*work[iRow];
2996        }
2997      }
2998      if (reference(sequenceIn)) 
2999        devex_+=1.0;
3000    }
3001  }
3002  double oldDevex = weights_[sequenceIn];
3003#ifdef CLP_DEBUG
3004  if ((model_->messageHandler()->logLevel()&32))
3005    printf("old weight %g, new %g\n",oldDevex,devex_);
3006#endif
3007  double check = max(devex_,oldDevex)+0.1;
3008  weights_[sequenceIn] = devex_;
3009  double testValue=0.1;
3010  if (mode_==4&&numberSwitched_==1)
3011    testValue = 0.5;
3012  if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
3013#ifdef CLP_DEBUG
3014    if ((model_->messageHandler()->logLevel()&48)==16)
3015      printf("old weight %g, new %g\n",oldDevex,devex_);
3016#endif
3017    //printf("old weight %g, new %g\n",oldDevex,devex_);
3018    testValue=0.99;
3019    if (mode_==1) 
3020      testValue=1.01e1; // make unlikely to do if steepest
3021    else if (mode_==4&&numberSwitched_==1)
3022      testValue=0.9;
3023    double difference = fabs(devex_-oldDevex);
3024    if ( difference > testValue * check ) {
3025      // need to redo
3026      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
3027                                        *model_->messagesPointer())
3028                                          <<oldDevex<<devex_
3029                                          <<CoinMessageEol;
3030      initializeWeights();
3031    }
3032  }
3033  if (pivotRow>=0) {
3034    // set outgoing weight here
3035    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
3036  }
3037}
3038// Checks accuracy - just for debug
3039void 
3040ClpPrimalColumnSteepest::checkAccuracy(int sequence,
3041                                       double relativeTolerance,
3042                                       CoinIndexedVector * rowArray1,
3043                                       CoinIndexedVector * rowArray2)
3044{
3045  if (mode_==4&&!numberSwitched_)
3046    return;
3047  model_->unpack(rowArray1,sequence);
3048  model_->factorization()->updateColumn(rowArray2,rowArray1);
3049  int number=rowArray1->getNumElements();
3050  int * which = rowArray1->getIndices();
3051  double * work = rowArray1->denseVector();
3052  const int * pivotVariable = model_->pivotVariable();
3053 
3054  double devex =0.0;
3055  int i;
3056
3057  if (mode_==1) {
3058    for (i=0;i<number;i++) {
3059      int iRow = which[i];
3060      devex += work[iRow]*work[iRow];
3061      work[iRow]=0.0;
3062    }
3063    devex += ADD_ONE;
3064  } else {
3065    for (i=0;i<number;i++) {
3066      int iRow = which[i];
3067      int iPivot = pivotVariable[iRow];
3068      if (reference(iPivot)) {
3069        devex += work[iRow]*work[iRow];
3070      }
3071      work[iRow]=0.0;
3072    }
3073    if (reference(sequence)) 
3074      devex+=1.0;
3075  }
3076
3077  double oldDevex = weights_[sequence];
3078  double check = max(devex,oldDevex);;
3079  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
3080    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
3081    // update so won't print again
3082    weights_[sequence]=devex;
3083  }
3084  rowArray1->setNumElements(0);
3085}
3086
3087// Initialize weights
3088void 
3089ClpPrimalColumnSteepest::initializeWeights()
3090{
3091  int numberRows = model_->numberRows();
3092  int numberColumns = model_->numberColumns();
3093  int number = numberRows + numberColumns;
3094  int iSequence;
3095  if (mode_!=1) {
3096    // initialize to 1.0
3097    // and set reference framework
3098    if (!reference_) {
3099      int nWords = (number+31)>>5;
3100      reference_ = new unsigned int[nWords];
3101      // tiny overhead to zero out (stops valgrind complaining)
3102      memset(reference_,0,nWords*sizeof(int));
3103    }
3104   
3105    for (iSequence=0;iSequence<number;iSequence++) {
3106      weights_[iSequence]=1.0;
3107      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
3108        setReference(iSequence,false);
3109      } else {
3110        setReference(iSequence,true);
3111      }
3112    }
3113  } else {
3114    CoinIndexedVector * temp = new CoinIndexedVector();
3115    temp->reserve(numberRows+
3116                  model_->factorization()->maximumPivots());
3117    double * array = alternateWeights_->denseVector();
3118    int * which = alternateWeights_->getIndices();
3119     
3120    for (iSequence=0;iSequence<number;iSequence++) {
3121      weights_[iSequence]=1.0+ADD_ONE;
3122      if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
3123          model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
3124        model_->unpack(alternateWeights_,iSequence);
3125        double value=ADD_ONE;
3126        model_->factorization()->updateColumn(temp,alternateWeights_);
3127        int number = alternateWeights_->getNumElements();       int j;
3128        for (j=0;j<number;j++) {
3129          int iRow=which[j];
3130          value+=array[iRow]*array[iRow];
3131          array[iRow]=0.0;
3132        }
3133        alternateWeights_->setNumElements(0);
3134        weights_[iSequence] = value;
3135      }
3136    }
3137    delete temp;
3138  }
3139}
3140// Gets rid of all arrays
3141void 
3142ClpPrimalColumnSteepest::clearArrays()
3143{
3144  delete [] weights_;
3145  weights_=NULL;
3146  delete infeasible_;
3147  infeasible_ = NULL;
3148  delete alternateWeights_;
3149  alternateWeights_ = NULL;
3150  delete [] savedWeights_;
3151  savedWeights_ = NULL;
3152  delete [] reference_;
3153  reference_ = NULL;
3154  pivotSequence_=-1;
3155  state_ = -1;
3156  savedPivotSequence_ = -1;
3157  savedSequenceOut_ = -1; 
3158  devex_ = 0.0;
3159}
3160// Returns true if would not find any column
3161bool 
3162ClpPrimalColumnSteepest::looksOptimal() const
3163{
3164  //**** THIS MUST MATCH the action coding above
3165  double tolerance=model_->currentDualTolerance();
3166  // we can't really trust infeasibilities if there is dual error
3167  // this coding has to mimic coding in checkDualSolution
3168  double error = min(1.0e-3,model_->largestDualError());
3169  // allow tolerance at least slightly bigger than standard
3170  tolerance = tolerance +  error;
3171  if(model_->numberIterations()<model_->lastBadIteration()+200) {
3172    // we can't really trust infeasibilities if there is dual error
3173    double checkTolerance = 1.0e-8;
3174    if (!model_->factorization()->pivots())
3175      checkTolerance = 1.0e-6;
3176    if (model_->largestDualError()>checkTolerance)
3177      tolerance *= model_->largestDualError()/checkTolerance;
3178    // But cap
3179    tolerance = min(1000.0,tolerance);
3180  }
3181  int number = model_->numberRows() + model_->numberColumns();
3182  int iSequence;
3183 
3184  double * reducedCost = model_->djRegion();
3185  int numberInfeasible=0;
3186  if (!model_->nonLinearCost()->lookBothWays()) {
3187    for (iSequence=0;iSequence<number;iSequence++) {
3188      double value = reducedCost[iSequence];
3189      ClpSimplex::Status status = model_->getStatus(iSequence);
3190     
3191      switch(status) {
3192
3193      case ClpSimplex::basic:
3194      case ClpSimplex::isFixed:
3195        break;
3196      case ClpSimplex::isFree:
3197      case ClpSimplex::superBasic:
3198        if (fabs(value)>FREE_ACCEPT*tolerance) 
3199          numberInfeasible++;
3200        break;
3201      case ClpSimplex::atUpperBound:
3202        if (value>tolerance) 
3203          numberInfeasible++;
3204        break;
3205      case ClpSimplex::atLowerBound:
3206        if (value<-tolerance) 
3207          numberInfeasible++;
3208      }
3209    }
3210  } else {
3211    ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
3212    // can go both ways
3213    for (iSequence=0;iSequence<number;iSequence++) {
3214      double value = reducedCost[iSequence];
3215      ClpSimplex::Status status = model_->getStatus(iSequence);
3216     
3217      switch(status) {
3218
3219      case ClpSimplex::basic:
3220      case ClpSimplex::isFixed:
3221        break;
3222      case ClpSimplex::isFree:
3223      case ClpSimplex::superBasic:
3224        if (fabs(value)>FREE_ACCEPT*tolerance) 
3225          numberInfeasible++;
3226        break;
3227      case ClpSimplex::atUpperBound:
3228        if (value>tolerance) {
3229          numberInfeasible++;
3230        } else {
3231          // look other way - change up should be negative
3232          value -= nonLinear->changeUpInCost(iSequence);
3233          if (value<-tolerance) 
3234            numberInfeasible++;
3235        }
3236        break;
3237      case ClpSimplex::atLowerBound:
3238        if (value<-tolerance) {
3239          numberInfeasible++;
3240        } else {
3241          // look other way - change down should be positive
3242          value -= nonLinear->changeDownInCost(iSequence);
3243          if (value>tolerance) 
3244            numberInfeasible++;
3245        }
3246      }
3247    }
3248  }
3249  return numberInfeasible==0;
3250}
Note: See TracBrowser for help on using the repository browser.