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

Last change on this file since 1732 was 1732, checked in by forrest, 8 years ago

various fixes plus slightly weighted pricing plus lagomory switches

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 42.0 KB
Line 
1/* $Id: ClpNetworkBasis.cpp 1732 2011-05-31 08:09:41Z forrest $ */
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}
Note: See TracBrowser for help on using the repository browser.