source: trunk/Clp/src/ClpPackedMatrix.cpp @ 1726

Last change on this file since 1726 was 1726, checked in by stefan, 8 years ago

fix compiler warnings, including one that pointed to a bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 274.1 KB
Line 
1/* $Id: ClpPackedMatrix.cpp 1726 2011-05-02 08:58:39Z stefan $ */
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//-------------------------------------------------------------------
54ClpPackedMatrix::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//-------------------------------------------------------------------
68ClpPackedMatrix::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//-------------------------------------------------------------------
100ClpPackedMatrix::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
112ClpPackedMatrix::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//-------------------------------------------------------------------
131ClpPackedMatrix::~ClpPackedMatrix ()
132{
133     delete matrix_;
134     delete rowCopy_;
135     delete columnCopy_;
136}
137
138//----------------------------------------------------------------
139// Assignment operator
140//-------------------------------------------------------------------
141ClpPackedMatrix &
142ClpPackedMatrix::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//-------------------------------------------------------------------
174ClpMatrixBase * ClpPackedMatrix::clone() const
175{
176     return new ClpPackedMatrix(*this);
177}
178// Copy contents - resizing if necessary - otherwise re-use memory
179void
180ClpPackedMatrix::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 */
189ClpMatrixBase *
190ClpPackedMatrix::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 */
199ClpPackedMatrix::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}
212ClpPackedMatrix::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 */
228ClpMatrixBase *
229ClpPackedMatrix::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
242void
243ClpPackedMatrix::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}
284void
285ClpPackedMatrix::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}
336void
337ClpPackedMatrix::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}
384void
385ClpPackedMatrix::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}
487void
488ClpPackedMatrix::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 */
554void
555ClpPackedMatrix::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 */
809void
810ClpPackedMatrix::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 */
1143void
1144ClpPackedMatrix::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
1310int
1311ClpPackedMatrix::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
1370int
1371ClpPackedMatrix::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
1414int
1415ClpPackedMatrix::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 */
1466int
1467ClpPackedMatrix::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
1622int
1623ClpPackedMatrix::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
1662int
1663ClpPackedMatrix::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
1715void
1716ClpPackedMatrix::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
1813void
1814ClpPackedMatrix::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 */
1842void
1843ClpPackedMatrix::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 */
1981bool
1982ClpPackedMatrix::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
2013void
2014ClpPackedMatrix::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
2308void
2309ClpPackedMatrix::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,
2403CoinBigIndex
2404ClpPackedMatrix::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}
2417void
2418ClpPackedMatrix::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
2526int
2527ClpPackedMatrix::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
3201static 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)
3209int
3210ClpPackedMatrix::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
3833void
3834ClpPackedMatrix::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 */
3873void
3874ClpPackedMatrix::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
3900Note that model is NOT const.  Bounds and objective could
3901be modified if doing column generation (just for this variable) */
3902void
3903ClpPackedMatrix::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 */
3946void
3947ClpPackedMatrix::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 */
3973void
3974ClpPackedMatrix::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*/
4009bool
4010ClpPackedMatrix::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                         return false;
4081                    }
4082                    if (mark[iRow] == -1) {
4083                         mark[iRow] = j;
4084                    } else {
4085                         // duplicate
4086                         numberDuplicate++;
4087                    }
4088                    //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
4089                    if (!value)
4090                         flags_ |= 1; // there are zero elements
4091                    if (value < smallest) {
4092                         numberSmall++;
4093                    } else if (!(value <= largest)) {
4094                         numberLarge++;
4095                         if (firstBadColumn < 0) {
4096                              firstBadColumn = iColumn;
4097                              firstBadRow = row[j];
4098                              firstBadElement = elementByColumn[j];
4099                         }
4100                    }
4101               }
4102               //clear mark
4103               for (j = columnStart[iColumn];
4104                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4105                    int iRow = row[j];
4106                    mark[iRow] = -1;
4107               }
4108          }
4109          delete [] mark;
4110     } else {
4111          // just check for out of range - not for duplicates
4112          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4113               CoinBigIndex j;
4114               CoinBigIndex start = columnStart[iColumn];
4115               CoinBigIndex end = start + columnLength[iColumn];
4116               if (end != columnStart[iColumn+1])
4117                    flags_ |= 2;
4118               for (j = start; j < end; j++) {
4119                    double value = fabs(elementByColumn[j]);
4120                    int iRow = row[j];
4121                    if (iRow < 0 || iRow >= numberRows) {
4122#ifndef COIN_BIG_INDEX
4123                         printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4124#elif COIN_BIG_INDEX==0
4125                         printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4126#elif COIN_BIG_INDEX==1
4127                         printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4128#else
4129                         printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4130#endif
4131                         return false;
4132                    }
4133                    if (!value)
4134                         flags_ |= 1; // there are zero elements
4135                    if (value < smallest) {
4136                         numberSmall++;
4137                    } else if (!(value <= largest)) {
4138                         numberLarge++;
4139                         if (firstBadColumn < 0) {
4140                              firstBadColumn = iColumn;
4141                              firstBadRow = iRow;
4142                              firstBadElement = value;
4143                         }
4144                    }
4145               }
4146          }
4147     }
4148     if (numberLarge) {
4149          model->messageHandler()->message(CLP_BAD_MATRIX, model->messages())
4150                    << numberLarge
4151                    << firstBadColumn << firstBadRow << firstBadElement
4152                    << CoinMessageEol;
4153          return false;
4154     }
4155     if (numberSmall)
4156          model->messageHandler()->message(CLP_SMALLELEMENTS, model->messages())
4157                    << numberSmall
4158                    << CoinMessageEol;
4159     if (numberDuplicate)
4160          model->messageHandler()->message(CLP_DUPLICATEELEMENTS, model->messages())
4161                    << numberDuplicate
4162                    << CoinMessageEol;
4163     if (numberDuplicate)
4164          matrix_->eliminateDuplicates(smallest);
4165     else if (numberSmall)
4166          matrix_->compress(smallest);
4167     // If smallest >0.0 then there can't be zero elements
4168     if (smallest > 0.0)
4169          flags_ &= ~1;;
4170     if (numberSmall || numberDuplicate)
4171          flags_ |= 2; // will have gaps
4172     return true;
4173}
4174int
4175ClpPackedMatrix::gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector * COIN_RESTRICT piVector,
4176          int * COIN_RESTRICT index,
4177          double * COIN_RESTRICT output,
4178          int * COIN_RESTRICT lookup,
4179          char * COIN_RESTRICT marked,
4180          const double tolerance,
4181          const double scalar) const
4182{
4183     const double * COIN_RESTRICT pi = piVector->denseVector();
4184     int numberNonZero = 0;
4185     int numberInRowArray = piVector->getNumElements();
4186     const int * COIN_RESTRICT column = matrix_->getIndices();
4187     const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4188     const double * COIN_RESTRICT element = matrix_->getElements();
4189     const int * COIN_RESTRICT whichRow = piVector->getIndices();
4190     int * fakeRow = const_cast<int *> (whichRow);
4191     fakeRow[numberInRowArray]=0; // so can touch
4192#ifndef NDEBUG
4193     int maxColumn = 0;
4194#endif
4195     // ** Row copy is already scaled
4196     int nextRow=whichRow[0];
4197     CoinBigIndex nextStart = rowStart[nextRow]; 
4198     CoinBigIndex nextEnd = rowStart[nextRow+1]; 
4199     for (int i = 0; i < numberInRowArray; i++) {
4200          double value = pi[i] * scalar;
4201          CoinBigIndex start=nextStart;
4202          CoinBigIndex end=nextEnd;
4203          nextRow=whichRow[i+1];
4204          nextStart = rowStart[nextRow]; 
4205          //coin_prefetch_const(column + nextStart);
4206          //coin_prefetch_const(element + nextStart);
4207          nextEnd = rowStart[nextRow+1]; 
4208          CoinBigIndex j;
4209          for (j = start; j < end; j++) {
4210               int iColumn = column[j];
4211#ifndef NDEBUG
4212               maxColumn = CoinMax(maxColumn, iColumn);
4213#endif
4214               double elValue = element[j];
4215               elValue *= value;
4216               if (marked[iColumn]) {
4217                    int k = lookup[iColumn];
4218                    output[k] += elValue;
4219               } else {
4220                    output[numberNonZero] = elValue;
4221                    marked[iColumn] = 1;
4222                    lookup[iColumn] = numberNonZero;
4223                    index[numberNonZero++] = iColumn;
4224               }
4225          }
4226     }
4227#ifndef NDEBUG
4228     int saveN = numberNonZero;
4229#endif
4230     // get rid of tiny values and zero out lookup
4231     for (int i = 0; i < numberNonZero; i++) {
4232          int iColumn = index[i];
4233          marked[iColumn] = 0;
4234          double value = output[i];
4235          if (fabs(value) <= tolerance) {
4236               while (fabs(value) <= tolerance) {
4237                    numberNonZero--;
4238                    value = output[numberNonZero];
4239                    iColumn = index[numberNonZero];
4240                    marked[iColumn] = 0;
4241                    if (i < numberNonZero) {
4242                         output[numberNonZero] = 0.0;
4243                         output[i] = value;
4244                         index[i] = iColumn;
4245                    } else {
4246                         output[i] = 0.0;
4247                         value = 1.0; // to force end of while
4248                    }
4249               }
4250          }
4251     }
4252#ifndef NDEBUG
4253     for (int i = numberNonZero; i < saveN; i++)
4254          assert(!output[i]);
4255     for (int i = 0; i <= maxColumn; i++)
4256          assert (!marked[i]);
4257#endif
4258     return numberNonZero;
4259}
4260int
4261ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
4262          int * COIN_RESTRICT index,
4263          double * COIN_RESTRICT output,
4264          double * COIN_RESTRICT array,
4265          const double tolerance,
4266          const double scalar) const
4267{
4268     const double * COIN_RESTRICT pi = piVector->denseVector();
4269     int numberNonZero = 0;
4270     int numberInRowArray = piVector->getNumElements();
4271     const int * COIN_RESTRICT column = matrix_->getIndices();
4272     const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4273     const double * COIN_RESTRICT element = matrix_->getElements();
4274     const int * COIN_RESTRICT whichRow = piVector->getIndices();
4275     // ** Row copy is already scaled
4276     for (int i = 0; i < numberInRowArray; i++) {
4277          int iRow = whichRow[i];
4278          double value = pi[i] * scalar;
4279          CoinBigIndex j;
4280          for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
4281               int iColumn = column[j];
4282               double inValue = array[iColumn];
4283               double elValue = element[j];
4284               elValue *= value;
4285               if (inValue) {
4286                 double outValue = inValue + elValue;
4287                 if (!outValue)
4288                   outValue = COIN_INDEXED_REALLY_TINY_ELEMENT;
4289                 array[iColumn] = outValue;
4290               } else {
4291                 array[iColumn] = elValue;
4292                 assert (elValue);
4293                 index[numberNonZero++] = iColumn;
4294               }
4295          }
4296     }
4297     int saveN = numberNonZero;
4298     // get rid of tiny values
4299     numberNonZero=0;
4300     for (int i = 0; i < saveN; i++) {
4301          int iColumn = index[i];
4302          double value = array[iColumn];
4303          array[iColumn] =0.0;
4304          if (fabs(value) > tolerance) {
4305            output[numberNonZero] = value;
4306            index[numberNonZero++] = iColumn;
4307          }
4308     }
4309     return numberNonZero;
4310}
4311/* Given positive integer weights for each row fills in sum of weights
4312   for each column (and slack).
4313   Returns weights vector
4314*/
4315CoinBigIndex *
4316ClpPackedMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
4317{
4318     int numberRows = model->numberRows();
4319     int numberColumns = matrix_->getNumCols();
4320     int number = numberRows + numberColumns;
4321     CoinBigIndex * weights = new CoinBigIndex[number];
4322     // get matrix data pointers
4323     const int * row = matrix_->getIndices();
4324     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4325     const int * columnLength = matrix_->getVectorLengths();
4326     int i;
4327     for (i = 0; i < numberColumns; i++) {
4328          CoinBigIndex j;
4329          CoinBigIndex count = 0;
4330          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4331               int iRow = row[j];
4332               count += inputWeights[iRow];
4333          }
4334          weights[i] = count;
4335     }
4336     for (i = 0; i < numberRows; i++) {
4337          weights[i+numberColumns] = inputWeights[i];
4338     }
4339     return weights;
4340}
4341/* Returns largest and smallest elements of both signs.
4342   Largest refers to largest absolute value.
4343*/
4344void
4345ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
4346                                 double & smallestPositive, double & largestPositive)
4347{
4348     smallestNegative = -COIN_DBL_MAX;
4349     largestNegative = 0.0;
4350     smallestPositive = COIN_DBL_MAX;
4351     largestPositive = 0.0;
4352     // get matrix data pointers
4353     const double * elementByColumn = matrix_->getElements();
4354     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4355     const int * columnLength = matrix_->getVectorLengths();
4356     int numberColumns = matrix_->getNumCols();
4357     int i;
4358     for (i = 0; i < numberColumns; i++) {
4359          CoinBigIndex j;
4360          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4361               double value = elementByColumn[j];
4362               if (value > 0.0) {
4363                    smallestPositive = CoinMin(smallestPositive, value);
4364                    largestPositive = CoinMax(largestPositive, value);
4365               } else if (value < 0.0) {
4366                    smallestNegative = CoinMax(smallestNegative, value);
4367                    largestNegative = CoinMin(largestNegative, value);
4368               }
4369          }
4370     }
4371}
4372// Says whether it can do partial pricing
4373bool
4374ClpPackedMatrix::canDoPartialPricing() const
4375{
4376     return true;
4377}
4378// Partial pricing
4379void
4380ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
4381                                int & bestSequence, int & numberWanted)
4382{
4383     numberWanted = currentWanted_;
4384     int start = static_cast<int> (startFraction * numberActiveColumns_);
4385     int end = CoinMin(static_cast<int> (endFraction * numberActiveColumns_ + 1), numberActiveColumns_);
4386     const double * element = matrix_->getElements();
4387     const int * row = matrix_->getIndices();
4388     const CoinBigIndex * startColumn = matrix_->getVectorStarts();
4389     const int * length = matrix_->getVectorLengths();
4390     const double * rowScale = model->rowScale();
4391     const double * columnScale = model->columnScale();
4392     int iSequence;
4393     CoinBigIndex j;
4394     double tolerance = model->currentDualTolerance();
4395     double * reducedCost = model->djRegion();
4396     const double * duals = model->dualRowSolution();
4397     const double * cost = model->costRegion();
4398     double bestDj;
4399     if (bestSequence >= 0)
4400          bestDj = fabs(model->clpMatrix()->reducedCost(model, bestSequence));
4401     else
4402          bestDj = tolerance;
4403     int sequenceOut = model->sequenceOut();
4404     int saveSequence = bestSequence;
4405     int lastScan = minimumObjectsScan_ < 0 ? end : start + minimumObjectsScan_;
4406     int minNeg = minimumGoodReducedCosts_ == -1 ? numberWanted : minimumGoodReducedCosts_;
4407     if (rowScale) {
4408          // scaled
4409          for (iSequence = start; iSequence < end; iSequence++) {
4410               if (iSequence != sequenceOut) {
4411                    double value;
4412                    ClpSimplex::Status status = model->getStatus(iSequence);
4413
4414                    switch(status) {
4415
4416                    case ClpSimplex::basic:
4417                    case ClpSimplex::isFixed:
4418                         break;
4419                    case ClpSimplex::isFree:
4420                    case ClpSimplex::superBasic:
4421                         value = 0.0;
4422                         // scaled
4423                         for (j = startColumn[iSequence];
4424                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4425                              int jRow = row[j];
4426                              value -= duals[jRow] * element[j] * rowScale[jRow];
4427                         }
4428                         value = fabs(cost[iSequence] + value * columnScale[iSequence]);
4429                         if (value > FREE_ACCEPT * tolerance) {
4430                              numberWanted--;
4431                              // we are going to bias towards free (but only if reasonable)
4432                              value *= FREE_BIAS;
4433                              if (value > bestDj) {
4434                                   // check flagged variable and correct dj
4435                                   if (!model->flagged(iSequence)) {
4436                                        bestDj = value;
4437                                        bestSequence = iSequence;
4438                                   } else {
4439                                        // just to make sure we don't exit before got something
4440                                        numberWanted++;
4441                                   }
4442                              }
4443                         }
4444                         break;
4445                    case ClpSimplex::atUpperBound:
4446                         value = 0.0;
4447                         // scaled
4448                         for (j = startColumn[iSequence];
4449                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4450                              int jRow = row[j];
4451                              value -= duals[jRow] * element[j] * rowScale[jRow];
4452                         }
4453                         value = cost[iSequence] + value * columnScale[iSequence];
4454                         if (value > tolerance) {
4455                              numberWanted--;
4456                              if (value > bestDj) {
4457                                   // check flagged variable and correct dj
4458                                   if (!model->flagged(iSequence)) {
4459                                        bestDj = value;
4460                                        bestSequence = iSequence;
4461                                   } else {
4462                                        // just to make sure we don't exit before got something
4463                                        numberWanted++;
4464                                   }
4465                              }
4466                         }
4467                         break;
4468                    case ClpSimplex::atLowerBound:
4469                         value = 0.0;
4470                         // scaled
4471                         for (j = startColumn[iSequence];
4472                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4473                              int jRow = row[j];
4474                              value -= duals[jRow] * element[j] * rowScale[jRow];
4475                         }
4476                         value = -(cost[iSequence] + value * columnScale[iSequence]);
4477                         if (value > tolerance) {
4478                              numberWanted--;
4479                              if (value > bestDj) {
4480                                   // check flagged variable and correct dj
4481                                   if (!model->flagged(iSequence)) {
4482                                        bestDj = value;
4483                                        bestSequence = iSequence;
4484                                   } else {
4485                                        // just to make sure we don't exit before got something
4486                                        numberWanted++;
4487                                   }
4488                              }
4489                         }
4490                         break;
4491                    }
4492               }
4493               if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4494                    // give up
4495                    break;
4496               }
4497               if (!numberWanted)
4498                    break;
4499          }
4500          if (bestSequence != saveSequence) {
4501               // recompute dj
4502               double value = 0.0;
4503               // scaled
4504               for (j = startColumn[bestSequence];
4505                         j < startColumn[bestSequence] + length[bestSequence]; j++) {
4506                    int jRow = row[j];
4507                    value -= duals[jRow] * element[j] * rowScale[jRow];
4508               }
4509               reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
4510               savedBestSequence_ = bestSequence;
4511               savedBestDj_ = reducedCost[savedBestSequence_];
4512          }
4513     } else {
4514          // not scaled
4515          for (iSequence = start; iSequence < end; iSequence++) {
4516               if (iSequence != sequenceOut) {
4517                    double value;
4518                    ClpSimplex::Status status = model->getStatus(iSequence);
4519
4520                    switch(status) {
4521
4522                    case ClpSimplex::basic:
4523                    case ClpSimplex::isFixed:
4524                         break;
4525                    case ClpSimplex::isFree:
4526                    case ClpSimplex::superBasic:
4527                         value = cost[iSequence];
4528                         for (j = startColumn[iSequence];
4529                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4530                              int jRow = row[j];
4531                              value -= duals[jRow] * element[j];
4532                         }
4533                         value = fabs(value);
4534                         if (value > FREE_ACCEPT * tolerance) {
4535                              numberWanted--;
4536                              // we are going to bias towards free (but only if reasonable)
4537                              value *= FREE_BIAS;
4538                              if (value > bestDj) {
4539                                   // check flagged variable and correct dj
4540                                   if (!model->flagged(iSequence)) {
4541                                        bestDj = value;
4542                                        bestSequence = iSequence;
4543                                   } else {
4544                                        // just to make sure we don't exit before got something
4545                                        numberWanted++;
4546                                   }
4547                              }
4548                         }
4549                         break;
4550                    case ClpSimplex::atUpperBound:
4551                         value = cost[iSequence];
4552                         // scaled
4553                         for (j = startColumn[iSequence];
4554                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4555                              int jRow = row[j];
4556                              value -= duals[jRow] * element[j];
4557                         }
4558                         if (value > tolerance) {
4559                              numberWanted--;
4560                              if (value > bestDj) {
4561                                   // check flagged variable and correct dj
4562                                   if (!model->flagged(iSequence)) {
4563                                        bestDj = value;
4564                                        bestSequence = iSequence;
4565                                   } else {
4566                                        // just to make sure we don't exit before got something
4567                                        numberWanted++;
4568                                   }
4569                              }
4570                         }
4571                         break;
4572                    case ClpSimplex::atLowerBound:
4573                         value = cost[iSequence];
4574                         for (j = startColumn[iSequence];
4575                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4576                              int jRow = row[j];
4577                              value -= duals[jRow] * element[j];
4578                         }
4579                         value = -value;
4580                         if (value > tolerance) {
4581                              numberWanted--;
4582                              if (value > bestDj) {
4583                                   // check flagged variable and correct dj
4584                                   if (!model->flagged(iSequence)) {
4585                                        bestDj = value;
4586                                        bestSequence = iSequence;
4587                                   } else {
4588                                        // just to make sure we don't exit before got something
4589                                        numberWanted++;
4590                                   }
4591                              }
4592                         }
4593                         break;
4594                    }
4595               }
4596               if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4597                    // give up
4598                    break;
4599               }
4600               if (!numberWanted)
4601                    break;
4602          }
4603          if (bestSequence != saveSequence) {
4604               // recompute dj
4605               double value = cost[bestSequence];
4606               for (j = startColumn[bestSequence];
4607                         j < startColumn[bestSequence] + length[bestSequence]; j++) {
4608                    int jRow = row[j];
4609                    value -= duals[jRow] * element[j];
4610               }
4611               reducedCost[bestSequence] = value;
4612               savedBestSequence_ = bestSequence;
4613               savedBestDj_ = reducedCost[savedBestSequence_];
4614          }
4615     }
4616     currentWanted_ = numberWanted;
4617}
4618// Sets up an effective RHS
4619void
4620ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
4621{
4622     delete [] rhsOffset_;
4623     int numberRows = model->numberRows();
4624     rhsOffset_ = new double[numberRows];
4625     rhsOffset(model, true);
4626}
4627// Gets rid of special copies
4628void
4629ClpPackedMatrix::clearCopies()
4630{
4631     delete rowCopy_;
4632     delete columnCopy_;
4633     rowCopy_ = NULL;
4634     columnCopy_ = NULL;
4635     flags_ &= ~(4 + 8);
4636     checkGaps();
4637}
4638// makes sure active columns correct
4639int
4640ClpPackedMatrix::refresh(ClpSimplex * )
4641{
4642     numberActiveColumns_ = matrix_->getNumCols();
4643#if 0
4644     ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
4645     ClpPackedMatrix* rowCopy =
4646          dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4647     // Make sure it is really a ClpPackedMatrix
4648     assert (rowCopy != NULL);
4649
4650     const int * column = rowCopy->matrix_->getIndices();
4651     const CoinBigIndex * rowStart = rowCopy->matrix_->getVectorStarts();
4652     const int * rowLength = rowCopy->matrix_->getVectorLengths();
4653     const double * element = rowCopy->matrix_->getElements();
4654     int numberRows = rowCopy->matrix_->getNumRows();
4655     for (int i = 0; i < numberRows; i++) {
4656          if (!rowLength[i])
4657               printf("zero row %d\n", i);
4658     }
4659     delete rowCopy;
4660#endif
4661     return 0;
4662}
4663
4664/* Scales rowCopy if column copy scaled
4665   Only called if scales already exist */
4666void
4667ClpPackedMatrix::scaleRowCopy(ClpModel * model) const
4668{
4669     if (model->rowCopy()) {
4670          // need to replace row by row
4671          int numberRows = model->numberRows();
4672#ifndef NDEBUG
4673          int numberColumns = matrix_->getNumCols();
4674#endif
4675          ClpMatrixBase * rowCopyBase = model->rowCopy();
4676#ifndef NDEBUG
4677          ClpPackedMatrix* rowCopy =
4678               dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4679          // Make sure it is really a ClpPackedMatrix
4680          assert (rowCopy != NULL);
4681#else
4682          ClpPackedMatrix* rowCopy =
4683               static_cast< ClpPackedMatrix*>(rowCopyBase);
4684#endif
4685
4686          const int * column = rowCopy->getIndices();
4687          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4688          double * element = rowCopy->getMutableElements();
4689          const double * rowScale = model->rowScale();
4690          const double * columnScale = model->columnScale();
4691          // scale row copy
4692          for (int iRow = 0; iRow < numberRows; iRow++) {
4693               CoinBigIndex j;
4694               double scale = rowScale[iRow];
4695               double * elementsInThisRow = element + rowStart[iRow];
4696               const int * columnsInThisRow = column + rowStart[iRow];
4697               int number = rowStart[iRow+1] - rowStart[iRow];
4698               assert (number <= numberColumns);
4699               for (j = 0; j < number; j++) {
4700                    int iColumn = columnsInThisRow[j];
4701                    elementsInThisRow[j] *= scale * columnScale[iColumn];
4702               }
4703          }
4704     }
4705}
4706/* Realy really scales column copy
4707   Only called if scales already exist.
4708   Up to user ro delete */
4709ClpMatrixBase *
4710ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const
4711{
4712     // need to replace column by column
4713#ifndef NDEBUG
4714     int numberRows = model->numberRows();
4715#endif
4716     int numberColumns = matrix_->getNumCols();
4717     ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
4718     const int * row = copy->getIndices();
4719     const CoinBigIndex * columnStart = copy->getVectorStarts();
4720     const int * length = copy->getVectorLengths();
4721     double * element = copy->getMutableElements();
4722     const double * rowScale = model->rowScale();
4723     const double * columnScale = model->columnScale();
4724     // scale column copy
4725     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4726          CoinBigIndex j;
4727          double scale = columnScale[iColumn];
4728          double * elementsInThisColumn = element + columnStart[iColumn];
4729          const int * rowsInThisColumn = row + columnStart[iColumn];
4730          int number = length[iColumn];
4731          assert (number <= numberRows);
4732          for (j = 0; j < number; j++) {
4733               int iRow = rowsInThisColumn[j];
4734               elementsInThisColumn[j] *= scale * rowScale[iRow];
4735          }
4736     }
4737     return copy;
4738}
4739// Really scale matrix
4740void
4741ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
4742{
4743     clearCopies();
4744     int numberColumns = matrix_->getNumCols();
4745     const int * row = matrix_->getIndices();
4746     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4747     const int * length = matrix_->getVectorLengths();
4748     double * element = matrix_->getMutableElements();
4749     // scale
4750     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4751          CoinBigIndex j;
4752          double scale = columnScale[iColumn];
4753          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length[iColumn]; j++) {
4754               int iRow = row[j];
4755               element[j] *= scale * rowScale[iRow];
4756          }
4757     }
4758}
4759/* Delete the columns whose indices are listed in <code>indDel</code>. */
4760void
4761ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
4762{
4763     if (matrix_->getNumCols())
4764          matrix_->deleteCols(numDel, indDel);
4765     clearCopies();
4766     numberActiveColumns_ = matrix_->getNumCols();
4767     // may now have gaps
4768     checkGaps();
4769     matrix_->setExtraGap(0.0);
4770}
4771/* Delete the rows whose indices are listed in <code>indDel</code>. */
4772void
4773ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
4774{
4775     if (matrix_->getNumRows())
4776          matrix_->deleteRows(numDel, indDel);
4777     clearCopies();
4778     numberActiveColumns_ = matrix_->getNumCols();
4779     // may now have gaps
4780     checkGaps();
4781     matrix_->setExtraGap(0.0);
4782}
4783#ifndef CLP_NO_VECTOR
4784// Append Columns
4785void
4786ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
4787{
4788     matrix_->appendCols(number, columns);
4789     numberActiveColumns_ = matrix_->getNumCols();
4790     clearCopies();
4791}
4792// Append Rows
4793void
4794ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
4795{
4796     matrix_->appendRows(number, rows);
4797     numberActiveColumns_ = matrix_->getNumCols();
4798     // may now have gaps
4799     checkGaps();
4800     clearCopies();
4801}
4802#endif
4803/* Set the dimensions of the matrix. In effect, append new empty
4804   columns/rows to the matrix. A negative number for either dimension
4805   means that that dimension doesn't change. Otherwise the new dimensions
4806   MUST be at least as large as the current ones otherwise an exception
4807   is thrown. */
4808void
4809ClpPackedMatrix::setDimensions(int numrows, int numcols)
4810{
4811     matrix_->setDimensions(numrows, numcols);
4812}
4813/* Append a set of rows/columns to the end of the matrix. Returns number of errors
4814   i.e. if any of the new rows/columns contain an index that's larger than the
4815   number of columns-1/rows-1 (if numberOther>0) or duplicates
4816   If 0 then rows, 1 if columns */
4817int
4818ClpPackedMatrix::appendMatrix(int number, int type,
4819                              const CoinBigIndex * starts, const int * index,
4820                              const double * element, int numberOther)
4821{
4822     int numberErrors = 0;
4823     // make sure other dimension is big enough
4824     if (type == 0) {
4825          // rows
4826          if (matrix_->isColOrdered() && numberOther > matrix_->getNumCols())
4827               matrix_->setDimensions(-1, numberOther);
4828          if (!matrix_->isColOrdered() || numberOther >= 0 || matrix_->getExtraGap()) {
4829               numberErrors = matrix_->appendRows(number, starts, index, element, numberOther);
4830          } else {
4831               //CoinPackedMatrix mm(*matrix_);
4832               matrix_->appendMinorFast(number, starts, index, element);
4833               //mm.appendRows(number,starts,index,element,numberOther);
4834               //if (!mm.isEquivalent(*matrix_)) {
4835               //printf("bad append\n");
4836               //abort();
4837               //}
4838          }
4839     } else {
4840          // columns
4841          if (!matrix_->isColOrdered() && numberOther > matrix_->getNumRows())
4842               matrix_->setDimensions(numberOther, -1);
4843          numberErrors = matrix_->appendCols(number, starts, index, element, numberOther);
4844     }
4845     clearCopies();
4846     numberActiveColumns_ = matrix_->getNumCols();
4847     return numberErrors;
4848}
4849void
4850ClpPackedMatrix::specialRowCopy(ClpSimplex * model, const ClpMatrixBase * rowCopy)
4851{
4852     delete rowCopy_;
4853     rowCopy_ = new ClpPackedMatrix2(model, rowCopy->getPackedMatrix());
4854     // See if anything in it
4855     if (!rowCopy_->usefulInfo()) {
4856          delete rowCopy_;
4857          rowCopy_ = NULL;
4858          flags_ &= ~4;
4859     } else {
4860          flags_ |= 4;
4861     }
4862}
4863void
4864ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
4865{
4866     delete columnCopy_;
4867     if ((flags_ & 16) != 0) {
4868          columnCopy_ = new ClpPackedMatrix3(model, matrix_);
4869          flags_ |= 8;
4870     } else {
4871          columnCopy_ = NULL;
4872     }
4873}
4874// Say we don't want special column copy
4875void
4876ClpPackedMatrix::releaseSpecialColumnCopy()
4877{
4878     flags_ &= ~(8 + 16);
4879     delete columnCopy_;
4880     columnCopy_ = NULL;
4881}
4882// Correct sequence in and out to give true value
4883void
4884ClpPackedMatrix::correctSequence(const ClpSimplex * model, int & sequenceIn, int & sequenceOut)
4885{
4886     if (columnCopy_) {
4887          if (sequenceIn != -999) {
4888               if (sequenceIn != sequenceOut) {
4889                    if (sequenceIn < numberActiveColumns_)
4890                         columnCopy_->swapOne(model, this, sequenceIn);
4891                    if (sequenceOut < numberActiveColumns_)
4892                         columnCopy_->swapOne(model, this, sequenceOut);
4893               }
4894          } else {
4895               // do all
4896               columnCopy_->sortBlocks(model);
4897          }
4898     }
4899}
4900// Check validity
4901void
4902ClpPackedMatrix::checkFlags(int type) const
4903{
4904     int iColumn;
4905     // get matrix data pointers
4906     //const int * row = matrix_->getIndices();
4907     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4908     const int * columnLength = matrix_->getVectorLengths();
4909     const double * elementByColumn = matrix_->getElements();
4910     if (!zeros()) {
4911          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4912               CoinBigIndex j;
4913               for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4914                    if (!elementByColumn[j])
4915                         abort();
4916               }
4917          }
4918     }
4919     if ((flags_ & 2) == 0) {
4920          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4921               if (columnStart[iColumn+1] != columnStart[iColumn] + columnLength[iColumn]) {
4922                    abort();
4923               }
4924          }
4925     }
4926     if (type) {
4927          if ((flags_ & 2) != 0) {
4928               bool ok = true;
4929               for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4930                    if (columnStart[iColumn+1] != columnStart[iColumn] + columnLength[iColumn]) {
4931                         ok = false;
4932                         break;
4933                    }
4934               }
4935               if (ok)
4936                 COIN_DETAIL_PRINT(printf("flags_ could be 0\n"));
4937          }
4938     }
4939}
4940//#############################################################################
4941// Constructors / Destructor / Assignment
4942//#############################################################################
4943
4944//-------------------------------------------------------------------
4945// Default Constructor
4946//-------------------------------------------------------------------
4947ClpPackedMatrix2::ClpPackedMatrix2 ()
4948     : numberBlocks_(0),
4949       numberRows_(0),
4950       offset_(NULL),
4951       count_(NULL),
4952       rowStart_(NULL),
4953       column_(NULL),
4954       work_(NULL)
4955{
4956#ifdef THREAD
4957     threadId_ = NULL;
4958     info_ = NULL;
4959#endif
4960}
4961//-------------------------------------------------------------------
4962// Useful Constructor
4963//-------------------------------------------------------------------
4964ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * , const CoinPackedMatrix * rowCopy)
4965     : numberBlocks_(0),
4966       numberRows_(0),
4967       offset_(NULL),
4968       count_(NULL),
4969       rowStart_(NULL),
4970       column_(NULL),
4971       work_(NULL)
4972{
4973#ifdef THREAD
4974     threadId_ = NULL;
4975     info_ = NULL;
4976#endif
4977     numberRows_ = rowCopy->getNumRows();
4978     if (!numberRows_)
4979          return;
4980     int numberColumns = rowCopy->getNumCols();
4981     const int * column = rowCopy->getIndices();
4982     const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4983     const int * length = rowCopy->getVectorLengths();
4984     const double * element = rowCopy->getElements();
4985     int chunk = 32768; // tune
4986     //chunk=100;
4987     // tune
4988#if 0
4989     int chunkY[7] = {1024, 2048, 4096, 8192, 16384, 32768, 65535};
4990     int its = model->maximumIterations();
4991     if (its >= 1000000 && its < 1000999) {
4992          its -= 1000000;
4993          its = its / 10;
4994          if (its >= 7) abort();
4995          chunk = chunkY[its];
4996          printf("chunk size %d\n", chunk);
4997#define cpuid(func,ax,bx,cx,dx)\
4998    __asm__ __volatile__ ("cpuid":\
4999    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
5000          unsigned int a, b, c, d;
5001          int func = 0;
5002          cpuid(func, a, b, c, d);
5003          {
5004               int i;
5005               unsigned int value;
5006               value = b;
5007               for (i = 0; i < 4; i++) {
5008                    printf("%c", (value & 0xff));
5009                    value = value >> 8;
5010               }
5011               value = d;
5012               for (i = 0; i < 4; i++) {
5013                    printf("%c", (value & 0xff));
5014                    value = value >> 8;
5015               }
5016               value = c;
5017               for (i = 0; i < 4; i++) {
5018                    printf("%c", (value & 0xff));
5019                    value = value >> 8;
5020               }
5021               printf("\n");
5022               int maxfunc = a;
5023               if (maxfunc > 10) {
5024                    printf("not intel?\n");
5025                    abort();
5026               }
5027               for (func = 1; func <= maxfunc; func++) {
5028                    cpuid(func, a, b, c, d);
5029                    printf("func %d, %x %x %x %x\n", func, a, b, c, d);
5030               }
5031          }
5032#else
5033     if (numberColumns > 10000 || chunk == 100) {
5034#endif
5035     } else {
5036          //printf("no chunk\n");
5037          return;
5038     }
5039     // Could also analyze matrix to get natural breaks
5040     numberBlocks_ = (numberColumns + chunk - 1) / chunk;
5041#ifdef THREAD
5042     // Get work areas
5043     threadId_ = new pthread_t [numberBlocks_];
5044     info_ = new dualColumn0Struct[numberBlocks_];
5045#endif
5046     // Even out
5047     chunk = (numberColumns + numberBlocks_ - 1) / numberBlocks_;
5048     offset_ = new int[numberBlocks_+1];
5049     offset_[numberBlocks_] = numberColumns;
5050     int nRow = numberBlocks_ * numberRows_;
5051     count_ = new unsigned short[nRow];
5052     memset(count_, 0, nRow * sizeof(unsigned short));
5053     rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
5054     CoinBigIndex nElement = rowStart[numberRows_];
5055     rowStart_[nRow+numberRows_] = nElement;
5056     column_ = new unsigned short [nElement];
5057     // assumes int <= double
5058     int sizeWork = 6 * numberBlocks_;
5059     work_ = new double[sizeWork];;
5060     int iBlock;
5061     int nZero = 0;
5062     for (iBlock = 0; iBlock < numberBlocks_; iBlock++) {
5063          int start = iBlock * chunk;
5064          offset_[iBlock] = start;
5065          int end = start + chunk;
5066          for (int iRow = 0; iRow < numberRows_; iRow++) {
5067               if (rowStart[iRow+1] != rowStart[iRow] + length[iRow]) {
5068                    printf("not packed correctly - gaps\n");
5069                    abort();
5070               }
5071               bool lastFound = false;
5072               int nFound = 0;
5073               for (CoinBigIndex j = rowStart[iRow];
5074                         j < rowStart[iRow] + length[iRow]; j++) {
5075                    int iColumn = column[j];
5076                    if (iColumn >= start) {
5077                         if (iColumn < end) {
5078                              if (!element[j]) {
5079                                   printf("not packed correctly - zero element\n");
5080                                   abort();
5081                              }
5082                              column_[j] = static_cast<unsigned short>(iColumn - start);
5083                              nFound++;
5084                              if (lastFound) {
5085                                   printf("not packed correctly - out of order\n");
5086                                   abort();
5087                              }
5088                         } else {
5089                              //can't find any more
5090                              lastFound = true;
5091                         }
5092                    }
5093               }
5094               count_[iRow*numberBlocks_+iBlock] = static_cast<unsigned short>(nFound);
5095               if (!nFound)
5096                    nZero++;
5097          }
5098     }
5099     //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
5100     //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
5101}
5102
5103//-------------------------------------------------------------------
5104// Copy constructor
5105//-------------------------------------------------------------------
5106ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs)
5107     : numberBlocks_(rhs.numberBlocks_),
5108       numberRows_(rhs.numberRows_)
5109{
5110     if (numberBlocks_) {
5111          offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + 1);
5112          int nRow = numberBlocks_ * numberRows_;
5113          count_ = CoinCopyOfArray(rhs.count_, nRow);
5114          rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + 1);
5115          CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5116          column_ = CoinCopyOfArray(rhs.column_, nElement);
5117          int sizeWork = 6 * numberBlocks_;
5118          work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5119#ifdef THREAD
5120          threadId_ = new pthread_t [numberBlocks_];
5121          info_ = new dualColumn0Struct[numberBlocks_];
5122#endif
5123     } else {
5124          offset_ = NULL;
5125          count_ = NULL;
5126          rowStart_ = NULL;
5127          column_ = NULL;
5128          work_ = NULL;
5129#ifdef THREAD
5130          threadId_ = NULL;
5131          info_ = NULL;
5132#endif
5133     }
5134}
5135//-------------------------------------------------------------------
5136// Destructor
5137//-------------------------------------------------------------------
5138ClpPackedMatrix2::~ClpPackedMatrix2 ()
5139{
5140     delete [] offset_;
5141     delete [] count_;
5142     delete [] rowStart_;
5143     delete [] column_;
5144     delete [] work_;
5145#ifdef THREAD
5146     delete [] threadId_;
5147     delete [] info_;
5148#endif
5149}
5150
5151//----------------------------------------------------------------
5152// Assignment operator
5153//-------------------------------------------------------------------
5154ClpPackedMatrix2 &
5155ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
5156{
5157     if (this != &rhs) {
5158          numberBlocks_ = rhs.numberBlocks_;
5159          numberRows_ = rhs.numberRows_;
5160          delete [] offset_;
5161          delete [] count_;
5162          delete [] rowStart_;
5163          delete [] column_;
5164          delete [] work_;
5165#ifdef THREAD
5166          delete [] threadId_;
5167          delete [] info_;
5168#endif
5169          if (numberBlocks_) {
5170               offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + 1);
5171               int nRow = numberBlocks_ * numberRows_;
5172               count_ = CoinCopyOfArray(rhs.count_, nRow);
5173               rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + 1);
5174               CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5175               column_ = CoinCopyOfArray(rhs.column_, nElement);
5176               int sizeWork = 6 * numberBlocks_;
5177               work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5178#ifdef THREAD
5179               threadId_ = new pthread_t [numberBlocks_];
5180               info_ = new dualColumn0Struct[numberBlocks_];
5181#endif
5182          } else {
5183               offset_ = NULL;
5184               count_ = NULL;
5185               rowStart_ = NULL;
5186               column_ = NULL;
5187               work_ = NULL;
5188#ifdef THREAD
5189               threadId_ = NULL;
5190               info_ = NULL;
5191#endif
5192          }
5193     }
5194     return *this;
5195}
5196static int dualColumn0(const ClpSimplex * model, double * spare,
5197                       int * spareIndex, const double * arrayTemp,
5198                       const int * indexTemp, int numberIn,
5199                       int offset, double acceptablePivot, double * bestPossiblePtr,
5200                       double * upperThetaPtr, int * posFreePtr, double * freePivotPtr)
5201{
5202     // do dualColumn0
5203     int i;
5204     int numberRemaining = 0;
5205     double bestPossible = 0.0;
5206     double upperTheta = 1.0e31;
5207     double freePivot = acceptablePivot;
5208     int posFree = -1;
5209     const double * reducedCost = model->djRegion(1);
5210     double dualTolerance = model->dualTolerance();
5211     // We can also see if infeasible or pivoting on free
5212     double tentativeTheta = 1.0e25;
5213     for (i = 0; i < numberIn; i++) {
5214          double alpha = arrayTemp[i];
5215          int iSequence = indexTemp[i] + offset;
5216          double oldValue;
5217          double value;
5218          bool keep;
5219
5220          switch(model->getStatus(iSequence)) {
5221
5222          case ClpSimplex::basic:
5223          case ClpSimplex::isFixed:
5224               break;
5225          case ClpSimplex::isFree:
5226          case ClpSimplex::superBasic:
5227               bestPossible = CoinMax(bestPossible, fabs(alpha));
5228               oldValue = reducedCost[iSequence];
5229               // If free has to be very large - should come in via dualRow
5230               if (model->getStatus(iSequence) == ClpSimplex::isFree && fabs(alpha) < 1.0e-3)
5231                    break;
5232               if (oldValue > dualTolerance) {
5233                    keep = true;
5234               } else if (oldValue < -dualTolerance) {
5235                    keep = true;
5236               } else {
5237                    if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5))
5238                         keep = true;
5239                    else
5240                         keep = false;
5241               }
5242               if (keep) {
5243                    // free - choose largest
5244                    if (fabs(alpha) > freePivot) {
5245                         freePivot = fabs(alpha);
5246                         posFree = i;
5247                    }
5248               }
5249               break;
5250          case ClpSimplex::atUpperBound:
5251               oldValue = reducedCost[iSequence];
5252               value = oldValue - tentativeTheta * alpha;
5253               //assert (oldValue<=dualTolerance*1.0001);
5254               if (value > dualTolerance) {
5255                    bestPossible = CoinMax(bestPossible, -alpha);
5256                    value = oldValue - upperTheta * alpha;
5257                    if (value > dualTolerance && -alpha >= acceptablePivot)
5258                         upperTheta = (oldValue - dualTolerance) / alpha;
5259                    // add to list
5260                    spare[numberRemaining] = alpha;
5261                    spareIndex[numberRemaining++] = iSequence;
5262               }
5263               break;
5264          case ClpSimplex::atLowerBound:
5265               oldValue = reducedCost[iSequence];
5266               value = oldValue - tentativeTheta * alpha;
5267               //assert (oldValue>=-dualTolerance*1.0001);
5268               if (value < -dualTolerance) {
5269                    bestPossible = CoinMax(bestPossible, alpha);
5270                    value = oldValue - upperTheta * alpha;
5271                    if (value < -dualTolerance && alpha >= acceptablePivot)
5272                         upperTheta = (oldValue + dualTolerance) / alpha;
5273                    // add to list
5274                    spare[numberRemaining] = alpha;
5275                    spareIndex[numberRemaining++] = iSequence;
5276               }
5277               break;
5278          }
5279     }
5280     *bestPossiblePtr = bestPossible;
5281     *upperThetaPtr = upperTheta;
5282     *freePivotPtr = freePivot;
5283     *posFreePtr = posFree;
5284     return numberRemaining;
5285}
5286static int doOneBlock(double * array, int * index,
5287                      const double * pi, const CoinBigIndex * rowStart, const double * element,
5288                      const unsigned short * column, int numberInRowArray, int numberLook)
5289{
5290     int iWhich = 0;
5291     int nextN = 0;
5292     CoinBigIndex nextStart = 0;
5293     double nextPi = 0.0;
5294     for (; iWhich < numberInRowArray; iWhich++) {
5295          nextStart = rowStart[0];
5296          nextN = rowStart[numberInRowArray] - nextStart;
5297          rowStart++;
5298          if (nextN) {
5299               nextPi = pi[iWhich];
5300               break;
5301          }
5302     }
5303     int i;
5304     while (iWhich < numberInRowArray) {
5305          double value = nextPi;
5306          CoinBigIndex j =  nextStart;
5307          int n = nextN;
5308          // get next
5309          iWhich++;
5310          for (; iWhich < numberInRowArray; iWhich++) {
5311               nextStart = rowStart[0];
5312               nextN = rowStart[numberInRowArray] - nextStart;
5313               rowStart++;
5314               if (nextN) {
5315                    //coin_prefetch_const(element + nextStart);
5316                    nextPi = pi[iWhich];
5317                    break;
5318               }
5319          }
5320          CoinBigIndex end = j + n;
5321          //coin_prefetch_const(element+rowStart_[i+1]);
5322          //coin_prefetch_const(column_+rowStart_[i+1]);
5323          if (n < 100) {
5324               if ((n & 1) != 0) {
5325                    unsigned int jColumn = column[j];
5326                    array[jColumn] -= value * element[j];
5327                    j++;
5328               }
5329               //coin_prefetch_const(column + nextStart);
5330               for (; j < end; j += 2) {
5331                    unsigned int jColumn0 = column[j];
5332                    double value0 = value * element[j];
5333                    unsigned int jColumn1 = column[j+1];
5334                    double value1 = value * element[j+1];
5335                    array[jColumn0] -= value0;
5336                    array[jColumn1] -= value1;
5337               }
5338          } else {
5339               if ((n & 1) != 0) {
5340                    unsigned int jColumn = column[j];
5341                    array[jColumn] -= value * element[j];
5342                    j++;
5343               }
5344               if ((n & 2) != 0) {
5345                    unsigned int jColumn0 = column[j];
5346                    double value0 = value * element[j];
5347                    unsigned int jColumn1 = column[j+1];
5348                    double value1 = value * element[j+1];
5349                    array[jColumn0] -= value0;
5350                    array[jColumn1] -= value1;
5351                    j += 2;
5352               }
5353               if ((n & 4) != 0) {
5354                    unsigned int jColumn0 = column[j];
5355                    double value0 = value * element[j];
5356                    unsigned int jColumn1 = column[j+1];