source: trunk/ClpPrimalColumnSteepest.cpp @ 114

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

Re-order variables

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