source: trunk/ClpDualRowSteepest.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: 17.1 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include "ClpSimplex.hpp"
6#include "ClpDualRowSteepest.hpp"
7#include "CoinIndexedVector.hpp"
8#include "ClpFactorization.hpp"
9#include "CoinHelperFunctions.hpp"
10#include <cstdio>
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpDualRowSteepest::ClpDualRowSteepest (int mode) 
20  : ClpDualRowPivot(),
21    state_(-1),
22    mode_(mode),
23    weights_(NULL),
24    infeasible_(NULL),
25    alternateWeights_(NULL),
26    savedWeights_(NULL)
27{
28  type_=2+64*mode;
29}
30
31//-------------------------------------------------------------------
32// Copy constructor
33//-------------------------------------------------------------------
34ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs) 
35: ClpDualRowPivot(rhs)
36{ 
37  state_=rhs.state_;
38  mode_ = rhs.mode_;
39  model_ = rhs.model_;
40  if (rhs.infeasible_) {
41    infeasible_= new CoinIndexedVector(rhs.infeasible_);
42  } else {
43    infeasible_=NULL;
44  }
45  if (rhs.weights_) {
46    assert(model_);
47    int number = model_->numberRows();
48    weights_= new double[number];
49    ClpDisjointCopyN(rhs.weights_,number,weights_);
50  } else {
51    weights_=NULL;
52  }
53  if (rhs.alternateWeights_) {
54    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
55  } else {
56    alternateWeights_=NULL;
57  }
58  if (rhs.savedWeights_) {
59    savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
60  } else {
61    savedWeights_=NULL;
62  }
63}
64
65//-------------------------------------------------------------------
66// Destructor
67//-------------------------------------------------------------------
68ClpDualRowSteepest::~ClpDualRowSteepest ()
69{
70  delete [] weights_;
71  delete infeasible_;
72  delete alternateWeights_;
73  delete savedWeights_;
74}
75
76//----------------------------------------------------------------
77// Assignment operator
78//-------------------------------------------------------------------
79ClpDualRowSteepest &
80ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
81{
82  if (this != &rhs) {
83    ClpDualRowPivot::operator=(rhs);
84    state_=rhs.state_;
85    mode_ = rhs.mode_;
86    model_ = rhs.model_;
87    delete [] weights_;
88    delete infeasible_;
89    delete alternateWeights_;
90    delete savedWeights_;
91    if (rhs.infeasible_!=NULL) {
92      infeasible_= new CoinIndexedVector(rhs.infeasible_);
93    } else {
94      infeasible_=NULL;
95    }
96    if (rhs.weights_!=NULL) {
97      assert(model_);
98      int number = model_->numberRows();
99      weights_= new double[number];
100      ClpDisjointCopyN(rhs.weights_,number,weights_);
101    } else {
102      weights_=NULL;
103    }
104    if (rhs.alternateWeights_!=NULL) {
105      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
106    } else {
107      alternateWeights_=NULL;
108    }
109    if (rhs.savedWeights_!=NULL) {
110      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
111    } else {
112      savedWeights_=NULL;
113    }
114  }
115  return *this;
116}
117
118// Returns pivot row, -1 if none
119int 
120ClpDualRowSteepest::pivotRow()
121{
122  assert(model_);
123  int i,iRow;
124  double * infeas = infeasible_->denseVector();
125  double largest=1.0e-50;
126  int * index = infeasible_->getIndices();
127  int number = infeasible_->getNumElements();
128  const int * pivotVariable =model_->pivotVariable();
129  int chosenRow=-1;
130  int lastPivotRow = model_->pivotRow();
131  double tolerance=model_->currentPrimalTolerance();
132  // we can't really trust infeasibilities if there is primal error
133  // this coding has to mimic coding in checkPrimalSolution
134  double error = min(1.0e-3,model_->largestPrimalError());
135  // allow tolerance at least slightly bigger than standard
136  tolerance = tolerance +  error;
137  tolerance *= tolerance; // as we are using squares
138  double * solution = model_->solutionRegion();
139  double * lower = model_->lowerRegion();
140  double * upper = model_->upperRegion();
141  // do last pivot row one here
142  //#define COLUMN_BIAS 4.0
143  //#define FIXED_BIAS 10.0
144  if (lastPivotRow>=0) {
145#ifdef COLUMN_BIAS
146    int numberColumns = model_->numberColumns();
147#endif
148    int iPivot=pivotVariable[lastPivotRow];
149    double value = solution[iPivot];
150    double lower = model_->lower(iPivot);
151    double upper = model_->upper(iPivot);
152    if (value>upper+tolerance) {
153      value -= upper;
154      value *= value;
155#ifdef COLUMN_BIAS
156      if (iPivot<numberColumns)
157        value *= COLUMN_BIAS; // bias towards columns
158#endif
159      // store square in list
160      if (infeas[lastPivotRow])
161        infeas[lastPivotRow] = value; // already there
162      else
163        infeasible_->quickAdd(lastPivotRow,value);
164    } else if (value<lower-tolerance) {
165      value -= lower;
166      value *= value;
167#ifdef COLUMN_BIAS
168      if (iPivot<numberColumns)
169        value *= COLUMN_BIAS; // bias towards columns
170#endif
171      // store square in list
172      if (infeas[lastPivotRow])
173        infeas[lastPivotRow] = value; // already there
174      else
175        infeasible_->add(lastPivotRow,value);
176    } else {
177      // feasible - was it infeasible - if so set tiny
178      if (infeas[lastPivotRow])
179        infeas[lastPivotRow] = 1.0e-100;
180    }
181    number = infeasible_->getNumElements();
182  }
183  for (i=0;i<number;i++) {
184    iRow = index[i];
185    double value = infeas[iRow];
186    if (value>largest*weights_[iRow]&&value>tolerance) {
187      // make last pivot row last resort choice
188      if (iRow==lastPivotRow) {
189        if (value*1.0e-10<largest*weights_[iRow]) 
190          continue;
191        else 
192          value *= 1.0e-10;
193      }
194      int iSequence = pivotVariable[iRow];
195      if (!model_->flagged(iSequence)) {
196        //#define CLP_DEBUG 1
197#ifdef CLP_DEBUG
198        double value2=0.0;
199        if (solution[iSequence]>upper[iSequence]+tolerance) 
200          value2=solution[iSequence]-upper[iSequence];
201        else if (solution[iSequence]<lower[iSequence]-tolerance) 
202          value2=solution[iSequence]-lower[iSequence];
203        assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
204#endif
205        if (solution[iSequence]>upper[iSequence]+tolerance||
206            solution[iSequence]<lower[iSequence]-tolerance) {
207          chosenRow=iRow;
208          largest=value/weights_[iRow];
209        }
210      }
211    }
212  }
213  return chosenRow;
214}
215#define TRY_NORM 1.0e-4
216// Updates weights
217void 
218ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
219                                  CoinIndexedVector * spare,
220                                  CoinIndexedVector * updatedColumn)
221{
222  // clear other region
223  alternateWeights_->clear();
224  double norm = 0.0;
225  double * work = input->denseVector();
226  int number = input->getNumElements();
227  int * which = input->getIndices();
228  double * work2 = spare->denseVector();
229  int * which2 = spare->getIndices();
230  double * work3 = alternateWeights_->denseVector();
231  int * which3 = alternateWeights_->getIndices();
232  int i;
233#if CLP_DEBUG>2
234  // Very expensive debug
235  {
236    int numberRows = model_->numberRows();
237    CoinIndexedVector * temp = new CoinIndexedVector();
238    temp->reserve(numberRows+
239                  model_->factorization()->maximumPivots());
240    double * array = alternateWeights_->denseVector();
241    int * which = alternateWeights_->getIndices();
242    for (i=0;i<numberRows;i++) {
243      double value=0.0;
244      array[i] = 1.0;
245      which[0] = i;
246      alternateWeights_->setNumElements(1);
247      model_->factorization()->updateColumnTranspose(temp,
248                                                     alternateWeights_);
249      int number = alternateWeights_->getNumElements();
250      int j;
251      for (j=0;j<number;j++) {
252        int iRow=which[j];
253        value+=array[iRow]*array[iRow];
254        array[iRow]=0.0;
255      }
256      alternateWeights_->setNumElements(0);
257      if (fabs(weights_[i]-value)>1.0e-4)
258        printf("%d old %g, true %g\n",i,weights_[i],value);
259    }
260    delete temp;
261  }
262#endif
263  for (i=0;i<number;i++) {
264    int iRow = which[i];
265    double value = work[iRow];
266    norm += value*value;
267    work2[iRow]=value;
268    which2[i]=iRow;
269  }
270  spare->setNumElements(number);
271  // ftran
272  model_->factorization()->updateColumn(alternateWeights_,spare);
273  // alternateWeights_ should still be empty
274  int pivotRow = model_->pivotRow();
275#ifdef CLP_DEBUG
276  if ( model_->logLevel (  ) >4  && 
277       fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
278    printf("on row %d, true weight %g, old %g\n",
279           pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
280#endif
281  // could re-initialize here (could be expensive)
282  norm /= model_->alpha() * model_->alpha();
283
284  assert(norm);
285  if (norm < TRY_NORM) 
286    norm = TRY_NORM;
287  if (norm != 0.) {
288    double multiplier = 2.0 / model_->alpha();
289    // look at updated column
290    work = updatedColumn->denseVector();
291    number = updatedColumn->getNumElements();
292    which = updatedColumn->getIndices();
293
294    int nSave=0;
295
296    for (i =0; i < number; i++) {
297      int iRow = which[i];
298      double theta = work[iRow];
299      if (theta) {
300        double devex = weights_[iRow];
301        work3[iRow]=devex; // save old
302        which3[nSave++]=iRow;
303        double value = work2[iRow];
304        devex +=  theta * (theta*norm+value * multiplier);
305        if (devex < TRY_NORM) 
306          devex = TRY_NORM;
307        weights_[iRow]=devex;
308      }
309    }
310#ifdef CLP_DEBUG
311    assert(work3[pivotRow]&&work[pivotRow]);
312#endif
313    alternateWeights_->setNumElements(nSave);
314    if (norm < TRY_NORM) 
315      norm = TRY_NORM;
316    weights_[pivotRow] = norm;
317  }
318  spare->clear();
319#ifdef CLP_DEBUG
320  spare->checkClear();
321#endif
322}
323 
324/* Updates primal solution (and maybe list of candidates)
325   Uses input vector which it deletes
326   Computes change in objective function
327*/
328void 
329ClpDualRowSteepest::updatePrimalSolution(
330                                        CoinIndexedVector * primalUpdate,
331                                        double primalRatio,
332                                        double & objectiveChange)
333{
334  double * work = primalUpdate->denseVector();
335  int number = primalUpdate->getNumElements();
336  int * which = primalUpdate->getIndices();
337  int i;
338  double changeObj=0.0;
339  double tolerance=model_->currentPrimalTolerance();
340  const int * pivotVariable = model_->pivotVariable();
341  double * infeas = infeasible_->denseVector();
342  int pivotRow = model_->pivotRow();
343  double * solution = model_->solutionRegion();
344#ifdef COLUMN_BIAS
345  int numberColumns = model_->numberColumns();
346#endif
347  for (i=0;i<number;i++) {
348    int iRow=which[i];
349    int iPivot=pivotVariable[iRow];
350    double value = solution[iPivot];
351    double cost = model_->cost(iPivot);
352    double change = primalRatio*work[iRow];
353    value -= change;
354    changeObj -= change*cost;
355    solution[iPivot] = value;
356    double lower = model_->lower(iPivot);
357    double upper = model_->upper(iPivot);
358    // But if pivot row then use value of incoming
359    // Although it is safer to recompute before next selection
360    // in case something odd happens
361    if (iRow==pivotRow) {
362      iPivot = model_->sequenceIn();
363      lower = model_->lower(iPivot);
364      upper = model_->upper(iPivot);
365      value = model_->valueIncomingDual();
366    }
367    if (value<lower-tolerance) {
368      value -= lower;
369      value *= value;
370#ifdef COLUMN_BIAS
371      if (iPivot<numberColumns)
372        value *= COLUMN_BIAS; // bias towards columns
373#endif
374#ifdef FIXED_BIAS
375      if (lower==upper)
376        value *= FIXED_BIAS; // bias towards taking out fixed variables
377#endif
378      // store square in list
379      if (infeas[iRow])
380        infeas[iRow] = value; // already there
381      else
382        infeasible_->quickAdd(iRow,value);
383    } else if (value>upper+tolerance) {
384      value -= upper;
385      value *= value;
386#ifdef COLUMN_BIAS
387      if (iPivot<numberColumns)
388        value *= COLUMN_BIAS; // bias towards columns
389#endif
390#ifdef FIXED_BIAS
391      if (lower==upper)
392        value *= FIXED_BIAS; // bias towards taking out fixed variables
393#endif
394      // store square in list
395      if (infeas[iRow])
396        infeas[iRow] = value; // already there
397      else
398        infeasible_->quickAdd(iRow,value);
399    } else {
400      // feasible - was it infeasible - if so set tiny
401      if (infeas[iRow])
402        infeas[iRow] = 1.0e-100;
403    }
404    work[iRow]=0.0;
405  }
406  primalUpdate->setNumElements(0);
407  objectiveChange += changeObj;
408}
409/* Saves any weights round factorization as pivot rows may change
410   1) before factorization
411   2) after factorization
412   3) just redo infeasibilities
413   4) restore weights
414*/
415void 
416ClpDualRowSteepest::saveWeights(ClpSimplex * model,int mode)
417{
418  // alternateWeights_ is defined as indexed but is treated oddly
419  model_ = model;
420  int numberRows = model_->numberRows();
421  int numberColumns = model_->numberColumns();
422  const int * pivotVariable = model_->pivotVariable();
423  int i;
424  if (mode==1&&weights_) {
425    alternateWeights_->clear();
426    // change from row numbers to sequence numbers
427    int * which = alternateWeights_->getIndices();
428    for (i=0;i<numberRows;i++) {
429      int iPivot=pivotVariable[i];
430      which[i]=iPivot;
431    }
432    state_=1;
433  } else if (mode==2||mode==4||mode==5) {
434    // restore
435    if (!weights_||state_==-1||mode==5) {
436      // initialize weights
437      delete [] weights_;
438      delete alternateWeights_;
439      weights_ = new double[numberRows];
440      alternateWeights_ = new CoinIndexedVector();
441      // enough space so can use it for factorization
442      alternateWeights_->reserve(numberRows+
443                                 model_->factorization()->maximumPivots());
444      if (!mode_||mode==5) {
445        // initialize to 1.0 (can we do better?)
446        for (i=0;i<numberRows;i++) {
447          weights_[i]=1.0;
448        }
449      } else {
450        CoinIndexedVector * temp = new CoinIndexedVector();
451        temp->reserve(numberRows+
452                      model_->factorization()->maximumPivots());
453        double * array = alternateWeights_->denseVector();
454        int * which = alternateWeights_->getIndices();
455        for (i=0;i<numberRows;i++) {
456          double value=0.0;
457          array[i] = 1.0;
458          which[0] = i;
459          alternateWeights_->setNumElements(1);
460          model_->factorization()->updateColumnTranspose(temp,
461                                                         alternateWeights_);
462          int number = alternateWeights_->getNumElements();
463          int j;
464          for (j=0;j<number;j++) {
465            int iRow=which[j];
466            value+=array[iRow]*array[iRow];
467            array[iRow]=0.0;
468          }
469          alternateWeights_->setNumElements(0);
470          weights_[i] = value;
471        }
472        delete temp;
473      }
474      // create saved weights (not really indexedvector)
475      savedWeights_ = new CoinIndexedVector();
476      savedWeights_->reserve(numberRows);
477     
478      double * array = savedWeights_->denseVector();
479      int * which = savedWeights_->getIndices();
480      for (i=0;i<numberRows;i++) {
481        array[i]=weights_[i];
482        which[i]=pivotVariable[i];
483      }
484    } else {
485      int * which = alternateWeights_->getIndices();
486      if (mode!=4) {
487        // save
488        memcpy(savedWeights_->getIndices(),which,
489               numberRows*sizeof(int));
490        memcpy(savedWeights_->denseVector(),weights_,
491               numberRows*sizeof(double));
492      } else {
493        // restore
494        memcpy(which,savedWeights_->getIndices(),
495               numberRows*sizeof(int));
496        memcpy(weights_,savedWeights_->denseVector(),
497               numberRows*sizeof(double));
498      }
499      // restore (a bit slow - but only every re-factorization)
500      double * array = new double[numberRows+numberColumns];
501      for (i=0;i<numberRows;i++) {
502        int iSeq=which[i];
503        array[iSeq]=weights_[i];
504      }
505      for (i=0;i<numberRows;i++) {
506        int iPivot=pivotVariable[i];
507        weights_[i]=array[iPivot];
508        if (weights_[i]<TRY_NORM)
509          weights_[i] = TRY_NORM; // may need to check more
510      }
511      delete [] array;
512    }
513    state_=0;
514    // set up infeasibilities
515    if (!infeasible_) {
516      infeasible_ = new CoinIndexedVector();
517      infeasible_->reserve(numberRows);
518    }
519  }
520  if (mode>=2) {
521    infeasible_->clear();
522    int iRow;
523    const int * pivotVariable = model_->pivotVariable();
524    double tolerance=model_->currentPrimalTolerance();
525    for (iRow=0;iRow<numberRows;iRow++) {
526      int iPivot=pivotVariable[iRow];
527      double value = model_->solution(iPivot);
528      double lower = model_->lower(iPivot);
529      double upper = model_->upper(iPivot);
530      if (value<lower-tolerance) {
531        value -= lower;
532        value *= value;
533#ifdef COLUMN_BIAS
534        if (iPivot<numberColumns)
535          value *= COLUMN_BIAS; // bias towards columns
536#endif
537#ifdef FIXED_BIAS
538        if (lower==upper)
539          value *= FIXED_BIAS; // bias towards taking out fixed variables
540#endif
541        // store square in list
542        infeasible_->quickAdd(iRow,value);
543      } else if (value>upper+tolerance) {
544        value -= upper;
545        value *= value;
546#ifdef COLUMN_BIAS
547        if (iPivot<numberColumns)
548          value *= COLUMN_BIAS; // bias towards columns
549#endif
550#ifdef FIXED_BIAS
551        if (lower==upper)
552          value *= FIXED_BIAS; // bias towards taking out fixed variables
553#endif
554        // store square in list
555        infeasible_->quickAdd(iRow,value);
556      }
557    }
558  }
559}
560// Gets rid of last update
561void 
562ClpDualRowSteepest::unrollWeights()
563{
564  double * saved = alternateWeights_->denseVector();
565  int number = alternateWeights_->getNumElements();
566  int * which = alternateWeights_->getIndices();
567  int i;
568  for (i=0;i<number;i++) {
569    int iRow = which[i];
570    weights_[iRow]=saved[iRow];
571    saved[iRow]=0.0;
572  }
573  alternateWeights_->setNumElements(0);
574}
575//-------------------------------------------------------------------
576// Clone
577//-------------------------------------------------------------------
578ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
579{
580  if (CopyData) {
581    return new ClpDualRowSteepest(*this);
582  } else {
583    return new ClpDualRowSteepest();
584  }
585}
586// Gets rid of all arrays
587void 
588ClpDualRowSteepest::clearArrays()
589{
590  delete [] weights_;
591  weights_=NULL;
592  delete infeasible_;
593  infeasible_ = NULL;
594  delete alternateWeights_;
595  alternateWeights_ = NULL;
596  delete savedWeights_;
597  savedWeights_ = NULL;
598  state_ =-1;
599}
600
Note: See TracBrowser for help on using the repository browser.