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

Last change on this file since 1141 was 1141, checked in by forrest, 12 years ago

for threaded random numbers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 148.0 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); 
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    if (scalingMethod==4) {
2177      // Aas auto
2178      scalingMethod=3;
2179    }
2180    // and see if there any empty rows
2181    for (iRow=0;iRow<numberRows;iRow++) {
2182      if (usefulRow[iRow]) {
2183        CoinBigIndex j;
2184        int useful=0;
2185        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
2186          int iColumn = column[j];
2187          if (usefulColumn[iColumn]) {
2188            useful=1;
2189            break;
2190          }
2191        }
2192        usefulRow[iRow]=useful;
2193      }
2194    }
2195    double savedOverallRatio=0.0;
2196    double tolerance = 5.0*model->primalTolerance();
2197    double overallLargest=-1.0e-20;
2198    double overallSmallest=1.0e20;
2199    bool finished=false;
2200    // if scalingMethod 3 then may change
2201    while (!finished) {
2202      ClpFillN ( rowScale, numberRows,1.0);
2203      ClpFillN ( columnScale, numberColumns,1.0);
2204      overallLargest=-1.0e-20;
2205      overallSmallest=1.0e20;
2206      if (scalingMethod==1||scalingMethod==3) {
2207        // Maximum in each row
2208        for (iRow=0;iRow<numberRows;iRow++) {
2209          if (usefulRow[iRow]) {
2210            CoinBigIndex j;
2211            largest=1.0e-10;
2212            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
2213              int iColumn = column[j];
2214              if (usefulColumn[iColumn]) {
2215                double value = fabs(element[j]);
2216                largest = CoinMax(largest,value);
2217                assert (largest<1.0e40);
2218              }
2219            }
2220            rowScale[iRow]=1.0/largest;
2221            overallLargest = CoinMax(overallLargest,largest);
2222            overallSmallest = CoinMin(overallSmallest,largest);
2223          }
2224        }
2225      } else {
2226        int numberPass=3;
2227#ifdef USE_OBJECTIVE
2228        // This will be used to help get scale factors
2229        double * objective = new double[numberColumns];
2230        memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
2231        double objScale=1.0;
2232#endif
2233        while (numberPass) {
2234          overallLargest=0.0;
2235          overallSmallest=1.0e50;
2236          numberPass--;
2237          // Geometric mean on row scales
2238          for (iRow=0;iRow<numberRows;iRow++) {
2239            if (usefulRow[iRow]) {
2240              CoinBigIndex j;
2241              largest=1.0e-20;
2242              smallest=1.0e50;
2243              for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
2244                int iColumn = column[j];
2245                if (usefulColumn[iColumn]) {
2246                  double value = fabs(element[j]);
2247                  // Don't bother with tiny elements
2248                  if (value>1.0e-20) {
2249                    value *= columnScale[iColumn];
2250                    largest = CoinMax(largest,value);
2251                    smallest = CoinMin(smallest,value);
2252                  }
2253                }
2254              }
2255              rowScale[iRow]=1.0/sqrt(smallest*largest);
2256              rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
2257              overallLargest = CoinMax(largest*rowScale[iRow],overallLargest);
2258              overallSmallest = CoinMin(smallest*rowScale[iRow],overallSmallest);
2259            }
2260          }
2261#ifdef USE_OBJECTIVE
2262          largest=1.0e-20;
2263          smallest=1.0e50;
2264          for (iColumn=0;iColumn<numberColumns;iColumn++) {
2265            if (usefulColumn[iColumn]) {
2266              double value = fabs(objective[iColumn]);
2267              // Don't bother with tiny elements
2268              if (value>1.0e-20) {
2269                value *= columnScale[iColumn];
2270                largest = CoinMax(largest,value);
2271                smallest = CoinMin(smallest,value);
2272              }
2273            }
2274          }
2275          objScale=1.0/sqrt(smallest*largest);
2276#endif
2277          model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
2278            <<overallSmallest
2279            <<overallLargest
2280            <<CoinMessageEol;
2281          // skip last column round
2282          if (numberPass==1)
2283            break;
2284          // Geometric mean on column scales
2285          for (iColumn=0;iColumn<numberColumns;iColumn++) {
2286            if (usefulColumn[iColumn]) {
2287              CoinBigIndex j;
2288              largest=1.0e-20;
2289              smallest=1.0e50;
2290              for (j=columnStart[iColumn];
2291                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
2292                iRow=row[j];
2293                double value = fabs(elementByColumn[j]);
2294                // Don't bother with tiny elements
2295                if (value>1.0e-20&&usefulRow[iRow]) {
2296                  value *= rowScale[iRow];
2297                  largest = CoinMax(largest,value);
2298                  smallest = CoinMin(smallest,value);
2299                }
2300              }
2301#ifdef USE_OBJECTIVE
2302              if (fabs(objective[iColumn])>1.0e-20) {
2303                double value = fabs(objective[iColumn])*objScale;
2304                largest = CoinMax(largest,value);
2305                smallest = CoinMin(smallest,value);
2306              }
2307#endif
2308              columnScale[iColumn]=1.0/sqrt(smallest*largest);
2309              columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
2310            }
2311          }
2312        }
2313#ifdef USE_OBJECTIVE
2314        delete [] objective;
2315        printf("obj scale %g - use it if you want\n",objScale);
2316#endif
2317      }
2318      // If ranges will make horrid then scale
2319      for (iRow=0;iRow<numberRows;iRow++) {
2320        if (usefulRow[iRow]) {
2321          double difference = rowUpper[iRow]-rowLower[iRow];
2322          double scaledDifference = difference*rowScale[iRow];
2323          if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
2324            // make gap larger
2325            rowScale[iRow] *= 1.0e-4/scaledDifference;
2326            rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
2327            //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
2328            // scaledDifference,difference*rowScale[iRow]);
2329          }
2330        }
2331      }
2332      // final pass to scale columns so largest is reasonable
2333      // See what smallest will be if largest is 1.0
2334      overallSmallest=1.0e50;
2335      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2336        if (usefulColumn[iColumn]) {
2337          CoinBigIndex j;
2338          largest=1.0e-20;
2339          smallest=1.0e50;
2340          for (j=columnStart[iColumn];
2341               j<columnStart[iColumn]+columnLength[iColumn];j++) {
2342            iRow=row[j];
2343            if(elementByColumn[j]&&usefulRow[iRow]) {
2344              double value = fabs(elementByColumn[j]*rowScale[iRow]);
2345              largest = CoinMax(largest,value);
2346              smallest = CoinMin(smallest,value);
2347            }
2348          }
2349          if (overallSmallest*largest>smallest)
2350            overallSmallest = smallest/largest;
2351        }
2352      }
2353      if (scalingMethod==1||scalingMethod==2) {
2354        finished=true;
2355      } else if (savedOverallRatio==0.0&&scalingMethod!=4) {
2356        savedOverallRatio=overallSmallest;
2357        scalingMethod=4;
2358      } else {
2359        assert (scalingMethod==4);
2360        if (overallSmallest>2.0*savedOverallRatio) {
2361          finished=true; // geometric was better
2362          if (model->scalingFlag()==4) {
2363            // If in Branch and bound change
2364            if ((model->specialOptions()&1024)!=0) {
2365              model->scaling(2);
2366            }
2367          }
2368        } else {
2369          scalingMethod=1; // redo equilibrium
2370          if (model->scalingFlag()==4) {
2371            // If in Branch and bound change
2372            if ((model->specialOptions()&1024)!=0) {
2373              model->scaling(1);
2374            }
2375          }
2376        }
2377#if 0
2378        if (model->logLevel()>2) {
2379          if (finished)
2380            printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
2381                   savedOverallRatio,overallSmallest);
2382          else
2383            printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
2384                   savedOverallRatio,overallSmallest);
2385        }
2386#endif
2387      }
2388    }
2389    //#define RANDOMIZE
2390#ifdef RANDOMIZE
2391    // randomize by up to 10%
2392    for (iRow=0;iRow<numberRows;iRow++) {
2393      double value = 0.5-randomNumberGenerator_.randomDouble();//between -0.5 to + 0.5
2394      rowScale[iRow] *= (1.0+0.1*value);
2395    }
2396#endif
2397    overallLargest=1.0;
2398    if (overallSmallest<1.0e-1)
2399      overallLargest = 1.0/sqrt(overallSmallest);
2400    overallLargest = CoinMin(100.0,overallLargest);
2401    overallSmallest=1.0e50;
2402    //printf("scaling %d\n",model->scalingFlag());
2403    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2404      if (columnUpper[iColumn]>
2405        columnLower[iColumn]+1.0e-12) {
2406        //if (usefulColumn[iColumn]) {
2407        CoinBigIndex j;
2408        largest=1.0e-20;
2409        smallest=1.0e50;
2410        for (j=columnStart[iColumn];
2411             j<columnStart[iColumn]+columnLength[iColumn];j++) {
2412          iRow=row[j];
2413          if(elementByColumn[j]&&usefulRow[iRow]) {
2414            double value = fabs(elementByColumn[j]*rowScale[iRow]);
2415            largest = CoinMax(largest,value);
2416            smallest = CoinMin(smallest,value);
2417          }
2418        }
2419        columnScale[iColumn]=overallLargest/largest;
2420        columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
2421#ifdef RANDOMIZE
2422        double value = 0.5-randomNumberGenerator_.randomDouble();//between -0.5 to + 0.5
2423        columnScale[iColumn] *= (1.0+0.1*value);
2424#endif
2425        double difference = columnUpper[iColumn]-columnLower[iColumn];
2426        if (difference<1.0e-5*columnScale[iColumn]) {
2427          // make gap larger
2428          columnScale[iColumn] = difference/1.0e-5;
2429          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
2430          // scaledDifference,difference*columnScale[iColumn]);
2431        }
2432        overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
2433      }
2434    }
2435    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
2436      <<overallSmallest
2437      <<overallLargest
2438      <<CoinMessageEol;
2439    delete [] usefulRow;
2440    delete [] usefulColumn;
2441#ifndef SLIM_CLP
2442    // If quadratic then make symmetric
2443    ClpObjective * obj = model->objectiveAsObject();
2444#ifndef NO_RTTI
2445    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2446#else
2447    ClpQuadraticObjective * quadraticObj = NULL;
2448    if (obj->type()==2)
2449      quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
2450#endif
2451    if (quadraticObj) {
2452      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2453      int numberXColumns = quadratic->getNumCols();
2454      if (numberXColumns<numberColumns) {
2455        // we assume symmetric
2456        int numberQuadraticColumns=0;
2457        int i;
2458        //const int * columnQuadratic = quadratic->getIndices();
2459        //const int * columnQuadraticStart = quadratic->getVectorStarts();
2460        const int * columnQuadraticLength = quadratic->getVectorLengths();
2461        for (i=0;i<numberXColumns;i++) {
2462          int length=columnQuadraticLength[i];
2463#ifndef CORRECT_COLUMN_COUNTS
2464          length=1;
2465#endif
2466          if (length)
2467            numberQuadraticColumns++;
2468        }
2469        int numberXRows = numberRows-numberQuadraticColumns;
2470        numberQuadraticColumns=0;
2471        for (i=0;i<numberXColumns;i++) { 
2472          int length=columnQuadraticLength[i];
2473#ifndef CORRECT_COLUMN_COUNTS
2474          length=1;
2475#endif
2476          if (length) {
2477            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
2478            numberQuadraticColumns++;
2479          }
2480        }   
2481        int numberQuadraticRows=0;
2482        for (i=0;i<numberXRows;i++) {
2483          // See if any in row quadratic
2484          CoinBigIndex j;
2485          int numberQ=0;
2486          for (j=rowStart[i];j<rowStart[i+1];j++) {
2487            int iColumn = column[j];
2488            if (columnQuadraticLength[iColumn])
2489              numberQ++;
2490          }
2491#ifndef CORRECT_ROW_COUNTS
2492          numberQ=1;
2493#endif
2494          if (numberQ) {
2495            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
2496            numberQuadraticRows++;
2497          }
2498        }
2499        // and make sure Sj okay
2500        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
2501          CoinBigIndex j=columnStart[iColumn];
2502          assert(columnLength[iColumn]==1);
2503          int iRow=row[j];
2504          double value = fabs(elementByColumn[j]*rowScale[iRow]);
2505          columnScale[iColumn]=1.0/value;
2506        }
2507      }
2508    }
2509#endif
2510    if (scaleLength>1) {
2511      // make copy
2512      double * inverseScale = rowScale + numberRows;
2513      for (iRow=0;iRow<numberRows;iRow++) 
2514        inverseScale[iRow] = 1.0/rowScale[iRow] ;
2515      inverseScale = columnScale + numberColumns;
2516      for (iColumn=0;iColumn<numberColumns;iColumn++) 
2517        inverseScale[iColumn] = 1.0/columnScale[iColumn] ;
2518    }
2519    model->setRowScale(rowScale);
2520    model->setColumnScale(columnScale);
2521    if (model->rowCopy()) {
2522      // need to replace row by row
2523      double * newElement = new double[numberColumns];
2524      // scale row copy
2525      for (iRow=0;iRow<numberRows;iRow++) {
2526        CoinBigIndex j;
2527        double scale = rowScale[iRow];
2528        const double * elementsInThisRow = element + rowStart[iRow];
2529        const int * columnsInThisRow = column + rowStart[iRow];
2530        int number = rowStart[iRow+1]-rowStart[iRow];
2531        assert (number<=numberColumns);
2532        for (j=0;j<number;j++) {
2533          int iColumn = columnsInThisRow[j];
2534          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
2535        }
2536        rowCopy->replaceVector(iRow,number,newElement);
2537      }
2538      delete [] newElement;
2539    } else {
2540      // no row copy
2541      delete rowCopyBase;
2542    }
2543    return 0;
2544  }
2545}
2546/* Unpacks a column into an CoinIndexedvector
2547 */
2548void 
2549ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
2550                   int iColumn) const 
2551{
2552  const double * rowScale = model->rowScale();
2553  const int * row = matrix_->getIndices();
2554  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2555  const int * columnLength = matrix_->getVectorLengths(); 
2556  const double * elementByColumn = matrix_->getElements();
2557  CoinBigIndex i;
2558  if (!rowScale) {
2559    for (i=columnStart[iColumn];
2560         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2561      rowArray->add(row[i],elementByColumn[i]);
2562    }
2563  } else {
2564    // apply scaling
2565    double scale = model->columnScale()[iColumn];
2566    for (i=columnStart[iColumn];
2567         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2568      int iRow = row[i];
2569      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
2570    }
2571  }
2572}
2573/* Unpacks a column into a CoinIndexedVector
2574** in packed format
2575Note that model is NOT const.  Bounds and objective could
2576be modified if doing column generation (just for this variable) */
2577void 
2578ClpPackedMatrix::unpackPacked(ClpSimplex * model,
2579                            CoinIndexedVector * rowArray,
2580                            int iColumn) const
2581{
2582  const double * rowScale = model->rowScale();
2583  const int * row = matrix_->getIndices();
2584  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2585  const int * columnLength = matrix_->getVectorLengths(); 
2586  const double * elementByColumn = matrix_->getElements();
2587  CoinBigIndex i;
2588  int * index = rowArray->getIndices();
2589  double * array = rowArray->denseVector();
2590  int number = 0;
2591  if (!rowScale) {
2592    for (i=columnStart[iColumn];
2593         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2594      int iRow = row[i];
2595      double value = elementByColumn[i];
2596      if (value) {
2597        array[number]=value;
2598        index[number++]=iRow;
2599      }
2600    }
2601    rowArray->setNumElements(number);
2602    rowArray->setPackedMode(true);
2603  } else {
2604    // apply scaling
2605    double scale = model->columnScale()[iColumn];
2606    for (i=columnStart[iColumn];
2607         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2608      int iRow = row[i];
2609      double value = elementByColumn[i]*scale*rowScale[iRow];
2610      if (value) {
2611        array[number]=value;
2612        index[number++]=iRow;
2613      }
2614    }
2615    rowArray->setNumElements(number);
2616    rowArray->setPackedMode(true);
2617  }
2618}
2619/* Adds multiple of a column into an CoinIndexedvector
2620      You can use quickAdd to add to vector */
2621void 
2622ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
2623                   int iColumn, double multiplier) const 
2624{
2625  const double * rowScale = model->rowScale();
2626  const int * row = matrix_->getIndices();
2627  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2628  const int * columnLength = matrix_->getVectorLengths(); 
2629  const double * elementByColumn = matrix_->getElements();
2630  CoinBigIndex i;
2631  if (!rowScale) {
2632    for (i=columnStart[iColumn];
2633         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2634      int iRow = row[i];
2635      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
2636    }
2637  } else {
2638    // apply scaling
2639    double scale = model->columnScale()[iColumn]*multiplier;
2640    for (i=columnStart[iColumn];
2641         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2642      int iRow = row[i];
2643      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
2644    }
2645  }
2646}
2647/* Adds multiple of a column into an array */
2648void 
2649ClpPackedMatrix::add(const ClpSimplex * model,double * array,
2650                    int iColumn, double multiplier) const
2651{
2652  const double * rowScale = model->rowScale();
2653  const int * row = matrix_->getIndices();
2654  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2655  const int * columnLength = matrix_->getVectorLengths(); 
2656  const double * elementByColumn = matrix_->getElements();
2657  CoinBigIndex i;
2658  if (!rowScale) {
2659    for (i=columnStart[iColumn];
2660         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2661      int iRow = row[i];
2662      array[iRow] += multiplier*elementByColumn[i];
2663    }
2664  } else {
2665    // apply scaling
2666    double scale = model->columnScale()[iColumn]*multiplier;
2667    for (i=columnStart[iColumn];
2668         i<columnStart[iColumn]+columnLength[iColumn];i++) {
2669      int iRow = row[i];
2670      array[iRow] += elementByColumn[i]*scale*rowScale[iRow];
2671    }
2672  }
2673}
2674/* Checks if all elements are in valid range.  Can just
2675   return true if you are not paranoid.  For Clp I will
2676   probably expect no zeros.  Code can modify matrix to get rid of
2677   small elements.
2678   check bits (can be turned off to save time) :
2679   1 - check if matrix has gaps
2680   2 - check if zero elements
2681   4 - check and compress duplicates
2682   8 - report on large and small
2683*/
2684bool 
2685ClpPackedMatrix::allElementsInRange(ClpModel * model,
2686                                    double smallest, double largest,
2687                                    int check)
2688{
2689  int iColumn;
2690  // make sure matrix correct size
2691  matrix_->setDimensions(model->numberRows(),model->numberColumns());
2692  CoinBigIndex numberLarge=0;;
2693  CoinBigIndex numberSmall=0;;
2694  CoinBigIndex numberDuplicate=0;;
2695  int firstBadColumn=-1;
2696  int firstBadRow=-1;
2697  double firstBadElement=0.0;
2698  // get matrix data pointers
2699  const int * row = matrix_->getIndices();
2700  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2701  const int * columnLength = matrix_->getVectorLengths(); 
2702  const double * elementByColumn = matrix_->getElements();
2703  int numberRows = matrix_->getNumRows();
2704  int numberColumns = matrix_->getNumCols();
2705  // Say no gaps
2706  flags_ &= ~2;
2707  if (check==14||check==10) {
2708    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2709      CoinBigIndex start = columnStart[iColumn];
2710      CoinBigIndex end = start + columnLength[iColumn];
2711      if (end!=columnStart[iColumn+1]) {
2712        flags_ |= 2;
2713        break;
2714      }
2715    }
2716    return true;
2717  }
2718  assert (check==15||check==11);
2719  if (check==15) {
2720    int * mark = new int [numberRows];
2721    int i;
2722    for (i=0;i<numberRows;i++)
2723      mark[i]=-1;
2724    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2725      CoinBigIndex j;
2726      CoinBigIndex start = columnStart[iColumn];
2727      CoinBigIndex end = start + columnLength[iColumn];
2728      if (end!=columnStart[iColumn+1])
2729        flags_ |= 2;
2730      for (j=start;j<end;j++) {
2731        double value = fabs(elementByColumn[j]);
2732        int iRow = row[j];
2733        if (iRow<0||iRow>=numberRows) {
2734#ifndef COIN_BIG_INDEX
2735          printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2736#elif COIN_BIG_INDEX==1
2737          printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2738#else
2739          printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2740#endif
2741          return false;
2742        }
2743        if (mark[iRow]==-1) {
2744          mark[iRow]=j;
2745        } else {
2746          // duplicate
2747          numberDuplicate++;
2748        }
2749        //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2750        if (!value)
2751          flags_ |= 1; // there are zero elements
2752        if (value<smallest) {
2753          numberSmall++;
2754        } else if (!(value<=largest)) {
2755          numberLarge++;
2756          if (firstBadColumn<0) {
2757            firstBadColumn=iColumn;
2758            firstBadRow=row[j];
2759            firstBadElement=elementByColumn[j];
2760          }
2761        }
2762      }
2763      //clear mark
2764      for (j=columnStart[iColumn];
2765           j<columnStart[iColumn]+columnLength[iColumn];j++) {
2766        int iRow = row[j];
2767        mark[iRow]=-1;
2768      }
2769    }
2770    delete [] mark;
2771  } else {
2772    // just check for out of range - not for duplicates
2773    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2774      CoinBigIndex j;
2775      CoinBigIndex start = columnStart[iColumn];
2776      CoinBigIndex end = start + columnLength[iColumn];
2777      if (end!=columnStart[iColumn+1])
2778        flags_ |= 2;
2779      for (j=start;j<end;j++) {
2780        double value = fabs(elementByColumn[j]);
2781        int iRow = row[j];
2782        if (iRow<0||iRow>=numberRows) {
2783#ifndef COIN_BIG_INDEX
2784          printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2785#elif COIN_BIG_INDEX==1
2786          printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2787#else
2788          printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
2789#endif
2790          return false;
2791        }
2792        if (!value)
2793          flags_ |= 1; // there are zero elements
2794        if (value<smallest) {
2795          numberSmall++;
2796        } else if (!(value<=largest)) {
2797          numberLarge++;
2798          if (firstBadColumn<0) {
2799            firstBadColumn=iColumn;
2800            firstBadRow=iRow;
2801            firstBadElement=value;
2802          }
2803        }
2804      }
2805    }
2806  }
2807  if (numberLarge) {
2808    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
2809      <<numberLarge
2810      <<firstBadColumn<<firstBadRow<<firstBadElement
2811      <<CoinMessageEol;
2812    return false;
2813  }
2814  if (numberSmall) 
2815    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
2816      <<numberSmall
2817      <<CoinMessageEol;
2818  if (numberDuplicate) 
2819    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
2820      <<numberDuplicate
2821      <<CoinMessageEol;
2822  if (numberDuplicate) 
2823    matrix_->eliminateDuplicates(smallest);
2824  else if (numberSmall) 
2825    matrix_->compress(smallest);
2826  // If smallest >0.0 then there can't be zero elements
2827  if (smallest>0.0)
2828    flags_ &= ~1;;
2829  return true;
2830}
2831/* Given positive integer weights for each row fills in sum of weights
2832   for each column (and slack).
2833   Returns weights vector
2834*/
2835CoinBigIndex * 
2836ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
2837{
2838  int numberRows = model->numberRows();
2839  int numberColumns = matrix_->getNumCols();
2840  int number = numberRows+numberColumns;
2841  CoinBigIndex * weights = new CoinBigIndex[number];
2842  // get matrix data pointers
2843  const int * row = matrix_->getIndices();
2844  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2845  const int * columnLength = matrix_->getVectorLengths(); 
2846  int i;
2847  for (i=0;i<numberColumns;i++) {
2848    CoinBigIndex j;
2849    CoinBigIndex count=0;
2850    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
2851      int iRow=row[j];
2852      count += inputWeights[iRow];
2853    }
2854    weights[i]=count;
2855  }
2856  for (i=0;i<numberRows;i++) {
2857    weights[i+numberColumns]=inputWeights[i];
2858  }
2859  return weights;
2860}
2861/* Returns largest and smallest elements of both signs.
2862   Largest refers to largest absolute value.
2863*/
2864void 
2865ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
2866                       double & smallestPositive, double & largestPositive)
2867{
2868  smallestNegative=-COIN_DBL_MAX;
2869  largestNegative=0.0;
2870  smallestPositive=COIN_DBL_MAX;
2871  largestPositive=0.0;
2872  // get matrix data pointers
2873  const double * elementByColumn = matrix_->getElements();
2874  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2875  const int * columnLength = matrix_->getVectorLengths(); 
2876  int numberColumns = matrix_->getNumCols();
2877  int i;
2878  for (i=0;i<numberColumns;i++) {
2879    CoinBigIndex j;
2880    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
2881      double value = elementByColumn[j];
2882      if (value>0.0) {
2883        smallestPositive = CoinMin(smallestPositive,value);
2884        largestPositive = CoinMax(largestPositive,value);
2885      } else if (value<0.0) {
2886        smallestNegative = CoinMax(smallestNegative,value);
2887        largestNegative = CoinMin(largestNegative,value);
2888      }
2889    }
2890  }
2891}
2892// Says whether it can do partial pricing
2893bool 
2894ClpPackedMatrix::canDoPartialPricing() const
2895{
2896  return true;
2897}
2898// Partial pricing
2899void 
2900ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
2901                              int & bestSequence, int & numberWanted)
2902{
2903  numberWanted=currentWanted_;
2904  int start = (int) (startFraction*numberActiveColumns_);
2905  int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
2906  const double * element =matrix_->getElements();
2907  const int * row = matrix_->getIndices();
2908  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
2909  const int * length = matrix_->getVectorLengths();
2910  const double * rowScale = model->rowScale();
2911  const double * columnScale = model->columnScale();
2912  int iSequence;
2913  CoinBigIndex j;
2914  double tolerance=model->currentDualTolerance();
2915  double * reducedCost = model->djRegion();
2916  const double * duals = model->dualRowSolution();
2917  const double * cost = model->costRegion();
2918  double bestDj;
2919  if (bestSequence>=0)
2920    bestDj = fabs(model->clpMatrix()->reducedCost(model,bestSequence));
2921  else
2922    bestDj=tolerance;
2923  int sequenceOut = model->sequenceOut();
2924  int saveSequence = bestSequence;
2925  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 
2926  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
2927  if (rowScale) {
2928    // scaled
2929    for (iSequence=start;iSequence<end;iSequence++) {
2930      if (iSequence!=sequenceOut) {
2931        double value;
2932        ClpSimplex::Status status = model->getStatus(iSequence);
2933       
2934        switch(status) {
2935         
2936        case ClpSimplex::basic:
2937        case ClpSimplex::isFixed:
2938          break;
2939        case ClpSimplex::isFree:
2940        case ClpSimplex::superBasic:
2941          value=0.0;
2942          // scaled
2943          for (j=startColumn[iSequence];
2944               j<startColumn[iSequence]+length[iSequence];j++) {
2945            int jRow=row[j];
2946            value -= duals[jRow]*element[j]*rowScale[jRow];
2947          }
2948          value = fabs(cost[iSequence] +value*columnScale[iSequence]);
2949          if (value>FREE_ACCEPT*tolerance) {
2950            numberWanted--;
2951            // we are going to bias towards free (but only if reasonable)
2952            value *= FREE_BIAS;
2953            if (value>bestDj) {
2954              // check flagged variable and correct dj
2955              if (!model->flagged(iSequence)) {
2956                bestDj=value;
2957                bestSequence = iSequence;
2958              } else {
2959                // just to make sure we don't exit before got something
2960                numberWanted++;
2961              }
2962            }
2963          }
2964          break;
2965        case ClpSimplex::atUpperBound:
2966          value=0.0;
2967          // scaled
2968          for (j=startColumn[iSequence];
2969               j<startColumn[iSequence]+length[iSequence];j++) {
2970            int jRow=row[j];
2971            value -= duals[jRow]*element[j]*rowScale[jRow];
2972          }
2973          value = cost[iSequence] +value*columnScale[iSequence];
2974          if (value>tolerance) {
2975            numberWanted--;
2976            if (value>bestDj) {
2977              // check flagged variable and correct dj
2978              if (!model->flagged(iSequence)) {
2979                bestDj=value;
2980                bestSequence = iSequence;
2981              } else {
2982                // just to make sure we don't exit before got something
2983                numberWanted++;
2984              }
2985            }
2986          }
2987          break;
2988        case ClpSimplex::atLowerBound:
2989          value=0.0;
2990          // scaled
2991          for (j=startColumn[iSequence];
2992               j<startColumn[iSequence]+length[iSequence];j++) {
2993            int jRow=row[j];
2994            value -= duals[jRow]*element[j]*rowScale[jRow];
2995          }
2996          value = -(cost[iSequence] +value*columnScale[iSequence]);
2997          if (value>tolerance) {
2998            numberWanted--;
2999            if (value>bestDj) {
3000              // check flagged variable and correct dj
3001              if (!model->flagged(iSequence)) {
3002                bestDj=value;
3003                bestSequence = iSequence;
3004              } else {
3005                // just to make sure we don't exit before got something
3006                numberWanted++;
3007              }
3008            }
3009          }
3010          break;
3011        }
3012      }
3013      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
3014        // give up
3015        break;
3016      }
3017      if (!numberWanted)
3018        break;
3019    }
3020    if (bestSequence!=saveSequence) {
3021      // recompute dj
3022      double value=0.0;
3023      // scaled
3024      for (j=startColumn[bestSequence];
3025           j<startColumn[bestSequence]+length[bestSequence];j++) {
3026        int jRow=row[j];
3027        value -= duals[jRow]*element[j]*rowScale[jRow];
3028      }
3029      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
3030      savedBestSequence_ = bestSequence;
3031      savedBestDj_ = reducedCost[savedBestSequence_];
3032    }
3033  } else {
3034    // not scaled
3035    for (iSequence=start;iSequence<end;iSequence++) {
3036      if (iSequence!=sequenceOut) {
3037        double value;
3038        ClpSimplex::Status status = model->getStatus(iSequence);
3039       
3040        switch(status) {
3041         
3042        case ClpSimplex::basic:
3043        case ClpSimplex::isFixed:
3044          break;
3045        case ClpSimplex::isFree:
3046        case ClpSimplex::superBasic:
3047          value=cost[iSequence];
3048          for (j=startColumn[iSequence];
3049               j<startColumn[iSequence]+length[iSequence];j++) {
3050            int jRow=row[j];
3051            value -= duals[jRow]*element[j];
3052          }
3053          value = fabs(value);
3054          if (value>FREE_ACCEPT*tolerance) {
3055            numberWanted--;
3056            // we are going to bias towards free (but only if reasonable)
3057            value *= FREE_BIAS;
3058            if (value>bestDj) {
3059              // check flagged variable and correct dj
3060              if (!model->flagged(iSequence)) {
3061                bestDj=value;
3062                bestSequence = iSequence;
3063              } else {
3064                // just to make sure we don't exit before got something
3065                numberWanted++;
3066              }
3067            }
3068          }
3069          break;
3070        case ClpSimplex::atUpperBound:
3071          value=cost[iSequence];
3072          // scaled
3073          for (j=startColumn[iSequence];
3074               j<startColumn[iSequence]+length[iSequence];j++) {
3075            int jRow=row[j];
3076            value -= duals[jRow]*element[j];
3077          }
3078          if (value>tolerance) {
3079            numberWanted--;
3080            if (value>bestDj) {
3081              // check flagged variable and correct dj
3082              if (!model->flagged(iSequence)) {
3083                bestDj=value;
3084                bestSequence = iSequence;
3085              } else {
3086                // just to make sure we don't exit before got something
3087                numberWanted++;
3088              }
3089            }
3090          }
3091          break;
3092        case ClpSimplex::atLowerBound:
3093          value=cost[iSequence];
3094          for (j=startColumn[iSequence];
3095               j<startColumn[iSequence]+length[iSequence];j++) {
3096            int jRow=row[j];
3097            value -= duals[jRow]*element[j];
3098          }
3099          value = -value;
3100          if (value>tolerance) {
3101            numberWanted--;
3102            if (value>bestDj) {
3103              // check flagged variable and correct dj
3104              if (!model->flagged(iSequence)) {
3105                bestDj=value;
3106                bestSequence = iSequence;
3107              } else {
3108                // just to make sure we don't exit before got something
3109                numberWanted++;
3110              }
3111            }
3112          }
3113          break;
3114        }
3115      }
3116      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
3117        // give up
3118        break;
3119      }
3120      if (!numberWanted)
3121        break;
3122    }
3123    if (bestSequence!=saveSequence) {
3124      // recompute dj
3125      double value=cost[bestSequence];
3126      for (j=startColumn[bestSequence];
3127           j<startColumn[bestSequence]+length[bestSequence];j++) {
3128        int jRow=row[j];
3129        value -= duals[jRow]*element[j];
3130      }
3131      reducedCost[bestSequence] = value;
3132      savedBestSequence_ = bestSequence;
3133      savedBestDj_ = reducedCost[savedBestSequence_];
3134    }
3135  }
3136  currentWanted_=numberWanted;
3137}
3138// Sets up an effective RHS
3139void 
3140ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
3141{
3142  delete [] rhsOffset_;
3143  int numberRows = model->numberRows();
3144  rhsOffset_ = new double[numberRows];
3145  rhsOffset(model,true);
3146}
3147// Gets rid of special copies
3148void 
3149ClpPackedMatrix::clearCopies()
3150{
3151  delete rowCopy_;
3152  delete columnCopy_;
3153  rowCopy_=NULL;
3154  columnCopy_=NULL;
3155  flags_ &= ~(4+8);
3156 }
3157// makes sure active columns correct
3158int 
3159ClpPackedMatrix::refresh(ClpSimplex * model)
3160{
3161  numberActiveColumns_ = matrix_->getNumCols();
3162#if 0
3163  ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
3164  ClpPackedMatrix* rowCopy =
3165    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3166  // Make sure it is really a ClpPackedMatrix
3167  assert (rowCopy!=NULL);
3168
3169  const int * column = rowCopy->getIndices();
3170  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3171  const int * rowLength = rowCopy->getVectorLengths();
3172  const double * element = rowCopy->getElements();
3173  for (int i=0;i<rowCopy->getNumRows();i++) {
3174    if (!rowLength[i])
3175      printf("zero row %d\n",i);
3176  }
3177  delete rowCopy;
3178#endif
3179  return 0;
3180}
3181
3182/* Scales rowCopy if column copy scaled
3183   Only called if scales already exist */
3184void 
3185ClpPackedMatrix::scaleRowCopy(ClpModel * model) const 
3186{
3187  if (model->rowCopy()) {
3188    // need to replace row by row
3189    int numberRows = model->numberRows();
3190    int numberColumns = matrix_->getNumCols();
3191    double * newElement = new double[numberColumns];
3192    ClpMatrixBase * rowCopyBase=model->rowCopy();
3193#ifndef NO_RTTI
3194    ClpPackedMatrix* rowCopy =
3195      dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3196    // Make sure it is really a ClpPackedMatrix
3197    assert (rowCopy!=NULL);
3198#else
3199    ClpPackedMatrix* rowCopy =
3200      static_cast< ClpPackedMatrix*>(rowCopyBase);
3201#endif
3202
3203    const int * column = rowCopy->getIndices();
3204    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3205    const double * element = rowCopy->getElements();
3206    const double * rowScale = model->rowScale();
3207    const double * columnScale = model->columnScale();
3208    // scale row copy
3209    for (int iRow=0;iRow<numberRows;iRow++) {
3210      CoinBigIndex j;
3211      double scale = rowScale[iRow];
3212      const double * elementsInThisRow = element + rowStart[iRow];
3213      const int * columnsInThisRow = column + rowStart[iRow];
3214      int number = rowStart[iRow+1]-rowStart[iRow];
3215      assert (number<=numberColumns);
3216      for (j=0;j<number;j++) {
3217        int iColumn = columnsInThisRow[j];
3218        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
3219      }
3220      rowCopy->replaceVector(iRow,number,newElement);
3221    }
3222    delete [] newElement;
3223  }
3224}
3225/* Realy really scales column copy
3226   Only called if scales already exist.
3227   Up to user ro delete */
3228ClpMatrixBase * 
3229ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const 
3230{
3231  // need to replace column by column
3232  int numberRows = model->numberRows();
3233  int numberColumns = matrix_->getNumCols();
3234  double * newElement = new double[numberRows];
3235  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
3236  const int * row = copy->getIndices();
3237  const CoinBigIndex * columnStart = copy->getVectorStarts();
3238  const int * length = copy->getVectorLengths();
3239  const double * element = copy->getElements();
3240  const double * rowScale = model->rowScale();
3241  const double * columnScale = model->columnScale();
3242  // scale column copy
3243  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3244    CoinBigIndex j;
3245    double scale = columnScale[iColumn];
3246    const double * elementsInThisColumn = element + columnStart[iColumn];
3247    const int * rowsInThisColumn = row + columnStart[iColumn];
3248    int number = length[iColumn];
3249    assert (number<=numberRows);
3250    for (j=0;j<number;j++) {
3251      int iRow = rowsInThisColumn[j];
3252      newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
3253    }
3254    copy->replaceVector(iColumn,number,newElement);
3255  }
3256  delete [] newElement;
3257  return copy;
3258}
3259// Really scale matrix
3260void 
3261ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
3262{
3263  clearCopies();
3264  int numberColumns = matrix_->getNumCols();
3265  const int * row = matrix_->getIndices();
3266  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3267  const int * length = matrix_->getVectorLengths();
3268  double * element = matrix_->getMutableElements();
3269  // scale
3270  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3271    CoinBigIndex j;
3272    double scale = columnScale[iColumn];
3273    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
3274      int iRow = row[j];
3275      element[j] *= scale*rowScale[iRow];
3276    }
3277  }
3278}
3279/* Delete the columns whose indices are listed in <code>indDel</code>. */
3280void 
3281ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
3282{ 
3283  if (matrix_->getNumCols())
3284    matrix_->deleteCols(numDel,indDel);
3285  clearCopies();
3286  numberActiveColumns_ = matrix_->getNumCols();
3287  // may now have gaps
3288  flags_ |= 2;
3289  matrix_->setExtraGap(1.0e-50);
3290}
3291/* Delete the rows whose indices are listed in <code>indDel</code>. */
3292void 
3293ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
3294{ 
3295  if (matrix_->getNumRows())
3296    matrix_->deleteRows(numDel,indDel);
3297  clearCopies();
3298  numberActiveColumns_ = matrix_->getNumCols();
3299  // may now have gaps
3300  flags_ |= 2;
3301  matrix_->setExtraGap(1.0e-50);
3302}
3303#ifndef CLP_NO_VECTOR
3304// Append Columns
3305void 
3306ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
3307{ 
3308  matrix_->appendCols(number,columns);
3309  numberActiveColumns_ = matrix_->getNumCols();
3310  clearCopies();
3311}
3312// Append Rows
3313void 
3314ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
3315{ 
3316  matrix_->appendRows(number,rows);
3317  numberActiveColumns_ = matrix_->getNumCols();
3318  // may now have gaps
3319  flags_ |= 2;
3320  clearCopies();
3321}
3322#endif
3323/* Set the dimensions of the matrix. In effect, append new empty
3324   columns/rows to the matrix. A negative number for either dimension
3325   means that that dimension doesn't change. Otherwise the new dimensions
3326   MUST be at least as large as the current ones otherwise an exception
3327   is thrown. */
3328void 
3329ClpPackedMatrix::setDimensions(int numrows, int numcols)
3330{
3331  matrix_->setDimensions(numrows,numcols);
3332}
3333/* Append a set of rows/columns to the end of the matrix. Returns number of errors
3334   i.e. if any of the new rows/columns contain an index that's larger than the
3335   number of columns-1/rows-1 (if numberOther>0) or duplicates
3336   If 0 then rows, 1 if columns */
3337int 
3338ClpPackedMatrix::appendMatrix(int number, int type,
3339                            const CoinBigIndex * starts, const int * index,
3340                            const double * element, int numberOther)
3341{
3342  int numberErrors=0;
3343  // make sure other dimension is big enough
3344  if (type==0) {
3345    // rows
3346    if (matrix_->isColOrdered()&&numberOther>matrix_->getNumCols())
3347      matrix_->setDimensions(-1,numberOther);
3348    numberErrors=matrix_->appendRows(number,starts,index,element,numberOther);
3349  } else {
3350    // columns
3351    if (!matrix_->isColOrdered()&&numberOther>matrix_->getNumRows())
3352      matrix_->setDimensions(numberOther,-1);
3353    numberErrors=matrix_->appendCols(number,starts,index,element,numberOther);
3354  }
3355  clearCopies();
3356  numberActiveColumns_ = matrix_->getNumCols();
3357  return numberErrors;
3358}
3359void
3360ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
3361{
3362  delete rowCopy_;
3363  rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
3364  // See if anything in it
3365  if (!rowCopy_->usefulInfo()) {
3366    delete rowCopy_;
3367    rowCopy_=NULL;
3368    flags_ &= ~4;
3369  } else {
3370    flags_ |= 4;
3371  }
3372}
3373void
3374ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
3375{
3376  delete columnCopy_;
3377  if ((flags_&8)!=0)
3378    columnCopy_ = new ClpPackedMatrix3(model,matrix_);
3379  else
3380    columnCopy_=NULL;
3381}
3382// Say we don't want special column copy
3383void 
3384ClpPackedMatrix::releaseSpecialColumnCopy()
3385{ 
3386  flags_ &= ~8;
3387  delete columnCopy_;
3388  columnCopy_=NULL;
3389}
3390// Correct sequence in and out to give true value
3391void 
3392ClpPackedMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) 
3393{
3394  if (columnCopy_) {
3395    if (sequenceIn!=-999) {
3396      if (sequenceIn!=sequenceOut) {
3397        if (sequenceIn<numberActiveColumns_)
3398          columnCopy_->swapOne(model,this,sequenceIn);
3399        if (sequenceOut<numberActiveColumns_)
3400          columnCopy_->swapOne(model,this,sequenceOut);
3401      }
3402    } else {
3403      // do all
3404      columnCopy_->sortBlocks(model);
3405    }
3406  }
3407}
3408// Check validity
3409void 
3410ClpPackedMatrix::checkFlags() const
3411{
3412  int iColumn;
3413  // get matrix data pointers
3414  //const int * row = matrix_->getIndices();
3415  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3416  const int * columnLength = matrix_->getVectorLengths(); 
3417  const double * elementByColumn = matrix_->getElements();
3418  if (!zeros()) {
3419    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3420      CoinBigIndex j;
3421      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
3422        if (!elementByColumn[j])
3423          abort();
3424      }
3425    }
3426  }
3427  if ((flags_&2)==0) {
3428    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3429      if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
3430        abort();
3431      }
3432    }
3433  }
3434}
3435//#############################################################################
3436// Constructors / Destructor / Assignment
3437//#############################################################################
3438
3439//-------------------------------------------------------------------
3440// Default Constructor
3441//-------------------------------------------------------------------
3442ClpPackedMatrix2::ClpPackedMatrix2 () 
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}
3456//-------------------------------------------------------------------
3457// Useful Constructor
3458//-------------------------------------------------------------------
3459ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy) 
3460  : numberBlocks_(0),
3461    numberRows_(0),
3462    offset_(NULL),
3463    count_(NULL),
3464    rowStart_(NULL),
3465    column_(NULL),
3466    work_(NULL)
3467{
3468#ifdef THREAD
3469  threadId_ = NULL;
3470  info_ = NULL;
3471#endif
3472  numberRows_ = rowCopy->getNumRows();
3473  if (!numberRows_)
3474    return;
3475  int numberColumns = rowCopy->getNumCols();
3476  const int * column = rowCopy->getIndices();
3477  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3478  const int * length = rowCopy->getVectorLengths();
3479  const double * element = rowCopy->getElements();
3480  int chunk = 32768; // tune
3481  //chunk=100;
3482  // tune
3483#if 0
3484  int chunkY[7]={1024,2048,4096,8192,16384,32768,65535};
3485  int its = model->maximumIterations();
3486  if (its>=1000000&&its<1000999) {
3487    its -= 1000000;
3488    its = its/10;
3489    if (its>=7) abort();
3490    chunk=chunkY[its];
3491    printf("chunk size %d\n",chunk);
3492#define cpuid(func,ax,bx,cx,dx)\
3493    __asm__ __volatile__ ("cpuid":\
3494    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
3495    unsigned int a,b,c,d;
3496    int func=0;
3497    cpuid(func,a,b,c,d);
3498    {
3499      int i;
3500      unsigned int value;
3501      value=b;
3502      for (i=0;i<4;i++) {
3503        printf("%c",(value&0xff));
3504        value = value >>8;
3505      }
3506      value=d;
3507      for (i=0;i<4;i++) {
3508        printf("%c",(value&0xff));
3509        value = value >>8;
3510      }
3511      value=c;
3512      for (i=0;i<4;i++) {
3513        printf("%c",(value&0xff));
3514        value = value >>8;
3515      }
3516      printf("\n");
3517      int maxfunc = a;
3518      if (maxfunc>10) {
3519        printf("not intel?\n");
3520        abort();
3521      }
3522      for (func=1;func<=maxfunc;func++) {
3523        cpuid(func,a,b,c,d);
3524        printf("func %d, %x %x %x %x\n",func,a,b,c,d);
3525      }
3526    }
3527#else
3528    if (numberColumns>10000||chunk==100) {
3529#endif
3530  } else {
3531    //printf("no chunk\n");
3532    return;
3533  }
3534  // Could also analyze matrix to get natural breaks
3535  numberBlocks_ = (numberColumns+chunk-1)/chunk;
3536#ifdef THREAD
3537  // Get work areas
3538  threadId_ = new pthread_t [numberBlocks_];
3539  info_ = new dualColumn0Struct[numberBlocks_];
3540#endif
3541  // Even out
3542  chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
3543  offset_ = new int[numberBlocks_+1];
3544  offset_[numberBlocks_]=numberColumns;
3545  int nRow = numberBlocks_*numberRows_;
3546  count_ = new unsigned short[nRow];
3547  memset(count_,0,nRow*sizeof(unsigned short));
3548  rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
3549  CoinBigIndex nElement = rowStart[numberRows_];
3550  rowStart_[nRow+numberRows_]=nElement;
3551  column_ = new unsigned short [nElement];
3552  // assumes int <= double
3553  int sizeWork = 6*numberBlocks_;
3554  work_ = new double[sizeWork];;
3555  int iBlock;
3556  int nZero=0;
3557  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3558    int start = iBlock*chunk;
3559    offset_[iBlock]=start;
3560    int end = start+chunk;
3561    for (int iRow=0;iRow<numberRows_;iRow++) {
3562      if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
3563        printf("not packed correctly - gaps\n");
3564        abort();
3565      }
3566      bool lastFound=false;
3567      int nFound=0;
3568      for (CoinBigIndex j = rowStart[iRow];
3569           j<rowStart[iRow]+length[iRow];j++) {
3570        int iColumn = column[j];
3571        if (iColumn>=start) {
3572          if (iColumn<end) {
3573            if (!element[j]) {
3574              printf("not packed correctly - zero element\n");
3575              abort();
3576            }
3577            column_[j]=iColumn-start;
3578            nFound++;
3579            if (lastFound) {
3580              printf("not packed correctly - out of order\n");
3581              abort();
3582            }
3583          } else {
3584            //can't find any more
3585            lastFound=true;
3586          }
3587        } 
3588      }
3589      count_[iRow*numberBlocks_+iBlock]=nFound;
3590      if (!nFound)
3591        nZero++;
3592    }
3593  }
3594  //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
3595  //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
3596}
3597
3598//-------------------------------------------------------------------
3599// Copy constructor
3600//-------------------------------------------------------------------
3601ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs) 
3602  : numberBlocks_(rhs.numberBlocks_),
3603    numberRows_(rhs.numberRows_)
3604{ 
3605  if (numberBlocks_) {
3606    offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3607    int nRow = numberBlocks_*numberRows_;
3608    count_ = CoinCopyOfArray(rhs.count_,nRow);
3609    rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3610    CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3611    column_ = CoinCopyOfArray(rhs.column_,nElement);
3612    int sizeWork = 6*numberBlocks_;
3613    work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3614#ifdef THREAD
3615    threadId_ = new pthread_t [numberBlocks_];
3616    info_ = new dualColumn0Struct[numberBlocks_];
3617#endif
3618  } else {
3619    offset_ = NULL;
3620    count_ = NULL;
3621    rowStart_ = NULL;
3622    column_ = NULL;
3623    work_ = NULL;
3624#ifdef THREAD
3625    threadId_ = NULL;
3626    info_ = NULL;
3627#endif
3628  }
3629}
3630//-------------------------------------------------------------------
3631// Destructor
3632//-------------------------------------------------------------------
3633ClpPackedMatrix2::~ClpPackedMatrix2 ()
3634{
3635  delete [] offset_;
3636  delete [] count_;
3637  delete [] rowStart_;
3638  delete [] column_;
3639  delete [] work_;
3640#ifdef THREAD
3641  delete [] threadId_;
3642  delete [] info_;
3643#endif
3644}
3645
3646//----------------------------------------------------------------
3647// Assignment operator
3648//-------------------------------------------------------------------
3649ClpPackedMatrix2 &
3650ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
3651{
3652  if (this != &rhs) {
3653    numberBlocks_ = rhs.numberBlocks_;
3654    numberRows_ = rhs.numberRows_;
3655    delete [] offset_;
3656    delete [] count_;
3657    delete [] rowStart_;
3658    delete [] column_;
3659    delete [] work_;
3660#ifdef THREAD
3661    delete [] threadId_;
3662    delete [] info_;
3663#endif
3664    if (numberBlocks_) {
3665      offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3666      int nRow = numberBlocks_*numberRows_;
3667      count_ = CoinCopyOfArray(rhs.count_,nRow);
3668      rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3669      CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3670      column_ = CoinCopyOfArray(rhs.column_,nElement);
3671      int sizeWork = 6*numberBlocks_;
3672      work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3673#ifdef THREAD
3674      threadId_ = new pthread_t [numberBlocks_];
3675      info_ = new dualColumn0Struct[numberBlocks_];
3676#endif
3677    } else {
3678      offset_ = NULL;
3679      count_ = NULL;
3680      rowStart_ = NULL;
3681      column_ = NULL;
3682      work_ = NULL;
3683#ifdef THREAD
3684      threadId_ = NULL;
3685      info_ = NULL;
3686#endif
3687    }
3688  }
3689  return *this;
3690}
3691static int dualColumn0(const ClpSimplex * model,double * spare,
3692                       int * spareIndex,const double * arrayTemp,
3693                       const int * indexTemp,int numberIn,
3694                       int offset, double acceptablePivot,double * bestPossiblePtr,
3695                       double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
3696{
3697  // do dualColumn0
3698  int i;
3699  int numberRemaining=0;
3700  double bestPossible=0.0;
3701  double upperTheta=1.0e31;
3702  double freePivot = acceptablePivot;
3703  int posFree=-1;
3704  const double * reducedCost=model->djRegion(1);
3705  double dualTolerance = model->dualTolerance();
3706  // We can also see if infeasible or pivoting on free
3707  double tentativeTheta = 1.0e25;
3708  for (i=0;i<numberIn;i++) {
3709    double alpha = arrayTemp[i];
3710    int iSequence = indexTemp[i]+offset;
3711    double oldValue;
3712    double value;
3713    bool keep;
3714   
3715    switch(model->getStatus(iSequence)) {
3716     
3717    case ClpSimplex::basic:
3718    case ClpSimplex::isFixed:
3719      break;
3720    case ClpSimplex::isFree:
3721    case ClpSimplex::superBasic:
3722      bestPossible = CoinMax(bestPossible,fabs(alpha));
3723      oldValue = reducedCost[iSequence];
3724      // If free has to be very large - should come in via dualRow
3725      if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3726        break;
3727      if (oldValue>dualTolerance) {
3728        keep = true;
3729      } else if (oldValue<-dualTolerance) {
3730        keep = true;
3731      } else {
3732        if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3733          keep = true;
3734        else
3735          keep=false;
3736      }
3737      if (keep) {
3738        // free - choose largest
3739        if (fabs(alpha)>freePivot) {
3740          freePivot=fabs(alpha);
3741          posFree = i;
3742        }
3743      }
3744      break;
3745    case ClpSimplex::atUpperBound:
3746      oldValue = reducedCost[iSequence];
3747      value = oldValue-tentativeTheta*alpha;
3748      //assert (oldValue<=dualTolerance*1.0001);
3749      if (value>dualTolerance) {
3750        bestPossible = CoinMax(bestPossible,-alpha);
3751        value = oldValue-upperTheta*alpha;
3752        if (value>dualTolerance && -alpha>=acceptablePivot) 
3753          upperTheta = (oldValue-dualTolerance)/alpha;
3754        // add to list
3755        spare[numberRemaining]=alpha;
3756        spareIndex[numberRemaining++]=iSequence;
3757      }
3758      break;
3759    case ClpSimplex::atLowerBound:
3760      oldValue = reducedCost[iSequence];
3761      value = oldValue-tentativeTheta*alpha;
3762      //assert (oldValue>=-dualTolerance*1.0001);
3763      if (value<-dualTolerance) {
3764        bestPossible = CoinMax(bestPossible,alpha);
3765        value = oldValue-upperTheta*alpha;
3766        if (value<-dualTolerance && alpha>=acceptablePivot) 
3767          upperTheta = (oldValue+dualTolerance)/alpha;
3768        // add to list
3769        spare[numberRemaining]=alpha;
3770        spareIndex[numberRemaining++]=iSequence;
3771      }
3772      break;
3773    }
3774  }
3775  *bestPossiblePtr = bestPossible;
3776  *upperThetaPtr = upperTheta;
3777  *freePivotPtr = freePivot;
3778  *posFreePtr = posFree;
3779  return numberRemaining;
3780}
3781static int doOneBlock(double * array, int * index, 
3782                      const double * pi,const CoinBigIndex * rowStart, const double * element,
3783                      const unsigned short * column, int numberInRowArray, int numberLook)
3784{
3785  int iWhich=0;
3786  int nextN=0;
3787  CoinBigIndex nextStart=0;
3788  double nextPi=0.0;
3789  for (;iWhich<numberInRowArray;iWhich++) {
3790    nextStart = rowStart[0];
3791    nextN = rowStart[numberInRowArray]-nextStart;
3792    rowStart++;
3793    if (nextN) {
3794      nextPi = pi[iWhich];
3795      break;
3796    }
3797  }
3798  int i;
3799  while (iWhich<numberInRowArray) {
3800    double value = nextPi;
3801    CoinBigIndex j=  nextStart;
3802    int n = nextN;
3803    // get next
3804    iWhich++;
3805    for (;iWhich<numberInRowArray;iWhich++) {
3806      nextStart = rowStart[0];
3807      nextN = rowStart[numberInRowArray]-nextStart;
3808      rowStart++;
3809      if (nextN) {
3810        coin_prefetch(element+nextStart);
3811        nextPi = pi[iWhich];
3812        break;
3813      }
3814    }
3815    CoinBigIndex end =j+n;
3816    //coin_prefetch(element+rowStart_[i+1]);
3817    //coin_prefetch(column_+rowStart_[i+1]);
3818    if (n<100) {
3819      if ((n&1)!=0) {
3820        unsigned int jColumn = column[j];
3821        array[jColumn] -= value*element[j];
3822        j++;
3823      }     
3824      coin_prefetch(column+nextStart);
3825      for (;j<end;j+=2) {
3826        unsigned int jColumn0 = column[j];
3827        double value0 = value*element[j];
3828        unsigned int jColumn1 = column[j+1];
3829        double value1 = value*element[j+1];
3830        array[jColumn0] -= value0;
3831        array[jColumn1] -= value1;
3832      }
3833    } else {
3834      if ((n&1)!=0) {
3835        unsigned int jColumn = column[j];
3836        array[jColumn] -= value*element[j];
3837        j++;
3838      }     
3839      if ((n&2)!=0) {
3840        unsigned int jColumn0 = column[j];
3841        double value0 = value*element[j];
3842        unsigned int jColumn1 = column[j+1];
3843        double value1 = value*element[j+1];
3844        array[jColumn0] -= value0;
3845        array[jColumn1] -= value1;
3846        j+=2;
3847      }     
3848      if ((n&4)!=0) {
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        j+=4;
3862      }
3863      //coin_prefetch(column+nextStart);
3864      for (;j<end;j+=8) {
3865        coin_prefetch(element+j+16);
3866        unsigned int jColumn0 = column[j];
3867        double value0 = value*element[j];
3868        unsigned int jColumn1 = column[j+1];
3869        double value1 = value*element[j+1];
3870        unsigned int jColumn2 = column[j+2];
3871        double value2 = value*element[j+2];
3872        unsigned int jColumn3 = column[j+3];
3873        double value3 = value*element[j+3];
3874        array[jColumn0] -= value0;
3875        array[jColumn1] -= value1;
3876        array[jColumn2] -= value2;
3877        array[jColumn3] -= value3;
3878        coin_prefetch(column+j+16);
3879        jColumn0 = column[j+4];
3880        value0 = value*element[j+4];
3881        jColumn1 = column[j+5];
3882        value1 = value*element[j+5];
3883        jColumn2 = column[j+6];
3884        value2 = value*element[j+6];
3885        jColumn3 = column[j+7];
3886        value3 = value*element[j+7];
3887        array[jColumn0] -= value0;
3888        array[jColumn1] -= value1;
3889        array[jColumn2] -= value2;
3890        array[jColumn3] -= value3;
3891      }
3892    }
3893  }
3894  // get rid of tiny values
3895  int nSmall = numberLook;
3896  int numberNonZero=0;
3897  for (i=0;i<nSmall;i++) {
3898    double value = array[i];
3899    array[i]=0.0; 
3900    if (fabs(value)>1.0e-12) {
3901      array[numberNonZero]=value;
3902      index[numberNonZero++]=i;
3903    }
3904  }
3905  for (;i<numberLook;i+=4) {
3906    double value0 = array[i+0];
3907    double value1 = array[i+1];
3908    double value2 = array[i+2];
3909    double value3 = array[i+3];
3910    array[i+0]=0.0; 
3911    array[i+1]=0.0; 
3912    array[i+2]=0.0; 
3913    array[i+3]=0.0; 
3914    if (fabs(value0)>1.0e-12) {
3915      array[numberNonZero]=value0;
3916      index[numberNonZero++]=i+0;
3917    }
3918    if (fabs(value1)>1.0e-12) {
3919      array[numberNonZero]=value1;
3920      index[numberNonZero++]=i+1;
3921    }
3922    if (fabs(value2)>1.0e-12) {
3923      array[numberNonZero]=value2;
3924      index[numberNonZero++]=i+2;
3925    }
3926    if (fabs(value3)>1.0e-12) {
3927      array[numberNonZero]=value3;
3928      index[numberNonZero++]=i+3;
3929    }
3930  }
3931  return numberNonZero;
3932}
3933#ifdef THREAD
3934static void * doOneBlockThread(void * voidInfo)
3935{
3936  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3937  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3938                                  info->rowStart,info->element,info->column,
3939                                  info->numberInRowArray,info->numberLook);
3940  return NULL;
3941}
3942static void * doOneBlockAnd0Thread(void * voidInfo)
3943{
3944  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3945  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3946                                  info->rowStart,info->element,info->column,
3947                                  info->numberInRowArray,info->numberLook);
3948  *(info->numberOutPtr) =  dualColumn0(info->model,info->spare,
3949                                    info->spareIndex,(const double *)info->arrayTemp,
3950                                    (const int *) info->indexTemp,*(info->numberInPtr),
3951                                    info->offset, info->acceptablePivot,info->bestPossiblePtr,
3952                                    info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
3953  return NULL;
3954}
3955#endif
3956/* Return <code>x * scalar * A in <code>z</code>.
3957   Note - x packed and z will be packed mode
3958   Squashes small elements and knows about ClpSimplex */
3959void 
3960ClpPackedMatrix2::transposeTimes(const ClpSimplex * model, 
3961                                 const CoinPackedMatrix * rowCopy,
3962                                 const CoinIndexedVector * rowArray,
3963                                 CoinIndexedVector * spareArray,
3964                                 CoinIndexedVector * columnArray) const
3965{
3966  // See if dualColumn0 coding wanted
3967  bool dualColumn = model->spareIntArray_[0]==1;
3968  double acceptablePivot = model->spareDoubleArray_[0];
3969  double bestPossible=0.0;
3970  double upperTheta=1.0e31;
3971  double freePivot = acceptablePivot;
3972  int posFree=-1;
3973  int numberRemaining=0;
3974  //if (model->numberIterations()>=200000) {
3975  //printf("time %g\n",CoinCpuTime()-startTime);
3976  //exit(77);
3977  //}
3978  double * pi = rowArray->denseVector();
3979  int numberNonZero=0;
3980  int * index = columnArray->getIndices();
3981  double * array = columnArray->denseVector();
3982  int numberInRowArray = rowArray->getNumElements();
3983  const int * whichRow = rowArray->getIndices();
3984  double * element = const_cast<double *>(rowCopy->getElements());
3985  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3986  int i;
3987  CoinBigIndex * rowStart2 = rowStart_;
3988  if (!dualColumn) {
3989    for (i=0;i<numberInRowArray;i++) {
3990      int iRow = whichRow[i]; 
3991      CoinBigIndex start = rowStart[iRow];
3992      *rowStart2 = start;
3993      unsigned short * count1 = count_+iRow*numberBlocks_;
3994      int put=0;
3995      for (int j=0;j<numberBlocks_;j++) {
3996        put += numberInRowArray;
3997        start += count1[j];
3998        rowStart2[put]=start;
3999      }
4000      rowStart2 ++;
4001    }
4002  } else {
4003    // also do dualColumn stuff
4004    double * spare = spareArray->denseVector();
4005    int * spareIndex = spareArray->getIndices();
4006    const double * reducedCost=model->djRegion(0);
4007    double dualTolerance = model->dualTolerance();
4008    // We can also see if infeasible or pivoting on free
4009    double tentativeTheta = 1.0e25;
4010    int addSequence=model->numberColumns();
4011    for (i=0;i<numberInRowArray;i++) {
4012      int iRow = whichRow[i]; 
4013      double alpha=pi[i];
4014      double oldValue;
4015      double value;
4016      bool keep;
4017
4018      switch(model->getStatus(iRow+addSequence)) {
4019         
4020      case ClpSimplex::basic:
4021      case ClpSimplex::isFixed:
4022        break;
4023      case ClpSimplex::isFree:
4024      case ClpSimplex::superBasic:
4025        bestPossible = CoinMax(bestPossible,fabs(alpha));
4026        oldValue = reducedCost[iRow];
4027        // If free has to be very large - should come in via dualRow
4028        if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
4029          break;
4030        if (oldValue>dualTolerance) {
4031          keep = true;
4032        } else if (oldValue<-dualTolerance) {
4033          keep = true;
4034        } else {
4035          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
4036            keep = true;
4037          else
4038            keep=false;
4039        }
4040        if (keep) {
4041          // free - choose largest
4042          if (fabs(alpha)>freePivot) {
4043            freePivot=fabs(alpha);
4044            posFree = i + addSequence;
4045          }
4046        }
4047        break;
4048      case ClpSimplex::atUpperBound:
4049        oldValue = reducedCost[iRow];
4050        value = oldValue-tentativeTheta*alpha;
4051        //assert (oldValue<=dualTolerance*1.0001);
4052        if (value>dualTolerance) {
4053          bestPossible = CoinMax(bestPossible,-alpha);
4054          value = oldValue-upperTheta*alpha;
4055          if (value>dualTolerance && -alpha>=acceptablePivot) 
4056            upperTheta = (oldValue-dualTolerance)/alpha;
4057          // add to list
4058          spare[numberRemaining]=alpha;
4059          spareIndex[numberRemaining++]=iRow+addSequence;
4060        }
4061        break;
4062      case ClpSimplex::atLowerBound:
4063        oldValue = reducedCost[iRow];
4064        value = oldValue-tentativeTheta*alpha;
4065        //assert (oldValue>=-dualTolerance*1.0001);
4066        if (value<-dualTolerance) {
4067          bestPossible = CoinMax(bestPossible,alpha);
4068          value = oldValue-upperTheta*alpha;
4069          if (value<-dualTolerance && alpha>=acceptablePivot) 
4070            upperTheta = (oldValue+dualTolerance)/alpha;
4071          // add to list
4072          spare[numberRemaining]=alpha;
4073          spareIndex[numberRemaining++]=iRow+addSequence;
4074        }
4075        break;
4076      }
4077      CoinBigIndex start = rowStart[iRow];
4078      *rowStart2 = start;
4079      unsigned short * count1 = count_+iRow*numberBlocks_;
4080      int put=0;
4081      for (int j=0;j<numberBlocks_;j++) {
4082        put += numberInRowArray;
4083        start += count1[j];
4084        rowStart2[put]=start;
4085      }
4086      rowStart2 ++;
4087    }
4088  }
4089  double * spare = spareArray->denseVector();
4090  int * spareIndex = spareArray->getIndices();
4091  int saveNumberRemaining=numberRemaining;
4092  int iBlock;
4093  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
4094    double * dwork = work_+6*iBlock;
4095    int * iwork = (int *) (dwork+3);
4096    if (!dualColumn) {
4097#ifndef THREAD
4098      int offset = offset_[iBlock];
4099      int offset3 = offset;
4100      offset = numberNonZero;
4101      double * arrayTemp = array+offset;
4102      int * indexTemp = index+offset; 
4103      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
4104                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
4105      int number = iwork[0];
4106      for (i=0;i<number;i++) {
4107        //double value = arrayTemp[i];
4108        //arrayTemp[i]=0.0;
4109        //array[numberNonZero]=value;
4110        index[numberNonZero++]=indexTemp[i]+offset3;
4111      }
4112#else
4113      int offset = offset_[iBlock];
4114      double * arrayTemp = array+offset;
4115      int * indexTemp = index+offset; 
4116      dualColumn0Struct * infoPtr = info_+iBlock;
4117      infoPtr->arrayTemp=arrayTemp;
4118      infoPtr->indexTemp=indexTemp;
4119      infoPtr->numberInPtr=&iwork[0];
4120      infoPtr->pi=pi;
4121      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
4122      infoPtr->element=element;
4123      infoPtr->column=column_;
4124      infoPtr->numberInRowArray=numberInRowArray;
4125      infoPtr->numberLook=offset_[iBlock+1]-offset;
4126      pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
4127#endif
4128    } else {
4129#ifndef THREAD
4130      int offset = offset_[iBlock];
4131      // allow for already saved
4132      int offset2 = offset + saveNumberRemaining;
4133      int offset3 = offset;
4134      offset = numberNonZero;
4135      offset2 = numberRemaining;
4136      double * arrayTemp = array+offset;
4137      int * indexTemp = index+offset; 
4138      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
4139                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
4140      iwork[1] = dualColumn0(model,spare+offset2,
4141                                   spareIndex+offset2,
4142                                   arrayTemp,indexTemp,
4143                                   iwork[0],offset3,acceptablePivot,
4144                                   &dwork[0],&dwork[1],&iwork[2],
4145                                   &dwork[2]);
4146      int number = iwork[0];
4147      int numberLook = iwork[1];
4148#if 1
4149      numberRemaining += numberLook;
4150#else
4151      double * spareTemp = spare+offset2;
4152      const int * spareIndexTemp = spareIndex+offset2; 
4153      for (i=0;i<numberLook;i++) {
4154        double value = spareTemp[i];
4155        spareTemp[i]=0.0;
4156        spare[numberRemaining] = value;
4157        spareIndex[numberRemaining++] = spareIndexTemp[i];
4158      }
4159#endif
4160      if (dwork[2]>freePivot) {
4161        freePivot=dwork[2];
4162        posFree = iwork[2]+numberNonZero;
4163      } 
4164      upperTheta =  CoinMin(dwork[1],upperTheta);
4165      bestPossible = CoinMax(dwork[0],bestPossible);
4166      for (i=0;i<number;i++) {
4167        // double value = arrayTemp[i];
4168        //arrayTemp[i]=0.0;
4169        //array[numberNonZero]=value;
4170        index[numberNonZero++]=indexTemp[i]+offset3;
4171      }
4172#else
4173      int offset = offset_[iBlock];
4174      // allow for already saved
4175      int offset2 = offset + saveNumberRemaining;
4176      double * arrayTemp = array+offset;
4177      int * indexTemp = index+offset; 
4178      dualColumn0Struct * infoPtr = info_+iBlock;
4179      infoPtr->model=model;
4180      infoPtr->spare=spare+offset2;
4181      infoPtr->spareIndex=spareIndex+offset2;
4182      infoPtr->arrayTemp=arrayTemp;
4183      infoPtr->indexTemp=indexTemp;
4184      infoPtr->numberInPtr=&iwork[0];
4185      infoPtr->offset=offset;
4186      infoPtr->acceptablePivot=acceptablePivot;
4187      infoPtr->bestPossiblePtr=&dwork[0];
4188      infoPtr->upperThetaPtr=&dwork[1];
4189      infoPtr->posFreePtr=&iwork[2];
4190      infoPtr->freePivotPtr=&dwork[2];
4191      infoPtr->numberOutPtr=&iwork[1];
4192      infoPtr->pi=pi;
4193      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
4194      infoPtr->element=element;
4195      infoPtr->column=column_;
4196      infoPtr->numberInRowArray=numberInRowArray;
4197      infoPtr->numberLook=offset_[iBlock+1]-offset;
4198      if (iBlock>=2)
4199        pthread_join(threadId_[iBlock-2],NULL);
4200      pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
4201      //pthread_join(threadId_[iBlock],NULL);
4202#endif
4203    }
4204  }
4205  for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
4206#ifdef THREAD
4207    pthread_join(threadId_[iBlock],NULL);
4208#endif
4209  }
4210#ifdef THREAD
4211  for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
4212    //pthread_join(threadId_[iBlock],NULL);
4213    int offset = offset_[iBlock];
4214    double * dwork = work_+6*iBlock;
4215    int * iwork = (int *) (dwork+3);
4216    int number = iwork[0];
4217    if (dualColumn) {
4218      // allow for already saved
4219      int offset2 = offset + saveNumberRemaining;
4220      int numberLook = iwork[1];
4221      double * spareTemp = spare+offset2;
4222      const int * spareIndexTemp = spareIndex+offset2; 
4223      for (i=0;i<numberLook;i++) {
4224        double value = spareTemp[i];
4225        spareTemp[i]=0.0;
4226        spare[numberRemaining] = value;
4227        spareIndex[numberRemaining++] = spareIndexTemp[i];
4228      }
4229      if (dwork[2]>freePivot) {
4230        freePivot=dwork[2];
4231        posFree = iwork[2]+numberNonZero;
4232      } 
4233      upperTheta =  CoinMin(dwork[1],upperTheta);
4234      bestPossible = CoinMax(dwork[0],bestPossible);
4235    }
4236    double * arrayTemp = array+offset;
4237    const int * indexTemp = index+offset; 
4238    for (i=0;i<number;i++) {
4239      double value = arrayTemp[i];
4240      arrayTemp[i]=0.0; 
4241      array[numberNonZero]=value;
4242      index[numberNonZero++]=indexTemp[i]+offset;
4243    }
4244  }
4245#endif
4246  columnArray->setNumElements(numberNonZero);
4247  columnArray->setPackedMode(true);
4248  if (dualColumn) {
4249    model->spareDoubleArray_[0]=upperTheta;
4250    model->spareDoubleArray_[1]=bestPossible;
4251    // and theta and alpha and sequence
4252    if (posFree<0) {
4253      model->spareIntArray_[1]=-1;
4254    } else {
4255      const double * reducedCost=model->djRegion(0);
4256      double alpha;
4257      int numberColumns=model->numberColumns();
4258      if (posFree<numberColumns) {
4259        alpha = columnArray->denseVector()[posFree];
4260        posFree = columnArray->getIndices()[posFree];
4261      } else {
4262        alpha = rowArray->denseVector()[posFree-numberColumns];
4263        posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
4264      }
4265      model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
4266      model->spareDoubleArray_[3]=alpha;
4267      model->spareIntArray_[1]=posFree;
4268    }
4269    spareArray->setNumElements(numberRemaining);
4270    // signal done
4271    model->spareIntArray_[0]=-1;
4272  }
4273}
4274/* Default constructor. */
4275ClpPackedMatrix3::ClpPackedMatrix3()
4276  : numberBlocks_(0),
4277    numberColumns_(0),
4278    column_(NULL),
4279    start_(NULL),
4280    row_(NULL),
4281    element_(NULL),
4282    block_(NULL)
4283{
4284}
4285/* Constructor from copy. */
4286ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy)
4287  : numberBlocks_(0),
4288    numberColumns_(0),
4289    column_(NULL),
4290    start_(NULL),
4291    row_(NULL),
4292    element_(NULL),
4293    block_(NULL)
4294{
4295#define MINBLOCK 6
4296#define MAXBLOCK 100
4297#define MAXUNROLL 10
4298  numberColumns_ = columnCopy->getNumCols();
4299  int numberRows = columnCopy->getNumRows();
4300  int * counts = new int[numberRows+1];
4301  CoinZeroN(counts,numberRows+1);
4302  CoinBigIndex nels=0;
4303  CoinBigIndex nZeroEl=0;
4304  int iColumn;
4305  // get matrix data pointers
4306  const int * row = columnCopy->getIndices();
4307  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4308  const int * columnLength = columnCopy->getVectorLengths(); 
4309  const double * elementByColumn = columnCopy->getElements();
4310  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4311    CoinBigIndex start = columnStart[iColumn];
4312    int n = columnLength[iColumn];
4313    CoinBigIndex end = start + n;
4314    int kZero=0;
4315    for (CoinBigIndex j=start;j<end;j++) {
4316      if(!elementByColumn[j])
4317        kZero++;
4318    }
4319    n -= kZero;
4320    nZeroEl += kZero;
4321    nels += n;
4322    counts[n]++;
4323  }
4324  int nZeroColumns = counts[0];
4325  counts[0]=-1;
4326  numberColumns_ -= nZeroColumns;
4327  column_ = new int [2*numberColumns_];
4328  int * lookup = column_ + numberColumns_;
4329  row_ = new int[nels];
4330  element_ = new double[nels];
4331  int nOdd=0;
4332  CoinBigIndex nInOdd=0;
4333  int i;
4334  for (i=1;i<=numberRows;i++) {
4335    int n = counts[i];
4336    if (n) {
4337      if (n<MINBLOCK||i>MAXBLOCK) {
4338        nOdd += n;
4339        counts[i]=-1;
4340        nInOdd += n*i;
4341      } else {
4342        numberBlocks_++;
4343      }
4344    } else {
4345      counts[i]=-1;
4346    }
4347  }
4348  start_ = new CoinBigIndex [nOdd+1];
4349  // even if no blocks do a dummy one
4350  numberBlocks_ = CoinMax(numberBlocks_,1);
4351  block_ = (blockStruct *) new char [numberBlocks_*sizeof(blockStruct)];
4352  memset(block_,0,numberBlocks_*sizeof(blockStruct));
4353  // Fill in what we can
4354  int nTotal=nOdd;
4355  // in case no blocks
4356  block_->startIndices_=nTotal;
4357  nels=nInOdd;
4358  int nBlock=0;
4359  for (i=0;i<=CoinMin(MAXBLOCK,numberRows);i++) {
4360    if (counts[i]>0) {
4361      blockStruct * block = block_ + nBlock;
4362      int n=counts[i];
4363      counts[i]= nBlock; // backward pointer
4364      nBlock++;
4365      block->startIndices_ = nTotal;
4366      block->startElements_ = nels;
4367      block->numberElements_ = i;
4368      // up counts
4369      nTotal += n;
4370      nels += n*i;
4371    }
4372  }
4373  // fill
4374  start_[0]=0;
4375  nOdd=0;
4376  nInOdd=0;
4377  const double * columnScale = model->columnScale();
4378  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4379    CoinBigIndex start = columnStart[iColumn];
4380    int n = columnLength[iColumn];
4381    CoinBigIndex end = start + n;
4382    int kZero=0;
4383    for (CoinBigIndex j=start;j<end;j++) {
4384      if(!elementByColumn[j])
4385        kZero++;
4386    }
4387    n -= kZero;
4388    if (n) {
4389      int iBlock = counts[n];
4390      if (iBlock>=0) {
4391        blockStruct * block = block_ + iBlock;
4392        int k = block->numberInBlock_;
4393        block->numberInBlock_ ++;
4394        column_[block->startIndices_+k]=iColumn;
4395        lookup[iColumn]=k;
4396        CoinBigIndex put = block->startElements_+k*n;
4397        for (CoinBigIndex j=start;j<end;j++) {
4398          double value = elementByColumn[j];
4399          if(value) {
4400            if (columnScale)
4401              value *= columnScale[iColumn];
4402            element_[put] = value;
4403            row_[put++]=row[j];
4404          }
4405        }
4406      } else {
4407        // odd ones
4408        for (CoinBigIndex j=start;j<end;j++) {
4409          double value = elementByColumn[j];
4410          if(value) {
4411            if (columnScale)
4412              value *= columnScale[iColumn];
4413            element_[nInOdd] = value;
4414            row_[nInOdd++]=row[j];
4415          }
4416        }
4417        column_[nOdd]=iColumn;
4418        lookup[iColumn]=-1;
4419        nOdd++;
4420        start_[nOdd]=nInOdd;
4421      }
4422    }
4423  }
4424  delete [] counts;
4425}
4426/* Destructor */
4427ClpPackedMatrix3::~ClpPackedMatrix3()
4428{
4429  delete [] column_;
4430  delete [] start_;
4431  delete [] row_;
4432  delete [] element_;
4433  delete [] block_;
4434}
4435/* The copy constructor. */
4436ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
4437  : numberBlocks_(rhs.numberBlocks_),
4438    numberColumns_(rhs.numberColumns_),
4439    column_(NULL),
4440    start_(NULL),
4441    row_(NULL),
4442    element_(NULL),
4443    block_(NULL)
4444{
4445  if (rhs.numberBlocks_) {
4446    block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
4447    column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
4448    int numberOdd = block_->startIndices_;
4449    start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
4450    blockStruct * lastBlock = block_ + (numberBlocks_-1);
4451    CoinBigIndex numberElements = lastBlock->startElements_+
4452      lastBlock->numberInBlock_*lastBlock->numberElements_;
4453    row_ = CoinCopyOfArray(rhs.row_,numberElements);
4454    element_ = CoinCopyOfArray(rhs.element_,numberElements);
4455  }
4456}
4457ClpPackedMatrix3& 
4458ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
4459{
4460  if (this != &rhs) {
4461    delete [] column_;
4462    delete [] start_;
4463    delete [] row_;
4464    delete [] element_;
4465    delete [] block_;
4466    numberBlocks_ = rhs.numberBlocks_;
4467    numberColumns_ = rhs.numberColumns_;
4468    if (rhs.numberBlocks_) {
4469      block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
4470      column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
4471      int numberOdd = block_->startIndices_;
4472      start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
4473      blockStruct * lastBlock = block_ + (numberBlocks_-1);
4474      CoinBigIndex numberElements = lastBlock->startElements_+
4475        lastBlock->numberInBlock_*lastBlock->numberElements_;
4476      row_ = CoinCopyOfArray(rhs.row_,numberElements);
4477      element_ = CoinCopyOfArray(rhs.element_,numberElements);
4478    } else {
4479      column_ = NULL;
4480      start_ = NULL;
4481      row_ = NULL;
4482      element_ = NULL;
4483      block_ = NULL;
4484    }
4485  }
4486  return *this;
4487}
4488/* Sort blocks */
4489void 
4490ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
4491{
4492  int * lookup = column_ + numberColumns_;
4493  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4494    blockStruct * block = block_ + iBlock;
4495    int numberInBlock = block->numberInBlock_; 
4496    int nel = block->numberElements_;
4497    int * row = row_ + block->startElements_;
4498    double * element = element_ + block->startElements_;
4499    int * column = column_ + block->startIndices_;
4500    int lastPrice=0;
4501    int firstNotPrice=numberInBlock-1;
4502    while (lastPrice<=firstNotPrice) {
4503      // find first basic or fixed
4504      int iColumn=numberInBlock;
4505      for (;lastPrice<=firstNotPrice;lastPrice++) {
4506        iColumn = column[lastPrice];
4507        if (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4508            model->getColumnStatus(iColumn)==ClpSimplex::isFixed)
4509          break;
4510      }
4511      // find last non basic or fixed
4512      int jColumn=-1;
4513      for (;firstNotPrice>lastPrice;firstNotPrice--) {
4514        jColumn = column[firstNotPrice];
4515        if (model->getColumnStatus(jColumn)!=ClpSimplex::basic&&
4516            model->getColumnStatus(jColumn)!=ClpSimplex::isFixed)
4517          break;
4518      }
4519      if (firstNotPrice>lastPrice) {
4520        assert (column[lastPrice]==iColumn);
4521        assert (column[firstNotPrice]==jColumn);
4522        // need to swap
4523        column[firstNotPrice]=iColumn;
4524        lookup[iColumn]=firstNotPrice;
4525        column[lastPrice]=jColumn;
4526        lookup[jColumn]=lastPrice;
4527        double * elementA = element + lastPrice*nel;
4528        int * rowA = row + lastPrice*nel;
4529        double * elementB = element + firstNotPrice*nel;
4530        int * rowB = row + firstNotPrice*nel;
4531        for (int i=0;i<nel;i++) {
4532          int temp = rowA[i];
4533          double tempE = elementA[i];
4534          rowA[i]=rowB[i];
4535          elementA[i]=elementB[i];
4536          rowB[i]=temp;
4537          elementB[i]=tempE;
4538        }
4539        firstNotPrice--;
4540        lastPrice++;
4541      } else if (lastPrice==firstNotPrice) {
4542        // make sure correct side
4543        iColumn = column[lastPrice];
4544        if (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4545            model->getColumnStatus(iColumn)!=ClpSimplex::isFixed)
4546          lastPrice++;
4547        break;
4548      }
4549    }
4550    block->numberPrice_=lastPrice;
4551#ifndef NDEBUG
4552    // check
4553    int i;
4554    for (i=0;i<lastPrice;i++) {
4555      int iColumn = column[i];
4556      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4557              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
4558      assert (lookup[iColumn]==i);
4559    }
4560    for (;i<numberInBlock;i++) {
4561      int iColumn = column[i];
4562      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4563              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4564      assert (lookup[iColumn]==i);
4565    }
4566#endif
4567  }
4568}
4569// Swap one variable
4570void 
4571ClpPackedMatrix3::swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
4572                          int iColumn)
4573{
4574  int * lookup = column_ + numberColumns_;
4575  // position in block
4576  int kA = lookup[iColumn];
4577  if (kA<0)
4578    return; // odd one
4579  // get matrix data pointers
4580  const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
4581  //const int * row = columnCopy->getIndices();
4582  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4583  const int * columnLength = columnCopy->getVectorLengths(); 
4584  const double * elementByColumn = columnCopy->getElements();
4585  CoinBigIndex start = columnStart[iColumn];
4586  int n = columnLength[iColumn];
4587  if (matrix->zeros()) {
4588    CoinBigIndex end = start + n;
4589    for (CoinBigIndex j=start;j<end;j++) {
4590      if(!elementByColumn[j])
4591        n--;
4592    }
4593  }
4594  // find block - could do binary search
4595  int iBlock=CoinMin(n,numberBlocks_)-1;
4596  while (block_[iBlock].numberElements_!=n)
4597    iBlock--;
4598  blockStruct * block = block_ + iBlock;
4599  int nel = block->numberElements_;
4600  int * row = row_ + block->startElements_;
4601  double * element = element_ + block->startElements_;
4602  int * column = column_ + block->startIndices_;
4603  assert (column[kA]==iColumn);
4604  bool moveUp = (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4605                 model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4606  int lastPrice=block->numberPrice_;
4607  int kB;
4608  if (moveUp) {
4609    kB=lastPrice-1;
4610    block->numberPrice_--;
4611  } else {
4612    kB=lastPrice;
4613    block->numberPrice_++;
4614  }
4615  int jColumn = column[kB];
4616  column[kA]=jColumn;
4617  lookup[jColumn]=kA;
4618  column[kB]=iColumn;
4619  lookup[iColumn]=kB;
4620  double * elementA = element + kB*nel;
4621  int * rowA = row + kB*nel;
4622  double * elementB = element + kA*nel;
4623  int * rowB = row + kA*nel;
4624  int i;
4625  for (i=0;i<nel;i++) {
4626    int temp = rowA[i];
4627    double tempE = elementA[i];
4628    rowA[i]=rowB[i];
4629    elementA[i]=elementB[i];
4630    rowB[i]=temp;
4631    elementB[i]=tempE;
4632  }
4633#if 1
4634#ifndef NDEBUG
4635  // check
4636  for (i=0;i<block->numberPrice_;i++) {
4637    int iColumn = column[i];
4638    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
4639      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4640              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
4641    assert (lookup[iColumn]==i);
4642  }
4643  int numberInBlock = block->numberInBlock_; 
4644  for (;i<numberInBlock;i++) {
4645    int iColumn = column[i];
4646    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
4647      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4648              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4649    assert (lookup[iColumn]==i);
4650  }
4651#endif
4652#endif
4653}
4654/* Return <code>x * -1 * A in <code>z</code>.
4655   Note - x packed and z will be packed mode
4656   Squashes small elements and knows about ClpSimplex */
4657void 
4658ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
4659                                 const double * pi,
4660                                 CoinIndexedVector * output) const
4661{
4662  int numberNonZero=0;
4663  int * index = output->getIndices();
4664  double * array = output->denseVector();
4665  double zeroTolerance = model->factorization()->zeroTolerance();
4666  double value = 0.0;
4667  CoinBigIndex j;
4668  int numberOdd = block_->startIndices_;
4669  if (numberOdd) {
4670    // A) as probably long may be worth unrolling
4671    CoinBigIndex end = start_[1];
4672    for (j=start_[0]; j<end;j++) {
4673      int iRow = row_[j];
4674      value += pi[iRow]*element_[j];
4675    }
4676    int iColumn;
4677    // int jColumn=column_[0];
4678   
4679    for (iColumn=0;iColumn<numberOdd-1;iColumn++) {
4680      CoinBigIndex start = end;
4681      end = start_[iColumn+2];
4682      if (fabs(value)>zeroTolerance) {
4683        array[numberNonZero]=value;
4684        index[numberNonZero++]=column_[iColumn];
4685        //index[numberNonZero++]=jColumn;
4686      }
4687      // jColumn = column_[iColumn+1];
4688      value = 0.0;
4689      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
4690        for (j=start; j<end;j++) {
4691          int iRow = row_[j];
4692          value += pi[iRow]*element_[j];
4693        }
4694        //}
4695    }
4696    if (fabs(value)>zeroTolerance) {
4697      array[numberNonZero]=value;
4698      index[numberNonZero++]=column_[iColumn];
4699      //index[numberNonZero++]=jColumn;
4700    }
4701  }
4702  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4703    // B) Can sort so just do nonbasic (and nonfixed)
4704    // C) Can do two at a time (if so put odd one into start_)
4705    // D) can use switch
4706    blockStruct * block = block_ + iBlock;
4707    //int numberPrice = block->numberInBlock_;
4708    int numberPrice = block->numberPrice_; 
4709    int nel = block->numberElements_;
4710    int * row = row_ + block->startElements_;
4711    double * element = element_ + block->startElements_;
4712    int * column = column_ + block->startIndices_;
4713#if 0
4714    // two at a time
4715    if ((numberPrice&1)!=0) {
4716      double value = 0.0;
4717      int nel2=nel;
4718      for (; nel2;nel2--) {
4719        int iRow = *row++;
4720        value += pi[iRow]*(*element++);
4721      }
4722      if (fabs(value)>zeroTolerance) {
4723        array[numberNonZero]=value;
4724        index[numberNonZero++]=*column;
4725      }
4726      column++;
4727    }
4728    numberPrice = numberPrice>>1;
4729    switch ((nel%2)) {
4730      // 2 k +0
4731    case 0:
4732      for (;numberPrice;numberPrice--) {
4733        double value0 = 0.0;
4734        double value1 = 0.0;
4735        int nel2=nel;
4736        for (; nel2;nel2--) {
4737          int iRow0 = *row;
4738          int iRow1 = *(row+nel);
4739          row++;
4740          double element0 = *element;
4741          double element1 = *(element+nel);
4742          element++;
4743          value0 += pi[iRow0]*element0;
4744          value1 += pi[iRow1]*element1;
4745        }
4746        row += nel;
4747        element += nel;
4748        if (fabs(value0)>zeroTolerance) {
4749          array[numberNonZero]=value0;
4750          index[numberNonZero++]=*column;
4751        }
4752        column++;
4753        if (fabs(value1)>zeroTolerance) {
4754          array[numberNonZero]=value1;
4755          index[numberNonZero++]=*column;
4756        }
4757        column++;
4758      } 
4759      break;
4760      // 2 k +1
4761    case 1:
4762      for (;numberPrice;numberPrice--) {
4763        double value0;
4764        double value1;
4765        int nel2=nel-1;
4766        {
4767          int iRow0 = row[0];
4768          int iRow1 = row[nel];
4769          double pi0 = pi[iRow0];
4770          double pi1 = pi[iRow1];
4771          value0 = pi0*element[0];
4772          value1 = pi1*element[nel];
4773          row++;
4774          element++;
4775        }
4776        for (; nel2;nel2--) {
4777          int iRow0 = *row;
4778          int iRow1 = *(row+nel);
4779          row++;
4780          double element0 = *element;
4781          double element1 = *(element+nel);
4782          element++;
4783          value0 += pi[iRow0]*element0;
4784          value1 += pi[iRow1]*element1;
4785        }
4786        row += nel;
4787        element += nel;
4788        if (fabs(value0)>zeroTolerance) {
4789          array[numberNonZero]=value0;
4790          index[numberNonZero++]=*column;
4791        }
4792        column++;
4793        if (fabs(value1)>zeroTolerance) {
4794          array[numberNonZero]=value1;
4795          index[numberNonZero++]=*column;
4796        }
4797        column++;
4798      }
4799      break;
4800    }
4801#else
4802    for (;numberPrice;numberPrice--) {
4803      double value = 0.0;
4804      int nel2=nel;
4805      for (; nel2;nel2--) {
4806        int iRow = *row++;
4807        value += pi[iRow]*(*element++);
4808      }
4809      if (fabs(value)>zeroTolerance) {
4810        array[numberNonZero]=value;
4811        index[numberNonZero++]=*column;
4812      }
4813      column++;
4814    }
4815#endif
4816  }
4817  output->setNumElements(numberNonZero);
4818}
4819// Updates two arrays for steepest
4820void 
4821ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
4822                                  const double * pi, CoinIndexedVector * output,
4823                                  const double * piWeight,
4824                                  double referenceIn, double devex,
4825                                  // Array for exact devex to say what is in reference framework
4826                                  unsigned int * reference,
4827                                  double * weights, double scaleFactor)
4828{
4829  int numberNonZero=0;
4830  int * index = output->getIndices();
4831  double * array = output->denseVector();
4832  double zeroTolerance = model->factorization()->zeroTolerance();
4833  double value = 0.0;
4834  bool killDjs = (scaleFactor==0.0);
4835  if (!scaleFactor)
4836    scaleFactor=1.0;
4837  int numberOdd = block_->startIndices_;
4838  int iColumn;
4839  CoinBigIndex end = start_[0];
4840  for (iColumn=0;iColumn<numberOdd;iColumn++) {
4841    CoinBigIndex start = end;
4842    CoinBigIndex j;
4843    int jColumn = column_[iColumn];
4844    end = start_[iColumn+1];
4845    value = 0.0;
4846    if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
4847      for (j=start; j<end;j++) {
4848        int iRow = row_[j];
4849        value -= pi[iRow]*element_[j];
4850      }
4851      if (fabs(value)>zeroTolerance) {
4852        // and do other array
4853        double modification = 0.0;
4854        for (j=start; j<end;j++) {
4855          int iRow = row_[j];
4856          modification += piWeight[iRow]*element_[j];
4857        }
4858        double thisWeight = weights[jColumn];
4859        double pivot = value*scaleFactor;
4860        double pivotSquared = pivot * pivot;
4861        thisWeight += pivotSquared * devex + pivot * modification;
4862        if (thisWeight<DEVEX_TRY_NORM) {
4863          if (referenceIn<0.0) {
4864            // steepest
4865            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
4866          } else {
4867            // exact
4868            thisWeight = referenceIn*pivotSquared;
4869            if (reference(jColumn))
4870              thisWeight += 1.0;
4871            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
4872          }
4873        }
4874        weights[jColumn] = thisWeight;
4875        if (!killDjs) {
4876          array[numberNonZero]=value;
4877          index[numberNonZero++]=jColumn;
4878        }
4879      }
4880    }
4881  }
4882  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4883    // B) Can sort so just do nonbasic (and nonfixed)
4884    // C) Can do two at a time (if so put odd one into start_)
4885    // D) can use switch
4886    blockStruct * block = block_ + iBlock;
4887    //int numberPrice = block->numberInBlock_;
4888    int numberPrice = block->numberPrice_; 
4889    int nel = block->numberElements_;
4890    int * row = row_ + block->startElements_;
4891    double * element = element_ + block->startElements_;
4892    int * column = column_ + block->startIndices_;
4893    for (;numberPrice;numberPrice--) {
4894      double value = 0.0;
4895      int nel2=nel;
4896      for (; nel2;nel2--) {
4897        int iRow = *row++;
4898        value -= pi[iRow]*(*element++);
4899      }
4900      if (fabs(value)>zeroTolerance) {
4901        int jColumn = *column;
4902        // back to beginning
4903        row -= nel;
4904        element -= nel;
4905        // and do other array
4906        double modification = 0.0;
4907        nel2=nel;
4908        for (; nel2;nel2--) {
4909          int iRow = *row++;
4910          modification += piWeight[iRow]*(*element++);
4911        }
4912        double thisWeight = weights[jColumn];
4913        double pivot = value*scaleFactor;
4914        double pivotSquared = pivot * pivot;
4915        thisWeight += pivotSquared * devex + pivot * modification;
4916        if (thisWeight<DEVEX_TRY_NORM) {
4917          if (referenceIn<0.0) {
4918            // steepest
4919            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
4920          } else {
4921            // exact
4922            thisWeight = referenceIn*pivotSquared;
4923            if (reference(jColumn))
4924              thisWeight += 1.0;
4925            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
4926          }
4927        }
4928        weights[jColumn] = thisWeight;
4929        if (!killDjs) {
4930          array[numberNonZero]=value;
4931          index[numberNonZero++]=jColumn;
4932        }
4933      }
4934      column++;
4935    }
4936  }
4937  output->setNumElements(numberNonZero);
4938  output->setPackedMode(true);
4939}
4940#ifdef CLP_ALL_ONE_FILE
4941#undef reference
4942#endif
Note: See TracBrowser for help on using the repository browser.