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

Last change on this file since 1197 was 1197, checked in by forrest, 11 years ago

many changes to try and improve performance

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