source: trunk/Clp/src/ClpPrimalColumnSteepest.cpp @ 1304

Last change on this file since 1304 was 1197, checked in by forrest, 12 years ago

many changes to try and improve performance

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