1/* $Id: ClpPlusMinusOneMatrix.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 "CoinPackedVector.hpp"
12#include "CoinHelperFunctions.hpp"
13
14#include "ClpSimplex.hpp"
15#include "ClpFactorization.hpp"
16// at end to get min/max!
17#include "ClpPlusMinusOneMatrix.hpp"
18#include "ClpMessage.hpp"
19
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix ()
28 : ClpMatrixBase()
29{
30 setType(12);
31 matrix_ = NULL;
32 startPositive_ = NULL;
33 startNegative_ = NULL;
34 lengths_ = NULL;
35 indices_ = NULL;
36 numberRows_ = 0;
37 numberColumns_ = 0;
38 columnOrdered_ = true;
39}
40
41//-------------------------------------------------------------------
42// Copy constructor
43//-------------------------------------------------------------------
44ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs)
45 : ClpMatrixBase(rhs)
46{
47 matrix_ = NULL;
48 startPositive_ = NULL;
49 startNegative_ = NULL;
50 lengths_ = NULL;
51 indices_ = NULL;
52 numberRows_ = rhs.numberRows_;
53 numberColumns_ = rhs.numberColumns_;
54 columnOrdered_ = rhs.columnOrdered_;
55 if (numberColumns_) {
56 CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
57 indices_ = new int [ numberElements];
58 CoinMemcpyN(rhs.indices_, numberElements, indices_);
59 startPositive_ = new CoinBigIndex [ numberColumns_+1];
60 CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
61 startNegative_ = new CoinBigIndex [ numberColumns_];
62 CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
63 }
64 int numberRows = getNumRows();
65 if (rhs.rhsOffset_ && numberRows) {
66 rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
67 } else {
68 rhsOffset_ = NULL;
69 }
70}
71// Constructor from arrays
72ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(int numberRows, int numberColumns,
73 bool columnOrdered, const int * indices,
74 const CoinBigIndex * startPositive,
75 const CoinBigIndex * startNegative)
76 : ClpMatrixBase()
77{
78 setType(12);
79 matrix_ = NULL;
80 lengths_ = NULL;
81 numberRows_ = numberRows;
82 numberColumns_ = numberColumns;
83 columnOrdered_ = columnOrdered;
84 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
85 int numberElements = startPositive[numberMajor];
86 startPositive_ = ClpCopyOfArray(startPositive, numberMajor + 1);
87 startNegative_ = ClpCopyOfArray(startNegative, numberMajor);
88 indices_ = ClpCopyOfArray(indices, numberElements);
89 // Check valid
90 checkValid(false);
91}
92
93ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs)
94 : ClpMatrixBase()
95{
96 setType(12);
97 matrix_ = NULL;
98 startPositive_ = NULL;
99 startNegative_ = NULL;
100 lengths_ = NULL;
101 indices_ = NULL;
102 int iColumn;
103 assert (rhs.isColOrdered());
104 // get matrix data pointers
105 const int * row = rhs.getIndices();
106 const CoinBigIndex * columnStart = rhs.getVectorStarts();
107 const int * columnLength = rhs.getVectorLengths();
108 const double * elementByColumn = rhs.getElements();
109 numberColumns_ = rhs.getNumCols();
110 numberRows_ = -1;
111 indices_ = new int[rhs.getNumElements()];
112 startPositive_ = new CoinBigIndex [numberColumns_+1];
113 startNegative_ = new CoinBigIndex [numberColumns_];
114 int * temp = new int [rhs.getNumRows()];
115 CoinBigIndex j = 0;
116 CoinBigIndex numberGoodP = 0;
117 CoinBigIndex numberGoodM = 0;
118 CoinBigIndex numberBad = 0;
119 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
120 CoinBigIndex k;
121 int iNeg = 0;
122 startPositive_[iColumn] = j;
123 for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn];
124 k++) {
125 int iRow;
126 if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
127 iRow = row[k];
128 numberRows_ = CoinMax(numberRows_, iRow);
129 indices_[j++] = iRow;
130 numberGoodP++;
131 } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
132 iRow = row[k];
133 numberRows_ = CoinMax(numberRows_, iRow);
134 temp[iNeg++] = iRow;
135 numberGoodM++;
136 } else {
137 numberBad++;
138 }
139 }
140 // move negative
141 startNegative_[iColumn] = j;
142 for (k = 0; k < iNeg; k++) {
143 indices_[j++] = temp[k];
144 }
145 }
146 startPositive_[numberColumns_] = j;
147 delete [] temp;
148 if (numberBad) {
149 delete [] indices_;
150 indices_ = NULL;
151 numberRows_ = 0;
152 numberColumns_ = 0;
153 delete [] startPositive_;
154 delete [] startNegative_;
155 // Put in statistics
156 startPositive_ = new CoinBigIndex [3];
157 startPositive_[0] = numberGoodP;
158 startPositive_[1] = numberGoodM;
159 startPositive_[2] = numberBad;
160 startNegative_ = NULL;
161 } else {
162 numberRows_ ++; // correct
163 // but number should be same as rhs
164 assert (numberRows_ <= rhs.getNumRows());
165 numberRows_ = rhs.getNumRows();
166 columnOrdered_ = true;
167 }
168 // Check valid
169 if (!numberBad)
170 checkValid(false);
171}
172
173//-------------------------------------------------------------------
174// Destructor
175//-------------------------------------------------------------------
176ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
177{
178 delete matrix_;
179 delete [] startPositive_;
180 delete [] startNegative_;
181 delete [] lengths_;
182 delete [] indices_;
183}
184
185//----------------------------------------------------------------
186// Assignment operator
187//-------------------------------------------------------------------
188ClpPlusMinusOneMatrix &
189ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
190{
191 if (this != &rhs) {
192 ClpMatrixBase::operator=(rhs);
193 delete matrix_;
194 delete [] startPositive_;
195 delete [] startNegative_;
196 delete [] lengths_;
197 delete [] indices_;
198 matrix_ = NULL;
199 startPositive_ = NULL;
200 lengths_ = NULL;
201 indices_ = NULL;
202 numberRows_ = rhs.numberRows_;
203 numberColumns_ = rhs.numberColumns_;
204 columnOrdered_ = rhs.columnOrdered_;
205 if (numberColumns_) {
206 CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
207 indices_ = new int [ numberElements];
208 CoinMemcpyN(rhs.indices_, numberElements, indices_);
209 startPositive_ = new CoinBigIndex [ numberColumns_+1];
210 CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
211 startNegative_ = new CoinBigIndex [ numberColumns_];
212 CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
213 }
214 }
215 return *this;
216}
217//-------------------------------------------------------------------
218// Clone
219//-------------------------------------------------------------------
220ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
221{
222 return new ClpPlusMinusOneMatrix(*this);
223}
224/* Subset clone (without gaps). Duplicates are allowed
225 and order is as given */
226ClpMatrixBase *
227ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
228 int numberColumns,
229 const int * whichColumns) const
230{
231 return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
232 numberColumns, whichColumns);
233}
234/* Subset constructor (without gaps). Duplicates are allowed
235 and order is as given */
236ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
237 const ClpPlusMinusOneMatrix & rhs,
238 int numberRows, const int * whichRow,
239 int numberColumns, const int * whichColumn)
240 : ClpMatrixBase(rhs)
241{
242 matrix_ = NULL;
243 startPositive_ = NULL;
244 startNegative_ = NULL;
245 lengths_ = NULL;
246 indices_ = NULL;
247 numberRows_ = 0;
248 numberColumns_ = 0;
249 columnOrdered_ = rhs.columnOrdered_;
250 if (numberRows <= 0 || numberColumns <= 0) {
251 startPositive_ = new CoinBigIndex [1];
252 startPositive_[0] = 0;
253 } else {
254 numberColumns_ = numberColumns;
255 numberRows_ = numberRows;
256 const int * index1 = rhs.indices_;
257 CoinBigIndex * startPositive1 = rhs.startPositive_;
258
259 int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
260 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
261 int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
262 int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
263 // Also swap incoming if not column ordered
264 if (!columnOrdered_) {
265 int temp1 = numberRows;
266 numberRows = numberColumns;
267 numberColumns = temp1;
268 const int * temp2;
269 temp2 = whichRow;
270 whichRow = whichColumn;
271 whichColumn = temp2;
272 }
273 // Throw exception if rhs empty
274 if (numberMajor1 <= 0 || numberMinor1 <= 0)
275 throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
276 // Array to say if an old row is in new copy
277 int * newRow = new int [numberMinor1];
278 int iRow;
279 for (iRow = 0; iRow < numberMinor1; iRow++)
280 newRow[iRow] = -1;
281 // and array for duplicating rows
282 int * duplicateRow = new int [numberMinor];
283 int numberBad = 0;
284 for (iRow = 0; iRow < numberMinor; iRow++) {
285 duplicateRow[iRow] = -1;
286 int kRow = whichRow[iRow];
287 if (kRow >= 0 && kRow < numberMinor1) {
288 if (newRow[kRow] < 0) {
289 // first time
290 newRow[kRow] = iRow;
291 } else {
292 // duplicate
293 int lastRow = newRow[kRow];
294 newRow[kRow] = iRow;
295 duplicateRow[iRow] = lastRow;
296 }
297 } else {
298 // bad row
299 numberBad++;
300 }
301 }
302
303 if (numberBad)
304 throw CoinError("bad minor entries",
305 "subset constructor", "ClpPlusMinusOneMatrix");
306 // now get size and check columns
307 CoinBigIndex size = 0;
308 int iColumn;
309 numberBad = 0;
310 for (iColumn = 0; iColumn < numberMajor; iColumn++) {
311 int kColumn = whichColumn[iColumn];
312 if (kColumn >= 0 && kColumn < numberMajor1) {
313 CoinBigIndex i;
314 for (i = startPositive1[kColumn]; i < startPositive1[kColumn+1]; i++) {
315 int kRow = index1[i];
316 kRow = newRow[kRow];
317 while (kRow >= 0) {
318 size++;
319 kRow = duplicateRow[kRow];
320 }
321 }
322 } else {
323 // bad column
324 numberBad++;
325 printf("%d %d %d %d\n", iColumn, numberMajor, numberMajor1, kColumn);
326 }
327 }
328 if (numberBad)
329 throw CoinError("bad major entries",
330 "subset constructor", "ClpPlusMinusOneMatrix");
331 // now create arrays
332 startPositive_ = new CoinBigIndex [numberMajor+1];
333 startNegative_ = new CoinBigIndex [numberMajor];
334 indices_ = new int[size];
335 // and fill them
336 size = 0;
337 startPositive_[0] = 0;
338 CoinBigIndex * startNegative1 = rhs.startNegative_;
339 for (iColumn = 0; iColumn < numberMajor; iColumn++) {
340 int kColumn = whichColumn[iColumn];
341 CoinBigIndex i;
342 for (i = startPositive1[kColumn]; i < startNegative1[kColumn]; i++) {
343 int kRow = index1[i];
344 kRow = newRow[kRow];
345 while (kRow >= 0) {
346 indices_[size++] = kRow;
347 kRow = duplicateRow[kRow];
348 }
349 }
350 startNegative_[iColumn] = size;
351 for (; i < startPositive1[kColumn+1]; i++) {
352 int kRow = index1[i];
353 kRow = newRow[kRow];
354 while (kRow >= 0) {
355 indices_[size++] = kRow;
356 kRow = duplicateRow[kRow];
357 }
358 }
359 startPositive_[iColumn+1] = size;
360 }
361 delete [] newRow;
362 delete [] duplicateRow;
363 }
364 // Check valid
365 checkValid(false);
366}
367
368
369/* Returns a new matrix in reverse order without gaps */
370ClpMatrixBase *
371ClpPlusMinusOneMatrix::reverseOrderedCopy() const
372{
373 int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
374 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
375 // count number in each row/column
376 CoinBigIndex * tempP = new CoinBigIndex [numberMinor];
377 CoinBigIndex * tempN = new CoinBigIndex [numberMinor];
378 memset(tempP, 0, numberMinor * sizeof(CoinBigIndex));
379 memset(tempN, 0, numberMinor * sizeof(CoinBigIndex));
380 CoinBigIndex j = 0;
381 int i;
382 for (i = 0; i < numberMajor; i++) {
383 for (; j < startNegative_[i]; j++) {
384 int iRow = indices_[j];
385 tempP[iRow]++;
386 }
387 for (; j < startPositive_[i+1]; j++) {
388 int iRow = indices_[j];
389 tempN[iRow]++;
390 }
391 }
392 int * newIndices = new int [startPositive_[numberMajor]];
393 CoinBigIndex * newP = new CoinBigIndex [numberMinor+1];
394 CoinBigIndex * newN = new CoinBigIndex[numberMinor];
395 int iRow;
396 j = 0;
397 // do starts
398 for (iRow = 0; iRow < numberMinor; iRow++) {
399 newP[iRow] = j;
400 j += tempP[iRow];
401 tempP[iRow] = newP[iRow];
402 newN[iRow] = j;
403 j += tempN[iRow];
404 tempN[iRow] = newN[iRow];
405 }
406 newP[numberMinor] = j;
407 j = 0;
408 for (i = 0; i < numberMajor; i++) {
409 for (; j < startNegative_[i]; j++) {
410 int iRow = indices_[j];
411 CoinBigIndex put = tempP[iRow];
412 newIndices[put++] = i;
413 tempP[iRow] = put;
414 }
415 for (; j < startPositive_[i+1]; j++) {
416 int iRow = indices_[j];
417 CoinBigIndex put = tempN[iRow];
418 newIndices[put++] = i;
419 tempN[iRow] = put;
420 }
421 }
422 delete [] tempP;
423 delete [] tempN;
424 ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
425 newCopy->passInCopy(numberMinor, numberMajor,
426 !columnOrdered_, newIndices, newP, newN);
427 return newCopy;
428}
429//unscaled versions
430void
431ClpPlusMinusOneMatrix::times(double scalar,
432 const double * x, double * y) const
433{
434 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
435 int i;
436 CoinBigIndex j;
437 assert (columnOrdered_);
438 for (i = 0; i < numberMajor; i++) {
439 double value = scalar * x[i];
440 if (value) {
441 for (j = startPositive_[i]; j < startNegative_[i]; j++) {
442 int iRow = indices_[j];
443 y[iRow] += value;
444 }
445 for (; j < startPositive_[i+1]; j++) {
446 int iRow = indices_[j];
447 y[iRow] -= value;
448 }
449 }
450 }
451}
452void
453ClpPlusMinusOneMatrix::transposeTimes(double scalar,
454 const double * x, double * y) const
455{
456 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
457 int i;
458 CoinBigIndex j = 0;
459 assert (columnOrdered_);
460 for (i = 0; i < numberMajor; i++) {
461 double value = 0.0;
462 for (; j < startNegative_[i]; j++) {
463 int iRow = indices_[j];
464 value += x[iRow];
465 }
466 for (; j < startPositive_[i+1]; j++) {
467 int iRow = indices_[j];
468 value -= x[iRow];
469 }
470 y[i] += scalar * value;
471 }
472}
473void
474ClpPlusMinusOneMatrix::times(double scalar,
475 const double * x, double * y,
476 const double * /*rowScale*/,
477 const double * /*columnScale*/) const
478{
479 // we know it is not scaled
480 times(scalar, x, y);
481}
482void
483ClpPlusMinusOneMatrix::transposeTimes( double scalar,
484 const double * x, double * y,
485 const double * /*rowScale*/,
486 const double * /*columnScale*/,
487 double * /*spare*/) const
488{
489 // we know it is not scaled
490 transposeTimes(scalar, x, y);
491}
492/* Return <code>x * A + y</code> in <code>z</code>.
493 Squashes small elements and knows about ClpSimplex */
494void
495ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
496 const CoinIndexedVector * rowArray,
497 CoinIndexedVector * y,
498 CoinIndexedVector * columnArray) const
499{
500 // we know it is not scaled
501 columnArray->clear();
502 double * pi = rowArray->denseVector();
503 int numberNonZero = 0;
504 int * index = columnArray->getIndices();
505 double * array = columnArray->denseVector();
506 int numberInRowArray = rowArray->getNumElements();
507 // maybe I need one in OsiSimplex
508 double zeroTolerance = model->zeroTolerance();
509 int numberRows = model->numberRows();
510 bool packed = rowArray->packedMode();
511#ifndef NO_RTTI
512 ClpPlusMinusOneMatrix* rowCopy =
513 dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
514#else
515 ClpPlusMinusOneMatrix* rowCopy =
516 static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
517#endif
518 double factor = 0.3;
519 // We may not want to do by row if there may be cache problems
520 int numberColumns = model->numberColumns();
521 // It would be nice to find L2 cache size - for moment 512K
522 // Be slightly optimistic
523 if (numberColumns * sizeof(double) > 1000000) {
524 if (numberRows * 10 < numberColumns)
525 factor = 0.1;
526 else if (numberRows * 4 < numberColumns)
527 factor = 0.15;
528 else if (numberRows * 2 < numberColumns)
529 factor = 0.2;
530 }
531 if (numberInRowArray > factor * numberRows || !rowCopy) {
532 assert (!y->getNumElements());
533 // do by column
534 // Need to expand if packed mode
535 int iColumn;
536 CoinBigIndex j = 0;
537 assert (columnOrdered_);
538 if (packed) {
539 // need to expand pi into y
540 assert(y->capacity() >= numberRows);
541 double * piOld = pi;
542 pi = y->denseVector();
543 const int * whichRow = rowArray->getIndices();
544 int i;
545 // modify pi so can collapse to one loop
546 for (i = 0; i < numberInRowArray; i++) {
547 int iRow = whichRow[i];
548 pi[iRow] = scalar * piOld[i];
549 }
550 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
551 double value = 0.0;
552 for (; j < startNegative_[iColumn]; j++) {
553 int iRow = indices_[j];
554 value += pi[iRow];
555 }
556 for (; j < startPositive_[iColumn+1]; j++) {
557 int iRow = indices_[j];
558 value -= pi[iRow];
559 }
560 if (fabs(value) > zeroTolerance) {
561 array[numberNonZero] = value;
562 index[numberNonZero++] = iColumn;
563 }
564 }
565 for (i = 0; i < numberInRowArray; i++) {
566 int iRow = whichRow[i];
567 pi[iRow] = 0.0;
568 }
569 } else {
570 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
571 double value = 0.0;
572 for (; j < startNegative_[iColumn]; j++) {
573 int iRow = indices_[j];
574 value += pi[iRow];
575 }
576 for (; j < startPositive_[iColumn+1]; j++) {
577 int iRow = indices_[j];
578 value -= pi[iRow];
579 }
580 value *= scalar;
581 if (fabs(value) > zeroTolerance) {
582 index[numberNonZero++] = iColumn;
583 array[iColumn] = value;
584 }
585 }
586 }
587 columnArray->setNumElements(numberNonZero);
588 } else {
589 // do by row
590 rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
591 }
592}
593/* Return <code>x * A + y</code> in <code>z</code>.
594 Squashes small elements and knows about ClpSimplex */
595void
596ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
597 const CoinIndexedVector * rowArray,
598 CoinIndexedVector * y,
599 CoinIndexedVector * columnArray) const
600{
601 columnArray->clear();
602 double * pi = rowArray->denseVector();
603 int numberNonZero = 0;
604 int * index = columnArray->getIndices();
605 double * array = columnArray->denseVector();
606 int numberInRowArray = rowArray->getNumElements();
607 // maybe I need one in OsiSimplex
608 double zeroTolerance = model->zeroTolerance();
609 const int * column = indices_;
610 const CoinBigIndex * startPositive = startPositive_;
611 const CoinBigIndex * startNegative = startNegative_;
612 const int * whichRow = rowArray->getIndices();
613 bool packed = rowArray->packedMode();
614 if (numberInRowArray > 2) {
615 // do by rows
616 int iRow;
617 double * markVector = y->denseVector(); // probably empty .. but
618 int numberOriginal = 0;
619 int i;
620 if (packed) {
621 numberNonZero = 0;
622 // and set up mark as char array
623 char * marked = reinterpret_cast<char *> (index + columnArray->capacity());
624 double * array2 = y->denseVector();
625#ifdef CLP_DEBUG
626 int numberColumns = model->numberColumns();
627 for (i = 0; i < numberColumns; i++) {
628 assert(!marked[i]);
629 assert(!array2[i]);
630 }
631#endif
632 for (i = 0; i < numberInRowArray; i++) {
633 iRow = whichRow[i];
634 double value = pi[i] * scalar;
635 CoinBigIndex j;
636 for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
637 int iColumn = column[j];
638 if (!marked[iColumn]) {
639 marked[iColumn] = 1;
640 index[numberNonZero++] = iColumn;
641 }
642 array2[iColumn] += value;
643 }
644 for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
645 int iColumn = column[j];
646 if (!marked[iColumn]) {
647 marked[iColumn] = 1;
648 index[numberNonZero++] = iColumn;
649 }
650 array2[iColumn] -= value;
651 }
652 }
653 // get rid of tiny values and zero out marked
654 numberOriginal = numberNonZero;
655 numberNonZero = 0;
656 for (i = 0; i < numberOriginal; i++) {
657 int iColumn = index[i];
658 if (marked[iColumn]) {
659 double value = array2[iColumn];
660 array2[iColumn] = 0.0;
661 marked[iColumn] = 0;
662 if (fabs(value) > zeroTolerance) {
663 array[numberNonZero] = value;
664 index[numberNonZero++] = iColumn;
665 }
666 }
667 }
668 } else {
669 numberNonZero = 0;
670 // and set up mark as char array
671 char * marked = reinterpret_cast<char *> (markVector);
672 for (i = 0; i < numberOriginal; i++) {
673 int iColumn = index[i];
674 marked[iColumn] = 0;
675 }
676 for (i = 0; i < numberInRowArray; i++) {
677 iRow = whichRow[i];
678 double value = pi[iRow] * scalar;
679 CoinBigIndex j;
680 for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
681 int iColumn = column[j];
682 if (!marked[iColumn]) {
683 marked[iColumn] = 1;
684 index[numberNonZero++] = iColumn;
685 }
686 array[iColumn] += value;
687 }
688 for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
689 int iColumn = column[j];
690 if (!marked[iColumn]) {
691 marked[iColumn] = 1;
692 index[numberNonZero++] = iColumn;
693 }
694 array[iColumn] -= value;
695 }
696 }
697 // get rid of tiny values and zero out marked
698 numberOriginal = numberNonZero;
699 numberNonZero = 0;
700 for (i = 0; i < numberOriginal; i++) {
701 int iColumn = index[i];
702 marked[iColumn] = 0;
703 if (fabs(array[iColumn]) > zeroTolerance) {
704 index[numberNonZero++] = iColumn;
705 } else {
706 array[iColumn] = 0.0;
707 }
708 }
709 }
710 } else if (numberInRowArray == 2) {
711 /* do by rows when two rows (do longer first when not packed
712 and shorter first if packed */
713 int iRow0 = whichRow[0];
714 int iRow1 = whichRow[1];
715 CoinBigIndex j;
716 if (packed) {
717 double pi0 = pi[0];
718 double pi1 = pi[1];
719 if (startPositive[iRow0+1] - startPositive[iRow0] >
720 startPositive[iRow1+1] - startPositive[iRow1]) {
721 int temp = iRow0;
722 iRow0 = iRow1;
723 iRow1 = temp;
724 pi0 = pi1;
725 pi1 = pi[0];
726 }
727 // and set up mark as char array
728 char * marked = reinterpret_cast<char *> (index + columnArray->capacity());
729 int * lookup = y->getIndices();
730 double value = pi0 * scalar;
731 for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
732 int iColumn = column[j];
733 array[numberNonZero] = value;
734 marked[iColumn] = 1;
735 lookup[iColumn] = numberNonZero;
736 index[numberNonZero++] = iColumn;
737 }
738 for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
739 int iColumn = column[j];
740 array[numberNonZero] = -value;
741 marked[iColumn] = 1;
742 lookup[iColumn] = numberNonZero;
743 index[numberNonZero++] = iColumn;
744 }
745 int numberOriginal = numberNonZero;
746 value = pi1 * scalar;
747 for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
748 int iColumn = column[j];
749 if (marked[iColumn]) {
750 int iLookup = lookup[iColumn];
751 array[iLookup] += value;
752 } else {
753 if (fabs(value) > zeroTolerance) {
754 array[numberNonZero] = value;
755 index[numberNonZero++] = iColumn;
756 }
757 }
758 }
759 for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
760 int iColumn = column[j];
761 if (marked[iColumn]) {
762 int iLookup = lookup[iColumn];
763 array[iLookup] -= value;
764 } else {
765 if (fabs(value) > zeroTolerance) {
766 array[numberNonZero] = -value;
767 index[numberNonZero++] = iColumn;
768 }
769 }
770 }
771 // get rid of tiny values and zero out marked
772 int nDelete = 0;
773 for (j = 0; j < numberOriginal; j++) {
774 int iColumn = index[j];
775 marked[iColumn] = 0;
776 if (fabs(array[j]) <= zeroTolerance)
777 nDelete++;
778 }
779 if (nDelete) {
780 numberOriginal = numberNonZero;
781 numberNonZero = 0;
782 for (j = 0; j < numberOriginal; j++) {
783 int iColumn = index[j];
784 double value = array[j];
785 array[j] = 0.0;
786 if (fabs(value) > zeroTolerance) {
787 array[numberNonZero] = value;
788 index[numberNonZero++] = iColumn;
789 }
790 }
791 }
792 } else {
793 if (startPositive[iRow0+1] - startPositive[iRow0] <
794 startPositive[iRow1+1] - startPositive[iRow1]) {
795 int temp = iRow0;
796 iRow0 = iRow1;
797 iRow1 = temp;
798 }
799 int numberOriginal;
800 int i;
801 numberNonZero = 0;
802 double value;
803 value = pi[iRow0] * scalar;
804 CoinBigIndex j;
805 for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
806 int iColumn = column[j];
807 index[numberNonZero++] = iColumn;
808 array[iColumn] = value;
809 }
810 for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
811 int iColumn = column[j];
812 index[numberNonZero++] = iColumn;
813 array[iColumn] = -value;
814 }
815 value = pi[iRow1] * scalar;
816 for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
817 int iColumn = column[j];
818 double value2 = array[iColumn];
819 if (value2) {
820 value2 += value;
821 } else {
822 value2 = value;
823 index[numberNonZero++] = iColumn;
824 }
825 array[iColumn] = value2;
826 }
827 for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
828 int iColumn = column[j];
829 double value2 = array[iColumn];
830 if (value2) {
831 value2 -= value;
832 } else {
833 value2 = -value;
834 index[numberNonZero++] = iColumn;
835 }
836 array[iColumn] = value2;
837 }
838 // get rid of tiny values and zero out marked
839 numberOriginal = numberNonZero;
840 numberNonZero = 0;
841 for (i = 0; i < numberOriginal; i++) {
842 int iColumn = index[i];
843 if (fabs(array[iColumn]) > zeroTolerance) {
844 index[numberNonZero++] = iColumn;
845 } else {
846 array[iColumn] = 0.0;
847 }
848 }
849 }
850 } else if (numberInRowArray == 1) {
851 // Just one row
852 int iRow = rowArray->getIndices()[0];
853 numberNonZero = 0;
854 double value;
855 iRow = whichRow[0];
856 CoinBigIndex j;
857 if (packed) {
858 value = pi[0] * scalar;
859 if (fabs(value) > zeroTolerance) {
860 for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
861 int iColumn = column[j];
862 array[numberNonZero] = value;
863 index[numberNonZero++] = iColumn;
864 }
865 for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
866 int iColumn = column[j];
867 array[numberNonZero] = -value;
868 index[numberNonZero++] = iColumn;
869 }
870 }
871 } else {
872 value = pi[iRow] * scalar;
873 if (fabs(value) > zeroTolerance) {
874 for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
875 int iColumn = column[j];
876 array[iColumn] = value;
877 index[numberNonZero++] = iColumn;
878 }
879 for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
880 int iColumn = column[j];
881 array[iColumn] = -value;
882 index[numberNonZero++] = iColumn;
883 }
884 }
885 }
886 }
887 columnArray->setNumElements(numberNonZero);
888 if (packed)
889 columnArray->setPacked();
890 y->setNumElements(0);
891}
892/* Return <code>x *A in <code>z</code> but
893 just for indices in y. */
894void
895ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * ,
896 const CoinIndexedVector * rowArray,
897 const CoinIndexedVector * y,
898 CoinIndexedVector * columnArray) const
899{
900 columnArray->clear();
901 double * pi = rowArray->denseVector();
902 double * array = columnArray->denseVector();
903 int jColumn;
904 int numberToDo = y->getNumElements();
905 const int * which = y->getIndices();
906 assert (!rowArray->packedMode());
907 columnArray->setPacked();
908 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
909 int iColumn = which[jColumn];
910 double value = 0.0;
911 CoinBigIndex j = startPositive_[iColumn];
912 for (; j < startNegative_[iColumn]; j++) {
913 int iRow = indices_[j];
914 value += pi[iRow];
915 }
916 for (; j < startPositive_[iColumn+1]; j++) {
917 int iRow = indices_[j];
918 value -= pi[iRow];
919 }
920 array[jColumn] = value;
921 }
922}
923/// returns number of elements in column part of basis,
924CoinBigIndex
925ClpPlusMinusOneMatrix::countBasis(const int * whichColumn,
926 int & numberColumnBasic)
927{
928 int i;
929 CoinBigIndex numberElements = 0;
930 for (i = 0; i < numberColumnBasic; i++) {
931 int iColumn = whichColumn[i];
932 numberElements += startPositive_[iColumn+1] - startPositive_[iColumn];
933 }
934 return numberElements;
935}
936void
937ClpPlusMinusOneMatrix::fillBasis(ClpSimplex * ,
938 const int * whichColumn,
939 int & numberColumnBasic,
940 int * indexRowU, int * start,
941 int * rowCount, int * columnCount,
942 CoinFactorizationDouble * elementU)
943{
944 int i;
945 CoinBigIndex numberElements = start[0];
946 assert (columnOrdered_);
947 for (i = 0; i < numberColumnBasic; i++) {
948 int iColumn = whichColumn[i];
949 CoinBigIndex j = startPositive_[iColumn];
950 for (; j < startNegative_[iColumn]; j++) {
951 int iRow = indices_[j];
952 indexRowU[numberElements] = iRow;
953 rowCount[iRow]++;
954 elementU[numberElements++] = 1.0;
955 }
956 for (; j < startPositive_[iColumn+1]; j++) {
957 int iRow = indices_[j];
958 indexRowU[numberElements] = iRow;
959 rowCount[iRow]++;
960 elementU[numberElements++] = -1.0;
961 }
962 start[i+1] = numberElements;
963 columnCount[i] = numberElements - start[i];
964 }
965}
966/* Unpacks a column into an CoinIndexedvector
967 */
968void
969ClpPlusMinusOneMatrix::unpack(const ClpSimplex * ,
970 CoinIndexedVector * rowArray,
971 int iColumn) const
972{
973 CoinBigIndex j = startPositive_[iColumn];
974 for (; j < startNegative_[iColumn]; j++) {
975 int iRow = indices_[j];
976 rowArray->add(iRow, 1.0);
977 }
978 for (; j < startPositive_[iColumn+1]; j++) {
979 int iRow = indices_[j];
980 rowArray->add(iRow, -1.0);
981 }
982}
983/* Unpacks a column into an CoinIndexedvector
984** in packed foramt
985Note that model is NOT const. Bounds and objective could
986be modified if doing column generation (just for this variable) */
987void
988ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * ,
989 CoinIndexedVector * rowArray,
990 int iColumn) const
991{
992 int * index = rowArray->getIndices();
993 double * array = rowArray->denseVector();
994 int number = 0;
995 CoinBigIndex j = startPositive_[iColumn];
996 for (; j < startNegative_[iColumn]; j++) {
997 int iRow = indices_[j];
998 array[number] = 1.0;
999 index[number++] = iRow;
1000 }
1001 for (; j < startPositive_[iColumn+1]; j++) {
1002 int iRow = indices_[j];
1003 array[number] = -1.0;
1004 index[number++] = iRow;
1005 }
1006 rowArray->setNumElements(number);
1007 rowArray->setPackedMode(true);
1008}
1009/* Adds multiple of a column into an CoinIndexedvector
1010 You can use quickAdd to add to vector */
1011void
1012ClpPlusMinusOneMatrix::add(const ClpSimplex * , CoinIndexedVector * rowArray,
1013 int iColumn, double multiplier) const
1014{
1015 CoinBigIndex j = startPositive_[iColumn];
1016 for (; j < startNegative_[iColumn]; j++) {
1017 int iRow = indices_[j];
1018 rowArray->quickAdd(iRow, multiplier);
1019 }
1020 for (; j < startPositive_[iColumn+1]; j++) {
1021 int iRow = indices_[j];
1022 rowArray->quickAdd(iRow, -multiplier);
1023 }
1024}
1025/* Adds multiple of a column into an array */
1026void
1027ClpPlusMinusOneMatrix::add(const ClpSimplex * , double * array,
1028 int iColumn, double multiplier) const
1029{
1030 CoinBigIndex j = startPositive_[iColumn];
1031 for (; j < startNegative_[iColumn]; j++) {
1032 int iRow = indices_[j];
1033 array[iRow] += multiplier;
1034 }
1035 for (; j < startPositive_[iColumn+1]; j++) {
1036 int iRow = indices_[j];
1037 array[iRow] -= multiplier;
1038 }
1039}
1040
1041// Return a complete CoinPackedMatrix
1042CoinPackedMatrix *
1043ClpPlusMinusOneMatrix::getPackedMatrix() const
1044{
1045 if (!matrix_) {
1046 int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
1047 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1048 int numberElements = startPositive_[numberMajor];
1049 double * elements = new double [numberElements];
1050 CoinBigIndex j = 0;
1051 int i;
1052 for (i = 0; i < numberMajor; i++) {
1053 for (; j < startNegative_[i]; j++) {
1054 elements[j] = 1.0;
1055 }
1056 for (; j < startPositive_[i+1]; j++) {
1057 elements[j] = -1.0;
1058 }
1059 }
1060 matrix_ = new CoinPackedMatrix(columnOrdered_, numberMinor, numberMajor,
1061 getNumElements(),
1062 elements, indices_,
1063 startPositive_, getVectorLengths());
1064 delete [] elements;
1065 delete [] lengths_;
1066 lengths_ = NULL;
1067 }
1068 return matrix_;
1069}
1070/* A vector containing the elements in the packed matrix. Note that there
1071 might be gaps in this list, entries that do not belong to any
1072 major-dimension vector. To get the actual elements one should look at
1073 this vector together with vectorStarts and vectorLengths. */
1074const double *
1075ClpPlusMinusOneMatrix::getElements() const
1076{
1077 if (!matrix_)
1078 getPackedMatrix();
1079 return matrix_->getElements();
1080}
1081
1082const CoinBigIndex *
1083ClpPlusMinusOneMatrix::getVectorStarts() const
1084{
1085 return startPositive_;
1086}
1087/* The lengths of the major-dimension vectors. */
1088const int *
1089ClpPlusMinusOneMatrix::getVectorLengths() const
1090{
1091 if (!lengths_) {
1092 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1093 lengths_ = new int [numberMajor];
1094 int i;
1095 for (i = 0; i < numberMajor; i++) {
1096 lengths_[i] = startPositive_[i+1] - startPositive_[i];
1097 }
1098 }
1099 return lengths_;
1100}
1101/* Delete the columns whose indices are listed in <code>indDel</code>. */
1102void
1103ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
1104{
1105 int iColumn;
1106 CoinBigIndex newSize = startPositive_[numberColumns_];
1107 int numberBad = 0;
1108 // Use array to make sure we can have duplicates
1109 int * which = new int[numberColumns_];
1110 memset(which, 0, numberColumns_ * sizeof(int));
1111 int nDuplicate = 0;
1112 for (iColumn = 0; iColumn < numDel; iColumn++) {
1113 int jColumn = indDel[iColumn];
1114 if (jColumn < 0 || jColumn >= numberColumns_) {
1115 numberBad++;
1116 } else {
1117 newSize -= startPositive_[jColumn+1] - startPositive_[jColumn];
1118 if (which[jColumn])
1119 nDuplicate++;
1120 else
1121 which[jColumn] = 1;
1122 }
1123 }
1124 if (numberBad)
1125 throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
1126 int newNumber = numberColumns_ - numDel + nDuplicate;
1127 // Get rid of temporary arrays
1128 delete [] lengths_;
1129 lengths_ = NULL;
1130 delete matrix_;
1131 matrix_ = NULL;
1132 CoinBigIndex * newPositive = new CoinBigIndex [newNumber+1];
1133 CoinBigIndex * newNegative = new CoinBigIndex [newNumber];
1134 int * newIndices = new int [newSize];
1135 newNumber = 0;
1136 newSize = 0;
1137 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1138 if (!which[iColumn]) {
1139 CoinBigIndex start, end;
1140 CoinBigIndex i;
1141 start = startPositive_[iColumn];
1142 end = startNegative_[iColumn];
1143 newPositive[newNumber] = newSize;
1144 for (i = start; i < end; i++)
1145 newIndices[newSize++] = indices_[i];
1146 start = startNegative_[iColumn];
1147 end = startPositive_[iColumn+1];
1148 newNegative[newNumber++] = newSize;
1149 for (i = start; i < end; i++)
1150 newIndices[newSize++] = indices_[i];
1151 }
1152 }
1153 newPositive[newNumber] = newSize;
1154 delete [] which;
1155 delete [] startPositive_;
1156 startPositive_ = newPositive;
1157 delete [] startNegative_;
1158 startNegative_ = newNegative;
1159 delete [] indices_;
1160 indices_ = newIndices;
1161 numberColumns_ = newNumber;
1162}
1163/* Delete the rows whose indices are listed in <code>indDel</code>. */
1164void
1165ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
1166{
1167 int iRow;
1168 int numberBad = 0;
1169 // Use array to make sure we can have duplicates
1170 int * which = new int[numberRows_];
1171 memset(which, 0, numberRows_ * sizeof(int));
1172 int nDuplicate = 0;
1173 for (iRow = 0; iRow < numDel; iRow++) {
1174 int jRow = indDel[iRow];
1175 if (jRow < 0 || jRow >= numberRows_) {
1176 numberBad++;
1177 } else {
1178 if (which[jRow])
1179 nDuplicate++;
1180 else
1181 which[jRow] = 1;
1182 }
1183 }
1184 if (numberBad)
1185 throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix");
1186 CoinBigIndex iElement;
1187 CoinBigIndex numberElements = startPositive_[numberColumns_];
1188 CoinBigIndex newSize = 0;
1189 for (iElement = 0; iElement < numberElements; iElement++) {
1190 iRow = indices_[iElement];
1191 if (!which[iRow])
1192 newSize++;
1193 }
1194 int newNumber = numberRows_ - numDel + nDuplicate;
1195 // Get rid of temporary arrays
1196 delete [] lengths_;
1197 lengths_ = NULL;
1198 delete matrix_;
1199 matrix_ = NULL;
1200 int * newIndices = new int [newSize];
1201 newSize = 0;
1202 int iColumn;
1203 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1204 CoinBigIndex start, end;
1205 CoinBigIndex i;
1206 start = startPositive_[iColumn];
1207 end = startNegative_[iColumn];
1208 startPositive_[newNumber] = newSize;
1209 for (i = start; i < end; i++) {
1210 iRow = indices_[i];
1211 if (!which[iRow])
1212 newIndices[newSize++] = iRow;
1213 }
1214 start = startNegative_[iColumn];
1215 end = startPositive_[iColumn+1];
1216 startNegative_[newNumber] = newSize;
1217 for (i = start; i < end; i++) {
1218 iRow = indices_[i];
1219 if (!which[iRow])
1220 newIndices[newSize++] = iRow;
1221 }
1222 }
1223 startPositive_[numberColumns_] = newSize;
1224 delete [] which;
1225 delete [] indices_;
1226 indices_ = newIndices;
1227 numberRows_ = newNumber;
1228}
1229bool
1230ClpPlusMinusOneMatrix::isColOrdered() const
1231{
1232 return columnOrdered_;
1233}
1234/* Number of entries in the packed matrix. */
1235CoinBigIndex
1236ClpPlusMinusOneMatrix::getNumElements() const
1237{
1238 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1239 if (startPositive_)
1240 return startPositive_[numberMajor];
1241 else
1242 return 0;
1243}
1244// pass in copy (object takes ownership)
1245void
1246ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
1247 bool columnOrdered, int * indices,
1248 CoinBigIndex * startPositive, CoinBigIndex * startNegative)
1249{
1250 columnOrdered_ = columnOrdered;
1251 startPositive_ = startPositive;
1252 startNegative_ = startNegative;
1253 indices_ = indices;
1254 numberRows_ = numberRows;
1255 numberColumns_ = numberColumns;
1256 // Check valid
1257 checkValid(false);
1258}
1259// Just checks matrix valid - will say if dimensions not quite right if detail
1260void
1261ClpPlusMinusOneMatrix::checkValid(bool detail) const
1262{
1263 int maxIndex = -1;
1264 int minIndex = columnOrdered_ ? numberRows_ : numberColumns_;
1265 int number = !columnOrdered_ ? numberRows_ : numberColumns_;
1266 int numberElements = getNumElements();
1267 CoinBigIndex last = -1;
1268 int bad = 0;
1269 for (int i = 0; i < number; i++) {
1270 if(startPositive_[i] < last)
1271 bad++;
1272 else
1273 last = startPositive_[i];
1274 if(startNegative_[i] < last)
1275 bad++;
1276 else
1277 last = startNegative_[i];
1278 }
1279 if(startPositive_[number] < last)
1280 bad++;
1281 CoinAssertHint(!bad, "starts are not monotonic");
1282 for (CoinBigIndex cbi = 0; cbi < numberElements; cbi++) {
1283 maxIndex = CoinMax(indices_[cbi], maxIndex);
1284 minIndex = CoinMin(indices_[cbi], minIndex);
1285 }
1286 CoinAssert(maxIndex < (columnOrdered_ ? numberRows_ : numberColumns_));
1287 CoinAssert(minIndex >= 0);
1288 if (detail) {
1289 if (minIndex > 0 || maxIndex + 1 < (columnOrdered_ ? numberRows_ : numberColumns_))
1290 printf("Not full range of indices - %d to %d\n", minIndex, maxIndex);
1291 }
1292}
1293/* Given positive integer weights for each row fills in sum of weights
1294 for each column (and slack).
1295 Returns weights vector
1296*/
1297CoinBigIndex *
1298ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
1299{
1300 int numberRows = model->numberRows();
1301 int numberColumns = model->numberColumns();
1302 int number = numberRows + numberColumns;
1303 CoinBigIndex * weights = new CoinBigIndex[number];
1304 int i;
1305 for (i = 0; i < numberColumns; i++) {
1306 CoinBigIndex j;
1307 CoinBigIndex count = 0;
1308 for (j = startPositive_[i]; j < startPositive_[i+1]; j++) {
1309 int iRow = indices_[j];
1310 count += inputWeights[iRow];
1311 }
1312 weights[i] = count;
1313 }
1314 for (i = 0; i < numberRows; i++) {
1315 weights[i+numberColumns] = inputWeights[i];
1316 }
1317 return weights;
1318}
1319// Append Columns
1320void
1321ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1322{
1323 int iColumn;
1324 CoinBigIndex size = 0;
1325 int numberBad = 0;
1326 for (iColumn = 0; iColumn < number; iColumn++) {
1327 int n = columns[iColumn]->getNumElements();
1328 const double * element = columns[iColumn]->getElements();
1329 size += n;
1330 int i;
1331 for (i = 0; i < n; i++) {
1332 if (fabs(element[i]) != 1.0)
1333 numberBad++;
1334 }
1335 }
1336 if (numberBad)
1337 throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
1338 // Get rid of temporary arrays
1339 delete [] lengths_;
1340 lengths_ = NULL;
1341 delete matrix_;
1342 matrix_ = NULL;
1343 int numberNow = startPositive_[numberColumns_];
1344 CoinBigIndex * temp;
1345 temp = new CoinBigIndex [numberColumns_+1+number];
1346 CoinMemcpyN(startPositive_, (numberColumns_ + 1), temp);
1347 delete [] startPositive_;
1348 startPositive_ = temp;
1349 temp = new CoinBigIndex [numberColumns_+number];
1350 CoinMemcpyN(startNegative_, numberColumns_, temp);
1351 delete [] startNegative_;
1352 startNegative_ = temp;
1353 int * temp2 = new int [numberNow+size];
1354 CoinMemcpyN(indices_, numberNow, temp2);
1355 delete [] indices_;
1356 indices_ = temp2;
1357 // now add
1358 size = numberNow;
1359 for (iColumn = 0; iColumn < number; iColumn++) {
1360 int n = columns[iColumn]->getNumElements();
1361 const int * row = columns[iColumn]->getIndices();
1362 const double * element = columns[iColumn]->getElements();
1363 int i;
1364 for (i = 0; i < n; i++) {
1365 if (element[i] == 1.0)
1366 indices_[size++] = row[i];
1367 }
1368 startNegative_[iColumn+numberColumns_] = size;
1369 for (i = 0; i < n; i++) {
1370 if (element[i] == -1.0)
1371 indices_[size++] = row[i];
1372 }
1373 startPositive_[iColumn+numberColumns_+1] = size;
1374 }
1375
1376 numberColumns_ += number;
1377}
1378// Append Rows
1379void
1380ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1381{
1382 // Allocate arrays to use for counting
1383 int * countPositive = new int [numberColumns_+1];
1384 memset(countPositive, 0, numberColumns_ * sizeof(int));
1385 int * countNegative = new int [numberColumns_];
1386 memset(countNegative, 0, numberColumns_ * sizeof(int));
1387 int iRow;
1388 CoinBigIndex size = 0;
1389 int numberBad = 0;
1390 for (iRow = 0; iRow < number; iRow++) {
1391 int n = rows[iRow]->getNumElements();
1392 const int * column = rows[iRow]->getIndices();
1393 const double * element = rows[iRow]->getElements();
1394 size += n;
1395 int i;
1396 for (i = 0; i < n; i++) {
1397 int iColumn = column[i];
1398 if (element[i] == 1.0)
1399 countPositive[iColumn]++;
1400 else if (element[i] == -1.0)
1401 countNegative[iColumn]++;
1402 else
1403 numberBad++;
1404 }
1405 }
1406 if (numberBad)
1407 throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
1408 // Get rid of temporary arrays
1409 delete [] lengths_;
1410 lengths_ = NULL;
1411 delete matrix_;
1412 matrix_ = NULL;
1413 int numberNow = startPositive_[numberColumns_];
1414 int * newIndices = new int [numberNow+size];
1415 // Update starts and turn counts into positions
1416 // also move current indices
1417 int iColumn;
1418 CoinBigIndex numberAdded = 0;
1419 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1420 int n, move;
1421 CoinBigIndex now;
1422 now = startPositive_[iColumn];
1423 move = startNegative_[iColumn] - now;
1424 n = countPositive[iColumn];
1425 startPositive_[iColumn] += numberAdded;
1426 CoinMemcpyN(newIndices + startPositive_[iColumn], move, indices_ + now);
1427 countPositive[iColumn] = startNegative_[iColumn] + numberAdded;
1428 numberAdded += n;
1429 now = startNegative_[iColumn];
1430 move = startPositive_[iColumn+1] - now;
1431 n = countNegative[iColumn];
1432 startNegative_[iColumn] += numberAdded;
1433 CoinMemcpyN(newIndices + startNegative_[iColumn], move, indices_ + now);
1434 countNegative[iColumn] = startPositive_[iColumn+1] + numberAdded;
1435 numberAdded += n;
1436 }
1437 delete [] indices_;
1438 indices_ = newIndices;
1439 startPositive_[numberColumns_] += numberAdded;
1440 // Now put in
1441 for (iRow = 0; iRow < number; iRow++) {
1442 int newRow = numberRows_ + iRow;
1443 int n = rows[iRow]->getNumElements();
1444 const int * column = rows[iRow]->getIndices();
1445 const double * element = rows[iRow]->getElements();
1446 int i;
1447 for (i = 0; i < n; i++) {
1448 int iColumn = column[i];
1449 int put;
1450 if (element[i] == 1.0) {
1451 put = countPositive[iColumn];
1452 countPositive[iColumn] = put + 1;
1453 } else {
1454 put = countNegative[iColumn];
1455 countNegative[iColumn] = put + 1;
1456 }
1457 indices_[put] = newRow;
1458 }
1459 }
1460 delete [] countPositive;
1461 delete [] countNegative;
1462 numberRows_ += number;
1463}
1464/* Returns largest and smallest elements of both signs.
1465 Largest refers to largest absolute value.
1466*/
1467void
1468ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1469 double & smallestPositive, double & largestPositive)
1470{
1471 int iColumn;
1472 bool plusOne = false;
1473 bool minusOne = false;
1474 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1475 if (startNegative_[iColumn] > startPositive_[iColumn])
1476 plusOne = true;
1477 if (startPositive_[iColumn+1] > startNegative_[iColumn])
1478 minusOne = true;
1479 }
1480 if (minusOne) {
1481 smallestNegative = -1.0;
1482 largestNegative = -1.0;
1483 } else {
1484 smallestNegative = 0.0;
1485 largestNegative = 0.0;
1486 }
1487 if (plusOne) {
1488 smallestPositive = 1.0;
1489 largestPositive = 1.0;
1490 } else {
1491 smallestPositive = 0.0;
1492 largestPositive = 0.0;
1493 }
1494}
1495// Says whether it can do partial pricing
1496bool
1497ClpPlusMinusOneMatrix::canDoPartialPricing() const
1498{
1499 return true;
1500}
1501// Partial pricing
1502void
1503ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1504 int & bestSequence, int & numberWanted)
1505{
1506 numberWanted = currentWanted_;
1507 int start = static_cast<int> (startFraction * numberColumns_);
1508 int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_);
1509 CoinBigIndex j;
1510 double tolerance = model->currentDualTolerance();
1511 double * reducedCost = model->djRegion();
1512 const double * duals = model->dualRowSolution();
1513 const double * cost = model->costRegion();
1514 double bestDj;
1515 if (bestSequence >= 0)
1516 bestDj = fabs(reducedCost[bestSequence]);
1517 else
1518 bestDj = tolerance;
1519 int sequenceOut = model->sequenceOut();
1520 int saveSequence = bestSequence;
1521 int iSequence;
1522 for (iSequence = start; iSequence < end; iSequence++) {
1523 if (iSequence != sequenceOut) {
1524 double value;
1525 ClpSimplex::Status status = model->getStatus(iSequence);
1526
1527 switch(status) {
1528
1529 case ClpSimplex::basic:
1530 case ClpSimplex::isFixed:
1531 break;
1532 case ClpSimplex::isFree:
1533 case ClpSimplex::superBasic:
1534 value = cost[iSequence];
1535 j = startPositive_[iSequence];
1536 for (; j < startNegative_[iSequence]; j++) {
1537 int iRow = indices_[j];
1538 value -= duals[iRow];
1539 }
1540 for (; j < startPositive_[iSequence+1]; j++) {
1541 int iRow = indices_[j];
1542 value += duals[iRow];
1543 }
1544 value = fabs(value);
1545 if (value > FREE_ACCEPT * tolerance) {
1546 numberWanted--;
1547 // we are going to bias towards free (but only if reasonable)
1548 value *= FREE_BIAS;
1549 if (value > bestDj) {
1550 // check flagged variable and correct dj
1551 if (!model->flagged(iSequence)) {
1552 bestDj = value;
1553 bestSequence = iSequence;
1554 } else {
1555 // just to make sure we don't exit before got something
1556 numberWanted++;
1557 }
1558 }
1559 }
1560 break;
1561 case ClpSimplex::atUpperBound:
1562 value = cost[iSequence];
1563 j = startPositive_[iSequence];
1564 for (; j < startNegative_[iSequence]; j++) {
1565 int iRow = indices_[j];
1566 value -= duals[iRow];
1567 }
1568 for (; j < startPositive_[iSequence+1]; j++) {
1569 int iRow = indices_[j];
1570 value += duals[iRow];
1571 }
1572 if (value > tolerance) {
1573 numberWanted--;
1574 if (value > bestDj) {
1575 // check flagged variable and correct dj
1576 if (!model->flagged(iSequence)) {
1577 bestDj = value;
1578 bestSequence = iSequence;
1579 } else {
1580 // just to make sure we don't exit before got something
1581 numberWanted++;
1582 }
1583 }
1584 }
1585 break;
1586 case ClpSimplex::atLowerBound:
1587 value = cost[iSequence];
1588 j = startPositive_[iSequence];
1589 for (; j < startNegative_[iSequence]; j++) {
1590 int iRow = indices_[j];
1591 value -= duals[iRow];
1592 }
1593 for (; j < startPositive_[iSequence+1]; j++) {
1594 int iRow = indices_[j];
1595 value += duals[iRow];
1596 }
1597 value = -value;
1598 if (value > tolerance) {
1599 numberWanted--;
1600 if (value > bestDj) {
1601 // check flagged variable and correct dj
1602 if (!model->flagged(iSequence)) {
1603 bestDj = value;
1604 bestSequence = iSequence;
1605 } else {
1606 // just to make sure we don't exit before got something
1607 numberWanted++;
1608 }
1609 }
1610 }
1611 break;
1612 }
1613 }
1614 if (!numberWanted)
1615 break;
1616 }
1617 if (bestSequence != saveSequence) {
1618 // recompute dj
1619 double value = cost[bestSequence];
1620 j = startPositive_[bestSequence];
1621 for (; j < startNegative_[bestSequence]; j++) {
1622 int iRow = indices_[j];
1623 value -= duals[iRow];
1624 }
1625 for (; j < startPositive_[bestSequence+1]; j++) {
1626 int iRow = indices_[j];
1627 value += duals[iRow];
1628 }
1629 reducedCost[bestSequence] = value;
1630 savedBestSequence_ = bestSequence;
1631 savedBestDj_ = reducedCost[savedBestSequence_];
1632 }
1633 currentWanted_ = numberWanted;
1634}
1635// Allow any parts of a created CoinMatrix to be deleted
1636void
1637ClpPlusMinusOneMatrix::releasePackedMatrix() const
1638{
1639 delete matrix_;
1640 delete [] lengths_;
1641 matrix_ = NULL;
1642 lengths_ = NULL;
1643}
1644/* Returns true if can combine transposeTimes and subsetTransposeTimes
1645 and if it would be faster */
1646bool
1647ClpPlusMinusOneMatrix::canCombine(const ClpSimplex * model,
1648 const CoinIndexedVector * pi) const
1649{
1650 int numberInRowArray = pi->getNumElements();
1651 int numberRows = model->numberRows();
1652 bool packed = pi->packedMode();
1653 // factor should be smaller if doing both with two pi vectors
1654 double factor = 0.27;
1655 // We may not want to do by row if there may be cache problems
1656 // It would be nice to find L2 cache size - for moment 512K
1657 // Be slightly optimistic
1658 if (numberColumns_ * sizeof(double) > 1000000) {
1659 if (numberRows * 10 < numberColumns_)
1660 factor *= 0.333333333;
1661 else if (numberRows * 4 < numberColumns_)
1662 factor *= 0.5;
1663 else if (numberRows * 2 < numberColumns_)
1664 factor *= 0.66666666667;
1665 //if (model->numberIterations()%50==0)
1666 //printf("%d nonzero\n",numberInRowArray);
1667 }
1668 // if not packed then bias a bit more towards by column
1669 if (!packed)
1670 factor *= 0.9;
1671 return (numberInRowArray > factor * numberRows || !model->rowCopy());
1672}
1673// These have to match ClpPrimalColumnSteepest version
1674#define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)
1675// Updates two arrays for steepest
1676void
1677ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex * model,
1678 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
1679 const CoinIndexedVector * pi2,
1680 CoinIndexedVector * spare,
1681 double referenceIn, double devex,
1682 // Array for exact devex to say what is in reference framework
1683 unsigned int * reference,
1684 double * weights, double scaleFactor)
1685{
1686 // put row of tableau in dj1
1687 double * pi = pi1->denseVector();
1688 int numberNonZero = 0;
1689 int * index = dj1->getIndices();
1690 double * array = dj1->denseVector();
1691 int numberInRowArray = pi1->getNumElements();
1692 double zeroTolerance = model->zeroTolerance();
1693 bool packed = pi1->packedMode();
1694 // do by column
1695 int iColumn;
1696 assert (!spare->getNumElements());
1697 double * piWeight = pi2->denseVector();
1698 assert (!pi2->packedMode());
1699 bool killDjs = (scaleFactor == 0.0);
1700 if (!scaleFactor)
1701 scaleFactor = 1.0;
1702 // Note scale factor was -1.0
1703 if (packed) {
1704 // need to expand pi into y
1705 assert(spare->capacity() >= model->numberRows());
1706 double * piOld = pi;
1707 pi = spare->denseVector();
1708 const int * whichRow = pi1->getIndices();
1709 int i;
1710 // modify pi so can collapse to one loop
1711 for (i = 0; i < numberInRowArray; i++) {
1712 int iRow = whichRow[i];
1713 pi[iRow] = piOld[i];
1714 }
1715 CoinBigIndex j;
1716 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1717 ClpSimplex::Status status = model->getStatus(iColumn);
1718 if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
1719 double value = 0.0;
1720 for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1721 int iRow = indices_[j];
1722 value -= pi[iRow];
1723 }
1724 for (; j < startPositive_[iColumn+1]; j++) {
1725 int iRow = indices_[j];
1726 value += pi[iRow];
1727 }
1728 if (fabs(value) > zeroTolerance) {
1729 // and do other array
1730 double modification = 0.0;
1731 for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1732 int iRow = indices_[j];
1733 modification += piWeight[iRow];
1734 }
1735 for (; j < startPositive_[iColumn+1]; j++) {
1736 int iRow = indices_[j];
1737 modification -= piWeight[iRow];
1738 }
1739 double thisWeight = weights[iColumn];
1740 double pivot = value * scaleFactor;
1741 double pivotSquared = pivot * pivot;
1742 thisWeight += pivotSquared * devex + pivot * modification;
1743 if (thisWeight < DEVEX_TRY_NORM) {
1744 if (referenceIn < 0.0) {
1745 // steepest
1746 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1747 } else {
1748 // exact
1749 thisWeight = referenceIn * pivotSquared;
1750 if (reference(iColumn))
1751 thisWeight += 1.0;
1752 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1753 }
1754 }
1755 weights[iColumn] = thisWeight;
1756 if (!killDjs) {
1757 array[numberNonZero] = value;
1758 index[numberNonZero++] = iColumn;
1759 }
1760 }
1761 }
1762 // zero out
1763 for (i = 0; i < numberInRowArray; i++) {
1764 int iRow = whichRow[i];
1765 pi[iRow] = 0.0;
1766 }
1767 } else {
1768 CoinBigIndex j;
1769 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1770 ClpSimplex::Status status = model->getStatus(iColumn);
1771 if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
1772 double value = 0.0;
1773 for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1774 int iRow = indices_[j];
1775 value -= pi[iRow];
1776 }
1777 for (; j < startPositive_[iColumn+1]; j++) {
1778 int iRow = indices_[j];
1779 value += pi[iRow];
1780 }
1781 if (fabs(value) > zeroTolerance) {
1782 // and do other array
1783 double modification = 0.0;
1784 for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1785 int iRow = indices_[j];
1786 modification += piWeight[iRow];
1787 }
1788 for (; j < startPositive_[iColumn+1]; j++) {
1789 int iRow = indices_[j];
1790 modification -= piWeight[iRow];
1791 }
1792 double thisWeight = weights[iColumn];
1793 double pivot = value * scaleFactor;
1794 double pivotSquared = pivot * pivot;
1795 thisWeight += pivotSquared * devex + pivot * modification;
1796 if (thisWeight < DEVEX_TRY_NORM) {
1797 if (referenceIn < 0.0) {
1798 // steepest
1799 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1800 } else {
1801 // exact
1802 thisWeight = referenceIn * pivotSquared;
1803 if (reference(iColumn))
1804 thisWeight += 1.0;
1805 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1806 }
1807 }
1808 weights[iColumn] = thisWeight;
1809 if (!killDjs) {
1810 array[iColumn] = value;
1811 index[numberNonZero++] = iColumn;
1812 }
1813 }
1814 }
1815 }
1816 dj1->setNumElements(numberNonZero);
1817 spare->setNumElements(0);
1818 if (packed)
1819 dj1->setPackedMode(true);
1820}
1821// Updates second array for steepest and does devex weights
1822void
1823ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex * ,
1824 CoinIndexedVector * dj1,
1825 const CoinIndexedVector * pi2, CoinIndexedVector *,
1826 double referenceIn, double devex,
1827 // Array for exact devex to say what is in reference framework
1828 unsigned int * reference,
1829 double * weights, double scaleFactor)
1830{
1831 int number = dj1->getNumElements();
1832 const int * index = dj1->getIndices();
1833 double * array = dj1->denseVector();
1834 assert( dj1->packedMode());
1835
1836 double * piWeight = pi2->denseVector();
1837 bool killDjs = (scaleFactor == 0.0);
1838 if (!scaleFactor)
1839 scaleFactor = 1.0;
1840 for (int k = 0; k < number; k++) {
1841 int iColumn = index[k];
1842 double pivot = array[k] * scaleFactor;
1843 if (killDjs)
1844 array[k] = 0.0;
1845 // and do other array
1846 double modification = 0.0;
1847 CoinBigIndex j;
1848 for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1849 int iRow = indices_[j];
1850 modification += piWeight[iRow];
1851 }
1852 for (; j < startPositive_[iColumn+1]; j++) {
1853 int iRow = indices_[j];
1854 modification -= piWeight[iRow];
1855 }
1856 double thisWeight = weights[iColumn];
1857 double pivotSquared = pivot * pivot;
1858 thisWeight += pivotSquared * devex + pivot * modification;
1859 if (thisWeight < DEVEX_TRY_NORM) {
1860 if (referenceIn < 0.0) {
1861 // steepest
1862 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1863 } else {
1864 // exact
1865 thisWeight = referenceIn * pivotSquared;
1866 if (reference(iColumn))
1867 thisWeight += 1.0;
1868 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1869 }
1870 }
1871 weights[iColumn] = thisWeight;
1872 }
1873}
1874/* Set the dimensions of the matrix. In effect, append new empty
1875 columns/rows to the matrix. A negative number for either dimension
1876 means that that dimension doesn't change. Otherwise the new dimensions
1877 MUST be at least as large as the current ones otherwise an exception
1878 is thrown. */
1879void
1880ClpPlusMinusOneMatrix::setDimensions(int newnumrows, int newnumcols)
1881{
1882 if (newnumrows < 0)
1883 newnumrows = numberRows_;
1884 if (newnumrows < numberRows_)
1885 throw CoinError("Bad new rownum (less than current)",
1886 "setDimensions", "CoinPackedMatrix");
1887
1888 if (newnumcols < 0)
1889 newnumcols = numberColumns_;
1890 if (newnumcols < numberColumns_)
1891 throw CoinError("Bad new colnum (less than current)",
1892 "setDimensions", "CoinPackedMatrix");
1893
1894 int number = 0;
1895 int length = 0;
1896 if (columnOrdered_) {
1897 length = numberColumns_;
1898 numberColumns_ = newnumcols;
1899 number = numberColumns_;
1900
1901 } else {
1902 length = numberRows_;
1903 numberRows_ = newnumrows;
1904 number = numberRows_;
1905 }
1906 if (number > length) {
1907 CoinBigIndex * temp;
1908 int i;
1909 CoinBigIndex end = startPositive_[length];
1910 temp = new CoinBigIndex [number+1];
1911 CoinMemcpyN(startPositive_, (length + 1), temp);
1912 delete [] startPositive_;
1913 for (i = length + 1; i < number + 1; i++)
1914 temp[i] = end;
1915 startPositive_ = temp;
1916 temp = new CoinBigIndex [number];
1917 CoinMemcpyN(startNegative_, length, temp);
1918 delete [] startNegative_;
1919 for (i = length; i < number; i++)
1920 temp[i] = end;
1921 startNegative_ = temp;
1922 }
1923}
1924#ifndef SLIM_CLP
1925/* Append a set of rows/columns to the end of the matrix. Returns number of errors
1926 i.e. if any of the new rows/columns contain an index that's larger than the
1927 number of columns-1/rows-1 (if numberOther>0) or duplicates
1928 If 0 then rows, 1 if columns */
1929int
1930ClpPlusMinusOneMatrix::appendMatrix(int number, int type,
1931 const CoinBigIndex * starts, const int * index,
1932 const double * element, int /*numberOther*/)
1933{
1934 int numberErrors = 0;
1935 // make into CoinPackedVector
1936 CoinPackedVectorBase ** vectors =
1937 new CoinPackedVectorBase * [number];
1938 int iVector;
1939 for (iVector = 0; iVector < number; iVector++) {
1940 int iStart = starts[iVector];
1941 vectors[iVector] =
1942 new CoinPackedVector(starts[iVector+1] - iStart,
1943 index + iStart, element + iStart);
1944 }
1945 if (type == 0) {
1946 // rows
1947 appendRows(number, vectors);
1948 } else {
1949 // columns
1950 appendCols(number, vectors);
1951 }
1952 for (iVector = 0; iVector < number; iVector++)
1953 delete vectors[iVector];
1954 delete [] vectors;
1955 return numberErrors;
1956}
1957#endif
1958