source: branches/devel/Clp/src/ClpPackedMatrix.cpp @ 989

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

this may be a mistake

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