source: branches/devel-1/ClpNetworkBasis.cpp @ 2351

Last change on this file since 2351 was 123, checked in by forrest, 17 years ago

What a mess

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.6 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include "ClpNetworkBasis.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "ClpSimplex.hpp"
8#include "ClpMatrixBase.hpp"
9#include "CoinIndexedVector.hpp"
10
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpNetworkBasis::ClpNetworkBasis () 
20{
21  slackValue_=-1.0;
22  numberRows_=0;
23  numberColumns_=0;
24  root_ = -1;
25  leaf_ = -1;
26  parent_ = NULL;
27  descendant_ = NULL;
28  pivot_ = NULL;
29  rightSibling_ = NULL;
30  leftSibling_ = NULL;
31  sign_ = NULL;
32  stack_ = NULL;
33  toLeaf_ = NULL;
34  toRoot_ = NULL;
35  mark_ = NULL;
36  model_=NULL;
37}
38// Constructor from CoinFactorization
39ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex * model,
40                                 int numberRows, const double * pivotRegion,
41                                 const int * permuteBack,
42                                 const int * startColumn, 
43                                 const int * numberInColumn,
44                                 const int * indexRow, const double * element)
45{
46  slackValue_=-1.0;
47  numberRows_=numberRows;
48  numberColumns_=numberRows;
49  parent_ = new int [ numberRows_+1];
50  descendant_ = new int [ numberRows_+1];
51  pivot_ = new int [ numberRows_+1];
52  rightSibling_ = new int [ numberRows_+1];
53  leftSibling_ = new int [ numberRows_+1];
54  sign_ = new double [ numberRows_+1];
55  stack_ = new int [ numberRows_+1];
56  toLeaf_ = new int [numberRows_+1];
57  toRoot_ = new int [numberRows_+1];
58  mark_ = new char[numberRows_+1];
59  int i;
60  for (i=0;i<numberRows_+1;i++) {
61    parent_[i]=-1;
62    descendant_[i]=-1;
63    pivot_[i]=-1;
64    rightSibling_[i]=-1;
65    leftSibling_[i]=-1;
66    sign_[i]=-1.0;
67    stack_[i]=-1;
68    mark_[i]=0;
69  }
70  // pivotColumnBack gives order of pivoting into basis
71  // so pivotColumnback[0] is first slack in basis and
72  // it pivots on row permuteBack[0]
73  // a known root is given by permuteBack[numberRows_-1]
74  root_ = numberRows_;
75  int lastPivot=numberRows_;
76  for (i=0;i<numberRows_;i++) {
77    int iPivot = permuteBack[i];
78    toRoot_[iPivot] = lastPivot;
79    toLeaf_[lastPivot]=iPivot;
80    lastPivot=iPivot;
81    double sign;
82    if (pivotRegion[i]>0.0)
83      sign = 1.0;
84    else
85      sign =-1.0;
86    int other;
87    if (numberInColumn[i]>0) {
88      int iRow = indexRow[startColumn[i]];
89      other = permuteBack[iRow];
90      assert (parent_[other]!=-1);
91    } else {
92      other = numberRows_;
93    }
94    sign_[iPivot] = sign;
95    int iParent = other;
96    parent_[iPivot] = other;
97    if (descendant_[iParent]>=0) {
98      // we have a sibling
99      int iRight = descendant_[iParent];
100      rightSibling_[iPivot]=iRight;
101      leftSibling_[iRight]=iPivot;
102    } else {
103      rightSibling_[iPivot]=-1;
104    }       
105    descendant_[iParent] = iPivot;
106    leftSibling_[iPivot]=-1;
107  }
108  toLeaf_[lastPivot]=numberRows_;
109  leaf_ = lastPivot;
110  model_=model;
111}
112
113//-------------------------------------------------------------------
114// Copy constructor
115//-------------------------------------------------------------------
116ClpNetworkBasis::ClpNetworkBasis (const ClpNetworkBasis & rhs) 
117{
118  slackValue_=rhs.slackValue_;
119  numberRows_=rhs.numberRows_;
120  numberColumns_=rhs.numberColumns_;
121  root_ = rhs.root_;
122  leaf_ = rhs.leaf_;
123  if (rhs.parent_) {
124    parent_ = new int [numberRows_+1];
125    memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
126  } else {
127    parent_ = NULL;
128  }
129  if (rhs.descendant_) {
130    descendant_ = new int [numberRows_+1];
131    memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
132  } else {
133    descendant_ = NULL;
134  }
135  if (rhs.pivot_) {
136    pivot_ = new int [numberRows_+1];
137    memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
138  } else {
139    pivot_ = NULL;
140  }
141  if (rhs.rightSibling_) {
142    rightSibling_ = new int [numberRows_+1];
143    memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
144  } else {
145    rightSibling_ = NULL;
146  }
147  if (rhs.leftSibling_) {
148    leftSibling_ = new int [numberRows_+1];
149    memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
150  } else {
151    leftSibling_ = NULL;
152  }
153  if (rhs.sign_) {
154    sign_ = new double [numberRows_+1];
155    memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
156  } else {
157    sign_ = NULL;
158  }
159  if (rhs.stack_) {
160    stack_ = new int [numberRows_+1];
161    memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
162  } else {
163    stack_ = NULL;
164  }
165  if (rhs.toLeaf_) {
166    toLeaf_ = new int [numberRows_+1];
167    memcpy(toLeaf_,rhs.toLeaf_,(numberRows_+1)*sizeof(int));
168  } else {
169    toLeaf_ = NULL;
170  }
171  if (rhs.toRoot_) {
172    toRoot_ = new int [numberRows_+1];
173    memcpy(toRoot_,rhs.toRoot_,(numberRows_+1)*sizeof(int));
174  } else {
175    toRoot_ = NULL;
176  }
177  if (rhs.mark_) {
178    mark_ = new char [numberRows_+1];
179    memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
180  } else {
181    mark_ = NULL;
182  }
183  model_=rhs.model_;
184}
185
186//-------------------------------------------------------------------
187// Destructor
188//-------------------------------------------------------------------
189ClpNetworkBasis::~ClpNetworkBasis () 
190{
191  delete [] parent_;
192  delete [] descendant_;
193  delete [] pivot_;
194  delete [] rightSibling_;
195  delete [] leftSibling_;
196  delete [] sign_;
197  delete [] stack_;
198  delete [] toLeaf_;
199  delete [] toRoot_;
200  delete [] mark_;
201}
202
203//----------------------------------------------------------------
204// Assignment operator
205//-------------------------------------------------------------------
206ClpNetworkBasis &
207ClpNetworkBasis::operator=(const ClpNetworkBasis& rhs)
208{
209  if (this != &rhs) {
210    delete [] parent_;
211    delete [] descendant_;
212    delete [] pivot_;
213    delete [] rightSibling_;
214    delete [] leftSibling_;
215    delete [] sign_;
216    delete [] stack_;
217    delete [] toLeaf_;
218    delete [] toRoot_;
219    delete [] mark_;
220    slackValue_=rhs.slackValue_;
221    numberRows_=rhs.numberRows_;
222    numberColumns_=rhs.numberColumns_;
223    root_ = rhs.root_;
224    leaf_ = rhs.leaf_;
225    if (rhs.parent_) {
226      parent_ = new int [numberRows_+1];
227      memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
228    } else {
229      parent_ = NULL;
230    }
231    if (rhs.descendant_) {
232      descendant_ = new int [numberRows_+1];
233      memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
234    } else {
235      descendant_ = NULL;
236    }
237    if (rhs.pivot_) {
238      pivot_ = new int [numberRows_+1];
239      memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
240    } else {
241      pivot_ = NULL;
242    }
243    if (rhs.rightSibling_) {
244      rightSibling_ = new int [numberRows_+1];
245      memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
246    } else {
247      rightSibling_ = NULL;
248    }
249    if (rhs.leftSibling_) {
250      leftSibling_ = new int [numberRows_+1];
251      memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
252    } else {
253      leftSibling_ = NULL;
254    }
255    if (rhs.sign_) {
256      sign_ = new double [numberRows_+1];
257      memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
258    } else {
259      sign_ = NULL;
260    }
261    if (rhs.stack_) {
262      stack_ = new int [numberRows_+1];
263      memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
264    } else {
265      stack_ = NULL;
266    }
267    if (rhs.toLeaf_) {
268      toLeaf_ = new int [numberRows_+1];
269      memcpy(toLeaf_,rhs.toLeaf_,(numberRows_+1)*sizeof(int));
270    } else {
271      toLeaf_ = NULL;
272    }
273    if (rhs.toRoot_) {
274      toRoot_ = new int [numberRows_+1];
275      memcpy(toRoot_,rhs.toRoot_,(numberRows_+1)*sizeof(int));
276    } else {
277      toRoot_ = NULL;
278    }
279    if (rhs.mark_) {
280      mark_ = new char [numberRows_+1];
281      memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
282    } else {
283    mark_ = NULL;
284    }
285    model_=rhs.model_;
286  }
287  return *this;
288}
289/* Replaces one Column to basis,
290   returns 0=OK
291*/
292int 
293ClpNetworkBasis::replaceColumn ( CoinIndexedVector * regionSparse,
294                                 int pivotRow)
295{
296  // regionSparse is empty
297  assert (!regionSparse->getNumElements());
298  model_->unpack(regionSparse, model_->sequenceIn());
299  // arc given by pivotRow is leaving basis
300  int iParent = parent_[pivotRow];
301  // arc coming in has these two nodes
302  int * indices = regionSparse->getIndices();
303  int iRow0 = indices[0];
304  int iRow1;
305  if (regionSparse->getNumElements()==2)
306    iRow1 = indices[1];
307  else
308    iRow1 = numberRows_;
309  double sign = -regionSparse->denseVector()[iRow0];
310  printf("In %d (%g) %d pivoting on %d\n",
311         iRow1, sign, iRow0,pivotRow);
312  regionSparse->clear();
313  // take out of tree
314  int iLeft = leftSibling_[pivotRow];
315  int iRight = rightSibling_[pivotRow];
316  if (iLeft>=0) {
317    rightSibling_[iLeft] = iRight;
318    if (iRight>=0) 
319      leftSibling_[iRight]=iLeft;
320  } else if (iRight>=0) {
321    leftSibling_[iRight]=-1;
322    descendant_[iParent]=iRight;
323  } else {
324    descendant_[iParent]=-1;;
325  }
326  // move to other end of chain
327  int descendant = descendant_[pivotRow];
328  if (descendant>=0) {
329    // make this descendant of that
330    if (descendant_[descendant]>=0) {
331      // we have a sibling
332      int iRight = descendant_[descendant];
333      rightSibling_[pivotRow]=iRight;
334      leftSibling_[iRight]=pivotRow;
335    } else {
336      rightSibling_[pivotRow]=-1;
337    }       
338    descendant_[descendant] = pivotRow;
339    leftSibling_[pivotRow]=-1;
340  }
341  // now insert new one
342  descendant = descendant_[iRow1];
343  parent_[iRow0] = iRow1;
344  sign_[iRow1]= sign;
345  if (descendant>=0) {
346    // we have a sibling
347    int iRight = descendant;
348    rightSibling_[iRow1]=iRight;
349    leftSibling_[iRight]=iRow1;
350  } else {
351    rightSibling_[iRow1]=-1;
352  }         
353  descendant_[descendant] = iRow1;
354  leftSibling_[iRow1]=-1;
355  return 0;
356}
357
358/* Updates one column (FTRAN) from region2 */
359int 
360ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse2)
361{
362  int iPivot=leaf_;
363  int numberNonZero=0;
364  double * array = regionSparse2->denseVector();
365  while (iPivot!=numberRows_) {
366    double pivotValue = array[iPivot];
367    if (pivotValue) {
368      numberNonZero++;
369      int otherRow = parent_[iPivot];
370      if (sign_[iPivot]<0) {
371        array[iPivot] = - pivotValue;
372        array[otherRow] += pivotValue;
373      } else {
374        array[otherRow] += pivotValue;
375      }
376    }
377    iPivot = toRoot_[iPivot];
378  }
379  array[numberRows_]=0.0;
380  return numberNonZero;
381}
382/* Updates one column (FTRAN) to/from array
383    ** For large problems you should ALWAYS know where the nonzeros
384   are, so please try and migrate to previous method after you
385   have got code working using this simple method - thank you!
386   (the only exception is if you know input is dense e.g. rhs) */
387int 
388ClpNetworkBasis::updateColumn ( double array[] ) const
389{
390  int iPivot=leaf_;
391  int numberNonZero=0;
392  while (iPivot!=numberRows_) {
393    double pivotValue = array[iPivot];
394    if (pivotValue) {
395      numberNonZero++;
396      int otherRow = parent_[iPivot];
397      if (sign_[iPivot]<0) {
398        array[iPivot] = - pivotValue;
399        array[otherRow] += pivotValue;
400      } else {
401        array[otherRow] += pivotValue;
402      }
403    }
404    iPivot = toRoot_[iPivot];
405  }
406  array[numberRows_]=0.0;
407  return numberNonZero;
408}
409/* Updates one column transpose (BTRAN)
410   For large problems you should ALWAYS know where the nonzeros
411   are, so please try and migrate to previous method after you
412   have got code working using this simple method - thank you!
413   (the only exception is if you know input is dense e.g. dense objective)
414   returns number of nonzeros */
415int 
416ClpNetworkBasis::updateColumnTranspose ( double array[] ) const
417{
418  abort();
419  return 1;
420}
421/* Updates one column (BTRAN) from region2 */
422int 
423ClpNetworkBasis::updateColumnTranspose ( 
424                                        CoinIndexedVector * regionSparse2) const
425{
426  abort();
427  return 1;
428}
Note: See TracBrowser for help on using the repository browser.