source: branches/devel-1/ClpPackedMatrix.cpp @ 34

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

Allow Presolve to work with gaps

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include <cstdio>
9
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12
13#include "ClpSimplex.hpp"
14#include "ClpFactorization.hpp"
15// at end to get min/max!
16#include "ClpPackedMatrix.hpp"
17#include "ClpMessage.hpp"
18
19//#############################################################################
20// Constructors / Destructor / Assignment
21//#############################################################################
22
23//-------------------------------------------------------------------
24// Default Constructor
25//-------------------------------------------------------------------
26ClpPackedMatrix::ClpPackedMatrix () 
27  : ClpMatrixBase(),
28    matrix_(NULL)
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 
41}
42
43//-------------------------------------------------------------------
44// assign matrix (for space reasons)
45//-------------------------------------------------------------------
46ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs) 
47: ClpMatrixBase()
48{ 
49  matrix_ = rhs;
50  setType(1);
51 
52}
53
54ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs) 
55: ClpMatrixBase()
56{ 
57  matrix_ = new CoinPackedMatrix(rhs);
58  setType(1);
59 
60}
61
62//-------------------------------------------------------------------
63// Destructor
64//-------------------------------------------------------------------
65ClpPackedMatrix::~ClpPackedMatrix ()
66{
67  delete matrix_;
68}
69
70//----------------------------------------------------------------
71// Assignment operator
72//-------------------------------------------------------------------
73ClpPackedMatrix &
74ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
75{
76  if (this != &rhs) {
77    ClpMatrixBase::operator=(rhs);
78    delete matrix_;
79    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
80  }
81  return *this;
82}
83//-------------------------------------------------------------------
84// Clone
85//-------------------------------------------------------------------
86ClpMatrixBase * ClpPackedMatrix::clone() const
87{
88  return new ClpPackedMatrix(*this);
89}
90
91/* Returns a new matrix in reverse order without gaps */
92ClpMatrixBase * 
93ClpPackedMatrix::reverseOrderedCopy() const
94{
95  ClpPackedMatrix * copy = new ClpPackedMatrix();
96  copy->matrix_= new CoinPackedMatrix();
97  copy->matrix_->reverseOrderedCopyOf(*matrix_);
98  copy->matrix_->removeGaps();
99  return copy;
100}
101//unscaled versions
102void 
103ClpPackedMatrix::times(double scalar,
104                   const double * x, double * y) const
105{
106  int iRow,iColumn;
107  // get matrix data pointers
108  const int * row = matrix_->getIndices();
109  const int * columnStart = matrix_->getVectorStarts();
110  const int * columnLength = matrix_->getVectorLengths(); 
111  const double * elementByColumn = matrix_->getElements();
112  int numberColumns = matrix_->getNumCols();
113  //memset(y,0,matrix_->getNumRows()*sizeof(double));
114  for (iColumn=0;iColumn<numberColumns;iColumn++) {
115    int j;
116    double value = scalar*x[iColumn];
117    if (value) {
118      for (j=columnStart[iColumn];
119           j<columnStart[iColumn]+columnLength[iColumn];j++) {
120        iRow=row[j];
121        y[iRow] += value*elementByColumn[j];
122      }
123    }
124  }
125}
126void 
127ClpPackedMatrix::transposeTimes(double scalar,
128                                const double * x, double * y) const
129{
130  int iColumn;
131  // get matrix data pointers
132  const int * row = matrix_->getIndices();
133  const int * columnStart = matrix_->getVectorStarts();
134  const int * columnLength = matrix_->getVectorLengths(); 
135  const double * elementByColumn = matrix_->getElements();
136  int numberColumns = matrix_->getNumCols();
137  //memset(y,0,numberColumns*sizeof(double));
138  for (iColumn=0;iColumn<numberColumns;iColumn++) {
139    int j;
140    double value=0.0;
141    for (j=columnStart[iColumn];
142         j<columnStart[iColumn]+columnLength[iColumn];j++) {
143      int jRow=row[j];
144      value += x[jRow]*elementByColumn[j];
145    }
146    y[iColumn] += value*scalar;
147  }
148}
149void 
150ClpPackedMatrix::times(double scalar,
151                       const double * x, double * y,
152                       const double * rowScale, 
153                       const double * columnScale) const
154{
155  int iRow,iColumn;
156  // get matrix data pointers
157  const int * row = matrix_->getIndices();
158  const int * columnStart = matrix_->getVectorStarts();
159  const int * columnLength = matrix_->getVectorLengths(); 
160  const double * elementByColumn = matrix_->getElements();
161  int numberColumns = matrix_->getNumCols();
162  //memset(y,0,matrix_->getNumRows()*sizeof(double));
163  for (iColumn=0;iColumn<numberColumns;iColumn++) {
164    int j;
165    double value = x[iColumn];
166    if (value) {
167      // scaled
168      value *= scalar*columnScale[iColumn];
169      for (j=columnStart[iColumn];
170           j<columnStart[iColumn]+columnLength[iColumn];j++) {
171        iRow=row[j];
172        y[iRow] += value*elementByColumn[j];
173      }
174    }
175  }
176  int numberRows = getNumRows();
177  for (iRow=0;iRow<numberRows;iRow++) {
178    y[iRow] *= rowScale[iRow];
179  }
180}
181void 
182ClpPackedMatrix::transposeTimes( double scalar,
183                                 const double * x, double * y,
184                                 const double * rowScale, 
185                                 const double * columnScale) const
186{
187  int iColumn;
188  // get matrix data pointers
189  const int * row = matrix_->getIndices();
190  const int * columnStart = matrix_->getVectorStarts();
191  const int * columnLength = matrix_->getVectorLengths(); 
192  const double * elementByColumn = matrix_->getElements();
193  int numberColumns = matrix_->getNumCols();
194  //memset(y,0,numberColumns*sizeof(double));
195  for (iColumn=0;iColumn<numberColumns;iColumn++) {
196    int j;
197    double value=0.0;
198    // scaled
199    for (j=columnStart[iColumn];
200         j<columnStart[iColumn]+columnLength[iColumn];j++) {
201      int jRow=row[j];
202      value += x[jRow]*elementByColumn[j]*rowScale[jRow];
203    }
204    y[iColumn] += value*scalar*columnScale[iColumn];
205  }
206}
207/* Return <code>x * A + y</code> in <code>z</code>.
208        Squashes small elements and knows about ClpSimplex */
209void 
210ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
211                              const CoinIndexedVector * rowArray,
212                              CoinIndexedVector * y,
213                              CoinIndexedVector * columnArray) const
214{
215  columnArray->clear();
216  double * pi = rowArray->denseVector();
217  int numberNonZero=0;
218  int * index = columnArray->getIndices();
219  double * array = columnArray->denseVector();
220  int numberInRowArray = rowArray->getNumElements();
221  // maybe I need one in OsiSimplex
222  double zeroTolerance = model->factorization()->zeroTolerance();
223  int numberRows = model->numberRows();
224  const ClpMatrixBase * rowCopy = model->rowCopy();
225  // make sure row copy is of correct type
226  // if not we would have to have another transposeTimes
227  assert (!rowCopy||rowCopy->type()==1);
228  if (numberInRowArray>0.3*numberRows||!rowCopy) {
229    // do by column
230    int iColumn;
231    // get matrix data pointers
232    const int * row = matrix_->getIndices();
233    const int * columnStart = matrix_->getVectorStarts();
234    const int * columnLength = matrix_->getVectorLengths(); 
235    const double * elementByColumn = matrix_->getElements();
236    double * markVector = y->denseVector(); // probably empty
237    const double * rowScale = model->rowScale();
238    int numberColumns = model->numberColumns();
239    if (!rowScale) {
240      for (iColumn=0;iColumn<numberColumns;iColumn++) {
241        double value = markVector[iColumn];
242        markVector[iColumn]=0.0;
243        double value2 = 0.0;
244        int j;
245        for (j=columnStart[iColumn];
246             j<columnStart[iColumn]+columnLength[iColumn];j++) {
247          int iRow = row[j];
248          value2 += pi[iRow]*elementByColumn[j];
249        }
250        value += value2*scalar;
251        if (fabs(value)>zeroTolerance) {
252          index[numberNonZero++]=iColumn;
253          array[iColumn]=value;
254        }
255      }
256    } else {
257      // scaled
258      for (iColumn=0;iColumn<numberColumns;iColumn++) {
259        double value = markVector[iColumn];
260        markVector[iColumn]=0.0;
261        int j;
262        const double * columnScale = model->columnScale();
263        double value2 = 0.0;
264        for (j=columnStart[iColumn];
265             j<columnStart[iColumn]+columnLength[iColumn];j++) {
266          int iRow = row[j];
267          value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
268        }
269        value += value2*scalar*columnScale[iColumn];
270        if (fabs(value)>zeroTolerance) {
271          index[numberNonZero++]=iColumn;
272          array[iColumn]=value;
273        }
274      }
275    }
276  } else if (numberInRowArray>2||y->getNumElements()) {
277    // do by rows
278    // ** Row copy is already scaled
279    int iRow;
280    double * markVector = y->denseVector(); // probably empty .. but
281    int * mark = y->getIndices();
282    int numberOriginal=y->getNumElements();
283    int i;
284    for (i=0;i<numberOriginal;i++) {
285      int iColumn = mark[i];
286      index[i]=iColumn;
287      array[iColumn]=markVector[iColumn];
288      markVector[iColumn]=0.0;
289    }
290    numberNonZero=numberOriginal;
291    // and set up mark as char array
292    char * marked = (char *) markVector;
293    for (i=0;i<numberOriginal;i++) {
294      int iColumn = index[i];
295      marked[iColumn]=0;
296    }
297
298    const int * column = rowCopy->getIndices();
299    const int * rowStart = rowCopy->getVectorStarts();
300    const double * element = rowCopy->getElements();
301    const int * whichRow = rowArray->getIndices();
302    for (i=0;i<numberInRowArray;i++) {
303      iRow = whichRow[i]; 
304      double value = pi[iRow]*scalar;
305      int j;
306      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
307        int iColumn = column[j];
308        if (!marked[iColumn]) {
309          marked[iColumn]=1;
310          index[numberNonZero++]=iColumn;
311        }
312        array[iColumn] += value*element[j];
313      }
314    }
315    // get rid of tiny values and zero out marked
316    numberOriginal=numberNonZero;
317    numberNonZero=0;
318    for (i=0;i<numberOriginal;i++) {
319      int iColumn = index[i];
320      marked[iColumn]=0;
321      if (fabs(array[iColumn])>zeroTolerance) {
322        index[numberNonZero++]=iColumn;
323      } else {
324        array[iColumn]=0.0;
325      }
326    }
327  } else if (numberInRowArray==2) {
328    // do by rows when two rows
329    int iRow;
330    int numberOriginal;
331    int i;
332    numberNonZero=0;
333
334    const int * column = rowCopy->getIndices();
335    const int * rowStart = rowCopy->getVectorStarts();
336    const double * element = rowCopy->getElements();
337    const int * whichRow = rowArray->getIndices();
338    double value;
339    iRow = whichRow[0]; 
340    value = pi[iRow]*scalar;
341    int j;
342    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
343      int iColumn = column[j];
344      double value2 = value*element[j];
345      index[numberNonZero++]=iColumn;
346      array[iColumn] = value2;
347    }
348    iRow = whichRow[1]; 
349    value = pi[iRow]*scalar;
350    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
351      int iColumn = column[j];
352      double value2 = value*element[j];
353      // I am assuming no zeros in matrix
354      if (array[iColumn])
355        value2 += array[iColumn];
356      else
357        index[numberNonZero++]=iColumn;
358      array[iColumn] = value2;
359    }
360    // get rid of tiny values and zero out marked
361    numberOriginal=numberNonZero;
362    numberNonZero=0;
363    for (i=0;i<numberOriginal;i++) {
364      int iColumn = index[i];
365      if (fabs(array[iColumn])>zeroTolerance) {
366        index[numberNonZero++]=iColumn;
367      } else {
368        array[iColumn]=0.0;
369      }
370    }
371  } else if (numberInRowArray==1) {
372    // Just one row
373    int iRow=rowArray->getIndices()[0];
374    numberNonZero=0;
375    const int * column = rowCopy->getIndices();
376    const int * rowStart = rowCopy->getVectorStarts();
377    const double * element = rowCopy->getElements();
378    double value = pi[iRow]*scalar;
379    int j;
380    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
381      int iColumn = column[j];
382      double value2 = value*element[j];
383      if (fabs(value2)>zeroTolerance) {
384        index[numberNonZero++]=iColumn;
385        array[iColumn] = value2;
386      }
387    }
388  }
389  columnArray->setNumElements(numberNonZero);
390  y->setNumElements(0);
391}
392/* Return <code>x *A in <code>z</code> but
393   just for indices in y.
394   Squashes small elements and knows about ClpSimplex */
395void 
396ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
397                              const CoinIndexedVector * rowArray,
398                              const CoinIndexedVector * y,
399                              CoinIndexedVector * columnArray) const
400{
401  columnArray->clear();
402  double * pi = rowArray->denseVector();
403  int numberNonZero=0;
404  int * index = columnArray->getIndices();
405  double * array = columnArray->denseVector();
406  // maybe I need one in OsiSimplex
407  double zeroTolerance = model->factorization()->zeroTolerance();
408  int jColumn;
409  // get matrix data pointers
410  const int * row = matrix_->getIndices();
411  const int * columnStart = matrix_->getVectorStarts();
412  const int * columnLength = matrix_->getVectorLengths(); 
413  const double * elementByColumn = matrix_->getElements();
414  const double * rowScale = model->rowScale();
415  int numberToDo = y->getNumElements();
416  const int * which = y->getIndices();
417  if (!rowScale) {
418    for (jColumn=0;jColumn<numberToDo;jColumn++) {
419      int iColumn = which[jColumn];
420      double value = 0.0;
421      int j;
422      for (j=columnStart[iColumn];
423           j<columnStart[iColumn]+columnLength[iColumn];j++) {
424        int iRow = row[j];
425        value += pi[iRow]*elementByColumn[j];
426      }
427      if (fabs(value)>zeroTolerance) {
428        index[numberNonZero++]=iColumn;
429        array[iColumn]=value;
430      }
431    }
432  } else {
433    // scaled
434    for (jColumn=0;jColumn<numberToDo;jColumn++) {
435      int iColumn = which[jColumn];
436      double value = 0.0;
437      int j;
438      const double * columnScale = model->columnScale();
439      for (j=columnStart[iColumn];
440           j<columnStart[iColumn]+columnLength[iColumn];j++) {
441        int iRow = row[j];
442        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
443      }
444      value *= columnScale[iColumn];
445      if (fabs(value)>zeroTolerance) {
446        index[numberNonZero++]=iColumn;
447        array[iColumn]=value;
448      }
449    }
450  }
451}
452/* Returns number of elements in basis
453   column is basic if entry >=0 */
454ClpBigIndex
455ClpPackedMatrix::numberInBasis(const int * columnIsBasic) const 
456{
457  int i;
458  int numberColumns = getNumCols();
459  const int * columnLength = matrix_->getVectorLengths(); 
460  ClpBigIndex numberElements=0;
461  for (i=0;i<numberColumns;i++) {
462    if (columnIsBasic[i]>=0) {
463      numberElements += columnLength[i];
464    }
465  }
466  return numberElements;
467}
468// Fills in basis (Returns number of elements and updates numberBasic)
469ClpBigIndex
470ClpPackedMatrix::fillBasis(const ClpSimplex * model,
471                                const int * columnIsBasic, int & numberBasic,
472                                int * indexRowU, int * indexColumnU,
473                                double * elementU) const 
474{
475  const double * rowScale = model->rowScale();
476  const int * row = matrix_->getIndices();
477  const int * columnStart = matrix_->getVectorStarts();
478  const int * columnLength = matrix_->getVectorLengths(); 
479  const double * elementByColumn = matrix_->getElements();
480  int i;
481  int numberColumns = getNumCols();
482  ClpBigIndex numberElements=0;
483  if (!rowScale) {
484    // no scaling
485    for (i=0;i<numberColumns;i++) {
486      if (columnIsBasic[i]>=0) {
487        int j;
488        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
489          indexRowU[numberElements]=row[j];
490          indexColumnU[numberElements]=numberBasic;
491          elementU[numberElements++]=elementByColumn[j];
492        }
493        numberBasic++;
494      }
495    }
496  } else {
497    // scaling
498    const double * columnScale = model->columnScale();
499    for (i=0;i<numberColumns;i++) {
500      if (columnIsBasic[i]>=0) {
501        int j;
502        double scale = columnScale[i];
503        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
504          int iRow = row[j];
505          indexRowU[numberElements]=iRow;
506          indexColumnU[numberElements]=numberBasic;
507          elementU[numberElements++]=elementByColumn[j]*scale*rowScale[iRow];
508        }
509        numberBasic++;
510      }
511    }
512  }
513  return numberElements;
514}
515// Creates scales for column copy (rowCopy in model may be modified)
516int 
517ClpPackedMatrix::scale(ClpSimplex * model) const 
518{
519  ClpMatrixBase * rowCopyBase=model->rowCopy();
520  if (!rowCopyBase) {
521    // temporary copy
522    rowCopyBase = reverseOrderedCopy();
523  }
524  ClpPackedMatrix* rowCopy =
525    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
526
527  // Make sure it is really a ClpPackedMatrix
528  assert (rowCopy!=NULL);
529  int numberRows = model->numberRows();
530  int numberColumns = model->numberColumns();
531  const int * column = rowCopy->getIndices();
532  const int * rowStart = rowCopy->getVectorStarts();
533  const double * element = rowCopy->getElements();
534  double * rowScale = new double [numberRows];
535  double * columnScale = new double [numberColumns];
536  // we are going to mark bits we are interested in
537  char * usefulRow = new char [numberRows];
538  char * usefulColumn = new char [numberColumns];
539  double * rowLower = model->rowLower();
540  double * rowUpper = model->rowUpper();
541  double * columnLower = model->columnLower();
542  double * columnUpper = model->columnUpper();
543  int iColumn, iRow;
544  // mark free rows
545  for (iRow=0;iRow<numberRows;iRow++) {
546    usefulRow[iRow]=0;
547    if (rowUpper[iRow]<1.0e20||
548        rowLower[iRow]>-1.0e20)
549      usefulRow[iRow]=1;
550  }
551  // mark empty and fixed columns
552  // also see if worth scaling
553  assert (model->scalingFlag()==1); // dynamic not implemented
554  double largest=0.0;
555  double smallest=1.0e50;
556  // get matrix data pointers
557  const int * row = matrix_->getIndices();
558  const int * columnStart = matrix_->getVectorStarts();
559  const int * columnLength = matrix_->getVectorLengths(); 
560  const double * elementByColumn = matrix_->getElements();
561  for (iColumn=0;iColumn<numberColumns;iColumn++) {
562    int j;
563    char useful=0;
564    if (columnUpper[iColumn]>
565        columnLower[iColumn]+1.0e-9) {
566      for (j=columnStart[iColumn];
567           j<columnStart[iColumn]+columnLength[iColumn];j++) {
568        iRow=row[j];
569        if(elementByColumn[j]&&usefulRow[iRow]) {
570          useful=1;
571          largest = max(largest,fabs(elementByColumn[j]));
572          smallest = min(smallest,fabs(elementByColumn[j]));
573        }
574      }
575    }
576    usefulColumn[iColumn]=useful;
577  }
578  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
579    <<smallest<<largest
580    <<CoinMessageEol;
581  if (smallest>=0.5&&largest<=2.0) {
582    // don't bother scaling
583    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
584      <<CoinMessageEol;
585    delete [] rowScale;
586    delete [] usefulRow;
587    delete [] columnScale;
588    delete [] usefulColumn;
589    return 1;
590  } else {
591    // and see if there any empty rows
592    for (iRow=0;iRow<numberRows;iRow++) {
593      if (usefulRow[iRow]) {
594        int j;
595        int useful=0;
596        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
597          int iColumn = column[j];
598          if (usefulColumn[iColumn]) {
599            useful=1;
600            break;
601          }
602        }
603        usefulRow[iRow]=useful;
604      }
605    }
606    CoinFillN ( rowScale, numberRows,1.0);
607    CoinFillN ( columnScale, numberColumns,1.0);
608    int numberPass=3;
609    double overallLargest;
610    double overallSmallest;
611    while (numberPass) {
612      overallLargest=0.0;
613      overallSmallest=1.0e50;
614      numberPass--;
615      // Geometric mean on row scales
616      for (iRow=0;iRow<numberRows;iRow++) {
617        if (usefulRow[iRow]) {
618          int j;
619          largest=1.0e-20;
620          smallest=1.0e50;
621          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
622            int iColumn = column[j];
623            if (usefulColumn[iColumn]) {
624              double value = fabs(element[j]*columnScale[iColumn]);
625              largest = max(largest,value);
626              smallest = min(smallest,value);
627            }
628          }
629          rowScale[iRow]=1.0/sqrt(smallest*largest);
630          overallLargest = max(largest*rowScale[iRow],overallLargest);
631          overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
632        }
633      }
634      model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
635        <<overallSmallest
636        <<overallLargest
637        <<CoinMessageEol;
638      // skip last column round
639      if (numberPass==1)
640        break;
641      // Geometric mean on column scales
642      for (iColumn=0;iColumn<numberColumns;iColumn++) {
643        if (usefulColumn[iColumn]) {
644          int j;
645          largest=1.0e-20;
646          smallest=1.0e50;
647          for (j=columnStart[iColumn];
648               j<columnStart[iColumn]+columnLength[iColumn];j++) {
649            iRow=row[j];
650            if(elementByColumn[j]&&usefulRow[iRow]) {
651              double value = fabs(elementByColumn[j]*rowScale[iRow]);
652              largest = max(largest,value);
653              smallest = min(smallest,value);
654            }
655          }
656          columnScale[iColumn]=1.0/sqrt(smallest*largest);
657        }
658      }
659    }
660    // final pass to scale columns so largest is 1.0
661    overallLargest=0.0;
662    overallSmallest=1.0e50;
663    for (iColumn=0;iColumn<numberColumns;iColumn++) {
664      if (usefulColumn[iColumn]) {
665        int j;
666        largest=1.0e-20;
667        smallest=1.0e50;
668        for (j=columnStart[iColumn];
669             j<columnStart[iColumn]+columnLength[iColumn];j++) {
670          iRow=row[j];
671          if(elementByColumn[j]&&usefulRow[iRow]) {
672            double value = fabs(elementByColumn[j]*rowScale[iRow]);
673            largest = max(largest,value);
674            smallest = min(smallest,value);
675          }
676        }
677        columnScale[iColumn]=1.0/largest;
678        overallLargest = max(overallLargest,largest*columnScale[iColumn]);
679        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
680      }
681    }
682    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
683      <<overallSmallest
684      <<overallLargest
685      <<CoinMessageEol;
686   
687    delete [] usefulRow;
688    delete [] usefulColumn;
689    model->setRowScale(rowScale);
690    model->setColumnScale(columnScale);
691    if (model->rowCopy()) {
692      // need to replace row by row
693      double * newElement = new double[numberColumns];
694      // scale row copy
695      for (iRow=0;iRow<numberRows;iRow++) {
696        int j;
697        double scale = rowScale[iRow];
698        const double * elementsInThisRow = element + rowStart[iRow];
699        const int * columnsInThisRow = column + rowStart[iRow];
700        int number = rowStart[iRow+1]-rowStart[iRow];
701        assert (number<=numberColumns);
702        for (j=0;j<number;j++) {
703          int iColumn = columnsInThisRow[j];
704          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
705        }
706        rowCopy->replaceVector(iRow,number,newElement);
707      }
708      delete [] newElement;
709    } else {
710      // no row copy
711      delete rowCopyBase;
712    }
713    return 0;
714  }
715}
716// Creates row copy and scales if necessary
717ClpMatrixBase * 
718ClpPackedMatrix::scaledRowCopy(ClpSimplex * model) const
719{
720  ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
721  ClpPackedMatrix* rowCopy =
722    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
723
724  // Make sure it is really a ClpPackedMatrix
725  assert (rowCopy!=NULL);
726
727  const double * rowScale = model->rowScale();
728  const double * columnScale = model->columnScale();
729
730  if (rowScale) {
731    // scale row copy
732    int numberRows = model->numberRows();
733    int numberColumns = model->numberColumns();
734    const int * column = rowCopy->getIndices();
735    const int * rowStart = rowCopy->getVectorStarts();
736    const double * element = rowCopy->getElements();
737    int iRow;
738    // need to replace row by row
739    double * newElement = new double[numberColumns];
740    // scale row copy
741    for (iRow=0;iRow<numberRows;iRow++) {
742      int j;
743      double scale = rowScale[iRow];
744      const double * elementsInThisRow = element + rowStart[iRow];
745      const int * columnsInThisRow = column + rowStart[iRow];
746      int number = rowStart[iRow+1]-rowStart[iRow];
747      assert (number<=numberColumns);
748      for (j=0;j<number;j++) {
749        int iColumn = columnsInThisRow[j];
750        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
751      }
752      rowCopy->replaceVector(iRow,number,newElement);
753    }
754    delete [] newElement;
755  }
756  return rowCopyBase;
757}
758/* Unpacks a column into an CoinIndexedvector
759      Note that model is NOT const.  Bounds and objective could
760      be modified if doing column generation */
761void 
762ClpPackedMatrix::unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
763                   int iColumn) const 
764{
765  const double * rowScale = model->rowScale();
766  const int * row = matrix_->getIndices();
767  const int * columnStart = matrix_->getVectorStarts();
768  const int * columnLength = matrix_->getVectorLengths(); 
769  const double * elementByColumn = matrix_->getElements();
770  int i;
771  if (!rowScale) {
772    for (i=columnStart[iColumn];
773         i<columnStart[iColumn]+columnLength[iColumn];i++) {
774      rowArray->add(row[i],elementByColumn[i]);
775    }
776  } else {
777    // apply scaling
778    double scale = model->columnScale()[iColumn];
779    for (i=columnStart[iColumn];
780         i<columnStart[iColumn]+columnLength[iColumn];i++) {
781      int iRow = row[i];
782      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
783    }
784  }
785}
786/* Adds multiple of a column into an CoinIndexedvector
787      You can use quickAdd to add to vector */
788void 
789ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
790                   int iColumn, double multiplier) const 
791{
792  const double * rowScale = model->rowScale();
793  const int * row = matrix_->getIndices();
794  const int * columnStart = matrix_->getVectorStarts();
795  const int * columnLength = matrix_->getVectorLengths(); 
796  const double * elementByColumn = matrix_->getElements();
797  int i;
798  if (!rowScale) {
799    for (i=columnStart[iColumn];
800         i<columnStart[iColumn]+columnLength[iColumn];i++) {
801      int iRow = row[i];
802      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
803    }
804  } else {
805    // apply scaling
806    double scale = model->columnScale()[iColumn]*multiplier;
807    for (i=columnStart[iColumn];
808         i<columnStart[iColumn]+columnLength[iColumn];i++) {
809      int iRow = row[i];
810      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
811    }
812  }
813}
814/* Checks if all elements are in valid range.  Can just
815   return true if you are not paranoid.  For Clp I will
816   probably expect no zeros.  Code can modify matrix to get rid of
817   small elements.
818*/
819bool 
820ClpPackedMatrix::allElementsInRange(ClpSimplex * model,
821                                    double smallest, double largest)
822{
823  int iColumn;
824  int numberLarge=0;;
825  int numberSmall=0;;
826  int firstBadColumn=-1;
827  int firstBadRow=-1;
828  double firstBadElement=0.0;
829  // get matrix data pointers
830  const int * row = matrix_->getIndices();
831  const int * columnStart = matrix_->getVectorStarts();
832  const int * columnLength = matrix_->getVectorLengths(); 
833  const double * elementByColumn = matrix_->getElements();
834  int numberColumns = matrix_->getNumCols();
835  for (iColumn=0;iColumn<numberColumns;iColumn++) {
836    int j;
837    for (j=columnStart[iColumn];
838         j<columnStart[iColumn]+columnLength[iColumn];j++) {
839      double value = fabs(elementByColumn[j]);
840      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
841      if (value<smallest) {
842        numberSmall++;
843      } else if (value>largest) {
844        numberLarge++;
845        if (firstBadColumn<0) {
846          firstBadColumn=iColumn;
847          firstBadRow=row[j];
848          firstBadElement=elementByColumn[j];
849        }
850      }
851    }
852  }
853  if (numberLarge) {
854    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
855      <<numberLarge
856      <<firstBadColumn<<firstBadRow<<firstBadElement
857      <<CoinMessageEol;
858    return false;
859  }
860  if (numberSmall) {
861    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
862      <<numberSmall
863      <<CoinMessageEol;
864    matrix_->compress(smallest);
865  }
866  return true;
867}
868
869
Note: See TracBrowser for help on using the repository browser.