source: trunk/ClpPackedMatrix.cpp @ 461

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

Trying to make faster

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