1 /* $Id: ClpPackedMatrix.cpp 1864 2012-06-28 10:27:20Z forrest $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7
8 #include <cstdio>
9
10 #include "CoinPragma.hpp"
11 #include "CoinIndexedVector.hpp"
12 #include "CoinHelperFunctions.hpp"
13 //#define THREAD
14
15 #include "ClpSimplex.hpp"
16 #include "ClpSimplexDual.hpp"
17 #include "ClpFactorization.hpp"
18 #ifndef SLIM_CLP
19 #include "ClpQuadraticObjective.hpp"
20 #endif
21 // at end to get min/max!
22 #include "ClpPackedMatrix.hpp"
23 #include "ClpMessage.hpp"
24 #ifdef INTEL_MKL
25 #include "mkl_spblas.h"
26 #endif
27
28 //=============================================================================
29 #ifdef COIN_PREFETCH
30 #if 1
31 #define coin_prefetch(mem) \
32 __asm__ __volatile__ ("prefetchnta %0" : : "m" (*(reinterpret_cast<char *>(mem))))
33 #define coin_prefetch_const(mem) \
34 __asm__ __volatile__ ("prefetchnta %0" : : "m" (*(reinterpret_cast<const char *>(mem))))
35 #else
36 #define coin_prefetch(mem) \
37 __asm__ __volatile__ ("prefetch %0" : : "m" (*(reinterpret_cast<char *>(mem))))
38 #define coin_prefetch_const(mem) \
39 __asm__ __volatile__ ("prefetch %0" : : "m" (*(reinterpret_cast<const char *>(mem))))
40 #endif
41 #else
42 // dummy
43 #define coin_prefetch(mem)
44 #define coin_prefetch_const(mem)
45 #endif
46
47 //#############################################################################
48 // Constructors / Destructor / Assignment
49 //#############################################################################
50
51 //-------------------------------------------------------------------
52 // Default Constructor
53 //-------------------------------------------------------------------
ClpPackedMatrix()54 ClpPackedMatrix::ClpPackedMatrix ()
55 : ClpMatrixBase(),
56 matrix_(NULL),
57 numberActiveColumns_(0),
58 flags_(2),
59 rowCopy_(NULL),
60 columnCopy_(NULL)
61 {
62 setType(1);
63 }
64
65 //-------------------------------------------------------------------
66 // Copy constructor
67 //-------------------------------------------------------------------
ClpPackedMatrix(const ClpPackedMatrix & rhs)68 ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs)
69 : ClpMatrixBase(rhs)
70 {
71 #ifndef COIN_SPARSE_MATRIX
72 matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -1, -1);
73 #else
74 matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
75 #endif
76 numberActiveColumns_ = rhs.numberActiveColumns_;
77 flags_ = rhs.flags_ & (~2);
78 int numberRows = matrix_->getNumRows();
79 if (rhs.rhsOffset_ && numberRows) {
80 rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
81 } else {
82 rhsOffset_ = NULL;
83 }
84 if (rhs.rowCopy_) {
85 assert ((flags_ & 4) != 0);
86 rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
87 } else {
88 rowCopy_ = NULL;
89 }
90 if (rhs.columnCopy_) {
91 assert ((flags_&(8 + 16)) == 8 + 16);
92 columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
93 } else {
94 columnCopy_ = NULL;
95 }
96 }
97 //-------------------------------------------------------------------
98 // assign matrix (for space reasons)
99 //-------------------------------------------------------------------
ClpPackedMatrix(CoinPackedMatrix * rhs)100 ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs)
101 : ClpMatrixBase()
102 {
103 matrix_ = rhs;
104 flags_ = matrix_->hasGaps() ? 2 : 0;
105 numberActiveColumns_ = matrix_->getNumCols();
106 rowCopy_ = NULL;
107 columnCopy_ = NULL;
108 setType(1);
109
110 }
111
ClpPackedMatrix(const CoinPackedMatrix & rhs)112 ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs)
113 : ClpMatrixBase()
114 {
115 #ifndef COIN_SPARSE_MATRIX
116 matrix_ = new CoinPackedMatrix(rhs, -1, -1);
117 #else
118 matrix_ = new CoinPackedMatrix(rhs, -0, -0);
119 #endif
120 numberActiveColumns_ = matrix_->getNumCols();
121 rowCopy_ = NULL;
122 flags_ = 0;
123 columnCopy_ = NULL;
124 setType(1);
125
126 }
127
128 //-------------------------------------------------------------------
129 // Destructor
130 //-------------------------------------------------------------------
~ClpPackedMatrix()131 ClpPackedMatrix::~ClpPackedMatrix ()
132 {
133 delete matrix_;
134 delete rowCopy_;
135 delete columnCopy_;
136 }
137
138 //----------------------------------------------------------------
139 // Assignment operator
140 //-------------------------------------------------------------------
141 ClpPackedMatrix &
operator =(const ClpPackedMatrix & rhs)142 ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
143 {
144 if (this != &rhs) {
145 ClpMatrixBase::operator=(rhs);
146 delete matrix_;
147 #ifndef COIN_SPARSE_MATRIX
148 matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
149 #else
150 matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
151 #endif
152 numberActiveColumns_ = rhs.numberActiveColumns_;
153 flags_ = rhs.flags_;
154 delete rowCopy_;
155 delete columnCopy_;
156 if (rhs.rowCopy_) {
157 assert ((flags_ & 4) != 0);
158 rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
159 } else {
160 rowCopy_ = NULL;
161 }
162 if (rhs.columnCopy_) {
163 assert ((flags_&(8 + 16)) == 8 + 16);
164 columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
165 } else {
166 columnCopy_ = NULL;
167 }
168 }
169 return *this;
170 }
171 //-------------------------------------------------------------------
172 // Clone
173 //-------------------------------------------------------------------
clone() const174 ClpMatrixBase * ClpPackedMatrix::clone() const
175 {
176 return new ClpPackedMatrix(*this);
177 }
178 // Copy contents - resizing if necessary - otherwise re-use memory
179 void
copy(const ClpPackedMatrix * rhs)180 ClpPackedMatrix::copy(const ClpPackedMatrix * rhs)
181 {
182 //*this = *rhs;
183 assert(numberActiveColumns_ == rhs->numberActiveColumns_);
184 assert (matrix_->isColOrdered() == rhs->matrix_->isColOrdered());
185 matrix_->copyReuseArrays(*rhs->matrix_);
186 }
187 /* Subset clone (without gaps). Duplicates are allowed
188 and order is as given */
189 ClpMatrixBase *
subsetClone(int numberRows,const int * whichRows,int numberColumns,const int * whichColumns) const190 ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
191 int numberColumns,
192 const int * whichColumns) const
193 {
194 return new ClpPackedMatrix(*this, numberRows, whichRows,
195 numberColumns, whichColumns);
196 }
197 /* Subset constructor (without gaps). Duplicates are allowed
198 and order is as given */
ClpPackedMatrix(const ClpPackedMatrix & rhs,int numberRows,const int * whichRows,int numberColumns,const int * whichColumns)199 ClpPackedMatrix::ClpPackedMatrix (
200 const ClpPackedMatrix & rhs,
201 int numberRows, const int * whichRows,
202 int numberColumns, const int * whichColumns)
203 : ClpMatrixBase(rhs)
204 {
205 matrix_ = new CoinPackedMatrix(*(rhs.matrix_), numberRows, whichRows,
206 numberColumns, whichColumns);
207 numberActiveColumns_ = matrix_->getNumCols();
208 rowCopy_ = NULL;
209 flags_ = rhs.flags_ & (~2); // no gaps
210 columnCopy_ = NULL;
211 }
ClpPackedMatrix(const CoinPackedMatrix & rhs,int numberRows,const int * whichRows,int numberColumns,const int * whichColumns)212 ClpPackedMatrix::ClpPackedMatrix (
213 const CoinPackedMatrix & rhs,
214 int numberRows, const int * whichRows,
215 int numberColumns, const int * whichColumns)
216 : ClpMatrixBase()
217 {
218 matrix_ = new CoinPackedMatrix(rhs, numberRows, whichRows,
219 numberColumns, whichColumns);
220 numberActiveColumns_ = matrix_->getNumCols();
221 rowCopy_ = NULL;
222 flags_ = 2;
223 columnCopy_ = NULL;
224 setType(1);
225 }
226
227 /* Returns a new matrix in reverse order without gaps */
228 ClpMatrixBase *
reverseOrderedCopy() const229 ClpPackedMatrix::reverseOrderedCopy() const
230 {
231 ClpPackedMatrix * copy = new ClpPackedMatrix();
232 copy->matrix_ = new CoinPackedMatrix();
233 copy->matrix_->setExtraGap(0.0);
234 copy->matrix_->setExtraMajor(0.0);
235 copy->matrix_->reverseOrderedCopyOf(*matrix_);
236 //copy->matrix_->removeGaps();
237 copy->numberActiveColumns_ = copy->matrix_->getNumCols();
238 copy->flags_ = flags_ & (~2); // no gaps
239 return copy;
240 }
241 //unscaled versions
242 void
times(double scalar,const double * x,double * y) const243 ClpPackedMatrix::times(double scalar,
244 const double * x, double * y) const
245 {
246 int iRow, iColumn;
247 // get matrix data pointers
248 const int * row = matrix_->getIndices();
249 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
250 const double * elementByColumn = matrix_->getElements();
251 //memset(y,0,matrix_->getNumRows()*sizeof(double));
252 assert (((flags_ & 2) != 0) == matrix_->hasGaps());
253 if (!(flags_ & 2)) {
254 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
255 CoinBigIndex j;
256 double value = x[iColumn];
257 if (value) {
258 CoinBigIndex start = columnStart[iColumn];
259 CoinBigIndex end = columnStart[iColumn+1];
260 value *= scalar;
261 for (j = start; j < end; j++) {
262 iRow = row[j];
263 y[iRow] += value * elementByColumn[j];
264 }
265 }
266 }
267 } else {
268 const int * columnLength = matrix_->getVectorLengths();
269 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
270 CoinBigIndex j;
271 double value = x[iColumn];
272 if (value) {
273 CoinBigIndex start = columnStart[iColumn];
274 CoinBigIndex end = start + columnLength[iColumn];
275 value *= scalar;
276 for (j = start; j < end; j++) {
277 iRow = row[j];
278 y[iRow] += value * elementByColumn[j];
279 }
280 }
281 }
282 }
283 }
284 void
transposeTimes(double scalar,const double * x,double * y) const285 ClpPackedMatrix::transposeTimes(double scalar,
286 const double * x, double * y) const
287 {
288 int iColumn;
289 // get matrix data pointers
290 const int * row = matrix_->getIndices();
291 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
292 const double * elementByColumn = matrix_->getElements();
293 if (!(flags_ & 2)) {
294 if (scalar == -1.0) {
295 CoinBigIndex start = columnStart[0];
296 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
297 CoinBigIndex j;
298 CoinBigIndex next = columnStart[iColumn+1];
299 double value = y[iColumn];
300 for (j = start; j < next; j++) {
301 int jRow = row[j];
302 value -= x[jRow] * elementByColumn[j];
303 }
304 start = next;
305 y[iColumn] = value;
306 }
307 } else {
308 CoinBigIndex start = columnStart[0];
309 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
310 CoinBigIndex j;
311 CoinBigIndex next = columnStart[iColumn+1];
312 double value = 0.0;
313 for (j = start; j < next; j++) {
314 int jRow = row[j];
315 value += x[jRow] * elementByColumn[j];
316 }
317 start = next;
318 y[iColumn] += value * scalar;
319 }
320 }
321 } else {
322 const int * columnLength = matrix_->getVectorLengths();
323 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
324 CoinBigIndex j;
325 double value = 0.0;
326 CoinBigIndex start = columnStart[iColumn];
327 CoinBigIndex end = start + columnLength[iColumn];
328 for (j = start; j < end; j++) {
329 int jRow = row[j];
330 value += x[jRow] * elementByColumn[j];
331 }
332 y[iColumn] += value * scalar;
333 }
334 }
335 }
336 void
times(double scalar,const double * COIN_RESTRICT x,double * COIN_RESTRICT y,const double * COIN_RESTRICT rowScale,const double * COIN_RESTRICT columnScale) const337 ClpPackedMatrix::times(double scalar,
338 const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
339 const double * COIN_RESTRICT rowScale,
340 const double * COIN_RESTRICT columnScale) const
341 {
342 if (rowScale) {
343 int iRow, iColumn;
344 // get matrix data pointers
345 const int * COIN_RESTRICT row = matrix_->getIndices();
346 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
347 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
348 if (!(flags_ & 2)) {
349 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
350 double value = x[iColumn];
351 if (value) {
352 // scaled
353 value *= scalar * columnScale[iColumn];
354 CoinBigIndex start = columnStart[iColumn];
355 CoinBigIndex end = columnStart[iColumn+1];
356 CoinBigIndex j;
357 for (j = start; j < end; j++) {
358 iRow = row[j];
359 y[iRow] += value * elementByColumn[j] * rowScale[iRow];
360 }
361 }
362 }
363 } else {
364 const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
365 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
366 double value = x[iColumn];
367 if (value) {
368 // scaled
369 value *= scalar * columnScale[iColumn];
370 CoinBigIndex start = columnStart[iColumn];
371 CoinBigIndex end = start + columnLength[iColumn];
372 CoinBigIndex j;
373 for (j = start; j < end; j++) {
374 iRow = row[j];
375 y[iRow] += value * elementByColumn[j] * rowScale[iRow];
376 }
377 }
378 }
379 }
380 } else {
381 times(scalar, x, y);
382 }
383 }
384 void
transposeTimes(double scalar,const double * COIN_RESTRICT x,double * COIN_RESTRICT y,const double * COIN_RESTRICT rowScale,const double * COIN_RESTRICT columnScale,double * COIN_RESTRICT spare) const385 ClpPackedMatrix::transposeTimes( double scalar,
386 const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
387 const double * COIN_RESTRICT rowScale,
388 const double * COIN_RESTRICT columnScale,
389 double * COIN_RESTRICT spare) const
390 {
391 if (rowScale) {
392 int iColumn;
393 // get matrix data pointers
394 const int * COIN_RESTRICT row = matrix_->getIndices();
395 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
396 const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
397 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
398 if (!spare) {
399 if (!(flags_ & 2)) {
400 CoinBigIndex start = columnStart[0];
401 if (scalar == -1.0) {
402 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
403 CoinBigIndex j;
404 CoinBigIndex next = columnStart[iColumn+1];
405 double value = 0.0;
406 // scaled
407 for (j = start; j < next; j++) {
408 int jRow = row[j];
409 value += x[jRow] * elementByColumn[j] * rowScale[jRow];
410 }
411 start = next;
412 y[iColumn] -= value * columnScale[iColumn];
413 }
414 } else {
415 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
416 CoinBigIndex j;
417 CoinBigIndex next = columnStart[iColumn+1];
418 double value = 0.0;
419 // scaled
420 for (j = start; j < next; j++) {
421 int jRow = row[j];
422 value += x[jRow] * elementByColumn[j] * rowScale[jRow];
423 }
424 start = next;
425 y[iColumn] += value * scalar * columnScale[iColumn];
426 }
427 }
428 } else {
429 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
430 CoinBigIndex j;
431 double value = 0.0;
432 // scaled
433 for (j = columnStart[iColumn];
434 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
435 int jRow = row[j];
436 value += x[jRow] * elementByColumn[j] * rowScale[jRow];
437 }
438 y[iColumn] += value * scalar * columnScale[iColumn];
439 }
440 }
441 } else {
442 // can use spare region
443 int iRow;
444 int numberRows = matrix_->getNumRows();
445 for (iRow = 0; iRow < numberRows; iRow++) {
446 double value = x[iRow];
447 if (value)
448 spare[iRow] = value * rowScale[iRow];
449 else
450 spare[iRow] = 0.0;
451 }
452 if (!(flags_ & 2)) {
453 CoinBigIndex start = columnStart[0];
454 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
455 CoinBigIndex j;
456 CoinBigIndex next = columnStart[iColumn+1];
457 double value = 0.0;
458 // scaled
459 for (j = start; j < next; j++) {
460 int jRow = row[j];
461 value += spare[jRow] * elementByColumn[j];
462 }
463 start = next;
464 y[iColumn] += value * scalar * columnScale[iColumn];
465 }
466 } else {
467 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
468 CoinBigIndex j;
469 double value = 0.0;
470 // scaled
471 for (j = columnStart[iColumn];
472 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
473 int jRow = row[j];
474 value += spare[jRow] * elementByColumn[j];
475 }
476 y[iColumn] += value * scalar * columnScale[iColumn];
477 }
478 }
479 // no need to zero out
480 //for (iRow=0;iRow<numberRows;iRow++)
481 //spare[iRow] = 0.0;
482 }
483 } else {
484 transposeTimes(scalar, x, y);
485 }
486 }
487 void
transposeTimesSubset(int number,const int * which,const double * COIN_RESTRICT x,double * COIN_RESTRICT y,const double * COIN_RESTRICT rowScale,const double * COIN_RESTRICT columnScale,double * COIN_RESTRICT spare) const488 ClpPackedMatrix::transposeTimesSubset( int number,
489 const int * which,
490 const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
491 const double * COIN_RESTRICT rowScale,
492 const double * COIN_RESTRICT columnScale,
493 double * COIN_RESTRICT spare) const
494 {
495 // get matrix data pointers
496 const int * COIN_RESTRICT row = matrix_->getIndices();
497 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
498 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
499 if (!spare || !rowScale) {
500 if (rowScale) {
501 for (int jColumn = 0; jColumn < number; jColumn++) {
502 int iColumn = which[jColumn];
503 CoinBigIndex j;
504 CoinBigIndex start = columnStart[iColumn];
505 CoinBigIndex next = columnStart[iColumn+1];
506 double value = 0.0;
507 for (j = start; j < next; j++) {
508 int jRow = row[j];
509 value += x[jRow] * elementByColumn[j] * rowScale[jRow];
510 }
511 y[iColumn] -= value * columnScale[iColumn];
512 }
513 } else {
514 for (int jColumn = 0; jColumn < number; jColumn++) {
515 int iColumn = which[jColumn];
516 CoinBigIndex j;
517 CoinBigIndex start = columnStart[iColumn];
518 CoinBigIndex next = columnStart[iColumn+1];
519 double value = 0.0;
520 for (j = start; j < next; j++) {
521 int jRow = row[j];
522 value += x[jRow] * elementByColumn[j];
523 }
524 y[iColumn] -= value;
525 }
526 }
527 } else {
528 // can use spare region
529 int iRow;
530 int numberRows = matrix_->getNumRows();
531 for (iRow = 0; iRow < numberRows; iRow++) {
532 double value = x[iRow];
533 if (value)
534 spare[iRow] = value * rowScale[iRow];
535 else
536 spare[iRow] = 0.0;
537 }
538 for (int jColumn = 0; jColumn < number; jColumn++) {
539 int iColumn = which[jColumn];
540 CoinBigIndex j;
541 CoinBigIndex start = columnStart[iColumn];
542 CoinBigIndex next = columnStart[iColumn+1];
543 double value = 0.0;
544 for (j = start; j < next; j++) {
545 int jRow = row[j];
546 value += spare[jRow] * elementByColumn[j];
547 }
548 y[iColumn] -= value * columnScale[iColumn];
549 }
550 }
551 }
552 /* Return <code>x * A + y</code> in <code>z</code>.
553 Squashes small elements and knows about ClpSimplex */
554 void
transposeTimes(const ClpSimplex * model,double scalar,const CoinIndexedVector * rowArray,CoinIndexedVector * y,CoinIndexedVector * columnArray) const555 ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
556 const CoinIndexedVector * rowArray,
557 CoinIndexedVector * y,
558 CoinIndexedVector * columnArray) const
559 {
560 columnArray->clear();
561 double * pi = rowArray->denseVector();
562 int numberNonZero = 0;
563 int * index = columnArray->getIndices();
564 double * array = columnArray->denseVector();
565 int numberInRowArray = rowArray->getNumElements();
566 // maybe I need one in OsiSimplex
567 double zeroTolerance = model->zeroTolerance();
568 #if 0 //def COIN_DEVELOP
569 if (zeroTolerance != 1.0e-13) {
570 printf("small element in matrix - zero tolerance %g\n", zeroTolerance);
571 }
572 #endif
573 int numberRows = model->numberRows();
574 ClpPackedMatrix* rowCopy =
575 static_cast< ClpPackedMatrix*>(model->rowCopy());
576 bool packed = rowArray->packedMode();
577 double factor = (numberRows < 100) ? 0.25 : 0.35;
578 factor = 0.5;
579 // We may not want to do by row if there may be cache problems
580 // It would be nice to find L2 cache size - for moment 512K
581 // Be slightly optimistic
582 if (numberActiveColumns_ * sizeof(double) > 1000000) {
583 if (numberRows * 10 < numberActiveColumns_)
584 factor *= 0.333333333;
585 else if (numberRows * 4 < numberActiveColumns_)
586 factor *= 0.5;
587 else if (numberRows * 2 < numberActiveColumns_)
588 factor *= 0.66666666667;
589 //if (model->numberIterations()%50==0)
590 //printf("%d nonzero\n",numberInRowArray);
591 }
592 // if not packed then bias a bit more towards by column
593 if (!packed)
594 factor *= 0.9;
595 assert (!y->getNumElements());
596 double multiplierX = 0.8;
597 double factor2 = factor * multiplierX;
598 if (packed && rowCopy_ && numberInRowArray > 2 && numberInRowArray > factor2 * numberRows &&
599 numberInRowArray < 0.9 * numberRows && 0) {
600 rowCopy_->transposeTimes(model, rowCopy->matrix_, rowArray, y, columnArray);
601 return;
602 }
603 if (numberInRowArray > factor * numberRows || !rowCopy) {
604 // do by column
605 // If no gaps - can do a bit faster
606 if (!(flags_ & 2) || columnCopy_) {
607 transposeTimesByColumn( model, scalar,
608 rowArray, y, columnArray);
609 return;
610 }
611 int iColumn;
612 // get matrix data pointers
613 const int * row = matrix_->getIndices();
614 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
615 const int * columnLength = matrix_->getVectorLengths();
616 const double * elementByColumn = matrix_->getElements();
617 const double * rowScale = model->rowScale();
618 #if 0
619 ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
620 if (rowScale && scaledMatrix) {
621 rowScale = NULL;
622 // get matrix data pointers
623 row = scaledMatrix->getIndices();
624 columnStart = scaledMatrix->getVectorStarts();
625 columnLength = scaledMatrix->getVectorLengths();
626 elementByColumn = scaledMatrix->getElements();
627 }
628 #endif
629 if (packed) {
630 // need to expand pi into y
631 assert(y->capacity() >= numberRows);
632 double * piOld = pi;
633 pi = y->denseVector();
634 const int * whichRow = rowArray->getIndices();
635 int i;
636 if (!rowScale) {
637 // modify pi so can collapse to one loop
638 if (scalar == -1.0) {
639 for (i = 0; i < numberInRowArray; i++) {
640 int iRow = whichRow[i];
641 pi[iRow] = -piOld[i];
642 }
643 } else {
644 for (i = 0; i < numberInRowArray; i++) {
645 int iRow = whichRow[i];
646 pi[iRow] = scalar * piOld[i];
647 }
648 }
649 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
650 double value = 0.0;
651 CoinBigIndex j;
652 for (j = columnStart[iColumn];
653 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
654 int iRow = row[j];
655 value += pi[iRow] * elementByColumn[j];
656 }
657 if (fabs(value) > zeroTolerance) {
658 array[numberNonZero] = value;
659 index[numberNonZero++] = iColumn;
660 }
661 }
662 } else {
663 #ifdef CLP_INVESTIGATE
664 if (model->clpScaledMatrix())
665 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
666 #endif
667 // scaled
668 // modify pi so can collapse to one loop
669 if (scalar == -1.0) {
670 for (i = 0; i < numberInRowArray; i++) {
671 int iRow = whichRow[i];
672 pi[iRow] = -piOld[i] * rowScale[iRow];
673 }
674 } else {
675 for (i = 0; i < numberInRowArray; i++) {
676 int iRow = whichRow[i];
677 pi[iRow] = scalar * piOld[i] * rowScale[iRow];
678 }
679 }
680 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
681 double value = 0.0;
682 CoinBigIndex j;
683 const double * columnScale = model->columnScale();
684 for (j = columnStart[iColumn];
685 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
686 int iRow = row[j];
687 value += pi[iRow] * elementByColumn[j];
688 }
689 value *= columnScale[iColumn];
690 if (fabs(value) > zeroTolerance) {
691 array[numberNonZero] = value;
692 index[numberNonZero++] = iColumn;
693 }
694 }
695 }
696 // zero out
697 for (i = 0; i < numberInRowArray; i++) {
698 int iRow = whichRow[i];
699 pi[iRow] = 0.0;
700 }
701 } else {
702 if (!rowScale) {
703 if (scalar == -1.0) {
704 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
705 double value = 0.0;
706 CoinBigIndex j;
707 for (j = columnStart[iColumn];
708 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
709 int iRow = row[j];
710 value += pi[iRow] * elementByColumn[j];
711 }
712 if (fabs(value) > zeroTolerance) {
713 index[numberNonZero++] = iColumn;
714 array[iColumn] = -value;
715 }
716 }
717 } else {
718 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
719 double value = 0.0;
720 CoinBigIndex j;
721 for (j = columnStart[iColumn];
722 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
723 int iRow = row[j];
724 value += pi[iRow] * elementByColumn[j];
725 }
726 value *= scalar;
727 if (fabs(value) > zeroTolerance) {
728 index[numberNonZero++] = iColumn;
729 array[iColumn] = value;
730 }
731 }
732 }
733 } else {
734 #ifdef CLP_INVESTIGATE
735 if (model->clpScaledMatrix())
736 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
737 #endif
738 // scaled
739 if (scalar == -1.0) {
740 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
741 double value = 0.0;
742 CoinBigIndex j;
743 const double * columnScale = model->columnScale();
744 for (j = columnStart[iColumn];
745 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
746 int iRow = row[j];
747 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
748 }
749 value *= columnScale[iColumn];
750 if (fabs(value) > zeroTolerance) {
751 index[numberNonZero++] = iColumn;
752 array[iColumn] = -value;
753 }
754 }
755 } else {
756 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
757 double value = 0.0;
758 CoinBigIndex j;
759 const double * columnScale = model->columnScale();
760 for (j = columnStart[iColumn];
761 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
762 int iRow = row[j];
763 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
764 }
765 value *= scalar * columnScale[iColumn];
766 if (fabs(value) > zeroTolerance) {
767 index[numberNonZero++] = iColumn;
768 array[iColumn] = value;
769 }
770 }
771 }
772 }
773 }
774 columnArray->setNumElements(numberNonZero);
775 y->setNumElements(0);
776 } else {
777 // do by row
778 rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
779 }
780 if (packed)
781 columnArray->setPackedMode(true);
782 if (0) {
783 columnArray->checkClean();
784 int numberNonZero = columnArray->getNumElements();
785 int * index = columnArray->getIndices();
786 double * array = columnArray->denseVector();
787 int i;
788 for (i = 0; i < numberNonZero; i++) {
789 int j = index[i];
790 double value;
791 if (packed)
792 value = array[i];
793 else
794 value = array[j];
795 printf("Ti %d %d %g\n", i, j, value);
796 }
797 }
798 }
799 //static int xA=0;
800 //static int xB=0;
801 //static int xC=0;
802 //static int xD=0;
803 //static double yA=0.0;
804 //static double yC=0.0;
805 /* Return <code>x * scalar * A + y</code> in <code>z</code>.
806 Note - If x packed mode - then z packed mode
807 This does by column and knows no gaps
808 Squashes small elements and knows about ClpSimplex */
809 void
transposeTimesByColumn(const ClpSimplex * model,double scalar,const CoinIndexedVector * rowArray,CoinIndexedVector * y,CoinIndexedVector * columnArray) const810 ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
811 const CoinIndexedVector * rowArray,
812 CoinIndexedVector * y,
813 CoinIndexedVector * columnArray) const
814 {
815 double * COIN_RESTRICT pi = rowArray->denseVector();
816 int numberNonZero = 0;
817 int * COIN_RESTRICT index = columnArray->getIndices();
818 double * COIN_RESTRICT array = columnArray->denseVector();
819 int numberInRowArray = rowArray->getNumElements();
820 // maybe I need one in OsiSimplex
821 double zeroTolerance = model->zeroTolerance();
822 bool packed = rowArray->packedMode();
823 // do by column
824 int iColumn;
825 // get matrix data pointers
826 const int * COIN_RESTRICT row = matrix_->getIndices();
827 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
828 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
829 const double * COIN_RESTRICT rowScale = model->rowScale();
830 assert (!y->getNumElements());
831 assert (numberActiveColumns_ > 0);
832 const ClpPackedMatrix * thisMatrix = this;
833 #if 0
834 ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
835 if (rowScale && scaledMatrix) {
836 rowScale = NULL;
837 // get matrix data pointers
838 row = scaledMatrix->getIndices();
839 columnStart = scaledMatrix->getVectorStarts();
840 elementByColumn = scaledMatrix->getElements();
841 thisMatrix = scaledMatrix;
842 //printf("scaledMatrix\n");
843 } else if (rowScale) {
844 //printf("no scaledMatrix\n");
845 } else {
846 //printf("no rowScale\n");
847 }
848 #endif
849 if (packed) {
850 // need to expand pi into y
851 assert(y->capacity() >= model->numberRows());
852 double * piOld = pi;
853 pi = y->denseVector();
854 const int * COIN_RESTRICT whichRow = rowArray->getIndices();
855 int i;
856 if (!rowScale) {
857 // modify pi so can collapse to one loop
858 if (scalar == -1.0) {
859 //yA += numberInRowArray;
860 for (i = 0; i < numberInRowArray; i++) {
861 int iRow = whichRow[i];
862 pi[iRow] = -piOld[i];
863 }
864 } else {
865 for (i = 0; i < numberInRowArray; i++) {
866 int iRow = whichRow[i];
867 pi[iRow] = scalar * piOld[i];
868 }
869 }
870 if (!columnCopy_) {
871 if ((model->specialOptions(), 131072) != 0) {
872 if(model->spareIntArray_[0] > 0) {
873 CoinIndexedVector * spareArray = model->rowArray(3);
874 // also do dualColumn stuff
875 double * spare = spareArray->denseVector();
876 int * spareIndex = spareArray->getIndices();
877 const double * reducedCost = model->djRegion(0);
878 double multiplier[] = { -1.0, 1.0};
879 double dualT = - model->currentDualTolerance();
880 double acceptablePivot = model->spareDoubleArray_[0];
881 // We can also see if infeasible or pivoting on free
882 double tentativeTheta = 1.0e15;
883 double upperTheta = 1.0e31;
884 double bestPossible = 0.0;
885 int addSequence = model->numberColumns();
886 const unsigned char * statusArray = model->statusArray() + addSequence;
887 int numberRemaining = 0;
888 assert (scalar == -1.0);
889 for (i = 0; i < numberInRowArray; i++) {
890 int iSequence = whichRow[i];
891 int iStatus = (statusArray[iSequence] & 3) - 1;
892 if (iStatus) {
893 double mult = multiplier[iStatus-1];
894 double alpha = piOld[i] * mult;
895 double oldValue;
896 double value;
897 if (alpha > 0.0) {
898 oldValue = reducedCost[iSequence] * mult;
899 value = oldValue - tentativeTheta * alpha;
900 if (value < dualT) {
901 bestPossible = CoinMax(bestPossible, alpha);
902 value = oldValue - upperTheta * alpha;
903 if (value < dualT && alpha >= acceptablePivot) {
904 upperTheta = (oldValue - dualT) / alpha;
905 //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
906 }
907 // add to list
908 spare[numberRemaining] = alpha * mult;
909 spareIndex[numberRemaining++] = iSequence + addSequence;
910 }
911 }
912 }
913 }
914 numberNonZero =
915 thisMatrix->gutsOfTransposeTimesUnscaled(pi,
916 columnArray->getIndices(),
917 columnArray->denseVector(),
918 model->statusArray(),
919 spareIndex,
920 spare,
921 model->djRegion(1),
922 upperTheta,
923 bestPossible,
924 acceptablePivot,
925 model->currentDualTolerance(),
926 numberRemaining,
927 zeroTolerance);
928 model->spareDoubleArray_[0] = upperTheta;
929 model->spareDoubleArray_[1] = bestPossible;
930 spareArray->setNumElements(numberRemaining);
931 // signal partially done
932 model->spareIntArray_[0] = -2;
933 } else {
934 numberNonZero =
935 thisMatrix->gutsOfTransposeTimesUnscaled(pi,
936 columnArray->getIndices(),
937 columnArray->denseVector(),
938 model->statusArray(),
939 zeroTolerance);
940 }
941 } else {
942 numberNonZero =
943 thisMatrix->gutsOfTransposeTimesUnscaled(pi,
944 columnArray->getIndices(),
945 columnArray->denseVector(),
946 zeroTolerance);
947 }
948 columnArray->setNumElements(numberNonZero);
949 //xA++;
950 } else {
951 columnCopy_->transposeTimes(model, pi, columnArray);
952 numberNonZero = columnArray->getNumElements();
953 //xB++;
954 }
955 } else {
956 #ifdef CLP_INVESTIGATE
957 if (model->clpScaledMatrix())
958 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
959 #endif
960 // scaled
961 // modify pi so can collapse to one loop
962 if (scalar == -1.0) {
963 //yC += numberInRowArray;
964 for (i = 0; i < numberInRowArray; i++) {
965 int iRow = whichRow[i];
966 pi[iRow] = -piOld[i] * rowScale[iRow];
967 }
968 } else {
969 for (i = 0; i < numberInRowArray; i++) {
970 int iRow = whichRow[i];
971 pi[iRow] = scalar * piOld[i] * rowScale[iRow];
972 }
973 }
974 const double * columnScale = model->columnScale();
975 if (!columnCopy_) {
976 if ((model->specialOptions(), 131072) != 0)
977 numberNonZero =
978 gutsOfTransposeTimesScaled(pi, columnScale,
979 columnArray->getIndices(),
980 columnArray->denseVector(),
981 model->statusArray(),
982 zeroTolerance);
983 else
984 numberNonZero =
985 gutsOfTransposeTimesScaled(pi, columnScale,
986 columnArray->getIndices(),
987 columnArray->denseVector(),
988 zeroTolerance);
989 columnArray->setNumElements(numberNonZero);
990 //xC++;
991 } else {
992 columnCopy_->transposeTimes(model, pi, columnArray);
993 numberNonZero = columnArray->getNumElements();
994 //xD++;
995 }
996 }
997 // zero out
998 int numberRows = model->numberRows();
999 if (numberInRowArray * 4 < numberRows) {
1000 for (i = 0; i < numberInRowArray; i++) {
1001 int iRow = whichRow[i];
1002 pi[iRow] = 0.0;
1003 }
1004 } else {
1005 CoinZeroN(pi, numberRows);
1006 }
1007 //int kA=xA+xB;
1008 //int kC=xC+xD;
1009 //if ((kA+kC)%10000==0)
1010 //printf("AA %d %d %g, CC %d %d %g\n",
1011 // xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0);
1012 } else {
1013 if (!rowScale) {
1014 if (scalar == -1.0) {
1015 double value = 0.0;
1016 CoinBigIndex j;
1017 CoinBigIndex end = columnStart[1];
1018 for (j = columnStart[0]; j < end; j++) {
1019 int iRow = row[j];
1020 value += pi[iRow] * elementByColumn[j];
1021 }
1022 for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
1023 CoinBigIndex start = end;
1024 end = columnStart[iColumn+2];
1025 if (fabs(value) > zeroTolerance) {
1026 array[iColumn] = -value;
1027 index[numberNonZero++] = iColumn;
1028 }
1029 value = 0.0;
1030 for (j = start; j < end; j++) {
1031 int iRow = row[j];
1032 value += pi[iRow] * elementByColumn[j];
1033 }
1034 }
1035 if (fabs(value) > zeroTolerance) {
1036 array[iColumn] = -value;
1037 index[numberNonZero++] = iColumn;
1038 }
1039 } else {
1040 double value = 0.0;
1041 CoinBigIndex j;
1042 CoinBigIndex end = columnStart[1];
1043 for (j = columnStart[0]; j < end; j++) {
1044 int iRow = row[j];
1045 value += pi[iRow] * elementByColumn[j];
1046 }
1047 for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
1048 value *= scalar;
1049 CoinBigIndex start = end;
1050 end = columnStart[iColumn+2];
1051 if (fabs(value) > zeroTolerance) {
1052 array[iColumn] = value;
1053 index[numberNonZero++] = iColumn;
1054 }
1055 value = 0.0;
1056 for (j = start; j < end; j++) {
1057 int iRow = row[j];
1058 value += pi[iRow] * elementByColumn[j];
1059 }
1060 }
1061 value *= scalar;
1062 if (fabs(value) > zeroTolerance) {
1063 array[iColumn] = value;
1064 index[numberNonZero++] = iColumn;
1065 }
1066 }
1067 } else {
1068 #ifdef CLP_INVESTIGATE
1069 if (model->clpScaledMatrix())
1070 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
1071 #endif
1072 // scaled
1073 if (scalar == -1.0) {
1074 const double * columnScale = model->columnScale();
1075 double value = 0.0;
1076 double scale = columnScale[0];
1077 CoinBigIndex j;
1078 CoinBigIndex end = columnStart[1];
1079 for (j = columnStart[0]; j < end; j++) {
1080 int iRow = row[j];
1081 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1082 }
1083 for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
1084 value *= scale;
1085 CoinBigIndex start = end;
1086 end = columnStart[iColumn+2];
1087 scale = columnScale[iColumn+1];
1088 if (fabs(value) > zeroTolerance) {
1089 array[iColumn] = -value;
1090 index[numberNonZero++] = iColumn;
1091 }
1092 value = 0.0;
1093 for (j = start; j < end; j++) {
1094 int iRow = row[j];
1095 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1096 }
1097 }
1098 value *= scale;
1099 if (fabs(value) > zeroTolerance) {
1100 array[iColumn] = -value;
1101 index[numberNonZero++] = iColumn;
1102 }
1103 } else {
1104 const double * columnScale = model->columnScale();
1105 double value = 0.0;
1106 double scale = columnScale[0] * scalar;
1107 CoinBigIndex j;
1108 CoinBigIndex end = columnStart[1];
1109 for (j = columnStart[0]; j < end; j++) {
1110 int iRow = row[j];
1111 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1112 }
1113 for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
1114 value *= scale;
1115 CoinBigIndex start = end;
1116 end = columnStart[iColumn+2];
1117 scale = columnScale[iColumn+1] * scalar;
1118 if (fabs(value) > zeroTolerance) {
1119 array[iColumn] = value;
1120 index[numberNonZero++] = iColumn;
1121 }
1122 value = 0.0;
1123 for (j = start; j < end; j++) {
1124 int iRow = row[j];
1125 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1126 }
1127 }
1128 value *= scale;
1129 if (fabs(value) > zeroTolerance) {
1130 array[iColumn] = value;
1131 index[numberNonZero++] = iColumn;
1132 }
1133 }
1134 }
1135 }
1136 columnArray->setNumElements(numberNonZero);
1137 y->setNumElements(0);
1138 if (packed)
1139 columnArray->setPackedMode(true);
1140 }
1141 /* Return <code>x * A + y</code> in <code>z</code>.
1142 Squashes small elements and knows about ClpSimplex */
1143 void
transposeTimesByRow(const ClpSimplex * model,double scalar,const CoinIndexedVector * rowArray,CoinIndexedVector * y,CoinIndexedVector * columnArray) const1144 ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
1145 const CoinIndexedVector * rowArray,
1146 CoinIndexedVector * y,
1147 CoinIndexedVector * columnArray) const
1148 {
1149 columnArray->clear();
1150 double * pi = rowArray->denseVector();
1151 int numberNonZero = 0;
1152 int * index = columnArray->getIndices();
1153 double * array = columnArray->denseVector();
1154 int numberInRowArray = rowArray->getNumElements();
1155 // maybe I need one in OsiSimplex
1156 double zeroTolerance = model->zeroTolerance();
1157 const int * column = matrix_->getIndices();
1158 const CoinBigIndex * rowStart = getVectorStarts();
1159 const double * element = getElements();
1160 const int * whichRow = rowArray->getIndices();
1161 bool packed = rowArray->packedMode();
1162 if (numberInRowArray > 2) {
1163 // do by rows
1164 // ** Row copy is already scaled
1165 int iRow;
1166 int i;
1167 int numberOriginal = 0;
1168 if (packed) {
1169 int * index = columnArray->getIndices();
1170 double * array = columnArray->denseVector();
1171 #if 0
1172 {
1173 double * array2 = y->denseVector();
1174 int numberColumns = matrix_->getNumCols();
1175 for (int i=0;i<numberColumns;i++) {
1176 assert(!array[i]);
1177 assert(!array2[i]);
1178 }
1179 }
1180 #endif
1181 //#define COIN_SPARSE_MATRIX 1
1182 #if COIN_SPARSE_MATRIX
1183 assert (!y->getNumElements());
1184 #if COIN_SPARSE_MATRIX != 2
1185 // and set up mark as char array
1186 char * marked = reinterpret_cast<char *> (index+columnArray->capacity());
1187 int * lookup = y->getIndices();
1188 #ifndef NDEBUG
1189 //int numberColumns = matrix_->getNumCols();
1190 //for (int i=0;i<numberColumns;i++)
1191 //assert(!marked[i]);
1192 #endif
1193 numberNonZero=gutsOfTransposeTimesByRowGE3a(rowArray,index,array,
1194 lookup,marked,zeroTolerance,scalar);
1195 #else
1196 double * array2 = y->denseVector();
1197 numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
1198 array2,zeroTolerance,scalar);
1199 #endif
1200 #else
1201 int numberColumns = matrix_->getNumCols();
1202 numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
1203 numberColumns, zeroTolerance, scalar);
1204 #endif
1205 columnArray->setNumElements(numberNonZero);
1206 } else {
1207 double * markVector = y->denseVector();
1208 numberNonZero = 0;
1209 // and set up mark as char array
1210 char * marked = reinterpret_cast<char *> (markVector);
1211 for (i = 0; i < numberOriginal; i++) {
1212 int iColumn = index[i];
1213 marked[iColumn] = 0;
1214 }
1215
1216 for (i = 0; i < numberInRowArray; i++) {
1217 iRow = whichRow[i];
1218 double value = pi[iRow] * scalar;
1219 CoinBigIndex j;
1220 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
1221 int iColumn = column[j];
1222 if (!marked[iColumn]) {
1223 marked[iColumn] = 1;
1224 index[numberNonZero++] = iColumn;
1225 }
1226 array[iColumn] += value * element[j];
1227 }
1228 }
1229 // get rid of tiny values and zero out marked
1230 numberOriginal = numberNonZero;
1231 numberNonZero = 0;
1232 for (i = 0; i < numberOriginal; i++) {
1233 int iColumn = index[i];
1234 marked[iColumn] = 0;
1235 if (fabs(array[iColumn]) > zeroTolerance) {
1236 index[numberNonZero++] = iColumn;
1237 } else {
1238 array[iColumn] = 0.0;
1239 }
1240 }
1241 }
1242 } else if (numberInRowArray == 2) {
1243 // do by rows when two rows
1244 int numberOriginal;
1245 int i;
1246 CoinBigIndex j;
1247 numberNonZero = 0;
1248
1249 double value;
1250 if (packed) {
1251 gutsOfTransposeTimesByRowEQ2(rowArray, columnArray, y, zeroTolerance, scalar);
1252 numberNonZero = columnArray->getNumElements();
1253 } else {
1254 int iRow = whichRow[0];
1255 value = pi[iRow] * scalar;
1256 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
1257 int iColumn = column[j];
1258 double value2 = value * element[j];
1259 index[numberNonZero++] = iColumn;
1260 array[iColumn] = value2;
1261 }
1262 iRow = whichRow[1];
1263 value = pi[iRow] * scalar;
1264 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
1265 int iColumn = column[j];
1266 double value2 = value * element[j];
1267 // I am assuming no zeros in matrix
1268 if (array[iColumn])
1269 value2 += array[iColumn];
1270 else
1271 index[numberNonZero++] = iColumn;
1272 array[iColumn] = value2;
1273 }
1274 // get rid of tiny values and zero out marked
1275 numberOriginal = numberNonZero;
1276 numberNonZero = 0;
1277 for (i = 0; i < numberOriginal; i++) {
1278 int iColumn = index[i];
1279 if (fabs(array[iColumn]) > zeroTolerance) {
1280 index[numberNonZero++] = iColumn;
1281 } else {
1282 array[iColumn] = 0.0;
1283 }
1284 }
1285 }
1286 } else if (numberInRowArray == 1) {
1287 // Just one row
1288 int iRow = rowArray->getIndices()[0];
1289 numberNonZero = 0;
1290 CoinBigIndex j;
1291 if (packed) {
1292 gutsOfTransposeTimesByRowEQ1(rowArray, columnArray, zeroTolerance, scalar);
1293 numberNonZero = columnArray->getNumElements();
1294 } else {
1295 double value = pi[iRow] * scalar;
1296 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
1297 int iColumn = column[j];
1298 double value2 = value * element[j];
1299 if (fabs(value2) > zeroTolerance) {
1300 index[numberNonZero++] = iColumn;
1301 array[iColumn] = value2;
1302 }
1303 }
1304 }
1305 }
1306 columnArray->setNumElements(numberNonZero);
1307 y->setNumElements(0);
1308 }
1309 // Meat of transposeTimes by column when not scaled
1310 int
gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,int * COIN_RESTRICT index,double * COIN_RESTRICT array,const double zeroTolerance) const1311 ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
1312 int * COIN_RESTRICT index,
1313 double * COIN_RESTRICT array,
1314 const double zeroTolerance) const
1315 {
1316 int numberNonZero = 0;
1317 // get matrix data pointers
1318 const int * COIN_RESTRICT row = matrix_->getIndices();
1319 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1320 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1321 #if 1 //ndef INTEL_MKL
1322 double value = 0.0;
1323 CoinBigIndex j;
1324 CoinBigIndex end = columnStart[1];
1325 for (j = columnStart[0]; j < end; j++) {
1326 int iRow = row[j];
1327 value += pi[iRow] * elementByColumn[j];
1328 }
1329 int iColumn;
1330 for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
1331 CoinBigIndex start = end;
1332 end = columnStart[iColumn+2];
1333 if (fabs(value) > zeroTolerance) {
1334 array[numberNonZero] = value;
1335 index[numberNonZero++] = iColumn;
1336 }
1337 value = 0.0;
1338 for (j = start; j < end; j++) {
1339 int iRow = row[j];
1340 value += pi[iRow] * elementByColumn[j];
1341 }
1342 }
1343 if (fabs(value) > zeroTolerance) {
1344 array[numberNonZero] = value;
1345 index[numberNonZero++] = iColumn;
1346 }
1347 #else
1348 char transA = 'N';
1349 //int numberRows = matrix_->getNumRows();
1350 mkl_cspblas_dcsrgemv(&transA, const_cast<int *>(&numberActiveColumns_),
1351 const_cast<double *>(elementByColumn),
1352 const_cast<int *>(columnStart),
1353 const_cast<int *>(row),
1354 const_cast<double *>(pi), array);
1355 int iColumn;
1356 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
1357 double value = array[iColumn];
1358 if (value) {
1359 array[iColumn] = 0.0;
1360 if (fabs(value) > zeroTolerance) {
1361 array[numberNonZero] = value;
1362 index[numberNonZero++] = iColumn;
1363 }
1364 }
1365 }
1366 #endif
1367 return numberNonZero;
1368 }
1369 // Meat of transposeTimes by column when scaled
1370 int
gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,const double * COIN_RESTRICT columnScale,int * COIN_RESTRICT index,double * COIN_RESTRICT array,const double zeroTolerance) const1371 ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
1372 const double * COIN_RESTRICT columnScale,
1373 int * COIN_RESTRICT index,
1374 double * COIN_RESTRICT array,
1375 const double zeroTolerance) const
1376 {
1377 int numberNonZero = 0;
1378 // get matrix data pointers
1379 const int * COIN_RESTRICT row = matrix_->getIndices();
1380 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1381 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1382 double value = 0.0;
1383 double scale = columnScale[0];
1384 CoinBigIndex j;
1385 CoinBigIndex end = columnStart[1];
1386 for (j = columnStart[0]; j < end; j++) {
1387 int iRow = row[j];
1388 value += pi[iRow] * elementByColumn[j];
1389 }
1390 int iColumn;
1391 for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
1392 value *= scale;
1393 CoinBigIndex start = end;
1394 scale = columnScale[iColumn+1];
1395 end = columnStart[iColumn+2];
1396 if (fabs(value) > zeroTolerance) {
1397 array[numberNonZero] = value;
1398 index[numberNonZero++] = iColumn;
1399 }
1400 value = 0.0;
1401 for (j = start; j < end; j++) {
1402 int iRow = row[j];
1403 value += pi[iRow] * elementByColumn[j];
1404 }
1405 }
1406 value *= scale;
1407 if (fabs(value) > zeroTolerance) {
1408 array[numberNonZero] = value;
1409 index[numberNonZero++] = iColumn;
1410 }
1411 return numberNonZero;
1412 }
1413 // Meat of transposeTimes by column when not scaled
1414 int
gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,int * COIN_RESTRICT index,double * COIN_RESTRICT array,const unsigned char * COIN_RESTRICT status,const double zeroTolerance) const1415 ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
1416 int * COIN_RESTRICT index,
1417 double * COIN_RESTRICT array,
1418 const unsigned char * COIN_RESTRICT status,
1419 const double zeroTolerance) const
1420 {
1421 int numberNonZero = 0;
1422 // get matrix data pointers
1423 const int * COIN_RESTRICT row = matrix_->getIndices();
1424 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1425 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1426 double value = 0.0;
1427 int jColumn = -1;
1428 for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
1429 bool wanted = ((status[iColumn] & 3) != 1);
1430 if (fabs(value) > zeroTolerance) {
1431 array[numberNonZero] = value;
1432 index[numberNonZero++] = jColumn;
1433 }
1434 value = 0.0;
1435 if (wanted) {
1436 CoinBigIndex start = columnStart[iColumn];
1437 CoinBigIndex end = columnStart[iColumn+1];
1438 jColumn = iColumn;
1439 int n = end - start;
1440 bool odd = (n & 1) != 0;
1441 n = n >> 1;
1442 const int * COIN_RESTRICT rowThis = row + start;
1443 const double * COIN_RESTRICT elementThis = elementByColumn + start;
1444 for (; n; n--) {
1445 int iRow0 = *rowThis;
1446 int iRow1 = *(rowThis + 1);
1447 rowThis += 2;
1448 value += pi[iRow0] * (*elementThis);
1449 value += pi[iRow1] * (*(elementThis + 1));
1450 elementThis += 2;
1451 }
1452 if (odd) {
1453 int iRow = *rowThis;
1454 value += pi[iRow] * (*elementThis);
1455 }
1456 }
1457 }
1458 if (fabs(value) > zeroTolerance) {
1459 array[numberNonZero] = value;
1460 index[numberNonZero++] = jColumn;
1461 }
1462 return numberNonZero;
1463 }
1464 /* Meat of transposeTimes by column when not scaled and skipping
1465 and doing part of dualColumn */
1466 int
gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,int * COIN_RESTRICT index,double * COIN_RESTRICT array,const unsigned char * COIN_RESTRICT status,int * COIN_RESTRICT spareIndex,double * COIN_RESTRICT spareArray,const double * COIN_RESTRICT reducedCost,double & upperThetaP,double & bestPossibleP,double acceptablePivot,double dualTolerance,int & numberRemainingP,const double zeroTolerance) const1467 ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
1468 int * COIN_RESTRICT index,
1469 double * COIN_RESTRICT array,
1470 const unsigned char * COIN_RESTRICT status,
1471 int * COIN_RESTRICT spareIndex,
1472 double * COIN_RESTRICT spareArray,
1473 const double * COIN_RESTRICT reducedCost,
1474 double & upperThetaP,
1475 double & bestPossibleP,
1476 double acceptablePivot,
1477 double dualTolerance,
1478 int & numberRemainingP,
1479 const double zeroTolerance) const
1480 {
1481 double tentativeTheta = 1.0e15;
1482 int numberRemaining = numberRemainingP;
1483 double upperTheta = upperThetaP;
1484 double bestPossible = bestPossibleP;
1485 int numberNonZero = 0;
1486 // get matrix data pointers
1487 const int * COIN_RESTRICT row = matrix_->getIndices();
1488 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1489 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1490 double multiplier[] = { -1.0, 1.0};
1491 double dualT = - dualTolerance;
1492 for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
1493 int wanted = (status[iColumn] & 3) - 1;
1494 if (wanted) {
1495 double value = 0.0;
1496 CoinBigIndex start = columnStart[iColumn];
1497 CoinBigIndex end = columnStart[iColumn+1];
1498 int n = end - start;
1499 #if 1
1500 bool odd = (n & 1) != 0;
1501 n = n >> 1;
1502 const int * COIN_RESTRICT rowThis = row + start;
1503 const double * COIN_RESTRICT elementThis = elementByColumn + start;
1504 for (; n; n--) {
1505 int iRow0 = *rowThis;
1506 int iRow1 = *(rowThis + 1);
1507 rowThis += 2;
1508 value += pi[iRow0] * (*elementThis);
1509 value += pi[iRow1] * (*(elementThis + 1));
1510 elementThis += 2;
1511 }
1512 if (odd) {
1513 int iRow = *rowThis;
1514 value += pi[iRow] * (*elementThis);
1515 }
1516 #else
1517 const int * COIN_RESTRICT rowThis = &row[end-16];
1518 const double * COIN_RESTRICT elementThis = &elementByColumn[end-16];
1519 bool odd = (n & 1) != 0;
1520 n = n >> 1;
1521 double value2 = 0.0;
1522 if (odd) {
1523 int iRow = row[start];
1524 value2 = pi[iRow] * elementByColumn[start];
1525 }
1526 switch (n) {
1527 default: {
1528 if (odd) {
1529 start++;
1530 }
1531 n -= 8;
1532 for (; n; n--) {
1533 int iRow0 = row[start];
1534 int iRow1 = row[start+1];
1535 value += pi[iRow0] * elementByColumn[start];
1536 value2 += pi[iRow1] * elementByColumn[start+1];
1537 start += 2;
1538 }
1539 case 8: {
1540 int iRow0 = rowThis[16-16];
1541 int iRow1 = rowThis[16-15];
1542 value += pi[iRow0] * elementThis[16-16];
1543 value2 += pi[iRow1] * elementThis[16-15];
1544 }
1545 case 7: {
1546 int iRow0 = rowThis[16-14];
1547 int iRow1 = rowThis[16-13];
1548 value += pi[iRow0] * elementThis[16-14];
1549 value2 += pi[iRow1] * elementThis[16-13];
1550 }
1551 case 6: {
1552 int iRow0 = rowThis[16-12];
1553 int iRow1 = rowThis[16-11];
1554 value += pi[iRow0] * elementThis[16-12];
1555 value2 += pi[iRow1] * elementThis[16-11];
1556 }
1557 case 5: {
1558 int iRow0 = rowThis[16-10];
1559 int iRow1 = rowThis[16-9];
1560 value += pi[iRow0] * elementThis[16-10];
1561 value2 += pi[iRow1] * elementThis[16-9];
1562 }
1563 case 4: {
1564 int iRow0 = rowThis[16-8];
1565 int iRow1 = rowThis[16-7];
1566 value += pi[iRow0] * elementThis[16-8];
1567 value2 += pi[iRow1] * elementThis[16-7];
1568 }
1569 case 3: {
1570 int iRow0 = rowThis[16-6];
1571 int iRow1 = rowThis[16-5];
1572 value += pi[iRow0] * elementThis[16-6];
1573 value2 += pi[iRow1] * elementThis[16-5];
1574 }
1575 case 2: {
1576 int iRow0 = rowThis[16-4];
1577 int iRow1 = rowThis[16-3];
1578 value += pi[iRow0] * elementThis[16-4];
1579 value2 += pi[iRow1] * elementThis[16-3];
1580 }
1581 case 1: {
1582 int iRow0 = rowThis[16-2];
1583 int iRow1 = rowThis[16-1];
1584 value += pi[iRow0] * elementThis[16-2];
1585 value2 += pi[iRow1] * elementThis[16-1];
1586 }
1587 case 0:
1588 ;
1589 }
1590 }
1591 value += value2;
1592 #endif
1593 if (fabs(value) > zeroTolerance) {
1594 double mult = multiplier[wanted-1];
1595 double alpha = value * mult;
1596 array[numberNonZero] = value;
1597 index[numberNonZero++] = iColumn;
1598 if (alpha > 0.0) {
1599 double oldValue = reducedCost[iColumn] * mult;
1600 double value = oldValue - tentativeTheta * alpha;
1601 if (value < dualT) {
1602 bestPossible = CoinMax(bestPossible, alpha);
1603 value = oldValue - upperTheta * alpha;
1604 if (value < dualT && alpha >= acceptablePivot) {
1605 upperTheta = (oldValue - dualT) / alpha;
1606 //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
1607 }
1608 // add to list
1609 spareArray[numberRemaining] = alpha * mult;
1610 spareIndex[numberRemaining++] = iColumn;
1611 }
1612 }
1613 }
1614 }
1615 }
1616 numberRemainingP = numberRemaining;
1617 upperThetaP = upperTheta;
1618 bestPossibleP = bestPossible;
1619 return numberNonZero;
1620 }
1621 // Meat of transposeTimes by column when scaled
1622 int
gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,const double * COIN_RESTRICT columnScale,int * COIN_RESTRICT index,double * COIN_RESTRICT array,const unsigned char * COIN_RESTRICT status,const double zeroTolerance) const1623 ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
1624 const double * COIN_RESTRICT columnScale,
1625 int * COIN_RESTRICT index,
1626 double * COIN_RESTRICT array,
1627 const unsigned char * COIN_RESTRICT status, const double zeroTolerance) const
1628 {
1629 int numberNonZero = 0;
1630 // get matrix data pointers
1631 const int * COIN_RESTRICT row = matrix_->getIndices();
1632 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1633 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1634 double value = 0.0;
1635 int jColumn = -1;
1636 for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
1637 bool wanted = ((status[iColumn] & 3) != 1);
1638 if (fabs(value) > zeroTolerance) {
1639 array[numberNonZero] = value;
1640 index[numberNonZero++] = jColumn;
1641 }
1642 value = 0.0;
1643 if (wanted) {
1644 double scale = columnScale[iColumn];
1645 CoinBigIndex start = columnStart[iColumn];
1646 CoinBigIndex end = columnStart[iColumn+1];
1647 jColumn = iColumn;
1648 for (CoinBigIndex j = start; j < end; j++) {
1649 int iRow = row[j];
1650 value += pi[iRow] * elementByColumn[j];
1651 }
1652 value *= scale;
1653 }
1654 }
1655 if (fabs(value) > zeroTolerance) {
1656 array[numberNonZero] = value;
1657 index[numberNonZero++] = jColumn;
1658 }
1659 return numberNonZero;
1660 }
1661 // Meat of transposeTimes by row n > K if packed - returns number nonzero
1662 int
gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,int * COIN_RESTRICT index,double * COIN_RESTRICT output,int numberColumns,const double tolerance,const double scalar) const1663 ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,
1664 int * COIN_RESTRICT index,
1665 double * COIN_RESTRICT output,
1666 int numberColumns,
1667 const double tolerance,
1668 const double scalar) const
1669 {
1670 const double * COIN_RESTRICT pi = piVector->denseVector();
1671 int numberInRowArray = piVector->getNumElements();
1672 const int * COIN_RESTRICT column = matrix_->getIndices();
1673 const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
1674 const double * COIN_RESTRICT element = matrix_->getElements();
1675 const int * COIN_RESTRICT whichRow = piVector->getIndices();
1676 // ** Row copy is already scaled
1677 for (int i = 0; i < numberInRowArray; i++) {
1678 int iRow = whichRow[i];
1679 double value = pi[i] * scalar;
1680 CoinBigIndex start = rowStart[iRow];
1681 CoinBigIndex end = rowStart[iRow+1];
1682 int n = end - start;
1683 const int * COIN_RESTRICT columnThis = column + start;
1684 const double * COIN_RESTRICT elementThis = element + start;
1685
1686 // could do by twos
1687 for (; n; n--) {
1688 int iColumn = *columnThis;
1689 columnThis++;
1690 double elValue = *elementThis;
1691 elementThis++;
1692 elValue *= value;
1693 output[iColumn] += elValue;
1694 }
1695 }
1696 // get rid of tiny values and count
1697 int numberNonZero = 0;
1698 for (int i = 0; i < numberColumns; i++) {
1699 double value = output[i];
1700 if (value) {
1701 output[i] = 0.0;
1702 if (fabs(value) > tolerance) {
1703 output[numberNonZero] = value;
1704 index[numberNonZero++] = i;
1705 }
1706 }
1707 }
1708 #ifndef NDEBUG
1709 for (int i = numberNonZero; i < numberColumns; i++)
1710 assert(!output[i]);
1711 #endif
1712 return numberNonZero;
1713 }
1714 // Meat of transposeTimes by row n == 2 if packed
1715 void
gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector,CoinIndexedVector * output,CoinIndexedVector * spareVector,const double tolerance,const double scalar) const1716 ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
1717 CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
1718 {
1719 double * pi = piVector->denseVector();
1720 int numberNonZero = 0;
1721 int * index = output->getIndices();
1722 double * array = output->denseVector();
1723 const int * column = matrix_->getIndices();
1724 const CoinBigIndex * rowStart = matrix_->getVectorStarts();
1725 const double * element = matrix_->getElements();
1726 const int * whichRow = piVector->getIndices();
1727 int iRow0 = whichRow[0];
1728 int iRow1 = whichRow[1];
1729 double pi0 = pi[0];
1730 double pi1 = pi[1];
1731 if (rowStart[iRow0+1] - rowStart[iRow0] >
1732 rowStart[iRow1+1] - rowStart[iRow1]) {
1733 // do one with fewer first
1734 iRow0 = iRow1;
1735 iRow1 = whichRow[0];
1736 pi0 = pi1;
1737 pi1 = pi[0];
1738 }
1739 // and set up mark as char array
1740 char * marked = reinterpret_cast<char *> (index + output->capacity());
1741 int * lookup = spareVector->getIndices();
1742 double value = pi0 * scalar;
1743 CoinBigIndex j;
1744 for (j = rowStart[iRow0]; j < rowStart[iRow0+1]; j++) {
1745 int iColumn = column[j];
1746 double elValue = element[j];
1747 double value2 = value * elValue;
1748 array[numberNonZero] = value2;
1749 marked[iColumn] = 1;
1750 lookup[iColumn] = numberNonZero;
1751 index[numberNonZero++] = iColumn;
1752 }
1753 int numberOriginal = numberNonZero;
1754 value = pi1 * scalar;
1755 for (j = rowStart[iRow1]; j < rowStart[iRow1+1]; j++) {
1756 int iColumn = column[j];
1757 double elValue = element[j];
1758 double value2 = value * elValue;
1759 // I am assuming no zeros in matrix
1760 if (marked[iColumn]) {
1761 int iLookup = lookup[iColumn];
1762 array[iLookup] += value2;
1763 } else {
1764 if (fabs(value2) > tolerance) {
1765 array[numberNonZero] = value2;
1766 index[numberNonZero++] = iColumn;
1767 }
1768 }
1769 }
1770 // get rid of tiny values and zero out marked
1771 int i;
1772 int iFirst = numberNonZero;
1773 for (i = 0; i < numberOriginal; i++) {
1774 int iColumn = index[i];
1775 marked[iColumn] = 0;
1776 if (fabs(array[i]) <= tolerance) {
1777 if (numberNonZero > numberOriginal) {
1778 numberNonZero--;
1779 double value = array[numberNonZero];
1780 array[numberNonZero] = 0.0;
1781 array[i] = value;
1782 index[i] = index[numberNonZero];
1783 } else {
1784 iFirst = i;
1785 }
1786 }
1787 }
1788
1789 if (iFirst < numberNonZero) {
1790 int n = iFirst;
1791 for (i = n; i < numberOriginal; i++) {
1792 int iColumn = index[i];
1793 double value = array[i];
1794 array[i] = 0.0;
1795 if (fabs(value) > tolerance) {
1796 array[n] = value;
1797 index[n++] = iColumn;
1798 }
1799 }
1800 for (; i < numberNonZero; i++) {
1801 int iColumn = index[i];
1802 double value = array[i];
1803 array[i] = 0.0;
1804 array[n] = value;
1805 index[n++] = iColumn;
1806 }
1807 numberNonZero = n;
1808 }
1809 output->setNumElements(numberNonZero);
1810 spareVector->setNumElements(0);
1811 }
1812 // Meat of transposeTimes by row n == 1 if packed
1813 void
gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector,CoinIndexedVector * output,const double tolerance,const double scalar) const1814 ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
1815 const double tolerance, const double scalar) const
1816 {
1817 double * pi = piVector->denseVector();
1818 int numberNonZero = 0;
1819 int * index = output->getIndices();
1820 double * array = output->denseVector();
1821 const int * column = matrix_->getIndices();
1822 const CoinBigIndex * rowStart = matrix_->getVectorStarts();
1823 const double * element = matrix_->getElements();
1824 int iRow = piVector->getIndices()[0];
1825 numberNonZero = 0;
1826 CoinBigIndex j;
1827 double value = pi[0] * scalar;
1828 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
1829 int iColumn = column[j];
1830 double elValue = element[j];
1831 double value2 = value * elValue;
1832 if (fabs(value2) > tolerance) {
1833 array[numberNonZero] = value2;
1834 index[numberNonZero++] = iColumn;
1835 }
1836 }
1837 output->setNumElements(numberNonZero);
1838 }
1839 /* Return <code>x *A in <code>z</code> but
1840 just for indices in y.
1841 Squashes small elements and knows about ClpSimplex */
1842 void
subsetTransposeTimes(const ClpSimplex * model,const CoinIndexedVector * rowArray,const CoinIndexedVector * y,CoinIndexedVector * columnArray) const1843 ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
1844 const CoinIndexedVector * rowArray,
1845 const CoinIndexedVector * y,
1846 CoinIndexedVector * columnArray) const
1847 {
1848 columnArray->clear();
1849 double * COIN_RESTRICT pi = rowArray->denseVector();
1850 double * COIN_RESTRICT array = columnArray->denseVector();
1851 int jColumn;
1852 // get matrix data pointers
1853 const int * COIN_RESTRICT row = matrix_->getIndices();
1854 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1855 const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
1856 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1857 const double * COIN_RESTRICT rowScale = model->rowScale();
1858 int numberToDo = y->getNumElements();
1859 const int * COIN_RESTRICT which = y->getIndices();
1860 assert (!rowArray->packedMode());
1861 columnArray->setPacked();
1862 ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
1863 int flags = flags_;
1864 if (rowScale && scaledMatrix && !(scaledMatrix->flags() & 2)) {
1865 flags = 0;
1866 rowScale = NULL;
1867 // get matrix data pointers
1868 row = scaledMatrix->getIndices();
1869 columnStart = scaledMatrix->getVectorStarts();
1870 elementByColumn = scaledMatrix->getElements();
1871 }
1872 if (!(flags & 2) && numberToDo>2) {
1873 // no gaps
1874 if (!rowScale) {
1875 int iColumn = which[0];
1876 double value = 0.0;
1877 CoinBigIndex j;
1878 int columnNext = which[1];
1879 CoinBigIndex startNext=columnStart[columnNext];
1880 //coin_prefetch_const(row+startNext);
1881 //coin_prefetch_const(elementByColumn+startNext);
1882 CoinBigIndex endNext=columnStart[columnNext+1];
1883 for (j = columnStart[iColumn];
1884 j < columnStart[iColumn+1]; j++) {
1885 int iRow = row[j];
1886 value += pi[iRow] * elementByColumn[j];
1887 }
1888 for (jColumn = 0; jColumn < numberToDo - 2; jColumn++) {
1889 CoinBigIndex start = startNext;
1890 CoinBigIndex end = endNext;
1891 columnNext = which[jColumn+2];
1892 startNext=columnStart[columnNext];
1893 //coin_prefetch_const(row+startNext);
1894 //coin_prefetch_const(elementByColumn+startNext);
1895 endNext=columnStart[columnNext+1];
1896 array[jColumn] = value;
1897 value = 0.0;
1898 for (j = start; j < end; j++) {
1899 int iRow = row[j];
1900 value += pi[iRow] * elementByColumn[j];
1901 }
1902 }
1903 array[jColumn++] = value;
1904 value = 0.0;
1905 for (j = startNext; j < endNext; j++) {
1906 int iRow = row[j];
1907 value += pi[iRow] * elementByColumn[j];
1908 }
1909 array[jColumn] = value;
1910 } else {
1911 #ifdef CLP_INVESTIGATE
1912 if (model->clpScaledMatrix())
1913 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
1914 #endif
1915 // scaled
1916 const double * columnScale = model->columnScale();
1917 int iColumn = which[0];
1918 double value = 0.0;
1919 double scale = columnScale[iColumn];
1920 CoinBigIndex j;
1921 for (j = columnStart[iColumn];
1922 j < columnStart[iColumn+1]; j++) {
1923 int iRow = row[j];
1924 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1925 }
1926 for (jColumn = 0; jColumn < numberToDo - 1; jColumn++) {
1927 int iColumn = which[jColumn+1];
1928 value *= scale;
1929 scale = columnScale[iColumn];
1930 CoinBigIndex start = columnStart[iColumn];
1931 CoinBigIndex end = columnStart[iColumn+1];
1932 array[jColumn] = value;
1933 value = 0.0;
1934 for (j = start; j < end; j++) {
1935 int iRow = row[j];
1936 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1937 }
1938 }
1939 value *= scale;
1940 array[jColumn] = value;
1941 }
1942 } else if (numberToDo) {
1943 // gaps
1944 if (!rowScale) {
1945 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
1946 int iColumn = which[jColumn];
1947 double value = 0.0;
1948 CoinBigIndex j;
1949 for (j = columnStart[iColumn];
1950 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1951 int iRow = row[j];
1952 value += pi[iRow] * elementByColumn[j];
1953 }
1954 array[jColumn] = value;
1955 }
1956 } else {
1957 #ifdef CLP_INVESTIGATE
1958 if (model->clpScaledMatrix())
1959 printf("scaledMatrix_ at %d of ClpPackedMatrix - flags %d (%d) n %d\n",
1960 __LINE__, flags_, model->clpScaledMatrix()->flags(), numberToDo);
1961 #endif
1962 // scaled
1963 const double * columnScale = model->columnScale();
1964 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
1965 int iColumn = which[jColumn];
1966 double value = 0.0;
1967 CoinBigIndex j;
1968 for (j = columnStart[iColumn];
1969 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1970 int iRow = row[j];
1971 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1972 }
1973 value *= columnScale[iColumn];
1974 array[jColumn] = value;
1975 }
1976 }
1977 }
1978 }
1979 /* Returns true if can combine transposeTimes and subsetTransposeTimes
1980 and if it would be faster */
1981 bool
canCombine(const ClpSimplex * model,const CoinIndexedVector * pi) const1982 ClpPackedMatrix::canCombine(const ClpSimplex * model,
1983 const CoinIndexedVector * pi) const
1984 {
1985 int numberInRowArray = pi->getNumElements();
1986 int numberRows = model->numberRows();
1987 bool packed = pi->packedMode();
1988 // factor should be smaller if doing both with two pi vectors
1989 double factor = 0.30;
1990 // We may not want to do by row if there may be cache problems
1991 // It would be nice to find L2 cache size - for moment 512K
1992 // Be slightly optimistic
1993 if (numberActiveColumns_ * sizeof(double) > 1000000) {
1994 if (numberRows * 10 < numberActiveColumns_)
1995 factor *= 0.333333333;
1996 else if (numberRows * 4 < numberActiveColumns_)
1997 factor *= 0.5;
1998 else if (numberRows * 2 < numberActiveColumns_)
1999 factor *= 0.66666666667;
2000 //if (model->numberIterations()%50==0)
2001 //printf("%d nonzero\n",numberInRowArray);
2002 }
2003 // if not packed then bias a bit more towards by column
2004 if (!packed)
2005 factor *= 0.9;
2006 return ((numberInRowArray > factor * numberRows || !model->rowCopy()) && !(flags_ & 2));
2007 }
2008 #ifndef CLP_ALL_ONE_FILE
2009 // These have to match ClpPrimalColumnSteepest version
2010 #define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)
2011 #endif
2012 // Updates two arrays for steepest
2013 void
transposeTimes2(const ClpSimplex * model,const CoinIndexedVector * pi1,CoinIndexedVector * dj1,const CoinIndexedVector * pi2,CoinIndexedVector * spare,double referenceIn,double devex,unsigned int * reference,double * weights,double scaleFactor)2014 ClpPackedMatrix::transposeTimes2(const ClpSimplex * model,
2015 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
2016 const CoinIndexedVector * pi2,
2017 CoinIndexedVector * spare,
2018 double referenceIn, double devex,
2019 // Array for exact devex to say what is in reference framework
2020 unsigned int * reference,
2021 double * weights, double scaleFactor)
2022 {
2023 // put row of tableau in dj1
2024 double * pi = pi1->denseVector();
2025 int numberNonZero = 0;
2026 int * index = dj1->getIndices();
2027 double * array = dj1->denseVector();
2028 int numberInRowArray = pi1->getNumElements();
2029 double zeroTolerance = model->zeroTolerance();
2030 bool packed = pi1->packedMode();
2031 // do by column
2032 int iColumn;
2033 // get matrix data pointers
2034 const int * row = matrix_->getIndices();
2035 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2036 const double * elementByColumn = matrix_->getElements();
2037 const double * rowScale = model->rowScale();
2038 assert (!spare->getNumElements());
2039 assert (numberActiveColumns_ > 0);
2040 double * piWeight = pi2->denseVector();
2041 assert (!pi2->packedMode());
2042 bool killDjs = (scaleFactor == 0.0);
2043 if (!scaleFactor)
2044 scaleFactor = 1.0;
2045 if (packed) {
2046 // need to expand pi into y
2047 assert(spare->capacity() >= model->numberRows());
2048 double * piOld = pi;
2049 pi = spare->denseVector();
2050 const int * whichRow = pi1->getIndices();
2051 int i;
2052 ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
2053 if (rowScale && scaledMatrix) {
2054 rowScale = NULL;
2055 // get matrix data pointers
2056 row = scaledMatrix->getIndices();
2057 columnStart = scaledMatrix->getVectorStarts();
2058 elementByColumn = scaledMatrix->getElements();
2059 }
2060 if (!rowScale) {
2061 // modify pi so can collapse to one loop
2062 for (i = 0; i < numberInRowArray; i++) {
2063 int iRow = whichRow[i];
2064 pi[iRow] = piOld[i];
2065 }
2066 if (!columnCopy_) {
2067 CoinBigIndex j;
2068 CoinBigIndex end = columnStart[0];
2069 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
2070 CoinBigIndex start = end;
2071 end = columnStart[iColumn+1];
2072 ClpSimplex::Status status = model->getStatus(iColumn);
2073 if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
2074 double value = 0.0;
2075 for (j = start; j < end; j++) {
2076 int iRow = row[j];
2077 value -= pi[iRow] * elementByColumn[j];
2078 }
2079 if (fabs(value) > zeroTolerance) {
2080 // and do other array
2081 double modification = 0.0;
2082 for (j = start; j < end; j++) {
2083 int iRow = row[j];
2084 modification += piWeight[iRow] * elementByColumn[j];
2085 }
2086 double thisWeight = weights[iColumn];
2087 double pivot = value * scaleFactor;
2088 double pivotSquared = pivot * pivot;
2089 thisWeight += pivotSquared * devex + pivot * modification;
2090 if (thisWeight < DEVEX_TRY_NORM) {
2091 if (referenceIn < 0.0) {
2092 // steepest
2093 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2094 } else {
2095 // exact
2096 thisWeight = referenceIn * pivotSquared;
2097 if (reference(iColumn))
2098 thisWeight += 1.0;
2099 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2100 }
2101 }
2102 weights[iColumn] = thisWeight;
2103 if (!killDjs) {
2104 array[numberNonZero] = value;
2105 index[numberNonZero++] = iColumn;
2106 }
2107 }
2108 }
2109 } else {
2110 // use special column copy
2111 // reset back
2112 if (killDjs)
2113 scaleFactor = 0.0;
2114 columnCopy_->transposeTimes2(model, pi, dj1, piWeight, referenceIn, devex,
2115 reference, weights, scaleFactor);
2116 numberNonZero = dj1->getNumElements();
2117 }
2118 } else {
2119 // scaled
2120 // modify pi so can collapse to one loop
2121 for (i = 0; i < numberInRowArray; i++) {
2122 int iRow = whichRow[i];
2123 pi[iRow] = piOld[i] * rowScale[iRow];
2124 }
2125 // can also scale piWeight as not used again
2126 int numberWeight = pi2->getNumElements();
2127 const int * indexWeight = pi2->getIndices();
2128 for (i = 0; i < numberWeight; i++) {
2129 int iRow = indexWeight[i];
2130 piWeight[iRow] *= rowScale[iRow];
2131 }
2132 if (!columnCopy_) {
2133 const double * columnScale = model->columnScale();
2134 CoinBigIndex j;
2135 CoinBigIndex end = columnStart[0];
2136 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
2137 CoinBigIndex start = end;
2138 end = columnStart[iColumn+1];
2139 ClpSimplex::Status status = model->getStatus(iColumn);
2140 if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
2141 double scale = columnScale[iColumn];
2142 double value = 0.0;
2143 for (j = start; j < end; j++) {
2144 int iRow = row[j];
2145 value -= pi[iRow] * elementByColumn[j];
2146 }
2147 value *= scale;
2148 if (fabs(value) > zeroTolerance) {
2149 double modification = 0.0;
2150 for (j = start; j < end; j++) {
2151 int iRow = row[j];
2152 modification += piWeight[iRow] * elementByColumn[j];
2153 }
2154 modification *= scale;
2155 double thisWeight = weights[iColumn];
2156 double pivot = value * scaleFactor;
2157 double pivotSquared = pivot * pivot;
2158 thisWeight += pivotSquared * devex + pivot * modification;
2159 if (thisWeight < DEVEX_TRY_NORM) {
2160 if (referenceIn < 0.0) {
2161 // steepest
2162 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2163 } else {
2164 // exact
2165 thisWeight = referenceIn * pivotSquared;
2166 if (reference(iColumn))
2167 thisWeight += 1.0;
2168 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2169 }
2170 }
2171 weights[iColumn] = thisWeight;
2172 if (!killDjs) {
2173 array[numberNonZero] = value;
2174 index[numberNonZero++] = iColumn;
2175 }
2176 }
2177 }
2178 } else {
2179 // use special column copy
2180 // reset back
2181 if (killDjs)
2182 scaleFactor = 0.0;
2183 columnCopy_->transposeTimes2(model, pi, dj1, piWeight, referenceIn, devex,
2184 reference, weights, scaleFactor);
2185 numberNonZero = dj1->getNumElements();
2186 }
2187 }
2188 // zero out
2189 int numberRows = model->numberRows();
2190 if (numberInRowArray * 4 < numberRows) {
2191 for (i = 0; i < numberInRowArray; i++) {
2192 int iRow = whichRow[i];
2193 pi[iRow] = 0.0;
2194 }
2195 } else {
2196 CoinZeroN(pi, numberRows);
2197 }
2198 } else {
2199 if (!rowScale) {
2200 CoinBigIndex j;
2201 CoinBigIndex end = columnStart[0];
2202 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
2203 CoinBigIndex start = end;
2204 end = columnStart[iColumn+1];
2205 ClpSimplex::Status status = model->getStatus(iColumn);
2206 if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
2207 double value = 0.0;
2208 for (j = start; j < end; j++) {
2209 int iRow = row[j];
2210 value -= pi[iRow] * elementByColumn[j];
2211 }
2212 if (fabs(value) > zeroTolerance) {
2213 // and do other array
2214 double modification = 0.0;
2215 for (j = start; j < end; j++) {
2216 int iRow = row[j];
2217 modification += piWeight[iRow] * elementByColumn[j];
2218 }
2219 double thisWeight = weights[iColumn];
2220 double pivot = value * scaleFactor;
2221 double pivotSquared = pivot * pivot;
2222 thisWeight += pivotSquared * devex + pivot * modification;
2223 if (thisWeight < DEVEX_TRY_NORM) {
2224 if (referenceIn < 0.0) {
2225 // steepest
2226 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2227 } else {
2228 // exact
2229 thisWeight = referenceIn * pivotSquared;
2230 if (reference(iColumn))
2231 thisWeight += 1.0;
2232 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2233 }
2234 }
2235 weights[iColumn] = thisWeight;
2236 if (!killDjs) {
2237 array[iColumn] = value;
2238 index[numberNonZero++] = iColumn;
2239 }
2240 }
2241 }
2242 } else {
2243 #ifdef CLP_INVESTIGATE
2244 if (model->clpScaledMatrix())
2245 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
2246 #endif
2247 // scaled
2248 // can also scale piWeight as not used again
2249 int numberWeight = pi2->getNumElements();
2250 const int * indexWeight = pi2->getIndices();
2251 for (int i = 0; i < numberWeight; i++) {
2252 int iRow = indexWeight[i];
2253 piWeight[iRow] *= rowScale[iRow];
2254 }
2255 const double * columnScale = model->columnScale();
2256 CoinBigIndex j;
2257 CoinBigIndex end = columnStart[0];
2258 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
2259 CoinBigIndex start = end;
2260 end = columnStart[iColumn+1];
2261 ClpSimplex::Status status = model->getStatus(iColumn);
2262 if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
2263 double scale = columnScale[iColumn];
2264 double value = 0.0;
2265 for (j = start; j < end; j++) {
2266 int iRow = row[j];
2267 value -= pi[iRow] * elementByColumn[j] * rowScale[iRow];
2268 }
2269 value *= scale;
2270 if (fabs(value) > zeroTolerance) {
2271 double modification = 0.0;
2272 for (j = start; j < end; j++) {
2273 int iRow = row[j];
2274 modification += piWeight[iRow] * elementByColumn[j];
2275 }
2276 modification *= scale;
2277 double thisWeight = weights[iColumn];
2278 double pivot = value * scaleFactor;
2279 double pivotSquared = pivot * pivot;
2280 thisWeight += pivotSquared * devex + pivot * modification;
2281 if (thisWeight < DEVEX_TRY_NORM) {
2282 if (referenceIn < 0.0) {
2283 // steepest
2284 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2285 } else {
2286 // exact
2287 thisWeight = referenceIn * pivotSquared;
2288 if (reference(iColumn))
2289 thisWeight += 1.0;
2290 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2291 }
2292 }
2293 weights[iColumn] = thisWeight;
2294 if (!killDjs) {
2295 array[iColumn] = value;
2296 index[numberNonZero++] = iColumn;
2297 }
2298 }
2299 }
2300 }
2301 }
2302 dj1->setNumElements(numberNonZero);
2303 spare->setNumElements(0);
2304 if (packed)
2305 dj1->setPackedMode(true);
2306 }
2307 // Updates second array for steepest and does devex weights
2308 void
subsetTimes2(const ClpSimplex * model,CoinIndexedVector * dj1,const CoinIndexedVector * pi2,CoinIndexedVector *,double referenceIn,double devex,unsigned int * reference,double * weights,double scaleFactor)2309 ClpPackedMatrix::subsetTimes2(const ClpSimplex * model,
2310 CoinIndexedVector * dj1,
2311 const CoinIndexedVector * pi2, CoinIndexedVector *,
2312 double referenceIn, double devex,
2313 // Array for exact devex to say what is in reference framework
2314 unsigned int * reference,
2315 double * weights, double scaleFactor)
2316 {
2317 int number = dj1->getNumElements();
2318 const int * index = dj1->getIndices();
2319 double * array = dj1->denseVector();
2320 assert( dj1->packedMode());
2321
2322 // get matrix data pointers
2323 const int * row = matrix_->getIndices();
2324 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2325 const int * columnLength = matrix_->getVectorLengths();
2326 const double * elementByColumn = matrix_->getElements();
2327 const double * rowScale = model->rowScale();
2328 double * piWeight = pi2->denseVector();
2329 bool killDjs = (scaleFactor == 0.0);
2330 if (!scaleFactor)
2331 scaleFactor = 1.0;
2332 if (!rowScale) {
2333 for (int k = 0; k < number; k++) {
2334 int iColumn = index[k];
2335 double pivot = array[k] * scaleFactor;
2336 if (killDjs)
2337 array[k] = 0.0;
2338 // and do other array
2339 double modification = 0.0;
2340 for (CoinBigIndex j = columnStart[iColumn];
2341 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2342 int iRow = row[j];
2343 modification += piWeight[iRow] * elementByColumn[j];
2344 }
2345 double thisWeight = weights[iColumn];
2346 double pivotSquared = pivot * pivot;
2347 thisWeight += pivotSquared * devex + pivot * modification;
2348 if (thisWeight < DEVEX_TRY_NORM) {
2349 if (referenceIn < 0.0) {
2350 // steepest
2351 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2352 } else {
2353 // exact
2354 thisWeight = referenceIn * pivotSquared;
2355 if (reference(iColumn))
2356 thisWeight += 1.0;
2357 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2358 }
2359 }
2360 weights[iColumn] = thisWeight;
2361 }
2362 } else {
2363 #ifdef CLP_INVESTIGATE
2364 if (model->clpScaledMatrix())
2365 printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
2366 #endif
2367 // scaled
2368 const double * columnScale = model->columnScale();
2369 for (int k = 0; k < number; k++) {
2370 int iColumn = index[k];
2371 double pivot = array[k] * scaleFactor;
2372 double scale = columnScale[iColumn];
2373 if (killDjs)
2374 array[k] = 0.0;
2375 // and do other array
2376 double modification = 0.0;
2377 for (CoinBigIndex j = columnStart[iColumn];
2378 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2379 int iRow = row[j];
2380 modification += piWeight[iRow] * elementByColumn[j] * rowScale[iRow];
2381 }
2382 double thisWeight = weights[iColumn];
2383 modification *= scale;
2384 double pivotSquared = pivot * pivot;
2385 thisWeight += pivotSquared * devex + pivot * modification;
2386 if (thisWeight < DEVEX_TRY_NORM) {
2387 if (referenceIn < 0.0) {
2388 // steepest
2389 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2390 } else {
2391 // exact
2392 thisWeight = referenceIn * pivotSquared;
2393 if (reference(iColumn))
2394 thisWeight += 1.0;
2395 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2396 }
2397 }
2398 weights[iColumn] = thisWeight;
2399 }
2400 }
2401 }
2402 /// returns number of elements in column part of basis,
2403 CoinBigIndex
countBasis(const int * whichColumn,int & numberColumnBasic)2404 ClpPackedMatrix::countBasis( const int * whichColumn,
2405 int & numberColumnBasic)
2406 {
2407 const int * columnLength = matrix_->getVectorLengths();
2408 int i;
2409 CoinBigIndex numberElements = 0;
2410 // just count - can be over so ignore zero problem
2411 for (i = 0; i < numberColumnBasic; i++) {
2412 int iColumn = whichColumn[i];
2413 numberElements += columnLength[iColumn];
2414 }
2415 return numberElements;
2416 }
2417 void
fillBasis(ClpSimplex * model,const int * COIN_RESTRICT whichColumn,int & numberColumnBasic,int * COIN_RESTRICT indexRowU,int * COIN_RESTRICT start,int * COIN_RESTRICT rowCount,int * COIN_RESTRICT columnCount,CoinFactorizationDouble * COIN_RESTRICT elementU)2418 ClpPackedMatrix::fillBasis(ClpSimplex * model,
2419 const int * COIN_RESTRICT whichColumn,
2420 int & numberColumnBasic,
2421 int * COIN_RESTRICT indexRowU,
2422 int * COIN_RESTRICT start,
2423 int * COIN_RESTRICT rowCount,
2424 int * COIN_RESTRICT columnCount,
2425 CoinFactorizationDouble * COIN_RESTRICT elementU)
2426 {
2427 const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
2428 int i;
2429 CoinBigIndex numberElements = start[0];
2430 // fill
2431 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
2432 const double * COIN_RESTRICT rowScale = model->rowScale();
2433 const int * COIN_RESTRICT row = matrix_->getIndices();
2434 const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2435 ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
2436 if (scaledMatrix && true) {
2437 columnLength = scaledMatrix->matrix_->getVectorLengths();
2438 columnStart = scaledMatrix->matrix_->getVectorStarts();
2439 rowScale = NULL;
2440 row = scaledMatrix->matrix_->getIndices();
2441 elementByColumn = scaledMatrix->matrix_->getElements();
2442 }
2443 if ((flags_ & 1) == 0) {
2444 if (!rowScale) {
2445 // no scaling
2446 for (i = 0; i < numberColumnBasic; i++) {
2447 int iColumn = whichColumn[i];
2448 int length = columnLength[iColumn];
2449 CoinBigIndex startThis = columnStart[iColumn];
2450 columnCount[i] = length;
2451 CoinBigIndex endThis = startThis + length;
2452 for (CoinBigIndex j = startThis; j < endThis; j++) {
2453 int iRow = row[j];
2454 indexRowU[numberElements] = iRow;
2455 rowCount[iRow]++;
2456 assert (elementByColumn[j]);
2457 elementU[numberElements++] = elementByColumn[j];
2458 }
2459 start[i+1] = numberElements;
2460 }
2461 } else {
2462 // scaling
2463 const double * COIN_RESTRICT columnScale = model->columnScale();
2464 for (i = 0; i < numberColumnBasic; i++) {
2465 int iColumn = whichColumn[i];
2466 double scale = columnScale[iColumn];
2467 int length = columnLength[iColumn];
2468 CoinBigIndex startThis = columnStart[iColumn];
2469 columnCount[i] = length;
2470 CoinBigIndex endThis = startThis + length;
2471 for (CoinBigIndex j = startThis; j < endThis; j++) {
2472 int iRow = row[j];
2473 indexRowU[numberElements] = iRow;
2474 rowCount[iRow]++;
2475 assert (elementByColumn[j]);
2476 elementU[numberElements++] =
2477 elementByColumn[j] * scale * rowScale[iRow];
2478 }
2479 start[i+1] = numberElements;
2480 }
2481 }
2482 } else {
2483 // there are zero elements so need to look more closely
2484 if (!rowScale) {
2485 // no scaling
2486 for (i = 0; i < numberColumnBasic; i++) {
2487 int iColumn = whichColumn[i];
2488 CoinBigIndex j;
2489 for (j = columnStart[iColumn];
2490 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2491 double value = elementByColumn[j];
2492 if (value) {
2493 int iRow = row[j];
2494 indexRowU[numberElements] = iRow;
2495 rowCount[iRow]++;
2496 elementU[numberElements++] = value;
2497 }
2498 }
2499 start[i+1] = numberElements;
2500 columnCount[i] = numberElements - start[i];
2501 }
2502 } else {
2503 // scaling
2504 const double * columnScale = model->columnScale();
2505 for (i = 0; i < numberColumnBasic; i++) {
2506 int iColumn = whichColumn[i];
2507 CoinBigIndex j;
2508 double scale = columnScale[iColumn];
2509 for (j = columnStart[iColumn];
2510 j < columnStart[iColumn] + columnLength[i]; j++) {
2511 double value = elementByColumn[j];
2512 if (value) {
2513 int iRow = row[j];
2514 indexRowU[numberElements] = iRow;
2515 rowCount[iRow]++;
2516 elementU[numberElements++] = value * scale * rowScale[iRow];
2517 }
2518 }
2519 start[i+1] = numberElements;
2520 columnCount[i] = numberElements - start[i];
2521 }
2522 }
2523 }
2524 }
2525 #if 0
2526 int
2527 ClpPackedMatrix::scale2(ClpModel * model) const
2528 {
2529 ClpSimplex * baseModel = NULL;
2530 #ifndef NDEBUG
2531 //checkFlags();
2532 #endif
2533 int numberRows = model->numberRows();
2534 int numberColumns = matrix_->getNumCols();
2535 model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
2536 // If empty - return as sanityCheck will trap
2537 if (!numberRows || !numberColumns) {
2538 model->setRowScale(NULL);
2539 model->setColumnScale(NULL);
2540 return 1;
2541 }
2542 ClpMatrixBase * rowCopyBase = model->rowCopy();
2543 double * rowScale;
2544 double * columnScale;
2545 //assert (!model->rowScale());
2546 bool arraysExist;
2547 double * inverseRowScale = NULL;
2548 double * inverseColumnScale = NULL;
2549 if (!model->rowScale()) {
2550 rowScale = new double [numberRows*2];
2551 columnScale = new double [numberColumns*2];
2552 inverseRowScale = rowScale + numberRows;
2553 inverseColumnScale = columnScale + numberColumns;
2554 arraysExist = false;
2555 } else {
2556 rowScale = model->mutableRowScale();
2557 columnScale = model->mutableColumnScale();
2558 inverseRowScale = model->mutableInverseRowScale();
2559 inverseColumnScale = model->mutableInverseColumnScale();
2560 arraysExist = true;
2561 }
2562 assert (inverseRowScale == rowScale + numberRows);
2563 assert (inverseColumnScale == columnScale + numberColumns);
2564 // we are going to mark bits we are interested in
2565 char * usefulRow = new char [numberRows];
2566 char * usefulColumn = new char [numberColumns];
2567 double * rowLower = model->rowLower();
2568 double * rowUpper = model->rowUpper();
2569 double * columnLower = model->columnLower();
2570 double * columnUpper = model->columnUpper();
2571 int iColumn, iRow;
2572 //#define LEAVE_FIXED
2573 // mark free rows
2574 for (iRow = 0; iRow < numberRows; iRow++) {
2575 #if 0 //ndef LEAVE_FIXED
2576 if (rowUpper[iRow] < 1.0e20 ||
2577 rowLower[iRow] > -1.0e20)
2578 usefulRow[iRow] = 1;
2579 else
2580 usefulRow[iRow] = 0;
2581 #else
2582 usefulRow[iRow] = 1;
2583 #endif
2584 }
2585 // mark empty and fixed columns
2586 // also see if worth scaling
2587 assert (model->scalingFlag() <= 4);
2588 // scale_stats[model->scalingFlag()]++;
2589 double largest = 0.0;
2590 double smallest = 1.0e50;
2591 // get matrix data pointers
2592 int * row = matrix_->getMutableIndices();
2593 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2594 int * columnLength = matrix_->getMutableVectorLengths();
2595 double * elementByColumn = matrix_->getMutableElements();
2596 bool deletedElements = false;
2597 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2598 CoinBigIndex j;
2599 char useful = 0;
2600 bool deleteSome = false;
2601 int start = columnStart[iColumn];
2602 int end = start + columnLength[iColumn];
2603 #ifndef LEAVE_FIXED
2604 if (columnUpper[iColumn] >
2605 columnLower[iColumn] + 1.0e-12) {
2606 #endif
2607 for (j = start; j < end; j++) {
2608 iRow = row[j];
2609 double value = fabs(elementByColumn[j]);
2610 if (value > 1.0e-20) {
2611 if(usefulRow[iRow]) {
2612 useful = 1;
2613 largest = CoinMax(largest, fabs(elementByColumn[j]));
2614 smallest = CoinMin(smallest, fabs(elementByColumn[j]));
2615 }
2616 } else {
2617 // small
2618 deleteSome = true;
2619 }
2620 }
2621 #ifndef LEAVE_FIXED
2622 } else {
2623 // just check values
2624 for (j = start; j < end; j++) {
2625 double value = fabs(elementByColumn[j]);
2626 if (value <= 1.0e-20) {
2627 // small
2628 deleteSome = true;
2629 }
2630 }
2631 }
2632 #endif
2633 usefulColumn[iColumn] = useful;
2634 if (deleteSome) {
2635 deletedElements = true;
2636 CoinBigIndex put = start;
2637 for (j = start; j < end; j++) {
2638 double value = elementByColumn[j];
2639 if (fabs(value) > 1.0e-20) {
2640 row[put] = row[j];
2641 elementByColumn[put++] = value;
2642 }
2643 }
2644 columnLength[iColumn] = put - start;
2645 }
2646 }
2647 model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL, *model->messagesPointer())
2648 << smallest << largest
2649 << CoinMessageEol;
2650 if (smallest >= 0.5 && largest <= 2.0 && !deletedElements) {
2651 // don't bother scaling
2652 model->messageHandler()->message(CLP_PACKEDSCALE_FORGET, *model->messagesPointer())
2653 << CoinMessageEol;
2654 if (!arraysExist) {
2655 delete [] rowScale;
2656 delete [] columnScale;
2657 } else {
2658 model->setRowScale(NULL);
2659 model->setColumnScale(NULL);
2660 }
2661 delete [] usefulRow;
2662 delete [] usefulColumn;
2663 return 1;
2664 } else {
2665 #ifdef CLP_INVESTIGATE
2666 if (deletedElements)
2667 printf("DEL_ELS\n");
2668 #endif
2669 if (!rowCopyBase) {
2670 // temporary copy
2671 rowCopyBase = reverseOrderedCopy();
2672 } else if (deletedElements) {
2673 rowCopyBase = reverseOrderedCopy();
2674 model->setNewRowCopy(rowCopyBase);
2675 }
2676 #ifndef NDEBUG
2677 ClpPackedMatrix* rowCopy =
2678 dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
2679 // Make sure it is really a ClpPackedMatrix
2680 assert (rowCopy != NULL);
2681 #else
2682 ClpPackedMatrix* rowCopy =
2683 static_cast< ClpPackedMatrix*>(rowCopyBase);
2684 #endif
2685
2686 const int * column = rowCopy->getIndices();
2687 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2688 const double * element = rowCopy->getElements();
2689 // need to scale
2690 if (largest > 1.0e13 * smallest) {
2691 // safer to have smaller zero tolerance
2692 double ratio = smallest / largest;
2693 ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
2694 double newTolerance = CoinMax(ratio * 0.5, 1.0e-18);
2695 if (simplex->zeroTolerance() > newTolerance)
2696 simplex->setZeroTolerance(newTolerance);
2697 }
2698 int scalingMethod = model->scalingFlag();
2699 if (scalingMethod == 4) {
2700 // As auto
2701 scalingMethod = 3;
2702 }
2703 // and see if there any empty rows
2704 for (iRow = 0; iRow < numberRows; iRow++) {
2705 if (usefulRow[iRow]) {
2706 CoinBigIndex j;
2707 int useful = 0;
2708 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
2709 int iColumn = column[j];
2710 if (usefulColumn[iColumn]) {
2711 useful = 1;
2712 break;
2713 }
2714 }
2715 usefulRow[iRow] = static_cast<char>(useful);
2716 }
2717 }
2718 double savedOverallRatio = 0.0;
2719 double tolerance = 5.0 * model->primalTolerance();
2720 double overallLargest = -1.0e-20;
2721 double overallSmallest = 1.0e20;
2722 bool finished = false;
2723 // if scalingMethod 3 then may change
2724 bool extraDetails = (model->logLevel() > 2);
2725 while (!finished) {
2726 int numberPass = 3;
2727 overallLargest = -1.0e-20;
2728 overallSmallest = 1.0e20;
2729 if (!baseModel) {
2730 ClpFillN ( rowScale, numberRows, 1.0);
2731 ClpFillN ( columnScale, numberColumns, 1.0);
2732 } else {
2733 // Copy scales and do quick scale on extra rows
2734 // Then just one? pass
2735 assert (numberColumns == baseModel->numberColumns());
2736 int numberRows2 = baseModel->numberRows();
2737 assert (numberRows >= numberRows2);
2738 assert (baseModel->rowScale());
2739 CoinMemcpyN(baseModel->rowScale(), numberRows2, rowScale);
2740 CoinMemcpyN(baseModel->columnScale(), numberColumns, columnScale);
2741 if (numberRows > numberRows2) {
2742 numberPass = 1;
2743 // do some scaling
2744 if (scalingMethod == 1 || scalingMethod == 3) {
2745 // Maximum in each row
2746 for (iRow = numberRows2; iRow < numberRows; iRow++) {
2747 if (usefulRow[iRow]) {
2748 CoinBigIndex j;
2749 largest = 1.0e-10;
2750 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
2751 int iColumn = column[j];
2752 if (usefulColumn[iColumn]) {
2753 double value = fabs(element[j] * columnScale[iColumn]);
2754 largest = CoinMax(largest, value);
2755 assert (largest < 1.0e40);
2756 }
2757 }
2758 rowScale[iRow] = 1.0 / largest;
2759 #ifdef COIN_DEVELOP
2760 if (extraDetails) {
2761 overallLargest = CoinMax(overallLargest, largest);
2762 overallSmallest = CoinMin(overallSmallest, largest);
2763 }
2764 #endif
2765 }
2766 }
2767 } else {
2768 overallLargest = 0.0;
2769 overallSmallest = 1.0e50;
2770 // Geometric mean on row scales
2771 for (iRow = 0; iRow < numberRows; iRow++) {
2772 if (usefulRow[iRow]) {
2773 CoinBigIndex j;
2774 largest = 1.0e-20;
2775 smallest = 1.0e50;
2776 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
2777 int iColumn = column[j];
2778 if (usefulColumn[iColumn]) {
2779 double value = fabs(element[j]);
2780 value *= columnScale[iColumn];
2781 largest = CoinMax(largest, value);
2782 smallest = CoinMin(smallest, value);
2783 }
2784 }
2785 if (iRow >= numberRows2) {
2786 rowScale[iRow] = 1.0 / sqrt(smallest * largest);
2787 //rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
2788 }
2789 #ifdef COIN_DEVELOP
2790 if (extraDetails) {
2791 overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
2792 overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
2793 }
2794 #endif
2795 }
2796 }
2797 }
2798 } else {
2799 // just use
2800 numberPass = 0;
2801 }
2802 }
2803 if (!baseModel && (scalingMethod == 1 || scalingMethod == 3)) {
2804 // Maximum in each row
2805 for (iRow = 0; iRow < numberRows; iRow++) {
2806 if (usefulRow[iRow]) {
2807 CoinBigIndex j;
2808 largest = 1.0e-10;
2809 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
2810 int iColumn = column[j];
2811 if (usefulColumn[iColumn]) {
2812 double value = fabs(element[j]);
2813 largest = CoinMax(largest, value);
2814 assert (largest < 1.0e40);
2815 }
2816 }
2817 rowScale[iRow] = 1.0 / largest;
2818 #ifdef COIN_DEVELOP
2819 if (extraDetails) {
2820 overallLargest = CoinMax(overallLargest, largest);
2821 overallSmallest = CoinMin(overallSmallest, largest);
2822 }
2823 #endif
2824 }
2825 }
2826 } else {
2827 #ifdef USE_OBJECTIVE
2828 // This will be used to help get scale factors
2829 double * objective = new double[numberColumns];
2830 CoinMemcpyN(model->costRegion(1), numberColumns, objective);
2831 double objScale = 1.0;
2832 #endif
2833 while (numberPass) {
2834 overallLargest = 0.0;
2835 overallSmallest = 1.0e50;
2836 numberPass--;
2837 // Geometric mean on row scales
2838 for (iRow = 0; iRow < numberRows; iRow++) {
2839 if (usefulRow[iRow]) {
2840 CoinBigIndex j;
2841 largest = 1.0e-20;
2842 smallest = 1.0e50;
2843 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
2844 int iColumn = column[j];
2845 if (usefulColumn[iColumn]) {
2846 double value = fabs(element[j]);
2847 value *= columnScale[iColumn];
2848 largest = CoinMax(largest, value);
2849 smallest = CoinMin(smallest, value);
2850 }
2851 }
2852
2853 rowScale[iRow] = 1.0 / sqrt(smallest * largest);
2854 //rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
2855 if (extraDetails) {
2856 overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
2857 overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
2858 }
2859 }
2860 }
2861 #ifdef USE_OBJECTIVE
2862 largest = 1.0e-20;
2863 smallest = 1.0e50;
2864 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2865 if (usefulColumn[iColumn]) {
2866 double value = fabs(objective[iColumn]);
2867 value *= columnScale[iColumn];
2868 largest = CoinMax(largest, value);
2869 smallest = CoinMin(smallest, value);
2870 }
2871 }
2872 objScale = 1.0 / sqrt(smallest * largest);
2873 #endif
2874 model->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model->messagesPointer())
2875 << overallSmallest
2876 << overallLargest
2877 << CoinMessageEol;
2878 // skip last column round
2879 if (numberPass == 1)
2880 break;
2881 // Geometric mean on column scales
2882 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2883 if (usefulColumn[iColumn]) {
2884 CoinBigIndex j;
2885 largest = 1.0e-20;
2886 smallest = 1.0e50;
2887 for (j = columnStart[iColumn];
2888 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2889 iRow = row[j];
2890 double value = fabs(elementByColumn[j]);
2891 if (usefulRow[iRow]) {
2892 value *= rowScale[iRow];
2893 largest = CoinMax(largest, value);
2894 smallest = CoinMin(smallest, value);
2895 }
2896 }
2897 #ifdef USE_OBJECTIVE
2898 if (fabs(objective[iColumn]) > 1.0e-20) {
2899 double value = fabs(objective[iColumn]) * objScale;
2900 largest = CoinMax(largest, value);
2901 smallest = CoinMin(smallest, value);
2902 }
2903 #endif
2904 columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
2905 //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
2906 }
2907 }
2908 }
2909 #ifdef USE_OBJECTIVE
2910 delete [] objective;
2911 printf("obj scale %g - use it if you want\n", objScale);
2912 #endif
2913 }
2914 // If ranges will make horrid then scale
2915 for (iRow = 0; iRow < numberRows; iRow++) {
2916 if (usefulRow[iRow]) {
2917 double difference = rowUpper[iRow] - rowLower[iRow];
2918 double scaledDifference = difference * rowScale[iRow];
2919 if (scaledDifference > tolerance && scaledDifference < 1.0e-4) {
2920 // make gap larger
2921 rowScale[iRow] *= 1.0e-4 / scaledDifference;
2922 rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
2923 //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
2924 // scaledDifference,difference*rowScale[iRow]);
2925 }
2926 }
2927 }
2928 // final pass to scale columns so largest is reasonable
2929 // See what smallest will be if largest is 1.0
2930 overallSmallest = 1.0e50;
2931 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2932 if (usefulColumn[iColumn]) {
2933 CoinBigIndex j;
2934 largest = 1.0e-20;
2935 smallest = 1.0e50;
2936 for (j = columnStart[iColumn];
2937 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2938 iRow = row[j];
2939 if(elementByColumn[j] && usefulRow[iRow]) {
2940 double value = fabs(elementByColumn[j] * rowScale[iRow]);
2941 largest = CoinMax(largest, value);
2942 smallest = CoinMin(smallest, value);
2943 }
2944 }
2945 if (overallSmallest * largest > smallest)
2946 overallSmallest = smallest / largest;
2947 }
2948 }
2949 if (scalingMethod == 1 || scalingMethod == 2) {
2950 finished = true;
2951 } else if (savedOverallRatio == 0.0 && scalingMethod != 4) {
2952 savedOverallRatio = overallSmallest;
2953 scalingMethod = 4;
2954 } else {
2955 assert (scalingMethod == 4);
2956 if (overallSmallest > 2.0 * savedOverallRatio) {
2957 finished = true; // geometric was better
2958 if (model->scalingFlag() == 4) {
2959 // If in Branch and bound change
2960 if ((model->specialOptions() & 1024) != 0) {
2961 model->scaling(2);
2962 }
2963 }
2964 } else {
2965 scalingMethod = 1; // redo equilibrium
2966 if (model->scalingFlag() == 4) {
2967 // If in Branch and bound change
2968 if ((model->specialOptions() & 1024) != 0) {
2969 model->scaling(1);
2970 }
2971 }
2972 }
2973 #if 0
2974 if (extraDetails) {
2975 if (finished)
2976 printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
2977 savedOverallRatio, overallSmallest);
2978 else
2979 printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
2980 savedOverallRatio, overallSmallest);
2981 }
2982 #endif
2983 }
2984 }
2985 //#define RANDOMIZE
2986 #ifdef RANDOMIZE
2987 // randomize by up to 10%
2988 for (iRow = 0; iRow < numberRows; iRow++) {
2989 double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
2990 rowScale[iRow] *= (1.0 + 0.1 * value);
2991 }
2992 #endif
2993 overallLargest = 1.0;
2994 if (overallSmallest < 1.0e-1)
2995 overallLargest = 1.0 / sqrt(overallSmallest);
2996 overallLargest = CoinMin(100.0, overallLargest);
2997 overallSmallest = 1.0e50;
2998 //printf("scaling %d\n",model->scalingFlag());
2999 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3000 if (columnUpper[iColumn] >
3001 columnLower[iColumn] + 1.0e-12) {
3002 //if (usefulColumn[iColumn]) {
3003 CoinBigIndex j;
3004 largest = 1.0e-20;
3005 smallest = 1.0e50;
3006 for (j = columnStart[iColumn];
3007 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3008 iRow = row[j];
3009 if(elementByColumn[j] && usefulRow[iRow]) {
3010 double value = fabs(elementByColumn[j] * rowScale[iRow]);
3011 largest = CoinMax(largest, value);
3012 smallest = CoinMin(smallest, value);
3013 }
3014 }
3015 columnScale[iColumn] = overallLargest / largest;
3016 //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3017 #ifdef RANDOMIZE
3018 double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3019 columnScale[iColumn] *= (1.0 + 0.1 * value);
3020 #endif
3021 double difference = columnUpper[iColumn] - columnLower[iColumn];
3022 if (difference < 1.0e-5 * columnScale[iColumn]) {
3023 // make gap larger
3024 columnScale[iColumn] = difference / 1.0e-5;
3025 //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
3026 // scaledDifference,difference*columnScale[iColumn]);
3027 }
3028 double value = smallest * columnScale[iColumn];
3029 if (overallSmallest > value)
3030 overallSmallest = value;
3031 //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
3032 }
3033 }
3034 model->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model->messagesPointer())
3035 << overallSmallest
3036 << overallLargest
3037 << CoinMessageEol;
3038 if (overallSmallest < 1.0e-13) {
3039 // Change factorization zero tolerance
3040 double newTolerance = CoinMax(1.0e-15 * (overallSmallest / 1.0e-13),
3041 1.0e-18);
3042 ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3043 if (simplex->factorization()->zeroTolerance() > newTolerance)
3044 simplex->factorization()->zeroTolerance(newTolerance);
3045 newTolerance = CoinMax(overallSmallest * 0.5, 1.0e-18);
3046 simplex->setZeroTolerance(newTolerance);
3047 }
3048 delete [] usefulRow;
3049 delete [] usefulColumn;
3050 #ifndef SLIM_CLP
3051 // If quadratic then make symmetric
3052 ClpObjective * obj = model->objectiveAsObject();
3053 #ifndef NO_RTTI
3054 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
3055 #else
3056 ClpQuadraticObjective * quadraticObj = NULL;
3057 if (obj->type() == 2)
3058 quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
3059 #endif
3060 if (quadraticObj) {
3061 if (!rowCopyBase) {
3062 // temporary copy
3063 rowCopyBase = reverseOrderedCopy();
3064 }
3065 #ifndef NDEBUG
3066 ClpPackedMatrix* rowCopy =
3067 dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3068 // Make sure it is really a ClpPackedMatrix
3069 assert (rowCopy != NULL);
3070 #else
3071 ClpPackedMatrix* rowCopy =
3072 static_cast< ClpPackedMatrix*>(rowCopyBase);
3073 #endif
3074 const int * column = rowCopy->getIndices();
3075 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3076 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
3077 int numberXColumns = quadratic->getNumCols();
3078 if (numberXColumns < numberColumns) {
3079 // we assume symmetric
3080 int numberQuadraticColumns = 0;
3081 int i;
3082 //const int * columnQuadratic = quadratic->getIndices();
3083 //const int * columnQuadraticStart = quadratic->getVectorStarts();
3084 const int * columnQuadraticLength = quadratic->getVectorLengths();
3085 for (i = 0; i < numberXColumns; i++) {
3086 int length = columnQuadraticLength[i];
3087 #ifndef CORRECT_COLUMN_COUNTS
3088 length = 1;
3089 #endif
3090 if (length)
3091 numberQuadraticColumns++;
3092 }
3093 int numberXRows = numberRows - numberQuadraticColumns;
3094 numberQuadraticColumns = 0;
3095 for (i = 0; i < numberXColumns; i++) {
3096 int length = columnQuadraticLength[i];
3097 #ifndef CORRECT_COLUMN_COUNTS
3098 length = 1;
3099 #endif
3100 if (length) {
3101 rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
3102 numberQuadraticColumns++;
3103 }
3104 }
3105 int numberQuadraticRows = 0;
3106 for (i = 0; i < numberXRows; i++) {
3107 // See if any in row quadratic
3108 CoinBigIndex j;
3109 int numberQ = 0;
3110 for (j = rowStart[i]; j < rowStart[i+1]; j++) {
3111 int iColumn = column[j];
3112 if (columnQuadraticLength[iColumn])
3113 numberQ++;
3114 }
3115 #ifndef CORRECT_ROW_COUNTS
3116 numberQ = 1;
3117 #endif
3118 if (numberQ) {
3119 columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
3120 numberQuadraticRows++;
3121 }
3122 }
3123 // and make sure Sj okay
3124 for (iColumn = numberQuadraticRows + numberXColumns; iColumn < numberColumns; iColumn++) {
3125 CoinBigIndex j = columnStart[iColumn];
3126 assert(columnLength[iColumn] == 1);
3127 int iRow = row[j];
3128 double value = fabs(elementByColumn[j] * rowScale[iRow]);
3129 columnScale[iColumn] = 1.0 / value;
3130 }
3131 }
3132 }
3133 #endif
3134 // make copy (could do faster by using previous values)
3135 // could just do partial
3136 for (iRow = 0; iRow < numberRows; iRow++)
3137 inverseRowScale[iRow] = 1.0 / rowScale[iRow] ;
3138 for (iColumn = 0; iColumn < numberColumns; iColumn++)
3139 inverseColumnScale[iColumn] = 1.0 / columnScale[iColumn] ;
3140 if (!arraysExist) {
3141 model->setRowScale(rowScale);
3142 model->setColumnScale(columnScale);
3143 }
3144 if (model->rowCopy()) {
3145 // need to replace row by row
3146 ClpPackedMatrix* rowCopy =
3147 static_cast< ClpPackedMatrix*>(model->rowCopy());
3148 double * element = rowCopy->getMutableElements();
3149 const int * column = rowCopy->getIndices();
3150 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3151 // scale row copy
3152 for (iRow = 0; iRow < numberRows; iRow++) {
3153 CoinBigIndex j;
3154 double scale = rowScale[iRow];
3155 double * elementsInThisRow = element + rowStart[iRow];
3156 const int * columnsInThisRow = column + rowStart[iRow];
3157 int number = rowStart[iRow+1] - rowStart[iRow];
3158 assert (number <= numberColumns);
3159 for (j = 0; j < number; j++) {
3160 int iColumn = columnsInThisRow[j];
3161 elementsInThisRow[j] *= scale * columnScale[iColumn];
3162 }
3163 }
3164 if ((model->specialOptions() & 262144) != 0) {
3165 //if ((model->specialOptions()&(COIN_CBC_USING_CLP|16384))!=0) {
3166 //if (model->inCbcBranchAndBound()&&false) {
3167 // copy without gaps
3168 CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, 0, 0);
3169 ClpPackedMatrix * scaled = new ClpPackedMatrix(scaledMatrix);
3170 model->setClpScaledMatrix(scaled);
3171 // get matrix data pointers
3172 const int * row = scaledMatrix->getIndices();
3173 const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3174 #ifndef NDEBUG
3175 const int * columnLength = scaledMatrix->getVectorLengths();
3176 #endif
3177 double * elementByColumn = scaledMatrix->getMutableElements();
3178 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3179 CoinBigIndex j;
3180 double scale = columnScale[iColumn];
3181 assert (columnStart[iColumn+1] == columnStart[iColumn] + columnLength[iColumn]);
3182 for (j = columnStart[iColumn];
3183 j < columnStart[iColumn+1]; j++) {
3184 int iRow = row[j];
3185 elementByColumn[j] *= scale * rowScale[iRow];
3186 }
3187 }
3188 } else {
3189 //printf("not in b&b\n");
3190 }
3191 } else {
3192 // no row copy
3193 delete rowCopyBase;
3194 }
3195 return 0;
3196 }
3197 }
3198 #endif
3199 //#define SQRT_ARRAY
3200 #ifdef SQRT_ARRAY
doSqrts(double * array,int n)3201 static void doSqrts(double * array, int n)
3202 {
3203 for (int i = 0; i < n; i++)
3204 array[i] = 1.0 / sqrt(array[i]);
3205 }
3206 #endif
3207 //static int scale_stats[5]={0,0,0,0,0};
3208 // Creates scales for column copy (rowCopy in model may be modified)
3209 int
scale(ClpModel * model,const ClpSimplex *) const3210 ClpPackedMatrix::scale(ClpModel * model, const ClpSimplex * /*baseModel*/) const
3211 {
3212 //const ClpSimplex * baseModel=NULL;
3213 //return scale2(model);
3214 #if 0
3215 ClpMatrixBase * rowClone = NULL;
3216 if (model->rowCopy())
3217 rowClone = model->rowCopy()->clone();
3218 assert (!model->rowScale());
3219 assert (!model->columnScale());
3220 int returnCode = scale2(model);
3221 if (returnCode)
3222 return returnCode;
3223 #endif
3224 #ifndef NDEBUG
3225 //checkFlags();
3226 #endif
3227 int numberRows = model->numberRows();
3228 int numberColumns = matrix_->getNumCols();
3229 model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
3230 // If empty - return as sanityCheck will trap
3231 if (!numberRows || !numberColumns) {
3232 model->setRowScale(NULL);
3233 model->setColumnScale(NULL);
3234 return 1;
3235 }
3236 #if 0
3237 // start fake
3238 double * rowScale2 = CoinCopyOfArray(model->rowScale(), numberRows);
3239 double * columnScale2 = CoinCopyOfArray(model->columnScale(), numberColumns);
3240 model->setRowScale(NULL);
3241 model->setColumnScale(NULL);
3242 model->setNewRowCopy(rowClone);
3243 #endif
3244 ClpMatrixBase * rowCopyBase = model->rowCopy();
3245 double * rowScale;
3246 double * columnScale;
3247 //assert (!model->rowScale());
3248 bool arraysExist;
3249 double * inverseRowScale = NULL;
3250 double * inverseColumnScale = NULL;
3251 if (!model->rowScale()) {
3252 rowScale = new double [numberRows*2];
3253 columnScale = new double [numberColumns*2];
3254 inverseRowScale = rowScale + numberRows;
3255 inverseColumnScale = columnScale + numberColumns;
3256 arraysExist = false;
3257 } else {
3258 rowScale = model->mutableRowScale();
3259 columnScale = model->mutableColumnScale();
3260 inverseRowScale = model->mutableInverseRowScale();
3261 inverseColumnScale = model->mutableInverseColumnScale();
3262 arraysExist = true;
3263 }
3264 assert (inverseRowScale == rowScale + numberRows);
3265 assert (inverseColumnScale == columnScale + numberColumns);
3266 // we are going to mark bits we are interested in
3267 char * usefulColumn = new char [numberColumns];
3268 double * rowLower = model->rowLower();
3269 double * rowUpper = model->rowUpper();
3270 double * columnLower = model->columnLower();
3271 double * columnUpper = model->columnUpper();
3272 int iColumn, iRow;
3273 //#define LEAVE_FIXED
3274 // mark empty and fixed columns
3275 // also see if worth scaling
3276 assert (model->scalingFlag() <= 5);
3277 // scale_stats[model->scalingFlag()]++;
3278 double largest = 0.0;
3279 double smallest = 1.0e50;
3280 // get matrix data pointers
3281 int * row = matrix_->getMutableIndices();
3282 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3283 int * columnLength = matrix_->getMutableVectorLengths();
3284 double * elementByColumn = matrix_->getMutableElements();
3285 int deletedElements = 0;
3286 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3287 CoinBigIndex j;
3288 char useful = 0;
3289 bool deleteSome = false;
3290 int start = columnStart[iColumn];
3291 int end = start + columnLength[iColumn];
3292 #ifndef LEAVE_FIXED
3293 if (columnUpper[iColumn] >
3294 columnLower[iColumn] + 1.0e-12) {
3295 #endif
3296 for (j = start; j < end; j++) {
3297 iRow = row[j];
3298 double value = fabs(elementByColumn[j]);
3299 if (value > 1.0e-20) {
3300 useful = 1;
3301 largest = CoinMax(largest, fabs(elementByColumn[j]));
3302 smallest = CoinMin(smallest, fabs(elementByColumn[j]));
3303 } else {
3304 // small
3305 deleteSome = true;
3306 }
3307 }
3308 #ifndef LEAVE_FIXED
3309 } else {
3310 // just check values
3311 for (j = start; j < end; j++) {
3312 double value = fabs(elementByColumn[j]);
3313 if (value <= 1.0e-20) {
3314 // small
3315 deleteSome = true;
3316 }
3317 }
3318 }
3319 #endif
3320 usefulColumn[iColumn] = useful;
3321 if (deleteSome) {
3322 CoinBigIndex put = start;
3323 for (j = start; j < end; j++) {
3324 double value = elementByColumn[j];
3325 if (fabs(value) > 1.0e-20) {
3326 row[put] = row[j];
3327 elementByColumn[put++] = value;
3328 }
3329 }
3330 deletedElements += end - put;
3331 columnLength[iColumn] = put - start;
3332 }
3333 }
3334 if (deletedElements)
3335 matrix_->setNumElements(matrix_->getNumElements()-deletedElements);
3336 model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL, *model->messagesPointer())
3337 << smallest << largest
3338 << CoinMessageEol;
3339 if (smallest >= 0.5 && largest <= 2.0 && !deletedElements) {
3340 // don't bother scaling
3341 model->messageHandler()->message(CLP_PACKEDSCALE_FORGET, *model->messagesPointer())
3342 << CoinMessageEol;
3343 if (!arraysExist) {
3344 delete [] rowScale;
3345 delete [] columnScale;
3346 } else {
3347 model->setRowScale(NULL);
3348 model->setColumnScale(NULL);
3349 }
3350 delete [] usefulColumn;
3351 return 1;
3352 } else {
3353 #ifdef CLP_INVESTIGATE
3354 if (deletedElements)
3355 printf("DEL_ELS\n");
3356 #endif
3357 if (!rowCopyBase) {
3358 // temporary copy
3359 rowCopyBase = reverseOrderedCopy();
3360 } else if (deletedElements) {
3361 rowCopyBase = reverseOrderedCopy();
3362 model->setNewRowCopy(rowCopyBase);
3363 }
3364 #ifndef NDEBUG
3365 ClpPackedMatrix* rowCopy =
3366 dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3367 // Make sure it is really a ClpPackedMatrix
3368 assert (rowCopy != NULL);
3369 #else
3370 ClpPackedMatrix* rowCopy =
3371 static_cast< ClpPackedMatrix*>(rowCopyBase);
3372 #endif
3373
3374 const int * column = rowCopy->getIndices();
3375 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3376 const double * element = rowCopy->getElements();
3377 // need to scale
3378 if (largest > 1.0e13 * smallest) {
3379 // safer to have smaller zero tolerance
3380 double ratio = smallest / largest;
3381 ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3382 double newTolerance = CoinMax(ratio * 0.5, 1.0e-18);
3383 if (simplex->zeroTolerance() > newTolerance)
3384 simplex->setZeroTolerance(newTolerance);
3385 }
3386 int scalingMethod = model->scalingFlag();
3387 if (scalingMethod == 4) {
3388 // As auto
3389 scalingMethod = 3;
3390 } else if (scalingMethod == 5) {
3391 // As geometric
3392 scalingMethod = 2;
3393 }
3394 double savedOverallRatio = 0.0;
3395 double tolerance = 5.0 * model->primalTolerance();
3396 double overallLargest = -1.0e-20;
3397 double overallSmallest = 1.0e20;
3398 bool finished = false;
3399 // if scalingMethod 3 then may change
3400 bool extraDetails = (model->logLevel() > 2);
3401 while (!finished) {
3402 int numberPass = 3;
3403 overallLargest = -1.0e-20;
3404 overallSmallest = 1.0e20;
3405 ClpFillN ( rowScale, numberRows, 1.0);
3406 ClpFillN ( columnScale, numberColumns, 1.0);
3407 if (scalingMethod == 1 || scalingMethod == 3) {
3408 // Maximum in each row
3409 for (iRow = 0; iRow < numberRows; iRow++) {
3410 CoinBigIndex j;
3411 largest = 1.0e-10;
3412 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
3413 int iColumn = column[j];
3414 if (usefulColumn[iColumn]) {
3415 double value = fabs(element[j]);
3416 largest = CoinMax(largest, value);
3417 assert (largest < 1.0e40);
3418 }
3419 }
3420 rowScale[iRow] = 1.0 / largest;
3421 #ifdef COIN_DEVELOP
3422 if (extraDetails) {
3423 overallLargest = CoinMax(overallLargest, largest);
3424 overallSmallest = CoinMin(overallSmallest, largest);
3425 }
3426 #endif
3427 }
3428 } else {
3429 #ifdef USE_OBJECTIVE
3430 // This will be used to help get scale factors
3431 double * objective = new double[numberColumns];
3432 CoinMemcpyN(model->costRegion(1), numberColumns, objective);
3433 double objScale = 1.0;
3434 #endif
3435 while (numberPass) {
3436 overallLargest = 0.0;
3437 overallSmallest = 1.0e50;
3438 numberPass--;
3439 // Geometric mean on row scales
3440 for (iRow = 0; iRow < numberRows; iRow++) {
3441 CoinBigIndex j;
3442 largest = 1.0e-50;
3443 smallest = 1.0e50;
3444 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
3445 int iColumn = column[j];
3446 if (usefulColumn[iColumn]) {
3447 double value = fabs(element[j]);
3448 value *= columnScale[iColumn];
3449 largest = CoinMax(largest, value);
3450 smallest = CoinMin(smallest, value);
3451 }
3452 }
3453
3454 #ifdef SQRT_ARRAY
3455 rowScale[iRow] = smallest * largest;
3456 #else
3457 rowScale[iRow] = 1.0 / sqrt(smallest * largest);
3458 #endif
3459 //rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
3460 if (extraDetails) {
3461 overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
3462 overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
3463 }
3464 }
3465 if (model->scalingFlag() == 5)
3466 break; // just scale rows
3467 #ifdef SQRT_ARRAY
3468 doSqrts(rowScale, numberRows);
3469 #endif
3470 #ifdef USE_OBJECTIVE
3471 largest = 1.0e-20;
3472 smallest = 1.0e50;
3473 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3474 if (usefulColumn[iColumn]) {
3475 double value = fabs(objective[iColumn]);
3476 value *= columnScale[iColumn];
3477 largest = CoinMax(largest, value);
3478 smallest = CoinMin(smallest, value);
3479 }
3480 }
3481 objScale = 1.0 / sqrt(smallest * largest);
3482 #endif
3483 model->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model->messagesPointer())
3484 << overallSmallest
3485 << overallLargest
3486 << CoinMessageEol;
3487 // skip last column round
3488 if (numberPass == 1)
3489 break;
3490 // Geometric mean on column scales
3491 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3492 if (usefulColumn[iColumn]) {
3493 CoinBigIndex j;
3494 largest = 1.0e-50;
3495 smallest = 1.0e50;
3496 for (j = columnStart[iColumn];
3497 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3498 iRow = row[j];
3499 double value = fabs(elementByColumn[j]);
3500 value *= rowScale[iRow];
3501 largest = CoinMax(largest, value);
3502 smallest = CoinMin(smallest, value);
3503 }
3504 #ifdef USE_OBJECTIVE
3505 if (fabs(objective[iColumn]) > 1.0e-20) {
3506 double value = fabs(objective[iColumn]) * objScale;
3507 largest = CoinMax(largest, value);
3508 smallest = CoinMin(smallest, value);
3509 }
3510 #endif
3511 #ifdef SQRT_ARRAY
3512 columnScale[iColumn] = smallest * largest;
3513 #else
3514 columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
3515 #endif
3516 //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3517 }
3518 }
3519 #ifdef SQRT_ARRAY
3520 doSqrts(columnScale, numberColumns);
3521 #endif
3522 }
3523 #ifdef USE_OBJECTIVE
3524 delete [] objective;
3525 printf("obj scale %g - use it if you want\n", objScale);
3526 #endif
3527 }
3528 // If ranges will make horrid then scale
3529 for (iRow = 0; iRow < numberRows; iRow++) {
3530 double difference = rowUpper[iRow] - rowLower[iRow];
3531 double scaledDifference = difference * rowScale[iRow];
3532 if (scaledDifference > tolerance && scaledDifference < 1.0e-4) {
3533 // make gap larger
3534 rowScale[iRow] *= 1.0e-4 / scaledDifference;
3535 rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
3536 //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
3537 // scaledDifference,difference*rowScale[iRow]);
3538 }
3539 }
3540 // final pass to scale columns so largest is reasonable
3541 // See what smallest will be if largest is 1.0
3542 if (model->scalingFlag() != 5) {
3543 overallSmallest = 1.0e50;
3544 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3545 if (usefulColumn[iColumn]) {
3546 CoinBigIndex j;
3547 largest = 1.0e-20;
3548 smallest = 1.0e50;
3549 for (j = columnStart[iColumn];
3550 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3551 iRow = row[j];
3552 double value = fabs(elementByColumn[j] * rowScale[iRow]);
3553 largest = CoinMax(largest, value);
3554 smallest = CoinMin(smallest, value);
3555 }
3556 if (overallSmallest * largest > smallest)
3557 overallSmallest = smallest / largest;
3558 }
3559 }
3560 }
3561 if (scalingMethod == 1 || scalingMethod == 2) {
3562 finished = true;
3563 } else if (savedOverallRatio == 0.0 && scalingMethod != 4) {
3564 savedOverallRatio = overallSmallest;
3565 scalingMethod = 4;
3566 } else {
3567 assert (scalingMethod == 4);
3568 if (overallSmallest > 2.0 * savedOverallRatio) {
3569 finished = true; // geometric was better
3570 if (model->scalingFlag() == 4) {
3571 // If in Branch and bound change
3572 if ((model->specialOptions() & 1024) != 0) {
3573 model->scaling(2);
3574 }
3575 }
3576 } else {
3577 scalingMethod = 1; // redo equilibrium
3578 if (model->scalingFlag() == 4) {
3579 // If in Branch and bound change
3580 if ((model->specialOptions() & 1024) != 0) {
3581 model->scaling(1);
3582 }
3583 }
3584 }
3585 #if 0
3586 if (extraDetails) {
3587 if (finished)
3588 printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
3589 savedOverallRatio, overallSmallest);
3590 else
3591 printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
3592 savedOverallRatio, overallSmallest);
3593 }
3594 #endif
3595 }
3596 }
3597 //#define RANDOMIZE
3598 #ifdef RANDOMIZE
3599 // randomize by up to 10%
3600 for (iRow = 0; iRow < numberRows; iRow++) {
3601 double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3602 rowScale[iRow] *= (1.0 + 0.1 * value);
3603 }
3604 #endif
3605 overallLargest = 1.0;
3606 if (overallSmallest < 1.0e-1)
3607 overallLargest = 1.0 / sqrt(overallSmallest);
3608 overallLargest = CoinMin(100.0, overallLargest);
3609 overallSmallest = 1.0e50;
3610 char * usedRow = reinterpret_cast<char *>(inverseRowScale);
3611 memset(usedRow, 0, numberRows);
3612 //printf("scaling %d\n",model->scalingFlag());
3613 if (model->scalingFlag() != 5) {
3614 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3615 if (columnUpper[iColumn] >
3616 columnLower[iColumn] + 1.0e-12) {
3617 //if (usefulColumn[iColumn]) {
3618 CoinBigIndex j;
3619 largest = 1.0e-20;
3620 smallest = 1.0e50;
3621 for (j = columnStart[iColumn];
3622 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3623 iRow = row[j];
3624 usedRow[iRow] = 1;
3625 double value = fabs(elementByColumn[j] * rowScale[iRow]);
3626 largest = CoinMax(largest, value);
3627 smallest = CoinMin(smallest, value);
3628 }
3629 columnScale[iColumn] = overallLargest / largest;
3630 //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3631 #ifdef RANDOMIZE
3632 double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3633 columnScale[iColumn] *= (1.0 + 0.1 * value);
3634 #endif
3635 double difference = columnUpper[iColumn] - columnLower[iColumn];
3636 if (difference < 1.0e-5 * columnScale[iColumn]) {
3637 // make gap larger
3638 columnScale[iColumn] = difference / 1.0e-5;
3639 //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
3640 // scaledDifference,difference*columnScale[iColumn]);
3641 }
3642 double value = smallest * columnScale[iColumn];
3643 if (overallSmallest > value)
3644 overallSmallest = value;
3645 //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
3646 } else {
3647 assert(columnScale[iColumn] == 1.0);
3648 //columnScale[iColumn]=1.0;
3649 }
3650 }
3651 for (iRow = 0; iRow < numberRows; iRow++) {
3652 if (!usedRow[iRow]) {
3653 rowScale[iRow] = 1.0;
3654 }
3655 }
3656 }
3657 model->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model->messagesPointer())
3658 << overallSmallest
3659 << overallLargest
3660 << CoinMessageEol;
3661 #if 0
3662 {
3663 for (iRow = 0; iRow < numberRows; iRow++) {
3664 assert (rowScale[iRow] == rowScale2[iRow]);
3665 }
3666 delete [] rowScale2;
3667 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3668 assert (columnScale[iColumn] == columnScale2[iColumn]);
3669 }
3670 delete [] columnScale2;
3671 }
3672 #endif
3673 if (overallSmallest < 1.0e-13) {
3674 // Change factorization zero tolerance
3675 double newTolerance = CoinMax(1.0e-15 * (overallSmallest / 1.0e-13),
3676 1.0e-18);
3677 ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3678 if (simplex->factorization()->zeroTolerance() > newTolerance)
3679 simplex->factorization()->zeroTolerance(newTolerance);
3680 newTolerance = CoinMax(overallSmallest * 0.5, 1.0e-18);
3681 simplex->setZeroTolerance(newTolerance);
3682 }
3683 delete [] usefulColumn;
3684 #ifndef SLIM_CLP
3685 // If quadratic then make symmetric
3686 ClpObjective * obj = model->objectiveAsObject();
3687 #ifndef NO_RTTI
3688 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
3689 #else
3690 ClpQuadraticObjective * quadraticObj = NULL;
3691 if (obj->type() == 2)
3692 quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
3693 #endif
3694 if (quadraticObj) {
3695 if (!rowCopyBase) {
3696 // temporary copy
3697 rowCopyBase = reverseOrderedCopy();
3698 }
3699 #ifndef NDEBUG
3700 ClpPackedMatrix* rowCopy =
3701 dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3702 // Make sure it is really a ClpPackedMatrix
3703 assert (rowCopy != NULL);
3704 #else
3705 ClpPackedMatrix* rowCopy =
3706 static_cast< ClpPackedMatrix*>(rowCopyBase);
3707 #endif
3708 const int * column = rowCopy->getIndices();
3709 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3710 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
3711 int numberXColumns = quadratic->getNumCols();
3712 if (numberXColumns < numberColumns) {
3713 // we assume symmetric
3714 int numberQuadraticColumns = 0;
3715 int i;
3716 //const int * columnQuadratic = quadratic->getIndices();
3717 //const int * columnQuadraticStart = quadratic->getVectorStarts();
3718 const int * columnQuadraticLength = quadratic->getVectorLengths();
3719 for (i = 0; i < numberXColumns; i++) {
3720 int length = columnQuadraticLength[i];
3721 #ifndef CORRECT_COLUMN_COUNTS
3722 length = 1;
3723 #endif
3724 if (length)
3725 numberQuadraticColumns++;
3726 }
3727 int numberXRows = numberRows - numberQuadraticColumns;
3728 numberQuadraticColumns = 0;
3729 for (i = 0; i < numberXColumns; i++) {
3730 int length = columnQuadraticLength[i];
3731 #ifndef CORRECT_COLUMN_COUNTS
3732 length = 1;
3733 #endif
3734 if (length) {
3735 rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
3736 numberQuadraticColumns++;
3737 }
3738 }
3739 int numberQuadraticRows = 0;
3740 for (i = 0; i < numberXRows; i++) {
3741 // See if any in row quadratic
3742 CoinBigIndex j;
3743 int numberQ = 0;
3744 for (j = rowStart[i]; j < rowStart[i+1]; j++) {
3745 int iColumn = column[j];
3746 if (columnQuadraticLength[iColumn])
3747 numberQ++;
3748 }
3749 #ifndef CORRECT_ROW_COUNTS
3750 numberQ = 1;
3751 #endif
3752 if (numberQ) {
3753 columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
3754 numberQuadraticRows++;
3755 }
3756 }
3757 // and make sure Sj okay
3758 for (iColumn = numberQuadraticRows + numberXColumns; iColumn < numberColumns; iColumn++) {
3759 CoinBigIndex j = columnStart[iColumn];
3760 assert(columnLength[iColumn] == 1);
3761 int iRow = row[j];
3762 double value = fabs(elementByColumn[j] * rowScale[iRow]);
3763 columnScale[iColumn] = 1.0 / value;
3764 }
3765 }
3766 }
3767 #endif
3768 // make copy (could do faster by using previous values)
3769 // could just do partial
3770 for (iRow = 0; iRow < numberRows; iRow++)
3771 inverseRowScale[iRow] = 1.0 / rowScale[iRow] ;
3772 for (iColumn = 0; iColumn < numberColumns; iColumn++)
3773 inverseColumnScale[iColumn] = 1.0 / columnScale[iColumn] ;
3774 if (!arraysExist) {
3775 model->setRowScale(rowScale);
3776 model->setColumnScale(columnScale);
3777 }
3778 if (model->rowCopy()) {
3779 // need to replace row by row
3780 ClpPackedMatrix* rowCopy =
3781 static_cast< ClpPackedMatrix*>(model->rowCopy());
3782 double * element = rowCopy->getMutableElements();
3783 const int * column = rowCopy->getIndices();
3784 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3785 // scale row copy
3786 for (iRow = 0; iRow < numberRows; iRow++) {
3787 CoinBigIndex j;
3788 double scale = rowScale[iRow];
3789 double * elementsInThisRow = element + rowStart[iRow];
3790 const int * columnsInThisRow = column + rowStart[iRow];
3791 int number = rowStart[iRow+1] - rowStart[iRow];
3792 assert (number <= numberColumns);
3793 for (j = 0; j < number; j++) {
3794 int iColumn = columnsInThisRow[j];
3795 elementsInThisRow[j] *= scale * columnScale[iColumn];
3796 }
3797 }
3798 if ((model->specialOptions() & 262144) != 0) {
3799 //if ((model->specialOptions()&(COIN_CBC_USING_CLP|16384))!=0) {
3800 //if (model->inCbcBranchAndBound()&&false) {
3801 // copy without gaps
3802 CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, 0, 0);
3803 ClpPackedMatrix * scaled = new ClpPackedMatrix(scaledMatrix);
3804 model->setClpScaledMatrix(scaled);
3805 // get matrix data pointers
3806 const int * row = scaledMatrix->getIndices();
3807 const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3808 #ifndef NDEBUG
3809 const int * columnLength = scaledMatrix->getVectorLengths();
3810 #endif
3811 double * elementByColumn = scaledMatrix->getMutableElements();
3812 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3813 CoinBigIndex j;
3814 double scale = columnScale[iColumn];
3815 assert (columnStart[iColumn+1] == columnStart[iColumn] + columnLength[iColumn]);
3816 for (j = columnStart[iColumn];
3817 j < columnStart[iColumn+1]; j++) {
3818 int iRow = row[j];
3819 elementByColumn[j] *= scale * rowScale[iRow];
3820 }
3821 }
3822 } else {
3823 //printf("not in b&b\n");
3824 }
3825 } else {
3826 // no row copy
3827 delete rowCopyBase;
3828 }
3829 return 0;
3830 }
3831 }
3832 // Creates scaled column copy if scales exist
3833 void
createScaledMatrix(ClpSimplex * model) const3834 ClpPackedMatrix::createScaledMatrix(ClpSimplex * model) const
3835 {
3836 int numberRows = model->numberRows();
3837 int numberColumns = matrix_->getNumCols();
3838 model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
3839 // If empty - return as sanityCheck will trap
3840 if (!numberRows || !numberColumns) {
3841 model->setRowScale(NULL);
3842 model->setColumnScale(NULL);
3843 return ;
3844 }
3845 if (!model->rowScale())
3846 return;
3847 double * rowScale = model->mutableRowScale();
3848 double * columnScale = model->mutableColumnScale();
3849 // copy without gaps
3850 CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, 0, 0);
3851 ClpPackedMatrix * scaled = new ClpPackedMatrix(scaledMatrix);
3852 model->setClpScaledMatrix(scaled);
3853 // get matrix data pointers
3854 const int * row = scaledMatrix->getIndices();
3855 const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3856 #ifndef NDEBUG
3857 const int * columnLength = scaledMatrix->getVectorLengths();
3858 #endif
3859 double * elementByColumn = scaledMatrix->getMutableElements();
3860 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3861 CoinBigIndex j;
3862 double scale = columnScale[iColumn];
3863 assert (columnStart[iColumn+1] == columnStart[iColumn] + columnLength[iColumn]);
3864 for (j = columnStart[iColumn];
3865 j < columnStart[iColumn+1]; j++) {
3866 int iRow = row[j];
3867 elementByColumn[j] *= scale * rowScale[iRow];
3868 }
3869 }
3870 }
3871 /* Unpacks a column into an CoinIndexedvector
3872 */
3873 void
unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,int iColumn) const3874 ClpPackedMatrix::unpack(const ClpSimplex * model, CoinIndexedVector * rowArray,
3875 int iColumn) const
3876 {
3877 const double * rowScale = model->rowScale();
3878 const int * row = matrix_->getIndices();
3879 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3880 const int * columnLength = matrix_->getVectorLengths();
3881 const double * elementByColumn = matrix_->getElements();
3882 CoinBigIndex i;
3883 if (!rowScale) {
3884 for (i = columnStart[iColumn];
3885 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3886 rowArray->add(row[i], elementByColumn[i]);
3887 }
3888 } else {
3889 // apply scaling
3890 double scale = model->columnScale()[iColumn];
3891 for (i = columnStart[iColumn];
3892 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3893 int iRow = row[i];
3894 rowArray->add(iRow, elementByColumn[i]*scale * rowScale[iRow]);
3895 }
3896 }
3897 }
3898 /* Unpacks a column into a CoinIndexedVector
3899 ** in packed format
3900 Note that model is NOT const. Bounds and objective could
3901 be modified if doing column generation (just for this variable) */
3902 void
unpackPacked(ClpSimplex * model,CoinIndexedVector * rowArray,int iColumn) const3903 ClpPackedMatrix::unpackPacked(ClpSimplex * model,
3904 CoinIndexedVector * rowArray,
3905 int iColumn) const
3906 {
3907 const double * rowScale = model->rowScale();
3908 const int * row = matrix_->getIndices();
3909 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3910 const int * columnLength = matrix_->getVectorLengths();
3911 const double * elementByColumn = matrix_->getElements();
3912 CoinBigIndex i;
3913 int * index = rowArray->getIndices();
3914 double * array = rowArray->denseVector();
3915 int number = 0;
3916 if (!rowScale) {
3917 for (i = columnStart[iColumn];
3918 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3919 int iRow = row[i];
3920 double value = elementByColumn[i];
3921 if (value) {
3922 array[number] = value;
3923 index[number++] = iRow;
3924 }
3925 }
3926 rowArray->setNumElements(number);
3927 rowArray->setPackedMode(true);
3928 } else {
3929 // apply scaling
3930 double scale = model->columnScale()[iColumn];
3931 for (i = columnStart[iColumn];
3932 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3933 int iRow = row[i];
3934 double value = elementByColumn[i] * scale * rowScale[iRow];
3935 if (value) {
3936 array[number] = value;
3937 index[number++] = iRow;
3938 }
3939 }
3940 rowArray->setNumElements(number);
3941 rowArray->setPackedMode(true);
3942 }
3943 }
3944 /* Adds multiple of a column into an CoinIndexedvector
3945 You can use quickAdd to add to vector */
3946 void
add(const ClpSimplex * model,CoinIndexedVector * rowArray,int iColumn,double multiplier) const3947 ClpPackedMatrix::add(const ClpSimplex * model, CoinIndexedVector * rowArray,
3948 int iColumn, double multiplier) const
3949 {
3950 const double * rowScale = model->rowScale();
3951 const int * row = matrix_->getIndices();
3952 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3953 const int * columnLength = matrix_->getVectorLengths();
3954 const double * elementByColumn = matrix_->getElements();
3955 CoinBigIndex i;
3956 if (!rowScale) {
3957 for (i = columnStart[iColumn];
3958 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3959 int iRow = row[i];
3960 rowArray->quickAdd(iRow, multiplier * elementByColumn[i]);
3961 }
3962 } else {
3963 // apply scaling
3964 double scale = model->columnScale()[iColumn] * multiplier;
3965 for (i = columnStart[iColumn];
3966 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3967 int iRow = row[i];
3968 rowArray->quickAdd(iRow, elementByColumn[i]*scale * rowScale[iRow]);
3969 }
3970 }
3971 }
3972 /* Adds multiple of a column into an array */
3973 void
add(const ClpSimplex * model,double * array,int iColumn,double multiplier) const3974 ClpPackedMatrix::add(const ClpSimplex * model, double * array,
3975 int iColumn, double multiplier) const
3976 {
3977 const double * rowScale = model->rowScale();
3978 const int * row = matrix_->getIndices();
3979 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3980 const int * columnLength = matrix_->getVectorLengths();
3981 const double * elementByColumn = matrix_->getElements();
3982 CoinBigIndex i;
3983 if (!rowScale) {
3984 for (i = columnStart[iColumn];
3985 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3986 int iRow = row[i];
3987 array[iRow] += multiplier * elementByColumn[i];
3988 }
3989 } else {
3990 // apply scaling
3991 double scale = model->columnScale()[iColumn] * multiplier;
3992 for (i = columnStart[iColumn];
3993 i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3994 int iRow = row[i];
3995 array[iRow] += elementByColumn[i] * scale * rowScale[iRow];
3996 }
3997 }
3998 }
3999 /* Checks if all elements are in valid range. Can just
4000 return true if you are not paranoid. For Clp I will
4001 probably expect no zeros. Code can modify matrix to get rid of
4002 small elements.
4003 check bits (can be turned off to save time) :
4004 1 - check if matrix has gaps
4005 2 - check if zero elements
4006 4 - check and compress duplicates
4007 8 - report on large and small
4008 */
4009 bool
allElementsInRange(ClpModel * model,double smallest,double largest,int check)4010 ClpPackedMatrix::allElementsInRange(ClpModel * model,
4011 double smallest, double largest,
4012 int check)
4013 {
4014 int iColumn;
4015 // make sure matrix correct size
4016 assert (matrix_->getNumRows() <= model->numberRows());
4017 if (model->clpScaledMatrix())
4018 assert (model->clpScaledMatrix()->getNumElements() == matrix_->getNumElements());
4019 assert (matrix_->getNumRows() <= model->numberRows());
4020 matrix_->setDimensions(model->numberRows(), model->numberColumns());
4021 CoinBigIndex numberLarge = 0;
4022 CoinBigIndex numberSmall = 0;
4023 CoinBigIndex numberDuplicate = 0;
4024 int firstBadColumn = -1;
4025 int firstBadRow = -1;
4026 double firstBadElement = 0.0;
4027 // get matrix data pointers
4028 const int * row = matrix_->getIndices();
4029 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4030 const int * columnLength = matrix_->getVectorLengths();
4031 const double * elementByColumn = matrix_->getElements();
4032 int numberRows = model->numberRows();
4033 int numberColumns = matrix_->getNumCols();
4034 // Say no gaps
4035 flags_ &= ~2;
4036 if (type_>=10)
4037 return true; // gub
4038 if (check == 14 || check == 10) {
4039 if (matrix_->getNumElements() < columnStart[numberColumns]) {
4040 // pack down!
4041 #if 0
4042 matrix_->removeGaps();
4043 #else
4044 checkGaps();
4045 #endif
4046 #ifdef COIN_DEVELOP
4047 //printf("flags set to 2\n");
4048 #endif
4049 } else if (numberColumns) {
4050 assert(columnStart[numberColumns] ==
4051 columnStart[numberColumns-1] + columnLength[numberColumns-1]);
4052 }
4053 return true;
4054 }
4055 assert (check == 15 || check == 11);
4056 if (check == 15) {
4057 int * mark = new int [numberRows];
4058 int i;
4059 for (i = 0; i < numberRows; i++)
4060 mark[i] = -1;
4061 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4062 CoinBigIndex j;
4063 CoinBigIndex start = columnStart[iColumn];
4064 CoinBigIndex end = start + columnLength[iColumn];
4065 if (end != columnStart[iColumn+1])
4066 flags_ |= 2;
4067 for (j = start; j < end; j++) {
4068 double value = fabs(elementByColumn[j]);
4069 int iRow = row[j];
4070 if (iRow < 0 || iRow >= numberRows) {
4071 #ifndef COIN_BIG_INDEX
4072 printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4073 #elif COIN_BIG_INDEX==0
4074 printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4075 #elif COIN_BIG_INDEX==1
4076 printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4077 #else
4078 printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4079 #endif
4080 delete [] mark;
4081 return false;
4082 }
4083 if (mark[iRow] == -1) {
4084 mark[iRow] = j;
4085 } else {
4086 // duplicate
4087 numberDuplicate++;
4088 }
4089 //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
4090 if (!value)
4091 flags_ |= 1; // there are zero elements
4092 if (value < smallest) {
4093 numberSmall++;
4094 } else if (!(value <= largest)) {
4095 numberLarge++;
4096 if (firstBadColumn < 0) {
4097 firstBadColumn = iColumn;
4098 firstBadRow = row[j];
4099 firstBadElement = elementByColumn[j];
4100 }
4101 }
4102 }
4103 //clear mark
4104 for (j = columnStart[iColumn];
4105 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4106 int iRow = row[j];
4107 mark[iRow] = -1;
4108 }
4109 }
4110 delete [] mark;
4111 } else {
4112 // just check for out of range - not for duplicates
4113 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4114 CoinBigIndex j;
4115 CoinBigIndex start = columnStart[iColumn];
4116 CoinBigIndex end = start + columnLength[iColumn];
4117 if (end != columnStart[iColumn+1])
4118 flags_ |= 2;
4119 for (j = start; j < end; j++) {
4120 double value = fabs(elementByColumn[j]);
4121 int iRow = row[j];
4122 if (iRow < 0 || iRow >= numberRows) {
4123 #ifndef COIN_BIG_INDEX
4124 printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4125 #elif COIN_BIG_INDEX==0
4126 printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4127 #elif COIN_BIG_INDEX==1
4128 printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4129 #else
4130 printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4131 #endif
4132 return false;
4133 }
4134 if (!value)
4135 flags_ |= 1; // there are zero elements
4136 if (value < smallest) {
4137 numberSmall++;
4138 } else if (!(value <= largest)) {
4139 numberLarge++;
4140 if (firstBadColumn < 0) {
4141 firstBadColumn = iColumn;
4142 firstBadRow = iRow;
4143 firstBadElement = value;
4144 }
4145 }
4146 }
4147 }
4148 }
4149 if (numberLarge) {
4150 model->messageHandler()->message(CLP_BAD_MATRIX, model->messages())
4151 << numberLarge
4152 << firstBadColumn << firstBadRow << firstBadElement
4153 << CoinMessageEol;
4154 return false;
4155 }
4156 if (numberSmall)
4157 model->messageHandler()->message(CLP_SMALLELEMENTS, model->messages())
4158 << numberSmall
4159 << CoinMessageEol;
4160 if (numberDuplicate)
4161 model->messageHandler()->message(CLP_DUPLICATEELEMENTS, model->messages())
4162 << numberDuplicate
4163 << CoinMessageEol;
4164 if (numberDuplicate)
4165 matrix_->eliminateDuplicates(smallest);
4166 else if (numberSmall)
4167 matrix_->compress(smallest);
4168 // If smallest >0.0 then there can't be zero elements
4169 if (smallest > 0.0)
4170 flags_ &= ~1;
4171 if (numberSmall || numberDuplicate)
4172 flags_ |= 2; // will have gaps
4173 return true;
4174 }
4175 int
gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector * COIN_RESTRICT piVector,int * COIN_RESTRICT index,double * COIN_RESTRICT output,int * COIN_RESTRICT lookup,char * COIN_RESTRICT marked,const double tolerance,const double scalar) const4176 ClpPackedMatrix::gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector * COIN_RESTRICT piVector,
4177 int * COIN_RESTRICT index,
4178 double * COIN_RESTRICT output,
4179 int * COIN_RESTRICT lookup,
4180 char * COIN_RESTRICT marked,
4181 const double tolerance,
4182 const double scalar) const
4183 {
4184 const double * COIN_RESTRICT pi = piVector->denseVector();
4185 int numberNonZero = 0;
4186 int numberInRowArray = piVector->getNumElements();
4187 const int * COIN_RESTRICT column = matrix_->getIndices();
4188 const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4189 const double * COIN_RESTRICT element = matrix_->getElements();
4190 const int * COIN_RESTRICT whichRow = piVector->getIndices();
4191 int * fakeRow = const_cast<int *> (whichRow);
4192 fakeRow[numberInRowArray]=0; // so can touch
4193 #ifndef NDEBUG
4194 int maxColumn = 0;
4195 #endif
4196 // ** Row copy is already scaled
4197 int nextRow=whichRow[0];
4198 CoinBigIndex nextStart = rowStart[nextRow];
4199 CoinBigIndex nextEnd = rowStart[nextRow+1];
4200 for (int i = 0; i < numberInRowArray; i++) {
4201 double value = pi[i] * scalar;
4202 CoinBigIndex start=nextStart;
4203 CoinBigIndex end=nextEnd;
4204 nextRow=whichRow[i+1];
4205 nextStart = rowStart[nextRow];
4206 //coin_prefetch_const(column + nextStart);
4207 //coin_prefetch_const(element + nextStart);
4208 nextEnd = rowStart[nextRow+1];
4209 CoinBigIndex j;
4210 for (j = start; j < end; j++) {
4211 int iColumn = column[j];
4212 #ifndef NDEBUG
4213 maxColumn = CoinMax(maxColumn, iColumn);
4214 #endif
4215 double elValue = element[j];
4216 elValue *= value;
4217 if (marked[iColumn]) {
4218 int k = lookup[iColumn];
4219 output[k] += elValue;
4220 } else {
4221 output[numberNonZero] = elValue;
4222 marked[iColumn] = 1;
4223 lookup[iColumn] = numberNonZero;
4224 index[numberNonZero++] = iColumn;
4225 }
4226 }
4227 }
4228 #ifndef NDEBUG
4229 int saveN = numberNonZero;
4230 #endif
4231 // get rid of tiny values and zero out lookup
4232 for (int i = 0; i < numberNonZero; i++) {
4233 int iColumn = index[i];
4234 marked[iColumn] = 0;
4235 double value = output[i];
4236 if (fabs(value) <= tolerance) {
4237 while (fabs(value) <= tolerance) {
4238 numberNonZero--;
4239 value = output[numberNonZero];
4240 iColumn = index[numberNonZero];
4241 marked[iColumn] = 0;
4242 if (i < numberNonZero) {
4243 output[numberNonZero] = 0.0;
4244 output[i] = value;
4245 index[i] = iColumn;
4246 } else {
4247 output[i] = 0.0;
4248 value = 1.0; // to force end of while
4249 }
4250 }
4251 }
4252 }
4253 #ifndef NDEBUG
4254 for (int i = numberNonZero; i < saveN; i++)
4255 assert(!output[i]);
4256 for (int i = 0; i <= maxColumn; i++)
4257 assert (!marked[i]);
4258 #endif
4259 return numberNonZero;
4260 }
4261 int
gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,int * COIN_RESTRICT index,double * COIN_RESTRICT output,double * COIN_RESTRICT array,const double tolerance,const double scalar) const4262 ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
4263 int * COIN_RESTRICT index,
4264 double * COIN_RESTRICT output,
4265 double * COIN_RESTRICT array,
4266 const double tolerance,
4267 const double scalar) const
4268 {
4269 const double * COIN_RESTRICT pi = piVector->denseVector();
4270 int numberNonZero = 0;
4271 int numberInRowArray = piVector->getNumElements();
4272 const int * COIN_RESTRICT column = matrix_->getIndices();
4273 const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4274 const double * COIN_RESTRICT element = matrix_->getElements();
4275 const int * COIN_RESTRICT whichRow = piVector->getIndices();
4276 // ** Row copy is already scaled
4277 for (int i = 0; i < numberInRowArray; i++) {
4278 int iRow = whichRow[i];
4279 double value = pi[i] * scalar;
4280 CoinBigIndex j;
4281 for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
4282 int iColumn = column[j];
4283 double inValue = array[iColumn];
4284 double elValue = element[j];
4285 elValue *= value;
4286 if (inValue) {
4287 double outValue = inValue + elValue;
4288 if (!outValue)
4289 outValue = COIN_INDEXED_REALLY_TINY_ELEMENT;
4290 array[iColumn] = outValue;
4291 } else {
4292 array[iColumn] = elValue;
4293 assert (elValue);
4294 index[numberNonZero++] = iColumn;
4295 }
4296 }
4297 }
4298 int saveN = numberNonZero;
4299 // get rid of tiny values
4300 numberNonZero=0;
4301 for (int i = 0; i < saveN; i++) {
4302 int iColumn = index[i];
4303 double value = array[iColumn];
4304 array[iColumn] =0.0;
4305 if (fabs(value) > tolerance) {
4306 output[numberNonZero] = value;
4307 index[numberNonZero++] = iColumn;
4308 }
4309 }
4310 return numberNonZero;
4311 }
4312 /* Given positive integer weights for each row fills in sum of weights
4313 for each column (and slack).
4314 Returns weights vector
4315 */
4316 CoinBigIndex *
dubiousWeights(const ClpSimplex * model,int * inputWeights) const4317 ClpPackedMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
4318 {
4319 int numberRows = model->numberRows();
4320 int numberColumns = matrix_->getNumCols();
4321 int number = numberRows + numberColumns;
4322 CoinBigIndex * weights = new CoinBigIndex[number];
4323 // get matrix data pointers
4324 const int * row = matrix_->getIndices();
4325 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4326 const int * columnLength = matrix_->getVectorLengths();
4327 int i;
4328 for (i = 0; i < numberColumns; i++) {
4329 CoinBigIndex j;
4330 CoinBigIndex count = 0;
4331 for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4332 int iRow = row[j];
4333 count += inputWeights[iRow];
4334 }
4335 weights[i] = count;
4336 }
4337 for (i = 0; i < numberRows; i++) {
4338 weights[i+numberColumns] = inputWeights[i];
4339 }
4340 return weights;
4341 }
4342 /* Returns largest and smallest elements of both signs.
4343 Largest refers to largest absolute value.
4344 */
4345 void
rangeOfElements(double & smallestNegative,double & largestNegative,double & smallestPositive,double & largestPositive)4346 ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
4347 double & smallestPositive, double & largestPositive)
4348 {
4349 smallestNegative = -COIN_DBL_MAX;
4350 largestNegative = 0.0;
4351 smallestPositive = COIN_DBL_MAX;
4352 largestPositive = 0.0;
4353 // get matrix data pointers
4354 const double * elementByColumn = matrix_->getElements();
4355 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4356 const int * columnLength = matrix_->getVectorLengths();
4357 int numberColumns = matrix_->getNumCols();
4358 int i;
4359 for (i = 0; i < numberColumns; i++) {
4360 CoinBigIndex j;
4361 for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4362 double value = elementByColumn[j];
4363 if (value > 0.0) {
4364 smallestPositive = CoinMin(smallestPositive, value);
4365 largestPositive = CoinMax(largestPositive, value);
4366 } else if (value < 0.0) {
4367 smallestNegative = CoinMax(smallestNegative, value);
4368 largestNegative = CoinMin(largestNegative, value);
4369 }
4370 }
4371 }
4372 }
4373 // Says whether it can do partial pricing
4374 bool
canDoPartialPricing() const4375 ClpPackedMatrix::canDoPartialPricing() const
4376 {
4377 return true;
4378 }
4379 // Partial pricing
4380 void
partialPricing(ClpSimplex * model,double startFraction,double endFraction,int & bestSequence,int & numberWanted)4381 ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
4382 int & bestSequence, int & numberWanted)
4383 {
4384 numberWanted = currentWanted_;
4385 int start = static_cast<int> (startFraction * numberActiveColumns_);
4386 int end = CoinMin(static_cast<int> (endFraction * numberActiveColumns_ + 1), numberActiveColumns_);
4387 const double * element = matrix_->getElements();
4388 const int * row = matrix_->getIndices();
4389 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
4390 const int * length = matrix_->getVectorLengths();
4391 const double * rowScale = model->rowScale();
4392 const double * columnScale = model->columnScale();
4393 int iSequence;
4394 CoinBigIndex j;
4395 double tolerance = model->currentDualTolerance();
4396 double * reducedCost = model->djRegion();
4397 const double * duals = model->dualRowSolution();
4398 const double * cost = model->costRegion();
4399 double bestDj;
4400 if (bestSequence >= 0)
4401 bestDj = fabs(model->clpMatrix()->reducedCost(model, bestSequence));
4402 else
4403 bestDj = tolerance;
4404 int sequenceOut = model->sequenceOut();
4405 int saveSequence = bestSequence;
4406 int lastScan = minimumObjectsScan_ < 0 ? end : start + minimumObjectsScan_;
4407 int minNeg = minimumGoodReducedCosts_ == -1 ? numberWanted : minimumGoodReducedCosts_;
4408 if (rowScale) {
4409 // scaled
4410 for (iSequence = start; iSequence < end; iSequence++) {
4411 if (iSequence != sequenceOut) {
4412 double value;
4413 ClpSimplex::Status status = model->getStatus(iSequence);
4414
4415 switch(status) {
4416
4417 case ClpSimplex::basic:
4418 case ClpSimplex::isFixed:
4419 break;
4420 case ClpSimplex::isFree:
4421 case ClpSimplex::superBasic:
4422 value = 0.0;
4423 // scaled
4424 for (j = startColumn[iSequence];
4425 j < startColumn[iSequence] + length[iSequence]; j++) {
4426 int jRow = row[j];
4427 value -= duals[jRow] * element[j] * rowScale[jRow];
4428 }
4429 value = fabs(cost[iSequence] + value * columnScale[iSequence]);
4430 if (value > FREE_ACCEPT * tolerance) {
4431 numberWanted--;
4432 // we are going to bias towards free (but only if reasonable)
4433 value *= FREE_BIAS;
4434 if (value > bestDj) {
4435 // check flagged variable and correct dj
4436 if (!model->flagged(iSequence)) {
4437 bestDj = value;
4438 bestSequence = iSequence;
4439 } else {
4440 // just to make sure we don't exit before got something
4441 numberWanted++;
4442 }
4443 }
4444 }
4445 break;
4446 case ClpSimplex::atUpperBound:
4447 value = 0.0;
4448 // scaled
4449 for (j = startColumn[iSequence];
4450 j < startColumn[iSequence] + length[iSequence]; j++) {
4451 int jRow = row[j];
4452 value -= duals[jRow] * element[j] * rowScale[jRow];
4453 }
4454 value = cost[iSequence] + value * columnScale[iSequence];
4455 if (value > tolerance) {
4456 numberWanted--;
4457 if (value > bestDj) {
4458 // check flagged variable and correct dj
4459 if (!model->flagged(iSequence)) {
4460 bestDj = value;
4461 bestSequence = iSequence;
4462 } else {
4463 // just to make sure we don't exit before got something
4464 numberWanted++;
4465 }
4466 }
4467 }
4468 break;
4469 case ClpSimplex::atLowerBound:
4470 value = 0.0;
4471 // scaled
4472 for (j = startColumn[iSequence];
4473 j < startColumn[iSequence] + length[iSequence]; j++) {
4474 int jRow = row[j];
4475 value -= duals[jRow] * element[j] * rowScale[jRow];
4476 }
4477 value = -(cost[iSequence] + value * columnScale[iSequence]);
4478 if (value > tolerance) {
4479 numberWanted--;
4480 if (value > bestDj) {
4481 // check flagged variable and correct dj
4482 if (!model->flagged(iSequence)) {
4483 bestDj = value;
4484 bestSequence = iSequence;
4485 } else {
4486 // just to make sure we don't exit before got something
4487 numberWanted++;
4488 }
4489 }
4490 }
4491 break;
4492 }
4493 }
4494 if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4495 // give up
4496 break;
4497 }
4498 if (!numberWanted)
4499 break;
4500 }
4501 if (bestSequence != saveSequence) {
4502 // recompute dj
4503 double value = 0.0;
4504 // scaled
4505 for (j = startColumn[bestSequence];
4506 j < startColumn[bestSequence] + length[bestSequence]; j++) {
4507 int jRow = row[j];
4508 value -= duals[jRow] * element[j] * rowScale[jRow];
4509 }
4510 reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
4511 savedBestSequence_ = bestSequence;
4512 savedBestDj_ = reducedCost[savedBestSequence_];
4513 }
4514 } else {
4515 // not scaled
4516 for (iSequence = start; iSequence < end; iSequence++) {
4517 if (iSequence != sequenceOut) {
4518 double value;
4519 ClpSimplex::Status status = model->getStatus(iSequence);
4520
4521 switch(status) {
4522
4523 case ClpSimplex::basic:
4524 case ClpSimplex::isFixed:
4525 break;
4526 case ClpSimplex::isFree:
4527 case ClpSimplex::superBasic:
4528 value = cost[iSequence];
4529 for (j = startColumn[iSequence];
4530 j < startColumn[iSequence] + length[iSequence]; j++) {
4531 int jRow = row[j];
4532 value -= duals[jRow] * element[j];
4533 }
4534 value = fabs(value);
4535 if (value > FREE_ACCEPT * tolerance) {
4536 numberWanted--;
4537 // we are going to bias towards free (but only if reasonable)
4538 value *= FREE_BIAS;
4539 if (value > bestDj) {
4540 // check flagged variable and correct dj
4541 if (!model->flagged(iSequence)) {
4542 bestDj = value;
4543 bestSequence = iSequence;
4544 } else {
4545 // just to make sure we don't exit before got something
4546 numberWanted++;
4547 }
4548 }
4549 }
4550 break;
4551 case ClpSimplex::atUpperBound:
4552 value = cost[iSequence];
4553 // scaled
4554 for (j = startColumn[iSequence];
4555 j < startColumn[iSequence] + length[iSequence]; j++) {
4556 int jRow = row[j];
4557 value -= duals[jRow] * element[j];
4558 }
4559 if (value > tolerance) {
4560 numberWanted--;
4561 if (value > bestDj) {
4562 // check flagged variable and correct dj
4563 if (!model->flagged(iSequence)) {
4564 bestDj = value;
4565 bestSequence = iSequence;
4566 } else {
4567 // just to make sure we don't exit before got something
4568 numberWanted++;
4569 }
4570 }
4571 }
4572 break;
4573 case ClpSimplex::atLowerBound:
4574 value = cost[iSequence];
4575 for (j = startColumn[iSequence];
4576 j < startColumn[iSequence] + length[iSequence]; j++) {
4577 int jRow = row[j];
4578 value -= duals[jRow] * element[j];
4579 }
4580 value = -value;
4581 if (value > tolerance) {
4582 numberWanted--;
4583 if (value > bestDj) {
4584 // check flagged variable and correct dj
4585 if (!model->flagged(iSequence)) {
4586 bestDj = value;
4587 bestSequence = iSequence;
4588 } else {
4589 // just to make sure we don't exit before got something
4590 numberWanted++;
4591 }
4592 }
4593 }
4594 break;
4595 }
4596 }
4597 if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4598 // give up
4599 break;
4600 }
4601 if (!numberWanted)
4602 break;
4603 }
4604 if (bestSequence != saveSequence) {
4605 // recompute dj
4606 double value = cost[bestSequence];
4607 for (j = startColumn[bestSequence];
4608 j < startColumn[bestSequence] + length[bestSequence]; j++) {
4609 int jRow = row[j];
4610 value -= duals[jRow] * element[j];
4611 }
4612 reducedCost[bestSequence] = value;
4613 savedBestSequence_ = bestSequence;
4614 savedBestDj_ = reducedCost[savedBestSequence_];
4615 }
4616 }
4617 currentWanted_ = numberWanted;
4618 }
4619 // Sets up an effective RHS
4620 void
useEffectiveRhs(ClpSimplex * model)4621 ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
4622 {
4623 delete [] rhsOffset_;
4624 int numberRows = model->numberRows();
4625 rhsOffset_ = new double[numberRows];
4626 rhsOffset(model, true);
4627 }
4628 // Gets rid of special copies
4629 void
clearCopies()4630 ClpPackedMatrix::clearCopies()
4631 {
4632 delete rowCopy_;
4633 delete columnCopy_;
4634 rowCopy_ = NULL;
4635 columnCopy_ = NULL;
4636 flags_ &= ~(4 + 8);
4637 checkGaps();
4638 }
4639 // makes sure active columns correct
4640 int
refresh(ClpSimplex *)4641 ClpPackedMatrix::refresh(ClpSimplex * )
4642 {
4643 numberActiveColumns_ = matrix_->getNumCols();
4644 #if 0
4645 ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
4646 ClpPackedMatrix* rowCopy =
4647 dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4648 // Make sure it is really a ClpPackedMatrix
4649 assert (rowCopy != NULL);
4650
4651 const int * column = rowCopy->matrix_->getIndices();
4652 const CoinBigIndex * rowStart = rowCopy->matrix_->getVectorStarts();
4653 const int * rowLength = rowCopy->matrix_->getVectorLengths();
4654 const double * element = rowCopy->matrix_->getElements();
4655 int numberRows = rowCopy->matrix_->getNumRows();
4656 for (int i = 0; i < numberRows; i++) {
4657 if (!rowLength[i])
4658 printf("zero row %d\n", i);
4659 }
4660 delete rowCopy;
4661 #endif
4662 return 0;
4663 }
4664
4665 /* Scales rowCopy if column copy scaled
4666 Only called if scales already exist */
4667 void
scaleRowCopy(ClpModel * model) const4668 ClpPackedMatrix::scaleRowCopy(ClpModel * model) const
4669 {
4670 if (model->rowCopy()) {
4671 // need to replace row by row
4672 int numberRows = model->numberRows();
4673 #ifndef NDEBUG
4674 int numberColumns = matrix_->getNumCols();
4675 #endif
4676 ClpMatrixBase * rowCopyBase = model->rowCopy();
4677 #ifndef NDEBUG
4678 ClpPackedMatrix* rowCopy =
4679 dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4680 // Make sure it is really a ClpPackedMatrix
4681 assert (rowCopy != NULL);
4682 #else
4683 ClpPackedMatrix* rowCopy =
4684 static_cast< ClpPackedMatrix*>(rowCopyBase);
4685 #endif
4686
4687 const int * column = rowCopy->getIndices();
4688 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4689 double * element = rowCopy->getMutableElements();
4690 const double * rowScale = model->rowScale();
4691 const double * columnScale = model->columnScale();
4692 // scale row copy
4693 for (int iRow = 0; iRow < numberRows; iRow++) {
4694 CoinBigIndex j;
4695 double scale = rowScale[iRow];
4696 double * elementsInThisRow = element + rowStart[iRow];
4697 const int * columnsInThisRow = column + rowStart[iRow];
4698 int number = rowStart[iRow+1] - rowStart[iRow];
4699 assert (number <= numberColumns);
4700 for (j = 0; j < number; j++) {
4701 int iColumn = columnsInThisRow[j];
4702 elementsInThisRow[j] *= scale * columnScale[iColumn];
4703 }
4704 }
4705 }
4706 }
4707 /* Realy really scales column copy
4708 Only called if scales already exist.
4709 Up to user ro delete */
4710 ClpMatrixBase *
scaledColumnCopy(ClpModel * model) const4711 ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const
4712 {
4713 // need to replace column by column
4714 #ifndef NDEBUG
4715 int numberRows = model->numberRows();
4716 #endif
4717 int numberColumns = matrix_->getNumCols();
4718 ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
4719 const int * row = copy->getIndices();
4720 const CoinBigIndex * columnStart = copy->getVectorStarts();
4721 const int * length = copy->getVectorLengths();
4722 double * element = copy->getMutableElements();
4723 const double * rowScale = model->rowScale();
4724 const double * columnScale = model->columnScale();
4725 // scale column copy
4726 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4727 CoinBigIndex j;
4728 double scale = columnScale[iColumn];
4729 double * elementsInThisColumn = element + columnStart[iColumn];
4730 const int * rowsInThisColumn = row + columnStart[iColumn];
4731 int number = length[iColumn];
4732 assert (number <= numberRows);
4733 for (j = 0; j < number; j++) {
4734 int iRow = rowsInThisColumn[j];
4735 elementsInThisColumn[j] *= scale * rowScale[iRow];
4736 }
4737 }
4738 return copy;
4739 }
4740 // Really scale matrix
4741 void
reallyScale(const double * rowScale,const double * columnScale)4742 ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
4743 {
4744 clearCopies();
4745 int numberColumns = matrix_->getNumCols();
4746 const int * row = matrix_->getIndices();
4747 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4748 const int * length = matrix_->getVectorLengths();
4749 double * element = matrix_->getMutableElements();
4750 // scale
4751 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4752 CoinBigIndex j;
4753 double scale = columnScale[iColumn];
4754 for (j = columnStart[iColumn]; j < columnStart[iColumn] + length[iColumn]; j++) {
4755 int iRow = row[j];
4756 element[j] *= scale * rowScale[iRow];
4757 }
4758 }
4759 }
4760 /* Delete the columns whose indices are listed in <code>indDel</code>. */
4761 void
deleteCols(const int numDel,const int * indDel)4762 ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
4763 {
4764 if (matrix_->getNumCols())
4765 matrix_->deleteCols(numDel, indDel);
4766 clearCopies();
4767 numberActiveColumns_ = matrix_->getNumCols();
4768 // may now have gaps
4769 checkGaps();
4770 matrix_->setExtraGap(0.0);
4771 }
4772 /* Delete the rows whose indices are listed in <code>indDel</code>. */
4773 void
deleteRows(const int numDel,const int * indDel)4774 ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
4775 {
4776 if (matrix_->getNumRows())
4777 matrix_->deleteRows(numDel, indDel);
4778 clearCopies();
4779 numberActiveColumns_ = matrix_->getNumCols();
4780 // may now have gaps
4781 checkGaps();
4782 matrix_->setExtraGap(0.0);
4783 }
4784 #ifndef CLP_NO_VECTOR
4785 // Append Columns
4786 void
appendCols(int number,const CoinPackedVectorBase * const * columns)4787 ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
4788 {
4789 matrix_->appendCols(number, columns);
4790 numberActiveColumns_ = matrix_->getNumCols();
4791 clearCopies();
4792 }
4793 // Append Rows
4794 void
appendRows(int number,const CoinPackedVectorBase * const * rows)4795 ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
4796 {
4797 matrix_->appendRows(number, rows);
4798 numberActiveColumns_ = matrix_->getNumCols();
4799 // may now have gaps
4800 checkGaps();
4801 clearCopies();
4802 }
4803 #endif
4804 /* Set the dimensions of the matrix. In effect, append new empty
4805 columns/rows to the matrix. A negative number for either dimension
4806 means that that dimension doesn't change. Otherwise the new dimensions
4807 MUST be at least as large as the current ones otherwise an exception
4808 is thrown. */
4809 void
setDimensions(int numrows,int numcols)4810 ClpPackedMatrix::setDimensions(int numrows, int numcols)
4811 {
4812 matrix_->setDimensions(numrows, numcols);
4813 }
4814 /* Append a set of rows/columns to the end of the matrix. Returns number of errors
4815 i.e. if any of the new rows/columns contain an index that's larger than the
4816 number of columns-1/rows-1 (if numberOther>0) or duplicates
4817 If 0 then rows, 1 if columns */
4818 int
appendMatrix(int number,int type,const CoinBigIndex * starts,const int * index,const double * element,int numberOther)4819 ClpPackedMatrix::appendMatrix(int number, int type,
4820 const CoinBigIndex * starts, const int * index,
4821 const double * element, int numberOther)
4822 {
4823 int numberErrors = 0;
4824 // make sure other dimension is big enough
4825 if (type == 0) {
4826 // rows
4827 if (matrix_->isColOrdered() && numberOther > matrix_->getNumCols())
4828 matrix_->setDimensions(-1, numberOther);
4829 if (!matrix_->isColOrdered() || numberOther >= 0 || matrix_->getExtraGap() || matrix_->hasGaps()) {
4830 numberErrors = matrix_->appendRows(number, starts, index, element, numberOther);
4831 } else {
4832 //CoinPackedMatrix mm(*matrix_);
4833 matrix_->appendMinorFast(number, starts, index, element);
4834 //mm.appendRows(number,starts,index,element,numberOther);
4835 //if (!mm.isEquivalent(*matrix_)) {
4836 //printf("bad append\n");
4837 //abort();
4838 //}
4839 }
4840 } else {
4841 // columns
4842 if (!matrix_->isColOrdered() && numberOther > matrix_->getNumRows())
4843 matrix_->setDimensions(numberOther, -1);
4844 numberErrors = matrix_->appendCols(number, starts, index, element, numberOther);
4845 }
4846 clearCopies();
4847 numberActiveColumns_ = matrix_->getNumCols();
4848 return numberErrors;
4849 }
4850 void
specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)4851 ClpPackedMatrix::specialRowCopy(ClpSimplex * model, const ClpMatrixBase * rowCopy)
4852 {
4853 delete rowCopy_;
4854 rowCopy_ = new ClpPackedMatrix2(model, rowCopy->getPackedMatrix());
4855 // See if anything in it
4856 if (!rowCopy_->usefulInfo()) {
4857 delete rowCopy_;
4858 rowCopy_ = NULL;
4859 flags_ &= ~4;
4860 } else {
4861 flags_ |= 4;
4862 }
4863 }
4864 void
specialColumnCopy(ClpSimplex * model)4865 ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
4866 {
4867 delete columnCopy_;
4868 if ((flags_ & 16) != 0) {
4869 columnCopy_ = new ClpPackedMatrix3(model, matrix_);
4870 flags_ |= 8;
4871 } else {
4872 columnCopy_ = NULL;
4873 }
4874 }
4875 // Say we don't want special column copy
4876 void
releaseSpecialColumnCopy()4877 ClpPackedMatrix::releaseSpecialColumnCopy()
4878 {
4879 flags_ &= ~(8 + 16);
4880 delete columnCopy_;
4881 columnCopy_ = NULL;
4882 }
4883 // Correct sequence in and out to give true value
4884 void
correctSequence(const ClpSimplex * model,int & sequenceIn,int & sequenceOut)4885 ClpPackedMatrix::correctSequence(const ClpSimplex * model, int & sequenceIn, int & sequenceOut)
4886 {
4887 if (columnCopy_) {
4888 if (sequenceIn != -999) {
4889 if (sequenceIn != sequenceOut) {
4890 if (sequenceIn < numberActiveColumns_)
4891 columnCopy_->swapOne(model, this, sequenceIn);
4892 if (sequenceOut < numberActiveColumns_)
4893 columnCopy_->swapOne(model, this, sequenceOut);
4894 }
4895 } else {
4896 // do all
4897 columnCopy_->sortBlocks(model);
4898 }
4899 }
4900 }
4901 // Check validity
4902 void
checkFlags(int type) const4903 ClpPackedMatrix::checkFlags(int type) const
4904 {
4905 int iColumn;
4906 // get matrix data pointers
4907 //const int * row = matrix_->getIndices();
4908 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4909 const int * columnLength = matrix_->getVectorLengths();
4910 const double * elementByColumn = matrix_->getElements();
4911 if (!zeros()) {
4912 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4913 CoinBigIndex j;
4914 for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4915 if (!elementByColumn[j])
4916 abort();
4917 }
4918 }
4919 }
4920 if ((flags_ & 2) == 0) {
4921 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4922 if (columnStart[iColumn+1] != columnStart[iColumn] + columnLength[iColumn]) {
4923 abort();
4924 }
4925 }
4926 }
4927 if (type) {
4928 if ((flags_ & 2) != 0) {
4929 bool ok = true;
4930 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4931 if (columnStart[iColumn+1] != columnStart[iColumn] + columnLength[iColumn]) {
4932 ok = false;
4933 break;
4934 }
4935 }
4936 if (ok)
4937 COIN_DETAIL_PRINT(printf("flags_ could be 0\n"));
4938 }
4939 }
4940 }
4941 //#############################################################################
4942 // Constructors / Destructor / Assignment
4943 //#############################################################################
4944
4945 //-------------------------------------------------------------------
4946 // Default Constructor
4947 //-------------------------------------------------------------------
ClpPackedMatrix2()4948 ClpPackedMatrix2::ClpPackedMatrix2 ()
4949 : numberBlocks_(0),
4950 numberRows_(0),
4951 offset_(NULL),
4952 count_(NULL),
4953 rowStart_(NULL),
4954 column_(NULL),
4955 work_(NULL)
4956 {
4957 #ifdef THREAD
4958 threadId_ = NULL;
4959 info_ = NULL;
4960 #endif
4961 }
4962 //-------------------------------------------------------------------
4963 // Useful Constructor
4964 //-------------------------------------------------------------------
ClpPackedMatrix2(ClpSimplex *,const CoinPackedMatrix * rowCopy)4965 ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * , const CoinPackedMatrix * rowCopy)
4966 : numberBlocks_(0),
4967 numberRows_(0),
4968 offset_(NULL),
4969 count_(NULL),
4970 rowStart_(NULL),
4971 column_(NULL),
4972 work_(NULL)
4973 {
4974 #ifdef THREAD
4975 threadId_ = NULL;
4976 info_ = NULL;
4977 #endif
4978 numberRows_ = rowCopy->getNumRows();
4979 if (!numberRows_)
4980 return;
4981 int numberColumns = rowCopy->getNumCols();
4982 const int * column = rowCopy->getIndices();
4983 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4984 const int * length = rowCopy->getVectorLengths();
4985 const double * element = rowCopy->getElements();
4986 int chunk = 32768; // tune
4987 //chunk=100;
4988 // tune
4989 #if 0
4990 int chunkY[7] = {1024, 2048, 4096, 8192, 16384, 32768, 65535};
4991 int its = model->maximumIterations();
4992 if (its >= 1000000 && its < 1000999) {
4993 its -= 1000000;
4994 its = its / 10;
4995 if (its >= 7) abort();
4996 chunk = chunkY[its];
4997 printf("chunk size %d\n", chunk);
4998 #define cpuid(func,ax,bx,cx,dx)\
4999 __asm__ __volatile__ ("cpuid":\
5000 "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
5001 unsigned int a, b, c, d;
5002 int func = 0;
5003 cpuid(func, a, b, c, d);
5004 {
5005 int i;
5006 unsigned int value;
5007 value = b;
5008 for (i = 0; i < 4; i++) {
5009 printf("%c", (value & 0xff));
5010 value = value >> 8;
5011 }
5012 value = d;
5013 for (i = 0; i < 4; i++) {
5014 printf("%c", (value & 0xff));
5015 value = value >> 8;
5016 }
5017 value = c;
5018 for (i = 0; i < 4; i++) {
5019 printf("%c", (value & 0xff));
5020 value = value >> 8;
5021 }
5022 printf("\n");
5023 int maxfunc = a;
5024 if (maxfunc > 10) {
5025 printf("not intel?\n");
5026 abort();
5027 }
5028 for (func = 1; func <= maxfunc; func++) {
5029 cpuid(func, a, b, c, d);
5030 printf("func %d, %x %x %x %x\n", func, a, b, c, d);
5031 }
5032 }
5033 #else
5034 if (numberColumns > 10000 || chunk == 100) {
5035 #endif
5036 } else {
5037 //printf("no chunk\n");
5038 return;
5039 }
5040 // Could also analyze matrix to get natural breaks
5041 numberBlocks_ = (numberColumns + chunk - 1) / chunk;
5042 #ifdef THREAD
5043 // Get work areas
5044 threadId_ = new pthread_t [numberBlocks_];
5045 info_ = new dualColumn0Struct[numberBlocks_];
5046 #endif
5047 // Even out
5048 chunk = (numberColumns + numberBlocks_ - 1) / numberBlocks_;
5049 offset_ = new int[numberBlocks_+1];
5050 offset_[numberBlocks_] = numberColumns;
5051 int nRow = numberBlocks_ * numberRows_;
5052 count_ = new unsigned short[nRow];
5053 memset(count_, 0, nRow * sizeof(unsigned short));
5054 rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
5055 CoinBigIndex nElement = rowStart[numberRows_];
5056 rowStart_[nRow+numberRows_] = nElement;
5057 column_ = new unsigned short [nElement];
5058 // assumes int <= double
5059 int sizeWork = 6 * numberBlocks_;
5060 work_ = new double[sizeWork];
5061 int iBlock;
5062 int nZero = 0;
5063 for (iBlock = 0; iBlock < numberBlocks_; iBlock++) {
5064 int start = iBlock * chunk;
5065 offset_[iBlock] = start;
5066 int end = start + chunk;
5067 for (int iRow = 0; iRow < numberRows_; iRow++) {
5068 if (rowStart[iRow+1] != rowStart[iRow] + length[iRow]) {
5069 printf("not packed correctly - gaps\n");
5070 abort();
5071 }
5072 bool lastFound = false;
5073 int nFound = 0;
5074 for (CoinBigIndex j = rowStart[iRow];
5075 j < rowStart[iRow] + length[iRow]; j++) {
5076 int iColumn = column[j];
5077 if (iColumn >= start) {
5078 if (iColumn < end) {
5079 if (!element[j]) {
5080 printf("not packed correctly - zero element\n");
5081 abort();
5082 }
5083 column_[j] = static_cast<unsigned short>(iColumn - start);
5084 nFound++;
5085 if (lastFound) {
5086 printf("not packed correctly - out of order\n");
5087 abort();
5088 }
5089 } else {
5090 //can't find any more
5091 lastFound = true;
5092 }
5093 }
5094 }
5095 count_[iRow*numberBlocks_+iBlock] = static_cast<unsigned short>(nFound);
5096 if (!nFound)
5097 nZero++;
5098 }
5099 }
5100 //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
5101 //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
5102 }
5103
5104 //-------------------------------------------------------------------
5105 // Copy constructor
5106 //-------------------------------------------------------------------
5107 ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs)
5108 : numberBlocks_(rhs.numberBlocks_),
5109 numberRows_(rhs.numberRows_)
5110 {
5111 if (numberBlocks_) {
5112 offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + 1);
5113 int nRow = numberBlocks_ * numberRows_;
5114 count_ = CoinCopyOfArray(rhs.count_, nRow);
5115 rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + 1);
5116 CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5117 column_ = CoinCopyOfArray(rhs.column_, nElement);
5118 int sizeWork = 6 * numberBlocks_;
5119 work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5120 #ifdef THREAD
5121 threadId_ = new pthread_t [numberBlocks_];
5122 info_ = new dualColumn0Struct[numberBlocks_];
5123 #endif
5124 } else {
5125 offset_ = NULL;
5126 count_ = NULL;
5127 rowStart_ = NULL;
5128 column_ = NULL;
5129 work_ = NULL;
5130 #ifdef THREAD
5131 threadId_ = NULL;
5132 info_ = NULL;
5133 #endif
5134 }
5135 }
5136 //-------------------------------------------------------------------
5137 // Destructor
5138 //-------------------------------------------------------------------
5139 ClpPackedMatrix2::~ClpPackedMatrix2 ()
5140 {
5141 delete [] offset_;
5142 delete [] count_;
5143 delete [] rowStart_;
5144 delete [] column_;
5145 delete [] work_;
5146 #ifdef THREAD
5147 delete [] threadId_;
5148 delete [] info_;
5149 #endif
5150 }
5151
5152 //----------------------------------------------------------------
5153 // Assignment operator
5154 //-------------------------------------------------------------------
5155 ClpPackedMatrix2 &
5156 ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
5157 {
5158 if (this != &rhs) {
5159 numberBlocks_ = rhs.numberBlocks_;
5160 numberRows_ = rhs.numberRows_;
5161 delete [] offset_;
5162 delete [] count_;
5163 delete [] rowStart_;
5164 delete [] column_;
5165 delete [] work_;
5166 #ifdef THREAD
5167 delete [] threadId_;
5168 delete [] info_;
5169 #endif
5170 if (numberBlocks_) {
5171 offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + 1);
5172 int nRow = numberBlocks_ * numberRows_;
5173 count_ = CoinCopyOfArray(rhs.count_, nRow);
5174 rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + 1);
5175 CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5176 column_ = CoinCopyOfArray(rhs.column_, nElement);
5177 int sizeWork = 6 * numberBlocks_;
5178 work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5179 #ifdef THREAD
5180 threadId_ = new pthread_t [numberBlocks_];
5181 info_ = new dualColumn0Struct[numberBlocks_];
5182 #endif
5183 } else {
5184 offset_ = NULL;
5185 count_ = NULL;
5186 rowStart_ = NULL;
5187 column_ = NULL;
5188 work_ = NULL;
5189 #ifdef THREAD
5190 threadId_ = NULL;
5191 info_ = NULL;
5192 #endif
5193 }
5194 }
5195 return *this;
5196 }
5197 static int dualColumn0(const ClpSimplex * model, double * spare,
5198 int * spareIndex, const double * arrayTemp,
5199 const int * indexTemp, int numberIn,
5200 int offset, double acceptablePivot, double * bestPossiblePtr,
5201 double * upperThetaPtr, int * posFreePtr, double * freePivotPtr)
5202 {
5203 // do dualColumn0
5204 int i;
5205 int numberRemaining = 0;
5206 double bestPossible = 0.0;
5207 double upperTheta = 1.0e31;
5208 double freePivot = acceptablePivot;
5209 int posFree = -1;
5210 const double * reducedCost = model->djRegion(1);
5211 double dualTolerance = model->dualTolerance();
5212 // We can also see if infeasible or pivoting on free
5213 double tentativeTheta = 1.0e25;
5214 for (i = 0; i < numberIn; i++) {
5215 double alpha = arrayTemp[i];
5216 int iSequence = indexTemp[i] + offset;
5217 double oldValue;
5218 double value;
5219 bool keep;
5220
5221 switch(model->getStatus(iSequence)) {
5222
5223 case ClpSimplex::basic:
5224 case ClpSimplex::isFixed:
5225 break;
5226 case ClpSimplex::isFree:
5227 case ClpSimplex::superBasic:
5228 bestPossible = CoinMax(bestPossible, fabs(alpha));
5229 oldValue = reducedCost[iSequence];
5230 // If free has to be very large - should come in via dualRow
5231 if (model->getStatus(iSequence) == ClpSimplex::isFree && fabs(alpha) < 1.0e-3)
5232 break;
5233 if (oldValue > dualTolerance) {
5234 keep = true;
5235 } else if (oldValue < -dualTolerance) {
5236 keep = true;
5237 } else {
5238 if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5))
5239 keep = true;
5240 else
5241 keep = false;
5242 }
5243 if (keep) {
5244 // free - choose largest
5245 if (fabs(alpha) > freePivot) {
5246 freePivot = fabs(alpha);
5247 posFree = i;
5248 }
5249 }
5250 break;
5251 case ClpSimplex::atUpperBound:
5252 oldValue = reducedCost[iSequence];
5253 value = oldValue - tentativeTheta * alpha;
5254 //assert (oldValue<=dualTolerance*1.0001);
5255 if (value > dualTolerance) {
5256 bestPossible = CoinMax(bestPossible, -alpha);
5257 value = oldValue - upperTheta * alpha;
5258 if (value > dualTolerance && -alpha >= acceptablePivot)
5259 upperTheta = (oldValue - dualTolerance) / alpha;
5260 // add to list
5261 spare[numberRemaining] = alpha;
5262 spareIndex[numberRemaining++] = iSequence;
5263 }
5264 break;
5265 case ClpSimplex::atLowerBound:
5266 oldValue = reducedCost[iSequence];
5267 value = oldValue - tentativeTheta * alpha;
5268 //assert (oldValue>=-dualTolerance*1.0001);
5269 if (value < -dualTolerance) {
5270 bestPossible = CoinMax(bestPossible, alpha);
5271 value = oldValue - upperTheta * alpha;
5272 if (value < -dualTolerance && alpha >= acceptablePivot)
5273 upperTheta = (oldValue + dualTolerance) / alpha;
5274 // add to list
5275 spare[numberRemaining] = alpha;
5276 spareIndex[numberRemaining++] = iSequence;
5277 }
5278 break;
5279 }
5280 }
5281 *bestPossiblePtr = bestPossible;
5282 *upperThetaPtr = upperTheta;
5283 *freePivotPtr = freePivot;
5284 *posFreePtr = posFree;
5285 return numberRemaining;
5286 }
5287 static int doOneBlock(double * array, int * index,
5288 const double * pi, const CoinBigIndex * rowStart, const double * element,
5289 const unsigned short * column, int numberInRowArray, int numberLook)
5290 {
5291 int iWhich = 0;
5292 int nextN = 0;
5293 CoinBigIndex nextStart = 0;
5294 double nextPi = 0.0;
5295 for (; iWhich < numberInRowArray; iWhich++) {
5296 nextStart = rowStart[0];
5297 nextN = rowStart[numberInRowArray] - nextStart;
5298 rowStart++;
5299 if (nextN) {
5300 nextPi = pi[iWhich];
5301 break;
5302 }
5303 }
5304 int i;
5305 while (iWhich < numberInRowArray) {
5306 double value = nextPi;
5307 CoinBigIndex j = nextStart;
5308 int n = nextN;
5309 // get next
5310 iWhich++;
5311 for (; iWhich < numberInRowArray; iWhich++) {
5312 nextStart = rowStart[0];
5313 nextN = rowStart[numberInRowArray] - nextStart;
5314 rowStart++;
5315 if (nextN) {
5316 //coin_prefetch_const(element + nextStart);
5317 nextPi = pi[iWhich];
5318 break;
5319 }
5320 }
5321 CoinBigIndex end = j + n;
5322 //coin_prefetch_const(element+rowStart_[i+1]);
5323 //coin_prefetch_const(column_+rowStart_[i+1]);
5324 if (n < 100) {
5325 if ((n & 1) != 0) {
5326 unsigned int jColumn = column[j];
5327 array[jColumn] -= value * element[j];
5328 j++;
5329 }
5330 //coin_prefetch_const(column + nextStart);
5331 for (; j < end; j += 2) {
5332 unsigned int jColumn0 = column[j];
5333 double value0 = value * element[j];
5334 unsigned int jColumn1 = column[j+1];
5335 double value1 = value * element[j+1];
5336 array[jColumn0] -= value0;
5337 array[jColumn1] -= value1;
5338 }
5339 } else {
5340 if ((n & 1) != 0) {
5341 unsigned int jColumn = column[j];
5342 array[jColumn] -= value * element[j];
5343 j++;
5344 }
5345 if ((n & 2) != 0) {
5346 unsigned int jColumn0 = column[j];
5347 double value0 = value * element[j];
5348 unsigned int jColumn1 = column[j+1];
5349 double value1 = value * element[j+1];
5350 array[jColumn0] -= value0;
5351 array[jColumn1] -= value1;
5352 j += 2;
5353 }
5354 if ((n & 4) != 0) {
5355 unsigned int jColumn0 = column[j];
5356 double value0 = value * element[j];
5357 unsigned int jColumn1 = column[j+1];
5358 double value1 = value * element[j+1];
5359 unsigned int jColumn2 = column[j+2];
5360 double value2 = value * element[j+2];
5361 unsigned int jColumn3 = column[j+3];
5362 double value3 = value * element[j+3];
5363 array[jColumn0] -= value0;
5364 array[jColumn1] -= value1;
5365 array[jColumn2] -= value2;
5366 array[jColumn3] -= value3;
5367 j += 4;
5368 }
5369 //coin_prefetch_const(column+nextStart);
5370 for (; j < end; j += 8) {
5371 //coin_prefetch_const(element + j + 16);
5372 unsigned int jColumn0 = column[j];
5373 double value0 = value * element[j];
5374 unsigned int jColumn1 = column[j+1];
5375 double value1 = value * element[j+1];
5376 unsigned int jColumn2 = column[j+2];
5377 double value2 = value * element[j+2];
5378 unsigned int jColumn3 = column[j+3];
5379 double value3 = value * element[j+3];
5380 array[jColumn0] -= value0;
5381 array[jColumn1] -= value1;
5382 array[jColumn2] -= value2;
5383 array[jColumn3] -= value3;
5384 //coin_prefetch_const(column + j + 16);
5385 jColumn0 = column[j+4];
5386 value0 = value * element[j+4];
5387 jColumn1 = column[j+5];
5388 value1 = value * element[j+5];
5389 jColumn2 = column[j+6];
5390 value2 = value * element[j+6];
5391 jColumn3 = column[j+7];
5392 value3 = value * element[j+7];
5393 array[jColumn0] -= value0;
5394 array[jColumn1] -= value1;
5395 array[jColumn2] -= value2;
5396 array[jColumn3] -= value3;
5397 }
5398 }
5399 }
5400 // get rid of tiny values
5401 int nSmall = numberLook;
5402 int numberNonZero = 0;
5403 for (i = 0; i < nSmall; i++) {
5404 double value = array[i];
5405 array[i] = 0.0;
5406 if (fabs(value) > 1.0e-12) {
5407 array[numberNonZero] = value;
5408 index[numberNonZero++] = i;
5409 }
5410 }
5411 for (; i < numberLook; i += 4) {
5412 double value0 = array[i+0];
5413 double value1 = array[i+1];
5414 double value2 = array[i+2];
5415 double value3 = array[i+3];
5416 array[i+0] = 0.0;
5417 array[i+1] = 0.0;
5418 array[i+2] = 0.0;
5419 array[i+3] = 0.0;
5420 if (fabs(value0) > 1.0e-12) {
5421 array[numberNonZero] = value0;
5422 index[numberNonZero++] = i + 0;
5423 }
5424 if (fabs(value1) > 1.0e-12) {
5425 array[numberNonZero] = value1;
5426 index[numberNonZero++] = i + 1;
5427 }
5428 if (fabs(value2) > 1.0e-12) {
5429 array[numberNonZero] = value2;
5430 index[numberNonZero++] = i + 2;
5431 }
5432 if (fabs(value3) > 1.0e-12) {
5433 array[numberNonZero] = value3;
5434 index[numberNonZero++] = i + 3;
5435 }
5436 }
5437 return numberNonZero;
5438 }
5439 #ifdef THREAD
5440 static void * doOneBlockThread(void * voidInfo)
5441 {
5442 dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
5443 *(info->numberInPtr) = doOneBlock(info->arrayTemp, info->indexTemp, info->pi,
5444 info->rowStart, info->element, info->column,
5445 info->numberInRowArray, info->numberLook);
5446 return NULL;
5447 }
5448 static void * doOneBlockAnd0Thread(void * voidInfo)
5449 {
5450 dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
5451 *(info->numberInPtr) = doOneBlock(info->arrayTemp, info->indexTemp, info->pi,
5452 info->rowStart, info->element, info->column,
5453 info->numberInRowArray, info->numberLook);
5454 *(info->numberOutPtr) = dualColumn0(info->model, info->spare,
5455 info->spareIndex, (const double *)info->arrayTemp,
5456 (const int *) info->indexTemp, *(info->numberInPtr),
5457 info->offset, info->acceptablePivot, info->bestPossiblePtr,
5458 info->upperThetaPtr, info->posFreePtr, info->freePivotPtr);
5459 return NULL;
5460 }
5461 #endif
5462 /* Return <code>x * scalar * A in <code>z</code>.
5463 Note - x packed and z will be packed mode
5464 Squashes small elements and knows about ClpSimplex */
5465 void
5466 ClpPackedMatrix2::transposeTimes(const ClpSimplex * model,
5467 const CoinPackedMatrix * rowCopy,
5468 const CoinIndexedVector * rowArray,
5469 CoinIndexedVector * spareArray,
5470 CoinIndexedVector * columnArray) const
5471 {
5472 // See if dualColumn0 coding wanted
5473 bool dualColumn = model->spareIntArray_[0] == 1;
5474 double acceptablePivot = model->spareDoubleArray_[0];
5475 double bestPossible = 0.0;
5476 double upperTheta = 1.0e31;
5477 double freePivot = acceptablePivot;
5478 int posFree = -1;
5479 int numberRemaining = 0;
5480 //if (model->numberIterations()>=200000) {
5481 //printf("time %g\n",CoinCpuTime()-startTime);
5482 //exit(77);
5483 //}
5484 double * pi = rowArray->denseVector();
5485 int numberNonZero = 0;
5486 int * index = columnArray->getIndices();
5487 double * array = columnArray->denseVector();
5488 int numberInRowArray = rowArray->getNumElements();
5489 const int * whichRow = rowArray->getIndices();
5490 double * element = const_cast<double *>(rowCopy->getElements());
5491 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
5492 int i;
5493 CoinBigIndex * rowStart2 = rowStart_;
5494 if (!dualColumn) {
5495 for (i = 0; i < numberInRowArray; i++) {
5496 int iRow = whichRow[i];
5497 CoinBigIndex start = rowStart[iRow];
5498 *rowStart2 = start;
5499 unsigned short * count1 = count_ + iRow * numberBlocks_;
5500 int put = 0;
5501 for (int j = 0; j < numberBlocks_; j++) {
5502 put += numberInRowArray;
5503 start += count1[j];
5504 rowStart2[put] = start;
5505 }
5506 rowStart2 ++;
5507 }
5508 } else {
5509 // also do dualColumn stuff
5510 double * spare = spareArray->denseVector();
5511 int * spareIndex = spareArray->getIndices();
5512 const double * reducedCost = model->djRegion(0);
5513 double dualTolerance = model->dualTolerance();
5514 // We can also see if infeasible or pivoting on free
5515 double tentativeTheta = 1.0e25;
5516 int addSequence = model->numberColumns();
5517 for (i = 0; i < numberInRowArray; i++) {
5518 int iRow = whichRow[i];
5519 double alpha = pi[i];
5520 double oldValue;
5521 double value;
5522 bool keep;
5523
5524 switch(model->getStatus(iRow + addSequence)) {
5525
5526 case ClpSimplex::basic:
5527 case ClpSimplex::isFixed:
5528 break;
5529 case ClpSimplex::isFree:
5530 case ClpSimplex::superBasic:
5531 bestPossible = CoinMax(bestPossible, fabs(alpha));
5532 oldValue = reducedCost[iRow];
5533 // If free has to be very large - should come in via dualRow
5534 if (model->getStatus(iRow + addSequence) == ClpSimplex::isFree && fabs(alpha) < 1.0e-3)
5535 break;
5536 if (oldValue > dualTolerance) {
5537 keep = true;
5538 } else if (oldValue < -dualTolerance) {
5539 keep = true;
5540 } else {
5541 if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5))
5542 keep = true;
5543 else
5544 keep = false;
5545 }
5546 if (keep) {
5547 // free - choose largest
5548 if (fabs(alpha) > freePivot) {
5549 freePivot = fabs(alpha);
5550 posFree = i + addSequence;
5551 }
5552 }
5553 break;
5554 case ClpSimplex::atUpperBound:
5555 oldValue = reducedCost[iRow];
5556 value = oldValue - tentativeTheta * alpha;
5557 //assert (oldValue<=dualTolerance*1.0001);
5558 if (value > dualTolerance) {
5559 bestPossible = CoinMax(bestPossible, -alpha);
5560 value = oldValue - upperTheta * alpha;
5561 if (value > dualTolerance && -alpha >= acceptablePivot)
5562 upperTheta = (oldValue - dualTolerance) / alpha;
5563 // add to list
5564 spare[numberRemaining] = alpha;
5565 spareIndex[numberRemaining++] = iRow + addSequence;
5566 }
5567 break;
5568 case ClpSimplex::atLowerBound:
5569 oldValue = reducedCost[iRow];
5570 value = oldValue - tentativeTheta * alpha;
5571 //assert (oldValue>=-dualTolerance*1.0001);
5572 if (value < -dualTolerance) {
5573 bestPossible = CoinMax(bestPossible, alpha);
5574 value = oldValue - upperTheta * alpha;
5575 if (value < -dualTolerance && alpha >= acceptablePivot)
5576 upperTheta = (oldValue + dualTolerance) / alpha;
5577 // add to list
5578 spare[numberRemaining] = alpha;
5579 spareIndex[numberRemaining++] = iRow + addSequence;
5580 }
5581 break;
5582 }
5583 CoinBigIndex start = rowStart[iRow];
5584 *rowStart2 = start;
5585 unsigned short * count1 = count_ + iRow * numberBlocks_;
5586 int put = 0;
5587 for (int j = 0; j < numberBlocks_; j++) {
5588 put += numberInRowArray;
5589 start += count1[j];
5590 rowStart2[put] = start;
5591 }
5592 rowStart2 ++;
5593 }
5594 }
5595 double * spare = spareArray->denseVector();
5596 int * spareIndex = spareArray->getIndices();
5597 int saveNumberRemaining = numberRemaining;
5598 int iBlock;
5599 for (iBlock = 0; iBlock < numberBlocks_; iBlock++) {
5600 double * dwork = work_ + 6 * iBlock;
5601 int * iwork = reinterpret_cast<int *> (dwork + 3);
5602 if (!dualColumn) {
5603 #ifndef THREAD
5604 int offset = offset_[iBlock];
5605 int offset3 = offset;
5606 offset = numberNonZero;
5607 double * arrayTemp = array + offset;
5608 int * indexTemp = index + offset;
5609 iwork[0] = doOneBlock(arrayTemp, indexTemp, pi, rowStart_ + numberInRowArray * iBlock,
5610 element, column_, numberInRowArray, offset_[iBlock+1] - offset);
5611 int number = iwork[0];
5612 for (i = 0; i < number; i++) {
5613 //double value = arrayTemp[i];
5614 //arrayTemp[i]=0.0;
5615 //array[numberNonZero]=value;
5616 index[numberNonZero++] = indexTemp[i] + offset3;
5617 }
5618 #else
5619 int offset = offset_[iBlock];
5620 double * arrayTemp = array + offset;
5621 int * indexTemp = index + offset;
5622 dualColumn0Struct * infoPtr = info_ + iBlock;
5623 infoPtr->arrayTemp = arrayTemp;
5624 infoPtr->indexTemp = indexTemp;
5625 infoPtr->numberInPtr = &iwork[0];
5626 infoPtr->pi = pi;
5627 infoPtr->rowStart = rowStart_ + numberInRowArray * iBlock;
5628 infoPtr->element = element;
5629 infoPtr->column = column_;
5630 infoPtr->numberInRowArray = numberInRowArray;
5631 infoPtr->numberLook = offset_[iBlock+1] - offset;
5632 pthread_create(&threadId_[iBlock], NULL, doOneBlockThread, infoPtr);
5633 #endif
5634 } else {
5635 #ifndef THREAD
5636 int offset = offset_[iBlock];
5637 // allow for already saved
5638 int offset2 = offset + saveNumberRemaining;
5639 int offset3 = offset;
5640 offset = numberNonZero;
5641 offset2 = numberRemaining;
5642 double * arrayTemp = array + offset;
5643 int * indexTemp = index + offset;
5644 iwork[0] = doOneBlock(arrayTemp, indexTemp, pi, rowStart_ + numberInRowArray * iBlock,
5645 element, column_, numberInRowArray, offset_[iBlock+1] - offset);
5646 iwork[1] = dualColumn0(model, spare + offset2,
5647 spareIndex + offset2,
5648 arrayTemp, indexTemp,
5649 iwork[0], offset3, acceptablePivot,
5650 &dwork[0], &dwork[1], &iwork[2],
5651 &dwork[2]);
5652 int number = iwork[0];
5653 int numberLook = iwork[1];
5654 #if 1
5655 numberRemaining += numberLook;
5656 #else
5657 double * spareTemp = spare + offset2;
5658 const int * spareIndexTemp = spareIndex + offset2;
5659 for (i = 0; i < numberLook; i++) {
5660 double value = spareTemp[i];
5661 spareTemp[i] = 0.0;
5662 spare[numberRemaining] = value;
5663 spareIndex[numberRemaining++] = spareIndexTemp[i];
5664 }
5665 #endif
5666 if (dwork[2] > freePivot) {
5667 freePivot = dwork[2];
5668 posFree = iwork[2] + numberNonZero;
5669 }
5670 upperTheta = CoinMin(dwork[1], upperTheta);
5671 bestPossible = CoinMax(dwork[0], bestPossible);
5672 for (i = 0; i < number; i++) {
5673 // double value = arrayTemp[i];
5674 //arrayTemp[i]=0.0;
5675 //array[numberNonZero]=value;
5676 index[numberNonZero++] = indexTemp[i] + offset3;
5677 }
5678 #else
5679 int offset = offset_[iBlock];
5680 // allow for already saved
5681 int offset2 = offset + saveNumberRemaining;
5682 double * arrayTemp = array + offset;
5683 int * indexTemp = index + offset;
5684 dualColumn0Struct * infoPtr = info_ + iBlock;
5685 infoPtr->model = model;
5686 infoPtr->spare = spare + offset2;
5687 infoPtr->spareIndex = spareIndex + offset2;
5688 infoPtr->arrayTemp = arrayTemp;
5689 infoPtr->indexTemp = indexTemp;
5690 infoPtr->numberInPtr = &iwork[0];
5691 infoPtr->offset = offset;
5692 infoPtr->acceptablePivot = acceptablePivot;
5693 infoPtr->bestPossiblePtr = &dwork[0];
5694 infoPtr->upperThetaPtr = &dwork[1];
5695 infoPtr->posFreePtr = &iwork[2];
5696 infoPtr->freePivotPtr = &dwork[2];
5697 infoPtr->numberOutPtr = &iwork[1];
5698 infoPtr->pi = pi;
5699 infoPtr->rowStart = rowStart_ + numberInRowArray * iBlock;
5700 infoPtr->element = element;
5701 infoPtr->column = column_;
5702 infoPtr->numberInRowArray = numberInRowArray;
5703 infoPtr->numberLook = offset_[iBlock+1] - offset;
5704 if (iBlock >= 2)
5705 pthread_join(threadId_[iBlock-2], NULL);
5706 pthread_create(threadId_ + iBlock, NULL, doOneBlockAnd0Thread, infoPtr);
5707 //pthread_join(threadId_[iBlock],NULL);
5708 #endif
5709 }
5710 }
5711 for ( iBlock = CoinMax(0, numberBlocks_ - 2); iBlock < numberBlocks_; iBlock++) {
5712 #ifdef THREAD
5713 pthread_join(threadId_[iBlock], NULL);
5714 #endif
5715 }
5716 #ifdef THREAD
5717 for ( iBlock = 0; iBlock < numberBlocks_; iBlock++) {
5718 //pthread_join(threadId_[iBlock],NULL);
5719 int offset = offset_[iBlock];
5720 double * dwork = work_ + 6 * iBlock;
5721 int * iwork = (int *) (dwork + 3);
5722 int number = iwork[0];
5723 if (dualColumn) {
5724 // allow for already saved
5725 int offset2 = offset + saveNumberRemaining;
5726 int numberLook = iwork[1];
5727 double * spareTemp = spare + offset2;
5728 const int * spareIndexTemp = spareIndex + offset2;
5729 for (i = 0; i < numberLook; i++) {
5730 double value = spareTemp[i];
5731 spareTemp[i] = 0.0;
5732 spare[numberRemaining] = value;
5733 spareIndex[numberRemaining++] = spareIndexTemp[i];
5734 }
5735 if (dwork[2] > freePivot) {
5736 freePivot = dwork[2];
5737 posFree = iwork[2] + numberNonZero;
5738 }
5739 upperTheta = CoinMin(dwork[1], upperTheta);
5740 bestPossible = CoinMax(dwork[0], bestPossible);
5741 }
5742 double * arrayTemp = array + offset;
5743 const int * indexTemp = index + offset;
5744 for (i = 0; i < number; i++) {
5745 double value = arrayTemp[i];
5746 arrayTemp[i] = 0.0;
5747 array[numberNonZero] = value;
5748 index[numberNonZero++] = indexTemp[i] + offset;
5749 }
5750 }
5751 #endif
5752 columnArray->setNumElements(numberNonZero);
5753 columnArray->setPackedMode(true);
5754 if (dualColumn) {
5755 model->spareDoubleArray_[0] = upperTheta;
5756 model->spareDoubleArray_[1] = bestPossible;
5757 // and theta and alpha and sequence
5758 if (posFree < 0) {
5759 model->spareIntArray_[1] = -1;
5760 } else {
5761 const double * reducedCost = model->djRegion(0);
5762 double alpha;
5763 int numberColumns = model->numberColumns();
5764 if (posFree < numberColumns) {
5765 alpha = columnArray->denseVector()[posFree];
5766 posFree = columnArray->getIndices()[posFree];
5767 } else {
5768 alpha = rowArray->denseVector()[posFree-numberColumns];
5769 posFree = rowArray->getIndices()[posFree-numberColumns] + numberColumns;
5770 }
5771 model->spareDoubleArray_[2] = fabs(reducedCost[posFree] / alpha);
5772 model->spareDoubleArray_[3] = alpha;
5773 model->spareIntArray_[1] = posFree;
5774 }
5775 spareArray->setNumElements(numberRemaining);
5776 // signal done
5777 model->spareIntArray_[0] = -1;
5778 }
5779 }
5780 /* Default constructor. */
5781 ClpPackedMatrix3::ClpPackedMatrix3()
5782 : numberBlocks_(0),
5783 numberColumns_(0),
5784 column_(NULL),
5785 start_(NULL),
5786 row_(NULL),
5787 element_(NULL),
5788 block_(NULL)
5789 {
5790 }
5791 /* Constructor from copy. */
5792 ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model, const CoinPackedMatrix * columnCopy)
5793 : numberBlocks_(0),
5794 numberColumns_(0),
5795 column_(NULL),
5796 start_(NULL),
5797 row_(NULL),
5798 element_(NULL),
5799 block_(NULL)
5800 {
5801 #define MINBLOCK 6
5802 #define MAXBLOCK 100
5803 #define MAXUNROLL 10
5804 numberColumns_ = model->getNumCols();
5805 int numberColumns = columnCopy->getNumCols();
5806 assert (numberColumns_ >= numberColumns);
5807 int numberRows = columnCopy->getNumRows();
5808 int * counts = new int[numberRows+1];
5809 CoinZeroN(counts, numberRows + 1);
5810 CoinBigIndex nels = 0;
5811 CoinBigIndex nZeroEl = 0;
5812 int iColumn;
5813 // get matrix data pointers
5814 const int * row = columnCopy->getIndices();
5815 const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
5816 const int * columnLength = columnCopy->getVectorLengths();
5817 const double * elementByColumn = columnCopy->getElements();
5818 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5819 CoinBigIndex start = columnStart[iColumn];
5820 int n = columnLength[iColumn];
5821 CoinBigIndex end = start + n;
5822 int kZero = 0;
5823 for (CoinBigIndex j = start; j < end; j++) {
5824 if(!elementByColumn[j])
5825 kZero++;
5826 }
5827 n -= kZero;
5828 nZeroEl += kZero;
5829 nels += n;
5830 counts[n]++;
5831 }
5832 counts[0] += numberColumns_ - numberColumns;
5833 int nZeroColumns = counts[0];
5834 counts[0] = -1;
5835 numberColumns_ -= nZeroColumns;
5836 column_ = new int [2*numberColumns_+nZeroColumns];
5837 int * lookup = column_ + numberColumns_;
5838 row_ = new int[nels];
5839 element_ = new double[nels];
5840 int nOdd = 0;
5841 CoinBigIndex nInOdd = 0;
5842 int i;
5843 for (i = 1; i <= numberRows; i++) {
5844 int n = counts[i];
5845 if (n) {
5846 if (n < MINBLOCK || i > MAXBLOCK) {
5847 nOdd += n;
5848 counts[i] = -1;
5849 nInOdd += n * i;
5850 } else {
5851 numberBlocks_++;
5852 }
5853 } else {
5854 counts[i] = -1;
5855 }
5856 }
5857 start_ = new CoinBigIndex [nOdd+1];
5858 // even if no blocks do a dummy one
5859 numberBlocks_ = CoinMax(numberBlocks_, 1);
5860 block_ = new blockStruct [numberBlocks_];
5861 memset(block_, 0, numberBlocks_ * sizeof(blockStruct));
5862 // Fill in what we can
5863 int nTotal = nOdd;
5864 // in case no blocks
5865 block_->startIndices_ = nTotal;
5866 nels = nInOdd;
5867 int nBlock = 0;
5868 for (i = 0; i <= CoinMin(MAXBLOCK, numberRows); i++) {
5869 if (counts[i] > 0) {
5870 blockStruct * block = block_ + nBlock;
5871 int n = counts[i];
5872 counts[i] = nBlock; // backward pointer
5873 nBlock++;
5874 block->startIndices_ = nTotal;
5875 block->startElements_ = nels;
5876 block->numberElements_ = i;
5877 // up counts
5878 nTotal += n;
5879 nels += n * i;
5880 }
5881 }
5882 for (iColumn = numberColumns; iColumn < numberColumns_; iColumn++)
5883 lookup[iColumn] = -1;
5884 // fill
5885 start_[0] = 0;
5886 nOdd = 0;
5887 nInOdd = 0;
5888 const double * columnScale = model->columnScale();
5889 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5890 CoinBigIndex start = columnStart[iColumn];
5891 int n = columnLength[iColumn];
5892 CoinBigIndex end = start + n;
5893 int kZero = 0;
5894 for (CoinBigIndex j = start; j < end; j++) {
5895 if(!elementByColumn[j])
5896 kZero++;
5897 }
5898 n -= kZero;
5899 if (n) {
5900 int iBlock = counts[n];
5901 if (iBlock >= 0) {
5902 blockStruct * block = block_ + iBlock;
5903 int k = block->numberInBlock_;
5904 block->numberInBlock_ ++;
5905 column_[block->startIndices_+k] = iColumn;
5906 lookup[iColumn] = k;
5907 CoinBigIndex put = block->startElements_ + k * n;
5908 for (CoinBigIndex j = start; j < end; j++) {
5909 double value = elementByColumn[j];
5910 if(value) {
5911 if (columnScale)
5912 value *= columnScale[iColumn];
5913 element_[put] = value;
5914 row_[put++] = row[j];
5915 }
5916 }
5917 } else {
5918 // odd ones
5919 for (CoinBigIndex j = start; j < end; j++) {
5920 double value = elementByColumn[j];
5921 if(value) {
5922 if (columnScale)
5923 value *= columnScale[iColumn];
5924 element_[nInOdd] = value;
5925 row_[nInOdd++] = row[j];
5926 }
5927 }
5928 column_[nOdd] = iColumn;
5929 lookup[iColumn] = -1;
5930 nOdd++;
5931 start_[nOdd] = nInOdd;
5932 }
5933 } else {
5934 // zero column
5935 lookup[iColumn] = -1;
5936 }
5937 }
5938 delete [] counts;
5939 }
5940 /* Destructor */
5941 ClpPackedMatrix3::~ClpPackedMatrix3()
5942 {
5943 delete [] column_;
5944 delete [] start_;
5945 delete [] row_;
5946 delete [] element_;
5947 delete [] block_;
5948 }
5949 /* The copy constructor. */
5950 ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
5951 : numberBlocks_(rhs.numberBlocks_),
5952 numberColumns_(rhs.numberColumns_),
5953 column_(NULL),
5954 start_(NULL),
5955 row_(NULL),
5956 element_(NULL),
5957 block_(NULL)
5958 {
5959 if (rhs.numberBlocks_) {
5960 block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
5961 column_ = CoinCopyOfArray(rhs.column_, 2 * numberColumns_);
5962 int numberOdd = block_->startIndices_;
5963 start_ = CoinCopyOfArray(rhs.start_, numberOdd + 1);
5964 blockStruct * lastBlock = block_ + (numberBlocks_ - 1);
5965 CoinBigIndex numberElements = lastBlock->startElements_ +
5966 lastBlock->numberInBlock_ * lastBlock->numberElements_;
5967 row_ = CoinCopyOfArray(rhs.row_, numberElements);
5968 element_ = CoinCopyOfArray(rhs.element_, numberElements);
5969 }
5970 }
5971 ClpPackedMatrix3&
5972 ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
5973 {
5974 if (this != &rhs) {
5975 delete [] column_;
5976 delete [] start_;
5977 delete [] row_;
5978 delete [] element_;
5979 delete [] block_;
5980 numberBlocks_ = rhs.numberBlocks_;
5981 numberColumns_ = rhs.numberColumns_;
5982 if (rhs.numberBlocks_) {
5983 block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
5984 column_ = CoinCopyOfArray(rhs.column_, 2 * numberColumns_);
5985 int numberOdd = block_->startIndices_;
5986 start_ = CoinCopyOfArray(rhs.start_, numberOdd + 1);
5987 blockStruct * lastBlock = block_ + (numberBlocks_ - 1);
5988 CoinBigIndex numberElements = lastBlock->startElements_ +
5989 lastBlock->numberInBlock_ * lastBlock->numberElements_;
5990 row_ = CoinCopyOfArray(rhs.row_, numberElements);
5991 element_ = CoinCopyOfArray(rhs.element_, numberElements);
5992 } else {
5993 column_ = NULL;
5994 start_ = NULL;
5995 row_ = NULL;
5996 element_ = NULL;
5997 block_ = NULL;
5998 }
5999 }
6000 return *this;
6001 }
6002 /* Sort blocks */
6003 void
6004 ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
6005 {
6006 int * lookup = column_ + numberColumns_;
6007 for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
6008 blockStruct * block = block_ + iBlock;
6009 int numberInBlock = block->numberInBlock_;
6010 int nel = block->numberElements_;
6011 int * row = row_ + block->startElements_;
6012 double * element = element_ + block->startElements_;
6013 int * column = column_ + block->startIndices_;
6014 int lastPrice = 0;
6015 int firstNotPrice = numberInBlock - 1;
6016 while (lastPrice <= firstNotPrice) {
6017 // find first basic or fixed
6018 int iColumn = numberInBlock;
6019 for (; lastPrice <= firstNotPrice; lastPrice++) {
6020 iColumn = column[lastPrice];
6021 if (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
6022 model->getColumnStatus(iColumn) == ClpSimplex::isFixed)
6023 break;
6024 }
6025 // find last non basic or fixed
6026 int jColumn = -1;
6027 for (; firstNotPrice > lastPrice; firstNotPrice--) {
6028 jColumn = column[firstNotPrice];
6029 if (model->getColumnStatus(jColumn) != ClpSimplex::basic &&
6030 model->getColumnStatus(jColumn) != ClpSimplex::isFixed)
6031 break;
6032 }
6033 if (firstNotPrice > lastPrice) {
6034 assert (column[lastPrice] == iColumn);
6035 assert (column[firstNotPrice] == jColumn);
6036 // need to swap
6037 column[firstNotPrice] = iColumn;
6038 lookup[iColumn] = firstNotPrice;
6039 column[lastPrice] = jColumn;
6040 lookup[jColumn] = lastPrice;
6041 double * elementA = element + lastPrice * nel;
6042 int * rowA = row + lastPrice * nel;
6043 double * elementB = element + firstNotPrice * nel;
6044 int * rowB = row + firstNotPrice * nel;
6045 for (int i = 0; i < nel; i++) {
6046 int temp = rowA[i];
6047 double tempE = elementA[i];
6048 rowA[i] = rowB[i];
6049 elementA[i] = elementB[i];
6050 rowB[i] = temp;
6051 elementB[i] = tempE;
6052 }
6053 firstNotPrice--;
6054 lastPrice++;
6055 } else if (lastPrice == firstNotPrice) {
6056 // make sure correct side
6057 iColumn = column[lastPrice];
6058 if (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
6059 model->getColumnStatus(iColumn) != ClpSimplex::isFixed)
6060 lastPrice++;
6061 break;
6062 }
6063 }
6064 block->numberPrice_ = lastPrice;
6065 #ifndef NDEBUG
6066 // check
6067 int i;
6068 for (i = 0; i < lastPrice; i++) {
6069 int iColumn = column[i];
6070 assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
6071 model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
6072 assert (lookup[iColumn] == i);
6073 }
6074 for (; i < numberInBlock; i++) {
6075 int iColumn = column[i];
6076 assert (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
6077 model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
6078 assert (lookup[iColumn] == i);
6079 }
6080 #endif
6081 }
6082 }
6083 // Swap one variable
6084 void
6085 ClpPackedMatrix3::swapOne(const ClpSimplex * model, const ClpPackedMatrix * matrix,
6086 int iColumn)
6087 {
6088 int * lookup = column_ + numberColumns_;
6089 // position in block
6090 int kA = lookup[iColumn];
6091 if (kA < 0)
6092 return; // odd one
6093 // get matrix data pointers
6094 const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
6095 //const int * row = columnCopy->getIndices();
6096 const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
6097 const int * columnLength = columnCopy->getVectorLengths();
6098 const double * elementByColumn = columnCopy->getElements();
6099 CoinBigIndex start = columnStart[iColumn];
6100 int n = columnLength[iColumn];
6101 if (matrix->zeros()) {
6102 CoinBigIndex end = start + n;
6103 for (CoinBigIndex j = start; j < end; j++) {
6104 if(!elementByColumn[j])
6105 n--;
6106 }
6107 }
6108 // find block - could do binary search
6109 int iBlock = CoinMin(n, numberBlocks_) - 1;
6110 while (block_[iBlock].numberElements_ != n)
6111 iBlock--;
6112 blockStruct * block = block_ + iBlock;
6113 int nel = block->numberElements_;
6114 int * row = row_ + block->startElements_;
6115 double * element = element_ + block->startElements_;
6116 int * column = column_ + block->startIndices_;
6117 assert (column[kA] == iColumn);
6118 bool moveUp = (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
6119 model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
6120 int lastPrice = block->numberPrice_;
6121 int kB;
6122 if (moveUp) {
6123 // May already be in correct place (e.g. fixed basic leaving basis)
6124 if (kA >= lastPrice)
6125 return;
6126 kB = lastPrice - 1;
6127 block->numberPrice_--;
6128 } else {
6129 assert (kA >= lastPrice);
6130 kB = lastPrice;
6131 block->numberPrice_++;
6132 }
6133 int jColumn = column[kB];
6134 column[kA] = jColumn;
6135 lookup[jColumn] = kA;
6136 column[kB] = iColumn;
6137 lookup[iColumn] = kB;
6138 double * elementA = element + kB * nel;
6139 int * rowA = row + kB * nel;
6140 double * elementB = element + kA * nel;
6141 int * rowB = row + kA * nel;
6142 int i;
6143 for (i = 0; i < nel; i++) {
6144 int temp = rowA[i];
6145 double tempE = elementA[i];
6146 rowA[i] = rowB[i];
6147 elementA[i] = elementB[i];
6148 rowB[i] = temp;
6149 elementB[i] = tempE;
6150 }
6151 #ifndef NDEBUG
6152 // check
6153 for (i = 0; i < block->numberPrice_; i++) {
6154 int iColumn = column[i];
6155 if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
6156 assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
6157 model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
6158 assert (lookup[iColumn] == i);
6159 }
6160 int numberInBlock = block->numberInBlock_;
6161 for (; i < numberInBlock; i++) {
6162 int iColumn = column[i];
6163 if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
6164 assert (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
6165 model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
6166 assert (lookup[iColumn] == i);
6167 }
6168 #endif
6169 }
6170 /* Return <code>x * -1 * A in <code>z</code>.
6171 Note - x packed and z will be packed mode
6172 Squashes small elements and knows about ClpSimplex */
6173 void
6174 ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
6175 const double * pi,
6176 CoinIndexedVector * output) const
6177 {
6178 int numberNonZero = 0;
6179 int * index = output->getIndices();
6180 double * array = output->denseVector();
6181 double zeroTolerance = model->zeroTolerance();
6182 double value = 0.0;
6183 CoinBigIndex j;
6184 int numberOdd = block_->startIndices_;
6185 if (numberOdd) {
6186 // A) as probably long may be worth unrolling
6187 CoinBigIndex end = start_[1];
6188 for (j = start_[0]; j < end; j++) {
6189 int iRow = row_[j];
6190 value += pi[iRow] * element_[j];
6191 }
6192 int iColumn;
6193 // int jColumn=column_[0];
6194
6195 for (iColumn = 0; iColumn < numberOdd - 1; iColumn++) {
6196 CoinBigIndex start = end;
6197 end = start_[iColumn+2];
6198 if (fabs(value) > zeroTolerance) {
6199 array[numberNonZero] = value;
6200 index[numberNonZero++] = column_[iColumn];
6201 //index[numberNonZero++]=jColumn;
6202 }
6203 // jColumn = column_[iColumn+1];
6204 value = 0.0;
6205 //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
6206 for (j = start; j < end; j++) {
6207 int iRow = row_[j];
6208 value += pi[iRow] * element_[j];
6209 }
6210 //}
6211 }
6212 if (fabs(value) > zeroTolerance) {
6213 array[numberNonZero] = value;
6214 index[numberNonZero++] = column_[iColumn];
6215 //index[numberNonZero++]=jColumn;
6216 }
6217 }
6218 for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
6219 // B) Can sort so just do nonbasic (and nonfixed)
6220 // C) Can do two at a time (if so put odd one into start_)
6221 // D) can use switch
6222 blockStruct * block = block_ + iBlock;
6223 //int numberPrice = block->numberInBlock_;
6224 int numberPrice = block->numberPrice_;
6225 int nel = block->numberElements_;
6226 int * row = row_ + block->startElements_;
6227 double * element = element_ + block->startElements_;
6228 int * column = column_ + block->startIndices_;
6229 #if 0
6230 // two at a time
6231 if ((numberPrice & 1) != 0) {
6232 double value = 0.0;
6233 int nel2 = nel;
6234 for (; nel2; nel2--) {
6235 int iRow = *row++;
6236 value += pi[iRow] * (*element++);
6237 }
6238 if (fabs(value) > zeroTolerance) {
6239 array[numberNonZero] = value;
6240 index[numberNonZero++] = *column;
6241 }
6242 column++;
6243 }
6244 numberPrice = numberPrice >> 1;
6245 switch ((nel % 2)) {
6246 // 2 k +0
6247 case 0:
6248 for (; numberPrice; numberPrice--) {
6249 double value0 = 0.0;
6250 double value1 = 0.0;
6251 int nel2 = nel;
6252 for (; nel2; nel2--) {
6253 int iRow0 = *row;
6254 int iRow1 = *(row + nel);
6255 row++;
6256 double element0 = *element;
6257 double element1 = *(element + nel);
6258 element++;
6259 value0 += pi[iRow0] * element0;
6260 value1 += pi[iRow1] * element1;
6261 }
6262 row += nel;
6263 element += nel;
6264 if (fabs(value0) > zeroTolerance) {
6265 array[numberNonZero] = value0;
6266 index[numberNonZero++] = *column;
6267 }
6268 column++;
6269 if (fabs(value1) > zeroTolerance) {
6270 array[numberNonZero] = value1;
6271 index[numberNonZero++] = *column;
6272 }
6273 column++;
6274 }
6275 break;
6276 // 2 k +1
6277 case 1:
6278 for (; numberPrice; numberPrice--) {
6279 double value0;
6280 double value1;
6281 int nel2 = nel - 1;
6282 {
6283 int iRow0 = row[0];
6284 int iRow1 = row[nel];
6285 double pi0 = pi[iRow0];
6286 double pi1 = pi[iRow1];
6287 value0 = pi0 * element[0];
6288 value1 = pi1 * element[nel];
6289 row++;
6290 element++;
6291 }
6292 for (; nel2; nel2--) {
6293 int iRow0 = *row;
6294 int iRow1 = *(row + nel);
6295 row++;
6296 double element0 = *element;
6297 double element1 = *(element + nel);
6298 element++;
6299 value0 += pi[iRow0] * element0;
6300 value1 += pi[iRow1] * element1;
6301 }
6302 row += nel;
6303 element += nel;
6304 if (fabs(value0) > zeroTolerance) {
6305 array[numberNonZero] = value0;
6306 index[numberNonZero++] = *column;
6307 }
6308 column++;
6309 if (fabs(value1) > zeroTolerance) {
6310 array[numberNonZero] = value1;
6311 index[numberNonZero++] = *column;
6312 }
6313 column++;
6314 }
6315 break;
6316 }
6317 #else
6318 for (; numberPrice; numberPrice--) {
6319 double value = 0.0;
6320 int nel2 = nel;
6321 for (; nel2; nel2--) {
6322 int iRow = *row++;
6323 value += pi[iRow] * (*element++);
6324 }
6325 if (fabs(value) > zeroTolerance) {
6326 array[numberNonZero] = value;
6327 index[numberNonZero++] = *column;
6328 }
6329 column++;
6330 }
6331 #endif
6332 }
6333 output->setNumElements(numberNonZero);
6334 }
6335 // Updates two arrays for steepest
6336 void
6337 ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
6338 const double * pi, CoinIndexedVector * output,
6339 const double * piWeight,
6340 double referenceIn, double devex,
6341 // Array for exact devex to say what is in reference framework
6342 unsigned int * reference,
6343 double * weights, double scaleFactor)
6344 {
6345 int numberNonZero = 0;
6346 int * index = output->getIndices();
6347 double * array = output->denseVector();
6348 double zeroTolerance = model->zeroTolerance();
6349 double value = 0.0;
6350 bool killDjs = (scaleFactor == 0.0);
6351 if (!scaleFactor)
6352 scaleFactor = 1.0;
6353 int numberOdd = block_->startIndices_;
6354 int iColumn;
6355 CoinBigIndex end = start_[0];
6356 for (iColumn = 0; iColumn < numberOdd; iColumn++) {
6357 CoinBigIndex start = end;
6358 CoinBigIndex j;
6359 int jColumn = column_[iColumn];
6360 end = start_[iColumn+1];
6361 value = 0.0;
6362 if (model->getColumnStatus(jColumn) != ClpSimplex::basic) {
6363 for (j = start; j < end; j++) {
6364 int iRow = row_[j];
6365 value -= pi[iRow] * element_[j];
6366 }
6367 if (fabs(value) > zeroTolerance) {
6368 // and do other array
6369 double modification = 0.0;
6370 for (j = start; j < end; j++) {
6371 int iRow = row_[j];
6372 modification += piWeight[iRow] * element_[j];
6373 }
6374 double thisWeight = weights[jColumn];
6375 double pivot = value * scaleFactor;
6376 double pivotSquared = pivot * pivot;
6377 thisWeight += pivotSquared * devex + pivot * modification;
6378 if (thisWeight < DEVEX_TRY_NORM) {
6379 if (referenceIn < 0.0) {
6380 // steepest
6381 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
6382 } else {
6383 // exact
6384 thisWeight = referenceIn * pivotSquared;
6385 if (reference(jColumn))
6386 thisWeight += 1.0;
6387 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
6388 }
6389 }
6390 weights[jColumn] = thisWeight;
6391 if (!killDjs) {
6392 array[numberNonZero] = value;
6393 index[numberNonZero++] = jColumn;
6394 }
6395 }
6396 }
6397 }
6398 for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
6399 // B) Can sort so just do nonbasic (and nonfixed)
6400 // C) Can do two at a time (if so put odd one into start_)
6401 // D) can use switch
6402 blockStruct * block = block_ + iBlock;
6403 //int numberPrice = block->numberInBlock_;
6404 int numberPrice = block->numberPrice_;
6405 int nel = block->numberElements_;
6406 int * row = row_ + block->startElements_;
6407 double * element = element_ + block->startElements_;
6408 int * column = column_ + block->startIndices_;
6409 for (; numberPrice; numberPrice--) {
6410 double value = 0.0;
6411 int nel2 = nel;
6412 for (; nel2; nel2--) {
6413 int iRow = *row++;
6414 value -= pi[iRow] * (*element++);
6415 }
6416 if (fabs(value) > zeroTolerance) {
6417 int jColumn = *column;
6418 // back to beginning
6419 row -= nel;
6420 element -= nel;
6421 // and do other array
6422 double modification = 0.0;
6423 nel2 = nel;
6424 for (; nel2; nel2--) {
6425 int iRow = *row++;
6426 modification += piWeight[iRow] * (*element++);
6427 }
6428 double thisWeight = weights[jColumn];
6429 double pivot = value * scaleFactor;
6430 double pivotSquared = pivot * pivot;
6431 thisWeight += pivotSquared * devex + pivot * modification;
6432 if (thisWeight < DEVEX_TRY_NORM) {
6433 if (referenceIn < 0.0) {
6434 // steepest
6435 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
6436 } else {
6437 // exact
6438 thisWeight = referenceIn * pivotSquared;
6439 if (reference(jColumn))
6440 thisWeight += 1.0;
6441 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
6442 }
6443 }
6444 weights[jColumn] = thisWeight;
6445 if (!killDjs) {
6446 array[numberNonZero] = value;
6447 index[numberNonZero++] = jColumn;
6448 }
6449 }
6450 column++;
6451 }
6452 }
6453 output->setNumElements(numberNonZero);
6454 output->setPackedMode(true);
6455 }
6456 #if COIN_LONG_WORK
6457 // For long double versions
6458 void
6459 ClpPackedMatrix::times(CoinWorkDouble scalar,
6460 const CoinWorkDouble * x, CoinWorkDouble * y) const
6461 {
6462 int iRow, iColumn;
6463 // get matrix data pointers
6464 const int * row = matrix_->getIndices();
6465 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
6466 const double * elementByColumn = matrix_->getElements();
6467 //memset(y,0,matrix_->getNumRows()*sizeof(double));
6468 assert (((flags_ & 2) != 0) == matrix_->hasGaps());
6469 if (!(flags_ & 2)) {
6470 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
6471 CoinBigIndex j;
6472 CoinWorkDouble value = x[iColumn];
6473 if (value) {
6474 CoinBigIndex start = columnStart[iColumn];
6475 CoinBigIndex end = columnStart[iColumn+1];
6476 value *= scalar;
6477 for (j = start; j < end; j++) {
6478 iRow = row[j];
6479 y[iRow] += value * elementByColumn[j];
6480 }
6481 }
6482 }
6483 } else {
6484 const int * columnLength = matrix_->getVectorLengths();
6485 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
6486 CoinBigIndex j;
6487 CoinWorkDouble value = x[iColumn];
6488 if (value) {
6489 CoinBigIndex start = columnStart[iColumn];
6490 CoinBigIndex end = start + columnLength[iColumn];
6491 value *= scalar;
6492 for (j = start; j < end; j++) {
6493 iRow = row[j];
6494 y[iRow] += value * elementByColumn[j];
6495 }
6496 }
6497 }
6498 }
6499 }
6500 void
6501 ClpPackedMatrix::transposeTimes(CoinWorkDouble scalar,
6502 const CoinWorkDouble * x, CoinWorkDouble * y) const
6503 {
6504 int iColumn;
6505 // get matrix data pointers
6506 const int * row = matrix_->getIndices();
6507 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
6508 const double * elementByColumn = matrix_->getElements();
6509 if (!(flags_ & 2)) {
6510 if (scalar == -1.0) {
6511 CoinBigIndex start = columnStart[0];
6512 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
6513 CoinBigIndex j;
6514 CoinBigIndex next = columnStart[iColumn+1];
6515 CoinWorkDouble value = y[iColumn];
6516 for (j = start; j < next; j++) {
6517 int jRow = row[j];
6518 value -= x[jRow] * elementByColumn[j];
6519 }
6520 start = next;
6521 y[iColumn] = value;
6522 }
6523 } else {
6524 CoinBigIndex start = columnStart[0];
6525 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
6526 CoinBigIndex j;
6527 CoinBigIndex next = columnStart[iColumn+1];
6528 CoinWorkDouble value = 0.0;
6529 for (j = start; j < next; j++) {
6530 int jRow = row[j];
6531 value += x[jRow] * elementByColumn[j];
6532 }
6533 start = next;
6534 y[iColumn] += value * scalar;
6535 }
6536 }
6537 } else {
6538 const int * columnLength = matrix_->getVectorLengths();
6539 for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
6540 CoinBigIndex j;
6541 CoinWorkDouble value = 0.0;
6542 CoinBigIndex start = columnStart[iColumn];
6543 CoinBigIndex end = start + columnLength[iColumn];
6544 for (j = start; j < end; j++) {
6545 int jRow = row[j];
6546 value += x[jRow] * elementByColumn[j];
6547 }
6548 y[iColumn] += value * scalar;
6549 }
6550 }
6551 }
6552 #endif
6553 #ifdef CLP_ALL_ONE_FILE
6554 #undef reference
6555 #endif
6556