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

Last change on this file since 2367 was 2367, checked in by forrest, 5 months ago

take out misleading comments

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