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

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

Synchronizing

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