source: trunk/Clp/src/ClpMatrixBase.cpp @ 1402

Last change on this file since 1402 was 1402, checked in by forrest, 10 years ago

get rid of compiler warnings

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