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

Last change on this file since 846 was 846, checked in by forrest, 15 years ago

fix obscure bug

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