source: stable/1.6/Clp/src/ClpMatrixBase.cpp @ 1609

Last change on this file since 1609 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.7 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,
127                              double * spare) const
128{
129  if (rowScale) {
130    std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
131    abort();
132  } else {
133    transposeTimes(scalar,x,y);
134  }
135}
136/* Subset clone (without gaps).  Duplicates are allowed
137   and order is as given.
138   Derived classes need not provide this as it may not always make
139   sense */
140ClpMatrixBase * 
141ClpMatrixBase::subsetClone (
142                            int numberRows, const int * whichRows,
143                            int numberColumns, const int * whichColumns) const
144 
145
146{
147  std::cerr<<"subsetClone not supported - ClpMatrixBase"<<std::endl;
148  abort();
149  return NULL;
150}
151/* Given positive integer weights for each row fills in sum of weights
152   for each column (and slack).
153   Returns weights vector
154   Default returns vector of ones
155*/
156CoinBigIndex * 
157ClpMatrixBase::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
158{
159  int number = model->numberRows()+model->numberColumns();
160  CoinBigIndex * weights = new CoinBigIndex[number];
161  int i;
162  for (i=0;i<number;i++)
163    weights[i]=1;
164  return weights;
165}
166#ifndef CLP_NO_VECTOR
167// Append Columns
168void 
169ClpMatrixBase::appendCols(int number, const CoinPackedVectorBase * const * columns)
170{
171  std::cerr<<"appendCols not supported - ClpMatrixBase"<<std::endl;
172  abort();
173}
174// Append Rows
175void 
176ClpMatrixBase::appendRows(int number, const CoinPackedVectorBase * const * rows)
177{
178  std::cerr<<"appendRows not supported - ClpMatrixBase"<<std::endl;
179  abort();
180}
181#endif
182/* Returns largest and smallest elements of both signs.
183   Largest refers to largest absolute value.
184*/
185void 
186ClpMatrixBase::rangeOfElements(double & smallestNegative, double & largestNegative,
187                       double & smallestPositive, double & largestPositive)
188{
189  smallestNegative=0.0;
190  largestNegative=0.0;
191  smallestPositive=0.0;
192  largestPositive=0.0;
193}
194/* The length of a major-dimension vector. */
195int
196ClpMatrixBase::getVectorLength(int index) const 
197{
198  return getVectorLengths()[index];
199}
200// Says whether it can do partial pricing
201bool 
202ClpMatrixBase::canDoPartialPricing() const
203{
204  return false; // default is no
205}
206/* Return <code>x *A</code> in <code>z</code> but
207   just for number indices in y.
208   Default cheats with fake CoinIndexedVector and
209   then calls subsetTransposeTimes */
210void 
211ClpMatrixBase::listTransposeTimes(const ClpSimplex * model,
212                                  double * x,
213                                  int * y,
214                                  int number,
215                                  double * z) const
216{
217  CoinIndexedVector pi;
218  CoinIndexedVector list;
219  CoinIndexedVector output;
220  int * saveIndices = list.getIndices();
221  list.setNumElements(number);
222  list.setIndexVector(y);
223  double * savePi = pi.denseVector();
224  pi.setDenseVector(x);
225  double * saveOutput = output.denseVector();
226  output.setDenseVector(z);
227  output.setPacked();
228  subsetTransposeTimes(model,&pi,&list,&output);
229  // restore settings
230  list.setIndexVector(saveIndices);
231  pi.setDenseVector(savePi);
232  output.setDenseVector(saveOutput);
233}
234// Partial pricing
235void 
236ClpMatrixBase::partialPricing(ClpSimplex * model, double start, double end,
237                              int & bestSequence, int & numberWanted)
238{
239  std::cerr<<"partialPricing not supported - ClpMatrixBase"<<std::endl;
240  abort();
241}
242/* expands an updated column to allow for extra rows which the main
243   solver does not know about and returns number added.  If the arrays are NULL
244   then returns number of extra entries needed.
245   
246   This will normally be a no-op - it is in for GUB!
247*/
248int 
249ClpMatrixBase::extendUpdated(ClpSimplex * model,CoinIndexedVector * update,int mode)
250{
251  return 0;
252}
253/*
254     utility primal function for dealing with dynamic constraints
255     mode=n see ClpGubMatrix.hpp for definition
256     Remember to update here when settled down
257*/
258void 
259ClpMatrixBase::primalExpanded(ClpSimplex * model,int mode)
260{
261}
262/*
263     utility dual function for dealing with dynamic constraints
264     mode=n see ClpGubMatrix.hpp for definition
265     Remember to update here when settled down
266*/
267void 
268ClpMatrixBase::dualExpanded(ClpSimplex * model,
269                            CoinIndexedVector * array,
270                            double * other,int mode)
271{
272}
273/*
274     general utility function for dealing with dynamic constraints
275     mode=n see ClpGubMatrix.hpp for definition
276     Remember to update here when settled down
277*/
278int
279ClpMatrixBase::generalExpanded(ClpSimplex * model,int mode, int &number)
280{
281  int returnCode=0;
282  switch (mode) {
283    // Fill in pivotVariable but not for key variables
284  case 0:
285    {
286      int i;
287      int numberBasic=number;
288      int numberColumns = model->numberColumns();
289      // Use different array so can build from true pivotVariable_
290      //int * pivotVariable = model->pivotVariable();
291      int * pivotVariable = model->rowArray(0)->getIndices();
292      for (i=0;i<numberColumns;i++) {
293        if (model->getColumnStatus(i) == ClpSimplex::basic) 
294          pivotVariable[numberBasic++]=i;
295      }
296      number = numberBasic;
297    }
298    break;
299    // Do initial extra rows + maximum basic
300  case 2:
301    {
302      number = model->numberRows();
303    }
304    break;
305    // To see if can dual or primal
306  case 4:
307    {
308      returnCode= 3;
309    }
310    break;
311  default:
312    break;
313  }
314  return returnCode;
315}
316// Sets up an effective RHS
317void 
318ClpMatrixBase::useEffectiveRhs(ClpSimplex * model)
319{
320  std::cerr<<"useEffectiveRhs not supported - ClpMatrixBase"<<std::endl;
321  abort();
322}
323/* Returns effective RHS if it is being used.  This is used for long problems
324   or big gub or anywhere where going through full columns is
325   expensive.  This may re-compute */
326double * 
327ClpMatrixBase::rhsOffset(ClpSimplex * model,bool forceRefresh,bool check)
328{
329  if (rhsOffset_) {
330#ifdef CLP_DEBUG
331    if (check) {
332      // no need - but check anyway
333      // zero out basic
334      int numberRows = model->numberRows();
335      int numberColumns = model->numberColumns();
336      double * solution = new double [numberColumns];
337      double * rhs = new double[numberRows];
338      const double * solutionSlack = model->solutionRegion(0);
339      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
340      int iRow;
341      for (iRow=0;iRow<numberRows;iRow++) {
342        if (model->getRowStatus(iRow)!=ClpSimplex::basic)
343          rhs[iRow]=solutionSlack[iRow];
344        else
345          rhs[iRow]=0.0;
346      }
347      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
348        if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
349          solution[iColumn]=0.0;
350      }
351      times(-1.0,solution,rhs);
352      delete [] solution;
353      for (iRow=0;iRow<numberRows;iRow++) {
354        if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
355          printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
356      }
357    }
358#endif
359    if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
360                       lastRefresh_+refreshFrequency_)) {
361      // zero out basic
362      int numberRows = model->numberRows();
363      int numberColumns = model->numberColumns();
364      double * solution = new double [numberColumns];
365      const double * solutionSlack = model->solutionRegion(0);
366      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
367      for (int iRow=0;iRow<numberRows;iRow++) {
368        if (model->getRowStatus(iRow)!=ClpSimplex::basic)
369          rhsOffset_[iRow]=solutionSlack[iRow];
370        else
371          rhsOffset_[iRow]=0.0;
372      }
373      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
374        if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
375          solution[iColumn]=0.0;
376      }
377      times(-1.0,solution,rhsOffset_);
378      delete [] solution;
379      lastRefresh_ = model->numberIterations();
380    }
381  }
382  return rhsOffset_;
383}
384/*
385   update information for a pivot (and effective rhs)
386*/
387int 
388ClpMatrixBase::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
389{
390  if (rhsOffset_) {
391    // update effective rhs
392    int sequenceIn = model->sequenceIn();
393    int sequenceOut = model->sequenceOut();
394    double * solution = model->solutionRegion();
395    int numberColumns = model->numberColumns();
396    if (sequenceIn==sequenceOut) {
397      if (sequenceIn<numberColumns)
398        add(model,rhsOffset_,sequenceIn,oldInValue-solution[sequenceIn]);
399    } else {
400      if (sequenceIn<numberColumns)
401        add(model,rhsOffset_,sequenceIn,oldInValue);
402      if (sequenceOut<numberColumns)
403        add(model,rhsOffset_,sequenceOut,-solution[sequenceOut]);
404    }
405  }
406  return 0;
407}
408int 
409ClpMatrixBase::hiddenRows() const
410{ 
411  return 0;
412}
413/* Creates a variable.  This is called after partial pricing and may modify matrix.
414   May update bestSequence.
415*/
416void 
417ClpMatrixBase::createVariable(ClpSimplex * model, int & bestSequence)
418{
419}
420// Returns reduced cost of a variable
421double 
422ClpMatrixBase::reducedCost(ClpSimplex * model,int sequence) const
423{
424  int numberRows = model->numberRows();
425  int numberColumns = model->numberColumns();
426  if (sequence<numberRows+numberColumns)
427    return model->djRegion()[sequence];
428  else
429    return savedBestDj_;
430}
431/* Just for debug if odd type matrix.
432   Returns number and sum of primal infeasibilities.
433*/
434int 
435ClpMatrixBase::checkFeasible(ClpSimplex * model, double & sum) const 
436{
437  int numberRows = model->numberRows();
438  double * rhs = new double[numberRows];
439  int numberColumns = model->numberColumns();
440  int iRow;
441  CoinZeroN(rhs,numberRows);
442  times(1.0,model->solutionRegion(),rhs,model->rowScale(),model->columnScale());
443  int iColumn;
444  int logLevel = model->messageHandler()->logLevel();
445  int numberInfeasible=0;
446  const double * rowLower = model->lowerRegion(0);
447  const double * rowUpper = model->upperRegion(0);
448  const double * solution;
449  solution = model->solutionRegion(0);
450  double tolerance = model->primalTolerance()*1.01;
451  sum=0.0;
452  for (iRow=0;iRow<numberRows;iRow++) {
453    double value=rhs[iRow];
454    double value2= solution[iRow];
455    if (logLevel>3) {
456      if (fabs(value-value2)>1.0e-8)
457        printf("Row %d stored %g, computed %g\n",iRow,value2,value);
458    }
459    if (value<rowLower[iRow]-tolerance||
460        value>rowUpper[iRow]+tolerance) {
461      numberInfeasible++;
462      sum += CoinMax(rowLower[iRow]-value,value-rowUpper[iRow]);
463    }
464    if (value2>rowLower[iRow]+tolerance&&
465        value2<rowUpper[iRow]-tolerance&&
466        model->getRowStatus(iRow)!=ClpSimplex::basic) {
467      assert (model->getRowStatus(iRow)==ClpSimplex::superBasic);
468    }
469  }
470  const double * columnLower = model->lowerRegion(1);
471  const double * columnUpper = model->upperRegion(1);
472  solution = model->solutionRegion(1);
473  for (iColumn=0;iColumn<numberColumns;iColumn++) {
474    double value= solution[iColumn];
475    if (value<columnLower[iColumn]-tolerance||
476        value>columnUpper[iColumn]+tolerance) {
477      numberInfeasible++;
478      sum += CoinMax(columnLower[iColumn]-value,value-columnUpper[iColumn]);
479    }
480    if (value>columnLower[iColumn]+tolerance&&
481        value<columnUpper[iColumn]-tolerance&&
482        model->getColumnStatus(iColumn)!=ClpSimplex::basic) {
483      assert (model->getColumnStatus(iColumn)==ClpSimplex::superBasic);
484    }
485  }
486  delete [] rhs;
487  return numberInfeasible;
488}
489// These have to match ClpPrimalColumnSteepest version
490#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
491// Updates second array for steepest and does devex weights (need not be coded)
492void 
493ClpMatrixBase::subsetTimes2(const ClpSimplex * model,
494                            CoinIndexedVector * dj1,
495                            const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
496                            double referenceIn, double devex,
497                            // Array for exact devex to say what is in reference framework
498                            unsigned int * reference,
499                            double * weights, double scaleFactor)
500{
501  // get subset which have nonzero tableau elements
502  subsetTransposeTimes(model,pi2,dj1,dj2);
503  bool killDjs = (scaleFactor==0.0);
504  if (!scaleFactor)
505    scaleFactor=1.0;
506  // columns
507 
508  int number = dj1->getNumElements();
509  const int * index = dj1->getIndices();
510  double * updateBy = dj1->denseVector();
511  double * updateBy2 = dj2->denseVector();
512 
513  for (int j=0;j<number;j++) {
514    double thisWeight;
515    double pivot;
516    double pivotSquared;
517    int iSequence = index[j];
518    double value2 = updateBy[j];
519    if (killDjs)
520      updateBy[j]=0.0;
521    double modification=updateBy2[j];
522    updateBy2[j]=0.0;
523    ClpSimplex::Status status = model->getStatus(iSequence);
524   
525    if (status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) {
526      thisWeight = weights[iSequence];
527      pivot = value2*scaleFactor;
528      pivotSquared = pivot * pivot;
529     
530      thisWeight += pivotSquared * devex + pivot * modification;
531      if (thisWeight<DEVEX_TRY_NORM) {
532        if (referenceIn<0.0) {
533          // steepest
534          thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
535        } else {
536          // exact
537          thisWeight = referenceIn*pivotSquared;
538          if (reference(iSequence))
539            thisWeight += 1.0;
540          thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
541        }
542      }
543      weights[iSequence] = thisWeight;
544    }
545  }
546  dj2->setNumElements(0);
547}
548// Correct sequence in and out to give true value
549void 
550ClpMatrixBase::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) 
551{
552}
553// Really scale matrix
554void 
555ClpMatrixBase::reallyScale(const double * rowScale, const double * columnScale)
556{
557  std::cerr<<"reallyScale not supported - ClpMatrixBase"<<std::endl;
558  abort();
559}
560// Updates two arrays for steepest
561void 
562ClpMatrixBase::transposeTimes2(const ClpSimplex * model,
563                               const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
564                               const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
565                               CoinIndexedVector * spare,
566                               double referenceIn, double devex,
567                               // Array for exact devex to say what is in reference framework
568                               unsigned int * reference,
569                               double * weights, double scaleFactor)
570{
571  std::cerr<<"transposeTimes2 not supported - ClpMatrixBase"<<std::endl;
572  abort();
573}
574/* Set the dimensions of the matrix. In effect, append new empty
575   columns/rows to the matrix. A negative number for either dimension
576   means that that dimension doesn't change. Otherwise the new dimensions
577   MUST be at least as large as the current ones otherwise an exception
578   is thrown. */
579void 
580ClpMatrixBase::setDimensions(int numrows, int numcols){
581  // If odd matrix assume user knows what they are doing
582}
583/* Append a set of rows/columns to the end of the matrix. Returns number of errors
584   i.e. if any of the new rows/columns contain an index that's larger than the
585   number of columns-1/rows-1 (if numberOther>0) or duplicates
586   If 0 then rows, 1 if columns */
587int 
588ClpMatrixBase::appendMatrix(int number, int type,
589                            const CoinBigIndex * starts, const int * index,
590                            const double * element, int numberOther)
591{
592  std::cerr<<"appendMatrix not supported - ClpMatrixBase"<<std::endl;
593  abort();
594  return -1;
595}
596
597/* Modify one element of packed matrix.  An element may be added.
598   This works for either ordering If the new element is zero it will be
599   deleted unless keepZero true */
600void 
601ClpMatrixBase::modifyCoefficient(int row, int column, double newElement,
602                        bool keepZero)
603{ 
604  std::cerr<<"modifyCoefficient not supported - ClpMatrixBase"<<std::endl;
605  abort();
606}
Note: See TracBrowser for help on using the repository browser.