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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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