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