source: branches/pre/ClpPackedMatrix.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago

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