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

Last change on this file since 1002 was 1002, checked in by forrest, 13 years ago

bug - active columns

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 114.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5#include <cstdio>
6
7#include "CoinPragma.hpp"
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10
11#include "ClpSimplex.hpp"
12#include "ClpSimplexDual.hpp"
13#include "ClpFactorization.hpp"
14#ifndef SLIM_CLP
15#include "ClpQuadraticObjective.hpp"
16#endif
17//#define THREAD
18// at end to get min/max!
19#include "ClpPackedMatrix.hpp"
20#include "ClpMessage.hpp"
21//#define COIN_PREFETCH
22#ifdef COIN_PREFETCH
23#if 1
24#define coin_prefetch(mem) \
25         __asm__ __volatile__ ("prefetchnta %0" : : "m" (*((char *)(mem))))
26#else
27#define coin_prefetch(mem) \
28         __asm__ __volatile__ ("prefetch %0" : : "m" (*((char *)(mem))))
29#endif
30#else
31// dummy
32#define coin_prefetch(mem)
33#endif
34
35//#############################################################################
36// Constructors / Destructor / Assignment
37//#############################################################################
38
39//-------------------------------------------------------------------
40// Default Constructor
41//-------------------------------------------------------------------
42ClpPackedMatrix::ClpPackedMatrix () 
43  : ClpMatrixBase(),
44    matrix_(NULL),
45    numberActiveColumns_(0),
46    zeroElements_(false),
47    hasGaps_(true),
48    rowCopy_(NULL)
49{
50  setType(1);
51}
52
53//-------------------------------------------------------------------
54// Copy constructor
55//-------------------------------------------------------------------
56ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs) 
57: ClpMatrixBase(rhs)
58{ 
59  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
60  numberActiveColumns_ = rhs.numberActiveColumns_;
61  zeroElements_ = rhs.zeroElements_;
62  hasGaps_ = rhs.hasGaps_;
63  int numberRows = getNumRows();
64  if (rhs.rhsOffset_&&numberRows) {
65    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
66  } else {
67    rhsOffset_=NULL;
68  }
69  if (rhs.rowCopy_) {
70    rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
71  } else {
72    rowCopy_ = NULL;
73  }
74}
75
76//-------------------------------------------------------------------
77// assign matrix (for space reasons)
78//-------------------------------------------------------------------
79ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs) 
80: ClpMatrixBase()
81{ 
82  matrix_ = rhs;
83  zeroElements_ = false;
84  hasGaps_ = true;
85  numberActiveColumns_ = matrix_->getNumCols();
86  rowCopy_ = NULL;
87  setType(1);
88 
89}
90
91ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs) 
92: ClpMatrixBase()
93{ 
94  matrix_ = new CoinPackedMatrix(rhs);
95  numberActiveColumns_ = matrix_->getNumCols();
96  rowCopy_ = NULL;
97  zeroElements_ = false;
98  hasGaps_ = true;
99  setType(1);
100 
101}
102
103//-------------------------------------------------------------------
104// Destructor
105//-------------------------------------------------------------------
106ClpPackedMatrix::~ClpPackedMatrix ()
107{
108  delete matrix_;
109  delete rowCopy_;
110}
111
112//----------------------------------------------------------------
113// Assignment operator
114//-------------------------------------------------------------------
115ClpPackedMatrix &
116ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
117{
118  if (this != &rhs) {
119    ClpMatrixBase::operator=(rhs);
120    delete matrix_;
121    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
122    numberActiveColumns_ = rhs.numberActiveColumns_;
123    zeroElements_ = rhs.zeroElements_;
124    hasGaps_ = rhs.hasGaps_;
125    delete rowCopy_;
126    if (rhs.rowCopy_) {
127      rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
128    } else {
129      rowCopy_ = NULL;
130    }
131  }
132  return *this;
133}
134//-------------------------------------------------------------------
135// Clone
136//-------------------------------------------------------------------
137ClpMatrixBase * ClpPackedMatrix::clone() const
138{
139  return new ClpPackedMatrix(*this);
140}
141/* Subset clone (without gaps).  Duplicates are allowed
142   and order is as given */
143ClpMatrixBase * 
144ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
145                              int numberColumns, 
146                              const int * whichColumns) const 
147{
148  return new ClpPackedMatrix(*this, numberRows, whichRows,
149                                   numberColumns, whichColumns);
150}
151/* Subset constructor (without gaps).  Duplicates are allowed
152   and order is as given */
153ClpPackedMatrix::ClpPackedMatrix (
154                       const ClpPackedMatrix & rhs,
155                       int numberRows, const int * whichRows,
156                       int numberColumns, const int * whichColumns)
157: ClpMatrixBase(rhs)
158{
159  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
160                                 numberColumns,whichColumns);
161  numberActiveColumns_ = matrix_->getNumCols();
162  zeroElements_ = rhs.zeroElements_;
163  hasGaps_ = rhs.hasGaps_;
164  rowCopy_ = NULL;
165}
166ClpPackedMatrix::ClpPackedMatrix (
167                       const CoinPackedMatrix & rhs,
168                       int numberRows, const int * whichRows,
169                       int numberColumns, const int * whichColumns)
170: ClpMatrixBase()
171{
172  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
173                                 numberColumns,whichColumns);
174  numberActiveColumns_ = matrix_->getNumCols();
175  zeroElements_ = false;
176  hasGaps_=true;
177  rowCopy_ = NULL;
178  setType(1);
179}
180
181/* Returns a new matrix in reverse order without gaps */
182ClpMatrixBase * 
183ClpPackedMatrix::reverseOrderedCopy() const
184{
185  ClpPackedMatrix * copy = new ClpPackedMatrix();
186  copy->matrix_= new CoinPackedMatrix();
187  copy->matrix_->setExtraGap(0.0);
188  copy->matrix_->setExtraMajor(0.0);
189  copy->matrix_->reverseOrderedCopyOf(*matrix_);
190  //copy->matrix_->removeGaps();
191  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
192  copy->hasGaps_=false;
193  return copy;
194}
195//unscaled versions
196void 
197ClpPackedMatrix::times(double scalar,
198                   const double * x, double * y) const
199{
200  int iRow,iColumn;
201  // get matrix data pointers
202  const int * row = matrix_->getIndices();
203  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
204  const int * columnLength = matrix_->getVectorLengths(); 
205  const double * elementByColumn = matrix_->getElements();
206  //memset(y,0,matrix_->getNumRows()*sizeof(double));
207  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
208    CoinBigIndex j;
209    double value = scalar*x[iColumn];
210    if (value) {
211      for (j=columnStart[iColumn];
212           j<columnStart[iColumn]+columnLength[iColumn];j++) {
213        iRow=row[j];
214        y[iRow] += value*elementByColumn[j];
215      }
216    }
217  }
218}
219void 
220ClpPackedMatrix::transposeTimes(double scalar,
221                                const double * x, double * y) const
222{
223  int iColumn;
224  // get matrix data pointers
225  const int * row = matrix_->getIndices();
226  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
227  const int * columnLength = matrix_->getVectorLengths(); 
228  const double * elementByColumn = matrix_->getElements();
229  if (!hasGaps_) {
230    if (scalar==1.0) {
231      CoinBigIndex start=columnStart[0];
232      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
233        CoinBigIndex j;
234        CoinBigIndex next=columnStart[iColumn+1];
235        double value=y[iColumn];
236        for (j=start;j<next;j++) {
237          int jRow=row[j];
238          value += x[jRow]*elementByColumn[j];
239        }
240        start=next;
241        y[iColumn] = value;
242      }
243    } else if (scalar==-1.0) {
244      CoinBigIndex start=columnStart[0];
245      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
246        CoinBigIndex j;
247        CoinBigIndex next=columnStart[iColumn+1];
248        double value=y[iColumn];
249        for (j=start;j<next;j++) {
250          int jRow=row[j];
251          value -= x[jRow]*elementByColumn[j];
252        }
253        start=next;
254        y[iColumn] = value;
255      }
256    } else {
257      CoinBigIndex start=columnStart[0];
258      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
259        CoinBigIndex j;
260        CoinBigIndex next=columnStart[iColumn+1];
261        double value=0.0;
262        for (j=start;j<next;j++) {
263          int jRow=row[j];
264          value += x[jRow]*elementByColumn[j];
265        }
266        start=next;
267        y[iColumn] += value*scalar;
268      }
269    }
270  } else {
271    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
272      CoinBigIndex j;
273      double value=0.0;
274      for (j=columnStart[iColumn];
275           j<columnStart[iColumn]+columnLength[iColumn];j++) {
276        int jRow=row[j];
277        value += x[jRow]*elementByColumn[j];
278      }
279      y[iColumn] += value*scalar;
280    }
281  }
282}
283void 
284ClpPackedMatrix::times(double scalar,
285                       const double * x, double * y,
286                       const double * rowScale, 
287                       const double * columnScale) const
288{
289  if (rowScale) {
290    int iRow,iColumn;
291    // get matrix data pointers
292    const int * row = matrix_->getIndices();
293    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
294    const int * columnLength = matrix_->getVectorLengths(); 
295    const double * elementByColumn = matrix_->getElements();
296    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
297      CoinBigIndex j;
298      double value = x[iColumn];
299      if (value) {
300        // scaled
301        value *= scalar*columnScale[iColumn];
302        for (j=columnStart[iColumn];
303             j<columnStart[iColumn]+columnLength[iColumn];j++) {
304          iRow=row[j];
305          y[iRow] += value*elementByColumn[j]*rowScale[iRow];
306        }
307      }
308    }
309  } else {
310    times(scalar,x,y);
311  }
312}
313void 
314ClpPackedMatrix::transposeTimes( double scalar,
315                                 const double * x, double * y,
316                                 const double * rowScale, 
317                                 const double * columnScale,
318                                 double * spare) const
319{
320  if (rowScale) {
321    int iColumn;
322    // get matrix data pointers
323    const int * row = matrix_->getIndices();
324    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
325    const int * columnLength = matrix_->getVectorLengths(); 
326    const double * elementByColumn = matrix_->getElements();
327    if (!spare) {
328      if (!hasGaps_) {
329        CoinBigIndex start=columnStart[0];
330        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
331          CoinBigIndex j;
332          CoinBigIndex next=columnStart[iColumn+1];
333          double value=0.0;
334          // scaled
335          for (j=start;j<next;j++) {
336            int jRow=row[j];
337            value += x[jRow]*elementByColumn[j]*rowScale[jRow];
338          }
339          start=next;
340          y[iColumn] += value*scalar*columnScale[iColumn];
341        }
342      } else {
343        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
344          CoinBigIndex j;
345          double value=0.0;
346          // scaled
347          for (j=columnStart[iColumn];
348               j<columnStart[iColumn]+columnLength[iColumn];j++) {
349            int jRow=row[j];
350            value += x[jRow]*elementByColumn[j]*rowScale[jRow];
351          }
352          y[iColumn] += value*scalar*columnScale[iColumn];
353        }
354      }
355    } else {
356      // can use spare region
357      int iRow;
358      int numberRows = getNumRows();
359      for (iRow=0;iRow<numberRows;iRow++)
360        spare[iRow] = x[iRow]*rowScale[iRow];
361      if (!hasGaps_) {
362        CoinBigIndex start=columnStart[0];
363        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
364          CoinBigIndex j;
365          CoinBigIndex next=columnStart[iColumn+1];
366          double value=0.0;
367          // scaled
368          for (j=start;j<next;j++) {
369            int jRow=row[j];
370            value += spare[jRow]*elementByColumn[j];
371          }
372          start=next;
373          y[iColumn] += value*scalar*columnScale[iColumn];
374        }
375      } else {
376        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
377          CoinBigIndex j;
378          double value=0.0;
379          // scaled
380          for (j=columnStart[iColumn];
381               j<columnStart[iColumn]+columnLength[iColumn];j++) {
382            int jRow=row[j];
383            value += spare[jRow]*elementByColumn[j];
384          }
385          y[iColumn] += value*scalar*columnScale[iColumn];
386        }
387      }
388      // no need to zero out
389      //for (iRow=0;iRow<numberRows;iRow++)
390      //spare[iRow] = 0.0;
391    }
392  } else {
393    transposeTimes(scalar,x,y);
394  }
395}
396/* Return <code>x * A + y</code> in <code>z</code>.
397        Squashes small elements and knows about ClpSimplex */
398void 
399ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
400                              const CoinIndexedVector * rowArray,
401                              CoinIndexedVector * y,
402                              CoinIndexedVector * columnArray) const
403{
404  columnArray->clear();
405  double * pi = rowArray->denseVector();
406  int numberNonZero=0;
407  int * index = columnArray->getIndices();
408  double * array = columnArray->denseVector();
409  int numberInRowArray = rowArray->getNumElements();
410  // maybe I need one in OsiSimplex
411  double zeroTolerance = model->factorization()->zeroTolerance();
412  int numberRows = model->numberRows();
413#ifndef NO_RTTI
414  ClpPackedMatrix* rowCopy =
415    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
416#else
417  ClpPackedMatrix* rowCopy =
418    static_cast< ClpPackedMatrix*>(model->rowCopy());
419#endif
420  bool packed = rowArray->packedMode();
421  double factor = 0.27;
422  // We may not want to do by row if there may be cache problems
423  // It would be nice to find L2 cache size - for moment 512K
424  // Be slightly optimistic
425  if (numberActiveColumns_*sizeof(double)>1000000) {
426    if (numberRows*10<numberActiveColumns_)
427      factor *= 0.333333333;
428    else if (numberRows*4<numberActiveColumns_)
429      factor *= 0.5;
430    else if (numberRows*2<numberActiveColumns_)
431      factor *= 0.66666666667;
432    //if (model->numberIterations()%50==0)
433    //printf("%d nonzero\n",numberInRowArray);
434  }
435  // if not packed then bias a bit more towards by column
436  if (!packed)
437    factor *= 0.9;
438  assert (!y->getNumElements());
439  double multiplierX=0.8;
440  double factor2 = factor*multiplierX;
441  if (packed&&rowCopy_&&numberInRowArray>2&&numberInRowArray>factor2*numberRows&&
442      numberInRowArray<0.9*numberRows) {
443    rowCopy_->transposeTimes(model,rowCopy->getPackedMatrix(),rowArray,y,columnArray);
444    return;
445  }
446  if (numberInRowArray>factor*numberRows||!rowCopy) {
447    // do by column
448    // If no gaps - can do a bit faster
449    if (!hasGaps_) {
450      transposeTimesByColumn( model,  scalar,
451                              rowArray, y, columnArray);
452      return;
453    }
454    int iColumn;
455    // get matrix data pointers
456    const int * row = matrix_->getIndices();
457    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
458    const int * columnLength = matrix_->getVectorLengths(); 
459    const double * elementByColumn = matrix_->getElements();
460    const double * rowScale = model->rowScale();
461    if (packed) {
462      // need to expand pi into y
463      assert(y->capacity()>=numberRows);
464      double * piOld = pi;
465      pi = y->denseVector();
466      const int * whichRow = rowArray->getIndices();
467      int i;
468      if (!rowScale) {
469        // modify pi so can collapse to one loop
470        if (scalar==-1.0) {
471          for (i=0;i<numberInRowArray;i++) {
472            int iRow = whichRow[i];
473            pi[iRow]=-piOld[i];
474          }
475        } else {
476          for (i=0;i<numberInRowArray;i++) {
477            int iRow = whichRow[i];
478            pi[iRow]=scalar*piOld[i];
479          }
480        }
481        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
482          double value = 0.0;
483          CoinBigIndex j;
484          for (j=columnStart[iColumn];
485               j<columnStart[iColumn]+columnLength[iColumn];j++) {
486            int iRow = row[j];
487            value += pi[iRow]*elementByColumn[j];
488          }
489          if (fabs(value)>zeroTolerance) {
490            array[numberNonZero]=value;
491            index[numberNonZero++]=iColumn;
492          }
493        }
494      } else {
495        // scaled
496        // modify pi so can collapse to one loop
497        if (scalar==-1.0) {
498          for (i=0;i<numberInRowArray;i++) {
499            int iRow = whichRow[i];
500            pi[iRow]=-piOld[i]*rowScale[iRow];
501          }
502        } else {
503          for (i=0;i<numberInRowArray;i++) {
504            int iRow = whichRow[i];
505            pi[iRow]=scalar*piOld[i]*rowScale[iRow];
506          }
507        }
508        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
509          double value = 0.0;
510          CoinBigIndex j;
511          const double * columnScale = model->columnScale();
512          for (j=columnStart[iColumn];
513               j<columnStart[iColumn]+columnLength[iColumn];j++) {
514            int iRow = row[j];
515            value += pi[iRow]*elementByColumn[j];
516          }
517          value *= columnScale[iColumn];
518          if (fabs(value)>zeroTolerance) {
519            array[numberNonZero]=value;
520            index[numberNonZero++]=iColumn;
521          }
522        }
523      }
524      // zero out
525      for (i=0;i<numberInRowArray;i++) {
526        int iRow = whichRow[i];
527        pi[iRow]=0.0;
528      }
529    } else {
530      if (!rowScale) {
531        if (scalar==-1.0) {
532          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
533            double value = 0.0;
534            CoinBigIndex j;
535            for (j=columnStart[iColumn];
536                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
537              int iRow = row[j];
538              value += pi[iRow]*elementByColumn[j];
539            }
540            if (fabs(value)>zeroTolerance) {
541              index[numberNonZero++]=iColumn;
542              array[iColumn]=-value;
543            }
544          }
545        } else {
546          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
547            double value = 0.0;
548            CoinBigIndex j;
549            for (j=columnStart[iColumn];
550                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
551              int iRow = row[j];
552              value += pi[iRow]*elementByColumn[j];
553            }
554            value *= scalar;
555            if (fabs(value)>zeroTolerance) {
556              index[numberNonZero++]=iColumn;
557              array[iColumn]=value;
558            }
559          }
560        }
561      } else {
562        // scaled
563        if (scalar==-1.0) {
564          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
565            double value = 0.0;
566            CoinBigIndex j;
567            const double * columnScale = model->columnScale();
568            for (j=columnStart[iColumn];
569                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
570              int iRow = row[j];
571              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
572            }
573            value *= columnScale[iColumn];
574            if (fabs(value)>zeroTolerance) {
575              index[numberNonZero++]=iColumn;
576              array[iColumn]=-value;
577            }
578          }
579        } else {
580          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
581            double value = 0.0;
582            CoinBigIndex j;
583            const double * columnScale = model->columnScale();
584            for (j=columnStart[iColumn];
585                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
586              int iRow = row[j];
587              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
588            }
589            value *= scalar*columnScale[iColumn];
590            if (fabs(value)>zeroTolerance) {
591              index[numberNonZero++]=iColumn;
592              array[iColumn]=value;
593            }
594          }
595        }
596      }
597    }
598    columnArray->setNumElements(numberNonZero);
599    y->setNumElements(0);
600  } else {
601    // do by row
602    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
603  }
604  if (packed)
605    columnArray->setPackedMode(true);
606  if (0) {
607    columnArray->checkClean();
608    int numberNonZero=columnArray->getNumElements();;
609    int * index = columnArray->getIndices();
610    double * array = columnArray->denseVector();
611    int i;
612    for (i=0;i<numberNonZero;i++) {
613      int j=index[i];
614      double value;
615      if (packed)
616        value=array[i];
617      else
618        value=array[j];
619      printf("Ti %d %d %g\n",i,j,value);
620    }
621  }
622}
623/* Return <code>x * scalar * A + y</code> in <code>z</code>.
624   Note - If x packed mode - then z packed mode
625   This does by column and knows no gaps
626   Squashes small elements and knows about ClpSimplex */
627void 
628ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
629                                        const CoinIndexedVector * rowArray,
630                                        CoinIndexedVector * y,
631                                        CoinIndexedVector * columnArray) const
632{
633  double * pi = rowArray->denseVector();
634  int numberNonZero=0;
635  int * index = columnArray->getIndices();
636  double * array = columnArray->denseVector();
637  int numberInRowArray = rowArray->getNumElements();
638  // maybe I need one in OsiSimplex
639  double zeroTolerance = model->factorization()->zeroTolerance();
640  bool packed = rowArray->packedMode();
641  // do by column
642  int iColumn;
643  // get matrix data pointers
644  const int * row = matrix_->getIndices();
645  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
646  const double * elementByColumn = matrix_->getElements();
647  const double * rowScale = model->rowScale();
648  assert (!y->getNumElements());
649  assert (numberActiveColumns_>0);
650  if (packed) {
651    // need to expand pi into y
652    assert(y->capacity()>=model->numberRows());
653    double * piOld = pi;
654    pi = y->denseVector();
655    const int * whichRow = rowArray->getIndices();
656    int i;
657    if (!rowScale) {
658      // modify pi so can collapse to one loop
659      if (scalar==-1.0) {
660        for (i=0;i<numberInRowArray;i++) {
661          int iRow = whichRow[i];
662          pi[iRow]=-piOld[i];
663        }
664      } else {
665        for (i=0;i<numberInRowArray;i++) {
666          int iRow = whichRow[i];
667          pi[iRow]=scalar*piOld[i];
668        }
669      }
670      double value = 0.0;
671      CoinBigIndex j;
672      CoinBigIndex end = columnStart[1];
673      for (j=columnStart[0]; j<end;j++) {
674        int iRow = row[j];
675        value += pi[iRow]*elementByColumn[j];
676      }
677      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
678        CoinBigIndex start = end;
679        end = columnStart[iColumn+2];
680        if (fabs(value)>zeroTolerance) {
681          array[numberNonZero]=value;
682          index[numberNonZero++]=iColumn;
683        }
684        value = 0.0;
685        for (j=start; j<end;j++) {
686          int iRow = row[j];
687          value += pi[iRow]*elementByColumn[j];
688        }
689      }
690      if (fabs(value)>zeroTolerance) {
691        array[numberNonZero]=value;
692        index[numberNonZero++]=iColumn;
693      }
694    } else {
695      // scaled
696      // modify pi so can collapse to one loop
697      if (scalar==-1.0) {
698        for (i=0;i<numberInRowArray;i++) {
699          int iRow = whichRow[i];
700          pi[iRow]=-piOld[i]*rowScale[iRow];
701        }
702      } else {
703        for (i=0;i<numberInRowArray;i++) {
704          int iRow = whichRow[i];
705          pi[iRow]=scalar*piOld[i]*rowScale[iRow];
706        }
707      }
708      const double * columnScale = model->columnScale();
709      double value = 0.0;
710      double scale=columnScale[0];
711      CoinBigIndex j;
712      CoinBigIndex end = columnStart[1];
713      for (j=columnStart[0]; j<end;j++) {
714        int iRow = row[j];
715        value += pi[iRow]*elementByColumn[j];
716      }
717      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
718        value *= scale;
719        CoinBigIndex start = end;
720        scale = columnScale[iColumn+1];
721        end = columnStart[iColumn+2];
722        if (fabs(value)>zeroTolerance) {
723          array[numberNonZero]=value;
724          index[numberNonZero++]=iColumn;
725        }
726        value = 0.0;
727        for (j=start; j<end;j++) {
728          int iRow = row[j];
729          value += pi[iRow]*elementByColumn[j];
730        }
731      }
732      value *= scale;
733      if (fabs(value)>zeroTolerance) {
734        array[numberNonZero]=value;
735        index[numberNonZero++]=iColumn;
736      }
737    }
738    // zero out
739    int numberRows = model->numberRows();
740    if (numberInRowArray*4<numberRows) {
741      for (i=0;i<numberInRowArray;i++) {
742        int iRow = whichRow[i];
743        pi[iRow]=0.0;
744      }
745    } else {
746      CoinZeroN(pi,numberRows);
747    }
748  } else {
749    if (!rowScale) {
750      if (scalar==-1.0) {
751        double value = 0.0;
752        CoinBigIndex j;
753        CoinBigIndex end = columnStart[1];
754        for (j=columnStart[0]; j<end;j++) {
755          int iRow = row[j];
756          value += pi[iRow]*elementByColumn[j];
757        }
758        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
759          CoinBigIndex start = end;
760          end = columnStart[iColumn+2];
761          if (fabs(value)>zeroTolerance) {
762            array[iColumn]=-value;
763            index[numberNonZero++]=iColumn;
764          }
765          value = 0.0;
766          for (j=start; j<end;j++) {
767            int iRow = row[j];
768            value += pi[iRow]*elementByColumn[j];
769          }
770        }
771        if (fabs(value)>zeroTolerance) {
772          array[iColumn]=-value;
773          index[numberNonZero++]=iColumn;
774        }
775      } else {
776        double value = 0.0;
777        CoinBigIndex j;
778        CoinBigIndex end = columnStart[1];
779        for (j=columnStart[0]; j<end;j++) {
780          int iRow = row[j];
781          value += pi[iRow]*elementByColumn[j];
782        }
783        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
784          value *= scalar;
785          CoinBigIndex start = end;
786          end = columnStart[iColumn+2];
787          if (fabs(value)>zeroTolerance) {
788            array[iColumn]=value;
789            index[numberNonZero++]=iColumn;
790          }
791          value = 0.0;
792          for (j=start; j<end;j++) {
793            int iRow = row[j];
794            value += pi[iRow]*elementByColumn[j];
795          }
796        }
797        value *= scalar;
798        if (fabs(value)>zeroTolerance) {
799          array[iColumn]=value;
800          index[numberNonZero++]=iColumn;
801        }
802      }
803    } else {
804      // scaled
805      if (scalar==-1.0) {
806        const double * columnScale = model->columnScale();
807        double value = 0.0;
808        double scale=columnScale[0];
809        CoinBigIndex j;
810        CoinBigIndex end = columnStart[1];
811        for (j=columnStart[0]; j<end;j++) {
812          int iRow = row[j];
813          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
814        }
815        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
816          value *= scale;
817          CoinBigIndex start = end;
818          end = columnStart[iColumn+2];
819          scale = columnScale[iColumn+1];
820          if (fabs(value)>zeroTolerance) {
821            array[iColumn]=-value;
822            index[numberNonZero++]=iColumn;
823          }
824          value = 0.0;
825          for (j=start; j<end;j++) {
826            int iRow = row[j];
827            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
828          }
829        }
830        value *= scale;
831        if (fabs(value)>zeroTolerance) {
832          array[iColumn]=-value;
833          index[numberNonZero++]=iColumn;
834        }
835      } else {
836        const double * columnScale = model->columnScale();
837        double value = 0.0;
838        double scale=columnScale[0]*scalar;
839        CoinBigIndex j;
840        CoinBigIndex end = columnStart[1];
841        for (j=columnStart[0]; j<end;j++) {
842          int iRow = row[j];
843          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
844        }
845        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
846          value *= scale;
847          CoinBigIndex start = end;
848          end = columnStart[iColumn+2];
849          scale = columnScale[iColumn+1]*scalar;
850          if (fabs(value)>zeroTolerance) {
851            array[iColumn]=value;
852            index[numberNonZero++]=iColumn;
853          }
854          value = 0.0;
855          for (j=start; j<end;j++) {
856            int iRow = row[j];
857            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
858          }
859        }
860        value *= scale;
861        if (fabs(value)>zeroTolerance) {
862          array[iColumn]=value;
863          index[numberNonZero++]=iColumn;
864        }
865      }
866    }
867  }
868  columnArray->setNumElements(numberNonZero);
869  y->setNumElements(0);
870  if (packed)
871    columnArray->setPackedMode(true);
872}
873/* Return <code>x * A + y</code> in <code>z</code>.
874        Squashes small elements and knows about ClpSimplex */
875void 
876ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
877                              const CoinIndexedVector * rowArray,
878                              CoinIndexedVector * y,
879                              CoinIndexedVector * columnArray) const
880{
881  columnArray->clear();
882  double * pi = rowArray->denseVector();
883  int numberNonZero=0;
884  int * index = columnArray->getIndices();
885  double * array = columnArray->denseVector();
886  int numberInRowArray = rowArray->getNumElements();
887  // maybe I need one in OsiSimplex
888  double zeroTolerance = model->factorization()->zeroTolerance();
889  const int * column = getIndices();
890  const CoinBigIndex * rowStart = getVectorStarts();
891  const double * element = getElements();
892  const int * whichRow = rowArray->getIndices();
893  bool packed = rowArray->packedMode();
894  if (numberInRowArray>2) {
895    // do by rows
896    // ** Row copy is already scaled
897    int iRow;
898    int i;
899    int numberOriginal = 0;
900    if (packed) {
901      numberNonZero=0;
902      // and set up mark as char array
903      char * marked = (char *) (index+columnArray->capacity());
904      double * array2 = y->denseVector();
905#ifdef CLP_DEBUG
906      for (i=0;i<numberActiveColumns_;i++) {
907        assert(!marked[i]);
908        assert(!array2[i]);
909      }
910#endif     
911      for (i=0;i<numberInRowArray;i++) {
912        iRow = whichRow[i]; 
913        double value = pi[i]*scalar;
914        CoinBigIndex j;
915        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
916          int iColumn = column[j];
917          if (!marked[iColumn]) {
918            marked[iColumn]=1;
919            index[numberNonZero++]=iColumn;
920          }
921          array2[iColumn] += value*element[j];
922        }
923      }
924      // get rid of tiny values and zero out marked
925      numberOriginal=numberNonZero;
926      numberNonZero=0;
927      for (i=0;i<numberOriginal;i++) {
928        int iColumn = index[i];
929        if (marked[iColumn]) {
930          double value = array2[iColumn];
931          array2[iColumn]=0.0; 
932          marked[iColumn]=0;
933          if (fabs(value)>zeroTolerance) {
934            array[numberNonZero]=value;
935            index[numberNonZero++]=iColumn;
936          }
937        }
938      }
939    } else {
940      double * markVector = y->denseVector();
941      numberNonZero=0;
942      // and set up mark as char array
943      char * marked = (char *) markVector;
944      for (i=0;i<numberOriginal;i++) {
945        int iColumn = index[i];
946        marked[iColumn]=0;
947      }
948     
949      for (i=0;i<numberInRowArray;i++) {
950        iRow = whichRow[i]; 
951        double value = pi[iRow]*scalar;
952        CoinBigIndex j;
953        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
954          int iColumn = column[j];
955          if (!marked[iColumn]) {
956            marked[iColumn]=1;
957            index[numberNonZero++]=iColumn;
958          }
959          array[iColumn] += value*element[j];
960        }
961      }
962      // get rid of tiny values and zero out marked
963      numberOriginal=numberNonZero;
964      numberNonZero=0;
965      for (i=0;i<numberOriginal;i++) {
966        int iColumn = index[i];
967        marked[iColumn]=0;
968        if (fabs(array[iColumn])>zeroTolerance) {
969          index[numberNonZero++]=iColumn;
970        } else {
971          array[iColumn]=0.0;
972        }
973      }
974    }
975  } else if (numberInRowArray==2) {
976    // do by rows when two rows
977    int numberOriginal;
978    int i;
979    CoinBigIndex j;
980    numberNonZero=0;
981
982    double value;
983    if (packed) {
984      int iRow0 = whichRow[0]; 
985      int iRow1 = whichRow[1]; 
986      double pi0 = pi[0];
987      double pi1 = pi[1];
988      if (rowStart[iRow0+1]-rowStart[iRow0]>
989          rowStart[iRow1+1]-rowStart[iRow1]) {
990        // do one with fewer first
991        iRow0=iRow1;
992        iRow1=whichRow[0];
993        pi0=pi1;
994        pi1=pi[0];
995      }
996      // and set up mark as char array
997      char * marked = (char *) (index+columnArray->capacity());
998      int * lookup = y->getIndices();
999      value = pi0*scalar;
1000      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
1001        int iColumn = column[j];
1002        double value2 = value*element[j];
1003        array[numberNonZero] = value2;
1004        marked[iColumn]=1;
1005        lookup[iColumn]=numberNonZero;
1006        index[numberNonZero++]=iColumn;
1007      }
1008      numberOriginal = numberNonZero;
1009      value = pi1*scalar;
1010      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
1011        int iColumn = column[j];
1012        double value2 = value*element[j];
1013        // I am assuming no zeros in matrix
1014        if (marked[iColumn]) {
1015          int iLookup = lookup[iColumn];
1016          array[iLookup] += value2;
1017        } else {
1018          if (fabs(value2)>zeroTolerance) {
1019            array[numberNonZero] = value2;
1020            index[numberNonZero++]=iColumn;
1021          }
1022        }
1023      }
1024      // get rid of tiny values and zero out marked
1025      int nDelete=0;
1026      for (i=0;i<numberOriginal;i++) {
1027        int iColumn = index[i];
1028        marked[iColumn]=0;
1029        if (fabs(array[i])<=zeroTolerance) 
1030          nDelete++;
1031      }
1032      if (nDelete) {
1033        numberOriginal=numberNonZero;
1034        numberNonZero=0;
1035        for (i=0;i<numberOriginal;i++) {
1036          int iColumn = index[i];
1037          double value = array[i];
1038          array[i]=0.0;
1039          if (fabs(value)>zeroTolerance) {
1040            array[numberNonZero]=value;
1041            index[numberNonZero++]=iColumn;
1042          }
1043        }
1044      }
1045    } else {
1046      int iRow = whichRow[0]; 
1047      value = pi[iRow]*scalar;
1048      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1049        int iColumn = column[j];
1050        double value2 = value*element[j];
1051        index[numberNonZero++]=iColumn;
1052        array[iColumn] = value2;
1053      }
1054      iRow = whichRow[1]; 
1055      value = pi[iRow]*scalar;
1056      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1057        int iColumn = column[j];
1058        double value2 = value*element[j];
1059        // I am assuming no zeros in matrix
1060        if (array[iColumn])
1061          value2 += array[iColumn];
1062        else
1063          index[numberNonZero++]=iColumn;
1064        array[iColumn] = value2;
1065      }
1066      // get rid of tiny values and zero out marked
1067      numberOriginal=numberNonZero;
1068      numberNonZero=0;
1069      for (i=0;i<numberOriginal;i++) {
1070        int iColumn = index[i];
1071        if (fabs(array[iColumn])>zeroTolerance) {
1072          index[numberNonZero++]=iColumn;
1073        } else {
1074          array[iColumn]=0.0;
1075        }
1076      }
1077    }
1078  } else if (numberInRowArray==1) {
1079    // Just one row
1080    int iRow=rowArray->getIndices()[0];
1081    numberNonZero=0;
1082    CoinBigIndex j;
1083    if (packed) {
1084      double value = pi[0]*scalar;
1085      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1086        int iColumn = column[j];
1087        double value2 = value*element[j];
1088        if (fabs(value2)>zeroTolerance) {
1089          array[numberNonZero] = value2;
1090          index[numberNonZero++]=iColumn;
1091        }
1092      }
1093    } else {
1094      double value = pi[iRow]*scalar;
1095      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1096        int iColumn = column[j];
1097        double value2 = value*element[j];
1098        if (fabs(value2)>zeroTolerance) {
1099          index[numberNonZero++]=iColumn;
1100          array[iColumn] = value2;
1101        }
1102      }
1103    }
1104  }
1105  columnArray->setNumElements(numberNonZero);
1106  y->setNumElements(0);
1107}
1108/* Return <code>x *A in <code>z</code> but
1109   just for indices in y.
1110   Squashes small elements and knows about ClpSimplex */
1111void 
1112ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
1113                              const CoinIndexedVector * rowArray,
1114                              const CoinIndexedVector * y,
1115                              CoinIndexedVector * columnArray) const
1116{
1117  columnArray->clear();
1118  double * pi = rowArray->denseVector();
1119  double * array = columnArray->denseVector();
1120  int jColumn;
1121  // get matrix data pointers
1122  const int * row = matrix_->getIndices();
1123  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1124  const int * columnLength = matrix_->getVectorLengths(); 
1125  const double * elementByColumn = matrix_->getElements();
1126  const double * rowScale = model->rowScale();
1127  int numberToDo = y->getNumElements();
1128  const int * which = y->getIndices();
1129  assert (!rowArray->packedMode());
1130  columnArray->setPacked();
1131  if (!hasGaps_&&numberToDo>5) {
1132    // no gaps and a reasonable number
1133    if (!rowScale) {
1134      int iColumn = which[0];
1135      double value = 0.0;
1136      CoinBigIndex j;
1137      for (j=columnStart[iColumn];
1138           j<columnStart[iColumn+1];j++) {
1139        int iRow = row[j];
1140        value += pi[iRow]*elementByColumn[j];
1141      }
1142      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
1143        int iColumn = which[jColumn+1];
1144        CoinBigIndex start = columnStart[iColumn];
1145        CoinBigIndex end = columnStart[iColumn+1];
1146        array[jColumn]=value;
1147        value = 0.0;
1148        for (j=start; j<end;j++) {
1149          int iRow = row[j];
1150          value += pi[iRow]*elementByColumn[j];
1151        }
1152      }
1153      array[jColumn]=value;
1154    } else {
1155      // scaled
1156      const double * columnScale = model->columnScale();
1157      int iColumn = which[0];
1158      double value = 0.0;
1159      double scale=columnScale[iColumn];
1160      CoinBigIndex j;
1161      for (j=columnStart[iColumn];
1162           j<columnStart[iColumn+1];j++) {
1163        int iRow = row[j];
1164        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1165      }
1166      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
1167        int iColumn = which[jColumn+1];
1168        value *= scale;
1169        scale = columnScale[iColumn];
1170        CoinBigIndex start = columnStart[iColumn];
1171        CoinBigIndex end = columnStart[iColumn+1];
1172        array[jColumn]=value;
1173        value = 0.0;
1174        for (j=start; j<end;j++) {
1175          int iRow = row[j];
1176          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1177        }
1178      }
1179      value *= scale;
1180      array[jColumn]=value;
1181    }
1182  } else {
1183    // gaps
1184    if (!rowScale) {
1185      for (jColumn=0;jColumn<numberToDo;jColumn++) {
1186        int iColumn = which[jColumn];
1187        double value = 0.0;
1188        CoinBigIndex j;
1189        for (j=columnStart[iColumn];
1190             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1191          int iRow = row[j];
1192          value += pi[iRow]*elementByColumn[j];
1193        }
1194        array[jColumn]=value;
1195      }
1196    } else {
1197      // scaled
1198      const double * columnScale = model->columnScale();
1199      for (jColumn=0;jColumn<numberToDo;jColumn++) {
1200        int iColumn = which[jColumn];
1201        double value = 0.0;
1202        CoinBigIndex j;
1203        for (j=columnStart[iColumn];
1204             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1205          int iRow = row[j];
1206          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1207        }
1208        value *= columnScale[iColumn];
1209        array[jColumn]=value;
1210      }
1211    }
1212  }
1213}
1214/* Returns true if can combine transposeTimes and subsetTransposeTimes
1215   and if it would be faster */
1216bool 
1217ClpPackedMatrix::canCombine(const ClpSimplex * model,
1218                            const CoinIndexedVector * pi) const
1219{
1220  int numberInRowArray = pi->getNumElements();
1221  int numberRows = model->numberRows();
1222  bool packed = pi->packedMode();
1223  // factor should be smaller if doing both with two pi vectors
1224  double factor = 0.27;
1225  // We may not want to do by row if there may be cache problems
1226  // It would be nice to find L2 cache size - for moment 512K
1227  // Be slightly optimistic
1228  if (numberActiveColumns_*sizeof(double)>1000000) {
1229    if (numberRows*10<numberActiveColumns_)
1230      factor *= 0.333333333;
1231    else if (numberRows*4<numberActiveColumns_)
1232      factor *= 0.5;
1233    else if (numberRows*2<numberActiveColumns_)
1234      factor *= 0.66666666667;
1235    //if (model->numberIterations()%50==0)
1236    //printf("%d nonzero\n",numberInRowArray);
1237  }
1238  // if not packed then bias a bit more towards by column
1239  if (!packed)
1240    factor *= 0.9;
1241  return ((numberInRowArray>factor*numberRows||!model->rowCopy())&&!hasGaps_);
1242}
1243#ifndef CLP_ALL_ONE_FILE
1244// These have to match ClpPrimalColumnSteepest version
1245#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
1246#endif
1247// Updates two arrays for steepest
1248void 
1249ClpPackedMatrix::transposeTimes2(const ClpSimplex * model,
1250                                 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
1251                                 const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
1252                                 CoinIndexedVector * spare,
1253                                 double referenceIn, double devex,
1254                                 // Array for exact devex to say what is in reference framework
1255                                 unsigned int * reference,
1256                                 double * weights, double scaleFactor)
1257{
1258  // put row of tableau in dj1
1259  double * pi = pi1->denseVector();
1260  int numberNonZero=0;
1261  int * index = dj1->getIndices();
1262  double * array = dj1->denseVector();
1263  int numberInRowArray = pi1->getNumElements();
1264  double zeroTolerance = model->factorization()->zeroTolerance();
1265  bool packed = pi1->packedMode();
1266  // do by column
1267  int iColumn;
1268  // get matrix data pointers
1269  const int * row = matrix_->getIndices();
1270  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1271  const double * elementByColumn = matrix_->getElements();
1272  const double * rowScale = model->rowScale();
1273  assert (!spare->getNumElements());
1274  assert (numberActiveColumns_>0);
1275  double * piWeight = pi2->denseVector();
1276  assert (!pi2->packedMode());
1277  bool killDjs = (scaleFactor==0.0);
1278  if (!scaleFactor)
1279    scaleFactor=1.0;
1280  if (packed) {
1281    // need to expand pi into y
1282    assert(spare->capacity()>=model->numberRows());
1283    double * piOld = pi;
1284    pi = spare->denseVector();
1285    const int * whichRow = pi1->getIndices();
1286    int i;
1287    if (!rowScale) {
1288      // modify pi so can collapse to one loop
1289      for (i=0;i<numberInRowArray;i++) {
1290        int iRow = whichRow[i];
1291        pi[iRow]=piOld[i];
1292      }
1293      CoinBigIndex j;
1294      CoinBigIndex end = columnStart[0];
1295      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
1296        CoinBigIndex start = end;
1297        end = columnStart[iColumn+1];
1298        ClpSimplex::Status status = model->getStatus(iColumn);
1299        if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
1300        double value = 0.0;
1301        for (j=start; j<end;j++) {
1302          int iRow = row[j];
1303          value -= pi[iRow]*elementByColumn[j];
1304        }
1305        if (fabs(value)>zeroTolerance) {
1306          // and do other array
1307          double modification = 0.0;
1308          for (j=start; j<end;j++) {
1309            int iRow = row[j];
1310            modification += piWeight[iRow]*elementByColumn[j];
1311          }
1312          double thisWeight = weights[iColumn];
1313          double pivot = value*scaleFactor;
1314          double pivotSquared = pivot * pivot;
1315          thisWeight += pivotSquared * devex + pivot * modification;
1316          if (thisWeight<DEVEX_TRY_NORM) {
1317            if (referenceIn<0.0) {
1318              // steepest
1319              thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
1320            } else {
1321              // exact
1322              thisWeight = referenceIn*pivotSquared;
1323              if (reference(iColumn))
1324                thisWeight += 1.0;
1325              thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
1326            }
1327          }
1328          weights[iColumn] = thisWeight;
1329          if (!killDjs) {
1330            array[numberNonZero]=value;
1331            index[numberNonZero++]=iColumn;
1332          }
1333        }
1334      }
1335    } else {
1336      // scaled
1337      // modify pi so can collapse to one loop
1338      for (i=0;i<numberInRowArray;i++) {
1339        int iRow = whichRow[i];
1340        pi[iRow]=piOld[i]*rowScale[iRow];
1341      }
1342      const double * columnScale = model->columnScale();
1343      CoinBigIndex j;
1344      CoinBigIndex end = columnStart[0];
1345      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
1346        CoinBigIndex start = end;
1347        end = columnStart[iColumn+1];
1348        ClpSimplex::Status status = model->getStatus(iColumn);
1349        if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
1350        double scale=columnScale[iColumn];
1351        double value = 0.0;
1352        for (j=start; j<end;j++) {
1353          int iRow = row[j];
1354          value -= pi[iRow]*elementByColumn[j];
1355        }
1356        value *= scale;
1357        if (fabs(value)>zeroTolerance) {
1358          double modification = 0.0;
1359          for (j=start; j<end;j++) {
1360            int iRow = row[j];
1361            modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
1362          }
1363          modification *= scale;
1364          double thisWeight = weights[iColumn];
1365          double pivot = value*scaleFactor;
1366          double pivotSquared = pivot * pivot;
1367          thisWeight += pivotSquared * devex + pivot * modification;
1368          if (thisWeight<DEVEX_TRY_NORM) {
1369            if (referenceIn<0.0) {
1370              // steepest
1371              thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
1372            } else {
1373              // exact
1374              thisWeight = referenceIn*pivotSquared;
1375              if (reference(iColumn))
1376                thisWeight += 1.0;
1377              thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
1378            }
1379          }
1380          weights[iColumn] = thisWeight;
1381          if (!killDjs) {
1382            array[numberNonZero]=value;
1383            index[numberNonZero++]=iColumn;
1384          }
1385        }
1386      }
1387    }
1388    // zero out
1389    int numberRows = model->numberRows();
1390    if (numberInRowArray*4<numberRows) {
1391      for (i=0;i<numberInRowArray;i++) {
1392        int iRow = whichRow[i];
1393        pi[iRow]=0.0;
1394      }
1395    } else {
1396      CoinZeroN(pi,numberRows);
1397    }
1398  } else {
1399    if (!rowScale) {
1400      CoinBigIndex j;
1401      CoinBigIndex end = columnStart[0];
1402      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
1403        CoinBigIndex start = end;
1404        end = columnStart[iColumn+1];
1405        ClpSimplex::Status status = model->getStatus(iColumn);
1406        if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
1407        double value = 0.0;
1408        for (j=start; j<end;j++) {
1409          int iRow = row[j];
1410          value -= pi[iRow]*elementByColumn[j];
1411        }
1412        if (fabs(value)>zeroTolerance) {
1413          // and do other array
1414          double modification = 0.0;
1415          for (j=start; j<end;j++) {
1416            int iRow = row[j];
1417            modification += piWeight[iRow]*elementByColumn[j];
1418          }
1419          double thisWeight = weights[iColumn];
1420          double pivot = value*scaleFactor;
1421          double pivotSquared = pivot * pivot;
1422          thisWeight += pivotSquared * devex + pivot * modification;
1423          if (thisWeight<DEVEX_TRY_NORM) {
1424            if (referenceIn<0.0) {
1425              // steepest
1426              thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
1427            } else {
1428              // exact
1429              thisWeight = referenceIn*pivotSquared;
1430              if (reference(iColumn))
1431                thisWeight += 1.0;
1432              thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
1433            }
1434          }
1435          weights[iColumn] = thisWeight;
1436          if (!killDjs) {
1437            array[iColumn]=value;
1438            index[numberNonZero++]=iColumn;
1439          }
1440        }
1441      }
1442    } else {
1443      // scaled
1444      const double * columnScale = model->columnScale();
1445      CoinBigIndex j;
1446      CoinBigIndex end = columnStart[0];
1447      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
1448        CoinBigIndex start = end;
1449        end = columnStart[iColumn+1];
1450        ClpSimplex::Status status = model->getStatus(iColumn);
1451        if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
1452        double scale=columnScale[iColumn];
1453        double value = 0.0;
1454        for (j=start; j<end;j++) {
1455          int iRow = row[j];
1456          value -= pi[iRow]*elementByColumn[j]*rowScale[iRow];
1457        }
1458        value *= scale;
1459        if (fabs(value)>zeroTolerance) {
1460          double modification = 0.0;
1461          for (j=start; j<end;j++) {
1462            int iRow = row[j];
1463            modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
1464          }
1465          modification *= scale;
1466          double thisWeight = weights[iColumn];
1467          double pivot = value*scaleFactor;
1468          double pivotSquared = pivot * pivot;
1469          thisWeight += pivotSquared * devex + pivot * modification;
1470          if (thisWeight<DEVEX_TRY_NORM) {
1471            if (referenceIn<0.0) {
1472              // steepest
1473              thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
1474            } else {
1475              // exact
1476              thisWeight = referenceIn*pivotSquared;
1477              if (reference(iColumn))
1478                thisWeight += 1.0;
1479              thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
1480            }
1481          }
1482          weights[iColumn] = thisWeight;
1483          if (!killDjs) {
1484            array[iColumn]=value;
1485            index[numberNonZero++]=iColumn;
1486          }
1487        }
1488      }
1489    }
1490  }
1491  dj1->setNumElements(numberNonZero);
1492  spare->setNumElements(0);
1493  if (packed)
1494    dj1->setPackedMode(true);
1495}
1496// Updates second array for steepest and does devex weights
1497void 
1498ClpPackedMatrix::subsetTimes2(const ClpSimplex * model,
1499                              CoinIndexedVector * dj1,
1500                            const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
1501                            double referenceIn, double devex,
1502                            // Array for exact devex to say what is in reference framework
1503                            unsigned int * reference,
1504                            double * weights, double scaleFactor)
1505{
1506  int number = dj1->getNumElements();
1507  const int * index = dj1->getIndices();
1508  double * array = dj1->denseVector();
1509  assert( dj1->packedMode());
1510
1511  // get matrix data pointers
1512  const int * row = matrix_->getIndices();
1513  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1514  const int * columnLength = matrix_->getVectorLengths(); 
1515  const double * elementByColumn = matrix_->getElements();
1516  const double * rowScale = model->rowScale();
1517  double * piWeight = pi2->denseVector();
1518  bool killDjs = (scaleFactor==0.0);
1519  if (!scaleFactor)
1520    scaleFactor=1.0;
1521  if (!rowScale) {
1522    for (int k=0;k<number;k++) {
1523      int iColumn = index[k];
1524      double pivot = array[k]*scaleFactor;
1525      if (killDjs)
1526        array[k]=0.0;
1527      // and do other array
1528      double modification = 0.0;
1529      for (CoinBigIndex j=columnStart[iColumn]; 
1530           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1531        int iRow = row[j];
1532        modification += piWeight[iRow]*elementByColumn[j];
1533      }
1534      double thisWeight = weights[iColumn];
1535      double pivotSquared = pivot * pivot;
1536      thisWeight += pivotSquared * devex + pivot * modification;
1537      if (thisWeight<DEVEX_TRY_NORM) {
1538        if (referenceIn<0.0) {
1539          // steepest
1540          thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
1541        } else {
1542          // exact
1543          thisWeight = referenceIn*pivotSquared;
1544          if (reference(iColumn))
1545            thisWeight += 1.0;
1546          thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
1547        }
1548      }
1549      weights[iColumn] = thisWeight;
1550    }
1551  } else {
1552    // scaled
1553    const double * columnScale = model->columnScale();
1554    for (int k=0;k<number;k++) {
1555      int iColumn = index[k];
1556      double pivot = array[k]*scaleFactor;
1557      double scale=columnScale[iColumn];
1558      if (killDjs)
1559        array[k]=0.0;
1560      // and do other array
1561      double modification = 0.0;
1562      for (CoinBigIndex j=columnStart[iColumn]; 
1563           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1564        int iRow = row[j];
1565        modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
1566      }
1567      double thisWeight = weights[iColumn];
1568      modification *= scale;
1569      double pivotSquared = pivot * pivot;
1570      thisWeight += pivotSquared * devex + pivot * modification;
1571      if (thisWeight<DEVEX_TRY_NORM) {
1572        if (referenceIn<0.0) {
1573          // steepest
1574          thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
1575        } else {
1576          // exact
1577          thisWeight = referenceIn*pivotSquared;
1578          if (reference(iColumn))
1579            thisWeight += 1.0;
1580          thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
1581        }
1582      }
1583      weights[iColumn] = thisWeight;
1584    }
1585  }
1586}
1587/// returns number of elements in column part of basis,
1588CoinBigIndex
1589ClpPackedMatrix::countBasis(ClpSimplex * model,
1590                           const int * whichColumn, 
1591                           int numberBasic,
1592                            int & numberColumnBasic)
1593{
1594  const int * columnLength = matrix_->getVectorLengths(); 
1595  int i;
1596  CoinBigIndex numberElements=0;
1597  // just count - can be over so ignore zero problem
1598  for (i=0;i<numberColumnBasic;i++) {
1599    int iColumn = whichColumn[i];
1600    numberElements += columnLength[iColumn];
1601  }
1602  return numberElements;
1603}
1604void
1605ClpPackedMatrix::fillBasis(ClpSimplex * model,
1606                         const int * whichColumn, 
1607                         int & numberColumnBasic,
1608                         int * indexRowU, int * start,
1609                         int * rowCount, int * columnCount,
1610                         double * elementU)
1611{
1612  const int * columnLength = matrix_->getVectorLengths(); 
1613  int i;
1614  CoinBigIndex numberElements=start[0];
1615  // fill
1616  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1617  const double * rowScale = model->rowScale();
1618  const int * row = matrix_->getIndices();
1619  const double * elementByColumn = matrix_->getElements();
1620  if (!zeroElements_) {
1621    if (!rowScale) {
1622      // no scaling
1623      for (i=0;i<numberColumnBasic;i++) {
1624        int iColumn = whichColumn[i];
1625        CoinBigIndex j;
1626        for (j=columnStart[iColumn];
1627             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1628          int iRow=row[j];
1629          indexRowU[numberElements]=iRow;
1630          rowCount[iRow]++;
1631          elementU[numberElements++]=elementByColumn[j];
1632        }
1633        start[i+1]=numberElements;
1634        columnCount[i]=columnLength[iColumn];
1635      }
1636    } else {
1637      // scaling
1638      const double * columnScale = model->columnScale();
1639      for (i=0;i<numberColumnBasic;i++) {
1640        int iColumn = whichColumn[i];
1641        CoinBigIndex j;
1642        double scale = columnScale[iColumn];
1643        for (j=columnStart[iColumn];
1644             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1645          int iRow = row[j];
1646          indexRowU[numberElements]=iRow;
1647          rowCount[iRow]++;
1648          elementU[numberElements++]=
1649            elementByColumn[j]*scale*rowScale[iRow];
1650        }
1651        start[i+1]=numberElements;
1652        columnCount[i]=columnLength[iColumn];
1653      }
1654    }
1655  } else {
1656    // there are zero elements so need to look more closely
1657    if (!rowScale) {
1658      // no scaling
1659      for (i=0;i<numberColumnBasic;i++) {
1660        int iColumn = whichColumn[i];
1661        CoinBigIndex j;
1662        for (j=columnStart[iColumn];
1663             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1664          double value = elementByColumn[j];
1665          if (value) {
1666            int iRow = row[j];
1667            indexRowU[numberElements]=iRow;
1668            rowCount[iRow]++;
1669            elementU[numberElements++]=value;
1670          }
1671        }
1672        start[i+1]=numberElements;
1673        columnCount[i]=numberElements-start[i];
1674      }
1675    } else {
1676      // scaling
1677      const double * columnScale = model->columnScale();
1678      for (i=0;i<numberColumnBasic;i++) {
1679        int iColumn = whichColumn[i];
1680        CoinBigIndex j;
1681        double scale = columnScale[iColumn];
1682        for (j=columnStart[iColumn];
1683             j<columnStart[iColumn]+columnLength[i];j++) {
1684          double value = elementByColumn[j];
1685          if (value) {
1686            int iRow = row[j];
1687            indexRowU[numberElements]=iRow;
1688            rowCount[iRow]++;
1689            elementU[numberElements++]=value*scale*rowScale[iRow];
1690          }
1691        }
1692        start[i+1]=numberElements;
1693        columnCount[i]=numberElements-start[i];
1694      }
1695    }
1696  }
1697}
1698// Creates scales for column copy (rowCopy in model may be modified)
1699int 
1700ClpPackedMatrix::scale(ClpModel * model) const 
1701{
1702  int numberRows = model->numberRows(); 
1703  int numberColumns = matrix_->getNumCols();
1704  // If empty - return as sanityCheck will trap
1705  if (!numberRows||!numberColumns)
1706    return 1;
1707  ClpMatrixBase * rowCopyBase=model->rowCopy();
1708  if (!rowCopyBase) {
1709    // temporary copy
1710    rowCopyBase = reverseOrderedCopy();
1711  }
1712#ifndef NO_RTTI
1713  ClpPackedMatrix* rowCopy =
1714    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
1715  // Make sure it is really a ClpPackedMatrix
1716  assert (rowCopy!=NULL);
1717#else
1718  ClpPackedMatrix* rowCopy =
1719    static_cast< ClpPackedMatrix*>(rowCopyBase);
1720#endif
1721
1722  const int * column = rowCopy->getIndices();
1723  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1724  const double * element = rowCopy->getElements();
1725  double * rowScale = new double [numberRows];
1726  double * columnScale = new double [numberColumns];
1727  // we are going to mark bits we are interested in
1728  char * usefulRow = new char [numberRows];
1729  char * usefulColumn = new char [numberColumns];
1730  double * rowLower = model->rowLower();
1731  double * rowUpper = model->rowUpper();
1732  double * columnLower = model->columnLower();
1733  double * columnUpper = model->columnUpper();
1734  int iColumn, iRow;
1735  //#define LEAVE_FIXED
1736  // mark free rows
1737  for (iRow=0;iRow<numberRows;iRow++) {
1738#ifndef LEAVE_FIXED
1739    usefulRow[iRow]=0;
1740    if (rowUpper[iRow]<1.0e20||
1741        rowLower[iRow]>-1.0e20)
1742#endif
1743      usefulRow[iRow]=1;
1744  }
1745  // mark empty and fixed columns
1746  // also see if worth scaling
1747  assert (model->scalingFlag()<4); // dynamic not implemented
1748  double largest=0.0;
1749  double smallest=1.0e50;
1750  // get matrix data pointers
1751  const int * row = matrix_->getIndices();
1752  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1753  const int * columnLength = matrix_->getVectorLengths(); 
1754  const double * elementByColumn = matrix_->getElements();
1755  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1756    CoinBigIndex j;
1757    char useful=0;
1758#ifndef LEAVE_FIXED
1759    if (columnUpper[iColumn]>
1760        columnLower[iColumn]+1.0e-12) {
1761#endif
1762      for (j=columnStart[iColumn];
1763           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1764        iRow=row[j];
1765        if(elementByColumn[j]&&usefulRow[iRow]) {
1766          useful=1;
1767          largest = CoinMax(largest,fabs(elementByColumn[j]));
1768          smallest = CoinMin(smallest,fabs(elementByColumn[j]));
1769        }
1770      }
1771#ifndef LEAVE_FIXED
1772    }
1773#endif
1774    usefulColumn[iColumn]=useful;
1775  }
1776  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
1777    <<smallest<<largest
1778    <<CoinMessageEol;
1779  if (smallest>=0.5&&largest<=2.0) {
1780    // don't bother scaling
1781    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
1782      <<CoinMessageEol;
1783    delete [] rowScale;
1784    delete [] usefulRow;
1785    delete [] columnScale;
1786    delete [] usefulColumn;
1787    if (!model->rowCopy()) 
1788      delete rowCopyBase; // get rid of temporary row copy
1789    return 1;
1790  } else {
1791      // need to scale
1792    int scalingMethod = model->scalingFlag();
1793    // and see if there any empty rows
1794    for (iRow=0;iRow<numberRows;iRow++) {
1795      if (usefulRow[iRow]) {
1796        CoinBigIndex j;
1797        int useful=0;
1798        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1799          int iColumn = column[j];
1800          if (usefulColumn[iColumn]) {
1801            useful=1;
1802            break;
1803          }
1804        }
1805        usefulRow[iRow]=useful;
1806      }
1807    }
1808    double savedOverallRatio=0.0;
1809    double tolerance = 5.0*model->primalTolerance();
1810    double overallLargest=-1.0e-20;
1811    double overallSmallest=1.0e20;
1812    bool finished=false;
1813    // if scalingMethod 3 then may change
1814    while (!finished) {
1815      ClpFillN ( rowScale, numberRows,1.0);
1816      ClpFillN ( columnScale, numberColumns,1.0);
1817      overallLargest=-1.0e-20;
1818      overallSmallest=1.0e20;
1819      if (scalingMethod==1||scalingMethod==3) {
1820        // Maximum in each row
1821        for (iRow=0;iRow<numberRows;iRow++) {
1822          if (usefulRow[iRow]) {
1823            CoinBigIndex j;
1824            largest=1.0e-10;
1825            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1826              int iColumn = column[j];
1827              if (usefulColumn[iColumn]) {
1828                double value = fabs(element[j]);
1829                largest = CoinMax(largest,value);
1830                assert (largest<1.0e40);
1831              }
1832            }
1833            rowScale[iRow]=1.0/largest;
1834            overallLargest = CoinMax(overallLargest,largest);
1835            overallSmallest = CoinMin(overallSmallest,largest);
1836          }
1837        }
1838      } else {
1839        int numberPass=3;
1840#ifdef USE_OBJECTIVE
1841        // This will be used to help get scale factors
1842        double * objective = new double[numberColumns];
1843        memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
1844        double objScale=1.0;
1845#endif
1846        while (numberPass) {
1847          overallLargest=0.0;
1848          overallSmallest=1.0e50;
1849          numberPass--;
1850          // Geometric mean on row scales
1851          for (iRow=0;iRow<numberRows;iRow++) {
1852            if (usefulRow[iRow]) {
1853              CoinBigIndex j;
1854              largest=1.0e-20;
1855              smallest=1.0e50;
1856              for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1857                int iColumn = column[j];
1858                if (usefulColumn[iColumn]) {
1859                  double value = fabs(element[j]);
1860                  // Don't bother with tiny elements
1861                  if (value>1.0e-20) {
1862                    value *= columnScale[iColumn];
1863                    largest = CoinMax(largest,value);
1864                    smallest = CoinMin(smallest,value);
1865                  }
1866                }
1867              }
1868              rowScale[iRow]=1.0/sqrt(smallest*largest);
1869              rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
1870              overallLargest = CoinMax(largest*rowScale[iRow],overallLargest);
1871              overallSmallest = CoinMin(smallest*rowScale[iRow],overallSmallest);
1872            }
1873          }
1874#ifdef USE_OBJECTIVE
1875          largest=1.0e-20;
1876          smallest=1.0e50;
1877          for (iColumn=0;iColumn<numberColumns;iColumn++) {
1878            if (usefulColumn[iColumn]) {
1879              double value = fabs(objective[iColumn]);
1880              // Don't bother with tiny elements
1881              if (value>1.0e-20) {
1882                value *= columnScale[iColumn];
1883                largest = CoinMax(largest,value);
1884                smallest = CoinMin(smallest,value);
1885              }
1886            }
1887          }
1888          objScale=1.0/sqrt(smallest*largest);
1889#endif
1890          model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
1891            <<overallSmallest
1892            <<overallLargest
1893            <<CoinMessageEol;
1894          // skip last column round
1895          if (numberPass==1)
1896            break;
1897          // Geometric mean on column scales
1898          for (iColumn=0;iColumn<numberColumns;iColumn++) {
1899            if (usefulColumn[iColumn]) {
1900              CoinBigIndex j;
1901              largest=1.0e-20;
1902              smallest=1.0e50;
1903              for (j=columnStart[iColumn];
1904                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
1905                iRow=row[j];
1906                double value = fabs(elementByColumn[j]);
1907                // Don't bother with tiny elements
1908                if (value>1.0e-20&&usefulRow[iRow]) {
1909                  value *= rowScale[iRow];
1910                  largest = CoinMax(largest,value);
1911                  smallest = CoinMin(smallest,value);
1912                }
1913              }
1914#ifdef USE_OBJECTIVE
1915              if (fabs(objective[iColumn])>1.0e-20) {
1916                double value = fabs(objective[iColumn])*objScale;
1917                largest = CoinMax(largest,value);
1918                smallest = CoinMin(smallest,value);
1919              }
1920#endif
1921              columnScale[iColumn]=1.0/sqrt(smallest*largest);
1922              columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
1923            }
1924          }
1925        }
1926#ifdef USE_OBJECTIVE
1927        delete [] objective;
1928        printf("obj scale %g - use it if you want\n",objScale);
1929#endif
1930      }
1931      // If ranges will make horrid then scale
1932      for (iRow=0;iRow<numberRows;iRow++) {
1933        if (usefulRow[iRow]) {
1934          double difference = rowUpper[iRow]-rowLower[iRow];
1935          double scaledDifference = difference*rowScale[iRow];
1936          if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1937            // make gap larger
1938            rowScale[iRow] *= 1.0e-4/scaledDifference;
1939            rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
1940            //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
1941            // scaledDifference,difference*rowScale[iRow]);
1942          }
1943        }
1944      }
1945      // final pass to scale columns so largest is reasonable
1946      // See what smallest will be if largest is 1.0
1947      overallSmallest=1.0e50;
1948      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1949        if (usefulColumn[iColumn]) {
1950          CoinBigIndex j;
1951          largest=1.0e-20;
1952          smallest=1.0e50;
1953          for (j=columnStart[iColumn];
1954               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1955            iRow=row[j];
1956            if(elementByColumn[j]&&usefulRow[iRow]) {
1957              double value = fabs(elementByColumn[j]*rowScale[iRow]);
1958              largest = CoinMax(largest,value);
1959              smallest = CoinMin(smallest,value);
1960            }
1961          }
1962          if (overallSmallest*largest>smallest)
1963            overallSmallest = smallest/largest;
1964        }
1965      }
1966      if (scalingMethod==1||scalingMethod==2) {
1967        finished=true;
1968      } else if (savedOverallRatio==0.0&&scalingMethod!=4) {
1969        savedOverallRatio=overallSmallest;
1970        scalingMethod=4;
1971      } else {
1972        assert (scalingMethod==4);
1973        if (overallSmallest>2.0*savedOverallRatio)
1974          finished=true; // geometric was better
1975        else
1976          scalingMethod=1; // redo equilibrium
1977#if 0
1978        if (model->logLevel()>2) {
1979          if (finished)
1980            printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
1981                   savedOverallRatio,overallSmallest);
1982          else
1983            printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
1984                   savedOverallRatio,overallSmallest);
1985        }
1986#endif
1987      }
1988    }
1989    //#define RANDOMIZE
1990#ifdef RANDOMIZE
1991    // randomize by up to 10%
1992    for (iRow=0;iRow<numberRows;iRow++) {
1993      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1994      rowScale[iRow] *= (1.0+0.1*value);
1995    }
1996#endif
1997    overallLargest=1.0;
1998    if (overallSmallest<1.0e-1)
1999      overallLargest = 1.0/sqrt(overallSmallest);
2000    overallLargest = CoinMin(100.0,overallLargest);
2001    overallSmallest=1.0e50;
2002    //printf("scaling %d\n",model->scalingFlag());
2003    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2004      if (columnUpper[iColumn]>
2005        columnLower[iColumn]+1.0e-12) {
2006        //if (usefulColumn[iColumn]) {
2007        CoinBigIndex j;
2008        largest=1.0e-20;
2009        smallest=1.0e50;
2010        for (j=columnStart[iColumn];
2011             j<columnStart[iColumn]+columnLength[iColumn];j++) {
2012          iRow=row[j];
2013          if(elementByColumn[j]&&usefulRow[iRow]) {
2014            double value = fabs(elementByColumn[j]*rowScale[iRow]);
2015            largest = CoinMax(largest,value);
2016            smallest = CoinMin(smallest,value);
2017          }
2018        }
2019        columnScale[iColumn]=overallLargest/largest;
2020        columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
2021#ifdef RANDOMIZE
2022        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
2023        columnScale[iColumn] *= (1.0+0.1*value);
2024#endif
2025        double difference = columnUpper[iColumn]-columnLower[iColumn];
2026        if (difference<1.0e-5*columnScale[iColumn]) {
2027          // make gap larger
2028          columnScale[iColumn] = difference/1.0e-5;
2029          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
2030          // scaledDifference,difference*columnScale[iColumn]);
2031        }
2032        overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
2033      }
2034    }
2035    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
2036      <<overallSmallest
2037      <<overallLargest
2038      <<CoinMessageEol;
2039    delete [] usefulRow;
2040    delete [] usefulColumn;
2041#ifndef SLIM_CLP
2042    // If quadratic then make symmetric
2043    ClpObjective * obj = model->objectiveAsObject();
2044#ifndef NO_RTTI
2045    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2046#else
2047    ClpQuadraticObjective * quadraticObj = NULL;
2048    if (obj->type()==2)
2049      quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
2050#endif
2051    if (quadraticObj) {
2052      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2053      int numberXColumns = quadratic->getNumCols();
2054      if (numberXColumns<numberColumns) {
2055        // we assume symmetric
2056        int numberQuadraticColumns=0;
2057        int i;
2058        //const int * columnQuadratic = quadratic->getIndices();
2059        //const int * columnQuadraticStart = quadratic->getVectorStarts();
2060        const int * columnQuadraticLength = quadratic->getVectorLengths();
2061        for (i=0;i<numberXColumns;i++) {
2062          int length=columnQuadraticLength[i];
2063#ifndef CORRECT_COLUMN_COUNTS
2064          length=1;
2065#endif
2066          if (length)
2067            numberQuadraticColumns++;
2068        }
2069        int numberXRows = numberRows-numberQuadraticColumns;
2070        numberQuadraticColumns=0;
2071        for (i=0;i<numberXColumns;i++) { 
2072          int length=columnQuadraticLength[i];
2073#ifndef CORRECT_COLUMN_COUNTS
2074          length=1;
2075#endif
2076          if (length) {
2077            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
2078            numberQuadraticColumns++;
2079          }
2080        }   
2081        int numberQuadraticRows=0;
2082        for (i=0;i<numberXRows;i++) {
2083          // See if any in row quadratic
2084          CoinBigIndex j;
2085          int numberQ=0;
2086          for (j=rowStart[i];j<rowStart[i+1];j++) {
2087            int iColumn = column[j];
2088            if (columnQuadraticLength[iColumn])
2089              numberQ++;
2090          }
2091#ifndef CORRECT_ROW_COUNTS
2092          numberQ=1;
2093#endif
2094          if (numberQ) {
2095            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
2096            numberQuadraticRows++;
2097          }
2098        }
2099        // and make sure Sj okay
2100        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
2101          CoinBigIndex j=columnStart[iColumn];
2102          assert(columnLength[iColumn]==1);
2103          int iRow=row[j];
2104          double value = fabs(elementByColumn[j]*rowScale[iRow]);
2105          columnScale[iColumn]=1.0/value;
2106        }
2107      }
2108    }
2109#endif
2110    model->setRowScale(rowScale);
2111    model->setColumnScale(columnScale);
2112    if (model->rowCopy()) {
2113      // need to replace row by row
2114      double * newElement = new double[numberColumns];
2115      // scale row copy
2116      for (iRow=0;iRow<numberRows;iRow++) {
2117        CoinBigIndex j;
2118        double scale = rowScale[iRow];
2119        const double * elementsInThisRow = element + rowStart[iRow];
2120        const int * columnsInThisRow = column + rowStart[iRow];
2121        int number = rowStart[iRow+1]-rowStart[iRow];
2122        assert (number<=numberColumns);
2123        for (j=0;j<number;j++) {
2124          int iColumn = columnsInThisRow[j];
2125          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
2126        }
2127        rowCopy->replaceVector(iRow,number,newElement);
2128      }
2129      delete [] newElement;
2130    } else {
2131      // no row copy
2132      delete rowCopyBase;
2133    }
2134    return 0;
2135  }
2136}
2137/* Unpacks a column into an CoinIndexedvector
2138 */
2139void 
2140ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
2141                   int iColumn) const 
2142{
2143  const double * rowScale = model->rowScale();
2144  const int * row = matrix_->getIndices();
2145  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2146  const int * columnLength = matrix_->getVectorLengths(); 
2147  const double * elementByColumn = matrix_->getElements();
2148  CoinBigIndex i;
2149  if (!rowScale) {
2150    for (i=columnStart[iColumn];
2151         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2152      rowArray->add(row[i],elementByColumn[i]);
2153    }
2154  } else {
2155    // apply scaling
2156    double scale = model->columnScale()[iColumn];
2157    for (i=columnStart[iColumn];
2158         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2159      int iRow = row[i];
2160      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
2161    }
2162  }
2163}
2164/* Unpacks a column into a CoinIndexedVector
2165** in packed format
2166Note that model is NOT const.  Bounds and objective could
2167be modified if doing column generation (just for this variable) */
2168void 
2169ClpPackedMatrix::unpackPacked(ClpSimplex * model,
2170                            CoinIndexedVector * rowArray,
2171                            int iColumn) const
2172{
2173  const double * rowScale = model->rowScale();
2174  const int * row = matrix_->getIndices();
2175  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2176  const int * columnLength = matrix_->getVectorLengths(); 
2177  const double * elementByColumn = matrix_->getElements();
2178  CoinBigIndex i;
2179  if (!rowScale) {
2180    CoinBigIndex j=columnStart[iColumn];
2181    rowArray->createPacked(columnLength[iColumn],
2182                           row+j,elementByColumn+j);
2183  } else {
2184    // apply scaling
2185    double scale = model->columnScale()[iColumn];
2186    int * index = rowArray->getIndices();
2187    double * array = rowArray->denseVector();
2188    int number = 0;
2189    for (i=columnStart[iColumn];
2190         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2191      int iRow = row[i];
2192      array[number]=elementByColumn[i]*scale*rowScale[iRow];
2193      index[number++]=iRow;
2194    }
2195    rowArray->setNumElements(number);
2196    rowArray->setPackedMode(true);
2197  }
2198}
2199/* Adds multiple of a column into an CoinIndexedvector
2200      You can use quickAdd to add to vector */
2201void 
2202ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
2203                   int iColumn, double multiplier) const 
2204{
2205  const double * rowScale = model->rowScale();
2206  const int * row = matrix_->getIndices();
2207  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2208  const int * columnLength = matrix_->getVectorLengths(); 
2209  const double * elementByColumn = matrix_->getElements();
2210  CoinBigIndex i;
2211  if (!rowScale) {
2212    for (i=columnStart[iColumn];
2213         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2214      int iRow = row[i];
2215      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
2216    }
2217  } else {
2218    // apply scaling
2219    double scale = model->columnScale()[iColumn]*multiplier;
2220    for (i=columnStart[iColumn];
2221         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2222      int iRow = row[i];
2223      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
2224    }
2225  }
2226}
2227/* Adds multiple of a column into an array */
2228void 
2229ClpPackedMatrix::add(const ClpSimplex * model,double * array,
2230                    int iColumn, double multiplier) const
2231{
2232  const double * rowScale = model->rowScale();
2233  const int * row = matrix_->getIndices();
2234  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2235  const int * columnLength = matrix_->getVectorLengths(); 
2236  const double * elementByColumn = matrix_->getElements();
2237  CoinBigIndex i;
2238  if (!rowScale) {
2239    for (i=columnStart[iColumn];
2240         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2241      int iRow = row[i];
2242      array[iRow] += multiplier*elementByColumn[i];
2243    }
2244  } else {
2245    // apply scaling
2246    double scale = model->columnScale()[iColumn]*multiplier;
2247    for (i=columnStart[iColumn];
2248         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2249      int iRow = row[i];
2250      array[iRow] += elementByColumn[i]*scale*rowScale[iRow];
2251    }
2252  }
2253}
2254/* Checks if all elements are in valid range.  Can just
2255   return true if you are not paranoid.  For Clp I will
2256   probably expect no zeros.  Code can modify matrix to get rid of
2257   small elements.
2258   check bits (can be turned off to save time) :
2259   1 - check if matrix has gaps
2260   2 - check if zero elements
2261   4 - check and compress duplicates
2262   8 - report on large and small
2263*/
2264bool 
2265ClpPackedMatrix::allElementsInRange(ClpModel * model,
2266                                    double smallest, double largest,
2267                                    int check)
2268{
2269  int iColumn;
2270  // make sure matrix correct size
2271  matrix_->setDimensions(model->numberRows(),model->numberColumns());
2272  CoinBigIndex numberLarge=0;;
2273  CoinBigIndex numberSmall=0;;
2274  CoinBigIndex numberDuplicate=0;;
2275  int firstBadColumn=-1;
2276  int firstBadRow=-1;
2277  double firstBadElement=0.0;
2278  // get matrix data pointers
2279  const int * row = matrix_->getIndices();
2280  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2281  const int * columnLength = matrix_->getVectorLengths(); 
2282  const double * elementByColumn = matrix_->getElements();
2283  int numberRows = matrix_->getNumRows();
2284  int numberColumns = matrix_->getNumCols();
2285  // Say no gaps
2286  hasGaps_=false;
2287  if (check==14) {
2288    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2289      CoinBigIndex start = columnStart[iColumn];
2290      CoinBigIndex end = start + columnLength[iColumn];
2291      if (end!=columnStart[iColumn+1]) {
2292        hasGaps_=true;
2293        break;
2294      }
2295    }
2296    return true;
2297  }
2298  assert (check==15);
2299  int * mark = new int [numberRows];
2300  int i;
2301  for (i=0;i<numberRows;i++)
2302    mark[i]=-1;
2303  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2304    CoinBigIndex j;
2305    CoinBigIndex start = columnStart[iColumn];
2306    CoinBigIndex end = start + columnLength[iColumn];
2307    if (end!=columnStart[iColumn+1])
2308      hasGaps_=true;
2309    for (j=start;j<end;j++) {
2310      double value = fabs(elementByColumn[j]);
2311      int iRow = row[j];
2312      if (iRow<0||iRow>=numberRows) {
2313#ifndef COIN_BIG_INDEX
2314        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2315#elif COIN_BIG_INDEX==1
2316        printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2317#else
2318        printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2319#endif
2320        return false;
2321      }
2322      if (mark[iRow]==-1) {
2323        mark[iRow]=j;
2324      } else {
2325        // duplicate
2326        numberDuplicate++;
2327      }
2328      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2329      if (!value)
2330        zeroElements_ = true; // there are zero elements
2331      if (value<smallest) {
2332        numberSmall++;
2333      } else if (!(value<=largest)) {
2334        numberLarge++;
2335        if (firstBadColumn<0) {
2336          firstBadColumn=iColumn;
2337          firstBadRow=row[j];
2338          firstBadElement=elementByColumn[j];
2339        }
2340      }
2341    }
2342    //clear mark
2343    for (j=columnStart[iColumn];
2344         j<columnStart[iColumn]+columnLength[iColumn];j++) {
2345      int iRow = row[j];
2346      mark[iRow]=-1;
2347    }
2348  }
2349  delete [] mark;
2350  if (numberLarge) {
2351    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
2352      <<numberLarge
2353      <<firstBadColumn<<firstBadRow<<firstBadElement
2354      <<CoinMessageEol;
2355    return false;
2356  }
2357  if (numberSmall) 
2358    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
2359      <<numberSmall
2360      <<CoinMessageEol;
2361  if (numberDuplicate) 
2362    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
2363      <<numberDuplicate
2364      <<CoinMessageEol;
2365  if (numberDuplicate) 
2366    matrix_->eliminateDuplicates(smallest);
2367  else if (numberSmall) 
2368    matrix_->compress(smallest);
2369  // If smallest >0.0 then there can't be zero elements
2370  if (smallest>0.0)
2371    zeroElements_=false;
2372  return true;
2373}
2374/* Given positive integer weights for each row fills in sum of weights
2375   for each column (and slack).
2376   Returns weights vector
2377*/
2378CoinBigIndex * 
2379ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
2380{
2381  int numberRows = model->numberRows();
2382  int numberColumns = matrix_->getNumCols();
2383  int number = numberRows+numberColumns;
2384  CoinBigIndex * weights = new CoinBigIndex[number];
2385  // get matrix data pointers
2386  const int * row = matrix_->getIndices();
2387  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2388  const int * columnLength = matrix_->getVectorLengths(); 
2389  int i;
2390  for (i=0;i<numberColumns;i++) {
2391    CoinBigIndex j;
2392    CoinBigIndex count=0;
2393    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
2394      int iRow=row[j];
2395      count += inputWeights[iRow];
2396    }
2397    weights[i]=count;
2398  }
2399  for (i=0;i<numberRows;i++) {
2400    weights[i+numberColumns]=inputWeights[i];
2401  }
2402  return weights;
2403}
2404/* Returns largest and smallest elements of both signs.
2405   Largest refers to largest absolute value.
2406*/
2407void 
2408ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
2409                       double & smallestPositive, double & largestPositive)
2410{
2411  smallestNegative=-COIN_DBL_MAX;
2412  largestNegative=0.0;
2413  smallestPositive=COIN_DBL_MAX;
2414  largestPositive=0.0;
2415  // get matrix data pointers
2416  const double * elementByColumn = matrix_->getElements();
2417  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2418  const int * columnLength = matrix_->getVectorLengths(); 
2419  int numberColumns = matrix_->getNumCols();
2420  int i;
2421  for (i=0;i<numberColumns;i++) {
2422    CoinBigIndex j;
2423    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
2424      double value = elementByColumn[j];
2425      if (value>0.0) {
2426        smallestPositive = CoinMin(smallestPositive,value);
2427        largestPositive = CoinMax(largestPositive,value);
2428      } else if (value<0.0) {
2429        smallestNegative = CoinMax(smallestNegative,value);
2430        largestNegative = CoinMin(largestNegative,value);
2431      }
2432    }
2433  }
2434}
2435// Says whether it can do partial pricing
2436bool 
2437ClpPackedMatrix::canDoPartialPricing() const
2438{
2439  return true;
2440}
2441// Partial pricing
2442void 
2443ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
2444                              int & bestSequence, int & numberWanted)
2445{
2446  numberWanted=currentWanted_;
2447  int start = (int) (startFraction*numberActiveColumns_);
2448  int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
2449  const double * element =matrix_->getElements();
2450  const int * row = matrix_->getIndices();
2451  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
2452  const int * length = matrix_->getVectorLengths();
2453  const double * rowScale = model->rowScale();
2454  const double * columnScale = model->columnScale();
2455  int iSequence;
2456  CoinBigIndex j;
2457  double tolerance=model->currentDualTolerance();
2458  double * reducedCost = model->djRegion();
2459  const double * duals = model->dualRowSolution();
2460  const double * cost = model->costRegion();
2461  double bestDj;
2462  if (bestSequence>=0)
2463    bestDj = fabs(model->clpMatrix()->reducedCost(model,bestSequence));
2464  else
2465    bestDj=tolerance;
2466  int sequenceOut = model->sequenceOut();
2467  int saveSequence = bestSequence;
2468  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 
2469  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
2470  if (rowScale) {
2471    // scaled
2472    for (iSequence=start;iSequence<end;iSequence++) {
2473      if (iSequence!=sequenceOut) {
2474        double value;
2475        ClpSimplex::Status status = model->getStatus(iSequence);
2476       
2477        switch(status) {
2478         
2479        case ClpSimplex::basic:
2480        case ClpSimplex::isFixed:
2481          break;
2482        case ClpSimplex::isFree:
2483        case ClpSimplex::superBasic:
2484          value=0.0;
2485          // scaled
2486          for (j=startColumn[iSequence];
2487               j<startColumn[iSequence]+length[iSequence];j++) {
2488            int jRow=row[j];
2489            value -= duals[jRow]*element[j]*rowScale[jRow];
2490          }
2491          value = fabs(cost[iSequence] +value*columnScale[iSequence]);
2492          if (value>FREE_ACCEPT*tolerance) {
2493            numberWanted--;
2494            // we are going to bias towards free (but only if reasonable)
2495            value *= FREE_BIAS;
2496            if (value>bestDj) {
2497              // check flagged variable and correct dj
2498              if (!model->flagged(iSequence)) {
2499                bestDj=value;
2500                bestSequence = iSequence;
2501              } else {
2502                // just to make sure we don't exit before got something
2503                numberWanted++;
2504              }
2505            }
2506          }
2507          break;
2508        case ClpSimplex::atUpperBound:
2509          value=0.0;
2510          // scaled
2511          for (j=startColumn[iSequence];
2512               j<startColumn[iSequence]+length[iSequence];j++) {
2513            int jRow=row[j];
2514            value -= duals[jRow]*element[j]*rowScale[jRow];
2515          }
2516          value = cost[iSequence] +value*columnScale[iSequence];
2517          if (value>tolerance) {
2518            numberWanted--;
2519            if (value>bestDj) {
2520              // check flagged variable and correct dj
2521              if (!model->flagged(iSequence)) {
2522                bestDj=value;
2523                bestSequence = iSequence;
2524              } else {
2525                // just to make sure we don't exit before got something
2526                numberWanted++;
2527              }
2528            }
2529          }
2530          break;
2531        case ClpSimplex::atLowerBound:
2532          value=0.0;
2533          // scaled
2534          for (j=startColumn[iSequence];
2535               j<startColumn[iSequence]+length[iSequence];j++) {
2536            int jRow=row[j];
2537            value -= duals[jRow]*element[j]*rowScale[jRow];
2538          }
2539          value = -(cost[iSequence] +value*columnScale[iSequence]);
2540          if (value>tolerance) {
2541            numberWanted--;
2542            if (value>bestDj) {
2543              // check flagged variable and correct dj
2544              if (!model->flagged(iSequence)) {
2545                bestDj=value;
2546                bestSequence = iSequence;
2547              } else {
2548                // just to make sure we don't exit before got something
2549                numberWanted++;
2550              }
2551            }
2552          }
2553          break;
2554        }
2555      }
2556      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2557        // give up
2558        break;
2559      }
2560      if (!numberWanted)
2561        break;
2562    }
2563    if (bestSequence!=saveSequence) {
2564      // recompute dj
2565      double value=0.0;
2566      // scaled
2567      for (j=startColumn[bestSequence];
2568           j<startColumn[bestSequence]+length[bestSequence];j++) {
2569        int jRow=row[j];
2570        value -= duals[jRow]*element[j]*rowScale[jRow];
2571      }
2572      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
2573      savedBestSequence_ = bestSequence;
2574      savedBestDj_ = reducedCost[savedBestSequence_];
2575    }
2576  } else {
2577    // not scaled
2578    for (iSequence=start;iSequence<end;iSequence++) {
2579      if (iSequence!=sequenceOut) {
2580        double value;
2581        ClpSimplex::Status status = model->getStatus(iSequence);
2582       
2583        switch(status) {
2584         
2585        case ClpSimplex::basic:
2586        case ClpSimplex::isFixed:
2587          break;
2588        case ClpSimplex::isFree:
2589        case ClpSimplex::superBasic:
2590          value=cost[iSequence];
2591          for (j=startColumn[iSequence];
2592               j<startColumn[iSequence]+length[iSequence];j++) {
2593            int jRow=row[j];
2594            value -= duals[jRow]*element[j];
2595          }
2596          value = fabs(value);
2597          if (value>FREE_ACCEPT*tolerance) {
2598            numberWanted--;
2599            // we are going to bias towards free (but only if reasonable)
2600            value *= FREE_BIAS;
2601            if (value>bestDj) {
2602              // check flagged variable and correct dj
2603              if (!model->flagged(iSequence)) {
2604                bestDj=value;
2605                bestSequence = iSequence;
2606              } else {
2607                // just to make sure we don't exit before got something
2608                numberWanted++;
2609              }
2610            }
2611          }
2612          break;
2613        case ClpSimplex::atUpperBound:
2614          value=cost[iSequence];
2615          // scaled
2616          for (j=startColumn[iSequence];
2617               j<startColumn[iSequence]+length[iSequence];j++) {
2618            int jRow=row[j];
2619            value -= duals[jRow]*element[j];
2620          }
2621          if (value>tolerance) {
2622            numberWanted--;
2623            if (value>bestDj) {
2624              // check flagged variable and correct dj
2625              if (!model->flagged(iSequence)) {
2626                bestDj=value;
2627                bestSequence = iSequence;
2628              } else {
2629                // just to make sure we don't exit before got something
2630                numberWanted++;
2631              }
2632            }
2633          }
2634          break;
2635        case ClpSimplex::atLowerBound:
2636          value=cost[iSequence];
2637          for (j=startColumn[iSequence];
2638               j<startColumn[iSequence]+length[iSequence];j++) {
2639            int jRow=row[j];
2640            value -= duals[jRow]*element[j];
2641          }
2642          value = -value;
2643          if (value>tolerance) {
2644            numberWanted--;
2645            if (value>bestDj) {
2646              // check flagged variable and correct dj
2647              if (!model->flagged(iSequence)) {
2648                bestDj=value;
2649                bestSequence = iSequence;
2650              } else {
2651                // just to make sure we don't exit before got something
2652                numberWanted++;
2653              }
2654            }
2655          }
2656          break;
2657        }
2658      }
2659      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2660        // give up
2661        break;
2662      }
2663      if (!numberWanted)
2664        break;
2665    }
2666    if (bestSequence!=saveSequence) {
2667      // recompute dj
2668      double value=cost[bestSequence];
2669      for (j=startColumn[bestSequence];
2670           j<startColumn[bestSequence]+length[bestSequence];j++) {
2671        int jRow=row[j];
2672        value -= duals[jRow]*element[j];
2673      }
2674      reducedCost[bestSequence] = value;
2675      savedBestSequence_ = bestSequence;
2676      savedBestDj_ = reducedCost[savedBestSequence_];
2677    }
2678  }
2679  currentWanted_=numberWanted;
2680}
2681// Sets up an effective RHS
2682void 
2683ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
2684{
2685  delete [] rhsOffset_;
2686  int numberRows = model->numberRows();
2687  rhsOffset_ = new double[numberRows];
2688  rhsOffset(model,true);
2689}
2690// makes sure active columns correct
2691int 
2692ClpPackedMatrix::refresh(ClpSimplex * model)
2693{
2694  numberActiveColumns_ = matrix_->getNumCols();
2695  return 0;
2696}
2697
2698/* Scales rowCopy if column copy scaled
2699   Only called if scales already exist */
2700void 
2701ClpPackedMatrix::scaleRowCopy(ClpModel * model) const 
2702{
2703  if (model->rowCopy()) {
2704    // need to replace row by row
2705    int numberRows = model->numberRows();
2706    int numberColumns = matrix_->getNumCols();
2707    double * newElement = new double[numberColumns];
2708    ClpMatrixBase * rowCopyBase=model->rowCopy();
2709#ifndef NO_RTTI
2710    ClpPackedMatrix* rowCopy =
2711      dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
2712    // Make sure it is really a ClpPackedMatrix
2713    assert (rowCopy!=NULL);
2714#else
2715    ClpPackedMatrix* rowCopy =
2716      static_cast< ClpPackedMatrix*>(rowCopyBase);
2717#endif
2718
2719    const int * column = rowCopy->getIndices();
2720    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2721    const double * element = rowCopy->getElements();
2722    const double * rowScale = model->rowScale();
2723    const double * columnScale = model->columnScale();
2724    // scale row copy
2725    for (int iRow=0;iRow<numberRows;iRow++) {
2726      CoinBigIndex j;
2727      double scale = rowScale[iRow];
2728      const double * elementsInThisRow = element + rowStart[iRow];
2729      const int * columnsInThisRow = column + rowStart[iRow];
2730      int number = rowStart[iRow+1]-rowStart[iRow];
2731      assert (number<=numberColumns);
2732      for (j=0;j<number;j++) {
2733        int iColumn = columnsInThisRow[j];
2734        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
2735      }
2736      rowCopy->replaceVector(iRow,number,newElement);
2737    }
2738    delete [] newElement;
2739  }
2740}
2741/* Realy really scales column copy
2742   Only called if scales already exist.
2743   Up to user ro delete */
2744ClpMatrixBase * 
2745ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const 
2746{
2747  // need to replace column by column
2748  int numberRows = model->numberRows();
2749  int numberColumns = matrix_->getNumCols();
2750  double * newElement = new double[numberRows];
2751  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
2752  const int * row = copy->getIndices();
2753  const CoinBigIndex * columnStart = copy->getVectorStarts();
2754  const int * length = copy->getVectorLengths();
2755  const double * element = copy->getElements();
2756  const double * rowScale = model->rowScale();
2757  const double * columnScale = model->columnScale();
2758  // scale column copy
2759  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2760    CoinBigIndex j;
2761    double scale = columnScale[iColumn];
2762    const double * elementsInThisColumn = element + columnStart[iColumn];
2763    const int * rowsInThisColumn = row + columnStart[iColumn];
2764    int number = length[iColumn];
2765    assert (number<=numberRows);
2766    for (j=0;j<number;j++) {
2767      int iRow = rowsInThisColumn[j];
2768      newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
2769    }
2770    copy->replaceVector(iColumn,number,newElement);
2771  }
2772  delete [] newElement;
2773  return copy;
2774}
2775// Really scale matrix
2776void 
2777ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
2778{
2779  int numberColumns = matrix_->getNumCols();
2780  const int * row = matrix_->getIndices();
2781  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2782  const int * length = matrix_->getVectorLengths();
2783  double * element = matrix_->getMutableElements();
2784  // scale
2785  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2786    CoinBigIndex j;
2787    double scale = columnScale[iColumn];
2788    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
2789      int iRow = row[j];
2790      element[j] *= scale*rowScale[iRow];
2791    }
2792  }
2793}
2794/* Delete the columns whose indices are listed in <code>indDel</code>. */
2795void 
2796ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
2797{ 
2798  if (matrix_->getNumCols())
2799    matrix_->deleteCols(numDel,indDel);
2800  numberActiveColumns_ = matrix_->getNumCols();
2801  // may now have gaps
2802  hasGaps_=true;
2803  matrix_->setExtraGap(1.0e-50);
2804}
2805/* Delete the rows whose indices are listed in <code>indDel</code>. */
2806void 
2807ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
2808{ 
2809  if (matrix_->getNumRows())
2810    matrix_->deleteRows(numDel,indDel);
2811  numberActiveColumns_ = matrix_->getNumCols();
2812  // may now have gaps
2813  hasGaps_=true;
2814  matrix_->setExtraGap(1.0e-50);
2815}
2816#ifndef CLP_NO_VECTOR
2817// Append Columns
2818void 
2819ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
2820{ 
2821  matrix_->appendCols(number,columns);
2822  numberActiveColumns_ = matrix_->getNumCols();
2823}
2824// Append Rows
2825void 
2826ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
2827{ 
2828  matrix_->appendRows(number,rows);
2829  numberActiveColumns_ = matrix_->getNumCols();
2830  // may now have gaps
2831  hasGaps_=true;
2832}
2833#endif
2834/* Set the dimensions of the matrix. In effect, append new empty
2835   columns/rows to the matrix. A negative number for either dimension
2836   means that that dimension doesn't change. Otherwise the new dimensions
2837   MUST be at least as large as the current ones otherwise an exception
2838   is thrown. */
2839void 
2840ClpPackedMatrix::setDimensions(int numrows, int numcols)
2841{
2842  matrix_->setDimensions(numrows,numcols);
2843}
2844/* Append a set of rows/columns to the end of the matrix. Returns number of errors
2845   i.e. if any of the new rows/columns contain an index that's larger than the
2846   number of columns-1/rows-1 (if numberOther>0) or duplicates
2847   If 0 then rows, 1 if columns */
2848int 
2849ClpPackedMatrix::appendMatrix(int number, int type,
2850                            const CoinBigIndex * starts, const int * index,
2851                            const double * element, int numberOther)
2852{
2853  int numberErrors=0;
2854  // make sure other dimension is big enough
2855  if (type==0) {
2856    // rows
2857    if (matrix_->isColOrdered()&&numberOther>matrix_->getNumCols())
2858      matrix_->setDimensions(-1,numberOther);
2859    numberErrors=matrix_->appendRows(number,starts,index,element,numberOther);
2860  } else {
2861    // columns
2862    if (!matrix_->isColOrdered()&&numberOther>matrix_->getNumRows())
2863      matrix_->setDimensions(numberOther,-1);
2864    numberErrors=matrix_->appendCols(number,starts,index,element,numberOther);
2865  }
2866  numberActiveColumns_ = matrix_->getNumCols();
2867  return numberErrors;
2868}
2869void
2870ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
2871{
2872  delete rowCopy_;
2873  rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
2874  // See if anything in it
2875  if (!rowCopy_->usefulInfo()) {
2876    delete rowCopy_;
2877    rowCopy_=NULL;
2878  }
2879}
2880//#############################################################################
2881// Constructors / Destructor / Assignment
2882//#############################################################################
2883
2884//-------------------------------------------------------------------
2885// Default Constructor
2886//-------------------------------------------------------------------
2887ClpPackedMatrix2::ClpPackedMatrix2 () 
2888  : numberBlocks_(0),
2889    numberRows_(0),
2890    offset_(NULL),
2891    count_(NULL),
2892    rowStart_(NULL),
2893    column_(NULL),
2894    work_(NULL)
2895{
2896#ifdef THREAD
2897  threadId_ = NULL;
2898  info_ = NULL;
2899#endif
2900}
2901//-------------------------------------------------------------------
2902// Useful Constructor
2903//-------------------------------------------------------------------
2904ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy) 
2905  : numberBlocks_(0),
2906    numberRows_(0),
2907    offset_(NULL),
2908    count_(NULL),
2909    rowStart_(NULL),
2910    column_(NULL),
2911    work_(NULL)
2912{
2913#ifdef THREAD
2914  threadId_ = NULL;
2915  info_ = NULL;
2916#endif
2917  numberRows_ = rowCopy->getNumRows();
2918  if (!numberRows_)
2919    return;
2920  int numberColumns = rowCopy->getNumCols();
2921  const int * column = rowCopy->getIndices();
2922  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2923  const int * length = rowCopy->getVectorLengths();
2924  const double * element = rowCopy->getElements();
2925  int chunk = 32768; // tune
2926  //chunk=100;
2927  // tune
2928#if 0
2929  int chunkY[5]={4096,8192,16384,32768,65535};
2930  int its = model->maximumIterations();
2931  if (its>=1000000&&its<1000999) {
2932    its -= 1000000;
2933    its = its/10;
2934    if (its>=5) abort();
2935    chunk=chunkY[its];
2936    printf("chunk size %d\n",chunk);
2937#define cpuid(func,ax,bx,cx,dx)\
2938    __asm__ __volatile__ ("cpuid":\
2939    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
2940    unsigned int a,b,c,d;
2941    int func=0;
2942    cpuid(func,a,b,c,d);
2943    {
2944      int i;
2945      unsigned int value;
2946      value=b;
2947      for (i=0;i<4;i++) {
2948        printf("%c",(value&0xff));
2949        value = value >>8;
2950      }
2951      value=d;
2952      for (i=0;i<4;i++) {
2953        printf("%c",(value&0xff));
2954        value = value >>8;
2955      }
2956      value=c;
2957      for (i=0;i<4;i++) {
2958        printf("%c",(value&0xff));
2959        value = value >>8;
2960      }
2961      printf("\n");
2962      int maxfunc = a;
2963      if (maxfunc>10) {
2964        printf("not intel?\n");
2965        abort();
2966      }
2967      for (func=1;func<=maxfunc;func++) {
2968        cpuid(func,a,b,c,d);
2969        printf("func %d, %x %x %x %x\n",func,a,b,c,d);
2970      }
2971    }
2972#else
2973    if (numberColumns>10000||chunk==100) {
2974#endif
2975  } else {
2976    //printf("no chunk\n");
2977    return;
2978  }
2979  // Could also analyze matrix to get natural breaks
2980  numberBlocks_ = (numberColumns+chunk-1)/chunk;
2981#ifdef THREAD
2982  // Get work areas
2983  threadId_ = new pthread_t [numberBlocks_];
2984  info_ = new dualColumn0Struct[numberBlocks_];
2985#endif
2986  // Even out
2987  chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
2988  offset_ = new int[numberBlocks_+1];
2989  offset_[numberBlocks_]=numberColumns;
2990  int nRow = numberBlocks_*numberRows_;
2991  count_ = new unsigned short[nRow];
2992  memset(count_,0,nRow*sizeof(unsigned short));
2993  rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
2994  CoinBigIndex nElement = rowStart[numberRows_];
2995  rowStart_[nRow+numberRows_]=nElement;
2996  column_ = new unsigned short [nElement];
2997  // assumes int <= double
2998  int sizeWork = 6*numberBlocks_;
2999  work_ = new double[sizeWork];;
3000  int iBlock;
3001  int nZero=0;
3002  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3003    int start = iBlock*chunk;
3004    offset_[iBlock]=start;
3005    int end = start+chunk;
3006    for (int iRow=0;iRow<numberRows_;iRow++) {
3007      if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
3008        printf("not packed correctly - gaps\n");
3009        abort();
3010      }
3011      bool lastFound=false;
3012      int nFound=0;
3013      for (CoinBigIndex j = rowStart[iRow];
3014           j<rowStart[iRow]+length[iRow];j++) {
3015        int iColumn = column[j];
3016        if (iColumn>=start) {
3017          if (iColumn<end) {
3018            if (!element[j]) {
3019              printf("not packed correctly - zero element\n");
3020              abort();
3021            }
3022            column_[j]=iColumn-start;
3023            nFound++;
3024            if (lastFound) {
3025              printf("not packed correctly - out of order\n");
3026              abort();
3027            }
3028          } else {
3029            //can't find any more
3030            lastFound=true;
3031          }
3032        } 
3033      }
3034      count_[iRow*numberBlocks_+iBlock]=nFound;
3035      if (!nFound)
3036        nZero++;
3037    }
3038  }
3039  //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
3040  //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
3041}
3042
3043//-------------------------------------------------------------------
3044// Copy constructor
3045//-------------------------------------------------------------------
3046ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs) 
3047  : numberBlocks_(rhs.numberBlocks_),
3048    numberRows_(rhs.numberRows_)
3049{ 
3050  if (numberBlocks_) {
3051    offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3052    int nRow = numberBlocks_*numberRows_;
3053    count_ = CoinCopyOfArray(rhs.count_,nRow);
3054    rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3055    CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3056    column_ = CoinCopyOfArray(rhs.column_,nElement);
3057    int sizeWork = 6*numberBlocks_;
3058    work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3059#ifdef THREAD
3060    threadId_ = new pthread_t [numberBlocks_];
3061    info_ = new dualColumn0Struct[numberBlocks_];
3062#endif
3063  } else {
3064    offset_ = NULL;
3065    count_ = NULL;
3066    rowStart_ = NULL;
3067    column_ = NULL;
3068    work_ = NULL;
3069#ifdef THREAD
3070    threadId_ = NULL;
3071    info_ = NULL;
3072#endif
3073  }
3074}
3075//-------------------------------------------------------------------
3076// Destructor
3077//-------------------------------------------------------------------
3078ClpPackedMatrix2::~ClpPackedMatrix2 ()
3079{
3080  delete [] offset_;
3081  delete [] count_;
3082  delete [] rowStart_;
3083  delete [] column_;
3084  delete [] work_;
3085#ifdef THREAD
3086  delete [] threadId_;
3087  delete [] info_;
3088#endif
3089}
3090
3091//----------------------------------------------------------------
3092// Assignment operator
3093//-------------------------------------------------------------------
3094ClpPackedMatrix2 &
3095ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
3096{
3097  if (this != &rhs) {
3098    numberBlocks_ = rhs.numberBlocks_;
3099    numberRows_ = rhs.numberRows_;
3100    delete [] offset_;
3101    delete [] count_;
3102    delete [] rowStart_;
3103    delete [] column_;
3104    delete [] work_;
3105#ifdef THREAD
3106    delete [] threadId_;
3107    delete [] info_;
3108#endif
3109    if (numberBlocks_) {
3110      offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3111      int nRow = numberBlocks_*numberRows_;
3112      count_ = CoinCopyOfArray(rhs.count_,nRow);
3113      rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3114      CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3115      column_ = CoinCopyOfArray(rhs.column_,nElement);
3116      int sizeWork = 6*numberBlocks_;
3117      work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3118#ifdef THREAD
3119      threadId_ = new pthread_t [numberBlocks_];
3120      info_ = new dualColumn0Struct[numberBlocks_];
3121#endif
3122    } else {
3123      offset_ = NULL;
3124      count_ = NULL;
3125      rowStart_ = NULL;
3126      column_ = NULL;
3127      work_ = NULL;
3128#ifdef THREAD
3129      threadId_ = NULL;
3130      info_ = NULL;
3131#endif
3132    }
3133  }
3134  return *this;
3135}
3136static int dualColumn0(const ClpSimplex * model,double * spare,
3137                       int * spareIndex,const double * arrayTemp,
3138                       const int * indexTemp,int numberIn,
3139                       int offset, double acceptablePivot,double * bestPossiblePtr,
3140                       double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
3141{
3142  // do dualColumn0
3143  int i;
3144  int numberRemaining=0;
3145  double bestPossible=0.0;
3146  double upperTheta=1.0e31;
3147  double freePivot = acceptablePivot;
3148  int posFree=-1;
3149  const double * reducedCost=model->djRegion(1);
3150  double dualTolerance = model->dualTolerance();
3151  // We can also see if infeasible or pivoting on free
3152  double tentativeTheta = 1.0e25;
3153  for (i=0;i<numberIn;i++) {
3154    double alpha = arrayTemp[i];
3155    int iSequence = indexTemp[i]+offset;
3156    double oldValue;
3157    double value;
3158    bool keep;
3159   
3160    switch(model->getStatus(iSequence)) {
3161     
3162    case ClpSimplex::basic:
3163    case ClpSimplex::isFixed:
3164      break;
3165    case ClpSimplex::isFree:
3166    case ClpSimplex::superBasic:
3167      bestPossible = CoinMax(bestPossible,fabs(alpha));
3168      oldValue = reducedCost[iSequence];
3169      // If free has to be very large - should come in via dualRow
3170      if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3171        break;
3172      if (oldValue>dualTolerance) {
3173        keep = true;
3174      } else if (oldValue<-dualTolerance) {
3175        keep = true;
3176      } else {
3177        if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3178          keep = true;
3179        else
3180          keep=false;
3181      }
3182      if (keep) {
3183        // free - choose largest
3184        if (fabs(alpha)>freePivot) {
3185          freePivot=fabs(alpha);
3186          posFree = i;
3187        }
3188      }
3189      break;
3190    case ClpSimplex::atUpperBound:
3191      oldValue = reducedCost[iSequence];
3192      value = oldValue-tentativeTheta*alpha;
3193      //assert (oldValue<=dualTolerance*1.0001);
3194      if (value>dualTolerance) {
3195        bestPossible = CoinMax(bestPossible,-alpha);
3196        value = oldValue-upperTheta*alpha;
3197        if (value>dualTolerance && -alpha>=acceptablePivot) 
3198          upperTheta = (oldValue-dualTolerance)/alpha;
3199        // add to list
3200        spare[numberRemaining]=alpha;
3201        spareIndex[numberRemaining++]=iSequence;
3202      }
3203      break;
3204    case ClpSimplex::atLowerBound:
3205      oldValue = reducedCost[iSequence];
3206      value = oldValue-tentativeTheta*alpha;
3207      //assert (oldValue>=-dualTolerance*1.0001);
3208      if (value<-dualTolerance) {
3209        bestPossible = CoinMax(bestPossible,alpha);
3210        value = oldValue-upperTheta*alpha;
3211        if (value<-dualTolerance && alpha>=acceptablePivot) 
3212          upperTheta = (oldValue+dualTolerance)/alpha;
3213        // add to list
3214        spare[numberRemaining]=alpha;
3215        spareIndex[numberRemaining++]=iSequence;
3216      }
3217      break;
3218    }
3219  }
3220  *bestPossiblePtr = bestPossible;
3221  *upperThetaPtr = upperTheta;
3222  *freePivotPtr = freePivot;
3223  *posFreePtr = posFree;
3224  return numberRemaining;
3225}
3226static int doOneBlock(double * array, int * index, 
3227                      const double * pi,const CoinBigIndex * rowStart, const double * element,
3228                      const unsigned short * column, int numberInRowArray, int numberLook)
3229{
3230  int iWhich=0;
3231  int nextN=0;
3232  CoinBigIndex nextStart=0;
3233  double nextPi=0.0;
3234  for (;iWhich<numberInRowArray;iWhich++) {
3235    nextStart = rowStart[0];
3236    nextN = rowStart[numberInRowArray]-nextStart;
3237    rowStart++;
3238    if (nextN) {
3239      nextPi = pi[iWhich];
3240      break;
3241    }
3242  }
3243  int i;
3244  while (iWhich<numberInRowArray) {
3245    double value = nextPi;
3246    CoinBigIndex j=  nextStart;
3247    int n = nextN;
3248    // get next
3249    iWhich++;
3250    for (;iWhich<numberInRowArray;iWhich++) {
3251      nextStart = rowStart[0];
3252      nextN = rowStart[numberInRowArray]-nextStart;
3253      rowStart++;
3254      if (nextN) {
3255        coin_prefetch(element+nextStart);
3256        nextPi = pi[iWhich];
3257        break;
3258      }
3259    }
3260    CoinBigIndex end =j+n;
3261    //coin_prefetch(element+rowStart_[i+1]);
3262    //coin_prefetch(column_+rowStart_[i+1]);
3263    if (n<100) {
3264      if ((n&1)!=0) {
3265        unsigned int jColumn = column[j];
3266        array[jColumn] -= value*element[j];
3267        j++;
3268      }     
3269      coin_prefetch(column+nextStart);
3270      for (;j<end;j+=2) {
3271        unsigned int jColumn0 = column[j];
3272        double value0 = value*element[j];
3273        unsigned int jColumn1 = column[j+1];
3274        double value1 = value*element[j+1];
3275        array[jColumn0] -= value0;
3276        array[jColumn1] -= value1;
3277      }
3278    } else {
3279      if ((n&1)!=0) {
3280        unsigned int jColumn = column[j];
3281        array[jColumn] -= value*element[j];
3282        j++;
3283      }     
3284      if ((n&2)!=0) {
3285        unsigned int jColumn0 = column[j];
3286        double value0 = value*element[j];
3287        unsigned int jColumn1 = column[j+1];
3288        double value1 = value*element[j+1];
3289        array[jColumn0] -= value0;
3290        array[jColumn1] -= value1;
3291        j+=2;
3292      }     
3293      if ((n&4)!=0) {
3294        unsigned int jColumn0 = column[j];
3295        double value0 = value*element[j];
3296        unsigned int jColumn1 = column[j+1];
3297        double value1 = value*element[j+1];
3298        unsigned int jColumn2 = column[j+2];
3299        double value2 = value*element[j+2];
3300        unsigned int jColumn3 = column[j+3];
3301        double value3 = value*element[j+3];
3302        array[jColumn0] -= value0;
3303        array[jColumn1] -= value1;
3304        array[jColumn2] -= value2;
3305        array[jColumn3] -= value3;
3306        j+=4;
3307      }
3308      //coin_prefetch(column+nextStart);
3309      for (;j<end;j+=8) {
3310        coin_prefetch(element+j+16);
3311        unsigned int jColumn0 = column[j];
3312        double value0 = value*element[j];
3313        unsigned int jColumn1 = column[j+1];
3314        double value1 = value*element[j+1];
3315        unsigned int jColumn2 = column[j+2];
3316        double value2 = value*element[j+2];
3317        unsigned int jColumn3 = column[j+3];
3318        double value3 = value*element[j+3];
3319        array[jColumn0] -= value0;
3320        array[jColumn1] -= value1;
3321        array[jColumn2] -= value2;
3322        array[jColumn3] -= value3;
3323        coin_prefetch(column+j+16);
3324        jColumn0 = column[j+4];
3325        value0 = value*element[j+4];
3326        jColumn1 = column[j+5];
3327        value1 = value*element[j+5];
3328        jColumn2 = column[j+6];
3329        value2 = value*element[j+6];
3330        jColumn3 = column[j+7];
3331        value3 = value*element[j+7];
3332        array[jColumn0] -= value0;
3333        array[jColumn1] -= value1;
3334        array[jColumn2] -= value2;
3335        array[jColumn3] -= value3;
3336      }
3337    }
3338  }
3339  // get rid of tiny values
3340  int nSmall = numberLook;
3341  int numberNonZero=0;
3342  for (i=0;i<nSmall;i++) {
3343    double value = array[i];
3344    array[i]=0.0; 
3345    if (fabs(value)>1.0e-12) {
3346      array[numberNonZero]=value;
3347      index[numberNonZero++]=i;
3348    }
3349  }
3350  for (;i<numberLook;i+=4) {
3351    double value0 = array[i+0];
3352    double value1 = array[i+1];
3353    double value2 = array[i+2];
3354    double value3 = array[i+3];
3355    array[i+0]=0.0; 
3356    array[i+1]=0.0; 
3357    array[i+2]=0.0; 
3358    array[i+3]=0.0; 
3359    if (fabs(value0)>1.0e-12) {
3360      array[numberNonZero]=value0;
3361      index[numberNonZero++]=i+0;
3362    }
3363    if (fabs(value1)>1.0e-12) {
3364      array[numberNonZero]=value1;
3365      index[numberNonZero++]=i+1;
3366    }
3367    if (fabs(value2)>1.0e-12) {
3368      array[numberNonZero]=value2;
3369      index[numberNonZero++]=i+2;
3370    }
3371    if (fabs(value3)>1.0e-12) {
3372      array[numberNonZero]=value3;
3373      index[numberNonZero++]=i+3;
3374    }
3375  }
3376  return numberNonZero;
3377}
3378#ifdef THREAD
3379static void * doOneBlockThread(void * voidInfo)
3380{
3381  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3382  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3383                                  info->rowStart,info->element,info->column,
3384                                  info->numberInRowArray,info->numberLook);
3385  return NULL;
3386}
3387static void * doOneBlockAnd0Thread(void * voidInfo)
3388{
3389  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3390  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3391                                  info->rowStart,info->element,info->column,
3392                                  info->numberInRowArray,info->numberLook);
3393  *(info->numberOutPtr) =  dualColumn0(info->model,info->spare,
3394                                    info->spareIndex,(const double *)info->arrayTemp,
3395                                    (const int *) info->indexTemp,*(info->numberInPtr),
3396                                    info->offset, info->acceptablePivot,info->bestPossiblePtr,
3397                                    info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
3398  return NULL;
3399}
3400#endif
3401/* Return <code>x * scalar * A in <code>z</code>.
3402   Note - x packed and z will be packed mode
3403   Squashes small elements and knows about ClpSimplex */
3404void 
3405ClpPackedMatrix2::transposeTimes(const ClpSimplex * model, 
3406                                 const CoinPackedMatrix * rowCopy,
3407                                 const CoinIndexedVector * rowArray,
3408                                 CoinIndexedVector * spareArray,
3409                                 CoinIndexedVector * columnArray) const
3410{
3411  // See if dualColumn0 coding wanted
3412  bool dualColumn = model->spareIntArray_[0]==1;
3413  double acceptablePivot = model->spareDoubleArray_[0];
3414  double bestPossible=0.0;
3415  double upperTheta=1.0e31;
3416  double freePivot = acceptablePivot;
3417  int posFree=-1;
3418  int numberRemaining=0;
3419  //if (model->numberIterations()>=200000) {
3420  //printf("time %g\n",CoinCpuTime()-startTime);
3421  //exit(77);
3422  //}
3423  double * pi = rowArray->denseVector();
3424  int numberNonZero=0;
3425  int * index = columnArray->getIndices();
3426  double * array = columnArray->denseVector();
3427  int numberInRowArray = rowArray->getNumElements();
3428  const int * whichRow = rowArray->getIndices();
3429  double * element = const_cast<double *>(rowCopy->getElements());
3430  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3431  int i;
3432  CoinBigIndex * rowStart2 = rowStart_;
3433  if (!dualColumn) {
3434    for (i=0;i<numberInRowArray;i++) {
3435      int iRow = whichRow[i]; 
3436      CoinBigIndex start = rowStart[iRow];
3437      *rowStart2 = start;
3438      unsigned short * count1 = count_+iRow*numberBlocks_;
3439      int put=0;
3440      for (int j=0;j<numberBlocks_;j++) {
3441        put += numberInRowArray;
3442        start += count1[j];
3443        rowStart2[put]=start;
3444      }
3445      rowStart2 ++;
3446    }
3447  } else {
3448    // also do dualColumn stuff
3449    double * spare = spareArray->denseVector();
3450    int * spareIndex = spareArray->getIndices();
3451    const double * reducedCost=model->djRegion(0);
3452    double dualTolerance = model->dualTolerance();
3453    // We can also see if infeasible or pivoting on free
3454    double tentativeTheta = 1.0e25;
3455    int addSequence=model->numberColumns();
3456    for (i=0;i<numberInRowArray;i++) {
3457      int iRow = whichRow[i]; 
3458      double alpha=pi[i];
3459      double oldValue;
3460      double value;
3461      bool keep;
3462
3463      switch(model->getStatus(iRow+addSequence)) {
3464         
3465      case ClpSimplex::basic:
3466      case ClpSimplex::isFixed:
3467        break;
3468      case ClpSimplex::isFree:
3469      case ClpSimplex::superBasic:
3470        bestPossible = CoinMax(bestPossible,fabs(alpha));
3471        oldValue = reducedCost[iRow];
3472        // If free has to be very large - should come in via dualRow
3473        if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3474          break;
3475        if (oldValue>dualTolerance) {
3476          keep = true;
3477        } else if (oldValue<-dualTolerance) {
3478          keep = true;
3479        } else {
3480          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3481            keep = true;
3482          else
3483            keep=false;
3484        }
3485        if (keep) {
3486          // free - choose largest
3487          if (fabs(alpha)>freePivot) {
3488            freePivot=fabs(alpha);
3489            posFree = i + addSequence;
3490          }
3491        }
3492        break;
3493      case ClpSimplex::atUpperBound:
3494        oldValue = reducedCost[iRow];
3495        value = oldValue-tentativeTheta*alpha;
3496        //assert (oldValue<=dualTolerance*1.0001);
3497        if (value>dualTolerance) {
3498          bestPossible = CoinMax(bestPossible,-alpha);
3499          value = oldValue-upperTheta*alpha;
3500          if (value>dualTolerance && -alpha>=acceptablePivot) 
3501            upperTheta = (oldValue-dualTolerance)/alpha;
3502          // add to list
3503          spare[numberRemaining]=alpha;
3504          spareIndex[numberRemaining++]=iRow+addSequence;
3505        }
3506        break;
3507      case ClpSimplex::atLowerBound:
3508        oldValue = reducedCost[iRow];
3509        value = oldValue-tentativeTheta*alpha;
3510        //assert (oldValue>=-dualTolerance*1.0001);
3511        if (value<-dualTolerance) {
3512          bestPossible = CoinMax(bestPossible,alpha);
3513          value = oldValue-upperTheta*alpha;
3514          if (value<-dualTolerance && alpha>=acceptablePivot) 
3515            upperTheta = (oldValue+dualTolerance)/alpha;
3516          // add to list
3517          spare[numberRemaining]=alpha;
3518          spareIndex[numberRemaining++]=iRow+addSequence;
3519        }
3520        break;
3521      }
3522      CoinBigIndex start = rowStart[iRow];
3523      *rowStart2 = start;
3524      unsigned short * count1 = count_+iRow*numberBlocks_;
3525      int put=0;
3526      for (int j=0;j<numberBlocks_;j++) {
3527        put += numberInRowArray;
3528        start += count1[j];
3529        rowStart2[put]=start;
3530      }
3531      rowStart2 ++;
3532    }
3533  }
3534  double * spare = spareArray->denseVector();
3535  int * spareIndex = spareArray->getIndices();
3536  int saveNumberRemaining=numberRemaining;
3537  int iBlock;
3538  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3539    double * dwork = work_+6*iBlock;
3540    int * iwork = (int *) (dwork+3);
3541    if (!dualColumn) {
3542#ifndef THREAD
3543      int offset = offset_[iBlock];
3544      int offset3 = offset;
3545      offset = numberNonZero;
3546      double * arrayTemp = array+offset;
3547      int * indexTemp = index+offset; 
3548      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
3549                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
3550      int number = iwork[0];
3551      for (i=0;i<number;i++) {
3552        //double value = arrayTemp[i];
3553        //arrayTemp[i]=0.0;
3554        //array[numberNonZero]=value;
3555        index[numberNonZero++]=indexTemp[i]+offset3;
3556      }
3557#else
3558      int offset = offset_[iBlock];
3559      double * arrayTemp = array+offset;
3560      int * indexTemp = index+offset; 
3561      dualColumn0Struct * infoPtr = info_+iBlock;
3562      infoPtr->arrayTemp=arrayTemp;
3563      infoPtr->indexTemp=indexTemp;
3564      infoPtr->numberInPtr=&iwork[0];
3565      infoPtr->pi=pi;
3566      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
3567      infoPtr->element=element;
3568      infoPtr->column=column_;
3569      infoPtr->numberInRowArray=numberInRowArray;
3570      infoPtr->numberLook=offset_[iBlock+1]-offset;
3571      pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
3572#endif
3573    } else {
3574#ifndef THREAD
3575      int offset = offset_[iBlock];
3576      // allow for already saved
3577      int offset2 = offset + saveNumberRemaining;
3578      int offset3 = offset;
3579      offset = numberNonZero;
3580      offset2 = numberRemaining;
3581      double * arrayTemp = array+offset;
3582      int * indexTemp = index+offset; 
3583      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
3584                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
3585      iwork[1] = dualColumn0(model,spare+offset2,
3586                                   spareIndex+offset2,
3587                                   arrayTemp,indexTemp,
3588                                   iwork[0],offset3,acceptablePivot,
3589                                   &dwork[0],&dwork[1],&iwork[2],
3590                                   &dwork[2]);
3591      int number = iwork[0];
3592      int numberLook = iwork[1];
3593#if 1
3594      numberRemaining += numberLook;
3595#else
3596      double * spareTemp = spare+offset2;
3597      const int * spareIndexTemp = spareIndex+offset2; 
3598      for (i=0;i<numberLook;i++) {
3599        double value = spareTemp[i];
3600        spareTemp[i]=0.0;
3601        spare[numberRemaining] = value;
3602        spareIndex[numberRemaining++] = spareIndexTemp[i];
3603      }
3604#endif
3605      if (dwork[2]>freePivot) {
3606        freePivot=dwork[2];
3607        posFree = iwork[2]+numberNonZero;
3608      } 
3609      upperTheta =  CoinMin(dwork[1],upperTheta);
3610      bestPossible = CoinMax(dwork[0],bestPossible);
3611      for (i=0;i<number;i++) {
3612        // double value = arrayTemp[i];
3613        //arrayTemp[i]=0.0;
3614        //array[numberNonZero]=value;
3615        index[numberNonZero++]=indexTemp[i]+offset3;
3616      }
3617#else
3618      int offset = offset_[iBlock];
3619      // allow for already saved
3620      int offset2 = offset + saveNumberRemaining;
3621      double * arrayTemp = array+offset;
3622      int * indexTemp = index+offset; 
3623      dualColumn0Struct * infoPtr = info_+iBlock;
3624      infoPtr->model=model;
3625      infoPtr->spare=spare+offset2;
3626      infoPtr->spareIndex=spareIndex+offset2;
3627      infoPtr->arrayTemp=arrayTemp;
3628      infoPtr->indexTemp=indexTemp;
3629      infoPtr->numberInPtr=&iwork[0];
3630      infoPtr->offset=offset;
3631      infoPtr->acceptablePivot=acceptablePivot;
3632      infoPtr->bestPossiblePtr=&dwork[0];
3633      infoPtr->upperThetaPtr=&dwork[1];
3634      infoPtr->posFreePtr=&iwork[2];
3635      infoPtr->freePivotPtr=&dwork[2];
3636      infoPtr->numberOutPtr=&iwork[1];
3637      infoPtr->pi=pi;
3638      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
3639      infoPtr->element=element;
3640      infoPtr->column=column_;
3641      infoPtr->numberInRowArray=numberInRowArray;
3642      infoPtr->numberLook=offset_[iBlock+1]-offset;
3643      if (iBlock>=2)
3644        pthread_join(threadId_[iBlock-2],NULL);
3645      pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
3646      //pthread_join(threadId_[iBlock],NULL);
3647#endif
3648    }
3649  }
3650  for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
3651#ifdef THREAD
3652    pthread_join(threadId_[iBlock],NULL);
3653#endif
3654  }
3655#ifdef THREAD
3656  for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
3657    //pthread_join(threadId_[iBlock],NULL);
3658    int offset = offset_[iBlock];
3659    double * dwork = work_+6*iBlock;
3660    int * iwork = (int *) (dwork+3);
3661    int number = iwork[0];
3662    if (dualColumn) {
3663      // allow for already saved
3664      int offset2 = offset + saveNumberRemaining;
3665      int numberLook = iwork[1];
3666      double * spareTemp = spare+offset2;
3667      const int * spareIndexTemp = spareIndex+offset2; 
3668      for (i=0;i<numberLook;i++) {
3669        double value = spareTemp[i];
3670        spareTemp[i]=0.0;
3671        spare[numberRemaining] = value;
3672        spareIndex[numberRemaining++] = spareIndexTemp[i];
3673      }
3674      if (dwork[2]>freePivot) {
3675        freePivot=dwork[2];
3676        posFree = iwork[2]+numberNonZero;
3677      } 
3678      upperTheta =  CoinMin(dwork[1],upperTheta);
3679      bestPossible = CoinMax(dwork[0],bestPossible);
3680    }
3681    double * arrayTemp = array+offset;
3682    const int * indexTemp = index+offset; 
3683    for (i=0;i<number;i++) {
3684      double value = arrayTemp[i];
3685      arrayTemp[i]=0.0; 
3686      array[numberNonZero]=value;
3687      index[numberNonZero++]=indexTemp[i]+offset;
3688    }
3689  }
3690#endif
3691  columnArray->setNumElements(numberNonZero);
3692  columnArray->setPackedMode(true);
3693  if (dualColumn) {
3694    model->spareDoubleArray_[0]=upperTheta;
3695    model->spareDoubleArray_[1]=bestPossible;
3696    // and theta and alpha and sequence
3697    if (posFree<0) {
3698      model->spareIntArray_[1]=-1;
3699    } else {
3700      const double * reducedCost=model->djRegion(0);
3701      double alpha;
3702      int numberColumns=model->numberColumns();
3703      if (posFree<numberColumns) {
3704        alpha = columnArray->denseVector()[posFree];
3705        posFree = columnArray->getIndices()[posFree];
3706      } else {
3707        alpha = rowArray->denseVector()[posFree-numberColumns];
3708        posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
3709      }
3710      model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
3711      model->spareDoubleArray_[3]=alpha;
3712      model->spareIntArray_[1]=posFree;
3713    }
3714    spareArray->setNumElements(numberRemaining);
3715    // signal done
3716    model->spareIntArray_[0]=-1;
3717  }
3718}
3719#ifdef CLP_ALL_ONE_FILE
3720#undef reference
3721#endif
Note: See TracBrowser for help on using the repository browser.