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

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

for threaded random numbers

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