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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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