source: trunk/Clp/src/ClpNetworkBasis.cpp

Last change on this file was 2385, checked in by unxusr, 8 months ago

formatting

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