source: trunk/ClpPackedMatrix.cpp @ 284

Last change on this file since 284 was 284, checked in by jpfasano, 16 years ago

Modified to compile with MS Visual Studio Version 6

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