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

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

out compiler warnings and stability improvements

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