1/* $Id: ClpMatrixBase.cpp 1665 2011-01-04 17:55:54Z lou $ */
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
566void
567ClpMatrixBase::transposeTimes2(const ClpSimplex * ,
568 const CoinIndexedVector * , CoinIndexedVector *,
569 const CoinIndexedVector * ,
570 CoinIndexedVector * ,
571 double , double ,
572 // Array for exact devex to say what is in reference framework
573 unsigned int * ,
574 double * , double )
575{
576 std::cerr << "transposeTimes2 not supported - ClpMatrixBase" << std::endl;
577 abort();
578}
579/* Set the dimensions of the matrix. In effect, append new empty
580 columns/rows to the matrix. A negative number for either dimension
581 means that that dimension doesn't change. Otherwise the new dimensions
582 MUST be at least as large as the current ones otherwise an exception
583 is thrown. */
584void
585ClpMatrixBase::setDimensions(int , int )
586{
587 // If odd matrix assume user knows what they are doing
588}
589/* Append a set of rows/columns to the end of the matrix. Returns number of errors
590 i.e. if any of the new rows/columns contain an index that's larger than the
591 number of columns-1/rows-1 (if numberOther>0) or duplicates
592 If 0 then rows, 1 if columns */
593int
594ClpMatrixBase::appendMatrix(int , int ,
595 const CoinBigIndex * , const int * ,
596 const double * , int )
597{
598 std::cerr << "appendMatrix not supported - ClpMatrixBase" << std::endl;
599 abort();
600 return -1;
601}
602
603/* Modify one element of packed matrix. An element may be added.
604 This works for either ordering If the new element is zero it will be
605 deleted unless keepZero true */
606void
607ClpMatrixBase::modifyCoefficient(int , int , double ,
608 bool )
609{
610 std::cerr << "modifyCoefficient not supported - ClpMatrixBase" << std::endl;
611 abort();
612}
613#if COIN_LONG_WORK
614// For long double versions (aborts if not supported)
615void
616ClpMatrixBase::times(CoinWorkDouble scalar,
617 const CoinWorkDouble * x, CoinWorkDouble * y) const
618{
619 std::cerr << "long times not supported - ClpMatrixBase" << std::endl;
620 abort();
621}
622void
623ClpMatrixBase::transposeTimes(CoinWorkDouble scalar,
624 const CoinWorkDouble * x, CoinWorkDouble * y) const
625{
626 std::cerr << "long transposeTimes not supported - ClpMatrixBase" << std::endl;
627 abort();
628}
629#endif
630