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

Last change on this file since 1836 was 1836, checked in by lou, 8 years ago

Fix bug of r1679 again, plus a few other changes to set ClpPackedMatrix::flags_
correctly.

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