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

Last change on this file since 2385 was 2385, checked in by unxusr, 3 months ago

formatting

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