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

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

for messages

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 147.1 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  return numberErrors;
3330}
3331void
3332ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
3333{
3334  delete rowCopy_;
3335  rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
3336  // See if anything in it
3337  if (!rowCopy_->usefulInfo()) {
3338    delete rowCopy_;
3339    rowCopy_=NULL;
3340    flags_ &= ~4;
3341  } else {
3342    flags_ |= 4;
3343  }
3344}
3345void
3346ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
3347{
3348  delete columnCopy_;
3349  if ((flags_&8)!=0)
3350    columnCopy_ = new ClpPackedMatrix3(model,matrix_);
3351  else
3352    columnCopy_=NULL;
3353}
3354// Say we don't want special column copy
3355void 
3356ClpPackedMatrix::releaseSpecialColumnCopy()
3357{ 
3358  flags_ &= ~8;
3359  delete columnCopy_;
3360  columnCopy_=NULL;
3361}
3362// Correct sequence in and out to give true value
3363void 
3364ClpPackedMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) 
3365{
3366  if (columnCopy_) {
3367    if (sequenceIn!=-999) {
3368      if (sequenceIn!=sequenceOut) {
3369        if (sequenceIn<numberActiveColumns_)
3370          columnCopy_->swapOne(model,this,sequenceIn);
3371        if (sequenceOut<numberActiveColumns_)
3372          columnCopy_->swapOne(model,this,sequenceOut);
3373      }
3374    } else {
3375      // do all
3376      columnCopy_->sortBlocks(model);
3377    }
3378  }
3379}
3380// Check validity
3381void 
3382ClpPackedMatrix::checkFlags() const
3383{
3384  int iColumn;
3385  // get matrix data pointers
3386  //const int * row = matrix_->getIndices();
3387  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3388  const int * columnLength = matrix_->getVectorLengths(); 
3389  const double * elementByColumn = matrix_->getElements();
3390  if (!zeros()) {
3391    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3392      CoinBigIndex j;
3393      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
3394        if (!elementByColumn[j])
3395          abort();
3396      }
3397    }
3398  }
3399  if ((flags_&2)==0) {
3400    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3401      if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
3402        abort();
3403      }
3404    }
3405  }
3406}
3407//#############################################################################
3408// Constructors / Destructor / Assignment
3409//#############################################################################
3410
3411//-------------------------------------------------------------------
3412// Default Constructor
3413//-------------------------------------------------------------------
3414ClpPackedMatrix2::ClpPackedMatrix2 () 
3415  : numberBlocks_(0),
3416    numberRows_(0),
3417    offset_(NULL),
3418    count_(NULL),
3419    rowStart_(NULL),
3420    column_(NULL),
3421    work_(NULL)
3422{
3423#ifdef THREAD
3424  threadId_ = NULL;
3425  info_ = NULL;
3426#endif
3427}
3428//-------------------------------------------------------------------
3429// Useful Constructor
3430//-------------------------------------------------------------------
3431ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy) 
3432  : numberBlocks_(0),
3433    numberRows_(0),
3434    offset_(NULL),
3435    count_(NULL),
3436    rowStart_(NULL),
3437    column_(NULL),
3438    work_(NULL)
3439{
3440#ifdef THREAD
3441  threadId_ = NULL;
3442  info_ = NULL;
3443#endif
3444  numberRows_ = rowCopy->getNumRows();
3445  if (!numberRows_)
3446    return;
3447  int numberColumns = rowCopy->getNumCols();
3448  const int * column = rowCopy->getIndices();
3449  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3450  const int * length = rowCopy->getVectorLengths();
3451  const double * element = rowCopy->getElements();
3452  int chunk = 32768; // tune
3453  //chunk=100;
3454  // tune
3455#if 0
3456  int chunkY[7]={1024,2048,4096,8192,16384,32768,65535};
3457  int its = model->maximumIterations();
3458  if (its>=1000000&&its<1000999) {
3459    its -= 1000000;
3460    its = its/10;
3461    if (its>=7) abort();
3462    chunk=chunkY[its];
3463    printf("chunk size %d\n",chunk);
3464#define cpuid(func,ax,bx,cx,dx)\
3465    __asm__ __volatile__ ("cpuid":\
3466    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
3467    unsigned int a,b,c,d;
3468    int func=0;
3469    cpuid(func,a,b,c,d);
3470    {
3471      int i;
3472      unsigned int value;
3473      value=b;
3474      for (i=0;i<4;i++) {
3475        printf("%c",(value&0xff));
3476        value = value >>8;
3477      }
3478      value=d;
3479      for (i=0;i<4;i++) {
3480        printf("%c",(value&0xff));
3481        value = value >>8;
3482      }
3483      value=c;
3484      for (i=0;i<4;i++) {
3485        printf("%c",(value&0xff));
3486        value = value >>8;
3487      }
3488      printf("\n");
3489      int maxfunc = a;
3490      if (maxfunc>10) {
3491        printf("not intel?\n");
3492        abort();
3493      }
3494      for (func=1;func<=maxfunc;func++) {
3495        cpuid(func,a,b,c,d);
3496        printf("func %d, %x %x %x %x\n",func,a,b,c,d);
3497      }
3498    }
3499#else
3500    if (numberColumns>10000||chunk==100) {
3501#endif
3502  } else {
3503    //printf("no chunk\n");
3504    return;
3505  }
3506  // Could also analyze matrix to get natural breaks
3507  numberBlocks_ = (numberColumns+chunk-1)/chunk;
3508#ifdef THREAD
3509  // Get work areas
3510  threadId_ = new pthread_t [numberBlocks_];
3511  info_ = new dualColumn0Struct[numberBlocks_];
3512#endif
3513  // Even out
3514  chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
3515  offset_ = new int[numberBlocks_+1];
3516  offset_[numberBlocks_]=numberColumns;
3517  int nRow = numberBlocks_*numberRows_;
3518  count_ = new unsigned short[nRow];
3519  memset(count_,0,nRow*sizeof(unsigned short));
3520  rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
3521  CoinBigIndex nElement = rowStart[numberRows_];
3522  rowStart_[nRow+numberRows_]=nElement;
3523  column_ = new unsigned short [nElement];
3524  // assumes int <= double
3525  int sizeWork = 6*numberBlocks_;
3526  work_ = new double[sizeWork];;
3527  int iBlock;
3528  int nZero=0;
3529  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3530    int start = iBlock*chunk;
3531    offset_[iBlock]=start;
3532    int end = start+chunk;
3533    for (int iRow=0;iRow<numberRows_;iRow++) {
3534      if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
3535        printf("not packed correctly - gaps\n");
3536        abort();
3537      }
3538      bool lastFound=false;
3539      int nFound=0;
3540      for (CoinBigIndex j = rowStart[iRow];
3541           j<rowStart[iRow]+length[iRow];j++) {
3542        int iColumn = column[j];
3543        if (iColumn>=start) {
3544          if (iColumn<end) {
3545            if (!element[j]) {
3546              printf("not packed correctly - zero element\n");
3547              abort();
3548            }
3549            column_[j]=iColumn-start;
3550            nFound++;
3551            if (lastFound) {
3552              printf("not packed correctly - out of order\n");
3553              abort();
3554            }
3555          } else {
3556            //can't find any more
3557            lastFound=true;
3558          }
3559        } 
3560      }
3561      count_[iRow*numberBlocks_+iBlock]=nFound;
3562      if (!nFound)
3563        nZero++;
3564    }
3565  }
3566  //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
3567  //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
3568}
3569
3570//-------------------------------------------------------------------
3571// Copy constructor
3572//-------------------------------------------------------------------
3573ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs) 
3574  : numberBlocks_(rhs.numberBlocks_),
3575    numberRows_(rhs.numberRows_)
3576{ 
3577  if (numberBlocks_) {
3578    offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3579    int nRow = numberBlocks_*numberRows_;
3580    count_ = CoinCopyOfArray(rhs.count_,nRow);
3581    rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3582    CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3583    column_ = CoinCopyOfArray(rhs.column_,nElement);
3584    int sizeWork = 6*numberBlocks_;
3585    work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3586#ifdef THREAD
3587    threadId_ = new pthread_t [numberBlocks_];
3588    info_ = new dualColumn0Struct[numberBlocks_];
3589#endif
3590  } else {
3591    offset_ = NULL;
3592    count_ = NULL;
3593    rowStart_ = NULL;
3594    column_ = NULL;
3595    work_ = NULL;
3596#ifdef THREAD
3597    threadId_ = NULL;
3598    info_ = NULL;
3599#endif
3600  }
3601}
3602//-------------------------------------------------------------------
3603// Destructor
3604//-------------------------------------------------------------------
3605ClpPackedMatrix2::~ClpPackedMatrix2 ()
3606{
3607  delete [] offset_;
3608  delete [] count_;
3609  delete [] rowStart_;
3610  delete [] column_;
3611  delete [] work_;
3612#ifdef THREAD
3613  delete [] threadId_;
3614  delete [] info_;
3615#endif
3616}
3617
3618//----------------------------------------------------------------
3619// Assignment operator
3620//-------------------------------------------------------------------
3621ClpPackedMatrix2 &
3622ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
3623{
3624  if (this != &rhs) {
3625    numberBlocks_ = rhs.numberBlocks_;
3626    numberRows_ = rhs.numberRows_;
3627    delete [] offset_;
3628    delete [] count_;
3629    delete [] rowStart_;
3630    delete [] column_;
3631    delete [] work_;
3632#ifdef THREAD
3633    delete [] threadId_;
3634    delete [] info_;
3635#endif
3636    if (numberBlocks_) {
3637      offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3638      int nRow = numberBlocks_*numberRows_;
3639      count_ = CoinCopyOfArray(rhs.count_,nRow);
3640      rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3641      CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3642      column_ = CoinCopyOfArray(rhs.column_,nElement);
3643      int sizeWork = 6*numberBlocks_;
3644      work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3645#ifdef THREAD
3646      threadId_ = new pthread_t [numberBlocks_];
3647      info_ = new dualColumn0Struct[numberBlocks_];
3648#endif
3649    } else {
3650      offset_ = NULL;
3651      count_ = NULL;
3652      rowStart_ = NULL;
3653      column_ = NULL;
3654      work_ = NULL;
3655#ifdef THREAD
3656      threadId_ = NULL;
3657      info_ = NULL;
3658#endif
3659    }
3660  }
3661  return *this;
3662}
3663static int dualColumn0(const ClpSimplex * model,double * spare,
3664                       int * spareIndex,const double * arrayTemp,
3665                       const int * indexTemp,int numberIn,
3666                       int offset, double acceptablePivot,double * bestPossiblePtr,
3667                       double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
3668{
3669  // do dualColumn0
3670  int i;
3671  int numberRemaining=0;
3672  double bestPossible=0.0;
3673  double upperTheta=1.0e31;
3674  double freePivot = acceptablePivot;
3675  int posFree=-1;
3676  const double * reducedCost=model->djRegion(1);
3677  double dualTolerance = model->dualTolerance();
3678  // We can also see if infeasible or pivoting on free
3679  double tentativeTheta = 1.0e25;
3680  for (i=0;i<numberIn;i++) {
3681    double alpha = arrayTemp[i];
3682    int iSequence = indexTemp[i]+offset;
3683    double oldValue;
3684    double value;
3685    bool keep;
3686   
3687    switch(model->getStatus(iSequence)) {
3688     
3689    case ClpSimplex::basic:
3690    case ClpSimplex::isFixed:
3691      break;
3692    case ClpSimplex::isFree:
3693    case ClpSimplex::superBasic:
3694      bestPossible = CoinMax(bestPossible,fabs(alpha));
3695      oldValue = reducedCost[iSequence];
3696      // If free has to be very large - should come in via dualRow
3697      if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3698        break;
3699      if (oldValue>dualTolerance) {
3700        keep = true;
3701      } else if (oldValue<-dualTolerance) {
3702        keep = true;
3703      } else {
3704        if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3705          keep = true;
3706        else
3707          keep=false;
3708      }
3709      if (keep) {
3710        // free - choose largest
3711        if (fabs(alpha)>freePivot) {
3712          freePivot=fabs(alpha);
3713          posFree = i;
3714        }
3715      }
3716      break;
3717    case ClpSimplex::atUpperBound:
3718      oldValue = reducedCost[iSequence];
3719      value = oldValue-tentativeTheta*alpha;
3720      //assert (oldValue<=dualTolerance*1.0001);
3721      if (value>dualTolerance) {
3722        bestPossible = CoinMax(bestPossible,-alpha);
3723        value = oldValue-upperTheta*alpha;
3724        if (value>dualTolerance && -alpha>=acceptablePivot) 
3725          upperTheta = (oldValue-dualTolerance)/alpha;
3726        // add to list
3727        spare[numberRemaining]=alpha;
3728        spareIndex[numberRemaining++]=iSequence;
3729      }
3730      break;
3731    case ClpSimplex::atLowerBound:
3732      oldValue = reducedCost[iSequence];
3733      value = oldValue-tentativeTheta*alpha;
3734      //assert (oldValue>=-dualTolerance*1.0001);
3735      if (value<-dualTolerance) {
3736        bestPossible = CoinMax(bestPossible,alpha);
3737        value = oldValue-upperTheta*alpha;
3738        if (value<-dualTolerance && alpha>=acceptablePivot) 
3739          upperTheta = (oldValue+dualTolerance)/alpha;
3740        // add to list
3741        spare[numberRemaining]=alpha;
3742        spareIndex[numberRemaining++]=iSequence;
3743      }
3744      break;
3745    }
3746  }
3747  *bestPossiblePtr = bestPossible;
3748  *upperThetaPtr = upperTheta;
3749  *freePivotPtr = freePivot;
3750  *posFreePtr = posFree;
3751  return numberRemaining;
3752}
3753static int doOneBlock(double * array, int * index, 
3754                      const double * pi,const CoinBigIndex * rowStart, const double * element,
3755                      const unsigned short * column, int numberInRowArray, int numberLook)
3756{
3757  int iWhich=0;
3758  int nextN=0;
3759  CoinBigIndex nextStart=0;
3760  double nextPi=0.0;
3761  for (;iWhich<numberInRowArray;iWhich++) {
3762    nextStart = rowStart[0];
3763    nextN = rowStart[numberInRowArray]-nextStart;
3764    rowStart++;
3765    if (nextN) {
3766      nextPi = pi[iWhich];
3767      break;
3768    }
3769  }
3770  int i;
3771  while (iWhich<numberInRowArray) {
3772    double value = nextPi;
3773    CoinBigIndex j=  nextStart;
3774    int n = nextN;
3775    // get next
3776    iWhich++;
3777    for (;iWhich<numberInRowArray;iWhich++) {
3778      nextStart = rowStart[0];
3779      nextN = rowStart[numberInRowArray]-nextStart;
3780      rowStart++;
3781      if (nextN) {
3782        coin_prefetch(element+nextStart);
3783        nextPi = pi[iWhich];
3784        break;
3785      }
3786    }
3787    CoinBigIndex end =j+n;
3788    //coin_prefetch(element+rowStart_[i+1]);
3789    //coin_prefetch(column_+rowStart_[i+1]);
3790    if (n<100) {
3791      if ((n&1)!=0) {
3792        unsigned int jColumn = column[j];
3793        array[jColumn] -= value*element[j];
3794        j++;
3795      }     
3796      coin_prefetch(column+nextStart);
3797      for (;j<end;j+=2) {
3798        unsigned int jColumn0 = column[j];
3799        double value0 = value*element[j];
3800        unsigned int jColumn1 = column[j+1];
3801        double value1 = value*element[j+1];
3802        array[jColumn0] -= value0;
3803        array[jColumn1] -= value1;
3804      }
3805    } else {
3806      if ((n&1)!=0) {
3807        unsigned int jColumn = column[j];
3808        array[jColumn] -= value*element[j];
3809        j++;
3810      }     
3811      if ((n&2)!=0) {
3812        unsigned int jColumn0 = column[j];
3813        double value0 = value*element[j];
3814        unsigned int jColumn1 = column[j+1];
3815        double value1 = value*element[j+1];
3816        array[jColumn0] -= value0;
3817        array[jColumn1] -= value1;
3818        j+=2;
3819      }     
3820      if ((n&4)!=0) {
3821        unsigned int jColumn0 = column[j];
3822        double value0 = value*element[j];
3823        unsigned int jColumn1 = column[j+1];
3824        double value1 = value*element[j+1];
3825        unsigned int jColumn2 = column[j+2];
3826        double value2 = value*element[j+2];
3827        unsigned int jColumn3 = column[j+3];
3828        double value3 = value*element[j+3];
3829        array[jColumn0] -= value0;
3830        array[jColumn1] -= value1;
3831        array[jColumn2] -= value2;
3832        array[jColumn3] -= value3;
3833        j+=4;
3834      }
3835      //coin_prefetch(column+nextStart);
3836      for (;j<end;j+=8) {
3837        coin_prefetch(element+j+16);
3838        unsigned int jColumn0 = column[j];
3839        double value0 = value*element[j];
3840        unsigned int jColumn1 = column[j+1];
3841        double value1 = value*element[j+1];
3842        unsigned int jColumn2 = column[j+2];
3843        double value2 = value*element[j+2];
3844        unsigned int jColumn3 = column[j+3];
3845        double value3 = value*element[j+3];
3846        array[jColumn0] -= value0;
3847        array[jColumn1] -= value1;
3848        array[jColumn2] -= value2;
3849        array[jColumn3] -= value3;
3850        coin_prefetch(column+j+16);
3851        jColumn0 = column[j+4];
3852        value0 = value*element[j+4];
3853        jColumn1 = column[j+5];
3854        value1 = value*element[j+5];
3855        jColumn2 = column[j+6];
3856        value2 = value*element[j+6];
3857        jColumn3 = column[j+7];
3858        value3 = value*element[j+7];
3859        array[jColumn0] -= value0;
3860        array[jColumn1] -= value1;
3861        array[jColumn2] -= value2;
3862        array[jColumn3] -= value3;
3863      }
3864    }
3865  }
3866  // get rid of tiny values
3867  int nSmall = numberLook;
3868  int numberNonZero=0;
3869  for (i=0;i<nSmall;i++) {
3870    double value = array[i];
3871    array[i]=0.0; 
3872    if (fabs(value)>1.0e-12) {
3873      array[numberNonZero]=value;
3874      index[numberNonZero++]=i;
3875    }
3876  }
3877  for (;i<numberLook;i+=4) {
3878    double value0 = array[i+0];
3879    double value1 = array[i+1];
3880    double value2 = array[i+2];
3881    double value3 = array[i+3];
3882    array[i+0]=0.0; 
3883    array[i+1]=0.0; 
3884    array[i+2]=0.0; 
3885    array[i+3]=0.0; 
3886    if (fabs(value0)>1.0e-12) {
3887      array[numberNonZero]=value0;
3888      index[numberNonZero++]=i+0;
3889    }
3890    if (fabs(value1)>1.0e-12) {
3891      array[numberNonZero]=value1;
3892      index[numberNonZero++]=i+1;
3893    }
3894    if (fabs(value2)>1.0e-12) {
3895      array[numberNonZero]=value2;
3896      index[numberNonZero++]=i+2;
3897    }
3898    if (fabs(value3)>1.0e-12) {
3899      array[numberNonZero]=value3;
3900      index[numberNonZero++]=i+3;
3901    }
3902  }
3903  return numberNonZero;
3904}
3905#ifdef THREAD
3906static void * doOneBlockThread(void * voidInfo)
3907{
3908  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3909  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3910                                  info->rowStart,info->element,info->column,
3911                                  info->numberInRowArray,info->numberLook);
3912  return NULL;
3913}
3914static void * doOneBlockAnd0Thread(void * voidInfo)
3915{
3916  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
3917  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
3918                                  info->rowStart,info->element,info->column,
3919                                  info->numberInRowArray,info->numberLook);
3920  *(info->numberOutPtr) =  dualColumn0(info->model,info->spare,
3921                                    info->spareIndex,(const double *)info->arrayTemp,
3922                                    (const int *) info->indexTemp,*(info->numberInPtr),
3923                                    info->offset, info->acceptablePivot,info->bestPossiblePtr,
3924                                    info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
3925  return NULL;
3926}
3927#endif
3928/* Return <code>x * scalar * A in <code>z</code>.
3929   Note - x packed and z will be packed mode
3930   Squashes small elements and knows about ClpSimplex */
3931void 
3932ClpPackedMatrix2::transposeTimes(const ClpSimplex * model, 
3933                                 const CoinPackedMatrix * rowCopy,
3934                                 const CoinIndexedVector * rowArray,
3935                                 CoinIndexedVector * spareArray,
3936                                 CoinIndexedVector * columnArray) const
3937{
3938  // See if dualColumn0 coding wanted
3939  bool dualColumn = model->spareIntArray_[0]==1;
3940  double acceptablePivot = model->spareDoubleArray_[0];
3941  double bestPossible=0.0;
3942  double upperTheta=1.0e31;
3943  double freePivot = acceptablePivot;
3944  int posFree=-1;
3945  int numberRemaining=0;
3946  //if (model->numberIterations()>=200000) {
3947  //printf("time %g\n",CoinCpuTime()-startTime);
3948  //exit(77);
3949  //}
3950  double * pi = rowArray->denseVector();
3951  int numberNonZero=0;
3952  int * index = columnArray->getIndices();
3953  double * array = columnArray->denseVector();
3954  int numberInRowArray = rowArray->getNumElements();
3955  const int * whichRow = rowArray->getIndices();
3956  double * element = const_cast<double *>(rowCopy->getElements());
3957  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3958  int i;
3959  CoinBigIndex * rowStart2 = rowStart_;
3960  if (!dualColumn) {
3961    for (i=0;i<numberInRowArray;i++) {
3962      int iRow = whichRow[i]; 
3963      CoinBigIndex start = rowStart[iRow];
3964      *rowStart2 = start;
3965      unsigned short * count1 = count_+iRow*numberBlocks_;
3966      int put=0;
3967      for (int j=0;j<numberBlocks_;j++) {
3968        put += numberInRowArray;
3969        start += count1[j];
3970        rowStart2[put]=start;
3971      }
3972      rowStart2 ++;
3973    }
3974  } else {
3975    // also do dualColumn stuff
3976    double * spare = spareArray->denseVector();
3977    int * spareIndex = spareArray->getIndices();
3978    const double * reducedCost=model->djRegion(0);
3979    double dualTolerance = model->dualTolerance();
3980    // We can also see if infeasible or pivoting on free
3981    double tentativeTheta = 1.0e25;
3982    int addSequence=model->numberColumns();
3983    for (i=0;i<numberInRowArray;i++) {
3984      int iRow = whichRow[i]; 
3985      double alpha=pi[i];
3986      double oldValue;
3987      double value;
3988      bool keep;
3989
3990      switch(model->getStatus(iRow+addSequence)) {
3991         
3992      case ClpSimplex::basic:
3993      case ClpSimplex::isFixed:
3994        break;
3995      case ClpSimplex::isFree:
3996      case ClpSimplex::superBasic:
3997        bestPossible = CoinMax(bestPossible,fabs(alpha));
3998        oldValue = reducedCost[iRow];
3999        // If free has to be very large - should come in via dualRow
4000        if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
4001          break;
4002        if (oldValue>dualTolerance) {
4003          keep = true;
4004        } else if (oldValue<-dualTolerance) {
4005          keep = true;
4006        } else {
4007          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
4008            keep = true;
4009          else
4010            keep=false;
4011        }
4012        if (keep) {
4013          // free - choose largest
4014          if (fabs(alpha)>freePivot) {
4015            freePivot=fabs(alpha);
4016            posFree = i + addSequence;
4017          }
4018        }
4019        break;
4020      case ClpSimplex::atUpperBound:
4021        oldValue = reducedCost[iRow];
4022        value = oldValue-tentativeTheta*alpha;
4023        //assert (oldValue<=dualTolerance*1.0001);
4024        if (value>dualTolerance) {
4025          bestPossible = CoinMax(bestPossible,-alpha);
4026          value = oldValue-upperTheta*alpha;
4027          if (value>dualTolerance && -alpha>=acceptablePivot) 
4028            upperTheta = (oldValue-dualTolerance)/alpha;
4029          // add to list
4030          spare[numberRemaining]=alpha;
4031          spareIndex[numberRemaining++]=iRow+addSequence;
4032        }
4033        break;
4034      case ClpSimplex::atLowerBound:
4035        oldValue = reducedCost[iRow];
4036        value = oldValue-tentativeTheta*alpha;
4037        //assert (oldValue>=-dualTolerance*1.0001);
4038        if (value<-dualTolerance) {
4039          bestPossible = CoinMax(bestPossible,alpha);
4040          value = oldValue-upperTheta*alpha;
4041          if (value<-dualTolerance && alpha>=acceptablePivot) 
4042            upperTheta = (oldValue+dualTolerance)/alpha;
4043          // add to list
4044          spare[numberRemaining]=alpha;
4045          spareIndex[numberRemaining++]=iRow+addSequence;
4046        }
4047        break;
4048      }
4049      CoinBigIndex start = rowStart[iRow];
4050      *rowStart2 = start;
4051      unsigned short * count1 = count_+iRow*numberBlocks_;
4052      int put=0;
4053      for (int j=0;j<numberBlocks_;j++) {
4054        put += numberInRowArray;
4055        start += count1[j];
4056        rowStart2[put]=start;
4057      }
4058      rowStart2 ++;
4059    }
4060  }
4061  double * spare = spareArray->denseVector();
4062  int * spareIndex = spareArray->getIndices();
4063  int saveNumberRemaining=numberRemaining;
4064  int iBlock;
4065  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
4066    double * dwork = work_+6*iBlock;
4067    int * iwork = (int *) (dwork+3);
4068    if (!dualColumn) {
4069#ifndef THREAD
4070      int offset = offset_[iBlock];
4071      int offset3 = offset;
4072      offset = numberNonZero;
4073      double * arrayTemp = array+offset;
4074      int * indexTemp = index+offset; 
4075      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
4076                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
4077      int number = iwork[0];
4078      for (i=0;i<number;i++) {
4079        //double value = arrayTemp[i];
4080        //arrayTemp[i]=0.0;
4081        //array[numberNonZero]=value;
4082        index[numberNonZero++]=indexTemp[i]+offset3;
4083      }
4084#else
4085      int offset = offset_[iBlock];
4086      double * arrayTemp = array+offset;
4087      int * indexTemp = index+offset; 
4088      dualColumn0Struct * infoPtr = info_+iBlock;
4089      infoPtr->arrayTemp=arrayTemp;
4090      infoPtr->indexTemp=indexTemp;
4091      infoPtr->numberInPtr=&iwork[0];
4092      infoPtr->pi=pi;
4093      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
4094      infoPtr->element=element;
4095      infoPtr->column=column_;
4096      infoPtr->numberInRowArray=numberInRowArray;
4097      infoPtr->numberLook=offset_[iBlock+1]-offset;
4098      pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
4099#endif
4100    } else {
4101#ifndef THREAD
4102      int offset = offset_[iBlock];
4103      // allow for already saved
4104      int offset2 = offset + saveNumberRemaining;
4105      int offset3 = offset;
4106      offset = numberNonZero;
4107      offset2 = numberRemaining;
4108      double * arrayTemp = array+offset;
4109      int * indexTemp = index+offset; 
4110      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
4111                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
4112      iwork[1] = dualColumn0(model,spare+offset2,
4113                                   spareIndex+offset2,
4114                                   arrayTemp,indexTemp,
4115                                   iwork[0],offset3,acceptablePivot,
4116                                   &dwork[0],&dwork[1],&iwork[2],
4117                                   &dwork[2]);
4118      int number = iwork[0];
4119      int numberLook = iwork[1];
4120#if 1
4121      numberRemaining += numberLook;
4122#else
4123      double * spareTemp = spare+offset2;
4124      const int * spareIndexTemp = spareIndex+offset2; 
4125      for (i=0;i<numberLook;i++) {
4126        double value = spareTemp[i];
4127        spareTemp[i]=0.0;
4128        spare[numberRemaining] = value;
4129        spareIndex[numberRemaining++] = spareIndexTemp[i];
4130      }
4131#endif
4132      if (dwork[2]>freePivot) {
4133        freePivot=dwork[2];
4134        posFree = iwork[2]+numberNonZero;
4135      } 
4136      upperTheta =  CoinMin(dwork[1],upperTheta);
4137      bestPossible = CoinMax(dwork[0],bestPossible);
4138      for (i=0;i<number;i++) {
4139        // double value = arrayTemp[i];
4140        //arrayTemp[i]=0.0;
4141        //array[numberNonZero]=value;
4142        index[numberNonZero++]=indexTemp[i]+offset3;
4143      }
4144#else
4145      int offset = offset_[iBlock];
4146      // allow for already saved
4147      int offset2 = offset + saveNumberRemaining;
4148      double * arrayTemp = array+offset;
4149      int * indexTemp = index+offset; 
4150      dualColumn0Struct * infoPtr = info_+iBlock;
4151      infoPtr->model=model;
4152      infoPtr->spare=spare+offset2;
4153      infoPtr->spareIndex=spareIndex+offset2;
4154      infoPtr->arrayTemp=arrayTemp;
4155      infoPtr->indexTemp=indexTemp;
4156      infoPtr->numberInPtr=&iwork[0];
4157      infoPtr->offset=offset;
4158      infoPtr->acceptablePivot=acceptablePivot;
4159      infoPtr->bestPossiblePtr=&dwork[0];
4160      infoPtr->upperThetaPtr=&dwork[1];
4161      infoPtr->posFreePtr=&iwork[2];
4162      infoPtr->freePivotPtr=&dwork[2];
4163      infoPtr->numberOutPtr=&iwork[1];
4164      infoPtr->pi=pi;
4165      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
4166      infoPtr->element=element;
4167      infoPtr->column=column_;
4168      infoPtr->numberInRowArray=numberInRowArray;
4169      infoPtr->numberLook=offset_[iBlock+1]-offset;
4170      if (iBlock>=2)
4171        pthread_join(threadId_[iBlock-2],NULL);
4172      pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
4173      //pthread_join(threadId_[iBlock],NULL);
4174#endif
4175    }
4176  }
4177  for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
4178#ifdef THREAD
4179    pthread_join(threadId_[iBlock],NULL);
4180#endif
4181  }
4182#ifdef THREAD
4183  for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
4184    //pthread_join(threadId_[iBlock],NULL);
4185    int offset = offset_[iBlock];
4186    double * dwork = work_+6*iBlock;
4187    int * iwork = (int *) (dwork+3);
4188    int number = iwork[0];
4189    if (dualColumn) {
4190      // allow for already saved
4191      int offset2 = offset + saveNumberRemaining;
4192      int numberLook = iwork[1];
4193      double * spareTemp = spare+offset2;
4194      const int * spareIndexTemp = spareIndex+offset2; 
4195      for (i=0;i<numberLook;i++) {
4196        double value = spareTemp[i];
4197        spareTemp[i]=0.0;
4198        spare[numberRemaining] = value;
4199        spareIndex[numberRemaining++] = spareIndexTemp[i];
4200      }
4201      if (dwork[2]>freePivot) {
4202        freePivot=dwork[2];
4203        posFree = iwork[2]+numberNonZero;
4204      } 
4205      upperTheta =  CoinMin(dwork[1],upperTheta);
4206      bestPossible = CoinMax(dwork[0],bestPossible);
4207    }
4208    double * arrayTemp = array+offset;
4209    const int * indexTemp = index+offset; 
4210    for (i=0;i<number;i++) {
4211      double value = arrayTemp[i];
4212      arrayTemp[i]=0.0; 
4213      array[numberNonZero]=value;
4214      index[numberNonZero++]=indexTemp[i]+offset;
4215    }
4216  }
4217#endif
4218  columnArray->setNumElements(numberNonZero);
4219  columnArray->setPackedMode(true);
4220  if (dualColumn) {
4221    model->spareDoubleArray_[0]=upperTheta;
4222    model->spareDoubleArray_[1]=bestPossible;
4223    // and theta and alpha and sequence
4224    if (posFree<0) {
4225      model->spareIntArray_[1]=-1;
4226    } else {
4227      const double * reducedCost=model->djRegion(0);
4228      double alpha;
4229      int numberColumns=model->numberColumns();
4230      if (posFree<numberColumns) {
4231        alpha = columnArray->denseVector()[posFree];
4232        posFree = columnArray->getIndices()[posFree];
4233      } else {
4234        alpha = rowArray->denseVector()[posFree-numberColumns];
4235        posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
4236      }
4237      model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
4238      model->spareDoubleArray_[3]=alpha;
4239      model->spareIntArray_[1]=posFree;
4240    }
4241    spareArray->setNumElements(numberRemaining);
4242    // signal done
4243    model->spareIntArray_[0]=-1;
4244  }
4245}
4246/* Default constructor. */
4247ClpPackedMatrix3::ClpPackedMatrix3()
4248  : numberBlocks_(0),
4249    numberColumns_(0),
4250    column_(NULL),
4251    start_(NULL),
4252    row_(NULL),
4253    element_(NULL),
4254    block_(NULL)
4255{
4256}
4257/* Constructor from copy. */
4258ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy)
4259  : numberBlocks_(0),
4260    numberColumns_(0),
4261    column_(NULL),
4262    start_(NULL),
4263    row_(NULL),
4264    element_(NULL),
4265    block_(NULL)
4266{
4267#define MINBLOCK 6
4268#define MAXBLOCK 100
4269#define MAXUNROLL 10
4270  numberColumns_ = columnCopy->getNumCols();
4271  int numberRows = columnCopy->getNumRows();
4272  int * counts = new int[numberRows+1];
4273  CoinZeroN(counts,numberRows+1);
4274  CoinBigIndex nels=0;
4275  CoinBigIndex nZeroEl=0;
4276  int iColumn;
4277  // get matrix data pointers
4278  const int * row = columnCopy->getIndices();
4279  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4280  const int * columnLength = columnCopy->getVectorLengths(); 
4281  const double * elementByColumn = columnCopy->getElements();
4282  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4283    CoinBigIndex start = columnStart[iColumn];
4284    int n = columnLength[iColumn];
4285    CoinBigIndex end = start + n;
4286    int kZero=0;
4287    for (CoinBigIndex j=start;j<end;j++) {
4288      if(!elementByColumn[j])
4289        kZero++;
4290    }
4291    n -= kZero;
4292    nZeroEl += kZero;
4293    nels += n;
4294    counts[n]++;
4295  }
4296  int nZeroColumns = counts[0];
4297  counts[0]=-1;
4298  numberColumns_ -= nZeroColumns;
4299  column_ = new int [2*numberColumns_];
4300  int * lookup = column_ + numberColumns_;
4301  row_ = new int[nels];
4302  element_ = new double[nels];
4303  int nOdd=0;
4304  CoinBigIndex nInOdd=0;
4305  int i;
4306  for (i=1;i<=numberRows;i++) {
4307    int n = counts[i];
4308    if (n) {
4309      if (n<MINBLOCK||i>MAXBLOCK) {
4310        nOdd += n;
4311        counts[i]=-1;
4312        nInOdd += n*i;
4313      } else {
4314        numberBlocks_++;
4315      }
4316    } else {
4317      counts[i]=-1;
4318    }
4319  }
4320  start_ = new CoinBigIndex [nOdd+1];
4321  // even if no blocks do a dummy one
4322  numberBlocks_ = CoinMax(numberBlocks_,1);
4323  block_ = (blockStruct *) new char [numberBlocks_*sizeof(blockStruct)];
4324  memset(block_,0,numberBlocks_*sizeof(blockStruct));
4325  // Fill in what we can
4326  int nTotal=nOdd;
4327  // in case no blocks
4328  block_->startIndices_=nTotal;
4329  nels=nInOdd;
4330  int nBlock=0;
4331  for (i=0;i<=CoinMin(MAXBLOCK,numberRows);i++) {
4332    if (counts[i]>0) {
4333      blockStruct * block = block_ + nBlock;
4334      int n=counts[i];
4335      counts[i]= nBlock; // backward pointer
4336      nBlock++;
4337      block->startIndices_ = nTotal;
4338      block->startElements_ = nels;
4339      block->numberElements_ = i;
4340      // up counts
4341      nTotal += n;
4342      nels += n*i;
4343    }
4344  }
4345  // fill
4346  start_[0]=0;
4347  nOdd=0;
4348  nInOdd=0;
4349  const double * columnScale = model->columnScale();
4350  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4351    CoinBigIndex start = columnStart[iColumn];
4352    int n = columnLength[iColumn];
4353    CoinBigIndex end = start + n;
4354    int kZero=0;
4355    for (CoinBigIndex j=start;j<end;j++) {
4356      if(!elementByColumn[j])
4357        kZero++;
4358    }
4359    n -= kZero;
4360    if (n) {
4361      int iBlock = counts[n];
4362      if (iBlock>=0) {
4363        blockStruct * block = block_ + iBlock;
4364        int k = block->numberInBlock_;
4365        block->numberInBlock_ ++;
4366        column_[block->startIndices_+k]=iColumn;
4367        lookup[iColumn]=k;
4368        CoinBigIndex put = block->startElements_+k*n;
4369        for (CoinBigIndex j=start;j<end;j++) {
4370          double value = elementByColumn[j];
4371          if(value) {
4372            if (columnScale)
4373              value *= columnScale[iColumn];
4374            element_[put] = value;
4375            row_[put++]=row[j];
4376          }
4377        }
4378      } else {
4379        // odd ones
4380        for (CoinBigIndex j=start;j<end;j++) {
4381          double value = elementByColumn[j];
4382          if(value) {
4383            if (columnScale)
4384              value *= columnScale[iColumn];
4385            element_[nInOdd] = value;
4386            row_[nInOdd++]=row[j];
4387          }
4388        }
4389        column_[nOdd]=iColumn;
4390        lookup[iColumn]=-1;
4391        nOdd++;
4392        start_[nOdd]=nInOdd;
4393      }
4394    }
4395  }
4396  delete [] counts;
4397}
4398/* Destructor */
4399ClpPackedMatrix3::~ClpPackedMatrix3()
4400{
4401  delete [] column_;
4402  delete [] start_;
4403  delete [] row_;
4404  delete [] element_;
4405  delete [] block_;
4406}
4407/* The copy constructor. */
4408ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
4409  : numberBlocks_(rhs.numberBlocks_),
4410    numberColumns_(rhs.numberColumns_),
4411    column_(NULL),
4412    start_(NULL),
4413    row_(NULL),
4414    element_(NULL),
4415    block_(NULL)
4416{
4417  if (rhs.numberBlocks_) {
4418    block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
4419    column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
4420    int numberOdd = block_->startIndices_;
4421    start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
4422    blockStruct * lastBlock = block_ + (numberBlocks_-1);
4423    CoinBigIndex numberElements = lastBlock->startElements_+
4424      lastBlock->numberInBlock_*lastBlock->numberElements_;
4425    row_ = CoinCopyOfArray(rhs.row_,numberElements);
4426    element_ = CoinCopyOfArray(rhs.element_,numberElements);
4427  }
4428}
4429ClpPackedMatrix3& 
4430ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
4431{
4432  if (this != &rhs) {
4433    delete [] column_;
4434    delete [] start_;
4435    delete [] row_;
4436    delete [] element_;
4437    delete [] block_;
4438    numberBlocks_ = rhs.numberBlocks_;
4439    numberColumns_ = rhs.numberColumns_;
4440    if (rhs.numberBlocks_) {
4441      block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
4442      column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
4443      int numberOdd = block_->startIndices_;
4444      start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
4445      blockStruct * lastBlock = block_ + (numberBlocks_-1);
4446      CoinBigIndex numberElements = lastBlock->startElements_+
4447        lastBlock->numberInBlock_*lastBlock->numberElements_;
4448      row_ = CoinCopyOfArray(rhs.row_,numberElements);
4449      element_ = CoinCopyOfArray(rhs.element_,numberElements);
4450    } else {
4451      column_ = NULL;
4452      start_ = NULL;
4453      row_ = NULL;
4454      element_ = NULL;
4455      block_ = NULL;
4456    }
4457  }
4458  return *this;
4459}
4460/* Sort blocks */
4461void 
4462ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
4463{
4464  int * lookup = column_ + numberColumns_;
4465  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4466    blockStruct * block = block_ + iBlock;
4467    int numberInBlock = block->numberInBlock_; 
4468    int nel = block->numberElements_;
4469    int * row = row_ + block->startElements_;
4470    double * element = element_ + block->startElements_;
4471    int * column = column_ + block->startIndices_;
4472    int lastPrice=0;
4473    int firstNotPrice=numberInBlock-1;
4474    while (lastPrice<=firstNotPrice) {
4475      // find first basic or fixed
4476      int iColumn=numberInBlock;
4477      for (;lastPrice<=firstNotPrice;lastPrice++) {
4478        iColumn = column[lastPrice];
4479        if (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4480            model->getColumnStatus(iColumn)==ClpSimplex::isFixed)
4481          break;
4482      }
4483      // find last non basic or fixed
4484      int jColumn=-1;
4485      for (;firstNotPrice>lastPrice;firstNotPrice--) {
4486        jColumn = column[firstNotPrice];
4487        if (model->getColumnStatus(jColumn)!=ClpSimplex::basic&&
4488            model->getColumnStatus(jColumn)!=ClpSimplex::isFixed)
4489          break;
4490      }
4491      if (firstNotPrice>lastPrice) {
4492        assert (column[lastPrice]==iColumn);
4493        assert (column[firstNotPrice]==jColumn);
4494        // need to swap
4495        column[firstNotPrice]=iColumn;
4496        lookup[iColumn]=firstNotPrice;
4497        column[lastPrice]=jColumn;
4498        lookup[jColumn]=lastPrice;
4499        double * elementA = element + lastPrice*nel;
4500        int * rowA = row + lastPrice*nel;
4501        double * elementB = element + firstNotPrice*nel;
4502        int * rowB = row + firstNotPrice*nel;
4503        for (int i=0;i<nel;i++) {
4504          int temp = rowA[i];
4505          double tempE = elementA[i];
4506          rowA[i]=rowB[i];
4507          elementA[i]=elementB[i];
4508          rowB[i]=temp;
4509          elementB[i]=tempE;
4510        }
4511        firstNotPrice--;
4512        lastPrice++;
4513      } else if (lastPrice==firstNotPrice) {
4514        // make sure correct side
4515        iColumn = column[lastPrice];
4516        if (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4517            model->getColumnStatus(iColumn)!=ClpSimplex::isFixed)
4518          lastPrice++;
4519        break;
4520      }
4521    }
4522    block->numberPrice_=lastPrice;
4523#ifndef NDEBUG
4524    // check
4525    int i;
4526    for (i=0;i<lastPrice;i++) {
4527      int iColumn = column[i];
4528      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4529              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
4530      assert (lookup[iColumn]==i);
4531    }
4532    for (;i<numberInBlock;i++) {
4533      int iColumn = column[i];
4534      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4535              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4536      assert (lookup[iColumn]==i);
4537    }
4538#endif
4539  }
4540}
4541// Swap one variable
4542void 
4543ClpPackedMatrix3::swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
4544                          int iColumn)
4545{
4546  int * lookup = column_ + numberColumns_;
4547  // position in block
4548  int kA = lookup[iColumn];
4549  if (kA<0)
4550    return; // odd one
4551  // get matrix data pointers
4552  const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
4553  //const int * row = columnCopy->getIndices();
4554  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4555  const int * columnLength = columnCopy->getVectorLengths(); 
4556  const double * elementByColumn = columnCopy->getElements();
4557  CoinBigIndex start = columnStart[iColumn];
4558  int n = columnLength[iColumn];
4559  if (matrix->zeros()) {
4560    CoinBigIndex end = start + n;
4561    for (CoinBigIndex j=start;j<end;j++) {
4562      if(!elementByColumn[j])
4563        n--;
4564    }
4565  }
4566  // find block - could do binary search
4567  int iBlock=CoinMin(n,numberBlocks_)-1;
4568  while (block_[iBlock].numberElements_!=n)
4569    iBlock--;
4570  blockStruct * block = block_ + iBlock;
4571  int nel = block->numberElements_;
4572  int * row = row_ + block->startElements_;
4573  double * element = element_ + block->startElements_;
4574  int * column = column_ + block->startIndices_;
4575  assert (column[kA]==iColumn);
4576  bool moveUp = (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4577                 model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4578  int lastPrice=block->numberPrice_;
4579  int kB;
4580  if (moveUp) {
4581    kB=lastPrice-1;
4582    block->numberPrice_--;
4583  } else {
4584    kB=lastPrice;
4585    block->numberPrice_++;
4586  }
4587  int jColumn = column[kB];
4588  column[kA]=jColumn;
4589  lookup[jColumn]=kA;
4590  column[kB]=iColumn;
4591  lookup[iColumn]=kB;
4592  double * elementA = element + kB*nel;
4593  int * rowA = row + kB*nel;
4594  double * elementB = element + kA*nel;
4595  int * rowB = row + kA*nel;
4596  for (int i=0;i<nel;i++) {
4597    int temp = rowA[i];
4598    double tempE = elementA[i];
4599    rowA[i]=rowB[i];
4600    elementA[i]=elementB[i];
4601    rowB[i]=temp;
4602    elementB[i]=tempE;
4603  }
4604#if 1
4605#ifndef NDEBUG
4606  // check
4607  int i;
4608  for (i=0;i<block->numberPrice_;i++) {
4609    int iColumn = column[i];
4610    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
4611      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4612              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
4613    assert (lookup[iColumn]==i);
4614  }
4615  int numberInBlock = block->numberInBlock_; 
4616  for (;i<numberInBlock;i++) {
4617    int iColumn = column[i];
4618    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
4619      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4620              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4621    assert (lookup[iColumn]==i);
4622  }
4623#endif
4624#endif
4625}
4626/* Return <code>x * -1 * A in <code>z</code>.
4627   Note - x packed and z will be packed mode
4628   Squashes small elements and knows about ClpSimplex */
4629void 
4630ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
4631                                 const double * pi,
4632                                 CoinIndexedVector * output) const
4633{
4634  int numberNonZero=0;
4635  int * index = output->getIndices();
4636  double * array = output->denseVector();
4637  double zeroTolerance = model->factorization()->zeroTolerance();
4638  double value = 0.0;
4639  CoinBigIndex j;
4640  int numberOdd = block_->startIndices_;
4641  if (numberOdd) {
4642    // A) as probably long may be worth unrolling
4643    CoinBigIndex end = start_[1];
4644    for (j=start_[0]; j<end;j++) {
4645      int iRow = row_[j];
4646      value += pi[iRow]*element_[j];
4647    }
4648    int iColumn;
4649    // int jColumn=column_[0];
4650   
4651    for (iColumn=0;iColumn<numberOdd-1;iColumn++) {
4652      CoinBigIndex start = end;
4653      end = start_[iColumn+2];
4654      if (fabs(value)>zeroTolerance) {
4655        array[numberNonZero]=value;
4656        index[numberNonZero++]=column_[iColumn];
4657        //index[numberNonZero++]=jColumn;
4658      }
4659      // jColumn = column_[iColumn+1];
4660      value = 0.0;
4661      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
4662        for (j=start; j<end;j++) {
4663          int iRow = row_[j];
4664          value += pi[iRow]*element_[j];
4665        }
4666        //}
4667    }
4668    if (fabs(value)>zeroTolerance) {
4669      array[numberNonZero]=value;
4670      index[numberNonZero++]=column_[iColumn];
4671      //index[numberNonZero++]=jColumn;
4672    }
4673  }
4674  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4675    // B) Can sort so just do nonbasic (and nonfixed)
4676    // C) Can do two at a time (if so put odd one into start_)
4677    // D) can use switch
4678    blockStruct * block = block_ + iBlock;
4679    //int numberPrice = block->numberInBlock_;
4680    int numberPrice = block->numberPrice_; 
4681    int nel = block->numberElements_;
4682    int * row = row_ + block->startElements_;
4683    double * element = element_ + block->startElements_;
4684    int * column = column_ + block->startIndices_;
4685#if 0
4686    // two at a time
4687    if ((numberPrice&1)!=0) {
4688      double value = 0.0;
4689      int nel2=nel;
4690      for (; nel2;nel2--) {
4691        int iRow = *row++;
4692        value += pi[iRow]*(*element++);
4693      }
4694      if (fabs(value)>zeroTolerance) {
4695        array[numberNonZero]=value;
4696        index[numberNonZero++]=*column;
4697      }
4698      column++;
4699    }
4700    numberPrice = numberPrice>>1;
4701    switch ((nel%2)) {
4702      // 2 k +0
4703    case 0:
4704      for (;numberPrice;numberPrice--) {
4705        double value0 = 0.0;
4706        double value1 = 0.0;
4707        int nel2=nel;
4708        for (; nel2;nel2--) {
4709          int iRow0 = *row;
4710          int iRow1 = *(row+nel);
4711          row++;
4712          double element0 = *element;
4713          double element1 = *(element+nel);
4714          element++;
4715          value0 += pi[iRow0]*element0;
4716          value1 += pi[iRow1]*element1;
4717        }
4718        row += nel;
4719        element += nel;
4720        if (fabs(value0)>zeroTolerance) {
4721          array[numberNonZero]=value0;
4722          index[numberNonZero++]=*column;
4723        }
4724        column++;
4725        if (fabs(value1)>zeroTolerance) {
4726          array[numberNonZero]=value1;
4727          index[numberNonZero++]=*column;
4728        }
4729        column++;
4730      } 
4731      break;
4732      // 2 k +1
4733    case 1:
4734      for (;numberPrice;numberPrice--) {
4735        double value0;
4736        double value1;
4737        int nel2=nel-1;
4738        {
4739          int iRow0 = row[0];
4740          int iRow1 = row[nel];
4741          double pi0 = pi[iRow0];
4742          double pi1 = pi[iRow1];
4743          value0 = pi0*element[0];
4744          value1 = pi1*element[nel];
4745          row++;
4746          element++;
4747        }
4748        for (; nel2;nel2--) {
4749          int iRow0 = *row;
4750          int iRow1 = *(row+nel);
4751          row++;
4752          double element0 = *element;
4753          double element1 = *(element+nel);
4754          element++;
4755          value0 += pi[iRow0]*element0;
4756          value1 += pi[iRow1]*element1;
4757        }
4758        row += nel;
4759        element += nel;
4760        if (fabs(value0)>zeroTolerance) {
4761          array[numberNonZero]=value0;
4762          index[numberNonZero++]=*column;
4763        }
4764        column++;
4765        if (fabs(value1)>zeroTolerance) {
4766          array[numberNonZero]=value1;
4767          index[numberNonZero++]=*column;
4768        }
4769        column++;
4770      }
4771      break;
4772    }
4773#else
4774    for (;numberPrice;numberPrice--) {
4775      double value = 0.0;
4776      int nel2=nel;
4777      for (; nel2;nel2--) {
4778        int iRow = *row++;
4779        value += pi[iRow]*(*element++);
4780      }
4781      if (fabs(value)>zeroTolerance) {
4782        array[numberNonZero]=value;
4783        index[numberNonZero++]=*column;
4784      }
4785      column++;
4786    }
4787#endif
4788  }
4789  output->setNumElements(numberNonZero);
4790}
4791// Updates two arrays for steepest
4792void 
4793ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
4794                                  const double * pi, CoinIndexedVector * output,
4795                                  const double * piWeight,
4796                                  double referenceIn, double devex,
4797                                  // Array for exact devex to say what is in reference framework
4798                                  unsigned int * reference,
4799                                  double * weights, double scaleFactor)
4800{
4801  int numberNonZero=0;
4802  int * index = output->getIndices();
4803  double * array = output->denseVector();
4804  double zeroTolerance = model->factorization()->zeroTolerance();
4805  double value = 0.0;
4806  bool killDjs = (scaleFactor==0.0);
4807  if (!scaleFactor)
4808    scaleFactor=1.0;
4809  int numberOdd = block_->startIndices_;
4810  int iColumn;
4811  CoinBigIndex end = start_[0];
4812  for (iColumn=0;iColumn<numberOdd;iColumn++) {
4813    CoinBigIndex start = end;
4814    CoinBigIndex j;
4815    int jColumn = column_[iColumn];
4816    end = start_[iColumn+1];
4817    value = 0.0;
4818    if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
4819      for (j=start; j<end;j++) {
4820        int iRow = row_[j];
4821        value -= pi[iRow]*element_[j];
4822      }
4823      if (fabs(value)>zeroTolerance) {
4824        // and do other array
4825        double modification = 0.0;
4826        for (j=start; j<end;j++) {
4827          int iRow = row_[j];
4828          modification += piWeight[iRow]*element_[j];
4829        }
4830        double thisWeight = weights[jColumn];
4831        double pivot = value*scaleFactor;
4832        double pivotSquared = pivot * pivot;
4833        thisWeight += pivotSquared * devex + pivot * modification;
4834        if (thisWeight<DEVEX_TRY_NORM) {
4835          if (referenceIn<0.0) {
4836            // steepest
4837            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
4838          } else {
4839            // exact
4840            thisWeight = referenceIn*pivotSquared;
4841            if (reference(jColumn))
4842              thisWeight += 1.0;
4843            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
4844          }
4845        }
4846        weights[jColumn] = thisWeight;
4847        if (!killDjs) {
4848          array[numberNonZero]=value;
4849          index[numberNonZero++]=jColumn;
4850        }
4851      }
4852    }
4853  }
4854  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4855    // B) Can sort so just do nonbasic (and nonfixed)
4856    // C) Can do two at a time (if so put odd one into start_)
4857    // D) can use switch
4858    blockStruct * block = block_ + iBlock;
4859    //int numberPrice = block->numberInBlock_;
4860    int numberPrice = block->numberPrice_; 
4861    int nel = block->numberElements_;
4862    int * row = row_ + block->startElements_;
4863    double * element = element_ + block->startElements_;
4864    int * column = column_ + block->startIndices_;
4865    for (;numberPrice;numberPrice--) {
4866      double value = 0.0;
4867      int nel2=nel;
4868      for (; nel2;nel2--) {
4869        int iRow = *row++;
4870        value -= pi[iRow]*(*element++);
4871      }
4872      if (fabs(value)>zeroTolerance) {
4873        int jColumn = *column;
4874        // back to beginning
4875        row -= nel;
4876        element -= nel;
4877        // and do other array
4878        double modification = 0.0;
4879        nel2=nel;
4880        for (; nel2;nel2--) {
4881          int iRow = *row++;
4882          modification += piWeight[iRow]*(*element++);
4883        }
4884        double thisWeight = weights[jColumn];
4885        double pivot = value*scaleFactor;
4886        double pivotSquared = pivot * pivot;
4887        thisWeight += pivotSquared * devex + pivot * modification;
4888        if (thisWeight<DEVEX_TRY_NORM) {
4889          if (referenceIn<0.0) {
4890            // steepest
4891            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
4892          } else {
4893            // exact
4894            thisWeight = referenceIn*pivotSquared;
4895            if (reference(jColumn))
4896              thisWeight += 1.0;
4897            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
4898          }
4899        }
4900        weights[jColumn] = thisWeight;
4901        if (!killDjs) {
4902          array[numberNonZero]=value;
4903          index[numberNonZero++]=jColumn;
4904        }
4905      }
4906      column++;
4907    }
4908  }
4909  output->setNumElements(numberNonZero);
4910  output->setPackedMode(true);
4911}
4912#ifdef CLP_ALL_ONE_FILE
4913#undef reference
4914#endif
Note: See TracBrowser for help on using the repository browser.