source: trunk/ClpPrimalColumnSteepest.cpp @ 137

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

faster Clp and fix memory leaks

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.6 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   
207    // put row of tableau in rowArray and columnArray
208    model_->clpMatrix()->transposeTimes(model_,-1.0,
209                                        updates,spareColumn2,spareColumn1);
210    for (iSection=0;iSection<2;iSection++) {
211     
212      reducedCost=model_->djRegion(iSection);
213      int addSequence;
214     
215      if (!iSection) {
216        number = updates->getNumElements();
217        index = updates->getIndices();
218        updateBy = updates->denseVector();
219        addSequence = model_->numberColumns();
220      } else {
221        number = spareColumn1->getNumElements();
222        index = spareColumn1->getIndices();
223        updateBy = spareColumn1->denseVector();
224        addSequence = 0;
225      }
226
227      for (j=0;j<number;j++) {
228        int iSequence = index[j];
229        double value = reducedCost[iSequence];
230        value -= updateBy[iSequence];
231        reducedCost[iSequence] = value;
232        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
233
234        switch(status) {
235         
236        case ClpSimplex::basic:
237          infeasible_->zero(iSequence+addSequence);
238        case ClpSimplex::isFixed:
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    ClpSimplex::Status status = model_->getStatus(sequenceOut);
293    double value = model_->reducedCost(sequenceOut);
294   
295    switch(status) {
296     
297    case ClpSimplex::basic:
298    case ClpSimplex::isFixed:
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#ifdef RECOMPUTE_ALL_DJS
339  for (iSection=0;iSection<2;iSection++) {
340     
341    reducedCost=model_->djRegion(iSection);
342    int addSequence;
343    int iSequence;
344     
345    if (!iSection) {
346      number = model_->numberRows();
347      addSequence = model_->numberColumns();
348    } else {
349      number = model_->numberColumns();
350      addSequence = 0;
351    }
352
353    for (iSequence=0;iSequence<number;iSequence++) {
354      double value = reducedCost[iSequence];
355      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
356     
357      switch(status) {
358
359      case ClpSimplex::basic:
360        infeasible_->zero(iSequence+addSequence);
361      case ClpSimplex::isFixed:
362        break;
363      case ClpSimplex::isFree:
364        if (fabs(value)>tolerance) {
365          // we are going to bias towards free (but only if reasonable)
366          value *= FREE_BIAS;
367          // store square in list
368          if (infeas[iSequence+addSequence])
369            infeas[iSequence+addSequence] = value*value; // already there
370          else
371            infeasible_->quickAdd(iSequence+addSequence,value*value);
372          } else {
373            infeasible_->zero(iSequence+addSequence);
374          }
375        break;
376      case ClpSimplex::atUpperBound:
377        if (value>tolerance) {
378            // store square in list
379          if (infeas[iSequence+addSequence])
380            infeas[iSequence+addSequence] = value*value; // already there
381          else
382            infeasible_->quickAdd(iSequence+addSequence,value*value);
383        } else {
384          infeasible_->zero(iSequence+addSequence);
385        }
386        break;
387      case ClpSimplex::atLowerBound:
388        if (value<-tolerance) {
389          // store square in list
390          if (infeas[iSequence+addSequence])
391            infeas[iSequence+addSequence] = value*value; // already there
392          else
393            infeasible_->quickAdd(iSequence+addSequence,value*value);
394        } else {
395          infeasible_->zero(iSequence+addSequence);
396        }
397      }
398    }
399  }
400#endif
401  // for weights update we use pivotSequence
402  pivotRow = pivotSequence_;
403  // unset in case sub flip
404  pivotSequence_=-1;
405  if (pivotRow>=0) {
406    // make sure infeasibility on incoming is 0.0
407    const int * pivotVariable = model_->pivotVariable();
408    int sequenceIn = pivotVariable[pivotRow];
409    infeasible_->zero(sequenceIn);
410    // and we can see if reference
411    double referenceIn=0.0;
412    if (!mode_&&reference(sequenceIn))
413      referenceIn=1.0;
414    // save outgoing weight round update
415    double outgoingWeight=weights_[sequenceOut];
416    // update weights
417    if (anyUpdates!=1) {
418      updates->setNumElements(0);
419      spareColumn1->setNumElements(0);
420      // might as well set dj to 1
421      dj = 1.0;
422      updates->insert(pivotRow,-dj);
423      model_->factorization()->updateColumnTranspose(spareRow2,updates);
424      // put row of tableau in rowArray and columnArray
425      model_->clpMatrix()->transposeTimes(model_,-1.0,
426                                          updates,spareColumn2,spareColumn1);
427    }
428    // now update weight update array
429    model_->factorization()->updateColumnTranspose(spareRow2,
430                                                   alternateWeights_);
431    double * other = alternateWeights_->denseVector();
432    double * weight;
433    int numberColumns = model_->numberColumns();
434    double scaleFactor = -1.0/dj; // as formula is with 1.0
435    // rows
436    number = updates->getNumElements();
437    index = updates->getIndices();
438    updateBy = updates->denseVector();
439    weight = weights_+numberColumns;
440   
441    for (j=0;j<number;j++) {
442      int iSequence = index[j];
443      double thisWeight = weight[iSequence];
444      // row has -1
445      double pivot = updateBy[iSequence]*scaleFactor;
446      updateBy[iSequence]=0.0;
447      double modification = other[iSequence];
448      double pivotSquared = pivot * pivot;
449     
450      thisWeight += pivotSquared * devex_ + pivot * modification;
451      if (thisWeight<TRY_NORM) {
452        if (mode_) {
453          // steepest
454          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
455        } else {
456          // exact
457          thisWeight = referenceIn*pivotSquared;
458          if (reference(iSequence+numberColumns))
459            thisWeight += 1.0;
460          thisWeight = max(thisWeight,TRY_NORM);
461        }
462      }
463      weight[iSequence] = thisWeight;
464    }
465   
466    // columns
467    weight = weights_;
468   
469    scaleFactor = -scaleFactor;
470   
471    number = spareColumn1->getNumElements();
472    index = spareColumn1->getIndices();
473    updateBy = spareColumn1->denseVector();
474    // get subset which have nonzero tableau elements
475    model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
476                                              spareColumn1,
477                                              spareColumn2);
478    double * updateBy2 = spareColumn2->denseVector();
479    for (j=0;j<number;j++) {
480      int iSequence = index[j];
481      double thisWeight = weight[iSequence];
482      double pivot = updateBy[iSequence]*scaleFactor;
483      updateBy[iSequence]=0.0;
484      double modification = updateBy2[iSequence];
485      updateBy2[iSequence]=0.0;
486      double pivotSquared = pivot * pivot;
487     
488      thisWeight += pivotSquared * devex_ + pivot * modification;
489      if (thisWeight<TRY_NORM) {
490        if (mode_) {
491          // steepest
492          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
493        } else {
494          // exact
495          thisWeight = referenceIn*pivotSquared;
496          if (reference(iSequence))
497            thisWeight += 1.0;
498          thisWeight = max(thisWeight,TRY_NORM);
499        }
500      }
501      weight[iSequence] = thisWeight;
502    }
503    // restore outgoing weight
504    weights_[sequenceOut]=outgoingWeight;
505    alternateWeights_->clear();
506    spareColumn2->setNumElements(0);
507#ifdef SOME_DEBUG_1
508    // check for accuracy
509    int iCheck=229;
510    printf("weight for iCheck is %g\n",weights_[iCheck]);
511    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
512        !model_->getStatus(iCheck)!=isFixed)
513      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
514#endif
515    updates->setNumElements(0);
516    spareColumn1->setNumElements(0);
517  }
518
519  // update of duals finished - now do pricing
520
521
522  double bestDj = 1.0e-30;
523  int bestSequence=-1;
524
525  int i,iSequence;
526  index = infeasible_->getIndices();
527  number = infeasible_->getNumElements();
528  if(model_->numberIterations()<model_->lastBadIteration()+200) {
529    // we can't really trust infeasibilities if there is dual error
530    double checkTolerance = 1.0e-8;
531    if (!model_->factorization()->pivots())
532      checkTolerance = 1.0e-6;
533    if (model_->largestDualError()>checkTolerance)
534      tolerance *= model_->largestDualError()/checkTolerance;
535  }
536  // stop last one coming immediately
537  double saveOutInfeasibility=0.0;
538  if (sequenceOut>=0) {
539    saveOutInfeasibility = infeas[sequenceOut];
540    infeas[sequenceOut]=0.0;
541  }
542  tolerance *= tolerance; // as we are using squares
543  for (i=0;i<number;i++) {
544    iSequence = index[i];
545    double value = infeas[iSequence];
546    double weight = weights_[iSequence];
547    /*if (model_->numberIterations()%100==0)
548      printf("%d inf %g wt %g\n",iSequence,value,weight);*/
549    //weight=1.0;
550    if (value>bestDj*weight&&value>tolerance) {
551      // check flagged variable
552      if (!model_->flagged(iSequence)) {
553        /*if (model_->numberIterations()%100==0)
554          printf("chosen %g => %g\n",bestDj,value/weight);*/
555        bestDj=value/weight;
556        bestSequence = iSequence;
557      }
558    }
559  }
560  if (sequenceOut>=0) {
561    infeas[sequenceOut]=saveOutInfeasibility;
562  }
563  /*if (model_->numberIterations()%100==0)
564    printf("%d best %g\n",bestSequence,bestDj);*/
565 
566#ifdef CLP_DEBUG
567  if (bestSequence>=0) {
568    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
569      assert(model_->reducedCost(bestSequence)<0.0);
570    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
571      assert(model_->reducedCost(bestSequence)>0.0);
572  }
573#endif
574  return bestSequence;
575}
576/*
577   1) before factorization
578   2) after factorization
579   3) just redo infeasibilities
580   4) restore weights
581*/
582void 
583ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
584{
585  // alternateWeights_ is defined as indexed but is treated oddly
586  // at times
587  model_ = model;
588  int numberRows = model_->numberRows();
589  int numberColumns = model_->numberColumns();
590  const int * pivotVariable = model_->pivotVariable();
591  if (mode==1&&weights_) {
592    if (pivotSequence_>=0) {
593      // save pivot order
594      memcpy(alternateWeights_->getIndices(),pivotVariable,
595             numberRows*sizeof(int));
596      // change from pivot row number to sequence number
597      pivotSequence_=pivotVariable[pivotSequence_];
598    } else {
599      alternateWeights_->clear();
600    }
601    state_=1;
602  } else if (mode==2||mode==4||mode==5) {
603    // restore
604    if (!weights_||state_==-1||mode==5) {
605      // initialize weights
606      delete [] weights_;
607      delete alternateWeights_;
608      weights_ = new double[numberRows+numberColumns];
609      alternateWeights_ = new CoinIndexedVector();
610      // enough space so can use it for factorization
611      alternateWeights_->reserve(numberRows+
612                                 model_->factorization()->maximumPivots());
613      initializeWeights();
614      // create saved weights
615      savedWeights_ = new double[numberRows+numberColumns];
616      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
617             sizeof(double));
618      savedPivotSequence_=-2;
619      savedSequenceOut_ = -2;
620     
621    } else {
622      if (mode!=4) {
623        // save
624        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
625               sizeof(double));
626        savedPivotSequence_= pivotSequence_;
627        savedSequenceOut_ = model_->sequenceOut();
628      } else {
629        // restore
630        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
631               sizeof(double));
632        // was - but I think should not be
633        //pivotSequence_= savedPivotSequence_;
634        //model_->setSequenceOut(savedSequenceOut_);
635        pivotSequence_= -1;
636        model_->setSequenceOut(-1); 
637      }
638    }
639    state_=0;
640    // set up infeasibilities
641    if (!infeasible_) {
642      infeasible_ = new CoinIndexedVector();
643      infeasible_->reserve(numberColumns+numberRows);
644    }
645  }
646  if (mode>=2&&mode!=5) {
647    if (mode!=3&&pivotSequence_>=0) {
648      // restore pivot row
649      int iRow;
650      // permute alternateWeights
651      double * temp = new double[numberRows+numberColumns];
652      double * work = alternateWeights_->denseVector();
653      int * oldPivotOrder = alternateWeights_->getIndices();
654      for (iRow=0;iRow<numberRows;iRow++) {
655        int iPivot=oldPivotOrder[iRow];
656        temp[iPivot]=work[iRow];
657      }
658      int number=0;
659      int found=-1;
660      int * which = oldPivotOrder;
661      // find pivot row and re-create alternateWeights
662      for (iRow=0;iRow<numberRows;iRow++) {
663        int iPivot=pivotVariable[iRow];
664        if (iPivot==pivotSequence_) 
665          found=iRow;
666        work[iRow]=temp[iPivot];
667        if (work[iRow])
668          which[number++]=iRow;
669      }
670      alternateWeights_->setNumElements(number);
671#ifdef CLP_DEBUG
672      // Can happen but I should clean up
673      assert(found>=0);
674#endif
675      pivotSequence_ = found;
676      delete [] temp;
677    }
678    infeasible_->clear();
679    double tolerance=model_->currentDualTolerance();
680    int number = model_->numberRows() + model_->numberColumns();
681    int iSequence;
682
683    double * reducedCost = model_->djRegion();
684     
685    for (iSequence=0;iSequence<number;iSequence++) {
686      double value = reducedCost[iSequence];
687      ClpSimplex::Status status = model_->getStatus(iSequence);
688     
689      switch(status) {
690
691      case ClpSimplex::basic:
692      case ClpSimplex::isFixed:
693        break;
694      case ClpSimplex::isFree:
695      case ClpSimplex::superBasic:
696        if (fabs(value)>1.0e2*tolerance) { 
697          // we are going to bias towards free (but only if reasonable)
698          value *= FREE_BIAS;
699          // store square in list
700          infeasible_->quickAdd(iSequence,value*value);
701        }
702        break;
703      case ClpSimplex::atUpperBound:
704        if (value>tolerance) {
705          infeasible_->quickAdd(iSequence,value*value);
706        }
707        break;
708      case ClpSimplex::atLowerBound:
709        if (value<-tolerance) {
710          infeasible_->quickAdd(iSequence,value*value);
711        }
712      }
713    }
714  }
715}
716// Gets rid of last update
717void 
718ClpPrimalColumnSteepest::unrollWeights()
719{
720  double * saved = alternateWeights_->denseVector();
721  int number = alternateWeights_->getNumElements();
722  int * which = alternateWeights_->getIndices();
723  int i;
724  for (i=0;i<number;i++) {
725    int iRow = which[i];
726    weights_[iRow]=saved[iRow];
727    saved[iRow]=0.0;
728  }
729  alternateWeights_->setNumElements(0);
730}
731 
732//-------------------------------------------------------------------
733// Clone
734//-------------------------------------------------------------------
735ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
736{
737  if (CopyData) {
738    return new ClpPrimalColumnSteepest(*this);
739  } else {
740    return new ClpPrimalColumnSteepest();
741  }
742}
743void
744ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
745{
746  int number=input->getNumElements();
747  int * which = input->getIndices();
748  double * work = input->denseVector();
749  int newNumber = 0;
750  int * newWhich = alternateWeights_->getIndices();
751  double * newWork = alternateWeights_->denseVector();
752  int i;
753  int sequenceIn = model_->sequenceIn();
754  int sequenceOut = model_->sequenceOut();
755  const int * pivotVariable = model_->pivotVariable();
756
757  int pivotRow = model_->pivotRow();
758  pivotSequence_ = pivotRow;
759
760  devex_ =0.0;
761
762  if (pivotRow>=0) {
763    if (mode_) {
764      for (i=0;i<number;i++) {
765        int iRow = which[i];
766        devex_ += work[iRow]*work[iRow];
767        newWork[iRow] = -2.0*work[iRow];
768      }
769      newWork[pivotRow] = -2.0*max(devex_,0.0);
770      devex_+=ADD_ONE;
771      weights_[sequenceOut]=1.0+ADD_ONE;
772      memcpy(newWhich,which,number*sizeof(int));
773      alternateWeights_->setNumElements(number);
774    } else {
775      for (i=0;i<number;i++) {
776        int iRow = which[i];
777        int iPivot = pivotVariable[iRow];
778        if (reference(iPivot)) {
779          devex_ += work[iRow]*work[iRow];
780          newWork[iRow] = -2.0*work[iRow];
781          newWhich[newNumber++]=iRow;
782        }
783      }
784      if (!newWork[pivotRow])
785        newWhich[newNumber++]=pivotRow; // add if not already in
786      newWork[pivotRow] = -2.0*max(devex_,0.0);
787      if (reference(sequenceIn)) {
788        devex_+=1.0;
789      } else {
790      }
791      if (reference(sequenceOut)) {
792        weights_[sequenceOut]=1.0+1.0;
793      } else {
794        weights_[sequenceOut]=1.0;
795      }
796      alternateWeights_->setNumElements(newNumber);
797    }
798  } else {
799    if (mode_) {
800      for (i=0;i<number;i++) {
801        int iRow = which[i];
802        devex_ += work[iRow]*work[iRow];
803      }
804      devex_ += ADD_ONE;
805    } else {
806      for (i=0;i<number;i++) {
807        int iRow = which[i];
808        int iPivot = pivotVariable[iRow];
809        if (reference(iPivot)) {
810          devex_ += work[iRow]*work[iRow];
811        }
812      }
813      if (reference(sequenceIn)) 
814        devex_+=1.0;
815    }
816  }
817  double oldDevex = weights_[sequenceIn];
818#ifdef CLP_DEBUG
819  if ((model_->messageHandler()->logLevel()&32))
820    printf("old weight %g, new %g\n",oldDevex,devex_);
821#endif
822  double check = max(devex_,oldDevex);;
823  weights_[sequenceIn] = devex_;
824  if ( fabs ( devex_ - oldDevex ) > 1.0e-1 * check ) {
825#ifdef CLP_DEBUG
826    if ((model_->messageHandler()->logLevel()&48)==16)
827      printf("old weight %g, new %g\n",oldDevex,devex_);
828#endif
829    if ( fabs ( devex_ - oldDevex ) > 1.0e5 * check ) {
830      // need to redo
831      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
832                                        *model_->messagesPointer())
833                                          <<oldDevex<<devex_
834                                          <<CoinMessageEol;
835      initializeWeights();
836    }
837  }
838  if (pivotRow>=0) {
839    // set outgoing weight here
840    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
841  }
842}
843// Checks accuracy - just for debug
844void 
845ClpPrimalColumnSteepest::checkAccuracy(int sequence,
846                                       double relativeTolerance,
847                                       CoinIndexedVector * rowArray1,
848                                       CoinIndexedVector * rowArray2)
849{
850  model_->unpack(rowArray1,sequence);
851  model_->factorization()->updateColumn(rowArray2,rowArray1);
852  int number=rowArray1->getNumElements();
853  int * which = rowArray1->getIndices();
854  double * work = rowArray1->denseVector();
855  const int * pivotVariable = model_->pivotVariable();
856 
857  double devex =0.0;
858  int i;
859
860  if (mode_) {
861    for (i=0;i<number;i++) {
862      int iRow = which[i];
863      devex += work[iRow]*work[iRow];
864      work[iRow]=0.0;
865    }
866    devex += ADD_ONE;
867  } else {
868    for (i=0;i<number;i++) {
869      int iRow = which[i];
870      int iPivot = pivotVariable[iRow];
871      if (reference(iPivot)) {
872        devex += work[iRow]*work[iRow];
873      }
874      work[iRow]=0.0;
875    }
876    if (reference(sequence)) 
877      devex+=1.0;
878  }
879
880  double oldDevex = weights_[sequence];
881  double check = max(devex,oldDevex);;
882  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
883    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
884    // update so won't print again
885    weights_[sequence]=devex;
886  }
887  rowArray1->setNumElements(0);
888}
889
890// Initialize weights
891void 
892ClpPrimalColumnSteepest::initializeWeights()
893{
894  int numberRows = model_->numberRows();
895  int numberColumns = model_->numberColumns();
896  int number = numberRows + numberColumns;
897  int iSequence;
898  if (!mode_) {
899    // initialize to 1.0
900    // and set reference framework
901    if (!reference_) {
902      int nWords = (number+31)>>5;
903      reference_ = new unsigned int[nWords];
904      assert (sizeof(unsigned int)==4);
905      // tiny overhead to zero out (stops valgrind complaining)
906      memset(reference_,0,nWords*sizeof(int));
907    }
908   
909    for (iSequence=0;iSequence<number;iSequence++) {
910      weights_[iSequence]=1.0;
911      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
912        setReference(iSequence,false);
913      } else {
914        setReference(iSequence,true);
915      }
916    }
917  } else {
918    CoinIndexedVector * temp = new CoinIndexedVector();
919    temp->reserve(numberRows+
920                  model_->factorization()->maximumPivots());
921    double * array = alternateWeights_->denseVector();
922    int * which = alternateWeights_->getIndices();
923     
924    for (iSequence=0;iSequence<number;iSequence++) {
925      weights_[iSequence]=1.0+ADD_ONE;
926      if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
927          model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
928        model_->unpack(alternateWeights_,iSequence);
929        double value=ADD_ONE;
930        model_->factorization()->updateColumn(temp,alternateWeights_);
931        int number = alternateWeights_->getNumElements();
932        int j;
933        for (j=0;j<number;j++) {
934          int iRow=which[j];
935          value+=array[iRow]*array[iRow];
936          array[iRow]=0.0;
937        }
938        alternateWeights_->setNumElements(0);
939        weights_[iSequence] = value;
940      }
941    }
942    delete temp;
943  }
944}
945// Gets rid of all arrays
946void 
947ClpPrimalColumnSteepest::clearArrays()
948{
949  delete [] weights_;
950  weights_=NULL;
951  delete infeasible_;
952  infeasible_ = NULL;
953  delete alternateWeights_;
954  alternateWeights_ = NULL;
955  delete [] savedWeights_;
956  savedWeights_ = NULL;
957  delete [] reference_;
958  reference_ = NULL;
959  pivotSequence_=-1;
960  state_ = -1;
961  savedPivotSequence_ = -1;
962  savedSequenceOut_ = -1; 
963  devex_ = 0.0;
964}
Note: See TracBrowser for help on using the repository browser.