source: tags/root-devel-1/ClpPrimalColumnSteepest.cpp @ 778

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

Adding Clp to development branch

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