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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 44.9 KB
Line 
1/* $Id: ClpNetworkMatrix.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5
6#include <cstdio>
7
8#include "CoinPragma.hpp"
9#include "CoinIndexedVector.hpp"
10#include "CoinHelperFunctions.hpp"
11#include "CoinPackedVector.hpp"
12
13#include "ClpSimplex.hpp"
14#include "ClpFactorization.hpp"
15// at end to get min/max!
16#include "ClpNetworkMatrix.hpp"
17#include "ClpPlusMinusOneMatrix.hpp"
18#include "ClpMessage.hpp"
19#include <iostream>
20#include <cassert>
21
22//#############################################################################
23// Constructors / Destructor / Assignment
24//#############################################################################
25
26//-------------------------------------------------------------------
27// Default Constructor
28//-------------------------------------------------------------------
29ClpNetworkMatrix::ClpNetworkMatrix ()
30     : ClpMatrixBase()
31{
32     setType(11);
33     matrix_ = NULL;
34     lengths_ = NULL;
35     indices_ = NULL;
36     numberRows_ = 0;
37     numberColumns_ = 0;
38     trueNetwork_ = false;
39}
40
41/* Constructor from two arrays */
42ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
43                                   const int * tail)
44     : ClpMatrixBase()
45{
46     setType(11);
47     matrix_ = NULL;
48     lengths_ = NULL;
49     indices_ = new int[2*numberColumns];;
50     numberRows_ = -1;
51     numberColumns_ = numberColumns;
52     trueNetwork_ = true;
53     int iColumn;
54     CoinBigIndex j = 0;
55     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
56          int iRow = head[iColumn];
57          numberRows_ = CoinMax(numberRows_, iRow);
58          indices_[j] = iRow;
59          iRow = tail[iColumn];
60          numberRows_ = CoinMax(numberRows_, iRow);
61          indices_[j+1] = iRow;
62     }
63     numberRows_++;
64}
65//-------------------------------------------------------------------
66// Copy constructor
67//-------------------------------------------------------------------
68ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs)
69     : ClpMatrixBase(rhs)
70{
71     matrix_ = 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          CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
80     }
81     int numberRows = getNumRows();
82     if (rhs.rhsOffset_ && numberRows) {
83          rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
84     } else {
85          rhsOffset_ = NULL;
86     }
87}
88
89ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs)
90     : ClpMatrixBase()
91{
92     setType(11);
93     matrix_ = NULL;
94     lengths_ = NULL;
95     indices_ = NULL;
96     int iColumn;
97     assert (rhs.isColOrdered());
98     // get matrix data pointers
99     const int * row = rhs.getIndices();
100     const CoinBigIndex * columnStart = rhs.getVectorStarts();
101     const int * columnLength = rhs.getVectorLengths();
102     const double * elementByColumn = rhs.getElements();
103     numberColumns_ = rhs.getNumCols();
104     int goodNetwork = 1;
105     numberRows_ = -1;
106     indices_ = new int[2*numberColumns_];
107     CoinBigIndex j = 0;
108     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
109          CoinBigIndex k = columnStart[iColumn];
110          int iRow;
111          switch (columnLength[iColumn]) {
112          case 0:
113               goodNetwork = -1; // not classic network
114               indices_[j] = -1;
115               indices_[j+1] = -1;
116               break;
117
118          case 1:
119               goodNetwork = -1; // not classic network
120               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
121                    indices_[j] = -1;
122                    iRow = row[k];
123                    numberRows_ = CoinMax(numberRows_, iRow);
124                    indices_[j+1] = iRow;
125               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
126                    indices_[j+1] = -1;
127                    iRow = row[k];
128                    numberRows_ = CoinMax(numberRows_, iRow);
129                    indices_[j] = iRow;
130               } else {
131                    goodNetwork = 0; // not a network
132               }
133               break;
134
135          case 2:
136               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
137                    if (fabs(elementByColumn[k+1] + 1.0) < 1.0e-10) {
138                         iRow = row[k];
139                         numberRows_ = CoinMax(numberRows_, iRow);
140                         indices_[j+1] = iRow;
141                         iRow = row[k+1];
142                         numberRows_ = CoinMax(numberRows_, iRow);
143                         indices_[j] = iRow;
144                    } else {
145                         goodNetwork = 0; // not a network
146                    }
147               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
148                    if (fabs(elementByColumn[k+1] - 1.0) < 1.0e-10) {
149                         iRow = row[k];
150                         numberRows_ = CoinMax(numberRows_, iRow);
151                         indices_[j] = iRow;
152                         iRow = row[k+1];
153                         numberRows_ = CoinMax(numberRows_, iRow);
154                         indices_[j+1] = iRow;
155                    } else {
156                         goodNetwork = 0; // not a network
157                    }
158               } else {
159                    goodNetwork = 0; // not a network
160               }
161               break;
162
163          default:
164               goodNetwork = 0; // not a network
165               break;
166          }
167          if (!goodNetwork)
168               break;
169     }
170     if (!goodNetwork) {
171          delete [] indices_;
172          // put in message
173          printf("Not a network - can test if indices_ null\n");
174          indices_ = NULL;
175          numberRows_ = 0;
176          numberColumns_ = 0;
177     } else {
178          numberRows_ ++; //  correct
179          trueNetwork_ = goodNetwork > 0;
180     }
181}
182
183//-------------------------------------------------------------------
184// Destructor
185//-------------------------------------------------------------------
186ClpNetworkMatrix::~ClpNetworkMatrix ()
187{
188     delete matrix_;
189     delete [] lengths_;
190     delete [] indices_;
191}
192
193//----------------------------------------------------------------
194// Assignment operator
195//-------------------------------------------------------------------
196ClpNetworkMatrix &
197ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
198{
199     if (this != &rhs) {
200          ClpMatrixBase::operator=(rhs);
201          delete matrix_;
202          delete [] lengths_;
203          delete [] indices_;
204          matrix_ = NULL;
205          lengths_ = NULL;
206          indices_ = NULL;
207          numberRows_ = rhs.numberRows_;
208          numberColumns_ = rhs.numberColumns_;
209          trueNetwork_ = rhs.trueNetwork_;
210          if (numberColumns_) {
211               indices_ = new int [ 2*numberColumns_];
212               CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
213          }
214     }
215     return *this;
216}
217//-------------------------------------------------------------------
218// Clone
219//-------------------------------------------------------------------
220ClpMatrixBase * ClpNetworkMatrix::clone() const
221{
222     return new ClpNetworkMatrix(*this);
223}
224
225/* Returns a new matrix in reverse order without gaps */
226ClpMatrixBase *
227ClpNetworkMatrix::reverseOrderedCopy() const
228{
229     // count number in each row
230     CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
231     CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
232     memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
233     memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
234     CoinBigIndex j = 0;
235     int i;
236     for (i = 0; i < numberColumns_; i++, j += 2) {
237          int iRow = indices_[j];
238          tempN[iRow]++;
239          iRow = indices_[j+1];
240          tempP[iRow]++;
241     }
242     int * newIndices = new int [2*numberColumns_];
243     CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
244     CoinBigIndex * newN = new CoinBigIndex[numberRows_];
245     int iRow;
246     j = 0;
247     // do starts
248     for (iRow = 0; iRow < numberRows_; iRow++) {
249          newP[iRow] = j;
250          j += tempP[iRow];
251          tempP[iRow] = newP[iRow];
252          newN[iRow] = j;
253          j += tempN[iRow];
254          tempN[iRow] = newN[iRow];
255     }
256     newP[numberRows_] = j;
257     j = 0;
258     for (i = 0; i < numberColumns_; i++, j += 2) {
259          int iRow = indices_[j];
260          CoinBigIndex put = tempN[iRow];
261          newIndices[put++] = i;
262          tempN[iRow] = put;
263          iRow = indices_[j+1];
264          put = tempP[iRow];
265          newIndices[put++] = i;
266          tempP[iRow] = put;
267     }
268     delete [] tempP;
269     delete [] tempN;
270     ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
271     newCopy->passInCopy(numberRows_, numberColumns_,
272                         false,  newIndices, newP, newN);
273     return newCopy;
274}
275//unscaled versions
276void
277ClpNetworkMatrix::times(double scalar,
278                        const double * x, double * y) const
279{
280     int iColumn;
281     CoinBigIndex j = 0;
282     if (trueNetwork_) {
283          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
284               double value = scalar * x[iColumn];
285               if (value) {
286                    int iRowM = indices_[j];
287                    int iRowP = indices_[j+1];
288                    y[iRowM] -= value;
289                    y[iRowP] += value;
290               }
291          }
292     } else {
293          // skip negative rows
294          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
295               double value = scalar * x[iColumn];
296               if (value) {
297                    int iRowM = indices_[j];
298                    int iRowP = indices_[j+1];
299                    if (iRowM >= 0)
300                         y[iRowM] -= value;
301                    if (iRowP >= 0)
302                         y[iRowP] += value;
303               }
304          }
305     }
306}
307void
308ClpNetworkMatrix::transposeTimes(double scalar,
309                                 const double * x, double * y) const
310{
311     int iColumn;
312     CoinBigIndex j = 0;
313     if (trueNetwork_) {
314          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
315               double value = y[iColumn];
316               int iRowM = indices_[j];
317               int iRowP = indices_[j+1];
318               value -= scalar * x[iRowM];
319               value += scalar * x[iRowP];
320               y[iColumn] = value;
321          }
322     } else {
323          // skip negative rows
324          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
325               double value = y[iColumn];
326               int iRowM = indices_[j];
327               int iRowP = indices_[j+1];
328               if (iRowM >= 0)
329                    value -= scalar * x[iRowM];
330               if (iRowP >= 0)
331                    value += scalar * x[iRowP];
332               y[iColumn] = value;
333          }
334     }
335}
336void
337ClpNetworkMatrix::times(double scalar,
338                        const double * x, double * y,
339                        const double * /*rowScale*/,
340                        const double * /*columnScale*/) const
341{
342     // we know it is not scaled
343     times(scalar, x, y);
344}
345void
346ClpNetworkMatrix::transposeTimes( double scalar,
347                                  const double * x, double * y,
348                                  const double * /*rowScale*/,
349                                  const double * /*columnScale*/,
350                                  double * /*spare*/) const
351{
352     // we know it is not scaled
353     transposeTimes(scalar, x, y);
354}
355/* Return <code>x * A + y</code> in <code>z</code>.
356        Squashes small elements and knows about ClpSimplex */
357void
358ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
359                                 const CoinIndexedVector * rowArray,
360                                 CoinIndexedVector * y,
361                                 CoinIndexedVector * columnArray) const
362{
363     // we know it is not scaled
364     columnArray->clear();
365     double * pi = rowArray->denseVector();
366     int numberNonZero = 0;
367     int * index = columnArray->getIndices();
368     double * array = columnArray->denseVector();
369     int numberInRowArray = rowArray->getNumElements();
370     // maybe I need one in OsiSimplex
371     double zeroTolerance = model->zeroTolerance();
372     int numberRows = model->numberRows();
373#ifndef NO_RTTI
374     ClpPlusMinusOneMatrix* rowCopy =
375          dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
376#else
377     ClpPlusMinusOneMatrix* rowCopy =
378          static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
379#endif
380     bool packed = rowArray->packedMode();
381     double factor = 0.3;
382     // We may not want to do by row if there may be cache problems
383     int numberColumns = model->numberColumns();
384     // It would be nice to find L2 cache size - for moment 512K
385     // Be slightly optimistic
386     if (numberColumns * sizeof(double) > 1000000) {
387          if (numberRows * 10 < numberColumns)
388               factor = 0.1;
389          else if (numberRows * 4 < numberColumns)
390               factor = 0.15;
391          else if (numberRows * 2 < numberColumns)
392               factor = 0.2;
393          //if (model->numberIterations()%50==0)
394          //printf("%d nonzero\n",numberInRowArray);
395     }
396     if (numberInRowArray > factor * numberRows || !rowCopy) {
397          // do by column
398          int iColumn;
399          assert (!y->getNumElements());
400          CoinBigIndex j = 0;
401          if (packed) {
402               // need to expand pi into y
403               assert(y->capacity() >= numberRows);
404               double * piOld = pi;
405               pi = y->denseVector();
406               const int * whichRow = rowArray->getIndices();
407               int i;
408               // modify pi so can collapse to one loop
409               for (i = 0; i < numberInRowArray; i++) {
410                    int iRow = whichRow[i];
411                    pi[iRow] = scalar * piOld[i];
412               }
413               if (trueNetwork_) {
414                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
415                         double value = 0.0;
416                         int iRowM = indices_[j];
417                         int iRowP = indices_[j+1];
418                         value -= pi[iRowM];
419                         value += pi[iRowP];
420                         if (fabs(value) > zeroTolerance) {
421                              array[numberNonZero] = value;
422                              index[numberNonZero++] = iColumn;
423                         }
424                    }
425               } else {
426                    // skip negative rows
427                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
428                         double value = 0.0;
429                         int iRowM = indices_[j];
430                         int iRowP = indices_[j+1];
431                         if (iRowM >= 0)
432                              value -= pi[iRowM];
433                         if (iRowP >= 0)
434                              value += pi[iRowP];
435                         if (fabs(value) > zeroTolerance) {
436                              array[numberNonZero] = value;
437                              index[numberNonZero++] = iColumn;
438                         }
439                    }
440               }
441               for (i = 0; i < numberInRowArray; i++) {
442                    int iRow = whichRow[i];
443                    pi[iRow] = 0.0;
444               }
445          } else {
446               if (trueNetwork_) {
447                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
448                         double value = 0.0;
449                         int iRowM = indices_[j];
450                         int iRowP = indices_[j+1];
451                         value -= scalar * pi[iRowM];
452                         value += scalar * pi[iRowP];
453                         if (fabs(value) > zeroTolerance) {
454                              index[numberNonZero++] = iColumn;
455                              array[iColumn] = value;
456                         }
457                    }
458               } else {
459                    // skip negative rows
460                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
461                         double value = 0.0;
462                         int iRowM = indices_[j];
463                         int iRowP = indices_[j+1];
464                         if (iRowM >= 0)
465                              value -= scalar * pi[iRowM];
466                         if (iRowP >= 0)
467                              value += scalar * pi[iRowP];
468                         if (fabs(value) > zeroTolerance) {
469                              index[numberNonZero++] = iColumn;
470                              array[iColumn] = value;
471                         }
472                    }
473               }
474          }
475          columnArray->setNumElements(numberNonZero);
476     } else {
477          // do by row
478          rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
479     }
480}
481/* Return <code>x *A in <code>z</code> but
482   just for indices in y. */
483void
484ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
485                                       const CoinIndexedVector * rowArray,
486                                       const CoinIndexedVector * y,
487                                       CoinIndexedVector * columnArray) const
488{
489     columnArray->clear();
490     double * pi = rowArray->denseVector();
491     double * array = columnArray->denseVector();
492     int jColumn;
493     int numberToDo = y->getNumElements();
494     const int * which = y->getIndices();
495     assert (!rowArray->packedMode());
496     columnArray->setPacked();
497     if (trueNetwork_) {
498          for (jColumn = 0; jColumn < numberToDo; jColumn++) {
499               int iColumn = which[jColumn];
500               double value = 0.0;
501               CoinBigIndex j = iColumn << 1;
502               int iRowM = indices_[j];
503               int iRowP = indices_[j+1];
504               value -= pi[iRowM];
505               value += pi[iRowP];
506               array[jColumn] = value;
507          }
508     } else {
509          // skip negative rows
510          for (jColumn = 0; jColumn < numberToDo; jColumn++) {
511               int iColumn = which[jColumn];
512               double value = 0.0;
513               CoinBigIndex j = iColumn << 1;
514               int iRowM = indices_[j];
515               int iRowP = indices_[j+1];
516               if (iRowM >= 0)
517                    value -= pi[iRowM];
518               if (iRowP >= 0)
519                    value += pi[iRowP];
520               array[jColumn] = value;
521          }
522     }
523}
524/// returns number of elements in column part of basis,
525CoinBigIndex
526ClpNetworkMatrix::countBasis( const int * whichColumn,
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 = static_cast<int> (startFraction * numberColumns_);
881     int end = CoinMin(static_cast<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.