source: trunk/ClpPackedMatrix.cpp @ 466

Last change on this file since 466 was 466, checked in by forrest, 15 years ago

deleting all and restoring all

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 67.4 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*/
1773bool 
1774ClpPackedMatrix::allElementsInRange(ClpModel * model,
1775                                    double smallest, double largest)
1776{
1777  int iColumn;
1778  // make sure matrix correct size
1779  matrix_->setDimensions(model->numberRows(),model->numberColumns());
1780  CoinBigIndex numberLarge=0;;
1781  CoinBigIndex numberSmall=0;;
1782  CoinBigIndex numberDuplicate=0;;
1783  int firstBadColumn=-1;
1784  int firstBadRow=-1;
1785  double firstBadElement=0.0;
1786  // get matrix data pointers
1787  const int * row = matrix_->getIndices();
1788  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1789  const int * columnLength = matrix_->getVectorLengths(); 
1790  const double * elementByColumn = matrix_->getElements();
1791  int numberRows = matrix_->getNumRows();
1792  int numberColumns = matrix_->getNumCols();
1793  int * mark = new int [numberRows];
1794  int i;
1795  // Say no gaps
1796  hasGaps_=false;
1797  for (i=0;i<numberRows;i++)
1798    mark[i]=-1;
1799  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1800    CoinBigIndex j;
1801    CoinBigIndex start = columnStart[iColumn];
1802    CoinBigIndex end = start + columnLength[iColumn];
1803    if (end!=columnStart[iColumn+1])
1804      hasGaps_=true;
1805    for (j=start;j<end;j++) {
1806      double value = fabs(elementByColumn[j]);
1807      int iRow = row[j];
1808      if (iRow<0||iRow>=numberRows) {
1809#ifndef COIN_BIG_INDEX
1810        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1811#elif COIN_BIG_INDEX==1
1812        printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1813#else
1814        printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1815#endif
1816        return false;
1817      }
1818      if (mark[iRow]==-1) {
1819        mark[iRow]=j;
1820      } else {
1821        // duplicate
1822        numberDuplicate++;
1823      }
1824      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1825      if (!value)
1826        zeroElements_ = true; // there are zero elements
1827      if (value<smallest) {
1828        numberSmall++;
1829      } else if (value>largest) {
1830        numberLarge++;
1831        if (firstBadColumn<0) {
1832          firstBadColumn=iColumn;
1833          firstBadRow=row[j];
1834          firstBadElement=elementByColumn[j];
1835        }
1836      }
1837    }
1838    //clear mark
1839    for (j=columnStart[iColumn];
1840         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1841      int iRow = row[j];
1842      mark[iRow]=-1;
1843    }
1844  }
1845  delete [] mark;
1846  if (numberLarge) {
1847    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
1848      <<numberLarge
1849      <<firstBadColumn<<firstBadRow<<firstBadElement
1850      <<CoinMessageEol;
1851    return false;
1852  }
1853  if (numberSmall) 
1854    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
1855      <<numberSmall
1856      <<CoinMessageEol;
1857  if (numberDuplicate) 
1858    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
1859      <<numberDuplicate
1860      <<CoinMessageEol;
1861  if (numberDuplicate) 
1862    matrix_->eliminateDuplicates(smallest);
1863  else if (numberSmall) 
1864    matrix_->compress(smallest);
1865  // If smallest >0.0 then there can't be zero elements
1866  if (smallest>0.0)
1867    zeroElements_=false;
1868  return true;
1869}
1870/* Given positive integer weights for each row fills in sum of weights
1871   for each column (and slack).
1872   Returns weights vector
1873*/
1874CoinBigIndex * 
1875ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
1876{
1877  int numberRows = model->numberRows();
1878  int numberColumns = matrix_->getNumCols();
1879  int number = numberRows+numberColumns;
1880  CoinBigIndex * weights = new CoinBigIndex[number];
1881  // get matrix data pointers
1882  const int * row = matrix_->getIndices();
1883  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1884  const int * columnLength = matrix_->getVectorLengths(); 
1885  int i;
1886  for (i=0;i<numberColumns;i++) {
1887    CoinBigIndex j;
1888    CoinBigIndex count=0;
1889    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1890      int iRow=row[j];
1891      count += inputWeights[iRow];
1892    }
1893    weights[i]=count;
1894  }
1895  for (i=0;i<numberRows;i++) {
1896    weights[i+numberColumns]=inputWeights[i];
1897  }
1898  return weights;
1899}
1900/* Returns largest and smallest elements of both signs.
1901   Largest refers to largest absolute value.
1902*/
1903void 
1904ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1905                       double & smallestPositive, double & largestPositive)
1906{
1907  smallestNegative=-COIN_DBL_MAX;
1908  largestNegative=0.0;
1909  smallestPositive=COIN_DBL_MAX;
1910  largestPositive=0.0;
1911  // get matrix data pointers
1912  const double * elementByColumn = matrix_->getElements();
1913  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1914  const int * columnLength = matrix_->getVectorLengths(); 
1915  int numberColumns = matrix_->getNumCols();
1916  int i;
1917  for (i=0;i<numberColumns;i++) {
1918    CoinBigIndex j;
1919    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1920      double value = elementByColumn[j];
1921      if (value>0.0) {
1922        smallestPositive = CoinMin(smallestPositive,value);
1923        largestPositive = CoinMax(largestPositive,value);
1924      } else if (value<0.0) {
1925        smallestNegative = CoinMax(smallestNegative,value);
1926        largestNegative = CoinMin(largestNegative,value);
1927      }
1928    }
1929  }
1930}
1931// Says whether it can do partial pricing
1932bool 
1933ClpPackedMatrix::canDoPartialPricing() const
1934{
1935  return true;
1936}
1937// Partial pricing
1938void 
1939ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1940                              int & bestSequence, int & numberWanted)
1941{
1942  numberWanted=currentWanted_;
1943  int start = (int) (startFraction*numberActiveColumns_);
1944  int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
1945  const double * element =matrix_->getElements();
1946  const int * row = matrix_->getIndices();
1947  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1948  const int * length = matrix_->getVectorLengths();
1949  const double * rowScale = model->rowScale();
1950  const double * columnScale = model->columnScale();
1951  int iSequence;
1952  CoinBigIndex j;
1953  double tolerance=model->currentDualTolerance();
1954  double * reducedCost = model->djRegion();
1955  const double * duals = model->dualRowSolution();
1956  const double * cost = model->costRegion();
1957  double bestDj;
1958  if (bestSequence>=0)
1959    bestDj = fabs(model->clpMatrix()->reducedCost(model,bestSequence));
1960  else
1961    bestDj=tolerance;
1962  int sequenceOut = model->sequenceOut();
1963  int saveSequence = bestSequence;
1964  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 
1965  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
1966  if (rowScale) {
1967    // scaled
1968    for (iSequence=start;iSequence<end;iSequence++) {
1969      if (iSequence!=sequenceOut) {
1970        double value;
1971        ClpSimplex::Status status = model->getStatus(iSequence);
1972       
1973        switch(status) {
1974         
1975        case ClpSimplex::basic:
1976        case ClpSimplex::isFixed:
1977          break;
1978        case ClpSimplex::isFree:
1979        case ClpSimplex::superBasic:
1980          value=0.0;
1981          // scaled
1982          for (j=startColumn[iSequence];
1983               j<startColumn[iSequence]+length[iSequence];j++) {
1984            int jRow=row[j];
1985            value -= duals[jRow]*element[j]*rowScale[jRow];
1986          }
1987          value = fabs(cost[iSequence] +value*columnScale[iSequence]);
1988          if (value>FREE_ACCEPT*tolerance) {
1989            numberWanted--;
1990            // we are going to bias towards free (but only if reasonable)
1991            value *= FREE_BIAS;
1992            if (value>bestDj) {
1993              // check flagged variable and correct dj
1994              if (!model->flagged(iSequence)) {
1995                bestDj=value;
1996                bestSequence = iSequence;
1997              } else {
1998                // just to make sure we don't exit before got something
1999                numberWanted++;
2000              }
2001            }
2002          }
2003          break;
2004        case ClpSimplex::atUpperBound:
2005          value=0.0;
2006          // scaled
2007          for (j=startColumn[iSequence];
2008               j<startColumn[iSequence]+length[iSequence];j++) {
2009            int jRow=row[j];
2010            value -= duals[jRow]*element[j]*rowScale[jRow];
2011          }
2012          value = cost[iSequence] +value*columnScale[iSequence];
2013          if (value>tolerance) {
2014            numberWanted--;
2015            if (value>bestDj) {
2016              // check flagged variable and correct dj
2017              if (!model->flagged(iSequence)) {
2018                bestDj=value;
2019                bestSequence = iSequence;
2020              } else {
2021                // just to make sure we don't exit before got something
2022                numberWanted++;
2023              }
2024            }
2025          }
2026          break;
2027        case ClpSimplex::atLowerBound:
2028          value=0.0;
2029          // scaled
2030          for (j=startColumn[iSequence];
2031               j<startColumn[iSequence]+length[iSequence];j++) {
2032            int jRow=row[j];
2033            value -= duals[jRow]*element[j]*rowScale[jRow];
2034          }
2035          value = -(cost[iSequence] +value*columnScale[iSequence]);
2036          if (value>tolerance) {
2037            numberWanted--;
2038            if (value>bestDj) {
2039              // check flagged variable and correct dj
2040              if (!model->flagged(iSequence)) {
2041                bestDj=value;
2042                bestSequence = iSequence;
2043              } else {
2044                // just to make sure we don't exit before got something
2045                numberWanted++;
2046              }
2047            }
2048          }
2049          break;
2050        }
2051      }
2052      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2053        // give up
2054        break;
2055      }
2056      if (!numberWanted)
2057        break;
2058    }
2059    if (bestSequence!=saveSequence) {
2060      // recompute dj
2061      double value=0.0;
2062      // scaled
2063      for (j=startColumn[bestSequence];
2064           j<startColumn[bestSequence]+length[bestSequence];j++) {
2065        int jRow=row[j];
2066        value -= duals[jRow]*element[j]*rowScale[jRow];
2067      }
2068      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
2069      savedBestSequence_ = bestSequence;
2070      savedBestDj_ = reducedCost[savedBestSequence_];
2071    }
2072  } else {
2073    // not scaled
2074    for (iSequence=start;iSequence<end;iSequence++) {
2075      if (iSequence!=sequenceOut) {
2076        double value;
2077        ClpSimplex::Status status = model->getStatus(iSequence);
2078       
2079        switch(status) {
2080         
2081        case ClpSimplex::basic:
2082        case ClpSimplex::isFixed:
2083          break;
2084        case ClpSimplex::isFree:
2085        case ClpSimplex::superBasic:
2086          value=cost[iSequence];
2087          for (j=startColumn[iSequence];
2088               j<startColumn[iSequence]+length[iSequence];j++) {
2089            int jRow=row[j];
2090            value -= duals[jRow]*element[j];
2091          }
2092          value = fabs(value);
2093          if (value>FREE_ACCEPT*tolerance) {
2094            numberWanted--;
2095            // we are going to bias towards free (but only if reasonable)
2096            value *= FREE_BIAS;
2097            if (value>bestDj) {
2098              // check flagged variable and correct dj
2099              if (!model->flagged(iSequence)) {
2100                bestDj=value;
2101                bestSequence = iSequence;
2102              } else {
2103                // just to make sure we don't exit before got something
2104                numberWanted++;
2105              }
2106            }
2107          }
2108          break;
2109        case ClpSimplex::atUpperBound:
2110          value=cost[iSequence];
2111          // scaled
2112          for (j=startColumn[iSequence];
2113               j<startColumn[iSequence]+length[iSequence];j++) {
2114            int jRow=row[j];
2115            value -= duals[jRow]*element[j];
2116          }
2117          if (value>tolerance) {
2118            numberWanted--;
2119            if (value>bestDj) {
2120              // check flagged variable and correct dj
2121              if (!model->flagged(iSequence)) {
2122                bestDj=value;
2123                bestSequence = iSequence;
2124              } else {
2125                // just to make sure we don't exit before got something
2126                numberWanted++;
2127              }
2128            }
2129          }
2130          break;
2131        case ClpSimplex::atLowerBound:
2132          value=cost[iSequence];
2133          for (j=startColumn[iSequence];
2134               j<startColumn[iSequence]+length[iSequence];j++) {
2135            int jRow=row[j];
2136            value -= duals[jRow]*element[j];
2137          }
2138          value = -value;
2139          if (value>tolerance) {
2140            numberWanted--;
2141            if (value>bestDj) {
2142              // check flagged variable and correct dj
2143              if (!model->flagged(iSequence)) {
2144                bestDj=value;
2145                bestSequence = iSequence;
2146              } else {
2147                // just to make sure we don't exit before got something
2148                numberWanted++;
2149              }
2150            }
2151          }
2152          break;
2153        }
2154      }
2155      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2156        // give up
2157        break;
2158      }
2159      if (!numberWanted)
2160        break;
2161    }
2162    if (bestSequence!=saveSequence) {
2163      // recompute dj
2164      double value=cost[bestSequence];
2165      for (j=startColumn[bestSequence];
2166           j<startColumn[bestSequence]+length[bestSequence];j++) {
2167        int jRow=row[j];
2168        value -= duals[jRow]*element[j];
2169      }
2170      reducedCost[bestSequence] = value;
2171      savedBestSequence_ = bestSequence;
2172      savedBestDj_ = reducedCost[savedBestSequence_];
2173    }
2174  }
2175  currentWanted_=numberWanted;
2176}
2177// Sets up an effective RHS
2178void 
2179ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
2180{
2181  delete [] rhsOffset_;
2182  int numberRows = model->numberRows();
2183  rhsOffset_ = new double[numberRows];
2184  rhsOffset(model,true);
2185}
2186// makes sure active columns correct
2187int 
2188ClpPackedMatrix::refresh(ClpSimplex * model)
2189{
2190  numberActiveColumns_ = matrix_->getNumCols();
2191  return 0;
2192}
2193
2194/* Scales rowCopy if column copy scaled
2195   Only called if scales already exist */
2196void 
2197ClpPackedMatrix::scaleRowCopy(ClpModel * model) const 
2198{
2199  if (model->rowCopy()) {
2200    // need to replace row by row
2201    int numberRows = model->numberRows();
2202    int numberColumns = matrix_->getNumCols();
2203    double * newElement = new double[numberColumns];
2204    ClpMatrixBase * rowCopyBase=model->rowCopy();
2205#ifndef NO_RTTI
2206    ClpPackedMatrix* rowCopy =
2207      dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
2208    // Make sure it is really a ClpPackedMatrix
2209    assert (rowCopy!=NULL);
2210#else
2211    ClpPackedMatrix* rowCopy =
2212      static_cast< ClpPackedMatrix*>(rowCopyBase);
2213#endif
2214
2215    const int * column = rowCopy->getIndices();
2216    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2217    const double * element = rowCopy->getElements();
2218    const double * rowScale = model->rowScale();
2219    const double * columnScale = model->columnScale();
2220    // scale row copy
2221    for (int iRow=0;iRow<numberRows;iRow++) {
2222      CoinBigIndex j;
2223      double scale = rowScale[iRow];
2224      const double * elementsInThisRow = element + rowStart[iRow];
2225      const int * columnsInThisRow = column + rowStart[iRow];
2226      int number = rowStart[iRow+1]-rowStart[iRow];
2227      assert (number<=numberColumns);
2228      for (j=0;j<number;j++) {
2229        int iColumn = columnsInThisRow[j];
2230        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
2231      }
2232      rowCopy->replaceVector(iRow,number,newElement);
2233    }
2234    delete [] newElement;
2235  }
2236}
2237/* Realy really scales column copy
2238   Only called if scales already exist.
2239   Up to user ro delete */
2240ClpMatrixBase * 
2241ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const 
2242{
2243  // need to replace column by column
2244  int numberRows = model->numberRows();
2245  int numberColumns = matrix_->getNumCols();
2246  double * newElement = new double[numberRows];
2247  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
2248  const int * row = copy->getIndices();
2249  const CoinBigIndex * columnStart = copy->getVectorStarts();
2250  const int * length = copy->getVectorLengths();
2251  const double * element = copy->getElements();
2252  const double * rowScale = model->rowScale();
2253  const double * columnScale = model->columnScale();
2254  // scale column copy
2255  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2256    CoinBigIndex j;
2257    double scale = columnScale[iColumn];
2258    const double * elementsInThisColumn = element + columnStart[iColumn];
2259    const int * rowsInThisColumn = row + columnStart[iColumn];
2260    int number = length[iColumn];
2261    assert (number<=numberRows);
2262    for (j=0;j<number;j++) {
2263      int iRow = rowsInThisColumn[j];
2264      newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
2265    }
2266    copy->replaceVector(iColumn,number,newElement);
2267  }
2268  delete [] newElement;
2269  return copy;
2270}
2271// Really scale matrix
2272void 
2273ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
2274{
2275  int numberColumns = matrix_->getNumCols();
2276  const int * row = matrix_->getIndices();
2277  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2278  const int * length = matrix_->getVectorLengths();
2279  double * element = matrix_->getMutableElements();
2280  // scale
2281  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2282    CoinBigIndex j;
2283    double scale = columnScale[iColumn];
2284    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
2285      int iRow = row[j];
2286      element[j] *= scale*rowScale[iRow];
2287    }
2288  }
2289}
Note: See TracBrowser for help on using the repository browser.