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

Last change on this file since 2235 was 2235, checked in by forrest, 4 years ago

stuff for vector matrix

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