source: trunk/ClpPackedMatrix.cpp @ 225

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

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 50.8 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    zeroElements_(false)
29{
30  setType(1);
31}
32
33//-------------------------------------------------------------------
34// Copy constructor
35//-------------------------------------------------------------------
36ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs) 
37: ClpMatrixBase(rhs)
38{ 
39  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
40  zeroElements_ = rhs.zeroElements_;
41 
42}
43
44//-------------------------------------------------------------------
45// assign matrix (for space reasons)
46//-------------------------------------------------------------------
47ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs) 
48: ClpMatrixBase()
49{ 
50  matrix_ = rhs;
51  zeroElements_ = false;
52  setType(1);
53 
54}
55
56ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs) 
57: ClpMatrixBase()
58{ 
59  matrix_ = new CoinPackedMatrix(rhs);
60  zeroElements_ = false;
61  setType(1);
62 
63}
64
65//-------------------------------------------------------------------
66// Destructor
67//-------------------------------------------------------------------
68ClpPackedMatrix::~ClpPackedMatrix ()
69{
70  delete matrix_;
71}
72
73//----------------------------------------------------------------
74// Assignment operator
75//-------------------------------------------------------------------
76ClpPackedMatrix &
77ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
78{
79  if (this != &rhs) {
80    ClpMatrixBase::operator=(rhs);
81    delete matrix_;
82    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
83    zeroElements_ = rhs.zeroElements_;
84  }
85  return *this;
86}
87//-------------------------------------------------------------------
88// Clone
89//-------------------------------------------------------------------
90ClpMatrixBase * ClpPackedMatrix::clone() const
91{
92  return new ClpPackedMatrix(*this);
93}
94/* Subset clone (without gaps).  Duplicates are allowed
95   and order is as given */
96ClpMatrixBase * 
97ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
98                              int numberColumns, 
99                              const int * whichColumns) const 
100{
101  return new ClpPackedMatrix(*this, numberRows, whichRows,
102                                   numberColumns, whichColumns);
103}
104/* Subset constructor (without gaps).  Duplicates are allowed
105   and order is as given */
106ClpPackedMatrix::ClpPackedMatrix (
107                       const ClpPackedMatrix & rhs,
108                       int numberRows, const int * whichRows,
109                       int numberColumns, const int * whichColumns)
110: ClpMatrixBase(rhs)
111{
112  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
113                                 numberColumns,whichColumns);
114  zeroElements_ = rhs.zeroElements_;
115}
116ClpPackedMatrix::ClpPackedMatrix (
117                       const CoinPackedMatrix & rhs,
118                       int numberRows, const int * whichRows,
119                       int numberColumns, const int * whichColumns)
120: ClpMatrixBase()
121{
122  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
123                                 numberColumns,whichColumns);
124  zeroElements_ = false;
125  setType(1);
126}
127
128/* Returns a new matrix in reverse order without gaps */
129ClpMatrixBase * 
130ClpPackedMatrix::reverseOrderedCopy() const
131{
132  ClpPackedMatrix * copy = new ClpPackedMatrix();
133  copy->matrix_= new CoinPackedMatrix();
134  copy->matrix_->reverseOrderedCopyOf(*matrix_);
135  copy->matrix_->removeGaps();
136  return copy;
137}
138//unscaled versions
139void 
140ClpPackedMatrix::times(double scalar,
141                   const double * x, double * y) const
142{
143  int iRow,iColumn;
144  // get matrix data pointers
145  const int * row = matrix_->getIndices();
146  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
147  const int * columnLength = matrix_->getVectorLengths(); 
148  const double * elementByColumn = matrix_->getElements();
149  int numberColumns = matrix_->getNumCols();
150  //memset(y,0,matrix_->getNumRows()*sizeof(double));
151  for (iColumn=0;iColumn<numberColumns;iColumn++) {
152    CoinBigIndex j;
153    double value = scalar*x[iColumn];
154    if (value) {
155      for (j=columnStart[iColumn];
156           j<columnStart[iColumn]+columnLength[iColumn];j++) {
157        iRow=row[j];
158        y[iRow] += value*elementByColumn[j];
159      }
160    }
161  }
162}
163void 
164ClpPackedMatrix::transposeTimes(double scalar,
165                                const double * x, double * y) const
166{
167  int 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  int numberColumns = matrix_->getNumCols();
174  //memset(y,0,numberColumns*sizeof(double));
175  for (iColumn=0;iColumn<numberColumns;iColumn++) {
176    CoinBigIndex j;
177    double value=0.0;
178    for (j=columnStart[iColumn];
179         j<columnStart[iColumn]+columnLength[iColumn];j++) {
180      int jRow=row[j];
181      value += x[jRow]*elementByColumn[j];
182    }
183    y[iColumn] += value*scalar;
184  }
185}
186void 
187ClpPackedMatrix::times(double scalar,
188                       const double * x, double * y,
189                       const double * rowScale, 
190                       const double * columnScale) const
191{
192  if (rowScale) {
193    int iRow,iColumn;
194    // get matrix data pointers
195    const int * row = matrix_->getIndices();
196    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
197    const int * columnLength = matrix_->getVectorLengths(); 
198    const double * elementByColumn = matrix_->getElements();
199    int numberColumns = matrix_->getNumCols();
200    //memset(y,0,matrix_->getNumRows()*sizeof(double));
201    for (iColumn=0;iColumn<numberColumns;iColumn++) {
202      CoinBigIndex j;
203      double value = x[iColumn];
204      if (value) {
205        // scaled
206        value *= scalar*columnScale[iColumn];
207        for (j=columnStart[iColumn];
208             j<columnStart[iColumn]+columnLength[iColumn];j++) {
209          iRow=row[j];
210          y[iRow] += value*elementByColumn[j];
211        }
212      }
213    }
214    int numberRows = getNumRows();
215    for (iRow=0;iRow<numberRows;iRow++) {
216      y[iRow] *= rowScale[iRow];
217    }
218  } else {
219    times(scalar,x,y);
220  }
221}
222void 
223ClpPackedMatrix::transposeTimes( double scalar,
224                                 const double * x, double * y,
225                                 const double * rowScale, 
226                                 const double * columnScale) const
227{
228  if (rowScale) {
229    int iColumn;
230    // get matrix data pointers
231    const int * row = matrix_->getIndices();
232    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
233    const int * columnLength = matrix_->getVectorLengths(); 
234    const double * elementByColumn = matrix_->getElements();
235    int numberColumns = matrix_->getNumCols();
236    //memset(y,0,numberColumns*sizeof(double));
237    for (iColumn=0;iColumn<numberColumns;iColumn++) {
238      CoinBigIndex j;
239      double value=0.0;
240      // scaled
241      for (j=columnStart[iColumn];
242           j<columnStart[iColumn]+columnLength[iColumn];j++) {
243        int jRow=row[j];
244      value += x[jRow]*elementByColumn[j]*rowScale[jRow];
245      }
246      y[iColumn] += value*scalar*columnScale[iColumn];
247    }
248  } else {
249    transposeTimes(scalar,x,y);
250  }
251}
252/* Return <code>x * A + y</code> in <code>z</code>.
253        Squashes small elements and knows about ClpSimplex */
254void 
255ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
256                              const CoinIndexedVector * rowArray,
257                              CoinIndexedVector * y,
258                              CoinIndexedVector * columnArray) const
259{
260  columnArray->clear();
261  double * pi = rowArray->denseVector();
262  int numberNonZero=0;
263  int * index = columnArray->getIndices();
264  double * array = columnArray->denseVector();
265  int numberInRowArray = rowArray->getNumElements();
266  // maybe I need one in OsiSimplex
267  double zeroTolerance = model->factorization()->zeroTolerance();
268  int numberRows = model->numberRows();
269  ClpPackedMatrix* rowCopy =
270    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
271  bool packed = rowArray->packedMode();
272  double factor = 0.3;
273  // We may not want to do by row if there may be cache problems
274  int numberColumns = model->numberColumns();
275  // It would be nice to find L2 cache size - for moment 512K
276  // Be slightly optimistic
277  if (numberColumns*sizeof(double)>1000000) {
278    if (numberRows*10<numberColumns)
279      factor=0.1;
280    else if (numberRows*4<numberColumns)
281      factor=0.15;
282    else if (numberRows*2<numberColumns)
283      factor=0.2;
284    //if (model->numberIterations()%50==0)
285    //printf("%d nonzero\n",numberInRowArray);
286  }
287  if (numberInRowArray>factor*numberRows||!rowCopy) {
288    // do by column
289    int iColumn;
290    // get matrix data pointers
291    const int * row = matrix_->getIndices();
292    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
293    const int * columnLength = matrix_->getVectorLengths(); 
294    const double * elementByColumn = matrix_->getElements();
295    const double * rowScale = model->rowScale();
296    int numberColumns = model->numberColumns();
297    if (!y->getNumElements()) {
298      if (packed) {
299        // need to expand pi into y
300        assert(y->capacity()>=numberRows);
301        double * piOld = pi;
302        pi = y->denseVector();
303        const int * whichRow = rowArray->getIndices();
304        int i;
305        if (!rowScale) {
306          // modify pi so can collapse to one loop
307          for (i=0;i<numberInRowArray;i++) {
308            int iRow = whichRow[i];
309            pi[iRow]=scalar*piOld[i];
310          }
311          for (iColumn=0;iColumn<numberColumns;iColumn++) {
312            double value = 0.0;
313            CoinBigIndex j;
314            for (j=columnStart[iColumn];
315                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
316              int iRow = row[j];
317              value += pi[iRow]*elementByColumn[j];
318            }
319            if (fabs(value)>zeroTolerance) {
320              array[numberNonZero]=value;
321              index[numberNonZero++]=iColumn;
322            }
323          }
324        } else {
325          // scaled
326          // modify pi so can collapse to one loop
327          for (i=0;i<numberInRowArray;i++) {
328            int iRow = whichRow[i];
329            pi[iRow]=scalar*piOld[i]*rowScale[iRow];
330          }
331          for (iColumn=0;iColumn<numberColumns;iColumn++) {
332            double value = 0.0;
333            CoinBigIndex j;
334            const double * columnScale = model->columnScale();
335            for (j=columnStart[iColumn];
336                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
337              int iRow = row[j];
338              value += pi[iRow]*elementByColumn[j];
339            }
340            value *= columnScale[iColumn];
341            if (fabs(value)>zeroTolerance) {
342              array[numberNonZero]=value;
343              index[numberNonZero++]=iColumn;
344            }
345          }
346        }
347        // zero out
348        for (i=0;i<numberInRowArray;i++) {
349          int iRow = whichRow[i];
350          pi[iRow]=0.0;
351        }
352      } else {
353        if (!rowScale) {
354          if (scalar==-1.0) {
355            for (iColumn=0;iColumn<numberColumns;iColumn++) {
356              double value = 0.0;
357              CoinBigIndex j;
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              if (fabs(value)>zeroTolerance) {
364                index[numberNonZero++]=iColumn;
365                array[iColumn]=-value;
366              }
367            }
368          } else if (scalar==1.0) {
369            for (iColumn=0;iColumn<numberColumns;iColumn++) {
370              double value = 0.0;
371              CoinBigIndex j;
372              for (j=columnStart[iColumn];
373                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
374                int iRow = row[j];
375                value += pi[iRow]*elementByColumn[j];
376              }
377              if (fabs(value)>zeroTolerance) {
378                index[numberNonZero++]=iColumn;
379                array[iColumn]=value;
380              }
381            }
382          } else {
383            for (iColumn=0;iColumn<numberColumns;iColumn++) {
384              double value = 0.0;
385              CoinBigIndex j;
386              for (j=columnStart[iColumn];
387                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
388                int iRow = row[j];
389                value += pi[iRow]*elementByColumn[j];
390              }
391              value *= scalar;
392              if (fabs(value)>zeroTolerance) {
393                index[numberNonZero++]=iColumn;
394                array[iColumn]=value;
395              }
396            }
397          }
398        } else {
399          // scaled
400          if (scalar==-1.0) {
401            for (iColumn=0;iColumn<numberColumns;iColumn++) {
402              double value = 0.0;
403              CoinBigIndex j;
404              const double * columnScale = model->columnScale();
405              for (j=columnStart[iColumn];
406                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
407                int iRow = row[j];
408                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
409              }
410              value *= columnScale[iColumn];
411              if (fabs(value)>zeroTolerance) {
412                index[numberNonZero++]=iColumn;
413                array[iColumn]=-value;
414              }
415            }
416          } else if (scalar==1.0) {
417            for (iColumn=0;iColumn<numberColumns;iColumn++) {
418              double value = 0.0;
419              CoinBigIndex j;
420              const double * columnScale = model->columnScale();
421              for (j=columnStart[iColumn];
422                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
423                int iRow = row[j];
424                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
425              }
426              value *= columnScale[iColumn];
427              if (fabs(value)>zeroTolerance) {
428                index[numberNonZero++]=iColumn;
429                array[iColumn]=value;
430              }
431            }
432          } else {
433            for (iColumn=0;iColumn<numberColumns;iColumn++) {
434              double value = 0.0;
435              CoinBigIndex j;
436              const double * columnScale = model->columnScale();
437              for (j=columnStart[iColumn];
438                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
439                int iRow = row[j];
440                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
441              }
442              value *= scalar*columnScale[iColumn];
443              if (fabs(value)>zeroTolerance) {
444                index[numberNonZero++]=iColumn;
445                array[iColumn]=value;
446              }
447            }
448          }
449        }
450      }
451    } else {
452      assert(!packed);
453      double * markVector = y->denseVector(); // not empty
454      if (!rowScale) {
455        for (iColumn=0;iColumn<numberColumns;iColumn++) {
456          double value = markVector[iColumn];
457          markVector[iColumn]=0.0;
458          double value2 = 0.0;
459          CoinBigIndex j;
460          for (j=columnStart[iColumn];
461               j<columnStart[iColumn]+columnLength[iColumn];j++) {
462            int iRow = row[j];
463            value2 += pi[iRow]*elementByColumn[j];
464          }
465          value += value2*scalar;
466          if (fabs(value)>zeroTolerance) {
467            index[numberNonZero++]=iColumn;
468            array[iColumn]=value;
469          }
470        }
471      } else {
472        // scaled
473        for (iColumn=0;iColumn<numberColumns;iColumn++) {
474          double value = markVector[iColumn];
475          markVector[iColumn]=0.0;
476          CoinBigIndex j;
477          const double * columnScale = model->columnScale();
478          double value2 = 0.0;
479          for (j=columnStart[iColumn];
480               j<columnStart[iColumn]+columnLength[iColumn];j++) {
481            int iRow = row[j];
482            value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
483          }
484          value += value2*scalar*columnScale[iColumn];
485          if (fabs(value)>zeroTolerance) {
486            index[numberNonZero++]=iColumn;
487            array[iColumn]=value;
488          }
489        }
490      }
491    }
492    columnArray->setNumElements(numberNonZero);
493    y->setNumElements(0);
494  } else {
495    // do by row
496    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
497  }
498  if (packed)
499    columnArray->setPackedMode(true);
500  if (0) {
501    columnArray->checkClean();
502    int numberNonZero=columnArray->getNumElements();;
503    int * index = columnArray->getIndices();
504    double * array = columnArray->denseVector();
505    int i;
506    for (i=0;i<numberNonZero;i++) {
507      int j=index[i];
508      double value;
509      if (packed)
510        value=array[i];
511      else
512        value=array[j];
513      printf("Ti %d %d %g\n",i,j,value);
514    }
515  }
516}
517/* Return <code>x * A + y</code> in <code>z</code>.
518        Squashes small elements and knows about ClpSimplex */
519void 
520ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
521                              const CoinIndexedVector * rowArray,
522                              CoinIndexedVector * y,
523                              CoinIndexedVector * columnArray) const
524{
525  columnArray->clear();
526  double * pi = rowArray->denseVector();
527  int numberNonZero=0;
528  int * index = columnArray->getIndices();
529  double * array = columnArray->denseVector();
530  int numberInRowArray = rowArray->getNumElements();
531  // maybe I need one in OsiSimplex
532  double zeroTolerance = model->factorization()->zeroTolerance();
533  const int * column = getIndices();
534  const CoinBigIndex * rowStart = getVectorStarts();
535  const double * element = getElements();
536  const int * whichRow = rowArray->getIndices();
537  bool packed = rowArray->packedMode();
538  if (numberInRowArray>2||y->getNumElements()) {
539    // do by rows
540    // ** Row copy is already scaled
541    int iRow;
542    int * mark = y->getIndices();
543    int numberOriginal=y->getNumElements();
544    int i;
545    if (packed) {
546      assert(!numberOriginal);
547      numberNonZero=0;
548      // and set up mark as char array
549      char * marked = (char *) (index+columnArray->capacity());
550      double * array2 = y->denseVector();
551#ifdef CLP_DEBUG
552      int numberColumns = model->numberColumns();
553      for (i=0;i<numberColumns;i++) {
554        assert(!marked[i]);
555        assert(!array2[i]);
556      }
557#endif     
558      for (i=0;i<numberInRowArray;i++) {
559        iRow = whichRow[i]; 
560        double value = pi[i]*scalar;
561        CoinBigIndex j;
562        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
563          int iColumn = column[j];
564          if (!marked[iColumn]) {
565            marked[iColumn]=1;
566            index[numberNonZero++]=iColumn;
567          }
568          array2[iColumn] += value*element[j];
569        }
570      }
571      // get rid of tiny values and zero out marked
572      numberOriginal=numberNonZero;
573      numberNonZero=0;
574      for (i=0;i<numberOriginal;i++) {
575        int iColumn = index[i];
576        if (marked[iColumn]) {
577          double value = array2[iColumn];
578          array2[iColumn]=0.0;
579          marked[iColumn]=0;
580          if (fabs(value)>zeroTolerance) {
581            array[numberNonZero]=value;
582            index[numberNonZero++]=iColumn;
583          }
584        }
585      }
586    } else {
587      double * markVector = y->denseVector(); // probably empty .. but
588      for (i=0;i<numberOriginal;i++) {
589        int iColumn = mark[i];
590        index[i]=iColumn;
591        array[iColumn]=markVector[iColumn];
592        markVector[iColumn]=0.0;
593      }
594      numberNonZero=numberOriginal;
595      // and set up mark as char array
596      char * marked = (char *) markVector;
597      for (i=0;i<numberOriginal;i++) {
598        int iColumn = index[i];
599        marked[iColumn]=0;
600      }
601     
602      for (i=0;i<numberInRowArray;i++) {
603        iRow = whichRow[i]; 
604        double value = pi[iRow]*scalar;
605        CoinBigIndex j;
606        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
607          int iColumn = column[j];
608          if (!marked[iColumn]) {
609            marked[iColumn]=1;
610            index[numberNonZero++]=iColumn;
611          }
612          array[iColumn] += value*element[j];
613        }
614      }
615      // get rid of tiny values and zero out marked
616      numberOriginal=numberNonZero;
617      numberNonZero=0;
618      for (i=0;i<numberOriginal;i++) {
619        int iColumn = index[i];
620        marked[iColumn]=0;
621        if (fabs(array[iColumn])>zeroTolerance) {
622          index[numberNonZero++]=iColumn;
623        } else {
624          array[iColumn]=0.0;
625        }
626      }
627    }
628  } else if (numberInRowArray==2) {
629    // do by rows when two rows
630    int numberOriginal;
631    int i;
632    CoinBigIndex j;
633    numberNonZero=0;
634
635    double value;
636    if (packed) {
637      int iRow0 = whichRow[0]; 
638      int iRow1 = whichRow[1]; 
639      double pi0 = pi[0];
640      double pi1 = pi[1];
641      if (rowStart[iRow0+1]-rowStart[iRow0]>
642          rowStart[iRow1+1]-rowStart[iRow1]) {
643        // do one with fewer first
644        iRow0=iRow1;
645        iRow1=whichRow[0];
646        pi0=pi1;
647        pi1=pi[0];
648      }
649      // and set up mark as char array
650      char * marked = (char *) (index+columnArray->capacity());
651      int * lookup = y->getIndices();
652      value = pi0*scalar;
653      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
654        int iColumn = column[j];
655        double value2 = value*element[j];
656        array[numberNonZero] = value2;
657        marked[iColumn]=1;
658        lookup[iColumn]=numberNonZero;
659        index[numberNonZero++]=iColumn;
660      }
661      numberOriginal = numberNonZero;
662      value = pi1*scalar;
663      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
664        int iColumn = column[j];
665        double value2 = value*element[j];
666        // I am assuming no zeros in matrix
667        if (marked[iColumn]) {
668          int iLookup = lookup[iColumn];
669          array[iLookup] += value2;
670        } else {
671          if (fabs(value2)>zeroTolerance) {
672            array[numberNonZero] = value2;
673            index[numberNonZero++]=iColumn;
674          }
675        }
676      }
677      // get rid of tiny values and zero out marked
678      int nDelete=0;
679      for (i=0;i<numberOriginal;i++) {
680        int iColumn = index[i];
681        marked[iColumn]=0;
682        if (fabs(array[i])<=zeroTolerance) 
683          nDelete++;
684      }
685      if (nDelete) {
686        numberOriginal=numberNonZero;
687        numberNonZero=0;
688        for (i=0;i<numberOriginal;i++) {
689          int iColumn = index[i];
690          double value = array[i];
691          array[i]=0.0;
692          if (fabs(value)>zeroTolerance) {
693            array[numberNonZero]=value;
694            index[numberNonZero++]=iColumn;
695          }
696        }
697      }
698    } else {
699      int iRow = whichRow[0]; 
700      value = pi[iRow]*scalar;
701      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
702        int iColumn = column[j];
703        double value2 = value*element[j];
704        index[numberNonZero++]=iColumn;
705        array[iColumn] = value2;
706      }
707      iRow = whichRow[1]; 
708      value = pi[iRow]*scalar;
709      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
710        int iColumn = column[j];
711        double value2 = value*element[j];
712        // I am assuming no zeros in matrix
713        if (array[iColumn])
714          value2 += array[iColumn];
715        else
716          index[numberNonZero++]=iColumn;
717        array[iColumn] = value2;
718      }
719      // get rid of tiny values and zero out marked
720      numberOriginal=numberNonZero;
721      numberNonZero=0;
722      for (i=0;i<numberOriginal;i++) {
723        int iColumn = index[i];
724        if (fabs(array[iColumn])>zeroTolerance) {
725          index[numberNonZero++]=iColumn;
726        } else {
727          array[iColumn]=0.0;
728        }
729      }
730    }
731  } else if (numberInRowArray==1) {
732    // Just one row
733    int iRow=rowArray->getIndices()[0];
734    numberNonZero=0;
735    CoinBigIndex j;
736    if (packed) {
737      double value = pi[0]*scalar;
738      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
739        int iColumn = column[j];
740        double value2 = value*element[j];
741        if (fabs(value2)>zeroTolerance) {
742          array[numberNonZero] = value2;
743          index[numberNonZero++]=iColumn;
744        }
745      }
746    } else {
747      double value = pi[iRow]*scalar;
748      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
749        int iColumn = column[j];
750        double value2 = value*element[j];
751        if (fabs(value2)>zeroTolerance) {
752          index[numberNonZero++]=iColumn;
753          array[iColumn] = value2;
754        }
755      }
756    }
757  }
758  columnArray->setNumElements(numberNonZero);
759  y->setNumElements(0);
760}
761/* Return <code>x *A in <code>z</code> but
762   just for indices in y.
763   Squashes small elements and knows about ClpSimplex */
764void 
765ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
766                              const CoinIndexedVector * rowArray,
767                              const CoinIndexedVector * y,
768                              CoinIndexedVector * columnArray) const
769{
770  columnArray->clear();
771  double * pi = rowArray->denseVector();
772  int numberNonZero=0;
773  int * index = columnArray->getIndices();
774  double * array = columnArray->denseVector();
775  // maybe I need one in OsiSimplex
776  double zeroTolerance = model->factorization()->zeroTolerance();
777  int jColumn;
778  // get matrix data pointers
779  const int * row = matrix_->getIndices();
780  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
781  const int * columnLength = matrix_->getVectorLengths(); 
782  const double * elementByColumn = matrix_->getElements();
783  const double * rowScale = model->rowScale();
784  int numberToDo = y->getNumElements();
785  const int * which = y->getIndices();
786  bool packed = rowArray->packedMode();
787  if (packed) {
788    // need to expand pi into y
789    int numberInRowArray = rowArray->getNumElements();
790    assert(y->capacity()>=model->numberRows());
791    double * piOld = pi;
792    pi = y->denseVector();
793    const int * whichRow = rowArray->getIndices();
794    int i;
795    // Do NOT squash small elements - must line up with  y
796    if (!rowScale) {
797      for (i=0;i<numberInRowArray;i++) {
798        int iRow = whichRow[i];
799        pi[iRow]=piOld[i];
800      }
801      for (jColumn=0;jColumn<numberToDo;jColumn++) {
802        int iColumn = which[jColumn];
803        double value = 0.0;
804        CoinBigIndex j;
805        for (j=columnStart[iColumn];
806             j<columnStart[iColumn]+columnLength[iColumn];j++) {
807          int iRow = row[j];
808          value += pi[iRow]*elementByColumn[j];
809        }
810        array[jColumn]=value;
811      }
812    } else {
813      // scaled
814      for (i=0;i<numberInRowArray;i++) {
815        int iRow = whichRow[i];
816        pi[iRow]=rowScale[iRow]*piOld[i];
817      }
818      for (jColumn=0;jColumn<numberToDo;jColumn++) {
819        int iColumn = which[jColumn];
820        double value = 0.0;
821        CoinBigIndex j;
822        const double * columnScale = model->columnScale();
823        for (j=columnStart[iColumn];
824             j<columnStart[iColumn]+columnLength[iColumn];j++) {
825          int iRow = row[j];
826          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
827        }
828        value *= columnScale[iColumn];
829        array[jColumn]=value;
830      }
831    }
832    for (i=0;i<numberInRowArray;i++) {
833      int iRow = whichRow[i];
834      pi[iRow]=0.0;
835    }
836  } else {
837    if (!rowScale) {
838      for (jColumn=0;jColumn<numberToDo;jColumn++) {
839        int iColumn = which[jColumn];
840        double value = 0.0;
841        CoinBigIndex j;
842        for (j=columnStart[iColumn];
843             j<columnStart[iColumn]+columnLength[iColumn];j++) {
844          int iRow = row[j];
845          value += pi[iRow]*elementByColumn[j];
846        }
847        if (fabs(value)>zeroTolerance) {
848          index[numberNonZero++]=iColumn;
849          array[iColumn]=value;
850        }
851      }
852    } else {
853      // scaled
854      for (jColumn=0;jColumn<numberToDo;jColumn++) {
855        int iColumn = which[jColumn];
856        double value = 0.0;
857        CoinBigIndex j;
858        const double * columnScale = model->columnScale();
859        for (j=columnStart[iColumn];
860             j<columnStart[iColumn]+columnLength[iColumn];j++) {
861          int iRow = row[j];
862          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
863        }
864        value *= columnScale[iColumn];
865        if (fabs(value)>zeroTolerance) {
866          index[numberNonZero++]=iColumn;
867          array[iColumn]=value;
868        }
869      }
870    }
871  }
872}
873/* Returns number of elements in basis
874   column is basic if entry >=0 */
875CoinBigIndex
876ClpPackedMatrix::numberInBasis(const int * columnIsBasic) const 
877{
878  int i;
879  int numberColumns = getNumCols();
880  const int * columnLength = matrix_->getVectorLengths(); 
881  CoinBigIndex numberElements=0;
882  if (!zeroElements_) {
883    for (i=0;i<numberColumns;i++) {
884      if (columnIsBasic[i]>=0) {
885        numberElements += columnLength[i];
886      }
887    }
888  } else {
889    // there are zero elements so need to look more closely
890    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
891    const double * elementByColumn = matrix_->getElements();
892    for (i=0;i<numberColumns;i++) {
893      if (columnIsBasic[i]>=0) {
894        CoinBigIndex j;
895        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
896          if (elementByColumn[j])
897            numberElements++;
898        }
899      }
900    }
901  }
902  return numberElements;
903}
904// Fills in basis (Returns number of elements and updates numberBasic)
905CoinBigIndex
906ClpPackedMatrix::fillBasis(const ClpSimplex * model,
907                                const int * columnIsBasic, int & numberBasic,
908                                int * indexRowU, int * indexColumnU,
909                                double * elementU) const 
910{
911  const double * rowScale = model->rowScale();
912  const int * row = matrix_->getIndices();
913  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
914  const int * columnLength = matrix_->getVectorLengths(); 
915  const double * elementByColumn = matrix_->getElements();
916  int i;
917  int numberColumns = getNumCols();
918  CoinBigIndex numberElements=0;
919  if (!zeroElements_) {
920    if (!rowScale) {
921      // no scaling
922      for (i=0;i<numberColumns;i++) {
923        if (columnIsBasic[i]>=0) {
924          CoinBigIndex j;
925          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
926            indexRowU[numberElements]=row[j];
927            indexColumnU[numberElements]=numberBasic;
928            elementU[numberElements++]=elementByColumn[j];
929          }
930          numberBasic++;
931        }
932      }
933    } else {
934      // scaling
935      const double * columnScale = model->columnScale();
936      for (i=0;i<numberColumns;i++) {
937        if (columnIsBasic[i]>=0) {
938          CoinBigIndex j;
939          double scale = columnScale[i];
940          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
941            int iRow = row[j];
942            indexRowU[numberElements]=iRow;
943            indexColumnU[numberElements]=numberBasic;
944            elementU[numberElements++]=elementByColumn[j]*scale*rowScale[iRow];
945          }
946          numberBasic++;
947        }
948      }
949    }
950  } else {
951    // there are zero elements so need to look more closely
952    if (!rowScale) {
953      // no scaling
954      for (i=0;i<numberColumns;i++) {
955        if (columnIsBasic[i]>=0) {
956          CoinBigIndex j;
957          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
958            double value = elementByColumn[j];
959            if (value) {
960              indexRowU[numberElements]=row[j];
961              indexColumnU[numberElements]=numberBasic;
962              elementU[numberElements++]=value;
963            }
964          }
965          numberBasic++;
966        }
967      }
968    } else {
969      // scaling
970      const double * columnScale = model->columnScale();
971      for (i=0;i<numberColumns;i++) {
972        if (columnIsBasic[i]>=0) {
973          CoinBigIndex j;
974          double scale = columnScale[i];
975          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
976            double value = elementByColumn[j];
977            if (value) {
978              int iRow = row[j];
979              indexRowU[numberElements]=iRow;
980              indexColumnU[numberElements]=numberBasic;
981              elementU[numberElements++]=value*scale*rowScale[iRow];
982            }
983          }
984          numberBasic++;
985        }
986      }
987    }
988  }
989  return numberElements;
990}
991/* If element NULL returns number of elements in column part of basis,
992   If not NULL fills in as well */
993CoinBigIndex
994ClpPackedMatrix::fillBasis(const ClpSimplex * model,
995                           const int * whichColumn, 
996                           int numberBasic,
997                           int numberColumnBasic,
998                           int * indexRowU, int * indexColumnU,
999                           double * elementU) const 
1000{
1001  const int * columnLength = matrix_->getVectorLengths(); 
1002  int i;
1003  CoinBigIndex numberElements=0;
1004  if (elementU!=NULL) {
1005    // fill
1006    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1007    const double * rowScale = model->rowScale();
1008    const int * row = matrix_->getIndices();
1009    const double * elementByColumn = matrix_->getElements();
1010    if (!zeroElements_) {
1011      if (!rowScale) {
1012        // no scaling
1013        for (i=0;i<numberColumnBasic;i++) {
1014          int iColumn = whichColumn[i];
1015          CoinBigIndex j;
1016          for (j=columnStart[iColumn];
1017               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1018            indexRowU[numberElements]=row[j];
1019            indexColumnU[numberElements]=numberBasic;
1020            elementU[numberElements++]=elementByColumn[j];
1021          }
1022          numberBasic++;
1023        }
1024      } else {
1025        // scaling
1026        const double * columnScale = model->columnScale();
1027        for (i=0;i<numberColumnBasic;i++) {
1028          int iColumn = whichColumn[i];
1029          CoinBigIndex j;
1030          double scale = columnScale[iColumn];
1031          for (j=columnStart[iColumn];
1032               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1033            int iRow = row[j];
1034            indexRowU[numberElements]=iRow;
1035            indexColumnU[numberElements]=numberBasic;
1036            elementU[numberElements++]=
1037              elementByColumn[j]*scale*rowScale[iRow];
1038          }
1039          numberBasic++;
1040        }
1041      }
1042    } else {
1043      // there are zero elements so need to look more closely
1044      if (!rowScale) {
1045        // no scaling
1046        for (i=0;i<numberColumnBasic;i++) {
1047          int iColumn = whichColumn[i];
1048          CoinBigIndex j;
1049          for (j=columnStart[iColumn];
1050               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1051            double value = elementByColumn[j];
1052            if (value) {
1053              indexRowU[numberElements]=row[j];
1054              indexColumnU[numberElements]=numberBasic;
1055              elementU[numberElements++]=value;
1056            }
1057          }
1058          numberBasic++;
1059        }
1060      } else {
1061        // scaling
1062        const double * columnScale = model->columnScale();
1063        for (i=0;i<numberColumnBasic;i++) {
1064          int iColumn = whichColumn[i];
1065          CoinBigIndex j;
1066          double scale = columnScale[iColumn];
1067          for (j=columnStart[iColumn];
1068               j<columnStart[iColumn]+columnLength[i];j++) {
1069            double value = elementByColumn[j];
1070            if (value) {
1071              int iRow = row[j];
1072              indexRowU[numberElements]=iRow;
1073              indexColumnU[numberElements]=numberBasic;
1074              elementU[numberElements++]=value*scale*rowScale[iRow];
1075            }
1076          }
1077          numberBasic++;
1078        }
1079      }
1080    }
1081  } else {
1082    // just count - can be over so ignore zero problem
1083    for (i=0;i<numberColumnBasic;i++) {
1084      int iColumn = whichColumn[i];
1085      numberElements += columnLength[iColumn];
1086    }
1087  }
1088  return numberElements;
1089}
1090// Creates scales for column copy (rowCopy in model may be modified)
1091int 
1092ClpPackedMatrix::scale(ClpSimplex * model) const 
1093{
1094  ClpMatrixBase * rowCopyBase=model->rowCopy();
1095  if (!rowCopyBase) {
1096    // temporary copy
1097    rowCopyBase = reverseOrderedCopy();
1098  }
1099  ClpPackedMatrix* rowCopy =
1100    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
1101
1102  // Make sure it is really a ClpPackedMatrix
1103  assert (rowCopy!=NULL);
1104  int numberRows = model->numberRows();
1105  int numberColumns = model->numberColumns();
1106  const int * column = rowCopy->getIndices();
1107  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1108  const double * element = rowCopy->getElements();
1109  double * rowScale = new double [numberRows];
1110  double * columnScale = new double [numberColumns];
1111  // we are going to mark bits we are interested in
1112  char * usefulRow = new char [numberRows];
1113  char * usefulColumn = new char [numberColumns];
1114  double * rowLower = model->rowLower();
1115  double * rowUpper = model->rowUpper();
1116  double * columnLower = model->columnLower();
1117  double * columnUpper = model->columnUpper();
1118  int iColumn, iRow;
1119  // mark free rows
1120  for (iRow=0;iRow<numberRows;iRow++) {
1121    usefulRow[iRow]=0;
1122    if (rowUpper[iRow]<1.0e20||
1123        rowLower[iRow]>-1.0e20)
1124      usefulRow[iRow]=1;
1125  }
1126  // mark empty and fixed columns
1127  // also see if worth scaling
1128  assert (model->scalingFlag()<4); // dynamic not implemented
1129  double largest=0.0;
1130  double smallest=1.0e50;
1131  // get matrix data pointers
1132  const int * row = matrix_->getIndices();
1133  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1134  const int * columnLength = matrix_->getVectorLengths(); 
1135  const double * elementByColumn = matrix_->getElements();
1136  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1137    CoinBigIndex j;
1138    char useful=0;
1139    if (columnUpper[iColumn]>
1140        columnLower[iColumn]+1.0e-9) {
1141      for (j=columnStart[iColumn];
1142           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1143        iRow=row[j];
1144        if(elementByColumn[j]&&usefulRow[iRow]) {
1145          useful=1;
1146          largest = max(largest,fabs(elementByColumn[j]));
1147          smallest = min(smallest,fabs(elementByColumn[j]));
1148        }
1149      }
1150    }
1151    usefulColumn[iColumn]=useful;
1152  }
1153  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
1154    <<smallest<<largest
1155    <<CoinMessageEol;
1156  if (smallest>=0.5&&largest<=2.0) {
1157    // don't bother scaling
1158    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
1159      <<CoinMessageEol;
1160    delete [] rowScale;
1161    delete [] usefulRow;
1162    delete [] columnScale;
1163    delete [] usefulColumn;
1164    return 1;
1165  } else {
1166      // need to scale
1167    int scalingMethod = model->scalingFlag();
1168    if (scalingMethod==3) {
1169      // Choose between 1 and 2
1170      if (smallest<1.0e-5||smallest*largest<1.0)
1171        scalingMethod=1;
1172      else
1173        scalingMethod=2;
1174    }
1175    // and see if there any empty rows
1176    for (iRow=0;iRow<numberRows;iRow++) {
1177      if (usefulRow[iRow]) {
1178        CoinBigIndex j;
1179        int useful=0;
1180        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1181          int iColumn = column[j];
1182          if (usefulColumn[iColumn]) {
1183            useful=1;
1184            break;
1185          }
1186        }
1187        usefulRow[iRow]=useful;
1188      }
1189    }
1190    ClpFillN ( rowScale, numberRows,1.0);
1191    ClpFillN ( columnScale, numberColumns,1.0);
1192    double overallLargest=-1.0e-30;
1193    double overallSmallest=1.0e30;
1194    if (scalingMethod==1) {
1195      // Maximum in each row
1196      for (iRow=0;iRow<numberRows;iRow++) {
1197        if (usefulRow[iRow]) {
1198          CoinBigIndex j;
1199          largest=1.0e-20;
1200          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1201            int iColumn = column[j];
1202            if (usefulColumn[iColumn]) {
1203              double value = fabs(element[j]);
1204              largest = max(largest,value);
1205            }
1206          }
1207          rowScale[iRow]=1.0/largest;
1208          overallLargest = max(overallLargest,largest);
1209          overallSmallest = min(overallSmallest,largest);
1210        }
1211      }
1212    } else {
1213      assert(scalingMethod==2);
1214      int numberPass=3;
1215#ifdef USE_OBJECTIVE
1216      // This will be used to help get scale factors
1217      double * objective = new double[numberColumns];
1218      memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
1219      double objScale=1.0;
1220#endif
1221      while (numberPass) {
1222        overallLargest=0.0;
1223        overallSmallest=1.0e50;
1224        numberPass--;
1225        // Geometric mean on row scales
1226        for (iRow=0;iRow<numberRows;iRow++) {
1227          if (usefulRow[iRow]) {
1228            CoinBigIndex j;
1229            largest=1.0e-20;
1230            smallest=1.0e50;
1231            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1232              int iColumn = column[j];
1233              if (usefulColumn[iColumn]) {
1234                double value = fabs(element[j]);
1235                // Don't bother with tiny elements
1236                if (value>1.0e-30) {
1237                  value *= columnScale[iColumn];
1238                  largest = max(largest,value);
1239                  smallest = min(smallest,value);
1240                }
1241              }
1242            }
1243            rowScale[iRow]=1.0/sqrt(smallest*largest);
1244            overallLargest = max(largest*rowScale[iRow],overallLargest);
1245            overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
1246          }
1247        }
1248#ifdef USE_OBJECTIVE
1249        largest=1.0e-20;
1250        smallest=1.0e50;
1251        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1252          if (usefulColumn[iColumn]) {
1253            double value = fabs(objective[iColumn]);
1254            // Don't bother with tiny elements
1255            if (value>1.0e-30) {
1256              value *= columnScale[iColumn];
1257              largest = max(largest,value);
1258              smallest = min(smallest,value);
1259            }
1260          }
1261        }
1262        objScale=1.0/sqrt(smallest*largest);
1263#endif
1264        model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
1265          <<overallSmallest
1266          <<overallLargest
1267          <<CoinMessageEol;
1268        // skip last column round
1269        if (numberPass==1)
1270          break;
1271        // Geometric mean on column scales
1272        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1273          if (usefulColumn[iColumn]) {
1274            CoinBigIndex j;
1275            largest=1.0e-20;
1276            smallest=1.0e50;
1277            for (j=columnStart[iColumn];
1278                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1279              iRow=row[j];
1280              double value = fabs(elementByColumn[j]);
1281              // Don't bother with tiny elements
1282              if (value>1.0e-30&&usefulRow[iRow]) {
1283                value *= rowScale[iRow];
1284                largest = max(largest,value);
1285                smallest = min(smallest,value);
1286              }
1287            }
1288#ifdef USE_OBJECTIVE
1289            if (fabs(objective[iColumn])>1.0e-30) {
1290              double value = fabs(objective[iColumn])*objScale;
1291              largest = max(largest,value);
1292              smallest = min(smallest,value);
1293            }
1294#endif
1295            columnScale[iColumn]=1.0/sqrt(smallest*largest);
1296          }
1297        }
1298      }
1299#ifdef USE_OBJECTIVE
1300      delete [] objective;
1301      printf("obj scale %g - use it if you want\n",objScale);
1302#endif
1303    }
1304    // If ranges will make horrid then scale
1305    double tolerance = 5.0*model->primalTolerance();
1306    for (iRow=0;iRow<numberRows;iRow++) {
1307      if (usefulRow[iRow]) {
1308        double difference = rowUpper[iRow]-rowLower[iRow];
1309        double scaledDifference = difference*rowScale[iRow];
1310        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1311          // make gap larger
1312          rowScale[iRow] *= 1.0e-4/scaledDifference;
1313          //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
1314          // scaledDifference,difference*rowScale[iRow]);
1315        }
1316      }
1317    }
1318    // final pass to scale columns so largest is reasonable
1319    // See what smallest will be if largest is 1.0
1320    overallSmallest=1.0e50;
1321    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1322      if (usefulColumn[iColumn]) {
1323        CoinBigIndex j;
1324        largest=1.0e-20;
1325        smallest=1.0e50;
1326        for (j=columnStart[iColumn];
1327             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1328          iRow=row[j];
1329          if(elementByColumn[j]&&usefulRow[iRow]) {
1330            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1331            largest = max(largest,value);
1332            smallest = min(smallest,value);
1333          }
1334        }
1335        if (overallSmallest*largest>smallest)
1336          overallSmallest = smallest/largest;
1337      }
1338    }
1339#define RANDOMIZE
1340#ifdef RANDOMIZE
1341    // randomize by up to 10%
1342    for (iRow=0;iRow<numberRows;iRow++) {
1343      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1344      rowScale[iRow] *= (1.0+0.1*value);
1345    }
1346#endif
1347    overallLargest=1.0;
1348    if (overallSmallest<1.0e-1)
1349      overallLargest = 1.0/sqrt(overallSmallest);
1350    overallLargest = min(1000.0,overallLargest);
1351    overallSmallest=1.0e50;
1352    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1353      if (usefulColumn[iColumn]) {
1354        CoinBigIndex j;
1355        largest=1.0e-20;
1356        smallest=1.0e50;
1357        for (j=columnStart[iColumn];
1358             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1359          iRow=row[j];
1360          if(elementByColumn[j]&&usefulRow[iRow]) {
1361            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1362            largest = max(largest,value);
1363            smallest = min(smallest,value);
1364          }
1365        }
1366        columnScale[iColumn]=overallLargest/largest;
1367#ifdef RANDOMIZE
1368        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1369        columnScale[iColumn] *= (1.0+0.1*value);
1370#endif
1371        double difference = columnUpper[iColumn]-columnLower[iColumn];
1372        double scaledDifference = difference*columnScale[iColumn];
1373        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1374          // make gap larger
1375          columnScale[iColumn] *= 1.0e-4/scaledDifference;
1376          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
1377          // scaledDifference,difference*columnScale[iColumn]);
1378        }
1379        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
1380      }
1381    }
1382    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
1383      <<overallSmallest
1384      <<overallLargest
1385      <<CoinMessageEol;
1386    delete [] usefulRow;
1387    delete [] usefulColumn;
1388    // If quadratic then make symmetric
1389    ClpObjective * obj = model->objectiveAsObject();
1390    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
1391    if (quadraticObj) {
1392      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1393      int numberXColumns = quadratic->getNumCols();
1394      if (numberXColumns<numberColumns) {
1395        // we assume symmetric
1396        int numberQuadraticColumns=0;
1397        int i;
1398        //const int * columnQuadratic = quadratic->getIndices();
1399        //const int * columnQuadraticStart = quadratic->getVectorStarts();
1400        const int * columnQuadraticLength = quadratic->getVectorLengths();
1401        for (i=0;i<numberXColumns;i++) {
1402          int length=columnQuadraticLength[i];
1403#ifndef CORRECT_COLUMN_COUNTS
1404          length=1;
1405#endif
1406          if (length)
1407            numberQuadraticColumns++;
1408        }
1409        int numberXRows = numberRows-numberQuadraticColumns;
1410        numberQuadraticColumns=0;
1411        for (i=0;i<numberXColumns;i++) { 
1412          int length=columnQuadraticLength[i];
1413#ifndef CORRECT_COLUMN_COUNTS
1414          length=1;
1415#endif
1416          if (length) {
1417            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
1418            numberQuadraticColumns++;
1419          }
1420        }   
1421        int numberQuadraticRows=0;
1422        for (i=0;i<numberXRows;i++) {
1423          // See if any in row quadratic
1424          int j;
1425          int numberQ=0;
1426          for (j=rowStart[i];j<rowStart[i+1];j++) {
1427            int iColumn = column[j];
1428            if (columnQuadraticLength[iColumn])
1429              numberQ++;
1430          }
1431#ifndef CORRECT_ROW_COUNTS
1432          numberQ=1;
1433#endif
1434          if (numberQ) {
1435            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
1436            numberQuadraticRows++;
1437          }
1438        }
1439        // and make sure Sj okay
1440        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
1441          CoinBigIndex j=columnStart[iColumn];
1442          assert(columnLength[iColumn]==1);
1443          int iRow=row[j];
1444          double value = fabs(elementByColumn[j]*rowScale[iRow]);
1445          columnScale[iColumn]=1.0/value;
1446        }
1447      }
1448    }
1449    model->setRowScale(rowScale);
1450    model->setColumnScale(columnScale);
1451    if (model->rowCopy()) {
1452      // need to replace row by row
1453      double * newElement = new double[numberColumns];
1454      // scale row copy
1455      for (iRow=0;iRow<numberRows;iRow++) {
1456        int j;
1457        double scale = rowScale[iRow];
1458        const double * elementsInThisRow = element + rowStart[iRow];
1459        const int * columnsInThisRow = column + rowStart[iRow];
1460        int number = rowStart[iRow+1]-rowStart[iRow];
1461        assert (number<=numberColumns);
1462        for (j=0;j<number;j++) {
1463          int iColumn = columnsInThisRow[j];
1464          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
1465        }
1466        rowCopy->replaceVector(iRow,number,newElement);
1467      }
1468      delete [] newElement;
1469    } else {
1470      // no row copy
1471      delete rowCopyBase;
1472    }
1473    return 0;
1474  }
1475}
1476/* Unpacks a column into an CoinIndexedvector
1477 */
1478void 
1479ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
1480                   int iColumn) const 
1481{
1482  const double * rowScale = model->rowScale();
1483  const int * row = matrix_->getIndices();
1484  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1485  const int * columnLength = matrix_->getVectorLengths(); 
1486  const double * elementByColumn = matrix_->getElements();
1487  CoinBigIndex i;
1488  if (!rowScale) {
1489    for (i=columnStart[iColumn];
1490         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1491      rowArray->add(row[i],elementByColumn[i]);
1492    }
1493  } else {
1494    // apply scaling
1495    double scale = model->columnScale()[iColumn];
1496    for (i=columnStart[iColumn];
1497         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1498      int iRow = row[i];
1499      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1500    }
1501  }
1502}
1503/* Unpacks a column into a CoinIndexedVector
1504** in packed format
1505Note that model is NOT const.  Bounds and objective could
1506be modified if doing column generation (just for this variable) */
1507void 
1508ClpPackedMatrix::unpackPacked(ClpSimplex * model,
1509                            CoinIndexedVector * rowArray,
1510                            int iColumn) const
1511{
1512  const double * rowScale = model->rowScale();
1513  const int * row = matrix_->getIndices();
1514  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1515  const int * columnLength = matrix_->getVectorLengths(); 
1516  const double * elementByColumn = matrix_->getElements();
1517  CoinBigIndex i;
1518  if (!rowScale) {
1519    int j=columnStart[iColumn];
1520    rowArray->createPacked(columnLength[iColumn],
1521                           row+j,elementByColumn+j);
1522  } else {
1523    // apply scaling
1524    double scale = model->columnScale()[iColumn];
1525    int * index = rowArray->getIndices();
1526    double * array = rowArray->denseVector();
1527    int number = 0;
1528    for (i=columnStart[iColumn];
1529         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1530      int iRow = row[i];
1531      array[number]=elementByColumn[i]*scale*rowScale[iRow];
1532      index[number++]=iRow;
1533    }
1534    rowArray->setNumElements(number);
1535    rowArray->setPackedMode(true);
1536  }
1537}
1538/* Adds multiple of a column into an CoinIndexedvector
1539      You can use quickAdd to add to vector */
1540void 
1541ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1542                   int iColumn, double multiplier) const 
1543{
1544  const double * rowScale = model->rowScale();
1545  const int * row = matrix_->getIndices();
1546  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1547  const int * columnLength = matrix_->getVectorLengths(); 
1548  const double * elementByColumn = matrix_->getElements();
1549  CoinBigIndex i;
1550  if (!rowScale) {
1551    for (i=columnStart[iColumn];
1552         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1553      int iRow = row[i];
1554      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
1555    }
1556  } else {
1557    // apply scaling
1558    double scale = model->columnScale()[iColumn]*multiplier;
1559    for (i=columnStart[iColumn];
1560         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1561      int iRow = row[i];
1562      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1563    }
1564  }
1565}
1566/* Checks if all elements are in valid range.  Can just
1567   return true if you are not paranoid.  For Clp I will
1568   probably expect no zeros.  Code can modify matrix to get rid of
1569   small elements.
1570*/
1571bool 
1572ClpPackedMatrix::allElementsInRange(ClpModel * model,
1573                                    double smallest, double largest)
1574{
1575  int iColumn;
1576  CoinBigIndex numberLarge=0;;
1577  CoinBigIndex numberSmall=0;;
1578  CoinBigIndex numberDuplicate=0;;
1579  int firstBadColumn=-1;
1580  int firstBadRow=-1;
1581  double firstBadElement=0.0;
1582  // get matrix data pointers
1583  const int * row = matrix_->getIndices();
1584  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1585  const int * columnLength = matrix_->getVectorLengths(); 
1586  const double * elementByColumn = matrix_->getElements();
1587  int numberColumns = matrix_->getNumCols();
1588  int numberRows = matrix_->getNumRows();
1589  int * mark = new int [numberRows];
1590  int i;
1591  for (i=0;i<numberRows;i++)
1592    mark[i]=-1;
1593  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1594    CoinBigIndex j;
1595    for (j=columnStart[iColumn];
1596         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1597      double value = fabs(elementByColumn[j]);
1598      int iRow = row[j];
1599      if (iRow<0||iRow>=numberRows) {
1600        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1601        return false;
1602      }
1603      if (mark[iRow]==-1) {
1604        mark[iRow]=j;
1605      } else {
1606        // duplicate
1607        numberDuplicate++;
1608      }
1609      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1610      if (!value)
1611        zeroElements_ = true; // there are zero elements
1612      if (value<smallest) {
1613        numberSmall++;
1614      } else if (value>largest) {
1615        numberLarge++;
1616        if (firstBadColumn<0) {
1617          firstBadColumn=iColumn;
1618          firstBadRow=row[j];
1619          firstBadElement=elementByColumn[j];
1620        }
1621      }
1622    }
1623    //clear mark
1624    for (j=columnStart[iColumn];
1625         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1626      int iRow = row[j];
1627      mark[iRow]=-1;
1628    }
1629  }
1630  delete [] mark;
1631  if (numberLarge) {
1632    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
1633      <<numberLarge
1634      <<firstBadColumn<<firstBadRow<<firstBadElement
1635      <<CoinMessageEol;
1636    return false;
1637  }
1638  if (numberSmall) 
1639    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
1640      <<numberSmall
1641      <<CoinMessageEol;
1642  if (numberDuplicate) 
1643    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
1644      <<numberDuplicate
1645      <<CoinMessageEol;
1646  if (numberDuplicate) 
1647    matrix_->eliminateDuplicates(smallest);
1648  else if (numberSmall) 
1649    matrix_->compress(smallest);
1650  // If smallest >0.0 then there can't be zero elements
1651  if (smallest>0.0)
1652    zeroElements_=false;
1653  return true;
1654}
1655/* Given positive integer weights for each row fills in sum of weights
1656   for each column (and slack).
1657   Returns weights vector
1658*/
1659CoinBigIndex * 
1660ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
1661{
1662  int numberRows = model->numberRows();
1663  int numberColumns =model->numberColumns();
1664  int number = numberRows+numberColumns;
1665  CoinBigIndex * weights = new CoinBigIndex[number];
1666  // get matrix data pointers
1667  const int * row = matrix_->getIndices();
1668  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1669  const int * columnLength = matrix_->getVectorLengths(); 
1670  int i;
1671  for (i=0;i<numberColumns;i++) {
1672    CoinBigIndex j;
1673    CoinBigIndex count=0;
1674    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1675      int iRow=row[j];
1676      count += inputWeights[iRow];
1677    }
1678    weights[i]=count;
1679  }
1680  for (i=0;i<numberRows;i++) {
1681    weights[i+numberColumns]=inputWeights[i];
1682  }
1683  return weights;
1684}
1685/* Returns largest and smallest elements of both signs.
1686   Largest refers to largest absolute value.
1687*/
1688void 
1689ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1690                       double & smallestPositive, double & largestPositive)
1691{
1692  smallestNegative=-COIN_DBL_MAX;
1693  largestNegative=0.0;
1694  smallestPositive=COIN_DBL_MAX;
1695  largestPositive=0.0;
1696  int numberColumns = matrix_->getNumCols();
1697  // get matrix data pointers
1698  const double * elementByColumn = matrix_->getElements();
1699  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1700  const int * columnLength = matrix_->getVectorLengths(); 
1701  int i;
1702  for (i=0;i<numberColumns;i++) {
1703    CoinBigIndex j;
1704    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1705      double value = elementByColumn[j];
1706      if (value>0.0) {
1707        smallestPositive = min(smallestPositive,value);
1708        largestPositive = max(largestPositive,value);
1709      } else if (value<0.0) {
1710        smallestNegative = max(smallestNegative,value);
1711        largestNegative = min(largestNegative,value);
1712      }
1713    }
1714  }
1715}
1716
1717
Note: See TracBrowser for help on using the repository browser.