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

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

adding possibility of long doubles

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