source: trunk/ClpPackedMatrix.cpp @ 468

Last change on this file since 468 was 468, checked in by forrest, 16 years ago

to make sbb faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 67.9 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
11#include "ClpSimplex.hpp"
12#include "ClpFactorization.hpp"
13#include "ClpQuadraticObjective.hpp"
14// at end to get min/max!
15#include "ClpPackedMatrix.hpp"
16#include "ClpMessage.hpp"
17
18//#############################################################################
19// Constructors / Destructor / Assignment
20//#############################################################################
21
22//-------------------------------------------------------------------
23// Default Constructor
24//-------------------------------------------------------------------
25ClpPackedMatrix::ClpPackedMatrix () 
26  : ClpMatrixBase(),
27    matrix_(NULL),
28    numberActiveColumns_(0),
29    zeroElements_(false),
30    hasGaps_(true)
31{
32  setType(1);
33}
34
35//-------------------------------------------------------------------
36// Copy constructor
37//-------------------------------------------------------------------
38ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs) 
39: ClpMatrixBase(rhs)
40{ 
41  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
42  numberActiveColumns_ = rhs.numberActiveColumns_;
43  zeroElements_ = rhs.zeroElements_;
44  hasGaps_ = rhs.hasGaps_;
45  int numberRows = getNumRows();
46  if (rhs.rhsOffset_&&numberRows) {
47    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
48  } else {
49    rhsOffset_=NULL;
50  }
51 
52}
53
54//-------------------------------------------------------------------
55// assign matrix (for space reasons)
56//-------------------------------------------------------------------
57ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs) 
58: ClpMatrixBase()
59{ 
60  matrix_ = rhs;
61  zeroElements_ = false;
62  hasGaps_ = true;
63  numberActiveColumns_ = matrix_->getNumCols();
64  setType(1);
65 
66}
67
68ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs) 
69: ClpMatrixBase()
70{ 
71  matrix_ = new CoinPackedMatrix(rhs);
72  numberActiveColumns_ = matrix_->getNumCols();
73  zeroElements_ = false;
74  hasGaps_ = true;
75  setType(1);
76 
77}
78
79//-------------------------------------------------------------------
80// Destructor
81//-------------------------------------------------------------------
82ClpPackedMatrix::~ClpPackedMatrix ()
83{
84  delete matrix_;
85}
86
87//----------------------------------------------------------------
88// Assignment operator
89//-------------------------------------------------------------------
90ClpPackedMatrix &
91ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
92{
93  if (this != &rhs) {
94    ClpMatrixBase::operator=(rhs);
95    delete matrix_;
96    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
97    numberActiveColumns_ = rhs.numberActiveColumns_;
98    zeroElements_ = rhs.zeroElements_;
99    hasGaps_ = rhs.hasGaps_;
100  }
101  return *this;
102}
103//-------------------------------------------------------------------
104// Clone
105//-------------------------------------------------------------------
106ClpMatrixBase * ClpPackedMatrix::clone() const
107{
108  return new ClpPackedMatrix(*this);
109}
110/* Subset clone (without gaps).  Duplicates are allowed
111   and order is as given */
112ClpMatrixBase * 
113ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
114                              int numberColumns, 
115                              const int * whichColumns) const 
116{
117  return new ClpPackedMatrix(*this, numberRows, whichRows,
118                                   numberColumns, whichColumns);
119}
120/* Subset constructor (without gaps).  Duplicates are allowed
121   and order is as given */
122ClpPackedMatrix::ClpPackedMatrix (
123                       const ClpPackedMatrix & rhs,
124                       int numberRows, const int * whichRows,
125                       int numberColumns, const int * whichColumns)
126: ClpMatrixBase(rhs)
127{
128  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
129                                 numberColumns,whichColumns);
130  numberActiveColumns_ = matrix_->getNumCols();
131  zeroElements_ = rhs.zeroElements_;
132  hasGaps_ = rhs.hasGaps_;
133}
134ClpPackedMatrix::ClpPackedMatrix (
135                       const CoinPackedMatrix & rhs,
136                       int numberRows, const int * whichRows,
137                       int numberColumns, const int * whichColumns)
138: ClpMatrixBase()
139{
140  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
141                                 numberColumns,whichColumns);
142  numberActiveColumns_ = matrix_->getNumCols();
143  zeroElements_ = false;
144  hasGaps_=true;
145  setType(1);
146}
147
148/* Returns a new matrix in reverse order without gaps */
149ClpMatrixBase * 
150ClpPackedMatrix::reverseOrderedCopy() const
151{
152  ClpPackedMatrix * copy = new ClpPackedMatrix();
153  copy->matrix_= new CoinPackedMatrix();
154  copy->matrix_->setExtraGap(0.0);
155  copy->matrix_->setExtraMajor(0.0);
156  copy->matrix_->reverseOrderedCopyOf(*matrix_);
157  //copy->matrix_->removeGaps();
158  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
159  copy->hasGaps_=false;
160  return copy;
161}
162//unscaled versions
163void 
164ClpPackedMatrix::times(double scalar,
165                   const double * x, double * y) const
166{
167  int iRow,iColumn;
168  // get matrix data pointers
169  const int * row = matrix_->getIndices();
170  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
171  const int * columnLength = matrix_->getVectorLengths(); 
172  const double * elementByColumn = matrix_->getElements();
173  //memset(y,0,matrix_->getNumRows()*sizeof(double));
174  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
175    CoinBigIndex j;
176    double value = scalar*x[iColumn];
177    if (value) {
178      for (j=columnStart[iColumn];
179           j<columnStart[iColumn]+columnLength[iColumn];j++) {
180        iRow=row[j];
181        y[iRow] += value*elementByColumn[j];
182      }
183    }
184  }
185}
186void 
187ClpPackedMatrix::transposeTimes(double scalar,
188                                const double * x, double * y) const
189{
190  int iColumn;
191  // get matrix data pointers
192  const int * row = matrix_->getIndices();
193  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
194  const int * columnLength = matrix_->getVectorLengths(); 
195  const double * elementByColumn = matrix_->getElements();
196  //memset(y,0,numberActiveColumns_*sizeof(double));
197  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
198    CoinBigIndex j;
199    double value=0.0;
200    for (j=columnStart[iColumn];
201         j<columnStart[iColumn]+columnLength[iColumn];j++) {
202      int jRow=row[j];
203      value += x[jRow]*elementByColumn[j];
204    }
205    y[iColumn] += value*scalar;
206  }
207}
208void 
209ClpPackedMatrix::times(double scalar,
210                       const double * x, double * y,
211                       const double * rowScale, 
212                       const double * columnScale) const
213{
214  if (rowScale) {
215    int iRow,iColumn;
216    // get matrix data pointers
217    const int * row = matrix_->getIndices();
218    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
219    const int * columnLength = matrix_->getVectorLengths(); 
220    const double * elementByColumn = matrix_->getElements();
221    //memset(y,0,matrix_->getNumRows()*sizeof(double));
222    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
223      CoinBigIndex j;
224      double value = x[iColumn];
225      if (value) {
226        // scaled
227        value *= scalar*columnScale[iColumn];
228        for (j=columnStart[iColumn];
229             j<columnStart[iColumn]+columnLength[iColumn];j++) {
230          iRow=row[j];
231          y[iRow] += value*elementByColumn[j]*rowScale[iRow];
232        }
233      }
234    }
235  } else {
236    times(scalar,x,y);
237  }
238}
239void 
240ClpPackedMatrix::transposeTimes( double scalar,
241                                 const double * x, double * y,
242                                 const double * rowScale, 
243                                 const double * columnScale) const
244{
245  if (rowScale) {
246    int iColumn;
247    // get matrix data pointers
248    const int * row = matrix_->getIndices();
249    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
250    const int * columnLength = matrix_->getVectorLengths(); 
251    const double * elementByColumn = matrix_->getElements();
252    //memset(y,0,numberActiveColumns_*sizeof(double));
253    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
254      CoinBigIndex j;
255      double value=0.0;
256      // scaled
257      for (j=columnStart[iColumn];
258           j<columnStart[iColumn]+columnLength[iColumn];j++) {
259        int jRow=row[j];
260      value += x[jRow]*elementByColumn[j]*rowScale[jRow];
261      }
262      y[iColumn] += value*scalar*columnScale[iColumn];
263    }
264  } else {
265    transposeTimes(scalar,x,y);
266  }
267}
268/* Return <code>x * A + y</code> in <code>z</code>.
269        Squashes small elements and knows about ClpSimplex */
270void 
271ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
272                              const CoinIndexedVector * rowArray,
273                              CoinIndexedVector * y,
274                              CoinIndexedVector * columnArray) const
275{
276  columnArray->clear();
277  double * pi = rowArray->denseVector();
278  int numberNonZero=0;
279  int * index = columnArray->getIndices();
280  double * array = columnArray->denseVector();
281  int numberInRowArray = rowArray->getNumElements();
282  // maybe I need one in OsiSimplex
283  double zeroTolerance = model->factorization()->zeroTolerance();
284  int numberRows = model->numberRows();
285#ifndef NO_RTTI
286  ClpPackedMatrix* rowCopy =
287    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
288#else
289  ClpPackedMatrix* rowCopy =
290    static_cast< ClpPackedMatrix*>(model->rowCopy());
291#endif
292  bool packed = rowArray->packedMode();
293  double factor = 0.3;
294  // We may not want to do by row if there may be cache problems
295  // It would be nice to find L2 cache size - for moment 512K
296  // Be slightly optimistic
297  if (numberActiveColumns_*sizeof(double)>1000000) {
298    if (numberRows*10<numberActiveColumns_)
299      factor=0.1;
300    else if (numberRows*4<numberActiveColumns_)
301      factor=0.15;
302    else if (numberRows*2<numberActiveColumns_)
303      factor=0.2;
304    //if (model->numberIterations()%50==0)
305    //printf("%d nonzero\n",numberInRowArray);
306  }
307  assert (!y->getNumElements());
308  if (numberInRowArray>factor*numberRows||!rowCopy) {
309    // do by column
310    // If no gaps - can do a bit faster
311    if (!hasGaps_&&true) {
312      transposeTimesByColumn( model,  scalar,
313                              rowArray, y, columnArray);
314      return;
315    }
316    int iColumn;
317    // get matrix data pointers
318    const int * row = matrix_->getIndices();
319    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
320    const int * columnLength = matrix_->getVectorLengths(); 
321    const double * elementByColumn = matrix_->getElements();
322    const double * rowScale = model->rowScale();
323    if (packed) {
324      // need to expand pi into y
325      assert(y->capacity()>=numberRows);
326      double * piOld = pi;
327      pi = y->denseVector();
328      const int * whichRow = rowArray->getIndices();
329      int i;
330      if (!rowScale) {
331        // modify pi so can collapse to one loop
332        for (i=0;i<numberInRowArray;i++) {
333          int iRow = whichRow[i];
334          pi[iRow]=scalar*piOld[i];
335        }
336        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
337          double value = 0.0;
338          CoinBigIndex j;
339          for (j=columnStart[iColumn];
340               j<columnStart[iColumn]+columnLength[iColumn];j++) {
341            int iRow = row[j];
342            value += pi[iRow]*elementByColumn[j];
343          }
344          if (fabs(value)>zeroTolerance) {
345            array[numberNonZero]=value;
346            index[numberNonZero++]=iColumn;
347          }
348        }
349      } else {
350        // scaled
351        // modify pi so can collapse to one loop
352        for (i=0;i<numberInRowArray;i++) {
353          int iRow = whichRow[i];
354          pi[iRow]=scalar*piOld[i]*rowScale[iRow];
355        }
356        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
357          double value = 0.0;
358          CoinBigIndex j;
359          const double * columnScale = model->columnScale();
360          for (j=columnStart[iColumn];
361               j<columnStart[iColumn]+columnLength[iColumn];j++) {
362            int iRow = row[j];
363            value += pi[iRow]*elementByColumn[j];
364          }
365          value *= columnScale[iColumn];
366          if (fabs(value)>zeroTolerance) {
367            array[numberNonZero]=value;
368            index[numberNonZero++]=iColumn;
369          }
370        }
371      }
372      // zero out
373      for (i=0;i<numberInRowArray;i++) {
374        int iRow = whichRow[i];
375        pi[iRow]=0.0;
376      }
377    } else {
378      if (!rowScale) {
379        if (scalar==-1.0) {
380          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
381            double value = 0.0;
382            CoinBigIndex j;
383            for (j=columnStart[iColumn];
384                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
385              int iRow = row[j];
386              value += pi[iRow]*elementByColumn[j];
387            }
388            if (fabs(value)>zeroTolerance) {
389              index[numberNonZero++]=iColumn;
390              array[iColumn]=-value;
391            }
392          }
393        } else if (scalar==1.0) {
394          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
395            double value = 0.0;
396            CoinBigIndex j;
397            for (j=columnStart[iColumn];
398                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
399              int iRow = row[j];
400              value += pi[iRow]*elementByColumn[j];
401            }
402            if (fabs(value)>zeroTolerance) {
403              index[numberNonZero++]=iColumn;
404              array[iColumn]=value;
405            }
406          }
407        } else {
408          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
409            double value = 0.0;
410            CoinBigIndex j;
411            for (j=columnStart[iColumn];
412                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
413              int iRow = row[j];
414              value += pi[iRow]*elementByColumn[j];
415            }
416            value *= scalar;
417            if (fabs(value)>zeroTolerance) {
418              index[numberNonZero++]=iColumn;
419              array[iColumn]=value;
420            }
421          }
422        }
423      } else {
424        // scaled
425        if (scalar==-1.0) {
426          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
427            double value = 0.0;
428            CoinBigIndex j;
429            const double * columnScale = model->columnScale();
430            for (j=columnStart[iColumn];
431                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
432              int iRow = row[j];
433              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
434            }
435            value *= columnScale[iColumn];
436            if (fabs(value)>zeroTolerance) {
437              index[numberNonZero++]=iColumn;
438              array[iColumn]=-value;
439            }
440          }
441        } else if (scalar==1.0) {
442          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
443            double value = 0.0;
444            CoinBigIndex j;
445            const double * columnScale = model->columnScale();
446            for (j=columnStart[iColumn];
447                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
448              int iRow = row[j];
449              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
450            }
451            value *= columnScale[iColumn];
452            if (fabs(value)>zeroTolerance) {
453              index[numberNonZero++]=iColumn;
454              array[iColumn]=value;
455            }
456          }
457        } else {
458          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
459            double value = 0.0;
460            CoinBigIndex j;
461            const double * columnScale = model->columnScale();
462            for (j=columnStart[iColumn];
463                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
464              int iRow = row[j];
465              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
466            }
467            value *= scalar*columnScale[iColumn];
468            if (fabs(value)>zeroTolerance) {
469              index[numberNonZero++]=iColumn;
470              array[iColumn]=value;
471            }
472          }
473        }
474      }
475    }
476    columnArray->setNumElements(numberNonZero);
477    y->setNumElements(0);
478  } else {
479    // do by row
480    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
481  }
482  if (packed)
483    columnArray->setPackedMode(true);
484  if (0) {
485    columnArray->checkClean();
486    int numberNonZero=columnArray->getNumElements();;
487    int * index = columnArray->getIndices();
488    double * array = columnArray->denseVector();
489    int i;
490    for (i=0;i<numberNonZero;i++) {
491      int j=index[i];
492      double value;
493      if (packed)
494        value=array[i];
495      else
496        value=array[j];
497      printf("Ti %d %d %g\n",i,j,value);
498    }
499  }
500}
501/* Return <code>x * scalar * A + y</code> in <code>z</code>.
502   Note - If x packed mode - then z packed mode
503   This does by column and knows no gaps
504   Squashes small elements and knows about ClpSimplex */
505void 
506ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
507                                        const CoinIndexedVector * rowArray,
508                                        CoinIndexedVector * y,
509                                        CoinIndexedVector * columnArray) const
510{
511  double * pi = rowArray->denseVector();
512  int numberNonZero=0;
513  int * index = columnArray->getIndices();
514  double * array = columnArray->denseVector();
515  int numberInRowArray = rowArray->getNumElements();
516  // maybe I need one in OsiSimplex
517  double zeroTolerance = model->factorization()->zeroTolerance();
518  bool packed = rowArray->packedMode();
519  // do by column
520  int iColumn;
521  // get matrix data pointers
522  const int * row = matrix_->getIndices();
523  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
524  const double * elementByColumn = matrix_->getElements();
525  const double * rowScale = model->rowScale();
526  assert (!y->getNumElements());
527  assert (numberActiveColumns_>0);
528  if (packed) {
529    // need to expand pi into y
530    assert(y->capacity()>=model->numberRows());
531    double * piOld = pi;
532    pi = y->denseVector();
533    const int * whichRow = rowArray->getIndices();
534    int i;
535    if (!rowScale) {
536      // modify pi so can collapse to one loop
537      for (i=0;i<numberInRowArray;i++) {
538        int iRow = whichRow[i];
539        pi[iRow]=scalar*piOld[i];
540      }
541      double value = 0.0;
542      CoinBigIndex j;
543      CoinBigIndex end = columnStart[1];
544      for (j=columnStart[0]; j<end;j++) {
545        int iRow = row[j];
546        value += pi[iRow]*elementByColumn[j];
547      }
548      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
549        CoinBigIndex start = end;
550        end = columnStart[iColumn+2];
551        if (fabs(value)>zeroTolerance) {
552          array[numberNonZero]=value;
553          index[numberNonZero++]=iColumn;
554        }
555        value = 0.0;
556        for (j=start; j<end;j++) {
557          int iRow = row[j];
558          value += pi[iRow]*elementByColumn[j];
559        }
560      }
561      if (fabs(value)>zeroTolerance) {
562        array[numberNonZero]=value;
563        index[numberNonZero++]=iColumn;
564      }
565    } else {
566      // scaled
567      // modify pi so can collapse to one loop
568      for (i=0;i<numberInRowArray;i++) {
569        int iRow = whichRow[i];
570        pi[iRow]=scalar*piOld[i]*rowScale[iRow];
571      }
572      const double * columnScale = model->columnScale();
573      double value = 0.0;
574      double scale=columnScale[0];
575      CoinBigIndex j;
576      CoinBigIndex end = columnStart[1];
577      for (j=columnStart[0]; j<end;j++) {
578        int iRow = row[j];
579        value += pi[iRow]*elementByColumn[j];
580      }
581      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
582        value *= scale;
583        CoinBigIndex start = end;
584        scale = columnScale[iColumn+1];
585        end = columnStart[iColumn+2];
586        if (fabs(value)>zeroTolerance) {
587          array[numberNonZero]=value;
588          index[numberNonZero++]=iColumn;
589        }
590        value = 0.0;
591        for (j=start; j<end;j++) {
592          int iRow = row[j];
593          value += pi[iRow]*elementByColumn[j];
594        }
595      }
596      value *= scale;
597      if (fabs(value)>zeroTolerance) {
598        array[numberNonZero]=value;
599        index[numberNonZero++]=iColumn;
600      }
601    }
602    // zero out
603    for (i=0;i<numberInRowArray;i++) {
604      int iRow = whichRow[i];
605      pi[iRow]=0.0;
606    }
607  } else {
608    if (!rowScale) {
609      if (scalar==-1.0) {
610        double value = 0.0;
611        CoinBigIndex j;
612        CoinBigIndex end = columnStart[1];
613        for (j=columnStart[0]; j<end;j++) {
614          int iRow = row[j];
615          value += pi[iRow]*elementByColumn[j];
616        }
617        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
618          CoinBigIndex start = end;
619          end = columnStart[iColumn+2];
620          if (fabs(value)>zeroTolerance) {
621            array[iColumn]=-value;
622            index[numberNonZero++]=iColumn;
623          }
624          value = 0.0;
625          for (j=start; j<end;j++) {
626            int iRow = row[j];
627            value += pi[iRow]*elementByColumn[j];
628          }
629        }
630        if (fabs(value)>zeroTolerance) {
631          array[iColumn]=-value;
632          index[numberNonZero++]=iColumn;
633        }
634      } else if (scalar==1.0) {
635        double value = 0.0;
636        CoinBigIndex j;
637        CoinBigIndex end = columnStart[1];
638        for (j=columnStart[0]; j<end;j++) {
639          int iRow = row[j];
640          value += pi[iRow]*elementByColumn[j];
641        }
642        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
643          CoinBigIndex start = end;
644          end = columnStart[iColumn+2];
645          if (fabs(value)>zeroTolerance) {
646            array[iColumn]=value;
647            index[numberNonZero++]=iColumn;
648          }
649          value = 0.0;
650          for (j=start; j<end;j++) {
651            int iRow = row[j];
652            value += pi[iRow]*elementByColumn[j];
653          }
654        }
655        if (fabs(value)>zeroTolerance) {
656          array[iColumn]=value;
657          index[numberNonZero++]=iColumn;
658        }
659      } else {
660        double value = 0.0;
661        CoinBigIndex j;
662        CoinBigIndex end = columnStart[1];
663        for (j=columnStart[0]; j<end;j++) {
664          int iRow = row[j];
665          value += pi[iRow]*elementByColumn[j];
666        }
667        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
668          value *= scalar;
669          CoinBigIndex start = end;
670          end = columnStart[iColumn+2];
671          if (fabs(value)>zeroTolerance) {
672            array[iColumn]=value;
673            index[numberNonZero++]=iColumn;
674          }
675          value = 0.0;
676          for (j=start; j<end;j++) {
677            int iRow = row[j];
678            value += pi[iRow]*elementByColumn[j];
679          }
680        }
681        value *= scalar;
682        if (fabs(value)>zeroTolerance) {
683          array[iColumn]=value;
684          index[numberNonZero++]=iColumn;
685        }
686      }
687    } else {
688      // unscaled
689      if (scalar==-1.0) {
690        const double * columnScale = model->columnScale();
691        double value = 0.0;
692        double scale=columnScale[0];
693        CoinBigIndex j;
694        CoinBigIndex end = columnStart[1];
695        for (j=columnStart[0]; j<end;j++) {
696          int iRow = row[j];
697          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
698        }
699        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
700          value *= scale;
701          CoinBigIndex start = end;
702          end = columnStart[iColumn+2];
703          scale = columnScale[iColumn+1];
704          if (fabs(value)>zeroTolerance) {
705            array[iColumn]=-value;
706            index[numberNonZero++]=iColumn;
707          }
708          value = 0.0;
709          for (j=start; j<end;j++) {
710            int iRow = row[j];
711            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
712          }
713        }
714        value *= scale;
715        if (fabs(value)>zeroTolerance) {
716          array[iColumn]=-value;
717          index[numberNonZero++]=iColumn;
718        }
719      } else if (scalar==1.0) {
720        const double * columnScale = model->columnScale();
721        double value = 0.0;
722        double scale=columnScale[0];
723        CoinBigIndex j;
724        CoinBigIndex end = columnStart[1];
725        for (j=columnStart[0]; j<end;j++) {
726          int iRow = row[j];
727          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
728        }
729        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
730          value *= scale;
731          CoinBigIndex start = end;
732          end = columnStart[iColumn+2];
733          scale = columnScale[iColumn+1];
734          if (fabs(value)>zeroTolerance) {
735            array[iColumn]=value;
736            index[numberNonZero++]=iColumn;
737          }
738          value = 0.0;
739          for (j=start; j<end;j++) {
740            int iRow = row[j];
741            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
742          }
743        }
744        value *= scale;
745        if (fabs(value)>zeroTolerance) {
746          array[iColumn]=value;
747          index[numberNonZero++]=iColumn;
748        }
749      } else {
750        const double * columnScale = model->columnScale();
751        double value = 0.0;
752        double scale=columnScale[0]*scalar;
753        CoinBigIndex j;
754        CoinBigIndex end = columnStart[1];
755        for (j=columnStart[0]; j<end;j++) {
756          int iRow = row[j];
757          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
758        }
759        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
760          value *= scale;
761          CoinBigIndex start = end;
762          end = columnStart[iColumn+2];
763          scale = columnScale[iColumn+1]*scalar;
764          if (fabs(value)>zeroTolerance) {
765            array[iColumn]=value;
766            index[numberNonZero++]=iColumn;
767          }
768          value = 0.0;
769          for (j=start; j<end;j++) {
770            int iRow = row[j];
771            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
772          }
773        }
774        value *= scale;
775        if (fabs(value)>zeroTolerance) {
776          array[iColumn]=value;
777          index[numberNonZero++]=iColumn;
778        }
779      }
780    }
781  }
782  columnArray->setNumElements(numberNonZero);
783  y->setNumElements(0);
784  if (packed)
785    columnArray->setPackedMode(true);
786}
787/* Return <code>x * scalar * A in <code>z</code>.
788   Note - this version when x packed mode and so will be z
789   This does by column and knows no gaps and knows y empty
790   Squashes small elements and knows about ClpSimplex */
791void 
792ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
793                                        const CoinIndexedVector * rowArray,
794                                        CoinIndexedVector * columnArray) const
795{
796}
797/* Return <code>x * A + y</code> in <code>z</code>.
798        Squashes small elements and knows about ClpSimplex */
799void 
800ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
801                              const CoinIndexedVector * rowArray,
802                              CoinIndexedVector * y,
803                              CoinIndexedVector * columnArray) const
804{
805  columnArray->clear();
806  double * pi = rowArray->denseVector();
807  int numberNonZero=0;
808  int * index = columnArray->getIndices();
809  double * array = columnArray->denseVector();
810  int numberInRowArray = rowArray->getNumElements();
811  // maybe I need one in OsiSimplex
812  double zeroTolerance = model->factorization()->zeroTolerance();
813  const int * column = getIndices();
814  const CoinBigIndex * rowStart = getVectorStarts();
815  const double * element = getElements();
816  const int * whichRow = rowArray->getIndices();
817  bool packed = rowArray->packedMode();
818  if (numberInRowArray>2) {
819    // do by rows
820    // ** Row copy is already scaled
821    int iRow;
822    int i;
823    int numberOriginal = 0;
824    if (packed) {
825      numberNonZero=0;
826      // and set up mark as char array
827      char * marked = (char *) (index+columnArray->capacity());
828      double * array2 = y->denseVector();
829#ifdef CLP_DEBUG
830      for (i=0;i<numberActiveColumns_;i++) {
831        assert(!marked[i]);
832        assert(!array2[i]);
833      }
834#endif     
835      for (i=0;i<numberInRowArray;i++) {
836        iRow = whichRow[i]; 
837        double value = pi[i]*scalar;
838        CoinBigIndex j;
839        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
840          int iColumn = column[j];
841          if (!marked[iColumn]) {
842            marked[iColumn]=1;
843            index[numberNonZero++]=iColumn;
844          }
845          array2[iColumn] += value*element[j];
846        }
847      }
848      // get rid of tiny values and zero out marked
849      numberOriginal=numberNonZero;
850      numberNonZero=0;
851      for (i=0;i<numberOriginal;i++) {
852        int iColumn = index[i];
853        if (marked[iColumn]) {
854          double value = array2[iColumn];
855          array2[iColumn]=0.0;
856          marked[iColumn]=0;
857          if (fabs(value)>zeroTolerance) {
858            array[numberNonZero]=value;
859            index[numberNonZero++]=iColumn;
860          }
861        }
862      }
863    } else {
864      double * markVector = y->denseVector();
865      numberNonZero=0;
866      // and set up mark as char array
867      char * marked = (char *) markVector;
868      for (i=0;i<numberOriginal;i++) {
869        int iColumn = index[i];
870        marked[iColumn]=0;
871      }
872     
873      for (i=0;i<numberInRowArray;i++) {
874        iRow = whichRow[i]; 
875        double value = pi[iRow]*scalar;
876        CoinBigIndex j;
877        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
878          int iColumn = column[j];
879          if (!marked[iColumn]) {
880            marked[iColumn]=1;
881            index[numberNonZero++]=iColumn;
882          }
883          array[iColumn] += value*element[j];
884        }
885      }
886      // get rid of tiny values and zero out marked
887      numberOriginal=numberNonZero;
888      numberNonZero=0;
889      for (i=0;i<numberOriginal;i++) {
890        int iColumn = index[i];
891        marked[iColumn]=0;
892        if (fabs(array[iColumn])>zeroTolerance) {
893          index[numberNonZero++]=iColumn;
894        } else {
895          array[iColumn]=0.0;
896        }
897      }
898    }
899  } else if (numberInRowArray==2) {
900    // do by rows when two rows
901    int numberOriginal;
902    int i;
903    CoinBigIndex j;
904    numberNonZero=0;
905
906    double value;
907    if (packed) {
908      int iRow0 = whichRow[0]; 
909      int iRow1 = whichRow[1]; 
910      double pi0 = pi[0];
911      double pi1 = pi[1];
912      if (rowStart[iRow0+1]-rowStart[iRow0]>
913          rowStart[iRow1+1]-rowStart[iRow1]) {
914        // do one with fewer first
915        iRow0=iRow1;
916        iRow1=whichRow[0];
917        pi0=pi1;
918        pi1=pi[0];
919      }
920      // and set up mark as char array
921      char * marked = (char *) (index+columnArray->capacity());
922      int * lookup = y->getIndices();
923      value = pi0*scalar;
924      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
925        int iColumn = column[j];
926        double value2 = value*element[j];
927        array[numberNonZero] = value2;
928        marked[iColumn]=1;
929        lookup[iColumn]=numberNonZero;
930        index[numberNonZero++]=iColumn;
931      }
932      numberOriginal = numberNonZero;
933      value = pi1*scalar;
934      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
935        int iColumn = column[j];
936        double value2 = value*element[j];
937        // I am assuming no zeros in matrix
938        if (marked[iColumn]) {
939          int iLookup = lookup[iColumn];
940          array[iLookup] += value2;
941        } else {
942          if (fabs(value2)>zeroTolerance) {
943            array[numberNonZero] = value2;
944            index[numberNonZero++]=iColumn;
945          }
946        }
947      }
948      // get rid of tiny values and zero out marked
949      int nDelete=0;
950      for (i=0;i<numberOriginal;i++) {
951        int iColumn = index[i];
952        marked[iColumn]=0;
953        if (fabs(array[i])<=zeroTolerance) 
954          nDelete++;
955      }
956      if (nDelete) {
957        numberOriginal=numberNonZero;
958        numberNonZero=0;
959        for (i=0;i<numberOriginal;i++) {
960          int iColumn = index[i];
961          double value = array[i];
962          array[i]=0.0;
963          if (fabs(value)>zeroTolerance) {
964            array[numberNonZero]=value;
965            index[numberNonZero++]=iColumn;
966          }
967        }
968      }
969    } else {
970      int iRow = whichRow[0]; 
971      value = pi[iRow]*scalar;
972      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
973        int iColumn = column[j];
974        double value2 = value*element[j];
975        index[numberNonZero++]=iColumn;
976        array[iColumn] = value2;
977      }
978      iRow = whichRow[1]; 
979      value = pi[iRow]*scalar;
980      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
981        int iColumn = column[j];
982        double value2 = value*element[j];
983        // I am assuming no zeros in matrix
984        if (array[iColumn])
985          value2 += array[iColumn];
986        else
987          index[numberNonZero++]=iColumn;
988        array[iColumn] = value2;
989      }
990      // get rid of tiny values and zero out marked
991      numberOriginal=numberNonZero;
992      numberNonZero=0;
993      for (i=0;i<numberOriginal;i++) {
994        int iColumn = index[i];
995        if (fabs(array[iColumn])>zeroTolerance) {
996          index[numberNonZero++]=iColumn;
997        } else {
998          array[iColumn]=0.0;
999        }
1000      }
1001    }
1002  } else if (numberInRowArray==1) {
1003    // Just one row
1004    int iRow=rowArray->getIndices()[0];
1005    numberNonZero=0;
1006    CoinBigIndex j;
1007    if (packed) {
1008      double value = pi[0]*scalar;
1009      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1010        int iColumn = column[j];
1011        double value2 = value*element[j];
1012        if (fabs(value2)>zeroTolerance) {
1013          array[numberNonZero] = value2;
1014          index[numberNonZero++]=iColumn;
1015        }
1016      }
1017    } else {
1018      double value = pi[iRow]*scalar;
1019      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1020        int iColumn = column[j];
1021        double value2 = value*element[j];
1022        if (fabs(value2)>zeroTolerance) {
1023          index[numberNonZero++]=iColumn;
1024          array[iColumn] = value2;
1025        }
1026      }
1027    }
1028  }
1029  columnArray->setNumElements(numberNonZero);
1030  y->setNumElements(0);
1031}
1032/* Return <code>x *A in <code>z</code> but
1033   just for indices in y.
1034   Squashes small elements and knows about ClpSimplex */
1035void 
1036ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
1037                              const CoinIndexedVector * rowArray,
1038                              const CoinIndexedVector * y,
1039                              CoinIndexedVector * columnArray) const
1040{
1041  columnArray->clear();
1042  double * pi = rowArray->denseVector();
1043  double * array = columnArray->denseVector();
1044  int jColumn;
1045  // get matrix data pointers
1046  const int * row = matrix_->getIndices();
1047  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1048  const int * columnLength = matrix_->getVectorLengths(); 
1049  const double * elementByColumn = matrix_->getElements();
1050  const double * rowScale = model->rowScale();
1051  int numberToDo = y->getNumElements();
1052  const int * which = y->getIndices();
1053  assert (!rowArray->packedMode());
1054  columnArray->setPacked();
1055  if (!hasGaps_&&numberToDo>5) {
1056    // no gaps and a reasonable number
1057    if (!rowScale) {
1058      int iColumn = which[0];
1059      double value = 0.0;
1060      CoinBigIndex j;
1061      for (j=columnStart[iColumn];
1062           j<columnStart[iColumn+1];j++) {
1063        int iRow = row[j];
1064        value += pi[iRow]*elementByColumn[j];
1065      }
1066      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
1067        int iColumn = which[jColumn+1];
1068        CoinBigIndex start = columnStart[iColumn];
1069        CoinBigIndex end = columnStart[iColumn+1];
1070        array[jColumn]=value;
1071        value = 0.0;
1072        for (j=start; j<end;j++) {
1073          int iRow = row[j];
1074          value += pi[iRow]*elementByColumn[j];
1075        }
1076      }
1077      array[jColumn]=value;
1078    } else {
1079      // scaled
1080      const double * columnScale = model->columnScale();
1081      int iColumn = which[0];
1082      double value = 0.0;
1083      double scale=columnScale[iColumn];
1084      CoinBigIndex j;
1085      for (j=columnStart[iColumn];
1086           j<columnStart[iColumn+1];j++) {
1087        int iRow = row[j];
1088        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1089      }
1090      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
1091        int iColumn = which[jColumn+1];
1092        value *= scale;
1093        scale = columnScale[iColumn];
1094        CoinBigIndex start = columnStart[iColumn];
1095        CoinBigIndex end = columnStart[iColumn+1];
1096        array[jColumn]=value;
1097        value = 0.0;
1098        for (j=start; j<end;j++) {
1099          int iRow = row[j];
1100          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1101        }
1102      }
1103      value *= scale;
1104      array[jColumn]=value;
1105    }
1106  } else {
1107    // gaps
1108    if (!rowScale) {
1109      for (jColumn=0;jColumn<numberToDo;jColumn++) {
1110        int iColumn = which[jColumn];
1111        double value = 0.0;
1112        CoinBigIndex j;
1113        for (j=columnStart[iColumn];
1114             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1115          int iRow = row[j];
1116          value += pi[iRow]*elementByColumn[j];
1117        }
1118        array[jColumn]=value;
1119      }
1120    } else {
1121      // scaled
1122      const double * columnScale = model->columnScale();
1123      for (jColumn=0;jColumn<numberToDo;jColumn++) {
1124        int iColumn = which[jColumn];
1125        double value = 0.0;
1126        CoinBigIndex j;
1127        for (j=columnStart[iColumn];
1128             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1129          int iRow = row[j];
1130          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1131        }
1132        value *= columnScale[iColumn];
1133        array[jColumn]=value;
1134      }
1135    }
1136  }
1137}
1138/// returns number of elements in column part of basis,
1139CoinBigIndex
1140ClpPackedMatrix::countBasis(ClpSimplex * model,
1141                           const int * whichColumn, 
1142                           int numberBasic,
1143                            int & numberColumnBasic)
1144{
1145  const int * columnLength = matrix_->getVectorLengths(); 
1146  int i;
1147  CoinBigIndex numberElements=0;
1148  // just count - can be over so ignore zero problem
1149  for (i=0;i<numberColumnBasic;i++) {
1150    int iColumn = whichColumn[i];
1151    numberElements += columnLength[iColumn];
1152  }
1153  return numberElements;
1154}
1155void
1156ClpPackedMatrix::fillBasis(ClpSimplex * model,
1157                         const int * whichColumn, 
1158                         int & numberColumnBasic,
1159                         int * indexRowU, int * start,
1160                         int * rowCount, int * columnCount,
1161                         double * elementU)
1162{
1163  const int * columnLength = matrix_->getVectorLengths(); 
1164  int i;
1165  CoinBigIndex numberElements=start[0];
1166  // fill
1167  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1168  const double * rowScale = model->rowScale();
1169  const int * row = matrix_->getIndices();
1170  const double * elementByColumn = matrix_->getElements();
1171  if (!zeroElements_) {
1172    if (!rowScale) {
1173      // no scaling
1174      for (i=0;i<numberColumnBasic;i++) {
1175        int iColumn = whichColumn[i];
1176        CoinBigIndex j;
1177        for (j=columnStart[iColumn];
1178             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1179          int iRow=row[j];
1180          indexRowU[numberElements]=iRow;
1181          rowCount[iRow]++;
1182          elementU[numberElements++]=elementByColumn[j];
1183        }
1184        start[i+1]=numberElements;
1185        columnCount[i]=columnLength[iColumn];
1186      }
1187    } else {
1188      // scaling
1189      const double * columnScale = model->columnScale();
1190      for (i=0;i<numberColumnBasic;i++) {
1191        int iColumn = whichColumn[i];
1192        CoinBigIndex j;
1193        double scale = columnScale[iColumn];
1194        for (j=columnStart[iColumn];
1195             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1196          int iRow = row[j];
1197          indexRowU[numberElements]=iRow;
1198          rowCount[iRow]++;
1199          elementU[numberElements++]=
1200            elementByColumn[j]*scale*rowScale[iRow];
1201        }
1202        start[i+1]=numberElements;
1203        columnCount[i]=columnLength[iColumn];
1204      }
1205    }
1206  } else {
1207    // there are zero elements so need to look more closely
1208    if (!rowScale) {
1209      // no scaling
1210      for (i=0;i<numberColumnBasic;i++) {
1211        int iColumn = whichColumn[i];
1212        CoinBigIndex j;
1213        for (j=columnStart[iColumn];
1214             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1215          double value = elementByColumn[j];
1216          if (value) {
1217            int iRow = row[j];
1218            indexRowU[numberElements]=iRow;
1219            rowCount[iRow]++;
1220            elementU[numberElements++]=value;
1221          }
1222        }
1223        start[i+1]=numberElements;
1224        columnCount[i]=numberElements-start[i];
1225      }
1226    } else {
1227      // scaling
1228      const double * columnScale = model->columnScale();
1229      for (i=0;i<numberColumnBasic;i++) {
1230        int iColumn = whichColumn[i];
1231        CoinBigIndex j;
1232        double scale = columnScale[iColumn];
1233        for (j=columnStart[iColumn];
1234             j<columnStart[iColumn]+columnLength[i];j++) {
1235          double value = elementByColumn[j];
1236          if (value) {
1237            int iRow = row[j];
1238            indexRowU[numberElements]=iRow;
1239            rowCount[iRow]++;
1240            elementU[numberElements++]=value*scale*rowScale[iRow];
1241          }
1242        }
1243        start[i+1]=numberElements;
1244        columnCount[i]=numberElements-start[i];
1245      }
1246    }
1247  }
1248}
1249// Creates scales for column copy (rowCopy in model may be modified)
1250int 
1251ClpPackedMatrix::scale(ClpModel * model) const 
1252{
1253  int numberRows = model->numberRows();
1254  int numberColumns = matrix_->getNumCols();
1255  // If empty - return as sanityCheck will trap
1256  if (!numberRows||!numberColumns)
1257    return 1;
1258  ClpMatrixBase * rowCopyBase=model->rowCopy();
1259  if (!rowCopyBase) {
1260    // temporary copy
1261    rowCopyBase = reverseOrderedCopy();
1262  }
1263#ifndef NO_RTTI
1264  ClpPackedMatrix* rowCopy =
1265    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
1266  // Make sure it is really a ClpPackedMatrix
1267  assert (rowCopy!=NULL);
1268#else
1269  ClpPackedMatrix* rowCopy =
1270    static_cast< ClpPackedMatrix*>(rowCopyBase);
1271#endif
1272
1273  const int * column = rowCopy->getIndices();
1274  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1275  const double * element = rowCopy->getElements();
1276  double * rowScale = new double [numberRows];
1277  double * columnScale = new double [numberColumns];
1278  // we are going to mark bits we are interested in
1279  char * usefulRow = new char [numberRows];
1280  char * usefulColumn = new char [numberColumns];
1281  double * rowLower = model->rowLower();
1282  double * rowUpper = model->rowUpper();
1283  double * columnLower = model->columnLower();
1284  double * columnUpper = model->columnUpper();
1285  int iColumn, iRow;
1286  // mark free rows
1287  for (iRow=0;iRow<numberRows;iRow++) {
1288    usefulRow[iRow]=0;
1289    if (rowUpper[iRow]<1.0e20||
1290        rowLower[iRow]>-1.0e20)
1291      usefulRow[iRow]=1;
1292  }
1293  // mark empty and fixed columns
1294  // also see if worth scaling
1295  assert (model->scalingFlag()<4); // dynamic not implemented
1296  double largest=0.0;
1297  double smallest=1.0e50;
1298  // get matrix data pointers
1299  const int * row = matrix_->getIndices();
1300  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1301  const int * columnLength = matrix_->getVectorLengths(); 
1302  const double * elementByColumn = matrix_->getElements();
1303  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1304    CoinBigIndex j;
1305    char useful=0;
1306    if (columnUpper[iColumn]>
1307        columnLower[iColumn]+1.0e-9) {
1308      for (j=columnStart[iColumn];
1309           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1310        iRow=row[j];
1311        if(elementByColumn[j]&&usefulRow[iRow]) {
1312          useful=1;
1313          largest = CoinMax(largest,fabs(elementByColumn[j]));
1314          smallest = CoinMin(smallest,fabs(elementByColumn[j]));
1315        }
1316      }
1317    }
1318    usefulColumn[iColumn]=useful;
1319  }
1320  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
1321    <<smallest<<largest
1322    <<CoinMessageEol;
1323  if (smallest>=0.5&&largest<=2.0) {
1324    // don't bother scaling
1325    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
1326      <<CoinMessageEol;
1327    delete [] rowScale;
1328    delete [] usefulRow;
1329    delete [] columnScale;
1330    delete [] usefulColumn;
1331    if (!model->rowCopy()) 
1332      delete rowCopyBase; // get rid of temporary row copy
1333    return 1;
1334  } else {
1335      // need to scale
1336    int scalingMethod = model->scalingFlag();
1337    if (scalingMethod==3) {
1338      // Choose between 1 and 2
1339      if (smallest<1.0e-5||smallest*largest<1.0)
1340        scalingMethod=1;
1341      else
1342        scalingMethod=2;
1343    }
1344    // and see if there any empty rows
1345    for (iRow=0;iRow<numberRows;iRow++) {
1346      if (usefulRow[iRow]) {
1347        CoinBigIndex j;
1348        int useful=0;
1349        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1350          int iColumn = column[j];
1351          if (usefulColumn[iColumn]) {
1352            useful=1;
1353            break;
1354          }
1355        }
1356        usefulRow[iRow]=useful;
1357      }
1358    }
1359    ClpFillN ( rowScale, numberRows,1.0);
1360    ClpFillN ( columnScale, numberColumns,1.0);
1361    double overallLargest=-1.0e-30;
1362    double overallSmallest=1.0e30;
1363    if (scalingMethod==1) {
1364      // Maximum in each row
1365      for (iRow=0;iRow<numberRows;iRow++) {
1366        if (usefulRow[iRow]) {
1367          CoinBigIndex j;
1368          largest=1.0e-20;
1369          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1370            int iColumn = column[j];
1371            if (usefulColumn[iColumn]) {
1372              double value = fabs(element[j]);
1373              largest = CoinMax(largest,value);
1374            }
1375          }
1376          rowScale[iRow]=1.0/largest;
1377          overallLargest = CoinMax(overallLargest,largest);
1378          overallSmallest = CoinMin(overallSmallest,largest);
1379        }
1380      }
1381    } else {
1382      assert(scalingMethod==2);
1383      int numberPass=3;
1384#ifdef USE_OBJECTIVE
1385      // This will be used to help get scale factors
1386      double * objective = new double[numberColumns];
1387      memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
1388      double objScale=1.0;
1389#endif
1390      while (numberPass) {
1391        overallLargest=0.0;
1392        overallSmallest=1.0e50;
1393        numberPass--;
1394        // Geometric mean on row scales
1395        for (iRow=0;iRow<numberRows;iRow++) {
1396          if (usefulRow[iRow]) {
1397            CoinBigIndex j;
1398            largest=1.0e-20;
1399            smallest=1.0e50;
1400            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1401              int iColumn = column[j];
1402              if (usefulColumn[iColumn]) {
1403                double value = fabs(element[j]);
1404                // Don't bother with tiny elements
1405                if (value>1.0e-30) {
1406                  value *= columnScale[iColumn];
1407                  largest = CoinMax(largest,value);
1408                  smallest = CoinMin(smallest,value);
1409                }
1410              }
1411            }
1412            rowScale[iRow]=1.0/sqrt(smallest*largest);
1413            overallLargest = CoinMax(largest*rowScale[iRow],overallLargest);
1414            overallSmallest = CoinMin(smallest*rowScale[iRow],overallSmallest);
1415          }
1416        }
1417#ifdef USE_OBJECTIVE
1418        largest=1.0e-20;
1419        smallest=1.0e50;
1420        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1421          if (usefulColumn[iColumn]) {
1422            double value = fabs(objective[iColumn]);
1423            // Don't bother with tiny elements
1424            if (value>1.0e-30) {
1425              value *= columnScale[iColumn];
1426              largest = CoinMax(largest,value);
1427              smallest = CoinMin(smallest,value);
1428            }
1429          }
1430        }
1431        objScale=1.0/sqrt(smallest*largest);
1432#endif
1433        model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
1434          <<overallSmallest
1435          <<overallLargest
1436          <<CoinMessageEol;
1437        // skip last column round
1438        if (numberPass==1)
1439          break;
1440        // Geometric mean on column scales
1441        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1442          if (usefulColumn[iColumn]) {
1443            CoinBigIndex j;
1444            largest=1.0e-20;
1445            smallest=1.0e50;
1446            for (j=columnStart[iColumn];
1447                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1448              iRow=row[j];
1449              double value = fabs(elementByColumn[j]);
1450              // Don't bother with tiny elements
1451              if (value>1.0e-30&&usefulRow[iRow]) {
1452                value *= rowScale[iRow];
1453                largest = CoinMax(largest,value);
1454                smallest = CoinMin(smallest,value);
1455              }
1456            }
1457#ifdef USE_OBJECTIVE
1458            if (fabs(objective[iColumn])>1.0e-30) {
1459              double value = fabs(objective[iColumn])*objScale;
1460              largest = CoinMax(largest,value);
1461              smallest = CoinMin(smallest,value);
1462            }
1463#endif
1464            columnScale[iColumn]=1.0/sqrt(smallest*largest);
1465          }
1466        }
1467      }
1468#ifdef USE_OBJECTIVE
1469      delete [] objective;
1470      printf("obj scale %g - use it if you want\n",objScale);
1471#endif
1472    }
1473    // If ranges will make horrid then scale
1474    double tolerance = 5.0*model->primalTolerance();
1475    for (iRow=0;iRow<numberRows;iRow++) {
1476      if (usefulRow[iRow]) {
1477        double difference = rowUpper[iRow]-rowLower[iRow];
1478        double scaledDifference = difference*rowScale[iRow];
1479        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1480          // make gap larger
1481          rowScale[iRow] *= 1.0e-4/scaledDifference;
1482          //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
1483          // scaledDifference,difference*rowScale[iRow]);
1484        }
1485      }
1486    }
1487    // final pass to scale columns so largest is reasonable
1488    // See what smallest will be if largest is 1.0
1489    overallSmallest=1.0e50;
1490    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1491      if (usefulColumn[iColumn]) {
1492        CoinBigIndex j;
1493        largest=1.0e-20;
1494        smallest=1.0e50;
1495        for (j=columnStart[iColumn];
1496             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1497          iRow=row[j];
1498          if(elementByColumn[j]&&usefulRow[iRow]) {
1499            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1500            largest = CoinMax(largest,value);
1501            smallest = CoinMin(smallest,value);
1502          }
1503        }
1504        if (overallSmallest*largest>smallest)
1505          overallSmallest = smallest/largest;
1506      }
1507    }
1508    //#define RANDOMIZE
1509#ifdef RANDOMIZE
1510    // randomize by up to 10%
1511    for (iRow=0;iRow<numberRows;iRow++) {
1512      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1513      rowScale[iRow] *= (1.0+0.1*value);
1514    }
1515#endif
1516    overallLargest=1.0;
1517    if (overallSmallest<1.0e-1)
1518      overallLargest = 1.0/sqrt(overallSmallest);
1519    overallLargest = CoinMin(100.0,overallLargest);
1520    overallSmallest=1.0e50;
1521    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1522      if (usefulColumn[iColumn]) {
1523        CoinBigIndex j;
1524        largest=1.0e-20;
1525        smallest=1.0e50;
1526        for (j=columnStart[iColumn];
1527             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1528          iRow=row[j];
1529          if(elementByColumn[j]&&usefulRow[iRow]) {
1530            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1531            largest = CoinMax(largest,value);
1532            smallest = CoinMin(smallest,value);
1533          }
1534        }
1535        columnScale[iColumn]=overallLargest/largest;
1536#ifdef RANDOMIZE
1537        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1538        columnScale[iColumn] *= (1.0+0.1*value);
1539#endif
1540        double difference = columnUpper[iColumn]-columnLower[iColumn];
1541        double scaledDifference = difference*columnScale[iColumn];
1542        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1543          // make gap larger
1544          columnScale[iColumn] *= 1.0e-4/scaledDifference;
1545          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
1546          // scaledDifference,difference*columnScale[iColumn]);
1547        }
1548        overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
1549      }
1550    }
1551    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
1552      <<overallSmallest
1553      <<overallLargest
1554      <<CoinMessageEol;
1555    delete [] usefulRow;
1556    delete [] usefulColumn;
1557    // If quadratic then make symmetric
1558    ClpObjective * obj = model->objectiveAsObject();
1559#ifndef NO_RTTI
1560    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
1561#else
1562    ClpQuadraticObjective * quadraticObj = NULL;
1563    if (obj->type()==2)
1564      quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
1565#endif
1566    if (quadraticObj) {
1567      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1568      int numberXColumns = quadratic->getNumCols();
1569      if (numberXColumns<numberColumns) {
1570        // we assume symmetric
1571        int numberQuadraticColumns=0;
1572        int i;
1573        //const int * columnQuadratic = quadratic->getIndices();
1574        //const int * columnQuadraticStart = quadratic->getVectorStarts();
1575        const int * columnQuadraticLength = quadratic->getVectorLengths();
1576        for (i=0;i<numberXColumns;i++) {
1577          int length=columnQuadraticLength[i];
1578#ifndef CORRECT_COLUMN_COUNTS
1579          length=1;
1580#endif
1581          if (length)
1582            numberQuadraticColumns++;
1583        }
1584        int numberXRows = numberRows-numberQuadraticColumns;
1585        numberQuadraticColumns=0;
1586        for (i=0;i<numberXColumns;i++) { 
1587          int length=columnQuadraticLength[i];
1588#ifndef CORRECT_COLUMN_COUNTS
1589          length=1;
1590#endif
1591          if (length) {
1592            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
1593            numberQuadraticColumns++;
1594          }
1595        }   
1596        int numberQuadraticRows=0;
1597        for (i=0;i<numberXRows;i++) {
1598          // See if any in row quadratic
1599          CoinBigIndex j;
1600          int numberQ=0;
1601          for (j=rowStart[i];j<rowStart[i+1];j++) {
1602            int iColumn = column[j];
1603            if (columnQuadraticLength[iColumn])
1604              numberQ++;
1605          }
1606#ifndef CORRECT_ROW_COUNTS
1607          numberQ=1;
1608#endif
1609          if (numberQ) {
1610            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
1611            numberQuadraticRows++;
1612          }
1613        }
1614        // and make sure Sj okay
1615        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
1616          CoinBigIndex j=columnStart[iColumn];
1617          assert(columnLength[iColumn]==1);
1618          int iRow=row[j];
1619          double value = fabs(elementByColumn[j]*rowScale[iRow]);
1620          columnScale[iColumn]=1.0/value;
1621        }
1622      }
1623    }
1624    model->setRowScale(rowScale);
1625    model->setColumnScale(columnScale);
1626    if (model->rowCopy()) {
1627      // need to replace row by row
1628      double * newElement = new double[numberColumns];
1629      // scale row copy
1630      for (iRow=0;iRow<numberRows;iRow++) {
1631        CoinBigIndex j;
1632        double scale = rowScale[iRow];
1633        const double * elementsInThisRow = element + rowStart[iRow];
1634        const int * columnsInThisRow = column + rowStart[iRow];
1635        int number = rowStart[iRow+1]-rowStart[iRow];
1636        assert (number<=numberColumns);
1637        for (j=0;j<number;j++) {
1638          int iColumn = columnsInThisRow[j];
1639          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
1640        }
1641        rowCopy->replaceVector(iRow,number,newElement);
1642      }
1643      delete [] newElement;
1644    } else {
1645      // no row copy
1646      delete rowCopyBase;
1647    }
1648    return 0;
1649  }
1650}
1651/* Unpacks a column into an CoinIndexedvector
1652 */
1653void 
1654ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
1655                   int iColumn) const 
1656{
1657  const double * rowScale = model->rowScale();
1658  const int * row = matrix_->getIndices();
1659  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1660  const int * columnLength = matrix_->getVectorLengths(); 
1661  const double * elementByColumn = matrix_->getElements();
1662  CoinBigIndex i;
1663  if (!rowScale) {
1664    for (i=columnStart[iColumn];
1665         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1666      rowArray->add(row[i],elementByColumn[i]);
1667    }
1668  } else {
1669    // apply scaling
1670    double scale = model->columnScale()[iColumn];
1671    for (i=columnStart[iColumn];
1672         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1673      int iRow = row[i];
1674      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1675    }
1676  }
1677}
1678/* Unpacks a column into a CoinIndexedVector
1679** in packed format
1680Note that model is NOT const.  Bounds and objective could
1681be modified if doing column generation (just for this variable) */
1682void 
1683ClpPackedMatrix::unpackPacked(ClpSimplex * model,
1684                            CoinIndexedVector * rowArray,
1685                            int iColumn) const
1686{
1687  const double * rowScale = model->rowScale();
1688  const int * row = matrix_->getIndices();
1689  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1690  const int * columnLength = matrix_->getVectorLengths(); 
1691  const double * elementByColumn = matrix_->getElements();
1692  CoinBigIndex i;
1693  if (!rowScale) {
1694    CoinBigIndex j=columnStart[iColumn];
1695    rowArray->createPacked(columnLength[iColumn],
1696                           row+j,elementByColumn+j);
1697  } else {
1698    // apply scaling
1699    double scale = model->columnScale()[iColumn];
1700    int * index = rowArray->getIndices();
1701    double * array = rowArray->denseVector();
1702    int number = 0;
1703    for (i=columnStart[iColumn];
1704         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1705      int iRow = row[i];
1706      array[number]=elementByColumn[i]*scale*rowScale[iRow];
1707      index[number++]=iRow;
1708    }
1709    rowArray->setNumElements(number);
1710    rowArray->setPackedMode(true);
1711  }
1712}
1713/* Adds multiple of a column into an CoinIndexedvector
1714      You can use quickAdd to add to vector */
1715void 
1716ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1717                   int iColumn, double multiplier) const 
1718{
1719  const double * rowScale = model->rowScale();
1720  const int * row = matrix_->getIndices();
1721  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1722  const int * columnLength = matrix_->getVectorLengths(); 
1723  const double * elementByColumn = matrix_->getElements();
1724  CoinBigIndex i;
1725  if (!rowScale) {
1726    for (i=columnStart[iColumn];
1727         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1728      int iRow = row[i];
1729      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
1730    }
1731  } else {
1732    // apply scaling
1733    double scale = model->columnScale()[iColumn]*multiplier;
1734    for (i=columnStart[iColumn];
1735         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1736      int iRow = row[i];
1737      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1738    }
1739  }
1740}
1741/* Adds multiple of a column into an array */
1742void 
1743ClpPackedMatrix::add(const ClpSimplex * model,double * array,
1744                    int iColumn, double multiplier) const
1745{
1746  const double * rowScale = model->rowScale();
1747  const int * row = matrix_->getIndices();
1748  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1749  const int * columnLength = matrix_->getVectorLengths(); 
1750  const double * elementByColumn = matrix_->getElements();
1751  CoinBigIndex i;
1752  if (!rowScale) {
1753    for (i=columnStart[iColumn];
1754         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1755      int iRow = row[i];
1756      array[iRow] += multiplier*elementByColumn[i];
1757    }
1758  } else {
1759    // apply scaling
1760    double scale = model->columnScale()[iColumn]*multiplier;
1761    for (i=columnStart[iColumn];
1762         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1763      int iRow = row[i];
1764      array[iRow] += elementByColumn[i]*scale*rowScale[iRow];
1765    }
1766  }
1767}
1768/* Checks if all elements are in valid range.  Can just
1769   return true if you are not paranoid.  For Clp I will
1770   probably expect no zeros.  Code can modify matrix to get rid of
1771   small elements.
1772   check bits (can be turned off to save time) :
1773   1 - check if matrix has gaps
1774   2 - check if zero elements
1775   4 - check and compress duplicates
1776   8 - report on large and small
1777*/
1778bool 
1779ClpPackedMatrix::allElementsInRange(ClpModel * model,
1780                                    double smallest, double largest,
1781                                    int check)
1782{
1783  int iColumn;
1784  // make sure matrix correct size
1785  matrix_->setDimensions(model->numberRows(),model->numberColumns());
1786  CoinBigIndex numberLarge=0;;
1787  CoinBigIndex numberSmall=0;;
1788  CoinBigIndex numberDuplicate=0;;
1789  int firstBadColumn=-1;
1790  int firstBadRow=-1;
1791  double firstBadElement=0.0;
1792  // get matrix data pointers
1793  const int * row = matrix_->getIndices();
1794  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1795  const int * columnLength = matrix_->getVectorLengths(); 
1796  const double * elementByColumn = matrix_->getElements();
1797  int numberRows = matrix_->getNumRows();
1798  int numberColumns = matrix_->getNumCols();
1799  // Say no gaps
1800  hasGaps_=false;
1801  if (check==14) {
1802    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1803      CoinBigIndex start = columnStart[iColumn];
1804      CoinBigIndex end = start + columnLength[iColumn];
1805      if (end!=columnStart[iColumn+1]) {
1806        hasGaps_=true;
1807        break;
1808      }
1809    }
1810    return true;
1811  }
1812  assert (check=15);
1813  int * mark = new int [numberRows];
1814  int i;
1815  for (i=0;i<numberRows;i++)
1816    mark[i]=-1;
1817  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1818    CoinBigIndex j;
1819    CoinBigIndex start = columnStart[iColumn];
1820    CoinBigIndex end = start + columnLength[iColumn];
1821    if (end!=columnStart[iColumn+1])
1822      hasGaps_=true;
1823    for (j=start;j<end;j++) {
1824      double value = fabs(elementByColumn[j]);
1825      int iRow = row[j];
1826      if (iRow<0||iRow>=numberRows) {
1827#ifndef COIN_BIG_INDEX
1828        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1829#elif COIN_BIG_INDEX==1
1830        printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1831#else
1832        printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1833#endif
1834        return false;
1835      }
1836      if (mark[iRow]==-1) {
1837        mark[iRow]=j;
1838      } else {
1839        // duplicate
1840        numberDuplicate++;
1841      }
1842      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1843      if (!value)
1844        zeroElements_ = true; // there are zero elements
1845      if (value<smallest) {
1846        numberSmall++;
1847      } else if (value>largest) {
1848        numberLarge++;
1849        if (firstBadColumn<0) {
1850          firstBadColumn=iColumn;
1851          firstBadRow=row[j];
1852          firstBadElement=elementByColumn[j];
1853        }
1854      }
1855    }
1856    //clear mark
1857    for (j=columnStart[iColumn];
1858         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1859      int iRow = row[j];
1860      mark[iRow]=-1;
1861    }
1862  }
1863  delete [] mark;
1864  if (numberLarge) {
1865    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
1866      <<numberLarge
1867      <<firstBadColumn<<firstBadRow<<firstBadElement
1868      <<CoinMessageEol;
1869    return false;
1870  }
1871  if (numberSmall) 
1872    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
1873      <<numberSmall
1874      <<CoinMessageEol;
1875  if (numberDuplicate) 
1876    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
1877      <<numberDuplicate
1878      <<CoinMessageEol;
1879  if (numberDuplicate) 
1880    matrix_->eliminateDuplicates(smallest);
1881  else if (numberSmall) 
1882    matrix_->compress(smallest);
1883  // If smallest >0.0 then there can't be zero elements
1884  if (smallest>0.0)
1885    zeroElements_=false;
1886  return true;
1887}
1888/* Given positive integer weights for each row fills in sum of weights
1889   for each column (and slack).
1890   Returns weights vector
1891*/
1892CoinBigIndex * 
1893ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
1894{
1895  int numberRows = model->numberRows();
1896  int numberColumns = matrix_->getNumCols();
1897  int number = numberRows+numberColumns;
1898  CoinBigIndex * weights = new CoinBigIndex[number];
1899  // get matrix data pointers
1900  const int * row = matrix_->getIndices();
1901  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1902  const int * columnLength = matrix_->getVectorLengths(); 
1903  int i;
1904  for (i=0;i<numberColumns;i++) {
1905    CoinBigIndex j;
1906    CoinBigIndex count=0;
1907    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1908      int iRow=row[j];
1909      count += inputWeights[iRow];
1910    }
1911    weights[i]=count;
1912  }
1913  for (i=0;i<numberRows;i++) {
1914    weights[i+numberColumns]=inputWeights[i];
1915  }
1916  return weights;
1917}
1918/* Returns largest and smallest elements of both signs.
1919   Largest refers to largest absolute value.
1920*/
1921void 
1922ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1923                       double & smallestPositive, double & largestPositive)
1924{
1925  smallestNegative=-COIN_DBL_MAX;
1926  largestNegative=0.0;
1927  smallestPositive=COIN_DBL_MAX;
1928  largestPositive=0.0;
1929  // get matrix data pointers
1930  const double * elementByColumn = matrix_->getElements();
1931  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1932  const int * columnLength = matrix_->getVectorLengths(); 
1933  int numberColumns = matrix_->getNumCols();
1934  int i;
1935  for (i=0;i<numberColumns;i++) {
1936    CoinBigIndex j;
1937    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1938      double value = elementByColumn[j];
1939      if (value>0.0) {
1940        smallestPositive = CoinMin(smallestPositive,value);
1941        largestPositive = CoinMax(largestPositive,value);
1942      } else if (value<0.0) {
1943        smallestNegative = CoinMax(smallestNegative,value);
1944        largestNegative = CoinMin(largestNegative,value);
1945      }
1946    }
1947  }
1948}
1949// Says whether it can do partial pricing
1950bool 
1951ClpPackedMatrix::canDoPartialPricing() const
1952{
1953  return true;
1954}
1955// Partial pricing
1956void 
1957ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1958                              int & bestSequence, int & numberWanted)
1959{
1960  numberWanted=currentWanted_;
1961  int start = (int) (startFraction*numberActiveColumns_);
1962  int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
1963  const double * element =matrix_->getElements();
1964  const int * row = matrix_->getIndices();
1965  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1966  const int * length = matrix_->getVectorLengths();
1967  const double * rowScale = model->rowScale();
1968  const double * columnScale = model->columnScale();
1969  int iSequence;
1970  CoinBigIndex j;
1971  double tolerance=model->currentDualTolerance();
1972  double * reducedCost = model->djRegion();
1973  const double * duals = model->dualRowSolution();
1974  const double * cost = model->costRegion();
1975  double bestDj;
1976  if (bestSequence>=0)
1977    bestDj = fabs(model->clpMatrix()->reducedCost(model,bestSequence));
1978  else
1979    bestDj=tolerance;
1980  int sequenceOut = model->sequenceOut();
1981  int saveSequence = bestSequence;
1982  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 
1983  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
1984  if (rowScale) {
1985    // scaled
1986    for (iSequence=start;iSequence<end;iSequence++) {
1987      if (iSequence!=sequenceOut) {
1988        double value;
1989        ClpSimplex::Status status = model->getStatus(iSequence);
1990       
1991        switch(status) {
1992         
1993        case ClpSimplex::basic:
1994        case ClpSimplex::isFixed:
1995          break;
1996        case ClpSimplex::isFree:
1997        case ClpSimplex::superBasic:
1998          value=0.0;
1999          // scaled
2000          for (j=startColumn[iSequence];
2001               j<startColumn[iSequence]+length[iSequence];j++) {
2002            int jRow=row[j];
2003            value -= duals[jRow]*element[j]*rowScale[jRow];
2004          }
2005          value = fabs(cost[iSequence] +value*columnScale[iSequence]);
2006          if (value>FREE_ACCEPT*tolerance) {
2007            numberWanted--;
2008            // we are going to bias towards free (but only if reasonable)
2009            value *= FREE_BIAS;
2010            if (value>bestDj) {
2011              // check flagged variable and correct dj
2012              if (!model->flagged(iSequence)) {
2013                bestDj=value;
2014                bestSequence = iSequence;
2015              } else {
2016                // just to make sure we don't exit before got something
2017                numberWanted++;
2018              }
2019            }
2020          }
2021          break;
2022        case ClpSimplex::atUpperBound:
2023          value=0.0;
2024          // scaled
2025          for (j=startColumn[iSequence];
2026               j<startColumn[iSequence]+length[iSequence];j++) {
2027            int jRow=row[j];
2028            value -= duals[jRow]*element[j]*rowScale[jRow];
2029          }
2030          value = cost[iSequence] +value*columnScale[iSequence];
2031          if (value>tolerance) {
2032            numberWanted--;
2033            if (value>bestDj) {
2034              // check flagged variable and correct dj
2035              if (!model->flagged(iSequence)) {
2036                bestDj=value;
2037                bestSequence = iSequence;
2038              } else {
2039                // just to make sure we don't exit before got something
2040                numberWanted++;
2041              }
2042            }
2043          }
2044          break;
2045        case ClpSimplex::atLowerBound:
2046          value=0.0;
2047          // scaled
2048          for (j=startColumn[iSequence];
2049               j<startColumn[iSequence]+length[iSequence];j++) {
2050            int jRow=row[j];
2051            value -= duals[jRow]*element[j]*rowScale[jRow];
2052          }
2053          value = -(cost[iSequence] +value*columnScale[iSequence]);
2054          if (value>tolerance) {
2055            numberWanted--;
2056            if (value>bestDj) {
2057              // check flagged variable and correct dj
2058              if (!model->flagged(iSequence)) {
2059                bestDj=value;
2060                bestSequence = iSequence;
2061              } else {
2062                // just to make sure we don't exit before got something
2063                numberWanted++;
2064              }
2065            }
2066          }
2067          break;
2068        }
2069      }
2070      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2071        // give up
2072        break;
2073      }
2074      if (!numberWanted)
2075        break;
2076    }
2077    if (bestSequence!=saveSequence) {
2078      // recompute dj
2079      double value=0.0;
2080      // scaled
2081      for (j=startColumn[bestSequence];
2082           j<startColumn[bestSequence]+length[bestSequence];j++) {
2083        int jRow=row[j];
2084        value -= duals[jRow]*element[j]*rowScale[jRow];
2085      }
2086      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
2087      savedBestSequence_ = bestSequence;
2088      savedBestDj_ = reducedCost[savedBestSequence_];
2089    }
2090  } else {
2091    // not scaled
2092    for (iSequence=start;iSequence<end;iSequence++) {
2093      if (iSequence!=sequenceOut) {
2094        double value;
2095        ClpSimplex::Status status = model->getStatus(iSequence);
2096       
2097        switch(status) {
2098         
2099        case ClpSimplex::basic:
2100        case ClpSimplex::isFixed:
2101          break;
2102        case ClpSimplex::isFree:
2103        case ClpSimplex::superBasic:
2104          value=cost[iSequence];
2105          for (j=startColumn[iSequence];
2106               j<startColumn[iSequence]+length[iSequence];j++) {
2107            int jRow=row[j];
2108            value -= duals[jRow]*element[j];
2109          }
2110          value = fabs(value);
2111          if (value>FREE_ACCEPT*tolerance) {
2112            numberWanted--;
2113            // we are going to bias towards free (but only if reasonable)
2114            value *= FREE_BIAS;
2115            if (value>bestDj) {
2116              // check flagged variable and correct dj
2117              if (!model->flagged(iSequence)) {
2118                bestDj=value;
2119                bestSequence = iSequence;
2120              } else {
2121                // just to make sure we don't exit before got something
2122                numberWanted++;
2123              }
2124            }
2125          }
2126          break;
2127        case ClpSimplex::atUpperBound:
2128          value=cost[iSequence];
2129          // scaled
2130          for (j=startColumn[iSequence];
2131               j<startColumn[iSequence]+length[iSequence];j++) {
2132            int jRow=row[j];
2133            value -= duals[jRow]*element[j];
2134          }
2135          if (value>tolerance) {
2136            numberWanted--;
2137            if (value>bestDj) {
2138              // check flagged variable and correct dj
2139              if (!model->flagged(iSequence)) {
2140                bestDj=value;
2141                bestSequence = iSequence;
2142              } else {
2143                // just to make sure we don't exit before got something
2144                numberWanted++;
2145              }
2146            }
2147          }
2148          break;
2149        case ClpSimplex::atLowerBound:
2150          value=cost[iSequence];
2151          for (j=startColumn[iSequence];
2152               j<startColumn[iSequence]+length[iSequence];j++) {
2153            int jRow=row[j];
2154            value -= duals[jRow]*element[j];
2155          }
2156          value = -value;
2157          if (value>tolerance) {
2158            numberWanted--;
2159            if (value>bestDj) {
2160              // check flagged variable and correct dj
2161              if (!model->flagged(iSequence)) {
2162                bestDj=value;
2163                bestSequence = iSequence;
2164              } else {
2165                // just to make sure we don't exit before got something
2166                numberWanted++;
2167              }
2168            }
2169          }
2170          break;
2171        }
2172      }
2173      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2174        // give up
2175        break;
2176      }
2177      if (!numberWanted)
2178        break;
2179    }
2180    if (bestSequence!=saveSequence) {
2181      // recompute dj
2182      double value=cost[bestSequence];
2183      for (j=startColumn[bestSequence];
2184           j<startColumn[bestSequence]+length[bestSequence];j++) {
2185        int jRow=row[j];
2186        value -= duals[jRow]*element[j];
2187      }
2188      reducedCost[bestSequence] = value;
2189      savedBestSequence_ = bestSequence;
2190      savedBestDj_ = reducedCost[savedBestSequence_];
2191    }
2192  }
2193  currentWanted_=numberWanted;
2194}
2195// Sets up an effective RHS
2196void 
2197ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
2198{
2199  delete [] rhsOffset_;
2200  int numberRows = model->numberRows();
2201  rhsOffset_ = new double[numberRows];
2202  rhsOffset(model,true);
2203}
2204// makes sure active columns correct
2205int 
2206ClpPackedMatrix::refresh(ClpSimplex * model)
2207{
2208  numberActiveColumns_ = matrix_->getNumCols();
2209  return 0;
2210}
2211
2212/* Scales rowCopy if column copy scaled
2213   Only called if scales already exist */
2214void 
2215ClpPackedMatrix::scaleRowCopy(ClpModel * model) const 
2216{
2217  if (model->rowCopy()) {
2218    // need to replace row by row
2219    int numberRows = model->numberRows();
2220    int numberColumns = matrix_->getNumCols();
2221    double * newElement = new double[numberColumns];
2222    ClpMatrixBase * rowCopyBase=model->rowCopy();
2223#ifndef NO_RTTI
2224    ClpPackedMatrix* rowCopy =
2225      dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
2226    // Make sure it is really a ClpPackedMatrix
2227    assert (rowCopy!=NULL);
2228#else
2229    ClpPackedMatrix* rowCopy =
2230      static_cast< ClpPackedMatrix*>(rowCopyBase);
2231#endif
2232
2233    const int * column = rowCopy->getIndices();
2234    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2235    const double * element = rowCopy->getElements();
2236    const double * rowScale = model->rowScale();
2237    const double * columnScale = model->columnScale();
2238    // scale row copy
2239    for (int iRow=0;iRow<numberRows;iRow++) {
2240      CoinBigIndex j;
2241      double scale = rowScale[iRow];
2242      const double * elementsInThisRow = element + rowStart[iRow];
2243      const int * columnsInThisRow = column + rowStart[iRow];
2244      int number = rowStart[iRow+1]-rowStart[iRow];
2245      assert (number<=numberColumns);
2246      for (j=0;j<number;j++) {
2247        int iColumn = columnsInThisRow[j];
2248        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
2249      }
2250      rowCopy->replaceVector(iRow,number,newElement);
2251    }
2252    delete [] newElement;
2253  }
2254}
2255/* Realy really scales column copy
2256   Only called if scales already exist.
2257   Up to user ro delete */
2258ClpMatrixBase * 
2259ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const 
2260{
2261  // need to replace column by column
2262  int numberRows = model->numberRows();
2263  int numberColumns = matrix_->getNumCols();
2264  double * newElement = new double[numberRows];
2265  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
2266  const int * row = copy->getIndices();
2267  const CoinBigIndex * columnStart = copy->getVectorStarts();
2268  const int * length = copy->getVectorLengths();
2269  const double * element = copy->getElements();
2270  const double * rowScale = model->rowScale();
2271  const double * columnScale = model->columnScale();
2272  // scale column copy
2273  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2274    CoinBigIndex j;
2275    double scale = columnScale[iColumn];
2276    const double * elementsInThisColumn = element + columnStart[iColumn];
2277    const int * rowsInThisColumn = row + columnStart[iColumn];
2278    int number = length[iColumn];
2279    assert (number<=numberRows);
2280    for (j=0;j<number;j++) {
2281      int iRow = rowsInThisColumn[j];
2282      newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
2283    }
2284    copy->replaceVector(iColumn,number,newElement);
2285  }
2286  delete [] newElement;
2287  return copy;
2288}
2289// Really scale matrix
2290void 
2291ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
2292{
2293  int numberColumns = matrix_->getNumCols();
2294  const int * row = matrix_->getIndices();
2295  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2296  const int * length = matrix_->getVectorLengths();
2297  double * element = matrix_->getMutableElements();
2298  // scale
2299  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2300    CoinBigIndex j;
2301    double scale = columnScale[iColumn];
2302    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
2303      int iRow = row[j];
2304      element[j] *= scale*rowScale[iRow];
2305    }
2306  }
2307}
Note: See TracBrowser for help on using the repository browser.