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

Last change on this file since 1722 was 1722, checked in by stefan, 10 years ago

adjust to changes in CoinUtils? header files

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