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

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

Improving reliability using IBM Burlington problems (Primal)

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