source: stable/1.6/Clp/src/ClpPrimalColumnSteepest.cpp @ 1166

Last change on this file since 1166 was 1166, checked in by forrest, 14 years ago

add bandaid

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