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

Last change on this file since 928 was 928, checked in by forrest, 14 years ago

active columns

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