source: branches/devel-1/ClpDualRowSteepest.cpp @ 2351

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

out CoinCopy? and CoinFill?

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