source: stable/1.15/Clp/src/ClpPackedMatrix.cpp @ 1995

Last change on this file since 1995 was 1957, checked in by forrest, 7 years ago

llow mor adding columns without elements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 275.7 KB
Line 
1/* $Id: ClpPackedMatrix.cpp 1957 2013-05-15 08:58:19Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7
8#include <cstdio>
9
10#include "CoinPragma.hpp"
11#include "CoinIndexedVector.hpp"
12#include "CoinHelperFunctions.hpp"
13//#define THREAD
14
15#include "ClpSimplex.hpp"
16#include "ClpSimplexDual.hpp"
17#include "ClpFactorization.hpp"
18#ifndef SLIM_CLP
19#include "ClpQuadraticObjective.hpp"
20#endif
21// at end to get min/max!
22#include "ClpPackedMatrix.hpp"
23#include "ClpMessage.hpp"
24#ifdef INTEL_MKL
25#include "mkl_spblas.h"
26#endif
27//#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     if (deletedElements) {
3369       matrix_->setNumElements(matrix_->getNumElements()-deletedElements);
3370       flags_ |= 0x02 ;
3371     }
3372     model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL, *model->messagesPointer())
3373               << smallest << largest
3374               << CoinMessageEol;
3375     if (smallest >= 0.5 && largest <= 2.0 && !deletedElements) {
3376          // don't bother scaling
3377          model->messageHandler()->message(CLP_PACKEDSCALE_FORGET, *model->messagesPointer())
3378                    << CoinMessageEol;
3379          if (!arraysExist) {
3380               delete [] rowScale;
3381               delete [] columnScale;
3382          } else {
3383               model->setRowScale(NULL);
3384               model->setColumnScale(NULL);
3385          }
3386          delete [] usefulColumn;
3387          return 1;
3388     } else {
3389#ifdef CLP_INVESTIGATE
3390          if (deletedElements)
3391               printf("DEL_ELS\n");
3392#endif
3393          if (!rowCopyBase) {
3394               // temporary copy
3395               rowCopyBase = reverseOrderedCopy();
3396          } else if (deletedElements) {
3397               rowCopyBase = reverseOrderedCopy();
3398               model->setNewRowCopy(rowCopyBase);
3399          }
3400#ifndef NDEBUG
3401          ClpPackedMatrix* rowCopy =
3402               dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3403          // Make sure it is really a ClpPackedMatrix
3404          assert (rowCopy != NULL);
3405#else
3406          ClpPackedMatrix* rowCopy =
3407               static_cast< ClpPackedMatrix*>(rowCopyBase);
3408#endif
3409
3410          const int * column = rowCopy->getIndices();
3411          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3412          const double * element = rowCopy->getElements();
3413          // need to scale
3414          if (largest > 1.0e13 * smallest) {
3415               // safer to have smaller zero tolerance
3416               double ratio = smallest / largest;
3417               ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3418               double newTolerance = CoinMax(ratio * 0.5, 1.0e-18);
3419               if (simplex->zeroTolerance() > newTolerance)
3420                    simplex->setZeroTolerance(newTolerance);
3421          }
3422          int scalingMethod = model->scalingFlag();
3423          if (scalingMethod == 4) {
3424               // As auto
3425               scalingMethod = 3;
3426          } else if (scalingMethod == 5) {
3427               // As geometric
3428               scalingMethod = 2;
3429          }
3430          double savedOverallRatio = 0.0;
3431          double tolerance = 5.0 * model->primalTolerance();
3432          double overallLargest = -1.0e-20;
3433          double overallSmallest = 1.0e20;
3434          bool finished = false;
3435          // if scalingMethod 3 then may change
3436          bool extraDetails = (model->logLevel() > 2);
3437#if 0
3438          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3439            if (columnUpper[iColumn] >
3440                columnLower[iColumn] + 1.0e-12 && columnLength[iColumn])
3441              assert(usefulColumn[iColumn]!=0);
3442            else
3443              assert(usefulColumn[iColumn]==0);
3444          }
3445#endif
3446          while (!finished) {
3447               int numberPass = 3;
3448               overallLargest = -1.0e-20;
3449               overallSmallest = 1.0e20;
3450               ClpFillN ( rowScale, numberRows, 1.0);
3451               ClpFillN ( columnScale, numberColumns, 1.0);
3452               if (scalingMethod == 1 || scalingMethod == 3) {
3453                    // Maximum in each row
3454                    for (iRow = 0; iRow < numberRows; iRow++) {
3455                         CoinBigIndex j;
3456                         largest = 1.0e-10;
3457                         for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
3458                              int iColumn = column[j];
3459                              if (usefulColumn[iColumn]) {
3460                                   double value = fabs(element[j]);
3461                                   largest = CoinMax(largest, value);
3462                                   assert (largest < 1.0e40);
3463                              }
3464                         }
3465                         rowScale[iRow] = 1.0 / largest;
3466#ifdef COIN_DEVELOP
3467                         if (extraDetails) {
3468                              overallLargest = CoinMax(overallLargest, largest);
3469                              overallSmallest = CoinMin(overallSmallest, largest);
3470                         }
3471#endif
3472                    }
3473               } else {
3474#ifdef USE_OBJECTIVE
3475                    // This will be used to help get scale factors
3476                    double * objective = new double[numberColumns];
3477                    CoinMemcpyN(model->costRegion(1), numberColumns, objective);
3478                    double objScale = 1.0;
3479#endif
3480                    while (numberPass) {
3481                         overallLargest = 0.0;
3482                         overallSmallest = 1.0e50;
3483                         numberPass--;
3484                         // Geometric mean on row scales
3485                         for (iRow = 0; iRow < numberRows; iRow++) {
3486                              CoinBigIndex j;
3487                              largest = 1.0e-50;
3488                              smallest = 1.0e50;
3489                              for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
3490                                   int iColumn = column[j];
3491                                   if (usefulColumn[iColumn]) {
3492                                        double value = fabs(element[j]);
3493                                        value *= columnScale[iColumn];
3494                                        largest = CoinMax(largest, value);
3495                                        smallest = CoinMin(smallest, value);
3496                                   }
3497                              }
3498
3499#ifdef SQRT_ARRAY
3500                              rowScale[iRow] = smallest * largest;
3501#else
3502                              rowScale[iRow] = 1.0 / sqrt(smallest * largest);
3503#endif
3504                              //rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
3505                              if (extraDetails) {
3506                                   overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
3507                                   overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
3508                              }
3509                         }
3510                         if (model->scalingFlag() == 5)
3511                              break; // just scale rows
3512#ifdef SQRT_ARRAY
3513                         doSqrts(rowScale, numberRows);
3514#endif
3515#ifdef USE_OBJECTIVE
3516                         largest = 1.0e-20;
3517                         smallest = 1.0e50;
3518                         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3519                              if (usefulColumn[iColumn]) {
3520                                   double value = fabs(objective[iColumn]);
3521                                   value *= columnScale[iColumn];
3522                                   largest = CoinMax(largest, value);
3523                                   smallest = CoinMin(smallest, value);
3524                              }
3525                         }
3526                         objScale = 1.0 / sqrt(smallest * largest);
3527#endif
3528                         model->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model->messagesPointer())
3529                                   << overallSmallest
3530                                   << overallLargest
3531                                   << CoinMessageEol;
3532                         // skip last column round
3533                         if (numberPass == 1)
3534                              break;
3535                         // Geometric mean on column scales
3536                         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3537                              if (usefulColumn[iColumn]) {
3538                                   CoinBigIndex j;
3539                                   largest = 1.0e-50;
3540                                   smallest = 1.0e50;
3541                                   for (j = columnStart[iColumn];
3542                                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3543                                        iRow = row[j];
3544                                        double value = fabs(elementByColumn[j]);
3545                                        value *= rowScale[iRow];
3546                                        largest = CoinMax(largest, value);
3547                                        smallest = CoinMin(smallest, value);
3548                                   }
3549#ifdef USE_OBJECTIVE
3550                                   if (fabs(objective[iColumn]) > 1.0e-20) {
3551                                        double value = fabs(objective[iColumn]) * objScale;
3552                                        largest = CoinMax(largest, value);
3553                                        smallest = CoinMin(smallest, value);
3554                                   }
3555#endif
3556#ifdef SQRT_ARRAY
3557                                   columnScale[iColumn] = smallest * largest;
3558#else
3559                                   columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
3560#endif
3561                                   //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3562                              }
3563                         }
3564#ifdef SQRT_ARRAY
3565                         doSqrts(columnScale, numberColumns);
3566#endif
3567                    }
3568#ifdef USE_OBJECTIVE
3569                    delete [] objective;
3570                    printf("obj scale %g - use it if you want\n", objScale);
3571#endif
3572               }
3573               // If ranges will make horrid then scale
3574               for (iRow = 0; iRow < numberRows; iRow++) {
3575                    double difference = rowUpper[iRow] - rowLower[iRow];
3576                    double scaledDifference = difference * rowScale[iRow];
3577                    if (scaledDifference > tolerance && scaledDifference < 1.0e-4) {
3578                         // make gap larger
3579                         rowScale[iRow] *= 1.0e-4 / scaledDifference;
3580                         rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
3581                         //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
3582                         // scaledDifference,difference*rowScale[iRow]);
3583                    }
3584               }
3585               // final pass to scale columns so largest is reasonable
3586               // See what smallest will be if largest is 1.0
3587               if (model->scalingFlag() != 5) {
3588                    overallSmallest = 1.0e50;
3589                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3590                         if (usefulColumn[iColumn]) {
3591                              CoinBigIndex j;
3592                              largest = 1.0e-20;
3593                              smallest = 1.0e50;
3594                              for (j = columnStart[iColumn];
3595                                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3596                                   iRow = row[j];
3597                                   double value = fabs(elementByColumn[j] * rowScale[iRow]);
3598                                   largest = CoinMax(largest, value);
3599                                   smallest = CoinMin(smallest, value);
3600                              }
3601                              if (overallSmallest * largest > smallest)
3602                                   overallSmallest = smallest / largest;
3603                         }
3604                    }
3605               }
3606               if (scalingMethod == 1 || scalingMethod == 2) {
3607                    finished = true;
3608               } else if (savedOverallRatio == 0.0 && scalingMethod != 4) {
3609                    savedOverallRatio = overallSmallest;
3610                    scalingMethod = 4;
3611               } else {
3612                    assert (scalingMethod == 4);
3613                    if (overallSmallest > 2.0 * savedOverallRatio) {
3614                         finished = true; // geometric was better
3615                         if (model->scalingFlag() == 4) {
3616                              // If in Branch and bound change
3617                              if ((model->specialOptions() & 1024) != 0) {
3618                                   model->scaling(2);
3619                              }
3620                         }
3621                    } else {
3622                         scalingMethod = 1; // redo equilibrium
3623                         if (model->scalingFlag() == 4) {
3624                              // If in Branch and bound change
3625                              if ((model->specialOptions() & 1024) != 0) {
3626                                   model->scaling(1);
3627                              }
3628                         }
3629                    }
3630#if 0
3631                    if (extraDetails) {
3632                         if (finished)
3633                              printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
3634                                     savedOverallRatio, overallSmallest);
3635                         else
3636                              printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
3637                                     savedOverallRatio, overallSmallest);
3638                    }
3639#endif
3640               }
3641          }
3642          //#define RANDOMIZE
3643#ifdef RANDOMIZE
3644          // randomize by up to 10%
3645          for (iRow = 0; iRow < numberRows; iRow++) {
3646               double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3647               rowScale[iRow] *= (1.0 + 0.1 * value);
3648          }
3649#endif
3650          overallLargest = 1.0;
3651          if (overallSmallest < 1.0e-1)
3652               overallLargest = 1.0 / sqrt(overallSmallest);
3653          overallLargest = CoinMin(100.0, overallLargest);
3654          overallSmallest = 1.0e50;
3655          char * usedRow = reinterpret_cast<char *>(inverseRowScale);
3656          memset(usedRow, 0, numberRows);
3657          //printf("scaling %d\n",model->scalingFlag());
3658          if (model->scalingFlag() != 5) {
3659               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3660                    if (columnUpper[iColumn] >
3661                              columnLower[iColumn] + 1.0e-12) {
3662                         //if (usefulColumn[iColumn]) {
3663                         CoinBigIndex j;
3664                         largest = 1.0e-20;
3665                         smallest = 1.0e50;
3666                         for (j = columnStart[iColumn];
3667                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3668                              iRow = row[j];
3669                              usedRow[iRow] = 1;
3670                              double value = fabs(elementByColumn[j] * rowScale[iRow]);
3671                              largest = CoinMax(largest, value);
3672                              smallest = CoinMin(smallest, value);
3673                         }
3674                         columnScale[iColumn] = overallLargest / largest;
3675                         //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3676#ifdef RANDOMIZE
3677                         double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3678                         columnScale[iColumn] *= (1.0 + 0.1 * value);
3679#endif
3680                         double difference = columnUpper[iColumn] - columnLower[iColumn];
3681                         if (difference < 1.0e-5 * columnScale[iColumn]) {
3682                              // make gap larger
3683                              columnScale[iColumn] = difference / 1.0e-5;
3684                              //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
3685                              // scaledDifference,difference*columnScale[iColumn]);
3686                         }
3687                         double value = smallest * columnScale[iColumn];
3688                         if (overallSmallest > value)
3689                              overallSmallest = value;
3690                         //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
3691                    } else {
3692                      assert(columnScale[iColumn] == 1.0);
3693                         //columnScale[iColumn]=1.0;
3694                    }
3695               }
3696               for (iRow = 0; iRow < numberRows; iRow++) {
3697                    if (!usedRow[iRow]) {
3698                         rowScale[iRow] = 1.0;
3699                    }
3700               }
3701          }
3702          model->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model->messagesPointer())
3703                    << overallSmallest
3704                    << overallLargest
3705                    << CoinMessageEol;
3706#if 0
3707          {
3708               for (iRow = 0; iRow < numberRows; iRow++) {
3709                    assert (rowScale[iRow] == rowScale2[iRow]);
3710               }
3711               delete [] rowScale2;
3712               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3713                    assert (columnScale[iColumn] == columnScale2[iColumn]);
3714               }
3715               delete [] columnScale2;
3716          }
3717#endif
3718          if (overallSmallest < 1.0e-13) {
3719               // Change factorization zero tolerance
3720               double newTolerance = CoinMax(1.0e-15 * (overallSmallest / 1.0e-13),
3721                                             1.0e-18);
3722               ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3723               if (simplex->factorization()->zeroTolerance() > newTolerance)
3724                    simplex->factorization()->zeroTolerance(newTolerance);
3725               newTolerance = CoinMax(overallSmallest * 0.5, 1.0e-18);
3726               simplex->setZeroTolerance(newTolerance);
3727          }
3728          delete [] usefulColumn;
3729#ifndef SLIM_CLP
3730          // If quadratic then make symmetric
3731          ClpObjective * obj = model->objectiveAsObject();
3732#ifndef NO_RTTI
3733          ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
3734#else
3735          ClpQuadraticObjective * quadraticObj = NULL;
3736          if (obj->type() == 2)
3737               quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
3738#endif
3739          if (quadraticObj) {
3740               if (!rowCopyBase) {
3741                    // temporary copy
3742                    rowCopyBase = reverseOrderedCopy();
3743               }
3744#ifndef NDEBUG
3745               ClpPackedMatrix* rowCopy =
3746                    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3747               // Make sure it is really a ClpPackedMatrix
3748               assert (rowCopy != NULL);
3749#else
3750               ClpPackedMatrix* rowCopy =
3751                    static_cast< ClpPackedMatrix*>(rowCopyBase);
3752#endif
3753               const int * column = rowCopy->getIndices();
3754               const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3755               CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
3756               int numberXColumns = quadratic->getNumCols();
3757               if (numberXColumns < numberColumns) {
3758                    // we assume symmetric
3759                    int numberQuadraticColumns = 0;
3760                    int i;
3761                    //const int * columnQuadratic = quadratic->getIndices();
3762                    //const int * columnQuadraticStart = quadratic->getVectorStarts();
3763                    const int * columnQuadraticLength = quadratic->getVectorLengths();
3764                    for (i = 0; i < numberXColumns; i++) {
3765                         int length = columnQuadraticLength[i];
3766#ifndef CORRECT_COLUMN_COUNTS
3767                         length = 1;
3768#endif
3769                         if (length)
3770                              numberQuadraticColumns++;
3771                    }
3772                    int numberXRows = numberRows - numberQuadraticColumns;
3773                    numberQuadraticColumns = 0;
3774                    for (i = 0; i < numberXColumns; i++) {
3775                         int length = columnQuadraticLength[i];
3776#ifndef CORRECT_COLUMN_COUNTS
3777                         length = 1;
3778#endif
3779                         if (length) {
3780                              rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
3781                              numberQuadraticColumns++;
3782                         }
3783                    }
3784                    int numberQuadraticRows = 0;
3785                    for (i = 0; i < numberXRows; i++) {
3786                         // See if any in row quadratic
3787                         CoinBigIndex j;
3788                         int numberQ = 0;
3789                         for (j = rowStart[i]; j < rowStart[i+1]; j++) {
3790                              int iColumn = column[j];
3791                              if (columnQuadraticLength[iColumn])
3792                                   numberQ++;
3793                         }
3794#ifndef CORRECT_ROW_COUNTS
3795                         numberQ = 1;
3796#endif
3797                         if (numberQ) {
3798                              columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
3799                              numberQuadraticRows++;
3800                         }
3801                    }
3802                    // and make sure Sj okay
3803                    for (iColumn = numberQuadraticRows + numberXColumns; iColumn < numberColumns; iColumn++) {
3804                         CoinBigIndex j = columnStart[iColumn];
3805                         assert(columnLength[iColumn] == 1);
3806                         int iRow = row[j];
3807                         double value = fabs(elementByColumn[j] * rowScale[iRow]);
3808                         columnScale[iColumn] = 1.0 / value;
3809                    }
3810               }
3811          }
3812#endif
3813          // make copy (could do faster by using previous values)
3814          // could just do partial
3815          for (iRow = 0; iRow < numberRows; iRow++)
3816               inverseRowScale[iRow] = 1.0 / rowScale[iRow] ;
3817          for (iColumn = 0; iColumn < numberColumns; iColumn++)
3818               inverseColumnScale[iColumn] = 1.0 / columnScale[iColumn] ;
3819          if (!arraysExist) {
3820               model->setRowScale(rowScale);
3821               model->setColumnScale(columnScale);
3822          }
3823          if (model->rowCopy()) {
3824               // need to replace row by row
3825               ClpPackedMatrix* rowCopy =
3826                    static_cast< ClpPackedMatrix*>(model->rowCopy());
3827               double * element = rowCopy->getMutableElements();
3828               const int * column = rowCopy->getIndices();
3829               const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3830               // scale row copy
3831               for (iRow = 0; iRow < numberRows; iRow++) {
3832                    CoinBigIndex j;
3833                    double scale = rowScale[iRow];
3834                    double * elementsInThisRow = element + rowStart[iRow];
3835                    const int * columnsInThisRow = column + rowStart[iRow];
3836                    int number = rowStart[iRow+1] - rowStart[iRow];
3837                    assert (number <= numberColumns);
3838                    for (j = 0; j < number; j++) {
3839                         int iColumn = columnsInThisRow[j];
3840                         elementsInThisRow[j] *= scale * columnScale[iColumn];
3841                    }
3842               }
3843               if ((model->specialOptions() & 262144) != 0) {
3844                    //if ((model->specialOptions()&(COIN_CBC_USING_CLP|16384))!=0) {
3845                    //if (model->inCbcBranchAndBound()&&false) {
3846                    // copy without gaps
3847                    CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, 0, 0);
3848                    ClpPackedMatrix * scaled = new ClpPackedMatrix(scaledMatrix);
3849                    model->setClpScaledMatrix(scaled);
3850                    // get matrix data pointers
3851                    const int * row = scaledMatrix->getIndices();
3852                    const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3853#ifndef NDEBUG
3854                    const int * columnLength = scaledMatrix->getVectorLengths();
3855#endif
3856                    double * elementByColumn = scaledMatrix->getMutableElements();
3857                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3858                         CoinBigIndex j;
3859                         double scale = columnScale[iColumn];
3860                         assert (columnStart[iColumn+1] == columnStart[iColumn] + columnLength[iColumn]);
3861                         for (j = columnStart[iColumn];
3862                                   j < columnStart[iColumn+1]; j++) {
3863                              int iRow = row[j];
3864                              elementByColumn[j] *= scale * rowScale[iRow];
3865                         }
3866                    }
3867               } else {
3868                    //printf("not in b&b\n");
3869               }
3870          } else {
3871               // no row copy
3872               delete rowCopyBase;
3873          }
3874          return 0;
3875     }
3876}
3877// Creates scaled column copy if scales exist
3878void
3879ClpPackedMatrix::createScaledMatrix(ClpSimplex * model) const
3880{
3881     int numberRows = model->numberRows();
3882     int numberColumns = matrix_->getNumCols();
3883     model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
3884     // If empty - return as sanityCheck will trap
3885     if (!numberRows || !numberColumns) {
3886          model->setRowScale(NULL);
3887          model->setColumnScale(NULL);
3888          return ;
3889     }
3890     if (!model->rowScale())
3891          return;
3892     double * rowScale = model->mutableRowScale();
3893     double * columnScale = model->mutableColumnScale();
3894     // copy without gaps
3895     CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, 0, 0);
3896     ClpPackedMatrix * scaled = new ClpPackedMatrix(scaledMatrix);
3897     model->setClpScaledMatrix(scaled);
3898     // get matrix data pointers
3899     const int * row = scaledMatrix->getIndices();
3900     const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3901#ifndef NDEBUG
3902     const int * columnLength = scaledMatrix->getVectorLengths();
3903#endif
3904     double * elementByColumn = scaledMatrix->getMutableElements();
3905     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3906          CoinBigIndex j;
3907          double scale = columnScale[iColumn];
3908          assert (columnStart[iColumn+1] == columnStart[iColumn] + columnLength[iColumn]);
3909          for (j = columnStart[iColumn];
3910                    j < columnStart[iColumn+1]; j++) {
3911               int iRow = row[j];
3912               elementByColumn[j] *= scale * rowScale[iRow];
3913          }
3914     }
3915#ifdef DO_CHECK_FLAGS
3916     checkFlags(0);
3917#endif
3918}
3919/* Unpacks a column into an CoinIndexedvector
3920 */
3921void
3922ClpPackedMatrix::unpack(const ClpSimplex * model, CoinIndexedVector * rowArray,
3923                        int iColumn) const
3924{
3925     const double * rowScale = model->rowScale();
3926     const int * row = matrix_->getIndices();
3927     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3928     const int * columnLength = matrix_->getVectorLengths();
3929     const double * elementByColumn = matrix_->getElements();
3930     CoinBigIndex i;
3931     if (!rowScale) {
3932          for (i = columnStart[iColumn];
3933                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3934               rowArray->quickAdd(row[i], elementByColumn[i]);
3935          }
3936     } else {
3937          // apply scaling
3938          double scale = model->columnScale()[iColumn];
3939          for (i = columnStart[iColumn];
3940                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3941               int iRow = row[i];
3942               rowArray->quickAdd(iRow, elementByColumn[i]*scale * rowScale[iRow]);
3943          }
3944     }
3945}
3946/* Unpacks a column into a CoinIndexedVector
3947** in packed format
3948Note that model is NOT const.  Bounds and objective could
3949be modified if doing column generation (just for this variable) */
3950void
3951ClpPackedMatrix::unpackPacked(ClpSimplex * model,
3952                              CoinIndexedVector * rowArray,
3953                              int iColumn) const
3954{
3955     const double * rowScale = model->rowScale();
3956     const int * row = matrix_->getIndices();
3957     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3958     const int * columnLength = matrix_->getVectorLengths();
3959     const double * elementByColumn = matrix_->getElements();
3960     CoinBigIndex i;
3961     int * index = rowArray->getIndices();
3962     double * array = rowArray->denseVector();
3963     int number = 0;
3964     if (!rowScale) {
3965          for (i = columnStart[iColumn];
3966                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3967               int iRow = row[i];
3968               double value = elementByColumn[i];
3969               if (value) {
3970                    array[number] = value;
3971                    index[number++] = iRow;
3972               }
3973          }
3974          rowArray->setNumElements(number);
3975          rowArray->setPackedMode(true);
3976     } else {
3977          // apply scaling
3978          double scale = model->columnScale()[iColumn];
3979          for (i = columnStart[iColumn];
3980                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3981               int iRow = row[i];
3982               double value = elementByColumn[i] * scale * rowScale[iRow];
3983               if (value) {
3984                    array[number] = value;
3985                    index[number++] = iRow;
3986               }
3987          }
3988          rowArray->setNumElements(number);
3989          rowArray->setPackedMode(true);
3990     }
3991}
3992/* Adds multiple of a column into an CoinIndexedvector
3993      You can use quickAdd to add to vector */
3994void
3995ClpPackedMatrix::add(const ClpSimplex * model, CoinIndexedVector * rowArray,
3996                     int iColumn, double multiplier) const
3997{
3998     const double * rowScale = model->rowScale();
3999     const int * row = matrix_->getIndices();
4000     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4001     const int * columnLength = matrix_->getVectorLengths();
4002     const double * elementByColumn = matrix_->getElements();
4003     CoinBigIndex i;
4004     if (!rowScale) {
4005          for (i = columnStart[iColumn];
4006                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
4007               int iRow = row[i];
4008               rowArray->quickAdd(iRow, multiplier * elementByColumn[i]);
4009          }
4010     } else {
4011          // apply scaling
4012          double scale = model->columnScale()[iColumn] * multiplier;
4013          for (i = columnStart[iColumn];
4014                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
4015               int iRow = row[i];
4016               rowArray->quickAdd(iRow, elementByColumn[i]*scale * rowScale[iRow]);
4017          }
4018     }
4019}
4020/* Adds multiple of a column into an array */
4021void
4022ClpPackedMatrix::add(const ClpSimplex * model, double * array,
4023                     int iColumn, double multiplier) const
4024{
4025     const double * rowScale = model->rowScale();
4026     const int * row = matrix_->getIndices();
4027     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4028     const int * columnLength = matrix_->getVectorLengths();
4029     const double * elementByColumn = matrix_->getElements();
4030     CoinBigIndex i;
4031     if (!rowScale) {
4032          for (i = columnStart[iColumn];
4033                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
4034               int iRow = row[i];
4035               array[iRow] += multiplier * elementByColumn[i];
4036          }
4037     } else {
4038          // apply scaling
4039          double scale = model->columnScale()[iColumn] * multiplier;
4040          for (i = columnStart[iColumn];
4041                    i < columnStart[iColumn] + columnLength[iColumn]; i++) {
4042               int iRow = row[i];
4043               array[iRow] += elementByColumn[i] * scale * rowScale[iRow];
4044          }
4045     }
4046}
4047/* Checks if all elements are in valid range.  Can just
4048   return true if you are not paranoid.  For Clp I will
4049   probably expect no zeros.  Code can modify matrix to get rid of
4050   small elements.
4051   check bits (can be turned off to save time) :
4052   1 - check if matrix has gaps
4053   2 - check if zero elements
4054   4 - check and compress duplicates
4055   8 - report on large and small
4056*/
4057bool
4058ClpPackedMatrix::allElementsInRange(ClpModel * model,
4059                                    double smallest, double largest,
4060                                    int check)
4061{
4062     int iColumn;
4063     // make sure matrix correct size
4064     assert (matrix_->getNumRows() <= model->numberRows());
4065     if (model->clpScaledMatrix())
4066          assert (model->clpScaledMatrix()->getNumElements() == matrix_->getNumElements());
4067     assert (matrix_->getNumRows() <= model->numberRows());
4068     matrix_->setDimensions(model->numberRows(), model->numberColumns());
4069     CoinBigIndex numberLarge = 0;;
4070     CoinBigIndex numberSmall = 0;;
4071     CoinBigIndex numberDuplicate = 0;;
4072     int firstBadColumn = -1;
4073     int firstBadRow = -1;
4074     double firstBadElement = 0.0;
4075     // get matrix data pointers
4076     const int * row = matrix_->getIndices();
4077     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4078     const int * columnLength = matrix_->getVectorLengths();
4079     const double * elementByColumn = matrix_->getElements();
4080     int numberRows = model->numberRows();
4081     int numberColumns = matrix_->getNumCols();
4082     // Say no gaps
4083     flags_ &= ~2;
4084     if (type_>=10)
4085       return true; // gub
4086     if (check == 14 || check == 10) {
4087          if (matrix_->getNumElements() < columnStart[numberColumns]) {
4088               // pack down!
4089#if 0
4090               matrix_->removeGaps();
4091#else
4092               checkGaps();
4093#ifdef DO_CHECK_FLAGS
4094               checkFlags(0);
4095#endif
4096#endif
4097#ifdef COIN_DEVELOP
4098               //printf("flags set to 2\n");
4099#endif
4100          } else if (numberColumns) {
4101               assert(columnStart[numberColumns] ==
4102                      columnStart[numberColumns-1] + columnLength[numberColumns-1]);
4103          }
4104          return true;
4105     }
4106     assert (check == 15 || check == 11);
4107     if (check == 15) {
4108          int * mark = new int [numberRows];
4109          int i;
4110          for (i = 0; i < numberRows; i++)
4111               mark[i] = -1;
4112          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4113               CoinBigIndex j;
4114               CoinBigIndex start = columnStart[iColumn];
4115               CoinBigIndex end = start + columnLength[iColumn];
4116               if (end != columnStart[iColumn+1])
4117                    flags_ |= 2;
4118               for (j = start; j < end; j++) {
4119                    double value = fabs(elementByColumn[j]);
4120                    int iRow = row[j];
4121                    if (iRow < 0 || iRow >= numberRows) {
4122#ifndef COIN_BIG_INDEX
4123                         printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4124#elif COIN_BIG_INDEX==0
4125                         printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4126#elif COIN_BIG_INDEX==1
4127                         printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4128#else
4129                         printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4130#endif
4131                         return false;
4132                    }
4133                    if (mark[iRow] == -1) {
4134                         mark[iRow] = j;
4135                    } else {
4136                         // duplicate
4137                         numberDuplicate++;
4138                    }
4139                    //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
4140                    if (!value)
4141                         flags_ |= 1; // there are zero elements
4142                    if (value < smallest) {
4143                         numberSmall++;
4144                    } else if (!(value <= largest)) {
4145                         numberLarge++;
4146                         if (firstBadColumn < 0) {
4147                              firstBadColumn = iColumn;
4148                              firstBadRow = row[j];
4149                              firstBadElement = elementByColumn[j];
4150                         }
4151                    }
4152               }
4153               //clear mark
4154               for (j = columnStart[iColumn];
4155                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4156                    int iRow = row[j];
4157                    mark[iRow] = -1;
4158               }
4159          }
4160          delete [] mark;
4161     } else {
4162          // just check for out of range - not for duplicates
4163          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4164               CoinBigIndex j;
4165               CoinBigIndex start = columnStart[iColumn];
4166               CoinBigIndex end = start + columnLength[iColumn];
4167               if (end != columnStart[iColumn+1])
4168                    flags_ |= 2;
4169               for (j = start; j < end; j++) {
4170                    double value = fabs(elementByColumn[j]);
4171                    int iRow = row[j];
4172                    if (iRow < 0 || iRow >= numberRows) {
4173#ifndef COIN_BIG_INDEX
4174                         printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4175#elif COIN_BIG_INDEX==0
4176                         printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4177#elif COIN_BIG_INDEX==1
4178                         printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4179#else
4180                         printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4181#endif
4182                         return false;
4183                    }
4184                    if (!value)
4185                         flags_ |= 1; // there are zero elements
4186                    if (value < smallest) {
4187                         numberSmall++;
4188                    } else if (!(value <= largest)) {
4189                         numberLarge++;
4190                         if (firstBadColumn < 0) {
4191                              firstBadColumn = iColumn;
4192                              firstBadRow = iRow;
4193                              firstBadElement = value;
4194                         }
4195                    }
4196               }
4197          }
4198     }
4199     if (numberLarge) {
4200          model->messageHandler()->message(CLP_BAD_MATRIX, model->messages())
4201                    << numberLarge
4202                    << firstBadColumn << firstBadRow << firstBadElement
4203                    << CoinMessageEol;
4204          return false;
4205     }
4206     if (numberSmall)
4207          model->messageHandler()->message(CLP_SMALLELEMENTS, model->messages())
4208                    << numberSmall
4209                    << CoinMessageEol;
4210     if (numberDuplicate)
4211          model->messageHandler()->message(CLP_DUPLICATEELEMENTS, model->messages())
4212                    << numberDuplicate
4213                    << CoinMessageEol;
4214     if (numberDuplicate)
4215          matrix_->eliminateDuplicates(smallest);
4216     else if (numberSmall)
4217          matrix_->compress(smallest);
4218     // If smallest >0.0 then there can't be zero elements
4219     if (smallest > 0.0)
4220          flags_ &= ~1;;
4221     if (numberSmall || numberDuplicate)
4222          flags_ |= 2; // will have gaps
4223     return true;
4224}
4225int
4226ClpPackedMatrix::gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector * COIN_RESTRICT piVector,
4227          int * COIN_RESTRICT index,
4228          double * COIN_RESTRICT output,
4229          int * COIN_RESTRICT lookup,
4230          char * COIN_RESTRICT marked,
4231          const double tolerance,
4232          const double scalar) const
4233{
4234     const double * COIN_RESTRICT pi = piVector->denseVector();
4235     int numberNonZero = 0;
4236     int numberInRowArray = piVector->getNumElements();
4237     const int * COIN_RESTRICT column = matrix_->getIndices();
4238     const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4239     const double * COIN_RESTRICT element = matrix_->getElements();
4240     const int * COIN_RESTRICT whichRow = piVector->getIndices();
4241     int * fakeRow = const_cast<int *> (whichRow);
4242     fakeRow[numberInRowArray]=0; // so can touch
4243#ifndef NDEBUG
4244     int maxColumn = 0;
4245#endif
4246     // ** Row copy is already scaled
4247     int nextRow=whichRow[0];
4248     CoinBigIndex nextStart = rowStart[nextRow]; 
4249     CoinBigIndex nextEnd = rowStart[nextRow+1]; 
4250     for (int i = 0; i < numberInRowArray; i++) {
4251          double value = pi[i] * scalar;
4252          CoinBigIndex start=nextStart;
4253          CoinBigIndex end=nextEnd;
4254          nextRow=whichRow[i+1];
4255          nextStart = rowStart[nextRow]; 
4256          //coin_prefetch_const(column + nextStart);
4257          //coin_prefetch_const(element + nextStart);
4258          nextEnd = rowStart[nextRow+1]; 
4259          CoinBigIndex j;
4260          for (j = start; j < end; j++) {
4261               int iColumn = column[j];
4262#ifndef NDEBUG
4263               maxColumn = CoinMax(maxColumn, iColumn);
4264#endif
4265               double elValue = element[j];
4266               elValue *= value;
4267               if (marked[iColumn]) {
4268                    int k = lookup[iColumn];
4269                    output[k] += elValue;
4270               } else {
4271                    output[numberNonZero] = elValue;
4272                    marked[iColumn] = 1;
4273                    lookup[iColumn] = numberNonZero;
4274                    index[numberNonZero++] = iColumn;
4275               }
4276          }
4277     }
4278#ifndef NDEBUG
4279     int saveN = numberNonZero;
4280#endif
4281     // get rid of tiny values and zero out lookup
4282     for (int i = 0; i < numberNonZero; i++) {
4283          int iColumn = index[i];
4284          marked[iColumn] = 0;
4285          double value = output[i];
4286          if (fabs(value) <= tolerance) {
4287               while (fabs(value) <= tolerance) {
4288                    numberNonZero--;
4289                    value = output[numberNonZero];
4290                    iColumn = index[numberNonZero];
4291                    marked[iColumn] = 0;
4292                    if (i < numberNonZero) {
4293                         output[numberNonZero] = 0.0;
4294                         output[i] = value;
4295                         index[i] = iColumn;
4296                    } else {
4297                         output[i] = 0.0;
4298                         value = 1.0; // to force end of while
4299                    }
4300               }
4301          }
4302     }
4303#ifndef NDEBUG
4304     for (int i = numberNonZero; i < saveN; i++)
4305          assert(!output[i]);
4306     for (int i = 0; i <= maxColumn; i++)
4307          assert (!marked[i]);
4308#endif
4309     return numberNonZero;
4310}
4311int
4312ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
4313          int * COIN_RESTRICT index,
4314          double * COIN_RESTRICT output,
4315          double * COIN_RESTRICT array,
4316          const double tolerance,
4317          const double scalar) const
4318{
4319     const double * COIN_RESTRICT pi = piVector->denseVector();
4320     int numberNonZero = 0;
4321     int numberInRowArray = piVector->getNumElements();
4322     const int * COIN_RESTRICT column = matrix_->getIndices();
4323     const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4324     const double * COIN_RESTRICT element = matrix_->getElements();
4325     const int * COIN_RESTRICT whichRow = piVector->getIndices();
4326     // ** Row copy is already scaled
4327     for (int i = 0; i < numberInRowArray; i++) {
4328          int iRow = whichRow[i];
4329          double value = pi[i] * scalar;
4330          CoinBigIndex j;
4331          for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
4332               int iColumn = column[j];
4333               double inValue = array[iColumn];
4334               double elValue = element[j];
4335               elValue *= value;
4336               if (inValue) {
4337                 double outValue = inValue + elValue;
4338                 if (!outValue)
4339                   outValue = COIN_INDEXED_REALLY_TINY_ELEMENT;
4340                 array[iColumn] = outValue;
4341               } else {
4342                 array[iColumn] = elValue;
4343                 assert (elValue);
4344                 index[numberNonZero++] = iColumn;
4345               }
4346          }
4347     }
4348     int saveN = numberNonZero;
4349     // get rid of tiny values
4350     numberNonZero=0;
4351     for (int i = 0; i < saveN; i++) {
4352          int iColumn = index[i];
4353          double value = array[iColumn];
4354          array[iColumn] =0.0;
4355          if (fabs(value) > tolerance) {
4356            output[numberNonZero] = value;
4357            index[numberNonZero++] = iColumn;
4358          }
4359     }
4360     return numberNonZero;
4361}
4362/* Given positive integer weights for each row fills in sum of weights
4363   for each column (and slack).
4364   Returns weights vector
4365*/
4366CoinBigIndex *
4367ClpPackedMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
4368{
4369     int numberRows = model->numberRows();
4370     int numberColumns = matrix_->getNumCols();
4371     int number = numberRows + numberColumns;
4372     CoinBigIndex * weights = new CoinBigIndex[number];
4373     // get matrix data pointers
4374     const int * row = matrix_->getIndices();
4375     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4376     const int * columnLength = matrix_->getVectorLengths();
4377     int i;
4378     for (i = 0; i < numberColumns; i++) {
4379          CoinBigIndex j;
4380          CoinBigIndex count = 0;
4381          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4382               int iRow = row[j];
4383               count += inputWeights[iRow];
4384          }
4385          weights[i] = count;
4386     }
4387     for (i = 0; i < numberRows; i++) {
4388          weights[i+numberColumns] = inputWeights[i];
4389     }
4390     return weights;
4391}
4392/* Returns largest and smallest elements of both signs.
4393   Largest refers to largest absolute value.
4394*/
4395void
4396ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
4397                                 double & smallestPositive, double & largestPositive)
4398{
4399     smallestNegative = -COIN_DBL_MAX;
4400     largestNegative = 0.0;
4401     smallestPositive = COIN_DBL_MAX;
4402     largestPositive = 0.0;
4403     // get matrix data pointers
4404     const double * elementByColumn = matrix_->getElements();
4405     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4406     const int * columnLength = matrix_->getVectorLengths();
4407     int numberColumns = matrix_->getNumCols();
4408     int i;
4409     for (i = 0; i < numberColumns; i++) {
4410          CoinBigIndex j;
4411          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4412               double value = elementByColumn[j];
4413               if (value > 0.0) {
4414                    smallestPositive = CoinMin(smallestPositive, value);
4415                    largestPositive = CoinMax(largestPositive, value);
4416               } else if (value < 0.0) {
4417                    smallestNegative = CoinMax(smallestNegative, value);
4418                    largestNegative = CoinMin(largestNegative, value);
4419               }
4420          }
4421     }
4422}
4423// Says whether it can do partial pricing
4424bool
4425ClpPackedMatrix::canDoPartialPricing() const
4426{
4427     return true;
4428}
4429// Partial pricing
4430void
4431ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
4432                                int & bestSequence, int & numberWanted)
4433{
4434     numberWanted = currentWanted_;
4435     int start = static_cast<int> (startFraction * numberActiveColumns_);
4436     int end = CoinMin(static_cast<int> (endFraction * numberActiveColumns_ + 1), numberActiveColumns_);
4437     const double * element = matrix_->getElements();
4438     const int * row = matrix_->getIndices();
4439     const CoinBigIndex * startColumn = matrix_->getVectorStarts();
4440     const int * length = matrix_->getVectorLengths();
4441     const double * rowScale = model->rowScale();
4442     const double * columnScale = model->columnScale();
4443     int iSequence;
4444     CoinBigIndex j;
4445     double tolerance = model->currentDualTolerance();
4446     double * reducedCost = model->djRegion();
4447     const double * duals = model->dualRowSolution();
4448     const double * cost = model->costRegion();
4449     double bestDj;
4450     if (bestSequence >= 0)
4451          bestDj = fabs(model->clpMatrix()->reducedCost(model, bestSequence));
4452     else
4453          bestDj = tolerance;
4454     int sequenceOut = model->sequenceOut();
4455     int saveSequence = bestSequence;
4456     int lastScan = minimumObjectsScan_ < 0 ? end : start + minimumObjectsScan_;
4457     int minNeg = minimumGoodReducedCosts_ == -1 ? numberWanted : minimumGoodReducedCosts_;
4458     if (rowScale) {
4459          // scaled
4460          for (iSequence = start; iSequence < end; iSequence++) {
4461               if (iSequence != sequenceOut) {
4462                    double value;
4463                    ClpSimplex::Status status = model->getStatus(iSequence);
4464
4465                    switch(status) {
4466
4467                    case ClpSimplex::basic:
4468                    case ClpSimplex::isFixed:
4469                         break;
4470                    case ClpSimplex::isFree:
4471                    case ClpSimplex::superBasic:
4472                         value = 0.0;
4473                         // scaled
4474                         for (j = startColumn[iSequence];
4475                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4476                              int jRow = row[j];
4477                              value -= duals[jRow] * element[j] * rowScale[jRow];
4478                         }
4479                         value = fabs(cost[iSequence] + value * columnScale[iSequence]);
4480                         if (value > FREE_ACCEPT * tolerance) {
4481                              numberWanted--;
4482                              // we are going to bias towards free (but only if reasonable)
4483                              value *= FREE_BIAS;
4484                              if (value > bestDj) {
4485                                   // check flagged variable and correct dj
4486                                   if (!model->flagged(iSequence)) {
4487                                        bestDj = value;
4488                                        bestSequence = iSequence;
4489                                   } else {
4490                                        // just to make sure we don't exit before got something
4491                                        numberWanted++;
4492                                   }
4493                              }
4494                         }
4495                         break;
4496                    case ClpSimplex::atUpperBound:
4497                         value = 0.0;
4498                         // scaled
4499                         for (j = startColumn[iSequence];
4500                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4501                              int jRow = row[j];
4502                              value -= duals[jRow] * element[j] * rowScale[jRow];
4503                         }
4504                         value = cost[iSequence] + value * columnScale[iSequence];
4505                         if (value > tolerance) {
4506                              numberWanted--;
4507                              if (value > bestDj) {
4508                                   // check flagged variable and correct dj
4509                                   if (!model->flagged(iSequence)) {
4510                                        bestDj = value;
4511                                        bestSequence = iSequence;
4512                                   } else {
4513                                        // just to make sure we don't exit before got something
4514                                        numberWanted++;
4515                                   }
4516                              }
4517                         }
4518                         break;
4519                    case ClpSimplex::atLowerBound:
4520                         value = 0.0;
4521                         // scaled
4522                         for (j = startColumn[iSequence];
4523                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4524                              int jRow = row[j];
4525                              value -= duals[jRow] * element[j] * rowScale[jRow];
4526                         }
4527                         value = -(cost[iSequence] + value * columnScale[iSequence]);
4528                         if (value > tolerance) {
4529                              numberWanted--;
4530                              if (value > bestDj) {
4531                                   // check flagged variable and correct dj
4532                                   if (!model->flagged(iSequence)) {
4533                                        bestDj = value;
4534                                        bestSequence = iSequence;
4535                                   } else {
4536                                        // just to make sure we don't exit before got something
4537                                        numberWanted++;
4538                                   }
4539                              }
4540                         }
4541                         break;
4542                    }
4543               }
4544               if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4545                    // give up
4546                    break;
4547               }
4548               if (!numberWanted)
4549                    break;
4550          }
4551          if (bestSequence != saveSequence) {
4552               // recompute dj
4553               double value = 0.0;
4554               // scaled
4555               for (j = startColumn[bestSequence];
4556                         j < startColumn[bestSequence] + length[bestSequence]; j++) {
4557                    int jRow = row[j];
4558                    value -= duals[jRow] * element[j] * rowScale[jRow];
4559               }
4560               reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
4561               savedBestSequence_ = bestSequence;
4562               savedBestDj_ = reducedCost[savedBestSequence_];
4563          }
4564     } else {
4565          // not scaled
4566          for (iSequence = start; iSequence < end; iSequence++) {
4567               if (iSequence != sequenceOut) {
4568                    double value;
4569                    ClpSimplex::Status status = model->getStatus(iSequence);
4570
4571                    switch(status) {
4572
4573                    case ClpSimplex::basic:
4574                    case ClpSimplex::isFixed:
4575                         break;
4576                    case ClpSimplex::isFree:
4577                    case ClpSimplex::superBasic:
4578                         value = cost[iSequence];
4579                         for (j = startColumn[iSequence];
4580                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4581                              int jRow = row[j];
4582                              value -= duals[jRow] * element[j];
4583                         }
4584                         value = fabs(value);
4585                         if (value > FREE_ACCEPT * tolerance) {
4586                              numberWanted--;
4587                              // we are going to bias towards free (but only if reasonable)
4588                              value *= FREE_BIAS;
4589                              if (value > bestDj) {
4590                                   // check flagged variable and correct dj
4591                                   if (!model->flagged(iSequence)) {
4592                                        bestDj = value;
4593                                        bestSequence = iSequence;
4594                                   } else {
4595                                        // just to make sure we don't exit before got something
4596                                        numberWanted++;
4597                                   }
4598                              }
4599                         }
4600                         break;
4601                    case ClpSimplex::atUpperBound:
4602                         value = cost[iSequence];
4603                         // scaled
4604                         for (j = startColumn[iSequence];
4605                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4606                              int jRow = row[j];
4607                              value -= duals[jRow] * element[j];
4608                         }
4609                         if (value > tolerance) {
4610                              numberWanted--;
4611                              if (value > bestDj) {
4612                                   // check flagged variable and correct dj
4613                                   if (!model->flagged(iSequence)) {
4614                                        bestDj = value;
4615                                        bestSequence = iSequence;
4616                                   } else {
4617                                        // just to make sure we don't exit before got something
4618                                        numberWanted++;
4619                                   }
4620                              }
4621                         }
4622                         break;
4623                    case ClpSimplex::atLowerBound:
4624                         value = cost[iSequence];
4625                         for (j = startColumn[iSequence];
4626                                   j < startColumn[iSequence] + length[iSequence]; j++) {
4627                              int jRow = row[j];
4628                              value -= duals[jRow] * element[j];
4629                         }
4630                         value = -value;
4631                         if (value > tolerance) {
4632                              numberWanted--;
4633                              if (value > bestDj) {
4634                                   // check flagged variable and correct dj
4635                                   if (!model->flagged(iSequence)) {
4636                                        bestDj = value;
4637                                        bestSequence = iSequence;
4638                                   } else {
4639                                        // just to make sure we don't exit before got something
4640                                        numberWanted++;
4641                                   }
4642                              }
4643                         }
4644                         break;
4645                    }
4646               }
4647               if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4648                    // give up
4649                    break;
4650               }
4651               if (!numberWanted)
4652                    break;
4653          }
4654          if (bestSequence != saveSequence) {
4655               // recompute dj
4656               double value = cost[bestSequence];
4657               for (j = startColumn[bestSequence];
4658                         j < startColumn[bestSequence] + length[bestSequence]; j++) {
4659                    int jRow = row[j];
4660                    value -= duals[jRow] * element[j];
4661               }
4662               reducedCost[bestSequence] = value;
4663               savedBestSequence_ = bestSequence;
4664               savedBestDj_ = reducedCost[savedBestSequence_];
4665          }
4666     }
4667     currentWanted_ = numberWanted;
4668}
4669// Sets up an effective RHS
4670void
4671ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
4672{
4673     delete [] rhsOffset_;
4674     int numberRows = model->numberRows();
4675     rhsOffset_ = new double[numberRows];
4676     rhsOffset(model, true);
4677}
4678// Gets rid of special copies
4679void
4680ClpPackedMatrix::clearCopies()
4681{
4682     delete rowCopy_;
4683     delete columnCopy_;
4684     rowCopy_ = NULL;
4685     columnCopy_ = NULL;
4686     flags_ &= ~(4 + 8);
4687     checkGaps();
4688#ifdef DO_CHECK_FLAGS
4689     checkFlags(0);
4690#endif
4691}
4692// makes sure active columns correct
4693int
4694ClpPackedMatrix::refresh(ClpSimplex * )
4695{
4696     numberActiveColumns_ = matrix_->getNumCols();
4697#if 0
4698     ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
4699     ClpPackedMatrix* rowCopy =
4700          dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4701     // Make sure it is really a ClpPackedMatrix
4702     assert (rowCopy != NULL);
4703
4704     const int * column = rowCopy->matrix_->getIndices();
4705     const CoinBigIndex * rowStart = rowCopy->matrix_->getVectorStarts();
4706     const int * rowLength = rowCopy->matrix_->getVectorLengths();
4707     const double * element = rowCopy->matrix_->getElements();
4708     int numberRows = rowCopy->matrix_->getNumRows();
4709     for (int i = 0; i < numberRows; i++) {
4710          if (!rowLength[i])
4711               printf("zero row %d\n", i);
4712     }
4713     delete rowCopy;
4714#endif
4715     checkGaps();
4716#ifdef DO_CHECK_FLAGS
4717     checkFlags(0);
4718#endif
4719     return 0;
4720}
4721
4722/* Scales rowCopy if column copy scaled
4723   Only called if scales already exist */
4724void
4725ClpPackedMatrix::scaleRowCopy(ClpModel * model) const
4726{
4727     if (model->rowCopy()) {
4728          // need to replace row by row
4729          int numberRows = model->numberRows();
4730#ifndef NDEBUG
4731          int numberColumns = matrix_->getNumCols();
4732#endif
4733          ClpMatrixBase * rowCopyBase = model->rowCopy();
4734#ifndef NDEBUG
4735          ClpPackedMatrix* rowCopy =
4736               dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4737          // Make sure it is really a ClpPackedMatrix
4738          assert (rowCopy != NULL);
4739#else
4740          ClpPackedMatrix* rowCopy =
4741               static_cast< ClpPackedMatrix*>(rowCopyBase);
4742#endif
4743
4744          const int * column = rowCopy->getIndices();
4745          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4746          double * element = rowCopy->getMutableElements();
4747          const double * rowScale = model->rowScale();
4748          const double * columnScale = model->columnScale();
4749          // scale row copy
4750          for (int iRow = 0; iRow < numberRows; iRow++) {
4751               CoinBigIndex j;
4752               double scale = rowScale[iRow];
4753               double * elementsInThisRow = element + rowStart[iRow];
4754               const int * columnsInThisRow = column + rowStart[iRow];
4755               int number = rowStart[iRow+1] - rowStart[iRow];
4756               assert (number <= numberColumns);
4757               for (j = 0; j < number; j++) {
4758                    int iColumn = columnsInThisRow[j];
4759                    elementsInThisRow[j] *= scale * columnScale[iColumn];
4760               }
4761          }
4762     }
4763}
4764/* Realy really scales column copy
4765   Only called if scales already exist.
4766   Up to user ro delete */
4767ClpMatrixBase *
4768ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const
4769{
4770     // need to replace column by column
4771#ifndef NDEBUG
4772     int numberRows = model->numberRows();
4773#endif
4774     int numberColumns = matrix_->getNumCols();
4775     ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
4776     const int * row = copy->getIndices();
4777     const CoinBigIndex * columnStart = copy->getVectorStarts();
4778     const int * length = copy->getVectorLengths();
4779     double * element = copy->getMutableElements();
4780     const double * rowScale = model->rowScale();
4781     const double * columnScale = model->columnScale();
4782     // scale column copy
4783     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4784          CoinBigIndex j;
4785          double scale = columnScale[iColumn];
4786          double * elementsInThisColumn = element + columnStart[iColumn];
4787          const int * rowsInThisColumn = row + columnStart[iColumn];
4788          int number = length[iColumn];
4789          assert (number <= numberRows);
4790          for (j = 0; j < number; j++) {
4791               int iRow = rowsInThisColumn[j];
4792               elementsInThisColumn[j] *= scale * rowScale[iRow];
4793          }
4794     }
4795     return copy;
4796}
4797// Really scale matrix
4798void
4799ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
4800{
4801     clearCopies();
4802     int numberColumns = matrix_->getNumCols();
4803     const int * row = matrix_->getIndices();
4804     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4805     const int * length = matrix_->getVectorLengths();
4806     double * element = matrix_->getMutableElements();
4807     // scale
4808     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4809          CoinBigIndex j;
4810          double scale = columnScale[iColumn];
4811          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length[iColumn]; j++) {
4812               int iRow = row[j];
4813               element[j] *= scale * rowScale[iRow];
4814          }
4815     }
4816}
4817/* Delete the columns whose indices are listed in <code>indDel</code>. */
4818void
4819ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
4820{
4821     if (matrix_->getNumCols())
4822          matrix_->deleteCols(numDel, indDel);
4823     clearCopies();
4824     numberActiveColumns_ = matrix_->getNumCols();
4825     // may now have gaps
4826     checkGaps();
4827#ifdef DO_CHECK_FLAGS
4828     checkFlags(0);
4829#endif
4830     matrix_->setExtraGap(0.0);
4831}
4832/* Delete the rows whose indices are listed in <code>indDel</code>. */
4833void
4834ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
4835{
4836     if (matrix_->getNumRows())
4837          matrix_->deleteRows(numDel, indDel);
4838     clearCopies();
4839     numberActiveColumns_ = matrix_->getNumCols();
4840     // may now have gaps
4841     checkGaps();
4842#ifdef DO_CHECK_FLAGS
4843     checkFlags(0);
4844#endif
4845     matrix_->setExtraGap(0.0);
4846}
4847#ifndef CLP_NO_VECTOR
4848// Append Columns
4849void
4850ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
4851{
4852     matrix_->appendCols(number, columns);
4853     numberActiveColumns_ = matrix_->getNumCols();
4854     clearCopies();
4855}
4856// Append Rows
4857void
4858ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
4859{
4860     matrix_->appendRows(number, rows);
4861     numberActiveColumns_ = matrix_->getNumCols();
4862     // may now have gaps
4863     checkGaps();
4864#ifdef DO_CHECK_FLAGS
4865     checkFlags(0);
4866#endif
4867     clearCopies();
4868}
4869#endif
4870/* Set the dimensions of the matrix. In effect, append new empty
4871   columns/rows to the matrix. A negative number for either dimension
4872   means that that dimension doesn't change. Otherwise the new dimensions
4873   MUST be at least as large as the current ones otherwise an exception
4874   is thrown. */
4875void
4876ClpPackedMatrix::setDimensions(int numrows, int numcols)
4877{
4878     matrix_->setDimensions(numrows, numcols);
4879#ifdef DO_CHECK_FLAGS
4880     checkFlags(0);
4881#endif
4882}
4883/* Append a set of rows/columns to the end of the matrix. Returns number of errors
4884   i.e. if any of the new rows/columns contain an index that's larger than the
4885   number of columns-1/rows-1 (if numberOther>0) or duplicates
4886   If 0 then rows, 1 if columns */
4887int
4888ClpPackedMatrix::appendMatrix(int number, int type,
4889                              const CoinBigIndex * starts, const int * index,
4890                              const double * element, int numberOther)
4891{
4892     int numberErrors = 0;
4893     // make sure other dimension is big enough
4894     if (type == 0) {
4895          // rows
4896          if (matrix_->isColOrdered() && numberOther > matrix_->getNumCols())
4897               matrix_->setDimensions(-1, numberOther);
4898          if (!matrix_->isColOrdered() || numberOther >= 0 || matrix_->getExtraGap()) {
4899               numberErrors = matrix_->appendRows(number, starts, index, element, numberOther);
4900          } else {
4901               //CoinPackedMatrix mm(*matrix_);
4902               matrix_->appendMinorFast(number, starts, index, element);
4903               //mm.appendRows(number,starts,index,element,numberOther);
4904               //if (!mm.isEquivalent(*matrix_)) {
4905               //printf("bad append\n");
4906               //abort();
4907               //}
4908          }
4909     } else {
4910          // columns
4911          if (!matrix_->isColOrdered() && numberOther > matrix_->getNumRows())
4912               matrix_->setDimensions(numberOther, -1);
4913          if (element)
4914            numberErrors = matrix_->appendCols(number, starts, index, element, numberOther);
4915          else
4916            matrix_->setDimensions(-1,matrix_->getNumCols()+number); // resize
4917     }
4918     clearCopies();
4919     numberActiveColumns_ = matrix_->getNumCols();
4920     return numberErrors;
4921}
4922void
4923ClpPackedMatrix::specialRowCopy(ClpSimplex * model, const ClpMatrixBase * rowCopy)
4924{
4925     delete rowCopy_;
4926     rowCopy_ = new ClpPackedMatrix2(model, rowCopy->getPackedMatrix());
4927     // See if anything in it
4928     if (!rowCopy_->usefulInfo()) {
4929          delete rowCopy_;
4930          rowCopy_ = NULL;
4931          flags_ &= ~4;
4932     } else {
4933          flags_ |= 4;
4934     }
4935}
4936void
4937ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
4938{
4939     delete columnCopy_;
4940     if ((flags_ & 16) != 0) {
4941          columnCopy_ = new ClpPackedMatrix3(model, matrix_);
4942          flags_ |= 8;
4943     } else {
4944          columnCopy_ = NULL;
4945     }
4946}
4947// Say we don't want special column copy
4948void
4949ClpPackedMatrix::releaseSpecialColumnCopy()
4950{
4951     flags_ &= ~(8 + 16);
4952     delete columnCopy_;
4953     columnCopy_ = NULL;
4954}
4955// Correct sequence in and out to give true value
4956void
4957ClpPackedMatrix::correctSequence(const ClpSimplex * model, int & sequenceIn, int & sequenceOut)
4958{
4959     if (columnCopy_) {
4960          if (sequenceIn != -999) {
4961               if (sequenceIn != sequenceOut) {
4962                    if (sequenceIn < numberActiveColumns_)
4963                         columnCopy_->swapOne(model, this, sequenceIn);
4964                    if (sequenceOut < numberActiveColumns_)
4965                         columnCopy_->swapOne(model, this, sequenceOut);
4966               }
4967          } else {
4968               // do all
4969               columnCopy_->sortBlocks(model);
4970          }
4971     }
4972}
4973// Check validity
4974void
4975ClpPackedMatrix::checkFlags(int type) const
4976{
4977     int iColumn;
4978     // get matrix data pointers
4979     //const int * row = matrix_->getIndices();
4980     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4981     const int * columnLength = matrix_->getVectorLengths();
4982     const double * elementByColumn = matrix_->getElements();
4983     if (!zeros()) {
4984          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4985               CoinBigIndex j;
4986               for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4987                    if (!elementByColumn[j])
4988                         abort();
4989               }
4990          }
4991     }
4992     if ((flags_ & 2) == 0) {
4993          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
4994               if (columnStart[iColumn+1] != columnStart[iColumn] + columnLength[iColumn]) {
4995                    abort();
4996               }
4997          }
4998     }
4999     if (type) {
5000          if ((flags_ & 2) != 0) {
5001               bool ok = true;
5002               for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
5003                    if (columnStart[iColumn+1] != columnStart[iColumn] + columnLength[iColumn]) {
5004                         ok = false;
5005                         break;
5006                    }
5007               }
5008               if (ok)
5009                 COIN_DETAIL_PRINT(printf("flags_ could be 0\n"));
5010          }
5011     }
5012}
5013//#############################################################################
5014// Constructors / Destructor / Assignment
5015//#############################################################################
5016
5017//-------------------------------------------------------------------
5018// Default Constructor
5019//-------------------------------------------------------------------
5020ClpPackedMatrix2::ClpPackedMatrix2 ()
5021     : numberBlocks_(0),
5022       numberRows_(0),
5023       offset_(NULL),
5024       count_(NULL),
5025       rowStart_(NULL),
5026       column_(NULL),
5027       work_(NULL)
5028{
5029#ifdef THREAD
5030     threadId_ = NULL;
5031     info_ = NULL;
5032#endif
5033}
5034//-------------------------------------------------------------------
5035// Useful Constructor
5036//-------------------------------------------------------------------
5037ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * , const CoinPackedMatrix * rowCopy)
5038     : numberBlocks_(0),
5039       numberRows_(0),
5040       offset_(NULL),
5041       count_(NULL),
5042       rowStart_(NULL),
5043       column_(NULL),
5044       work_(NULL)
5045{
5046#ifdef THREAD
5047     threadId_ = NULL;
5048     info_ = NULL;
5049#endif
5050     numberRows_ = rowCopy->getNumRows();
5051     if (!numberRows_)
5052          return;
5053     int numberColumns = rowCopy->getNumCols();
5054     const int * column = rowCopy->getIndices();
5055     const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
5056     const int * length = rowCopy->getVectorLengths();
5057     const double * element = rowCopy->getElements();
5058     int chunk = 32768; // tune
5059     //chunk=100;
5060     // tune
5061#if 0
5062     int chunkY[7] = {1024, 2048, 4096, 8192, 16384, 32768, 65535};
5063     int its = model->maximumIterations();
5064     if (its >= 1000000 && its < 1000999) {
5065          its -= 1000000;
5066          its = its / 10;
5067          if (its >= 7) abort();
5068          chunk = chunkY[its];
5069          printf("chunk size %d\n", chunk);
5070#define cpuid(func,ax,bx,cx,dx)\
5071    __asm__ __volatile__ ("cpuid":\
5072    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
5073          unsigned int a, b, c, d;
5074          int func = 0;
5075          cpuid(func, a, b, c, d);
5076          {
5077               int i;
5078               unsigned int value;
5079               value = b;
5080               for (i = 0; i < 4; i++) {
5081                    printf("%c", (value & 0xff));
5082                    value = value >> 8;
5083               }
5084               value = d;
5085               for (i = 0; i < 4; i++) {
5086                    printf("%c", (value & 0xff));
5087                    value = value >> 8;
5088               }
5089               value = c;
5090               for (i = 0; i < 4; i++) {
5091                    printf("%c", (value & 0xff));
5092                    value = value >> 8;
5093               }
5094               printf("\n");
5095               int maxfunc = a;
5096               if (maxfunc > 10) {
5097                    printf("not intel?\n");
5098                    abort();
5099               }
5100               for (func = 1; func <= maxfunc; func++) {
5101                    cpuid(func, a, b, c, d);
5102                    printf("func %d, %x %x %x %x\n", func, a, b, c, d);
5103               }
5104          }
5105#else
5106     if (numberColumns > 10000 || chunk == 100) {
5107#endif
5108     } else {
5109          //printf("no chunk\n");
5110          return;
5111     }
5112     // Could also analyze matrix to get natural breaks
5113     numberBlocks_ = (numberColumns + chunk - 1) / chunk;
5114#ifdef THREAD
5115     // Get work areas
5116     threadId_ = new pthread_t [numberBlocks_];
5117     info_ = new dualColumn0Struct[numberBlocks_];
5118#endif
5119     // Even out
5120     chunk = (numberColumns + numberBlocks_ - 1) / numberBlocks_;
5121     offset_ = new int[numberBlocks_+1];
5122     offset_[numberBlocks_] = numberColumns;
5123     int nRow = numberBlocks_ * numberRows_;
5124     count_ = new unsigned short[nRow];
5125     memset(count_, 0, nRow * sizeof(unsigned short));
5126     rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
5127     CoinBigIndex nElement = rowStart[numberRows_];
5128     rowStart_[nRow+numberRows_] = nElement;
5129     column_ = new unsigned short [nElement];
5130     // assumes int <= double
5131     int sizeWork = 6 * numberBlocks_;
5132     work_ = new double[sizeWork];;
5133     int iBlock;
5134     int nZero = 0;
5135     for (iBlock = 0; iBlock < numberBlocks_; iBlock++) {
5136          int start = iBlock * chunk;
5137          offset_[iBlock] = start;
5138          int end = start + chunk;
5139          for (int iRow = 0; iRow < numberRows_; iRow++) {
5140               if (rowStart[iRow+1] != rowStart[iRow] + length[iRow]) {
5141                    printf("not packed correctly - gaps\n");
5142                    abort();
5143               }
5144               bool lastFound = false;
5145               int nFound = 0;
5146               for (CoinBigIndex j = rowStart[iRow];
5147                         j < rowStart[iRow] + length[iRow]; j++) {
5148                    int iColumn = column[j];
5149                    if (iColumn >= start) {
5150                         if (iColumn < end) {
5151                              if (!element[j]) {
5152                                   printf("not packed correctly - zero element\n");
5153                                   abort();
5154                              }
5155                              column_[j] = static_cast<unsigned short>(iColumn - start);
5156                              nFound++;
5157                              if (lastFound) {
5158                                   printf("not packed correctly - out of order\n");
5159                                   abort();
5160                              }
5161                         } else {
5162                              //can't find any more
5163                              lastFound = true;
5164                         }
5165                    }
5166               }
5167               count_[iRow*numberBlocks_+iBlock] = static_cast<unsigned short>(nFound);
5168               if (!nFound)
5169                    nZero++;
5170          }
5171     }
5172     //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
5173     //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
5174}
5175
5176//-------------------------------------------------------------------
5177// Copy constructor
5178//-------------------------------------------------------------------
5179ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs)
5180     : numberBlocks_(rhs.numberBlocks_),
5181       numberRows_(rhs.numberRows_)
5182{
5183     if (numberBlocks_) {
5184          offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + 1);
5185          int nRow = numberBlocks_ * numberRows_;
5186          count_ = CoinCopyOfArray(rhs.count_, nRow);
5187          rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + 1);
5188          CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5189          column_ = CoinCopyOfArray(rhs.column_, nElement);
5190          int sizeWork = 6 * numberBlocks_;
5191          work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5192#ifdef THREAD
5193          threadId_ = new pthread_t [numberBlocks_];
5194          info_ = new dualColumn0Struct[numberBlocks_];
5195#endif
5196     } else {
5197          offset_ = NULL;
5198          count_ = NULL;
5199          rowStart_ = NULL;
5200          column_ = NULL;
5201          work_ = NULL;
5202#ifdef THREAD
5203          threadId_ = NULL;
5204          info_ = NULL;
5205#endif
5206     }
5207}
5208//-------------------------------------------------------------------
5209// Destructor
5210//-------------------------------------------------------------------
5211ClpPackedMatrix2::~ClpPackedMatrix2 ()
5212{
5213     delete [] offset_;
5214     delete [] count_;
5215     delete [] rowStart_;
5216     delete [] column_;
5217     delete [] work_;
5218#ifdef THREAD
5219     delete [] threadId_;
5220     delete [] info_;
5221#endif
5222}
5223
5224//----------------------------------------------------------------
5225// Assignment operator
5226//-------------------------------------------------------------------
5227ClpPackedMatrix2 &
5228ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
5229{
5230     if (this != &rhs) {
5231          numberBlocks_ = rhs.numberBlocks_;
5232          numberRows_ = rhs.numberRows_;
5233          delete [] offset_;
5234          delete [] count_;
5235          delete [] rowStart_;
5236          delete [] column_;
5237          delete [] work_;
5238#ifdef THREAD
5239          delete [] threadId_;
5240          delete [] info_;
5241#endif
5242          if (numberBlocks_) {
5243               offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + 1);
5244               int nRow = numberBlocks_ * numberRows_;
5245               count_ = CoinCopyOfArray(rhs.count_, nRow);
5246               rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + 1);
5247               CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5248               column_ = CoinCopyOfArray(rhs.column_, nElement);
5249               int sizeWork = 6 * numberBlocks_;
5250               work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5251#ifdef THREAD
5252               threadId_ = new pthread_t [numberBlocks_];
5253               info_ = new dualColumn0Struct[numberBlocks_];
5254#endif
5255          } else {
5256               offset_ = NULL;
5257               count_ = NULL;
5258               rowStart_ = NULL;
5259               column_ = NULL;
5260               work_ = NULL;
5261#ifdef THREAD
5262               threadId_ = NULL;
5263               info_ = NULL;
5264#endif
5265          }
5266     }
5267     return *this;
5268}
5269static int dualColumn0(const ClpSimplex * model, double * spare,
5270                       int * spareIndex, const double * arrayTemp,
5271                       const int * indexTemp, int numberIn,
5272                       int offset, double acceptablePivot, double * bestPossiblePtr,
5273                       double * upperThetaPtr, int * posFreePtr, double * freePivotPtr)
5274{
5275     // do dualColumn0
5276     int i;
5277     int numberRemaining = 0;
5278     double bestPossible = 0.0;
5279     double upperTheta = 1.0e31;
5280     double freePivot = acceptablePivot;
5281     int posFree = -1;
5282     const double * reducedCost = model->djRegion(1);
5283     double dualTolerance = model->dualTolerance();
5284     // We can also see if infeasible or pivoting on free
5285     double tentativeTheta = 1.0e25;
5286     for (i = 0; i < numberIn; i++) {
5287          double alpha = arrayTemp[i];
5288          int iSequence = indexTemp[i] + offset;
5289          double oldValue;
5290          double value;
5291          bool keep;
5292
5293          switch(model->getStatus(iSequence)) {
5294
5295          case ClpSimplex::basic:
5296          case ClpSimplex::isFixed:
5297               break;
5298          case ClpSimplex::isFree:
5299          case ClpSimplex::superBasic:
5300               bestPossible = CoinMax(bestPossible, fabs(alpha));
5301               oldValue = reducedCost[iSequence];
5302               // If free has to be very large - should come in via dualRow
5303               if (model->getStatus(iSequence) == ClpSimplex::isFree && fabs(alpha) < 1.0e-3)
5304                    break;
5305               if (oldValue > dualTolerance) {
5306                    keep = true;
5307               } else if (oldValue < -dualTolerance) {
5308                    keep = true;
5309               } else {
5310                    if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5))
5311                         keep = true;
5312                    else
5313                         keep = false;
5314               }
5315               if (keep) {
5316                    // free - choose largest
5317                    if (fabs(alpha) > freePivot) {
5318                         freePivot = fabs(alpha);
5319                         posFree = i;
5320                    }
5321               }
5322               break;
5323          case ClpSimplex::atUpperBound:
5324               oldValue = reducedCost[iSequence];
5325               value = oldValue - tentativeTheta * alpha;
5326               //assert (oldValue<=dualTolerance*1.0001);
5327               if (value > dualTolerance) {
5328                    bestPossible = CoinMax(bestPossible, -alpha);
5329                    value = oldValue - upperTheta * alpha;
5330                    if (value > dualTolerance && -alpha >= acceptablePivot)
5331                         upperTheta = (oldValue - dualTolerance) / alpha;
5332                    // add to list
5333                    spare[numberRemaining] = alpha;
5334                    spareIndex[numberRemaining++] = iSequence;
5335               }
5336               break;
5337          case ClpSimplex::atLowerBound:
5338               oldValue = reducedCost[iSequence];
5339               value = oldValue - tentativeTheta * alpha;
5340               //assert (oldValue>=-dualTolerance*1.0001);
5341               if (value < -dualTolerance) {
5342                    bestPossible = CoinMax(bestPossible, alpha);
5343                    value = oldValue - upperTheta * alpha;
5344                    if (value < -dualTolerance && alpha >= acceptablePivot)
5345                         upperTheta = (oldValue + dualTolerance) / alpha;
5346                    // add to list
5347                    spare[numberRemaining] = alpha;
5348                    spareIndex[numberRemaining++] = iSequence;
5349               }
5350               break;
5351          }
5352     }
5353     *bestPossiblePtr = bestPossible;
5354     *upperThetaPtr = upperTheta;
5355     *freePivotPtr = freePivot;
5356     *posFreePtr = posFree;
5357     return numberRemaining;
5358}
5359static int doOneBlock(double * array, int * index,
5360                      const double * pi, const CoinBigIndex * rowStart, const double * element,
5361                      const unsigned short * column, int numberInRowArray, int numberLook)
5362{
5363     int iWhich = 0;
5364     int nextN = 0;
5365     CoinBigIndex nextStart = 0;
5366     double nextPi = 0.0;
5367     for (; iWhich < numberInRowArray; iWhich++) {
5368          nextStart = rowStart[0];
5369          nextN = rowStart[numberInRowArray] - nextStart;
5370          rowStart++;
5371          if (nextN) {
5372               nextPi = pi[iWhich];
5373               break;
5374          }
5375     }
5376     int i;
5377     while (iW