source: tags/move-to-subversion/ClpMatrixBase.cpp @ 1355

Last change on this file since 1355 was 744, checked in by forrest, 14 years ago

for quadratic

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