source: trunk/Clp/src/ClpNetworkBasis.cpp @ 1402

Last change on this file since 1402 was 1402, checked in by forrest, 10 years ago

get rid of compiler warnings

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 30.9 KB
RevLine 
[1370]1/* $Id: ClpNetworkBasis.cpp 1402 2009-07-25 08:39:55Z forrest $ */
[124]2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6#include "ClpNetworkBasis.hpp"
7#include "CoinHelperFunctions.hpp"
8#include "ClpSimplex.hpp"
9#include "ClpMatrixBase.hpp"
10#include "CoinIndexedVector.hpp"
11
12
13//#############################################################################
14// Constructors / Destructor / Assignment
15//#############################################################################
16
17//-------------------------------------------------------------------
18// Default Constructor
19//-------------------------------------------------------------------
20ClpNetworkBasis::ClpNetworkBasis () 
21{
[1197]22#ifndef COIN_FAST_CODE
[124]23  slackValue_=-1.0;
[1197]24#endif
[124]25  numberRows_=0;
26  numberColumns_=0;
27  parent_ = NULL;
28  descendant_ = NULL;
29  pivot_ = NULL;
30  rightSibling_ = NULL;
31  leftSibling_ = NULL;
32  sign_ = NULL;
33  stack_ = NULL;
[130]34  permute_ = NULL;
35  permuteBack_ = NULL;
36  stack2_=NULL;
37  depth_ = NULL;
[124]38  mark_ = NULL;
39  model_=NULL;
40}
41// Constructor from CoinFactorization
42ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex * model,
[1302]43                                 int numberRows, const CoinFactorizationDouble * pivotRegion,
[124]44                                 const int * permuteBack,
[277]45                                 const CoinBigIndex * startColumn, 
[124]46                                 const int * numberInColumn,
[1402]47                                 const int * indexRow, const CoinFactorizationDouble * /*element*/)
[124]48{
[1197]49#ifndef COIN_FAST_CODE
[124]50  slackValue_=-1.0;
[1197]51#endif
[124]52  numberRows_=numberRows;
53  numberColumns_=numberRows;
54  parent_ = new int [ numberRows_+1];
55  descendant_ = new int [ numberRows_+1];
56  pivot_ = new int [ numberRows_+1];
57  rightSibling_ = new int [ numberRows_+1];
58  leftSibling_ = new int [ numberRows_+1];
59  sign_ = new double [ numberRows_+1];
60  stack_ = new int [ numberRows_+1];
[130]61  stack2_ = new int[numberRows_+1];
62  depth_ = new int[numberRows_+1];
[124]63  mark_ = new char[numberRows_+1];
[130]64  permute_ = new int [numberRows_ + 1];
65  permuteBack_ = new int [numberRows_ + 1];
[124]66  int i;
67  for (i=0;i<numberRows_+1;i++) {
68    parent_[i]=-1;
69    descendant_[i]=-1;
70    pivot_[i]=-1;
71    rightSibling_[i]=-1;
72    leftSibling_[i]=-1;
73    sign_[i]=-1.0;
74    stack_[i]=-1;
[130]75    permute_[i]=i;
76    permuteBack_[i]=i;
77    stack2_[i]=-1;
78    depth_[i]=-1;
[124]79    mark_[i]=0;
80  }
[130]81  mark_[numberRows_]=1;
[124]82  // pivotColumnBack gives order of pivoting into basis
83  // so pivotColumnback[0] is first slack in basis and
84  // it pivots on row permuteBack[0]
85  // a known root is given by permuteBack[numberRows_-1]
86  int lastPivot=numberRows_;
87  for (i=0;i<numberRows_;i++) {
88    int iPivot = permuteBack[i];
89    lastPivot=iPivot;
90    double sign;
91    if (pivotRegion[i]>0.0)
92      sign = 1.0;
93    else
94      sign =-1.0;
95    int other;
96    if (numberInColumn[i]>0) {
97      int iRow = indexRow[startColumn[i]];
98      other = permuteBack[iRow];
[131]99      //assert (parent_[other]!=-1);
[124]100    } else {
101      other = numberRows_;
102    }
103    sign_[iPivot] = sign;
104    int iParent = other;
105    parent_[iPivot] = other;
106    if (descendant_[iParent]>=0) {
107      // we have a sibling
108      int iRight = descendant_[iParent];
109      rightSibling_[iPivot]=iRight;
110      leftSibling_[iRight]=iPivot;
111    } else {
112      rightSibling_[iPivot]=-1;
113    }       
114    descendant_[iParent] = iPivot;
115    leftSibling_[iPivot]=-1;
116  }
[130]117  // do depth
118  int iPivot = numberRows_;
119  int nStack = 1;
120  stack_[0]=descendant_[numberRows_];
121  depth_[numberRows_]=-1; // root
122  while (nStack) {
123    // take off
124    int iNext = stack_[--nStack];
125    if (iNext>=0) {
126      depth_[iNext] = nStack;
127      iPivot = iNext;
128      int iRight = rightSibling_[iNext];
129      stack_[nStack++] = iRight;
130      if (descendant_[iNext]>=0)
131        stack_[nStack++] = descendant_[iNext];
132    }
133  }
[124]134  model_=model;
[130]135  check();
[124]136}
137
138//-------------------------------------------------------------------
139// Copy constructor
140//-------------------------------------------------------------------
141ClpNetworkBasis::ClpNetworkBasis (const ClpNetworkBasis & rhs) 
142{
[1197]143#ifndef COIN_FAST_CODE
[124]144  slackValue_=rhs.slackValue_;
[1197]145#endif
[124]146  numberRows_=rhs.numberRows_;
147  numberColumns_=rhs.numberColumns_;
148  if (rhs.parent_) {
149    parent_ = new int [numberRows_+1];
[1197]150    CoinMemcpyN(rhs.parent_,(numberRows_+1),parent_);
[124]151  } else {
152    parent_ = NULL;
153  }
154  if (rhs.descendant_) {
155    descendant_ = new int [numberRows_+1];
[1197]156    CoinMemcpyN(rhs.descendant_,(numberRows_+1),descendant_);
[124]157  } else {
158    descendant_ = NULL;
159  }
160  if (rhs.pivot_) {
161    pivot_ = new int [numberRows_+1];
[1197]162    CoinMemcpyN(rhs.pivot_,(numberRows_+1),pivot_);
[124]163  } else {
164    pivot_ = NULL;
165  }
166  if (rhs.rightSibling_) {
167    rightSibling_ = new int [numberRows_+1];
[1197]168    CoinMemcpyN(rhs.rightSibling_,(numberRows_+1),rightSibling_);
[124]169  } else {
170    rightSibling_ = NULL;
171  }
172  if (rhs.leftSibling_) {
173    leftSibling_ = new int [numberRows_+1];
[1197]174    CoinMemcpyN(rhs.leftSibling_,(numberRows_+1),leftSibling_);
[124]175  } else {
176    leftSibling_ = NULL;
177  }
178  if (rhs.sign_) {
179    sign_ = new double [numberRows_+1];
[1197]180    CoinMemcpyN(rhs.sign_,(numberRows_+1),sign_);
[124]181  } else {
182    sign_ = NULL;
183  }
184  if (rhs.stack_) {
185    stack_ = new int [numberRows_+1];
[1197]186    CoinMemcpyN(rhs.stack_,(numberRows_+1),stack_);
[124]187  } else {
188    stack_ = NULL;
189  }
[130]190  if (rhs.permute_) {
191    permute_ = new int [numberRows_+1];
[1197]192    CoinMemcpyN(rhs.permute_,(numberRows_+1),permute_);
[124]193  } else {
[130]194    permute_ = NULL;
[124]195  }
[130]196  if (rhs.permuteBack_) {
197    permuteBack_ = new int [numberRows_+1];
[1197]198    CoinMemcpyN(rhs.permuteBack_,(numberRows_+1),permuteBack_);
[124]199  } else {
[130]200    permuteBack_ = NULL;
[124]201  }
[130]202  if (rhs.stack2_) {
203    stack2_ = new int [numberRows_+1];
[1197]204    CoinMemcpyN(rhs.stack2_,(numberRows_+1),stack2_);
[130]205  } else {
206    stack2_ = NULL;
207  }
208  if (rhs.depth_) {
209    depth_ = new int [numberRows_+1];
[1197]210    CoinMemcpyN(rhs.depth_,(numberRows_+1),depth_);
[130]211  } else {
212    depth_ = NULL;
213  }
[124]214  if (rhs.mark_) {
215    mark_ = new char [numberRows_+1];
[1197]216    CoinMemcpyN(rhs.mark_,(numberRows_+1),mark_);
[124]217  } else {
218    mark_ = NULL;
219  }
220  model_=rhs.model_;
221}
222
223//-------------------------------------------------------------------
224// Destructor
225//-------------------------------------------------------------------
226ClpNetworkBasis::~ClpNetworkBasis () 
227{
228  delete [] parent_;
229  delete [] descendant_;
230  delete [] pivot_;
231  delete [] rightSibling_;
232  delete [] leftSibling_;
233  delete [] sign_;
234  delete [] stack_;
[130]235  delete [] permute_;
236  delete [] permuteBack_;
237  delete [] stack2_;
238  delete [] depth_;
[124]239  delete [] mark_;
240}
241
242//----------------------------------------------------------------
243// Assignment operator
244//-------------------------------------------------------------------
245ClpNetworkBasis &
246ClpNetworkBasis::operator=(const ClpNetworkBasis& rhs)
247{
248  if (this != &rhs) {
249    delete [] parent_;
250    delete [] descendant_;
251    delete [] pivot_;
252    delete [] rightSibling_;
253    delete [] leftSibling_;
254    delete [] sign_;
255    delete [] stack_;
[130]256    delete [] permute_;
257    delete [] permuteBack_;
258    delete [] stack2_;
259    delete [] depth_;
[124]260    delete [] mark_;
[1197]261#ifndef COIN_FAST_CODE
[124]262    slackValue_=rhs.slackValue_;
[1197]263#endif
[124]264    numberRows_=rhs.numberRows_;
265    numberColumns_=rhs.numberColumns_;
266    if (rhs.parent_) {
267      parent_ = new int [numberRows_+1];
[1197]268      CoinMemcpyN(rhs.parent_,(numberRows_+1),parent_);
[124]269    } else {
270      parent_ = NULL;
271    }
272    if (rhs.descendant_) {
273      descendant_ = new int [numberRows_+1];
[1197]274      CoinMemcpyN(rhs.descendant_,(numberRows_+1),descendant_);
[124]275    } else {
276      descendant_ = NULL;
277    }
278    if (rhs.pivot_) {
279      pivot_ = new int [numberRows_+1];
[1197]280      CoinMemcpyN(rhs.pivot_,(numberRows_+1),pivot_);
[124]281    } else {
282      pivot_ = NULL;
283    }
284    if (rhs.rightSibling_) {
285      rightSibling_ = new int [numberRows_+1];
[1197]286      CoinMemcpyN(rhs.rightSibling_,(numberRows_+1),rightSibling_);
[124]287    } else {
288      rightSibling_ = NULL;
289    }
290    if (rhs.leftSibling_) {
291      leftSibling_ = new int [numberRows_+1];
[1197]292      CoinMemcpyN(rhs.leftSibling_,(numberRows_+1),leftSibling_);
[124]293    } else {
294      leftSibling_ = NULL;
295    }
296    if (rhs.sign_) {
297      sign_ = new double [numberRows_+1];
[1197]298      CoinMemcpyN(rhs.sign_,(numberRows_+1),sign_);
[124]299    } else {
300      sign_ = NULL;
301    }
302    if (rhs.stack_) {
303      stack_ = new int [numberRows_+1];
[1197]304      CoinMemcpyN(rhs.stack_,(numberRows_+1),stack_);
[124]305    } else {
306      stack_ = NULL;
307    }
[130]308    if (rhs.permute_) {
309      permute_ = new int [numberRows_+1];
[1197]310      CoinMemcpyN(rhs.permute_,(numberRows_+1),permute_);
[124]311    } else {
[130]312      permute_ = NULL;
[124]313    }
[130]314    if (rhs.permuteBack_) {
315      permuteBack_ = new int [numberRows_+1];
[1197]316      CoinMemcpyN(rhs.permuteBack_,(numberRows_+1),permuteBack_);
[124]317    } else {
[130]318      permuteBack_ = NULL;
[124]319    }
[130]320    if (rhs.stack2_) {
321      stack2_ = new int [numberRows_+1];
[1197]322      CoinMemcpyN(rhs.stack2_,(numberRows_+1),stack2_);
[130]323    } else {
324      stack2_ = NULL;
325    }
326    if (rhs.depth_) {
327      depth_ = new int [numberRows_+1];
[1197]328      CoinMemcpyN(rhs.depth_,(numberRows_+1),depth_);
[130]329    } else {
330      depth_ = NULL;
331    }
[124]332    if (rhs.mark_) {
333      mark_ = new char [numberRows_+1];
[1197]334      CoinMemcpyN(rhs.mark_,(numberRows_+1),mark_);
[124]335    } else {
[130]336      mark_ = NULL;
[124]337    }
338  }
339  return *this;
340}
[130]341// checks looks okay
342void ClpNetworkBasis::check()
343{
344  // check depth
345  int iPivot = numberRows_;
346  int nStack = 1;
347  stack_[0]=descendant_[numberRows_];
348  depth_[numberRows_]=-1; // root
349  while (nStack) {
350    // take off
351    int iNext = stack_[--nStack];
352    if (iNext>=0) {
[131]353      //assert (depth_[iNext]==nStack);
[130]354      depth_[iNext] = nStack;
355      iPivot = iNext;
356      int iRight = rightSibling_[iNext];
357      stack_[nStack++] = iRight;
358      if (descendant_[iNext]>=0)
359        stack_[nStack++] = descendant_[iNext];
360    }
361  }
362}
363// prints
364void ClpNetworkBasis::print()
365{
366  int i;
367  printf("       parent descendant     left    right   sign    depth\n");
368  for (i=0;i<numberRows_+1;i++)
369    printf("%4d  %7d   %8d  %7d  %7d  %5g  %7d\n",
370           i,parent_[i],descendant_[i],leftSibling_[i],rightSibling_[i],
371           sign_[i],depth_[i]);
372}
[124]373/* Replaces one Column to basis,
374   returns 0=OK
375*/
376int 
377ClpNetworkBasis::replaceColumn ( CoinIndexedVector * regionSparse,
378                                 int pivotRow)
379{
[130]380  // When things have settled down then redo this to make more elegant
381  // I am sure lots of loops can be combined
[124]382  // regionSparse is empty
383  assert (!regionSparse->getNumElements());
384  model_->unpack(regionSparse, model_->sequenceIn());
385  // arc given by pivotRow is leaving basis
[130]386  //int kParent = parent_[pivotRow];
[124]387  // arc coming in has these two nodes
388  int * indices = regionSparse->getIndices();
389  int iRow0 = indices[0];
390  int iRow1;
391  if (regionSparse->getNumElements()==2)
392    iRow1 = indices[1];
393  else
394    iRow1 = numberRows_;
395  double sign = -regionSparse->denseVector()[iRow0];
396  regionSparse->clear();
[130]397  // and outgoing
398  model_->unpack(regionSparse,model_->pivotVariable()[pivotRow]);
399  int jRow0 = indices[0];
400  int jRow1;
401  if (regionSparse->getNumElements()==2)
402    jRow1 = indices[1];
403  else
404    jRow1 = numberRows_;
405  regionSparse->clear();
406  // get correct pivotRow
407  //#define FULL_DEBUG
408#ifdef FULL_DEBUG
409  printf ("irow %d %d, jrow %d %d\n",
410          iRow0,iRow1,jRow0,jRow1);
411#endif
412  if (parent_[jRow0]==jRow1) {
413    int newPivot = jRow0;
414    if (newPivot!=pivotRow) {
415#ifdef FULL_DEBUG
416      printf("pivot row of %d permuted to %d\n",pivotRow,newPivot);
417#endif
418      pivotRow = newPivot;
419    }
[124]420  } else {
[131]421    //assert (parent_[jRow1]==jRow0);
[130]422    int newPivot = jRow1;
423    if (newPivot!=pivotRow) {
424#ifdef FULL_DEBUG
425      printf("pivot row of %d permuted to %d\n",pivotRow,newPivot);
426#endif
427      pivotRow = newPivot;
428    }
[124]429  }
[130]430  bool extraPrint = (model_->numberIterations()>-3)&&
431    (model_->logLevel()>10);
432  if (extraPrint)
433    print();
434#ifdef FULL_DEBUG
435  printf("In %d (region= %g, stored %g) %d (%g) pivoting on %d (%g)\n",
436         iRow1, sign, sign_[iRow1], iRow0, sign_[iRow0] ,pivotRow,sign_[pivotRow]);
437#endif
438  // see which path outgoing pivot is on
439  int kRow = -1;
440  int jRow = iRow1;
441  while (jRow!=numberRows_) {
442    if (jRow==pivotRow) {
443      kRow = iRow1;
444      break;
445    } else {
446      jRow = parent_[jRow];
447    }
448  }
449  if (kRow<0) {
450    jRow = iRow0;
451    while (jRow!=numberRows_) {
452      if (jRow==pivotRow) {
453        kRow = iRow0;
454        break;
455      } else {
456        jRow = parent_[jRow];
457      }
458    }
459  }
[131]460  //assert (kRow>=0);
[130]461  if (iRow0==kRow) {
462    iRow0 = iRow1;
463    iRow1 = kRow;
464    sign = -sign;
465  }
466  // pivot row is on path from iRow1 back to root
467  // get stack of nodes to change
468  // Also get precursors for cleaning order
469  int nStack=1;
470  stack_[0]=iRow0;
471  while (kRow!=pivotRow) {
472    stack_[nStack++] = kRow;
473    if (sign*sign_[kRow]<0.0) {
474      sign_[kRow]= -sign_[kRow];
475    } else {
476      sign = -sign;
477    }
478    kRow = parent_[kRow];
479    //sign *= sign_[kRow];
480  }
481  stack_[nStack++]=pivotRow;
482  if (sign*sign_[pivotRow]<0.0) {
483    sign_[pivotRow]= -sign_[pivotRow];
484  } else {
485    sign = -sign;
486  }
487  int iParent = parent_[pivotRow];
488  while (nStack>1) {
489    int iLeft;
490    int iRight;
491    kRow = stack_[--nStack];
492    int newParent = stack_[nStack-1];
493#ifdef FULL_DEBUG
494    printf("row %d, old parent %d, new parent %d, pivotrow %d\n",kRow,
495           iParent,newParent,pivotRow);
496#endif
497    int i1 = permuteBack_[pivotRow];
498    int i2 = permuteBack_[kRow];
499    permuteBack_[pivotRow]=i2;
500    permuteBack_[kRow]=i1;
501    // do Btran permutation
502    permute_[i1]=kRow;
503    permute_[i2]=pivotRow;
504    pivotRow = kRow;
505    // Take out of old parent
506    iLeft = leftSibling_[kRow];
507    iRight = rightSibling_[kRow];
508    if (iLeft<0) {
509      // take out of tree
510      if (iRight>=0) {
511        leftSibling_[iRight]=iLeft;
512        descendant_[iParent]=iRight;
513      } else {
514#ifdef FULL_DEBUG
515        printf("Saying %d (old parent of %d) has no descendants\n",
516               iParent, kRow);
517#endif
518        descendant_[iParent]=-1;
519      }
520    } else {
521      // take out of tree
522      rightSibling_[iLeft] = iRight;
523      if (iRight>=0) 
524        leftSibling_[iRight]=iLeft;
525    }
526    leftSibling_[kRow]=-1;
527    rightSibling_[kRow]=-1;
528   
529    // now insert new one
[124]530    // make this descendant of that
[130]531    if (descendant_[newParent]>=0) {
532      // we will have a sibling
533      int jRight = descendant_[newParent];
534      rightSibling_[kRow]=jRight;
535      leftSibling_[jRight]=kRow;
[124]536    } else {
[130]537      rightSibling_[kRow]=-1;
[124]538    }       
[130]539    descendant_[newParent] = kRow;
540    leftSibling_[kRow]=-1;
541    parent_[kRow]=newParent;
542     
543    iParent = kRow;
[124]544  }
[130]545  // now redo all depths from stack_[1]
546  // This must be possible to combine - but later
547  {
548    int iPivot  = stack_[1];
549    int iDepth=depth_[parent_[iPivot]]; //depth of parent
550    iDepth ++; //correct for this one
551    int nStack = 1;
552    stack_[0]=iPivot;
553    while (nStack) {
554      // take off
555      int iNext = stack_[--nStack];
556      if (iNext>=0) {
557        // add stack level
558        depth_[iNext]=nStack+iDepth;
559        stack_[nStack++] = rightSibling_[iNext];
560        if (descendant_[iNext]>=0)
561          stack_[nStack++] = descendant_[iNext];
562      }
563    }
564  }
565  if (extraPrint)
566    print();
567  //check();
[124]568  return 0;
569}
570
571/* Updates one column (FTRAN) from region2 */
[225]572double
[130]573ClpNetworkBasis::updateColumn (  CoinIndexedVector * regionSparse,
[225]574                                 CoinIndexedVector * regionSparse2,
575                                 int pivotRow)
[124]576{
[130]577  regionSparse->clear (  );
578  double *region = regionSparse->denseVector (  );
579  double *region2 = regionSparse2->denseVector (  );
580  int *regionIndex2 = regionSparse2->getIndices (  );
581  int numberNonZero = regionSparse2->getNumElements (  );
582  int *regionIndex = regionSparse->getIndices (  );
583  int i;
584  bool doTwo = (numberNonZero==2);
[225]585  int i0=-1;
586  int i1=-1;
[130]587  if (doTwo) {
588    i0 = regionIndex2[0];
589    i1 = regionIndex2[1];
590  }
[225]591  double returnValue=0.0;
592  bool packed = regionSparse2->packedMode();
593  if (packed) {
594    if (doTwo&&region2[0]*region2[1]<0.0) {
595      region[i0] = region2[0];
596      region2[0]=0.0;
597      region[i1] = region2[1];
598      region2[1]=0.0;
599      int iDepth0 = depth_[i0];
600      int iDepth1 = depth_[i1];
601      if (iDepth1>iDepth0) {
602        int temp = i0;
603        i0 = i1;
604        i1 = temp;
605        temp = iDepth0;
606        iDepth0 = iDepth1;
607        iDepth1 = temp;
[124]608      }
[225]609      numberNonZero=0;
610      if (pivotRow<0) {
611        while (iDepth0>iDepth1) {
612          double pivotValue = region[i0];
[130]613          // put back now ?
[225]614          int iBack = permuteBack_[i0];
615          region2[numberNonZero] = pivotValue*sign_[i0];
[130]616          regionIndex2[numberNonZero++] = iBack;
[225]617          int otherRow = parent_[i0];
618          region[i0] =0.0;
[130]619          region[otherRow] += pivotValue;
[225]620          iDepth0--;
621          i0 = otherRow;
[130]622        }
[225]623        while (i0!=i1) {
624          double pivotValue = region[i0];
625          // put back now ?
626          int iBack = permuteBack_[i0];
627          region2[numberNonZero] = pivotValue*sign_[i0];
628          regionIndex2[numberNonZero++] = iBack;
629          int otherRow = parent_[i0];
630          region[i0] =0.0;
631          region[otherRow] += pivotValue;
632          i0 = otherRow;
633          double pivotValue1 = region[i1];
634          // put back now ?
635          int iBack1 = permuteBack_[i1];
636          region2[numberNonZero] = pivotValue1*sign_[i1];
637          regionIndex2[numberNonZero++] = iBack1;
638          int otherRow1 = parent_[i1];
639          region[i1] =0.0;
640          region[otherRow1] += pivotValue1;
641          i1 = otherRow1;
642        }
643      } else {
644        while (iDepth0>iDepth1) {
645          double pivotValue = region[i0];
646          // put back now ?
647          int iBack = permuteBack_[i0];
648          double value = pivotValue*sign_[i0];
649          region2[numberNonZero] = value;
650          regionIndex2[numberNonZero++] = iBack;
651          if (iBack==pivotRow)
652            returnValue = value;
653          int otherRow = parent_[i0];
654          region[i0] =0.0;
655          region[otherRow] += pivotValue;
656          iDepth0--;
657          i0 = otherRow;
658        }
659        while (i0!=i1) {
660          double pivotValue = region[i0];
661          // put back now ?
662          int iBack = permuteBack_[i0];
663          double value = pivotValue*sign_[i0];
664          region2[numberNonZero] = value;
665          regionIndex2[numberNonZero++] = iBack;
666          if (iBack==pivotRow)
667            returnValue = value;
668          int otherRow = parent_[i0];
669          region[i0] =0.0;
670          region[otherRow] += pivotValue;
671          i0 = otherRow;
672          double pivotValue1 = region[i1];
673          // put back now ?
674          int iBack1 = permuteBack_[i1];
675          value = pivotValue1*sign_[i1];
676          region2[numberNonZero] = value;
677          regionIndex2[numberNonZero++] = iBack1;
678          if (iBack1==pivotRow)
679            returnValue = value;
680          int otherRow1 = parent_[i1];
681          region[i1] =0.0;
682          region[otherRow1] += pivotValue1;
683          i1 = otherRow1;
684        }
[130]685      }
[225]686    } else {
687      // set up linked lists at each depth
688      // stack2 is start, stack is next
689      int greatestDepth=-1;
690      //mark_[numberRows_]=1;
691      for (i=0;i<numberNonZero;i++) {
692        int j = regionIndex2[i];
693        double value = region2[i];
694        region2[i]=0.0;
695        region[j]=value;
696        regionIndex[i]=j;
697        int iDepth = depth_[j];
698        if (iDepth>greatestDepth) 
699          greatestDepth = iDepth;
700        // and back until marked
701        while (!mark_[j]) {
702          int iNext = stack2_[iDepth];
703          stack2_[iDepth]=j;
704          stack_[j]=iNext;
705          mark_[j]=1;
706          iDepth--;
707          j=parent_[j];
708        }
709      }
710      numberNonZero=0;
711      if (pivotRow<0) {
712        for (;greatestDepth>=0; greatestDepth--) {
713          int iPivot = stack2_[greatestDepth];
714          stack2_[greatestDepth]=-1;
715          while (iPivot>=0) {
716            mark_[iPivot]=0;
717            double pivotValue = region[iPivot];
718            if (pivotValue) {
719              // put back now ?
720              int iBack = permuteBack_[iPivot];
721              region2[numberNonZero] = pivotValue*sign_[iPivot];
722              regionIndex2[numberNonZero++] = iBack;
723              int otherRow = parent_[iPivot];
724              region[iPivot] =0.0;
725            region[otherRow] += pivotValue;
726            }
727            iPivot = stack_[iPivot];
728          }
729        }
730      } else {
731        for (;greatestDepth>=0; greatestDepth--) {
732          int iPivot = stack2_[greatestDepth];
733          stack2_[greatestDepth]=-1;
734          while (iPivot>=0) {
735            mark_[iPivot]=0;
736            double pivotValue = region[iPivot];
737            if (pivotValue) {
738              // put back now ?
739              int iBack = permuteBack_[iPivot];
740              double value = pivotValue*sign_[iPivot];
741              region2[numberNonZero] = value;
742              regionIndex2[numberNonZero++] = iBack;
743              if (iBack==pivotRow)
744                returnValue = value;
745              int otherRow = parent_[iPivot];
746              region[iPivot] =0.0;
747              region[otherRow] += pivotValue;
748            }
749            iPivot = stack_[iPivot];
750          }
751        }
752      }
[130]753    }
[225]754  } else {
755    if (doTwo&&region2[i0]*region2[i1]<0.0) {
756      // If just +- 1 then could go backwards on depth until join
757      region[i0] = region2[i0];
758      region2[i0]=0.0;
759      region[i1] = region2[i1];
760      region2[i1]=0.0;
761      int iDepth0 = depth_[i0];
762      int iDepth1 = depth_[i1];
763      if (iDepth1>iDepth0) {
764        int temp = i0;
765        i0 = i1;
766        i1 = temp;
767        temp = iDepth0;
768        iDepth0 = iDepth1;
769        iDepth1 = temp;
770      }
771      numberNonZero=0;
772      while (iDepth0>iDepth1) {
773        double pivotValue = region[i0];
774        // put back now ?
775        int iBack = permuteBack_[i0];
776        regionIndex2[numberNonZero++] = iBack;
777        int otherRow = parent_[i0];
778        region2[iBack] = pivotValue*sign_[i0];
779        region[i0] =0.0;
780        region[otherRow] += pivotValue;
781        iDepth0--;
782        i0 = otherRow;
783      }
784      while (i0!=i1) {
785        double pivotValue = region[i0];
786        // put back now ?
787        int iBack = permuteBack_[i0];
788        regionIndex2[numberNonZero++] = iBack;
789        int otherRow = parent_[i0];
790        region2[iBack] = pivotValue*sign_[i0];
791        region[i0] =0.0;
792        region[otherRow] += pivotValue;
793        i0 = otherRow;
794        double pivotValue1 = region[i1];
795        // put back now ?
796        int iBack1 = permuteBack_[i1];
797        regionIndex2[numberNonZero++] = iBack1;
798        int otherRow1 = parent_[i1];
799        region2[iBack1] = pivotValue1*sign_[i1];
800        region[i1] =0.0;
801        region[otherRow1] += pivotValue1;
802        i1 = otherRow1;
803      }
804    } else {
805      // set up linked lists at each depth
806      // stack2 is start, stack is next
807      int greatestDepth=-1;
808      //mark_[numberRows_]=1;
809      for (i=0;i<numberNonZero;i++) {
810        int j = regionIndex2[i];
811        double value = region2[j];
812        region2[j]=0.0;
813        region[j]=value;
814        regionIndex[i]=j;
815        int iDepth = depth_[j];
816        if (iDepth>greatestDepth) 
817          greatestDepth = iDepth;
818        // and back until marked
819        while (!mark_[j]) {
820          int iNext = stack2_[iDepth];
821          stack2_[iDepth]=j;
822          stack_[j]=iNext;
823          mark_[j]=1;
824          iDepth--;
825          j=parent_[j];
826        }
827      }
828      numberNonZero=0;
829      for (;greatestDepth>=0; greatestDepth--) {
830        int iPivot = stack2_[greatestDepth];
831        stack2_[greatestDepth]=-1;
832        while (iPivot>=0) {
833          mark_[iPivot]=0;
834          double pivotValue = region[iPivot];
835          if (pivotValue) {
836            // put back now ?
837            int iBack = permuteBack_[iPivot];
838            regionIndex2[numberNonZero++] = iBack;
839            int otherRow = parent_[iPivot];
840            region2[iBack] = pivotValue*sign_[iPivot];
841            region[iPivot] =0.0;
842            region[otherRow] += pivotValue;
843          }
844          iPivot = stack_[iPivot];
845        }
846      }
847    }
848    if (pivotRow>=0)
849      returnValue = region2[pivotRow];
[124]850  }
[130]851  region[numberRows_]=0.0;
852  regionSparse2->setNumElements(numberNonZero);
853#ifdef FULL_DEBUG
854 {
855   int i;
856   for (i=0;i<numberRows_;i++) {
857     assert(!mark_[i]);
858     assert (stack2_[i]==-1);
859   }
860 }
861#endif
[225]862  return returnValue;
[124]863}
864/* Updates one column (FTRAN) to/from array
865    ** For large problems you should ALWAYS know where the nonzeros
866   are, so please try and migrate to previous method after you
867   have got code working using this simple method - thank you!
868   (the only exception is if you know input is dense e.g. rhs) */
869int 
[130]870ClpNetworkBasis::updateColumn (  CoinIndexedVector * regionSparse,
871                                 double region2[] ) const
[124]872{
[130]873  regionSparse->clear (  );
874  double *region = regionSparse->denseVector (  );
875  int numberNonZero = 0;
876  int *regionIndex = regionSparse->getIndices (  );
877  int i;
878  // set up linked lists at each depth
879  // stack2 is start, stack is next
880  int greatestDepth=-1;
881  for (i=0;i<numberRows_;i++) {
882    double value = region2[i];
883    if (value) {
884      region2[i]=0.0;
885      region[i]=value;
886      regionIndex[numberNonZero++]=i;
887      int j=i;
888      int iDepth = depth_[j];
889      if (iDepth>greatestDepth) 
890        greatestDepth = iDepth;
891      // and back until marked
892      while (!mark_[j]) {
893        int iNext = stack2_[iDepth];
894        stack2_[iDepth]=j;
895        stack_[j]=iNext;
896        mark_[j]=1;
897        iDepth--;
898        j=parent_[j];
[124]899      }
900    }
901  }
[130]902  numberNonZero=0;
903  for (;greatestDepth>=0; greatestDepth--) {
904    int iPivot = stack2_[greatestDepth];
905    stack2_[greatestDepth]=-1;
906    while (iPivot>=0) {
907      mark_[iPivot]=0;
908      double pivotValue = region[iPivot];
909      if (pivotValue) {
910        // put back now ?
911        int iBack = permuteBack_[iPivot];
912        numberNonZero++;
913        int otherRow = parent_[iPivot];
914        region2[iBack] = pivotValue*sign_[iPivot];
915        region[iPivot] =0.0;
916        region[otherRow] += pivotValue;
917      }
918      iPivot = stack_[iPivot];
919    }
920  }
921  region[numberRows_]=0.0;
[124]922  return numberNonZero;
923}
924/* Updates one column transpose (BTRAN)
925   For large problems you should ALWAYS know where the nonzeros
926   are, so please try and migrate to previous method after you
927   have got code working using this simple method - thank you!
928   (the only exception is if you know input is dense e.g. dense objective)
929   returns number of nonzeros */
930int 
[130]931ClpNetworkBasis::updateColumnTranspose (  CoinIndexedVector * regionSparse,
932                                          double region2[] ) const
[124]933{
[130]934  // permute in after copying
935  // so will end up in right place
936  double *region = regionSparse->denseVector (  );
937  int *regionIndex = regionSparse->getIndices (  );
938  int i;
939  int numberNonZero=0;
[1197]940  CoinMemcpyN(region2,numberRows_,region);
[130]941  for (i=0;i<numberRows_;i++) {
942    double value = region[i];
943    if (value) {
944      int k = permute_[i];
945      region[i]=0.0;
946      region2[k]=value;
947      regionIndex[numberNonZero++]=k;
948      mark_[k]=1;
949    }
950  }
951  // copy back
952  // set up linked lists at each depth
953  // stack2 is start, stack is next
954  int greatestDepth=-1;
955  int smallestDepth=numberRows_;
956  for (i=0;i<numberNonZero;i++) {
957    int j = regionIndex[i];
958    // add in
959    int iDepth = depth_[j];
[399]960    smallestDepth = CoinMin(iDepth,smallestDepth) ;
961    greatestDepth = CoinMax(iDepth,greatestDepth) ;
[130]962    int jNext = stack2_[iDepth];
963    stack2_[iDepth]=j;
964    stack_[j]=jNext;
965    // and put all descendants on list
966    int iChild = descendant_[j];
967    while (iChild>=0) {
968      if (!mark_[iChild]) {
969        regionIndex[numberNonZero++] = iChild;
970        mark_[iChild]=1;
971      }
972      iChild = rightSibling_[iChild];
973    }
974  }
975  numberNonZero=0;
976  region2[numberRows_]=0.0;
977  int iDepth;
978  for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
979    int iPivot = stack2_[iDepth];
980    stack2_[iDepth]=-1;
981    while (iPivot>=0) {
982      mark_[iPivot]=0;
983      double pivotValue = region2[iPivot];
984      int otherRow = parent_[iPivot];
985      double otherValue = region2[otherRow];
986      pivotValue = sign_[iPivot]*pivotValue+otherValue;
987      region2[iPivot]=pivotValue;
988      if (pivotValue) 
989        numberNonZero++;
990      iPivot = stack_[iPivot];
991    }
992  }
993  return numberNonZero;
[124]994}
995/* Updates one column (BTRAN) from region2 */
996int 
[130]997ClpNetworkBasis::updateColumnTranspose (  CoinIndexedVector * regionSparse,
[124]998                                        CoinIndexedVector * regionSparse2) const
999{
[130]1000  // permute in - presume small number so copy back
1001  // so will end up in right place
1002  regionSparse->clear (  );
1003  double *region = regionSparse->denseVector (  );
1004  double *region2 = regionSparse2->denseVector (  );
1005  int *regionIndex2 = regionSparse2->getIndices (  );
1006  int numberNonZero2 = regionSparse2->getNumElements (  );
1007  int *regionIndex = regionSparse->getIndices (  );
1008  int i;
1009  int numberNonZero=0;
[225]1010  bool packed = regionSparse2->packedMode();
1011  if (packed) {
1012    for (i=0;i<numberNonZero2;i++) {
1013      int k = regionIndex2[i];
1014      int j = permute_[k];
1015      double value = region2[i];
1016      region2[i]=0.0;
1017      region[j]=value;
1018      mark_[j]=1;
1019      regionIndex[numberNonZero++]=j;
1020    }
1021    // set up linked lists at each depth
1022    // stack2 is start, stack is next
1023    int greatestDepth=-1;
1024    int smallestDepth=numberRows_;
1025    //mark_[numberRows_]=1;
1026    for (i=0;i<numberNonZero2;i++) {
1027      int j = regionIndex[i];
1028      regionIndex2[i]=j;
1029      // add in
1030      int iDepth = depth_[j];
[399]1031      smallestDepth = CoinMin(iDepth,smallestDepth) ;
1032      greatestDepth = CoinMax(iDepth,greatestDepth) ;
[225]1033      int jNext = stack2_[iDepth];
1034      stack2_[iDepth]=j;
1035      stack_[j]=jNext;
1036      // and put all descendants on list
1037      int iChild = descendant_[j];
1038      while (iChild>=0) {
1039        if (!mark_[iChild]) {
1040          regionIndex2[numberNonZero++] = iChild;
1041          mark_[iChild]=1;
1042        }
1043        iChild = rightSibling_[iChild];
[130]1044      }
1045    }
[225]1046    for (;i<numberNonZero;i++) {
1047      int j = regionIndex2[i];
1048      // add in
1049      int iDepth = depth_[j];
[399]1050      smallestDepth = CoinMin(iDepth,smallestDepth) ;
1051      greatestDepth = CoinMax(iDepth,greatestDepth) ;
[225]1052      int jNext = stack2_[iDepth];
1053      stack2_[iDepth]=j;
1054      stack_[j]=jNext;
1055      // and put all descendants on list
1056      int iChild = descendant_[j];
1057      while (iChild>=0) {
1058        if (!mark_[iChild]) {
1059          regionIndex2[numberNonZero++] = iChild;
1060          mark_[iChild]=1;
1061        }
1062        iChild = rightSibling_[iChild];
[130]1063      }
1064    }
[225]1065    numberNonZero2=0;
1066    region[numberRows_]=0.0;
1067    int iDepth;
1068    for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
1069      int iPivot = stack2_[iDepth];
1070      stack2_[iDepth]=-1;
1071      while (iPivot>=0) {
1072        mark_[iPivot]=0;
1073        double pivotValue = region[iPivot];
1074        int otherRow = parent_[iPivot];
1075        double otherValue = region[otherRow];
1076        pivotValue = sign_[iPivot]*pivotValue+otherValue;
1077        region[iPivot]=pivotValue;
1078        if (pivotValue) {
1079          region2[numberNonZero2]=pivotValue;
1080          regionIndex2[numberNonZero2++]=iPivot;
1081        }
1082        iPivot = stack_[iPivot];
1083      }
[130]1084    }
[225]1085    // zero out
1086    for (i=0;i<numberNonZero2;i++) {
1087      int k = regionIndex2[i];
1088      region[k]=0.0;
1089    }
1090  } else {
1091    for (i=0;i<numberNonZero2;i++) {
1092      int k = regionIndex2[i];
1093      int j = permute_[k];
1094      double value = region2[k];
1095      region2[k]=0.0;
1096      region[j]=value;
1097      mark_[j]=1;
1098      regionIndex[numberNonZero++]=j;
1099    }
1100    // copy back
1101    // set up linked lists at each depth
1102    // stack2 is start, stack is next
1103    int greatestDepth=-1;
1104    int smallestDepth=numberRows_;
1105    //mark_[numberRows_]=1;
1106    for (i=0;i<numberNonZero2;i++) {
1107      int j = regionIndex[i];
1108      double value = region[j];
1109      region[j]=0.0;
1110      region2[j]=value;
1111      regionIndex2[i]=j;
1112      // add in
1113      int iDepth = depth_[j];
[399]1114      smallestDepth = CoinMin(iDepth,smallestDepth) ;
1115      greatestDepth = CoinMax(iDepth,greatestDepth) ;
[225]1116      int jNext = stack2_[iDepth];
1117      stack2_[iDepth]=j;
1118      stack_[j]=jNext;
1119      // and put all descendants on list
1120      int iChild = descendant_[j];
1121      while (iChild>=0) {
1122        if (!mark_[iChild]) {
1123          regionIndex2[numberNonZero++] = iChild;
1124          mark_[iChild]=1;
1125        }
1126        iChild = rightSibling_[iChild];
1127      }
1128    }
1129    for (;i<numberNonZero;i++) {
1130      int j = regionIndex2[i];
1131      // add in
1132      int iDepth = depth_[j];
[399]1133      smallestDepth = CoinMin(iDepth,smallestDepth) ;
1134      greatestDepth = CoinMax(iDepth,greatestDepth) ;
[225]1135      int jNext = stack2_[iDepth];
1136      stack2_[iDepth]=j;
1137      stack_[j]=jNext;
1138      // and put all descendants on list
1139      int iChild = descendant_[j];
1140      while (iChild>=0) {
1141        if (!mark_[iChild]) {
1142          regionIndex2[numberNonZero++] = iChild;
1143          mark_[iChild]=1;
1144        }
1145        iChild = rightSibling_[iChild];
1146      }
1147    }
1148    numberNonZero2=0;
1149    region2[numberRows_]=0.0;
1150    int iDepth;
1151    for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
1152      int iPivot = stack2_[iDepth];
1153      stack2_[iDepth]=-1;
1154      while (iPivot>=0) {
1155        mark_[iPivot]=0;
1156        double pivotValue = region2[iPivot];
1157        int otherRow = parent_[iPivot];
1158        double otherValue = region2[otherRow];
1159        pivotValue = sign_[iPivot]*pivotValue+otherValue;
1160        region2[iPivot]=pivotValue;
1161        if (pivotValue) 
1162          regionIndex2[numberNonZero2++]=iPivot;
1163        iPivot = stack_[iPivot];
1164      }
1165    }
[130]1166  }
1167  regionSparse2->setNumElements(numberNonZero2);
1168#ifdef FULL_DEBUG
1169 {
1170   int i;
1171   for (i=0;i<numberRows_;i++) {
1172     assert(!mark_[i]);
1173     assert (stack2_[i]==-1);
1174   }
1175 }
1176#endif
1177  return numberNonZero2;
[124]1178}
Note: See TracBrowser for help on using the repository browser.