1/* $Id: ClpNetworkMatrix.cpp 1665 2011-01-04 17:55:54Z lou $ */
2// Copyright (C) 2003, 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
7#include <cstdio>
8
9#include "CoinPragma.hpp"
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12#include "CoinPackedVector.hpp"
13
14#include "ClpSimplex.hpp"
15#include "ClpFactorization.hpp"
16// at end to get min/max!
17#include "ClpNetworkMatrix.hpp"
18#include "ClpPlusMinusOneMatrix.hpp"
19#include "ClpMessage.hpp"
20#include <iostream>
21#include <cassert>
22
23//#############################################################################
24// Constructors / Destructor / Assignment
25//#############################################################################
26
27//-------------------------------------------------------------------
28// Default Constructor
29//-------------------------------------------------------------------
30ClpNetworkMatrix::ClpNetworkMatrix ()
31 : ClpMatrixBase()
32{
33 setType(11);
34 matrix_ = NULL;
35 lengths_ = NULL;
36 indices_ = NULL;
37 numberRows_ = 0;
38 numberColumns_ = 0;
39 trueNetwork_ = false;
40}
41
42/* Constructor from two arrays */
43ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
44 const int * tail)
45 : ClpMatrixBase()
46{
47 setType(11);
48 matrix_ = NULL;
49 lengths_ = NULL;
50 indices_ = new int[2*numberColumns];
51 numberRows_ = -1;
52 numberColumns_ = numberColumns;
53 trueNetwork_ = true;
54 int iColumn;
55 CoinBigIndex j = 0;
56 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
57 int iRow = head[iColumn];
58 numberRows_ = CoinMax(numberRows_, iRow);
59 indices_[j] = iRow;
60 iRow = tail[iColumn];
61 numberRows_ = CoinMax(numberRows_, iRow);
62 indices_[j+1] = iRow;
63 }
64 numberRows_++;
65}
66//-------------------------------------------------------------------
67// Copy constructor
68//-------------------------------------------------------------------
69ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs)
70 : ClpMatrixBase(rhs)
71{
72 matrix_ = NULL;
73 lengths_ = NULL;
74 indices_ = NULL;
75 numberRows_ = rhs.numberRows_;
76 numberColumns_ = rhs.numberColumns_;
77 trueNetwork_ = rhs.trueNetwork_;
78 if (numberColumns_) {
79 indices_ = new int [ 2*numberColumns_];
80 CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
81 }
82 int numberRows = getNumRows();
83 if (rhs.rhsOffset_ && numberRows) {
84 rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
85 } else {
86 rhsOffset_ = NULL;
87 }
88}
89
90ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs)
91 : ClpMatrixBase()
92{
93 setType(11);
94 matrix_ = NULL;
95 lengths_ = NULL;
96 indices_ = NULL;
97 int iColumn;
98 assert (rhs.isColOrdered());
99 // get matrix data pointers
100 const int * row = rhs.getIndices();
101 const CoinBigIndex * columnStart = rhs.getVectorStarts();
102 const int * columnLength = rhs.getVectorLengths();
103 const double * elementByColumn = rhs.getElements();
104 numberColumns_ = rhs.getNumCols();
105 int goodNetwork = 1;
106 numberRows_ = -1;
107 indices_ = new int[2*numberColumns_];
108 CoinBigIndex j = 0;
109 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
110 CoinBigIndex k = columnStart[iColumn];
111 int iRow;
112 switch (columnLength[iColumn]) {
113 case 0:
114 goodNetwork = -1; // not classic network
115 indices_[j] = -1;
116 indices_[j+1] = -1;
117 break;
118
119 case 1:
120 goodNetwork = -1; // not classic network
121 if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
122 indices_[j] = -1;
123 iRow = row[k];
124 numberRows_ = CoinMax(numberRows_, iRow);
125 indices_[j+1] = iRow;
126 } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
127 indices_[j+1] = -1;
128 iRow = row[k];
129 numberRows_ = CoinMax(numberRows_, iRow);
130 indices_[j] = iRow;
131 } else {
132 goodNetwork = 0; // not a network
133 }
134 break;
135
136 case 2:
137 if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
138 if (fabs(elementByColumn[k+1] + 1.0) < 1.0e-10) {
139 iRow = row[k];
140 numberRows_ = CoinMax(numberRows_, iRow);
141 indices_[j+1] = iRow;
142 iRow = row[k+1];
143 numberRows_ = CoinMax(numberRows_, iRow);
144 indices_[j] = iRow;
145 } else {
146 goodNetwork = 0; // not a network
147 }
148 } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
149 if (fabs(elementByColumn[k+1] - 1.0) < 1.0e-10) {
150 iRow = row[k];
151 numberRows_ = CoinMax(numberRows_, iRow);
152 indices_[j] = iRow;
153 iRow = row[k+1];
154 numberRows_ = CoinMax(numberRows_, iRow);
155 indices_[j+1] = iRow;
156 } else {
157 goodNetwork = 0; // not a network
158 }
159 } else {
160 goodNetwork = 0; // not a network
161 }
162 break;
163
164 default:
165 goodNetwork = 0; // not a network
166 break;
167 }
168 if (!goodNetwork)
169 break;
170 }
171 if (!goodNetwork) {
172 delete [] indices_;
173 // put in message
174 printf("Not a network - can test if indices_ null\n");
175 indices_ = NULL;
176 numberRows_ = 0;
177 numberColumns_ = 0;
178 } else {
179 numberRows_ ++; // correct
180 trueNetwork_ = goodNetwork > 0;
181 }
182}
183
184//-------------------------------------------------------------------
185// Destructor
186//-------------------------------------------------------------------
187ClpNetworkMatrix::~ClpNetworkMatrix ()
188{
189 delete matrix_;
190 delete [] lengths_;
191 delete [] indices_;
192}
193
194//----------------------------------------------------------------
195// Assignment operator
196//-------------------------------------------------------------------
197ClpNetworkMatrix &
198ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
199{
200 if (this != &rhs) {
201 ClpMatrixBase::operator=(rhs);
202 delete matrix_;
203 delete [] lengths_;
204 delete [] indices_;
205 matrix_ = NULL;
206 lengths_ = NULL;
207 indices_ = NULL;
208 numberRows_ = rhs.numberRows_;
209 numberColumns_ = rhs.numberColumns_;
210 trueNetwork_ = rhs.trueNetwork_;
211 if (numberColumns_) {
212 indices_ = new int [ 2*numberColumns_];
213 CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
214 }
215 }
216 return *this;
217}
218//-------------------------------------------------------------------
219// Clone
220//-------------------------------------------------------------------
221ClpMatrixBase * ClpNetworkMatrix::clone() const
222{
223 return new ClpNetworkMatrix(*this);
224}
225
226/* Returns a new matrix in reverse order without gaps */
227ClpMatrixBase *
228ClpNetworkMatrix::reverseOrderedCopy() const
229{
230 // count number in each row
231 CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
232 CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
233 memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
234 memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
235 CoinBigIndex j = 0;
236 int i;
237 for (i = 0; i < numberColumns_; i++, j += 2) {
238 int iRow = indices_[j];
239 tempN[iRow]++;
240 iRow = indices_[j+1];
241 tempP[iRow]++;
242 }
243 int * newIndices = new int [2*numberColumns_];
244 CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
245 CoinBigIndex * newN = new CoinBigIndex[numberRows_];
246 int iRow;
247 j = 0;
248 // do starts
249 for (iRow = 0; iRow < numberRows_; iRow++) {
250 newP[iRow] = j;
251 j += tempP[iRow];
252 tempP[iRow] = newP[iRow];
253 newN[iRow] = j;
254 j += tempN[iRow];
255 tempN[iRow] = newN[iRow];
256 }
257 newP[numberRows_] = j;
258 j = 0;
259 for (i = 0; i < numberColumns_; i++, j += 2) {
260 int iRow = indices_[j];
261 CoinBigIndex put = tempN[iRow];
262 newIndices[put++] = i;
263 tempN[iRow] = put;
264 iRow = indices_[j+1];
265 put = tempP[iRow];
266 newIndices[put++] = i;
267 tempP[iRow] = put;
268 }
269 delete [] tempP;
270 delete [] tempN;
271 ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
272 newCopy->passInCopy(numberRows_, numberColumns_,
273 false, newIndices, newP, newN);
274 return newCopy;
275}
276//unscaled versions
277void
278ClpNetworkMatrix::times(double scalar,
279 const double * x, double * y) const
280{
281 int iColumn;
282 CoinBigIndex j = 0;
283 if (trueNetwork_) {
284 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
285 double value = scalar * x[iColumn];
286 if (value) {
287 int iRowM = indices_[j];
288 int iRowP = indices_[j+1];
289 y[iRowM] -= value;
290 y[iRowP] += value;
291 }
292 }
293 } else {
294 // skip negative rows
295 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
296 double value = scalar * x[iColumn];
297 if (value) {
298 int iRowM = indices_[j];
299 int iRowP = indices_[j+1];
300 if (iRowM >= 0)
301 y[iRowM] -= value;
302 if (iRowP >= 0)
303 y[iRowP] += value;
304 }
305 }
306 }
307}
308void
309ClpNetworkMatrix::transposeTimes(double scalar,
310 const double * x, double * y) const
311{
312 int iColumn;
313 CoinBigIndex j = 0;
314 if (trueNetwork_) {
315 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
316 double value = y[iColumn];
317 int iRowM = indices_[j];
318 int iRowP = indices_[j+1];
319 value -= scalar * x[iRowM];
320 value += scalar * x[iRowP];
321 y[iColumn] = value;
322 }
323 } else {
324 // skip negative rows
325 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
326 double value = y[iColumn];
327 int iRowM = indices_[j];
328 int iRowP = indices_[j+1];
329 if (iRowM >= 0)
330 value -= scalar * x[iRowM];
331 if (iRowP >= 0)
332 value += scalar * x[iRowP];
333 y[iColumn] = value;
334 }
335 }
336}
337void
338ClpNetworkMatrix::times(double scalar,
339 const double * x, double * y,
340 const double * /*rowScale*/,
341 const double * /*columnScale*/) const
342{
343 // we know it is not scaled
344 times(scalar, x, y);
345}
346void
347ClpNetworkMatrix::transposeTimes( double scalar,
348 const double * x, double * y,
349 const double * /*rowScale*/,
350 const double * /*columnScale*/,
351 double * /*spare*/) const
352{
353 // we know it is not scaled
354 transposeTimes(scalar, x, y);
355}
356/* Return <code>x * A + y</code> in <code>z</code>.
357 Squashes small elements and knows about ClpSimplex */
358void
359ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
360 const CoinIndexedVector * rowArray,
361 CoinIndexedVector * y,
362 CoinIndexedVector * columnArray) const
363{
364 // we know it is not scaled
365 columnArray->clear();
366 double * pi = rowArray->denseVector();
367 int numberNonZero = 0;
368 int * index = columnArray->getIndices();
369 double * array = columnArray->denseVector();
370 int numberInRowArray = rowArray->getNumElements();
371 // maybe I need one in OsiSimplex
372 double zeroTolerance = model->zeroTolerance();
373 int numberRows = model->numberRows();
374#ifndef NO_RTTI
375 ClpPlusMinusOneMatrix* rowCopy =
376 dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
377#else
378 ClpPlusMinusOneMatrix* rowCopy =
379 static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
380#endif
381 bool packed = rowArray->packedMode();
382 double factor = 0.3;
383 // We may not want to do by row if there may be cache problems
384 int numberColumns = model->numberColumns();
385 // It would be nice to find L2 cache size - for moment 512K
386 // Be slightly optimistic
387 if (numberColumns * sizeof(double) > 1000000) {
388 if (numberRows * 10 < numberColumns)
389 factor = 0.1;
390 else if (numberRows * 4 < numberColumns)
391 factor = 0.15;
392 else if (numberRows * 2 < numberColumns)
393 factor = 0.2;
394 //if (model->numberIterations()%50==0)
395 //printf("%d nonzero\n",numberInRowArray);
396 }
397 if (numberInRowArray > factor * numberRows || !rowCopy) {
398 // do by column
399 int iColumn;
400 assert (!y->getNumElements());
401 CoinBigIndex j = 0;
402 if (packed) {
403 // need to expand pi into y
404 assert(y->capacity() >= numberRows);
405 double * piOld = pi;
406 pi = y->denseVector();
407 const int * whichRow = rowArray->getIndices();
408 int i;
409 // modify pi so can collapse to one loop
410 for (i = 0; i < numberInRowArray; i++) {
411 int iRow = whichRow[i];
412 pi[iRow] = scalar * piOld[i];
413 }
414 if (trueNetwork_) {
415 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
416 double value = 0.0;
417 int iRowM = indices_[j];
418 int iRowP = indices_[j+1];
419 value -= pi[iRowM];
420 value += pi[iRowP];
421 if (fabs(value) > zeroTolerance) {
422 array[numberNonZero] = value;
423 index[numberNonZero++] = iColumn;
424 }
425 }
426 } else {
427 // skip negative rows
428 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
429 double value = 0.0;
430 int iRowM = indices_[j];
431 int iRowP = indices_[j+1];
432 if (iRowM >= 0)
433 value -= pi[iRowM];
434 if (iRowP >= 0)
435 value += pi[iRowP];
436 if (fabs(value) > zeroTolerance) {
437 array[numberNonZero] = value;
438 index[numberNonZero++] = iColumn;
439 }
440 }
441 }
442 for (i = 0; i < numberInRowArray; i++) {
443 int iRow = whichRow[i];
444 pi[iRow] = 0.0;
445 }
446 } else {
447 if (trueNetwork_) {
448 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
449 double value = 0.0;
450 int iRowM = indices_[j];
451 int iRowP = indices_[j+1];
452 value -= scalar * pi[iRowM];
453 value += scalar * pi[iRowP];
454 if (fabs(value) > zeroTolerance) {
455 index[numberNonZero++] = iColumn;
456 array[iColumn] = value;
457 }
458 }
459 } else {
460 // skip negative rows
461 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
462 double value = 0.0;
463 int iRowM = indices_[j];
464 int iRowP = indices_[j+1];
465 if (iRowM >= 0)
466 value -= scalar * pi[iRowM];
467 if (iRowP >= 0)
468 value += scalar * pi[iRowP];
469 if (fabs(value) > zeroTolerance) {
470 index[numberNonZero++] = iColumn;
471 array[iColumn] = value;
472 }
473 }
474 }
475 }
476 columnArray->setNumElements(numberNonZero);
477 } else {
478 // do by row
479 rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
480 }
481}
482/* Return <code>x *A in <code>z</code> but
483 just for indices in y. */
484void
485ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
486 const CoinIndexedVector * rowArray,
487 const CoinIndexedVector * y,
488 CoinIndexedVector * columnArray) const
489{
490 columnArray->clear();
491 double * pi = rowArray->denseVector();
492 double * array = columnArray->denseVector();
493 int jColumn;
494 int numberToDo = y->getNumElements();
495 const int * which = y->getIndices();
496 assert (!rowArray->packedMode());
497 columnArray->setPacked();
498 if (trueNetwork_) {
499 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
500 int iColumn = which[jColumn];
501 double value = 0.0;
502 CoinBigIndex j = iColumn << 1;
503 int iRowM = indices_[j];
504 int iRowP = indices_[j+1];
505 value -= pi[iRowM];
506 value += pi[iRowP];
507 array[jColumn] = value;
508 }
509 } else {
510 // skip negative rows
511 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
512 int iColumn = which[jColumn];
513 double value = 0.0;
514 CoinBigIndex j = iColumn << 1;
515 int iRowM = indices_[j];
516 int iRowP = indices_[j+1];
517 if (iRowM >= 0)
518 value -= pi[iRowM];
519 if (iRowP >= 0)
520 value += pi[iRowP];
521 array[jColumn] = value;
522 }
523 }
524}
525/// returns number of elements in column part of basis,
526CoinBigIndex
527ClpNetworkMatrix::countBasis( const int * whichColumn,
528 int & numberColumnBasic)
529{
530 int i;
531 CoinBigIndex numberElements = 0;
532 if (trueNetwork_) {
533 numberElements = 2 * numberColumnBasic;
534 } else {
535 for (i = 0; i < numberColumnBasic; i++) {
536 int iColumn = whichColumn[i];
537 CoinBigIndex j = iColumn << 1;
538 int iRowM = indices_[j];
539 int iRowP = indices_[j+1];
540 if (iRowM >= 0)
541 numberElements ++;
542 if (iRowP >= 0)
543 numberElements ++;
544 }
545 }
546 return numberElements;
547}
548void
549ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/,
550 const int * whichColumn,
551 int & numberColumnBasic,
552 int * indexRowU, int * start,
553 int * rowCount, int * columnCount,
554 CoinFactorizationDouble * elementU)
555{
556 int i;
557 CoinBigIndex numberElements = start[0];
558 if (trueNetwork_) {
559 for (i = 0; i < numberColumnBasic; i++) {
560 int iColumn = whichColumn[i];
561 CoinBigIndex j = iColumn << 1;
562 int iRowM = indices_[j];
563 int iRowP = indices_[j+1];
564 indexRowU[numberElements] = iRowM;
565 rowCount[iRowM]++;
566 elementU[numberElements] = -1.0;
567 indexRowU[numberElements+1] = iRowP;
568 rowCount[iRowP]++;
569 elementU[numberElements+1] = 1.0;
570 numberElements += 2;
571 start[i+1] = numberElements;
572 columnCount[i] = 2;
573 }
574 } else {
575 for (i = 0; i < numberColumnBasic; i++) {
576 int iColumn = whichColumn[i];
577 CoinBigIndex j = iColumn << 1;
578 int iRowM = indices_[j];
579 int iRowP = indices_[j+1];
580 if (iRowM >= 0) {
581 indexRowU[numberElements] = iRowM;
582 rowCount[iRowM]++;
583 elementU[numberElements++] = -1.0;
584 }
585 if (iRowP >= 0) {
586 indexRowU[numberElements] = iRowP;
587 rowCount[iRowP]++;
588 elementU[numberElements++] = 1.0;
589 }
590 start[i+1] = numberElements;
591 columnCount[i] = numberElements - start[i];
592 }
593 }
594}
595/* Unpacks a column into an CoinIndexedvector
596 */
597void
598ClpNetworkMatrix::unpack(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray,
599 int iColumn) const
600{
601 CoinBigIndex j = iColumn << 1;
602 int iRowM = indices_[j];
603 int iRowP = indices_[j+1];
604 if (iRowM >= 0)
605 rowArray->add(iRowM, -1.0);
606 if (iRowP >= 0)
607 rowArray->add(iRowP, 1.0);
608}
609/* Unpacks a column into an CoinIndexedvector
610** in packed foramt
611Note that model is NOT const. Bounds and objective could
612be modified if doing column generation (just for this variable) */
613void
614ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/,
615 CoinIndexedVector * rowArray,
616 int iColumn) const
617{
618 int * index = rowArray->getIndices();
619 double * array = rowArray->denseVector();
620 int number = 0;
621 CoinBigIndex j = iColumn << 1;
622 int iRowM = indices_[j];
623 int iRowP = indices_[j+1];
624 if (iRowM >= 0) {
625 array[number] = -1.0;
626 index[number++] = iRowM;
627 }
628 if (iRowP >= 0) {
629 array[number] = 1.0;
630 index[number++] = iRowP;
631 }
632 rowArray->setNumElements(number);
633 rowArray->setPackedMode(true);
634}
635/* Adds multiple of a column into an CoinIndexedvector
636 You can use quickAdd to add to vector */
637void
638ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray,
639 int iColumn, double multiplier) const
640{
641 CoinBigIndex j = iColumn << 1;
642 int iRowM = indices_[j];
643 int iRowP = indices_[j+1];
644 if (iRowM >= 0)
645 rowArray->quickAdd(iRowM, -multiplier);
646 if (iRowP >= 0)
647 rowArray->quickAdd(iRowP, multiplier);
648}
649/* Adds multiple of a column into an array */
650void
651ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double * array,
652 int iColumn, double multiplier) const
653{
654 CoinBigIndex j = iColumn << 1;
655 int iRowM = indices_[j];
656 int iRowP = indices_[j+1];
657 if (iRowM >= 0)
658 array[iRowM] -= multiplier;
659 if (iRowP >= 0)
660 array[iRowP] += multiplier;
661}
662
663// Return a complete CoinPackedMatrix
664CoinPackedMatrix *
665ClpNetworkMatrix::getPackedMatrix() const
666{
667 if (!matrix_) {
668 assert (trueNetwork_); // fix later
669 int numberElements = 2 * numberColumns_;
670 double * elements = new double [numberElements];
671 CoinBigIndex i;
672 for (i = 0; i < 2 * numberColumns_; i += 2) {
673 elements[i] = -1.0;
674 elements[i+1] = 1.0;
675 }
676 CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
677 for (i = 0; i < numberColumns_ + 1; i++) {
678 starts[i] = 2 * i;
679 }
680 // use assignMatrix to save space
681 delete [] lengths_;
682 lengths_ = NULL;
683 matrix_ = new CoinPackedMatrix();
684 int * indices = CoinCopyOfArray(indices_, 2 * numberColumns_);
685 matrix_->assignMatrix(true, numberRows_, numberColumns_,
686 getNumElements(),
687 elements, indices,
688 starts, lengths_);
689 assert(!elements);
690 assert(!starts);
691 assert (!indices);
692 assert (!lengths_);
693 }
694 return matrix_;
695}
696/* A vector containing the elements in the packed matrix. Note that there
697 might be gaps in this list, entries that do not belong to any
698 major-dimension vector. To get the actual elements one should look at
699 this vector together with vectorStarts and vectorLengths. */
700const double *
701ClpNetworkMatrix::getElements() const
702{
703 if (!matrix_)
704 getPackedMatrix();
705 return matrix_->getElements();
706}
707
708const CoinBigIndex *
709ClpNetworkMatrix::getVectorStarts() const
710{
711 if (!matrix_)
712 getPackedMatrix();
713 return matrix_->getVectorStarts();
714}
715/* The lengths of the major-dimension vectors. */
716const int *
717ClpNetworkMatrix::getVectorLengths() const
718{
719 assert (trueNetwork_); // fix later
720 if (!lengths_) {
721 lengths_ = new int [numberColumns_];
722 int i;
723 for (i = 0; i < numberColumns_; i++) {
724 lengths_[i] = 2;
725 }
726 }
727 return lengths_;
728}
729/* Delete the columns whose indices are listed in <code>indDel</code>. */
730void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel)
731{
732 assert (trueNetwork_);
733 int iColumn;
734 int numberBad = 0;
735 // Use array to make sure we can have duplicates
736 char * which = new char[numberColumns_];
737 memset(which, 0, numberColumns_);
738 int nDuplicate = 0;
739 for (iColumn = 0; iColumn < numDel; iColumn++) {
740 int jColumn = indDel[iColumn];
741 if (jColumn < 0 || jColumn >= numberColumns_) {
742 numberBad++;
743 } else {
744 if (which[jColumn])
745 nDuplicate++;
746 else
747 which[jColumn] = 1;
748 }
749 }
750 if (numberBad)
751 throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
752 int newNumber = numberColumns_ - numDel + nDuplicate;
753 // Get rid of temporary arrays
754 delete [] lengths_;
755 lengths_ = NULL;
756 delete matrix_;
757 matrix_ = NULL;
758 int newSize = 2 * newNumber;
759 int * newIndices = new int [newSize];
760 newSize = 0;
761 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
762 if (!which[iColumn]) {
763 CoinBigIndex start;
764 CoinBigIndex i;
765 start = 2 * iColumn;
766 for (i = start; i < start + 2; i++)
767 newIndices[newSize++] = indices_[i];
768 }
769 }
770 delete [] which;
771 delete [] indices_;
772 indices_ = newIndices;
773 numberColumns_ = newNumber;
774}
775/* Delete the rows whose indices are listed in <code>indDel</code>. */
776void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel)
777{
778 int iRow;
779 int numberBad = 0;
780 // Use array to make sure we can have duplicates
781 int * which = new int [numberRows_];
782 memset(which, 0, numberRows_ * sizeof(int));
783 for (iRow = 0; iRow < numDel; iRow++) {
784 int jRow = indDel[iRow];
785 if (jRow < 0 || jRow >= numberRows_) {
786 numberBad++;
787 } else {
788 which[jRow] = 1;
789 }
790 }
791 if (numberBad)
792 throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
793 // Only valid of all columns have 0 entries
794 int iColumn;
795 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
796 CoinBigIndex start;
797 CoinBigIndex i;
798 start = 2 * iColumn;
799 for (i = start; i < start + 2; i++) {
800 int iRow = indices_[i];
801 if (which[iRow])
802 numberBad++;
803 }
804 }
805 if (numberBad)
806 throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
807 int newNumber = 0;
808 for (iRow = 0; iRow < numberRows_; iRow++) {
809 if (!which[iRow])
810 which[iRow] = newNumber++;
811 else
812 which[iRow] = -1;
813 }
814 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
815 CoinBigIndex start;
816 CoinBigIndex i;
817 start = 2 * iColumn;
818 for (i = start; i < start + 2; i++) {
819 int iRow = indices_[i];
820 indices_[i] = which[iRow];
821 }
822 }
823 delete [] which;
824 numberRows_ = newNumber;
825}
826/* Given positive integer weights for each row fills in sum of weights
827 for each column (and slack).
828 Returns weights vector
829*/
830CoinBigIndex *
831ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
832{
833 int numberRows = model->numberRows();
834 int numberColumns = model->numberColumns();
835 int number = numberRows + numberColumns;
836 CoinBigIndex * weights = new CoinBigIndex[number];
837 int i;
838 for (i = 0; i < numberColumns; i++) {
839 CoinBigIndex j = i << 1;
840 CoinBigIndex count = 0;
841 int iRowM = indices_[j];
842 int iRowP = indices_[j+1];
843 if (iRowM >= 0) {
844 count += inputWeights[iRowM];
845 }
846 if (iRowP >= 0) {
847 count += inputWeights[iRowP];
848 }
849 weights[i] = count;
850 }
851 for (i = 0; i < numberRows; i++) {
852 weights[i+numberColumns] = inputWeights[i];
853 }
854 return weights;
855}
856/* Returns largest and smallest elements of both signs.
857 Largest refers to largest absolute value.
858*/
859void
860ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
861 double & smallestPositive, double & largestPositive)
862{
863 smallestNegative = -1.0;
864 largestNegative = -1.0;
865 smallestPositive = 1.0;
866 largestPositive = 1.0;
867}
868// Says whether it can do partial pricing
869bool
870ClpNetworkMatrix::canDoPartialPricing() const
871{
872 return true;
873}
874// Partial pricing
875void
876ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
877 int & bestSequence, int & numberWanted)
878{
879 numberWanted = currentWanted_;
880 int j;
881 int start = static_cast<int> (startFraction * numberColumns_);
882 int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_);
883 double tolerance = model->currentDualTolerance();
884 double * reducedCost = model->djRegion();
885 const double * duals = model->dualRowSolution();
886 const double * cost = model->costRegion();
887 double bestDj;
888 if (bestSequence >= 0)
889 bestDj = fabs(reducedCost[bestSequence]);
890 else
891 bestDj = tolerance;
892 int sequenceOut = model->sequenceOut();
893 int saveSequence = bestSequence;
894 if (!trueNetwork_) {
895 // Not true network
896 int iSequence;
897 for (iSequence = start; iSequence < end; iSequence++) {
898 if (iSequence != sequenceOut) {
899 double value;
900 int iRowM, iRowP;
901 ClpSimplex::Status status = model->getStatus(iSequence);
902
903 switch(status) {
904
905 case ClpSimplex::basic:
906 case ClpSimplex::isFixed:
907 break;
908 case ClpSimplex::isFree:
909 case ClpSimplex::superBasic:
910 value = cost[iSequence];
911 j = iSequence << 1;
912 // skip negative rows
913 iRowM = indices_[j];
914 iRowP = indices_[j+1];
915 if (iRowM >= 0)
916 value += duals[iRowM];
917 if (iRowP >= 0)
918 value -= duals[iRowP];
919 value = fabs(value);
920 if (value > FREE_ACCEPT * tolerance) {
921 numberWanted--;
922 // we are going to bias towards free (but only if reasonable)
923 value *= FREE_BIAS;
924 if (value > bestDj) {
925 // check flagged variable and correct dj
926 if (!model->flagged(iSequence)) {
927 bestDj = value;
928 bestSequence = iSequence;
929 } else {
930 // just to make sure we don't exit before got something
931 numberWanted++;
932 }
933 }
934 }
935 break;
936 case ClpSimplex::atUpperBound:
937 value = cost[iSequence];
938 j = iSequence << 1;
939 // skip negative rows
940 iRowM = indices_[j];
941 iRowP = indices_[j+1];
942 if (iRowM >= 0)
943 value += duals[iRowM];
944 if (iRowP >= 0)
945 value -= duals[iRowP];
946 if (value > tolerance) {
947 numberWanted--;
948 if (value > bestDj) {
949 // check flagged variable and correct dj
950 if (!model->flagged(iSequence)) {
951 bestDj = value;
952 bestSequence = iSequence;
953 } else {
954 // just to make sure we don't exit before got something
955 numberWanted++;
956 }
957 }
958 }
959 break;
960 case ClpSimplex::atLowerBound:
961 value = cost[iSequence];
962 j = iSequence << 1;
963 // skip negative rows
964 iRowM = indices_[j];
965 iRowP = indices_[j+1];
966 if (iRowM >= 0)
967 value += duals[iRowM];
968 if (iRowP >= 0)
969 value -= duals[iRowP];
970 value = -value;
971 if (value > tolerance) {
972 numberWanted--;
973 if (value > bestDj) {
974 // check flagged variable and correct dj
975 if (!model->flagged(iSequence)) {
976 bestDj = value;
977 bestSequence = iSequence;
978 } else {
979 // just to make sure we don't exit before got something
980 numberWanted++;
981 }
982 }
983 }
984 break;
985 }
986 }
987 if (!numberWanted)
988 break;
989 }
990 if (bestSequence != saveSequence) {
991 // recompute dj
992 double value = cost[bestSequence];
993 j = bestSequence << 1;
994 // skip negative rows
995 int iRowM = indices_[j];
996 int iRowP = indices_[j+1];
997 if (iRowM >= 0)
998 value += duals[iRowM];
999 if (iRowP >= 0)
1000 value -= duals[iRowP];
1001 reducedCost[bestSequence] = value;
1002 savedBestSequence_ = bestSequence;
1003 savedBestDj_ = reducedCost[savedBestSequence_];
1004 }
1005 } else {
1006 // true network
1007 int iSequence;
1008 for (iSequence = start; iSequence < end; iSequence++) {
1009 if (iSequence != sequenceOut) {
1010 double value;
1011 int iRowM, iRowP;
1012 ClpSimplex::Status status = model->getStatus(iSequence);
1013
1014 switch(status) {
1015
1016 case ClpSimplex::basic:
1017 case ClpSimplex::isFixed:
1018 break;
1019 case ClpSimplex::isFree:
1020 case ClpSimplex::superBasic:
1021 value = cost[iSequence];
1022 j = iSequence << 1;
1023 iRowM = indices_[j];
1024 iRowP = indices_[j+1];
1025 value += duals[iRowM];
1026 value -= duals[iRowP];
1027 value = fabs(value);
1028 if (value > FREE_ACCEPT * tolerance) {
1029 numberWanted--;
1030 // we are going to bias towards free (but only if reasonable)
1031 value *= FREE_BIAS;
1032 if (value > bestDj) {
1033 // check flagged variable and correct dj
1034 if (!model->flagged(iSequence)) {
1035 bestDj = value;
1036 bestSequence = iSequence;
1037 } else {
1038 // just to make sure we don't exit before got something
1039 numberWanted++;
1040 }
1041 }
1042 }
1043 break;
1044 case ClpSimplex::atUpperBound:
1045 value = cost[iSequence];
1046 j = iSequence << 1;
1047 iRowM = indices_[j];
1048 iRowP = indices_[j+1];
1049 value += duals[iRowM];
1050 value -= duals[iRowP];
1051 if (value > tolerance) {
1052 numberWanted--;
1053 if (value > bestDj) {
1054 // check flagged variable and correct dj
1055 if (!model->flagged(iSequence)) {
1056 bestDj = value;
1057 bestSequence = iSequence;
1058 } else {
1059 // just to make sure we don't exit before got something
1060 numberWanted++;
1061 }
1062 }
1063 }
1064 break;
1065 case ClpSimplex::atLowerBound:
1066 value = cost[iSequence];
1067 j = iSequence << 1;
1068 iRowM = indices_[j];
1069 iRowP = indices_[j+1];
1070 value += duals[iRowM];
1071 value -= duals[iRowP];
1072 value = -value;
1073 if (value > tolerance) {
1074 numberWanted--;
1075 if (value > bestDj) {
1076 // check flagged variable and correct dj
1077 if (!model->flagged(iSequence)) {
1078 bestDj = value;
1079 bestSequence = iSequence;
1080 } else {
1081 // just to make sure we don't exit before got something
1082 numberWanted++;
1083 }
1084 }
1085 }
1086 break;
1087 }
1088 }
1089 if (!numberWanted)
1090 break;
1091 }
1092 if (bestSequence != saveSequence) {
1093 // recompute dj
1094 double value = cost[bestSequence];
1095 j = bestSequence << 1;
1096 int iRowM = indices_[j];
1097 int iRowP = indices_[j+1];
1098 value += duals[iRowM];
1099 value -= duals[iRowP];
1100 reducedCost[bestSequence] = value;
1101 savedBestSequence_ = bestSequence;
1102 savedBestDj_ = reducedCost[savedBestSequence_];
1103 }
1104 }
1105 currentWanted_ = numberWanted;
1106}
1107// Allow any parts of a created CoinMatrix to be deleted
1108void
1109ClpNetworkMatrix::releasePackedMatrix() const
1110{
1111 delete matrix_;
1112 delete [] lengths_;
1113 matrix_ = NULL;
1114 lengths_ = NULL;
1115}
1116// Append Columns
1117void
1118ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1119{
1120 int iColumn;
1121 int numberBad = 0;
1122 for (iColumn = 0; iColumn < number; iColumn++) {
1123 int n = columns[iColumn]->getNumElements();
1124 const double * element = columns[iColumn]->getElements();
1125 if (n != 2)
1126 numberBad++;
1127 if (fabs(element[0]) != 1.0 || fabs(element[1]) != 1.0)
1128 numberBad++;
1129 else if (element[0]*element[1] != -1.0)
1130 numberBad++;
1131 }
1132 if (numberBad)
1133 throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
1134 // Get rid of temporary arrays
1135 delete [] lengths_;
1136 lengths_ = NULL;
1137 delete matrix_;
1138 matrix_ = NULL;
1139 CoinBigIndex size = 2 * number;
1140 int * temp2 = new int [numberColumns_*2+size];
1141 CoinMemcpyN(indices_, numberColumns_ * 2, temp2);
1142 delete [] indices_;
1143 indices_ = temp2;
1144 // now add
1145 size = 2 * numberColumns_;
1146 for (iColumn = 0; iColumn < number; iColumn++) {
1147 const int * row = columns[iColumn]->getIndices();
1148 const double * element = columns[iColumn]->getElements();
1149 if (element[0] == -1.0) {
1150 indices_[size++] = row[0];
1151 indices_[size++] = row[1];
1152 } else {
1153 indices_[size++] = row[1];
1154 indices_[size++] = row[0];
1155 }
1156 }
1157
1158 numberColumns_ += number;
1159}
1160// Append Rows
1161void
1162ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1163{
1164 // must be zero arrays
1165 int numberBad = 0;
1166 int iRow;
1167 for (iRow = 0; iRow < number; iRow++) {
1168 numberBad += rows[iRow]->getNumElements();
1169 }
1170 if (numberBad)
1171 throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
1172 numberRows_ += number;
1173}
1174#ifndef SLIM_CLP
1175/* Append a set of rows/columns to the end of the matrix. Returns number of errors
1176 i.e. if any of the new rows/columns contain an index that's larger than the
1177 number of columns-1/rows-1 (if numberOther>0) or duplicates
1178 If 0 then rows, 1 if columns */
1179int
1180ClpNetworkMatrix::appendMatrix(int number, int type,
1181 const CoinBigIndex * starts, const int * index,
1182 const double * element, int /*numberOther*/)
1183{
1184 int numberErrors = 0;
1185 // make into CoinPackedVector
1186 CoinPackedVectorBase ** vectors =
1187 new CoinPackedVectorBase * [number];
1188 int iVector;
1189 for (iVector = 0; iVector < number; iVector++) {
1190 int iStart = starts[iVector];
1191 vectors[iVector] =
1192 new CoinPackedVector(starts[iVector+1] - iStart,
1193 index + iStart, element + iStart);
1194 }
1195 if (type == 0) {
1196 // rows
1197 appendRows(number, vectors);
1198 } else {
1199 // columns
1200 appendCols(number, vectors);
1201 }
1202 for (iVector = 0; iVector < number; iVector++)
1203 delete vectors[iVector];
1204 delete [] vectors;
1205 return numberErrors;
1206}
1207#endif
1208/* Subset clone (without gaps). Duplicates are allowed
1209 and order is as given */
1210ClpMatrixBase *
1211ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
1212 int numberColumns,
1213 const int * whichColumns) const
1214{
1215 return new ClpNetworkMatrix(*this, numberRows, whichRows,
1216 numberColumns, whichColumns);
1217}
1218/* Subset constructor (without gaps). Duplicates are allowed
1219 and order is as given */
1220ClpNetworkMatrix::ClpNetworkMatrix (
1221 const ClpNetworkMatrix & rhs,
1222 int numberRows, const int * whichRow,
1223 int numberColumns, const int * whichColumn)
1224 : ClpMatrixBase(rhs)
1225{
1226 setType(11);
1227 matrix_ = NULL;
1228 lengths_ = NULL;
1229 indices_ = new int[2*numberColumns];
1230 numberRows_ = numberRows;
1231 numberColumns_ = numberColumns;
1232 trueNetwork_ = true;
1233 int iColumn;
1234 int numberBad = 0;
1235 int * which = new int [rhs.numberRows_];
1236 int iRow;
1237 for (iRow = 0; iRow < rhs.numberRows_; iRow++)
1238 which[iRow] = -1;
1239 int n = 0;
1240 for (iRow = 0; iRow < numberRows; iRow++) {
1241 int jRow = whichRow[iRow];
1242 assert (jRow >= 0 && jRow < rhs.numberRows_);
1243 which[jRow] = n++;
1244 }
1245 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1246 CoinBigIndex start;
1247 CoinBigIndex i;
1248 start = 2 * iColumn;
1249 CoinBigIndex offset = 2 * whichColumn[iColumn] - start;
1250 for (i = start; i < start + 2; i++) {
1251 int iRow = rhs.indices_[i+offset];
1252 iRow = which[iRow];
1253 if (iRow < 0)
1254 numberBad++;
1255 else
1256 indices_[i] = iRow;
1257 }
1258 }
1259 if (numberBad)
1260 throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
1261}
1262