source: branches/pre/ClpPlusMinusOneMatrix.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago

new stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.9 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5#include <cstdio>
6
7#include "CoinPragma.hpp"
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10
11#include "ClpSimplex.hpp"
12#include "ClpFactorization.hpp"
13// at end to get min/max!
14#include "ClpPlusMinusOneMatrix.hpp"
15#include "ClpMessage.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix () 
25  : ClpMatrixBase()
26{
27  setType(12);
28  elements_ = NULL;
29  startPositive_ = NULL;
30  startNegative_ = NULL;
31  lengths_=NULL;
32  indices_=NULL;
33  numberRows_=0;
34  numberColumns_=0;
35  columnOrdered_=true;
36}
37
38//-------------------------------------------------------------------
39// Copy constructor
40//-------------------------------------------------------------------
41ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs) 
42: ClpMatrixBase(rhs)
43{ 
44  elements_ = NULL;
45  startPositive_ = NULL;
46  startNegative_ = NULL;
47  lengths_=NULL;
48  indices_=NULL;
49  numberRows_=rhs.numberRows_;
50  numberColumns_=rhs.numberColumns_;
51  columnOrdered_=rhs.columnOrdered_;
52  if (numberColumns_) {
53    int numberElements = rhs.startPositive_[numberColumns_];
54    indices_ = new int [ numberElements];
55    memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
56    startPositive_ = new int [ numberColumns_+1];
57    memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
58    startNegative_ = new int [ numberColumns_];
59    memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int));
60  }
61}
62
63ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs) 
64  : ClpMatrixBase()
65{ 
66  setType(12);
67  elements_ = NULL;
68  startPositive_ = NULL;
69  startNegative_ = NULL;
70  lengths_=NULL;
71  indices_=NULL;
72  int iColumn;
73  assert (rhs.isColOrdered());
74  // get matrix data pointers
75  const int * row = rhs.getIndices();
76  const CoinBigIndex * columnStart = rhs.getVectorStarts();
77  const int * columnLength = rhs.getVectorLengths(); 
78  const double * elementByColumn = rhs.getElements();
79  numberColumns_ = rhs.getNumCols();
80  bool goodPlusMinusOne=true;
81  numberRows_=-1;
82  indices_ = new int[rhs.getNumElements()];
83  startPositive_ = new int [numberColumns_+1];
84  startNegative_ = new int [numberColumns_];
85  int * temp = new int [rhs.getNumRows()];
86  CoinBigIndex j=0;
87  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
88    CoinBigIndex k;
89    int iNeg=0;
90    startPositive_[iColumn]=j;
91    for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
92         k++) {
93      int iRow;
94      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
95        iRow = row[k];
96        numberRows_ = max(numberRows_,iRow);
97        indices_[j++]=iRow;
98      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
99        iRow = row[k];
100        numberRows_ = max(numberRows_,iRow);
101        temp[iNeg++]=iRow;
102      } else {
103        goodPlusMinusOne = false; // not a network
104      }
105    }
106    if (goodPlusMinusOne) {
107      // move negative
108      startNegative_[iColumn]=j;
109      for (k=0;k<iNeg;k++) {
110        indices_[j++] = temp[k];
111      }
112    } else {
113      break;
114    }
115  }
116  startPositive_[numberColumns_]=j;
117  delete [] temp;
118  if (!goodPlusMinusOne) {
119    delete [] indices_;
120    // put in message
121    printf("Not all +-1 - can test if indices_ null\n");
122    indices_=NULL;
123    numberRows_=0;
124    numberColumns_=0;
125    delete [] startPositive_;
126    delete [] startNegative_;
127    startPositive_ = NULL;
128    startNegative_ = NULL;
129  } else {
130    numberRows_ ++; //  correct
131    columnOrdered_ = true;
132  }
133}
134
135//-------------------------------------------------------------------
136// Destructor
137//-------------------------------------------------------------------
138ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
139{
140  delete [] elements_;
141  delete [] startPositive_;
142  delete [] startNegative_;
143  delete [] lengths_;
144  delete [] indices_;
145}
146
147//----------------------------------------------------------------
148// Assignment operator
149//-------------------------------------------------------------------
150ClpPlusMinusOneMatrix &
151ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
152{
153  if (this != &rhs) {
154    ClpMatrixBase::operator=(rhs);
155    delete [] elements_;
156    delete [] startPositive_;
157    delete [] startNegative_;
158    delete [] lengths_;
159    delete [] indices_;
160    elements_ = NULL;
161    startPositive_ = NULL;
162    lengths_=NULL;
163    indices_=NULL;
164    numberRows_=rhs.numberRows_;
165    numberColumns_=rhs.numberColumns_;
166    columnOrdered_=rhs.columnOrdered_;
167    if (numberColumns_) {
168      int numberElements = rhs.startPositive_[numberColumns_];
169      indices_ = new int [ numberElements];
170      memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
171      startPositive_ = new int [ numberColumns_+1];
172      memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
173      startNegative_ = new int [ numberColumns_];
174      memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int));
175    }
176  }
177  return *this;
178}
179//-------------------------------------------------------------------
180// Clone
181//-------------------------------------------------------------------
182ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
183{
184  return new ClpPlusMinusOneMatrix(*this);
185}
186/* Subset clone (without gaps).  Duplicates are allowed
187   and order is as given */
188ClpMatrixBase * 
189ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
190                              int numberColumns, 
191                              const int * whichColumns) const 
192{
193  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
194                                   numberColumns, whichColumns);
195}
196/* Subset constructor (without gaps).  Duplicates are allowed
197   and order is as given */
198ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
199                       const ClpPlusMinusOneMatrix & rhs,
200                       int numberRows, const int * whichRow,
201                       int numberColumns, const int * whichColumn)
202: ClpMatrixBase(rhs)
203{ 
204  elements_ = NULL;
205  startPositive_ = NULL;
206  startNegative_ = NULL;
207  lengths_=NULL;
208  indices_=NULL;
209  numberRows_=0;
210  numberColumns_=0;
211  columnOrdered_=rhs.columnOrdered_;
212  if (numberRows<=0||numberColumns<=0) {
213    startPositive_ = new int[1];
214    startPositive_[0] = 0;
215  } else {
216    numberColumns_ = numberColumns;
217    numberRows_ = numberRows;
218    const int * index1 = rhs.indices_;
219    int * startPositive1 = rhs.startPositive_;
220
221    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
222    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
223    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
224    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
225    // Throw exception if rhs empty
226    if (numberMajor1 <= 0 || numberMinor1 <= 0)
227      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
228    // Array to say if an old row is in new copy
229    int * newRow = new int [numberMinor1];
230    int iRow;
231    for (iRow=0;iRow<numberMinor1;iRow++) 
232      newRow[iRow] = -1;
233    // and array for duplicating rows
234    int * duplicateRow = new int [numberMinor];
235    int numberBad=0;
236    for (iRow=0;iRow<numberMinor;iRow++) {
237      duplicateRow[iRow] = -1;
238      int kRow = whichRow[iRow];
239      if (kRow>=0  && kRow < numberMinor1) {
240        if (newRow[kRow]<0) {
241          // first time
242          newRow[kRow]=iRow;
243        } else {
244          // duplicate
245          int lastRow = newRow[kRow];
246          newRow[kRow]=iRow;
247          duplicateRow[iRow] = lastRow;
248        }
249      } else {
250        // bad row
251        numberBad++;
252      }
253    }
254
255    if (numberBad)
256      throw CoinError("bad minor entries", 
257                      "subset constructor", "ClpPlusMinusOneMatrix");
258    // now get size and check columns
259    int size = 0;
260    int iColumn;
261    numberBad=0;
262    for (iColumn=0;iColumn<numberMajor;iColumn++) {
263      int kColumn = whichColumn[iColumn];
264      if (kColumn>=0  && kColumn <numberMajor1) {
265        int i;
266        for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
267          int kRow = index1[i];
268          kRow = newRow[kRow];
269          while (kRow>=0) {
270            size++;
271            kRow = duplicateRow[kRow];
272          }
273        }
274      } else {
275        // bad column
276        numberBad++;
277      }
278    }
279    if (numberBad)
280      throw CoinError("bad major entries", 
281                      "subset constructor", "ClpPlusMinusOneMatrix");
282    // now create arrays
283    startPositive_ = new int [numberMajor+1];
284    startNegative_ = new int [numberMajor];
285    indices_ = new int[size];
286    // and fill them
287    size = 0;
288    startPositive_[0]=0;
289    int * startNegative1 = rhs.startNegative_;
290    for (iColumn=0;iColumn<numberMajor;iColumn++) {
291      int kColumn = whichColumn[iColumn];
292      int i;
293      for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
294        int kRow = index1[i];
295        kRow = newRow[kRow];
296        while (kRow>=0) {
297          indices_[size++] = kRow;
298          kRow = duplicateRow[kRow];
299        }
300      }
301      startNegative_[iColumn] = size;
302      for (;i<startPositive1[kColumn+1];i++) {
303        int kRow = index1[i];
304        kRow = newRow[kRow];
305        while (kRow>=0) {
306          indices_[size++] = kRow;
307          kRow = duplicateRow[kRow];
308        }
309      }
310      startPositive_[iColumn+1] = size;
311    }
312    delete [] newRow;
313    delete [] duplicateRow;
314  }
315}
316
317
318/* Returns a new matrix in reverse order without gaps */
319ClpMatrixBase * 
320ClpPlusMinusOneMatrix::reverseOrderedCopy() const
321{
322  int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
323  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
324  // count number in each row/column
325  int * tempP = new int [numberMinor];
326  int * tempN = new int [numberMinor];
327  memset(tempP,0,numberMinor*sizeof(int));
328  memset(tempN,0,numberMinor*sizeof(int));
329  CoinBigIndex j=0;
330  int i;
331  for (i=0;i<numberMajor;i++) {
332    for (;j<startNegative_[i];j++) {
333      int iRow = indices_[j];
334      tempP[iRow]++;
335    }
336    for (;j<startPositive_[i+1];j++) {
337      int iRow = indices_[j];
338      tempN[iRow]++;
339    }
340  }
341  int * newIndices = new int [startPositive_[numberMajor]];
342  int * newP = new int [numberMinor+1];
343  int * newN = new int[numberMinor];
344  int iRow;
345  j=0;
346  // do starts
347  for (iRow=0;iRow<numberMinor;iRow++) {
348    newP[iRow]=j;
349    j += tempP[iRow];
350    tempP[iRow]=newP[iRow];
351    newN[iRow] = j;
352    j += tempN[iRow];
353    tempN[iRow]=newN[iRow];
354  }
355  newP[numberMinor]=j;
356  j=0;
357  for (i=0;i<numberMajor;i++) {
358    for (;j<startNegative_[i];j++) {
359      int iRow = indices_[j];
360      int put = tempP[iRow];
361      newIndices[put++] = i;
362      tempP[iRow] = put;
363    }
364    for (;j<startPositive_[i+1];j++) {
365      int iRow = indices_[j];
366      int put = tempN[iRow];
367      newIndices[put++] = i;
368      tempN[iRow] = put;
369    }
370  }
371  delete [] tempP;
372  delete [] tempN;
373  ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
374  newCopy->passInCopy(numberMinor, numberMajor,
375                      !columnOrdered_,  newIndices, newP, newN);
376  return newCopy;
377}
378//unscaled versions
379void 
380ClpPlusMinusOneMatrix::times(double scalar,
381                   const double * x, double * y) const
382{
383  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
384  int i;
385  CoinBigIndex j;
386  assert (columnOrdered_);
387  for (i=0;i<numberMajor;i++) {
388    double value = scalar*x[i];
389    if (value) {
390      for (j=startPositive_[i];j<startNegative_[i];j++) {
391        int iRow = indices_[j];
392        y[iRow] += value;
393      }
394      for (;j<startPositive_[i+1];j++) {
395        int iRow = indices_[j];
396        y[iRow] -= value;
397      }
398    }
399  }
400}
401void 
402ClpPlusMinusOneMatrix::transposeTimes(double scalar,
403                                const double * x, double * y) const
404{
405  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
406  int i;
407  CoinBigIndex j=0;
408  assert (columnOrdered_);
409  for (i=0;i<numberMajor;i++) {
410    double value = 0.0;
411    for (;j<startNegative_[i];j++) {
412      int iRow = indices_[j];
413      value += x[iRow];
414    }
415    for (;j<startPositive_[i+1];j++) {
416      int iRow = indices_[j];
417      value -= x[iRow];
418    }
419    y[i] += scalar*value;
420  }
421}
422void 
423ClpPlusMinusOneMatrix::times(double scalar,
424                       const double * x, double * y,
425                       const double * rowScale, 
426                       const double * columnScale) const
427{
428  // we know it is not scaled
429  times(scalar, x, y);
430}
431void 
432ClpPlusMinusOneMatrix::transposeTimes( double scalar,
433                                 const double * x, double * y,
434                                 const double * rowScale, 
435                                 const double * columnScale) const
436{
437  // we know it is not scaled
438  transposeTimes(scalar, x, y);
439}
440/* Return <code>x * A + y</code> in <code>z</code>.
441        Squashes small elements and knows about ClpSimplex */
442void 
443ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
444                              const CoinIndexedVector * rowArray,
445                              CoinIndexedVector * y,
446                              CoinIndexedVector * columnArray) const
447{
448  // we know it is not scaled
449  columnArray->clear();
450  double * pi = rowArray->denseVector();
451  int numberNonZero=0;
452  int * index = columnArray->getIndices();
453  double * array = columnArray->denseVector();
454  int numberInRowArray = rowArray->getNumElements();
455  // maybe I need one in OsiSimplex
456  double zeroTolerance = model->factorization()->zeroTolerance();
457  int numberRows = model->numberRows();
458  bool packed = rowArray->packedMode();
459  ClpPlusMinusOneMatrix* rowCopy =
460    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
461  if (numberInRowArray>0.3*numberRows||!rowCopy) {
462    assert (!y->getNumElements());
463    // do by column
464    // Need to expand if packed mode
465    int iColumn;
466    CoinBigIndex j=0;
467    assert (columnOrdered_);
468    if (packed) {
469      // need to expand pi into y
470      assert(y->capacity()>=numberRows);
471      double * piOld = pi;
472      pi = y->denseVector();
473      const int * whichRow = rowArray->getIndices();
474      int i;
475      // modify pi so can collapse to one loop
476      for (i=0;i<numberInRowArray;i++) {
477        int iRow = whichRow[i];
478        pi[iRow]=scalar*piOld[i];
479      }
480      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
481        double value = 0.0;
482        for (;j<startNegative_[iColumn];j++) {
483          int iRow = indices_[j];
484          value += pi[iRow];
485        }
486        for (;j<startPositive_[iColumn+1];j++) {
487          int iRow = indices_[j];
488          value -= pi[iRow];
489        }
490        if (fabs(value)>zeroTolerance) {
491          array[numberNonZero]=value;
492          index[numberNonZero++]=iColumn;
493        }
494      }
495      for (i=0;i<numberInRowArray;i++) {
496        int iRow = whichRow[i];
497        pi[iRow]=0.0;
498      }
499    } else {
500      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
501        double value = 0.0;
502        for (;j<startNegative_[iColumn];j++) {
503          int iRow = indices_[j];
504          value += pi[iRow];
505        }
506        for (;j<startPositive_[iColumn+1];j++) {
507          int iRow = indices_[j];
508          value -= pi[iRow];
509        }
510        value += scalar;
511        if (fabs(value)>zeroTolerance) {
512          index[numberNonZero++]=iColumn;
513          array[iColumn]=value;
514        }
515      }
516    }
517    columnArray->setNumElements(numberNonZero);
518  } else {
519    // do by row
520    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
521  }
522}
523/* Return <code>x * A + y</code> in <code>z</code>.
524        Squashes small elements and knows about ClpSimplex */
525void 
526ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
527                              const CoinIndexedVector * rowArray,
528                              CoinIndexedVector * y,
529                              CoinIndexedVector * columnArray) const
530{
531  columnArray->clear();
532  double * pi = rowArray->denseVector();
533  int numberNonZero=0;
534  int * index = columnArray->getIndices();
535  double * array = columnArray->denseVector();
536  int numberInRowArray = rowArray->getNumElements();
537  // maybe I need one in OsiSimplex
538  double zeroTolerance = model->factorization()->zeroTolerance();
539  const int * column = indices_;
540  const CoinBigIndex * startPositive = startPositive_;
541  const CoinBigIndex * startNegative = startNegative_;
542  const int * whichRow = rowArray->getIndices();
543  bool packed = rowArray->packedMode();
544  if (numberInRowArray>2||y->getNumElements()) {
545    // do by rows
546    int iRow;
547    double * markVector = y->denseVector(); // probably empty .. but
548    int * mark = y->getIndices();
549    int numberOriginal=y->getNumElements();
550    int i;
551    if (packed) {
552      assert(!numberOriginal);
553      numberNonZero=0;
554      // and set up mark as char array
555      char * marked = (char *) (index+columnArray->capacity());
556      double * array2 = y->denseVector();
557#ifdef CLP_DEBUG
558      int numberColumns = model->numberColumns();
559      for (i=0;i<numberColumns;i++) {
560        assert(!marked[i]);
561        assert(!array2[i]);
562      }
563#endif
564      for (i=0;i<numberInRowArray;i++) {
565        iRow = whichRow[i]; 
566        double value = pi[i]*scalar;
567        CoinBigIndex j;
568        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
569          int iColumn = column[j];
570          if (!marked[iColumn]) {
571            marked[iColumn]=1;
572            index[numberNonZero++]=iColumn;
573          }
574          array2[iColumn] += value;
575        }
576        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
577          int iColumn = column[j];
578          if (!marked[iColumn]) {
579            marked[iColumn]=1;
580            index[numberNonZero++]=iColumn;
581          }
582          array2[iColumn] -= value;
583        }
584      }
585      // get rid of tiny values and zero out marked
586      numberOriginal=numberNonZero;
587      numberNonZero=0;
588      for (i=0;i<numberOriginal;i++) {
589        int iColumn = index[i];
590        if (marked[iColumn]) {
591          double value = array2[iColumn];
592          array2[iColumn]=0.0;
593          marked[iColumn]=0;
594          if (fabs(value)>zeroTolerance) {
595            array[numberNonZero]=value;
596            index[numberNonZero++]=iColumn;
597          }
598        }
599      }
600    } else {     
601      for (i=0;i<numberOriginal;i++) {
602        int iColumn = mark[i];
603        index[i]=iColumn;
604        array[iColumn]=markVector[iColumn];
605        markVector[iColumn]=0.0;
606      }
607      numberNonZero=numberOriginal;
608      // and set up mark as char array
609      char * marked = (char *) markVector;
610      for (i=0;i<numberOriginal;i++) {
611        int iColumn = index[i];
612        marked[iColumn]=0;
613      }
614      for (i=0;i<numberInRowArray;i++) {
615        iRow = whichRow[i]; 
616        double value = pi[iRow]*scalar;
617        CoinBigIndex j;
618        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
619          int iColumn = column[j];
620          if (!marked[iColumn]) {
621            marked[iColumn]=1;
622            index[numberNonZero++]=iColumn;
623          }
624          array[iColumn] += value;
625        }
626        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
627          int iColumn = column[j];
628          if (!marked[iColumn]) {
629            marked[iColumn]=1;
630            index[numberNonZero++]=iColumn;
631          }
632          array[iColumn] -= value;
633        }
634      }
635      // get rid of tiny values and zero out marked
636      numberOriginal=numberNonZero;
637      numberNonZero=0;
638      for (i=0;i<numberOriginal;i++) {
639        int iColumn = index[i];
640        marked[iColumn]=0;
641        if (fabs(array[iColumn])>zeroTolerance) {
642          index[numberNonZero++]=iColumn;
643        } else {
644          array[iColumn]=0.0;
645        }
646      }
647    }
648  } else if (numberInRowArray==2) {
649    /* do by rows when two rows (do longer first when not packed
650       and shorter first if packed */
651    int iRow0 = whichRow[0];
652    int iRow1 = whichRow[1];
653    int j;
654    if (packed) {
655      double pi0 = pi[0];
656      double pi1 = pi[1];
657      if (startPositive[iRow0+1]-startPositive[iRow0]>
658          startPositive[iRow1+1]-startPositive[iRow1]) {
659        int temp = iRow0;
660        iRow0 = iRow1;
661        iRow1 = temp;
662        pi0=pi1;
663        pi1=pi[0];
664      }
665      // and set up mark as char array
666      char * marked = (char *) (index+columnArray->capacity());
667      int * lookup = y->getIndices();
668      double value = pi0*scalar;
669      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
670        int iColumn = column[j];
671        array[numberNonZero] = value;
672        marked[iColumn]=1;
673        lookup[iColumn]=numberNonZero;
674        index[numberNonZero++]=iColumn;
675      }
676      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
677        int iColumn = column[j];
678        array[numberNonZero] = -value;
679        marked[iColumn]=1;
680        lookup[iColumn]=numberNonZero;
681        index[numberNonZero++]=iColumn;
682      }
683      int numberOriginal = numberNonZero;
684      value = pi1*scalar;
685      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
686        int iColumn = column[j];
687        if (marked[iColumn]) {
688          int iLookup = lookup[iColumn];
689          array[iLookup] += value;
690        } else {
691          if (fabs(value)>zeroTolerance) {
692            array[numberNonZero] = value;
693            index[numberNonZero++]=iColumn;
694          }
695        }
696      }
697      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
698        int iColumn = column[j];
699        if (marked[iColumn]) {
700          int iLookup = lookup[iColumn];
701          array[iLookup] -= value;
702        } else {
703          if (fabs(value)>zeroTolerance) {
704            array[numberNonZero] = -value;
705            index[numberNonZero++]=iColumn;
706          }
707        }
708      }
709      // get rid of tiny values and zero out marked
710      int nDelete=0;
711      for (j=0;j<numberOriginal;j++) {
712        int iColumn = index[j];
713        marked[iColumn]=0;
714        if (fabs(array[j])<=zeroTolerance) 
715          nDelete++;
716      }
717      if (nDelete) {
718        numberOriginal=numberNonZero;
719        numberNonZero=0;
720        for (j=0;j<numberOriginal;j++) {
721          int iColumn = index[j];
722          double value = array[j];
723          array[j]=0.0;
724          if (fabs(value)>zeroTolerance) {
725            array[numberNonZero]=value;
726            index[numberNonZero++]=iColumn;
727          }
728        }
729      }
730    } else {
731      if (startPositive[iRow0+1]-startPositive[iRow0]<
732          startPositive[iRow1+1]-startPositive[iRow1]) {
733        int temp = iRow0;
734        iRow0 = iRow1;
735        iRow1 = temp;
736      }
737      int numberOriginal;
738      int i;
739      numberNonZero=0;
740      double value;
741      value = pi[iRow0]*scalar;
742      CoinBigIndex j;
743      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
744        int iColumn = column[j];
745        index[numberNonZero++]=iColumn;
746        array[iColumn] = value;
747      }
748      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
749        int iColumn = column[j];
750        index[numberNonZero++]=iColumn;
751        array[iColumn] = -value;
752      }
753      value = pi[iRow1]*scalar;
754      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
755        int iColumn = column[j];
756        double value2= array[iColumn];
757        if (value2) {
758          value2 += value;
759        } else {
760          value2 = value;
761          index[numberNonZero++]=iColumn;
762        }
763        array[iColumn] = value2;
764      }
765      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
766        int iColumn = column[j];
767        double value2= array[iColumn];
768        if (value2) {
769          value2 -= value;
770        } else {
771          value2 = -value;
772          index[numberNonZero++]=iColumn;
773        }
774        array[iColumn] = value2;
775      }
776      // get rid of tiny values and zero out marked
777      numberOriginal=numberNonZero;
778      numberNonZero=0;
779      for (i=0;i<numberOriginal;i++) {
780        int iColumn = index[i];
781        if (fabs(array[iColumn])>zeroTolerance) {
782          index[numberNonZero++]=iColumn;
783        } else {
784          array[iColumn]=0.0;
785        }
786      }
787    }
788  } else if (numberInRowArray==1) {
789    // Just one row
790    int iRow=rowArray->getIndices()[0];
791    numberNonZero=0;
792    double value;
793    iRow = whichRow[0]; 
794    CoinBigIndex j;
795    if (packed) {
796      value = pi[0]*scalar;
797      if (fabs(value)>zeroTolerance) {
798        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
799          int iColumn = column[j];
800          array[numberNonZero] = value;
801          index[numberNonZero++]=iColumn;
802        }
803        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
804          int iColumn = column[j];
805          array[numberNonZero] = -value;
806          index[numberNonZero++]=iColumn;
807        }
808      }
809    } else {
810      value = pi[iRow]*scalar;
811      if (fabs(value)>zeroTolerance) {
812        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
813          int iColumn = column[j];
814          array[iColumn] = value;
815          index[numberNonZero++]=iColumn;
816        }
817        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
818          int iColumn = column[j];
819          array[iColumn] = -value;
820          index[numberNonZero++]=iColumn;
821        }
822      }
823    }
824  }
825  columnArray->setNumElements(numberNonZero);
826  y->setNumElements(0);
827}
828/* Return <code>x *A in <code>z</code> but
829   just for indices in y.
830   Squashes small elements and knows about ClpSimplex */
831void 
832ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * model,
833                              const CoinIndexedVector * rowArray,
834                              const CoinIndexedVector * y,
835                              CoinIndexedVector * columnArray) const
836{
837  columnArray->clear();
838  double * pi = rowArray->denseVector();
839  int numberNonZero=0;
840  int * index = columnArray->getIndices();
841  double * array = columnArray->denseVector();
842  // maybe I need one in OsiSimplex
843  double zeroTolerance = model->factorization()->zeroTolerance();
844  int jColumn;
845  int numberToDo = y->getNumElements();
846  const int * which = y->getIndices();
847  bool packed = rowArray->packedMode();
848  if (packed) {
849    // need to expand pi into y
850    int numberInRowArray = rowArray->getNumElements();
851    assert(y->capacity()>=model->numberRows());
852    double * piOld = pi;
853    pi = y->denseVector();
854    const int * whichRow = rowArray->getIndices();
855    int i;
856    for (i=0;i<numberInRowArray;i++) {
857      int iRow = whichRow[i];
858      pi[iRow]=piOld[i];
859    }
860    // Must line up with y
861    for (jColumn=0;jColumn<numberToDo;jColumn++) {
862      int iColumn = which[jColumn];
863      double value = 0.0;
864      CoinBigIndex j=startPositive_[iColumn];
865      for (;j<startNegative_[iColumn];j++) {
866        int iRow = indices_[j];
867        value += pi[iRow];
868      }
869      for (;j<startPositive_[iColumn+1];j++) {
870        int iRow = indices_[j];
871        value -= pi[iRow];
872      }
873      array[jColumn]=value;
874    }
875    for (i=0;i<numberInRowArray;i++) {
876      int iRow = whichRow[i];
877      pi[iRow]=0.0;
878    }
879  } else {
880    for (jColumn=0;jColumn<numberToDo;jColumn++) {
881      int iColumn = which[jColumn];
882      double value = 0.0;
883      CoinBigIndex j=startPositive_[iColumn];
884      for (;j<startNegative_[iColumn];j++) {
885        int iRow = indices_[j];
886        value += pi[iRow];
887      }
888      for (;j<startPositive_[iColumn+1];j++) {
889        int iRow = indices_[j];
890        value -= pi[iRow];
891      }
892      if (fabs(value)>zeroTolerance) {
893        index[numberNonZero++]=iColumn;
894        array[iColumn]=value;
895      }
896    }
897  }
898}
899/* Returns number of elements in basis
900   column is basic if entry >=0 */
901CoinBigIndex
902ClpPlusMinusOneMatrix::numberInBasis(const int * columnIsBasic) const 
903{
904  int i;
905  CoinBigIndex numberElements=0;
906  assert (columnOrdered_);
907  for (i=0;i<numberColumns_;i++) {
908    if (columnIsBasic[i]>=0) 
909      numberElements += startPositive_[i+1]-startPositive_[i];
910  }
911  return numberElements;
912}
913// Fills in basis (Returns number of elements and updates numberBasic)
914CoinBigIndex
915ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
916                                const int * columnIsBasic, int & numberBasic,
917                                int * indexRowU, int * indexColumnU,
918                                double * elementU) const 
919{
920#ifdef CLPDEBUG
921  const double * rowScale = model->rowScale();
922  assert (!rowScale);
923#endif
924  int i;
925  CoinBigIndex numberElements=0;
926  assert (columnOrdered_);
927  for (i=0;i<numberColumns_;i++) {
928    if (columnIsBasic[i]>=0) {
929      CoinBigIndex j=startPositive_[i];
930      for (;j<startNegative_[i];j++) {
931        int iRow = indices_[j];
932        indexRowU[numberElements]=iRow;
933        indexColumnU[numberElements]=numberBasic;
934        elementU[numberElements++]=1.0;
935      }
936      for (;j<startPositive_[i+1];j++) {
937        int iRow = indices_[j];
938        indexRowU[numberElements]=iRow;
939        indexColumnU[numberElements]=numberBasic;
940        elementU[numberElements++]=-1.0;
941      }
942      numberBasic++;
943    }
944  }
945  return numberElements;
946}
947/* If element NULL returns number of elements in column part of basis,
948   If not NULL fills in as well */
949CoinBigIndex
950ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
951                                 const int * whichColumn, 
952                                 int numberBasic,
953                                 int numberColumnBasic,
954                                 int * indexRowU, int * indexColumnU,
955                                 double * elementU) const 
956{
957  int i;
958  CoinBigIndex numberElements=0;
959  if (elementU!=NULL) {
960    assert (columnOrdered_);
961    for (i=0;i<numberColumnBasic;i++) {
962      int iColumn = whichColumn[i];
963      CoinBigIndex j=startPositive_[iColumn];
964      for (;j<startNegative_[iColumn];j++) {
965        int iRow = indices_[j];
966        indexRowU[numberElements]=iRow;
967        indexColumnU[numberElements]=numberBasic;
968        elementU[numberElements++]=1.0;
969      }
970      for (;j<startPositive_[iColumn+1];j++) {
971        int iRow = indices_[j];
972        indexRowU[numberElements]=iRow;
973        indexColumnU[numberElements]=numberBasic;
974        elementU[numberElements++]=-1.0;
975      }
976      numberBasic++;
977    }
978  } else {
979    for (i=0;i<numberColumnBasic;i++) {
980      int iColumn = whichColumn[i];
981      numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
982    }
983  }
984  return numberElements;
985}
986/* Unpacks a column into an CoinIndexedvector
987 */
988void 
989ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
990                              CoinIndexedVector * rowArray,
991                              int iColumn) const 
992{
993  CoinBigIndex j=startPositive_[iColumn];
994  for (;j<startNegative_[iColumn];j++) {
995    int iRow = indices_[j];
996    rowArray->add(iRow,1.0);
997  }
998  for (;j<startPositive_[iColumn+1];j++) {
999    int iRow = indices_[j];
1000    rowArray->add(iRow,-1.0);
1001  }
1002}
1003/* Unpacks a column into an CoinIndexedvector
1004** in packed foramt
1005Note that model is NOT const.  Bounds and objective could
1006be modified if doing column generation (just for this variable) */
1007void 
1008ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
1009                            CoinIndexedVector * rowArray,
1010                            int iColumn) const
1011{
1012  int * index = rowArray->getIndices();
1013  double * array = rowArray->denseVector();
1014  int number = 0;
1015  CoinBigIndex j=startPositive_[iColumn];
1016  for (;j<startNegative_[iColumn];j++) {
1017    int iRow = indices_[j];
1018    array[number]=1.0;
1019    index[number++]=iRow;
1020  }
1021  for (;j<startPositive_[iColumn+1];j++) {
1022    int iRow = indices_[j];
1023    array[number]=-1.0;
1024    index[number++]=iRow;
1025  }
1026  rowArray->setNumElements(number);
1027  rowArray->setPackedMode(true);
1028}
1029/* Adds multiple of a column into an CoinIndexedvector
1030      You can use quickAdd to add to vector */
1031void 
1032ClpPlusMinusOneMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1033                   int iColumn, double multiplier) const 
1034{
1035  CoinBigIndex j=startPositive_[iColumn];
1036  for (;j<startNegative_[iColumn];j++) {
1037    int iRow = indices_[j];
1038    rowArray->quickAdd(iRow,multiplier);
1039  }
1040  for (;j<startPositive_[iColumn+1];j++) {
1041    int iRow = indices_[j];
1042    rowArray->quickAdd(iRow,-multiplier);
1043  }
1044}
1045
1046// Return a complete CoinPackedMatrix
1047CoinPackedMatrix * 
1048ClpPlusMinusOneMatrix::getPackedMatrix() const 
1049{
1050  int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
1051  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1052  return new CoinPackedMatrix(columnOrdered_,numberMinor,numberMajor,
1053                              getNumElements(),
1054                              getElements(),indices_,
1055                              startPositive_,getVectorLengths());
1056
1057}
1058/* A vector containing the elements in the packed matrix. Note that there
1059   might be gaps in this list, entries that do not belong to any
1060   major-dimension vector. To get the actual elements one should look at
1061   this vector together with vectorStarts and vectorLengths. */
1062const double * 
1063ClpPlusMinusOneMatrix::getElements() const 
1064{
1065  if (!elements_) {
1066    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1067    int numberElements = startPositive_[numberMajor];
1068    elements_ = new double [numberElements];
1069    CoinBigIndex j=0;
1070    int i;
1071    for (i=0;i<numberMajor;i++) {
1072      for (;j<startNegative_[i];j++) {
1073        elements_[j]=1.0;
1074      }
1075      for (;j<startPositive_[i+1];j++) {
1076        elements_[j]=-1.0;
1077      }
1078    }
1079  }
1080  return elements_;
1081}
1082
1083const CoinBigIndex * 
1084ClpPlusMinusOneMatrix::getVectorStarts() const 
1085{
1086  return startPositive_;
1087}
1088/* The lengths of the major-dimension vectors. */
1089const int * 
1090ClpPlusMinusOneMatrix::getVectorLengths() const
1091{
1092  if (!lengths_) {
1093    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1094    lengths_ = new int [numberMajor];
1095    int i;
1096    for (i=0;i<numberMajor;i++) {
1097      lengths_[i]=startPositive_[i+1]-startPositive_[i];
1098    }
1099  }
1100  return lengths_;
1101}
1102/* Delete the columns whose indices are listed in <code>indDel</code>. */
1103void 
1104ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel) 
1105{
1106  abort();
1107}
1108/* Delete the rows whose indices are listed in <code>indDel</code>. */
1109void 
1110ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel) 
1111{
1112  abort();
1113}
1114bool 
1115ClpPlusMinusOneMatrix::isColOrdered() const 
1116{ 
1117  return columnOrdered_;
1118}
1119/* Number of entries in the packed matrix. */
1120CoinBigIndex
1121ClpPlusMinusOneMatrix::getNumElements() const 
1122{
1123  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1124  if (startPositive_) 
1125    return startPositive_[numberMajor];
1126  else
1127    return 0;
1128}
1129// pass in copy (object takes ownership)
1130void 
1131ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
1132                  bool columnOrdered, int * indices,
1133                  int * startPositive, int * startNegative)
1134{
1135  columnOrdered_=columnOrdered;
1136  startPositive_ = startPositive;
1137  startNegative_ = startNegative;
1138  indices_ = indices;
1139  numberRows_=numberRows;
1140  numberColumns_=numberColumns;
1141}
Note: See TracBrowser for help on using the repository browser.