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

Last change on this file since 2235 was 2235, checked in by forrest, 3 years ago

stuff for vector matrix

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