source: trunk/Clp/src/ClpNetworkMatrix.cpp @ 1304

Last change on this file since 1304 was 1302, checked in by forrest, 12 years ago

adding possibility of long doubles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.2 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#include "CoinPackedVector.hpp"
11
12#include "ClpSimplex.hpp"
13#include "ClpFactorization.hpp"
14// at end to get min/max!
15#include "ClpNetworkMatrix.hpp"
16#include "ClpPlusMinusOneMatrix.hpp"
17#include "ClpMessage.hpp"
18#include <iostream>
19#include <cassert>
20
21//#############################################################################
22// Constructors / Destructor / Assignment
23//#############################################################################
24
25//-------------------------------------------------------------------
26// Default Constructor
27//-------------------------------------------------------------------
28ClpNetworkMatrix::ClpNetworkMatrix () 
29  : ClpMatrixBase()
30{
31  setType(11);
32  matrix_ = NULL;
33  lengths_=NULL;
34  indices_=NULL;
35  numberRows_=0;
36  numberColumns_=0;
37  trueNetwork_=false;
38}
39
40/* Constructor from two arrays */
41ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
42                                   const int * tail)
43  : ClpMatrixBase()
44{
45  setType(11);
46  matrix_ = 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_ = CoinMax(numberRows_,iRow);
57    indices_[j]=iRow;
58    iRow = tail[iColumn];
59    numberRows_ = CoinMax(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  matrix_ = 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    CoinMemcpyN(rhs.indices_,2*numberColumns_,indices_);
79  }
80  int numberRows = getNumRows();
81  if (rhs.rhsOffset_&&numberRows) {
82    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
83  } else {
84    rhsOffset_=NULL;
85  }
86}
87
88ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs) 
89  : ClpMatrixBase()
90{ 
91  setType(11);
92  matrix_ = NULL;
93  lengths_=NULL;
94  indices_=NULL;
95  int iColumn;
96  assert (rhs.isColOrdered());
97  // get matrix data pointers
98  const int * row = rhs.getIndices();
99  const CoinBigIndex * columnStart = rhs.getVectorStarts();
100  const int * columnLength = rhs.getVectorLengths(); 
101  const double * elementByColumn = rhs.getElements();
102  numberColumns_ = rhs.getNumCols();
103  int goodNetwork=1;
104  numberRows_=-1;
105  indices_ = new int[2*numberColumns_];
106  CoinBigIndex j=0;
107  for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
108    CoinBigIndex k=columnStart[iColumn];
109    int iRow;
110    switch (columnLength[iColumn]) {
111    case 0:
112      goodNetwork=-1; // not classic network
113      indices_[j]=-1;
114      indices_[j+1]=-1;
115      break;
116     
117    case 1:
118      goodNetwork=-1; // not classic network
119      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
120        indices_[j] = -1;
121        iRow = row[k];
122        numberRows_ = CoinMax(numberRows_,iRow);
123        indices_[j+1]=iRow;
124      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
125        indices_[j+1] = -1;
126        iRow = row[k];
127        numberRows_ = CoinMax(numberRows_,iRow);
128        indices_[j]=iRow;
129      } else {
130        goodNetwork = 0; // not a network
131      }
132      break;
133     
134    case 2:
135      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
136        if (fabs(elementByColumn[k+1]+1.0)<1.0e-10) {
137          iRow = row[k];
138          numberRows_ = CoinMax(numberRows_,iRow);
139          indices_[j+1]=iRow;
140          iRow = row[k+1];
141          numberRows_ = CoinMax(numberRows_,iRow);
142          indices_[j] = iRow;
143        } else {
144          goodNetwork = 0; // not a network
145        }
146      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
147        if (fabs(elementByColumn[k+1]-1.0)<1.0e-10) {
148          iRow = row[k];
149          numberRows_ = CoinMax(numberRows_,iRow);
150          indices_[j]=iRow;
151          iRow = row[k+1];
152          numberRows_ = CoinMax(numberRows_,iRow);
153          indices_[j+1] = iRow;
154        } else {
155          goodNetwork = 0; // not a network
156        }
157      } else {
158        goodNetwork = 0; // not a network
159      }
160      break;
161
162    default:
163      goodNetwork = 0; // not a network
164      break;
165    }
166    if (!goodNetwork)
167      break;
168  }
169  if (!goodNetwork) {
170    delete [] indices_;
171    // put in message
172    printf("Not a network - can test if indices_ null\n");
173    indices_=NULL;
174    numberRows_=0;
175    numberColumns_=0;
176  } else {
177    numberRows_ ++; //  correct
178    trueNetwork_ = goodNetwork>0;
179  }
180}
181
182//-------------------------------------------------------------------
183// Destructor
184//-------------------------------------------------------------------
185ClpNetworkMatrix::~ClpNetworkMatrix ()
186{
187  delete matrix_;
188  delete [] lengths_;
189  delete [] indices_;
190}
191
192//----------------------------------------------------------------
193// Assignment operator
194//-------------------------------------------------------------------
195ClpNetworkMatrix &
196ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
197{
198  if (this != &rhs) {
199    ClpMatrixBase::operator=(rhs);
200    delete matrix_;
201    delete [] lengths_;
202    delete [] indices_;
203    matrix_ = NULL;
204    lengths_=NULL;
205    indices_=NULL;
206    numberRows_=rhs.numberRows_;
207    numberColumns_=rhs.numberColumns_;
208    trueNetwork_=rhs.trueNetwork_;
209    if (numberColumns_) {
210      indices_ = new int [ 2*numberColumns_];
211      CoinMemcpyN(rhs.indices_,2*numberColumns_,indices_);
212    }
213  }
214  return *this;
215}
216//-------------------------------------------------------------------
217// Clone
218//-------------------------------------------------------------------
219ClpMatrixBase * ClpNetworkMatrix::clone() const
220{
221  return new ClpNetworkMatrix(*this);
222}
223
224/* Returns a new matrix in reverse order without gaps */
225ClpMatrixBase * 
226ClpNetworkMatrix::reverseOrderedCopy() const
227{
228  // count number in each row
229  CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
230  CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
231  memset(tempP,0,numberRows_*sizeof(CoinBigIndex));
232  memset(tempN,0,numberRows_*sizeof(CoinBigIndex));
233  CoinBigIndex j=0;
234  int i;
235  for (i=0;i<numberColumns_;i++,j+=2) {
236    int iRow = indices_[j];
237    tempN[iRow]++;
238    iRow = indices_[j+1];
239    tempP[iRow]++;
240  }
241  int * newIndices = new int [2*numberColumns_];
242  CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
243  CoinBigIndex * newN = new CoinBigIndex[numberRows_];
244  int iRow;
245  j=0;
246  // do starts
247  for (iRow=0;iRow<numberRows_;iRow++) {
248    newP[iRow]=j;
249    j += tempP[iRow];
250    tempP[iRow]=newP[iRow];
251    newN[iRow] = j;
252    j += tempN[iRow];
253    tempN[iRow]=newN[iRow];
254  }
255  newP[numberRows_]=j;
256  j=0;
257  for (i=0;i<numberColumns_;i++,j+=2) {
258    int iRow = indices_[j];
259    CoinBigIndex put = tempN[iRow];
260    newIndices[put++] = i;
261    tempN[iRow] = put;
262    iRow = indices_[j+1];
263    put = tempP[iRow];
264    newIndices[put++] = i;
265    tempP[iRow] = put;
266  }
267  delete [] tempP;
268  delete [] tempN;
269  ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
270  newCopy->passInCopy(numberRows_, numberColumns_,
271                      false,  newIndices, newP, newN);
272  return newCopy;
273}
274//unscaled versions
275void 
276ClpNetworkMatrix::times(double scalar,
277                   const double * x, double * y) const
278{
279  int iColumn;
280  CoinBigIndex j=0;
281  if (trueNetwork_) {
282    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
283      double value = scalar*x[iColumn];
284      if (value) {
285        int iRowM = indices_[j];
286        int iRowP = indices_[j+1];
287        y[iRowM] -= value;
288        y[iRowP] += value;
289      }
290    }
291  } else {
292    // skip negative rows
293    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
294      double value = scalar*x[iColumn];
295      if (value) {
296        int iRowM = indices_[j];
297        int iRowP = indices_[j+1];
298        if (iRowM>=0)
299          y[iRowM] -= value;
300        if (iRowP>=0)
301          y[iRowP] += value;
302      }
303    }
304  }
305}
306void 
307ClpNetworkMatrix::transposeTimes(double scalar,
308                                const double * x, double * y) const
309{
310  int iColumn;
311  CoinBigIndex j=0;
312  if (trueNetwork_) {
313    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
314      double value = y[iColumn];
315      int iRowM = indices_[j];
316      int iRowP = indices_[j+1];
317      value -= scalar*x[iRowM];
318      value += scalar*x[iRowP];
319      y[iColumn] = value;
320    }
321  } else {
322    // skip negative rows
323    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
324      double value = y[iColumn];
325      int iRowM = indices_[j];
326      int iRowP = indices_[j+1];
327      if (iRowM>=0)
328        value -= scalar*x[iRowM];
329      if (iRowP>=0)
330        value += scalar*x[iRowP];
331      y[iColumn] = value;
332    }
333  }
334}
335void 
336ClpNetworkMatrix::times(double scalar,
337                       const double * x, double * y,
338                       const double * rowScale, 
339                       const double * columnScale) const
340{
341  // we know it is not scaled
342  times(scalar, x, y);
343}
344void 
345ClpNetworkMatrix::transposeTimes( double scalar,
346                                 const double * x, double * y,
347                                 const double * rowScale, 
348                                 const double * columnScale, double * spare) const
349{
350  // we know it is not scaled
351  transposeTimes(scalar, x, y);
352}
353/* Return <code>x * A + y</code> in <code>z</code>.
354        Squashes small elements and knows about ClpSimplex */
355void 
356ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
357                              const CoinIndexedVector * rowArray,
358                              CoinIndexedVector * y,
359                              CoinIndexedVector * columnArray) const
360{
361  // we know it is not scaled
362  columnArray->clear();
363  double * pi = rowArray->denseVector();
364  int numberNonZero=0;
365  int * index = columnArray->getIndices();
366  double * array = columnArray->denseVector();
367  int numberInRowArray = rowArray->getNumElements();
368  // maybe I need one in OsiSimplex
369  double zeroTolerance = model->factorization()->zeroTolerance();
370  int numberRows = model->numberRows();
371#ifndef NO_RTTI
372  ClpPlusMinusOneMatrix* rowCopy =
373    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
374#else
375  ClpPlusMinusOneMatrix* rowCopy =
376    static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
377#endif
378  bool packed = rowArray->packedMode();
379  double factor = 0.3;
380  // We may not want to do by row if there may be cache problems
381  int numberColumns = model->numberColumns();
382  // It would be nice to find L2 cache size - for moment 512K
383  // Be slightly optimistic
384  if (numberColumns*sizeof(double)>1000000) {
385    if (numberRows*10<numberColumns)
386      factor=0.1;
387    else if (numberRows*4<numberColumns)
388      factor=0.15;
389    else if (numberRows*2<numberColumns)
390      factor=0.2;
391    //if (model->numberIterations()%50==0)
392    //printf("%d nonzero\n",numberInRowArray);
393  }
394  if (numberInRowArray>factor*numberRows||!rowCopy) {
395    // do by column
396    int iColumn;
397    assert (!y->getNumElements());
398    CoinBigIndex j=0;
399    if (packed) {
400      // need to expand pi into y
401      assert(y->capacity()>=numberRows);
402      double * piOld = pi;
403      pi = y->denseVector();
404      const int * whichRow = rowArray->getIndices();
405      int i;
406      // modify pi so can collapse to one loop
407      for (i=0;i<numberInRowArray;i++) {
408        int iRow = whichRow[i];
409        pi[iRow]=scalar*piOld[i];
410      }
411      if (trueNetwork_) {
412        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
413          double value = 0.0;
414          int iRowM = indices_[j];
415          int iRowP = indices_[j+1];
416          value -= pi[iRowM];
417          value += pi[iRowP];
418          if (fabs(value)>zeroTolerance) {
419            array[numberNonZero]=value;
420            index[numberNonZero++]=iColumn;
421          }
422        }
423      } else {
424        // skip negative rows
425        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
426          double value = 0.0;
427          int iRowM = indices_[j];
428          int iRowP = indices_[j+1];
429          if (iRowM>=0)
430            value -= pi[iRowM];
431          if (iRowP>=0)
432            value += pi[iRowP];
433          if (fabs(value)>zeroTolerance) {
434            array[numberNonZero]=value;
435            index[numberNonZero++]=iColumn;
436          }
437        }
438      }
439      for (i=0;i<numberInRowArray;i++) {
440        int iRow = whichRow[i];
441        pi[iRow]=0.0;
442      }
443    } else {
444      if (trueNetwork_) {
445        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
446          double value = 0.0;
447          int iRowM = indices_[j];
448          int iRowP = indices_[j+1];
449          value -= scalar*pi[iRowM];
450          value += scalar*pi[iRowP];
451          if (fabs(value)>zeroTolerance) {
452            index[numberNonZero++]=iColumn;
453            array[iColumn]=value;
454          }
455        }
456      } else {
457        // skip negative rows
458        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
459          double value = 0.0;
460          int iRowM = indices_[j];
461          int iRowP = indices_[j+1];
462          if (iRowM>=0)
463            value -= scalar*pi[iRowM];
464          if (iRowP>=0)
465            value += scalar*pi[iRowP];
466          if (fabs(value)>zeroTolerance) {
467            index[numberNonZero++]=iColumn;
468            array[iColumn]=value;
469          }
470        }
471      }
472    }
473    columnArray->setNumElements(numberNonZero);
474  } else {
475    // do by row
476    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
477  }
478}
479/* Return <code>x *A in <code>z</code> but
480   just for indices in y. */
481void 
482ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * model,
483                              const CoinIndexedVector * rowArray,
484                              const CoinIndexedVector * y,
485                              CoinIndexedVector * columnArray) const
486{
487  columnArray->clear();
488  double * pi = rowArray->denseVector();
489  double * array = columnArray->denseVector();
490  int jColumn;
491  int numberToDo = y->getNumElements();
492  const int * which = y->getIndices();
493  assert (!rowArray->packedMode());
494  columnArray->setPacked();
495  if (trueNetwork_) {
496    for (jColumn=0;jColumn<numberToDo;jColumn++) {
497      int iColumn = which[jColumn];
498      double value = 0.0;
499      CoinBigIndex j=iColumn<<1;
500      int iRowM = indices_[j];
501      int iRowP = indices_[j+1];
502      value -= pi[iRowM];
503      value += pi[iRowP];
504      array[jColumn]=value;
505    }
506  } else {
507    // skip negative rows
508    for (jColumn=0;jColumn<numberToDo;jColumn++) {
509      int iColumn = which[jColumn];
510      double value = 0.0;
511      CoinBigIndex j=iColumn<<1;
512      int iRowM = indices_[j];
513      int iRowP = indices_[j+1];
514      if (iRowM>=0)
515        value -= pi[iRowM];
516      if (iRowP>=0)
517        value += pi[iRowP];
518      array[jColumn]=value;
519    }
520  }
521}
522/// returns number of elements in column part of basis,
523CoinBigIndex
524ClpNetworkMatrix::countBasis(ClpSimplex * model,
525                                 const int * whichColumn, 
526                                 int numberBasic,
527                             int & numberColumnBasic)
528{
529  int i;
530  CoinBigIndex numberElements=0;
531  if (trueNetwork_) {
532    numberElements = 2*numberColumnBasic;
533  } else {
534    for (i=0;i<numberColumnBasic;i++) {
535      int iColumn = whichColumn[i];
536      CoinBigIndex j=iColumn<<1;
537      int iRowM = indices_[j];
538      int iRowP = indices_[j+1];
539      if (iRowM>=0)
540        numberElements ++;
541      if (iRowP>=0)
542        numberElements ++;
543    }
544  }
545  return numberElements;
546}
547void
548ClpNetworkMatrix::fillBasis(ClpSimplex * model,
549                         const int * whichColumn, 
550                         int & numberColumnBasic,
551                         int * indexRowU, int * start,
552                         int * rowCount, int * columnCount,
553                         CoinFactorizationDouble * elementU)
554{
555  int i;
556  CoinBigIndex numberElements=start[0];
557  if (trueNetwork_) {
558    for (i=0;i<numberColumnBasic;i++) {
559      int iColumn = whichColumn[i];
560      CoinBigIndex j=iColumn<<1;
561      int iRowM = indices_[j];
562      int iRowP = indices_[j+1];
563      indexRowU[numberElements]=iRowM;
564      rowCount[iRowM]++;
565      elementU[numberElements]=-1.0;
566      indexRowU[numberElements+1]=iRowP;
567      rowCount[iRowP]++;
568      elementU[numberElements+1]=1.0;
569      numberElements+=2;
570      start[i+1]=numberElements;
571      columnCount[i]=2;
572    }
573  } else {
574    for (i=0;i<numberColumnBasic;i++) {
575      int iColumn = whichColumn[i];
576      CoinBigIndex j=iColumn<<1;
577      int iRowM = indices_[j];
578      int iRowP = indices_[j+1];
579      if (iRowM>=0) {
580        indexRowU[numberElements]=iRowM;
581        rowCount[iRowM]++;
582        elementU[numberElements++]=-1.0;
583      }
584      if (iRowP>=0) {
585        indexRowU[numberElements]=iRowP;
586        rowCount[iRowP]++;
587        elementU[numberElements++]=1.0;
588      }
589      start[i+1]=numberElements;
590      columnCount[i]=numberElements-start[i];
591    }
592  }
593}
594/* Unpacks a column into an CoinIndexedvector
595 */
596void 
597ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
598                   int iColumn) const 
599{
600  CoinBigIndex j=iColumn<<1;
601  int iRowM = indices_[j];
602  int iRowP = indices_[j+1];
603  if (iRowM>=0) 
604    rowArray->add(iRowM,-1.0);
605  if (iRowP>=0) 
606    rowArray->add(iRowP,1.0);
607}
608/* Unpacks a column into an CoinIndexedvector
609** in packed foramt
610Note that model is NOT const.  Bounds and objective could
611be modified if doing column generation (just for this variable) */
612void 
613ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
614                            CoinIndexedVector * rowArray,
615                            int iColumn) const
616{
617  int * index = rowArray->getIndices();
618  double * array = rowArray->denseVector();
619  int number = 0;
620  CoinBigIndex j=iColumn<<1;
621  int iRowM = indices_[j];
622  int iRowP = indices_[j+1];
623  if (iRowM>=0) {
624    array[number]=-1.0;
625    index[number++]=iRowM;
626  }
627  if (iRowP>=0) {
628    array[number]=1.0;
629    index[number++]=iRowP;
630  }
631  rowArray->setNumElements(number);
632  rowArray->setPackedMode(true);
633}
634/* Adds multiple of a column into an CoinIndexedvector
635      You can use quickAdd to add to vector */
636void 
637ClpNetworkMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
638                   int iColumn, double multiplier) const 
639{
640  CoinBigIndex j=iColumn<<1;
641  int iRowM = indices_[j];
642  int iRowP = indices_[j+1];
643  if (iRowM>=0) 
644    rowArray->quickAdd(iRowM,-multiplier);
645  if (iRowP>=0) 
646    rowArray->quickAdd(iRowP,multiplier);
647}
648/* Adds multiple of a column into an array */
649void 
650ClpNetworkMatrix::add(const ClpSimplex * model,double * array,
651                    int iColumn, double multiplier) const
652{
653  CoinBigIndex j=iColumn<<1;
654  int iRowM = indices_[j];
655  int iRowP = indices_[j+1];
656  if (iRowM>=0) 
657    array[iRowM] -= multiplier;
658  if (iRowP>=0) 
659    array[iRowP] += multiplier;
660}
661
662// Return a complete CoinPackedMatrix
663CoinPackedMatrix * 
664ClpNetworkMatrix::getPackedMatrix() const 
665{
666  if (!matrix_) {
667    assert (trueNetwork_); // fix later
668    int numberElements = 2*numberColumns_;
669    double * elements = new double [numberElements];
670    CoinBigIndex i;
671    for (i=0;i<2*numberColumns_;i+=2) {
672      elements[i]=-1.0;
673      elements[i+1]=1.0;
674    }
675    CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
676    for (i=0;i<numberColumns_+1;i++) {
677      starts[i]=2*i;
678    }
679    // use assignMatrix to save space
680    delete [] lengths_;
681    lengths_ = NULL;
682    matrix_ = new CoinPackedMatrix();
683    int * indices = CoinCopyOfArray(indices_,2*numberColumns_);
684    matrix_->assignMatrix(true,numberRows_,numberColumns_,
685                          getNumElements(),
686                          elements,indices,
687                          starts,lengths_);
688    assert(!elements);
689    assert(!starts);
690    assert (!indices);
691    assert (!lengths_);
692  }
693  return matrix_;
694}
695/* A vector containing the elements in the packed matrix. Note that there
696   might be gaps in this list, entries that do not belong to any
697   major-dimension vector. To get the actual elements one should look at
698   this vector together with vectorStarts and vectorLengths. */
699const double * 
700ClpNetworkMatrix::getElements() const 
701{
702  if (!matrix_) 
703    getPackedMatrix();
704  return matrix_->getElements();
705}
706
707const CoinBigIndex * 
708ClpNetworkMatrix::getVectorStarts() const 
709{
710  if (!matrix_) 
711    getPackedMatrix();
712  return matrix_->getVectorStarts();
713}
714/* The lengths of the major-dimension vectors. */
715const int * 
716ClpNetworkMatrix::getVectorLengths() const
717{
718  assert (trueNetwork_); // fix later
719  if (!lengths_) {
720    lengths_ = new int [numberColumns_];
721    int i;
722    for (i=0;i<numberColumns_;i++) {
723      lengths_[i]=2;
724    }
725  }
726  return lengths_;
727}
728/* Delete the columns whose indices are listed in <code>indDel</code>. */
729void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
730{
731  assert (trueNetwork_);
732  int iColumn;
733  int numberBad=0;
734  // Use array to make sure we can have duplicates
735  char * which = new char[numberColumns_];
736  memset(which,0,numberColumns_);
737  int nDuplicate=0;
738  for (iColumn=0;iColumn<numDel;iColumn++) {
739    int jColumn = indDel[iColumn];
740    if (jColumn<0||jColumn>=numberColumns_) {
741      numberBad++;
742    } else {
743      if (which[jColumn])
744        nDuplicate++;
745      else
746        which[jColumn]=1;
747    }
748  }
749  if (numberBad)
750    throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
751  int newNumber = numberColumns_-numDel+nDuplicate;
752  // Get rid of temporary arrays
753  delete [] lengths_;
754  lengths_=NULL;
755  delete matrix_;
756  matrix_= NULL;
757  int newSize=2*newNumber;
758  int * newIndices = new int [newSize];
759  newSize=0;
760  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
761    if (!which[iColumn]) {
762      CoinBigIndex start;
763      CoinBigIndex i;
764      start = 2*iColumn;
765      for (i=start;i<start+2;i++) 
766        newIndices[newSize++]=indices_[i];
767    }
768  }
769  delete [] which;
770  delete [] indices_;
771  indices_= newIndices;
772  numberColumns_ = newNumber;
773}
774/* Delete the rows whose indices are listed in <code>indDel</code>. */
775void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
776{
777  int iRow;
778  int numberBad=0;
779  // Use array to make sure we can have duplicates
780  int * which = new int [numberRows_];
781  memset(which,0,numberRows_*sizeof(int));
782  for (iRow=0;iRow<numDel;iRow++) {
783    int jRow = indDel[iRow];
784    if (jRow<0||jRow>=numberRows_) {
785      numberBad++;
786    } else {
787      which[jRow]=1;
788    }
789  }
790  if (numberBad)
791    throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
792  // Only valid of all columns have 0 entries
793  int iColumn;
794  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
795    CoinBigIndex start;
796    CoinBigIndex i;
797    start = 2*iColumn;
798    for (i=start;i<start+2;i++) {
799      int iRow = indices_[i];
800      if (which[iRow])
801        numberBad++;
802    }
803  }
804  if (numberBad)
805    throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
806  int newNumber=0;
807  for (iRow=0;iRow<numberRows_;iRow++) {
808    if (!which[iRow])
809      which[iRow]=newNumber++;
810    else
811      which[iRow]=-1;
812  }
813  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
814    CoinBigIndex start;
815    CoinBigIndex i;
816    start = 2*iColumn;
817    for (i=start;i<start+2;i++) {
818      int iRow = indices_[i];
819      indices_[i]=which[iRow];
820    }
821  }
822  delete [] which;
823  numberRows_ = newNumber;
824}
825/* Given positive integer weights for each row fills in sum of weights
826   for each column (and slack).
827   Returns weights vector
828*/
829CoinBigIndex * 
830ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
831{
832  int numberRows = model->numberRows();
833  int numberColumns =model->numberColumns();
834  int number = numberRows+numberColumns;
835  CoinBigIndex * weights = new CoinBigIndex[number];
836  int i;
837  for (i=0;i<numberColumns;i++) {
838    CoinBigIndex j=i<<1;
839    CoinBigIndex count=0;
840    int iRowM = indices_[j];
841    int iRowP = indices_[j+1];
842    if (iRowM>=0) {
843      count += inputWeights[iRowM];
844    }
845    if (iRowP>=0) {
846      count += inputWeights[iRowP];
847    }
848    weights[i]=count;
849  }
850  for (i=0;i<numberRows;i++) {
851    weights[i+numberColumns]=inputWeights[i];
852  }
853  return weights;
854}
855/* Returns largest and smallest elements of both signs.
856   Largest refers to largest absolute value.
857*/
858void 
859ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
860                       double & smallestPositive, double & largestPositive)
861{
862  smallestNegative=-1.0;
863  largestNegative=-1.0;
864  smallestPositive=1.0;
865  largestPositive=1.0;
866}
867// Says whether it can do partial pricing
868bool 
869ClpNetworkMatrix::canDoPartialPricing() const
870{
871  return true; 
872}
873// Partial pricing
874void 
875ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
876                              int & bestSequence, int & numberWanted)
877{
878  numberWanted=currentWanted_;
879  int j;
880  int start = (int) (startFraction*numberColumns_);
881  int end = CoinMin((int) (endFraction*numberColumns_+1),numberColumns_);
882  double tolerance=model->currentDualTolerance();
883  double * reducedCost = model->djRegion();
884  const double * duals = model->dualRowSolution();
885  const double * cost = model->costRegion();
886  double bestDj;
887  if (bestSequence>=0)
888    bestDj = fabs(reducedCost[bestSequence]);
889  else
890    bestDj=tolerance;
891  int sequenceOut = model->sequenceOut();
892  int saveSequence = bestSequence;
893  if (!trueNetwork_) {
894    // Not true network
895    int iSequence;
896    for (iSequence=start;iSequence<end;iSequence++) {
897      if (iSequence!=sequenceOut) {
898        double value;
899        int iRowM,iRowP;
900        ClpSimplex::Status status = model->getStatus(iSequence);
901
902        switch(status) {
903         
904        case ClpSimplex::basic:
905        case ClpSimplex::isFixed:
906          break;
907        case ClpSimplex::isFree:
908        case ClpSimplex::superBasic:
909          value=cost[iSequence];
910          j = iSequence<<1;
911          // skip negative rows
912          iRowM = indices_[j];
913          iRowP = indices_[j+1];
914          if (iRowM>=0)
915            value += duals[iRowM];
916          if (iRowP>=0)
917            value -= duals[iRowP];
918          value = fabs(value);
919          if (value>FREE_ACCEPT*tolerance) {
920            numberWanted--;
921            // we are going to bias towards free (but only if reasonable)
922            value *= FREE_BIAS;
923            if (value>bestDj) {
924              // check flagged variable and correct dj
925              if (!model->flagged(iSequence)) {
926                bestDj=value;
927                bestSequence = iSequence;
928              } else {
929                // just to make sure we don't exit before got something
930                numberWanted++;
931              }
932            }
933          }
934          break;
935        case ClpSimplex::atUpperBound:
936          value=cost[iSequence];
937          j = iSequence<<1;
938          // skip negative rows
939          iRowM = indices_[j];
940          iRowP = indices_[j+1];
941          if (iRowM>=0)
942            value += duals[iRowM];
943          if (iRowP>=0)
944            value -= duals[iRowP];
945          if (value>tolerance) {
946            numberWanted--;
947            if (value>bestDj) {
948              // check flagged variable and correct dj
949              if (!model->flagged(iSequence)) {
950                bestDj=value;
951                bestSequence = iSequence;
952              } else {
953                // just to make sure we don't exit before got something
954                numberWanted++;
955              }
956            }
957          }
958          break;
959        case ClpSimplex::atLowerBound:
960          value=cost[iSequence];
961          j = iSequence<<1;
962          // skip negative rows
963          iRowM = indices_[j];
964          iRowP = indices_[j+1];
965          if (iRowM>=0)
966            value += duals[iRowM];
967          if (iRowP>=0)
968            value -= duals[iRowP];
969          value = -value;
970          if (value>tolerance) {
971            numberWanted--;
972            if (value>bestDj) {
973              // check flagged variable and correct dj
974              if (!model->flagged(iSequence)) {
975                bestDj=value;
976                bestSequence = iSequence;
977              } else {
978                // just to make sure we don't exit before got something
979                numberWanted++;
980              }
981            }
982          }
983          break;
984        }
985      }
986      if (!numberWanted)
987        break;
988    }
989    if (bestSequence!=saveSequence) {
990      // recompute dj
991      double value=cost[bestSequence];
992      j = bestSequence<<1;
993      // skip negative rows
994      int iRowM = indices_[j];
995      int iRowP = indices_[j+1];
996      if (iRowM>=0)
997        value += duals[iRowM];
998      if (iRowP>=0)
999        value -= duals[iRowP];
1000      reducedCost[bestSequence] = value;
1001      savedBestSequence_ = bestSequence;
1002      savedBestDj_ = reducedCost[savedBestSequence_];
1003    }
1004  } else {
1005    // true network
1006    int iSequence;
1007    for (iSequence=start;iSequence<end;iSequence++) {
1008      if (iSequence!=sequenceOut) {
1009        double value;
1010        int iRowM,iRowP;
1011        ClpSimplex::Status status = model->getStatus(iSequence);
1012
1013        switch(status) {
1014         
1015        case ClpSimplex::basic:
1016        case ClpSimplex::isFixed:
1017          break;
1018        case ClpSimplex::isFree:
1019        case ClpSimplex::superBasic:
1020          value=cost[iSequence];
1021          j = iSequence<<1;
1022          iRowM = indices_[j];
1023          iRowP = indices_[j+1];
1024          value += duals[iRowM];
1025          value -= duals[iRowP];
1026          value = fabs(value);
1027          if (value>FREE_ACCEPT*tolerance) {
1028            numberWanted--;
1029            // we are going to bias towards free (but only if reasonable)
1030            value *= FREE_BIAS;
1031            if (value>bestDj) {
1032              // check flagged variable and correct dj
1033              if (!model->flagged(iSequence)) {
1034                bestDj=value;
1035                bestSequence = iSequence;
1036              } else {
1037                // just to make sure we don't exit before got something
1038                numberWanted++;
1039              }
1040            }
1041          }
1042          break;
1043        case ClpSimplex::atUpperBound:
1044          value=cost[iSequence];
1045          j = iSequence<<1;
1046          iRowM = indices_[j];
1047          iRowP = indices_[j+1];
1048          value += duals[iRowM];
1049          value -= duals[iRowP];
1050          if (value>tolerance) {
1051            numberWanted--;
1052            if (value>bestDj) {
1053              // check flagged variable and correct dj
1054              if (!model->flagged(iSequence)) {
1055                bestDj=value;
1056                bestSequence = iSequence;
1057              } else {
1058                // just to make sure we don't exit before got something
1059                numberWanted++;
1060              }
1061            }
1062          }
1063          break;
1064        case ClpSimplex::atLowerBound:
1065          value=cost[iSequence];
1066          j = iSequence<<1;
1067          iRowM = indices_[j];
1068          iRowP = indices_[j+1];
1069          value += duals[iRowM];
1070          value -= duals[iRowP];
1071          value = -value;
1072          if (value>tolerance) {
1073            numberWanted--;
1074            if (value>bestDj) {
1075              // check flagged variable and correct dj
1076              if (!model->flagged(iSequence)) {
1077                bestDj=value;
1078                bestSequence = iSequence;
1079              } else {
1080                // just to make sure we don't exit before got something
1081                numberWanted++;
1082              }
1083            }
1084          }
1085          break;
1086        }
1087      }
1088      if (!numberWanted)
1089        break;
1090    }
1091    if (bestSequence!=saveSequence) {
1092      // recompute dj
1093      double value=cost[bestSequence];
1094      j = bestSequence<<1;
1095      int iRowM = indices_[j];
1096      int iRowP = indices_[j+1];
1097      value += duals[iRowM];
1098      value -= duals[iRowP];
1099      reducedCost[bestSequence] = value;
1100      savedBestSequence_ = bestSequence;
1101      savedBestDj_ = reducedCost[savedBestSequence_];
1102    }
1103  }
1104  currentWanted_=numberWanted;
1105}
1106// Allow any parts of a created CoinMatrix to be deleted
1107void 
1108ClpNetworkMatrix::releasePackedMatrix() const 
1109{
1110  delete matrix_;
1111  delete [] lengths_;
1112  matrix_=NULL;
1113  lengths_=NULL;
1114}
1115// Append Columns
1116void 
1117ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1118{
1119  int iColumn;
1120  int numberBad=0;
1121  for (iColumn=0;iColumn<number;iColumn++) {
1122    int n=columns[iColumn]->getNumElements();
1123    const double * element = columns[iColumn]->getElements();
1124    if (n!=2)
1125      numberBad++;
1126    if (fabs(element[0])!=1.0||fabs(element[1])!=1.0)
1127        numberBad++;
1128    else if (element[0]*element[1]!=-1.0)
1129        numberBad++;
1130  }
1131  if (numberBad)
1132    throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
1133  // Get rid of temporary arrays
1134  delete [] lengths_;
1135  lengths_=NULL;
1136  delete matrix_;
1137  matrix_= NULL;
1138  CoinBigIndex size = 2*number;
1139  int * temp2 = new int [numberColumns_*2+size];
1140  CoinMemcpyN(indices_,numberColumns_*2,temp2);
1141  delete [] indices_;
1142  indices_= temp2;
1143  // now add
1144  size=2*numberColumns_;
1145  for (iColumn=0;iColumn<number;iColumn++) {
1146    const int * row = columns[iColumn]->getIndices();
1147    const double * element = columns[iColumn]->getElements();
1148    if (element[0]==-1.0) {
1149      indices_[size++]=row[0];
1150      indices_[size++]=row[1];
1151    } else {
1152      indices_[size++]=row[1];
1153      indices_[size++]=row[0];
1154    }
1155  }
1156 
1157  numberColumns_ += number;
1158}
1159// Append Rows
1160void 
1161ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1162{
1163  // must be zero arrays
1164  int numberBad=0;
1165  int iRow;
1166  for (iRow=0;iRow<number;iRow++) {
1167    numberBad += rows[iRow]->getNumElements();
1168  }
1169  if (numberBad)
1170    throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
1171  numberRows_ += number;
1172}
1173#ifndef SLIM_CLP
1174/* Append a set of rows/columns to the end of the matrix. Returns number of errors
1175   i.e. if any of the new rows/columns contain an index that's larger than the
1176   number of columns-1/rows-1 (if numberOther>0) or duplicates
1177   If 0 then rows, 1 if columns */
1178int 
1179ClpNetworkMatrix::appendMatrix(int number, int type,
1180                                    const CoinBigIndex * starts, const int * index,
1181                                    const double * element, int numberOther)
1182{
1183  int numberErrors=0;
1184  // make into CoinPackedVector
1185  CoinPackedVectorBase ** vectors=
1186    new CoinPackedVectorBase * [number];
1187  int iVector;
1188  for (iVector=0;iVector<number;iVector++) {
1189    int iStart = starts[iVector];
1190    vectors[iVector] = 
1191      new CoinPackedVector(starts[iVector+1]-iStart,
1192                           index+iStart,element+iStart);
1193  }
1194  if (type==0) {
1195    // rows
1196    appendRows(number,vectors);
1197  } else {
1198    // columns
1199    appendCols(number,vectors);
1200  }
1201  for (iVector=0;iVector<number;iVector++) 
1202    delete vectors[iVector];
1203  delete [] vectors;
1204  return numberErrors;
1205}
1206#endif
1207/* Subset clone (without gaps).  Duplicates are allowed
1208   and order is as given */
1209ClpMatrixBase * 
1210ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
1211                              int numberColumns, 
1212                              const int * whichColumns) const 
1213{
1214  return new ClpNetworkMatrix(*this, numberRows, whichRows,
1215                                   numberColumns, whichColumns);
1216}
1217/* Subset constructor (without gaps).  Duplicates are allowed
1218   and order is as given */
1219ClpNetworkMatrix::ClpNetworkMatrix (
1220                       const ClpNetworkMatrix & rhs,
1221                       int numberRows, const int * whichRow,
1222                       int numberColumns, const int * whichColumn)
1223: ClpMatrixBase(rhs)
1224{ 
1225  setType(11);
1226  matrix_ = NULL;
1227  lengths_=NULL;
1228  indices_=new int[2*numberColumns];;
1229  numberRows_=numberRows;
1230  numberColumns_=numberColumns;
1231  trueNetwork_=true;
1232  int iColumn;
1233  int numberBad=0;
1234  int * which = new int [rhs.numberRows_];
1235  int iRow;
1236  for (iRow=0;iRow<rhs.numberRows_;iRow++) 
1237    which[iRow]=-1;
1238  int n=0;
1239  for (iRow=0;iRow<numberRows;iRow++) {
1240    int jRow=whichRow[iRow];
1241    assert (jRow>=0&&jRow<rhs.numberRows_);
1242    which[jRow]=n++;
1243  }
1244  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1245    CoinBigIndex start;
1246    CoinBigIndex i;
1247    start = 2*iColumn;
1248    CoinBigIndex offset = 2*whichColumn[iColumn]-start;
1249    for (i=start;i<start+2;i++) {
1250      int iRow = rhs.indices_[i+offset];
1251      iRow = which[iRow];
1252      if (iRow<0)
1253        numberBad++;
1254      else
1255        indices_[i]=iRow;
1256    }
1257  }
1258  if (numberBad)
1259    throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
1260}
Note: See TracBrowser for help on using the repository browser.