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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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