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

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

changes to allow more options and stop on feasible (and a few other things)

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