source: branches/devel-1/ClpPrimalColumnSteepest.cpp @ 117

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

Some tidying up

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8
9#include "ClpSimplex.hpp"
10#include "ClpPrimalColumnSteepest.hpp"
11#include "CoinIndexedVector.hpp"
12#include "ClpFactorization.hpp"
13#include "ClpMessage.hpp"
14#include "CoinHelperFunctions.hpp"
15#include <stdio.h>
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
23ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode) 
24  : ClpPrimalColumnPivot(),
25    state_(-1),
26    mode_(mode),
27    weights_(NULL),
28    infeasible_(NULL),
29    alternateWeights_(NULL),
30    savedWeights_(NULL),
31    pivotSequence_(-1),
32    savedPivotSequence_(-1),
33    savedSequenceOut_(-1), 
34    reference_(NULL),
35    devex_(0.0)
36{
37  type_=2+64*mode;
38}
39
40//-------------------------------------------------------------------
41// Copy constructor
42//-------------------------------------------------------------------
43ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs) 
44: ClpPrimalColumnPivot(rhs)
45{ 
46  state_=rhs.state_;
47  mode_ = rhs.mode_;
48  model_ = rhs.model_;
49  pivotSequence_ = rhs.pivotSequence_;
50  savedPivotSequence_ = rhs.savedPivotSequence_;
51  savedSequenceOut_ = rhs.savedSequenceOut_;
52  devex_ = rhs.devex_;
53  if (rhs.infeasible_) {
54    infeasible_= new CoinIndexedVector(rhs.infeasible_);
55  } else {
56    infeasible_=NULL;
57  }
58  reference_=NULL;
59  if (rhs.weights_) {
60    assert(model_);
61    int number = model_->numberRows()+model_->numberColumns();
62    weights_= new double[number];
63    ClpDisjointCopyN(rhs.weights_,number,weights_);
64    savedWeights_= new double[number];
65    ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
66    if (!mode_) {
67      reference_ = new unsigned int[(number+31)>>5];
68      memcpy(reference_,rhs.reference_,((number+31)>>5)*sizeof(unsigned int));
69    }
70  } else {
71    weights_=NULL;
72    savedWeights_=NULL;
73  }
74  if (rhs.alternateWeights_) {
75    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
76  } else {
77    alternateWeights_=NULL;
78  }
79}
80
81//-------------------------------------------------------------------
82// Destructor
83//-------------------------------------------------------------------
84ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
85{
86  delete [] weights_;
87  delete infeasible_;
88  delete alternateWeights_;
89  delete [] savedWeights_;
90  delete [] reference_;
91}
92
93//----------------------------------------------------------------
94// Assignment operator
95//-------------------------------------------------------------------
96ClpPrimalColumnSteepest &
97ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
98{
99  if (this != &rhs) {
100    ClpPrimalColumnPivot::operator=(rhs);
101    state_=rhs.state_;
102    mode_ = rhs.mode_;
103    model_ = rhs.model_;
104    pivotSequence_ = rhs.pivotSequence_;
105    savedPivotSequence_ = rhs.savedPivotSequence_;
106    savedSequenceOut_ = rhs.savedSequenceOut_;
107    devex_ = rhs.devex_;
108    delete [] weights_;
109    delete [] reference_;
110    reference_=NULL;
111    delete infeasible_;
112    delete alternateWeights_;
113    delete [] savedWeights_;
114    savedWeights_ = NULL;
115    if (rhs.infeasible_!=NULL) {
116      infeasible_= new CoinIndexedVector(rhs.infeasible_);
117    } else {
118      infeasible_=NULL;
119    }
120    if (rhs.weights_!=NULL) {
121      assert(model_);
122      int number = model_->numberRows()+model_->numberColumns();
123      weights_= new double[number];
124      ClpDisjointCopyN(rhs.weights_,number,weights_);
125      savedWeights_= new double[number];
126      ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
127      if (!mode_) {
128        reference_ = new unsigned int[(number+31)>>5];
129        memcpy(reference_,rhs.reference_,
130               ((number+31)>>5)*sizeof(unsigned int));
131      }
132    } else {
133      weights_=NULL;
134    }
135    if (rhs.alternateWeights_!=NULL) {
136      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
137    } else {
138      alternateWeights_=NULL;
139    }
140  }
141  return *this;
142}
143
144#define TRY_NORM 1.0e-4
145#define ADD_ONE 1.0
146// Returns pivot column, -1 if none
147int 
148ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
149                                    CoinIndexedVector * spareRow1,
150                                    CoinIndexedVector * spareRow2,
151                                    CoinIndexedVector * spareColumn1,
152                                    CoinIndexedVector * spareColumn2)
153{
154  assert(model_);
155  int iSection,j;
156  int number;
157  int * index;
158  double * updateBy;
159  double * reducedCost;
160  // dj could be very small (or even zero - take care)
161  double dj = model_->dualIn();
162  double * infeas = infeasible_->denseVector();
163  double tolerance=model_->currentDualTolerance();
164  int pivotRow = model_->pivotRow();
165
166  int anyUpdates;
167
168  if (updates->getNumElements()) {
169    // would have to have two goes for devex, three for steepest
170    anyUpdates=2;
171    // add in pivot contribution
172    if (pivotRow>=0) 
173      updates->add(pivotRow,-dj);
174  } else if (pivotRow>=0) {
175    if (fabs(dj)>1.0e-15) {
176      // some dj
177      updates->insert(pivotRow,-dj);
178      if (fabs(dj)>1.0e-6) {
179        // reasonable size
180        anyUpdates=1;
181      } else {
182        // too small
183        anyUpdates=2;
184      }
185    } else {
186      // zero dj
187      anyUpdates=-1;
188    }
189  } else if (pivotSequence_>=0){
190    // just after re-factorization
191    anyUpdates=-1;
192  } else {
193    // sub flip - nothing to do
194    anyUpdates=0;
195  }
196
197  if (anyUpdates>0) {
198    model_->factorization()->updateColumnTranspose(spareRow2,updates);
199    // put row of tableau in rowArray and columnArray
200    model_->clpMatrix()->transposeTimes(model_,-1.0,
201                                        updates,spareColumn2,spareColumn1);
202    for (iSection=0;iSection<2;iSection++) {
203     
204      reducedCost=model_->djRegion(iSection);
205      int addSequence;
206     
207      if (!iSection) {
208        number = updates->getNumElements();
209        index = updates->getIndices();
210        updateBy = updates->denseVector();
211        addSequence = model_->numberColumns();
212      } else {
213        number = spareColumn1->getNumElements();
214        index = spareColumn1->getIndices();
215        updateBy = spareColumn1->denseVector();
216        addSequence = 0;
217      }
218
219      for (j=0;j<number;j++) {
220        int iSequence = index[j];
221        double value = reducedCost[iSequence];
222        value -= updateBy[iSequence];
223        reducedCost[iSequence] = value;
224        if (model_->fixed(iSequence+addSequence))
225          continue;
226        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
227
228        switch(status) {
229         
230        case ClpSimplex::basic:
231          infeasible_->zero(iSequence+addSequence);
232          break;
233        case ClpSimplex::isFree:
234        case ClpSimplex::superBasic:
235          if (fabs(value)>tolerance) {
236            // we are going to bias towards free
237            value *= 10.0;
238            // store square in list
239            if (infeas[iSequence+addSequence])
240              infeas[iSequence+addSequence] = value*value; // already there
241            else
242              infeasible_->quickAdd(iSequence+addSequence,value*value);
243          } else {
244            infeasible_->zero(iSequence+addSequence);
245          }
246          break;
247        case ClpSimplex::atUpperBound:
248          if (value>tolerance) {
249            // store square in list
250            if (infeas[iSequence+addSequence])
251              infeas[iSequence+addSequence] = value*value; // already there
252            else
253              infeasible_->quickAdd(iSequence+addSequence,value*value);
254          } else {
255            infeasible_->zero(iSequence+addSequence);
256          }
257          break;
258        case ClpSimplex::atLowerBound:
259          if (value<-tolerance) {
260            // store square in list
261            if (infeas[iSequence+addSequence])
262              infeas[iSequence+addSequence] = value*value; // already there
263            else
264              infeasible_->quickAdd(iSequence+addSequence,value*value);
265          } else {
266            infeasible_->zero(iSequence+addSequence);
267          }
268        }
269      }
270    }
271    if (anyUpdates==2) {
272      // we can zero out as will have to get pivot row
273      updates->clear();
274      spareColumn1->clear();
275    }
276    if (pivotRow>=0) {
277      // make sure infeasibility on incoming is 0.0
278      int sequenceIn = model_->sequenceIn();
279      infeasible_->zero(sequenceIn);
280    }
281  }
282  // make sure outgoing from last iteration okay
283  int sequenceOut = model_->sequenceOut();
284  if (sequenceOut>=0) {
285    if (!model_->fixed(sequenceOut)) {
286      ClpSimplex::Status status = model_->getStatus(sequenceOut);
287      double value = model_->reducedCost(sequenceOut);
288     
289      switch(status) {
290
291      case ClpSimplex::basic:
292        break;
293      case ClpSimplex::isFree:
294      case ClpSimplex::superBasic:
295        if (fabs(value)>tolerance) { 
296          // we are going to bias towards free
297          value *= 10.0;
298          // store square in list
299          if (infeas[sequenceOut])
300            infeas[sequenceOut] = value*value; // already there
301          else
302            infeasible_->quickAdd(sequenceOut,value*value);
303        } else {
304          infeasible_->zero(sequenceOut);
305        }
306        break;
307      case ClpSimplex::atUpperBound:
308        if (value>tolerance) {
309          // store square in list
310          if (infeas[sequenceOut])
311            infeas[sequenceOut] = value*value; // already there
312          else
313            infeasible_->quickAdd(sequenceOut,value*value);
314        } else {
315          infeasible_->zero(sequenceOut);
316        }
317        break;
318      case ClpSimplex::atLowerBound:
319        if (value<-tolerance) {
320          // store square in list
321          if (infeas[sequenceOut])
322            infeas[sequenceOut] = value*value; // already there
323          else
324            infeasible_->quickAdd(sequenceOut,value*value);
325        } else {
326          infeasible_->zero(sequenceOut);
327        }
328      }
329    }
330  }
331
332#ifdef RECOMPUTE_ALL_DJS
333  for (iSection=0;iSection<2;iSection++) {
334     
335    reducedCost=model_->djRegion(iSection);
336    int addSequence;
337    int iSequence;
338     
339    if (!iSection) {
340      number = model_->numberRows();
341      addSequence = model_->numberColumns();
342    } else {
343      number = model_->numberColumns();
344      addSequence = 0;
345    }
346
347    for (iSequence=0;iSequence<number;iSequence++) {
348      double value = reducedCost[iSequence];
349      if (model_->fixed(iSequence+addSequence))
350        continue;
351      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
352     
353      switch(status) {
354
355      case ClpSimplex::basic:
356        infeasible_->zero(iSequence+addSequence);
357        break;
358      case ClpSimplex::isFree:
359        if (fabs(value)>tolerance) {
360          // we are going to bias towards free
361          value *= 10.0;
362          // store square in list
363          if (infeas[iSequence+addSequence])
364            infeas[iSequence+addSequence] = value*value; // already there
365          else
366            infeasible_->quickAdd(iSequence+addSequence,value*value);
367          } else {
368            infeasible_->zero(iSequence+addSequence);
369          }
370        break;
371      case ClpSimplex::atUpperBound:
372        if (value>tolerance) {
373            // store square in list
374          if (infeas[iSequence+addSequence])
375            infeas[iSequence+addSequence] = value*value; // already there
376          else
377            infeasible_->quickAdd(iSequence+addSequence,value*value);
378        } else {
379          infeasible_->zero(iSequence+addSequence);
380        }
381        break;
382      case ClpSimplex::atLowerBound:
383        if (value<-tolerance) {
384          // store square in list
385          if (infeas[iSequence+addSequence])
386            infeas[iSequence+addSequence] = value*value; // already there
387          else
388            infeasible_->quickAdd(iSequence+addSequence,value*value);
389        } else {
390          infeasible_->zero(iSequence+addSequence);
391        }
392      }
393    }
394  }
395#endif
396  // for weights update we use pivotSequence
397  pivotRow = pivotSequence_;
398  // unset in case sub flip
399  pivotSequence_=-1;
400  if (pivotRow>=0) {
401    // make sure infeasibility on incoming is 0.0
402    const int * pivotVariable = model_->pivotVariable();
403    int sequenceIn = pivotVariable[pivotRow];
404    infeasible_->zero(sequenceIn);
405    // and we can see if reference
406    double referenceIn=0.0;
407    if (!mode_&&reference(sequenceIn))
408      referenceIn=1.0;
409    // save outgoing weight round update
410    double outgoingWeight=weights_[sequenceOut];
411    // update weights
412    if (anyUpdates!=1) {
413      updates->setNumElements(0);
414      spareColumn1->setNumElements(0);
415      // might as well set dj to 1
416      dj = 1.0;
417      updates->insert(pivotRow,-dj);
418      model_->factorization()->updateColumnTranspose(spareRow2,updates);
419      // put row of tableau in rowArray and columnArray
420      model_->clpMatrix()->transposeTimes(model_,-1.0,
421                                          updates,spareColumn2,spareColumn1);
422    }
423    // now update weight update array
424    model_->factorization()->updateColumnTranspose(spareRow2,
425                                                   alternateWeights_);
426    double * other = alternateWeights_->denseVector();
427    double * weight;
428    int numberColumns = model_->numberColumns();
429    double scaleFactor = -1.0/dj; // as formula is with 1.0
430    // rows
431    number = updates->getNumElements();
432    index = updates->getIndices();
433    updateBy = updates->denseVector();
434    weight = weights_+numberColumns;
435   
436    for (j=0;j<number;j++) {
437      int iSequence = index[j];
438      double thisWeight = weight[iSequence];
439      // row has -1
440      double pivot = updateBy[iSequence]*scaleFactor;
441      updateBy[iSequence]=0.0;
442      double modification = other[iSequence];
443      double pivotSquared = pivot * pivot;
444     
445      thisWeight += pivotSquared * devex_ + pivot * modification;
446      if (thisWeight<TRY_NORM) {
447        if (mode_) {
448          // steepest
449          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
450        } else {
451          // exact
452          thisWeight = referenceIn*pivotSquared;
453          if (reference(iSequence+numberColumns))
454            thisWeight += 1.0;
455          thisWeight = max(thisWeight,TRY_NORM);
456        }
457      }
458      weight[iSequence] = thisWeight;
459    }
460   
461    // columns
462    weight = weights_;
463   
464    scaleFactor = -scaleFactor;
465   
466    number = spareColumn1->getNumElements();
467    index = spareColumn1->getIndices();
468    updateBy = spareColumn1->denseVector();
469    // get subset which have nonzero tableau elements
470    model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
471                                              spareColumn1,
472                                              spareColumn2);
473    double * updateBy2 = spareColumn2->denseVector();
474    for (j=0;j<number;j++) {
475      int iSequence = index[j];
476      double thisWeight = weight[iSequence];
477      double pivot = updateBy[iSequence]*scaleFactor;
478      updateBy[iSequence]=0.0;
479      double modification = updateBy2[iSequence];
480      updateBy2[iSequence]=0.0;
481      double pivotSquared = pivot * pivot;
482     
483      thisWeight += pivotSquared * devex_ + pivot * modification;
484      if (thisWeight<TRY_NORM) {
485        if (mode_) {
486          // steepest
487          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
488        } else {
489          // exact
490          thisWeight = referenceIn*pivotSquared;
491          if (reference(iSequence))
492            thisWeight += 1.0;
493          thisWeight = max(thisWeight,TRY_NORM);
494        }
495      }
496      weight[iSequence] = thisWeight;
497    }
498    // restore outgoing weight
499    weights_[sequenceOut]=outgoingWeight;
500    alternateWeights_->clear();
501    spareColumn2->setNumElements(0);
502#ifdef SOME_DEBUG_1
503#if 1
504    // check for accuracy
505    int iCheck=229;
506    printf("weight for iCheck is %g\n",weights_[iCheck]);
507    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
508        !model_->fixed(iCheck))
509      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
510#else
511    // check for accuracy
512    number = updates->getNumElements();
513    index = updates->getIndices();
514    for (j=0;j<number;j++) {
515      if (model_->getStatus(index[j])!=ClpSimplex::basic&&
516          !model_->fixed(index[j]))
517        checkAccuracy(index[j],1.0e-1,updates,spareRow2);
518    }
519#endif
520#endif
521    updates->setNumElements(0);
522    spareColumn1->setNumElements(0);
523  }
524  infeasible_->stopQuickAdd();
525
526
527  // update of duals finished - now do pricing
528
529
530  double bestDj = 1.0e-30;
531  int bestSequence=-1;
532
533  int i,iSequence;
534  index = infeasible_->getIndices();
535  number = infeasible_->getNumElements();
536  if(model_->numberIterations()<model_->lastBadIteration()+200) {
537    // we can't really trust infeasibilities if there is dual error
538    double checkTolerance = 1.0e-8;
539    if (!model_->factorization()->pivots())
540      checkTolerance = 1.0e-6;
541    if (model_->largestDualError()>checkTolerance)
542      tolerance *= model_->largestDualError()/checkTolerance;
543  }
544  // stop last one coming immediately
545  double saveOutInfeasibility=0.0;
546  if (sequenceOut>=0) {
547    saveOutInfeasibility = infeas[sequenceOut];
548    infeas[sequenceOut]=0.0;
549  }
550  tolerance *= tolerance; // as we are using squares
551  for (i=0;i<number;i++) {
552    iSequence = index[i];
553    double value = infeas[iSequence];
554    double weight = weights_[iSequence];
555    /*if (model_->numberIterations()%100==0)
556      printf("%d inf %g wt %g\n",iSequence,value,weight);*/
557    //weight=1.0;
558    if (value>bestDj*weight&&value>tolerance) {
559      // check flagged variable
560      if (!model_->flagged(iSequence)) {
561        /*if (model_->numberIterations()%100==0)
562          printf("chosen %g => %g\n",bestDj,value/weight);*/
563        bestDj=value/weight;
564        bestSequence = iSequence;
565      }
566    }
567  }
568  if (sequenceOut>=0) {
569    infeas[sequenceOut]=saveOutInfeasibility;
570  }
571  /*if (model_->numberIterations()%100==0)
572    printf("%d best %g\n",bestSequence,bestDj);*/
573 
574#ifdef CLP_DEBUG
575  if (bestSequence>=0) {
576    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
577      assert(model_->reducedCost(bestSequence)<0.0);
578    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
579      assert(model_->reducedCost(bestSequence)>0.0);
580  }
581#endif
582  return bestSequence;
583}
584/*
585   1) before factorization
586   2) after factorization
587   3) just redo infeasibilities
588   4) restore weights
589*/
590void 
591ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
592{
593  // alternateWeights_ is defined as indexed but is treated oddly
594  // at times
595  model_ = model;
596  int numberRows = model_->numberRows();
597  int numberColumns = model_->numberColumns();
598  const int * pivotVariable = model_->pivotVariable();
599  if (mode==1&&weights_) {
600    if (pivotSequence_>=0) {
601      // save pivot order
602      memcpy(alternateWeights_->getIndices(),pivotVariable,
603             numberRows*sizeof(int));
604      // change from pivot row number to sequence number
605      pivotSequence_=pivotVariable[pivotSequence_];
606    } else {
607      alternateWeights_->clear();
608    }
609    state_=1;
610  } else if (mode==2||mode==4||mode==5) {
611    // restore
612    if (!weights_||state_==-1||mode==5) {
613      // initialize weights
614      delete [] weights_;
615      delete alternateWeights_;
616      weights_ = new double[numberRows+numberColumns];
617      alternateWeights_ = new CoinIndexedVector();
618      // enough space so can use it for factorization
619      alternateWeights_->reserve(numberRows+
620                                 model_->factorization()->maximumPivots());
621      initializeWeights();
622      // create saved weights
623      savedWeights_ = new double[numberRows+numberColumns];
624      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
625             sizeof(double));
626      savedPivotSequence_=-2;
627      savedSequenceOut_ = -2;
628     
629    } else {
630      if (mode!=4) {
631        // save
632        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
633               sizeof(double));
634        savedPivotSequence_= pivotSequence_;
635        savedSequenceOut_ = model_->sequenceOut();
636      } else {
637        // restore
638        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
639               sizeof(double));
640        // was - but I think should not be
641        //pivotSequence_= savedPivotSequence_;
642        //model_->setSequenceOut(savedSequenceOut_);
643        pivotSequence_= -1;
644        model_->setSequenceOut(-1); 
645      }
646    }
647    state_=0;
648    // set up infeasibilities
649    if (!infeasible_) {
650      infeasible_ = new CoinIndexedVector();
651      infeasible_->reserve(numberColumns+numberRows);
652    }
653  }
654  if (mode>=2&&mode!=5) {
655    if (mode!=3&&pivotSequence_>=0) {
656      // restore pivot row
657      int iRow;
658      // permute alternateWeights
659      double * temp = new double[numberRows+numberColumns];
660      double * work = alternateWeights_->denseVector();
661      int * oldPivotOrder = alternateWeights_->getIndices();
662      for (iRow=0;iRow<numberRows;iRow++) {
663        int iPivot=oldPivotOrder[iRow];
664        temp[iPivot]=work[iRow];
665      }
666      int number=0;
667      int found=-1;
668      int * which = oldPivotOrder;
669      // find pivot row and re-create alternateWeights
670      for (iRow=0;iRow<numberRows;iRow++) {
671        int iPivot=pivotVariable[iRow];
672        if (iPivot==pivotSequence_) 
673          found=iRow;
674        work[iRow]=temp[iPivot];
675        if (work[iRow])
676          which[number++]=iRow;
677      }
678      alternateWeights_->setNumElements(number);
679#ifdef CLP_DEBUG
680      // Can happen but I should clean up
681      assert(found>=0);
682#endif
683      pivotSequence_ = found;
684      delete [] temp;
685    }
686    infeasible_->clear();
687    double tolerance=model_->currentDualTolerance();
688    int number = model_->numberRows() + model_->numberColumns();
689    int iSequence;
690
691    double * reducedCost = model_->djRegion();
692     
693    for (iSequence=0;iSequence<number;iSequence++) {
694      if (model_->fixed(iSequence))
695        continue;
696      double value = reducedCost[iSequence];
697      ClpSimplex::Status status = model_->getStatus(iSequence);
698     
699      switch(status) {
700
701      case ClpSimplex::basic:
702        break;
703      case ClpSimplex::isFree:
704      case ClpSimplex::superBasic:
705        if (fabs(value)>tolerance) { 
706          // store square in list
707          infeasible_->quickAdd(iSequence,value*value);
708        }
709        break;
710      case ClpSimplex::atUpperBound:
711        if (value>tolerance) {
712          infeasible_->quickAdd(iSequence,value*value);
713        }
714        break;
715      case ClpSimplex::atLowerBound:
716        if (value<-tolerance) {
717          infeasible_->quickAdd(iSequence,value*value);
718        }
719      }
720    }
721    infeasible_->stopQuickAdd();
722  }
723}
724// Gets rid of last update
725void 
726ClpPrimalColumnSteepest::unrollWeights()
727{
728  double * saved = alternateWeights_->denseVector();
729  int number = alternateWeights_->getNumElements();
730  int * which = alternateWeights_->getIndices();
731  int i;
732  for (i=0;i<number;i++) {
733    int iRow = which[i];
734    weights_[iRow]=saved[iRow];
735    saved[iRow]=0.0;
736  }
737  alternateWeights_->setNumElements(0);
738}
739 
740//-------------------------------------------------------------------
741// Clone
742//-------------------------------------------------------------------
743ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
744{
745  if (CopyData) {
746    return new ClpPrimalColumnSteepest(*this);
747  } else {
748    return new ClpPrimalColumnSteepest();
749  }
750}
751void
752ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
753{
754  int number=input->getNumElements();
755  int * which = input->getIndices();
756  double * work = input->denseVector();
757  int newNumber = 0;
758  int * newWhich = alternateWeights_->getIndices();
759  double * newWork = alternateWeights_->denseVector();
760  int i;
761  int sequenceIn = model_->sequenceIn();
762  int sequenceOut = model_->sequenceOut();
763  const int * pivotVariable = model_->pivotVariable();
764
765  int pivotRow = model_->pivotRow();
766  pivotSequence_ = pivotRow;
767
768  devex_ =0.0;
769
770  if (pivotRow>=0) {
771    if (mode_) {
772      for (i=0;i<number;i++) {
773        int iRow = which[i];
774        devex_ += work[iRow]*work[iRow];
775        newWork[iRow] = -2.0*work[iRow];
776      }
777      newWork[pivotRow] = -2.0*max(devex_,0.0);
778      devex_+=ADD_ONE;
779      weights_[sequenceOut]=1.0+ADD_ONE;
780      memcpy(newWhich,which,number*sizeof(int));
781      alternateWeights_->setNumElements(number);
782    } else {
783      for (i=0;i<number;i++) {
784        int iRow = which[i];
785        int iPivot = pivotVariable[iRow];
786        if (reference(iPivot)) {
787          devex_ += work[iRow]*work[iRow];
788          newWork[iRow] = -2.0*work[iRow];
789          newWhich[newNumber++]=iRow;
790        }
791      }
792      if (!newWork[pivotRow])
793        newWhich[newNumber++]=pivotRow; // add if not already in
794      newWork[pivotRow] = -2.0*max(devex_,0.0);
795      if (reference(sequenceIn)) {
796        devex_+=1.0;
797      } else {
798      }
799      if (reference(sequenceOut)) {
800        weights_[sequenceOut]=1.0+1.0;
801      } else {
802        weights_[sequenceOut]=1.0;
803      }
804      alternateWeights_->setNumElements(newNumber);
805    }
806  } else {
807    if (mode_) {
808      for (i=0;i<number;i++) {
809        int iRow = which[i];
810        devex_ += work[iRow]*work[iRow];
811      }
812      devex_ += ADD_ONE;
813    } else {
814      for (i=0;i<number;i++) {
815        int iRow = which[i];
816        int iPivot = pivotVariable[iRow];
817        if (reference(iPivot)) {
818          devex_ += work[iRow]*work[iRow];
819        }
820      }
821      if (reference(sequenceIn)) 
822        devex_+=1.0;
823    }
824  }
825  double oldDevex = weights_[sequenceIn];
826#ifdef CLP_DEBUG
827  if ((model_->messageHandler()->logLevel()&32))
828    printf("old weight %g, new %g\n",oldDevex,devex_);
829#endif
830  double check = max(devex_,oldDevex);;
831  weights_[sequenceIn] = devex_;
832  if ( fabs ( devex_ - oldDevex ) > 1.0e-1 * check ) {
833#ifdef CLP_DEBUG
834    if ((model_->messageHandler()->logLevel()&48)==16)
835      printf("old weight %g, new %g\n",oldDevex,devex_);
836#endif
837    if ( fabs ( devex_ - oldDevex ) > 1.0e5 * check ) {
838      // need to redo
839      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
840                                        *model_->messagesPointer())
841                                          <<oldDevex<<devex_
842                                          <<CoinMessageEol;
843      initializeWeights();
844    }
845  }
846  if (pivotRow>=0) {
847    // set outgoing weight here
848    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
849  }
850}
851// Checks accuracy - just for debug
852void 
853ClpPrimalColumnSteepest::checkAccuracy(int sequence,
854                                       double relativeTolerance,
855                                       CoinIndexedVector * rowArray1,
856                                       CoinIndexedVector * rowArray2)
857{
858  model_->unpack(rowArray1,sequence);
859  model_->factorization()->updateColumn(rowArray2,rowArray1);
860  int number=rowArray1->getNumElements();
861  int * which = rowArray1->getIndices();
862  double * work = rowArray1->denseVector();
863  const int * pivotVariable = model_->pivotVariable();
864 
865  double devex =0.0;
866  int i;
867
868  if (mode_) {
869    for (i=0;i<number;i++) {
870      int iRow = which[i];
871      devex += work[iRow]*work[iRow];
872      work[iRow]=0.0;
873    }
874    devex += ADD_ONE;
875  } else {
876    for (i=0;i<number;i++) {
877      int iRow = which[i];
878      int iPivot = pivotVariable[iRow];
879      if (reference(iPivot)) {
880        devex += work[iRow]*work[iRow];
881      }
882      work[iRow]=0.0;
883    }
884    if (reference(sequence)) 
885      devex+=1.0;
886  }
887
888  double oldDevex = weights_[sequence];
889  double check = max(devex,oldDevex);;
890  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
891    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
892    // update so won't print again
893    weights_[sequence]=devex;
894  }
895  rowArray1->setNumElements(0);
896}
897
898// Initialize weights
899void 
900ClpPrimalColumnSteepest::initializeWeights()
901{
902  int numberRows = model_->numberRows();
903  int numberColumns = model_->numberColumns();
904  int number = numberRows + numberColumns;
905  int iSequence;
906  if (!mode_) {
907    // initialize to 1.0
908    // and set reference framework
909    if (!reference_) {
910      int nWords = (number+31)>>5;
911      reference_ = new unsigned int[nWords];
912      assert (sizeof(unsigned int)==4);
913      // tiny overhead to zero out (stops valgrind complaining)
914      memset(reference_,0,nWords*sizeof(int));
915    }
916   
917    for (iSequence=0;iSequence<number;iSequence++) {
918      weights_[iSequence]=1.0;
919      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
920        setReference(iSequence,false);
921      } else {
922        setReference(iSequence,true);
923      }
924    }
925  } else {
926    CoinIndexedVector * temp = new CoinIndexedVector();
927    temp->reserve(numberRows+
928                  model_->factorization()->maximumPivots());
929    double * array = alternateWeights_->denseVector();
930    int * which = alternateWeights_->getIndices();
931     
932    for (iSequence=0;iSequence<number;iSequence++) {
933      weights_[iSequence]=1.0+ADD_ONE;
934      if (model_->fixed(iSequence))
935        continue;
936      if (model_->getStatus(iSequence)!=ClpSimplex::basic) {
937        model_->unpack(alternateWeights_,iSequence);
938        double value=ADD_ONE;
939        model_->factorization()->updateColumn(temp,alternateWeights_);
940        int number = alternateWeights_->getNumElements();
941        int j;
942        for (j=0;j<number;j++) {
943          int iRow=which[j];
944          value+=array[iRow]*array[iRow];
945          array[iRow]=0.0;
946        }
947        alternateWeights_->setNumElements(0);
948        weights_[iSequence] = value;
949      }
950    }
951    delete temp;
952  }
953}
954// Gets rid of all arrays
955void 
956ClpPrimalColumnSteepest::clearArrays()
957{
958  delete [] weights_;
959  weights_=NULL;
960  delete infeasible_;
961  infeasible_ = NULL;
962  delete alternateWeights_;
963  alternateWeights_ = NULL;
964  delete [] savedWeights_;
965  savedWeights_ = NULL;
966  delete [] reference_;
967  reference_ = NULL;
968  pivotSequence_=-1;
969  state_ = -1;
970  savedPivotSequence_ = -1;
971  savedSequenceOut_ = -1; 
972  devex_ = 0.0;
973}
Note: See TracBrowser for help on using the repository browser.