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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 21.7 KB
Line 
1/* $Id: ClpMatrixBase.cpp 1525 2010-02-26 17:27:59Z mjs $ */
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          int i;
289          int numberBasic = number;
290          int numberColumns = model->numberColumns();
291          // Use different array so can build from true pivotVariable_
292          //int * pivotVariable = model->pivotVariable();
293          int * pivotVariable = model->rowArray(0)->getIndices();
294          for (i = 0; i < numberColumns; i++) {
295               if (model->getColumnStatus(i) == ClpSimplex::basic)
296                    pivotVariable[numberBasic++] = i;
297          }
298          number = numberBasic;
299     }
300     break;
301     // Do initial extra rows + maximum basic
302     case 2: {
303          number = model->numberRows();
304     }
305     break;
306     // To see if can dual or primal
307     case 4: {
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 * )
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
328#ifdef CLP_DEBUG
329                         check
330#endif
331                        )
332{
333     if (rhsOffset_) {
334#ifdef CLP_DEBUG
335          if (check) {
336               // no need - but check anyway
337               // zero out basic
338               int numberRows = model->numberRows();
339               int numberColumns = model->numberColumns();
340               double * solution = new double [numberColumns];
341               double * rhs = new double[numberRows];
342               const double * solutionSlack = model->solutionRegion(0);
343               CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
344               int iRow;
345               for (iRow = 0; iRow < numberRows; iRow++) {
346                    if (model->getRowStatus(iRow) != ClpSimplex::basic)
347                         rhs[iRow] = solutionSlack[iRow];
348                    else
349                         rhs[iRow] = 0.0;
350               }
351               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
352                    if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
353                         solution[iColumn] = 0.0;
354               }
355               times(-1.0, solution, rhs);
356               delete [] solution;
357               for (iRow = 0; iRow < numberRows; iRow++) {
358                    if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
359                         printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
360               }
361          }
362#endif
363          if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
364                               lastRefresh_ + refreshFrequency_)) {
365               // zero out basic
366               int numberRows = model->numberRows();
367               int numberColumns = model->numberColumns();
368               double * solution = new double [numberColumns];
369               const double * solutionSlack = model->solutionRegion(0);
370               CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
371               for (int iRow = 0; iRow < numberRows; iRow++) {
372                    if (model->getRowStatus(iRow) != ClpSimplex::basic)
373                         rhsOffset_[iRow] = solutionSlack[iRow];
374                    else
375                         rhsOffset_[iRow] = 0.0;
376               }
377               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
378                    if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
379                         solution[iColumn] = 0.0;
380               }
381               times(-1.0, solution, rhsOffset_);
382               delete [] solution;
383               lastRefresh_ = model->numberIterations();
384          }
385     }
386     return rhsOffset_;
387}
388/*
389   update information for a pivot (and effective rhs)
390*/
391int
392ClpMatrixBase::updatePivot(ClpSimplex * model, double oldInValue, double )
393{
394     if (rhsOffset_) {
395          // update effective rhs
396          int sequenceIn = model->sequenceIn();
397          int sequenceOut = model->sequenceOut();
398          double * solution = model->solutionRegion();
399          int numberColumns = model->numberColumns();
400          if (sequenceIn == sequenceOut) {
401               if (sequenceIn < numberColumns)
402                    add(model, rhsOffset_, sequenceIn, oldInValue - solution[sequenceIn]);
403          } else {
404               if (sequenceIn < numberColumns)
405                    add(model, rhsOffset_, sequenceIn, oldInValue);
406               if (sequenceOut < numberColumns)
407                    add(model, rhsOffset_, sequenceOut, -solution[sequenceOut]);
408          }
409     }
410     return 0;
411}
412int
413ClpMatrixBase::hiddenRows() const
414{
415     return 0;
416}
417/* Creates a variable.  This is called after partial pricing and may modify matrix.
418   May update bestSequence.
419*/
420void
421ClpMatrixBase::createVariable(ClpSimplex *, int &)
422{
423}
424// Returns reduced cost of a variable
425double
426ClpMatrixBase::reducedCost(ClpSimplex * model, int sequence) const
427{
428     int numberRows = model->numberRows();
429     int numberColumns = model->numberColumns();
430     if (sequence < numberRows + numberColumns)
431          return model->djRegion()[sequence];
432     else
433          return savedBestDj_;
434}
435/* Just for debug if odd type matrix.
436   Returns number and sum of primal infeasibilities.
437*/
438int
439ClpMatrixBase::checkFeasible(ClpSimplex * model, double & sum) const
440{
441     int numberRows = model->numberRows();
442     double * rhs = new double[numberRows];
443     int numberColumns = model->numberColumns();
444     int iRow;
445     CoinZeroN(rhs, numberRows);
446     times(1.0, model->solutionRegion(), rhs, model->rowScale(), model->columnScale());
447     int iColumn;
448     int logLevel = model->messageHandler()->logLevel();
449     int numberInfeasible = 0;
450     const double * rowLower = model->lowerRegion(0);
451     const double * rowUpper = model->upperRegion(0);
452     const double * solution;
453     solution = model->solutionRegion(0);
454     double tolerance = model->primalTolerance() * 1.01;
455     sum = 0.0;
456     for (iRow = 0; iRow < numberRows; iRow++) {
457          double value = rhs[iRow];
458          double value2 = solution[iRow];
459          if (logLevel > 3) {
460               if (fabs(value - value2) > 1.0e-8)
461                    printf("Row %d stored %g, computed %g\n", iRow, value2, value);
462          }
463          if (value < rowLower[iRow] - tolerance ||
464                    value > rowUpper[iRow] + tolerance) {
465               numberInfeasible++;
466               sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
467          }
468          if (value2 > rowLower[iRow] + tolerance &&
469                    value2 < rowUpper[iRow] - tolerance &&
470                    model->getRowStatus(iRow) != ClpSimplex::basic) {
471               assert (model->getRowStatus(iRow) == ClpSimplex::superBasic);
472          }
473     }
474     const double * columnLower = model->lowerRegion(1);
475     const double * columnUpper = model->upperRegion(1);
476     solution = model->solutionRegion(1);
477     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
478          double value = solution[iColumn];
479          if (value < columnLower[iColumn] - tolerance ||
480                    value > columnUpper[iColumn] + tolerance) {
481               numberInfeasible++;
482               sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
483          }
484          if (value > columnLower[iColumn] + tolerance &&
485                    value < columnUpper[iColumn] - tolerance &&
486                    model->getColumnStatus(iColumn) != ClpSimplex::basic) {
487               assert (model->getColumnStatus(iColumn) == ClpSimplex::superBasic);
488          }
489     }
490     delete [] rhs;
491     return numberInfeasible;
492}
493// These have to match ClpPrimalColumnSteepest version
494#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
495// Updates second array for steepest and does devex weights (need not be coded)
496void
497ClpMatrixBase::subsetTimes2(const ClpSimplex * model,
498                            CoinIndexedVector * dj1,
499                            const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
500                            double referenceIn, double devex,
501                            // Array for exact devex to say what is in reference framework
502                            unsigned int * reference,
503                            double * weights, double scaleFactor)
504{
505     // get subset which have nonzero tableau elements
506     subsetTransposeTimes(model, pi2, dj1, dj2);
507     bool killDjs = (scaleFactor == 0.0);
508     if (!scaleFactor)
509          scaleFactor = 1.0;
510     // columns
511
512     int number = dj1->getNumElements();
513     const int * index = dj1->getIndices();
514     double * updateBy = dj1->denseVector();
515     double * updateBy2 = dj2->denseVector();
516
517     for (int j = 0; j < number; j++) {
518          double thisWeight;
519          double pivot;
520          double pivotSquared;
521          int iSequence = index[j];
522          double value2 = updateBy[j];
523          if (killDjs)
524               updateBy[j] = 0.0;
525          double modification = updateBy2[j];
526          updateBy2[j] = 0.0;
527          ClpSimplex::Status status = model->getStatus(iSequence);
528
529          if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
530               thisWeight = weights[iSequence];
531               pivot = value2 * scaleFactor;
532               pivotSquared = pivot * pivot;
533
534               thisWeight += pivotSquared * devex + pivot * modification;
535               if (thisWeight < DEVEX_TRY_NORM) {
536                    if (referenceIn < 0.0) {
537                         // steepest
538                         thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
539                    } else {
540                         // exact
541                         thisWeight = referenceIn * pivotSquared;
542                         if (reference(iSequence))
543                              thisWeight += 1.0;
544                         thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
545                    }
546               }
547               weights[iSequence] = thisWeight;
548          }
549     }
550     dj2->setNumElements(0);
551}
552// Correct sequence in and out to give true value
553void
554ClpMatrixBase::correctSequence(const ClpSimplex * , int & , int & )
555{
556}
557// Really scale matrix
558void
559ClpMatrixBase::reallyScale(const double * , const double * )
560{
561     std::cerr << "reallyScale not supported - ClpMatrixBase" << std::endl;
562     abort();
563}
564// Updates two arrays for steepest
565void
566ClpMatrixBase::transposeTimes2(const ClpSimplex * ,
567                               const CoinIndexedVector * , CoinIndexedVector *,
568                               const CoinIndexedVector * ,
569                               CoinIndexedVector * ,
570                               double , double ,
571                               // Array for exact devex to say what is in reference framework
572                               unsigned int * ,
573                               double * , double )
574{
575     std::cerr << "transposeTimes2 not supported - ClpMatrixBase" << std::endl;
576     abort();
577}
578/* Set the dimensions of the matrix. In effect, append new empty
579   columns/rows to the matrix. A negative number for either dimension
580   means that that dimension doesn't change. Otherwise the new dimensions
581   MUST be at least as large as the current ones otherwise an exception
582   is thrown. */
583void
584ClpMatrixBase::setDimensions(int , int )
585{
586     // If odd matrix assume user knows what they are doing
587}
588/* Append a set of rows/columns to the end of the matrix. Returns number of errors
589   i.e. if any of the new rows/columns contain an index that's larger than the
590   number of columns-1/rows-1 (if numberOther>0) or duplicates
591   If 0 then rows, 1 if columns */
592int
593ClpMatrixBase::appendMatrix(int , int ,
594                            const CoinBigIndex * , const int * ,
595                            const double * , int )
596{
597     std::cerr << "appendMatrix not supported - ClpMatrixBase" << std::endl;
598     abort();
599     return -1;
600}
601
602/* Modify one element of packed matrix.  An element may be added.
603   This works for either ordering If the new element is zero it will be
604   deleted unless keepZero true */
605void
606ClpMatrixBase::modifyCoefficient(int , int , double ,
607                                 bool )
608{
609     std::cerr << "modifyCoefficient not supported - ClpMatrixBase" << std::endl;
610     abort();
611}
612#if COIN_LONG_WORK
613// For long double versions (aborts if not supported)
614void
615ClpMatrixBase::times(CoinWorkDouble scalar,
616                     const CoinWorkDouble * x, CoinWorkDouble * y) const
617{
618     std::cerr << "long times not supported - ClpMatrixBase" << std::endl;
619     abort();
620}
621void
622ClpMatrixBase::transposeTimes(CoinWorkDouble scalar,
623                              const CoinWorkDouble * x, CoinWorkDouble * y) const
624{
625     std::cerr << "long transposeTimes not supported - ClpMatrixBase" << std::endl;
626     abort();
627}
628#endif
Note: See TracBrowser for help on using the repository browser.