source: trunk/ClpNetworkMatrix.cpp @ 225

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

This should break everything

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