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

Last change on this file since 2271 was 2271, checked in by forrest, 2 years ago

change some ints to CoinBigIndex?

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