source: tags/preroot/ClpNetworkMatrix.cpp @ 778

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

Hopefully this will fix tag problem

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.7 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  if (numberInRowArray>0.3*numberRows||!rowCopy) {
372    // do by column
373    int iColumn;
374    double * markVector = y->denseVector(); // probably empty
375    CoinBigIndex j=0;
376    if (trueNetwork_) {
377      for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
378        double value = markVector[iColumn];
379        markVector[iColumn]=0.0;
380        int iRowM = indices_[j];
381        int iRowP = indices_[j+1];
382        value -= scalar*pi[iRowM];
383        value += scalar*pi[iRowP];
384        if (fabs(value)>zeroTolerance) {
385          index[numberNonZero++]=iColumn;
386          array[iColumn]=value;
387        }
388      }
389    } else {
390      // skip negative rows
391      for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
392        double value = markVector[iColumn];
393        markVector[iColumn]=0.0;
394        int iRowM = indices_[j];
395        int iRowP = indices_[j+1];
396        if (iRowM>=0)
397          value -= scalar*pi[iRowM];
398        if (iRowP>=0)
399          value += scalar*pi[iRowP];
400        if (fabs(value)>zeroTolerance) {
401          index[numberNonZero++]=iColumn;
402          array[iColumn]=value;
403        }
404      }
405    }
406    columnArray->setNumElements(numberNonZero);
407    y->setNumElements(0);
408  } else {
409    // do by row
410    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
411  }
412}
413/* Return <code>x *A in <code>z</code> but
414   just for indices in y.
415   Squashes small elements and knows about ClpSimplex */
416void 
417ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * model,
418                              const CoinIndexedVector * rowArray,
419                              const CoinIndexedVector * y,
420                              CoinIndexedVector * columnArray) const
421{
422  columnArray->clear();
423  double * pi = rowArray->denseVector();
424  int numberNonZero=0;
425  int * index = columnArray->getIndices();
426  double * array = columnArray->denseVector();
427  // maybe I need one in OsiSimplex
428  double zeroTolerance = model->factorization()->zeroTolerance();
429  int jColumn;
430  int numberToDo = y->getNumElements();
431  const int * which = y->getIndices();
432  if (trueNetwork_) {
433    for (jColumn=0;jColumn<numberToDo;jColumn++) {
434      int iColumn = which[jColumn];
435      double value = 0.0;
436      CoinBigIndex j=iColumn<<1;
437      int iRowM = indices_[j];
438      int iRowP = indices_[j+1];
439      value -= pi[iRowM];
440      value += pi[iRowP];
441      if (fabs(value)>zeroTolerance) {
442        index[numberNonZero++]=iColumn;
443        array[iColumn]=value;
444      }
445    }
446  } else {
447    // skip negative rows
448    for (jColumn=0;jColumn<numberToDo;jColumn++) {
449      int iColumn = which[jColumn];
450      double value = 0.0;
451      CoinBigIndex j=iColumn<<1;
452      int iRowM = indices_[j];
453      int iRowP = indices_[j+1];
454      if (iRowM>=0)
455        value -= pi[iRowM];
456      if (iRowP>=0)
457        value += pi[iRowP];
458      if (fabs(value)>zeroTolerance) {
459        index[numberNonZero++]=iColumn;
460        array[iColumn]=value;
461      }
462    }
463  }
464}
465/* Returns number of elements in basis
466   column is basic if entry >=0 */
467CoinBigIndex
468ClpNetworkMatrix::numberInBasis(const int * columnIsBasic) const 
469{
470  int i;
471  CoinBigIndex numberElements=0;
472  if (trueNetwork_) {
473    for (i=0;i<numberColumns_;i++) {
474      if (columnIsBasic[i]>=0) 
475        numberElements ++;
476    }
477    numberElements *= 2;
478  } else {
479    for (i=0;i<numberColumns_;i++) {
480      if (columnIsBasic[i]>=0) {
481        CoinBigIndex j=i<<1;
482        int iRowM = indices_[j];
483        int iRowP = indices_[j+1];
484        if (iRowM>=0)
485          numberElements ++;
486        if (iRowP>=0)
487          numberElements ++;
488      }
489    }
490  }
491  return numberElements;
492}
493// Fills in basis (Returns number of elements and updates numberBasic)
494CoinBigIndex
495ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
496                                const int * columnIsBasic, int & numberBasic,
497                                int * indexRowU, int * indexColumnU,
498                                double * elementU) const 
499{
500#ifdef CLPDEBUG
501  const double * rowScale = model->rowScale();
502  assert (!rowScale);
503#endif
504  int i;
505  CoinBigIndex numberElements=0;
506  if (trueNetwork_) {
507    for (i=0;i<numberColumns_;i++) {
508      if (columnIsBasic[i]>=0) {
509        CoinBigIndex j=i<<1;
510        int iRowM = indices_[j];
511        int iRowP = indices_[j+1];
512        indexRowU[numberElements]=iRowM;
513        indexColumnU[numberElements]=numberBasic;
514        elementU[numberElements]=-1.0;
515        indexRowU[numberElements+1]=iRowP;
516        indexColumnU[numberElements+1]=numberBasic;
517        elementU[numberElements+1]=1.0;
518        numberElements+=2;
519        numberBasic++;
520      }
521    }
522  } else {
523    for (i=0;i<numberColumns_;i++) {
524      if (columnIsBasic[i]>=0) {
525        CoinBigIndex j=i<<1;
526        int iRowM = indices_[j];
527        int iRowP = indices_[j+1];
528        if (iRowM>=0) {
529          indexRowU[numberElements]=iRowM;
530          indexColumnU[numberElements]=numberBasic;
531          elementU[numberElements++]=-1.0;
532        }
533        if (iRowP>=0) {
534          indexRowU[numberElements]=iRowP;
535          indexColumnU[numberElements]=numberBasic;
536          elementU[numberElements++]=1.0;
537        }
538        numberBasic++;
539      }
540    }
541  }
542  return numberElements;
543}
544/* Unpacks a column into an CoinIndexedvector
545      Note that model is NOT const.  Bounds and objective could
546      be modified if doing column generation */
547void 
548ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
549                   int iColumn) const 
550{
551  CoinBigIndex j=iColumn<<1;
552  int iRowM = indices_[j];
553  int iRowP = indices_[j+1];
554  if (iRowM>=0) 
555    rowArray->add(iRowM,-1.0);
556  if (iRowP>=0) 
557    rowArray->add(iRowP,1.0);
558}
559/* Adds multiple of a column into an CoinIndexedvector
560      You can use quickAdd to add to vector */
561void 
562ClpNetworkMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
563                   int iColumn, double multiplier) const 
564{
565  CoinBigIndex j=iColumn<<1;
566  int iRowM = indices_[j];
567  int iRowP = indices_[j+1];
568  if (iRowM>=0) 
569    rowArray->quickAdd(iRowM,-multiplier);
570  if (iRowP>=0) 
571    rowArray->quickAdd(iRowP,multiplier);
572}
573
574// Return a complete CoinPackedMatrix
575CoinPackedMatrix * 
576ClpNetworkMatrix::getPackedMatrix() const 
577{
578  return new CoinPackedMatrix(true,numberRows_,numberColumns_,
579                              2*numberColumns_,
580                              getElements(),indices_,
581                              getVectorStarts(),getVectorLengths());
582
583}
584/* A vector containing the elements in the packed matrix. Note that there
585   might be gaps in this list, entries that do not belong to any
586   major-dimension vector. To get the actual elements one should look at
587   this vector together with vectorStarts and vectorLengths. */
588const double * 
589ClpNetworkMatrix::getElements() const 
590{
591  assert (trueNetwork_); // fix later
592  if (!elements_) {
593    elements_ = new double [2*numberColumns_];
594    int i;
595    for (i=0;i<2*numberColumns_;i+=2) {
596      elements_[i]=-1.0;
597      elements_[i+1]=1.0;
598    }
599  }
600  return elements_;
601}
602
603const CoinBigIndex * 
604ClpNetworkMatrix::getVectorStarts() const 
605{
606  assert (trueNetwork_); // fix later
607  if (!starts_) {
608    starts_ = new int [numberColumns_+1];
609    int i;
610    for (i=0;i<numberColumns_+1;i++) {
611      starts_[i]=i;
612    }
613  }
614  return starts_;
615}
616/* The lengths of the major-dimension vectors. */
617const int * 
618ClpNetworkMatrix::getVectorLengths() const
619{
620  assert (trueNetwork_); // fix later
621  if (!lengths_) {
622    lengths_ = new int [numberColumns_];
623    int i;
624    for (i=0;i<numberColumns_;i++) {
625      lengths_[i]=2;
626    }
627  }
628  return lengths_;
629}
630/* Delete the columns whose indices are listed in <code>indDel</code>. */
631void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
632{
633  abort();
634}
635/* Delete the rows whose indices are listed in <code>indDel</code>. */
636void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
637{
638  abort();
639}
Note: See TracBrowser for help on using the repository browser.