source: trunk/ClpPackedMatrix.cpp @ 118

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

Adding Network matrix and PlusMinusOne?

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