source: branches/pre/ClpPrimalColumnSteepest.cpp @ 222

Last change on this file since 222 was 222, checked in by forrest, 16 years ago

Start of mini sprint

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