source: branches/pre/ClpPrimalColumnSteepest.cpp @ 210

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

Trying

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include "ClpSimplex.hpp"
7#include "ClpPrimalColumnSteepest.hpp"
8#include "CoinIndexedVector.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpNonLinearCost.hpp"
11#include "ClpMessage.hpp"
12#include "CoinHelperFunctions.hpp"
13#include <stdio.h>
14
15
16// bias for free variables
17#define FREE_BIAS 1.0e1
18// Acceptance criteria for free variables
19#define FREE_ACCEPT 1.0e2
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode) 
28  : ClpPrimalColumnPivot(),
29    devex_(0.0),
30    weights_(NULL),
31    infeasible_(NULL),
32    alternateWeights_(NULL),
33    savedWeights_(NULL),
34    reference_(NULL),
35    state_(-1),
36    mode_(mode),
37    numberSwitched_(0),
38    pivotSequence_(-1),
39    savedPivotSequence_(-1),
40    savedSequenceOut_(-1)
41{
42  type_=2+64*mode;
43}
44
45//-------------------------------------------------------------------
46// Copy constructor
47//-------------------------------------------------------------------
48ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs) 
49: ClpPrimalColumnPivot(rhs)
50{ 
51  state_=rhs.state_;
52  mode_ = rhs.mode_;
53  numberSwitched_ = rhs.numberSwitched_;
54  model_ = rhs.model_;
55  pivotSequence_ = rhs.pivotSequence_;
56  savedPivotSequence_ = rhs.savedPivotSequence_;
57  savedSequenceOut_ = rhs.savedSequenceOut_;
58  devex_ = rhs.devex_;
59  if (rhs.infeasible_) {
60    infeasible_= new CoinIndexedVector(rhs.infeasible_);
61  } else {
62    infeasible_=NULL;
63  }
64  reference_=NULL;
65  if (rhs.weights_) {
66    assert(model_);
67    int number = model_->numberRows()+model_->numberColumns();
68    weights_= new double[number];
69    ClpDisjointCopyN(rhs.weights_,number,weights_);
70    savedWeights_= new double[number];
71    ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
72    if (mode_!=1) {
73      reference_ = new unsigned int[(number+31)>>5];
74      memcpy(reference_,rhs.reference_,((number+31)>>5)*sizeof(unsigned int));
75    }
76  } else {
77    weights_=NULL;
78    savedWeights_=NULL;
79  }
80  if (rhs.alternateWeights_) {
81    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
82  } else {
83    alternateWeights_=NULL;
84  }
85}
86
87//-------------------------------------------------------------------
88// Destructor
89//-------------------------------------------------------------------
90ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
91{
92  delete [] weights_;
93  delete infeasible_;
94  delete alternateWeights_;
95  delete [] savedWeights_;
96  delete [] reference_;
97}
98
99//----------------------------------------------------------------
100// Assignment operator
101//-------------------------------------------------------------------
102ClpPrimalColumnSteepest &
103ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
104{
105  if (this != &rhs) {
106    ClpPrimalColumnPivot::operator=(rhs);
107    state_=rhs.state_;
108    mode_ = rhs.mode_;
109    numberSwitched_ = rhs.numberSwitched_;
110    model_ = rhs.model_;
111    pivotSequence_ = rhs.pivotSequence_;
112    savedPivotSequence_ = rhs.savedPivotSequence_;
113    savedSequenceOut_ = rhs.savedSequenceOut_;
114    devex_ = rhs.devex_;
115    delete [] weights_;
116    delete [] reference_;
117    reference_=NULL;
118    delete infeasible_;
119    delete alternateWeights_;
120    delete [] savedWeights_;
121    savedWeights_ = NULL;
122    if (rhs.infeasible_!=NULL) {
123      infeasible_= new CoinIndexedVector(rhs.infeasible_);
124    } else {
125      infeasible_=NULL;
126    }
127    if (rhs.weights_!=NULL) {
128      assert(model_);
129      int number = model_->numberRows()+model_->numberColumns();
130      weights_= new double[number];
131      ClpDisjointCopyN(rhs.weights_,number,weights_);
132      savedWeights_= new double[number];
133      ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
134      if (mode_!=1) {
135        reference_ = new unsigned int[(number+31)>>5];
136        memcpy(reference_,rhs.reference_,
137               ((number+31)>>5)*sizeof(unsigned int));
138      }
139    } else {
140      weights_=NULL;
141    }
142    if (rhs.alternateWeights_!=NULL) {
143      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
144    } else {
145      alternateWeights_=NULL;
146    }
147  }
148  return *this;
149}
150
151#define TRY_NORM 1.0e-4
152#define ADD_ONE 1.0
153// Returns pivot column, -1 if none
154int 
155ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
156                                    CoinIndexedVector * spareRow1,
157                                    CoinIndexedVector * spareRow2,
158                                    CoinIndexedVector * spareColumn1,
159                                    CoinIndexedVector * spareColumn2)
160{
161  assert(model_);
162  int iSection,j;
163  int number;
164  int * index;
165  double * updateBy;
166  double * reducedCost;
167  // dj could be very small (or even zero - take care)
168  double dj = model_->dualIn();
169  double tolerance=model_->currentDualTolerance();
170  // we can't really trust infeasibilities if there is dual error
171  // this coding has to mimic coding in checkDualSolution
172  double error = min(1.0e-3,model_->largestDualError());
173  // allow tolerance at least slightly bigger than standard
174  tolerance = tolerance +  error;
175  int pivotRow = model_->pivotRow();
176  // Local copy of mode so can decide what to do
177  int switchType = mode_;
178  if (mode_==4) {
179    if (numberSwitched_) {
180      switchType=3;
181    } else {
182      // See if to switch
183      int numberElements = model_->factorization()->numberElements();
184      int numberRows = model_->numberRows();
185      // ratio is done on number of columns here
186      int numberColumns = model_->numberColumns();
187      double ratio = (double) numberElements/(double) numberColumns;
188      //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
189      int numberWanted;
190      bool doPartial=true;
191      int numberTotal = numberColumns+model_->numberRows();
192      if (ratio<0.1) {
193        numberWanted = max(100,numberTotal/200);
194      } else if (ratio<1.0) {
195        numberWanted = max(100,numberTotal/20);
196      } else if (ratio<2.0) {
197        numberWanted = max(200,numberTotal/10);
198      } else {
199        doPartial=false; // switch off
200        numberWanted=0; // just for compliler message
201      }
202      if (doPartial) {
203        if (model_->numberIterations()%100==0)
204          printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
205        // No do partial
206        // Always do slacks
207        int anyUpdates;
208        double * dual = model_->dualRowSolution();
209        if (updates->getNumElements()) {
210          // would have to have two goes for devex, three for steepest
211          anyUpdates=2;
212          // add in pivot contribution
213          if (pivotRow>=0) 
214            updates->add(pivotRow,-dj);
215        } else if (pivotRow>=0) {
216          if (fabs(dj)>1.0e-15) {
217            // some dj
218            updates->insert(pivotRow,-dj);
219            if (fabs(dj)>1.0e-6) {
220              // reasonable size
221              anyUpdates=1;
222            } else {
223              // too small
224              anyUpdates=2;
225            }
226          } else {
227            // zero dj
228            anyUpdates=-1;
229          }
230        } else if (pivotSequence_>=0){
231          // just after re-factorization
232          anyUpdates=-1;
233        } else {
234          // sub flip - nothing to do
235          anyUpdates=0;
236        }
237
238        // For now (as we can always switch off)
239        assert (!model_->nonLinearCost()->lookBothWays());
240        double * slackDj =model_->djRegion(0);
241        double * dj =model_->djRegion(1);
242        if (anyUpdates>0) {
243          model_->factorization()->updateColumnTranspose(spareRow2,updates);
244   
245          number = updates->getNumElements();
246          index = updates->getIndices();
247          updateBy = updates->denseVector();
248         
249
250          for (j=0;j<number;j++) {
251            int iSequence = index[j];
252            double value = updateBy[iSequence];
253            slackDj[iSequence] -= value;
254            dual[iSequence] -= value;
255            updateBy[iSequence]=0.0;
256          }
257          updates->setNumElements(0);
258        }
259       
260        int iPass;
261        // Setup two passes AND treat slacks differently
262        int start[2][4];
263        // slacks
264        start[0][1]=numberRows;
265        start[0][2]=0;
266        double dstart = ((double) numberRows) * CoinDrand48();
267        start[0][0]=(int) dstart;
268        start[0][3]=start[0][0];
269        for (iPass=0;iPass<4;iPass++)
270          start[0][iPass] += numberColumns;
271        start[1][1]=numberColumns;
272        start[1][2]=0;
273        dstart = ((double) numberColumns) * CoinDrand48();
274        start[1][0]=(int) dstart;
275        start[1][3]=start[1][0];
276        int bestSequence=-1;
277        double bestDj=tolerance;
278        // We are going to fake it so pi looks sparse
279        double * saveDense = updates->denseVector();
280        updates->setDenseVector(dual);
281        int kChunk = max(2000,numberWanted);
282        for (iPass=0;iPass<2;iPass++) {
283          // Slacks
284          int end = start[0][2*iPass+1];
285          for (int iSequence=start[0][2*iPass];iSequence<end;iSequence++) {
286            double value;
287            ClpSimplex::Status status = model_->getStatus(iSequence);
288         
289            switch(status) {
290             
291            case ClpSimplex::basic:
292            case ClpSimplex::isFixed:
293              break;
294            case ClpSimplex::isFree:
295            case ClpSimplex::superBasic:
296              value = dj[iSequence];
297              if (fabs(value)>FREE_ACCEPT*tolerance) {
298                numberWanted--;
299                // we are going to bias towards free (but only if reasonable)
300                value = fabs(value*FREE_BIAS);
301                if (value>bestDj) {
302                  // check flagged variable and correct dj
303                  if (!model_->flagged(iSequence)) {
304                    bestDj=value;
305                    bestSequence = iSequence;
306                  } else {
307                    // just to make sure we don't exit before got something
308                    numberWanted++;
309                  }
310                }
311              }
312              break;
313            case ClpSimplex::atUpperBound:
314              value = dj[iSequence];
315              if (value>tolerance) {
316                numberWanted--;
317                if (value>bestDj) {
318                  // check flagged variable and correct dj
319                  if (!model_->flagged(iSequence)) {
320                    bestDj=value;
321                    bestSequence = iSequence;
322                  } else {
323                    // just to make sure we don't exit before got something
324                    numberWanted++;
325                  }
326                }
327              }
328              break;
329            case ClpSimplex::atLowerBound:
330              value = dj[iSequence];
331              if (value<-tolerance) {
332                numberWanted--;
333                if (-value>bestDj) {
334                  // check flagged variable and correct dj
335                  if (!model_->flagged(iSequence)) {
336                    bestDj=-value;
337                    bestSequence = iSequence;
338                  } else {
339                    // just to make sure we don't exit before got something
340                    numberWanted++;
341                  }
342                }
343              }
344            }
345            if (!numberWanted)
346              break;
347          }
348          if (!numberWanted)
349            break;
350          // For Columns break up into manageable chunks
351          end = start[1][2*iPass+1];
352          index = spareColumn1->getIndices();
353          updateBy = spareColumn1->denseVector();
354          for (int jSequence=start[1][2*iPass];jSequence<end;jSequence+=kChunk) {
355            int endThis = min(end,jSequence+kChunk);
356            number=0;
357            for (int iSequence=jSequence;iSequence<endThis;iSequence++) {
358              ClpSimplex::Status status = model_->getStatus(iSequence);
359         
360              if(status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) {
361                index[number++]=iSequence;
362              }
363            }
364            spareColumn1->setNumElements(number);
365            // get subset which have nonzero djs
366            // updates is actually slack djs
367            model_->clpMatrix()->subsetTransposeTimes(model_,updates,
368                                                      spareColumn1,
369                                                      spareColumn2);
370            double * updateBy2 = spareColumn2->denseVector();
371            const double * cost = model_->costRegion();
372            for (int j=0;j<number;j++) {
373              int iSequence = index[j];
374              double value = cost[iSequence]-updateBy2[iSequence];
375              updateBy2[iSequence]=0.0;
376              ClpSimplex::Status status = model_->getStatus(iSequence);
377         
378              switch(status) {
379               
380              case ClpSimplex::basic:
381              case ClpSimplex::isFixed:
382                break;
383              case ClpSimplex::isFree:
384              case ClpSimplex::superBasic:
385                if (fabs(value)>FREE_ACCEPT*tolerance) {
386                  numberWanted--;
387                  // we are going to bias towards free (but only if reasonable)
388                  if (fabs(value*FREE_BIAS)>bestDj) {
389                    // check flagged variable and correct dj
390                    if (!model_->flagged(iSequence)) {
391                      bestDj=value*FREE_BIAS;
392                      dj[iSequence]=value;
393                      bestSequence = iSequence;
394                    } else {
395                      // just to make sure we don't exit before got something
396                      numberWanted++;
397                    }
398                  }
399                }
400                break;
401              case ClpSimplex::atUpperBound:
402                if (value>tolerance) {
403                  numberWanted--;
404                  if (value>bestDj) {
405                    // check flagged variable and correct dj
406                    if (!model_->flagged(iSequence)) {
407                      bestDj=value;
408                      dj[iSequence]=value;
409                      bestSequence = iSequence;
410                    } else {
411                      // just to make sure we don't exit before got something
412                      numberWanted++;
413                    }
414                  }
415                }
416                break;
417              case ClpSimplex::atLowerBound:
418                if (value<-tolerance) {
419                  numberWanted--;
420                  if (-value>bestDj) {
421                    // check flagged variable and correct dj
422                    if (!model_->flagged(iSequence)) {
423                      bestDj=-value;
424                      dj[iSequence]=value;
425                      bestSequence = iSequence;
426                    } else {
427                      // just to make sure we don't exit before got something
428                      numberWanted++;
429                    }
430                  }
431                }
432                break;
433              }
434              if (!numberWanted) {
435                // erase rest
436                for (;j<number;j++) {
437                  int iSequence = index[j];
438                  updateBy2[iSequence]=0.0;
439                }
440                break;
441              }
442            }
443            if (!numberWanted)
444              break;
445          }
446          if (!numberWanted)
447            break;
448        }
449        spareColumn1->setNumElements(0);
450        spareColumn2->setNumElements(0);
451        updates->setDenseVector(saveDense);
452        /*if (model_->numberIterations()%100==0)
453          printf("%d best %g\n",bestSequence,bestDj);*/
454       
455#ifdef CLP_DEBUG
456        if (bestSequence>=0) {
457          if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
458            assert(model_->reducedCost(bestSequence)<0.0);
459          if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
460            assert(model_->reducedCost(bestSequence)>0.0);
461        }
462#endif
463        return bestSequence;
464      } else {
465        // Yes initialize
466        pivotSequence_=-1;
467        pivotRow=-1;
468        numberSwitched_++;
469        // Redo dualsolution
470        model_->computeDuals(NULL);
471        saveWeights(model_,4);
472        printf("switching %d nel ratio %g\n",numberElements,ratio);
473        updates->clear();
474      }
475    }
476  }
477  int anyUpdates;
478  double * infeas = infeasible_->denseVector();
479
480  if (updates->getNumElements()) {
481    // would have to have two goes for devex, three for steepest
482    anyUpdates=2;
483    // add in pivot contribution
484    if (pivotRow>=0) 
485      updates->add(pivotRow,-dj);
486  } else if (pivotRow>=0) {
487    if (fabs(dj)>1.0e-15) {
488      // some dj
489      updates->insert(pivotRow,-dj);
490      if (fabs(dj)>1.0e-6) {
491        // reasonable size
492        anyUpdates=1;
493      } else {
494        // too small
495        anyUpdates=2;
496      }
497    } else {
498      // zero dj
499      anyUpdates=-1;
500    }
501  } else if (pivotSequence_>=0){
502    // just after re-factorization
503    anyUpdates=-1;
504  } else {
505    // sub flip - nothing to do
506    anyUpdates=0;
507  }
508
509  if (anyUpdates>0) {
510    model_->factorization()->updateColumnTranspose(spareRow2,updates);
511   
512    // put row of tableau in rowArray and columnArray
513    model_->clpMatrix()->transposeTimes(model_,-1.0,
514                                        updates,spareColumn2,spareColumn1);
515    // normal
516    for (iSection=0;iSection<2;iSection++) {
517     
518      reducedCost=model_->djRegion(iSection);
519      int addSequence;
520     
521      if (!iSection) {
522        number = updates->getNumElements();
523        index = updates->getIndices();
524        updateBy = updates->denseVector();
525        addSequence = model_->numberColumns();
526      } else {
527        number = spareColumn1->getNumElements();
528        index = spareColumn1->getIndices();
529        updateBy = spareColumn1->denseVector();
530        addSequence = 0;
531      }
532      if (!model_->nonLinearCost()->lookBothWays()) {
533
534        for (j=0;j<number;j++) {
535          int iSequence = index[j];
536          double value = reducedCost[iSequence];
537          value -= updateBy[iSequence];
538          reducedCost[iSequence] = value;
539          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
540         
541          switch(status) {
542           
543          case ClpSimplex::basic:
544            infeasible_->zero(iSequence+addSequence);
545          case ClpSimplex::isFixed:
546            break;
547          case ClpSimplex::isFree:
548          case ClpSimplex::superBasic:
549            if (fabs(value)>FREE_ACCEPT*tolerance) {
550              // we are going to bias towards free (but only if reasonable)
551              value *= FREE_BIAS;
552              // store square in list
553              if (infeas[iSequence+addSequence])
554                infeas[iSequence+addSequence] = value*value; // already there
555              else
556                infeasible_->quickAdd(iSequence+addSequence,value*value);
557            } else {
558              infeasible_->zero(iSequence+addSequence);
559            }
560            break;
561          case ClpSimplex::atUpperBound:
562            if (value>tolerance) {
563              // store square in list
564              if (infeas[iSequence+addSequence])
565                infeas[iSequence+addSequence] = value*value; // already there
566              else
567                infeasible_->quickAdd(iSequence+addSequence,value*value);
568            } else {
569              infeasible_->zero(iSequence+addSequence);
570            }
571            break;
572          case ClpSimplex::atLowerBound:
573            if (value<-tolerance) {
574              // store square in list
575              if (infeas[iSequence+addSequence])
576                infeas[iSequence+addSequence] = value*value; // already there
577              else
578                infeasible_->quickAdd(iSequence+addSequence,value*value);
579            } else {
580              infeasible_->zero(iSequence+addSequence);
581            }
582          }
583        }
584      } else {
585        ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
586        // We can go up OR down
587        for (j=0;j<number;j++) {
588          int iSequence = index[j];
589          double value = reducedCost[iSequence];
590          value -= updateBy[iSequence];
591          reducedCost[iSequence] = value;
592          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
593         
594          switch(status) {
595           
596          case ClpSimplex::basic:
597            infeasible_->zero(iSequence+addSequence);
598          case ClpSimplex::isFixed:
599            break;
600          case ClpSimplex::isFree:
601          case ClpSimplex::superBasic:
602            if (fabs(value)>FREE_ACCEPT*tolerance) {
603              // we are going to bias towards free (but only if reasonable)
604              value *= FREE_BIAS;
605              // store square in list
606              if (infeas[iSequence+addSequence])
607                infeas[iSequence+addSequence] = value*value; // already there
608              else
609                infeasible_->quickAdd(iSequence+addSequence,value*value);
610            } else {
611              infeasible_->zero(iSequence+addSequence);
612            }
613            break;
614          case ClpSimplex::atUpperBound:
615            if (value>tolerance) {
616              // store square in list
617              if (infeas[iSequence+addSequence])
618                infeas[iSequence+addSequence] = value*value; // already there
619              else
620                infeasible_->quickAdd(iSequence+addSequence,value*value);
621            } else {
622              // look other way - change up should be negative
623              value -= nonLinear->changeUpInCost(iSequence+addSequence);
624              if (value>-tolerance) {
625                infeasible_->zero(iSequence+addSequence);
626              } else {
627                // store square in list
628                if (infeas[iSequence+addSequence])
629                  infeas[iSequence+addSequence] = value*value; // already there
630                else
631                  infeasible_->quickAdd(iSequence+addSequence,value*value);
632              }
633            }
634            break;
635          case ClpSimplex::atLowerBound:
636            if (value<-tolerance) {
637              // store square in list
638              if (infeas[iSequence+addSequence])
639                infeas[iSequence+addSequence] = value*value; // already there
640              else
641                infeasible_->quickAdd(iSequence+addSequence,value*value);
642            } else {
643              // look other way - change down should be positive
644              value -= nonLinear->changeDownInCost(iSequence+addSequence);
645              if (value<tolerance) {
646                infeasible_->zero(iSequence+addSequence);
647              } else {
648                // store square in list
649                if (infeas[iSequence+addSequence])
650                  infeas[iSequence+addSequence] = value*value; // already there
651                else
652                  infeasible_->quickAdd(iSequence+addSequence,value*value);
653              }
654            }
655          }
656        }
657      }
658    }
659    if (anyUpdates==2) {
660      // we can zero out as will have to get pivot row
661      updates->clear();
662      spareColumn1->clear();
663    }
664    if (pivotRow>=0) {
665      // make sure infeasibility on incoming is 0.0
666      int sequenceIn = model_->sequenceIn();
667      infeasible_->zero(sequenceIn);
668    }
669  }
670  // make sure outgoing from last iteration okay
671  int sequenceOut = model_->sequenceOut();
672  if (sequenceOut>=0) {
673    ClpSimplex::Status status = model_->getStatus(sequenceOut);
674    double value = model_->reducedCost(sequenceOut);
675   
676    switch(status) {
677     
678    case ClpSimplex::basic:
679    case ClpSimplex::isFixed:
680      break;
681    case ClpSimplex::isFree:
682    case ClpSimplex::superBasic:
683      if (fabs(value)>FREE_ACCEPT*tolerance) { 
684        // we are going to bias towards free (but only if reasonable)
685        value *= FREE_BIAS;
686        // store square in list
687        if (infeas[sequenceOut])
688          infeas[sequenceOut] = value*value; // already there
689        else
690          infeasible_->quickAdd(sequenceOut,value*value);
691      } else {
692        infeasible_->zero(sequenceOut);
693      }
694      break;
695    case ClpSimplex::atUpperBound:
696      if (value>tolerance) {
697        // store square in list
698        if (infeas[sequenceOut])
699          infeas[sequenceOut] = value*value; // already there
700        else
701          infeasible_->quickAdd(sequenceOut,value*value);
702      } else {
703        infeasible_->zero(sequenceOut);
704      }
705      break;
706    case ClpSimplex::atLowerBound:
707      if (value<-tolerance) {
708        // store square in list
709        if (infeas[sequenceOut])
710          infeas[sequenceOut] = value*value; // already there
711        else
712          infeasible_->quickAdd(sequenceOut,value*value);
713      } else {
714        infeasible_->zero(sequenceOut);
715      }
716    }
717  }
718
719  // If in quadratic re-compute all
720  if (model_->algorithm()==2) {
721    for (iSection=0;iSection<2;iSection++) {
722     
723      reducedCost=model_->djRegion(iSection);
724      int addSequence;
725      int iSequence;
726     
727      if (!iSection) {
728        number = model_->numberRows();
729        addSequence = model_->numberColumns();
730      } else {
731        number = model_->numberColumns();
732        addSequence = 0;
733      }
734     
735      if (!model_->nonLinearCost()->lookBothWays()) {
736        for (iSequence=0;iSequence<number;iSequence++) {
737          double value = reducedCost[iSequence];
738          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
739         
740          switch(status) {
741           
742          case ClpSimplex::basic:
743            infeasible_->zero(iSequence+addSequence);
744          case ClpSimplex::isFixed:
745            break;
746          case ClpSimplex::isFree:
747          case ClpSimplex::superBasic:
748            if (fabs(value)>tolerance) {
749              // we are going to bias towards free (but only if reasonable)
750              value *= FREE_BIAS;
751              // store square in list
752              if (infeas[iSequence+addSequence])
753                infeas[iSequence+addSequence] = value*value; // already there
754              else
755                infeasible_->quickAdd(iSequence+addSequence,value*value);
756            } else {
757              infeasible_->zero(iSequence+addSequence);
758            }
759            break;
760          case ClpSimplex::atUpperBound:
761            if (value>tolerance) {
762              // store square in list
763              if (infeas[iSequence+addSequence])
764                infeas[iSequence+addSequence] = value*value; // already there
765              else
766                infeasible_->quickAdd(iSequence+addSequence,value*value);
767            } else {
768              infeasible_->zero(iSequence+addSequence);
769            }
770            break;
771          case ClpSimplex::atLowerBound:
772            if (value<-tolerance) {
773              // store square in list
774              if (infeas[iSequence+addSequence])
775                infeas[iSequence+addSequence] = value*value; // already there
776              else
777                infeasible_->quickAdd(iSequence+addSequence,value*value);
778            } else {
779              infeasible_->zero(iSequence+addSequence);
780            }
781          }
782        }
783      } else {
784        // we can go both ways
785        ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
786        for (iSequence=0;iSequence<number;iSequence++) {
787          double value = reducedCost[iSequence];
788          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
789         
790          switch(status) {
791           
792          case ClpSimplex::basic:
793            infeasible_->zero(iSequence+addSequence);
794          case ClpSimplex::isFixed:
795            break;
796          case ClpSimplex::isFree:
797          case ClpSimplex::superBasic:
798            if (fabs(value)>tolerance) {
799              // we are going to bias towards free (but only if reasonable)
800              value *= FREE_BIAS;
801              // store square in list
802              if (infeas[iSequence+addSequence])
803                infeas[iSequence+addSequence] = value*value; // already there
804              else
805                infeasible_->quickAdd(iSequence+addSequence,value*value);
806            } else {
807              infeasible_->zero(iSequence+addSequence);
808            }
809            break;
810          case ClpSimplex::atUpperBound:
811            if (value>tolerance) {
812              // store square in list
813              if (infeas[iSequence+addSequence])
814                infeas[iSequence+addSequence] = value*value; // already there
815              else
816                infeasible_->quickAdd(iSequence+addSequence,value*value);
817            } else {
818              // look other way - change up should be negative
819              value -= nonLinear->changeUpInCost(iSequence+addSequence);
820              if (value>-tolerance) {
821                infeasible_->zero(iSequence+addSequence);
822              } else {
823                // store square in list
824                if (infeas[iSequence+addSequence])
825                  infeas[iSequence+addSequence] = value*value; // already there
826                else
827                  infeasible_->quickAdd(iSequence+addSequence,value*value);
828              }
829            }
830            break;
831          case ClpSimplex::atLowerBound:
832            if (value<-tolerance) {
833              // store square in list
834              if (infeas[iSequence+addSequence])
835                infeas[iSequence+addSequence] = value*value; // already there
836              else
837                infeasible_->quickAdd(iSequence+addSequence,value*value);
838            } else {
839              // look other way - change down should be positive
840              value -= nonLinear->changeDownInCost(iSequence+addSequence);
841              if (value<tolerance) {
842                infeasible_->zero(iSequence+addSequence);
843              } else {
844                // store square in list
845                if (infeas[iSequence+addSequence])
846                  infeas[iSequence+addSequence] = value*value; // already there
847                else
848                  infeasible_->quickAdd(iSequence+addSequence,value*value);
849              }
850            }
851          }
852        }
853      }
854    }
855  }
856  // for weights update we use pivotSequence
857  pivotRow = pivotSequence_;
858  // unset in case sub flip
859  pivotSequence_=-1;
860  if (pivotRow>=0) {
861    // make sure infeasibility on incoming is 0.0
862    const int * pivotVariable = model_->pivotVariable();
863    int sequenceIn = pivotVariable[pivotRow];
864    infeasible_->zero(sequenceIn);
865    // and we can see if reference
866    double referenceIn=0.0;
867    if (switchType!=1&&reference(sequenceIn))
868      referenceIn=1.0;
869    // save outgoing weight round update
870    double outgoingWeight=0.0;
871    if (sequenceOut>=0)
872      outgoingWeight=weights_[sequenceOut];
873    // update weights
874    if (anyUpdates!=1) {
875      updates->setNumElements(0);
876      spareColumn1->setNumElements(0);
877      // might as well set dj to 1
878      dj = 1.0;
879      updates->insert(pivotRow,-dj);
880      model_->factorization()->updateColumnTranspose(spareRow2,updates);
881      // put row of tableau in rowArray and columnArray
882      model_->clpMatrix()->transposeTimes(model_,-1.0,
883                                          updates,spareColumn2,spareColumn1);
884    }
885    // now update weight update array
886    model_->factorization()->updateColumnTranspose(spareRow2,
887                                                   alternateWeights_);
888    double * other = alternateWeights_->denseVector();
889    double * weight;
890    int numberColumns = model_->numberColumns();
891    double scaleFactor = -1.0/dj; // as formula is with 1.0
892    // rows
893    number = updates->getNumElements();
894    index = updates->getIndices();
895    updateBy = updates->denseVector();
896    weight = weights_+numberColumns;
897   
898    for (j=0;j<number;j++) {
899      int iSequence = index[j];
900      double thisWeight = weight[iSequence];
901      // row has -1
902      double pivot = updateBy[iSequence]*scaleFactor;
903      updateBy[iSequence]=0.0;
904      double modification = other[iSequence];
905      double pivotSquared = pivot * pivot;
906     
907      thisWeight += pivotSquared * devex_ + pivot * modification;
908      if (thisWeight<TRY_NORM) {
909        if (switchType==1) {
910          // steepest
911          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
912        } else {
913          // exact
914          thisWeight = referenceIn*pivotSquared;
915          if (reference(iSequence+numberColumns))
916            thisWeight += 1.0;
917          thisWeight = max(thisWeight,TRY_NORM);
918        }
919      }
920      weight[iSequence] = thisWeight;
921    }
922   
923    // columns
924    weight = weights_;
925   
926    scaleFactor = -scaleFactor;
927   
928    number = spareColumn1->getNumElements();
929    index = spareColumn1->getIndices();
930    updateBy = spareColumn1->denseVector();
931    // get subset which have nonzero tableau elements
932    model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
933                                              spareColumn1,
934                                              spareColumn2);
935    double * updateBy2 = spareColumn2->denseVector();
936    for (j=0;j<number;j++) {
937      int iSequence = index[j];
938      double thisWeight = weight[iSequence];
939      double pivot = updateBy[iSequence]*scaleFactor;
940      updateBy[iSequence]=0.0;
941      double modification = updateBy2[iSequence];
942      updateBy2[iSequence]=0.0;
943      double pivotSquared = pivot * pivot;
944     
945      thisWeight += pivotSquared * devex_ + pivot * modification;
946      if (thisWeight<TRY_NORM) {
947        if (switchType==1) {
948          // steepest
949          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
950        } else {
951          // exact
952          thisWeight = referenceIn*pivotSquared;
953          if (reference(iSequence))
954            thisWeight += 1.0;
955          thisWeight = max(thisWeight,TRY_NORM);
956        }
957      }
958      weight[iSequence] = thisWeight;
959    }
960    // restore outgoing weight
961    if (sequenceOut>=0)
962      weights_[sequenceOut]=outgoingWeight;
963    alternateWeights_->clear();
964    spareColumn2->setNumElements(0);
965#ifdef SOME_DEBUG_1
966    // check for accuracy
967    int iCheck=229;
968    printf("weight for iCheck is %g\n",weights_[iCheck]);
969    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
970        !model_->getStatus(iCheck)!=isFixed)
971      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
972#endif
973    updates->setNumElements(0);
974    spareColumn1->setNumElements(0);
975  }
976
977  // update of duals finished - now do pricing
978
979
980  double bestDj = 1.0e-30;
981  int bestSequence=-1;
982
983  int i,iSequence;
984  index = infeasible_->getIndices();
985  number = infeasible_->getNumElements();
986  if(model_->numberIterations()<model_->lastBadIteration()+200) {
987    // we can't really trust infeasibilities if there is dual error
988    double checkTolerance = 1.0e-8;
989    if (!model_->factorization()->pivots())
990      checkTolerance = 1.0e-6;
991    if (model_->largestDualError()>checkTolerance)
992      tolerance *= model_->largestDualError()/checkTolerance;
993    // But cap
994    tolerance = min(1000.0,tolerance);
995  }
996#ifdef CLP_DEBUG
997  if (model_->numberDualInfeasibilities()==1) 
998    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
999           model_->largestDualError(),model_,model_->messageHandler(),
1000           number);
1001#endif
1002  // stop last one coming immediately
1003  double saveOutInfeasibility=0.0;
1004  if (sequenceOut>=0) {
1005    saveOutInfeasibility = infeas[sequenceOut];
1006    infeas[sequenceOut]=0.0;
1007  }
1008  tolerance *= tolerance; // as we are using squares
1009
1010  int numberWanted;
1011  if (switchType<2 ) {
1012    numberWanted = number+1;
1013  } else if (switchType==2) {
1014    numberWanted = max(2000,number/8);
1015  } else {
1016    int numberElements = model_->factorization()->numberElements();
1017    double ratio = (double) numberElements/(double) model_->numberRows();
1018    numberWanted = max(2000,number/8);
1019    if (ratio<1.0) {
1020      numberWanted = max(2000,number/20);
1021    } else if (ratio>10.0) {
1022      ratio = number * (ratio/80.0);
1023      if (ratio>number)
1024        numberWanted=number+1;
1025      else
1026        numberWanted = max(2000,(int) ratio);
1027    }
1028  }
1029  int iPass;
1030  // Setup two passes
1031  int start[4];
1032  start[1]=number;
1033  start[2]=0;
1034  double dstart = ((double) number) * CoinDrand48();
1035  start[0]=(int) dstart;
1036  start[3]=start[0];
1037  //double largestWeight=0.0;
1038  //double smallestWeight=1.0e100;
1039  for (iPass=0;iPass<2;iPass++) {
1040    int end = start[2*iPass+1];
1041    for (i=start[2*iPass];i<end;i++) {
1042      iSequence = index[i];
1043      double value = infeas[iSequence];
1044      if (value>tolerance) {
1045        double weight = weights_[iSequence];
1046        //weight=1.0;
1047        if (value>bestDj*weight) {
1048          // check flagged variable and correct dj
1049          if (!model_->flagged(iSequence)) {
1050            bestDj=value/weight;
1051            bestSequence = iSequence;
1052          } else {
1053            // just to make sure we don't exit before got something
1054            numberWanted++;
1055          }
1056        }
1057      }
1058      numberWanted--;
1059      if (!numberWanted)
1060        break;
1061    }
1062    if (!numberWanted)
1063      break;
1064  }
1065  if (sequenceOut>=0) {
1066    infeas[sequenceOut]=saveOutInfeasibility;
1067  }
1068  /*if (model_->numberIterations()%100==0)
1069    printf("%d best %g\n",bestSequence,bestDj);*/
1070 
1071#ifdef CLP_DEBUG
1072  if (bestSequence>=0) {
1073    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
1074      assert(model_->reducedCost(bestSequence)<0.0);
1075    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
1076      assert(model_->reducedCost(bestSequence)>0.0);
1077  }
1078#endif
1079  return bestSequence;
1080}
1081/*
1082   1) before factorization
1083   2) after factorization
1084   3) just redo infeasibilities
1085   4) restore weights
1086*/
1087void 
1088ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
1089{
1090  model_ = model;
1091  if (mode_==4) {
1092    if (mode==1&&!weights_) 
1093      numberSwitched_=0; // Reset
1094    if(!numberSwitched_)
1095      return;
1096  }
1097  // alternateWeights_ is defined as indexed but is treated oddly
1098  // at times
1099  int numberRows = model_->numberRows();
1100  int numberColumns = model_->numberColumns();
1101  const int * pivotVariable = model_->pivotVariable();
1102  if (mode==1&&weights_) {
1103    if (pivotSequence_>=0) {
1104      // save pivot order
1105      memcpy(alternateWeights_->getIndices(),pivotVariable,
1106             numberRows*sizeof(int));
1107      // change from pivot row number to sequence number
1108      pivotSequence_=pivotVariable[pivotSequence_];
1109    } else {
1110      alternateWeights_->clear();
1111    }
1112    state_=1;
1113  } else if (mode==2||mode==4||mode==5) {
1114    // restore
1115    if (!weights_||state_==-1||mode==5) {
1116      // initialize weights
1117      delete [] weights_;
1118      delete alternateWeights_;
1119      weights_ = new double[numberRows+numberColumns];
1120      alternateWeights_ = new CoinIndexedVector();
1121      // enough space so can use it for factorization
1122      alternateWeights_->reserve(numberRows+
1123                                 model_->factorization()->maximumPivots());
1124      initializeWeights();
1125      // create saved weights
1126      delete [] savedWeights_;
1127      savedWeights_ = new double[numberRows+numberColumns];
1128      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
1129             sizeof(double));
1130      savedPivotSequence_=-2;
1131      savedSequenceOut_ = -2;
1132     
1133    } else {
1134      if (mode!=4) {
1135        // save
1136        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
1137               sizeof(double));
1138        savedPivotSequence_= pivotSequence_;
1139        savedSequenceOut_ = model_->sequenceOut();
1140      } else {
1141        // restore
1142        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
1143               sizeof(double));
1144        // was - but I think should not be
1145        //pivotSequence_= savedPivotSequence_;
1146        //model_->setSequenceOut(savedSequenceOut_);
1147        pivotSequence_= -1;
1148        model_->setSequenceOut(-1); 
1149      }
1150    }
1151    state_=0;
1152    // set up infeasibilities
1153    if (!infeasible_) {
1154      infeasible_ = new CoinIndexedVector();
1155      infeasible_->reserve(numberColumns+numberRows);
1156    }
1157  }
1158  if (mode>=2&&mode!=5) {
1159    if (mode!=3&&pivotSequence_>=0) {
1160      // restore pivot row
1161      int iRow;
1162      // permute alternateWeights
1163      double * temp = new double[numberRows+numberColumns];
1164      double * work = alternateWeights_->denseVector();
1165      int * oldPivotOrder = alternateWeights_->getIndices();
1166      for (iRow=0;iRow<numberRows;iRow++) {
1167        int iPivot=oldPivotOrder[iRow];
1168        temp[iPivot]=work[iRow];
1169      }
1170      int number=0;
1171      int found=-1;
1172      int * which = oldPivotOrder;
1173      // find pivot row and re-create alternateWeights
1174      for (iRow=0;iRow<numberRows;iRow++) {
1175        int iPivot=pivotVariable[iRow];
1176        if (iPivot==pivotSequence_) 
1177          found=iRow;
1178        work[iRow]=temp[iPivot];
1179        if (work[iRow])
1180          which[number++]=iRow;
1181      }
1182      alternateWeights_->setNumElements(number);
1183#ifdef CLP_DEBUG
1184      // Can happen but I should clean up
1185      assert(found>=0);
1186#endif
1187      pivotSequence_ = found;
1188      delete [] temp;
1189    }
1190    infeasible_->clear();
1191    double tolerance=model_->currentDualTolerance();
1192    int number = model_->numberRows() + model_->numberColumns();
1193    int iSequence;
1194   
1195    double * reducedCost = model_->djRegion();
1196     
1197    if (!model_->nonLinearCost()->lookBothWays()) {
1198      for (iSequence=0;iSequence<number;iSequence++) {
1199        double value = reducedCost[iSequence];
1200        ClpSimplex::Status status = model_->getStatus(iSequence);
1201
1202        switch(status) {
1203         
1204        case ClpSimplex::basic:
1205        case ClpSimplex::isFixed:
1206          break;
1207        case ClpSimplex::isFree:
1208        case ClpSimplex::superBasic:
1209          if (fabs(value)>FREE_ACCEPT*tolerance) { 
1210            // we are going to bias towards free (but only if reasonable)
1211            value *= FREE_BIAS;
1212            // store square in list
1213            infeasible_->quickAdd(iSequence,value*value);
1214          }
1215          break;
1216        case ClpSimplex::atUpperBound:
1217          if (value>tolerance) {
1218            infeasible_->quickAdd(iSequence,value*value);
1219          }
1220          break;
1221        case ClpSimplex::atLowerBound:
1222          if (value<-tolerance) {
1223            infeasible_->quickAdd(iSequence,value*value);
1224          }
1225        }
1226      }
1227    } else {
1228      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
1229      // can go both ways
1230      for (iSequence=0;iSequence<number;iSequence++) {
1231        double value = reducedCost[iSequence];
1232        ClpSimplex::Status status = model_->getStatus(iSequence);
1233
1234        switch(status) {
1235         
1236        case ClpSimplex::basic:
1237        case ClpSimplex::isFixed:
1238          break;
1239        case ClpSimplex::isFree:
1240        case ClpSimplex::superBasic:
1241          if (fabs(value)>FREE_ACCEPT*tolerance) { 
1242            // we are going to bias towards free (but only if reasonable)
1243            value *= FREE_BIAS;
1244            // store square in list
1245            infeasible_->quickAdd(iSequence,value*value);
1246          }
1247          break;
1248        case ClpSimplex::atUpperBound:
1249          if (value>tolerance) {
1250            infeasible_->quickAdd(iSequence,value*value);
1251          } else {
1252            // look other way - change up should be negative
1253            value -= nonLinear->changeUpInCost(iSequence);
1254            if (value<-tolerance) {
1255              // store square in list
1256              infeasible_->quickAdd(iSequence,value*value);
1257            }
1258          }
1259          break;
1260        case ClpSimplex::atLowerBound:
1261          if (value<-tolerance) {
1262            infeasible_->quickAdd(iSequence,value*value);
1263          } else {
1264            // look other way - change down should be positive
1265            value -= nonLinear->changeDownInCost(iSequence);
1266            if (value>tolerance) {
1267              // store square in list
1268              infeasible_->quickAdd(iSequence,value*value);
1269            }
1270          }
1271        }
1272      }
1273    }
1274  }
1275}
1276// Gets rid of last update
1277void 
1278ClpPrimalColumnSteepest::unrollWeights()
1279{
1280  if (mode_==4&&!numberSwitched_)
1281    return;
1282  double * saved = alternateWeights_->denseVector();
1283  int number = alternateWeights_->getNumElements();
1284  int * which = alternateWeights_->getIndices();
1285  int i;
1286  for (i=0;i<number;i++) {
1287    int iRow = which[i];
1288    weights_[iRow]=saved[iRow];
1289    saved[iRow]=0.0;
1290  }
1291  alternateWeights_->setNumElements(0);
1292}
1293 
1294//-------------------------------------------------------------------
1295// Clone
1296//-------------------------------------------------------------------
1297ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
1298{
1299  if (CopyData) {
1300    return new ClpPrimalColumnSteepest(*this);
1301  } else {
1302    return new ClpPrimalColumnSteepest();
1303  }
1304}
1305void
1306ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
1307{
1308  // Local copy of mode so can decide what to do
1309  int switchType = mode_;
1310  if (mode_==4&&numberSwitched_)
1311    switchType=3;
1312  else if (mode_==4)
1313    return;
1314  int number=input->getNumElements();
1315  int * which = input->getIndices();
1316  double * work = input->denseVector();
1317  int newNumber = 0;
1318  int * newWhich = alternateWeights_->getIndices();
1319  double * newWork = alternateWeights_->denseVector();
1320  int i;
1321  int sequenceIn = model_->sequenceIn();
1322  int sequenceOut = model_->sequenceOut();
1323  const int * pivotVariable = model_->pivotVariable();
1324
1325  int pivotRow = model_->pivotRow();
1326  pivotSequence_ = pivotRow;
1327
1328  devex_ =0.0;
1329
1330  if (pivotRow>=0) {
1331    if (switchType==1) {
1332      for (i=0;i<number;i++) {
1333        int iRow = which[i];
1334        devex_ += work[iRow]*work[iRow];
1335        newWork[iRow] = -2.0*work[iRow];
1336      }
1337      newWork[pivotRow] = -2.0*max(devex_,0.0);
1338      devex_+=ADD_ONE;
1339      weights_[sequenceOut]=1.0+ADD_ONE;
1340      memcpy(newWhich,which,number*sizeof(int));
1341      alternateWeights_->setNumElements(number);
1342    } else {
1343      for (i=0;i<number;i++) {
1344        int iRow = which[i];
1345        int iPivot = pivotVariable[iRow];
1346        if (reference(iPivot)) {
1347          devex_ += work[iRow]*work[iRow];
1348          newWork[iRow] = -2.0*work[iRow];
1349          newWhich[newNumber++]=iRow;
1350        }
1351      }
1352      if (!newWork[pivotRow])
1353        newWhich[newNumber++]=pivotRow; // add if not already in
1354      newWork[pivotRow] = -2.0*max(devex_,0.0);
1355      if (reference(sequenceIn)) {
1356        devex_+=1.0;
1357      } else {
1358      }
1359      if (reference(sequenceOut)) {
1360        weights_[sequenceOut]=1.0+1.0;
1361      } else {
1362        weights_[sequenceOut]=1.0;
1363      }
1364      alternateWeights_->setNumElements(newNumber);
1365    }
1366  } else {
1367    if (switchType==1) {
1368      for (i=0;i<number;i++) {
1369        int iRow = which[i];
1370        devex_ += work[iRow]*work[iRow];
1371      }
1372      devex_ += ADD_ONE;
1373    } else {
1374      for (i=0;i<number;i++) {
1375        int iRow = which[i];
1376        int iPivot = pivotVariable[iRow];
1377        if (reference(iPivot)) {
1378          devex_ += work[iRow]*work[iRow];
1379        }
1380      }
1381      if (reference(sequenceIn)) 
1382        devex_+=1.0;
1383    }
1384  }
1385  double oldDevex = weights_[sequenceIn];
1386#ifdef CLP_DEBUG
1387  if ((model_->messageHandler()->logLevel()&32))
1388    printf("old weight %g, new %g\n",oldDevex,devex_);
1389#endif
1390  double check = max(devex_,oldDevex);;
1391  weights_[sequenceIn] = devex_;
1392  if ( fabs ( devex_ - oldDevex ) > 1.0e-1 * check ) {
1393#ifdef CLP_DEBUG
1394    if ((model_->messageHandler()->logLevel()&48)==16)
1395      printf("old weight %g, new %g\n",oldDevex,devex_);
1396#endif
1397    if ( fabs ( devex_ - oldDevex ) > 1.0e5 * check ) {
1398      // need to redo
1399      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
1400                                        *model_->messagesPointer())
1401                                          <<oldDevex<<devex_
1402                                          <<CoinMessageEol;
1403      initializeWeights();
1404    }
1405  }
1406  if (pivotRow>=0) {
1407    // set outgoing weight here
1408    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
1409  }
1410}
1411// Checks accuracy - just for debug
1412void 
1413ClpPrimalColumnSteepest::checkAccuracy(int sequence,
1414                                       double relativeTolerance,
1415                                       CoinIndexedVector * rowArray1,
1416                                       CoinIndexedVector * rowArray2)
1417{
1418  if (mode_==4&&!numberSwitched_)
1419    return;
1420  model_->unpack(rowArray1,sequence);
1421  model_->factorization()->updateColumn(rowArray2,rowArray1);
1422  int number=rowArray1->getNumElements();
1423  int * which = rowArray1->getIndices();
1424  double * work = rowArray1->denseVector();
1425  const int * pivotVariable = model_->pivotVariable();
1426 
1427  double devex =0.0;
1428  int i;
1429
1430  if (mode_==1) {
1431    for (i=0;i<number;i++) {
1432      int iRow = which[i];
1433      devex += work[iRow]*work[iRow];
1434      work[iRow]=0.0;
1435    }
1436    devex += ADD_ONE;
1437  } else {
1438    for (i=0;i<number;i++) {
1439      int iRow = which[i];
1440      int iPivot = pivotVariable[iRow];
1441      if (reference(iPivot)) {
1442        devex += work[iRow]*work[iRow];
1443      }
1444      work[iRow]=0.0;
1445    }
1446    if (reference(sequence)) 
1447      devex+=1.0;
1448  }
1449
1450  double oldDevex = weights_[sequence];
1451  double check = max(devex,oldDevex);;
1452  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
1453    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
1454    // update so won't print again
1455    weights_[sequence]=devex;
1456  }
1457  rowArray1->setNumElements(0);
1458}
1459
1460// Initialize weights
1461void 
1462ClpPrimalColumnSteepest::initializeWeights()
1463{
1464  int numberRows = model_->numberRows();
1465  int numberColumns = model_->numberColumns();
1466  int number = numberRows + numberColumns;
1467  int iSequence;
1468  if (mode_!=1) {
1469    // initialize to 1.0
1470    // and set reference framework
1471    if (!reference_) {
1472      int nWords = (number+31)>>5;
1473      reference_ = new unsigned int[nWords];
1474      // tiny overhead to zero out (stops valgrind complaining)
1475      memset(reference_,0,nWords*sizeof(int));
1476    }
1477   
1478    for (iSequence=0;iSequence<number;iSequence++) {
1479      weights_[iSequence]=1.0;
1480      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
1481        setReference(iSequence,false);
1482      } else {
1483        setReference(iSequence,true);
1484      }
1485    }
1486  } else {
1487    CoinIndexedVector * temp = new CoinIndexedVector();
1488    temp->reserve(numberRows+
1489                  model_->factorization()->maximumPivots());
1490    double * array = alternateWeights_->denseVector();
1491    int * which = alternateWeights_->getIndices();
1492     
1493    for (iSequence=0;iSequence<number;iSequence++) {
1494      weights_[iSequence]=1.0+ADD_ONE;
1495      if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
1496          model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
1497        model_->unpack(alternateWeights_,iSequence);
1498        double value=ADD_ONE;
1499        model_->factorization()->updateColumn(temp,alternateWeights_);
1500        int number = alternateWeights_->getNumElements();       int j;
1501        for (j=0;j<number;j++) {
1502          int iRow=which[j];
1503          value+=array[iRow]*array[iRow];
1504          array[iRow]=0.0;
1505        }
1506        alternateWeights_->setNumElements(0);
1507        weights_[iSequence] = value;
1508      }
1509    }
1510    delete temp;
1511  }
1512}
1513// Gets rid of all arrays
1514void 
1515ClpPrimalColumnSteepest::clearArrays()
1516{
1517  delete [] weights_;
1518  weights_=NULL;
1519  delete infeasible_;
1520  infeasible_ = NULL;
1521  delete alternateWeights_;
1522  alternateWeights_ = NULL;
1523  delete [] savedWeights_;
1524  savedWeights_ = NULL;
1525  delete [] reference_;
1526  reference_ = NULL;
1527  pivotSequence_=-1;
1528  state_ = -1;
1529  savedPivotSequence_ = -1;
1530  savedSequenceOut_ = -1; 
1531  devex_ = 0.0;
1532}
1533// Returns true if would not find any column
1534bool 
1535ClpPrimalColumnSteepest::looksOptimal() const
1536{
1537  //**** THIS MUST MATCH the action coding above
1538  double tolerance=model_->currentDualTolerance();
1539  // we can't really trust infeasibilities if there is dual error
1540  // this coding has to mimic coding in checkDualSolution
1541  double error = min(1.0e-3,model_->largestDualError());
1542  // allow tolerance at least slightly bigger than standard
1543  tolerance = tolerance +  error;
1544  if(model_->numberIterations()<model_->lastBadIteration()+200) {
1545    // we can't really trust infeasibilities if there is dual error
1546    double checkTolerance = 1.0e-8;
1547    if (!model_->factorization()->pivots())
1548      checkTolerance = 1.0e-6;
1549    if (model_->largestDualError()>checkTolerance)
1550      tolerance *= model_->largestDualError()/checkTolerance;
1551    // But cap
1552    tolerance = min(1000.0,tolerance);
1553  }
1554  int number = model_->numberRows() + model_->numberColumns();
1555  int iSequence;
1556 
1557  double * reducedCost = model_->djRegion();
1558  int numberInfeasible=0;
1559  if (!model_->nonLinearCost()->lookBothWays()) {
1560    for (iSequence=0;iSequence<number;iSequence++) {
1561      double value = reducedCost[iSequence];
1562      ClpSimplex::Status status = model_->getStatus(iSequence);
1563     
1564      switch(status) {
1565
1566      case ClpSimplex::basic:
1567      case ClpSimplex::isFixed:
1568        break;
1569      case ClpSimplex::isFree:
1570      case ClpSimplex::superBasic:
1571        if (fabs(value)>FREE_ACCEPT*tolerance) 
1572          numberInfeasible++;
1573        break;
1574      case ClpSimplex::atUpperBound:
1575        if (value>tolerance) 
1576          numberInfeasible++;
1577        break;
1578      case ClpSimplex::atLowerBound:
1579        if (value<-tolerance) 
1580          numberInfeasible++;
1581      }
1582    }
1583  } else {
1584    ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
1585    // can go both ways
1586    for (iSequence=0;iSequence<number;iSequence++) {
1587      double value = reducedCost[iSequence];
1588      ClpSimplex::Status status = model_->getStatus(iSequence);
1589     
1590      switch(status) {
1591
1592      case ClpSimplex::basic:
1593      case ClpSimplex::isFixed:
1594        break;
1595      case ClpSimplex::isFree:
1596      case ClpSimplex::superBasic:
1597        if (fabs(value)>FREE_ACCEPT*tolerance) 
1598          numberInfeasible++;
1599        break;
1600      case ClpSimplex::atUpperBound:
1601        if (value>tolerance) {
1602          numberInfeasible++;
1603        } else {
1604          // look other way - change up should be negative
1605          value -= nonLinear->changeUpInCost(iSequence);
1606          if (value<-tolerance) 
1607            numberInfeasible++;
1608        }
1609        break;
1610      case ClpSimplex::atLowerBound:
1611        if (value<-tolerance) {
1612          numberInfeasible++;
1613        } else {
1614          // look other way - change down should be positive
1615          value -= nonLinear->changeDownInCost(iSequence);
1616          if (value>tolerance) 
1617            numberInfeasible++;
1618        }
1619      }
1620    }
1621  }
1622  return numberInfeasible==0;
1623}
Note: See TracBrowser for help on using the repository browser.