source: trunk/ClpMatrixBase.cpp @ 437

Last change on this file since 437 was 426, checked in by forrest, 16 years ago

correcting errors

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <iostream>
7
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10#include "ClpMatrixBase.hpp"
11#include "ClpSimplex.hpp"
12
13//#############################################################################
14// Constructors / Destructor / Assignment
15//#############################################################################
16
17//-------------------------------------------------------------------
18// Default Constructor
19//-------------------------------------------------------------------
20ClpMatrixBase::ClpMatrixBase () :
21  rhsOffset_(NULL),
22  startFraction_(0.0),
23  endFraction_(1.0),
24  savedBestDj_(0.0),
25  originalWanted_(0),
26  currentWanted_(0),
27  savedBestSequence_(-1),
28  type_(-1),
29  lastRefresh_(-1),
30  refreshFrequency_(0),
31  minimumObjectsScan_(-1),
32  minimumGoodReducedCosts_(-1),
33  trueSequenceIn_(-1),
34  trueSequenceOut_(-1),
35  skipDualCheck_(false)
36{
37
38}
39
40//-------------------------------------------------------------------
41// Copy constructor
42//-------------------------------------------------------------------
43ClpMatrixBase::ClpMatrixBase (const ClpMatrixBase & rhs) :
44  type_(rhs.type_),
45  skipDualCheck_(rhs.skipDualCheck_)
46{ 
47  startFraction_ = rhs.startFraction_;
48  endFraction_ = rhs.endFraction_;
49  savedBestDj_ = rhs.savedBestDj_;
50  originalWanted_ = rhs.originalWanted_;
51  currentWanted_ = rhs.currentWanted_;
52  savedBestSequence_ = rhs.savedBestSequence_;
53  lastRefresh_ = rhs.lastRefresh_;
54  refreshFrequency_ = rhs.refreshFrequency_;
55  minimumObjectsScan_ = rhs.minimumObjectsScan_;
56  minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
57  trueSequenceIn_ = rhs.trueSequenceIn_;
58  trueSequenceOut_ = rhs.trueSequenceOut_;
59  skipDualCheck_ = rhs.skipDualCheck_;
60  int numberRows = rhs.getNumRows();
61  if (rhs.rhsOffset_&&numberRows) {
62    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
63  } else {
64    rhsOffset_=NULL;
65  }
66}
67
68//-------------------------------------------------------------------
69// Destructor
70//-------------------------------------------------------------------
71ClpMatrixBase::~ClpMatrixBase ()
72{
73  delete [] rhsOffset_;
74}
75
76//----------------------------------------------------------------
77// Assignment operator
78//-------------------------------------------------------------------
79ClpMatrixBase &
80ClpMatrixBase::operator=(const ClpMatrixBase& rhs)
81{
82  if (this != &rhs) {
83    type_ = rhs.type_;
84    delete [] rhsOffset_;
85    int numberRows = rhs.getNumRows();
86    if (rhs.rhsOffset_&&numberRows) {
87      rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
88    } else {
89      rhsOffset_=NULL;
90    }
91    startFraction_ = rhs.startFraction_;
92    endFraction_ = rhs.endFraction_;
93    savedBestDj_ = rhs.savedBestDj_;
94    originalWanted_ = rhs.originalWanted_;
95    currentWanted_ = rhs.currentWanted_;
96    savedBestSequence_ = rhs.savedBestSequence_;
97    lastRefresh_ = rhs.lastRefresh_;
98    refreshFrequency_ = rhs.refreshFrequency_;
99    minimumObjectsScan_ = rhs.minimumObjectsScan_;
100    minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
101    trueSequenceIn_ = rhs.trueSequenceIn_;
102    trueSequenceOut_ = rhs.trueSequenceOut_;
103    skipDualCheck_ = rhs.skipDualCheck_;
104  }
105  return *this;
106}
107// And for scaling - default aborts for when scaling not supported
108void 
109ClpMatrixBase::times(double scalar,
110                     const double * x, double * y,
111                     const double * rowScale, 
112                     const double * columnScale) const
113{
114  if (rowScale) {
115    std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
116    abort();
117  } else {
118    times(scalar,x,y);
119  }
120}
121// And for scaling - default aborts for when scaling not supported
122void 
123ClpMatrixBase::transposeTimes(double scalar,
124                                const double * x, double * y,
125                                const double * rowScale, 
126                                const double * columnScale) const
127{
128  if (rowScale) {
129    std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
130    abort();
131  } else {
132    transposeTimes(scalar,x,y);
133  }
134}
135/* Subset clone (without gaps).  Duplicates are allowed
136   and order is as given.
137   Derived classes need not provide this as it may not always make
138   sense */
139ClpMatrixBase * 
140ClpMatrixBase::subsetClone (
141                            int numberRows, const int * whichRows,
142                            int numberColumns, const int * whichColumns) const
143 
144
145{
146  std::cerr<<"subsetClone not supported - ClpMatrixBase"<<std::endl;
147  abort();
148  return NULL;
149}
150/* Given positive integer weights for each row fills in sum of weights
151   for each column (and slack).
152   Returns weights vector
153   Default returns vector of ones
154*/
155CoinBigIndex * 
156ClpMatrixBase::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
157{
158  int number = model->numberRows()+model->numberColumns();
159  CoinBigIndex * weights = new CoinBigIndex[number];
160  int i;
161  for (i=0;i<number;i++)
162    weights[i]=1;
163  return weights;
164}
165// Append Columns
166void 
167ClpMatrixBase::appendCols(int number, const CoinPackedVectorBase * const * columns)
168{
169  std::cerr<<"appendCols not supported - ClpMatrixBase"<<std::endl;
170  abort();
171}
172// Append Rows
173void 
174ClpMatrixBase::appendRows(int number, const CoinPackedVectorBase * const * rows)
175{
176  std::cerr<<"appendRows not supported - ClpMatrixBase"<<std::endl;
177  abort();
178}
179/* Returns largest and smallest elements of both signs.
180   Largest refers to largest absolute value.
181*/
182void 
183ClpMatrixBase::rangeOfElements(double & smallestNegative, double & largestNegative,
184                       double & smallestPositive, double & largestPositive)
185{
186  smallestNegative=0.0;
187  largestNegative=0.0;
188  smallestPositive=0.0;
189  largestPositive=0.0;
190}
191// Says whether it can do partial pricing
192bool 
193ClpMatrixBase::canDoPartialPricing() const
194{
195  return false; // default is no
196}
197// Partial pricing
198void 
199ClpMatrixBase::partialPricing(ClpSimplex * model, double start, double end,
200                              int & bestSequence, int & numberWanted)
201{
202  std::cerr<<"partialPricing not supported - ClpMatrixBase"<<std::endl;
203  abort();
204}
205/* expands an updated column to allow for extra rows which the main
206   solver does not know about and returns number added.  If the arrays are NULL
207   then returns number of extra entries needed.
208   
209   This will normally be a no-op - it is in for GUB!
210*/
211int 
212ClpMatrixBase::extendUpdated(ClpSimplex * model,CoinIndexedVector * update,int mode)
213{
214  return 0;
215}
216/*
217     utility primal function for dealing with dynamic constraints
218     mode=n see ClpGubMatrix.hpp for definition
219     Remember to update here when settled down
220*/
221void 
222ClpMatrixBase::primalExpanded(ClpSimplex * model,int mode)
223{
224}
225/*
226     utility dual function for dealing with dynamic constraints
227     mode=n see ClpGubMatrix.hpp for definition
228     Remember to update here when settled down
229*/
230void 
231ClpMatrixBase::dualExpanded(ClpSimplex * model,
232                            CoinIndexedVector * array,
233                            double * other,int mode)
234{
235}
236/*
237     general utility function for dealing with dynamic constraints
238     mode=n see ClpGubMatrix.hpp for definition
239     Remember to update here when settled down
240*/
241int
242ClpMatrixBase::generalExpanded(ClpSimplex * model,int mode, int &number)
243{
244  int returnCode=0;
245  switch (mode) {
246    // Fill in pivotVariable but not for key variables
247  case 0:
248    {
249      int i;
250      int numberBasic=number;
251      int numberColumns = model->numberColumns();
252      // Use different array so can build from true pivotVariable_
253      //int * pivotVariable = model->pivotVariable();
254      int * pivotVariable = model->rowArray(0)->getIndices();
255      for (i=0;i<numberColumns;i++) {
256        if (model->getColumnStatus(i) == ClpSimplex::basic) 
257          pivotVariable[numberBasic++]=i;
258      }
259      number = numberBasic;
260    }
261    break;
262    // Do initial extra rows + maximum basic
263  case 2:
264    {
265      number = model->numberRows();
266    }
267    break;
268    // To see if can dual or primal
269  case 4:
270    {
271      returnCode= 3;
272    }
273    break;
274  default:
275    break;
276  }
277  return returnCode;
278}
279// Sets up an effective RHS
280void 
281ClpMatrixBase::useEffectiveRhs(ClpSimplex * model)
282{
283  std::cerr<<"useEffectiveRhs not supported - ClpMatrixBase"<<std::endl;
284  abort();
285}
286/* Returns effective RHS if it is being used.  This is used for long problems
287   or big gub or anywhere where going through full columns is
288   expensive.  This may re-compute */
289double * 
290ClpMatrixBase::rhsOffset(ClpSimplex * model,bool forceRefresh,bool check)
291{
292  if (rhsOffset_) {
293#ifdef CLP_DEBUG
294    if (check) {
295      // no need - but check anyway
296      // zero out basic
297      int numberRows = model->numberRows();
298      int numberColumns = model->numberColumns();
299      double * solution = new double [numberColumns];
300      double * rhs = new double[numberRows];
301      const double * solutionSlack = model->solutionRegion(0);
302      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
303      int iRow;
304      for (iRow=0;iRow<numberRows;iRow++) {
305        if (model->getRowStatus(iRow)!=ClpSimplex::basic)
306          rhs[iRow]=solutionSlack[iRow];
307        else
308          rhs[iRow]=0.0;
309      }
310      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
311        if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
312          solution[iColumn]=0.0;
313      }
314      times(-1.0,solution,rhs);
315      delete [] solution;
316      for (iRow=0;iRow<numberRows;iRow++) {
317        if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
318          printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
319      }
320    }
321#endif
322    if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
323                       lastRefresh_+refreshFrequency_)) {
324      // zero out basic
325      int numberRows = model->numberRows();
326      int numberColumns = model->numberColumns();
327      double * solution = new double [numberColumns];
328      const double * solutionSlack = model->solutionRegion(0);
329      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
330      for (int iRow=0;iRow<numberRows;iRow++) {
331        if (model->getRowStatus(iRow)!=ClpSimplex::basic)
332          rhsOffset_[iRow]=solutionSlack[iRow];
333        else
334          rhsOffset_[iRow]=0.0;
335      }
336      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
337        if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
338          solution[iColumn]=0.0;
339      }
340      times(-1.0,solution,rhsOffset_);
341      delete [] solution;
342      lastRefresh_ = model->numberIterations();
343    }
344  }
345  return rhsOffset_;
346}
347/*
348   update information for a pivot (and effective rhs)
349*/
350int 
351ClpMatrixBase::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
352{
353  if (rhsOffset_) {
354    // update effective rhs
355    int sequenceIn = model->sequenceIn();
356    int sequenceOut = model->sequenceOut();
357    double * solution = model->solutionRegion();
358    int numberColumns = model->numberColumns();
359    if (sequenceIn==sequenceOut) {
360      if (sequenceIn<numberColumns)
361        add(model,rhsOffset_,sequenceIn,oldInValue-solution[sequenceIn]);
362    } else {
363      if (sequenceIn<numberColumns)
364        add(model,rhsOffset_,sequenceIn,oldInValue);
365      if (sequenceOut<numberColumns)
366        add(model,rhsOffset_,sequenceOut,-solution[sequenceOut]);
367    }
368  }
369  return 0;
370}
371int 
372ClpMatrixBase::hiddenRows() const
373{ 
374  return 0;
375}
376/* Creates a variable.  This is called after partial pricing and may modify matrix.
377   May update bestSequence.
378*/
379void 
380ClpMatrixBase::createVariable(ClpSimplex * model, int & bestSequence)
381{
382}
383// Returns reduced cost of a variable
384double 
385ClpMatrixBase::reducedCost(ClpSimplex * model,int sequence) const
386{
387  int numberRows = model->numberRows();
388  int numberColumns = model->numberColumns();
389  if (sequence<numberRows+numberColumns)
390    return model->djRegion()[sequence];
391  else
392    return savedBestDj_;
393}
394/* Just for debug if odd type matrix.
395   Returns number of primal infeasibilities.
396*/
397int 
398ClpMatrixBase::checkFeasible(ClpSimplex * model) const 
399{
400  int numberRows = model->numberRows();
401  double * rhs = new double[numberRows];
402  int numberColumns = model->numberColumns();
403  int iRow;
404  CoinZeroN(rhs,numberRows);
405  times(1.0,model->solutionRegion(),rhs,model->rowScale(),model->columnScale());
406  int iColumn;
407  int logLevel = model->messageHandler()->logLevel();
408  int numberInfeasible=0;
409  const double * rowLower = model->lowerRegion(0);
410  const double * rowUpper = model->upperRegion(0);
411  const double * solution;
412  solution = model->solutionRegion(0);
413  double tolerance = model->primalTolerance()*1.01;
414  for (iRow=0;iRow<numberRows;iRow++) {
415    double value=rhs[iRow];
416    double value2= solution[iRow];
417    if (logLevel>3) {
418      if (fabs(value-value2)>1.0e-8)
419        printf("Row %d stored %g, computed %g\n",iRow,value2,value);
420    }
421    if (value<rowLower[iRow]-tolerance||
422        value>rowUpper[iRow]+tolerance) {
423      numberInfeasible++;
424    }
425    if (value2>rowLower[iRow]+tolerance&&
426        value2<rowUpper[iRow]-tolerance&&
427        model->getRowStatus(iRow)!=ClpSimplex::basic) {
428      assert (model->getRowStatus(iRow)==ClpSimplex::superBasic);
429    }
430  }
431  const double * columnLower = model->lowerRegion(1);
432  const double * columnUpper = model->upperRegion(1);
433  solution = model->solutionRegion(1);
434  for (iColumn=0;iColumn<numberColumns;iColumn++) {
435    double value= solution[iColumn];
436    if (value<columnLower[iColumn]-tolerance||
437        value>columnUpper[iColumn]+tolerance) {
438      numberInfeasible++;
439    }
440    if (value>columnLower[iColumn]+tolerance&&
441        value<columnUpper[iColumn]-tolerance&&
442        model->getColumnStatus(iColumn)!=ClpSimplex::basic) {
443      assert (model->getColumnStatus(iColumn)==ClpSimplex::superBasic);
444    }
445  }
446  delete [] rhs;
447  return numberInfeasible;
448}
449// Correct sequence in and out to give true value
450void 
451ClpMatrixBase::correctSequence(int & sequenceIn, int & sequenceOut) const
452{
453}
454// Really scale matrix
455void 
456ClpMatrixBase::reallyScale(const double * rowScale, const double * columnScale)
457{
458  std::cerr<<"reallyScale not supported - ClpMatrixBase"<<std::endl;
459  abort();
460}
461
462
Note: See TracBrowser for help on using the repository browser.