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

Last change on this file since 1679 was 1679, checked in by forrest, 10 years ago

Add code to help when bad cuts are added by user i.e. ones with tiny
or zero elements. May not fix all problems as Cbc assumes users check
the cuts they generate and so does not do all checks.

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