source: branches/pre/ClpPrimalColumnSteepest.cpp @ 1926

Last change on this file since 1926 was 223, checked in by forrest, 17 years ago

Stuff

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