source: trunk/ClpPackedMatrix.cpp @ 2

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

Adding Clp to development branch

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