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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 114.4 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) throw(CoinError)
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  return numberErrors;
2867}
2868void
2869ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
2870{
2871  delete rowCopy_;
2872  rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
2873  // See if anything in it
2874  if (!rowCopy_->usefulInfo()) {
2875    delete rowCopy_;
2876    rowCopy_=NULL;
2877  }
2878}
2879//#############################################################################
2880// Constructors / Destructor / Assignment
2881//#############################################################################
2882
2883//-------------------------------------------------------------------
2884// Default Constructor
2885//-------------------------------------------------------------------
2886ClpPackedMatrix2::ClpPackedMatrix2 () 
2887  : numberBlocks_(0),
2888    numberRows_(0),
2889    offset_(NULL),
2890    count_(NULL),
2891    rowStart_(NULL),
2892    column_(NULL),
2893    work_(NULL)
2894{
2895#ifdef THREAD
2896  threadId_ = NULL;
2897  info_ = NULL;
2898#endif
2899}
2900//-------------------------------------------------------------------
2901// Useful Constructor
2902//-------------------------------------------------------------------
2903ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy) 
2904  : numberBlocks_(0),
2905    numberRows_(0),
2906    offset_(NULL),
2907    count_(NULL),
2908    rowStart_(NULL),
2909    column_(NULL),
2910    work_(NULL)
2911{
2912#ifdef THREAD
2913  threadId_ = NULL;
2914  info_ = NULL;
2915#endif
2916  numberRows_ = rowCopy->getNumRows();
2917  if (!numberRows_)
2918    return;
2919  int numberColumns = rowCopy->getNumCols();
2920  const int * column = rowCopy->getIndices();
2921  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2922  const int * length = rowCopy->getVectorLengths();
2923  const double * element = rowCopy->getElements();
2924  int chunk = 32768; // tune
2925  //chunk=100;
2926  // tune
2927#if 0
2928  int chunkY[5]={4096,8192,16384,32768,65535};
2929  int its = model->maximumIterations();
2930  if (its>=1000000&&its<1000999) {
2931    its -= 1000000;
2932    its = its/10;
2933    if (its>=5) abort();
2934    chunk=chunkY[its];
2935    printf("chunk size %d\n",chunk);
2936#define cpuid(func,ax,bx,cx,dx)\
2937    __asm__ __volatile__ ("cpuid":\
2938    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
2939    unsigned int a,b,c,d;
2940    int func=0;
2941    cpuid(func,a,b,c,d);
2942    {
2943      int i;
2944      unsigned int value;
2945      value=b;
2946      for (i=0;i<4;i++) {
2947        printf("%c",(value&0xff));
2948        value = value >>8;
2949      }
2950      value=d;
2951      for (i=0;i<4;i++) {
2952        printf("%c",(value&0xff));
2953        value = value >>8;
2954      }
2955      value=c;
2956      for (i=0;i<4;i++) {
2957        printf("%c",(value&0xff));
2958        value = value >>8;
2959      }
2960      printf("\n");
2961      int maxfunc = a;
2962      if (maxfunc>10) {
2963        printf("not intel?\n");
2964        abort();
2965      }
2966      for (func=1;func<=maxfunc;func++) {
2967        cpuid(func,a,b,c,d);
2968        printf("func %d, %x %x %x %x\n",func,a,b,c,d);
2969      }
2970    }
2971#else
2972    if (numberColumns>10000||chunk==100) {
2973#endif
2974  } else {
2975    //printf("no chunk\n");
2976    return;
2977  }
2978  // Could also analyze matrix to get natural breaks
2979  numberBlocks_ = (numberColumns+chunk-1)/chunk;
2980#ifdef THREAD
2981  // Get work areas
2982  threadId_ = new pthread_t [numberBlocks_];
2983  info_ = new dualColumn0Struct[numberBlocks_];
2984#endif
2985  // Even out
2986  chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
2987  offset_ = new int[numberBlocks_+1];
2988  offset_[numberBlocks_]=numberColumns;
2989  int nRow = numberBlocks_*numberRows_;
2990  count_ = new unsigned short[nRow];
2991  memset(count_,0,nRow*sizeof(unsigned short));
2992  rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
2993  CoinBigIndex nElement = rowStart[numberRows_];
2994  rowStart_[nRow+numberRows_]=nElement;
2995  column_ = new unsigned short [nElement];
2996  // assumes int <= double
2997  int sizeWork = 6*numberBlocks_;
2998  work_ = new double[sizeWork];;
2999  int iBlock;
3000  int nZero=0;
3001  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3002    int start = iBlock*chunk;
3003    offset_[iBlock]=start;
3004    int end = start+chunk;
3005    for (int iRow=0;iRow<numberRows_;iRow++) {
3006      if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
3007        printf("not packed correctly - gaps\n");
3008        abort();
3009      }
3010      bool lastFound=false;
3011      int nFound=0;
3012      for (CoinBigIndex j = rowStart[iRow];
3013           j<rowStart[iRow]+length[iRow];j++) {
3014        int iColumn = column[j];
3015        if (iColumn>=start) {
3016          if (iColumn<end) {
3017            if (!element[j]) {
3018              printf("not packed correctly - zero element\n");
3019              abort();
3020            }
3021            column_[j]=iColumn-start;
3022            nFound++;
3023            if (lastFound) {
3024              printf("not packed correctly - out of order\n");
3025              abort();
3026            }
3027          } else {
3028            //can't find any more
3029            lastFound=true;
3030          }
3031        } 
3032      }
3033      count_[iRow*numberBlocks_+iBlock]=nFound;
3034      if (!nFound)
3035        nZero++;
3036    }
3037  }
3038  //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
3039  //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
3040}
3041
3042//-------------------------------------------------------------------
3043// Copy constructor
3044//-------------------------------------------------------------------
3045ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs) 
3046  : numberBlocks_(rhs.numberBlocks_),
3047    numberRows_(rhs.numberRows_)
3048{ 
3049  if (numberBlocks_) {
3050    offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3051    int nRow = numberBlocks_*numberRows_;
3052    count_ = CoinCopyOfArray(rhs.count_,nRow);
3053    rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3054    CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3055    column_ = CoinCopyOfArray(rhs.column_,nElement);
3056    int sizeWork = 6*numberBlocks_;
3057    work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3058#ifdef THREAD
3059    threadId_ = new pthread_t [numberBlocks_];
3060    info_ = new dualColumn0Struct[numberBlocks_];
3061#endif
3062  } else {
3063    offset_ = NULL;
3064    count_ = NULL;
3065    rowStart_ = NULL;
3066    column_ = NULL;
3067    work_ = NULL;
3068#ifdef THREAD
3069    threadId_ = NULL;
3070    info_ = NULL;
3071#endif
3072  }
3073}
3074//-------------------------------------------------------------------
3075// Destructor
3076//-------------------------------------------------------------------
3077ClpPackedMatrix2::~ClpPackedMatrix2 ()
3078{
3079  delete [] offset_;
3080  delete [] count_;
3081  delete [] rowStart_;
3082  delete [] column_;
3083  delete [] work_;
3084#ifdef THREAD
3085  delete [] threadId_;
3086  delete [] info_;
3087#endif
3088}
3089
3090//----------------------------------------------------------------
3091// Assignment operator
3092//-------------------------------------------------------------------
3093ClpPackedMatrix2 &
3094ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
3095{
3096  if (this != &rhs) {
3097    numberBlocks_ = rhs.numberBlocks_;
3098    numberRows_ = rhs.numberRows_;
3099    delete [] offset_;
3100    delete [] count_;
3101    delete [] rowStart_;
3102    delete [] column_;
3103    delete [] work_;
3104#ifdef THREAD
3105    delete [] threadId_;
3106    delete [] info_;
3107#endif
3108    if (numberBlocks_) {
3109      offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3110      int nRow = numberBlocks_*numberRows_;
3111      count_ = CoinCopyOfArray(rhs.count_,nRow);
3112      rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3113      CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3114      column_ = CoinCopyOfArray(rhs.column_,nElement);
3115      int sizeWork = 6*numberBlocks_;
3116      work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3117#ifdef THREAD
3118      threadId_ = new pthread_t [numberBlocks_];
3119      info_ = new dualColumn0Struct[numberBlocks_];
3120#endif
3121    } else {
3122      offset_ = NULL;
3123      count_ = NULL;
3124      rowStart_ = NULL;
3125      column_ = NULL;
3126      work_ = NULL;
3127#ifdef THREAD
3128      threadId_ = NULL;
3129      info_ = NULL;
3130#endif
3131    }
3132  }
3133  return *this;
3134}
3135static int dualColumn0(const ClpSimplex * model,double * spare,
3136                       int * spareIndex,const double * arrayTemp,
3137                       const int * indexTemp,int numberIn,
3138                       int offset, double acceptablePivot,double * bestPossiblePtr,
3139                       double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
3140{
3141  // do dualColumn0
3142  int i;
3143  int numberRemaining=0;
3144  double bestPossible=0.0;
3145  double upperTheta=1.0e31;
3146  double freePivot = acceptablePivot;
3147  int posFree=-1;
3148  const double * reducedCost=model->djRegion(1);
3149  double dualTolerance = model->dualTolerance();
3150  // We can also see if infeasible or pivoting on free
3151  double tentativeTheta = 1.0e25;
3152  for (i=0;i<numberIn;i++) {
3153    double alpha = arrayTemp[i];
3154    int iSequence = indexTemp[i]+offset;
3155    double oldValue;
3156    double value;
3157    bool keep;
3158   
3159    switch(model->getStatus(iSequence)) {
3160     
3161    case ClpSimplex::basic:
3162    case ClpSimplex::isFixed:
3163      break;
3164    case ClpSimplex::isFree:
3165    case ClpSimplex::superBasic:
3166      bestPossible = CoinMax(bestPossible,fabs(alpha));
3167      oldValue = reducedCost[iSequence];
3168      // If free has to be very large - should come in via dualRow
3169      if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3170        break;
3171      if (oldValue>dualTolerance) {
3172        keep = true;
3173      } else if (oldValue<-dualTolerance) {
3174        keep = true;
3175      } else {
3176        if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3177          keep = true;
3178        else
3179          keep=false;
3180      }
3181      if (keep) {
3182        // free - choose largest
3183        if (fabs(alpha)>freePivot) {
3184          freePivot=fabs(alpha);
3185          posFree = i;
3186        }
3187      }
3188      break;
3189    case ClpSimplex::atUpperBound:
3190      oldValue = reducedCost[iSequence];
3191      value = oldValue-tentativeTheta*alpha;
3192      //assert (oldValue<=dualTolerance*1.0001);
3193      if (value>dualTolerance) {
3194        bestPossible = CoinMax(bestPossible,-alpha);
3195        value = oldValue-upperTheta*alpha;
3196        if (value>dualTolerance && -alpha>=acceptablePivot) 
3197          upperTheta = (oldValue-dualTolerance)/alpha;
3198        // add to list
3199        spare[numberRemaining]=alpha;
3200        spareIndex[numberRemaining++]=iSequence;
3201      }
3202      break;
3203    case ClpSimplex::atLowerBound:
3204      oldValue = reducedCost[iSequence];
3205      value = oldValue-tentativeTheta*alpha;
3206      //assert (oldValue>=-dualTolerance*1.0001);
3207      if (value<-dualTolerance) {
3208        bestPossible = CoinMax(bestPossible,alpha);
3209        value = oldValue-upperTheta*alpha;
3210        if (value<-dualTolerance && alpha>=acceptablePivot) 
3211          upperTheta = (oldValue+dualTolerance)/alpha;
3212        // add to list
3213        spare[numberRemaining]=alpha;
3214        spareIndex[numberRemaining++]=iSequence;
3215      }
3216      break;
3217    }
3218  }
3219  *bestPossiblePtr = bestPossible;
3220  *upperThetaPtr = upperTheta;
3221  *freePivotPtr = freePivot;
3222  *posFreePtr = posFree;
3223  return numberRemaining;
3224}
3225static int doOneBlock(double * array, int * index, 
3226                      const double * pi,const CoinBigIndex * rowStart, const double * element,
3227                      const unsigned short * column, int numberInRowArray, int numberLook)
3228{
3229  int iWhich=0;
3230  int nextN=0;
3231  CoinBigIndex nextStart=0;
3232  double nextPi=0.0;
3233  for (;iWhich<numberInRowArray;iWhich++) {
3234    nextStart = rowStart[0];
3235    nextN = rowStart[numberInRowArray]-nextStart;
3236    rowStart++;
3237    if (nextN) {
3238      nextPi = pi[iWhich];
3239      break;
3240    }
3241  }
3242  int i;
3243  while (iWhich<numberInRowArray) {
3244    double value = nextPi;
3245    CoinBigIndex j=  nextStart;
3246    int n = nextN;
3247    // get next
3248    iWhich++;
3249    for (;iWhich<numberInRowArray;iWhich++) {
3250      nextStart = rowStart[0];
3251      nextN = rowStart[numberInRowArray]-nextStart;
3252      rowStart++;
3253      if (nextN) {
3254        coin_prefetch(element+nextStart);
3255        nextPi = pi[iWhich];
3256        break;
3257      }
3258    }
3259    CoinBigIndex end =j+n;
3260    //coin_prefetch(element+rowStart_[i+1]);
3261    //coin_prefetch(column_+rowStart_[i+1]);
3262    if (n<100) {
3263      if ((n&1)!=0) {
3264        unsigned int jColumn = column[j];
3265        array[jColumn] -= value*element[j];
3266        j++;
3267      }     
3268      coin_prefetch(column+nextStart);
3269      for (;j<end;j+=2) {
3270        unsigned int jColumn0 = column[j];
3271        double value0 = value*element[j];
3272        unsigned int jColumn1 = column[j+1];
3273        double value1 = value*element[j+1];
3274        array[jColumn0] -= value0;
3275        array[jColumn1] -= value1;
3276      }
3277    } else {
3278      if ((n&1)!=0) {
3279        unsigned int jColumn = column[j];
3280        array[jColumn] -= value*element[j];
3281        j++;
3282      }     
3283      if ((n&2)!=0) {
3284        unsigned int jColumn0 = column[j];
3285        double value0 = value*element[j];
3286        unsigned int jColumn1 = column[j+1];
3287        double value1 = value*element[j+1];
3288        array[jColumn0] -= value0;
3289        array[jColumn1] -= value1;
3290        j+=2;
3291      }     
3292      if ((n&4)!=0) {
3293        unsigned int jColumn0 = column[j];
3294        double value0 = value*element[j];
3295        unsigned int jColumn1 = column[j+1];
3296        double value1 = value*element[j+1];
3297        unsigned int jColumn2 = column[j+2];
3298        double value2 = value*element[j+2];
3299        unsigned int jColumn3 = column[j+3];
3300        double value3 = value*element[j+3];
3301        array[jColumn0] -= value0;
3302        array[jColumn1] -= value1;
3303        array[jColumn2] -= value2;
3304        array[jColumn3] -= value3;
3305        j+=4;
3306      }
3307      //coin_prefetch(column+nextStart);
3308      for (;j<end;j+=8) {
3309        coin_prefetch(element+j+16);
3310        unsigned int jColumn0 = column[j];
3311        double value0 = value*element[j];
3312        unsigned int jColumn1 = column[j+1];
3313        double value1 = value*element[j+1];
3314        unsigned int jColumn2 = column[j+2];
3315        double value2 = value*element[j+2];
3316        unsigned int jColumn3 = column[j+3];
3317        double value3 = value*element[j+3];
3318        array[jColumn0] -= value0;
3319        array[jColumn1] -= value1;
3320        array[jColumn2] -= value2;
3321        array[jColumn3] -= value3;
3322        coin_prefetch(column+j+16);
3323        jColumn0 = column[j+4];
3324        value0 = value*element[j+4];
3325        jColumn1 = column[j+5];
3326        value1 = value*element[j+5];
3327        jColumn2 = column[j+6];
3328        value2 = value*element[j+6];
3329        jColumn3 = column[j+7];
3330        value3 = value*element[j+7];
3331        array[jColumn0] -= value0;
3332        array[jColumn1] -= value1;
3333        array[jColumn2] -= value2;
3334        array[jColumn3] -= value3;
3335      }
3336    }
3337  }
3338  // get rid of tiny values
3339  int nSmall = numberLook;
3340  int numberNonZero=0;
3341  for (i=0;i<nSmall;i++) {
3342    double value = array[i];
3343    array[i]=0.0; 
3344    if (fabs(value)>1.0e-12) {
3345      array[numberNonZero]=value;
3346      index[numberNonZero++]=i;
3347    }
3348  }
3349  for (;i<numberLook;i+=4) {
3350    double value0 = array[i+0];
3351    double value1 = array[i+1];
3352    double value2 = array[i+2];
3353    double value3 = array[i+3];
3354    array[i+0]=0.0; 
3355    array[i+1]=0.0; 
3356    array[i+2]=0.0; 
3357    array[i+3]=0.0; 
3358    if (fabs(value0)>1.0e-12) {
3359      array[numberNonZero]=value0;
3360      index[numberNonZero++]=i+0;
3361    }
3362    if (fabs(value1)>1.0e-12) {
3363      array[numberNonZero]=value1;
3364      index[numberNonZero++]=i+1;
3365    }
3366    if (fabs(value2)>1.0e-12) {
3367      array[numberNonZero]=value2;
3368      index[numberNonZero++]=i+2;
3369    }
3370    if (fabs(value3)>1.0e-12) {
3371      array[numberNonZero]=value3;
3372      index[numberNonZero++]=i+3;
3373    }
3374  }
3375  return numberNonZero;
3376}
3377#ifdef THREAD
3378static void * doOneBlockThread(void * voidInfo)
3379{
3380  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3381  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3382                                  info->rowStart,info->element,info->column,
3383                                  info->numberInRowArray,info->numberLook);
3384  return NULL;
3385}
3386static void * doOneBlockAnd0Thread(void * voidInfo)
3387{
3388  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3389  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3390                                  info->rowStart,info->element,info->column,
3391                                  info->numberInRowArray,info->numberLook);
3392  *(info->numberOutPtr) =  dualColumn0(info->model,info->spare,
3393                                    info->spareIndex,(const double *)info->arrayTemp,
3394                                    (const int *) info->indexTemp,*(info->numberInPtr),
3395                                    info->offset, info->acceptablePivot,info->bestPossiblePtr,
3396                                    info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
3397  return NULL;
3398}
3399#endif
3400/* Return <code>x * scalar * A in <code>z</code>.
3401   Note - x packed and z will be packed mode
3402   Squashes small elements and knows about ClpSimplex */
3403void 
3404ClpPackedMatrix2::transposeTimes(const ClpSimplex * model, 
3405                                 const CoinPackedMatrix * rowCopy,
3406                                 const CoinIndexedVector * rowArray,
3407                                 CoinIndexedVector * spareArray,
3408                                 CoinIndexedVector * columnArray) const
3409{
3410  // See if dualColumn0 coding wanted
3411  bool dualColumn = model->spareIntArray_[0]==1;
3412  double acceptablePivot = model->spareDoubleArray_[0];
3413  double bestPossible=0.0;
3414  double upperTheta=1.0e31;
3415  double freePivot = acceptablePivot;
3416  int posFree=-1;
3417  int numberRemaining=0;
3418  //if (model->numberIterations()>=200000) {
3419  //printf("time %g\n",CoinCpuTime()-startTime);
3420  //exit(77);
3421  //}
3422  double * pi = rowArray->denseVector();
3423  int numberNonZero=0;
3424  int * index = columnArray->getIndices();
3425  double * array = columnArray->denseVector();
3426  int numberInRowArray = rowArray->getNumElements();
3427  const int * whichRow = rowArray->getIndices();
3428  double * element = const_cast<double *>(rowCopy->getElements());
3429  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3430  int i;
3431  CoinBigIndex * rowStart2 = rowStart_;
3432  if (!dualColumn) {
3433    for (i=0;i<numberInRowArray;i++) {
3434      int iRow = whichRow[i]; 
3435      CoinBigIndex start = rowStart[iRow];
3436      *rowStart2 = start;
3437      unsigned short * count1 = count_+iRow*numberBlocks_;
3438      int put=0;
3439      for (int j=0;j<numberBlocks_;j++) {
3440        put += numberInRowArray;
3441        start += count1[j];
3442        rowStart2[put]=start;
3443      }
3444      rowStart2 ++;
3445    }
3446  } else {
3447    // also do dualColumn stuff
3448    double * spare = spareArray->denseVector();
3449    int * spareIndex = spareArray->getIndices();
3450    const double * reducedCost=model->djRegion(0);
3451    double dualTolerance = model->dualTolerance();
3452    // We can also see if infeasible or pivoting on free
3453    double tentativeTheta = 1.0e25;
3454    int addSequence=model->numberColumns();
3455    for (i=0;i<numberInRowArray;i++) {
3456      int iRow = whichRow[i]; 
3457      double alpha=pi[i];
3458      double oldValue;
3459      double value;
3460      bool keep;
3461
3462      switch(model->getStatus(iRow+addSequence)) {
3463         
3464      case ClpSimplex::basic:
3465      case ClpSimplex::isFixed:
3466        break;
3467      case ClpSimplex::isFree:
3468      case ClpSimplex::superBasic:
3469        bestPossible = CoinMax(bestPossible,fabs(alpha));
3470        oldValue = reducedCost[iRow];
3471        // If free has to be very large - should come in via dualRow
3472        if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3473          break;
3474        if (oldValue>dualTolerance) {
3475          keep = true;
3476        } else if (oldValue<-dualTolerance) {
3477          keep = true;
3478        } else {
3479          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3480            keep = true;
3481          else
3482            keep=false;
3483        }
3484        if (keep) {
3485          // free - choose largest
3486          if (fabs(alpha)>freePivot) {
3487            freePivot=fabs(alpha);
3488            posFree = i + addSequence;
3489          }
3490        }
3491        break;
3492      case ClpSimplex::atUpperBound:
3493        oldValue = reducedCost[iRow];
3494        value = oldValue-tentativeTheta*alpha;
3495        //assert (oldValue<=dualTolerance*1.0001);
3496        if (value>dualTolerance) {
3497          bestPossible = CoinMax(bestPossible,-alpha);
3498          value = oldValue-upperTheta*alpha;
3499          if (value>dualTolerance && -alpha>=acceptablePivot) 
3500            upperTheta = (oldValue-dualTolerance)/alpha;
3501          // add to list
3502          spare[numberRemaining]=alpha;
3503          spareIndex[numberRemaining++]=iRow+addSequence;
3504        }
3505        break;
3506      case ClpSimplex::atLowerBound:
3507        oldValue = reducedCost[iRow];
3508        value = oldValue-tentativeTheta*alpha;
3509        //assert (oldValue>=-dualTolerance*1.0001);
3510        if (value<-dualTolerance) {
3511          bestPossible = CoinMax(bestPossible,alpha);
3512          value = oldValue-upperTheta*alpha;
3513          if (value<-dualTolerance && alpha>=acceptablePivot) 
3514            upperTheta = (oldValue+dualTolerance)/alpha;
3515          // add to list
3516          spare[numberRemaining]=alpha;
3517          spareIndex[numberRemaining++]=iRow+addSequence;
3518        }
3519        break;
3520      }
3521      CoinBigIndex start = rowStart[iRow];
3522      *rowStart2 = start;
3523      unsigned short * count1 = count_+iRow*numberBlocks_;
3524      int put=0;
3525      for (int j=0;j<numberBlocks_;j++) {
3526        put += numberInRowArray;
3527        start += count1[j];
3528        rowStart2[put]=start;
3529      }
3530      rowStart2 ++;
3531    }
3532  }
3533  double * spare = spareArray->denseVector();
3534  int * spareIndex = spareArray->getIndices();
3535  int saveNumberRemaining=numberRemaining;
3536  int iBlock;
3537  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3538    double * dwork = work_+6*iBlock;
3539    int * iwork = (int *) (dwork+3);
3540    if (!dualColumn) {
3541#ifndef THREAD
3542      int offset = offset_[iBlock];
3543      int offset3 = offset;
3544      offset = numberNonZero;
3545      double * arrayTemp = array+offset;
3546      int * indexTemp = index+offset; 
3547      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
3548                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
3549      int number = iwork[0];
3550      for (i=0;i<number;i++) {
3551        //double value = arrayTemp[i];
3552        //arrayTemp[i]=0.0;
3553        //array[numberNonZero]=value;
3554        index[numberNonZero++]=indexTemp[i]+offset3;
3555      }
3556#else
3557      int offset = offset_[iBlock];
3558      double * arrayTemp = array+offset;
3559      int * indexTemp = index+offset; 
3560      dualColumn0Struct * infoPtr = info_+iBlock;
3561      infoPtr->arrayTemp=arrayTemp;
3562      infoPtr->indexTemp=indexTemp;
3563      infoPtr->numberInPtr=&iwork[0];
3564      infoPtr->pi=pi;
3565      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
3566      infoPtr->element=element;
3567      infoPtr->column=column_;
3568      infoPtr->numberInRowArray=numberInRowArray;
3569      infoPtr->numberLook=offset_[iBlock+1]-offset;
3570      pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
3571#endif
3572    } else {
3573#ifndef THREAD
3574      int offset = offset_[iBlock];
3575      // allow for already saved
3576      int offset2 = offset + saveNumberRemaining;
3577      int offset3 = offset;
3578      offset = numberNonZero;
3579      offset2 = numberRemaining;
3580      double * arrayTemp = array+offset;
3581      int * indexTemp = index+offset; 
3582      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
3583                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
3584      iwork[1] = dualColumn0(model,spare+offset2,
3585                                   spareIndex+offset2,
3586                                   arrayTemp,indexTemp,
3587                                   iwork[0],offset3,acceptablePivot,
3588                                   &dwork[0],&dwork[1],&iwork[2],
3589                                   &dwork[2]);
3590      int number = iwork[0];
3591      int numberLook = iwork[1];
3592#if 1
3593      numberRemaining += numberLook;
3594#else
3595      double * spareTemp = spare+offset2;
3596      const int * spareIndexTemp = spareIndex+offset2; 
3597      for (i=0;i<numberLook;i++) {
3598        double value = spareTemp[i];
3599        spareTemp[i]=0.0;
3600        spare[numberRemaining] = value;
3601        spareIndex[numberRemaining++] = spareIndexTemp[i];
3602      }
3603#endif
3604      if (dwork[2]>freePivot) {
3605        freePivot=dwork[2];
3606        posFree = iwork[2]+numberNonZero;
3607      } 
3608      upperTheta =  CoinMin(dwork[1],upperTheta);
3609      bestPossible = CoinMax(dwork[0],bestPossible);
3610      for (i=0;i<number;i++) {
3611        // double value = arrayTemp[i];
3612        //arrayTemp[i]=0.0;
3613        //array[numberNonZero]=value;
3614        index[numberNonZero++]=indexTemp[i]+offset3;
3615      }
3616#else
3617      int offset = offset_[iBlock];
3618      // allow for already saved
3619      int offset2 = offset + saveNumberRemaining;
3620      double * arrayTemp = array+offset;
3621      int * indexTemp = index+offset; 
3622      dualColumn0Struct * infoPtr = info_+iBlock;
3623      infoPtr->model=model;
3624      infoPtr->spare=spare+offset2;
3625      infoPtr->spareIndex=spareIndex+offset2;
3626      infoPtr->arrayTemp=arrayTemp;
3627      infoPtr->indexTemp=indexTemp;
3628      infoPtr->numberInPtr=&iwork[0];
3629      infoPtr->offset=offset;
3630      infoPtr->acceptablePivot=acceptablePivot;
3631      infoPtr->bestPossiblePtr=&dwork[0];
3632      infoPtr->upperThetaPtr=&dwork[1];
3633      infoPtr->posFreePtr=&iwork[2];
3634      infoPtr->freePivotPtr=&dwork[2];
3635      infoPtr->numberOutPtr=&iwork[1];
3636      infoPtr->pi=pi;
3637      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
3638      infoPtr->element=element;
3639      infoPtr->column=column_;
3640      infoPtr->numberInRowArray=numberInRowArray;
3641      infoPtr->numberLook=offset_[iBlock+1]-offset;
3642      if (iBlock>=2)
3643        pthread_join(threadId_[iBlock-2],NULL);
3644      pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
3645      //pthread_join(threadId_[iBlock],NULL);
3646#endif
3647    }
3648  }
3649  for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
3650#ifdef THREAD
3651    pthread_join(threadId_[iBlock],NULL);
3652#endif
3653  }
3654#ifdef THREAD
3655  for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
3656    //pthread_join(threadId_[iBlock],NULL);
3657    int offset = offset_[iBlock];
3658    double * dwork = work_+6*iBlock;
3659    int * iwork = (int *) (dwork+3);
3660    int number = iwork[0];
3661    if (dualColumn) {
3662      // allow for already saved
3663      int offset2 = offset + saveNumberRemaining;
3664      int numberLook = iwork[1];
3665      double * spareTemp = spare+offset2;
3666      const int * spareIndexTemp = spareIndex+offset2; 
3667      for (i=0;i<numberLook;i++) {
3668        double value = spareTemp[i];
3669        spareTemp[i]=0.0;
3670        spare[numberRemaining] = value;
3671        spareIndex[numberRemaining++] = spareIndexTemp[i];
3672      }
3673      if (dwork[2]>freePivot) {
3674        freePivot=dwork[2];
3675        posFree = iwork[2]+numberNonZero;
3676      } 
3677      upperTheta =  CoinMin(dwork[1],upperTheta);
3678      bestPossible = CoinMax(dwork[0],bestPossible);
3679    }
3680    double * arrayTemp = array+offset;
3681    const int * indexTemp = index+offset; 
3682    for (i=0;i<number;i++) {
3683      double value = arrayTemp[i];
3684      arrayTemp[i]=0.0; 
3685      array[numberNonZero]=value;
3686      index[numberNonZero++]=indexTemp[i]+offset;
3687    }
3688  }
3689#endif
3690  columnArray->setNumElements(numberNonZero);
3691  columnArray->setPackedMode(true);
3692  if (dualColumn) {
3693    model->spareDoubleArray_[0]=upperTheta;
3694    model->spareDoubleArray_[1]=bestPossible;
3695    // and theta and alpha and sequence
3696    if (posFree<0) {
3697      model->spareIntArray_[1]=-1;
3698    } else {
3699      const double * reducedCost=model->djRegion(0);
3700      double alpha;
3701      int numberColumns=model->numberColumns();
3702      if (posFree<numberColumns) {
3703        alpha = columnArray->denseVector()[posFree];
3704        posFree = columnArray->getIndices()[posFree];
3705      } else {
3706        alpha = rowArray->denseVector()[posFree-numberColumns];
3707        posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
3708      }
3709      model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
3710      model->spareDoubleArray_[3]=alpha;
3711      model->spareIntArray_[1]=posFree;
3712    }
3713    spareArray->setNumElements(numberRemaining);
3714    // signal done
3715    model->spareIntArray_[0]=-1;
3716  }
3717}
3718#ifdef CLP_ALL_ONE_FILE
3719#undef reference
3720#endif
Note: See TracBrowser for help on using the repository browser.