source: branches/pre/ClpNetworkMatrix.cpp @ 213

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

Trying to make faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.4 KB
Line 
1// Copyright (C) 2003, 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 "ClpNetworkMatrix.hpp"
15#include "ClpPlusMinusOneMatrix.hpp"
16#include "ClpMessage.hpp"
17
18//#############################################################################
19// Constructors / Destructor / Assignment
20//#############################################################################
21
22//-------------------------------------------------------------------
23// Default Constructor
24//-------------------------------------------------------------------
25ClpNetworkMatrix::ClpNetworkMatrix () 
26  : ClpMatrixBase()
27{
28  setType(11);
29  elements_ = NULL;
30  starts_ = NULL;
31  lengths_=NULL;
32  indices_=NULL;
33  numberRows_=0;
34  numberColumns_=0;
35  trueNetwork_=false;
36}
37
38/* Constructor from two arrays */
39ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
40                                   const int * tail)
41  : ClpMatrixBase()
42{
43  setType(11);
44  elements_ = NULL;
45  starts_ = NULL;
46  lengths_=NULL;
47  indices_=new int[2*numberColumns];;
48  numberRows_=-1;
49  numberColumns_=numberColumns;
50  trueNetwork_=true;
51  int iColumn;
52  CoinBigIndex j=0;
53  for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
54    int iRow = head[iColumn];
55    numberRows_ = max(numberRows_,iRow);
56    indices_[j]=iRow;
57    iRow = tail[iColumn];
58    numberRows_ = max(numberRows_,iRow);
59    indices_[j+1]=iRow;
60  }
61  numberRows_++;
62}
63//-------------------------------------------------------------------
64// Copy constructor
65//-------------------------------------------------------------------
66ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs) 
67: ClpMatrixBase(rhs)
68{ 
69  elements_ = NULL;
70  starts_ = NULL;
71  lengths_=NULL;
72  indices_=NULL;
73  numberRows_=rhs.numberRows_;
74  numberColumns_=rhs.numberColumns_;
75  trueNetwork_=rhs.trueNetwork_;
76  if (numberColumns_) {
77    indices_ = new int [ 2*numberColumns_];
78    memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
79  }
80}
81
82ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs) 
83  : ClpMatrixBase()
84{ 
85  setType(11);
86  elements_ = NULL;
87  starts_ = NULL;
88  lengths_=NULL;
89  indices_=NULL;
90  int iColumn;
91  assert (rhs.isColOrdered());
92  // get matrix data pointers
93  const int * row = rhs.getIndices();
94  const CoinBigIndex * columnStart = rhs.getVectorStarts();
95  const int * columnLength = rhs.getVectorLengths(); 
96  const double * elementByColumn = rhs.getElements();
97  numberColumns_ = rhs.getNumCols();
98  int goodNetwork=1;
99  numberRows_=-1;
100  indices_ = new int[2*numberColumns_];
101  CoinBigIndex j=0;
102  for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
103    CoinBigIndex k=columnStart[iColumn];
104    int iRow;
105    switch (columnLength[iColumn]) {
106    case 0:
107      goodNetwork=-1; // not classic network
108      indices_[j]=-1;
109      indices_[j+1]=-1;
110      break;
111     
112    case 1:
113      goodNetwork=-1; // not classic network
114      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
115        indices_[j] = -1;
116        iRow = row[k];
117        numberRows_ = max(numberRows_,iRow);
118        indices_[j+1]=iRow;
119      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
120        indices_[j+1] = -1;
121        iRow = row[k];
122        numberRows_ = max(numberRows_,iRow);
123        indices_[j]=iRow;
124      } else {
125        goodNetwork = 0; // not a network
126      }
127      break;
128     
129    case 2:
130      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
131        if (fabs(elementByColumn[k+1]+1.0)<1.0e-10) {
132          iRow = row[k];
133          numberRows_ = max(numberRows_,iRow);
134          indices_[j+1]=iRow;
135          iRow = row[k+1];
136          numberRows_ = max(numberRows_,iRow);
137          indices_[j] = iRow;
138        } else {
139          goodNetwork = 0; // not a network
140        }
141      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
142        if (fabs(elementByColumn[k+1]-1.0)<1.0e-10) {
143          iRow = row[k];
144          numberRows_ = max(numberRows_,iRow);
145          indices_[j]=iRow;
146          iRow = row[k+1];
147          numberRows_ = max(numberRows_,iRow);
148          indices_[j+1] = iRow;
149        } else {
150          goodNetwork = 0; // not a network
151        }
152      } else {
153        goodNetwork = 0; // not a network
154      }
155      break;
156
157    default:
158      goodNetwork = 0; // not a network
159      break;
160    }
161    if (!goodNetwork)
162      break;
163  }
164  if (!goodNetwork) {
165    delete [] indices_;
166    // put in message
167    printf("Not a network - can test if indices_ null\n");
168    indices_=NULL;
169    numberRows_=0;
170    numberColumns_=0;
171  } else {
172    numberRows_ ++; //  correct
173    trueNetwork_ = goodNetwork>0;
174  }
175}
176
177//-------------------------------------------------------------------
178// Destructor
179//-------------------------------------------------------------------
180ClpNetworkMatrix::~ClpNetworkMatrix ()
181{
182  delete [] elements_;
183  delete [] starts_;
184  delete [] lengths_;
185  delete [] indices_;
186}
187
188//----------------------------------------------------------------
189// Assignment operator
190//-------------------------------------------------------------------
191ClpNetworkMatrix &
192ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
193{
194  if (this != &rhs) {
195    ClpMatrixBase::operator=(rhs);
196    delete [] elements_;
197    delete [] starts_;
198    delete [] lengths_;
199    delete [] indices_;
200    elements_ = NULL;
201    starts_ = NULL;
202    lengths_=NULL;
203    indices_=NULL;
204    numberRows_=rhs.numberRows_;
205    numberColumns_=rhs.numberColumns_;
206    trueNetwork_=rhs.trueNetwork_;
207    if (numberColumns_) {
208      indices_ = new int [ 2*numberColumns_];
209      memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
210    }
211  }
212  return *this;
213}
214//-------------------------------------------------------------------
215// Clone
216//-------------------------------------------------------------------
217ClpMatrixBase * ClpNetworkMatrix::clone() const
218{
219  return new ClpNetworkMatrix(*this);
220}
221
222/* Returns a new matrix in reverse order without gaps */
223ClpMatrixBase * 
224ClpNetworkMatrix::reverseOrderedCopy() const
225{
226  // count number in each row
227  int * tempP = new int [numberRows_];
228  int * tempN = new int [numberRows_];
229  memset(tempP,0,numberRows_*sizeof(int));
230  memset(tempN,0,numberRows_*sizeof(int));
231  CoinBigIndex j=0;
232  int i;
233  for (i=0;i<numberColumns_;i++,j+=2) {
234    int iRow = indices_[j];
235    tempN[iRow]++;
236    iRow = indices_[j+1];
237    tempP[iRow]++;
238  }
239  int * newIndices = new int [2*numberColumns_];
240  int * newP = new int [numberRows_+1];
241  int * newN = new int[numberRows_];
242  int iRow;
243  j=0;
244  // do starts
245  for (iRow=0;iRow<numberRows_;iRow++) {
246    newP[iRow]=j;
247    j += tempP[iRow];
248    tempP[iRow]=newP[iRow];
249    newN[iRow] = j;
250    j += tempN[iRow];
251    tempN[iRow]=newN[iRow];
252  }
253  newP[numberRows_]=j;
254  j=0;
255  for (i=0;i<numberColumns_;i++,j+=2) {
256    int iRow = indices_[j];
257    int put = tempN[iRow];
258    newIndices[put++] = i;
259    tempN[iRow] = put;
260    iRow = indices_[j+1];
261    put = tempP[iRow];
262    newIndices[put++] = i;
263    tempP[iRow] = put;
264  }
265  delete [] tempP;
266  delete [] tempN;
267  ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
268  newCopy->passInCopy(numberRows_, numberColumns_,
269                      false,  newIndices, newP, newN);
270  return newCopy;
271}
272//unscaled versions
273void 
274ClpNetworkMatrix::times(double scalar,
275                   const double * x, double * y) const
276{
277  int iColumn;
278  CoinBigIndex j=0;
279  if (trueNetwork_) {
280    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
281      double value = scalar*x[iColumn];
282      if (value) {
283        int iRowM = indices_[j];
284        int iRowP = indices_[j+1];
285        y[iRowM] -= value;
286        y[iRowP] += value;
287      }
288    }
289  } else {
290    // skip negative rows
291    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
292      double value = scalar*x[iColumn];
293      if (value) {
294        int iRowM = indices_[j];
295        int iRowP = indices_[j+1];
296        if (iRowM>=0)
297          y[iRowM] -= value;
298        if (iRowP>=0)
299          y[iRowP] += value;
300      }
301    }
302  }
303}
304void 
305ClpNetworkMatrix::transposeTimes(double scalar,
306                                const double * x, double * y) const
307{
308  int iColumn;
309  CoinBigIndex j=0;
310  if (trueNetwork_) {
311    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
312      double value = y[iColumn];
313      int iRowM = indices_[j];
314      int iRowP = indices_[j+1];
315      value -= scalar*x[iRowM];
316      value += scalar*x[iRowP];
317      y[iColumn] = value;
318    }
319  } else {
320    // skip negative rows
321    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
322      double value = y[iColumn];
323      int iRowM = indices_[j];
324      int iRowP = indices_[j+1];
325      if (iRowM>=0)
326        value -= scalar*x[iRowM];
327      if (iRowP>=0)
328        value += scalar*x[iRowP];
329      y[iColumn] = value;
330    }
331  }
332}
333void 
334ClpNetworkMatrix::times(double scalar,
335                       const double * x, double * y,
336                       const double * rowScale, 
337                       const double * columnScale) const
338{
339  // we know it is not scaled
340  times(scalar, x, y);
341}
342void 
343ClpNetworkMatrix::transposeTimes( double scalar,
344                                 const double * x, double * y,
345                                 const double * rowScale, 
346                                 const double * columnScale) const
347{
348  // we know it is not scaled
349  transposeTimes(scalar, x, y);
350}
351/* Return <code>x * A + y</code> in <code>z</code>.
352        Squashes small elements and knows about ClpSimplex */
353void 
354ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
355                              const CoinIndexedVector * rowArray,
356                              CoinIndexedVector * y,
357                              CoinIndexedVector * columnArray) const
358{
359  // we know it is not scaled
360  columnArray->clear();
361  double * pi = rowArray->denseVector();
362  int numberNonZero=0;
363  int * index = columnArray->getIndices();
364  double * array = columnArray->denseVector();
365  int numberInRowArray = rowArray->getNumElements();
366  // maybe I need one in OsiSimplex
367  double zeroTolerance = model->factorization()->zeroTolerance();
368  int numberRows = model->numberRows();
369  ClpPlusMinusOneMatrix* rowCopy =
370    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
371  bool packed = rowArray->packedMode();
372  double factor = 0.3;
373  // We may not want to do by row if there may be cache problems
374  int numberColumns = model->numberColumns();
375  // It would be nice to find L2 cache size - for moment 512K
376  // Be slightly optimistic
377  if (numberColumns*sizeof(double)>1000000) {
378    if (numberRows*10<numberColumns)
379      factor=0.1;
380    else if (numberRows*4<numberColumns)
381      factor=0.15;
382    else if (numberRows*2<numberColumns)
383      factor=0.2;
384    if (model->numberIterations()%50==0)
385      printf("%d nonzero\n",numberInRowArray);
386  }
387  if (numberInRowArray>factor*numberRows||!rowCopy) {
388    // do by column
389    int iColumn;
390    assert (!y->getNumElements());
391    CoinBigIndex j=0;
392    if (packed) {
393      // need to expand pi into y
394      assert(y->capacity()>=numberRows);
395      double * piOld = pi;
396      pi = y->denseVector();
397      const int * whichRow = rowArray->getIndices();
398      int i;
399      // modify pi so can collapse to one loop
400      for (i=0;i<numberInRowArray;i++) {
401        int iRow = whichRow[i];
402        pi[iRow]=scalar*piOld[i];
403      }
404      if (trueNetwork_) {
405        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
406          double value = 0.0;
407          int iRowM = indices_[j];
408          int iRowP = indices_[j+1];
409          value -= pi[iRowM];
410          value += pi[iRowP];
411          if (fabs(value)>zeroTolerance) {
412            array[numberNonZero]=value;
413            index[numberNonZero++]=iColumn;
414          }
415        }
416      } else {
417        // skip negative rows
418        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
419          double value = 0.0;
420          int iRowM = indices_[j];
421          int iRowP = indices_[j+1];
422          if (iRowM>=0)
423            value -= pi[iRowM];
424          if (iRowP>=0)
425            value += pi[iRowP];
426          if (fabs(value)>zeroTolerance) {
427            array[numberNonZero]=value;
428            index[numberNonZero++]=iColumn;
429          }
430        }
431      }
432      for (i=0;i<numberInRowArray;i++) {
433        int iRow = whichRow[i];
434        pi[iRow]=0.0;
435      }
436    } else {
437      if (trueNetwork_) {
438        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
439          double value = 0.0;
440          int iRowM = indices_[j];
441          int iRowP = indices_[j+1];
442          value -= scalar*pi[iRowM];
443          value += scalar*pi[iRowP];
444          if (fabs(value)>zeroTolerance) {
445            index[numberNonZero++]=iColumn;
446            array[iColumn]=value;
447          }
448        }
449      } else {
450        // skip negative rows
451        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
452          double value = 0.0;
453          int iRowM = indices_[j];
454          int iRowP = indices_[j+1];
455          if (iRowM>=0)
456            value -= scalar*pi[iRowM];
457          if (iRowP>=0)
458            value += scalar*pi[iRowP];
459          if (fabs(value)>zeroTolerance) {
460            index[numberNonZero++]=iColumn;
461            array[iColumn]=value;
462          }
463        }
464      }
465    }
466    columnArray->setNumElements(numberNonZero);
467  } else {
468    // do by row
469    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
470  }
471}
472/* Return <code>x *A in <code>z</code> but
473   just for indices in y.
474   Squashes small elements and knows about ClpSimplex */
475void 
476ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * model,
477                              const CoinIndexedVector * rowArray,
478                              const CoinIndexedVector * y,
479                              CoinIndexedVector * columnArray) const
480{
481  columnArray->clear();
482  double * pi = rowArray->denseVector();
483  int numberNonZero=0;
484  int * index = columnArray->getIndices();
485  double * array = columnArray->denseVector();
486  // maybe I need one in OsiSimplex
487  double zeroTolerance = model->factorization()->zeroTolerance();
488  int jColumn;
489  int numberToDo = y->getNumElements();
490  const int * which = y->getIndices();
491  bool packed = rowArray->packedMode();
492  if (packed) {
493    // Must line up with y
494    // need to expand pi into y
495    int numberInRowArray = rowArray->getNumElements();
496    assert(y->capacity()>=model->numberRows());
497    double * piOld = pi;
498    pi = y->denseVector();
499    const int * whichRow = rowArray->getIndices();
500    int i;
501    for (i=0;i<numberInRowArray;i++) {
502      int iRow = whichRow[i];
503      pi[iRow]=piOld[i];
504    }
505    if (trueNetwork_) {
506      for (jColumn=0;jColumn<numberToDo;jColumn++) {
507        int iColumn = which[jColumn];
508        double value = 0.0;
509        CoinBigIndex j=iColumn<<1;
510        int iRowM = indices_[j];
511        int iRowP = indices_[j+1];
512        value -= pi[iRowM];
513        value += pi[iRowP];
514        array[jColumn]=value;
515      }
516    } else {
517      // skip negative rows
518      for (jColumn=0;jColumn<numberToDo;jColumn++) {
519        int iColumn = which[jColumn];
520        double value = 0.0;
521        CoinBigIndex j=iColumn<<1;
522        int iRowM = indices_[j];
523        int iRowP = indices_[j+1];
524        if (iRowM>=0)
525          value -= pi[iRowM];
526        if (iRowP>=0)
527          value += pi[iRowP];
528        array[jColumn]=value;
529      }
530    }
531    for (i=0;i<numberInRowArray;i++) {
532      int iRow = whichRow[i];
533      pi[iRow]=0.0;
534    }
535  } else {
536    if (trueNetwork_) {
537      for (jColumn=0;jColumn<numberToDo;jColumn++) {
538        int iColumn = which[jColumn];
539        double value = 0.0;
540        CoinBigIndex j=iColumn<<1;
541        int iRowM = indices_[j];
542        int iRowP = indices_[j+1];
543        value -= pi[iRowM];
544        value += pi[iRowP];
545        if (fabs(value)>zeroTolerance) {
546          index[numberNonZero++]=iColumn;
547          array[iColumn]=value;
548        }
549      }
550    } else {
551      // skip negative rows
552      for (jColumn=0;jColumn<numberToDo;jColumn++) {
553        int iColumn = which[jColumn];
554        double value = 0.0;
555        CoinBigIndex j=iColumn<<1;
556        int iRowM = indices_[j];
557        int iRowP = indices_[j+1];
558        if (iRowM>=0)
559          value -= pi[iRowM];
560        if (iRowP>=0)
561          value += pi[iRowP];
562        if (fabs(value)>zeroTolerance) {
563          index[numberNonZero++]=iColumn;
564          array[iColumn]=value;
565        }
566      }
567    }
568  }
569}
570/* Returns number of elements in basis
571   column is basic if entry >=0 */
572CoinBigIndex
573ClpNetworkMatrix::numberInBasis(const int * columnIsBasic) const 
574{
575  int i;
576  CoinBigIndex numberElements=0;
577  if (trueNetwork_) {
578    for (i=0;i<numberColumns_;i++) {
579      if (columnIsBasic[i]>=0) 
580        numberElements ++;
581    }
582    numberElements *= 2;
583  } else {
584    for (i=0;i<numberColumns_;i++) {
585      if (columnIsBasic[i]>=0) {
586        CoinBigIndex j=i<<1;
587        int iRowM = indices_[j];
588        int iRowP = indices_[j+1];
589        if (iRowM>=0)
590          numberElements ++;
591        if (iRowP>=0)
592          numberElements ++;
593      }
594    }
595  }
596  return numberElements;
597}
598// Fills in basis (Returns number of elements and updates numberBasic)
599CoinBigIndex
600ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
601                                const int * columnIsBasic, int & numberBasic,
602                                int * indexRowU, int * indexColumnU,
603                                double * elementU) const 
604{
605#ifdef CLPDEBUG
606  const double * rowScale = model->rowScale();
607  assert (!rowScale);
608#endif
609  int i;
610  CoinBigIndex numberElements=0;
611  if (trueNetwork_) {
612    for (i=0;i<numberColumns_;i++) {
613      if (columnIsBasic[i]>=0) {
614        CoinBigIndex j=i<<1;
615        int iRowM = indices_[j];
616        int iRowP = indices_[j+1];
617        indexRowU[numberElements]=iRowM;
618        indexColumnU[numberElements]=numberBasic;
619        elementU[numberElements]=-1.0;
620        indexRowU[numberElements+1]=iRowP;
621        indexColumnU[numberElements+1]=numberBasic;
622        elementU[numberElements+1]=1.0;
623        numberElements+=2;
624        numberBasic++;
625      }
626    }
627  } else {
628    for (i=0;i<numberColumns_;i++) {
629      if (columnIsBasic[i]>=0) {
630        CoinBigIndex j=i<<1;
631        int iRowM = indices_[j];
632        int iRowP = indices_[j+1];
633        if (iRowM>=0) {
634          indexRowU[numberElements]=iRowM;
635          indexColumnU[numberElements]=numberBasic;
636          elementU[numberElements++]=-1.0;
637        }
638        if (iRowP>=0) {
639          indexRowU[numberElements]=iRowP;
640          indexColumnU[numberElements]=numberBasic;
641          elementU[numberElements++]=1.0;
642        }
643        numberBasic++;
644      }
645    }
646  }
647  return numberElements;
648}
649/* If element NULL returns number of elements in column part of basis,
650   If not NULL fills in as well */
651CoinBigIndex
652ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
653                                 const int * whichColumn, 
654                                 int numberBasic,
655                                 int numberColumnBasic,
656                                 int * indexRowU, int * indexColumnU,
657                                 double * elementU) const 
658{
659  int i;
660  CoinBigIndex numberElements=0;
661  if (elementU!=NULL) {
662    if (trueNetwork_) {
663      for (i=0;i<numberColumnBasic;i++) {
664        int iColumn = whichColumn[i];
665        CoinBigIndex j=iColumn<<1;
666        int iRowM = indices_[j];
667        int iRowP = indices_[j+1];
668        indexRowU[numberElements]=iRowM;
669        indexColumnU[numberElements]=numberBasic;
670        elementU[numberElements]=-1.0;
671        indexRowU[numberElements+1]=iRowP;
672        indexColumnU[numberElements+1]=numberBasic;
673        elementU[numberElements+1]=1.0;
674        numberElements+=2;
675        numberBasic++;
676      }
677    } else {
678      for (i=0;i<numberColumnBasic;i++) {
679        int iColumn = whichColumn[i];
680        CoinBigIndex j=iColumn<<1;
681        int iRowM = indices_[j];
682        int iRowP = indices_[j+1];
683        if (iRowM>=0) {
684          indexRowU[numberElements]=iRowM;
685          indexColumnU[numberElements]=numberBasic;
686          elementU[numberElements++]=-1.0;
687        }
688        if (iRowP>=0) {
689          indexRowU[numberElements]=iRowP;
690          indexColumnU[numberElements]=numberBasic;
691          elementU[numberElements++]=1.0;
692        }
693        numberBasic++;
694      }
695    }
696  } else {
697    if (trueNetwork_) {
698      numberElements = 2*numberColumnBasic;
699    } else {
700      for (i=0;i<numberColumnBasic;i++) {
701        int iColumn = whichColumn[i];
702        CoinBigIndex j=iColumn<<1;
703        int iRowM = indices_[j];
704        int iRowP = indices_[j+1];
705        if (iRowM>=0)
706          numberElements ++;
707        if (iRowP>=0)
708          numberElements ++;
709      }
710    }
711  }
712  return numberElements;
713}
714/* Unpacks a column into an CoinIndexedvector
715 */
716void 
717ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
718                   int iColumn) const 
719{
720  CoinBigIndex j=iColumn<<1;
721  int iRowM = indices_[j];
722  int iRowP = indices_[j+1];
723  if (iRowM>=0) 
724    rowArray->add(iRowM,-1.0);
725  if (iRowP>=0) 
726    rowArray->add(iRowP,1.0);
727}
728/* Unpacks a column into an CoinIndexedvector
729** in packed foramt
730Note that model is NOT const.  Bounds and objective could
731be modified if doing column generation (just for this variable) */
732void 
733ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
734                            CoinIndexedVector * rowArray,
735                            int iColumn) const
736{
737  int * index = rowArray->getIndices();
738  double * array = rowArray->denseVector();
739  int number = 0;
740  CoinBigIndex j=iColumn<<1;
741  int iRowM = indices_[j];
742  int iRowP = indices_[j+1];
743  if (iRowM>=0) {
744    array[number]=-1.0;
745    index[number++]=iRowM;
746  }
747  if (iRowP>=0) {
748    array[number]=1.0;
749    index[number++]=iRowP;
750  }
751  rowArray->setNumElements(number);
752  rowArray->setPackedMode(true);
753}
754/* Adds multiple of a column into an CoinIndexedvector
755      You can use quickAdd to add to vector */
756void 
757ClpNetworkMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
758                   int iColumn, double multiplier) const 
759{
760  CoinBigIndex j=iColumn<<1;
761  int iRowM = indices_[j];
762  int iRowP = indices_[j+1];
763  if (iRowM>=0) 
764    rowArray->quickAdd(iRowM,-multiplier);
765  if (iRowP>=0) 
766    rowArray->quickAdd(iRowP,multiplier);
767}
768
769// Return a complete CoinPackedMatrix
770CoinPackedMatrix * 
771ClpNetworkMatrix::getPackedMatrix() const 
772{
773  return new CoinPackedMatrix(true,numberRows_,numberColumns_,
774                              2*numberColumns_,
775                              getElements(),indices_,
776                              getVectorStarts(),getVectorLengths());
777
778}
779/* A vector containing the elements in the packed matrix. Note that there
780   might be gaps in this list, entries that do not belong to any
781   major-dimension vector. To get the actual elements one should look at
782   this vector together with vectorStarts and vectorLengths. */
783const double * 
784ClpNetworkMatrix::getElements() const 
785{
786  assert (trueNetwork_); // fix later
787  if (!elements_) {
788    elements_ = new double [2*numberColumns_];
789    int i;
790    for (i=0;i<2*numberColumns_;i+=2) {
791      elements_[i]=-1.0;
792      elements_[i+1]=1.0;
793    }
794  }
795  return elements_;
796}
797
798const CoinBigIndex * 
799ClpNetworkMatrix::getVectorStarts() const 
800{
801  assert (trueNetwork_); // fix later
802  if (!starts_) {
803    starts_ = new int [numberColumns_+1];
804    int i;
805    for (i=0;i<numberColumns_+1;i++) {
806      starts_[i]=i;
807    }
808  }
809  return starts_;
810}
811/* The lengths of the major-dimension vectors. */
812const int * 
813ClpNetworkMatrix::getVectorLengths() const
814{
815  assert (trueNetwork_); // fix later
816  if (!lengths_) {
817    lengths_ = new int [numberColumns_];
818    int i;
819    for (i=0;i<numberColumns_;i++) {
820      lengths_[i]=2;
821    }
822  }
823  return lengths_;
824}
825/* Delete the columns whose indices are listed in <code>indDel</code>. */
826void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
827{
828  abort();
829}
830/* Delete the rows whose indices are listed in <code>indDel</code>. */
831void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
832{
833  abort();
834}
835/* Given positive integer weights for each row fills in sum of weights
836   for each column (and slack).
837   Returns weights vector
838*/
839CoinBigIndex * 
840ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
841{
842  int numberRows = model->numberRows();
843  int numberColumns =model->numberColumns();
844  int number = numberRows+numberColumns;
845  CoinBigIndex * weights = new CoinBigIndex[number];
846  int i;
847  for (i=0;i<numberColumns;i++) {
848    CoinBigIndex j=i<<1;
849    CoinBigIndex count=0;
850    int iRowM = indices_[j];
851    int iRowP = indices_[j+1];
852    if (iRowM>=0) {
853      count += inputWeights[iRowM];
854    }
855    if (iRowP>=0) {
856      count += inputWeights[iRowP];
857    }
858    weights[i]=count;
859  }
860  for (i=0;i<numberRows;i++) {
861    weights[i+numberColumns]=inputWeights[i];
862  }
863  return weights;
864}
Note: See TracBrowser for help on using the repository browser.