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

Last change on this file since 2078 was 2078, checked in by forrest, 6 years ago

make faster on difficult problems

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