source: trunk/ClpPackedMatrix.cpp @ 426

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

correcting errors

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 67.0 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#include "ClpQuadraticObjective.hpp"
14// at end to get min/max!
15#include "ClpPackedMatrix.hpp"
16#include "ClpMessage.hpp"
17
18//#############################################################################
19// Constructors / Destructor / Assignment
20//#############################################################################
21
22//-------------------------------------------------------------------
23// Default Constructor
24//-------------------------------------------------------------------
25ClpPackedMatrix::ClpPackedMatrix () 
26  : ClpMatrixBase(),
27    matrix_(NULL),
28    numberActiveColumns_(0),
29    zeroElements_(false),
30    hasGaps_(true)
31{
32  setType(1);
33}
34
35//-------------------------------------------------------------------
36// Copy constructor
37//-------------------------------------------------------------------
38ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs) 
39: ClpMatrixBase(rhs)
40{ 
41  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
42  numberActiveColumns_ = rhs.numberActiveColumns_;
43  zeroElements_ = rhs.zeroElements_;
44  hasGaps_ = rhs.hasGaps_;
45  int numberRows = getNumRows();
46  if (rhs.rhsOffset_&&numberRows) {
47    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
48  } else {
49    rhsOffset_=NULL;
50  }
51 
52}
53
54//-------------------------------------------------------------------
55// assign matrix (for space reasons)
56//-------------------------------------------------------------------
57ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs) 
58: ClpMatrixBase()
59{ 
60  matrix_ = rhs;
61  zeroElements_ = false;
62  hasGaps_ = true;
63  numberActiveColumns_ = matrix_->getNumCols();
64  setType(1);
65 
66}
67
68ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs) 
69: ClpMatrixBase()
70{ 
71  matrix_ = new CoinPackedMatrix(rhs);
72  numberActiveColumns_ = matrix_->getNumCols();
73  zeroElements_ = false;
74  hasGaps_ = true;
75  setType(1);
76 
77}
78
79//-------------------------------------------------------------------
80// Destructor
81//-------------------------------------------------------------------
82ClpPackedMatrix::~ClpPackedMatrix ()
83{
84  delete matrix_;
85}
86
87//----------------------------------------------------------------
88// Assignment operator
89//-------------------------------------------------------------------
90ClpPackedMatrix &
91ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
92{
93  if (this != &rhs) {
94    ClpMatrixBase::operator=(rhs);
95    delete matrix_;
96    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
97    numberActiveColumns_ = rhs.numberActiveColumns_;
98    zeroElements_ = rhs.zeroElements_;
99    hasGaps_ = rhs.hasGaps_;
100  }
101  return *this;
102}
103//-------------------------------------------------------------------
104// Clone
105//-------------------------------------------------------------------
106ClpMatrixBase * ClpPackedMatrix::clone() const
107{
108  return new ClpPackedMatrix(*this);
109}
110/* Subset clone (without gaps).  Duplicates are allowed
111   and order is as given */
112ClpMatrixBase * 
113ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
114                              int numberColumns, 
115                              const int * whichColumns) const 
116{
117  return new ClpPackedMatrix(*this, numberRows, whichRows,
118                                   numberColumns, whichColumns);
119}
120/* Subset constructor (without gaps).  Duplicates are allowed
121   and order is as given */
122ClpPackedMatrix::ClpPackedMatrix (
123                       const ClpPackedMatrix & rhs,
124                       int numberRows, const int * whichRows,
125                       int numberColumns, const int * whichColumns)
126: ClpMatrixBase(rhs)
127{
128  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
129                                 numberColumns,whichColumns);
130  numberActiveColumns_ = matrix_->getNumCols();
131  zeroElements_ = rhs.zeroElements_;
132  hasGaps_ = rhs.hasGaps_;
133}
134ClpPackedMatrix::ClpPackedMatrix (
135                       const CoinPackedMatrix & rhs,
136                       int numberRows, const int * whichRows,
137                       int numberColumns, const int * whichColumns)
138: ClpMatrixBase()
139{
140  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
141                                 numberColumns,whichColumns);
142  numberActiveColumns_ = matrix_->getNumCols();
143  zeroElements_ = false;
144  hasGaps_=true;
145  setType(1);
146}
147
148/* Returns a new matrix in reverse order without gaps */
149ClpMatrixBase * 
150ClpPackedMatrix::reverseOrderedCopy() const
151{
152  ClpPackedMatrix * copy = new ClpPackedMatrix();
153  copy->matrix_= new CoinPackedMatrix();
154  copy->matrix_->reverseOrderedCopyOf(*matrix_);
155  copy->matrix_->removeGaps();
156  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
157  copy->hasGaps_=false;
158  return copy;
159}
160//unscaled versions
161void 
162ClpPackedMatrix::times(double scalar,
163                   const double * x, double * y) const
164{
165  int iRow,iColumn;
166  // get matrix data pointers
167  const int * row = matrix_->getIndices();
168  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
169  const int * columnLength = matrix_->getVectorLengths(); 
170  const double * elementByColumn = matrix_->getElements();
171  //memset(y,0,matrix_->getNumRows()*sizeof(double));
172  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
173    CoinBigIndex j;
174    double value = scalar*x[iColumn];
175    if (value) {
176      for (j=columnStart[iColumn];
177           j<columnStart[iColumn]+columnLength[iColumn];j++) {
178        iRow=row[j];
179        y[iRow] += value*elementByColumn[j];
180      }
181    }
182  }
183}
184void 
185ClpPackedMatrix::transposeTimes(double scalar,
186                                const double * x, double * y) const
187{
188  int iColumn;
189  // get matrix data pointers
190  const int * row = matrix_->getIndices();
191  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
192  const int * columnLength = matrix_->getVectorLengths(); 
193  const double * elementByColumn = matrix_->getElements();
194  //memset(y,0,numberActiveColumns_*sizeof(double));
195  for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
196    CoinBigIndex j;
197    double value=0.0;
198    for (j=columnStart[iColumn];
199         j<columnStart[iColumn]+columnLength[iColumn];j++) {
200      int jRow=row[j];
201      value += x[jRow]*elementByColumn[j];
202    }
203    y[iColumn] += value*scalar;
204  }
205}
206void 
207ClpPackedMatrix::times(double scalar,
208                       const double * x, double * y,
209                       const double * rowScale, 
210                       const double * columnScale) const
211{
212  if (rowScale) {
213    int iRow,iColumn;
214    // get matrix data pointers
215    const int * row = matrix_->getIndices();
216    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
217    const int * columnLength = matrix_->getVectorLengths(); 
218    const double * elementByColumn = matrix_->getElements();
219    //memset(y,0,matrix_->getNumRows()*sizeof(double));
220    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
221      CoinBigIndex j;
222      double value = x[iColumn];
223      if (value) {
224        // scaled
225        value *= scalar*columnScale[iColumn];
226        for (j=columnStart[iColumn];
227             j<columnStart[iColumn]+columnLength[iColumn];j++) {
228          iRow=row[j];
229          y[iRow] += value*elementByColumn[j]*rowScale[iRow];
230        }
231      }
232    }
233  } else {
234    times(scalar,x,y);
235  }
236}
237void 
238ClpPackedMatrix::transposeTimes( double scalar,
239                                 const double * x, double * y,
240                                 const double * rowScale, 
241                                 const double * columnScale) const
242{
243  if (rowScale) {
244    int iColumn;
245    // get matrix data pointers
246    const int * row = matrix_->getIndices();
247    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
248    const int * columnLength = matrix_->getVectorLengths(); 
249    const double * elementByColumn = matrix_->getElements();
250    //memset(y,0,numberActiveColumns_*sizeof(double));
251    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
252      CoinBigIndex j;
253      double value=0.0;
254      // scaled
255      for (j=columnStart[iColumn];
256           j<columnStart[iColumn]+columnLength[iColumn];j++) {
257        int jRow=row[j];
258      value += x[jRow]*elementByColumn[j]*rowScale[jRow];
259      }
260      y[iColumn] += value*scalar*columnScale[iColumn];
261    }
262  } else {
263    transposeTimes(scalar,x,y);
264  }
265}
266/* Return <code>x * A + y</code> in <code>z</code>.
267        Squashes small elements and knows about ClpSimplex */
268void 
269ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
270                              const CoinIndexedVector * rowArray,
271                              CoinIndexedVector * y,
272                              CoinIndexedVector * columnArray) const
273{
274  columnArray->clear();
275  double * pi = rowArray->denseVector();
276  int numberNonZero=0;
277  int * index = columnArray->getIndices();
278  double * array = columnArray->denseVector();
279  int numberInRowArray = rowArray->getNumElements();
280  // maybe I need one in OsiSimplex
281  double zeroTolerance = model->factorization()->zeroTolerance();
282  int numberRows = model->numberRows();
283#ifndef NO_RTTI
284  ClpPackedMatrix* rowCopy =
285    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
286#else
287  ClpPackedMatrix* rowCopy =
288    static_cast< ClpPackedMatrix*>(model->rowCopy());
289#endif
290  bool packed = rowArray->packedMode();
291  double factor = 0.3;
292  // We may not want to do by row if there may be cache problems
293  // It would be nice to find L2 cache size - for moment 512K
294  // Be slightly optimistic
295  if (numberActiveColumns_*sizeof(double)>1000000) {
296    if (numberRows*10<numberActiveColumns_)
297      factor=0.1;
298    else if (numberRows*4<numberActiveColumns_)
299      factor=0.15;
300    else if (numberRows*2<numberActiveColumns_)
301      factor=0.2;
302    //if (model->numberIterations()%50==0)
303    //printf("%d nonzero\n",numberInRowArray);
304  }
305  assert (!y->getNumElements());
306  if (numberInRowArray>factor*numberRows||!rowCopy) {
307    // do by column
308    // If no gaps - can do a bit faster
309    if (!hasGaps_&&true) {
310      transposeTimesByColumn( model,  scalar,
311                              rowArray, y, columnArray);
312      return;
313    }
314    int iColumn;
315    // get matrix data pointers
316    const int * row = matrix_->getIndices();
317    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
318    const int * columnLength = matrix_->getVectorLengths(); 
319    const double * elementByColumn = matrix_->getElements();
320    const double * rowScale = model->rowScale();
321    if (packed) {
322      // need to expand pi into y
323      assert(y->capacity()>=numberRows);
324      double * piOld = pi;
325      pi = y->denseVector();
326      const int * whichRow = rowArray->getIndices();
327      int i;
328      if (!rowScale) {
329        // modify pi so can collapse to one loop
330        for (i=0;i<numberInRowArray;i++) {
331          int iRow = whichRow[i];
332          pi[iRow]=scalar*piOld[i];
333        }
334        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
335          double value = 0.0;
336          CoinBigIndex j;
337          for (j=columnStart[iColumn];
338               j<columnStart[iColumn]+columnLength[iColumn];j++) {
339            int iRow = row[j];
340            value += pi[iRow]*elementByColumn[j];
341          }
342          if (fabs(value)>zeroTolerance) {
343            array[numberNonZero]=value;
344            index[numberNonZero++]=iColumn;
345          }
346        }
347      } else {
348        // scaled
349        // modify pi so can collapse to one loop
350        for (i=0;i<numberInRowArray;i++) {
351          int iRow = whichRow[i];
352          pi[iRow]=scalar*piOld[i]*rowScale[iRow];
353        }
354        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
355          double value = 0.0;
356          CoinBigIndex j;
357          const double * columnScale = model->columnScale();
358          for (j=columnStart[iColumn];
359               j<columnStart[iColumn]+columnLength[iColumn];j++) {
360            int iRow = row[j];
361            value += pi[iRow]*elementByColumn[j];
362          }
363          value *= columnScale[iColumn];
364          if (fabs(value)>zeroTolerance) {
365            array[numberNonZero]=value;
366            index[numberNonZero++]=iColumn;
367          }
368        }
369      }
370      // zero out
371      for (i=0;i<numberInRowArray;i++) {
372        int iRow = whichRow[i];
373        pi[iRow]=0.0;
374      }
375    } else {
376      if (!rowScale) {
377        if (scalar==-1.0) {
378          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
379            double value = 0.0;
380            CoinBigIndex j;
381            for (j=columnStart[iColumn];
382                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
383              int iRow = row[j];
384              value += pi[iRow]*elementByColumn[j];
385            }
386            if (fabs(value)>zeroTolerance) {
387              index[numberNonZero++]=iColumn;
388              array[iColumn]=-value;
389            }
390          }
391        } else if (scalar==1.0) {
392          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
393            double value = 0.0;
394            CoinBigIndex j;
395            for (j=columnStart[iColumn];
396                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
397              int iRow = row[j];
398              value += pi[iRow]*elementByColumn[j];
399            }
400            if (fabs(value)>zeroTolerance) {
401              index[numberNonZero++]=iColumn;
402              array[iColumn]=value;
403            }
404          }
405        } else {
406          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
407            double value = 0.0;
408            CoinBigIndex 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            value *= scalar;
415            if (fabs(value)>zeroTolerance) {
416              index[numberNonZero++]=iColumn;
417              array[iColumn]=value;
418            }
419          }
420        }
421      } else {
422        // scaled
423        if (scalar==-1.0) {
424          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
425            double value = 0.0;
426            CoinBigIndex j;
427            const double * columnScale = model->columnScale();
428            for (j=columnStart[iColumn];
429                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
430              int iRow = row[j];
431              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
432            }
433            value *= columnScale[iColumn];
434            if (fabs(value)>zeroTolerance) {
435              index[numberNonZero++]=iColumn;
436              array[iColumn]=-value;
437            }
438          }
439        } else if (scalar==1.0) {
440          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
441            double value = 0.0;
442            CoinBigIndex j;
443            const double * columnScale = model->columnScale();
444            for (j=columnStart[iColumn];
445                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
446              int iRow = row[j];
447              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
448            }
449            value *= columnScale[iColumn];
450            if (fabs(value)>zeroTolerance) {
451              index[numberNonZero++]=iColumn;
452              array[iColumn]=value;
453            }
454          }
455        } else {
456          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
457            double value = 0.0;
458            CoinBigIndex j;
459            const double * columnScale = model->columnScale();
460            for (j=columnStart[iColumn];
461                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
462              int iRow = row[j];
463              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
464            }
465            value *= scalar*columnScale[iColumn];
466            if (fabs(value)>zeroTolerance) {
467              index[numberNonZero++]=iColumn;
468              array[iColumn]=value;
469            }
470          }
471        }
472      }
473    }
474    columnArray->setNumElements(numberNonZero);
475    y->setNumElements(0);
476  } else {
477    // do by row
478    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
479  }
480  if (packed)
481    columnArray->setPackedMode(true);
482  if (0) {
483    columnArray->checkClean();
484    int numberNonZero=columnArray->getNumElements();;
485    int * index = columnArray->getIndices();
486    double * array = columnArray->denseVector();
487    int i;
488    for (i=0;i<numberNonZero;i++) {
489      int j=index[i];
490      double value;
491      if (packed)
492        value=array[i];
493      else
494        value=array[j];
495      printf("Ti %d %d %g\n",i,j,value);
496    }
497  }
498}
499/* Return <code>x * scalar * A + y</code> in <code>z</code>.
500   Note - If x packed mode - then z packed mode
501   This does by column and knows no gaps
502   Squashes small elements and knows about ClpSimplex */
503void 
504ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
505                                        const CoinIndexedVector * rowArray,
506                                        CoinIndexedVector * y,
507                                        CoinIndexedVector * columnArray) const
508{
509  double * pi = rowArray->denseVector();
510  int numberNonZero=0;
511  int * index = columnArray->getIndices();
512  double * array = columnArray->denseVector();
513  int numberInRowArray = rowArray->getNumElements();
514  // maybe I need one in OsiSimplex
515  double zeroTolerance = model->factorization()->zeroTolerance();
516  bool packed = rowArray->packedMode();
517  // do by column
518  int iColumn;
519  // get matrix data pointers
520  const int * row = matrix_->getIndices();
521  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
522  const double * elementByColumn = matrix_->getElements();
523  const double * rowScale = model->rowScale();
524  assert (!y->getNumElements());
525  assert (numberActiveColumns_>0);
526  if (packed) {
527    // need to expand pi into y
528    assert(y->capacity()>=model->numberRows());
529    double * piOld = pi;
530    pi = y->denseVector();
531    const int * whichRow = rowArray->getIndices();
532    int i;
533    if (!rowScale) {
534      // modify pi so can collapse to one loop
535      for (i=0;i<numberInRowArray;i++) {
536        int iRow = whichRow[i];
537        pi[iRow]=scalar*piOld[i];
538      }
539      double value = 0.0;
540      CoinBigIndex j;
541      CoinBigIndex end = columnStart[1];
542      for (j=columnStart[0]; j<end;j++) {
543        int iRow = row[j];
544        value += pi[iRow]*elementByColumn[j];
545      }
546      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
547        CoinBigIndex start = end;
548        end = columnStart[iColumn+2];
549        if (fabs(value)>zeroTolerance) {
550          array[numberNonZero]=value;
551          index[numberNonZero++]=iColumn;
552        }
553        value = 0.0;
554        for (j=start; j<end;j++) {
555          int iRow = row[j];
556          value += pi[iRow]*elementByColumn[j];
557        }
558      }
559      if (fabs(value)>zeroTolerance) {
560        array[numberNonZero]=value;
561        index[numberNonZero++]=iColumn;
562      }
563    } else {
564      // scaled
565      // modify pi so can collapse to one loop
566      for (i=0;i<numberInRowArray;i++) {
567        int iRow = whichRow[i];
568        pi[iRow]=scalar*piOld[i]*rowScale[iRow];
569      }
570      const double * columnScale = model->columnScale();
571      double value = 0.0;
572      double scale=columnScale[0];
573      CoinBigIndex j;
574      CoinBigIndex end = columnStart[1];
575      for (j=columnStart[0]; j<end;j++) {
576        int iRow = row[j];
577        value += pi[iRow]*elementByColumn[j];
578      }
579      for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
580        value *= scale;
581        CoinBigIndex start = end;
582        scale = columnScale[iColumn+1];
583        end = columnStart[iColumn+2];
584        if (fabs(value)>zeroTolerance) {
585          array[numberNonZero]=value;
586          index[numberNonZero++]=iColumn;
587        }
588        value = 0.0;
589        for (j=start; j<end;j++) {
590          int iRow = row[j];
591          value += pi[iRow]*elementByColumn[j];
592        }
593      }
594      value *= scale;
595      if (fabs(value)>zeroTolerance) {
596        array[numberNonZero]=value;
597        index[numberNonZero++]=iColumn;
598      }
599    }
600    // zero out
601    for (i=0;i<numberInRowArray;i++) {
602      int iRow = whichRow[i];
603      pi[iRow]=0.0;
604    }
605  } else {
606    if (!rowScale) {
607      if (scalar==-1.0) {
608        double value = 0.0;
609        CoinBigIndex j;
610        CoinBigIndex end = columnStart[1];
611        for (j=columnStart[0]; j<end;j++) {
612          int iRow = row[j];
613          value += pi[iRow]*elementByColumn[j];
614        }
615        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
616          CoinBigIndex start = end;
617          end = columnStart[iColumn+2];
618          if (fabs(value)>zeroTolerance) {
619            array[iColumn]=-value;
620            index[numberNonZero++]=iColumn;
621          }
622          value = 0.0;
623          for (j=start; j<end;j++) {
624            int iRow = row[j];
625            value += pi[iRow]*elementByColumn[j];
626          }
627        }
628        if (fabs(value)>zeroTolerance) {
629          array[iColumn]=-value;
630          index[numberNonZero++]=iColumn;
631        }
632      } else if (scalar==1.0) {
633        double value = 0.0;
634        CoinBigIndex j;
635        CoinBigIndex end = columnStart[1];
636        for (j=columnStart[0]; j<end;j++) {
637          int iRow = row[j];
638          value += pi[iRow]*elementByColumn[j];
639        }
640        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
641          CoinBigIndex start = end;
642          end = columnStart[iColumn+2];
643          if (fabs(value)>zeroTolerance) {
644            array[iColumn]=value;
645            index[numberNonZero++]=iColumn;
646          }
647          value = 0.0;
648          for (j=start; j<end;j++) {
649            int iRow = row[j];
650            value += pi[iRow]*elementByColumn[j];
651          }
652        }
653        if (fabs(value)>zeroTolerance) {
654          array[iColumn]=value;
655          index[numberNonZero++]=iColumn;
656        }
657      } else {
658        double value = 0.0;
659        CoinBigIndex j;
660        CoinBigIndex end = columnStart[1];
661        for (j=columnStart[0]; j<end;j++) {
662          int iRow = row[j];
663          value += pi[iRow]*elementByColumn[j];
664        }
665        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
666          value *= scalar;
667          CoinBigIndex start = end;
668          end = columnStart[iColumn+2];
669          if (fabs(value)>zeroTolerance) {
670            array[iColumn]=value;
671            index[numberNonZero++]=iColumn;
672          }
673          value = 0.0;
674          for (j=start; j<end;j++) {
675            int iRow = row[j];
676            value += pi[iRow]*elementByColumn[j];
677          }
678        }
679        value *= scalar;
680        if (fabs(value)>zeroTolerance) {
681          array[iColumn]=value;
682          index[numberNonZero++]=iColumn;
683        }
684      }
685    } else {
686      // unscaled
687      if (scalar==-1.0) {
688        const double * columnScale = model->columnScale();
689        double value = 0.0;
690        double scale=columnScale[0];
691        CoinBigIndex j;
692        CoinBigIndex end = columnStart[1];
693        for (j=columnStart[0]; j<end;j++) {
694          int iRow = row[j];
695          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
696        }
697        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
698          value *= scale;
699          CoinBigIndex start = end;
700          end = columnStart[iColumn+2];
701          scale = columnScale[iColumn+1];
702          if (fabs(value)>zeroTolerance) {
703            array[iColumn]=-value;
704            index[numberNonZero++]=iColumn;
705          }
706          value = 0.0;
707          for (j=start; j<end;j++) {
708            int iRow = row[j];
709            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
710          }
711        }
712        value *= scale;
713        if (fabs(value)>zeroTolerance) {
714          array[iColumn]=-value;
715          index[numberNonZero++]=iColumn;
716        }
717      } else if (scalar==1.0) {
718        const double * columnScale = model->columnScale();
719        double value = 0.0;
720        double scale=columnScale[0];
721        CoinBigIndex j;
722        CoinBigIndex end = columnStart[1];
723        for (j=columnStart[0]; j<end;j++) {
724          int iRow = row[j];
725          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
726        }
727        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
728          value *= scale;
729          CoinBigIndex start = end;
730          end = columnStart[iColumn+2];
731          scale = columnScale[iColumn+1];
732          if (fabs(value)>zeroTolerance) {
733            array[iColumn]=value;
734            index[numberNonZero++]=iColumn;
735          }
736          value = 0.0;
737          for (j=start; j<end;j++) {
738            int iRow = row[j];
739            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
740          }
741        }
742        value *= scale;
743        if (fabs(value)>zeroTolerance) {
744          array[iColumn]=value;
745          index[numberNonZero++]=iColumn;
746        }
747      } else {
748        const double * columnScale = model->columnScale();
749        double value = 0.0;
750        double scale=columnScale[0]*scalar;
751        CoinBigIndex j;
752        CoinBigIndex end = columnStart[1];
753        for (j=columnStart[0]; j<end;j++) {
754          int iRow = row[j];
755          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
756        }
757        for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
758          value *= scale;
759          CoinBigIndex start = end;
760          end = columnStart[iColumn+2];
761          scale = columnScale[iColumn+1]*scalar;
762          if (fabs(value)>zeroTolerance) {
763            array[iColumn]=value;
764            index[numberNonZero++]=iColumn;
765          }
766          value = 0.0;
767          for (j=start; j<end;j++) {
768            int iRow = row[j];
769            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
770          }
771        }
772        value *= scale;
773        if (fabs(value)>zeroTolerance) {
774          array[iColumn]=value;
775          index[numberNonZero++]=iColumn;
776        }
777      }
778    }
779  }
780  columnArray->setNumElements(numberNonZero);
781  y->setNumElements(0);
782  if (packed)
783    columnArray->setPackedMode(true);
784}
785/* Return <code>x * scalar * A in <code>z</code>.
786   Note - this version when x packed mode and so will be z
787   This does by column and knows no gaps and knows y empty
788   Squashes small elements and knows about ClpSimplex */
789void 
790ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
791                                        const CoinIndexedVector * rowArray,
792                                        CoinIndexedVector * columnArray) const
793{
794}
795/* Return <code>x * A + y</code> in <code>z</code>.
796        Squashes small elements and knows about ClpSimplex */
797void 
798ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
799                              const CoinIndexedVector * rowArray,
800                              CoinIndexedVector * y,
801                              CoinIndexedVector * columnArray) const
802{
803  columnArray->clear();
804  double * pi = rowArray->denseVector();
805  int numberNonZero=0;
806  int * index = columnArray->getIndices();
807  double * array = columnArray->denseVector();
808  int numberInRowArray = rowArray->getNumElements();
809  // maybe I need one in OsiSimplex
810  double zeroTolerance = model->factorization()->zeroTolerance();
811  const int * column = getIndices();
812  const CoinBigIndex * rowStart = getVectorStarts();
813  const double * element = getElements();
814  const int * whichRow = rowArray->getIndices();
815  bool packed = rowArray->packedMode();
816  if (numberInRowArray>2) {
817    // do by rows
818    // ** Row copy is already scaled
819    int iRow;
820    int i;
821    int numberOriginal = 0;
822    if (packed) {
823      numberNonZero=0;
824      // and set up mark as char array
825      char * marked = (char *) (index+columnArray->capacity());
826      double * array2 = y->denseVector();
827#ifdef CLP_DEBUG
828      for (i=0;i<numberActiveColumns_;i++) {
829        assert(!marked[i]);
830        assert(!array2[i]);
831      }
832#endif     
833      for (i=0;i<numberInRowArray;i++) {
834        iRow = whichRow[i]; 
835        double value = pi[i]*scalar;
836        CoinBigIndex j;
837        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
838          int iColumn = column[j];
839          if (!marked[iColumn]) {
840            marked[iColumn]=1;
841            index[numberNonZero++]=iColumn;
842          }
843          array2[iColumn] += value*element[j];
844        }
845      }
846      // get rid of tiny values and zero out marked
847      numberOriginal=numberNonZero;
848      numberNonZero=0;
849      for (i=0;i<numberOriginal;i++) {
850        int iColumn = index[i];
851        if (marked[iColumn]) {
852          double value = array2[iColumn];
853          array2[iColumn]=0.0;
854          marked[iColumn]=0;
855          if (fabs(value)>zeroTolerance) {
856            array[numberNonZero]=value;
857            index[numberNonZero++]=iColumn;
858          }
859        }
860      }
861    } else {
862      double * markVector = y->denseVector();
863      numberNonZero=0;
864      // and set up mark as char array
865      char * marked = (char *) markVector;
866      for (i=0;i<numberOriginal;i++) {
867        int iColumn = index[i];
868        marked[iColumn]=0;
869      }
870     
871      for (i=0;i<numberInRowArray;i++) {
872        iRow = whichRow[i]; 
873        double value = pi[iRow]*scalar;
874        CoinBigIndex j;
875        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
876          int iColumn = column[j];
877          if (!marked[iColumn]) {
878            marked[iColumn]=1;
879            index[numberNonZero++]=iColumn;
880          }
881          array[iColumn] += value*element[j];
882        }
883      }
884      // get rid of tiny values and zero out marked
885      numberOriginal=numberNonZero;
886      numberNonZero=0;
887      for (i=0;i<numberOriginal;i++) {
888        int iColumn = index[i];
889        marked[iColumn]=0;
890        if (fabs(array[iColumn])>zeroTolerance) {
891          index[numberNonZero++]=iColumn;
892        } else {
893          array[iColumn]=0.0;
894        }
895      }
896    }
897  } else if (numberInRowArray==2) {
898    // do by rows when two rows
899    int numberOriginal;
900    int i;
901    CoinBigIndex j;
902    numberNonZero=0;
903
904    double value;
905    if (packed) {
906      int iRow0 = whichRow[0]; 
907      int iRow1 = whichRow[1]; 
908      double pi0 = pi[0];
909      double pi1 = pi[1];
910      if (rowStart[iRow0+1]-rowStart[iRow0]>
911          rowStart[iRow1+1]-rowStart[iRow1]) {
912        // do one with fewer first
913        iRow0=iRow1;
914        iRow1=whichRow[0];
915        pi0=pi1;
916        pi1=pi[0];
917      }
918      // and set up mark as char array
919      char * marked = (char *) (index+columnArray->capacity());
920      int * lookup = y->getIndices();
921      value = pi0*scalar;
922      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
923        int iColumn = column[j];
924        double value2 = value*element[j];
925        array[numberNonZero] = value2;
926        marked[iColumn]=1;
927        lookup[iColumn]=numberNonZero;
928        index[numberNonZero++]=iColumn;
929      }
930      numberOriginal = numberNonZero;
931      value = pi1*scalar;
932      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
933        int iColumn = column[j];
934        double value2 = value*element[j];
935        // I am assuming no zeros in matrix
936        if (marked[iColumn]) {
937          int iLookup = lookup[iColumn];
938          array[iLookup] += value2;
939        } else {
940          if (fabs(value2)>zeroTolerance) {
941            array[numberNonZero] = value2;
942            index[numberNonZero++]=iColumn;
943          }
944        }
945      }
946      // get rid of tiny values and zero out marked
947      int nDelete=0;
948      for (i=0;i<numberOriginal;i++) {
949        int iColumn = index[i];
950        marked[iColumn]=0;
951        if (fabs(array[i])<=zeroTolerance) 
952          nDelete++;
953      }
954      if (nDelete) {
955        numberOriginal=numberNonZero;
956        numberNonZero=0;
957        for (i=0;i<numberOriginal;i++) {
958          int iColumn = index[i];
959          double value = array[i];
960          array[i]=0.0;
961          if (fabs(value)>zeroTolerance) {
962            array[numberNonZero]=value;
963            index[numberNonZero++]=iColumn;
964          }
965        }
966      }
967    } else {
968      int iRow = whichRow[0]; 
969      value = pi[iRow]*scalar;
970      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
971        int iColumn = column[j];
972        double value2 = value*element[j];
973        index[numberNonZero++]=iColumn;
974        array[iColumn] = value2;
975      }
976      iRow = whichRow[1]; 
977      value = pi[iRow]*scalar;
978      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
979        int iColumn = column[j];
980        double value2 = value*element[j];
981        // I am assuming no zeros in matrix
982        if (array[iColumn])
983          value2 += array[iColumn];
984        else
985          index[numberNonZero++]=iColumn;
986        array[iColumn] = value2;
987      }
988      // get rid of tiny values and zero out marked
989      numberOriginal=numberNonZero;
990      numberNonZero=0;
991      for (i=0;i<numberOriginal;i++) {
992        int iColumn = index[i];
993        if (fabs(array[iColumn])>zeroTolerance) {
994          index[numberNonZero++]=iColumn;
995        } else {
996          array[iColumn]=0.0;
997        }
998      }
999    }
1000  } else if (numberInRowArray==1) {
1001    // Just one row
1002    int iRow=rowArray->getIndices()[0];
1003    numberNonZero=0;
1004    CoinBigIndex j;
1005    if (packed) {
1006      double value = pi[0]*scalar;
1007      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1008        int iColumn = column[j];
1009        double value2 = value*element[j];
1010        if (fabs(value2)>zeroTolerance) {
1011          array[numberNonZero] = value2;
1012          index[numberNonZero++]=iColumn;
1013        }
1014      }
1015    } else {
1016      double value = pi[iRow]*scalar;
1017      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1018        int iColumn = column[j];
1019        double value2 = value*element[j];
1020        if (fabs(value2)>zeroTolerance) {
1021          index[numberNonZero++]=iColumn;
1022          array[iColumn] = value2;
1023        }
1024      }
1025    }
1026  }
1027  columnArray->setNumElements(numberNonZero);
1028  y->setNumElements(0);
1029}
1030/* Return <code>x *A in <code>z</code> but
1031   just for indices in y.
1032   Squashes small elements and knows about ClpSimplex */
1033void 
1034ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
1035                              const CoinIndexedVector * rowArray,
1036                              const CoinIndexedVector * y,
1037                              CoinIndexedVector * columnArray) const
1038{
1039  columnArray->clear();
1040  double * pi = rowArray->denseVector();
1041  double * array = columnArray->denseVector();
1042  int jColumn;
1043  // get matrix data pointers
1044  const int * row = matrix_->getIndices();
1045  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1046  const int * columnLength = matrix_->getVectorLengths(); 
1047  const double * elementByColumn = matrix_->getElements();
1048  const double * rowScale = model->rowScale();
1049  int numberToDo = y->getNumElements();
1050  const int * which = y->getIndices();
1051  assert (!rowArray->packedMode());
1052  columnArray->setPacked();
1053  if (!hasGaps_&&numberToDo>5) {
1054    // no gaps and a reasonable number
1055    if (!rowScale) {
1056      int iColumn = which[0];
1057      double value = 0.0;
1058      CoinBigIndex j;
1059      for (j=columnStart[iColumn];
1060           j<columnStart[iColumn+1];j++) {
1061        int iRow = row[j];
1062        value += pi[iRow]*elementByColumn[j];
1063      }
1064      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
1065        int iColumn = which[jColumn+1];
1066        CoinBigIndex start = columnStart[iColumn];
1067        CoinBigIndex end = columnStart[iColumn+1];
1068        array[jColumn]=value;
1069        value = 0.0;
1070        for (j=start; j<end;j++) {
1071          int iRow = row[j];
1072          value += pi[iRow]*elementByColumn[j];
1073        }
1074      }
1075      array[jColumn]=value;
1076    } else {
1077      // scaled
1078      const double * columnScale = model->columnScale();
1079      int iColumn = which[0];
1080      double value = 0.0;
1081      double scale=columnScale[iColumn];
1082      CoinBigIndex j;
1083      for (j=columnStart[iColumn];
1084           j<columnStart[iColumn+1];j++) {
1085        int iRow = row[j];
1086        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1087      }
1088      for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
1089        int iColumn = which[jColumn+1];
1090        value *= scale;
1091        scale = columnScale[iColumn];
1092        CoinBigIndex start = columnStart[iColumn];
1093        CoinBigIndex end = columnStart[iColumn+1];
1094        array[jColumn]=value;
1095        value = 0.0;
1096        for (j=start; j<end;j++) {
1097          int iRow = row[j];
1098          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1099        }
1100      }
1101      value *= scale;
1102      array[jColumn]=value;
1103    }
1104  } else {
1105    // gaps
1106    if (!rowScale) {
1107      for (jColumn=0;jColumn<numberToDo;jColumn++) {
1108        int iColumn = which[jColumn];
1109        double value = 0.0;
1110        CoinBigIndex j;
1111        for (j=columnStart[iColumn];
1112             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1113          int iRow = row[j];
1114          value += pi[iRow]*elementByColumn[j];
1115        }
1116        array[jColumn]=value;
1117      }
1118    } else {
1119      // scaled
1120      const double * columnScale = model->columnScale();
1121      for (jColumn=0;jColumn<numberToDo;jColumn++) {
1122        int iColumn = which[jColumn];
1123        double value = 0.0;
1124        CoinBigIndex j;
1125        for (j=columnStart[iColumn];
1126             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1127          int iRow = row[j];
1128          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
1129        }
1130        value *= columnScale[iColumn];
1131        array[jColumn]=value;
1132      }
1133    }
1134  }
1135}
1136/* If element NULL returns number of elements in column part of basis,
1137   If not NULL fills in as well */
1138CoinBigIndex
1139ClpPackedMatrix::fillBasis(ClpSimplex * model,
1140                           const int * whichColumn, 
1141                           int numberBasic,
1142                           int & numberColumnBasic,
1143                           int * indexRowU, int * indexColumnU,
1144                           double * elementU)
1145{
1146  const int * columnLength = matrix_->getVectorLengths(); 
1147  int i;
1148  CoinBigIndex numberElements=0;
1149  if (elementU!=NULL) {
1150    // fill
1151    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1152    const double * rowScale = model->rowScale();
1153    const int * row = matrix_->getIndices();
1154    const double * elementByColumn = matrix_->getElements();
1155    if (!zeroElements_) {
1156      if (!rowScale) {
1157        // no scaling
1158        for (i=0;i<numberColumnBasic;i++) {
1159          int iColumn = whichColumn[i];
1160          CoinBigIndex j;
1161          for (j=columnStart[iColumn];
1162               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1163            indexRowU[numberElements]=row[j];
1164            indexColumnU[numberElements]=numberBasic;
1165            elementU[numberElements++]=elementByColumn[j];
1166          }
1167          numberBasic++;
1168        }
1169      } else {
1170        // scaling
1171        const double * columnScale = model->columnScale();
1172        for (i=0;i<numberColumnBasic;i++) {
1173          int iColumn = whichColumn[i];
1174          CoinBigIndex j;
1175          double scale = columnScale[iColumn];
1176          for (j=columnStart[iColumn];
1177               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1178            int iRow = row[j];
1179            indexRowU[numberElements]=iRow;
1180            indexColumnU[numberElements]=numberBasic;
1181            elementU[numberElements++]=
1182              elementByColumn[j]*scale*rowScale[iRow];
1183          }
1184          numberBasic++;
1185        }
1186      }
1187    } else {
1188      // there are zero elements so need to look more closely
1189      if (!rowScale) {
1190        // no scaling
1191        for (i=0;i<numberColumnBasic;i++) {
1192          int iColumn = whichColumn[i];
1193          CoinBigIndex j;
1194          for (j=columnStart[iColumn];
1195               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1196            double value = elementByColumn[j];
1197            if (value) {
1198              indexRowU[numberElements]=row[j];
1199              indexColumnU[numberElements]=numberBasic;
1200              elementU[numberElements++]=value;
1201            }
1202          }
1203          numberBasic++;
1204        }
1205      } else {
1206        // scaling
1207        const double * columnScale = model->columnScale();
1208        for (i=0;i<numberColumnBasic;i++) {
1209          int iColumn = whichColumn[i];
1210          CoinBigIndex j;
1211          double scale = columnScale[iColumn];
1212          for (j=columnStart[iColumn];
1213               j<columnStart[iColumn]+columnLength[i];j++) {
1214            double value = elementByColumn[j];
1215            if (value) {
1216              int iRow = row[j];
1217              indexRowU[numberElements]=iRow;
1218              indexColumnU[numberElements]=numberBasic;
1219              elementU[numberElements++]=value*scale*rowScale[iRow];
1220            }
1221          }
1222          numberBasic++;
1223        }
1224      }
1225    }
1226  } else {
1227    // just count - can be over so ignore zero problem
1228    for (i=0;i<numberColumnBasic;i++) {
1229      int iColumn = whichColumn[i];
1230      numberElements += columnLength[iColumn];
1231    }
1232  }
1233  return numberElements;
1234}
1235// Creates scales for column copy (rowCopy in model may be modified)
1236int 
1237ClpPackedMatrix::scale(ClpModel * model) const 
1238{
1239  int numberRows = model->numberRows();
1240  int numberColumns = matrix_->getNumCols();
1241  // If empty - return as sanityCheck will trap
1242  if (!numberRows||!numberColumns)
1243    return 1;
1244  ClpMatrixBase * rowCopyBase=model->rowCopy();
1245  if (!rowCopyBase) {
1246    // temporary copy
1247    rowCopyBase = reverseOrderedCopy();
1248  }
1249#ifndef NO_RTTI
1250  ClpPackedMatrix* rowCopy =
1251    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
1252  // Make sure it is really a ClpPackedMatrix
1253  assert (rowCopy!=NULL);
1254#else
1255  ClpPackedMatrix* rowCopy =
1256    static_cast< ClpPackedMatrix*>(rowCopyBase);
1257#endif
1258
1259  const int * column = rowCopy->getIndices();
1260  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1261  const double * element = rowCopy->getElements();
1262  double * rowScale = new double [numberRows];
1263  double * columnScale = new double [numberColumns];
1264  // we are going to mark bits we are interested in
1265  char * usefulRow = new char [numberRows];
1266  char * usefulColumn = new char [numberColumns];
1267  double * rowLower = model->rowLower();
1268  double * rowUpper = model->rowUpper();
1269  double * columnLower = model->columnLower();
1270  double * columnUpper = model->columnUpper();
1271  int iColumn, iRow;
1272  // mark free rows
1273  for (iRow=0;iRow<numberRows;iRow++) {
1274    usefulRow[iRow]=0;
1275    if (rowUpper[iRow]<1.0e20||
1276        rowLower[iRow]>-1.0e20)
1277      usefulRow[iRow]=1;
1278  }
1279  // mark empty and fixed columns
1280  // also see if worth scaling
1281  assert (model->scalingFlag()<4); // dynamic not implemented
1282  double largest=0.0;
1283  double smallest=1.0e50;
1284  // get matrix data pointers
1285  const int * row = matrix_->getIndices();
1286  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1287  const int * columnLength = matrix_->getVectorLengths(); 
1288  const double * elementByColumn = matrix_->getElements();
1289  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1290    CoinBigIndex j;
1291    char useful=0;
1292    if (columnUpper[iColumn]>
1293        columnLower[iColumn]+1.0e-9) {
1294      for (j=columnStart[iColumn];
1295           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1296        iRow=row[j];
1297        if(elementByColumn[j]&&usefulRow[iRow]) {
1298          useful=1;
1299          largest = CoinMax(largest,fabs(elementByColumn[j]));
1300          smallest = CoinMin(smallest,fabs(elementByColumn[j]));
1301        }
1302      }
1303    }
1304    usefulColumn[iColumn]=useful;
1305  }
1306  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
1307    <<smallest<<largest
1308    <<CoinMessageEol;
1309  if (smallest>=0.5&&largest<=2.0) {
1310    // don't bother scaling
1311    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
1312      <<CoinMessageEol;
1313    delete [] rowScale;
1314    delete [] usefulRow;
1315    delete [] columnScale;
1316    delete [] usefulColumn;
1317    if (!model->rowCopy()) 
1318      delete rowCopyBase; // get rid of temporary row copy
1319    return 1;
1320  } else {
1321      // need to scale
1322    int scalingMethod = model->scalingFlag();
1323    if (scalingMethod==3) {
1324      // Choose between 1 and 2
1325      if (smallest<1.0e-5||smallest*largest<1.0)
1326        scalingMethod=1;
1327      else
1328        scalingMethod=2;
1329    }
1330    // and see if there any empty rows
1331    for (iRow=0;iRow<numberRows;iRow++) {
1332      if (usefulRow[iRow]) {
1333        CoinBigIndex j;
1334        int useful=0;
1335        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1336          int iColumn = column[j];
1337          if (usefulColumn[iColumn]) {
1338            useful=1;
1339            break;
1340          }
1341        }
1342        usefulRow[iRow]=useful;
1343      }
1344    }
1345    ClpFillN ( rowScale, numberRows,1.0);
1346    ClpFillN ( columnScale, numberColumns,1.0);
1347    double overallLargest=-1.0e-30;
1348    double overallSmallest=1.0e30;
1349    if (scalingMethod==1) {
1350      // Maximum in each row
1351      for (iRow=0;iRow<numberRows;iRow++) {
1352        if (usefulRow[iRow]) {
1353          CoinBigIndex j;
1354          largest=1.0e-20;
1355          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1356            int iColumn = column[j];
1357            if (usefulColumn[iColumn]) {
1358              double value = fabs(element[j]);
1359              largest = CoinMax(largest,value);
1360            }
1361          }
1362          rowScale[iRow]=1.0/largest;
1363          overallLargest = CoinMax(overallLargest,largest);
1364          overallSmallest = CoinMin(overallSmallest,largest);
1365        }
1366      }
1367    } else {
1368      assert(scalingMethod==2);
1369      int numberPass=3;
1370#ifdef USE_OBJECTIVE
1371      // This will be used to help get scale factors
1372      double * objective = new double[numberColumns];
1373      memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
1374      double objScale=1.0;
1375#endif
1376      while (numberPass) {
1377        overallLargest=0.0;
1378        overallSmallest=1.0e50;
1379        numberPass--;
1380        // Geometric mean on row scales
1381        for (iRow=0;iRow<numberRows;iRow++) {
1382          if (usefulRow[iRow]) {
1383            CoinBigIndex j;
1384            largest=1.0e-20;
1385            smallest=1.0e50;
1386            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
1387              int iColumn = column[j];
1388              if (usefulColumn[iColumn]) {
1389                double value = fabs(element[j]);
1390                // Don't bother with tiny elements
1391                if (value>1.0e-30) {
1392                  value *= columnScale[iColumn];
1393                  largest = CoinMax(largest,value);
1394                  smallest = CoinMin(smallest,value);
1395                }
1396              }
1397            }
1398            rowScale[iRow]=1.0/sqrt(smallest*largest);
1399            overallLargest = CoinMax(largest*rowScale[iRow],overallLargest);
1400            overallSmallest = CoinMin(smallest*rowScale[iRow],overallSmallest);
1401          }
1402        }
1403#ifdef USE_OBJECTIVE
1404        largest=1.0e-20;
1405        smallest=1.0e50;
1406        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1407          if (usefulColumn[iColumn]) {
1408            double value = fabs(objective[iColumn]);
1409            // Don't bother with tiny elements
1410            if (value>1.0e-30) {
1411              value *= columnScale[iColumn];
1412              largest = CoinMax(largest,value);
1413              smallest = CoinMin(smallest,value);
1414            }
1415          }
1416        }
1417        objScale=1.0/sqrt(smallest*largest);
1418#endif
1419        model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
1420          <<overallSmallest
1421          <<overallLargest
1422          <<CoinMessageEol;
1423        // skip last column round
1424        if (numberPass==1)
1425          break;
1426        // Geometric mean on column scales
1427        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1428          if (usefulColumn[iColumn]) {
1429            CoinBigIndex j;
1430            largest=1.0e-20;
1431            smallest=1.0e50;
1432            for (j=columnStart[iColumn];
1433                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1434              iRow=row[j];
1435              double value = fabs(elementByColumn[j]);
1436              // Don't bother with tiny elements
1437              if (value>1.0e-30&&usefulRow[iRow]) {
1438                value *= rowScale[iRow];
1439                largest = CoinMax(largest,value);
1440                smallest = CoinMin(smallest,value);
1441              }
1442            }
1443#ifdef USE_OBJECTIVE
1444            if (fabs(objective[iColumn])>1.0e-30) {
1445              double value = fabs(objective[iColumn])*objScale;
1446              largest = CoinMax(largest,value);
1447              smallest = CoinMin(smallest,value);
1448            }
1449#endif
1450            columnScale[iColumn]=1.0/sqrt(smallest*largest);
1451          }
1452        }
1453      }
1454#ifdef USE_OBJECTIVE
1455      delete [] objective;
1456      printf("obj scale %g - use it if you want\n",objScale);
1457#endif
1458    }
1459    // If ranges will make horrid then scale
1460    double tolerance = 5.0*model->primalTolerance();
1461    for (iRow=0;iRow<numberRows;iRow++) {
1462      if (usefulRow[iRow]) {
1463        double difference = rowUpper[iRow]-rowLower[iRow];
1464        double scaledDifference = difference*rowScale[iRow];
1465        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1466          // make gap larger
1467          rowScale[iRow] *= 1.0e-4/scaledDifference;
1468          //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
1469          // scaledDifference,difference*rowScale[iRow]);
1470        }
1471      }
1472    }
1473    // final pass to scale columns so largest is reasonable
1474    // See what smallest will be if largest is 1.0
1475    overallSmallest=1.0e50;
1476    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1477      if (usefulColumn[iColumn]) {
1478        CoinBigIndex j;
1479        largest=1.0e-20;
1480        smallest=1.0e50;
1481        for (j=columnStart[iColumn];
1482             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1483          iRow=row[j];
1484          if(elementByColumn[j]&&usefulRow[iRow]) {
1485            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1486            largest = CoinMax(largest,value);
1487            smallest = CoinMin(smallest,value);
1488          }
1489        }
1490        if (overallSmallest*largest>smallest)
1491          overallSmallest = smallest/largest;
1492      }
1493    }
1494    //#define RANDOMIZE
1495#ifdef RANDOMIZE
1496    // randomize by up to 10%
1497    for (iRow=0;iRow<numberRows;iRow++) {
1498      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1499      rowScale[iRow] *= (1.0+0.1*value);
1500    }
1501#endif
1502    overallLargest=1.0;
1503    if (overallSmallest<1.0e-1)
1504      overallLargest = 1.0/sqrt(overallSmallest);
1505    overallLargest = CoinMin(100.0,overallLargest);
1506    overallSmallest=1.0e50;
1507    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1508      if (usefulColumn[iColumn]) {
1509        CoinBigIndex j;
1510        largest=1.0e-20;
1511        smallest=1.0e50;
1512        for (j=columnStart[iColumn];
1513             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1514          iRow=row[j];
1515          if(elementByColumn[j]&&usefulRow[iRow]) {
1516            double value = fabs(elementByColumn[j]*rowScale[iRow]);
1517            largest = CoinMax(largest,value);
1518            smallest = CoinMin(smallest,value);
1519          }
1520        }
1521        columnScale[iColumn]=overallLargest/largest;
1522#ifdef RANDOMIZE
1523        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
1524        columnScale[iColumn] *= (1.0+0.1*value);
1525#endif
1526        double difference = columnUpper[iColumn]-columnLower[iColumn];
1527        double scaledDifference = difference*columnScale[iColumn];
1528        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
1529          // make gap larger
1530          columnScale[iColumn] *= 1.0e-4/scaledDifference;
1531          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
1532          // scaledDifference,difference*columnScale[iColumn]);
1533        }
1534        overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
1535      }
1536    }
1537    model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
1538      <<overallSmallest
1539      <<overallLargest
1540      <<CoinMessageEol;
1541    delete [] usefulRow;
1542    delete [] usefulColumn;
1543    // If quadratic then make symmetric
1544    ClpObjective * obj = model->objectiveAsObject();
1545#ifndef NO_RTTI
1546    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
1547#else
1548    ClpQuadraticObjective * quadraticObj = NULL;
1549    if (obj->type()==2)
1550      quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
1551#endif
1552    if (quadraticObj) {
1553      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1554      int numberXColumns = quadratic->getNumCols();
1555      if (numberXColumns<numberColumns) {
1556        // we assume symmetric
1557        int numberQuadraticColumns=0;
1558        int i;
1559        //const int * columnQuadratic = quadratic->getIndices();
1560        //const int * columnQuadraticStart = quadratic->getVectorStarts();
1561        const int * columnQuadraticLength = quadratic->getVectorLengths();
1562        for (i=0;i<numberXColumns;i++) {
1563          int length=columnQuadraticLength[i];
1564#ifndef CORRECT_COLUMN_COUNTS
1565          length=1;
1566#endif
1567          if (length)
1568            numberQuadraticColumns++;
1569        }
1570        int numberXRows = numberRows-numberQuadraticColumns;
1571        numberQuadraticColumns=0;
1572        for (i=0;i<numberXColumns;i++) { 
1573          int length=columnQuadraticLength[i];
1574#ifndef CORRECT_COLUMN_COUNTS
1575          length=1;
1576#endif
1577          if (length) {
1578            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
1579            numberQuadraticColumns++;
1580          }
1581        }   
1582        int numberQuadraticRows=0;
1583        for (i=0;i<numberXRows;i++) {
1584          // See if any in row quadratic
1585          CoinBigIndex j;
1586          int numberQ=0;
1587          for (j=rowStart[i];j<rowStart[i+1];j++) {
1588            int iColumn = column[j];
1589            if (columnQuadraticLength[iColumn])
1590              numberQ++;
1591          }
1592#ifndef CORRECT_ROW_COUNTS
1593          numberQ=1;
1594#endif
1595          if (numberQ) {
1596            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
1597            numberQuadraticRows++;
1598          }
1599        }
1600        // and make sure Sj okay
1601        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
1602          CoinBigIndex j=columnStart[iColumn];
1603          assert(columnLength[iColumn]==1);
1604          int iRow=row[j];
1605          double value = fabs(elementByColumn[j]*rowScale[iRow]);
1606          columnScale[iColumn]=1.0/value;
1607        }
1608      }
1609    }
1610    model->setRowScale(rowScale);
1611    model->setColumnScale(columnScale);
1612    if (model->rowCopy()) {
1613      // need to replace row by row
1614      double * newElement = new double[numberColumns];
1615      // scale row copy
1616      for (iRow=0;iRow<numberRows;iRow++) {
1617        CoinBigIndex j;
1618        double scale = rowScale[iRow];
1619        const double * elementsInThisRow = element + rowStart[iRow];
1620        const int * columnsInThisRow = column + rowStart[iRow];
1621        int number = rowStart[iRow+1]-rowStart[iRow];
1622        assert (number<=numberColumns);
1623        for (j=0;j<number;j++) {
1624          int iColumn = columnsInThisRow[j];
1625          newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
1626        }
1627        rowCopy->replaceVector(iRow,number,newElement);
1628      }
1629      delete [] newElement;
1630    } else {
1631      // no row copy
1632      delete rowCopyBase;
1633    }
1634    return 0;
1635  }
1636}
1637/* Unpacks a column into an CoinIndexedvector
1638 */
1639void 
1640ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
1641                   int iColumn) const 
1642{
1643  const double * rowScale = model->rowScale();
1644  const int * row = matrix_->getIndices();
1645  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1646  const int * columnLength = matrix_->getVectorLengths(); 
1647  const double * elementByColumn = matrix_->getElements();
1648  CoinBigIndex i;
1649  if (!rowScale) {
1650    for (i=columnStart[iColumn];
1651         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1652      rowArray->add(row[i],elementByColumn[i]);
1653    }
1654  } else {
1655    // apply scaling
1656    double scale = model->columnScale()[iColumn];
1657    for (i=columnStart[iColumn];
1658         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1659      int iRow = row[i];
1660      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1661    }
1662  }
1663}
1664/* Unpacks a column into a CoinIndexedVector
1665** in packed format
1666Note that model is NOT const.  Bounds and objective could
1667be modified if doing column generation (just for this variable) */
1668void 
1669ClpPackedMatrix::unpackPacked(ClpSimplex * model,
1670                            CoinIndexedVector * rowArray,
1671                            int iColumn) const
1672{
1673  const double * rowScale = model->rowScale();
1674  const int * row = matrix_->getIndices();
1675  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1676  const int * columnLength = matrix_->getVectorLengths(); 
1677  const double * elementByColumn = matrix_->getElements();
1678  CoinBigIndex i;
1679  if (!rowScale) {
1680    CoinBigIndex j=columnStart[iColumn];
1681    rowArray->createPacked(columnLength[iColumn],
1682                           row+j,elementByColumn+j);
1683  } else {
1684    // apply scaling
1685    double scale = model->columnScale()[iColumn];
1686    int * index = rowArray->getIndices();
1687    double * array = rowArray->denseVector();
1688    int number = 0;
1689    for (i=columnStart[iColumn];
1690         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1691      int iRow = row[i];
1692      array[number]=elementByColumn[i]*scale*rowScale[iRow];
1693      index[number++]=iRow;
1694    }
1695    rowArray->setNumElements(number);
1696    rowArray->setPackedMode(true);
1697  }
1698}
1699/* Adds multiple of a column into an CoinIndexedvector
1700      You can use quickAdd to add to vector */
1701void 
1702ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1703                   int iColumn, double multiplier) const 
1704{
1705  const double * rowScale = model->rowScale();
1706  const int * row = matrix_->getIndices();
1707  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1708  const int * columnLength = matrix_->getVectorLengths(); 
1709  const double * elementByColumn = matrix_->getElements();
1710  CoinBigIndex i;
1711  if (!rowScale) {
1712    for (i=columnStart[iColumn];
1713         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1714      int iRow = row[i];
1715      rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
1716    }
1717  } else {
1718    // apply scaling
1719    double scale = model->columnScale()[iColumn]*multiplier;
1720    for (i=columnStart[iColumn];
1721         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1722      int iRow = row[i];
1723      rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
1724    }
1725  }
1726}
1727/* Adds multiple of a column into an array */
1728void 
1729ClpPackedMatrix::add(const ClpSimplex * model,double * array,
1730                    int iColumn, double multiplier) const
1731{
1732  const double * rowScale = model->rowScale();
1733  const int * row = matrix_->getIndices();
1734  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1735  const int * columnLength = matrix_->getVectorLengths(); 
1736  const double * elementByColumn = matrix_->getElements();
1737  CoinBigIndex i;
1738  if (!rowScale) {
1739    for (i=columnStart[iColumn];
1740         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1741      int iRow = row[i];
1742      array[iRow] += multiplier*elementByColumn[i];
1743    }
1744  } else {
1745    // apply scaling
1746    double scale = model->columnScale()[iColumn]*multiplier;
1747    for (i=columnStart[iColumn];
1748         i<columnStart[iColumn]+columnLength[iColumn];i++) {
1749      int iRow = row[i];
1750      array[iRow] += elementByColumn[i]*scale*rowScale[iRow];
1751    }
1752  }
1753}
1754/* Checks if all elements are in valid range.  Can just
1755   return true if you are not paranoid.  For Clp I will
1756   probably expect no zeros.  Code can modify matrix to get rid of
1757   small elements.
1758*/
1759bool 
1760ClpPackedMatrix::allElementsInRange(ClpModel * model,
1761                                    double smallest, double largest)
1762{
1763  int iColumn;
1764  CoinBigIndex numberLarge=0;;
1765  CoinBigIndex numberSmall=0;;
1766  CoinBigIndex numberDuplicate=0;;
1767  int firstBadColumn=-1;
1768  int firstBadRow=-1;
1769  double firstBadElement=0.0;
1770  // get matrix data pointers
1771  const int * row = matrix_->getIndices();
1772  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1773  const int * columnLength = matrix_->getVectorLengths(); 
1774  const double * elementByColumn = matrix_->getElements();
1775  int numberRows = matrix_->getNumRows();
1776  int numberColumns = matrix_->getNumCols();
1777  int * mark = new int [numberRows];
1778  int i;
1779  // Say no gaps
1780  hasGaps_=false;
1781  for (i=0;i<numberRows;i++)
1782    mark[i]=-1;
1783  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1784    CoinBigIndex j;
1785    CoinBigIndex start = columnStart[iColumn];
1786    CoinBigIndex end = start + columnLength[iColumn];
1787    if (end!=columnStart[iColumn+1])
1788      hasGaps_=true;
1789    for (j=start;j<end;j++) {
1790      double value = fabs(elementByColumn[j]);
1791      int iRow = row[j];
1792      if (iRow<0||iRow>=numberRows) {
1793#ifndef COIN_BIG_INDEX
1794        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1795#elif COIN_BIG_INDEX==1
1796        printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1797#else
1798        printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1799#endif
1800        return false;
1801      }
1802      if (mark[iRow]==-1) {
1803        mark[iRow]=j;
1804      } else {
1805        // duplicate
1806        numberDuplicate++;
1807      }
1808      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
1809      if (!value)
1810        zeroElements_ = true; // there are zero elements
1811      if (value<smallest) {
1812        numberSmall++;
1813      } else if (value>largest) {
1814        numberLarge++;
1815        if (firstBadColumn<0) {
1816          firstBadColumn=iColumn;
1817          firstBadRow=row[j];
1818          firstBadElement=elementByColumn[j];
1819        }
1820      }
1821    }
1822    //clear mark
1823    for (j=columnStart[iColumn];
1824         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1825      int iRow = row[j];
1826      mark[iRow]=-1;
1827    }
1828  }
1829  delete [] mark;
1830  if (numberLarge) {
1831    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
1832      <<numberLarge
1833      <<firstBadColumn<<firstBadRow<<firstBadElement
1834      <<CoinMessageEol;
1835    return false;
1836  }
1837  if (numberSmall) 
1838    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
1839      <<numberSmall
1840      <<CoinMessageEol;
1841  if (numberDuplicate) 
1842    model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
1843      <<numberDuplicate
1844      <<CoinMessageEol;
1845  if (numberDuplicate) 
1846    matrix_->eliminateDuplicates(smallest);
1847  else if (numberSmall) 
1848    matrix_->compress(smallest);
1849  // If smallest >0.0 then there can't be zero elements
1850  if (smallest>0.0)
1851    zeroElements_=false;
1852  return true;
1853}
1854/* Given positive integer weights for each row fills in sum of weights
1855   for each column (and slack).
1856   Returns weights vector
1857*/
1858CoinBigIndex * 
1859ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
1860{
1861  int numberRows = model->numberRows();
1862  int numberColumns = matrix_->getNumCols();
1863  int number = numberRows+numberColumns;
1864  CoinBigIndex * weights = new CoinBigIndex[number];
1865  // get matrix data pointers
1866  const int * row = matrix_->getIndices();
1867  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1868  const int * columnLength = matrix_->getVectorLengths(); 
1869  int i;
1870  for (i=0;i<numberColumns;i++) {
1871    CoinBigIndex j;
1872    CoinBigIndex count=0;
1873    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1874      int iRow=row[j];
1875      count += inputWeights[iRow];
1876    }
1877    weights[i]=count;
1878  }
1879  for (i=0;i<numberRows;i++) {
1880    weights[i+numberColumns]=inputWeights[i];
1881  }
1882  return weights;
1883}
1884/* Returns largest and smallest elements of both signs.
1885   Largest refers to largest absolute value.
1886*/
1887void 
1888ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1889                       double & smallestPositive, double & largestPositive)
1890{
1891  smallestNegative=-COIN_DBL_MAX;
1892  largestNegative=0.0;
1893  smallestPositive=COIN_DBL_MAX;
1894  largestPositive=0.0;
1895  // get matrix data pointers
1896  const double * elementByColumn = matrix_->getElements();
1897  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1898  const int * columnLength = matrix_->getVectorLengths(); 
1899  int numberColumns = matrix_->getNumCols();
1900  int i;
1901  for (i=0;i<numberColumns;i++) {
1902    CoinBigIndex j;
1903    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1904      double value = elementByColumn[j];
1905      if (value>0.0) {
1906        smallestPositive = CoinMin(smallestPositive,value);
1907        largestPositive = CoinMax(largestPositive,value);
1908      } else if (value<0.0) {
1909        smallestNegative = CoinMax(smallestNegative,value);
1910        largestNegative = CoinMin(largestNegative,value);
1911      }
1912    }
1913  }
1914}
1915// Says whether it can do partial pricing
1916bool 
1917ClpPackedMatrix::canDoPartialPricing() const
1918{
1919  return true;
1920}
1921// Partial pricing
1922void 
1923ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1924                              int & bestSequence, int & numberWanted)
1925{
1926  numberWanted=currentWanted_;
1927  int start = (int) (startFraction*numberActiveColumns_);
1928  int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
1929  const double * element =matrix_->getElements();
1930  const int * row = matrix_->getIndices();
1931  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1932  const int * length = matrix_->getVectorLengths();
1933  const double * rowScale = model->rowScale();
1934  const double * columnScale = model->columnScale();
1935  int iSequence;
1936  CoinBigIndex j;
1937  double tolerance=model->currentDualTolerance();
1938  double * reducedCost = model->djRegion();
1939  const double * duals = model->dualRowSolution();
1940  const double * cost = model->costRegion();
1941  double bestDj;
1942  if (bestSequence>=0)
1943    bestDj = fabs(model->clpMatrix()->reducedCost(model,bestSequence));
1944  else
1945    bestDj=tolerance;
1946  int sequenceOut = model->sequenceOut();
1947  int saveSequence = bestSequence;
1948  int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_; 
1949  int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
1950  if (rowScale) {
1951    // scaled
1952    for (iSequence=start;iSequence<end;iSequence++) {
1953      if (iSequence!=sequenceOut) {
1954        double value;
1955        ClpSimplex::Status status = model->getStatus(iSequence);
1956       
1957        switch(status) {
1958         
1959        case ClpSimplex::basic:
1960        case ClpSimplex::isFixed:
1961          break;
1962        case ClpSimplex::isFree:
1963        case ClpSimplex::superBasic:
1964          value=0.0;
1965          // scaled
1966          for (j=startColumn[iSequence];
1967               j<startColumn[iSequence]+length[iSequence];j++) {
1968            int jRow=row[j];
1969            value -= duals[jRow]*element[j]*rowScale[jRow];
1970          }
1971          value = fabs(cost[iSequence] +value*columnScale[iSequence]);
1972          if (value>FREE_ACCEPT*tolerance) {
1973            numberWanted--;
1974            // we are going to bias towards free (but only if reasonable)
1975            value *= FREE_BIAS;
1976            if (value>bestDj) {
1977              // check flagged variable and correct dj
1978              if (!model->flagged(iSequence)) {
1979                bestDj=value;
1980                bestSequence = iSequence;
1981              } else {
1982                // just to make sure we don't exit before got something
1983                numberWanted++;
1984              }
1985            }
1986          }
1987          break;
1988        case ClpSimplex::atUpperBound:
1989          value=0.0;
1990          // scaled
1991          for (j=startColumn[iSequence];
1992               j<startColumn[iSequence]+length[iSequence];j++) {
1993            int jRow=row[j];
1994            value -= duals[jRow]*element[j]*rowScale[jRow];
1995          }
1996          value = cost[iSequence] +value*columnScale[iSequence];
1997          if (value>tolerance) {
1998            numberWanted--;
1999            if (value>bestDj) {
2000              // check flagged variable and correct dj
2001              if (!model->flagged(iSequence)) {
2002                bestDj=value;
2003                bestSequence = iSequence;
2004              } else {
2005                // just to make sure we don't exit before got something
2006                numberWanted++;
2007              }
2008            }
2009          }
2010          break;
2011        case ClpSimplex::atLowerBound:
2012          value=0.0;
2013          // scaled
2014          for (j=startColumn[iSequence];
2015               j<startColumn[iSequence]+length[iSequence];j++) {
2016            int jRow=row[j];
2017            value -= duals[jRow]*element[j]*rowScale[jRow];
2018          }
2019          value = -(cost[iSequence] +value*columnScale[iSequence]);
2020          if (value>tolerance) {
2021            numberWanted--;
2022            if (value>bestDj) {
2023              // check flagged variable and correct dj
2024              if (!model->flagged(iSequence)) {
2025                bestDj=value;
2026                bestSequence = iSequence;
2027              } else {
2028                // just to make sure we don't exit before got something
2029                numberWanted++;
2030              }
2031            }
2032          }
2033          break;
2034        }
2035      }
2036      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2037        // give up
2038        break;
2039      }
2040      if (!numberWanted)
2041        break;
2042    }
2043    if (bestSequence!=saveSequence) {
2044      // recompute dj
2045      double value=0.0;
2046      // scaled
2047      for (j=startColumn[bestSequence];
2048           j<startColumn[bestSequence]+length[bestSequence];j++) {
2049        int jRow=row[j];
2050        value -= duals[jRow]*element[j]*rowScale[jRow];
2051      }
2052      reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
2053      savedBestSequence_ = bestSequence;
2054      savedBestDj_ = reducedCost[savedBestSequence_];
2055    }
2056  } else {
2057    // not scaled
2058    for (iSequence=start;iSequence<end;iSequence++) {
2059      if (iSequence!=sequenceOut) {
2060        double value;
2061        ClpSimplex::Status status = model->getStatus(iSequence);
2062       
2063        switch(status) {
2064         
2065        case ClpSimplex::basic:
2066        case ClpSimplex::isFixed:
2067          break;
2068        case ClpSimplex::isFree:
2069        case ClpSimplex::superBasic:
2070          value=cost[iSequence];
2071          for (j=startColumn[iSequence];
2072               j<startColumn[iSequence]+length[iSequence];j++) {
2073            int jRow=row[j];
2074            value -= duals[jRow]*element[j];
2075          }
2076          value = fabs(value);
2077          if (value>FREE_ACCEPT*tolerance) {
2078            numberWanted--;
2079            // we are going to bias towards free (but only if reasonable)
2080            value *= FREE_BIAS;
2081            if (value>bestDj) {
2082              // check flagged variable and correct dj
2083              if (!model->flagged(iSequence)) {
2084                bestDj=value;
2085                bestSequence = iSequence;
2086              } else {
2087                // just to make sure we don't exit before got something
2088                numberWanted++;
2089              }
2090            }
2091          }
2092          break;
2093        case ClpSimplex::atUpperBound:
2094          value=cost[iSequence];
2095          // scaled
2096          for (j=startColumn[iSequence];
2097               j<startColumn[iSequence]+length[iSequence];j++) {
2098            int jRow=row[j];
2099            value -= duals[jRow]*element[j];
2100          }
2101          if (value>tolerance) {
2102            numberWanted--;
2103            if (value>bestDj) {
2104              // check flagged variable and correct dj
2105              if (!model->flagged(iSequence)) {
2106                bestDj=value;
2107                bestSequence = iSequence;
2108              } else {
2109                // just to make sure we don't exit before got something
2110                numberWanted++;
2111              }
2112            }
2113          }
2114          break;
2115        case ClpSimplex::atLowerBound:
2116          value=cost[iSequence];
2117          for (j=startColumn[iSequence];
2118               j<startColumn[iSequence]+length[iSequence];j++) {
2119            int jRow=row[j];
2120            value -= duals[jRow]*element[j];
2121          }
2122          value = -value;
2123          if (value>tolerance) {
2124            numberWanted--;
2125            if (value>bestDj) {
2126              // check flagged variable and correct dj
2127              if (!model->flagged(iSequence)) {
2128                bestDj=value;
2129                bestSequence = iSequence;
2130              } else {
2131                // just to make sure we don't exit before got something
2132                numberWanted++;
2133              }
2134            }
2135          }
2136          break;
2137        }
2138      }
2139      if (numberWanted+minNeg<originalWanted_&&iSequence>lastScan) {
2140        // give up
2141        break;
2142      }
2143      if (!numberWanted)
2144        break;
2145    }
2146    if (bestSequence!=saveSequence) {
2147      // recompute dj
2148      double value=cost[bestSequence];
2149      for (j=startColumn[bestSequence];
2150           j<startColumn[bestSequence]+length[bestSequence];j++) {
2151        int jRow=row[j];
2152        value -= duals[jRow]*element[j];
2153      }
2154      reducedCost[bestSequence] = value;
2155      savedBestSequence_ = bestSequence;
2156      savedBestDj_ = reducedCost[savedBestSequence_];
2157    }
2158  }
2159  currentWanted_=numberWanted;
2160}
2161// Sets up an effective RHS
2162void 
2163ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
2164{
2165  delete [] rhsOffset_;
2166  int numberRows = model->numberRows();
2167  rhsOffset_ = new double[numberRows];
2168  rhsOffset(model,true);
2169}
2170// makes sure active columns correct
2171int 
2172ClpPackedMatrix::refresh(ClpSimplex * model)
2173{
2174  numberActiveColumns_ = matrix_->getNumCols();
2175  return 0;
2176}
2177
2178/* Scales rowCopy if column copy scaled
2179   Only called if scales already exist */
2180void 
2181ClpPackedMatrix::scaleRowCopy(ClpModel * model) const 
2182{
2183  if (model->rowCopy()) {
2184    // need to replace row by row
2185    int numberRows = model->numberRows();
2186    int numberColumns = matrix_->getNumCols();
2187    double * newElement = new double[numberColumns];
2188    ClpMatrixBase * rowCopyBase=model->rowCopy();
2189#ifndef NO_RTTI
2190    ClpPackedMatrix* rowCopy =
2191      dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
2192    // Make sure it is really a ClpPackedMatrix
2193    assert (rowCopy!=NULL);
2194#else
2195    ClpPackedMatrix* rowCopy =
2196      static_cast< ClpPackedMatrix*>(rowCopyBase);
2197#endif
2198
2199    const int * column = rowCopy->getIndices();
2200    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2201    const double * element = rowCopy->getElements();
2202    const double * rowScale = model->rowScale();
2203    const double * columnScale = model->columnScale();
2204    // scale row copy
2205    for (int iRow=0;iRow<numberRows;iRow++) {
2206      CoinBigIndex j;
2207      double scale = rowScale[iRow];
2208      const double * elementsInThisRow = element + rowStart[iRow];
2209      const int * columnsInThisRow = column + rowStart[iRow];
2210      int number = rowStart[iRow+1]-rowStart[iRow];
2211      assert (number<=numberColumns);
2212      for (j=0;j<number;j++) {
2213        int iColumn = columnsInThisRow[j];
2214        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
2215      }
2216      rowCopy->replaceVector(iRow,number,newElement);
2217    }
2218    delete [] newElement;
2219  }
2220}
2221/* Realy really scales column copy
2222   Only called if scales already exist.
2223   Up to user ro delete */
2224ClpMatrixBase * 
2225ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const 
2226{
2227  // need to replace column by column
2228  int numberRows = model->numberRows();
2229  int numberColumns = matrix_->getNumCols();
2230  double * newElement = new double[numberRows];
2231  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
2232  const int * row = copy->getIndices();
2233  const CoinBigIndex * columnStart = copy->getVectorStarts();
2234  const int * length = copy->getVectorLengths();
2235  const double * element = copy->getElements();
2236  const double * rowScale = model->rowScale();
2237  const double * columnScale = model->columnScale();
2238  // scale column copy
2239  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2240    CoinBigIndex j;
2241    double scale = columnScale[iColumn];
2242    const double * elementsInThisColumn = element + columnStart[iColumn];
2243    const int * rowsInThisColumn = row + columnStart[iColumn];
2244    int number = length[iColumn];
2245    assert (number<=numberRows);
2246    for (j=0;j<number;j++) {
2247      int iRow = rowsInThisColumn[j];
2248      newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
2249    }
2250    copy->replaceVector(iColumn,number,newElement);
2251  }
2252  delete [] newElement;
2253  return copy;
2254}
2255// Really scale matrix
2256void 
2257ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
2258{
2259  int numberColumns = matrix_->getNumCols();
2260  const int * row = matrix_->getIndices();
2261  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2262  const int * length = matrix_->getVectorLengths();
2263  double * element = matrix_->getMutableElements();
2264  // scale
2265  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2266    CoinBigIndex j;
2267    double scale = columnScale[iColumn];
2268    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
2269      int iRow = row[j];
2270      element[j] *= scale*rowScale[iRow];
2271    }
2272  }
2273}
Note: See TracBrowser for help on using the repository browser.