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

Last change on this file since 1310 was 1310, checked in by forrest, 11 years ago

second try at fixing duplicate problem

  • 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  if (numberSmall||numberDuplicate)
2955    flags_ |= 2; // will have gaps
2956  return true;
2957}
2958int
2959ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector, 
2960                                              int * COIN_RESTRICT index, 
2961                                              double * COIN_RESTRICT output,
2962                                              int * COIN_RESTRICT lookup,
2963                                              char * COIN_RESTRICT marked,
2964                                              const double tolerance, 
2965                                              const double scalar) const
2966{
2967  const double * COIN_RESTRICT pi = piVector->denseVector();
2968  int numberNonZero=0;
2969  int numberInRowArray = piVector->getNumElements();
2970  const int * COIN_RESTRICT column = matrix_->getIndices();
2971  const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
2972  const double * COIN_RESTRICT element = matrix_->getElements();
2973  const int * COIN_RESTRICT whichRow = piVector->getIndices();
2974#ifndef NDEBUG
2975  int maxColumn=0;
2976#endif
2977  // ** Row copy is already scaled
2978  int iRow;
2979  int i;
2980  for (i=0;i<numberInRowArray;i++) {
2981    iRow = whichRow[i]; 
2982    double value = pi[i]*scalar;
2983    CoinBigIndex j;
2984    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
2985      int iColumn = column[j];
2986#ifndef NDEBUG
2987      maxColumn=CoinMax(maxColumn,iColumn);
2988#endif
2989      double elValue = element[j];
2990      elValue *= value;
2991      if (marked[iColumn]) {
2992        int k = lookup[iColumn];
2993        output[k] += elValue;
2994      } else {
2995        output[numberNonZero] = elValue;
2996        marked[iColumn]=1;
2997        lookup[iColumn]=numberNonZero;
2998        index[numberNonZero++]=iColumn;
2999      }
3000    }
3001  }
3002#ifndef NDEBUG
3003  int saveN = numberNonZero;
3004#endif
3005  // get rid of tiny values and zero out lookup
3006  for (i=0;i<numberNonZero;i++) {
3007    int iColumn = index[i];
3008    marked[iColumn]=0;
3009    double value = output[i];
3010    if (fabs(value)<=tolerance) {
3011      while (fabs(value)<=tolerance) {
3012        numberNonZero--;
3013        value = output[numberNonZero];
3014        iColumn = index[numberNonZero];
3015        marked[iColumn]=0;
3016        if (i<numberNonZero) {
3017          output[numberNonZero]=0.0;
3018          output[i] = value;
3019          index[i] = iColumn;
3020        } else {
3021          output[i]=0.0;
3022          value =1.0; // to force end of while
3023        }
3024      }
3025    }
3026  }
3027#ifndef NDEBUG
3028  for (i=numberNonZero;i<saveN;i++)
3029    assert(!output[i]);
3030  for (i=0;i<=maxColumn;i++)
3031    assert (!marked[i]);
3032#endif
3033  return numberNonZero;
3034}
3035/* Given positive integer weights for each row fills in sum of weights
3036   for each column (and slack).
3037   Returns weights vector
3038*/
3039CoinBigIndex * 
3040ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
3041{
3042  int numberRows = model->numberRows();
3043  int numberColumns = matrix_->getNumCols();
3044  int number = numberRows+numberColumns;
3045  CoinBigIndex * weights = new CoinBigIndex[number];
3046  // get matrix data pointers
3047  const int * row = matrix_->getIndices();
3048  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3049  const int * columnLength = matrix_->getVectorLengths(); 
3050  int i;
3051  for (i=0;i<numberColumns;i++) {
3052    CoinBigIndex j;
3053    CoinBigIndex count=0;
3054    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
3055      int iRow=row[j];
3056      count += inputWeights[iRow];
3057    }
3058    weights[i]=count;
3059  }
3060  for (i=0;i<numberRows;i++) {
3061    weights[i+numberColumns]=inputWeights[i];
3062  }
3063  return weights;
3064}
3065/* Returns largest and smallest elements of both signs.
3066   Largest refers to largest absolute value.
3067*/
3068void 
3069ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
3070                       double & smallestPositive, double & largestPositive)
3071{
3072  smallestNegative=-COIN_DBL_MAX;
3073  largestNegative=0.0;
3074  smallestPositive=COIN_DBL_MAX;
3075  largestPositive=0.0;
3076  // get matrix data pointers
3077  const double * elementByColumn = matrix_->getElements();
3078  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3079  const int * columnLength = matrix_->getVectorLengths(); 
3080  int numberColumns = matrix_->getNumCols();
3081  int i;
3082  for (i=0;i<numberColumns;i++) {
3083    CoinBigIndex j;
3084    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
3085      double value = elementByColumn[j];
3086      if (value>0.0) {
3087        smallestPositive = CoinMin(smallestPositive,value);
3088        largestPositive = CoinMax(largestPositive,value);
3089      } else if (value<0.0) {
3090        smallestNegative = CoinMax(smallestNegative,value);
3091        largestNegative = CoinMin(largestNegative,value);
3092      }
3093    }
3094  }
3095}
3096// Says whether it can do partial pricing
3097bool 
3098ClpPackedMatrix::canDoPartialPricing() const
3099{
3100  return true;
3101}
3102// Partial pricing
3103void 
3104ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
3105                              int & bestSequence, int & numberWanted)
3106{
3107  numberWanted=currentWanted_;
3108  int start = (int) (startFraction*numberActiveColumns_);
3109  int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
3110  const double * element =matrix_->getElements();
3111  const int * row = matrix_->getIndices();
3112  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
3113  const int * length = matrix_->getVectorLengths();
3114  const double * rowScale = model->rowScale();
3115  const double * columnScale = model->columnScale();
3116  int iSequence;
3117  CoinBigIndex j;
3118  double tolerance=model->currentDualTolerance();
3119  double * reducedCost = model->djRegion();
3120  const double * duals = model->dualRowSolution();
3121  const double * cost = model->costRegion();
3122  double bestDj;
3123  if (bestSequence>=0)
3124    bestDj = fabs(model->clpMatrix()->reducedCost(model,bestSequence));
3125  else
3126    bestDj=tolerance;
3127  int sequenceOut = model->sequenceOut();
3128  int saveSequence = bestSequence;
3129  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 
3130  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
3131  if (rowScale) {
3132    // scaled
3133    for (iSequence=start;iSequence<end;iSequence++) {
3134      if (iSequence!=sequenceOut) {
3135        double value;
3136        ClpSimplex::Status status = model->getStatus(iSequence);
3137
3138        switch(status) {
3139         
3140        case ClpSimplex::basic:
3141        case ClpSimplex::isFixed:
3142          break;
3143        case ClpSimplex::isFree:
3144        case ClpSimplex::superBasic:
3145          value=0.0;
3146          // scaled
3147          for (j=startColumn[iSequence];
3148               j<startColumn[iSequence]+length[iSequence];j++) {
3149            int jRow=row[j];
3150            value -= duals[jRow]*element[j]*rowScale[jRow];
3151          }
3152          value = fabs(cost[iSequence] +value*columnScale[iSequence]);
3153          if (value>FREE_ACCEPT*tolerance) {
3154            numberWanted--;
3155            // we are going to bias towards free (but only if reasonable)
3156            value *= FREE_BIAS;
3157            if (value>bestDj) {
3158              // check flagged variable and correct dj
3159              if (!model->flagged(iSequence)) {
3160                bestDj=value;
3161                bestSequence = iSequence;
3162              } else {
3163                // just to make sure we don't exit before got something
3164                numberWanted++;
3165              }
3166            }
3167          }
3168          break;
3169        case ClpSimplex::atUpperBound:
3170          value=0.0;
3171          // scaled
3172          for (j=startColumn[iSequence];
3173               j<startColumn[iSequence]+length[iSequence];j++) {
3174            int jRow=row[j];
3175            value -= duals[jRow]*element[j]*rowScale[jRow];
3176          }
3177          value = cost[iSequence] +value*columnScale[iSequence];
3178          if (value>tolerance) {
3179            numberWanted--;
3180            if (value>bestDj) {
3181              // check flagged variable and correct dj
3182              if (!model->flagged(iSequence)) {
3183                bestDj=value;
3184                bestSequence = iSequence;
3185              } else {
3186                // just to make sure we don't exit before got something
3187                numberWanted++;
3188              }
3189            }
3190          }
3191          break;
3192        case ClpSimplex::atLowerBound:
3193          value=0.0;
3194          // scaled
3195          for (j=startColumn[iSequence];
3196               j<startColumn[iSequence]+length[iSequence];j++) {
3197            int jRow=row[j];
3198            value -= duals[jRow]*element[j]*rowScale[jRow];
3199          }
3200          value = -(cost[iSequence] +value*columnScale[iSequence]);
3201          if (value>tolerance) {
3202            numberWanted--;
3203            if (value>bestDj) {
3204              // check flagged variable and correct dj
3205              if (!model->flagged(iSequence)) {
3206                bestDj=value;
3207                bestSequence = iSequence;
3208              } else {
3209                // just to make sure we don't exit before got something
3210                numberWanted++;
3211              }
3212            }
3213          }
3214          break;
3215        }
3216      }
3217      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
3218        // give up
3219        break;
3220      }
3221      if (!numberWanted)
3222        break;
3223    }
3224    if (bestSequence!=saveSequence) {
3225      // recompute dj
3226      double value=0.0;
3227      // scaled
3228      for (j=startColumn[bestSequence];
3229           j<startColumn[bestSequence]+length[bestSequence];j++) {
3230        int jRow=row[j];
3231        value -= duals[jRow]*element[j]*rowScale[jRow];
3232      }
3233      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
3234      savedBestSequence_ = bestSequence;
3235      savedBestDj_ = reducedCost[savedBestSequence_];
3236    }
3237  } else {
3238    // not scaled
3239    for (iSequence=start;iSequence<end;iSequence++) {
3240      if (iSequence!=sequenceOut) {
3241        double value;
3242        ClpSimplex::Status status = model->getStatus(iSequence);
3243
3244        switch(status) {
3245         
3246        case ClpSimplex::basic:
3247        case ClpSimplex::isFixed:
3248          break;
3249        case ClpSimplex::isFree:
3250        case ClpSimplex::superBasic:
3251          value=cost[iSequence];
3252          for (j=startColumn[iSequence];
3253               j<startColumn[iSequence]+length[iSequence];j++) {
3254            int jRow=row[j];
3255            value -= duals[jRow]*element[j];
3256          }
3257          value = fabs(value);
3258          if (value>FREE_ACCEPT*tolerance) {
3259            numberWanted--;
3260            // we are going to bias towards free (but only if reasonable)
3261            value *= FREE_BIAS;
3262            if (value>bestDj) {
3263              // check flagged variable and correct dj
3264              if (!model->flagged(iSequence)) {
3265                bestDj=value;
3266                bestSequence = iSequence;
3267              } else {
3268                // just to make sure we don't exit before got something
3269                numberWanted++;
3270              }
3271            }
3272          }
3273          break;
3274        case ClpSimplex::atUpperBound:
3275          value=cost[iSequence];
3276          // scaled
3277          for (j=startColumn[iSequence];
3278               j<startColumn[iSequence]+length[iSequence];j++) {
3279            int jRow=row[j];
3280            value -= duals[jRow]*element[j];
3281          }
3282          if (value>tolerance) {
3283            numberWanted--;
3284            if (value>bestDj) {
3285              // check flagged variable and correct dj
3286              if (!model->flagged(iSequence)) {
3287                bestDj=value;
3288                bestSequence = iSequence;
3289              } else {
3290                // just to make sure we don't exit before got something
3291                numberWanted++;
3292              }
3293            }
3294          }
3295          break;
3296        case ClpSimplex::atLowerBound:
3297          value=cost[iSequence];
3298          for (j=startColumn[iSequence];
3299               j<startColumn[iSequence]+length[iSequence];j++) {
3300            int jRow=row[j];
3301            value -= duals[jRow]*element[j];
3302          }
3303          value = -value;
3304          if (value>tolerance) {
3305            numberWanted--;
3306            if (value>bestDj) {
3307              // check flagged variable and correct dj
3308              if (!model->flagged(iSequence)) {
3309                bestDj=value;
3310                bestSequence = iSequence;
3311              } else {
3312                // just to make sure we don't exit before got something
3313                numberWanted++;
3314              }
3315            }
3316          }
3317          break;
3318        }
3319      }
3320      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
3321        // give up
3322        break;
3323      }
3324      if (!numberWanted)
3325        break;
3326    }
3327    if (bestSequence!=saveSequence) {
3328      // recompute dj
3329      double value=cost[bestSequence];
3330      for (j=startColumn[bestSequence];
3331           j<startColumn[bestSequence]+length[bestSequence];j++) {
3332        int jRow=row[j];
3333        value -= duals[jRow]*element[j];
3334      }
3335      reducedCost[bestSequence] = value;
3336      savedBestSequence_ = bestSequence;
3337      savedBestDj_ = reducedCost[savedBestSequence_];
3338    }
3339  }
3340  currentWanted_=numberWanted;
3341}
3342// Sets up an effective RHS
3343void 
3344ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
3345{
3346  delete [] rhsOffset_;
3347  int numberRows = model->numberRows();
3348  rhsOffset_ = new double[numberRows];
3349  rhsOffset(model,true);
3350}
3351// Gets rid of special copies
3352void 
3353ClpPackedMatrix::clearCopies()
3354{
3355  delete rowCopy_;
3356  delete columnCopy_;
3357  rowCopy_=NULL;
3358  columnCopy_=NULL;
3359  flags_ &= ~(4+8);
3360 }
3361// makes sure active columns correct
3362int 
3363ClpPackedMatrix::refresh(ClpSimplex * model)
3364{
3365  numberActiveColumns_ = matrix_->getNumCols();
3366#if 0
3367  ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
3368  ClpPackedMatrix* rowCopy =
3369    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3370  // Make sure it is really a ClpPackedMatrix
3371  assert (rowCopy!=NULL);
3372
3373  const int * column = rowCopy->matrix_->getIndices();
3374  const CoinBigIndex * rowStart = rowCopy->matrix_->getVectorStarts();
3375  const int * rowLength = rowCopy->matrix_->getVectorLengths();
3376  const double * element = rowCopy->matrix_->getElements();
3377  int numberRows = rowCopy->matrix_->getNumRows();
3378  for (int i=0;i<numberRows;i++) {
3379    if (!rowLength[i])
3380      printf("zero row %d\n",i);
3381  }
3382  delete rowCopy;
3383#endif
3384  return 0;
3385}
3386
3387/* Scales rowCopy if column copy scaled
3388   Only called if scales already exist */
3389void 
3390ClpPackedMatrix::scaleRowCopy(ClpModel * model) const 
3391{
3392  if (model->rowCopy()) {
3393    // need to replace row by row
3394    int numberRows = model->numberRows();
3395#ifndef NDEBUG
3396    int numberColumns = matrix_->getNumCols();
3397#endif
3398    ClpMatrixBase * rowCopyBase=model->rowCopy();
3399#ifndef NO_RTTI
3400    ClpPackedMatrix* rowCopy =
3401      dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3402    // Make sure it is really a ClpPackedMatrix
3403    assert (rowCopy!=NULL);
3404#else
3405    ClpPackedMatrix* rowCopy =
3406      static_cast< ClpPackedMatrix*>(rowCopyBase);
3407#endif
3408
3409    const int * column = rowCopy->getIndices();
3410    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3411    double * element = rowCopy->getMutableElements();
3412    const double * rowScale = model->rowScale();
3413    const double * columnScale = model->columnScale();
3414    // scale row copy
3415    for (int iRow=0;iRow<numberRows;iRow++) {
3416      CoinBigIndex j;
3417      double scale = rowScale[iRow];
3418      double * elementsInThisRow = element + rowStart[iRow];
3419      const int * columnsInThisRow = column + rowStart[iRow];
3420      int number = rowStart[iRow+1]-rowStart[iRow];
3421      assert (number<=numberColumns);
3422      for (j=0;j<number;j++) {
3423        int iColumn = columnsInThisRow[j];
3424        elementsInThisRow[j] *= scale*columnScale[iColumn];
3425      }
3426    }
3427  }
3428}
3429/* Realy really scales column copy
3430   Only called if scales already exist.
3431   Up to user ro delete */
3432ClpMatrixBase * 
3433ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const 
3434{
3435  // need to replace column by column
3436#ifndef NDEBUG
3437  int numberRows = model->numberRows();
3438#endif
3439  int numberColumns = matrix_->getNumCols();
3440  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
3441  const int * row = copy->getIndices();
3442  const CoinBigIndex * columnStart = copy->getVectorStarts();
3443  const int * length = copy->getVectorLengths();
3444  double * element = copy->getMutableElements();
3445  const double * rowScale = model->rowScale();
3446  const double * columnScale = model->columnScale();
3447  // scale column copy
3448  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3449    CoinBigIndex j;
3450    double scale = columnScale[iColumn];
3451    double * elementsInThisColumn = element + columnStart[iColumn];
3452    const int * rowsInThisColumn = row + columnStart[iColumn];
3453    int number = length[iColumn];
3454    assert (number<=numberRows);
3455    for (j=0;j<number;j++) {
3456      int iRow = rowsInThisColumn[j];
3457      elementsInThisColumn[j] *= scale*rowScale[iRow];
3458    }
3459  }
3460  return copy;
3461}
3462// Really scale matrix
3463void 
3464ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
3465{
3466  clearCopies();
3467  int numberColumns = matrix_->getNumCols();
3468  const int * row = matrix_->getIndices();
3469  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3470  const int * length = matrix_->getVectorLengths();
3471  double * element = matrix_->getMutableElements();
3472  // scale
3473  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3474    CoinBigIndex j;
3475    double scale = columnScale[iColumn];
3476    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
3477      int iRow = row[j];
3478      element[j] *= scale*rowScale[iRow];
3479    }
3480  }
3481}
3482/* Delete the columns whose indices are listed in <code>indDel</code>. */
3483void 
3484ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
3485{ 
3486  if (matrix_->getNumCols())
3487    matrix_->deleteCols(numDel,indDel);
3488  clearCopies();
3489  numberActiveColumns_ = matrix_->getNumCols();
3490  // may now have gaps
3491  flags_ |= 2;
3492  matrix_->setExtraGap(0.0);
3493}
3494/* Delete the rows whose indices are listed in <code>indDel</code>. */
3495void 
3496ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
3497{ 
3498  if (matrix_->getNumRows())
3499    matrix_->deleteRows(numDel,indDel);
3500  clearCopies();
3501  numberActiveColumns_ = matrix_->getNumCols();
3502  // may now have gaps
3503  flags_ |= 2;
3504  matrix_->setExtraGap(0.0);
3505}
3506#ifndef CLP_NO_VECTOR
3507// Append Columns
3508void 
3509ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
3510{ 
3511  matrix_->appendCols(number,columns);
3512  numberActiveColumns_ = matrix_->getNumCols();
3513  clearCopies();
3514}
3515// Append Rows
3516void 
3517ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
3518{ 
3519  matrix_->appendRows(number,rows);
3520  numberActiveColumns_ = matrix_->getNumCols();
3521  // may now have gaps
3522  flags_ |= 2;
3523  clearCopies();
3524}
3525#endif
3526/* Set the dimensions of the matrix. In effect, append new empty
3527   columns/rows to the matrix. A negative number for either dimension
3528   means that that dimension doesn't change. Otherwise the new dimensions
3529   MUST be at least as large as the current ones otherwise an exception
3530   is thrown. */
3531void 
3532ClpPackedMatrix::setDimensions(int numrows, int numcols)
3533{
3534  matrix_->setDimensions(numrows,numcols);
3535}
3536/* Append a set of rows/columns to the end of the matrix. Returns number of errors
3537   i.e. if any of the new rows/columns contain an index that's larger than the
3538   number of columns-1/rows-1 (if numberOther>0) or duplicates
3539   If 0 then rows, 1 if columns */
3540int 
3541ClpPackedMatrix::appendMatrix(int number, int type,
3542                            const CoinBigIndex * starts, const int * index,
3543                            const double * element, int numberOther)
3544{
3545  int numberErrors=0;
3546  // make sure other dimension is big enough
3547  if (type==0) {
3548    // rows
3549    if (matrix_->isColOrdered()&&numberOther>matrix_->getNumCols())
3550      matrix_->setDimensions(-1,numberOther);
3551    if (!matrix_->isColOrdered()||numberOther>=0||matrix_->getExtraGap()) {
3552      numberErrors=matrix_->appendRows(number,starts,index,element,numberOther);
3553    } else {
3554      //CoinPackedMatrix mm(*matrix_);
3555      matrix_->appendMinorFast(number,starts,index,element);
3556      //mm.appendRows(number,starts,index,element,numberOther);
3557      //if (!mm.isEquivalent(*matrix_)) {
3558      //printf("bad append\n");
3559      //abort();
3560      //}
3561    } 
3562  } else {
3563    // columns
3564    if (!matrix_->isColOrdered()&&numberOther>matrix_->getNumRows())
3565      matrix_->setDimensions(numberOther,-1);
3566    numberErrors=matrix_->appendCols(number,starts,index,element,numberOther);
3567  }
3568  clearCopies();
3569  numberActiveColumns_ = matrix_->getNumCols();
3570  return numberErrors;
3571}
3572void
3573ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
3574{
3575  delete rowCopy_;
3576  rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
3577  // See if anything in it
3578  if (!rowCopy_->usefulInfo()) {
3579    delete rowCopy_;
3580    rowCopy_=NULL;
3581    flags_ &= ~4;
3582  } else {
3583    flags_ |= 4;
3584  }
3585}
3586void
3587ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
3588{
3589  delete columnCopy_;
3590  if ((flags_&16)!=0) {
3591    columnCopy_ = new ClpPackedMatrix3(model,matrix_);
3592    flags_ |= 8;
3593  } else {
3594    columnCopy_=NULL;
3595  }
3596}
3597// Say we don't want special column copy
3598void 
3599ClpPackedMatrix::releaseSpecialColumnCopy()
3600{ 
3601  flags_ &= ~(8+16);
3602  delete columnCopy_;
3603  columnCopy_=NULL;
3604}
3605// Correct sequence in and out to give true value
3606void 
3607ClpPackedMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) 
3608{
3609  if (columnCopy_) {
3610    if (sequenceIn!=-999) {
3611      if (sequenceIn!=sequenceOut) {
3612        if (sequenceIn<numberActiveColumns_)
3613          columnCopy_->swapOne(model,this,sequenceIn);
3614        if (sequenceOut<numberActiveColumns_)
3615          columnCopy_->swapOne(model,this,sequenceOut);
3616      }
3617    } else {
3618      // do all
3619      columnCopy_->sortBlocks(model);
3620    }
3621  }
3622}
3623// Check validity
3624void 
3625ClpPackedMatrix::checkFlags(int type) const
3626{
3627  int iColumn;
3628  // get matrix data pointers
3629  //const int * row = matrix_->getIndices();
3630  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3631  const int * columnLength = matrix_->getVectorLengths(); 
3632  const double * elementByColumn = matrix_->getElements();
3633  if (!zeros()) {
3634    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3635      CoinBigIndex j;
3636      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
3637        if (!elementByColumn[j])
3638          abort();
3639      }
3640    }
3641  }
3642  if ((flags_&2)==0) {
3643    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3644      if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
3645        abort();
3646      }
3647    }
3648  }
3649  if (type) {
3650    if ((flags_&2)!=0) {
3651      bool ok=true;
3652      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
3653        if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
3654          ok=false;
3655          break;
3656        }
3657      }
3658      if (ok)
3659        printf("flags_ could be 0\n");
3660    }
3661  }
3662}
3663//#############################################################################
3664// Constructors / Destructor / Assignment
3665//#############################################################################
3666
3667//-------------------------------------------------------------------
3668// Default Constructor
3669//-------------------------------------------------------------------
3670ClpPackedMatrix2::ClpPackedMatrix2 () 
3671  : numberBlocks_(0),
3672    numberRows_(0),
3673    offset_(NULL),
3674    count_(NULL),
3675    rowStart_(NULL),
3676    column_(NULL),
3677    work_(NULL)
3678{
3679#ifdef THREAD
3680  threadId_ = NULL;
3681  info_ = NULL;
3682#endif
3683}
3684//-------------------------------------------------------------------
3685// Useful Constructor
3686//-------------------------------------------------------------------
3687ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy) 
3688  : numberBlocks_(0),
3689    numberRows_(0),
3690    offset_(NULL),
3691    count_(NULL),
3692    rowStart_(NULL),
3693    column_(NULL),
3694    work_(NULL)
3695{
3696#ifdef THREAD
3697  threadId_ = NULL;
3698  info_ = NULL;
3699#endif
3700  numberRows_ = rowCopy->getNumRows();
3701  if (!numberRows_)
3702    return;
3703  int numberColumns = rowCopy->getNumCols();
3704  const int * column = rowCopy->getIndices();
3705  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3706  const int * length = rowCopy->getVectorLengths();
3707  const double * element = rowCopy->getElements();
3708  int chunk = 32768; // tune
3709  //chunk=100;
3710  // tune
3711#if 0
3712  int chunkY[7]={1024,2048,4096,8192,16384,32768,65535};
3713  int its = model->maximumIterations();
3714  if (its>=1000000&&its<1000999) {
3715    its -= 1000000;
3716    its = its/10;
3717    if (its>=7) abort();
3718    chunk=chunkY[its];
3719    printf("chunk size %d\n",chunk);
3720#define cpuid(func,ax,bx,cx,dx)\
3721    __asm__ __volatile__ ("cpuid":\
3722    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
3723    unsigned int a,b,c,d;
3724    int func=0;
3725    cpuid(func,a,b,c,d);
3726    {
3727      int i;
3728      unsigned int value;
3729      value=b;
3730      for (i=0;i<4;i++) {
3731        printf("%c",(value&0xff));
3732        value = value >>8;
3733      }
3734      value=d;
3735      for (i=0;i<4;i++) {
3736        printf("%c",(value&0xff));
3737        value = value >>8;
3738      }
3739      value=c;
3740      for (i=0;i<4;i++) {
3741        printf("%c",(value&0xff));
3742        value = value >>8;
3743      }
3744      printf("\n");
3745      int maxfunc = a;
3746      if (maxfunc>10) {
3747        printf("not intel?\n");
3748        abort();
3749      }
3750      for (func=1;func<=maxfunc;func++) {
3751        cpuid(func,a,b,c,d);
3752        printf("func %d, %x %x %x %x\n",func,a,b,c,d);
3753      }
3754    }
3755#else
3756    if (numberColumns>10000||chunk==100) {
3757#endif
3758  } else {
3759    //printf("no chunk\n");
3760    return;
3761  }
3762  // Could also analyze matrix to get natural breaks
3763  numberBlocks_ = (numberColumns+chunk-1)/chunk;
3764#ifdef THREAD
3765  // Get work areas
3766  threadId_ = new pthread_t [numberBlocks_];
3767  info_ = new dualColumn0Struct[numberBlocks_];
3768#endif
3769  // Even out
3770  chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
3771  offset_ = new int[numberBlocks_+1];
3772  offset_[numberBlocks_]=numberColumns;
3773  int nRow = numberBlocks_*numberRows_;
3774  count_ = new unsigned short[nRow];
3775  memset(count_,0,nRow*sizeof(unsigned short));
3776  rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
3777  CoinBigIndex nElement = rowStart[numberRows_];
3778  rowStart_[nRow+numberRows_]=nElement;
3779  column_ = new unsigned short [nElement];
3780  // assumes int <= double
3781  int sizeWork = 6*numberBlocks_;
3782  work_ = new double[sizeWork];;
3783  int iBlock;
3784  int nZero=0;
3785  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
3786    int start = iBlock*chunk;
3787    offset_[iBlock]=start;
3788    int end = start+chunk;
3789    for (int iRow=0;iRow<numberRows_;iRow++) {
3790      if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
3791        printf("not packed correctly - gaps\n");
3792        abort();
3793      }
3794      bool lastFound=false;
3795      int nFound=0;
3796      for (CoinBigIndex j = rowStart[iRow];
3797           j<rowStart[iRow]+length[iRow];j++) {
3798        int iColumn = column[j];
3799        if (iColumn>=start) {
3800          if (iColumn<end) {
3801            if (!element[j]) {
3802              printf("not packed correctly - zero element\n");
3803              abort();
3804            }
3805            column_[j]=static_cast<unsigned short>(iColumn-start);
3806            nFound++;
3807            if (lastFound) {
3808              printf("not packed correctly - out of order\n");
3809              abort();
3810            }
3811          } else {
3812            //can't find any more
3813            lastFound=true;
3814          }
3815        } 
3816      }
3817      count_[iRow*numberBlocks_+iBlock]=static_cast<unsigned short>(nFound);
3818      if (!nFound)
3819        nZero++;
3820    }
3821  }
3822  //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
3823  //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
3824}
3825
3826//-------------------------------------------------------------------
3827// Copy constructor
3828//-------------------------------------------------------------------
3829ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs) 
3830  : numberBlocks_(rhs.numberBlocks_),
3831    numberRows_(rhs.numberRows_)
3832{ 
3833  if (numberBlocks_) {
3834    offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3835    int nRow = numberBlocks_*numberRows_;
3836    count_ = CoinCopyOfArray(rhs.count_,nRow);
3837    rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3838    CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3839    column_ = CoinCopyOfArray(rhs.column_,nElement);
3840    int sizeWork = 6*numberBlocks_;
3841    work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3842#ifdef THREAD
3843    threadId_ = new pthread_t [numberBlocks_];
3844    info_ = new dualColumn0Struct[numberBlocks_];
3845#endif
3846  } else {
3847    offset_ = NULL;
3848    count_ = NULL;
3849    rowStart_ = NULL;
3850    column_ = NULL;
3851    work_ = NULL;
3852#ifdef THREAD
3853    threadId_ = NULL;
3854    info_ = NULL;
3855#endif
3856  }
3857}
3858//-------------------------------------------------------------------
3859// Destructor
3860//-------------------------------------------------------------------
3861ClpPackedMatrix2::~ClpPackedMatrix2 ()
3862{
3863  delete [] offset_;
3864  delete [] count_;
3865  delete [] rowStart_;
3866  delete [] column_;
3867  delete [] work_;
3868#ifdef THREAD
3869  delete [] threadId_;
3870  delete [] info_;
3871#endif
3872}
3873
3874//----------------------------------------------------------------
3875// Assignment operator
3876//-------------------------------------------------------------------
3877ClpPackedMatrix2 &
3878ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
3879{
3880  if (this != &rhs) {
3881    numberBlocks_ = rhs.numberBlocks_;
3882    numberRows_ = rhs.numberRows_;
3883    delete [] offset_;
3884    delete [] count_;
3885    delete [] rowStart_;
3886    delete [] column_;
3887    delete [] work_;
3888#ifdef THREAD
3889    delete [] threadId_;
3890    delete [] info_;
3891#endif
3892    if (numberBlocks_) {
3893      offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
3894      int nRow = numberBlocks_*numberRows_;
3895      count_ = CoinCopyOfArray(rhs.count_,nRow);
3896      rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
3897      CoinBigIndex nElement = rowStart_[nRow+numberRows_];
3898      column_ = CoinCopyOfArray(rhs.column_,nElement);
3899      int sizeWork = 6*numberBlocks_;
3900      work_ = CoinCopyOfArray(rhs.work_,sizeWork);
3901#ifdef THREAD
3902      threadId_ = new pthread_t [numberBlocks_];
3903      info_ = new dualColumn0Struct[numberBlocks_];
3904#endif
3905    } else {
3906      offset_ = NULL;
3907      count_ = NULL;
3908      rowStart_ = NULL;
3909      column_ = NULL;
3910      work_ = NULL;
3911#ifdef THREAD
3912      threadId_ = NULL;
3913      info_ = NULL;
3914#endif
3915    }
3916  }
3917  return *this;
3918}
3919static int dualColumn0(const ClpSimplex * model,double * spare,
3920                       int * spareIndex,const double * arrayTemp,
3921                       const int * indexTemp,int numberIn,
3922                       int offset, double acceptablePivot,double * bestPossiblePtr,
3923                       double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
3924{
3925  // do dualColumn0
3926  int i;
3927  int numberRemaining=0;
3928  double bestPossible=0.0;
3929  double upperTheta=1.0e31;
3930  double freePivot = acceptablePivot;
3931  int posFree=-1;
3932  const double * reducedCost=model->djRegion(1);
3933  double dualTolerance = model->dualTolerance();
3934  // We can also see if infeasible or pivoting on free
3935  double tentativeTheta = 1.0e25;
3936  for (i=0;i<numberIn;i++) {
3937    double alpha = arrayTemp[i];
3938    int iSequence = indexTemp[i]+offset;
3939    double oldValue;
3940    double value;
3941    bool keep;
3942   
3943    switch(model->getStatus(iSequence)) {
3944     
3945    case ClpSimplex::basic:
3946    case ClpSimplex::isFixed:
3947      break;
3948    case ClpSimplex::isFree:
3949    case ClpSimplex::superBasic:
3950      bestPossible = CoinMax(bestPossible,fabs(alpha));
3951      oldValue = reducedCost[iSequence];
3952      // If free has to be very large - should come in via dualRow
3953      if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
3954        break;
3955      if (oldValue>dualTolerance) {
3956        keep = true;
3957      } else if (oldValue<-dualTolerance) {
3958        keep = true;
3959      } else {
3960        if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
3961          keep = true;
3962        else
3963          keep=false;
3964      }
3965      if (keep) {
3966        // free - choose largest
3967        if (fabs(alpha)>freePivot) {
3968          freePivot=fabs(alpha);
3969          posFree = i;
3970        }
3971      }
3972      break;
3973    case ClpSimplex::atUpperBound:
3974      oldValue = reducedCost[iSequence];
3975      value = oldValue-tentativeTheta*alpha;
3976      //assert (oldValue<=dualTolerance*1.0001);
3977      if (value>dualTolerance) {
3978        bestPossible = CoinMax(bestPossible,-alpha);
3979        value = oldValue-upperTheta*alpha;
3980        if (value>dualTolerance && -alpha>=acceptablePivot) 
3981          upperTheta = (oldValue-dualTolerance)/alpha;
3982        // add to list
3983        spare[numberRemaining]=alpha;
3984        spareIndex[numberRemaining++]=iSequence;
3985      }
3986      break;
3987    case ClpSimplex::atLowerBound:
3988      oldValue = reducedCost[iSequence];
3989      value = oldValue-tentativeTheta*alpha;
3990      //assert (oldValue>=-dualTolerance*1.0001);
3991      if (value<-dualTolerance) {
3992        bestPossible = CoinMax(bestPossible,alpha);
3993        value = oldValue-upperTheta*alpha;
3994        if (value<-dualTolerance && alpha>=acceptablePivot) 
3995          upperTheta = (oldValue+dualTolerance)/alpha;
3996        // add to list
3997        spare[numberRemaining]=alpha;
3998        spareIndex[numberRemaining++]=iSequence;
3999      }
4000      break;
4001    }
4002  }
4003  *bestPossiblePtr = bestPossible;
4004  *upperThetaPtr = upperTheta;
4005  *freePivotPtr = freePivot;
4006  *posFreePtr = posFree;
4007  return numberRemaining;
4008}
4009static int doOneBlock(double * array, int * index, 
4010                      const double * pi,const CoinBigIndex * rowStart, const double * element,
4011                      const unsigned short * column, int numberInRowArray, int numberLook)
4012{
4013  int iWhich=0;
4014  int nextN=0;
4015  CoinBigIndex nextStart=0;
4016  double nextPi=0.0;
4017  for (;iWhich<numberInRowArray;iWhich++) {
4018    nextStart = rowStart[0];
4019    nextN = rowStart[numberInRowArray]-nextStart;
4020    rowStart++;
4021    if (nextN) {
4022      nextPi = pi[iWhich];
4023      break;
4024    }
4025  }
4026  int i;
4027  while (iWhich<numberInRowArray) {
4028    double value = nextPi;
4029    CoinBigIndex j=  nextStart;
4030    int n = nextN;
4031    // get next
4032    iWhich++;
4033    for (;iWhich<numberInRowArray;iWhich++) {
4034      nextStart = rowStart[0];
4035      nextN = rowStart[numberInRowArray]-nextStart;
4036      rowStart++;
4037      if (nextN) {
4038        coin_prefetch(element+nextStart);
4039        nextPi = pi[iWhich];
4040        break;
4041      }
4042    }
4043    CoinBigIndex end =j+n;
4044    //coin_prefetch(element+rowStart_[i+1]);
4045    //coin_prefetch(column_+rowStart_[i+1]);
4046    if (n<100) {
4047      if ((n&1)!=0) {
4048        unsigned int jColumn = column[j];
4049        array[jColumn] -= value*element[j];
4050        j++;
4051      }     
4052      coin_prefetch(column+nextStart);
4053      for (;j<end;j+=2) {
4054        unsigned int jColumn0 = column[j];
4055        double value0 = value*element[j];
4056        unsigned int jColumn1 = column[j+1];
4057        double value1 = value*element[j+1];
4058        array[jColumn0] -= value0;
4059        array[jColumn1] -= value1;
4060      }
4061    } else {
4062      if ((n&1)!=0) {
4063        unsigned int jColumn = column[j];
4064        array[jColumn] -= value*element[j];
4065        j++;
4066      }     
4067      if ((n&2)!=0) {
4068        unsigned int jColumn0 = column[j];
4069        double value0 = value*element[j];
4070        unsigned int jColumn1 = column[j+1];
4071        double value1 = value*element[j+1];
4072        array[jColumn0] -= value0;
4073        array[jColumn1] -= value1;
4074        j+=2;
4075      }     
4076      if ((n&4)!=0) {
4077        unsigned int jColumn0 = column[j];
4078        double value0 = value*element[j];
4079        unsigned int jColumn1 = column[j+1];
4080        double value1 = value*element[j+1];
4081        unsigned int jColumn2 = column[j+2];
4082        double value2 = value*element[j+2];
4083        unsigned int jColumn3 = column[j+3];
4084        double value3 = value*element[j+3];
4085        array[jColumn0] -= value0;
4086        array[jColumn1] -= value1;
4087        array[jColumn2] -= value2;
4088        array[jColumn3] -= value3;
4089        j+=4;
4090      }
4091      //coin_prefetch(column+nextStart);
4092      for (;j<end;j+=8) {
4093        coin_prefetch(element+j+16);
4094        unsigned int jColumn0 = column[j];
4095        double value0 = value*element[j];
4096        unsigned int jColumn1 = column[j+1];
4097        double value1 = value*element[j+1];
4098        unsigned int jColumn2 = column[j+2];
4099        double value2 = value*element[j+2];
4100        unsigned int jColumn3 = column[j+3];
4101        double value3 = value*element[j+3];
4102        array[jColumn0] -= value0;
4103        array[jColumn1] -= value1;
4104        array[jColumn2] -= value2;
4105        array[jColumn3] -= value3;
4106        coin_prefetch(column+j+16);
4107        jColumn0 = column[j+4];
4108        value0 = value*element[j+4];
4109        jColumn1 = column[j+5];
4110        value1 = value*element[j+5];
4111        jColumn2 = column[j+6];
4112        value2 = value*element[j+6];
4113        jColumn3 = column[j+7];
4114        value3 = value*element[j+7];
4115        array[jColumn0] -= value0;
4116        array[jColumn1] -= value1;
4117        array[jColumn2] -= value2;
4118        array[jColumn3] -= value3;
4119      }
4120    }
4121  }
4122  // get rid of tiny values
4123  int nSmall = numberLook;
4124  int numberNonZero=0;
4125  for (i=0;i<nSmall;i++) {
4126    double value = array[i];
4127    array[i]=0.0; 
4128    if (fabs(value)>1.0e-12) {
4129      array[numberNonZero]=value;
4130      index[numberNonZero++]=i;
4131    }
4132  }
4133  for (;i<numberLook;i+=4) {
4134    double value0 = array[i+0];
4135    double value1 = array[i+1];
4136    double value2 = array[i+2];
4137    double value3 = array[i+3];
4138    array[i+0]=0.0; 
4139    array[i+1]=0.0; 
4140    array[i+2]=0.0; 
4141    array[i+3]=0.0; 
4142    if (fabs(value0)>1.0e-12) {
4143      array[numberNonZero]=value0;
4144      index[numberNonZero++]=i+0;
4145    }
4146    if (fabs(value1)>1.0e-12) {
4147      array[numberNonZero]=value1;
4148      index[numberNonZero++]=i+1;
4149    }
4150    if (fabs(value2)>1.0e-12) {
4151      array[numberNonZero]=value2;
4152      index[numberNonZero++]=i+2;
4153    }
4154    if (fabs(value3)>1.0e-12) {
4155      array[numberNonZero]=value3;
4156      index[numberNonZero++]=i+3;
4157    }
4158  }
4159  return numberNonZero;
4160}
4161#ifdef THREAD
4162static void * doOneBlockThread(void * voidInfo)
4163{
4164  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
4165  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
4166                                  info->rowStart,info->element,info->column,
4167                                  info->numberInRowArray,info->numberLook);
4168  return NULL;
4169}
4170static void * doOneBlockAnd0Thread(void * voidInfo)
4171{
4172  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
4173  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
4174                                  info->rowStart,info->element,info->column,
4175                                  info->numberInRowArray,info->numberLook);
4176  *(info->numberOutPtr) =  dualColumn0(info->model,info->spare,
4177                                    info->spareIndex,(const double *)info->arrayTemp,
4178                                    (const int *) info->indexTemp,*(info->numberInPtr),
4179                                    info->offset, info->acceptablePivot,info->bestPossiblePtr,
4180                                    info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
4181  return NULL;
4182}
4183#endif
4184/* Return <code>x * scalar * A in <code>z</code>.
4185   Note - x packed and z will be packed mode
4186   Squashes small elements and knows about ClpSimplex */
4187void 
4188ClpPackedMatrix2::transposeTimes(const ClpSimplex * model, 
4189                                 const CoinPackedMatrix * rowCopy,
4190                                 const CoinIndexedVector * rowArray,
4191                                 CoinIndexedVector * spareArray,
4192                                 CoinIndexedVector * columnArray) const
4193{
4194  // See if dualColumn0 coding wanted
4195  bool dualColumn = model->spareIntArray_[0]==1;
4196  double acceptablePivot = model->spareDoubleArray_[0];
4197  double bestPossible=0.0;
4198  double upperTheta=1.0e31;
4199  double freePivot = acceptablePivot;
4200  int posFree=-1;
4201  int numberRemaining=0;
4202  //if (model->numberIterations()>=200000) {
4203  //printf("time %g\n",CoinCpuTime()-startTime);
4204  //exit(77);
4205  //}
4206  double * pi = rowArray->denseVector();
4207  int numberNonZero=0;
4208  int * index = columnArray->getIndices();
4209  double * array = columnArray->denseVector();
4210  int numberInRowArray = rowArray->getNumElements();
4211  const int * whichRow = rowArray->getIndices();
4212  double * element = const_cast<double *>(rowCopy->getElements());
4213  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4214  int i;
4215  CoinBigIndex * rowStart2 = rowStart_;
4216  if (!dualColumn) {
4217    for (i=0;i<numberInRowArray;i++) {
4218      int iRow = whichRow[i]; 
4219      CoinBigIndex start = rowStart[iRow];
4220      *rowStart2 = start;
4221      unsigned short * count1 = count_+iRow*numberBlocks_;
4222      int put=0;
4223      for (int j=0;j<numberBlocks_;j++) {
4224        put += numberInRowArray;
4225        start += count1[j];
4226        rowStart2[put]=start;
4227      }
4228      rowStart2 ++;
4229    }
4230  } else {
4231    // also do dualColumn stuff
4232    double * spare = spareArray->denseVector();
4233    int * spareIndex = spareArray->getIndices();
4234    const double * reducedCost=model->djRegion(0);
4235    double dualTolerance = model->dualTolerance();
4236    // We can also see if infeasible or pivoting on free
4237    double tentativeTheta = 1.0e25;
4238    int addSequence=model->numberColumns();
4239    for (i=0;i<numberInRowArray;i++) {
4240      int iRow = whichRow[i]; 
4241      double alpha=pi[i];
4242      double oldValue;
4243      double value;
4244      bool keep;
4245
4246      switch(model->getStatus(iRow+addSequence)) {
4247         
4248      case ClpSimplex::basic:
4249      case ClpSimplex::isFixed:
4250        break;
4251      case ClpSimplex::isFree:
4252      case ClpSimplex::superBasic:
4253        bestPossible = CoinMax(bestPossible,fabs(alpha));
4254        oldValue = reducedCost[iRow];
4255        // If free has to be very large - should come in via dualRow
4256        if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
4257          break;
4258        if (oldValue>dualTolerance) {
4259          keep = true;
4260        } else if (oldValue<-dualTolerance) {
4261          keep = true;
4262        } else {
4263          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) 
4264            keep = true;
4265          else
4266            keep=false;
4267        }
4268        if (keep) {
4269          // free - choose largest
4270          if (fabs(alpha)>freePivot) {
4271            freePivot=fabs(alpha);
4272            posFree = i + addSequence;
4273          }
4274        }
4275        break;
4276      case ClpSimplex::atUpperBound:
4277        oldValue = reducedCost[iRow];
4278        value = oldValue-tentativeTheta*alpha;
4279        //assert (oldValue<=dualTolerance*1.0001);
4280        if (value>dualTolerance) {
4281          bestPossible = CoinMax(bestPossible,-alpha);
4282          value = oldValue-upperTheta*alpha;
4283          if (value>dualTolerance && -alpha>=acceptablePivot) 
4284            upperTheta = (oldValue-dualTolerance)/alpha;
4285          // add to list
4286          spare[numberRemaining]=alpha;
4287          spareIndex[numberRemaining++]=iRow+addSequence;
4288        }
4289        break;
4290      case ClpSimplex::atLowerBound:
4291        oldValue = reducedCost[iRow];
4292        value = oldValue-tentativeTheta*alpha;
4293        //assert (oldValue>=-dualTolerance*1.0001);
4294        if (value<-dualTolerance) {
4295          bestPossible = CoinMax(bestPossible,alpha);
4296          value = oldValue-upperTheta*alpha;
4297          if (value<-dualTolerance && alpha>=acceptablePivot) 
4298            upperTheta = (oldValue+dualTolerance)/alpha;
4299          // add to list
4300          spare[numberRemaining]=alpha;
4301          spareIndex[numberRemaining++]=iRow+addSequence;
4302        }
4303        break;
4304      }
4305      CoinBigIndex start = rowStart[iRow];
4306      *rowStart2 = start;
4307      unsigned short * count1 = count_+iRow*numberBlocks_;
4308      int put=0;
4309      for (int j=0;j<numberBlocks_;j++) {
4310        put += numberInRowArray;
4311        start += count1[j];
4312        rowStart2[put]=start;
4313      }
4314      rowStart2 ++;
4315    }
4316  }
4317  double * spare = spareArray->denseVector();
4318  int * spareIndex = spareArray->getIndices();
4319  int saveNumberRemaining=numberRemaining;
4320  int iBlock;
4321  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
4322    double * dwork = work_+6*iBlock;
4323    int * iwork = (int *) (dwork+3);
4324    if (!dualColumn) {
4325#ifndef THREAD
4326      int offset = offset_[iBlock];
4327      int offset3 = offset;
4328      offset = numberNonZero;
4329      double * arrayTemp = array+offset;
4330      int * indexTemp = index+offset; 
4331      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
4332                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
4333      int number = iwork[0];
4334      for (i=0;i<number;i++) {
4335        //double value = arrayTemp[i];
4336        //arrayTemp[i]=0.0;
4337        //array[numberNonZero]=value;
4338        index[numberNonZero++]=indexTemp[i]+offset3;
4339      }
4340#else
4341      int offset = offset_[iBlock];
4342      double * arrayTemp = array+offset;
4343      int * indexTemp = index+offset; 
4344      dualColumn0Struct * infoPtr = info_+iBlock;
4345      infoPtr->arrayTemp=arrayTemp;
4346      infoPtr->indexTemp=indexTemp;
4347      infoPtr->numberInPtr=&iwork[0];
4348      infoPtr->pi=pi;
4349      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
4350      infoPtr->element=element;
4351      infoPtr->column=column_;
4352      infoPtr->numberInRowArray=numberInRowArray;
4353      infoPtr->numberLook=offset_[iBlock+1]-offset;
4354      pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
4355#endif
4356    } else {
4357#ifndef THREAD
4358      int offset = offset_[iBlock];
4359      // allow for already saved
4360      int offset2 = offset + saveNumberRemaining;
4361      int offset3 = offset;
4362      offset = numberNonZero;
4363      offset2 = numberRemaining;
4364      double * arrayTemp = array+offset;
4365      int * indexTemp = index+offset; 
4366      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
4367                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
4368      iwork[1] = dualColumn0(model,spare+offset2,
4369                                   spareIndex+offset2,
4370                                   arrayTemp,indexTemp,
4371                                   iwork[0],offset3,acceptablePivot,
4372                                   &dwork[0],&dwork[1],&iwork[2],
4373                                   &dwork[2]);
4374      int number = iwork[0];
4375      int numberLook = iwork[1];
4376#if 1
4377      numberRemaining += numberLook;
4378#else
4379      double * spareTemp = spare+offset2;
4380      const int * spareIndexTemp = spareIndex+offset2; 
4381      for (i=0;i<numberLook;i++) {
4382        double value = spareTemp[i];
4383        spareTemp[i]=0.0;
4384        spare[numberRemaining] = value;
4385        spareIndex[numberRemaining++] = spareIndexTemp[i];
4386      }
4387#endif
4388      if (dwork[2]>freePivot) {
4389        freePivot=dwork[2];
4390        posFree = iwork[2]+numberNonZero;
4391      } 
4392      upperTheta =  CoinMin(dwork[1],upperTheta);
4393      bestPossible = CoinMax(dwork[0],bestPossible);
4394      for (i=0;i<number;i++) {
4395        // double value = arrayTemp[i];
4396        //arrayTemp[i]=0.0;
4397        //array[numberNonZero]=value;
4398        index[numberNonZero++]=indexTemp[i]+offset3;
4399      }
4400#else
4401      int offset = offset_[iBlock];
4402      // allow for already saved
4403      int offset2 = offset + saveNumberRemaining;
4404      double * arrayTemp = array+offset;
4405      int * indexTemp = index+offset; 
4406      dualColumn0Struct * infoPtr = info_+iBlock;
4407      infoPtr->model=model;
4408      infoPtr->spare=spare+offset2;
4409      infoPtr->spareIndex=spareIndex+offset2;
4410      infoPtr->arrayTemp=arrayTemp;
4411      infoPtr->indexTemp=indexTemp;
4412      infoPtr->numberInPtr=&iwork[0];
4413      infoPtr->offset=offset;
4414      infoPtr->acceptablePivot=acceptablePivot;
4415      infoPtr->bestPossiblePtr=&dwork[0];
4416      infoPtr->upperThetaPtr=&dwork[1];
4417      infoPtr->posFreePtr=&iwork[2];
4418      infoPtr->freePivotPtr=&dwork[2];
4419      infoPtr->numberOutPtr=&iwork[1];
4420      infoPtr->pi=pi;
4421      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
4422      infoPtr->element=element;
4423      infoPtr->column=column_;
4424      infoPtr->numberInRowArray=numberInRowArray;
4425      infoPtr->numberLook=offset_[iBlock+1]-offset;
4426      if (iBlock>=2)
4427        pthread_join(threadId_[iBlock-2],NULL);
4428      pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
4429      //pthread_join(threadId_[iBlock],NULL);
4430#endif
4431    }
4432  }
4433  for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
4434#ifdef THREAD
4435    pthread_join(threadId_[iBlock],NULL);
4436#endif
4437  }
4438#ifdef THREAD
4439  for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
4440    //pthread_join(threadId_[iBlock],NULL);
4441    int offset = offset_[iBlock];
4442    double * dwork = work_+6*iBlock;
4443    int * iwork = (int *) (dwork+3);
4444    int number = iwork[0];
4445    if (dualColumn) {
4446      // allow for already saved
4447      int offset2 = offset + saveNumberRemaining;
4448      int numberLook = iwork[1];
4449      double * spareTemp = spare+offset2;
4450      const int * spareIndexTemp = spareIndex+offset2; 
4451      for (i=0;i<numberLook;i++) {
4452        double value = spareTemp[i];
4453        spareTemp[i]=0.0;
4454        spare[numberRemaining] = value;
4455        spareIndex[numberRemaining++] = spareIndexTemp[i];
4456      }
4457      if (dwork[2]>freePivot) {
4458        freePivot=dwork[2];
4459        posFree = iwork[2]+numberNonZero;
4460      } 
4461      upperTheta =  CoinMin(dwork[1],upperTheta);
4462      bestPossible = CoinMax(dwork[0],bestPossible);
4463    }
4464    double * arrayTemp = array+offset;
4465    const int * indexTemp = index+offset; 
4466    for (i=0;i<number;i++) {
4467      double value = arrayTemp[i];
4468      arrayTemp[i]=0.0; 
4469      array[numberNonZero]=value;
4470      index[numberNonZero++]=indexTemp[i]+offset;
4471    }
4472  }
4473#endif
4474  columnArray->setNumElements(numberNonZero);
4475  columnArray->setPackedMode(true);
4476  if (dualColumn) {
4477    model->spareDoubleArray_[0]=upperTheta;
4478    model->spareDoubleArray_[1]=bestPossible;
4479    // and theta and alpha and sequence
4480    if (posFree<0) {
4481      model->spareIntArray_[1]=-1;
4482    } else {
4483      const double * reducedCost=model->djRegion(0);
4484      double alpha;
4485      int numberColumns=model->numberColumns();
4486      if (posFree<numberColumns) {
4487        alpha = columnArray->denseVector()[posFree];
4488        posFree = columnArray->getIndices()[posFree];
4489      } else {
4490        alpha = rowArray->denseVector()[posFree-numberColumns];
4491        posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
4492      }
4493      model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
4494      model->spareDoubleArray_[3]=alpha;
4495      model->spareIntArray_[1]=posFree;
4496    }
4497    spareArray->setNumElements(numberRemaining);
4498    // signal done
4499    model->spareIntArray_[0]=-1;
4500  }
4501}
4502/* Default constructor. */
4503ClpPackedMatrix3::ClpPackedMatrix3()
4504  : numberBlocks_(0),
4505    numberColumns_(0),
4506    column_(NULL),
4507    start_(NULL),
4508    row_(NULL),
4509    element_(NULL),
4510    block_(NULL)
4511{
4512}
4513/* Constructor from copy. */
4514ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy)
4515  : numberBlocks_(0),
4516    numberColumns_(0),
4517    column_(NULL),
4518    start_(NULL),
4519    row_(NULL),
4520    element_(NULL),
4521    block_(NULL)
4522{
4523#define MINBLOCK 6
4524#define MAXBLOCK 100
4525#define MAXUNROLL 10
4526  numberColumns_ = model->getNumCols();
4527  int numberColumns = columnCopy->getNumCols();
4528  assert (numberColumns_ >= numberColumns);
4529  int numberRows = columnCopy->getNumRows();
4530  int * counts = new int[numberRows+1];
4531  CoinZeroN(counts,numberRows+1);
4532  CoinBigIndex nels=0;
4533  CoinBigIndex nZeroEl=0;
4534  int iColumn;
4535  // get matrix data pointers
4536  const int * row = columnCopy->getIndices();
4537  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4538  const int * columnLength = columnCopy->getVectorLengths(); 
4539  const double * elementByColumn = columnCopy->getElements();
4540  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4541    CoinBigIndex start = columnStart[iColumn];
4542    int n = columnLength[iColumn];
4543    CoinBigIndex end = start + n;
4544    int kZero=0;
4545    for (CoinBigIndex j=start;j<end;j++) {
4546      if(!elementByColumn[j])
4547        kZero++;
4548    }
4549    n -= kZero;
4550    nZeroEl += kZero;
4551    nels += n;
4552    counts[n]++;
4553  }
4554  counts[0] += numberColumns_-numberColumns;
4555  int nZeroColumns = counts[0];
4556  counts[0]=-1;
4557  numberColumns_ -= nZeroColumns;
4558  column_ = new int [2*numberColumns_+nZeroColumns];
4559  int * lookup = column_ + numberColumns_;
4560  row_ = new int[nels];
4561  element_ = new double[nels];
4562  int nOdd=0;
4563  CoinBigIndex nInOdd=0;
4564  int i;
4565  for (i=1;i<=numberRows;i++) {
4566    int n = counts[i];
4567    if (n) {
4568      if (n<MINBLOCK||i>MAXBLOCK) {
4569        nOdd += n;
4570        counts[i]=-1;
4571        nInOdd += n*i;
4572      } else {
4573        numberBlocks_++;
4574      }
4575    } else {
4576      counts[i]=-1;
4577    }
4578  }
4579  start_ = new CoinBigIndex [nOdd+1];
4580  // even if no blocks do a dummy one
4581  numberBlocks_ = CoinMax(numberBlocks_,1);
4582  block_ = (blockStruct *) new char [numberBlocks_*sizeof(blockStruct)];
4583  memset(block_,0,numberBlocks_*sizeof(blockStruct));
4584  // Fill in what we can
4585  int nTotal=nOdd;
4586  // in case no blocks
4587  block_->startIndices_=nTotal;
4588  nels=nInOdd;
4589  int nBlock=0;
4590  for (i=0;i<=CoinMin(MAXBLOCK,numberRows);i++) {
4591    if (counts[i]>0) {
4592      blockStruct * block = block_ + nBlock;
4593      int n=counts[i];
4594      counts[i]= nBlock; // backward pointer
4595      nBlock++;
4596      block->startIndices_ = nTotal;
4597      block->startElements_ = nels;
4598      block->numberElements_ = i;
4599      // up counts
4600      nTotal += n;
4601      nels += n*i;
4602    }
4603  }
4604  for (iColumn=numberColumns;iColumn<numberColumns_;iColumn++) 
4605    lookup[iColumn]=-1;
4606  // fill
4607  start_[0]=0;
4608  nOdd=0;
4609  nInOdd=0;
4610  const double * columnScale = model->columnScale();
4611  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4612    CoinBigIndex start = columnStart[iColumn];
4613    int n = columnLength[iColumn];
4614    CoinBigIndex end = start + n;
4615    int kZero=0;
4616    for (CoinBigIndex j=start;j<end;j++) {
4617      if(!elementByColumn[j])
4618        kZero++;
4619    }
4620    n -= kZero;
4621    if (n) {
4622      int iBlock = counts[n];
4623      if (iBlock>=0) {
4624        blockStruct * block = block_ + iBlock;
4625        int k = block->numberInBlock_;
4626        block->numberInBlock_ ++;
4627        column_[block->startIndices_+k]=iColumn;
4628        lookup[iColumn]=k;
4629        CoinBigIndex put = block->startElements_+k*n;
4630        for (CoinBigIndex j=start;j<end;j++) {
4631          double value = elementByColumn[j];
4632          if(value) {
4633            if (columnScale)
4634              value *= columnScale[iColumn];
4635            element_[put] = value;
4636            row_[put++]=row[j];
4637          }
4638        }
4639      } else {
4640        // odd ones
4641        for (CoinBigIndex j=start;j<end;j++) {
4642          double value = elementByColumn[j];
4643          if(value) {
4644            if (columnScale)
4645              value *= columnScale[iColumn];
4646            element_[nInOdd] = value;
4647            row_[nInOdd++]=row[j];
4648          }
4649        }
4650        column_[nOdd]=iColumn;
4651        lookup[iColumn]=-1;
4652        nOdd++;
4653        start_[nOdd]=nInOdd;
4654      }
4655    } else {
4656      // zero column
4657      lookup[iColumn]=-1;
4658    }
4659  }
4660  delete [] counts;
4661}
4662/* Destructor */
4663ClpPackedMatrix3::~ClpPackedMatrix3()
4664{
4665  delete [] column_;
4666  delete [] start_;
4667  delete [] row_;
4668  delete [] element_;
4669  delete [] block_;
4670}
4671/* The copy constructor. */
4672ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
4673  : numberBlocks_(rhs.numberBlocks_),
4674    numberColumns_(rhs.numberColumns_),
4675    column_(NULL),
4676    start_(NULL),
4677    row_(NULL),
4678    element_(NULL),
4679    block_(NULL)
4680{
4681  if (rhs.numberBlocks_) {
4682    block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
4683    column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
4684    int numberOdd = block_->startIndices_;
4685    start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
4686    blockStruct * lastBlock = block_ + (numberBlocks_-1);
4687    CoinBigIndex numberElements = lastBlock->startElements_+
4688      lastBlock->numberInBlock_*lastBlock->numberElements_;
4689    row_ = CoinCopyOfArray(rhs.row_,numberElements);
4690    element_ = CoinCopyOfArray(rhs.element_,numberElements);
4691  }
4692}
4693ClpPackedMatrix3& 
4694ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
4695{
4696  if (this != &rhs) {
4697    delete [] column_;
4698    delete [] start_;
4699    delete [] row_;
4700    delete [] element_;
4701    delete [] block_;
4702    numberBlocks_ = rhs.numberBlocks_;
4703    numberColumns_ = rhs.numberColumns_;
4704    if (rhs.numberBlocks_) {
4705      block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
4706      column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
4707      int numberOdd = block_->startIndices_;
4708      start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
4709      blockStruct * lastBlock = block_ + (numberBlocks_-1);
4710      CoinBigIndex numberElements = lastBlock->startElements_+
4711        lastBlock->numberInBlock_*lastBlock->numberElements_;
4712      row_ = CoinCopyOfArray(rhs.row_,numberElements);
4713      element_ = CoinCopyOfArray(rhs.element_,numberElements);
4714    } else {
4715      column_ = NULL;
4716      start_ = NULL;
4717      row_ = NULL;
4718      element_ = NULL;
4719      block_ = NULL;
4720    }
4721  }
4722  return *this;
4723}
4724/* Sort blocks */
4725void 
4726ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
4727{
4728  int * lookup = column_ + numberColumns_;
4729  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4730    blockStruct * block = block_ + iBlock;
4731    int numberInBlock = block->numberInBlock_; 
4732    int nel = block->numberElements_;
4733    int * row = row_ + block->startElements_;
4734    double * element = element_ + block->startElements_;
4735    int * column = column_ + block->startIndices_;
4736    int lastPrice=0;
4737    int firstNotPrice=numberInBlock-1;
4738    while (lastPrice<=firstNotPrice) {
4739      // find first basic or fixed
4740      int iColumn=numberInBlock;
4741      for (;lastPrice<=firstNotPrice;lastPrice++) {
4742        iColumn = column[lastPrice];
4743        if (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4744            model->getColumnStatus(iColumn)==ClpSimplex::isFixed)
4745          break;
4746      }
4747      // find last non basic or fixed
4748      int jColumn=-1;
4749      for (;firstNotPrice>lastPrice;firstNotPrice--) {
4750        jColumn = column[firstNotPrice];
4751        if (model->getColumnStatus(jColumn)!=ClpSimplex::basic&&
4752            model->getColumnStatus(jColumn)!=ClpSimplex::isFixed)
4753          break;
4754      }
4755      if (firstNotPrice>lastPrice) {
4756        assert (column[lastPrice]==iColumn);
4757        assert (column[firstNotPrice]==jColumn);
4758        // need to swap
4759        column[firstNotPrice]=iColumn;
4760        lookup[iColumn]=firstNotPrice;
4761        column[lastPrice]=jColumn;
4762        lookup[jColumn]=lastPrice;
4763        double * elementA = element + lastPrice*nel;
4764        int * rowA = row + lastPrice*nel;
4765        double * elementB = element + firstNotPrice*nel;
4766        int * rowB = row + firstNotPrice*nel;
4767        for (int i=0;i<nel;i++) {
4768          int temp = rowA[i];
4769          double tempE = elementA[i];
4770          rowA[i]=rowB[i];
4771          elementA[i]=elementB[i];
4772          rowB[i]=temp;
4773          elementB[i]=tempE;
4774        }
4775        firstNotPrice--;
4776        lastPrice++;
4777      } else if (lastPrice==firstNotPrice) {
4778        // make sure correct side
4779        iColumn = column[lastPrice];
4780        if (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4781            model->getColumnStatus(iColumn)!=ClpSimplex::isFixed)
4782          lastPrice++;
4783        break;
4784      }
4785    }
4786    block->numberPrice_=lastPrice;
4787#ifndef NDEBUG
4788    // check
4789    int i;
4790    for (i=0;i<lastPrice;i++) {
4791      int iColumn = column[i];
4792      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4793              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
4794      assert (lookup[iColumn]==i);
4795    }
4796    for (;i<numberInBlock;i++) {
4797      int iColumn = column[i];
4798      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4799              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4800      assert (lookup[iColumn]==i);
4801    }
4802#endif
4803  }
4804}
4805// Swap one variable
4806void 
4807ClpPackedMatrix3::swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
4808                          int iColumn)
4809{
4810  int * lookup = column_ + numberColumns_;
4811  // position in block
4812  int kA = lookup[iColumn];
4813  if (kA<0)
4814    return; // odd one
4815  // get matrix data pointers
4816  const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
4817  //const int * row = columnCopy->getIndices();
4818  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4819  const int * columnLength = columnCopy->getVectorLengths(); 
4820  const double * elementByColumn = columnCopy->getElements();
4821  CoinBigIndex start = columnStart[iColumn];
4822  int n = columnLength[iColumn];
4823  if (matrix->zeros()) {
4824    CoinBigIndex end = start + n;
4825    for (CoinBigIndex j=start;j<end;j++) {
4826      if(!elementByColumn[j])
4827        n--;
4828    }
4829  }
4830  // find block - could do binary search
4831  int iBlock=CoinMin(n,numberBlocks_)-1;
4832  while (block_[iBlock].numberElements_!=n)
4833    iBlock--;
4834  blockStruct * block = block_ + iBlock;
4835  int nel = block->numberElements_;
4836  int * row = row_ + block->startElements_;
4837  double * element = element_ + block->startElements_;
4838  int * column = column_ + block->startIndices_;
4839  assert (column[kA]==iColumn);
4840  bool moveUp = (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4841                 model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4842  int lastPrice=block->numberPrice_;
4843  int kB;
4844  if (moveUp) {
4845    // May already be in correct place (e.g. fixed basic leaving basis)
4846    if (kA>=lastPrice)
4847      return;
4848    kB=lastPrice-1;
4849    block->numberPrice_--;
4850  } else {
4851    assert (kA>=lastPrice);
4852    kB=lastPrice;
4853    block->numberPrice_++;
4854  }
4855  int jColumn = column[kB];
4856  column[kA]=jColumn;
4857  lookup[jColumn]=kA;
4858  column[kB]=iColumn;
4859  lookup[iColumn]=kB;
4860  double * elementA = element + kB*nel;
4861  int * rowA = row + kB*nel;
4862  double * elementB = element + kA*nel;
4863  int * rowB = row + kA*nel;
4864  int i;
4865  for (i=0;i<nel;i++) {
4866    int temp = rowA[i];
4867    double tempE = elementA[i];
4868    rowA[i]=rowB[i];
4869    elementA[i]=elementB[i];
4870    rowB[i]=temp;
4871    elementB[i]=tempE;
4872  }
4873#ifndef NDEBUG
4874  // check
4875  for (i=0;i<block->numberPrice_;i++) {
4876    int iColumn = column[i];
4877    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
4878      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
4879              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
4880    assert (lookup[iColumn]==i);
4881  }
4882  int numberInBlock = block->numberInBlock_; 
4883  for (;i<numberInBlock;i++) {
4884    int iColumn = column[i];
4885    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
4886      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
4887              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4888    assert (lookup[iColumn]==i);
4889  }
4890#endif
4891}
4892/* Return <code>x * -1 * A in <code>z</code>.
4893   Note - x packed and z will be packed mode
4894   Squashes small elements and knows about ClpSimplex */
4895void 
4896ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
4897                                 const double * pi,
4898                                 CoinIndexedVector * output) const
4899{
4900  int numberNonZero=0;
4901  int * index = output->getIndices();
4902  double * array = output->denseVector();
4903  double zeroTolerance = model->factorization()->zeroTolerance();
4904  double value = 0.0;
4905  CoinBigIndex j;
4906  int numberOdd = block_->startIndices_;
4907  if (numberOdd) {
4908    // A) as probably long may be worth unrolling
4909    CoinBigIndex end = start_[1];
4910    for (j=start_[0]; j<end;j++) {
4911      int iRow = row_[j];
4912      value += pi[iRow]*element_[j];
4913    }
4914    int iColumn;
4915    // int jColumn=column_[0];
4916   
4917    for (iColumn=0;iColumn<numberOdd-1;iColumn++) {
4918      CoinBigIndex start = end;
4919      end = start_[iColumn+2];
4920      if (fabs(value)>zeroTolerance) {
4921        array[numberNonZero]=value;
4922        index[numberNonZero++]=column_[iColumn];
4923        //index[numberNonZero++]=jColumn;
4924      }
4925      // jColumn = column_[iColumn+1];
4926      value = 0.0;
4927      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
4928        for (j=start; j<end;j++) {
4929          int iRow = row_[j];
4930          value += pi[iRow]*element_[j];
4931        }
4932        //}
4933    }
4934    if (fabs(value)>zeroTolerance) {
4935      array[numberNonZero]=value;
4936      index[numberNonZero++]=column_[iColumn];
4937      //index[numberNonZero++]=jColumn;
4938    }
4939  }
4940  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
4941    // B) Can sort so just do nonbasic (and nonfixed)
4942    // C) Can do two at a time (if so put odd one into start_)
4943    // D) can use switch
4944    blockStruct * block = block_ + iBlock;
4945    //int numberPrice = block->numberInBlock_;
4946    int numberPrice = block->numberPrice_; 
4947    int nel = block->numberElements_;
4948    int * row = row_ + block->startElements_;
4949    double * element = element_ + block->startElements_;
4950    int * column = column_ + block->startIndices_;
4951#if 0
4952    // two at a time
4953    if ((numberPrice&1)!=0) {
4954      double value = 0.0;
4955      int nel2=nel;
4956      for (; nel2;nel2--) {
4957        int iRow = *row++;
4958        value += pi[iRow]*(*element++);
4959      }
4960      if (fabs(value)>zeroTolerance) {
4961        array[numberNonZero]=value;
4962        index[numberNonZero++]=*column;
4963      }
4964      column++;
4965    }
4966    numberPrice = numberPrice>>1;
4967    switch ((nel%2)) {
4968      // 2 k +0
4969    case 0:
4970      for (;numberPrice;numberPrice--) {
4971        double value0 = 0.0;
4972        double value1 = 0.0;
4973        int nel2=nel;
4974        for (; nel2;nel2--) {
4975          int iRow0 = *row;
4976          int iRow1 = *(row+nel);
4977          row++;
4978          double element0 = *element;
4979          double element1 = *(element+nel);
4980          element++;
4981          value0 += pi[iRow0]*element0;
4982          value1 += pi[iRow1]*element1;
4983        }
4984        row += nel;
4985        element += nel;
4986        if (fabs(value0)>zeroTolerance) {
4987          array[numberNonZero]=value0;
4988          index[numberNonZero++]=*column;
4989        }
4990        column++;
4991        if (fabs(value1)>zeroTolerance) {
4992          array[numberNonZero]=value1;
4993          index[numberNonZero++]=*column;
4994        }
4995        column++;
4996      } 
4997      break;
4998      // 2 k +1
4999    case 1:
5000      for (;numberPrice;numberPrice--) {
5001        double value0;
5002        double value1;
5003        int nel2=nel-1;
5004        {
5005          int iRow0 = row[0];
5006          int iRow1 = row[nel];
5007          double pi0 = pi[iRow0];
5008          double pi1 = pi[iRow1];
5009          value0 = pi0*element[0];
5010          value1 = pi1*element[nel];
5011          row++;
5012          element++;
5013        }
5014        for (; nel2;nel2--) {
5015          int iRow0 = *row;
5016          int iRow1 = *(row+nel);
5017          row++;
5018          double element0 = *element;
5019          double element1 = *(element+nel);
5020          element++;
5021          value0 += pi[iRow0]*element0;
5022          value1 += pi[iRow1]*element1;
5023        }
5024        row += nel;
5025        element += nel;
5026        if (fabs(value0)>zeroTolerance) {
5027          array[numberNonZero]=value0;
5028          index[numberNonZero++]=*column;
5029        }
5030        column++;
5031        if (fabs(value1)>zeroTolerance) {
5032          array[numberNonZero]=value1;
5033          index[numberNonZero++]=*column;
5034        }
5035        column++;
5036      }
5037      break;
5038    }
5039#else
5040    for (;numberPrice;numberPrice--) {
5041      double value = 0.0;
5042      int nel2=nel;
5043      for (; nel2;nel2--) {
5044        int iRow = *row++;
5045        value += pi[iRow]*(*element++);
5046      }
5047      if (fabs(value)>zeroTolerance) {
5048        array[numberNonZero]=value;
5049        index[numberNonZero++]=*column;
5050      }
5051      column++;
5052    }
5053#endif
5054  }
5055  output->setNumElements(numberNonZero);
5056}
5057// Updates two arrays for steepest
5058void 
5059ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
5060                                  const double * pi, CoinIndexedVector * output,
5061                                  const double * piWeight,
5062                                  double referenceIn, double devex,
5063                                  // Array for exact devex to say what is in reference framework
5064                                  unsigned int * reference,
5065                                  double * weights, double scaleFactor)
5066{
5067  int numberNonZero=0;
5068  int * index = output->getIndices();
5069  double * array = output->denseVector();
5070  double zeroTolerance = model->factorization()->zeroTolerance();
5071  double value = 0.0;
5072  bool killDjs = (scaleFactor==0.0);
5073  if (!scaleFactor)
5074    scaleFactor=1.0;
5075  int numberOdd = block_->startIndices_;
5076  int iColumn;
5077  CoinBigIndex end = start_[0];
5078  for (iColumn=0;iColumn<numberOdd;iColumn++) {
5079    CoinBigIndex start = end;
5080    CoinBigIndex j;
5081    int jColumn = column_[iColumn];
5082    end = start_[iColumn+1];
5083    value = 0.0;
5084    if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
5085      for (j=start; j<end;j++) {
5086        int iRow = row_[j];
5087        value -= pi[iRow]*element_[j];
5088      }
5089      if (fabs(value)>zeroTolerance) {
5090        // and do other array
5091        double modification = 0.0;
5092        for (j=start; j<end;j++) {
5093          int iRow = row_[j];
5094          modification += piWeight[iRow]*element_[j];
5095        }
5096        double thisWeight = weights[jColumn];
5097        double pivot = value*scaleFactor;
5098        double pivotSquared = pivot * pivot;
5099        thisWeight += pivotSquared * devex + pivot * modification;
5100        if (thisWeight<DEVEX_TRY_NORM) {
5101          if (referenceIn<0.0) {
5102            // steepest
5103            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
5104          } else {
5105            // exact
5106            thisWeight = referenceIn*pivotSquared;
5107            if (reference(jColumn))
5108              thisWeight += 1.0;
5109            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
5110          }
5111        }
5112        weights[jColumn] = thisWeight;
5113        if (!killDjs) {
5114          array[numberNonZero]=value;
5115          index[numberNonZero++]=jColumn;
5116        }
5117      }
5118    }
5119  }
5120  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
5121    // B) Can sort so just do nonbasic (and nonfixed)
5122    // C) Can do two at a time (if so put odd one into start_)
5123    // D) can use switch
5124    blockStruct * block = block_ + iBlock;
5125    //int numberPrice = block->numberInBlock_;
5126    int numberPrice = block->numberPrice_; 
5127    int nel = block->numberElements_;
5128    int * row = row_ + block->startElements_;
5129    double * element = element_ + block->startElements_;
5130    int * column = column_ + block->startIndices_;
5131    for (;numberPrice;numberPrice--) {
5132      double value = 0.0;
5133      int nel2=nel;
5134      for (; nel2;nel2--) {
5135        int iRow = *row++;
5136        value -= pi[iRow]*(*element++);
5137      }
5138      if (fabs(value)>zeroTolerance) {
5139        int jColumn = *column;
5140        // back to beginning
5141        row -= nel;
5142        element -= nel;
5143        // and do other array
5144        double modification = 0.0;
5145        nel2=nel;
5146        for (; nel2;nel2--) {
5147          int iRow = *row++;
5148          modification += piWeight[iRow]*(*element++);
5149        }
5150        double thisWeight = weights[jColumn];
5151        double pivot = value*scaleFactor;
5152        double pivotSquared = pivot * pivot;
5153        thisWeight += pivotSquared * devex + pivot * modification;
5154        if (thisWeight<DEVEX_TRY_NORM) {
5155          if (referenceIn<0.0) {
5156            // steepest
5157            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
5158          } else {
5159            // exact
5160            thisWeight = referenceIn*pivotSquared;
5161            if (reference(jColumn))
5162              thisWeight += 1.0;
5163            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
5164          }
5165        }
5166        weights[jColumn] = thisWeight;
5167        if (!killDjs) {
5168          array[numberNonZero]=value;
5169          index[numberNonZero++]=jColumn;
5170        }
5171      }
5172      column++;
5173    }
5174  }
5175  output->setNumElements(numberNonZero);
5176  output->setPackedMode(true);
5177}
5178#ifdef CLP_ALL_ONE_FILE
5179#undef reference
5180#endif
Note: See TracBrowser for help on using the repository browser.