source: branches/pre/ClpPackedMatrix.cpp @ 212

Last change on this file since 212 was 212, checked in by forrest, 17 years ago

lots of stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.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    int scalingMethod = model->scalingFlag();
1167    if (scalingMethod==3) {
1168      // Choose between 1 and 2
1169      if (smallest<1.0e-5||smallest*largest<1.0)
1170        scalingMethod=1;
1171      else
1172        scalingMethod=2;
1173    }
1174    // and see if there any empty rows
1175    for (iRow=0;iRow<numberRows;iRow++) {
1176      if (usefulRow[iRow]) {
1177        CoinBigIndex j;
1178        int useful=0;
1179        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1180          int iColumn = column[j];
1181          if (usefulColumn[iColumn]) {
1182            useful=1;
1183            break;
1184          }
1185        }
1186        usefulRow[iRow]=useful;
1187      }
1188    }
1189    ClpFillN ( rowScale, numberRows,1.0);
1190    ClpFillN ( columnScale, numberColumns,1.0);
1191    double overallLargest=-1.0e-30;
1192    double overallSmallest=1.0e30;
1193    if (scalingMethod==1) {
1194      // Maximum in each row
1195      for (iRow=0;iRow<numberRows;iRow++) {
1196        if (usefulRow[iRow]) {
1197          CoinBigIndex j;
1198          largest=1.0e-20;
1199          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1200            int iColumn = column[j];
1201            if (usefulColumn[iColumn]) {
1202              double value = fabs(element[j]);
1203              largest = max(largest,value);
1204            }
1205          }
1206          rowScale[iRow]=1.0/largest;
1207          overallLargest = max(overallLargest,largest);
1208          overallSmallest = min(overallSmallest,largest);
1209        }
1210      }
1211    } else {
1212      assert(scalingMethod==2);
1213      int numberPass=3;
1214#ifdef USE_OBJECTIVE
1215      // This will be used to help get scale factors
1216      double * objective = new double[numberColumns];
1217      memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
1218      double objScale=1.0;
1219#endif
1220      while (numberPass) {
1221        overallLargest=0.0;
1222        overallSmallest=1.0e50;
1223        numberPass--;
1224        // Geometric mean on row scales
1225        for (iRow=0;iRow<numberRows;iRow++) {
1226          if (usefulRow[iRow]) {
1227            CoinBigIndex j;
1228            largest=1.0e-20;
1229            smallest=1.0e50;
1230            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1231              int iColumn = column[j];
1232              if (usefulColumn[iColumn]) {
1233                double value = fabs(element[j]);
1234                // Don't bother with tiny elements
1235                if (value>1.0e-30) {
1236                  value *= columnScale[iColumn];
1237                  largest = max(largest,value);
1238                  smallest = min(smallest,value);
1239                }
1240              }
1241            }
1242            rowScale[iRow]=1.0/sqrt(smallest*largest);
1243            overallLargest = max(largest*rowScale[iRow],overallLargest);
1244            overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
1245          }
1246        }
1247#ifdef USE_OBJECTIVE
1248        largest=1.0e-20;
1249        smallest=1.0e50;
1250        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1251          if (usefulColumn[iColumn]) {
1252            double value = fabs(objective[iColumn]);
1253            // Don't bother with tiny elements
1254            if (value>1.0e-30) {
1255              value *= columnScale[iColumn];
1256              largest = max(largest,value);
1257              smallest = min(smallest,value);
1258            }
1259          }
1260        }
1261        objScale=1.0/sqrt(smallest*largest);
1262#endif
1263        model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
1264          <<overallSmallest
1265          <<overallLargest
1266          <<CoinMessageEol;
1267        // skip last column round
1268        if (numberPass==1)
1269          break;
1270        // Geometric mean on column scales
1271        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1272          if (usefulColumn[iColumn]) {
1273            CoinBigIndex j;
1274            largest=1.0e-20;
1275            smallest=1.0e50;
1276            for (j=columnStart[iColumn];
1277                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1278              iRow=row[j];
1279              double value = fabs(elementByColumn[j]);
1280              // Don't bother with tiny elements
1281              if (value>1.0e-30&&usefulRow[iRow]) {
1282                value *= rowScale[iRow];
1283                largest = max(largest,value);
1284                smallest = min(smallest,value);
1285              }
1286            }
1287#ifdef USE_OBJECTIVE
1288            if (fabs(objective[iColumn])>1.0e-30) {
1289              double value = fabs(objective[iColumn])*objScale;
1290              largest = max(largest,value);
1291              smallest = min(smallest,value);
1292            }
1293#endif
1294            columnScale[iColumn]=1.0/sqrt(smallest*largest);
1295          }
1296        }
1297      }
1298#ifdef USE_OBJECTIVE
1299      delete [] objective;
1300      printf("obj scale %g - use it if you want\n",objScale);
1301#endif
1302    }
1303    // final pass to scale columns so largest is reasonable
1304    // See what smallest will be if largest is 1.0
1305    overallSmallest=1.0e50;
1306    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1307      if (usefulColumn[iColumn]) {
1308        CoinBigIndex j;
1309        largest=1.0e-20;
1310        smallest=1.0e50;
1311        for (j=columnStart[iColumn];
1312             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1313          iRow=row[j];
1314          if(elementByColumn[j]&&usefulRow[iRow]) {
1315            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1316            largest = max(largest,value);
1317            smallest = min(smallest,value);
1318          }
1319        }
1320        if (overallSmallest*largest>smallest)
1321          overallSmallest = smallest/largest;
1322      }
1323    }
1324#define RANDOMIZE
1325#ifdef RANDOMIZE
1326    // randomize by up to 10%
1327    for (iRow=0;iRow<numberRows;iRow++) {
1328      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1329      rowScale[iRow] *= (1.0+0.1*value);
1330    }
1331#endif
1332    overallLargest=1.0;
1333    if (overallSmallest<1.0e-1)
1334      overallLargest = 1.0/sqrt(overallSmallest);
1335    overallLargest = min(1000.0,overallLargest);
1336    overallSmallest=1.0e50;
1337    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1338      if (usefulColumn[iColumn]) {
1339        CoinBigIndex j;
1340        largest=1.0e-20;
1341        smallest=1.0e50;
1342        for (j=columnStart[iColumn];
1343             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1344          iRow=row[j];
1345          if(elementByColumn[j]&&usefulRow[iRow]) {
1346            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1347            largest = max(largest,value);
1348            smallest = min(smallest,value);
1349          }
1350        }
1351        columnScale[iColumn]=overallLargest/largest;
1352#ifdef RANDOMIZE
1353        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1354        columnScale[iColumn] *= (1.0+0.1*value);
1355#endif
1356        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
1357      }
1358    }
1359    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
1360      <<overallSmallest
1361      <<overallLargest
1362      <<CoinMessageEol;
1363    delete [] usefulRow;
1364    delete [] usefulColumn;
1365    // If quadratic then make symmetric
1366    ClpObjective * obj = model->objectiveAsObject();
1367    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
1368    if (quadraticObj) {
1369      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1370      int numberXColumns = quadratic->getNumCols();
1371      if (numberXColumns<numberColumns) {
1372        // we assume symmetric
1373        int numberQuadraticColumns=0;
1374        int i;
1375        //const int * columnQuadratic = quadratic->getIndices();
1376        //const int * columnQuadraticStart = quadratic->getVectorStarts();
1377        const int * columnQuadraticLength = quadratic->getVectorLengths();
1378        for (i=0;i<numberXColumns;i++) {
1379          int length=columnQuadraticLength[i];
1380#ifndef CORRECT_COLUMN_COUNTS
1381          length=1;
1382#endif
1383          if (length)
1384            numberQuadraticColumns++;
1385        }
1386        int numberXRows = numberRows-numberQuadraticColumns;
1387        numberQuadraticColumns=0;
1388        for (i=0;i<numberXColumns;i++) { 
1389          int length=columnQuadraticLength[i];
1390#ifndef CORRECT_COLUMN_COUNTS
1391          length=1;
1392#endif
1393          if (length) {
1394            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
1395            numberQuadraticColumns++;
1396          }
1397        }   
1398        int numberQuadraticRows=0;
1399        for (i=0;i<numberXRows;i++) {
1400          // See if any in row quadratic
1401          int j;
1402          int numberQ=0;
1403          for (j=rowStart[i];j<rowStart[i+1];j++) {
1404            int iColumn = column[j];
1405            if (columnQuadraticLength[iColumn])
1406              numberQ++;
1407          }
1408#ifndef CORRECT_ROW_COUNTS
1409          numberQ=1;
1410#endif
1411          if (numberQ) {
1412            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
1413            numberQuadraticRows++;
1414          }
1415        }
1416        // and make sure Sj okay
1417        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
1418          CoinBigIndex j=columnStart[iColumn];
1419          assert(columnLength[iColumn]==1);
1420          int iRow=row[j];
1421          double value = fabs(elementByColumn[j]*rowScale[iRow]);
1422          columnScale[iColumn]=1.0/value;
1423        }
1424      }
1425    }
1426    model->setRowScale(rowScale);
1427    model->setColumnScale(columnScale);
1428    if (model->rowCopy()) {
1429      // need to replace row by row
1430      double * newElement = new double[numberColumns];
1431      // scale row copy
1432      for (iRow=0;iRow<numberRows;iRow++) {
1433        int j;
1434        double scale = rowScale[iRow];
1435        const double * elementsInThisRow = element + rowStart[iRow];
1436        const int * columnsInThisRow = column + rowStart[iRow];
1437        int number = rowStart[iRow+1]-rowStart[iRow];
1438        assert (number<=numberColumns);
1439        for (j=0;j<number;j++) {
1440          int iColumn = columnsInThisRow[j];
1441          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
1442        }
1443        rowCopy->replaceVector(iRow,number,newElement);
1444      }
1445      delete [] newElement;
1446    } else {
1447      // no row copy
1448      delete rowCopyBase;
1449    }
1450    return 0;
1451  }
1452}
1453/* Unpacks a column into an CoinIndexedvector
1454 */
1455void 
1456ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
1457                   int iColumn) const 
1458{
1459  const double * rowScale = model->rowScale();
1460  const int * row = matrix_->getIndices();
1461  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1462  const int * columnLength = matrix_->getVectorLengths(); 
1463  const double * elementByColumn = matrix_->getElements();
1464  CoinBigIndex i;
1465  if (!rowScale) {
1466    for (i=columnStart[iColumn];
1467         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1468      rowArray->add(row[i],elementByColumn[i]);
1469    }
1470  } else {
1471    // apply scaling
1472    double scale = model->columnScale()[iColumn];
1473    for (i=columnStart[iColumn];
1474         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1475      int iRow = row[i];
1476      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1477    }
1478  }
1479}
1480/* Unpacks a column into a CoinIndexedVector
1481** in packed format
1482Note that model is NOT const.  Bounds and objective could
1483be modified if doing column generation (just for this variable) */
1484void 
1485ClpPackedMatrix::unpackPacked(ClpSimplex * model,
1486                            CoinIndexedVector * rowArray,
1487                            int iColumn) const
1488{
1489  const double * rowScale = model->rowScale();
1490  const int * row = matrix_->getIndices();
1491  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1492  const int * columnLength = matrix_->getVectorLengths(); 
1493  const double * elementByColumn = matrix_->getElements();
1494  CoinBigIndex i;
1495  if (!rowScale) {
1496    int j=columnStart[iColumn];
1497    rowArray->createPacked(columnLength[iColumn],
1498                           row+j,elementByColumn+j);
1499  } else {
1500    // apply scaling
1501    double scale = model->columnScale()[iColumn];
1502    int * index = rowArray->getIndices();
1503    double * array = rowArray->denseVector();
1504    int number = 0;
1505    for (i=columnStart[iColumn];
1506         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1507      int iRow = row[i];
1508      array[number]=elementByColumn[i]*scale*rowScale[iRow];
1509      index[number++]=iRow;
1510    }
1511    rowArray->setNumElements(number);
1512    rowArray->setPackedMode(true);
1513  }
1514}
1515/* Adds multiple of a column into an CoinIndexedvector
1516      You can use quickAdd to add to vector */
1517void 
1518ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1519                   int iColumn, double multiplier) const 
1520{
1521  const double * rowScale = model->rowScale();
1522  const int * row = matrix_->getIndices();
1523  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1524  const int * columnLength = matrix_->getVectorLengths(); 
1525  const double * elementByColumn = matrix_->getElements();
1526  CoinBigIndex i;
1527  if (!rowScale) {
1528    for (i=columnStart[iColumn];
1529         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1530      int iRow = row[i];
1531      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
1532    }
1533  } else {
1534    // apply scaling
1535    double scale = model->columnScale()[iColumn]*multiplier;
1536    for (i=columnStart[iColumn];
1537         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1538      int iRow = row[i];
1539      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1540    }
1541  }
1542}
1543/* Checks if all elements are in valid range.  Can just
1544   return true if you are not paranoid.  For Clp I will
1545   probably expect no zeros.  Code can modify matrix to get rid of
1546   small elements.
1547*/
1548bool 
1549ClpPackedMatrix::allElementsInRange(ClpSimplex * model,
1550                                    double smallest, double largest)
1551{
1552  int iColumn;
1553  CoinBigIndex numberLarge=0;;
1554  CoinBigIndex numberSmall=0;;
1555  CoinBigIndex numberDuplicate=0;;
1556  int firstBadColumn=-1;
1557  int firstBadRow=-1;
1558  double firstBadElement=0.0;
1559  // get matrix data pointers
1560  const int * row = matrix_->getIndices();
1561  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1562  const int * columnLength = matrix_->getVectorLengths(); 
1563  const double * elementByColumn = matrix_->getElements();
1564  int numberColumns = matrix_->getNumCols();
1565  int numberRows = matrix_->getNumRows();
1566  int * mark = new int [numberRows];
1567  int i;
1568  for (i=0;i<numberRows;i++)
1569    mark[i]=-1;
1570  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1571    CoinBigIndex j;
1572    for (j=columnStart[iColumn];
1573         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1574      double value = fabs(elementByColumn[j]);
1575      int iRow = row[j];
1576      if (iRow<0||iRow>=numberRows) {
1577        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1578        return false;
1579      }
1580      if (mark[iRow]==-1) {
1581        mark[iRow]=j;
1582      } else {
1583        // duplicate
1584        numberDuplicate++;
1585      }
1586      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1587      if (!value)
1588        zeroElements_ = true; // there are zero elements
1589      if (value<smallest) {
1590        numberSmall++;
1591      } else if (value>largest) {
1592        numberLarge++;
1593        if (firstBadColumn<0) {
1594          firstBadColumn=iColumn;
1595          firstBadRow=row[j];
1596          firstBadElement=elementByColumn[j];
1597        }
1598      }
1599    }
1600    //clear mark
1601    for (j=columnStart[iColumn];
1602         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1603      int iRow = row[j];
1604      mark[iRow]=-1;
1605    }
1606  }
1607  delete [] mark;
1608  if (numberLarge) {
1609    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
1610      <<numberLarge
1611      <<firstBadColumn<<firstBadRow<<firstBadElement
1612      <<CoinMessageEol;
1613    return false;
1614  }
1615  if (numberSmall) 
1616    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
1617      <<numberSmall
1618      <<CoinMessageEol;
1619  if (numberDuplicate) 
1620    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
1621      <<numberDuplicate
1622      <<CoinMessageEol;
1623  if (numberDuplicate) 
1624    matrix_->eliminateDuplicates(smallest);
1625  else if (numberSmall) 
1626    matrix_->compress(smallest);
1627  // If smallest >0.0 then there can't be zero elements
1628  if (smallest>0.0)
1629    zeroElements_=false;
1630  return true;
1631}
1632/* Given positive integer weights for each row fills in sum of weights
1633   for each column (and slack).
1634   Returns weights vector
1635*/
1636CoinBigIndex * 
1637ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
1638{
1639  int numberRows = model->numberRows();
1640  int numberColumns =model->numberColumns();
1641  int number = numberRows+numberColumns;
1642  CoinBigIndex * weights = new CoinBigIndex[number];
1643  // get matrix data pointers
1644  const int * row = matrix_->getIndices();
1645  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1646  const int * columnLength = matrix_->getVectorLengths(); 
1647  int i;
1648  for (i=0;i<numberColumns;i++) {
1649    CoinBigIndex j;
1650    CoinBigIndex count=0;
1651    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1652      int iRow=row[j];
1653      count += inputWeights[iRow];
1654    }
1655    weights[i]=count;
1656  }
1657  for (i=0;i<numberRows;i++) {
1658    weights[i+numberColumns]=inputWeights[i];
1659  }
1660  return weights;
1661}
1662
1663
Note: See TracBrowser for help on using the repository browser.