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

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

Changes to make more reliable if problem size changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.1 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    CoinDisjointCopyN(rhs.weights_,number,weights_);
64    savedWeights_= new double[number];
65    CoinDisjointCopyN(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      CoinDisjointCopyN(rhs.weights_,number,weights_);
125      savedWeights_= new double[number];
126      CoinDisjointCopyN(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  tolerance *= tolerance; // as we are using squares
545  for (i=0;i<number;i++) {
546    iSequence = index[i];
547    double value = infeas[iSequence];
548    double weight = weights_[iSequence];
549    /*if (model_->numberIterations()%100==0)
550      printf("%d inf %g wt %g\n",iSequence,value,weight);*/
551    //weight=1.0;
552    if (value>bestDj*weight&&value>tolerance) {
553      // check flagged variable
554      if (!model_->flagged(iSequence)) {
555        /*if (model_->numberIterations()%100==0)
556          printf("chosen %g => %g\n",bestDj,value/weight);*/
557        bestDj=value/weight;
558        bestSequence = iSequence;
559      }
560    }
561  }
562  /*if (model_->numberIterations()%100==0)
563    printf("%d best %g\n",bestSequence,bestDj);*/
564 
565#ifdef CLP_DEBUG
566  if (bestSequence>=0) {
567    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
568      assert(model_->reducedCost(bestSequence)<0.0);
569    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
570      assert(model_->reducedCost(bestSequence)>0.0);
571  }
572#endif
573  return bestSequence;
574}
575/*
576   1) before factorization
577   2) after factorization
578   3) just redo infeasibilities
579   4) restore weights
580*/
581void 
582ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
583{
584  // alternateWeights_ is defined as indexed but is treated oddly
585  // at times
586  model_ = model;
587  int numberRows = model_->numberRows();
588  int numberColumns = model_->numberColumns();
589  const int * pivotVariable = model_->pivotVariable();
590  if (mode==1&&weights_) {
591    if (pivotSequence_>=0) {
592      // save pivot order
593      memcpy(alternateWeights_->getIndices(),pivotVariable,
594             numberRows*sizeof(int));
595      // change from pivot row number to sequence number
596      pivotSequence_=pivotVariable[pivotSequence_];
597    } else {
598      alternateWeights_->clear();
599    }
600    state_=1;
601  } else if (mode==2||mode==4||mode==5) {
602    // restore
603    if (!weights_||state_==-1||mode==5) {
604      // initialize weights
605      delete [] weights_;
606      delete alternateWeights_;
607      weights_ = new double[numberRows+numberColumns];
608      alternateWeights_ = new CoinIndexedVector();
609      // enough space so can use it for factorization
610      alternateWeights_->reserve(numberRows+
611                                 model_->factorization()->maximumPivots());
612      initializeWeights();
613      // create saved weights
614      savedWeights_ = new double[numberRows+numberColumns];
615      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
616             sizeof(double));
617      savedPivotSequence_=-2;
618      savedSequenceOut_ = -2;
619     
620    } else {
621      if (mode!=4) {
622        // save
623        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
624               sizeof(double));
625        savedPivotSequence_= pivotSequence_;
626        savedSequenceOut_ = model_->sequenceOut();
627      } else {
628        // restore
629        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
630               sizeof(double));
631        pivotSequence_= savedPivotSequence_;
632        model_->setSequenceOut(savedSequenceOut_); 
633      }
634    }
635    state_=0;
636    // set up infeasibilities
637    if (!infeasible_) {
638      infeasible_ = new CoinIndexedVector();
639      infeasible_->reserve(numberColumns+numberRows);
640    }
641  }
642  if (mode>=2&&mode!=5) {
643    if (mode!=3&&pivotSequence_>=0) {
644      // restore pivot row
645      int iRow;
646      // permute alternateWeights
647      double * temp = new double[numberRows+numberColumns];
648      double * work = alternateWeights_->denseVector();
649      int * oldPivotOrder = alternateWeights_->getIndices();
650      for (iRow=0;iRow<numberRows;iRow++) {
651        int iPivot=oldPivotOrder[iRow];
652        temp[iPivot]=work[iRow];
653      }
654      int number=0;
655      int found=-1;
656      int * which = oldPivotOrder;
657      // find pivot row and re-create alternateWeights
658      for (iRow=0;iRow<numberRows;iRow++) {
659        int iPivot=pivotVariable[iRow];
660        if (iPivot==pivotSequence_) 
661          found=iRow;
662        work[iRow]=temp[iPivot];
663        if (work[iRow])
664          which[number++]=iRow;
665      }
666      alternateWeights_->setNumElements(number);
667#ifdef CLP_DEBUG
668      // Can happen but I should clean up
669      assert(found>=0);
670#endif
671      pivotSequence_ = found;
672      delete [] temp;
673    }
674    infeasible_->clear();
675    double tolerance=model_->currentDualTolerance();
676    int number = model_->numberRows() + model_->numberColumns();
677    int iSequence;
678
679    double * reducedCost = model_->djRegion();
680     
681    for (iSequence=0;iSequence<number;iSequence++) {
682      if (model_->fixed(iSequence))
683        continue;
684      double value = reducedCost[iSequence];
685      ClpSimplex::Status status = model_->getStatus(iSequence);
686     
687      switch(status) {
688
689      case ClpSimplex::basic:
690        break;
691      case ClpSimplex::isFree:
692      case ClpSimplex::superBasic:
693        if (fabs(value)>tolerance) { 
694          // store square in list
695          infeasible_->quickAdd(iSequence,value*value);
696        }
697        break;
698      case ClpSimplex::atUpperBound:
699        if (value>tolerance) {
700          infeasible_->quickAdd(iSequence,value*value);
701        }
702        break;
703      case ClpSimplex::atLowerBound:
704        if (value<-tolerance) {
705          infeasible_->quickAdd(iSequence,value*value);
706        }
707      }
708    }
709    infeasible_->stopQuickAdd();
710  }
711}
712// Gets rid of last update
713void 
714ClpPrimalColumnSteepest::unrollWeights()
715{
716  double * saved = alternateWeights_->denseVector();
717  int number = alternateWeights_->getNumElements();
718  int * which = alternateWeights_->getIndices();
719  int i;
720  for (i=0;i<number;i++) {
721    int iRow = which[i];
722    weights_[iRow]=saved[iRow];
723    saved[iRow]=0.0;
724  }
725  alternateWeights_->setNumElements(0);
726}
727 
728//-------------------------------------------------------------------
729// Clone
730//-------------------------------------------------------------------
731ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
732{
733  if (CopyData) {
734    return new ClpPrimalColumnSteepest(*this);
735  } else {
736    return new ClpPrimalColumnSteepest();
737  }
738}
739void
740ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
741{
742  int number=input->getNumElements();
743  int * which = input->getIndices();
744  double * work = input->denseVector();
745  int newNumber = 0;
746  int * newWhich = alternateWeights_->getIndices();
747  double * newWork = alternateWeights_->denseVector();
748  int i;
749  int sequenceIn = model_->sequenceIn();
750  int sequenceOut = model_->sequenceOut();
751  const int * pivotVariable = model_->pivotVariable();
752
753  int pivotRow = model_->pivotRow();
754  pivotSequence_ = pivotRow;
755
756  devex_ =0.0;
757
758  if (pivotRow>=0) {
759    if (mode_) {
760      for (i=0;i<number;i++) {
761        int iRow = which[i];
762        devex_ += work[iRow]*work[iRow];
763        newWork[iRow] = -2.0*work[iRow];
764      }
765      newWork[pivotRow] = -2.0*max(devex_,0.0);
766      devex_+=ADD_ONE;
767      weights_[sequenceOut]=1.0+ADD_ONE;
768      memcpy(newWhich,which,number*sizeof(int));
769      alternateWeights_->setNumElements(number);
770    } else {
771      for (i=0;i<number;i++) {
772        int iRow = which[i];
773        int iPivot = pivotVariable[iRow];
774        if (reference(iPivot)) {
775          devex_ += work[iRow]*work[iRow];
776          newWork[iRow] = -2.0*work[iRow];
777          newWhich[newNumber++]=iRow;
778        }
779      }
780      if (!newWork[pivotRow])
781        newWhich[newNumber++]=pivotRow; // add if not already in
782      newWork[pivotRow] = -2.0*max(devex_,0.0);
783      if (reference(sequenceIn)) {
784        devex_+=1.0;
785      } else {
786      }
787      if (reference(sequenceOut)) {
788        weights_[sequenceOut]=1.0+1.0;
789      } else {
790        weights_[sequenceOut]=1.0;
791      }
792      alternateWeights_->setNumElements(newNumber);
793    }
794  } else {
795    if (mode_) {
796      for (i=0;i<number;i++) {
797        int iRow = which[i];
798        devex_ += work[iRow]*work[iRow];
799      }
800      devex_ += ADD_ONE;
801    } else {
802      for (i=0;i<number;i++) {
803        int iRow = which[i];
804        int iPivot = pivotVariable[iRow];
805        if (reference(iPivot)) {
806          devex_ += work[iRow]*work[iRow];
807        }
808      }
809      if (reference(sequenceIn)) 
810        devex_+=1.0;
811    }
812  }
813  double oldDevex = weights_[sequenceIn];
814#ifdef CLP_DEBUG
815  if ((model_->messageHandler()->logLevel()&32))
816    printf("old weight %g, new %g\n",oldDevex,devex_);
817#endif
818  double check = max(devex_,oldDevex);;
819  weights_[sequenceIn] = devex_;
820  if ( fabs ( devex_ - oldDevex ) > 1.0e-1 * check ) {
821#ifdef CLP_DEBUG
822    if ((model_->messageHandler()->logLevel()&48)==16)
823      printf("old weight %g, new %g\n",oldDevex,devex_);
824#endif
825    if ( fabs ( devex_ - oldDevex ) > 1.0e5 * check ) {
826      // need to redo
827      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
828                                        *model_->messagesPointer())
829                                          <<oldDevex<<devex_
830                                          <<CoinMessageEol;
831      initializeWeights();
832    }
833  }
834  if (pivotRow>=0) {
835    // set outgoing weight here
836    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
837  }
838}
839// Checks accuracy - just for debug
840void 
841ClpPrimalColumnSteepest::checkAccuracy(int sequence,
842                                       double relativeTolerance,
843                                       CoinIndexedVector * rowArray1,
844                                       CoinIndexedVector * rowArray2)
845{
846  model_->unpack(rowArray1,sequence);
847  model_->factorization()->updateColumn(rowArray2,rowArray1);
848  int number=rowArray1->getNumElements();
849  int * which = rowArray1->getIndices();
850  double * work = rowArray1->denseVector();
851  const int * pivotVariable = model_->pivotVariable();
852 
853  double devex =0.0;
854  int i;
855
856  if (mode_) {
857    for (i=0;i<number;i++) {
858      int iRow = which[i];
859      devex += work[iRow]*work[iRow];
860      work[iRow]=0.0;
861    }
862    devex += ADD_ONE;
863  } else {
864    for (i=0;i<number;i++) {
865      int iRow = which[i];
866      int iPivot = pivotVariable[iRow];
867      if (reference(iPivot)) {
868        devex += work[iRow]*work[iRow];
869      }
870      work[iRow]=0.0;
871    }
872    if (reference(sequence)) 
873      devex+=1.0;
874  }
875
876  double oldDevex = weights_[sequence];
877  double check = max(devex,oldDevex);;
878  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
879    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
880    // update so won't print again
881    weights_[sequence]=devex;
882  }
883  rowArray1->setNumElements(0);
884}
885
886// Initialize weights
887void 
888ClpPrimalColumnSteepest::initializeWeights()
889{
890  int numberRows = model_->numberRows();
891  int numberColumns = model_->numberColumns();
892  int number = numberRows + numberColumns;
893  int iSequence;
894  if (!mode_) {
895    // initialize to 1.0
896    // and set reference framework
897    if (!reference_) {
898      int nWords = (number+31)>>5;
899      reference_ = new unsigned int[nWords];
900      assert (sizeof(unsigned int)==4);
901      // tiny overhead to zero out (stops valgrind complaining)
902      memset(reference_,0,nWords*sizeof(int));
903    }
904   
905    for (iSequence=0;iSequence<number;iSequence++) {
906      weights_[iSequence]=1.0;
907      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
908        setReference(iSequence,false);
909      } else {
910        setReference(iSequence,true);
911      }
912    }
913  } else {
914    CoinIndexedVector * temp = new CoinIndexedVector();
915    temp->reserve(numberRows+
916                  model_->factorization()->maximumPivots());
917    double * array = alternateWeights_->denseVector();
918    int * which = alternateWeights_->getIndices();
919     
920    for (iSequence=0;iSequence<number;iSequence++) {
921      weights_[iSequence]=1.0+ADD_ONE;
922      if (model_->fixed(iSequence))
923        continue;
924      if (model_->getStatus(iSequence)!=ClpSimplex::basic) {
925        model_->unpack(alternateWeights_,iSequence);
926        double value=ADD_ONE;
927        model_->factorization()->updateColumn(temp,alternateWeights_);
928        int number = alternateWeights_->getNumElements();
929        int j;
930        for (j=0;j<number;j++) {
931          int iRow=which[j];
932          value+=array[iRow]*array[iRow];
933          array[iRow]=0.0;
934        }
935        alternateWeights_->setNumElements(0);
936        weights_[iSequence] = value;
937      }
938    }
939    delete temp;
940  }
941}
942// Gets rid of all arrays
943void 
944ClpPrimalColumnSteepest::clearArrays()
945{
946  delete [] weights_;
947  weights_=NULL;
948  delete infeasible_;
949  infeasible_ = NULL;
950  delete alternateWeights_;
951  alternateWeights_ = NULL;
952  delete [] savedWeights_;
953  savedWeights_ = NULL;
954  delete [] reference_;
955  reference_ = NULL;
956}
Note: See TracBrowser for help on using the repository browser.