source: tags/preroot/ClpNetworkBasis.cpp @ 778

Last change on this file since 778 was 131, checked in by forrest, 17 years ago

Network and faster updates

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.4 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 */
563int 
564ClpNetworkBasis::updateColumn (  CoinIndexedVector * regionSparse,
565                                 CoinIndexedVector * regionSparse2)
566{
567  regionSparse->clear (  );
568  double *region = regionSparse->denseVector (  );
569  double *region2 = regionSparse2->denseVector (  );
570  int *regionIndex2 = regionSparse2->getIndices (  );
571  int numberNonZero = regionSparse2->getNumElements (  );
572  int *regionIndex = regionSparse->getIndices (  );
573  int i;
574  bool doTwo = (numberNonZero==2);
575  int i0;
576  int i1;
577  if (doTwo) {
578    i0 = regionIndex2[0];
579    i1 = regionIndex2[1];
580    doTwo =  (region2[i0]*region2[i1]<0.0);
581  }
582  if (doTwo) {
583    // If just +- 1 then could go backwards on depth until join
584    region[i0] = region2[i0];
585    region2[i0]=0.0;
586    region[i1] = region2[i1];
587    region2[i1]=0.0;
588    int iDepth0 = depth_[i0];
589    int iDepth1 = depth_[i1];
590    if (iDepth1>iDepth0) {
591      int temp = i0;
592      i0 = i1;
593      i1 = temp;
594      temp = iDepth0;
595      iDepth0 = iDepth1;
596      iDepth1 = temp;
597    }
598    numberNonZero=0;
599    while (iDepth0>iDepth1) {
600      double pivotValue = region[i0];
601      // put back now ?
602      int iBack = permuteBack_[i0];
603      regionIndex2[numberNonZero++] = iBack;
604      int otherRow = parent_[i0];
605      region2[iBack] = pivotValue*sign_[i0];
606      region[i0] =0.0;
607      region[otherRow] += pivotValue;
608      iDepth0--;
609      i0 = otherRow;
610    }
611    while (i0!=i1) {
612      double pivotValue = region[i0];
613      // put back now ?
614      int iBack = permuteBack_[i0];
615      regionIndex2[numberNonZero++] = iBack;
616      int otherRow = parent_[i0];
617      region2[iBack] = pivotValue*sign_[i0];
618      region[i0] =0.0;
619      region[otherRow] += pivotValue;
620      i0 = otherRow;
621      double pivotValue1 = region[i1];
622      // put back now ?
623      int iBack1 = permuteBack_[i1];
624      regionIndex2[numberNonZero++] = iBack1;
625      int otherRow1 = parent_[i1];
626      region2[iBack1] = pivotValue1*sign_[i1];
627      region[i1] =0.0;
628      region[otherRow1] += pivotValue1;
629      i1 = otherRow1;
630    }
631  } else {
632    // set up linked lists at each depth
633    // stack2 is start, stack is next
634    int greatestDepth=-1;
635    //mark_[numberRows_]=1;
636    for (i=0;i<numberNonZero;i++) {
637      int j = regionIndex2[i];
638      double value = region2[j];
639      region2[j]=0.0;
640      region[j]=value;
641      regionIndex[i]=j;
642      int iDepth = depth_[j];
643      if (iDepth>greatestDepth) 
644        greatestDepth = iDepth;
645      // and back until marked
646      while (!mark_[j]) {
647        int iNext = stack2_[iDepth];
648        stack2_[iDepth]=j;
649        stack_[j]=iNext;
650        mark_[j]=1;
651        iDepth--;
652        j=parent_[j];
653      }
654    }
655    numberNonZero=0;
656    for (;greatestDepth>=0; greatestDepth--) {
657      int iPivot = stack2_[greatestDepth];
658      stack2_[greatestDepth]=-1;
659      while (iPivot>=0) {
660        mark_[iPivot]=0;
661        double pivotValue = region[iPivot];
662        if (pivotValue) {
663          // put back now ?
664          int iBack = permuteBack_[iPivot];
665          regionIndex2[numberNonZero++] = iBack;
666          int otherRow = parent_[iPivot];
667          region2[iBack] = pivotValue*sign_[iPivot];
668          region[iPivot] =0.0;
669          region[otherRow] += pivotValue;
670        }
671        iPivot = stack_[iPivot];
672      }
673    }
674  }
675  region[numberRows_]=0.0;
676  regionSparse2->setNumElements(numberNonZero);
677#ifdef FULL_DEBUG
678 {
679   int i;
680   for (i=0;i<numberRows_;i++) {
681     assert(!mark_[i]);
682     assert (stack2_[i]==-1);
683   }
684 }
685#endif
686  return numberNonZero;
687}
688/* Updates one column (FTRAN) to/from array
689    ** For large problems you should ALWAYS know where the nonzeros
690   are, so please try and migrate to previous method after you
691   have got code working using this simple method - thank you!
692   (the only exception is if you know input is dense e.g. rhs) */
693int 
694ClpNetworkBasis::updateColumn (  CoinIndexedVector * regionSparse,
695                                 double region2[] ) const
696{
697  regionSparse->clear (  );
698  double *region = regionSparse->denseVector (  );
699  int numberNonZero = 0;
700  int *regionIndex = regionSparse->getIndices (  );
701  int i;
702  // set up linked lists at each depth
703  // stack2 is start, stack is next
704  int greatestDepth=-1;
705  for (i=0;i<numberRows_;i++) {
706    double value = region2[i];
707    if (value) {
708      region2[i]=0.0;
709      region[i]=value;
710      regionIndex[numberNonZero++]=i;
711      int j=i;
712      int iDepth = depth_[j];
713      if (iDepth>greatestDepth) 
714        greatestDepth = iDepth;
715      // and back until marked
716      while (!mark_[j]) {
717        int iNext = stack2_[iDepth];
718        stack2_[iDepth]=j;
719        stack_[j]=iNext;
720        mark_[j]=1;
721        iDepth--;
722        j=parent_[j];
723      }
724    }
725  }
726  numberNonZero=0;
727  for (;greatestDepth>=0; greatestDepth--) {
728    int iPivot = stack2_[greatestDepth];
729    stack2_[greatestDepth]=-1;
730    while (iPivot>=0) {
731      mark_[iPivot]=0;
732      double pivotValue = region[iPivot];
733      if (pivotValue) {
734        // put back now ?
735        int iBack = permuteBack_[iPivot];
736        numberNonZero++;
737        int otherRow = parent_[iPivot];
738        region2[iBack] = pivotValue*sign_[iPivot];
739        region[iPivot] =0.0;
740        region[otherRow] += pivotValue;
741      }
742      iPivot = stack_[iPivot];
743    }
744  }
745  region[numberRows_]=0.0;
746  return numberNonZero;
747}
748/* Updates one column transpose (BTRAN)
749   For large problems you should ALWAYS know where the nonzeros
750   are, so please try and migrate to previous method after you
751   have got code working using this simple method - thank you!
752   (the only exception is if you know input is dense e.g. dense objective)
753   returns number of nonzeros */
754int 
755ClpNetworkBasis::updateColumnTranspose (  CoinIndexedVector * regionSparse,
756                                          double region2[] ) const
757{
758  // permute in after copying
759  // so will end up in right place
760  double *region = regionSparse->denseVector (  );
761  int *regionIndex = regionSparse->getIndices (  );
762  int i;
763  int numberNonZero=0;
764  memcpy(region,region2,numberRows_*sizeof(double));
765  for (i=0;i<numberRows_;i++) {
766    double value = region[i];
767    if (value) {
768      int k = permute_[i];
769      region[i]=0.0;
770      region2[k]=value;
771      regionIndex[numberNonZero++]=k;
772      mark_[k]=1;
773    }
774  }
775  // copy back
776  // set up linked lists at each depth
777  // stack2 is start, stack is next
778  int greatestDepth=-1;
779  int smallestDepth=numberRows_;
780  for (i=0;i<numberNonZero;i++) {
781    int j = regionIndex[i];
782    // add in
783    int iDepth = depth_[j];
784    smallestDepth = min(iDepth,smallestDepth) ;
785    greatestDepth = max(iDepth,greatestDepth) ;
786    int jNext = stack2_[iDepth];
787    stack2_[iDepth]=j;
788    stack_[j]=jNext;
789    // and put all descendants on list
790    int iChild = descendant_[j];
791    while (iChild>=0) {
792      if (!mark_[iChild]) {
793        regionIndex[numberNonZero++] = iChild;
794        mark_[iChild]=1;
795      }
796      iChild = rightSibling_[iChild];
797    }
798  }
799  numberNonZero=0;
800  region2[numberRows_]=0.0;
801  int iDepth;
802  for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
803    int iPivot = stack2_[iDepth];
804    stack2_[iDepth]=-1;
805    while (iPivot>=0) {
806      mark_[iPivot]=0;
807      double pivotValue = region2[iPivot];
808      int otherRow = parent_[iPivot];
809      double otherValue = region2[otherRow];
810      pivotValue = sign_[iPivot]*pivotValue+otherValue;
811      region2[iPivot]=pivotValue;
812      if (pivotValue) 
813        numberNonZero++;
814      iPivot = stack_[iPivot];
815    }
816  }
817  return numberNonZero;
818}
819/* Updates one column (BTRAN) from region2 */
820int 
821ClpNetworkBasis::updateColumnTranspose (  CoinIndexedVector * regionSparse,
822                                        CoinIndexedVector * regionSparse2) const
823{
824  // permute in - presume small number so copy back
825  // so will end up in right place
826  regionSparse->clear (  );
827  double *region = regionSparse->denseVector (  );
828  double *region2 = regionSparse2->denseVector (  );
829  int *regionIndex2 = regionSparse2->getIndices (  );
830  int numberNonZero2 = regionSparse2->getNumElements (  );
831  int *regionIndex = regionSparse->getIndices (  );
832  int i;
833  int numberNonZero=0;
834  for (i=0;i<numberNonZero2;i++) {
835    int k = regionIndex2[i];
836    int j = permute_[k];
837    double value = region2[k];
838    region2[k]=0.0;
839    region[j]=value;
840    mark_[j]=1;
841    regionIndex[numberNonZero++]=j;
842  }
843  // copy back
844  // set up linked lists at each depth
845  // stack2 is start, stack is next
846  int greatestDepth=-1;
847  int smallestDepth=numberRows_;
848  //mark_[numberRows_]=1;
849  for (i=0;i<numberNonZero2;i++) {
850    int j = regionIndex[i];
851    double value = region[j];
852    region[j]=0.0;
853    region2[j]=value;
854    regionIndex2[i]=j;
855    // add in
856    int iDepth = depth_[j];
857    smallestDepth = min(iDepth,smallestDepth) ;
858    greatestDepth = max(iDepth,greatestDepth) ;
859    int jNext = stack2_[iDepth];
860    stack2_[iDepth]=j;
861    stack_[j]=jNext;
862    // and put all descendants on list
863    int iChild = descendant_[j];
864    while (iChild>=0) {
865      if (!mark_[iChild]) {
866        regionIndex2[numberNonZero++] = iChild;
867        mark_[iChild]=1;
868      }
869      iChild = rightSibling_[iChild];
870    }
871  }
872  for (;i<numberNonZero;i++) {
873    int j = regionIndex2[i];
874    // add in
875    int iDepth = depth_[j];
876    smallestDepth = min(iDepth,smallestDepth) ;
877    greatestDepth = max(iDepth,greatestDepth) ;
878    int jNext = stack2_[iDepth];
879    stack2_[iDepth]=j;
880    stack_[j]=jNext;
881    // and put all descendants on list
882    int iChild = descendant_[j];
883    while (iChild>=0) {
884      if (!mark_[iChild]) {
885        regionIndex2[numberNonZero++] = iChild;
886        mark_[iChild]=1;
887      }
888      iChild = rightSibling_[iChild];
889    }
890  }
891  numberNonZero2=0;
892  region2[numberRows_]=0.0;
893  int iDepth;
894  for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
895    int iPivot = stack2_[iDepth];
896    stack2_[iDepth]=-1;
897    while (iPivot>=0) {
898      mark_[iPivot]=0;
899      double pivotValue = region2[iPivot];
900      int otherRow = parent_[iPivot];
901      double otherValue = region2[otherRow];
902      pivotValue = sign_[iPivot]*pivotValue+otherValue;
903      region2[iPivot]=pivotValue;
904      if (pivotValue) 
905        regionIndex2[numberNonZero2++]=iPivot;
906      iPivot = stack_[iPivot];
907    }
908  }
909  regionSparse2->setNumElements(numberNonZero2);
910#ifdef FULL_DEBUG
911 {
912   int i;
913   for (i=0;i<numberRows_;i++) {
914     assert(!mark_[i]);
915     assert (stack2_[i]==-1);
916   }
917 }
918#endif
919  return numberNonZero2;
920}
Note: See TracBrowser for help on using the repository browser.