source: trunk/ClpPackedMatrix.cpp @ 145

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

Zeros in matrix, and minor stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.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// 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
94/* Returns a new matrix in reverse order without gaps */
95ClpMatrixBase * 
96ClpPackedMatrix::reverseOrderedCopy() const
97{
98  ClpPackedMatrix * copy = new ClpPackedMatrix();
99  copy->matrix_= new CoinPackedMatrix();
100  copy->matrix_->reverseOrderedCopyOf(*matrix_);
101  copy->matrix_->removeGaps();
102  return copy;
103}
104//unscaled versions
105void 
106ClpPackedMatrix::times(double scalar,
107                   const double * x, double * y) const
108{
109  int iRow,iColumn;
110  // get matrix data pointers
111  const int * row = matrix_->getIndices();
112  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
113  const int * columnLength = matrix_->getVectorLengths(); 
114  const double * elementByColumn = matrix_->getElements();
115  int numberColumns = matrix_->getNumCols();
116  //memset(y,0,matrix_->getNumRows()*sizeof(double));
117  for (iColumn=0;iColumn<numberColumns;iColumn++) {
118    CoinBigIndex j;
119    double value = scalar*x[iColumn];
120    if (value) {
121      for (j=columnStart[iColumn];
122           j<columnStart[iColumn]+columnLength[iColumn];j++) {
123        iRow=row[j];
124        y[iRow] += value*elementByColumn[j];
125      }
126    }
127  }
128}
129void 
130ClpPackedMatrix::transposeTimes(double scalar,
131                                const double * x, double * y) const
132{
133  int iColumn;
134  // get matrix data pointers
135  const int * row = matrix_->getIndices();
136  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
137  const int * columnLength = matrix_->getVectorLengths(); 
138  const double * elementByColumn = matrix_->getElements();
139  int numberColumns = matrix_->getNumCols();
140  //memset(y,0,numberColumns*sizeof(double));
141  for (iColumn=0;iColumn<numberColumns;iColumn++) {
142    CoinBigIndex j;
143    double value=0.0;
144    for (j=columnStart[iColumn];
145         j<columnStart[iColumn]+columnLength[iColumn];j++) {
146      int jRow=row[j];
147      value += x[jRow]*elementByColumn[j];
148    }
149    y[iColumn] += value*scalar;
150  }
151}
152void 
153ClpPackedMatrix::times(double scalar,
154                       const double * x, double * y,
155                       const double * rowScale, 
156                       const double * columnScale) const
157{
158  int iRow,iColumn;
159  // get matrix data pointers
160  const int * row = matrix_->getIndices();
161  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
162  const int * columnLength = matrix_->getVectorLengths(); 
163  const double * elementByColumn = matrix_->getElements();
164  int numberColumns = matrix_->getNumCols();
165  //memset(y,0,matrix_->getNumRows()*sizeof(double));
166  for (iColumn=0;iColumn<numberColumns;iColumn++) {
167    CoinBigIndex j;
168    double value = x[iColumn];
169    if (value) {
170      // scaled
171      value *= scalar*columnScale[iColumn];
172      for (j=columnStart[iColumn];
173           j<columnStart[iColumn]+columnLength[iColumn];j++) {
174        iRow=row[j];
175        y[iRow] += value*elementByColumn[j];
176      }
177    }
178  }
179  int numberRows = getNumRows();
180  for (iRow=0;iRow<numberRows;iRow++) {
181    y[iRow] *= rowScale[iRow];
182  }
183}
184void 
185ClpPackedMatrix::transposeTimes( double scalar,
186                                 const double * x, double * y,
187                                 const double * rowScale, 
188                                 const double * columnScale) const
189{
190  int iColumn;
191  // get matrix data pointers
192  const int * row = matrix_->getIndices();
193  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
194  const int * columnLength = matrix_->getVectorLengths(); 
195  const double * elementByColumn = matrix_->getElements();
196  int numberColumns = matrix_->getNumCols();
197  //memset(y,0,numberColumns*sizeof(double));
198  for (iColumn=0;iColumn<numberColumns;iColumn++) {
199    CoinBigIndex j;
200    double value=0.0;
201    // scaled
202    for (j=columnStart[iColumn];
203         j<columnStart[iColumn]+columnLength[iColumn];j++) {
204      int jRow=row[j];
205      value += x[jRow]*elementByColumn[j]*rowScale[jRow];
206    }
207    y[iColumn] += value*scalar*columnScale[iColumn];
208  }
209}
210/* Return <code>x * A + y</code> in <code>z</code>.
211        Squashes small elements and knows about ClpSimplex */
212void 
213ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
214                              const CoinIndexedVector * rowArray,
215                              CoinIndexedVector * y,
216                              CoinIndexedVector * columnArray) const
217{
218  columnArray->clear();
219  double * pi = rowArray->denseVector();
220  int numberNonZero=0;
221  int * index = columnArray->getIndices();
222  double * array = columnArray->denseVector();
223  int numberInRowArray = rowArray->getNumElements();
224  // maybe I need one in OsiSimplex
225  double zeroTolerance = model->factorization()->zeroTolerance();
226  int numberRows = model->numberRows();
227  ClpPackedMatrix* rowCopy =
228    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
229  if (numberInRowArray>0.333*numberRows||!rowCopy) {
230    // do by column
231    int iColumn;
232    // get matrix data pointers
233    const int * row = matrix_->getIndices();
234    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
235    const int * columnLength = matrix_->getVectorLengths(); 
236    const double * elementByColumn = matrix_->getElements();
237    const double * rowScale = model->rowScale();
238    int numberColumns = model->numberColumns();
239    if (!y->getNumElements()) {
240      if (!rowScale) {
241        if (scalar==-1.0) {
242          for (iColumn=0;iColumn<numberColumns;iColumn++) {
243            double value = 0.0;
244            CoinBigIndex j;
245            for (j=columnStart[iColumn];
246                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
247              int iRow = row[j];
248              value += pi[iRow]*elementByColumn[j];
249            }
250            if (fabs(value)>zeroTolerance) {
251              index[numberNonZero++]=iColumn;
252              array[iColumn]=-value;
253            }
254          }
255        } else if (scalar==1.0) {
256          for (iColumn=0;iColumn<numberColumns;iColumn++) {
257            double value = 0.0;
258            CoinBigIndex j;
259            for (j=columnStart[iColumn];
260                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
261              int iRow = row[j];
262              value += pi[iRow]*elementByColumn[j];
263            }
264            if (fabs(value)>zeroTolerance) {
265              index[numberNonZero++]=iColumn;
266              array[iColumn]=value;
267            }
268          }
269        } else {
270          for (iColumn=0;iColumn<numberColumns;iColumn++) {
271            double value = 0.0;
272            CoinBigIndex j;
273            for (j=columnStart[iColumn];
274                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
275              int iRow = row[j];
276              value += pi[iRow]*elementByColumn[j];
277            }
278            value *= scalar;
279            if (fabs(value)>zeroTolerance) {
280              index[numberNonZero++]=iColumn;
281              array[iColumn]=value;
282            }
283          }
284        }
285      } else {
286        // scaled
287        if (scalar==-1.0) {
288          for (iColumn=0;iColumn<numberColumns;iColumn++) {
289            double value = 0.0;
290            CoinBigIndex j;
291            const double * columnScale = model->columnScale();
292            for (j=columnStart[iColumn];
293                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
294              int iRow = row[j];
295              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
296            }
297            value *= columnScale[iColumn];
298            if (fabs(value)>zeroTolerance) {
299              index[numberNonZero++]=iColumn;
300              array[iColumn]=-value;
301            }
302          }
303        } else if (scalar==1.0) {
304          for (iColumn=0;iColumn<numberColumns;iColumn++) {
305            double value = 0.0;
306            CoinBigIndex j;
307            const double * columnScale = model->columnScale();
308            for (j=columnStart[iColumn];
309                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
310              int iRow = row[j];
311              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
312            }
313            value *= columnScale[iColumn];
314            if (fabs(value)>zeroTolerance) {
315              index[numberNonZero++]=iColumn;
316              array[iColumn]=value;
317            }
318          }
319        } else {
320          for (iColumn=0;iColumn<numberColumns;iColumn++) {
321            double value = 0.0;
322            CoinBigIndex j;
323            const double * columnScale = model->columnScale();
324            for (j=columnStart[iColumn];
325                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
326              int iRow = row[j];
327              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
328            }
329            value *= scalar*columnScale[iColumn];
330            if (fabs(value)>zeroTolerance) {
331              index[numberNonZero++]=iColumn;
332              array[iColumn]=value;
333            }
334          }
335        }
336      }
337    } else {
338      double * markVector = y->denseVector(); // not empty
339      if (!rowScale) {
340        for (iColumn=0;iColumn<numberColumns;iColumn++) {
341          double value = markVector[iColumn];
342          markVector[iColumn]=0.0;
343          double value2 = 0.0;
344          CoinBigIndex j;
345          for (j=columnStart[iColumn];
346               j<columnStart[iColumn]+columnLength[iColumn];j++) {
347            int iRow = row[j];
348            value2 += pi[iRow]*elementByColumn[j];
349          }
350          value += value2*scalar;
351          if (fabs(value)>zeroTolerance) {
352            index[numberNonZero++]=iColumn;
353            array[iColumn]=value;
354          }
355        }
356      } else {
357        // scaled
358        for (iColumn=0;iColumn<numberColumns;iColumn++) {
359          double value = markVector[iColumn];
360          markVector[iColumn]=0.0;
361          CoinBigIndex j;
362          const double * columnScale = model->columnScale();
363          double value2 = 0.0;
364          for (j=columnStart[iColumn];
365               j<columnStart[iColumn]+columnLength[iColumn];j++) {
366            int iRow = row[j];
367            value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
368          }
369          value += value2*scalar*columnScale[iColumn];
370          if (fabs(value)>zeroTolerance) {
371            index[numberNonZero++]=iColumn;
372            array[iColumn]=value;
373          }
374        }
375      }
376    }
377    columnArray->setNumElements(numberNonZero);
378    y->setNumElements(0);
379  } else {
380    // do by row
381    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
382  }
383}
384/* Return <code>x * A + y</code> in <code>z</code>.
385        Squashes small elements and knows about ClpSimplex */
386void 
387ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
388                              const CoinIndexedVector * rowArray,
389                              CoinIndexedVector * y,
390                              CoinIndexedVector * columnArray) const
391{
392  columnArray->clear();
393  double * pi = rowArray->denseVector();
394  int numberNonZero=0;
395  int * index = columnArray->getIndices();
396  double * array = columnArray->denseVector();
397  int numberInRowArray = rowArray->getNumElements();
398  // maybe I need one in OsiSimplex
399  double zeroTolerance = model->factorization()->zeroTolerance();
400  const int * column = getIndices();
401  const CoinBigIndex * rowStart = getVectorStarts();
402  const double * element = getElements();
403  const int * whichRow = rowArray->getIndices();
404  if (numberInRowArray>2||y->getNumElements()) {
405    // do by rows
406    // ** Row copy is already scaled
407    int iRow;
408    double * markVector = y->denseVector(); // probably empty .. but
409    int * mark = y->getIndices();
410    int numberOriginal=y->getNumElements();
411    int i;
412    for (i=0;i<numberOriginal;i++) {
413      int iColumn = mark[i];
414      index[i]=iColumn;
415      array[iColumn]=markVector[iColumn];
416      markVector[iColumn]=0.0;
417    }
418    numberNonZero=numberOriginal;
419    // and set up mark as char array
420    char * marked = (char *) markVector;
421    for (i=0;i<numberOriginal;i++) {
422      int iColumn = index[i];
423      marked[iColumn]=0;
424    }
425
426    for (i=0;i<numberInRowArray;i++) {
427      iRow = whichRow[i]; 
428      double value = pi[iRow]*scalar;
429      CoinBigIndex j;
430      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
431        int iColumn = column[j];
432        if (!marked[iColumn]) {
433          marked[iColumn]=1;
434          index[numberNonZero++]=iColumn;
435        }
436        array[iColumn] += value*element[j];
437      }
438    }
439    // get rid of tiny values and zero out marked
440    numberOriginal=numberNonZero;
441    numberNonZero=0;
442    for (i=0;i<numberOriginal;i++) {
443      int iColumn = index[i];
444      marked[iColumn]=0;
445      if (fabs(array[iColumn])>zeroTolerance) {
446        index[numberNonZero++]=iColumn;
447      } else {
448        array[iColumn]=0.0;
449      }
450    }
451  } else if (numberInRowArray==2) {
452    // do by rows when two rows
453    int iRow;
454    int numberOriginal;
455    int i;
456    numberNonZero=0;
457
458    double value;
459    iRow = whichRow[0]; 
460    value = pi[iRow]*scalar;
461    CoinBigIndex j;
462    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
463      int iColumn = column[j];
464      double value2 = value*element[j];
465      index[numberNonZero++]=iColumn;
466      array[iColumn] = value2;
467    }
468    iRow = whichRow[1]; 
469    value = pi[iRow]*scalar;
470    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
471      int iColumn = column[j];
472      double value2 = value*element[j];
473      // I am assuming no zeros in matrix
474      if (array[iColumn])
475        value2 += array[iColumn];
476      else
477        index[numberNonZero++]=iColumn;
478      array[iColumn] = value2;
479    }
480    // get rid of tiny values and zero out marked
481    numberOriginal=numberNonZero;
482    numberNonZero=0;
483    for (i=0;i<numberOriginal;i++) {
484      int iColumn = index[i];
485      if (fabs(array[iColumn])>zeroTolerance) {
486        index[numberNonZero++]=iColumn;
487      } else {
488        array[iColumn]=0.0;
489      }
490    }
491  } else if (numberInRowArray==1) {
492    // Just one row
493    int iRow=rowArray->getIndices()[0];
494    numberNonZero=0;
495    double value = pi[iRow]*scalar;
496    CoinBigIndex j;
497    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
498      int iColumn = column[j];
499      double value2 = value*element[j];
500      if (fabs(value2)>zeroTolerance) {
501        index[numberNonZero++]=iColumn;
502        array[iColumn] = value2;
503      }
504    }
505  }
506  columnArray->setNumElements(numberNonZero);
507  y->setNumElements(0);
508}
509/* Return <code>x *A in <code>z</code> but
510   just for indices in y.
511   Squashes small elements and knows about ClpSimplex */
512void 
513ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
514                              const CoinIndexedVector * rowArray,
515                              const CoinIndexedVector * y,
516                              CoinIndexedVector * columnArray) const
517{
518  columnArray->clear();
519  double * pi = rowArray->denseVector();
520  int numberNonZero=0;
521  int * index = columnArray->getIndices();
522  double * array = columnArray->denseVector();
523  // maybe I need one in OsiSimplex
524  double zeroTolerance = model->factorization()->zeroTolerance();
525  int jColumn;
526  // get matrix data pointers
527  const int * row = matrix_->getIndices();
528  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
529  const int * columnLength = matrix_->getVectorLengths(); 
530  const double * elementByColumn = matrix_->getElements();
531  const double * rowScale = model->rowScale();
532  int numberToDo = y->getNumElements();
533  const int * which = y->getIndices();
534  if (!rowScale) {
535    for (jColumn=0;jColumn<numberToDo;jColumn++) {
536      int iColumn = which[jColumn];
537      double value = 0.0;
538      CoinBigIndex j;
539      for (j=columnStart[iColumn];
540           j<columnStart[iColumn]+columnLength[iColumn];j++) {
541        int iRow = row[j];
542        value += pi[iRow]*elementByColumn[j];
543      }
544      if (fabs(value)>zeroTolerance) {
545        index[numberNonZero++]=iColumn;
546        array[iColumn]=value;
547      }
548    }
549  } else {
550    // scaled
551    for (jColumn=0;jColumn<numberToDo;jColumn++) {
552      int iColumn = which[jColumn];
553      double value = 0.0;
554      CoinBigIndex j;
555      const double * columnScale = model->columnScale();
556      for (j=columnStart[iColumn];
557           j<columnStart[iColumn]+columnLength[iColumn];j++) {
558        int iRow = row[j];
559        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
560      }
561      value *= columnScale[iColumn];
562      if (fabs(value)>zeroTolerance) {
563        index[numberNonZero++]=iColumn;
564        array[iColumn]=value;
565      }
566    }
567  }
568}
569/* Returns number of elements in basis
570   column is basic if entry >=0 */
571CoinBigIndex
572ClpPackedMatrix::numberInBasis(const int * columnIsBasic) const 
573{
574  int i;
575  int numberColumns = getNumCols();
576  const int * columnLength = matrix_->getVectorLengths(); 
577  CoinBigIndex numberElements=0;
578  if (!zeroElements_) {
579    for (i=0;i<numberColumns;i++) {
580      if (columnIsBasic[i]>=0) {
581        numberElements += columnLength[i];
582      }
583    }
584  } else {
585    // there are zero elements so need to look more closely
586    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
587    const double * elementByColumn = matrix_->getElements();
588    for (i=0;i<numberColumns;i++) {
589      if (columnIsBasic[i]>=0) {
590        CoinBigIndex j;
591        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
592          if (elementByColumn[j])
593            numberElements++;
594        }
595      }
596    }
597  }
598  return numberElements;
599}
600// Fills in basis (Returns number of elements and updates numberBasic)
601CoinBigIndex
602ClpPackedMatrix::fillBasis(const ClpSimplex * model,
603                                const int * columnIsBasic, int & numberBasic,
604                                int * indexRowU, int * indexColumnU,
605                                double * elementU) const 
606{
607  const double * rowScale = model->rowScale();
608  const int * row = matrix_->getIndices();
609  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
610  const int * columnLength = matrix_->getVectorLengths(); 
611  const double * elementByColumn = matrix_->getElements();
612  int i;
613  int numberColumns = getNumCols();
614  CoinBigIndex numberElements=0;
615  if (!zeroElements_) {
616    if (!rowScale) {
617      // no scaling
618      for (i=0;i<numberColumns;i++) {
619        if (columnIsBasic[i]>=0) {
620          CoinBigIndex j;
621          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
622            indexRowU[numberElements]=row[j];
623            indexColumnU[numberElements]=numberBasic;
624            elementU[numberElements++]=elementByColumn[j];
625          }
626          numberBasic++;
627        }
628      }
629    } else {
630      // scaling
631      const double * columnScale = model->columnScale();
632      for (i=0;i<numberColumns;i++) {
633        if (columnIsBasic[i]>=0) {
634          CoinBigIndex j;
635          double scale = columnScale[i];
636          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
637            int iRow = row[j];
638            indexRowU[numberElements]=iRow;
639            indexColumnU[numberElements]=numberBasic;
640            elementU[numberElements++]=elementByColumn[j]*scale*rowScale[iRow];
641          }
642          numberBasic++;
643        }
644      }
645    }
646  } else {
647    // there are zero elements so need to look more closely
648    if (!rowScale) {
649      // no scaling
650      for (i=0;i<numberColumns;i++) {
651        if (columnIsBasic[i]>=0) {
652          CoinBigIndex j;
653          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
654            double value = elementByColumn[j];
655            if (value) {
656              indexRowU[numberElements]=row[j];
657              indexColumnU[numberElements]=numberBasic;
658              elementU[numberElements++]=value;
659            }
660          }
661          numberBasic++;
662        }
663      }
664    } else {
665      // scaling
666      const double * columnScale = model->columnScale();
667      for (i=0;i<numberColumns;i++) {
668        if (columnIsBasic[i]>=0) {
669          CoinBigIndex j;
670          double scale = columnScale[i];
671          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
672            double value = elementByColumn[j];
673            if (value) {
674              int iRow = row[j];
675              indexRowU[numberElements]=iRow;
676              indexColumnU[numberElements]=numberBasic;
677              elementU[numberElements++]=value*scale*rowScale[iRow];
678            }
679          }
680          numberBasic++;
681        }
682      }
683    }
684  }
685  return numberElements;
686}
687// Creates scales for column copy (rowCopy in model may be modified)
688int 
689ClpPackedMatrix::scale(ClpSimplex * model) const 
690{
691  ClpMatrixBase * rowCopyBase=model->rowCopy();
692  if (!rowCopyBase) {
693    // temporary copy
694    rowCopyBase = reverseOrderedCopy();
695  }
696  ClpPackedMatrix* rowCopy =
697    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
698
699  // Make sure it is really a ClpPackedMatrix
700  assert (rowCopy!=NULL);
701  int numberRows = model->numberRows();
702  int numberColumns = model->numberColumns();
703  const int * column = rowCopy->getIndices();
704  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
705  const double * element = rowCopy->getElements();
706  double * rowScale = new double [numberRows];
707  double * columnScale = new double [numberColumns];
708  // we are going to mark bits we are interested in
709  char * usefulRow = new char [numberRows];
710  char * usefulColumn = new char [numberColumns];
711  double * rowLower = model->rowLower();
712  double * rowUpper = model->rowUpper();
713  double * columnLower = model->columnLower();
714  double * columnUpper = model->columnUpper();
715  int iColumn, iRow;
716  // mark free rows
717  for (iRow=0;iRow<numberRows;iRow++) {
718    usefulRow[iRow]=0;
719    if (rowUpper[iRow]<1.0e20||
720        rowLower[iRow]>-1.0e20)
721      usefulRow[iRow]=1;
722  }
723  // mark empty and fixed columns
724  // also see if worth scaling
725  assert (model->scalingFlag()==1); // dynamic not implemented
726  double largest=0.0;
727  double smallest=1.0e50;
728  // get matrix data pointers
729  const int * row = matrix_->getIndices();
730  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
731  const int * columnLength = matrix_->getVectorLengths(); 
732  const double * elementByColumn = matrix_->getElements();
733  for (iColumn=0;iColumn<numberColumns;iColumn++) {
734    CoinBigIndex j;
735    char useful=0;
736    if (columnUpper[iColumn]>
737        columnLower[iColumn]+1.0e-9) {
738      for (j=columnStart[iColumn];
739           j<columnStart[iColumn]+columnLength[iColumn];j++) {
740        iRow=row[j];
741        if(elementByColumn[j]&&usefulRow[iRow]) {
742          useful=1;
743          largest = max(largest,fabs(elementByColumn[j]));
744          smallest = min(smallest,fabs(elementByColumn[j]));
745        }
746      }
747    }
748    usefulColumn[iColumn]=useful;
749  }
750  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
751    <<smallest<<largest
752    <<CoinMessageEol;
753  if (smallest>=0.5&&largest<=2.0) {
754    // don't bother scaling
755    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
756      <<CoinMessageEol;
757    delete [] rowScale;
758    delete [] usefulRow;
759    delete [] columnScale;
760    delete [] usefulColumn;
761    return 1;
762  } else {
763    // and see if there any empty rows
764    for (iRow=0;iRow<numberRows;iRow++) {
765      if (usefulRow[iRow]) {
766        CoinBigIndex j;
767        int useful=0;
768        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
769          int iColumn = column[j];
770          if (usefulColumn[iColumn]) {
771            useful=1;
772            break;
773          }
774        }
775        usefulRow[iRow]=useful;
776      }
777    }
778    ClpFillN ( rowScale, numberRows,1.0);
779    ClpFillN ( columnScale, numberColumns,1.0);
780    int numberPass=3;
781    double overallLargest;
782    double overallSmallest;
783    while (numberPass) {
784      overallLargest=0.0;
785      overallSmallest=1.0e50;
786      numberPass--;
787      // Geometric mean on row scales
788      for (iRow=0;iRow<numberRows;iRow++) {
789        if (usefulRow[iRow]) {
790          CoinBigIndex j;
791          largest=1.0e-20;
792          smallest=1.0e50;
793          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
794            int iColumn = column[j];
795            if (usefulColumn[iColumn]) {
796              double value = fabs(element[j]);
797              // Don't bother with tiny elements
798              if (value>1.0e-30) {
799                value *= columnScale[iColumn];
800                largest = max(largest,value);
801                smallest = min(smallest,value);
802              }
803            }
804          }
805          rowScale[iRow]=1.0/sqrt(smallest*largest);
806          overallLargest = max(largest*rowScale[iRow],overallLargest);
807          overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
808        }
809      }
810      model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
811        <<overallSmallest
812        <<overallLargest
813        <<CoinMessageEol;
814      // skip last column round
815      if (numberPass==1)
816        break;
817      // Geometric mean on column scales
818      for (iColumn=0;iColumn<numberColumns;iColumn++) {
819        if (usefulColumn[iColumn]) {
820          CoinBigIndex j;
821          largest=1.0e-20;
822          smallest=1.0e50;
823          for (j=columnStart[iColumn];
824               j<columnStart[iColumn]+columnLength[iColumn];j++) {
825            iRow=row[j];
826            double value = fabs(elementByColumn[j]);
827            // Don't bother with tiny elements
828            if (value>1.0e-30&&usefulRow[iRow]) {
829              value *= rowScale[iRow];
830              largest = max(largest,value);
831              smallest = min(smallest,value);
832            }
833          }
834          columnScale[iColumn]=1.0/sqrt(smallest*largest);
835        }
836      }
837    }
838    // final pass to scale columns so largest is 1.0
839    overallLargest=0.0;
840    overallSmallest=1.0e50;
841    for (iColumn=0;iColumn<numberColumns;iColumn++) {
842      if (usefulColumn[iColumn]) {
843        CoinBigIndex j;
844        largest=1.0e-20;
845        smallest=1.0e50;
846        for (j=columnStart[iColumn];
847             j<columnStart[iColumn]+columnLength[iColumn];j++) {
848          iRow=row[j];
849          if(elementByColumn[j]&&usefulRow[iRow]) {
850            double value = fabs(elementByColumn[j]*rowScale[iRow]);
851            largest = max(largest,value);
852            smallest = min(smallest,value);
853          }
854        }
855        columnScale[iColumn]=1.0/largest;
856        overallLargest = max(overallLargest,largest*columnScale[iColumn]);
857        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
858      }
859    }
860    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
861      <<overallSmallest
862      <<overallLargest
863      <<CoinMessageEol;
864   
865    delete [] usefulRow;
866    delete [] usefulColumn;
867    model->setRowScale(rowScale);
868    model->setColumnScale(columnScale);
869    if (model->rowCopy()) {
870      // need to replace row by row
871      double * newElement = new double[numberColumns];
872      // scale row copy
873      for (iRow=0;iRow<numberRows;iRow++) {
874        int j;
875        double scale = rowScale[iRow];
876        const double * elementsInThisRow = element + rowStart[iRow];
877        const int * columnsInThisRow = column + rowStart[iRow];
878        int number = rowStart[iRow+1]-rowStart[iRow];
879        assert (number<=numberColumns);
880        for (j=0;j<number;j++) {
881          int iColumn = columnsInThisRow[j];
882          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
883        }
884        rowCopy->replaceVector(iRow,number,newElement);
885      }
886      delete [] newElement;
887    } else {
888      // no row copy
889      delete rowCopyBase;
890    }
891    return 0;
892  }
893}
894/* Unpacks a column into an CoinIndexedvector
895      Note that model is NOT const.  Bounds and objective could
896      be modified if doing column generation */
897void 
898ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
899                   int iColumn) const 
900{
901  const double * rowScale = model->rowScale();
902  const int * row = matrix_->getIndices();
903  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
904  const int * columnLength = matrix_->getVectorLengths(); 
905  const double * elementByColumn = matrix_->getElements();
906  CoinBigIndex i;
907  if (!rowScale) {
908    for (i=columnStart[iColumn];
909         i<columnStart[iColumn]+columnLength[iColumn];i++) {
910      rowArray->add(row[i],elementByColumn[i]);
911    }
912  } else {
913    // apply scaling
914    double scale = model->columnScale()[iColumn];
915    for (i=columnStart[iColumn];
916         i<columnStart[iColumn]+columnLength[iColumn];i++) {
917      int iRow = row[i];
918      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
919    }
920  }
921}
922/* Adds multiple of a column into an CoinIndexedvector
923      You can use quickAdd to add to vector */
924void 
925ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
926                   int iColumn, double multiplier) const 
927{
928  const double * rowScale = model->rowScale();
929  const int * row = matrix_->getIndices();
930  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
931  const int * columnLength = matrix_->getVectorLengths(); 
932  const double * elementByColumn = matrix_->getElements();
933  CoinBigIndex i;
934  if (!rowScale) {
935    for (i=columnStart[iColumn];
936         i<columnStart[iColumn]+columnLength[iColumn];i++) {
937      int iRow = row[i];
938      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
939    }
940  } else {
941    // apply scaling
942    double scale = model->columnScale()[iColumn]*multiplier;
943    for (i=columnStart[iColumn];
944         i<columnStart[iColumn]+columnLength[iColumn];i++) {
945      int iRow = row[i];
946      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
947    }
948  }
949}
950/* Checks if all elements are in valid range.  Can just
951   return true if you are not paranoid.  For Clp I will
952   probably expect no zeros.  Code can modify matrix to get rid of
953   small elements.
954*/
955bool 
956ClpPackedMatrix::allElementsInRange(ClpSimplex * model,
957                                    double smallest, double largest)
958{
959  int iColumn;
960  CoinBigIndex numberLarge=0;;
961  CoinBigIndex numberSmall=0;;
962  CoinBigIndex numberDuplicate=0;;
963  int firstBadColumn=-1;
964  int firstBadRow=-1;
965  double firstBadElement=0.0;
966  // get matrix data pointers
967  const int * row = matrix_->getIndices();
968  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
969  const int * columnLength = matrix_->getVectorLengths(); 
970  const double * elementByColumn = matrix_->getElements();
971  int numberColumns = matrix_->getNumCols();
972  int numberRows = matrix_->getNumRows();
973  int * mark = new int [numberRows];
974  int i;
975  for (i=0;i<numberRows;i++)
976    mark[i]=-1;
977  for (iColumn=0;iColumn<numberColumns;iColumn++) {
978    CoinBigIndex j;
979    for (j=columnStart[iColumn];
980         j<columnStart[iColumn]+columnLength[iColumn];j++) {
981      double value = fabs(elementByColumn[j]);
982      int iRow = row[j];
983      if (iRow<0||iRow>=numberRows) {
984        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
985        return false;
986      }
987      if (mark[iRow]==-1) {
988        mark[iRow]=j;
989      } else {
990        // duplicate
991        numberDuplicate++;
992      }
993      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
994      if (!value)
995        zeroElements_ = true; // there are zero elements
996      if (value<smallest) {
997        numberSmall++;
998      } else if (value>largest) {
999        numberLarge++;
1000        if (firstBadColumn<0) {
1001          firstBadColumn=iColumn;
1002          firstBadRow=row[j];
1003          firstBadElement=elementByColumn[j];
1004        }
1005      }
1006    }
1007    //clear mark
1008    for (j=columnStart[iColumn];
1009         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1010      int iRow = row[j];
1011      mark[iRow]=-1;
1012    }
1013  }
1014  delete [] mark;
1015  if (numberLarge) {
1016    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
1017      <<numberLarge
1018      <<firstBadColumn<<firstBadRow<<firstBadElement
1019      <<CoinMessageEol;
1020    return false;
1021  }
1022  if (numberSmall) 
1023    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
1024      <<numberSmall
1025      <<CoinMessageEol;
1026  if (numberDuplicate) 
1027    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
1028      <<numberDuplicate
1029      <<CoinMessageEol;
1030  if (numberDuplicate) 
1031    matrix_->eliminateDuplicates(smallest);
1032  else if (numberSmall) 
1033    matrix_->compress(smallest);
1034  // If smallest >0.0 then there can't be zero elements
1035  if (smallest>0.0)
1036    zeroElements_=false;
1037  return true;
1038}
1039
1040
Note: See TracBrowser for help on using the repository browser.