CLP: ClpPrimalColumnSteepest.cpp Source File

00001 // Copyright (C) 2002, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 
00004 #include "CoinPragma.hpp"
00005 
00006 #include "ClpSimplex.hpp"
00007 #include "ClpPrimalColumnSteepest.hpp"
00008 #include "CoinIndexedVector.hpp"
00009 #include "ClpFactorization.hpp"
00010 #include "ClpNonLinearCost.hpp"
00011 #include "ClpMessage.hpp"
00012 #include "CoinHelperFunctions.hpp"
00013 #include <stdio.h>
00014 
00015 
00016 //#############################################################################
00017 // Constructors / Destructor / Assignment
00018 //#############################################################################
00019 
00020 //-------------------------------------------------------------------
00021 // Default Constructor 
00022 //-------------------------------------------------------------------
00023 ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode) 
00024   : ClpPrimalColumnPivot(),
00025     devex_(0.0),
00026     weights_(NULL),
00027     infeasible_(NULL),
00028     alternateWeights_(NULL),
00029     savedWeights_(NULL),
00030     reference_(NULL),
00031     state_(-1),
00032     mode_(mode),
00033     numberSwitched_(0),
00034     pivotSequence_(-1),
00035     savedPivotSequence_(-1),
00036     savedSequenceOut_(-1),
00037     sizeFactorization_(0)
00038 {
00039   type_=2+64*mode;
00040 }
00041 
00042 //-------------------------------------------------------------------
00043 // Copy constructor 
00044 //-------------------------------------------------------------------
00045 ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs) 
00046 : ClpPrimalColumnPivot(rhs)
00047 {  
00048   state_=rhs.state_;
00049   mode_ = rhs.mode_;
00050   numberSwitched_ = rhs.numberSwitched_;
00051   model_ = rhs.model_;
00052   pivotSequence_ = rhs.pivotSequence_;
00053   savedPivotSequence_ = rhs.savedPivotSequence_;
00054   savedSequenceOut_ = rhs.savedSequenceOut_;
00055   sizeFactorization_ = rhs.sizeFactorization_;
00056   devex_ = rhs.devex_;
00057   if (rhs.infeasible_) {
00058     infeasible_= new CoinIndexedVector(rhs.infeasible_);
00059   } else {
00060     infeasible_=NULL;
00061   }
00062   reference_=NULL;
00063   if (rhs.weights_) {
00064     assert(model_);
00065     int number = model_->numberRows()+model_->numberColumns();
00066     weights_= new double[number];
00067     ClpDisjointCopyN(rhs.weights_,number,weights_);
00068     savedWeights_= new double[number];
00069     ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
00070     if (mode_!=1) {
00071       reference_ = new unsigned int[(number+31)>>5];
00072       memcpy(reference_,rhs.reference_,((number+31)>>5)*sizeof(unsigned int));
00073     }
00074   } else {
00075     weights_=NULL;
00076     savedWeights_=NULL;
00077   }
00078   if (rhs.alternateWeights_) {
00079     alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
00080   } else {
00081     alternateWeights_=NULL;
00082   }
00083 }
00084 
00085 //-------------------------------------------------------------------
00086 // Destructor 
00087 //-------------------------------------------------------------------
00088 ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
00089 {
00090   delete [] weights_;
00091   delete infeasible_;
00092   delete alternateWeights_;
00093   delete [] savedWeights_;
00094   delete [] reference_;
00095 }
00096 
00097 //----------------------------------------------------------------
00098 // Assignment operator 
00099 //-------------------------------------------------------------------
00100 ClpPrimalColumnSteepest &
00101 ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
00102 {
00103   if (this != &rhs) {
00104     ClpPrimalColumnPivot::operator=(rhs);
00105     state_=rhs.state_;
00106     mode_ = rhs.mode_;
00107     numberSwitched_ = rhs.numberSwitched_;
00108     model_ = rhs.model_;
00109     pivotSequence_ = rhs.pivotSequence_;
00110     savedPivotSequence_ = rhs.savedPivotSequence_;
00111     savedSequenceOut_ = rhs.savedSequenceOut_;
00112     sizeFactorization_ = rhs.sizeFactorization_;
00113     devex_ = rhs.devex_;
00114     delete [] weights_;
00115     delete [] reference_;
00116     reference_=NULL;
00117     delete infeasible_;
00118     delete alternateWeights_;
00119     delete [] savedWeights_;
00120     savedWeights_ = NULL;
00121     if (rhs.infeasible_!=NULL) {
00122       infeasible_= new CoinIndexedVector(rhs.infeasible_);
00123     } else {
00124       infeasible_=NULL;
00125     }
00126     if (rhs.weights_!=NULL) {
00127       assert(model_);
00128       int number = model_->numberRows()+model_->numberColumns();
00129       weights_= new double[number];
00130       ClpDisjointCopyN(rhs.weights_,number,weights_);
00131       savedWeights_= new double[number];
00132       ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
00133       if (mode_!=1) {
00134         reference_ = new unsigned int[(number+31)>>5];
00135         memcpy(reference_,rhs.reference_,
00136                ((number+31)>>5)*sizeof(unsigned int));
00137       }
00138     } else {
00139       weights_=NULL;
00140     }
00141     if (rhs.alternateWeights_!=NULL) {
00142       alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
00143     } else {
00144       alternateWeights_=NULL;
00145     }
00146   }
00147   return *this;
00148 }
00149 
00150 #define TRY_NORM 1.0e-4
00151 #define ADD_ONE 1.0
00152 // Returns pivot column, -1 if none
00153 /*      The Packed CoinIndexedVector updates has cost updates - for normal LP
00154         that is just +-weight where a feasibility changed.  It also has 
00155         reduced cost from last iteration in pivot row*/
00156 int 
00157 ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
00158                                     CoinIndexedVector * spareRow1,
00159                                     CoinIndexedVector * spareRow2,
00160                                     CoinIndexedVector * spareColumn1,
00161                                     CoinIndexedVector * spareColumn2)
00162 {
00163   assert(model_);
00164   if (model_->nonLinearCost()->lookBothWays()||model_->algorithm()==2) {
00165     // Do old way
00166     updates->expand();
00167     return pivotColumnOldMethod(updates,spareRow1,spareRow2,
00168                                 spareColumn1,spareColumn2);
00169   }
00170   int number=0;
00171   int * index;
00172   double tolerance=model_->currentDualTolerance();
00173   // we can't really trust infeasibilities if there is dual error
00174   // this coding has to mimic coding in checkDualSolution
00175   double error = min(1.0e-3,model_->largestDualError());
00176   // allow tolerance at least slightly bigger than standard
00177   tolerance = tolerance +  error;
00178   int pivotRow = model_->pivotRow();
00179   int anyUpdates;
00180   double * infeas = infeasible_->denseVector();
00181 
00182   // Local copy of mode so can decide what to do
00183   int switchType;
00184   if (mode_==4) 
00185     switchType = 5-numberSwitched_;
00186   else if (mode_>=10)
00187     switchType=3;
00188   else
00189     switchType = mode_;
00190   /* switchType - 
00191      0 - all exact devex
00192      1 - all steepest
00193      2 - some exact devex
00194      3 - auto some exact devex
00195      4 - devex
00196      5 - dantzig
00197      10 - can go to mini-sprint
00198   */
00199   if (updates->getNumElements()>1) {
00200     // would have to have two goes for devex, three for steepest
00201     anyUpdates=2;
00202   } else if (updates->getNumElements()) {
00203     if (updates->getIndices()[0]==pivotRow&&fabs(updates->denseVector()[0])>1.0e-6) {
00204       // reasonable size
00205       anyUpdates=1;
00206       //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
00207       //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
00208     } else {
00209       // too small
00210       anyUpdates=2;
00211     }
00212   } else if (pivotSequence_>=0){
00213     // just after re-factorization
00214     anyUpdates=-1;
00215   } else {
00216     // sub flip - nothing to do
00217     anyUpdates=0;
00218   }
00219   int sequenceOut = model_->sequenceOut();
00220   if (switchType==5) {
00221     // If known matrix then we will do partial pricing
00222     if (model_->clpMatrix()->canDoPartialPricing()) {
00223       pivotSequence_=-1;
00224       pivotRow=-1;
00225       // See if to switch
00226       int numberRows = model_->numberRows();
00227       int numberWanted=0;
00228       int numberColumns = model_->numberColumns();
00229       double ratio = (double) sizeFactorization_/(double) numberRows;
00230       // Number of dual infeasibilities at last invert
00231       int numberDual = model_->numberDualInfeasibilities();
00232       int numberLook = min(numberDual,numberColumns/10);
00233       if (ratio<1.0) {
00234         numberWanted = 100;
00235         numberWanted = max(numberWanted,numberLook/20);
00236       } else if (ratio<3.0) {
00237         numberWanted = 500;
00238         numberWanted = max(numberWanted,numberLook/15);
00239       } else if (ratio<4.0) {
00240         numberWanted = 1000;
00241         numberWanted = max(numberWanted,numberLook/10);
00242       } else {
00243         switchType=4;
00244         // initialize
00245         numberSwitched_++;
00246         // Make sure will re-do
00247         delete [] weights_;
00248         weights_=NULL;
00249         model_->computeDuals(NULL);
00250         saveWeights(model_,4);
00251         anyUpdates=0;
00252         printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
00253       }
00254       if (switchType==5) {
00255         if (model_->numberIterations()%1000==0)
00256           printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
00257         // Update duals and row djs
00258         // Do partial pricing
00259         return partialPricing(updates,spareRow2,
00260                               numberWanted);
00261       }
00262     }
00263   }
00264   if (switchType==5) {
00265     if (anyUpdates>0) {
00266       justDjs(updates,spareRow1,spareRow2,
00267         spareColumn1,spareColumn2);
00268     }
00269   } else if (anyUpdates==1) {
00270     if (switchType<4) {
00271       // exact etc when can use dj 
00272       djsAndSteepest(updates,spareRow1,spareRow2,
00273         spareColumn1,spareColumn2);
00274     } else {
00275       // devex etc when can use dj 
00276       djsAndDevex(updates,spareRow1,spareRow2,
00277         spareColumn1,spareColumn2);
00278     }
00279   } else if (anyUpdates==-1) {
00280     if (switchType<4) {
00281       // exact etc when djs okay 
00282       justSteepest(updates,spareRow1,spareRow2,
00283         spareColumn1,spareColumn2);
00284     } else {
00285       // devex etc when djs okay 
00286       justDevex(updates,spareRow1,spareRow2,
00287         spareColumn1,spareColumn2);
00288     }
00289   } else if (anyUpdates==2) {
00290     if (switchType<4) {
00291       // exact etc when have to use pivot
00292       djsAndSteepest2(updates,spareRow1,spareRow2,
00293         spareColumn1,spareColumn2);
00294     } else {
00295       // devex etc when have to use pivot
00296       djsAndDevex2(updates,spareRow1,spareRow2,
00297         spareColumn1,spareColumn2);
00298     }
00299   } 
00300   alternateWeights_->checkClear();
00301   // make sure outgoing from last iteration okay
00302   if (sequenceOut>=0) {
00303     ClpSimplex::Status status = model_->getStatus(sequenceOut);
00304     double value = model_->reducedCost(sequenceOut);
00305     
00306     switch(status) {
00307       
00308     case ClpSimplex::basic:
00309     case ClpSimplex::isFixed:
00310       break;
00311     case ClpSimplex::isFree:
00312     case ClpSimplex::superBasic:
00313       if (fabs(value)>FREE_ACCEPT*tolerance) { 
00314         // we are going to bias towards free (but only if reasonable)
00315         value *= FREE_BIAS;
00316         // store square in list
00317         if (infeas[sequenceOut])
00318           infeas[sequenceOut] = value*value; // already there
00319         else
00320           infeasible_->quickAdd(sequenceOut,value*value);
00321       } else {
00322         infeasible_->zero(sequenceOut);
00323       }
00324       break;
00325     case ClpSimplex::atUpperBound:
00326       if (value>tolerance) {
00327         // store square in list
00328         if (infeas[sequenceOut])
00329           infeas[sequenceOut] = value*value; // already there
00330         else
00331           infeasible_->quickAdd(sequenceOut,value*value);
00332       } else {
00333         infeasible_->zero(sequenceOut);
00334       }
00335       break;
00336     case ClpSimplex::atLowerBound:
00337       if (value<-tolerance) {
00338         // store square in list
00339         if (infeas[sequenceOut])
00340           infeas[sequenceOut] = value*value; // already there
00341         else
00342           infeasible_->quickAdd(sequenceOut,value*value);
00343       } else {
00344         infeasible_->zero(sequenceOut);
00345       }
00346     }
00347   }
00348   // update of duals finished - now do pricing
00349   // See what sort of pricing
00350   int numberWanted=0;
00351   number = infeasible_->getNumElements();
00352   int numberColumns = model_->numberColumns();
00353   if (switchType==5) {
00354     pivotSequence_=-1;
00355     pivotRow=-1;
00356     // See if to switch
00357     int numberRows = model_->numberRows();
00358     // ratio is done on number of columns here
00359     //double ratio = (double) sizeFactorization_/(double) numberColumns;
00360     double ratio = (double) sizeFactorization_/(double) numberRows;
00361     //double ratio = (double) sizeFactorization_/(double) model_->clpMatrix()->getNumElements();
00362     if (ratio<1.0) {
00363       numberWanted = max(100,number/200);
00364     } else if (ratio<2.0-1.0) {
00365       numberWanted = max(500,number/40);
00366     } else if (ratio<4.0-3.0) {
00367       numberWanted = max(2000,number/10);
00368       numberWanted = max(numberWanted,numberColumns/30);
00369     } else {
00370       switchType=4;
00371       // initialize
00372       numberSwitched_++;
00373       // Make sure will re-do
00374       delete [] weights_;
00375       weights_=NULL;
00376       saveWeights(model_,4);
00377       printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
00378     }
00379     //if (model_->numberIterations()%1000==0)
00380     //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
00381   }
00382   int numberRows = model_->numberRows();
00383   // ratio is done on number of rows here
00384   double ratio = (double) sizeFactorization_/(double) numberRows;
00385   if(switchType==4) {
00386     // Still in devex mode
00387     // Go to steepest if lot of iterations?
00388     if (ratio<5.0) {
00389       numberWanted = max(2000,number/10);
00390       numberWanted = max(numberWanted,numberColumns/20);
00391     } else if (ratio<7.0) {
00392       numberWanted = max(2000,number/5);
00393       numberWanted = max(numberWanted,numberColumns/10);
00394     } else {
00395       // we can zero out
00396       updates->clear();
00397       spareColumn1->clear();
00398       switchType=3;
00399       // initialize
00400       pivotSequence_=-1;
00401       pivotRow=-1;
00402       numberSwitched_++;
00403       // Make sure will re-do
00404       delete [] weights_;
00405       weights_=NULL;
00406       saveWeights(model_,4);
00407       printf("switching to exact %d nel ratio %g\n",sizeFactorization_,ratio);
00408       updates->clear();
00409     }
00410     if (model_->numberIterations()%1000==0)
00411     printf("numels %d ratio %g wanted %d type x\n",sizeFactorization_,ratio,numberWanted);
00412   } 
00413   if (switchType<4) {
00414     if (switchType<2 ) {
00415       numberWanted = number+1;
00416     } else if (switchType==2) {
00417       numberWanted = max(2000,number/8);
00418     } else {
00419       if (ratio<1.0) {
00420         numberWanted = max(2000,number/20);
00421       } else if (ratio<5.0) {
00422         numberWanted = max(2000,number/10);
00423         numberWanted = max(numberWanted,numberColumns/40);
00424       } else if (ratio<10.0) {
00425         numberWanted = max(2000,number/8);
00426         numberWanted = max(numberWanted,numberColumns/20);
00427       } else {
00428         ratio = number * (ratio/80.0);
00429         if (ratio>number) {
00430           numberWanted=number+1;
00431         } else {
00432           numberWanted = max(2000,(int) ratio);
00433           numberWanted = max(numberWanted,numberColumns/10);
00434         }
00435       }
00436     }
00437     //if (model_->numberIterations()%1000==0)
00438     //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
00439     //switchType);
00440   }
00441 
00442 
00443   double bestDj = 1.0e-30;
00444   int bestSequence=-1;
00445 
00446   int i,iSequence;
00447   index = infeasible_->getIndices();
00448   number = infeasible_->getNumElements();
00449   // Re-sort infeasible every 100 pivots
00450   // Not a good idea
00451   if (0&&model_->factorization()->pivots()>0&&
00452       (model_->factorization()->pivots()%100)==0) {
00453     int nLook = model_->numberRows()+numberColumns;
00454     number=0;
00455     for (i=0;i<nLook;i++) {
00456       if (infeas[i]) { 
00457         if (fabs(infeas[i])>COIN_INDEXED_TINY_ELEMENT) 
00458           index[number++]=i;
00459         else
00460           infeas[i]=0.0;
00461       }
00462     }
00463     infeasible_->setNumElements(number);
00464   }
00465   if(model_->numberIterations()<model_->lastBadIteration()+200) {
00466     // we can't really trust infeasibilities if there is dual error
00467     double checkTolerance = 1.0e-8;
00468     if (!model_->factorization()->pivots())
00469       checkTolerance = 1.0e-6;
00470     if (model_->largestDualError()>checkTolerance)
00471       tolerance *= model_->largestDualError()/checkTolerance;
00472     // But cap
00473     tolerance = min(1000.0,tolerance);
00474   }
00475 #ifdef CLP_DEBUG
00476   if (model_->numberDualInfeasibilities()==1) 
00477     printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
00478            model_->largestDualError(),model_,model_->messageHandler(),
00479            number);
00480 #endif
00481   // stop last one coming immediately
00482   double saveOutInfeasibility=0.0;
00483   if (sequenceOut>=0) {
00484     saveOutInfeasibility = infeas[sequenceOut];
00485     infeas[sequenceOut]=0.0;
00486   }
00487   if (model_->factorization()->pivots()&&model_->numberPrimalInfeasibilities())
00488     tolerance = max(tolerance,1.0e-10*model_->infeasibilityCost());
00489   tolerance *= tolerance; // as we are using squares
00490 
00491   int iPass;
00492   // Setup two passes
00493   int start[4];
00494   start[1]=number;
00495   start[2]=0;
00496   double dstart = ((double) number) * CoinDrand48();
00497   start[0]=(int) dstart;
00498   start[3]=start[0];
00499   //double largestWeight=0.0;
00500   //double smallestWeight=1.0e100;
00501   for (iPass=0;iPass<2;iPass++) {
00502     int end = start[2*iPass+1];
00503     if (switchType<5) {
00504       for (i=start[2*iPass];i<end;i++) {
00505         iSequence = index[i];
00506         double value = infeas[iSequence];
00507         if (value>tolerance) {
00508           double weight = weights_[iSequence];
00509           //weight=1.0;
00510           if (value>bestDj*weight) {
00511             // check flagged variable and correct dj
00512             if (!model_->flagged(iSequence)) {
00513               bestDj=value/weight;
00514               bestSequence = iSequence;
00515             } else {
00516               // just to make sure we don't exit before got something
00517               numberWanted++;
00518             }
00519           }
00520         }
00521         numberWanted--;
00522         if (!numberWanted)
00523           break;
00524       }
00525     } else {
00526       // Dantzig
00527       for (i=start[2*iPass];i<end;i++) {
00528         iSequence = index[i];
00529         double value = infeas[iSequence];
00530         if (value>tolerance) {
00531           if (value>bestDj) {
00532             // check flagged variable and correct dj
00533             if (!model_->flagged(iSequence)) {
00534               bestDj=value;
00535               bestSequence = iSequence;
00536             } else {
00537               // just to make sure we don't exit before got something
00538               numberWanted++;
00539             }
00540           }
00541         }
00542         numberWanted--;
00543         if (!numberWanted)
00544           break;
00545       }
00546     }
00547     if (!numberWanted)
00548       break;
00549   }
00550   if (sequenceOut>=0) {
00551     infeas[sequenceOut]=saveOutInfeasibility;
00552   }
00553   /*if (model_->numberIterations()%100==0)
00554     printf("%d best %g\n",bestSequence,bestDj);*/
00555   
00556 #ifndef NDEBUG
00557   if (bestSequence>=0) {
00558     if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
00559       assert(model_->reducedCost(bestSequence)<0.0);
00560     if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
00561       assert(model_->reducedCost(bestSequence)>0.0);
00562   }
00563 #endif
00564   return bestSequence;
00565 }
00566 // Just update djs
00567 void 
00568 ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
00569                CoinIndexedVector * spareRow1,
00570                CoinIndexedVector * spareRow2,
00571                CoinIndexedVector * spareColumn1,
00572                CoinIndexedVector * spareColumn2)
00573 {
00574   int iSection,j;
00575   int number=0;
00576   int * index;
00577   double * updateBy;
00578   double * reducedCost;
00579   double tolerance=model_->currentDualTolerance();
00580   // we can't really trust infeasibilities if there is dual error
00581   // this coding has to mimic coding in checkDualSolution
00582   double error = min(1.0e-3,model_->largestDualError());
00583   // allow tolerance at least slightly bigger than standard
00584   tolerance = tolerance +  error;
00585   int pivotRow = model_->pivotRow();
00586   double * infeas = infeasible_->denseVector();
00587   //updates->scanAndPack();
00588   model_->factorization()->updateColumnTranspose(spareRow2,updates);
00589   
00590   // put row of tableau in rowArray and columnArray (packed mode)
00591   model_->clpMatrix()->transposeTimes(model_,-1.0,
00592                                       updates,spareColumn2,spareColumn1);
00593   // normal
00594   for (iSection=0;iSection<2;iSection++) {
00595     
00596     reducedCost=model_->djRegion(iSection);
00597     int addSequence;
00598     
00599     if (!iSection) {
00600       number = updates->getNumElements();
00601       index = updates->getIndices();
00602       updateBy = updates->denseVector();
00603       addSequence = model_->numberColumns();
00604     } else {
00605       number = spareColumn1->getNumElements();
00606       index = spareColumn1->getIndices();
00607       updateBy = spareColumn1->denseVector();
00608       addSequence = 0;
00609     }
00610     
00611     for (j=0;j<number;j++) {
00612       int iSequence = index[j];
00613       double value = reducedCost[iSequence];
00614       value -= updateBy[j];
00615       updateBy[j]=0.0;
00616       reducedCost[iSequence] = value;
00617       ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
00618       
00619       switch(status) {
00620         
00621       case ClpSimplex::basic:
00622         infeasible_->zero(iSequence+addSequence);
00623       case ClpSimplex::isFixed:
00624         break;
00625       case ClpSimplex::isFree:
00626       case ClpSimplex::superBasic:
00627         if (fabs(value)>FREE_ACCEPT*tolerance) {
00628           // we are going to bias towards free (but only if reasonable)
00629           value *= FREE_BIAS;
00630           // store square in list
00631           if (infeas[iSequence+addSequence])
00632             infeas[iSequence+addSequence] = value*value; // already there
00633           else
00634             infeasible_->quickAdd(iSequence+addSequence,value*value);
00635         } else {
00636           infeasible_->zero(iSequence+addSequence);
00637         }
00638         break;
00639       case ClpSimplex::atUpperBound:
00640         if (value>tolerance) {
00641           // store square in list
00642           if (infeas[iSequence+addSequence])
00643             infeas[iSequence+addSequence] = value*value; // already there
00644           else
00645             infeasible_->quickAdd(iSequence+addSequence,value*value);
00646         } else {
00647           infeasible_->zero(iSequence+addSequence);
00648         }
00649         break;
00650       case ClpSimplex::atLowerBound:
00651         if (value<-tolerance) {
00652           // store square in list
00653           if (infeas[iSequence+addSequence])
00654             infeas[iSequence+addSequence] = value*value; // already there
00655           else
00656             infeasible_->quickAdd(iSequence+addSequence,value*value);
00657         } else {
00658           infeasible_->zero(iSequence+addSequence);
00659         }
00660       }
00661     }
00662   }
00663   updates->setNumElements(0);
00664   spareColumn1->setNumElements(0);
00665   if (pivotRow>=0) {
00666     // make sure infeasibility on incoming is 0.0
00667     int sequenceIn = model_->sequenceIn();
00668     infeasible_->zero(sequenceIn);
00669   }
00670 }
00671 // Update djs, weights for Devex
00672 void 
00673 ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
00674                CoinIndexedVector * spareRow1,
00675                CoinIndexedVector * spareRow2,
00676                CoinIndexedVector * spareColumn1,
00677                CoinIndexedVector * spareColumn2)
00678 {
00679   int j;
00680   int number=0;
00681   int * index;
00682   double * updateBy;
00683   double * reducedCost;
00684   double tolerance=model_->currentDualTolerance();
00685   // we can't really trust infeasibilities if there is dual error
00686   // this coding has to mimic coding in checkDualSolution
00687   double error = min(1.0e-3,model_->largestDualError());
00688   // allow tolerance at least slightly bigger than standard
00689   tolerance = tolerance +  error;
00690   // for weights update we use pivotSequence
00691   // unset in case sub flip
00692   assert (pivotSequence_>=0);
00693   pivotSequence_=-1;
00694   double * infeas = infeasible_->denseVector();
00695   //updates->scanAndPack();
00696   model_->factorization()->updateColumnTranspose(spareRow2,updates);
00697   // and we can see if reference
00698   double referenceIn=0.0;
00699   int sequenceIn = model_->sequenceIn();
00700   if (mode_!=1&&reference(sequenceIn))
00701     referenceIn=1.0;
00702   // save outgoing weight round update
00703   double outgoingWeight=0.0;
00704   int sequenceOut = model_->sequenceOut();
00705   if (sequenceOut>=0)
00706     outgoingWeight=weights_[sequenceOut];
00707     
00708   double scaleFactor = 1.0/updates->denseVector()[0]; // as formula is with 1.0
00709   // put row of tableau in rowArray and columnArray (packed mode)
00710   model_->clpMatrix()->transposeTimes(model_,-1.0,
00711                                       updates,spareColumn2,spareColumn1);
00712   // update weights
00713   double * weight;
00714   int numberColumns = model_->numberColumns();
00715   // rows
00716   reducedCost=model_->djRegion(0);
00717   int addSequence=model_->numberColumns();;
00718     
00719   number = updates->getNumElements();
00720   index = updates->getIndices();
00721   updateBy = updates->denseVector();
00722   weight = weights_+numberColumns;
00723   // Devex
00724   for (j=0;j<number;j++) {
00725     double thisWeight;
00726     double pivot;
00727     double value3;
00728     int iSequence = index[j];
00729     double value = reducedCost[iSequence];
00730     double value2 = updateBy[j];
00731     updateBy[j]=0.0;
00732     value -= value2;
00733     reducedCost[iSequence] = value;
00734     ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
00735     
00736     switch(status) {
00737       
00738     case ClpSimplex::basic:
00739       infeasible_->zero(iSequence+addSequence);
00740     case ClpSimplex::isFixed:
00741       break;
00742     case ClpSimplex::isFree:
00743     case ClpSimplex::superBasic:
00744       thisWeight = weight[iSequence];
00745       // row has -1 
00746       pivot = value2*scaleFactor;
00747       value3 = pivot * pivot*devex_;
00748       if (reference(iSequence+numberColumns))
00749         value3 += 1.0;
00750       weight[iSequence] = max(0.99*thisWeight,value3);
00751       if (fabs(value)>FREE_ACCEPT*tolerance) {
00752         // we are going to bias towards free (but only if reasonable)
00753         value *= FREE_BIAS;
00754         // store square in list
00755         if (infeas[iSequence+addSequence])
00756           infeas[iSequence+addSequence] = value*value; // already there
00757         else
00758           infeasible_->quickAdd(iSequence+addSequence,value*value);
00759       } else {
00760         infeasible_->zero(iSequence+addSequence);
00761       }
00762       break;
00763     case ClpSimplex::atUpperBound:
00764       thisWeight = weight[iSequence];
00765       // row has -1 
00766       pivot = value2*scaleFactor;
00767       value3 = pivot * pivot*devex_;
00768       if (reference(iSequence+numberColumns))
00769         value3 += 1.0;
00770       weight[iSequence] = max(0.99*thisWeight,value3);
00771       if (value>tolerance) {
00772         // store square in list
00773         if (infeas[iSequence+addSequence])
00774           infeas[iSequence+addSequence] = value*value; // already there
00775         else
00776           infeasible_->quickAdd(iSequence+addSequence,value*value);
00777       } else {
00778         infeasible_->zero(iSequence+addSequence);
00779       }
00780       break;
00781     case ClpSimplex::atLowerBound:
00782       thisWeight = weight[iSequence];
00783       // row has -1 
00784       pivot = value2*scaleFactor;
00785       value3 = pivot * pivot*devex_;
00786       if (reference(iSequence+numberColumns))
00787         value3 += 1.0;
00788       weight[iSequence] = max(0.99*thisWeight,value3);
00789       if (value<-tolerance) {
00790         // store square in list
00791         if (infeas[iSequence+addSequence])
00792           infeas[iSequence+addSequence] = value*value; // already there
00793         else
00794           infeasible_->quickAdd(iSequence+addSequence,value*value);
00795       } else {
00796         infeasible_->zero(iSequence+addSequence);
00797       }
00798     }
00799   }
00800   
00801   // columns
00802   weight = weights_;
00803   
00804   scaleFactor = -scaleFactor;
00805   reducedCost=model_->djRegion(1);
00806   number = spareColumn1->getNumElements();
00807   index = spareColumn1->getIndices();
00808   updateBy = spareColumn1->denseVector();
00809 
00810   // Devex
00811   
00812   for (j=0;j<number;j++) {
00813     double thisWeight;
00814     double pivot;
00815     double value3;
00816     int iSequence = index[j];
00817     double value = reducedCost[iSequence];
00818     double value2 = updateBy[j];
00819     value -= value2;
00820     updateBy[j]=0.0;
00821     reducedCost[iSequence] = value;
00822     ClpSimplex::Status status = model_->getStatus(iSequence);
00823     
00824     switch(status) {
00825       
00826     case ClpSimplex::basic:
00827       infeasible_->zero(iSequence);
00828     case ClpSimplex::isFixed:
00829       break;
00830     case ClpSimplex::isFree:
00831     case ClpSimplex::superBasic:
00832       thisWeight = weight[iSequence];
00833       // row has -1 
00834       pivot = value2*scaleFactor;
00835       value3 = pivot * pivot*devex_;
00836       if (reference(iSequence))
00837         value3 += 1.0;
00838       weight[iSequence] = max(0.99*thisWeight,value3);
00839       if (fabs(value)>FREE_ACCEPT*tolerance) {
00840         // we are going to bias towards free (but only if reasonable)
00841         value *= FREE_BIAS;
00842         // store square in list
00843         if (infeas[iSequence])
00844           infeas[iSequence] = value*value; // already there
00845         else
00846           infeasible_->quickAdd(iSequence,value*value);
00847       } else {
00848         infeasible_->zero(iSequence);
00849       }
00850       break;
00851     case ClpSimplex::atUpperBound:
00852       thisWeight = weight[iSequence];
00853       // row has -1 
00854       pivot = value2*scaleFactor;
00855       value3 = pivot * pivot*devex_;
00856       if (reference(iSequence))
00857         value3 += 1.0;
00858       weight[iSequence] = max(0.99*thisWeight,value3);
00859       if (value>tolerance) {
00860         // store square in list
00861         if (infeas[iSequence])
00862           infeas[iSequence] = value*value; // already there
00863         else
00864           infeasible_->quickAdd(iSequence,value*value);
00865       } else {
00866         infeasible_->zero(iSequence);
00867       }
00868       break;
00869     case ClpSimplex::atLowerBound:
00870       thisWeight = weight[iSequence];
00871       // row has -1 
00872       pivot = value2*scaleFactor;
00873       value3 = pivot * pivot*devex_;
00874       if (reference(iSequence))
00875         value3 += 1.0;
00876       weight[iSequence] = max(0.99*thisWeight,value3);
00877       if (value<-tolerance) {
00878         // store square in list
00879         if (infeas[iSequence])
00880           infeas[iSequence] = value*value; // already there
00881         else
00882           infeasible_->quickAdd(iSequence,value*value);
00883       } else {
00884         infeasible_->zero(iSequence);
00885       }
00886     }
00887   }
00888   // restore outgoing weight
00889   if (sequenceOut>=0)
00890     weights_[sequenceOut]=outgoingWeight;
00891   // make sure infeasibility on incoming is 0.0
00892   infeasible_->zero(sequenceIn);
00893   spareRow2->setNumElements(0);
00894   //#define SOME_DEBUG_1
00895 #ifdef SOME_DEBUG_1
00896   // check for accuracy
00897   int iCheck=229;
00898   //printf("weight for iCheck is %g\n",weights_[iCheck]);
00899   int numberRows = model_->numberRows();
00900   //int numberColumns = model_->numberColumns();
00901   for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
00902     if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
00903         !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
00904       checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
00905   }
00906 #endif
00907   updates->setNumElements(0);
00908   spareColumn1->setNumElements(0);
00909 }
00910 // Update djs, weights for Steepest
00911 void 
00912 ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates,
00913                CoinIndexedVector * spareRow1,
00914                CoinIndexedVector * spareRow2,
00915                CoinIndexedVector * spareColumn1,
00916                CoinIndexedVector * spareColumn2)
00917 {
00918   int j;
00919   int number=0;
00920   int * index;
00921   double * updateBy;
00922   double * reducedCost;
00923   double tolerance=model_->currentDualTolerance();
00924   // we can't really trust infeasibilities if there is dual error
00925   // this coding has to mimic coding in checkDualSolution
00926   double error = min(1.0e-3,model_->largestDualError());
00927   // allow tolerance at least slightly bigger than standard
00928   tolerance = tolerance +  error;
00929   // for weights update we use pivotSequence
00930   // unset in case sub flip
00931   assert (pivotSequence_>=0);
00932   pivotSequence_=-1;
00933   double * infeas = infeasible_->denseVector();
00934   double scaleFactor = 1.0/updates->denseVector()[0]; // as formula is with 1.0
00935   //updates->scanAndPack();
00936   model_->factorization()->updateColumnTranspose(spareRow2,updates);
00937   //alternateWeights_->scanAndPack();
00938   model_->factorization()->updateColumnTranspose(spareRow2,
00939                                                  alternateWeights_);
00940   // and we can see if reference
00941   double referenceIn=0.0;
00942   int sequenceIn = model_->sequenceIn();
00943   if (mode_!=1&&reference(sequenceIn))
00944     referenceIn=1.0;
00945   // save outgoing weight round update
00946   double outgoingWeight=0.0;
00947   int sequenceOut = model_->sequenceOut();
00948   if (sequenceOut>=0)
00949     outgoingWeight=weights_[sequenceOut];
00950     
00951   // put row of tableau in rowArray and columnArray (packed)
00952   model_->clpMatrix()->transposeTimes(model_,-1.0,
00953                                       updates,spareColumn2,spareColumn1);
00954   // get subset which have nonzero tableau elements
00955   // Luckily spareRow2 is long enough (rowArray_[3])
00956   model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
00957                                               spareColumn1,
00958                                               spareRow2);
00959   // update weights
00960   double * weight;
00961   double * other = alternateWeights_->denseVector();
00962   int numberColumns = model_->numberColumns();
00963   // rows
00964   reducedCost=model_->djRegion(0);
00965   int addSequence=model_->numberColumns();;
00966     
00967   number = updates->getNumElements();
00968   index = updates->getIndices();
00969   updateBy = updates->denseVector();
00970   weight = weights_+numberColumns;
00971   
00972   for (j=0;j<number;j++) {
00973     double thisWeight;
00974     double pivot;
00975     double modification;
00976     double pivotSquared;
00977     int iSequence = index[j];
00978     double value = reducedCost[iSequence];
00979     double value2 = updateBy[j];
00980     value -= value2;
00981     reducedCost[iSequence] = value;
00982     ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
00983     modification = other[iSequence];
00984     updateBy[j]=0.0;
00985     
00986     switch(status) {
00987       
00988     case ClpSimplex::basic:
00989       infeasible_->zero(iSequence+addSequence);
00990     case ClpSimplex::isFixed:
00991       break;
00992     case ClpSimplex::isFree:
00993     case ClpSimplex::superBasic:
00994       thisWeight = weight[iSequence];
00995       // row has -1 
00996       pivot = value2*scaleFactor;
00997       pivotSquared = pivot * pivot;
00998       
00999       thisWeight += pivotSquared * devex_ + pivot * modification;
01000       if (thisWeight<TRY_NORM) {
01001         if (mode_==1) {
01002           // steepest
01003           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01004         } else {
01005           // exact
01006           thisWeight = referenceIn*pivotSquared;
01007           if (reference(iSequence+numberColumns))
01008             thisWeight += 1.0;
01009           thisWeight = max(thisWeight,TRY_NORM);
01010         }
01011       }
01012       weight[iSequence] = thisWeight;
01013       if (fabs(value)>FREE_ACCEPT*tolerance) {
01014         // we are going to bias towards free (but only if reasonable)
01015         value *= FREE_BIAS;
01016         // store square in list
01017         if (infeas[iSequence+addSequence])
01018           infeas[iSequence+addSequence] = value*value; // already there
01019         else
01020           infeasible_->quickAdd(iSequence+addSequence,value*value);
01021       } else {
01022         infeasible_->zero(iSequence+addSequence);
01023       }
01024       break;
01025     case ClpSimplex::atUpperBound:
01026       thisWeight = weight[iSequence];
01027       // row has -1 
01028       pivot = value2*scaleFactor;
01029       pivotSquared = pivot * pivot;
01030       
01031       thisWeight += pivotSquared * devex_ + pivot * modification;
01032       if (thisWeight<TRY_NORM) {
01033         if (mode_==1) {
01034           // steepest
01035           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01036         } else {
01037           // exact
01038           thisWeight = referenceIn*pivotSquared;
01039           if (reference(iSequence+numberColumns))
01040             thisWeight += 1.0;
01041           thisWeight = max(thisWeight,TRY_NORM);
01042         }
01043       }
01044       weight[iSequence] = thisWeight;
01045       if (value>tolerance) {
01046         // store square in list
01047         if (infeas[iSequence+addSequence])
01048           infeas[iSequence+addSequence] = value*value; // already there
01049         else
01050           infeasible_->quickAdd(iSequence+addSequence,value*value);
01051       } else {
01052         infeasible_->zero(iSequence+addSequence);
01053       }
01054       break;
01055     case ClpSimplex::atLowerBound:
01056       thisWeight = weight[iSequence];
01057       // row has -1 
01058       pivot = value2*scaleFactor;
01059       pivotSquared = pivot * pivot;
01060       
01061       thisWeight += pivotSquared * devex_ + pivot * modification;
01062       if (thisWeight<TRY_NORM) {
01063         if (mode_==1) {
01064           // steepest
01065           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01066         } else {
01067           // exact
01068           thisWeight = referenceIn*pivotSquared;
01069           if (reference(iSequence+numberColumns))
01070             thisWeight += 1.0;
01071           thisWeight = max(thisWeight,TRY_NORM);
01072         }
01073       }
01074       weight[iSequence] = thisWeight;
01075       if (value<-tolerance) {
01076         // store square in list
01077         if (infeas[iSequence+addSequence])
01078           infeas[iSequence+addSequence] = value*value; // already there
01079         else
01080           infeasible_->quickAdd(iSequence+addSequence,value*value);
01081       } else {
01082         infeasible_->zero(iSequence+addSequence);
01083       }
01084     }
01085   }
01086   alternateWeights_->clear();
01087   // columns
01088   weight = weights_;
01089   
01090   scaleFactor = -scaleFactor;
01091   reducedCost=model_->djRegion(1);
01092   number = spareColumn1->getNumElements();
01093   index = spareColumn1->getIndices();
01094   updateBy = spareColumn1->denseVector();
01095   double * updateBy2 = spareRow2->denseVector();
01096 
01097   for (j=0;j<number;j++) {
01098     double thisWeight;
01099     double pivot;
01100     double modification;
01101     double pivotSquared;
01102     int iSequence = index[j];
01103     double value = reducedCost[iSequence];
01104     double value2 = updateBy[j];
01105     updateBy[j]=0.0;
01106     value -= value2;
01107     reducedCost[iSequence] = value;
01108     ClpSimplex::Status status = model_->getStatus(iSequence);
01109     
01110     switch(status) {
01111       
01112     case ClpSimplex::basic:
01113       infeasible_->zero(iSequence);
01114       updateBy2[j]=0.0;
01115     case ClpSimplex::isFixed:
01116       updateBy2[j]=0.0;
01117       break;
01118     case ClpSimplex::isFree:
01119     case ClpSimplex::superBasic:
01120       thisWeight = weight[iSequence];
01121       pivot = value2*scaleFactor;
01122       modification = updateBy2[j];
01123       updateBy2[j]=0.0;
01124       pivotSquared = pivot * pivot;
01125       
01126       thisWeight += pivotSquared * devex_ + pivot * modification;
01127       if (thisWeight<TRY_NORM) {
01128         if (mode_==1) {
01129           // steepest
01130           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01131         } else {
01132           // exact
01133           thisWeight = referenceIn*pivotSquared;
01134           if (reference(iSequence))
01135             thisWeight += 1.0;
01136           thisWeight = max(thisWeight,TRY_NORM);
01137         }
01138       }
01139       weight[iSequence] = thisWeight;
01140       if (fabs(value)>FREE_ACCEPT*tolerance) {
01141         // we are going to bias towards free (but only if reasonable)
01142         value *= FREE_BIAS;
01143         // store square in list
01144         if (infeas[iSequence])
01145           infeas[iSequence] = value*value; // already there
01146         else
01147           infeasible_->quickAdd(iSequence,value*value);
01148       } else {
01149         infeasible_->zero(iSequence);
01150       }
01151       break;
01152     case ClpSimplex::atUpperBound:
01153       thisWeight = weight[iSequence];
01154       pivot = value2*scaleFactor;
01155       modification = updateBy2[j];
01156       updateBy2[j]=0.0;
01157       pivotSquared = pivot * pivot;
01158       
01159       thisWeight += pivotSquared * devex_ + pivot * modification;
01160       if (thisWeight<TRY_NORM) {
01161         if (mode_==1) {
01162           // steepest
01163           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01164         } else {
01165           // exact
01166           thisWeight = referenceIn*pivotSquared;
01167           if (reference(iSequence))
01168             thisWeight += 1.0;
01169           thisWeight = max(thisWeight,TRY_NORM);
01170         }
01171       }
01172       weight[iSequence] = thisWeight;
01173       if (value>tolerance) {
01174         // store square in list
01175         if (infeas[iSequence])
01176           infeas[iSequence] = value*value; // already there
01177         else
01178           infeasible_->quickAdd(iSequence,value*value);
01179       } else {
01180         infeasible_->zero(iSequence);
01181       }
01182       break;
01183     case ClpSimplex::atLowerBound:
01184       thisWeight = weight[iSequence];
01185       pivot = value2*scaleFactor;
01186       modification = updateBy2[j];
01187       updateBy2[j]=0.0;
01188       pivotSquared = pivot * pivot;
01189       
01190       thisWeight += pivotSquared * devex_ + pivot * modification;
01191       if (thisWeight<TRY_NORM) {
01192         if (mode_==1) {
01193           // steepest
01194           thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01195         } else {
01196           // exact
01197           thisWeight = referenceIn*pivotSquared;
01198           if (reference(iSequence))
01199             thisWeight += 1.0;
01200           thisWeight = max(thisWeight,TRY_NORM);
01201         }
01202       }
01203       weight[iSequence] = thisWeight;
01204       if (value<-tolerance) {
01205         // store square in list
01206         if (infeas[iSequence])
01207           infeas[iSequence] = value*value; // already there
01208         else
01209           infeasible_->quickAdd(iSequence,value*value);
01210       } else {
01211         infeasible_->zero(iSequence);
01212       }
01213     }
01214   }
01215   // restore outgoing weight
01216   if (sequenceOut>=0)
01217     weights_[sequenceOut]=outgoingWeight;
01218   // make sure infeasibility on incoming is 0.0
01219   infeasible_->zero(sequenceIn);
01220   spareRow2->setNumElements(0);
01221   //#define SOME_DEBUG_1
01222 #ifdef SOME_DEBUG_1
01223   // check for accuracy
01224   int iCheck=229;
01225   //printf("weight for iCheck is %g\n",weights_[iCheck]);
01226   int numberRows = model_->numberRows();
01227   //int numberColumns = model_->numberColumns();
01228   for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
01229     if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
01230         !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
01231       checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
01232   }
01233 #endif
01234   updates->setNumElements(0);
01235   spareColumn1->setNumElements(0);
01236 }
01237 // Update djs, weights for Devex
01238 void 
01239 ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
01240                CoinIndexedVector * spareRow1,
01241                CoinIndexedVector * spareRow2,
01242                CoinIndexedVector * spareColumn1,
01243                CoinIndexedVector * spareColumn2)
01244 {
01245   int iSection,j;
01246   int number=0;
01247   int * index;
01248   double * updateBy;
01249   double * reducedCost;
01250   // dj could be very small (or even zero - take care)
01251   double dj = model_->dualIn();
01252   double tolerance=model_->currentDualTolerance();
01253   // we can't really trust infeasibilities if there is dual error
01254   // this coding has to mimic coding in checkDualSolution
01255   double error = min(1.0e-3,model_->largestDualError());
01256   // allow tolerance at least slightly bigger than standard
01257   tolerance = tolerance +  error;
01258   int pivotRow = model_->pivotRow();
01259   double * infeas = infeasible_->denseVector();
01260   //updates->scanAndPack();
01261   model_->factorization()->updateColumnTranspose(spareRow2,updates);
01262   
01263   // put row of tableau in rowArray and columnArray
01264   model_->clpMatrix()->transposeTimes(model_,-1.0,
01265                                       updates,spareColumn2,spareColumn1);
01266   // normal
01267   for (iSection=0;iSection<2;iSection++) {
01268     
01269     reducedCost=model_->djRegion(iSection);
01270     int addSequence;
01271     
01272     if (!iSection) {
01273       number = updates->getNumElements();
01274       index = updates->getIndices();
01275       updateBy = updates->denseVector();
01276       addSequence = model_->numberColumns();
01277     } else {
01278       number = spareColumn1->getNumElements();
01279       index = spareColumn1->getIndices();
01280       updateBy = spareColumn1->denseVector();
01281       addSequence = 0;
01282     }
01283     
01284     for (j=0;j<number;j++) {
01285       int iSequence = index[j];
01286       double value = reducedCost[iSequence];
01287       value -= updateBy[j];
01288       updateBy[j]=0.0;
01289       reducedCost[iSequence] = value;
01290       ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
01291       
01292       switch(status) {
01293         
01294       case ClpSimplex::basic:
01295         infeasible_->zero(iSequence+addSequence);
01296       case ClpSimplex::isFixed:
01297         break;
01298       case ClpSimplex::isFree:
01299       case ClpSimplex::superBasic:
01300         if (fabs(value)>FREE_ACCEPT*tolerance) {
01301           // we are going to bias towards free (but only if reasonable)
01302           value *= FREE_BIAS;
01303           // store square in list
01304           if (infeas[iSequence+addSequence])
01305             infeas[iSequence+addSequence] = value*value; // already there
01306           else
01307             infeasible_->quickAdd(iSequence+addSequence,value*value);
01308         } else {
01309           infeasible_->zero(iSequence+addSequence);
01310         }
01311         break;
01312       case ClpSimplex::atUpperBound:
01313         if (value>tolerance) {
01314           // store square in list
01315           if (infeas[iSequence+addSequence])
01316             infeas[iSequence+addSequence] = value*value; // already there
01317           else
01318             infeasible_->quickAdd(iSequence+addSequence,value*value);
01319         } else {
01320           infeasible_->zero(iSequence+addSequence);
01321         }
01322         break;
01323       case ClpSimplex::atLowerBound:
01324         if (value<-tolerance) {
01325           // store square in list
01326           if (infeas[iSequence+addSequence])
01327             infeas[iSequence+addSequence] = value*value; // already there
01328           else
01329             infeasible_->quickAdd(iSequence+addSequence,value*value);
01330         } else {
01331           infeasible_->zero(iSequence+addSequence);
01332         }
01333       }
01334     }
01335   }
01336   // They are empty
01337   updates->setNumElements(0);
01338   spareColumn1->setNumElements(0);
01339   // make sure infeasibility on incoming is 0.0
01340   int sequenceIn = model_->sequenceIn();
01341   infeasible_->zero(sequenceIn);
01342   // for weights update we use pivotSequence
01343   if (pivotSequence_>=0) {
01344     pivotRow = pivotSequence_;
01345     // unset in case sub flip
01346     pivotSequence_=-1;
01347     // make sure infeasibility on incoming is 0.0
01348     const int * pivotVariable = model_->pivotVariable();
01349     sequenceIn = pivotVariable[pivotRow];
01350     infeasible_->zero(sequenceIn);
01351     // and we can see if reference
01352     double referenceIn=0.0;
01353     if (mode_!=1&&reference(sequenceIn))
01354       referenceIn=1.0;
01355     // save outgoing weight round update
01356     double outgoingWeight=0.0;
01357     int sequenceOut = model_->sequenceOut();
01358     if (sequenceOut>=0)
01359       outgoingWeight=weights_[sequenceOut];
01360     // update weights
01361     updates->setNumElements(0);
01362     spareColumn1->setNumElements(0);
01363     // might as well set dj to 1
01364     dj = 1.0;
01365     updates->insert(pivotRow,-dj);
01366     model_->factorization()->updateColumnTranspose(spareRow2,updates);
01367     // put row of tableau in rowArray and columnArray
01368     model_->clpMatrix()->transposeTimes(model_,-1.0,
01369                                         updates,spareColumn2,spareColumn1);
01370     double * weight;
01371     int numberColumns = model_->numberColumns();
01372     // rows
01373     number = updates->getNumElements();
01374     index = updates->getIndices();
01375     updateBy = updates->denseVector();
01376     weight = weights_+numberColumns;
01377     
01378     assert (devex_>0.0);
01379     for (j=0;j<number;j++) {
01380       int iSequence = index[j];
01381       double thisWeight = weight[iSequence];
01382       // row has -1 
01383       double pivot = - updateBy[iSequence];
01384       updateBy[iSequence]=0.0;
01385       double value = pivot * pivot*devex_;
01386       if (reference(iSequence+numberColumns))
01387         value += 1.0;
01388       weight[iSequence] = max(0.99*thisWeight,value);
01389     }
01390     
01391     // columns
01392     weight = weights_;
01393     
01394     number = spareColumn1->getNumElements();
01395     index = spareColumn1->getIndices();
01396     updateBy = spareColumn1->denseVector();
01397     for (j=0;j<number;j++) {
01398       int iSequence = index[j];
01399       double thisWeight = weight[iSequence];
01400       // row has -1 
01401       double pivot = updateBy[iSequence];
01402       updateBy[iSequence]=0.0;
01403       double value = pivot * pivot*devex_;
01404       if (reference(iSequence))
01405         value += 1.0;
01406       weight[iSequence] = max(0.99*thisWeight,value);
01407     }
01408     // restore outgoing weight
01409     if (sequenceOut>=0)
01410       weights_[sequenceOut]=outgoingWeight;
01411     spareColumn2->setNumElements(0);
01412     //#define SOME_DEBUG_1
01413 #ifdef SOME_DEBUG_1
01414     // check for accuracy
01415     int iCheck=229;
01416     //printf("weight for iCheck is %g\n",weights_[iCheck]);
01417     int numberRows = model_->numberRows();
01418     //int numberColumns = model_->numberColumns();
01419     for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
01420       if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
01421           !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
01422         checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
01423     }
01424 #endif
01425     updates->setNumElements(0);
01426     spareColumn1->setNumElements(0);
01427   }
01428 }
01429 // Update djs, weights for Steepest
01430 void 
01431 ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates,
01432                CoinIndexedVector * spareRow1,
01433                CoinIndexedVector * spareRow2,
01434                CoinIndexedVector * spareColumn1,
01435                CoinIndexedVector * spareColumn2)
01436 {
01437   int iSection,j;
01438   int number=0;
01439   int * index;
01440   double * updateBy;
01441   double * reducedCost;
01442   // dj could be very small (or even zero - take care)
01443   double dj = model_->dualIn();
01444   double tolerance=model_->currentDualTolerance();
01445   // we can't really trust infeasibilities if there is dual error
01446   // this coding has to mimic coding in checkDualSolution
01447   double error = min(1.0e-3,model_->largestDualError());
01448   // allow tolerance at least slightly bigger than standard
01449   tolerance = tolerance +  error;
01450   int pivotRow = model_->pivotRow();
01451   double * infeas = infeasible_->denseVector();
01452   //updates->scanAndPack();
01453   model_->factorization()->updateColumnTranspose(spareRow2,updates);
01454   
01455   // put row of tableau in rowArray and columnArray
01456   model_->clpMatrix()->transposeTimes(model_,-1.0,
01457                                       updates,spareColumn2,spareColumn1);
01458   // normal
01459   for (iSection=0;iSection<2;iSection++) {
01460     
01461     reducedCost=model_->djRegion(iSection);
01462     int addSequence;
01463     
01464     if (!iSection) {
01465       number = updates->getNumElements();
01466       index = updates->getIndices();
01467       updateBy = updates->denseVector();
01468       addSequence = model_->numberColumns();
01469     } else {
01470       number = spareColumn1->getNumElements();
01471       index = spareColumn1->getIndices();
01472       updateBy = spareColumn1->denseVector();
01473       addSequence = 0;
01474     }
01475     
01476     for (j=0;j<number;j++) {
01477       int iSequence = index[j];
01478       double value = reducedCost[iSequence];
01479       value -= updateBy[j];
01480       updateBy[j]=0.0;
01481       reducedCost[iSequence] = value;
01482       ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
01483       
01484       switch(status) {
01485         
01486       case ClpSimplex::basic:
01487         infeasible_->zero(iSequence+addSequence);
01488       case ClpSimplex::isFixed:
01489         break;
01490       case ClpSimplex::isFree:
01491       case ClpSimplex::superBasic:
01492         if (fabs(value)>FREE_ACCEPT*tolerance) {
01493           // we are going to bias towards free (but only if reasonable)
01494           value *= FREE_BIAS;
01495           // store square in list
01496           if (infeas[iSequence+addSequence])
01497             infeas[iSequence+addSequence] = value*value; // already there
01498           else
01499             infeasible_->quickAdd(iSequence+addSequence,value*value);
01500         } else {
01501           infeasible_->zero(iSequence+addSequence);
01502         }
01503         break;
01504       case ClpSimplex::atUpperBound:
01505         if (value>tolerance) {
01506           // store square in list
01507           if (infeas[iSequence+addSequence])
01508             infeas[iSequence+addSequence] = value*value; // already there
01509           else
01510             infeasible_->quickAdd(iSequence+addSequence,value*value);
01511         } else {
01512           infeasible_->zero(iSequence+addSequence);
01513         }
01514         break;
01515       case ClpSimplex::atLowerBound:
01516         if (value<-tolerance) {
01517           // store square in list
01518           if (infeas[iSequence+addSequence])
01519             infeas[iSequence+addSequence] = value*value; // already there
01520           else
01521             infeasible_->quickAdd(iSequence+addSequence,value*value);
01522         } else {
01523           infeasible_->zero(iSequence+addSequence);
01524         }
01525       }
01526     }
01527   }
01528   // we can zero out as will have to get pivot row
01529   // ***** move
01530   updates->setNumElements(0);
01531   spareColumn1->setNumElements(0);
01532   if (pivotRow>=0) {
01533     // make sure infeasibility on incoming is 0.0
01534     int sequenceIn = model_->sequenceIn();
01535     infeasible_->zero(sequenceIn);
01536   }
01537   // for weights update we use pivotSequence
01538   pivotRow = pivotSequence_;
01539   // unset in case sub flip
01540   pivotSequence_=-1;
01541   if (pivotRow>=0) {
01542     // make sure infeasibility on incoming is 0.0
01543     const int * pivotVariable = model_->pivotVariable();
01544     int sequenceIn = pivotVariable[pivotRow];
01545     infeasible_->zero(sequenceIn);
01546     // and we can see if reference
01547     double referenceIn=0.0;
01548     if (mode_!=1&&reference(sequenceIn))
01549       referenceIn=1.0;
01550     // save outgoing weight round update
01551     double outgoingWeight=0.0;
01552     int sequenceOut = model_->sequenceOut();
01553     if (sequenceOut>=0)
01554       outgoingWeight=weights_[sequenceOut];
01555     // update weights
01556     updates->setNumElements(0);
01557     spareColumn1->setNumElements(0);
01558     // might as well set dj to 1
01559     dj = -1.0;
01560     updates->createPacked(1,&pivotRow,&dj);
01561     model_->factorization()->updateColumnTranspose(spareRow2,updates);
01562     // put row of tableau in rowArray and columnArray
01563     model_->clpMatrix()->transposeTimes(model_,-1.0,
01564                                         updates,spareColumn2,spareColumn1);
01565     double * weight;
01566     double * other = alternateWeights_->denseVector();
01567     int numberColumns = model_->numberColumns();
01568     // rows
01569     number = updates->getNumElements();
01570     index = updates->getIndices();
01571     updateBy = updates->denseVector();
01572     weight = weights_+numberColumns;
01573     
01574     if (mode_<4||numberSwitched_>1||mode_>=10) {
01575       // Exact
01576       // now update weight update array
01577       //alternateWeights_->scanAndPack();
01578       model_->factorization()->updateColumnTranspose(spareRow2,
01579                                                      alternateWeights_);
01580       // Exact
01581       // get subset which have nonzero tableau elements
01582       model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
01583                                                 spareColumn1,
01584                                                 spareColumn2);
01585       for (j=0;j<number;j++) {
01586         int iSequence = index[j];
01587         double thisWeight = weight[iSequence];
01588         // row has -1 
01589         double pivot = - updateBy[j];
01590         updateBy[j]=0.0;
01591         double modification = other[iSequence];
01592         double pivotSquared = pivot * pivot;
01593         
01594         thisWeight += pivotSquared * devex_ + pivot * modification;
01595         if (thisWeight<TRY_NORM) {
01596           if (mode_==1) {
01597             // steepest
01598             thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01599           } else {
01600             // exact
01601             thisWeight = referenceIn*pivotSquared;
01602             if (reference(iSequence+numberColumns))
01603               thisWeight += 1.0;
01604             thisWeight = max(thisWeight,TRY_NORM);
01605           }
01606         }
01607         weight[iSequence] = thisWeight;
01608       }
01609     } else if (mode_==4) {
01610       // Devex
01611       assert (devex_>0.0);
01612       for (j=0;j<number;j++) {
01613         int iSequence = index[j];
01614         double thisWeight = weight[iSequence];
01615         // row has -1 
01616         double pivot = -updateBy[j];
01617         updateBy[j]=0.0;
01618         double value = pivot * pivot*devex_;
01619         if (reference(iSequence+numberColumns))
01620           value += 1.0;
01621         weight[iSequence] = max(0.99*thisWeight,value);
01622       }
01623     }
01624     
01625     // columns
01626     weight = weights_;
01627     
01628     number = spareColumn1->getNumElements();
01629     index = spareColumn1->getIndices();
01630     updateBy = spareColumn1->denseVector();
01631     if (mode_<4||numberSwitched_>1||mode_>=10) {
01632       // Exact
01633       double * updateBy2 = spareColumn2->denseVector();
01634       for (j=0;j<number;j++) {
01635         int iSequence = index[j];
01636         double thisWeight = weight[iSequence];
01637         double pivot = updateBy[j];
01638         updateBy[j]=0.0;
01639         double modification = updateBy2[j];
01640         updateBy2[j]=0.0;
01641         double pivotSquared = pivot * pivot;
01642         
01643         thisWeight += pivotSquared * devex_ + pivot * modification;
01644         if (thisWeight<TRY_NORM) {
01645           if (mode_==1) {
01646             // steepest
01647             thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01648           } else {
01649             // exact
01650             thisWeight = referenceIn*pivotSquared;
01651             if (reference(iSequence))
01652               thisWeight += 1.0;
01653             thisWeight = max(thisWeight,TRY_NORM);
01654           }
01655         }
01656         weight[iSequence] = thisWeight;
01657       }
01658     } else if (mode_==4) {
01659       // Devex
01660       for (j=0;j<number;j++) {
01661         int iSequence = index[j];
01662         double thisWeight = weight[iSequence];
01663         // row has -1 
01664         double pivot = updateBy[j];
01665         updateBy[j]=0.0;
01666         double value = pivot * pivot*devex_;
01667         if (reference(iSequence))
01668           value += 1.0;
01669         weight[iSequence] = max(0.99*thisWeight,value);
01670       }
01671     }
01672     // restore outgoing weight
01673     if (sequenceOut>=0)
01674       weights_[sequenceOut]=outgoingWeight;
01675     alternateWeights_->clear();
01676     spareColumn2->setNumElements(0);
01677     //#define SOME_DEBUG_1
01678 #ifdef SOME_DEBUG_1
01679     // check for accuracy
01680     int iCheck=229;
01681     //printf("weight for iCheck is %g\n",weights_[iCheck]);
01682     int numberRows = model_->numberRows();
01683     //int numberColumns = model_->numberColumns();
01684     for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
01685       if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
01686           !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
01687         checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
01688     }
01689 #endif
01690   }
01691   updates->setNumElements(0);
01692   spareColumn1->setNumElements(0);
01693 }
01694 // Update weights for Devex
01695 void 
01696 ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
01697                CoinIndexedVector * spareRow1,
01698                CoinIndexedVector * spareRow2,
01699                CoinIndexedVector * spareColumn1,
01700                CoinIndexedVector * spareColumn2)
01701 {
01702   int j;
01703   int number=0;
01704   int * index;
01705   double * updateBy;
01706   // dj could be very small (or even zero - take care)
01707   double dj = model_->dualIn();
01708   double tolerance=model_->currentDualTolerance();
01709   // we can't really trust infeasibilities if there is dual error
01710   // this coding has to mimic coding in checkDualSolution
01711   double error = min(1.0e-3,model_->largestDualError());
01712   // allow tolerance at least slightly bigger than standard
01713   tolerance = tolerance +  error;
01714   int pivotRow = model_->pivotRow();
01715   // for weights update we use pivotSequence
01716   pivotRow = pivotSequence_;
01717   assert (pivotRow>=0);
01718   // make sure infeasibility on incoming is 0.0
01719   const int * pivotVariable = model_->pivotVariable();
01720   int sequenceIn = pivotVariable[pivotRow];
01721   infeasible_->zero(sequenceIn);
01722   // and we can see if reference
01723   double referenceIn=0.0;
01724   if (mode_!=1&&reference(sequenceIn))
01725     referenceIn=1.0;
01726   // save outgoing weight round update
01727   double outgoingWeight=0.0;
01728   int sequenceOut = model_->sequenceOut();
01729   if (sequenceOut>=0)
01730     outgoingWeight=weights_[sequenceOut];
01731   assert (!updates->getNumElements());
01732   assert (!spareColumn1->getNumElements());
01733   // unset in case sub flip
01734   pivotSequence_=-1;
01735   // might as well set dj to 1
01736   dj = -1.0;
01737   updates->createPacked(1,&pivotRow,&dj);
01738   model_->factorization()->updateColumnTranspose(spareRow2,updates);
01739   // put row of tableau in rowArray and columnArray
01740   model_->clpMatrix()->transposeTimes(model_,-1.0,
01741                                       updates,spareColumn2,spareColumn1);
01742   double * weight;
01743   int numberColumns = model_->numberColumns();
01744   // rows
01745   number = updates->getNumElements();
01746   index = updates->getIndices();
01747   updateBy = updates->denseVector();
01748   weight = weights_+numberColumns;
01749   
01750   // Devex
01751   assert (devex_>0.0);
01752   for (j=0;j<number;j++) {
01753     int iSequence = index[j];
01754     double thisWeight = weight[iSequence];
01755     // row has -1 
01756     double pivot = - updateBy[j];
01757     updateBy[j]=0.0;
01758     double value = pivot * pivot*devex_;
01759     if (reference(iSequence+numberColumns))
01760       value += 1.0;
01761     weight[iSequence] = max(0.99*thisWeight,value);
01762   }
01763   
01764   // columns
01765   weight = weights_;
01766   
01767   number = spareColumn1->getNumElements();
01768   index = spareColumn1->getIndices();
01769   updateBy = spareColumn1->denseVector();
01770   // Devex
01771   for (j=0;j<number;j++) {
01772     int iSequence = index[j];
01773     double thisWeight = weight[iSequence];
01774     // row has -1 
01775     double pivot = updateBy[j];
01776     updateBy[j]=0.0;
01777     double value = pivot * pivot*devex_;
01778     if (reference(iSequence))
01779       value += 1.0;
01780     weight[iSequence] = max(0.99*thisWeight,value);
01781   }
01782   // restore outgoing weight
01783   if (sequenceOut>=0)
01784     weights_[sequenceOut]=outgoingWeight;
01785   spareColumn2->setNumElements(0);
01786   //#define SOME_DEBUG_1
01787 #ifdef SOME_DEBUG_1
01788   // check for accuracy
01789   int iCheck=229;
01790   //printf("weight for iCheck is %g\n",weights_[iCheck]);
01791   int numberRows = model_->numberRows();
01792   //int numberColumns = model_->numberColumns();
01793   for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
01794     if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
01795         !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
01796       checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
01797   }
01798 #endif
01799   updates->setNumElements(0);
01800   spareColumn1->setNumElements(0);
01801 }
01802 // Update weights for Steepest
01803 void 
01804 ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates,
01805                CoinIndexedVector * spareRow1,
01806                CoinIndexedVector * spareRow2,
01807                CoinIndexedVector * spareColumn1,
01808                CoinIndexedVector * spareColumn2)
01809 {
01810   int j;
01811   int number=0;
01812   int * index;
01813   double * updateBy;
01814   // dj could be very small (or even zero - take care)
01815   double dj = model_->dualIn();
01816   double tolerance=model_->currentDualTolerance();
01817   // we can't really trust infeasibilities if there is dual error
01818   // this coding has to mimic coding in checkDualSolution
01819   double error = min(1.0e-3,model_->largestDualError());
01820   // allow tolerance at least slightly bigger than standard
01821   tolerance = tolerance +  error;
01822   int pivotRow = model_->pivotRow();
01823   // for weights update we use pivotSequence
01824   pivotRow = pivotSequence_;
01825   // unset in case sub flip
01826   pivotSequence_=-1;
01827   assert (pivotRow>=0);
01828   // make sure infeasibility on incoming is 0.0
01829   const int * pivotVariable = model_->pivotVariable();
01830   int sequenceIn = pivotVariable[pivotRow];
01831   infeasible_->zero(sequenceIn);
01832   // and we can see if reference
01833   double referenceIn=0.0;
01834   if (mode_!=1&&reference(sequenceIn))
01835     referenceIn=1.0;
01836   // save outgoing weight round update
01837   double outgoingWeight=0.0;
01838   int sequenceOut = model_->sequenceOut();
01839   if (sequenceOut>=0)
01840     outgoingWeight=weights_[sequenceOut];
01841   assert (!updates->getNumElements());
01842   assert (!spareColumn1->getNumElements());
01843   // update weights
01844   //updates->setNumElements(0);
01845   //spareColumn1->setNumElements(0);
01846   // might as well set dj to 1
01847   dj = -1.0;
01848   updates->createPacked(1,&pivotRow,&dj);
01849   model_->factorization()->updateColumnTranspose(spareRow2,updates);
01850   // put row of tableau in rowArray and columnArray
01851   model_->clpMatrix()->transposeTimes(model_,-1.0,
01852                                       updates,spareColumn2,spareColumn1);
01853   double * weight;
01854   double * other = alternateWeights_->denseVector();
01855   int numberColumns = model_->numberColumns();
01856   // rows
01857   number = updates->getNumElements();
01858   index = updates->getIndices();
01859   updateBy = updates->denseVector();
01860   weight = weights_+numberColumns;
01861   
01862   // Exact
01863   // now update weight update array
01864   //alternateWeights_->scanAndPack();
01865   model_->factorization()->updateColumnTranspose(spareRow2,
01866                                                  alternateWeights_);
01867   // get subset which have nonzero tableau elements
01868   model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
01869                                             spareColumn1,
01870                                             spareColumn2);
01871   for (j=0;j<number;j++) {
01872     int iSequence = index[j];
01873     double thisWeight = weight[iSequence];
01874     // row has -1 
01875     double pivot = -updateBy[j];
01876     updateBy[j]=0.0;
01877     double modification = other[iSequence];
01878     double pivotSquared = pivot * pivot;
01879     
01880     thisWeight += pivotSquared * devex_ + pivot * modification;
01881     if (thisWeight<TRY_NORM) {
01882       if (mode_==1) {
01883         // steepest
01884         thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01885       } else {
01886         // exact
01887         thisWeight = referenceIn*pivotSquared;
01888         if (reference(iSequence+numberColumns))
01889           thisWeight += 1.0;
01890         thisWeight = max(thisWeight,TRY_NORM);
01891       }
01892     }
01893     weight[iSequence] = thisWeight;
01894   }
01895   
01896   // columns
01897   weight = weights_;
01898   
01899   number = spareColumn1->getNumElements();
01900   index = spareColumn1->getIndices();
01901   updateBy = spareColumn1->denseVector();
01902   // Exact
01903   double * updateBy2 = spareColumn2->denseVector();
01904   for (j=0;j<number;j++) {
01905     int iSequence = index[j];
01906     double thisWeight = weight[iSequence];
01907     double pivot = updateBy[j];
01908     updateBy[j]=0.0;
01909     double modification = updateBy2[j];
01910     updateBy2[j]=0.0;
01911     double pivotSquared = pivot * pivot;
01912     
01913     thisWeight += pivotSquared * devex_ + pivot * modification;
01914     if (thisWeight<TRY_NORM) {
01915       if (mode_==1) {
01916         // steepest
01917         thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
01918       } else {
01919         // exact
01920         thisWeight = referenceIn*pivotSquared;
01921         if (reference(iSequence))
01922           thisWeight += 1.0;
01923         thisWeight = max(thisWeight,TRY_NORM);
01924       }
01925     }
01926     weight[iSequence] = thisWeight;
01927   }
01928   // restore outgoing weight
01929   if (sequenceOut>=0)
01930     weights_[sequenceOut]=outgoingWeight;
01931   alternateWeights_->clear();
01932   spareColumn2->setNumElements(0);
01933   //#define SOME_DEBUG_1
01934 #ifdef SOME_DEBUG_1
01935   // check for accuracy
01936   int iCheck=229;
01937   //printf("weight for iCheck is %g\n",weights_[iCheck]);
01938   int numberRows = model_->numberRows();
01939   //int numberColumns = model_->numberColumns();
01940   for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
01941     if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
01942         !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
01943       checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
01944   }
01945 #endif
01946   updates->setNumElements(0);
01947   spareColumn1->setNumElements(0);
01948 }
01949 // Returns pivot column, -1 if none
01950 int 
01951 ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates,
01952                                     CoinIndexedVector * spareRow1,
01953                                     CoinIndexedVector * spareRow2,
01954                                     CoinIndexedVector * spareColumn1,
01955                                     CoinIndexedVector * spareColumn2)
01956 {
01957   assert(model_);
01958   int iSection,j;
01959   int number=0;
01960   int * index;
01961   double * updateBy;
01962   double * reducedCost;
01963   // dj could be very small (or even zero - take care)
01964   double dj = model_->dualIn();
01965   double tolerance=model_->currentDualTolerance();
01966   // we can't really trust infeasibilities if there is dual error
01967   // this coding has to mimic coding in checkDualSolution
01968   double error = min(1.0e-3,model_->largestDualError());
01969   // allow tolerance at least slightly bigger than standard
01970   tolerance = tolerance +  error;
01971   int pivotRow = model_->pivotRow();
01972   int anyUpdates;
01973   double * infeas = infeasible_->denseVector();
01974 
01975   // Local copy of mode so can decide what to do
01976   int switchType;
01977   if (mode_==4) 
01978     switchType = 5-numberSwitched_;
01979   else if (mode_>=10)
01980     switchType=3;
01981   else
01982     switchType = mode_;
01983   /* switchType - 
01984      0 - all exact devex
01985      1 - all steepest
01986      2 - some exact devex
01987      3 - auto some exact devex
01988      4 - devex
01989      5 - dantzig
01990   */
01991   if (updates->getNumElements()) {
01992     // would have to have two goes for devex, three for steepest
01993     anyUpdates=2;
01994     // add in pivot contribution
01995     if (pivotRow>=0) 
01996       updates->add(pivotRow,-dj);
01997   } else if (pivotRow>=0) {
01998     if (fabs(dj)>1.0e-15) {
01999       // some dj
02000       updates->insert(pivotRow,-dj);
02001       if (fabs(dj)>1.0e-6) {
02002         // reasonable size
02003         anyUpdates=1;
02004       } else {
02005         // too small
02006         anyUpdates=2;
02007       }
02008     } else {
02009       // zero dj
02010       anyUpdates=-1;
02011     }
02012   } else if (pivotSequence_>=0){
02013     // just after re-factorization
02014     anyUpdates=-1;
02015   } else {
02016     // sub flip - nothing to do
02017     anyUpdates=0;
02018   }
02019 
02020   if (anyUpdates>0) {
02021     model_->factorization()->updateColumnTranspose(spareRow2,updates);
02022     
02023     // put row of tableau in rowArray and columnArray
02024     model_->clpMatrix()->transposeTimes(model_,-1.0,
02025                                         updates,spareColumn2,spareColumn1);
02026     // normal
02027     for (iSection=0;iSection<2;iSection++) {
02028       
02029       reducedCost=model_->djRegion(iSection);
02030       int addSequence;
02031       
02032       if (!iSection) {
02033         number = updates->getNumElements();
02034         index = updates->getIndices();
02035         updateBy = updates->denseVector();
02036         addSequence = model_->numberColumns();
02037       } else {
02038         number = spareColumn1->getNumElements();
02039         index = spareColumn1->getIndices();
02040         updateBy = spareColumn1->denseVector();
02041         addSequence = 0;
02042       }
02043       if (!model_->nonLinearCost()->lookBothWays()) {
02044 
02045         for (j=0;j<number;j++) {
02046           int iSequence = index[j];
02047           double value = reducedCost[iSequence];
02048           value -= updateBy[iSequence];
02049           reducedCost[iSequence] = value;
02050           ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
02051           
02052           switch(status) {
02053             
02054           case ClpSimplex::basic:
02055             infeasible_->zero(iSequence+addSequence);
02056           case ClpSimplex::isFixed:
02057             break;
02058           case ClpSimplex::isFree:
02059           case ClpSimplex::superBasic:
02060             if (fabs(value)>FREE_ACCEPT*tolerance) {
02061               // we are going to bias towards free (but only if reasonable)
02062               value *= FREE_BIAS;
02063               // store square in list
02064               if (infeas[iSequence+addSequence])
02065                 infeas[iSequence+addSequence] = value*value; // already there
02066               else
02067                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02068             } else {
02069               infeasible_->zero(iSequence+addSequence);
02070             }
02071             break;
02072           case ClpSimplex::atUpperBound:
02073             if (value>tolerance) {
02074               // store square in list
02075               if (infeas[iSequence+addSequence])
02076                 infeas[iSequence+addSequence] = value*value; // already there
02077               else
02078                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02079             } else {
02080               infeasible_->zero(iSequence+addSequence);
02081             }
02082             break;
02083           case ClpSimplex::atLowerBound:
02084             if (value<-tolerance) {
02085               // store square in list
02086               if (infeas[iSequence+addSequence])
02087                 infeas[iSequence+addSequence] = value*value; // already there
02088               else
02089                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02090             } else {
02091               infeasible_->zero(iSequence+addSequence);
02092             }
02093           }
02094         }
02095       } else {
02096         ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
02097         // We can go up OR down
02098         for (j=0;j<number;j++) {
02099           int iSequence = index[j];
02100           double value = reducedCost[iSequence];
02101           value -= updateBy[iSequence];
02102           reducedCost[iSequence] = value;
02103           ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
02104           
02105           switch(status) {
02106             
02107           case ClpSimplex::basic:
02108             infeasible_->zero(iSequence+addSequence);
02109           case ClpSimplex::isFixed:
02110             break;
02111           case ClpSimplex::isFree:
02112           case ClpSimplex::superBasic:
02113             if (fabs(value)>FREE_ACCEPT*tolerance) {
02114               // we are going to bias towards free (but only if reasonable)
02115               value *= FREE_BIAS;
02116               // store square in list
02117               if (infeas[iSequence+addSequence])
02118                 infeas[iSequence+addSequence] = value*value; // already there
02119               else
02120                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02121             } else {
02122               infeasible_->zero(iSequence+addSequence);
02123             }
02124             break;
02125           case ClpSimplex::atUpperBound:
02126             if (value>tolerance) {
02127               // store square in list
02128               if (infeas[iSequence+addSequence])
02129                 infeas[iSequence+addSequence] = value*value; // already there
02130               else
02131                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02132             } else {
02133               // look other way - change up should be negative
02134               value -= nonLinear->changeUpInCost(iSequence+addSequence);
02135               if (value>-tolerance) {
02136                 infeasible_->zero(iSequence+addSequence);
02137               } else {
02138                 // store square in list
02139                 if (infeas[iSequence+addSequence])
02140                   infeas[iSequence+addSequence] = value*value; // already there
02141                 else
02142                   infeasible_->quickAdd(iSequence+addSequence,value*value);
02143               }
02144             }
02145             break;
02146           case ClpSimplex::atLowerBound:
02147             if (value<-tolerance) {
02148               // store square in list
02149               if (infeas[iSequence+addSequence])
02150                 infeas[iSequence+addSequence] = value*value; // already there
02151               else
02152                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02153             } else {
02154               // look other way - change down should be positive
02155               value -= nonLinear->changeDownInCost(iSequence+addSequence);
02156               if (value<tolerance) {
02157                 infeasible_->zero(iSequence+addSequence);
02158               } else {
02159                 // store square in list
02160                 if (infeas[iSequence+addSequence])
02161                   infeas[iSequence+addSequence] = value*value; // already there
02162                 else
02163                   infeasible_->quickAdd(iSequence+addSequence,value*value);
02164               }
02165             }
02166           }
02167         }
02168       }
02169     }
02170     if (anyUpdates==2) {
02171       // we can zero out as will have to get pivot row
02172       updates->clear();
02173       spareColumn1->clear();
02174     }
02175     if (pivotRow>=0) {
02176       // make sure infeasibility on incoming is 0.0
02177       int sequenceIn = model_->sequenceIn();
02178       infeasible_->zero(sequenceIn);
02179     }
02180   }
02181   // make sure outgoing from last iteration okay
02182   int sequenceOut = model_->sequenceOut();
02183   if (sequenceOut>=0) {
02184     ClpSimplex::Status status = model_->getStatus(sequenceOut);
02185     double value = model_->reducedCost(sequenceOut);
02186     
02187     switch(status) {
02188       
02189     case ClpSimplex::basic:
02190     case ClpSimplex::isFixed:
02191       break;
02192     case ClpSimplex::isFree:
02193     case ClpSimplex::superBasic:
02194       if (fabs(value)>FREE_ACCEPT*tolerance) { 
02195         // we are going to bias towards free (but only if reasonable)
02196         value *= FREE_BIAS;
02197         // store square in list
02198         if (infeas[sequenceOut])
02199           infeas[sequenceOut] = value*value; // already there
02200         else
02201           infeasible_->quickAdd(sequenceOut,value*value);
02202       } else {
02203         infeasible_->zero(sequenceOut);
02204       }
02205       break;
02206     case ClpSimplex::atUpperBound:
02207       if (value>tolerance) {
02208         // store square in list
02209         if (infeas[sequenceOut])
02210           infeas[sequenceOut] = value*value; // already there
02211         else
02212           infeasible_->quickAdd(sequenceOut,value*value);
02213       } else {
02214         infeasible_->zero(sequenceOut);
02215       }
02216       break;
02217     case ClpSimplex::atLowerBound:
02218       if (value<-tolerance) {
02219         // store square in list
02220         if (infeas[sequenceOut])
02221           infeas[sequenceOut] = value*value; // already there
02222         else
02223           infeasible_->quickAdd(sequenceOut,value*value);
02224       } else {
02225         infeasible_->zero(sequenceOut);
02226       }
02227     }
02228   }
02229 
02230   // If in quadratic re-compute all
02231   if (model_->algorithm()==2) {
02232     for (iSection=0;iSection<2;iSection++) {
02233       
02234       reducedCost=model_->djRegion(iSection);
02235       int addSequence;
02236       int iSequence;
02237       
02238       if (!iSection) {
02239         number = model_->numberRows();
02240         addSequence = model_->numberColumns();
02241       } else {
02242         number = model_->numberColumns();
02243         addSequence = 0;
02244       }
02245       
02246       if (!model_->nonLinearCost()->lookBothWays()) {
02247         for (iSequence=0;iSequence<number;iSequence++) {
02248           double value = reducedCost[iSequence];
02249           ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
02250           
02251           switch(status) {
02252             
02253           case ClpSimplex::basic:
02254             infeasible_->zero(iSequence+addSequence);
02255           case ClpSimplex::isFixed:
02256             break;
02257           case ClpSimplex::isFree:
02258           case ClpSimplex::superBasic:
02259             if (fabs(value)>tolerance) {
02260               // we are going to bias towards free (but only if reasonable)
02261               value *= FREE_BIAS;
02262               // store square in list
02263               if (infeas[iSequence+addSequence])
02264                 infeas[iSequence+addSequence] = value*value; // already there
02265               else
02266                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02267             } else {
02268               infeasible_->zero(iSequence+addSequence);
02269             }
02270             break;
02271           case ClpSimplex::atUpperBound:
02272             if (value>tolerance) {
02273               // store square in list
02274               if (infeas[iSequence+addSequence])
02275                 infeas[iSequence+addSequence] = value*value; // already there
02276               else
02277                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02278             } else {
02279               infeasible_->zero(iSequence+addSequence);
02280             }
02281             break;
02282           case ClpSimplex::atLowerBound:
02283             if (value<-tolerance) {
02284               // store square in list
02285               if (infeas[iSequence+addSequence])
02286                 infeas[iSequence+addSequence] = value*value; // already there
02287               else
02288                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02289             } else {
02290               infeasible_->zero(iSequence+addSequence);
02291             }
02292           }
02293         }
02294       } else {
02295         // we can go both ways
02296         ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
02297         for (iSequence=0;iSequence<number;iSequence++) {
02298           double value = reducedCost[iSequence];
02299           ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
02300           
02301           switch(status) {
02302             
02303           case ClpSimplex::basic:
02304             infeasible_->zero(iSequence+addSequence);
02305           case ClpSimplex::isFixed:
02306             break;
02307           case ClpSimplex::isFree:
02308           case ClpSimplex::superBasic:
02309             if (fabs(value)>tolerance) {
02310               // we are going to bias towards free (but only if reasonable)
02311               value *= FREE_BIAS;
02312               // store square in list
02313               if (infeas[iSequence+addSequence])
02314                 infeas[iSequence+addSequence] = value*value; // already there
02315               else
02316                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02317             } else {
02318               infeasible_->zero(iSequence+addSequence);
02319             }
02320             break;
02321           case ClpSimplex::atUpperBound:
02322             if (value>tolerance) {
02323               // store square in list
02324               if (infeas[iSequence+addSequence])
02325                 infeas[iSequence+addSequence] = value*value; // already there
02326               else
02327                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02328             } else {
02329               // look other way - change up should be negative
02330               value -= nonLinear->changeUpInCost(iSequence+addSequence);
02331               if (value>-tolerance) {
02332                 infeasible_->zero(iSequence+addSequence);
02333               } else {
02334                 // store square in list
02335                 if (infeas[iSequence+addSequence])
02336                   infeas[iSequence+addSequence] = value*value; // already there
02337                 else
02338                   infeasible_->quickAdd(iSequence+addSequence,value*value);
02339               }
02340             }
02341             break;
02342           case ClpSimplex::atLowerBound:
02343             if (value<-tolerance) {
02344               // store square in list
02345               if (infeas[iSequence+addSequence])
02346                 infeas[iSequence+addSequence] = value*value; // already there
02347               else
02348                 infeasible_->quickAdd(iSequence+addSequence,value*value);
02349             } else {
02350               // look other way - change down should be positive
02351               value -= nonLinear->changeDownInCost(iSequence+addSequence);
02352               if (value<tolerance) {
02353                 infeasible_->zero(iSequence+addSequence);
02354               } else {
02355                 // store square in list
02356                 if (infeas[iSequence+addSequence])
02357                   infeas[iSequence+addSequence] = value*value; // already there
02358                 else
02359                   infeasible_->quickAdd(iSequence+addSequence,value*value);
02360               }
02361             }
02362           }
02363         }
02364       }
02365     }
02366   }
02367   // See what sort of pricing
02368   int numberWanted=0;
02369   number = infeasible_->getNumElements();
02370   int numberColumns = model_->numberColumns();
02371   if (switchType==5) {
02372     // we can zero out
02373     updates->clear();
02374     spareColumn1->clear();
02375     pivotSequence_=-1;
02376     pivotRow=-1;
02377     // See if to switch
02378     int numberRows = model_->numberRows();
02379     // ratio is done on number of columns here
02380     //double ratio = (double) sizeFactorization_/(double) numberColumns;
02381     double ratio = (double) sizeFactorization_/(double) numberRows;
02382     //double ratio = (double) sizeFactorization_/(double) model_->clpMatrix()->getNumElements();
02383     if (ratio<0.1) {
02384       numberWanted = max(100,number/200);
02385     } else if (ratio<0.3) {
02386       numberWanted = max(500,number/40);
02387     } else if (ratio<0.5) {
02388       numberWanted = max(2000,number/10);
02389       numberWanted = max(numberWanted,numberColumns/30);
02390     } else {
02391       switchType=4;
02392       // initialize
02393       numberSwitched_++;
02394       // Make sure will re-do
02395       delete [] weights_;
02396       weights_=NULL;
02397       saveWeights(model_,4);
02398       printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
02399       updates->clear();
02400     }
02401     if (model_->numberIterations()%1000==0)
02402       printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
02403   }
02404   if(switchType==4) {
02405     // Still in devex mode
02406     int numberRows = model_->numberRows();
02407     // ratio is done on number of rows here
02408     double ratio = (double) sizeFactorization_/(double) numberRows;
02409     // Go to steepest if lot of iterations?
02410     if (ratio<1.0) {
02411       numberWanted = max(2000,number/20);
02412     } else if (ratio<5.0) {
02413       numberWanted = max(2000,number/10);
02414       numberWanted = max(numberWanted,numberColumns/20);
02415     } else {
02416       // we can zero out
02417       updates->clear();
02418       spareColumn1->clear();
02419       switchType=3;
02420       // initialize
02421       pivotSequence_=-1;
02422       pivotRow=-1;
02423       numberSwitched_++;
02424       // Make sure will re-do
02425       delete [] weights_;
02426       weights_=NULL;
02427       saveWeights(model_,4);
02428       printf("switching to exact %d nel ratio %g\n",sizeFactorization_,ratio);
02429       updates->clear();
02430     }
02431     if (model_->numberIterations()%1000==0)
02432       printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
02433   } 
02434   if (switchType<4) {
02435     if (switchType<2 ) {
02436       numberWanted = number+1;
02437     } else if (switchType==2) {
02438       numberWanted = max(2000,number/8);
02439     } else {
02440       double ratio = (double) sizeFactorization_/(double) model_->numberRows();
02441       if (ratio<1.0) {
02442         numberWanted = max(2000,number/20);
02443       } else if (ratio<5.0) {
02444         numberWanted = max(2000,number/10);
02445         numberWanted = max(numberWanted,numberColumns/20);
02446       } else if (ratio<10.0) {
02447         numberWanted = max(2000,number/8);
02448         numberWanted = max(numberWanted,numberColumns/20);
02449       } else {
02450         ratio = number * (ratio/80.0);
02451         if (ratio>number) {
02452           numberWanted=number+1;
02453         } else {
02454           numberWanted = max(2000,(int) ratio);
02455           numberWanted = max(numberWanted,numberColumns/10);
02456         }
02457       }
02458     }
02459   }
02460   // for weights update we use pivotSequence
02461   pivotRow = pivotSequence_;
02462   // unset in case sub flip
02463   pivotSequence_=-1;
02464   if (pivotRow>=0) {
02465     // make sure infeasibility on incoming is 0.0
02466     const int * pivotVariable = model_->pivotVariable();
02467     int sequenceIn = pivotVariable[pivotRow];
02468     infeasible_->zero(sequenceIn);
02469     // and we can see if reference
02470     double referenceIn=0.0;
02471     if (switchType!=1&&reference(sequenceIn))
02472       referenceIn=1.0;
02473     // save outgoing weight round update
02474     double outgoingWeight=0.0;
02475     if (sequenceOut>=0)
02476       outgoingWeight=weights_[sequenceOut];
02477     // update weights
02478     if (anyUpdates!=1) {
02479       updates->setNumElements(0);
02480       spareColumn1->setNumElements(0);
02481       // might as well set dj to 1
02482       dj = 1.0;
02483       updates->insert(pivotRow,-dj);
02484       model_->factorization()->updateColumnTranspose(spareRow2,updates);
02485       // put row of tableau in rowArray and columnArray
02486       model_->clpMatrix()->transposeTimes(model_,-1.0,
02487                                           updates,spareColumn2,spareColumn1);
02488     }
02489     double * weight;
02490     double * other = alternateWeights_->denseVector();
02491     int numberColumns = model_->numberColumns();
02492     double scaleFactor = -1.0/dj; // as formula is with 1.0
02493     // rows
02494     number = updates->getNumElements();
02495     index = updates->getIndices();
02496     updateBy = updates->denseVector();
02497     weight = weights_+numberColumns;
02498       
02499     if (switchType<4) {
02500       // Exact
02501       // now update weight update array
02502       model_->factorization()->updateColumnTranspose(spareRow2,
02503                                                      alternateWeights_);
02504       for (j=0;j<number;j++) {
02505         int iSequence = index[j];
02506         double thisWeight = weight[iSequence];
02507         // row has -1 
02508         double pivot = updateBy[iSequence]*scaleFactor;
02509         updateBy[iSequence]=0.0;
02510         double modification = other[iSequence];
02511         double pivotSquared = pivot * pivot;
02512         
02513         thisWeight += pivotSquared * devex_ + pivot * modification;
02514         if (thisWeight<TRY_NORM) {
02515           if (switchType==1) {
02516             // steepest
02517             thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
02518           } else {
02519             // exact
02520             thisWeight = referenceIn*pivotSquared;
02521             if (reference(iSequence+numberColumns))
02522               thisWeight += 1.0;
02523             thisWeight = max(thisWeight,TRY_NORM);
02524           }
02525         }
02526         weight[iSequence] = thisWeight;
02527       }
02528     } else if (switchType==4) {
02529       // Devex
02530       assert (devex_>0.0);
02531       for (j=0;j<number;j++) {
02532         int iSequence = index[j];
02533         double thisWeight = weight[iSequence];
02534         // row has -1 
02535         double pivot = updateBy[iSequence]*scaleFactor;
02536         updateBy[iSequence]=0.0;
02537         double value = pivot * pivot*devex_;
02538         if (reference(iSequence+numberColumns))
02539           value += 1.0;
02540         weight[iSequence] = max(0.99*thisWeight,value);
02541       }
02542     }
02543     
02544     // columns
02545     weight = weights_;
02546     
02547     scaleFactor = -scaleFactor;
02548     
02549     number = spareColumn1->getNumElements();
02550     index = spareColumn1->getIndices();
02551     updateBy = spareColumn1->denseVector();
02552     if (switchType<4) {
02553       // Exact
02554       // get subset which have nonzero tableau elements
02555       model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
02556                                                 spareColumn1,
02557                                                 spareColumn2);
02558       double * updateBy2 = spareColumn2->denseVector();
02559       for (j=0;j<number;j++) {
02560         int iSequence = index[j];
02561         double thisWeight = weight[iSequence];
02562         double pivot = updateBy[iSequence]*scaleFactor;
02563         updateBy[iSequence]=0.0;
02564         double modification = updateBy2[j];
02565         updateBy2[j]=0.0;
02566         double pivotSquared = pivot * pivot;
02567         
02568         thisWeight += pivotSquared * devex_ + pivot * modification;
02569         if (thisWeight<TRY_NORM) {
02570           if (switchType==1) {
02571             // steepest
02572             thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
02573           } else {
02574             // exact
02575             thisWeight = referenceIn*pivotSquared;
02576             if (reference(iSequence))
02577               thisWeight += 1.0;
02578             thisWeight = max(thisWeight,TRY_NORM);
02579           }
02580         }
02581         weight[iSequence] = thisWeight;
02582       }
02583     } else if (switchType==4) {
02584       // Devex
02585       for (j=0;j<number;j++) {
02586         int iSequence = index[j];
02587         double thisWeight = weight[iSequence];
02588         // row has -1 
02589         double pivot = updateBy[iSequence]*scaleFactor;
02590         updateBy[iSequence]=0.0;
02591         double value = pivot * pivot*devex_;
02592         if (reference(iSequence))
02593           value += 1.0;
02594         weight[iSequence] = max(0.99*thisWeight,value);
02595       }
02596     }
02597     // restore outgoing weight
02598     if (sequenceOut>=0)
02599       weights_[sequenceOut]=outgoingWeight;
02600     alternateWeights_->clear();
02601     spareColumn2->setNumElements(0);
02602     //#define SOME_DEBUG_1
02603 #ifdef SOME_DEBUG_1
02604     // check for accuracy
02605     int iCheck=229;
02606     //printf("weight for iCheck is %g\n",weights_[iCheck]);
02607     int numberRows = model_->numberRows();
02608     //int numberColumns = model_->numberColumns();
02609     for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
02610       if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
02611           !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
02612         checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
02613     }
02614 #endif
02615     updates->setNumElements(0);
02616     spareColumn1->setNumElements(0);
02617   }
02618 
02619   // update of duals finished - now do pricing
02620 
02621 
02622   double bestDj = 1.0e-30;
02623   int bestSequence=-1;
02624 
02625   int i,iSequence;
02626   index = infeasible_->getIndices();
02627   number = infeasible_->getNumElements();
02628   if(model_->numberIterations()<model_->lastBadIteration()+200) {
02629     // we can't really trust infeasibilities if there is dual error
02630     double checkTolerance = 1.0e-8;
02631     if (!model_->factorization()->pivots())
02632       checkTolerance = 1.0e-6;
02633     if (model_->largestDualError()>checkTolerance)
02634       tolerance *= model_->largestDualError()/checkTolerance;
02635     // But cap
02636     tolerance = min(1000.0,tolerance);
02637   }
02638 #ifdef CLP_DEBUG
02639   if (model_->numberDualInfeasibilities()==1) 
02640     printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
02641            model_->largestDualError(),model_,model_->messageHandler(),
02642            number);
02643 #endif
02644   // stop last one coming immediately
02645   double saveOutInfeasibility=0.0;
02646   if (sequenceOut>=0) {
02647     saveOutInfeasibility = infeas[sequenceOut];
02648     infeas[sequenceOut]=0.0;
02649   }
02650   tolerance *= tolerance; // as we are using squares
02651 
02652   int iPass;
02653   // Setup two passes
02654   int start[4];
02655   start[1]=number;
02656   start[2]=0;
02657   double dstart = ((double) number) * CoinDrand48();
02658   start[0]=(int) dstart;
02659   start[3]=start[0];
02660   //double largestWeight=0.0;
02661   //double smallestWeight=1.0e100;
02662   for (iPass=0;iPass<2;iPass++) {
02663     int end = start[2*iPass+1];
02664     if (switchType<5) {
02665       for (i=start[2*iPass];i<end;i++) {
02666         iSequence = index[i];
02667         double value = infeas[iSequence];
02668         if (value>tolerance) {
02669           double weight = weights_[iSequence];
02670           //weight=1.0;
02671           if (value>bestDj*weight) {
02672             // check flagged variable and correct dj
02673             if (!model_->flagged(iSequence)) {
02674               bestDj=value/weight;
02675               bestSequence = iSequence;
02676             } else {
02677               // just to make sure we don't exit before got something
02678               numberWanted++;
02679             }
02680           }
02681         }
02682         numberWanted--;
02683         if (!numberWanted)
02684           break;
02685       }
02686     } else {
02687       // Dantzig
02688       for (i=start[2*iPass];i<end;i++) {
02689         iSequence = index[i];
02690         double value = infeas[iSequence];
02691         if (value>tolerance) {
02692           if (value>bestDj) {
02693             // check flagged variable and correct dj
02694             if (!model_->flagged(iSequence)) {
02695               bestDj=value;
02696               bestSequence = iSequence;
02697             } else {
02698               // just to make sure we don't exit before got something
02699               numberWanted++;
02700             }
02701           }
02702         }
02703         numberWanted--;
02704         if (!numberWanted)
02705           break;
02706       }
02707     }
02708     if (!numberWanted)
02709       break;
02710   }
02711   if (sequenceOut>=0) {
02712     infeas[sequenceOut]=saveOutInfeasibility;
02713   }
02714   /*if (model_->numberIterations()%100==0)
02715     printf("%d best %g\n",bestSequence,bestDj);*/
02716   
02717 #ifdef CLP_DEBUG
02718   if (bestSequence>=0) {
02719     if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
02720       assert(model_->reducedCost(bestSequence)<0.0);
02721     if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
02722       assert(model_->reducedCost(bestSequence)>0.0);
02723   }
02724 #endif
02725   return bestSequence;
02726 }
02727 /* 
02728    1) before factorization
02729    2) after factorization
02730    3) just redo infeasibilities
02731    4) restore weights
02732 */
02733 void 
02734 ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
02735 {
02736   model_ = model;
02737   if (mode_==4) {
02738     if (mode==1&&!weights_) 
02739       numberSwitched_=0; // Reset
02740   }
02741   // alternateWeights_ is defined as indexed but is treated oddly
02742   // at times
02743   int numberRows = model_->numberRows();
02744   int numberColumns = model_->numberColumns();
02745   const int * pivotVariable = model_->pivotVariable();
02746   bool doInfeasibilities=true;
02747   if (mode==1&&weights_) {
02748     if (pivotSequence_>=0) {
02749       // save pivot order
02750       memcpy(alternateWeights_->getIndices(),pivotVariable,
02751              numberRows*sizeof(int));
02752       // change from pivot row number to sequence number
02753       pivotSequence_=pivotVariable[pivotSequence_];
02754     } else {
02755       alternateWeights_->clear();
02756     }
02757     state_=1;
02758   } else if (mode==2||mode==4||mode==5) {
02759     // restore
02760     if (!weights_||state_==-1||mode==5) {
02761       // Partial is only allowed with certain types of matrix
02762       if (mode_!=4||numberSwitched_||!model_->clpMatrix()->canDoPartialPricing()) {
02763         // initialize weights
02764         delete [] weights_;
02765         delete alternateWeights_;
02766         weights_ = new double[numberRows+numberColumns];
02767         alternateWeights_ = new CoinIndexedVector();
02768         // enough space so can use it for factorization
02769         alternateWeights_->reserve(numberRows+
02770                                    model_->factorization()->maximumPivots());
02771         initializeWeights();
02772         // create saved weights 
02773         delete [] savedWeights_;
02774         savedWeights_ = new double[numberRows+numberColumns];
02775         memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
02776                sizeof(double));
02777       } else {
02778         // Partial pricing
02779         // use region as somewhere to save non-fixed slacks
02780         // set up infeasibilities
02781         if (!infeasible_) {
02782           infeasible_ = new CoinIndexedVector();
02783           infeasible_->reserve(numberColumns+numberRows);
02784         }
02785         infeasible_->clear();
02786         int number = model_->numberRows() + model_->numberColumns();
02787         int iSequence;
02788         int numberLook=0;
02789         int * which = infeasible_->getIndices();
02790         for (iSequence=model_->numberColumns();iSequence<number;iSequence++) {
02791           ClpSimplex::Status status = model_->getStatus(iSequence);
02792           if (status!=ClpSimplex::isFixed)
02793             which[numberLook++]=iSequence;
02794         }
02795         infeasible_->setNumElements(numberLook);
02796         doInfeasibilities=false;
02797       }
02798       savedPivotSequence_=-2;
02799       savedSequenceOut_ = -2;
02800       
02801     } else {
02802       if (mode!=4) {
02803         // save
02804         memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
02805                sizeof(double));
02806         savedPivotSequence_= pivotSequence_;
02807         savedSequenceOut_ = model_->sequenceOut();
02808       } else {
02809         // restore
02810         memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
02811                sizeof(double));
02812         // was - but I think should not be
02813         //pivotSequence_= savedPivotSequence_;
02814         //model_->setSequenceOut(savedSequenceOut_); 
02815         pivotSequence_= -1;
02816         model_->setSequenceOut(-1); 
02817       }
02818     }
02819     state_=0;
02820     // set up infeasibilities
02821     if (!infeasible_) {
02822       infeasible_ = new CoinIndexedVector();
02823       infeasible_->reserve(numberColumns+numberRows);
02824     }
02825   }
02826   if (mode>=2&&mode!=5) {
02827     if (mode!=3&&pivotSequence_>=0) {
02828       // restore pivot row
02829       int iRow;
02830       // permute alternateWeights
02831       double * temp = new double[numberRows+numberColumns];
02832       double * work = alternateWeights_->denseVector();
02833       int * oldPivotOrder = alternateWeights_->getIndices();
02834       for (iRow=0;iRow<numberRows;iRow++) {
02835         int iPivot=oldPivotOrder[iRow];
02836         temp[iPivot]=work[iRow];
02837       }
02838       int number=0;
02839       int found=-1;
02840       int * which = oldPivotOrder;
02841       // find pivot row and re-create alternateWeights
02842       for (iRow=0;iRow<numberRows;iRow++) {
02843         int iPivot=pivotVariable[iRow];
02844         if (iPivot==pivotSequence_) 
02845           found=iRow;
02846         work[iRow]=temp[iPivot];
02847         if (work[iRow])
02848           which[number++]=iRow;
02849       }
02850       alternateWeights_->setNumElements(number);
02851 #ifdef CLP_DEBUG
02852       // Can happen but I should clean up
02853       assert(found>=0);
02854 #endif
02855       pivotSequence_ = found;
02856       delete [] temp;
02857     }
02858     // Save size of factorization
02859     if (!model->factorization()->pivots())
02860       sizeFactorization_ = model_->factorization()->numberElements();
02861     if(!doInfeasibilities)
02862       return; // don't disturb infeasibilities
02863     infeasible_->clear();
02864     double tolerance=model_->currentDualTolerance();
02865     int number = model_->numberRows() + model_->numberColumns();
02866     int iSequence;
02867    
02868     double * reducedCost = model_->djRegion();
02869       
02870     if (!model_->nonLinearCost()->lookBothWays()) {
02871       for (iSequence=0;iSequence<number;iSequence++) {
02872         double value = reducedCost[iSequence];
02873         ClpSimplex::Status status = model_->getStatus(iSequence);
02874 
02875         switch(status) {
02876           
02877         case ClpSimplex::basic:
02878         case ClpSimplex::isFixed:
02879           break;
02880         case ClpSimplex::isFree:
02881         case ClpSimplex::superBasic:
02882           if (fabs(value)>FREE_ACCEPT*tolerance) { 
02883             // we are going to bias towards free (but only if reasonable)
02884             value *= FREE_BIAS;
02885             // store square in list
02886             infeasible_->quickAdd(iSequence,value*value);
02887           }
02888           break;
02889         case ClpSimplex::atUpperBound:
02890           if (value>tolerance) {
02891             infeasible_->quickAdd(iSequence,value*value);
02892           }
02893           break;
02894         case ClpSimplex::atLowerBound:
02895           if (value<-tolerance) {
02896             infeasible_->quickAdd(iSequence,value*value);
02897           }
02898         }
02899       }
02900     } else {
02901       ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
02902       // can go both ways
02903       for (iSequence=0;iSequence<number;iSequence++) {
02904         double value = reducedCost[iSequence];
02905         ClpSimplex::Status status = model_->getStatus(iSequence);
02906 
02907         switch(status) {
02908           
02909         case ClpSimplex::basic:
02910         case ClpSimplex::isFixed:
02911           break;
02912         case ClpSimplex::isFree:
02913         case ClpSimplex::superBasic:
02914           if (fabs(value)>FREE_ACCEPT*tolerance) { 
02915             // we are going to bias towards free (but only if reasonable)
02916             value *= FREE_BIAS;
02917             // store square in list
02918             infeasible_->quickAdd(iSequence,value*value);
02919           }
02920           break;
02921         case ClpSimplex::atUpperBound:
02922           if (value>tolerance) {
02923             infeasible_->quickAdd(iSequence,value*value);
02924           } else {
02925             // look other way - change up should be negative
02926             value -= nonLinear->changeUpInCost(iSequence);
02927             if (value<-tolerance) {
02928               // store square in list
02929               infeasible_->quickAdd(iSequence,value*value);
02930             }
02931           }
02932           break;
02933         case ClpSimplex::atLowerBound:
02934           if (value<-tolerance) {
02935             infeasible_->quickAdd(iSequence,value*value);
02936           } else {
02937             // look other way - change down should be positive
02938             value -= nonLinear->changeDownInCost(iSequence);
02939             if (value>tolerance) {
02940               // store square in list
02941               infeasible_->quickAdd(iSequence,value*value);
02942             }
02943           }
02944         }
02945       }
02946     }
02947   }
02948 }
02949 // Gets rid of last update
02950 void 
02951 ClpPrimalColumnSteepest::unrollWeights()
02952 {
02953   if (mode_==4&&!numberSwitched_)
02954     return;
02955   double * saved = alternateWeights_->denseVector();
02956   int number = alternateWeights_->getNumElements();
02957   int * which = alternateWeights_->getIndices();
02958   int i;
02959   for (i=0;i<number;i++) {
02960     int iRow = which[i];
02961     weights_[iRow]=saved[iRow];
02962     saved[iRow]=0.0;
02963   }
02964   alternateWeights_->setNumElements(0);
02965 }
02966   
02967 //-------------------------------------------------------------------
02968 // Clone
02969 //-------------------------------------------------------------------
02970 ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
02971 {
02972   if (CopyData) {
02973     return new ClpPrimalColumnSteepest(*this);
02974   } else {
02975     return new ClpPrimalColumnSteepest();
02976   }
02977 }
02978 void
02979 ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
02980 {
02981   // Local copy of mode so can decide what to do
02982   int switchType = mode_;
02983   if (mode_==4&&numberSwitched_)
02984     switchType=3;
02985   else if (mode_==4)
02986     return;
02987   int number=input->getNumElements();
02988   int * which = input->getIndices();
02989   double * work = input->denseVector();
02990   int newNumber = 0;
02991   int * newWhich = alternateWeights_->getIndices();
02992   double * newWork = alternateWeights_->denseVector();
02993   int i;
02994   int sequenceIn = model_->sequenceIn();
02995   int sequenceOut = model_->sequenceOut();
02996   const int * pivotVariable = model_->pivotVariable();
02997 
02998   int pivotRow = model_->pivotRow();
02999   pivotSequence_ = pivotRow;
03000 
03001   devex_ =0.0;
03002   // Can't create alternateWeights_ as packed as needed unpacked
03003   if (!input->packedMode()) {
03004     if (pivotRow>=0) {
03005       if (switchType==1) {
03006         for (i=0;i<number;i++) {
03007           int iRow = which[i];
03008           devex_ += work[iRow]*work[iRow];
03009           newWork[iRow] = -2.0*work[iRow];
03010         }
03011         newWork[pivotRow] = -2.0*max(devex_,0.0);
03012         devex_+=ADD_ONE;
03013         weights_[sequenceOut]=1.0+ADD_ONE;
03014         memcpy(newWhich,which,number*sizeof(int));
03015         alternateWeights_->setNumElements(number);
03016       } else {
03017         if (mode_!=4||numberSwitched_>1) {
03018           for (i=0;i<number;i++) {
03019             int iRow = which[i];
03020             int iPivot = pivotVariable[iRow];
03021             if (reference(iPivot)) {
03022               devex_ += work[iRow]*work[iRow];
03023               newWork[iRow] = -2.0*work[iRow];
03024               newWhich[newNumber++]=iRow;
03025             }
03026           }
03027           if (!newWork[pivotRow]&&devex_>0.0)
03028             newWhich[newNumber++]=pivotRow; // add if not already in
03029           newWork[pivotRow] = -2.0*max(devex_,0.0);
03030         } else {
03031           for (i=0;i<number;i++) {
03032             int iRow = which[i];
03033             int iPivot = pivotVariable[iRow];
03034             if (reference(iPivot)) 
03035               devex_ += work[iRow]*work[iRow];
03036           }
03037         }
03038         if (reference(sequenceIn)) {
03039           devex_+=1.0;
03040         } else {
03041         }
03042         if (reference(sequenceOut)) {
03043           weights_[sequenceOut]=1.0+1.0;
03044         } else {
03045           weights_[sequenceOut]=1.0;
03046         }
03047         alternateWeights_->setNumElements(newNumber);
03048       }
03049     } else {
03050       if (switchType==1) {
03051         for (i=0;i<number;i++) {
03052           int iRow = which[i];
03053           devex_ += work[iRow]*work[iRow];
03054         }
03055         devex_ += ADD_ONE;
03056       } else {
03057         for (i=0;i<number;i++) {
03058           int iRow = which[i];
03059           int iPivot = pivotVariable[iRow];
03060           if (reference(iPivot)) {
03061             devex_ += work[iRow]*work[iRow];
03062           }
03063         }
03064         if (reference(sequenceIn)) 
03065           devex_+=1.0;
03066       }
03067     }
03068   } else {
03069     // packed input
03070     if (pivotRow>=0) {
03071       if (switchType==1) {
03072         for (i=0;i<number;i++) {
03073           int iRow = which[i];
03074           devex_ += work[i]*work[i];
03075           newWork[iRow] = -2.0*work[i];
03076         }
03077         newWork[pivotRow] = -2.0*max(devex_,0.0);
03078         devex_+=ADD_ONE;
03079         weights_[sequenceOut]=1.0+ADD_ONE;
03080         memcpy(newWhich,which,number*sizeof(int));
03081         alternateWeights_->setNumElements(number);
03082       } else {
03083         if (mode_!=4||numberSwitched_>1) {
03084           for (i=0;i<number;i++) {
03085             int iRow = which[i];
03086             int iPivot = pivotVariable[iRow];
03087             if (reference(iPivot)) {
03088               devex_ += work[i]*work[i];
03089               newWork[iRow] = -2.0*work[i];
03090               newWhich[newNumber++]=iRow;
03091             }
03092           }
03093           if (!newWork[pivotRow]&&devex_>0.0)
03094             newWhich[newNumber++]=pivotRow; // add if not already in
03095           newWork[pivotRow] = -2.0*max(devex_,0.0);
03096         } else {
03097           for (i=0;i<number;i++) {
03098             int iRow = which[i];
03099             int iPivot = pivotVariable[iRow];
03100             if (reference(iPivot)) 
03101               devex_ += work[i]*work[i];
03102           }
03103         }
03104         if (reference(sequenceIn)) {
03105           devex_+=1.0;
03106         } else {
03107         }
03108         if (reference(sequenceOut)) {
03109           weights_[sequenceOut]=1.0+1.0;
03110         } else {
03111           weights_[sequenceOut]=1.0;
03112         }
03113         alternateWeights_->setNumElements(newNumber);
03114       }
03115     } else {
03116       if (switchType==1) {
03117         for (i=0;i<number;i++) {
03118           devex_ += work[i]*work[i];
03119         }
03120         devex_ += ADD_ONE;
03121       } else {
03122         for (i=0;i<number;i++) {
03123           int iRow = which[i];
03124           int iPivot = pivotVariable[iRow];
03125           if (reference(iPivot)) {
03126             devex_ += work[i]*work[i];
03127           }
03128         }
03129         if (reference(sequenceIn)) 
03130           devex_+=1.0;
03131       }
03132     }
03133   }
03134   double oldDevex = weights_[sequenceIn];
03135 #ifdef CLP_DEBUG
03136   if ((model_->messageHandler()->logLevel()&32))
03137     printf("old weight %g, new %g\n",oldDevex,devex_);
03138 #endif
03139   double check = max(devex_,oldDevex)+0.1;
03140   weights_[sequenceIn] = devex_;
03141   double testValue=0.1;
03142   if (mode_==4&&numberSwitched_==1)
03143     testValue = 0.5;
03144   if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
03145 #ifdef CLP_DEBUG
03146     if ((model_->messageHandler()->logLevel()&48)==16)
03147       printf("old weight %g, new %g\n",oldDevex,devex_);
03148 #endif
03149     //printf("old weight %g, new %g\n",oldDevex,devex_);
03150     testValue=0.99;
03151     if (mode_==1) 
03152       testValue=1.01e1; // make unlikely to do if steepest
03153     else if (mode_==4&&numberSwitched_==1)
03154       testValue=0.9;
03155     double difference = fabs(devex_-oldDevex);
03156     if ( difference > testValue * check ) {
03157       // need to redo
03158       model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
03159                                         *model_->messagesPointer())
03160                                           <<oldDevex<<devex_
03161                                           <<CoinMessageEol;
03162       initializeWeights();
03163     }
03164   }
03165   if (pivotRow>=0) {
03166     // set outgoing weight here
03167     weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
03168   }
03169 }
03170 // Checks accuracy - just for debug
03171 void 
03172 ClpPrimalColumnSteepest::checkAccuracy(int sequence,
03173                                        double relativeTolerance,
03174                                        CoinIndexedVector * rowArray1,
03175                                        CoinIndexedVector * rowArray2)
03176 {
03177   if (mode_==4&&!numberSwitched_)
03178     return;
03179   model_->unpack(rowArray1,sequence);
03180   model_->factorization()->updateColumn(rowArray2,rowArray1);
03181   int number=rowArray1->getNumElements();
03182   int * which = rowArray1->getIndices();
03183   double * work = rowArray1->denseVector();
03184   const int * pivotVariable = model_->pivotVariable();
03185   
03186   double devex =0.0;
03187   int i;
03188 
03189   if (mode_==1) {
03190     for (i=0;i<number;i++) {
03191       int iRow = which[i];
03192       devex += work[iRow]*work[iRow];
03193       work[iRow]=0.0;
03194     }
03195     devex += ADD_ONE;
03196   } else {
03197     for (i=0;i<number;i++) {
03198       int iRow = which[i];
03199       int iPivot = pivotVariable[iRow];
03200       if (reference(iPivot)) {
03201         devex += work[iRow]*work[iRow];
03202       }
03203       work[iRow]=0.0;
03204     }
03205     if (reference(sequence)) 
03206       devex+=1.0;
03207   }
03208 
03209   double oldDevex = weights_[sequence];
03210   double check = max(devex,oldDevex);;
03211   if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
03212     printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
03213     // update so won't print again
03214     weights_[sequence]=devex;
03215   }
03216   rowArray1->setNumElements(0);
03217 }
03218 
03219 // Initialize weights
03220 void 
03221 ClpPrimalColumnSteepest::initializeWeights()
03222 {
03223   int numberRows = model_->numberRows();
03224   int numberColumns = model_->numberColumns();
03225   int number = numberRows + numberColumns;
03226   int iSequence;
03227   if (mode_!=1) {
03228     // initialize to 1.0 
03229     // and set reference framework
03230     if (!reference_) {
03231       int nWords = (number+31)>>5;
03232       reference_ = new unsigned int[nWords];
03233       // tiny overhead to zero out (stops valgrind complaining)
03234       memset(reference_,0,nWords*sizeof(int));
03235     }
03236     
03237     for (iSequence=0;iSequence<number;iSequence++) {
03238       weights_[iSequence]=1.0;
03239       if (model_->getStatus(iSequence)==ClpSimplex::basic) {
03240         setReference(iSequence,false);
03241       } else {
03242         setReference(iSequence,true);
03243       }
03244     }
03245   } else {
03246     CoinIndexedVector * temp = new CoinIndexedVector();
03247     temp->reserve(numberRows+
03248                   model_->factorization()->maximumPivots());
03249     double * array = alternateWeights_->denseVector();
03250     int * which = alternateWeights_->getIndices();
03251       
03252     for (iSequence=0;iSequence<number;iSequence++) {
03253       weights_[iSequence]=1.0+ADD_ONE;
03254       if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
03255           model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
03256         model_->unpack(alternateWeights_,iSequence);
03257         double value=ADD_ONE;
03258         model_->factorization()->updateColumn(temp,alternateWeights_);
03259         int number = alternateWeights_->getNumElements();       int j;
03260         for (j=0;j<number;j++) {
03261           int iRow=which[j];
03262           value+=array[iRow]*array[iRow];
03263           array[iRow]=0.0;
03264         }
03265         alternateWeights_->setNumElements(0);
03266         weights_[iSequence] = value;
03267       }
03268     }
03269     delete temp;
03270   }
03271 }
03272 // Gets rid of all arrays
03273 void 
03274 ClpPrimalColumnSteepest::clearArrays()
03275 {
03276   delete [] weights_;
03277   weights_=NULL;
03278   delete infeasible_;
03279   infeasible_ = NULL;
03280   delete alternateWeights_;
03281   alternateWeights_ = NULL;
03282   delete [] savedWeights_;
03283   savedWeights_ = NULL;
03284   delete [] reference_;
03285   reference_ = NULL;
03286   pivotSequence_=-1;
03287   state_ = -1;
03288   savedPivotSequence_ = -1;
03289   savedSequenceOut_ = -1;  
03290   devex_ = 0.0;
03291 }
03292 // Returns true if would not find any column
03293 bool 
03294 ClpPrimalColumnSteepest::looksOptimal() const
03295 {
03296   //**** THIS MUST MATCH the action coding above
03297   double tolerance=model_->currentDualTolerance();
03298   // we can't really trust infeasibilities if there is dual error
03299   // this coding has to mimic coding in checkDualSolution
03300   double error = min(1.0e-3,model_->largestDualError());
03301   // allow tolerance at least slightly bigger than standard
03302   tolerance = tolerance +  error;
03303   if(model_->numberIterations()<model_->lastBadIteration()+200) {
03304     // we can't really trust infeasibilities if there is dual error
03305     double checkTolerance = 1.0e-8;
03306     if (!model_->factorization()->pivots())
03307       checkTolerance = 1.0e-6;
03308     if (model_->largestDualError()>checkTolerance)
03309       tolerance *= model_->largestDualError()/checkTolerance;
03310     // But cap
03311     tolerance = min(1000.0,tolerance);
03312   }
03313   int number = model_->numberRows() + model_->numberColumns();
03314   int iSequence;
03315   
03316   double * reducedCost = model_->djRegion();
03317   int numberInfeasible=0;
03318   if (!model_->nonLinearCost()->lookBothWays()) {
03319     for (iSequence=0;iSequence<number;iSequence++) {
03320       double value = reducedCost[iSequence];
03321       ClpSimplex::Status status = model_->getStatus(iSequence);
03322       
03323       switch(status) {
03324 
03325       case ClpSimplex::basic:
03326       case ClpSimplex::isFixed:
03327         break;
03328       case ClpSimplex::isFree:
03329       case ClpSimplex::superBasic:
03330         if (fabs(value)>FREE_ACCEPT*tolerance) 
03331           numberInfeasible++;
03332         break;
03333       case ClpSimplex::atUpperBound:
03334         if (value>tolerance) 
03335           numberInfeasible++;
03336         break;
03337       case ClpSimplex::atLowerBound:
03338         if (value<-tolerance) 
03339           numberInfeasible++;
03340       }
03341     }
03342   } else {
03343     ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
03344     // can go both ways
03345     for (iSequence=0;iSequence<number;iSequence++) {
03346       double value = reducedCost[iSequence];
03347       ClpSimplex::Status status = model_->getStatus(iSequence);
03348       
03349       switch(status) {
03350 
03351       case ClpSimplex::basic:
03352       case ClpSimplex::isFixed:
03353         break;
03354       case ClpSimplex::isFree:
03355       case ClpSimplex::superBasic:
03356         if (fabs(value)>FREE_ACCEPT*tolerance) 
03357           numberInfeasible++;
03358         break;
03359       case ClpSimplex::atUpperBound:
03360         if (value>tolerance) {
03361           numberInfeasible++;
03362         } else {
03363           // look other way - change up should be negative
03364           value -= nonLinear->changeUpInCost(iSequence);
03365           if (value<-tolerance) 
03366             numberInfeasible++;
03367         }
03368         break;
03369       case ClpSimplex::atLowerBound:
03370         if (value<-tolerance) {
03371           numberInfeasible++;
03372         } else {
03373           // look other way - change down should be positive
03374           value -= nonLinear->changeDownInCost(iSequence);
03375           if (value>tolerance) 
03376             numberInfeasible++;
03377         }
03378       }
03379     }
03380   }
03381   return numberInfeasible==0;
03382 }
03383 /* Returns number of extra columns for sprint algorithm - 0 means off.
03384    Also number of iterations before recompute
03385 */
03386 int 
03387 ClpPrimalColumnSteepest::numberSprintColumns(int & numberIterations) const
03388 {
03389   numberIterations=0;
03390   int numberAdd=0;
03391   if (!numberSwitched_&&mode_>=10) {
03392     numberIterations = min(2000,model_->numberRows()/5);
03393     numberIterations = max(numberIterations,model_->factorizationFrequency());
03394     numberIterations = max(numberIterations,500);
03395     if (mode_==10) {
03396       numberAdd=max(300,model_->numberColumns()/10);
03397       numberAdd=max(numberAdd,model_->numberRows()/5);
03398       // fake all
03399       //numberAdd=1000000;
03400       numberAdd = min(numberAdd,model_->numberColumns());
03401     } else {
03402       abort();
03403     }
03404   }
03405   return numberAdd;
03406 }
03407 // Switch off sprint idea
03408 void 
03409 ClpPrimalColumnSteepest::switchOffSprint()
03410 {
03411   numberSwitched_=10;
03412 }
03413 // Update djs doing partial pricing (dantzig)
03414 int 
03415 ClpPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
03416                                         CoinIndexedVector * spareRow2,
03417                                         int numberWanted)
03418 {
03419   int number=0;
03420   int * index;
03421   double * updateBy;
03422   double * reducedCost;
03423   double saveTolerance = model_->currentDualTolerance();
03424   double tolerance=model_->currentDualTolerance();
03425   // we can't really trust infeasibilities if there is dual error
03426   // this coding has to mimic coding in checkDualSolution
03427   double error = min(1.0e-3,model_->largestDualError());
03428   // allow tolerance at least slightly bigger than standard
03429   tolerance = tolerance +  error;
03430   if(model_->numberIterations()<model_->lastBadIteration()+200) {
03431     // we can't really trust infeasibilities if there is dual error
03432     double checkTolerance = 1.0e-8;
03433     if (!model_->factorization()->pivots())
03434       checkTolerance = 1.0e-6;
03435     if (model_->largestDualError()>checkTolerance)
03436       tolerance *= model_->largestDualError()/checkTolerance;
03437     // But cap
03438     tolerance = min(1000.0,tolerance);
03439   }
03440   if (model_->factorization()->pivots()&&model_->numberPrimalInfeasibilities())
03441     tolerance = max(tolerance,1.0e-10*model_->infeasibilityCost());
03442   // So partial pricing can use
03443   model_->setCurrentDualTolerance(tolerance);
03444   model_->factorization()->updateColumnTranspose(spareRow2,updates);
03445   int numberColumns = model_->numberColumns();
03446 
03447   // Rows
03448   reducedCost=model_->djRegion(0);
03449   int addSequence;
03450     
03451   number = updates->getNumElements();
03452   index = updates->getIndices();
03453   updateBy = updates->denseVector();
03454   addSequence = numberColumns;
03455   int j;  
03456   double * duals = model_->dualRowSolution();
03457   for (j=0;j<number;j++) {
03458     int iSequence = index[j];
03459     double value = duals[iSequence];
03460     value -= updateBy[j];
03461     updateBy[j]=0.0;
03462     duals[iSequence] = value;
03463   }
03464   double bestDj = tolerance;
03465   int bestSequence=-1;
03466 
03467   const double * cost = model_->costRegion(1);
03468 
03469   int iPassR=0,iPassC=0;
03470   // Setup two passes
03471   // This biases towards picking row variables
03472   // This probably should be fixed
03473   int startR[4];
03474   const int * which = infeasible_->getIndices();
03475   int nSlacks = infeasible_->getNumElements();
03476   startR[1]=nSlacks;
03477   startR[2]=0;
03478   double randomR = CoinDrand48();
03479   double dstart = ((double) nSlacks) * randomR;
03480   startR[0]=(int) dstart;
03481   startR[3]=startR[0];
03482   int startC[4];
03483   startC[1]=numberColumns;
03484   startC[2]=0;
03485   double randomC = CoinDrand48();
03486   dstart = ((double) numberColumns) * randomC;
03487   startC[0]=(int) dstart;
03488   startC[3]=startC[0];
03489   reducedCost = model_->djRegion(1);
03490   int sequenceOut = model_->sequenceOut();
03491   double * duals2 = duals-numberColumns;
03492   int chunk = min(1024,(numberColumns+nSlacks)/32);
03493   if (model_->numberIterations()%1000==0) {
03494     printf("%d wanted, nSlacks %d, chunk %d\n",numberWanted,nSlacks,chunk);
03495     int i;
03496     for (i=0;i<4;i++)
03497       printf("start R %d C %d ",startR[i],startC[i]);
03498     printf("\n");
03499   }
03500   chunk = max(chunk,256);
03501   bool finishedR=false,finishedC=false;
03502   bool doingR = randomR>randomC;
03503   doingR=false;
03504   while (!finishedR||!finishedC) {
03505     if (finishedR)
03506       doingR=false;
03507     if (doingR) {
03508       int saveSequence = bestSequence;
03509       int start = startR[iPassR];
03510       int end = min(startR[iPassR+1],start+chunk/2);
03511       int jSequence;
03512       for (jSequence=start;jSequence<end;jSequence++) {
03513         int iSequence=which[jSequence];
03514         if (iSequence!=sequenceOut) {
03515           double value;
03516           ClpSimplex::Status status = model_->getStatus(iSequence);
03517           
03518           switch(status) {
03519             
03520           case ClpSimplex::basic:
03521           case ClpSimplex::isFixed:
03522             break;
03523           case ClpSimplex::isFree:
03524           case ClpSimplex::superBasic:
03525             value = fabs(cost[iSequence]+duals2[iSequence]);
03526             if (value>FREE_ACCEPT*tolerance) {
03527               numberWanted--;
03528               // we are going to bias towards free (but only if reasonable)
03529               value *= FREE_BIAS;
03530               if (value>bestDj) {
03531                 // check flagged variable and correct dj
03532                 if (!model_->flagged(iSequence)) {
03533                   bestDj=value;
03534                   bestSequence = iSequence;
03535                 } else {
03536                   // just to make sure we don't exit before got something
03537                   numberWanted++;
03538                 }
03539               }
03540             }
03541             break;
03542           case ClpSimplex::atUpperBound:
03543             value = cost[iSequence]+duals2[iSequence];
03544             if (value>tolerance) {
03545               numberWanted--;
03546               if (value>bestDj) {
03547                 // check flagged variable and correct dj
03548                 if (!model_->flagged(iSequence)) {
03549                   bestDj=value;
03550                   bestSequence = iSequence;
03551                 } else {
03552                   // just to make sure we don't exit before got something
03553                   numberWanted++;
03554                 }
03555               }
03556             }
03557             break;
03558           case ClpSimplex::atLowerBound:
03559             value = -(cost[iSequence]+duals2[iSequence]);
03560             if (value>tolerance) {
03561               numberWanted--;
03562               if (value>bestDj) {
03563                 // check flagged variable and correct dj
03564                 if (!model_->flagged(iSequence)) {
03565                   bestDj=value;
03566                   bestSequence = iSequence;
03567                 } else {
03568                   // just to make sure we don't exit before got something
03569                   numberWanted++;
03570                 }
03571               }
03572             }
03573             break;
03574           }
03575         }
03576         if (!numberWanted)
03577           break;
03578       }
03579       if (saveSequence!=bestSequence) {
03580         // dj
03581         reducedCost[bestSequence] = cost[bestSequence] +duals[bestSequence-numberColumns];
03582         bestDj=reducedCost[bestSequence];
03583       }
03584       if (!numberWanted) 
03585         break;
03586       doingR=false;
03587       // update start
03588       startR[iPassR]=jSequence;
03589       if (jSequence>=startR[iPassR+1]) {
03590         if (iPassR)
03591           finishedR=true;
03592         else
03593           iPassR=2;
03594       }
03595     }
03596     if (finishedC)
03597       doingR=true;
03598     if (!doingR) {
03599       int saveSequence = bestSequence;
03600       // Columns
03601       int start = startC[iPassC];
03602       int end = min(startC[iPassC+1],start+chunk);;
03603       model_->clpMatrix()->partialPricing(model_,start,end,bestSequence,numberWanted);
03604       if (saveSequence!=bestSequence) {
03605         // dj
03606         bestDj=reducedCost[bestSequence];
03607       }
03608       if (!numberWanted)
03609         break;
03610       doingR=true;
03611       // update start
03612       int iSequence = end;
03613       startC[iPassC]=iSequence;
03614       if (iSequence>=startC[iPassC+1]) {
03615         if (iPassC)
03616           finishedC=true;
03617         else
03618           iPassC=2;
03619       }
03620     }
03621   }
03622   updates->setNumElements(0);
03623 
03624   // Restore tolerance
03625   model_->setCurrentDualTolerance(saveTolerance);
03626 #ifndef NDEBUG
03627   if (bestSequence>=0) {
03628     if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
03629       assert(model_->reducedCost(bestSequence)<0.0);
03630     if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
03631       assert(model_->reducedCost(bestSequence)>0.0);
03632   }
03633 #endif
03634   return bestSequence;
03635 }