source: trunk/ClpNetworkBasis.cpp @ 225

Last change on this file since 225 was 225, checked in by forrest, 16 years ago

This should break everything

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