source: branches/pre/ClpPackedMatrix.cpp @ 208

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

Trying to make faster

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