1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkPKdTree.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 /*----------------------------------------------------------------------------
16 Copyright (c) Sandia Corporation
17 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
18 ----------------------------------------------------------------------------*/
19
20 #include "vtkPKdTree.h"
21 #include "vtkKdNode.h"
22 #include "vtkDataSet.h"
23 #include "vtkCellCenters.h"
24 #include "vtkPoints.h"
25 #include "vtkUnstructuredGrid.h"
26 #include "vtkObjectFactory.h"
27 #include "vtkMultiProcessController.h"
28 #include "vtkSocketController.h"
29 #include "vtkTimerLog.h"
30 #include "vtkCellData.h"
31 #include "vtkPointData.h"
32 #include "vtkIntArray.h"
33 #include "vtkIdList.h"
34 #include "vtkSubGroup.h"
35 #include "vtkCommand.h"
36
37 #include <queue>
38 #include <algorithm>
39 #include <cassert>
40
41 // Timing data ---------------------------------------------
42
43 #if 0
44 #define MSGSIZE 60
45
46 static char dots[MSGSIZE] = "...........................................................";
47 static char msg[MSGSIZE];
48
49 static char * makeEntry(const char *s)
50 {
51 memcpy(msg, dots, MSGSIZE);
52 int len = strlen(s);
53 len = (len >= MSGSIZE) ? MSGSIZE-1 : len;
54
55 memcpy(msg, s, len);
56
57 return msg;
58 }
59
60 #define TIMER(s) \
61 if (this->GetTiming()) \
62 { \
63 char *s2 = makeEntry(s); \
64 if (this->TimerLog == NULL) \
65 { \
66 this->TimerLog = vtkTimerLog::New(); \
67 } \
68 this->TimerLog->MarkStartEvent(s2); \
69 }
70
71 #define TIMERDONE(s) \
72 if (this->GetTiming())\
73 { char *s2 = makeEntry(s); this->TimerLog->MarkEndEvent(s2); }
74
75 #else
76 #define TIMER(s)
77 #define TIMERDONE(s)
78 #endif
79
80 // Timing data ---------------------------------------------
81
82 vtkStandardNewMacro(vtkPKdTree);
83
84 const int vtkPKdTree::NoRegionAssignment = 0; // default
85 const int vtkPKdTree::ContiguousAssignment = 1; // default if RegionAssignmentOn
86 const int vtkPKdTree::UserDefinedAssignment = 2;
87 const int vtkPKdTree::RoundRobinAssignment = 3;
88
89 #define FreeList(list) if (list) {delete [] list; list = NULL;}
90 #define FreeObject(item) if (item) {item->Delete(); item = NULL;}
91
92
93 #define VTKERROR(s) \
94 { \
95 vtkErrorMacro(<< "(process " << this->MyId << ") " << s);\
96 }
97 #define VTKWARNING(s) \
98 { \
99 vtkWarningMacro(<< "(process " << this->MyId << ") " << s);\
100 }
101
vtkPKdTree()102 vtkPKdTree::vtkPKdTree()
103 {
104 this->RegionAssignment = ContiguousAssignment;
105
106 this->Controller = NULL;
107 this->SubGroup = NULL;
108
109 this->NumProcesses = 1;
110 this->MyId = 0;
111
112 this->InitializeRegionAssignmentLists();
113 this->InitializeProcessDataLists();
114 this->InitializeFieldArrayMinMax();
115 this->InitializeGlobalIndexLists();
116
117 this->TotalNumCells = 0;
118
119 this->PtArray = NULL;
120 this->PtArray2 = NULL;
121 this->CurrentPtArray = NULL;
122 this->NextPtArray = NULL;
123
124 this->SelectBuffer = NULL;
125 }
~vtkPKdTree()126 vtkPKdTree::~vtkPKdTree()
127 {
128 this->SetController(NULL);
129 this->FreeSelectBuffer();
130 this->FreeDoubleBuffer();
131
132 this->FreeGlobalIndexLists();
133 this->FreeRegionAssignmentLists();
134 this->FreeProcessDataLists();
135 this->FreeFieldArrayMinMax();
136 }
SetController(vtkMultiProcessController * c)137 void vtkPKdTree::SetController(vtkMultiProcessController *c)
138 {
139 if (this->Controller == c)
140 {
141 return;
142 }
143
144 if ((c == NULL) || (c->GetNumberOfProcesses() == 0))
145 {
146 this->NumProcesses = 1;
147 this->MyId = 0;
148 }
149
150 this->Modified();
151
152 if (this->Controller != NULL)
153 {
154 this->Controller->UnRegister(this);
155 this->Controller = NULL;
156 }
157
158 if (c == NULL)
159 {
160 return;
161 }
162
163 vtkSocketController *sc = vtkSocketController::SafeDownCast(c);
164
165 if (sc)
166 {
167 vtkErrorMacro(<<
168 "vtkPKdTree communication will fail with a socket controller");
169
170 return;
171 }
172
173 this->NumProcesses = c->GetNumberOfProcesses();
174
175 this->Controller = c;
176 this->MyId = c->GetLocalProcessId();
177 c->Register(this);
178 }
179 //--------------------------------------------------------------------
180 // Parallel k-d tree build, Floyd and Rivest (1975) select algorithm
181 // for median finding.
182 //--------------------------------------------------------------------
183
AllCheckForFailure(int rc,const char * where,const char * how)184 int vtkPKdTree::AllCheckForFailure(int rc, const char *where, const char *how)
185 {
186 int vote;
187 char errmsg[256];
188
189 if (this->NumProcesses > 1){
190 this->SubGroup->ReduceSum(&rc, &vote, 1, 0);
191 this->SubGroup->Broadcast(&vote, 1, 0);
192 }
193 else{
194 vote = rc;
195 }
196
197 if (vote)
198 {
199 if (rc)
200 {
201 sprintf(errmsg,"%s on my node (%s)",how, where);
202 }
203 else
204 {
205 sprintf(errmsg,"%s on a remote node (%s)",how, where);
206 }
207 VTKWARNING(errmsg);
208
209 return 1;
210 }
211 return 0;
212 }
213
AllCheckParameters()214 void vtkPKdTree::AllCheckParameters()
215 {
216 int param[10];
217 int param0[10];
218
219 // All the parameters that determine how k-d tree is built and
220 // what tables get created afterward - there's no point in
221 // trying to build unless these match on all processes.
222
223 param[0] = this->ValidDirections;
224 param[1] = this->GetMinCells();
225 param[2] = this->GetNumberOfRegionsOrLess();
226 param[3] = this->GetNumberOfRegionsOrMore();
227 param[4] = this->RegionAssignment;
228 param[5] = 0;
229 param[6] = 0;
230 param[7] = 0;
231 param[8] = 0;
232 param[9] = 0;
233
234 if (this->MyId == 0)
235 {
236 this->SubGroup->Broadcast(param, 10, 0);
237 return;
238 }
239
240 this->SubGroup->Broadcast(param0, 10, 0);
241
242 int diff = 0;
243
244 for (int i=0; i < 10; i++)
245 {
246 if (param0[i] != param[i])
247 {
248 diff = 1;
249 break;
250 }
251 }
252 if (diff)
253 {
254 VTKWARNING("Changing my runtime parameters to match process 0");
255
256 this->ValidDirections = param0[0];
257 this->SetMinCells(param0[1]);
258 this->SetNumberOfRegionsOrLess(param0[2]);
259 this->SetNumberOfRegionsOrMore(param0[3]);
260 this->RegionAssignment = param0[4];
261 }
262 return;
263 }
264
265 #define BoundsToMinMax(bounds,min,max) \
266 { \
267 min[0] = bounds[0]; min[1] = bounds[2]; min[2] = bounds[4]; \
268 max[0] = bounds[1]; max[1] = bounds[3]; max[2] = bounds[5]; \
269 }
270 #define MinMaxToBounds(bounds,min,max) \
271 { \
272 bounds[0] = min[0]; bounds[2] = min[1]; bounds[4] = min[2]; \
273 bounds[1] = max[0]; bounds[3] = max[1]; bounds[5] = max[2]; \
274 }
275 #define BoundsToMinMaxUpdate(bounds,min,max) \
276 { \
277 min[0] = (bounds[0] < min[0] ? bounds[0] : min[0]); \
278 min[1] = (bounds[2] < min[1] ? bounds[2] : min[1]); \
279 min[2] = (bounds[4] < min[2] ? bounds[4] : min[2]); \
280 max[0] = (bounds[1] > max[0] ? bounds[1] : max[0]); \
281 max[1] = (bounds[3] > max[1] ? bounds[3] : max[1]); \
282 max[2] = (bounds[5] > max[2] ? bounds[5] : max[2]); \
283 }
284
VolumeBounds(double * volBounds)285 bool vtkPKdTree::VolumeBounds(double* volBounds)
286 {
287 int i;
288
289 // Get the spatial bounds of the whole volume
290 double localMin[3], localMax[3], globalMin[3], globalMax[3];
291
292 int number_of_datasets = this->GetNumberOfDataSets();
293 if (number_of_datasets == 0)
294 {
295 VTKERROR("NumberOfDatasets = 0, cannot determine volume bounds.");
296 return false;
297 }
298
299 for (i=0; i < number_of_datasets; i++)
300 {
301 this->GetDataSet(i)->GetBounds(volBounds);
302
303 if (i==0)
304 {
305 BoundsToMinMax(volBounds, localMin, localMax);
306 }
307 else
308 {
309 BoundsToMinMaxUpdate(volBounds, localMin, localMax);
310 }
311 }
312
313 // trick to reduce the number of global communications for getting both
314 // min and max
315 double localReduce[6], globalReduce[6];
316 for(i=0;i<3;i++)
317 {
318 localReduce[i] = localMin[i];
319 localReduce[i+3] = -localMax[i];
320 }
321 this->SubGroup->ReduceMin(localReduce, globalReduce, 6, 0);
322 this->SubGroup->Broadcast(globalReduce, 6, 0);
323
324 for(i=0;i<3;i++)
325 {
326 globalMin[i] = globalReduce[i];
327 globalMax[i] = -globalReduce[i+3];
328 }
329
330
331 MinMaxToBounds(volBounds, globalMin, globalMax);
332
333 // push out a little if flat
334
335 double diff[3], aLittle = 0.0;
336
337 for (i=0; i<3; i++)
338 {
339 diff[i] = volBounds[2*i+1] - volBounds[2*i];
340 aLittle = (diff[i] > aLittle) ? diff[i] : aLittle;
341 }
342 if ((aLittle /= 100.0) <= 0.0)
343 {
344 VTKERROR("VolumeBounds - degenerate volume");
345 return false;
346 }
347
348 this->FudgeFactor = aLittle * 10e-4;
349
350 for (i=0; i<3; i++)
351 {
352 if (diff[i] <= 0)
353 {
354 volBounds[2*i] -= aLittle;
355 volBounds[2*i+1] += aLittle;
356 }
357 else
358 {
359 volBounds[2*i] -= this->GetFudgeFactor();
360 volBounds[2*i+1] += this->GetFudgeFactor();
361 }
362 }
363 return true;
364 }
365
366 // BuildLocator must be called by all processes in the parallel application
367
BuildLocator()368 void vtkPKdTree::BuildLocator()
369 {
370 int fail = 0;
371 int rebuildLocator = 0;
372
373 if ((this->Top == NULL) ||
374 (this->BuildTime < this->GetMTime()) ||
375 this->NewGeometry())
376 {
377 // We don't have a k-d tree, or parameters that affect the
378 // build of the tree have changed, or input geometry has changed.
379
380 rebuildLocator = 1;
381 }
382
383 if (this->NumProcesses == 1)
384 {
385 if (rebuildLocator)
386 {
387 this->SingleProcessBuildLocator();
388 }
389 return;
390 }
391 this->UpdateProgress(0);
392
393 TIMER("Determine if we need to rebuild");
394
395 this->SubGroup = vtkSubGroup::New();
396 this->SubGroup->Initialize(0, this->NumProcesses-1,
397 this->MyId, 0x00001000, this->Controller->GetCommunicator());
398
399 int vote;
400 this->SubGroup->ReduceSum(&rebuildLocator, &vote, 1, 0);
401 this->SubGroup->Broadcast(&vote, 1, 0);
402
403 rebuildLocator = (vote > 0);
404
405 TIMERDONE("Determine if we need to rebuild");
406
407 if (rebuildLocator)
408 {
409 TIMER("Build k-d tree");
410 this->InvokeEvent(vtkCommand::StartEvent);
411
412 this->FreeSearchStructure();
413 this->ReleaseTables();
414
415 this->AllCheckParameters(); // global operation to ensure same parameters
416
417 double volBounds[6];
418 if(this->VolumeBounds(volBounds) == false) // global operation to get bounds
419 {
420 goto doneError;
421 }
422 this->UpdateProgress(0.1);
423
424 if (this->UserDefinedCuts)
425 {
426 fail = this->ProcessUserDefinedCuts(volBounds);
427 }
428 else
429 {
430 fail = this->MultiProcessBuildLocator(volBounds);
431 }
432
433 if (fail) goto doneError;
434
435 this->SetActualLevel();
436 this->BuildRegionList();
437
438 TIMERDONE("Build k-d tree");
439 this->InvokeEvent(vtkCommand::EndEvent);
440 }
441
442 // Even if locator is not rebuilt, we should update
443 // region assignments since they may have changed.
444
445 this->UpdateRegionAssignment();
446
447 goto done;
448
449 doneError:
450
451 this->FreeRegionAssignmentLists();
452 this->FreeSearchStructure();
453
454 done:
455
456 FreeObject(this->SubGroup);
457
458 this->SetCalculator(this->Top);
459
460 this->UpdateBuildTime();
461 this->UpdateProgress(1.0);
462 return;
463 }
MultiProcessBuildLocator(double * volBounds)464 int vtkPKdTree::MultiProcessBuildLocator(double *volBounds)
465 {
466 int retVal = 0;
467
468 vtkDebugMacro( << "Creating Kdtree in parallel" );
469
470 if (this->GetTiming())
471 {
472 if (this->TimerLog == NULL) this->TimerLog = vtkTimerLog::New();
473 }
474
475 // Locally, create a single list of the coordinates of the centers of the
476 // cells of my data sets
477
478 TIMER("Compute cell centers");
479
480 this->PtArray = NULL;
481
482 this->ProgressOffset = 0.1;
483 this->ProgressScale = 0.5;
484
485 this->PtArray = this->ComputeCellCenters();
486 vtkIdType totalPts = this->GetNumberOfCells(); // total on local node
487 this->CurrentPtArray = this->PtArray;
488
489 // int fail = (this->PtArray == NULL);
490 int fail = ((this->PtArray == NULL) && (totalPts > 0));
491
492 if (this->AllCheckForFailure(fail,
493 "MultiProcessBuildLocator", "memory allocation"))
494 {
495 goto doneError6;
496 }
497
498 TIMERDONE("Compute cell centers");
499
500 // Get total number of cells across all processes, assign global indices
501 // for select operation
502
503 TIMER("Build index lists");
504
505 fail = this->BuildGlobalIndexLists(totalPts);
506 this->UpdateProgress(0.7);
507
508 TIMERDONE("Build index lists");
509
510 if (fail)
511 {
512 goto doneError6;
513 }
514
515 // In parallel, build the k-d tree structure, partitioning all
516 // the points into spatial regions. Sub-groups of processors
517 // will form vtkSubGroups to divide sub-regions of space.
518
519 FreeObject(this->SubGroup);
520
521 TIMER("Compute tree");
522
523 fail = this->BreadthFirstDivide(volBounds);
524 this->UpdateProgress(0.9);
525
526 TIMERDONE("Compute tree");
527
528 this->SubGroup = vtkSubGroup::New();
529 this->SubGroup->Initialize(0, this->NumProcesses-1,
530 this->MyId, 0x00002000, this->Controller->GetCommunicator());
531
532 if (this->AllCheckForFailure(fail, "BreadthFirstDivide", "memory allocation"))
533 {
534 goto doneError6;
535 }
536
537 FreeObject(this->SubGroup);
538
539 // I only have a partial tree at this point, the regions in which
540 // I participated. Now collect the entire tree.
541
542 this->SubGroup = vtkSubGroup::New();
543 this->SubGroup->Initialize(0, this->NumProcesses-1,
544 this->MyId, 0x00003000, this->Controller->GetCommunicator());
545
546 TIMER("Complete tree");
547
548 fail = this->CompleteTree();
549
550 TIMERDONE("Complete tree");
551
552 if (fail)
553 {
554 goto doneError6;
555 }
556
557 goto done6;
558
559 doneError6:
560
561 this->FreeSearchStructure();
562 retVal = 1;
563
564 done6:
565 // no longer valid, we overwrote them during k-d tree parallel build
566 delete [] this->PtArray;
567 this->CurrentPtArray = this->PtArray = NULL;
568
569 FreeObject(this->SubGroup);
570
571 this->FreeGlobalIndexLists();
572
573 return retVal;
574 }
575
SingleProcessBuildLocator()576 void vtkPKdTree::SingleProcessBuildLocator()
577 {
578 vtkKdTree::BuildLocator();
579
580 this->TotalNumCells = this->GetNumberOfCells();
581
582 if (this->RegionAssignment != vtkPKdTree::NoRegionAssignment)
583 {
584 this->UpdateRegionAssignment();
585 }
586
587 return;
588 }
589 typedef struct _vtkNodeInfo{
590 vtkKdNode *kd;
591 int L;
592 int level;
593 int tag;
594 } *vtkNodeInfo;
595
596 #define ENQUEUE(a, b, c, d) \
597 { \
598 vtkNodeInfo rec = new struct _vtkNodeInfo; \
599 rec->kd = a; \
600 rec->L = b; \
601 rec->level = c; \
602 rec->tag = d; \
603 Queue.push(rec); \
604 }
605
BreadthFirstDivide(double * volBounds)606 int vtkPKdTree::BreadthFirstDivide(double *volBounds)
607 {
608 int returnVal = 0;
609
610 std::queue <vtkNodeInfo> Queue;
611
612 if (this->AllocateDoubleBuffer())
613 {
614 VTKERROR("memory allocation for double buffering");
615 return 1;
616 }
617
618 if (this->AllocateSelectBuffer())
619 {
620 this->FreeDoubleBuffer();
621
622 VTKERROR("memory allocation for select buffers");
623 return 1;
624 }
625
626 vtkKdNode *kd = this->Top = vtkKdNode::New();
627
628 kd->SetBounds(volBounds[0], volBounds[1],
629 volBounds[2], volBounds[3],
630 volBounds[4], volBounds[5]);
631
632 kd->SetNumberOfPoints(this->TotalNumCells);
633
634 kd->SetDataBounds(volBounds[0], volBounds[1],
635 volBounds[2], volBounds[3],
636 volBounds[4], volBounds[5]);
637
638 int midpt = this->DivideRegion(kd, 0, 0, 0x00000001);
639
640 if (midpt >= 0)
641 {
642 ENQUEUE(kd->GetLeft(), 0, 1, 0x00000002);
643 ENQUEUE(kd->GetRight(), midpt, 1, 0x00000003);
644 }
645 else if (midpt < -1)
646 {
647 this->FreeSelectBuffer();
648 this->FreeDoubleBuffer();
649
650 return 1;
651 }
652
653 while (!Queue.empty())
654 {
655 vtkNodeInfo info = Queue.front();
656 Queue.pop();
657
658 kd = info->kd;
659 int L = info->L;
660 int level = info->level;
661 int tag = info->tag;
662
663 midpt = this->DivideRegion(kd, L, level, tag);
664
665 if (midpt >= 0)
666 {
667 ENQUEUE(kd->GetLeft(), L, level+1, tag << 1);
668
669 ENQUEUE(kd->GetRight(), midpt, level+1, (tag << 1) | 1);
670 }
671 else if (midpt < -1)
672 {
673 returnVal = 1; // have to keep going, or remote ops may hang
674 }
675 delete info;
676 }
677
678 this->FreeSelectBuffer();
679
680 if (this->CurrentPtArray == this->PtArray2)
681 {
682 memcpy(this->PtArray, this->PtArray2, this->PtArraySize * sizeof(float));
683 }
684
685 this->FreeDoubleBuffer();
686
687 return returnVal;
688 }
DivideRegion(vtkKdNode * kd,int L,int level,int tag)689 int vtkPKdTree::DivideRegion(vtkKdNode *kd, int L, int level, int tag)
690 {
691 if (!this->DivideTest(kd->GetNumberOfPoints(), level)) return -1;
692
693 int numpoints = kd->GetNumberOfPoints();
694 int R = L + numpoints - 1;
695
696 if (numpoints < 2)
697 {
698 // Special case: not enough points to go around.
699 int p = this->WhoHas(L);
700 if (this->MyId != p) return -1;
701
702 int maxdim = this->SelectCutDirection(kd);
703 kd->SetDim(maxdim);
704
705 vtkKdNode *left = vtkKdNode::New();
706 vtkKdNode *right = vtkKdNode::New();
707 kd->AddChildNodes(left, right);
708
709 double bounds[6];
710 kd->GetBounds(bounds);
711
712 float *val = this->GetLocalVal(L);
713
714 double coord;
715 if (numpoints > 0)
716 {
717 coord = val[maxdim];
718 }
719 else
720 {
721 coord = (bounds[maxdim*2] + bounds[maxdim*2+1])*0.5;
722 }
723
724 left->SetBounds(bounds[0], ((maxdim == XDIM) ? coord : bounds[1]),
725 bounds[2], ((maxdim == YDIM) ? coord : bounds[3]),
726 bounds[4], ((maxdim == ZDIM) ? coord : bounds[5]));
727
728 left->SetNumberOfPoints(numpoints);
729
730 right->SetBounds(((maxdim == XDIM) ? coord : bounds[0]), bounds[1],
731 ((maxdim == YDIM) ? coord : bounds[2]), bounds[3],
732 ((maxdim == ZDIM) ? coord : bounds[4]), bounds[5]);
733
734 right->SetNumberOfPoints(0);
735
736 // Set the data bounds tightly around L. This will inevitably mean some
737 // regions that are empty will have their data bounds outside of them.
738 // Hopefully, that will not screw up anything down the road.
739 left ->SetDataBounds(val[0], val[0], val[1], val[1], val[2], val[2]);
740 right->SetDataBounds(val[0], val[0], val[1], val[1], val[2], val[2]);
741
742 // Return L as the midpoint to guarantee that both left and right trees
743 // are "owned" by the same process as the parent. This is important
744 // because only one process has not culled this node in the tree.
745 return L;
746 }
747
748 int p1 = this->WhoHas(L);
749 int p2 = this->WhoHas(R);
750
751 if ((this->MyId < p1) || (this->MyId > p2))
752 {
753 return -1;
754 }
755
756 this->SubGroup = vtkSubGroup::New();
757 this->SubGroup->Initialize(p1, p2, this->MyId, tag,
758 this->Controller->GetCommunicator());
759
760 int maxdim = this->SelectCutDirection(kd);
761
762 kd->SetDim(maxdim);
763
764 int midpt = this->Select(maxdim, L, R);
765
766 if (midpt < L + 1)
767 {
768 // Couldn't divide. Try a different direction.
769 int newdim = vtkKdTree::XDIM - 1;
770 vtkDebugMacro(<< "Could not divide along maxdim"
771 << " maxdim " << maxdim
772 << " L " << L << " R " << R << " midpt " << midpt);
773 while (midpt < L + 1)
774 {
775 do
776 {
777 newdim++;
778 if (newdim > vtkKdTree::ZDIM)
779 {
780 // Exhausted all possible divisions. All points must be at same
781 // location. Just split in the middle and hope for the best.
782 vtkDebugMacro(<< "Must have coincident points.");
783 newdim = maxdim;
784 kd->SetDim(maxdim);
785 // Add one to make sure there is always something to the left.
786 midpt = (L+R)/2 + 1;
787 goto FindMidptBreakout;
788 }
789 } while ( (newdim == maxdim)
790 || ((this->ValidDirections & (1 << newdim)) == 0) );
791 kd->SetDim(newdim);
792 midpt = this->Select(newdim, L, R);
793 vtkDebugMacro(<< " newdim " << newdim
794 << " L " << L << " R " << R << " midpt " << midpt);
795 }
796 FindMidptBreakout:
797 // Pretend the dimension we used was the minimum.
798 maxdim = newdim;
799 }
800
801 float *newDataBounds = this->DataBounds(L, midpt, R);
802 vtkKdNode *left = vtkKdNode::New();
803 vtkKdNode *right = vtkKdNode::New();
804
805 int fail = ( (newDataBounds == NULL) || (left == NULL) || (right == NULL) );
806
807 if (this->AllCheckForFailure(fail, "Divide Region", "memory allocation"))
808 {
809 FreeList(newDataBounds);
810 left->Delete();
811 right->Delete();
812 FreeObject(this->SubGroup);
813 return -3;
814 }
815
816 double coord = (newDataBounds[maxdim*2 + 1] + // max on left side
817 newDataBounds[6 + maxdim*2] )* // min on right side
818 0.5;
819
820 kd->AddChildNodes(left, right);
821
822 double bounds[6];
823 kd->GetBounds(bounds);
824
825 left->SetBounds(
826 bounds[0], ((maxdim == XDIM) ? coord : bounds[1]),
827 bounds[2], ((maxdim == YDIM) ? coord : bounds[3]),
828 bounds[4], ((maxdim == ZDIM) ? coord : bounds[5]));
829
830 left->SetNumberOfPoints(midpt - L);
831
832 right->SetBounds(
833 ((maxdim == XDIM) ? coord : bounds[0]), bounds[1],
834 ((maxdim == YDIM) ? coord : bounds[2]), bounds[3],
835 ((maxdim == ZDIM) ? coord : bounds[4]), bounds[5]);
836
837 right->SetNumberOfPoints(R - midpt + 1);
838
839 left->SetDataBounds(newDataBounds[0], newDataBounds[1],
840 newDataBounds[2], newDataBounds[3],
841 newDataBounds[4], newDataBounds[5]);
842
843 right->SetDataBounds(newDataBounds[6], newDataBounds[7],
844 newDataBounds[8], newDataBounds[9],
845 newDataBounds[10], newDataBounds[11]);
846
847 delete [] newDataBounds;
848
849 FreeObject(this->SubGroup);
850
851 return midpt;
852 }
853
ExchangeVals(int pos1,int pos2)854 void vtkPKdTree::ExchangeVals(int pos1, int pos2)
855 {
856 vtkCommunicator *comm = this->Controller->GetCommunicator();
857
858 float *myval;
859 float otherval[3];
860
861 int player1 = this->WhoHas(pos1);
862 int player2 = this->WhoHas(pos2);
863
864 if ((player1 == this->MyId) && (player2 == this->MyId))
865 {
866 this->ExchangeLocalVals(pos1, pos2);
867 }
868
869 else if (player1 == this->MyId)
870 {
871 myval = this->GetLocalVal(pos1);
872
873 comm->Send(myval, 3, player2, this->SubGroup->tag);
874
875 comm->Receive(otherval, 3, player2, this->SubGroup->tag);
876
877 this->SetLocalVal(pos1, otherval);
878 }
879 else if (player2 == this->MyId)
880 {
881 myval = this->GetLocalVal(pos2);
882
883 comm->Receive(otherval, 3, player1, this->SubGroup->tag);
884
885 comm->Send(myval, 3, player1, this->SubGroup->tag);
886
887 this->SetLocalVal(pos2, otherval);
888 }
889 return;
890 }
891
892 // Given an array X with element indices ranging from L to R, and
893 // a K such that L <= K <= R, rearrange the elements such that
894 // X[K] contains the ith sorted element, where i = K - L + 1, and
895 // all the elements X[j], j < k satisfy X[j] <= X[K], and all the
896 // elements X[j], j > k satisfy X[j] >= X[K].
897
898 #define sign(x) ((x<0) ? (-1) : (1))
899 #ifndef max
900 #define max(x,y) ((x>y) ? (x) : (y))
901 #endif
902 #ifndef min
903 #define min(x,y) ((x<y) ? (x) : (y))
904 #endif
905
_select(int L,int R,int K,int dim)906 void vtkPKdTree::_select(int L, int R, int K, int dim)
907 {
908 int N, I, J, S, SD, LL, RR;
909 float Z;
910
911 while (R > L)
912 {
913 if ( R - L > 600)
914 {
915 // "Recurse on a sample of size S to get an estimate for the
916 // (K-L+1)-th smallest element into X[K], biased slightly so
917 // that the (K-L+1)-th element is expected to lie in the
918 // smaller set after partitioning"
919
920 N = R - L + 1;
921 I = K - L + 1;
922 Z = static_cast<float>(log(float(N)));
923 S = static_cast<int>(.5 * exp(2*Z/3));
924 SD = static_cast<int>(.5 * sqrt(Z*S*((float)(N-S)/N)) * sign(I - N/2));
925 LL = max(L, K - static_cast<int>((I*((float)S/N))) + SD);
926 RR = min(R, K + static_cast<int>((N-I) * ((float)S/N)) + SD);
927 this->_select(LL, RR, K, dim);
928 }
929
930 int p1 = this->WhoHas(L);
931 int p2 = this->WhoHas(R);
932
933 // "now adjust L,R so they surround the subset containing
934 // the (K-L+1)-th smallest element"
935
936 // Due to very severe worst case behavior when the
937 // value at K (call it "T") is repeated many times in the array, we
938 // rearrange the array into three intervals: the leftmost being values
939 // less than T, the center being values equal to T, and the rightmost
940 // being values greater than T. Two integers are returned. This first
941 // is the global index of the start of the second interval. The second
942 // is the global index of the start of the third interval. (If there
943 // are no values greater than "T", the second integer will be R+1.)
944 //
945 // The original Floyd&Rivest arranged the array into two intervals,
946 // one less than "T", one greater than (or equal to) "T".
947
948 int *idx = this->PartitionSubArray(L, R, K, dim, p1, p2);
949
950 I = idx[0];
951 J = idx[1];
952
953 if (K >= J)
954 {
955 L = J;
956 }
957 else if (K >= I)
958 {
959 L = R; // partitioning is done, K is in the interval of T's
960 }
961 else
962 {
963 R = I-1;
964 }
965 }
966 }
Select(int dim,int L,int R)967 int vtkPKdTree::Select(int dim, int L, int R)
968 {
969 int K = ((R + L) / 2) + 1;
970
971 this->_select(L, R, K, dim);
972
973 if (K == L) return K;
974
975 // The global array is now re-ordered, partitioned around X[K].
976 // (In particular, for all i, i<K, X[i] <= X[K] and for all i,
977 // i > K, X[i] >= X[K].)
978 // However the value at X[K] may occur more than once, and by
979 // construction of the reordered array, there is a J <= K such that
980 // for all i < J, X[i] < X[K] and for all J <= i < K X[i] = X[K].
981 //
982 // We want to roll K back to this value J, so that all points are
983 // unambiguously assigned to one region or the other.
984
985 int hasK = this->WhoHas(K);
986 int hasKrank = this->SubGroup->getLocalRank(hasK);
987
988 int hasKleft = this->WhoHas(K-1);
989 int hasKleftrank = this->SubGroup->getLocalRank(hasKleft);
990
991 float Kval;
992 float Kleftval;
993 float *pt;
994
995 if (hasK == this->MyId)
996 {
997 pt = this->GetLocalVal(K) + dim;
998 Kval = *pt;
999 }
1000
1001 this->SubGroup->Broadcast(&Kval, 1, hasKrank);
1002
1003 if (hasKleft == this->MyId)
1004 {
1005 pt = this->GetLocalVal(K-1) + dim;
1006 Kleftval = *pt;
1007 }
1008
1009 this->SubGroup->Broadcast(&Kleftval, 1, hasKleftrank);
1010
1011 if (Kleftval != Kval) return K;
1012
1013 int firstKval = this->TotalNumCells; // greater than any valid index
1014
1015 if ((this->MyId <= hasKleft) && (this->NumCells[this->MyId] > 0))
1016 {
1017 int start = this->EndVal[this->MyId];
1018 if (start > K-1) start = K-1;
1019
1020 pt = this->GetLocalVal(start) + dim;
1021
1022 if (*pt == Kval)
1023 {
1024 firstKval = start;
1025
1026 int finish = this->StartVal[this->MyId];
1027
1028 for (int idx=start-1; idx >= finish; idx--)
1029 {
1030 pt -= 3;
1031 if (*pt < Kval) break;
1032
1033 firstKval--;
1034 }
1035 }
1036 }
1037
1038 int newK;
1039
1040 this->SubGroup->ReduceMin(&firstKval, &newK, 1, hasKrank);
1041 this->SubGroup->Broadcast(&newK, 1, hasKrank);
1042
1043 return newK;
1044 }
1045
_whoHas(int L,int R,int pos)1046 int vtkPKdTree::_whoHas(int L, int R, int pos)
1047 {
1048 if (L == R)
1049 {
1050 return L;
1051 }
1052
1053 int M = (L + R) >> 1;
1054
1055 if ( pos < this->StartVal[M])
1056 {
1057 return _whoHas(L, M-1, pos);
1058 }
1059 else if (pos < this->StartVal[M+1])
1060 {
1061 return M;
1062 }
1063 else
1064 {
1065 return _whoHas(M+1, R, pos);
1066 }
1067 }
WhoHas(int pos)1068 int vtkPKdTree::WhoHas(int pos)
1069 {
1070 if ( (pos < 0) || (pos >= this->TotalNumCells))
1071 {
1072 return -1;
1073 }
1074 return _whoHas(0, this->NumProcesses-1, pos);
1075 }
GetLocalVal(int pos)1076 float *vtkPKdTree::GetLocalVal(int pos)
1077 {
1078 if ( (pos < this->StartVal[this->MyId]) || (pos > this->EndVal[this->MyId]))
1079 {
1080 return NULL;
1081 }
1082 int localPos = pos - this->StartVal[this->MyId];
1083
1084 return this->CurrentPtArray + (3*localPos);
1085 }
GetLocalValNext(int pos)1086 float *vtkPKdTree::GetLocalValNext(int pos)
1087 {
1088 if ( (pos < this->StartVal[this->MyId]) || (pos > this->EndVal[this->MyId]))
1089 {
1090 return NULL;
1091 }
1092 int localPos = pos - this->StartVal[this->MyId];
1093
1094 return this->NextPtArray + (3*localPos);
1095 }
SetLocalVal(int pos,float * val)1096 void vtkPKdTree::SetLocalVal(int pos, float *val)
1097 {
1098 if ( (pos < this->StartVal[this->MyId]) || (pos > this->EndVal[this->MyId]))
1099 {
1100 VTKERROR("SetLocalVal - bad index");
1101 return;
1102 }
1103
1104 int localOffset = (pos - this->StartVal[this->MyId]) * 3;
1105
1106 this->CurrentPtArray[localOffset] = val[0];
1107 this->CurrentPtArray[localOffset+1] = val[1];
1108 this->CurrentPtArray[localOffset+2] = val[2];
1109
1110 return;
1111 }
ExchangeLocalVals(int pos1,int pos2)1112 void vtkPKdTree::ExchangeLocalVals(int pos1, int pos2)
1113 {
1114 float temp[3];
1115
1116 float *pt1 = this->GetLocalVal(pos1);
1117 float *pt2 = this->GetLocalVal(pos2);
1118
1119 if (!pt1 || !pt2)
1120 {
1121 VTKERROR("ExchangeLocalVal - bad index");
1122 return;
1123 }
1124
1125 temp[0] = pt1[0];
1126 temp[1] = pt1[1];
1127 temp[2] = pt1[2];
1128
1129 pt1[0] = pt2[0];
1130 pt1[1] = pt2[1];
1131 pt1[2] = pt2[2];
1132
1133 pt2[0] = temp[0];
1134 pt2[1] = temp[1];
1135 pt2[2] = temp[2];
1136
1137 return;
1138 }
1139
DoTransfer(int from,int to,int fromIndex,int toIndex,int count)1140 void vtkPKdTree::DoTransfer(int from, int to, int fromIndex, int toIndex, int count)
1141 {
1142 float *fromPt, *toPt;
1143
1144 vtkCommunicator *comm = this->Controller->GetCommunicator();
1145
1146 int nitems = count * 3;
1147
1148 int me = this->MyId;
1149
1150 int tag = this->SubGroup->tag;
1151
1152 if ( (from==me) && (to==me))
1153 {
1154 fromPt = this->GetLocalVal(fromIndex);
1155 toPt = this->GetLocalValNext(toIndex);
1156
1157 memcpy(toPt, fromPt, nitems * sizeof(float));
1158 }
1159 else if (from == me)
1160 {
1161 fromPt = this->GetLocalVal(fromIndex);
1162
1163 comm->Send(fromPt, nitems, to, tag);
1164 }
1165 else if (to == me)
1166 {
1167 toPt = this->GetLocalValNext(toIndex);
1168
1169 comm->Receive(toPt, nitems, from, tag);
1170 }
1171 }
1172
1173 // Partition global array into three intervals, the first all values < T,
1174 // the second all values = T, the third all values > T. Return two
1175 // global indices: The index to the beginning of the second interval, and
1176 // the index to the beginning of the third interval. "T" is the value
1177 // at array index K.
1178 //
1179 // If there is no third interval, the second index returned will be R+1.
1180
PartitionSubArray(int L,int R,int K,int dim,int p1,int p2)1181 int *vtkPKdTree::PartitionSubArray(int L, int R, int K, int dim, int p1, int p2)
1182 {
1183 int rootrank = this->SubGroup->getLocalRank(p1);
1184
1185 int me = this->MyId;
1186
1187 if ( (me < p1) || (me > p2))
1188 {
1189 this->SubGroup->Broadcast(this->SelectBuffer, 2, rootrank);
1190 return this->SelectBuffer;
1191 }
1192
1193 if (p1 == p2)
1194 {
1195 int *idx = this->PartitionAboutMyValue(L, R, K, dim);
1196
1197 this->SubGroup->Broadcast(idx, 2, rootrank);
1198
1199 return idx;
1200 }
1201
1202 // Each process will rearrange their subarray myL-myR into a left region
1203 // of values less than X[K], a center region of values equal to X[K], and
1204 // a right region of values greater than X[K]. "I" will be the index
1205 // of the first value in the center region, or it will equal "J" if there
1206 // is no center region. "J" will be the index to the start of the
1207 // right region, or it will be R+1 if there is no right region.
1208
1209 int tag = this->SubGroup->tag;
1210
1211 vtkSubGroup *sg = vtkSubGroup::New();
1212 sg->Initialize(p1, p2, me, tag, this->Controller->GetCommunicator());
1213
1214 int hasK = this->WhoHas(K);
1215
1216 int Krank = sg->getLocalRank(hasK);
1217
1218 int myL = this->StartVal[me];
1219 int myR = this->EndVal[me];
1220
1221 if (myL < L) myL = L;
1222 if (myR > R) myR = R;
1223
1224 // Get Kth element
1225
1226 float T;
1227
1228 if (hasK == me)
1229 {
1230 T = this->GetLocalVal(K)[dim];
1231 }
1232
1233 sg->Broadcast(&T, 1, Krank);
1234
1235 int *idx; // dividing points in rearranged sub array
1236
1237 if (hasK == me)
1238 {
1239 idx = this->PartitionAboutMyValue(myL, myR, K, dim);
1240 }
1241 else
1242 {
1243 idx = this->PartitionAboutOtherValue(myL, myR, T, dim);
1244 }
1245
1246 // Copy these right away. Implementation uses SelectBuffer
1247 // which is about to be overwritten.
1248
1249 int I = idx[0];
1250 int J = idx[1];
1251
1252 // Now the ugly part. The processes redistribute the array so that
1253 // globally the interval [L:R] is partitioned into an interval of values
1254 // less than T, and interval of values equal to T, and an interval of
1255 // values greater than T.
1256
1257 int nprocs = p2 - p1 + 1;
1258
1259 int *buf = this->SelectBuffer;
1260
1261 int *left = buf; buf += nprocs; // global index of my leftmost
1262 int *right = buf; buf += nprocs; // global index of my rightmost
1263 int *Ival = buf; buf += nprocs; // global index of my first val = T
1264 int *Jval = buf; buf += nprocs; // global index of my first val > T
1265
1266 int *leftArray = buf; buf += nprocs; // number of my vals < T
1267 int *leftUsed = buf; buf += nprocs; // how many scheduled to be sent so far
1268
1269 int *centerArray = buf; buf += nprocs; // number of my vals = T
1270 int *centerUsed = buf; buf += nprocs; // how many scheduled to be sent so far
1271
1272 int *rightArray = buf; buf += nprocs; // number of my vals > T
1273 int *rightUsed = buf; buf += nprocs; // how many scheduled to be sent so far
1274
1275 rootrank = sg->getLocalRank(p1);
1276
1277 sg->Gather(&myL, left, 1, rootrank);
1278 sg->Broadcast(left, nprocs, rootrank);
1279
1280 sg->Gather(&myR, right, 1, rootrank);
1281 sg->Broadcast(right, nprocs, rootrank);
1282
1283 sg->Gather(&I, Ival, 1, rootrank);
1284 sg->Broadcast(Ival, nprocs, rootrank);
1285
1286 sg->Gather(&J, Jval, 1, rootrank);
1287 sg->Broadcast(Jval, nprocs, rootrank);
1288
1289 sg->Delete();
1290
1291 int leftRemaining = 0;
1292 int centerRemaining = 0;
1293
1294 int p, sndr, recvr;
1295
1296 for (p = 0; p < nprocs; p++)
1297 {
1298 leftArray[p] = Ival[p] - left[p];
1299 centerArray[p] = Jval[p] - Ival[p];
1300 rightArray[p] = right[p] - Jval[p] + 1;
1301
1302 leftRemaining += leftArray[p];
1303 centerRemaining += centerArray[p];
1304
1305 leftUsed[p] = 0;
1306 centerUsed[p] = 0;
1307 rightUsed[p] = 0;
1308 }
1309
1310 int FirstCenter = left[0] + leftRemaining;
1311 int FirstRight = FirstCenter + centerRemaining;
1312
1313 int nextLeftProc = 0;
1314 int nextCenterProc = 0;
1315 int nextRightProc = 0;
1316
1317 int need, have, take;
1318
1319 if ( (myL > this->StartVal[me]) || (myR < this->EndVal[me]))
1320 {
1321 memcpy(this->NextPtArray,
1322 this->CurrentPtArray,
1323 this->PtArraySize * sizeof(float));
1324 }
1325
1326 for (recvr = 0; recvr < nprocs; recvr++)
1327 {
1328 need = leftArray[recvr] + centerArray[recvr] + rightArray[recvr];
1329 have = 0;
1330
1331 if (leftRemaining >= 0)
1332 {
1333 for (sndr = nextLeftProc; sndr < nprocs; sndr++)
1334 {
1335 take = leftArray[sndr] - leftUsed[sndr];
1336
1337 if (take == 0) continue;
1338
1339 take = (take > need) ? need : take;
1340
1341 this->DoTransfer(sndr + p1, recvr + p1,
1342 left[sndr] + leftUsed[sndr], left[recvr] + have, take);
1343
1344 have += take;
1345 need -= take;
1346 leftRemaining -= take;
1347
1348 leftUsed[sndr] += take;
1349
1350 if (need == 0) break;
1351 }
1352
1353 if (leftUsed[sndr] == leftArray[sndr])
1354 {
1355 nextLeftProc = sndr+1;
1356 }
1357 else
1358 {
1359 nextLeftProc = sndr;
1360 }
1361 }
1362
1363 if (need == 0) continue;
1364
1365 if (centerRemaining >= 0)
1366 {
1367 for (sndr = nextCenterProc; sndr < nprocs; sndr++)
1368 {
1369 take = centerArray[sndr] - centerUsed[sndr];
1370
1371 if (take == 0) continue;
1372
1373 take = (take > need) ? need : take;
1374
1375 // Just copy the values, since we know what they are
1376 this->DoTransfer(sndr + p1, recvr + p1,
1377 left[sndr] + leftArray[sndr] + centerUsed[sndr],
1378 left[recvr] + have, take);
1379
1380 have += take;
1381 need -= take;
1382 centerRemaining -= take;
1383
1384 centerUsed[sndr] += take;
1385
1386 if (need == 0) break;
1387 }
1388
1389 if (centerUsed[sndr] == centerArray[sndr])
1390 {
1391 nextCenterProc = sndr+1;
1392 }
1393 else
1394 {
1395 nextCenterProc = sndr;
1396 }
1397 }
1398
1399 if (need == 0) continue;
1400
1401 for (sndr = nextRightProc; sndr < nprocs; sndr++)
1402 {
1403 take = rightArray[sndr] - rightUsed[sndr];
1404
1405 if (take == 0) continue;
1406
1407 take = (take > need) ? need : take;
1408
1409 this->DoTransfer(
1410 sndr + p1, recvr + p1,
1411 left[sndr] + leftArray[sndr] + centerArray[sndr] + rightUsed[sndr],
1412 left[recvr] + have, take);
1413
1414 have += take;
1415 need -= take;
1416
1417 rightUsed[sndr] += take;
1418
1419 if (need == 0) break;
1420 }
1421
1422 if (rightUsed[sndr] == rightArray[sndr])
1423 {
1424 nextRightProc = sndr+1;
1425 }
1426 else
1427 {
1428 nextRightProc = sndr;
1429 }
1430 }
1431
1432 this->SwitchDoubleBuffer();
1433
1434 this->SelectBuffer[0] = FirstCenter;
1435 this->SelectBuffer[1] = FirstRight;
1436
1437 rootrank = this->SubGroup->getLocalRank(p1);
1438
1439 this->SubGroup->Broadcast(this->SelectBuffer, 2, rootrank);
1440
1441 return this->SelectBuffer;
1442 }
1443
1444 // This routine partitions the array from element L through element
1445 // R into three segments. This first contains values less than T, the
1446 // next contains values equal to T, the last has values greater than T.
1447 //
1448 // This routine returns two values. The first is the index of the
1449 // first value equal to T, the second is the index of the first value
1450 // greater than T. If there is no value equal to T, the first value
1451 // will equal the second value. If there is no value greater than T,
1452 // the second value returned will be R+1.
1453 //
1454 // This function is different than PartitionAboutMyValue, because in
1455 // that functin we know that "T" appears in the array. In this
1456 // function, "T" may or may not appear in the array.
1457
PartitionAboutOtherValue(int L,int R,float T,int dim)1458 int *vtkPKdTree::PartitionAboutOtherValue(int L, int R, float T, int dim)
1459 {
1460 float *Ipt, *Jpt, Lval, Rval;
1461 int *vals = this->SelectBuffer;
1462 int numTValues = 0;
1463 int numGreater = 0;
1464 int numLess = 0;
1465 int totalVals = R - L + 1;
1466
1467 if (totalVals == 0)
1468 {
1469 // Special case: no values.
1470 vals[0] = vals[1] = L;
1471 return vals;
1472 }
1473
1474 Ipt = this->GetLocalVal(L) + dim;
1475 Lval = *Ipt;
1476
1477 if (Lval == T)
1478 {
1479 numTValues++;
1480 }
1481 else if (Lval > T)
1482 {
1483 numGreater++;
1484 }
1485 else
1486 {
1487 numLess++;
1488 }
1489
1490 Jpt = this->GetLocalVal(R) + dim;
1491 Rval = *Jpt;
1492
1493 if (Rval == T) numTValues++;
1494 else if (Rval > T) numGreater++;
1495 else numLess++;
1496
1497 int I = L;
1498 int J = R;
1499
1500 if ((Lval >= T) && (Rval >= T))
1501 {
1502 while (--J > I)
1503 {
1504 Jpt -= 3;
1505 if (*Jpt < T)
1506 {
1507 break;
1508 }
1509 if (*Jpt == T)
1510 {
1511 numTValues++;
1512 }
1513 else
1514 {
1515 numGreater++;
1516 }
1517 }
1518 }
1519 else if ((Lval < T) && (Rval < T))
1520 {
1521 Ipt = this->GetLocalVal(I) + dim;
1522
1523 while (++I < J)
1524 {
1525 Ipt += 3;
1526 if (*Ipt >= T)
1527 {
1528 if (*Ipt == T)
1529 {
1530 numTValues++;
1531 }
1532 break;
1533 }
1534 numLess++;
1535 }
1536 }
1537 else if ((Lval < T) && (Rval >= T))
1538 {
1539 this->ExchangeLocalVals(I, J);
1540 }
1541 else if ((Lval >= T) && (Rval < T))
1542 {
1543 // first loop will fix this
1544 }
1545
1546 if (numLess == totalVals)
1547 {
1548 vals[0] = vals[1] = R+1; // special case - all less than T
1549 return vals;
1550 }
1551 else if (numTValues == totalVals)
1552 {
1553 vals[0] = L; // special case - all equal to T
1554 vals[1] = R+1;
1555 return vals;
1556 }
1557 else if (numGreater == totalVals)
1558 {
1559 vals[0] = vals[1] = L; // special case - all greater than T
1560 return vals;
1561 }
1562
1563 while (I < J)
1564 {
1565 // By design, I < J and value at I is >= T, and value
1566 // at J is < T, hence the exchange.
1567
1568 this->ExchangeLocalVals(I, J);
1569
1570 while (++I < J)
1571 {
1572 Ipt += 3;
1573 if (*Ipt >= T)
1574 {
1575 if (*Ipt == T)
1576 {
1577 numTValues++;
1578 }
1579 break;
1580 }
1581 }
1582 if (I == J) break;
1583
1584 while (--J > I)
1585 {
1586 Jpt -= 3;
1587 if (*Jpt < T)
1588 {
1589 break;
1590 }
1591 if (*Jpt == T)
1592 {
1593 numTValues++;
1594 }
1595 }
1596 }
1597
1598 // I and J are at the first value that is >= T.
1599
1600 if (numTValues == 0)
1601 {
1602 vals[0] = I;
1603 vals[1] = I;
1604 return vals;
1605 }
1606
1607 // Move all T's to the center interval
1608
1609 vals[0] = I; // the first T will be here when we're done
1610
1611 Ipt = this->GetLocalVal(I) + dim;
1612 I = I-1;
1613 Ipt -= 3;
1614
1615 J = R+1;
1616 Jpt = this->GetLocalVal(R) + dim;
1617 Jpt += 3;
1618
1619 while (I < J)
1620 {
1621 while (++I < J)
1622 {
1623 Ipt += 3;
1624 if (*Ipt != T)
1625 {
1626 break;
1627 }
1628 }
1629 if (I == J)
1630 {
1631 break;
1632 }
1633
1634 while (--J > I)
1635 {
1636 Jpt -= 3;
1637 if (*Jpt == T)
1638 {
1639 break;
1640 }
1641 }
1642
1643 if (I < J)
1644 {
1645 this->ExchangeLocalVals(I, J);
1646 }
1647 }
1648
1649 // Now I and J are at the first value that is > T, and the T's are
1650 // to the left.
1651
1652 vals[1] = I; // the first > T
1653
1654 return vals;
1655 }
1656
1657 // This routine partitions the array from element L through element
1658 // R into three segments. This first contains values less than T, the
1659 // next contains values equal to T, the last has values greater than T.
1660 // T is the element at K, where L <= K <= R.
1661 //
1662 // This routine returns two integers. The first is the index of the
1663 // first value equal to T, the second is the index of the first value
1664 // greater than T. If there is no value greater than T, the second
1665 // value returned will be R+1.
1666
PartitionAboutMyValue(int L,int R,int K,int dim)1667 int *vtkPKdTree::PartitionAboutMyValue(int L, int R, int K, int dim)
1668 {
1669 float *Ipt, *Jpt;
1670 float T;
1671 int I, J;
1672 int manyTValues = 0;
1673 int *vals = this->SelectBuffer;
1674
1675 // Set up so after first exchange in the loop, we have either
1676 // X[L] = T and X[R] >= T
1677 // or
1678 // X[L] < T and X[R] = T
1679 //
1680
1681 float *pt = this->GetLocalVal(K);
1682
1683 T = pt[dim];
1684
1685 this->ExchangeLocalVals(L, K);
1686
1687 pt = this->GetLocalVal(R);
1688
1689 if (pt[dim] >= T)
1690 {
1691 if (pt[dim] == T)
1692 {
1693 manyTValues = 1;
1694 }
1695 else
1696 {
1697 this->ExchangeLocalVals(R, L);
1698 }
1699 }
1700
1701 I = L;
1702 J = R;
1703
1704 Ipt = this->GetLocalVal(I) + dim;
1705 Jpt = this->GetLocalVal(J) + dim;
1706
1707 while (I < J)
1708 {
1709 this->ExchangeLocalVals(I, J);
1710
1711 while (--J > I)
1712 {
1713 Jpt -= 3;
1714 if (*Jpt < T)
1715 {
1716 break;
1717 }
1718 if (!manyTValues && (J > L) && (*Jpt == T))
1719 {
1720 manyTValues = 1;
1721 }
1722 }
1723
1724 if (I == J)
1725 {
1726 break;
1727 }
1728
1729 while (++I < J)
1730 {
1731 Ipt += 3;
1732
1733 if (*Ipt >= T)
1734 {
1735 if (!manyTValues && (*Ipt == T))
1736 {
1737 manyTValues = 1;
1738 }
1739 break;
1740 }
1741 }
1742 }
1743
1744 // I and J are at the rightmost value < T ( or at L if all values
1745 // are >= T)
1746
1747 pt = this->GetLocalVal(L);
1748
1749 float Lval = pt[dim];
1750
1751 if (Lval == T)
1752 {
1753 this->ExchangeLocalVals(L, J);
1754 }
1755 else
1756 {
1757 this->ExchangeLocalVals(++J, R);
1758 }
1759
1760 // Now J is at the leftmost value >= T. (It is at a T value.)
1761
1762 vals[0] = J;
1763 vals[1] = J+1;
1764
1765 // Arrange array so all values equal to T are together
1766
1767 if (manyTValues)
1768 {
1769 I = J;
1770 Ipt = this->GetLocalVal(I) + dim;
1771
1772 J = R+1;
1773 Jpt = this->GetLocalVal(R) + dim;
1774 Jpt += 3;
1775
1776 while (I < J)
1777 {
1778 while (++I < J)
1779 {
1780 Ipt += 3;
1781 if (*Ipt != T)
1782 {
1783 break;
1784 }
1785 }
1786 if (I == J)
1787 {
1788 break;
1789 }
1790
1791 while (--J > I)
1792 {
1793 Jpt -= 3;
1794 if (*Jpt == T)
1795 {
1796 break;
1797 }
1798 }
1799
1800 if (I < J)
1801 {
1802 this->ExchangeLocalVals(I, J);
1803 }
1804 }
1805 // I and J are at the first value that is > T
1806
1807 vals[1] = I;
1808 }
1809
1810 return vals;
1811 }
1812
1813 //--------------------------------------------------------------------
1814 // Compute the bounds for the data in a region
1815 //--------------------------------------------------------------------
1816
GetLocalMinMax(int L,int R,int me,float * min,float * max)1817 void vtkPKdTree::GetLocalMinMax(int L, int R, int me,
1818 float *min, float *max)
1819 {
1820 int i, d;
1821 int from = this->StartVal[me];
1822 int to = this->EndVal[me];
1823
1824 if (L > from)
1825 {
1826 from = L;
1827 }
1828 if (R < to)
1829 {
1830 to = R;
1831 }
1832
1833 if (from <= to)
1834 {
1835 from -= this->StartVal[me];
1836 to -= this->StartVal[me];
1837
1838 float *val = this->CurrentPtArray + from*3;
1839
1840 for (d=0; d<3; d++)
1841 {
1842 min[d] = max[d] = val[d];
1843 }
1844
1845 for (i= from+1; i<=to; i++)
1846 {
1847 val += 3;
1848
1849 for (d=0; d<3; d++)
1850 {
1851 if (val[d] < min[d])
1852 {
1853 min[d] = val[d];
1854 }
1855 else if (val[d] > max[d])
1856 {
1857 max[d] = val[d];
1858 }
1859 }
1860 }
1861 }
1862 else
1863 {
1864 // this guy has none of the data, but still must participate
1865 // in ReduceMax and ReduceMin
1866
1867 double *regionMin = this->Top->GetMinBounds();
1868 double *regionMax = this->Top->GetMaxBounds();
1869
1870 for (d=0; d<3; d++)
1871 {
1872 min[d] = (float)regionMax[d];
1873 max[d] = (float)regionMin[d];
1874 }
1875 }
1876 }
DataBounds(int L,int K,int R)1877 float *vtkPKdTree::DataBounds(int L, int K, int R)
1878 {
1879 float localMinLeft[3]; // Left region is L through K-1
1880 float localMaxLeft[3];
1881 float globalMinLeft[3];
1882 float globalMaxLeft[3];
1883 float localMinRight[3]; // Right region is K through R
1884 float localMaxRight[3];
1885 float globalMinRight[3];
1886 float globalMaxRight[3];
1887
1888 float *globalBounds = new float [12];
1889
1890 int fail = (globalBounds == NULL);
1891
1892 if (this->AllCheckForFailure(fail, "DataBounds", "memory allocation"))
1893 {
1894 delete [] globalBounds;
1895 return NULL;
1896 }
1897
1898 this->GetLocalMinMax(L, K-1, this->MyId, localMinLeft, localMaxLeft);
1899
1900 this->GetLocalMinMax(K, R, this->MyId, localMinRight, localMaxRight);
1901
1902 this->SubGroup->ReduceMin(localMinLeft, globalMinLeft, 3, 0);
1903 this->SubGroup->Broadcast(globalMinLeft, 3, 0);
1904
1905 this->SubGroup->ReduceMax(localMaxLeft, globalMaxLeft, 3, 0);
1906 this->SubGroup->Broadcast(globalMaxLeft, 3, 0);
1907
1908 this->SubGroup->ReduceMin(localMinRight, globalMinRight, 3, 0);
1909 this->SubGroup->Broadcast(globalMinRight, 3, 0);
1910
1911 this->SubGroup->ReduceMax(localMaxRight, globalMaxRight, 3, 0);
1912 this->SubGroup->Broadcast(globalMaxRight, 3, 0);
1913
1914 float *left = globalBounds;
1915 float *right = globalBounds + 6;
1916
1917 MinMaxToBounds(left, globalMinLeft, globalMaxLeft);
1918
1919 MinMaxToBounds(right, globalMinRight, globalMaxRight);
1920
1921 return globalBounds;
1922 }
1923
1924 //--------------------------------------------------------------------
1925 // Complete the tree - Different nodes of tree were computed by different
1926 // processors. Now put it together.
1927 //--------------------------------------------------------------------
1928
CompleteTree()1929 int vtkPKdTree::CompleteTree()
1930 {
1931 // calculate depth of entire tree
1932
1933 int depth;
1934 int myDepth = vtkPKdTree::ComputeDepth(this->Top);
1935
1936 this->SubGroup->ReduceMax(&myDepth, &depth, 1, 0);
1937 this->SubGroup->Broadcast(&depth, 1, 0);
1938
1939 // fill out nodes of tree
1940
1941 int fail = vtkPKdTree::FillOutTree(this->Top, depth);
1942
1943 if (this->AllCheckForFailure(fail, "CompleteTree", "memory allocation"))
1944 {
1945 return 1;
1946 }
1947
1948 // Processor 0 collects all the nodes of the k-d tree, and then
1949 // processes the tree to ensure region boundaries are
1950 // consistent. The completed tree is then broadcast.
1951
1952 int *buf = new int [this->NumProcesses];
1953
1954 fail = (buf == NULL);
1955
1956 if (this->AllCheckForFailure(fail, "CompleteTree", "memory allocation"))
1957 {
1958 if (buf) delete [] buf;
1959 return 1;
1960 }
1961
1962 #ifdef YIELDS_INCONSISTENT_REGION_BOUNDARIES
1963
1964 this->RetrieveData(this->Top, buf);
1965
1966 #else
1967
1968 this->ReduceData(this->Top, buf);
1969
1970 if (this->MyId == 0)
1971 {
1972 CheckFixRegionBoundaries(this->Top);
1973 }
1974
1975 this->BroadcastData(this->Top);
1976 #endif
1977
1978 delete [] buf;
1979
1980 return 0;
1981 }
1982
PackData(vtkKdNode * kd,double * data)1983 void vtkPKdTree::PackData(vtkKdNode *kd, double *data)
1984 {
1985 int i, v;
1986
1987 data[0] = (double)kd->GetDim();
1988 data[1] = (double)kd->GetLeft()->GetNumberOfPoints();
1989 data[2] = (double)kd->GetRight()->GetNumberOfPoints();
1990
1991 double *lmin = kd->GetLeft()->GetMinBounds();
1992 double *lmax = kd->GetLeft()->GetMaxBounds();
1993 double *lminData = kd->GetLeft()->GetMinDataBounds();
1994 double *lmaxData = kd->GetLeft()->GetMaxDataBounds();
1995 double *rmin = kd->GetRight()->GetMinBounds();
1996 double *rmax = kd->GetRight()->GetMaxBounds();
1997 double *rminData = kd->GetRight()->GetMinDataBounds();
1998 double *rmaxData = kd->GetRight()->GetMaxDataBounds();
1999
2000 v = 3;
2001 for (i=0; i<3; i++)
2002 {
2003 data[v++] = lmin[i];
2004 data[v++] = lmax[i];
2005 data[v++] = lminData[i];
2006 data[v++] = lmaxData[i];
2007 data[v++] = rmin[i];
2008 data[v++] = rmax[i];
2009 data[v++] = rminData[i];
2010 data[v++] = rmaxData[i];
2011 }
2012 }
UnpackData(vtkKdNode * kd,double * data)2013 void vtkPKdTree::UnpackData(vtkKdNode *kd, double *data)
2014 {
2015 int i, v;
2016
2017 kd->SetDim((int)data[0]);
2018 kd->GetLeft()->SetNumberOfPoints((int)data[1]);
2019 kd->GetRight()->SetNumberOfPoints((int)data[2]);
2020
2021 double lmin[3], rmin[3], lmax[3], rmax[3];
2022 double lminData[3], rminData[3], lmaxData[3], rmaxData[3];
2023
2024 v = 3;
2025 for (i=0; i<3; i++)
2026 {
2027 lmin[i] = data[v++];
2028 lmax[i] = data[v++];
2029 lminData[i] = data[v++];
2030 lmaxData[i] = data[v++];
2031 rmin[i] = data[v++];
2032 rmax[i] = data[v++];
2033 rminData[i] = data[v++];
2034 rmaxData[i] = data[v++];
2035 }
2036
2037 kd->GetLeft()->SetBounds(lmin[0], lmax[0],
2038 lmin[1], lmax[1],
2039 lmin[2], lmax[2]);
2040 kd->GetLeft()->SetDataBounds(lminData[0], lmaxData[0],
2041 lminData[1], lmaxData[1],
2042 lminData[2], lmaxData[2]);
2043
2044 kd->GetRight()->SetBounds(rmin[0], rmax[0],
2045 rmin[1], rmax[1],
2046 rmin[2], rmax[2]);
2047 kd->GetRight()->SetDataBounds(rminData[0], rmaxData[0],
2048 rminData[1], rmaxData[1],
2049 rminData[2], rmaxData[2]);
2050 }
ReduceData(vtkKdNode * kd,int * sources)2051 void vtkPKdTree::ReduceData(vtkKdNode *kd, int *sources)
2052 {
2053 int i;
2054 double data[27];
2055 vtkCommunicator *comm = this->Controller->GetCommunicator();
2056
2057 if (kd->GetLeft() == NULL) return;
2058
2059 int ihave = (kd->GetDim() < 3);
2060
2061 this->SubGroup->Gather(&ihave, sources, 1, 0);
2062 this->SubGroup->Broadcast(sources, this->NumProcesses, 0);
2063
2064 // a contiguous group of process IDs built this node, the first
2065 // in the group sends it to node 0 if node 0 doesn't have it
2066
2067 if (sources[0] == 0)
2068 {
2069 int root = -1;
2070
2071 for (i=1; i<this->NumProcesses; i++)
2072 {
2073 if (sources[i])
2074 {
2075 root = i;
2076 break;
2077 }
2078 }
2079 if (root == -1)
2080 {
2081
2082 // Normally BuildLocator will create a complete tree, but
2083 // it may refuse to divide a region if all the data is at
2084 // the same point along the axis it wishes to divide. In
2085 // that case, this region was not divided, so just return.
2086
2087 vtkKdTree::DeleteAllDescendants(kd);
2088
2089 return;
2090 }
2091
2092 if (root == this->MyId)
2093 {
2094 vtkPKdTree::PackData(kd, data);
2095
2096 comm->Send(data, 27, 0, 0x1111);
2097 }
2098 else if (0 == this->MyId)
2099 {
2100 comm->Receive(data, 27, root, 0x1111);
2101
2102 vtkPKdTree::UnpackData(kd, data);
2103 }
2104 }
2105
2106 this->ReduceData(kd->GetLeft(), sources);
2107
2108 this->ReduceData(kd->GetRight(), sources);
2109
2110 return;
2111 }
BroadcastData(vtkKdNode * kd)2112 void vtkPKdTree::BroadcastData(vtkKdNode *kd)
2113 {
2114 double data[27];
2115
2116 if (kd->GetLeft() == NULL) return;
2117
2118 if (0 == this->MyId)
2119 {
2120 vtkPKdTree::PackData(kd, data);
2121 }
2122
2123 this->SubGroup->Broadcast(data, 27, 0);
2124
2125 if (this->MyId > 0)
2126 {
2127 vtkPKdTree::UnpackData(kd, data);
2128 }
2129
2130 this->BroadcastData(kd->GetLeft());
2131
2132 this->BroadcastData(kd->GetRight());
2133
2134 return;
2135 }
CheckFixRegionBoundaries(vtkKdNode * tree)2136 void vtkPKdTree::CheckFixRegionBoundaries(vtkKdNode *tree)
2137 {
2138 if (tree->GetLeft() == NULL) return;
2139
2140 int nextDim = tree->GetDim();
2141
2142 vtkKdNode *left = tree->GetLeft();
2143 vtkKdNode *right = tree->GetRight();
2144
2145 double *min = tree->GetMinBounds();
2146 double *max = tree->GetMaxBounds();
2147 double *lmin = left->GetMinBounds();
2148 double *lmax = left->GetMaxBounds();
2149 double *rmin = right->GetMinBounds();
2150 double *rmax = right->GetMaxBounds();
2151
2152 for (int dim=0; dim<3; dim++)
2153 {
2154 if ((lmin[dim] - min[dim]) != 0.0) lmin[dim] = min[dim];
2155 if ((rmax[dim] - max[dim]) != 0.0) rmax[dim] = max[dim];
2156
2157 if (dim != nextDim) // the dimension I did *not* divide along
2158 {
2159 if ((lmax[dim] - max[dim]) != 0.0) lmax[dim] = max[dim];
2160 if ((rmin[dim] - min[dim]) != 0.0) rmin[dim] = min[dim];
2161 }
2162 else
2163 {
2164 if ((lmax[dim] - rmin[dim]) != 0.0) lmax[dim] = rmin[dim];
2165 }
2166 }
2167
2168 CheckFixRegionBoundaries(left);
2169 CheckFixRegionBoundaries(right);
2170 }
2171 #ifdef YIELDS_INCONSISTENT_REGION_BOUNDARIES
2172
RetrieveData(vtkKdNode * kd,int * sources)2173 void vtkPKdTree::RetrieveData(vtkKdNode *kd, int *sources)
2174 {
2175 int i;
2176 double data[27];
2177
2178 if (kd->GetLeft() == NULL) return;
2179
2180 int ihave = (kd->GetDim() < 3);
2181
2182 this->SubGroup->Gather(&ihave, sources, 1, 0);
2183 this->SubGroup->Broadcast(sources, this->NumProcesses, 0);
2184
2185 // a contiguous group of process IDs built this node, the first
2186 // in the group broadcasts the results to everyone else
2187
2188 int root = -1;
2189
2190 for (i=0; i<this->NumProcesses; i++)
2191 {
2192 if (sources[i])
2193 {
2194 root = i;
2195 break;
2196 }
2197 }
2198 if (root == -1)
2199 {
2200 // Normally BuildLocator will create a complete tree, but
2201 // it may refuse to divide a region if all the data is at
2202 // the same point along the axis it wishes to divide. In
2203 // that case, this region was not divided, so just return.
2204
2205 vtkKdTree::DeleteAllDescendants(kd);
2206
2207 return;
2208 }
2209
2210 if (root == this->MyId)
2211 {
2212 vtkPKdTree::PackData(kd, data);
2213 }
2214
2215 this->SubGroup->Broadcast(data, 27, root);
2216
2217 if (!ihave)
2218 {
2219 vtkPKdTree::UnpackData(kd, data);
2220 }
2221
2222 this->RetrieveData(kd->GetLeft(), sources);
2223
2224 this->RetrieveData(kd->GetRight(), sources);
2225
2226 return;
2227 }
2228 #endif
2229
FillOutTree(vtkKdNode * kd,int level)2230 int vtkPKdTree::FillOutTree(vtkKdNode *kd, int level)
2231 {
2232 if (level == 0) return 0;
2233
2234 if (kd->GetLeft() == NULL)
2235 {
2236 vtkKdNode *left = vtkKdNode::New();
2237
2238 if (!left) goto doneError2;
2239
2240 left->SetBounds(-1,-1,-1,-1,-1,-1);
2241 left->SetDataBounds(-1,-1,-1,-1,-1,-1);
2242 left->SetNumberOfPoints(-1);
2243
2244 vtkKdNode *right = vtkKdNode::New();
2245
2246 if (!right) goto doneError2;
2247
2248 right->SetBounds(-1,-1,-1,-1,-1,-1);
2249 right->SetDataBounds(-1,-1,-1,-1,-1,-1);
2250 right->SetNumberOfPoints(-1);
2251
2252 kd->AddChildNodes(left, right);
2253 }
2254
2255 if (vtkPKdTree::FillOutTree(kd->GetLeft(), level-1)) goto doneError2;
2256
2257 if (vtkPKdTree::FillOutTree(kd->GetRight(), level-1)) goto doneError2;
2258
2259 return 0;
2260
2261 doneError2:
2262
2263 return 1;
2264 }
2265
ComputeDepth(vtkKdNode * kd)2266 int vtkPKdTree::ComputeDepth(vtkKdNode *kd)
2267 {
2268 int leftDepth = 0;
2269 int rightDepth = 0;
2270
2271 if ((kd->GetLeft() == NULL) && (kd->GetRight() == NULL)) return 0;
2272
2273 if (kd->GetLeft())
2274 {
2275 leftDepth = vtkPKdTree::ComputeDepth(kd->GetLeft());
2276 }
2277 if (kd->GetRight())
2278 {
2279 rightDepth = vtkPKdTree::ComputeDepth(kd->GetRight());
2280 }
2281
2282 if (leftDepth > rightDepth) return leftDepth + 1;
2283 else return rightDepth + 1;
2284 }
2285
2286 //--------------------------------------------------------------------
2287 // lists, lists, lists
2288 //--------------------------------------------------------------------
2289
AllocateDoubleBuffer()2290 int vtkPKdTree::AllocateDoubleBuffer()
2291 {
2292 this->FreeDoubleBuffer();
2293
2294 this->PtArraySize = this->NumCells[this->MyId] * 3;
2295
2296 this->PtArray2 = new float [this->PtArraySize];
2297
2298 this->CurrentPtArray = this->PtArray;
2299 this->NextPtArray = this->PtArray2;
2300
2301 return (this->PtArray2 == NULL);
2302 }
SwitchDoubleBuffer()2303 void vtkPKdTree::SwitchDoubleBuffer()
2304 {
2305 float *temp = this->CurrentPtArray;
2306
2307 this->CurrentPtArray = this->NextPtArray;
2308 this->NextPtArray = temp;
2309 }
FreeDoubleBuffer()2310 void vtkPKdTree::FreeDoubleBuffer()
2311 {
2312 FreeList(this->PtArray2);
2313 this->CurrentPtArray = this->PtArray;
2314 this->NextPtArray = NULL;
2315 }
2316
AllocateSelectBuffer()2317 int vtkPKdTree::AllocateSelectBuffer()
2318 {
2319 this->FreeSelectBuffer();
2320
2321 this->SelectBuffer = new int [this->NumProcesses * 10];
2322
2323 return (this->SelectBuffer == NULL);
2324 }
FreeSelectBuffer()2325 void vtkPKdTree::FreeSelectBuffer()
2326 {
2327 delete [] this->SelectBuffer;
2328 this->SelectBuffer = NULL;
2329 }
2330
2331 #define FreeListOfLists(list, len) \
2332 { \
2333 int i; \
2334 if (list) \
2335 { \
2336 for (i=0; i<len; i++) \
2337 { \
2338 if (list[i]) delete [] list[i]; \
2339 } \
2340 delete [] list; \
2341 list = NULL; \
2342 } \
2343 }
2344
2345 #define MakeList(field, type, len) \
2346 { \
2347 if (len>0) \
2348 { \
2349 field = new type [len]; \
2350 if (field) \
2351 { \
2352 memset(field, 0, (len) * sizeof(type)); \
2353 } \
2354 } \
2355 }
2356
2357 // global index lists -----------------------------------------------
2358
InitializeGlobalIndexLists()2359 void vtkPKdTree::InitializeGlobalIndexLists()
2360 {
2361 this->StartVal = NULL;
2362 this->EndVal = NULL;
2363 this->NumCells = NULL;
2364 }
AllocateAndZeroGlobalIndexLists()2365 int vtkPKdTree::AllocateAndZeroGlobalIndexLists()
2366 {
2367 this->FreeGlobalIndexLists();
2368
2369 MakeList(this->StartVal, vtkIdType, this->NumProcesses);
2370 MakeList(this->EndVal, vtkIdType, this->NumProcesses);
2371 MakeList(this->NumCells, vtkIdType, this->NumProcesses);
2372
2373 int defined = ((this->StartVal != NULL) &&
2374 (this->EndVal != NULL) &&
2375 (this->NumCells != NULL));
2376
2377 if (!defined) this->FreeGlobalIndexLists();
2378
2379 return !defined;
2380 }
FreeGlobalIndexLists()2381 void vtkPKdTree::FreeGlobalIndexLists()
2382 {
2383 FreeList(StartVal);
2384 FreeList(EndVal);
2385 FreeList(NumCells);
2386 }
BuildGlobalIndexLists(vtkIdType numMyCells)2387 int vtkPKdTree::BuildGlobalIndexLists(vtkIdType numMyCells)
2388 {
2389 int fail = this->AllocateAndZeroGlobalIndexLists();
2390
2391 if (this->AllCheckForFailure(fail,
2392 "BuildGlobalIndexLists",
2393 "memory allocation"))
2394 {
2395 this->FreeGlobalIndexLists();
2396 return 1;
2397 }
2398
2399 this->SubGroup->Gather(&numMyCells, this->NumCells, 1, 0);
2400
2401 this->SubGroup->Broadcast(this->NumCells, this->NumProcesses, 0);
2402
2403 this->StartVal[0] = 0;
2404 this->EndVal[0] = this->NumCells[0] - 1;
2405
2406 this->TotalNumCells = this->NumCells[0];
2407
2408 for (int i=1; i<this->NumProcesses; i++)
2409 {
2410 this->StartVal[i] = this->EndVal[i-1] + 1;
2411 this->EndVal[i] = this->EndVal[i-1] + this->NumCells[i];
2412
2413 this->TotalNumCells += this->NumCells[i];
2414 }
2415
2416 return 0;
2417 }
2418
2419 // Region assignment lists ---------------------------------------------
2420
InitializeRegionAssignmentLists()2421 void vtkPKdTree::InitializeRegionAssignmentLists()
2422 {
2423 this->RegionAssignmentMap = NULL;
2424 this->RegionAssignmentMapLength = 0;
2425 this->ProcessAssignmentMap = NULL;
2426 this->NumRegionsAssigned = NULL;
2427 }
AllocateAndZeroRegionAssignmentLists()2428 int vtkPKdTree::AllocateAndZeroRegionAssignmentLists()
2429 {
2430 this->FreeRegionAssignmentLists();
2431
2432 this->RegionAssignmentMapLength = this->GetNumberOfRegions();
2433
2434 MakeList(this->RegionAssignmentMap, int , this->GetNumberOfRegions());
2435 MakeList(this->NumRegionsAssigned, int, this->NumProcesses);
2436
2437 MakeList(this->ProcessAssignmentMap, int *, this->NumProcesses);
2438
2439 int defined = ((this->RegionAssignmentMap != NULL) &&
2440 (this->ProcessAssignmentMap != NULL) &&
2441 (this->NumRegionsAssigned != NULL) );
2442
2443 if (!defined) this->FreeRegionAssignmentLists();
2444
2445 return !defined;
2446 }
FreeRegionAssignmentLists()2447 void vtkPKdTree::FreeRegionAssignmentLists()
2448 {
2449 FreeList(this->RegionAssignmentMap);
2450 FreeList(this->NumRegionsAssigned);
2451 FreeListOfLists(this->ProcessAssignmentMap, this->NumProcesses);
2452
2453 this->RegionAssignmentMapLength = 0;
2454 }
2455
2456 // Process data tables ------------------------------------------------
2457
InitializeProcessDataLists()2458 void vtkPKdTree::InitializeProcessDataLists()
2459 {
2460 this->DataLocationMap = NULL;
2461
2462 this->NumProcessesInRegion = NULL;
2463 this->ProcessList = NULL;
2464
2465 this->NumRegionsInProcess = NULL;
2466 this->RegionList = NULL;
2467
2468 this->CellCountList = NULL;
2469 }
2470
AllocateAndZeroProcessDataLists()2471 int vtkPKdTree::AllocateAndZeroProcessDataLists()
2472 {
2473 int nRegions = this->GetNumberOfRegions();
2474 int nProcesses = this->NumProcesses;
2475
2476 this->FreeProcessDataLists();
2477
2478 MakeList(this->DataLocationMap, char, nRegions * nProcesses);
2479
2480 if (this->DataLocationMap == NULL) goto doneError3;
2481
2482 MakeList(this->NumProcessesInRegion, int ,nRegions);
2483
2484 if (this->NumProcessesInRegion == NULL) goto doneError3;
2485
2486 MakeList(this->ProcessList, int * ,nRegions);
2487
2488 if (this->ProcessList == NULL) goto doneError3;
2489
2490 MakeList(this->NumRegionsInProcess, int ,nProcesses);
2491
2492 if (this->NumRegionsInProcess == NULL) goto doneError3;
2493
2494 MakeList(this->RegionList, int * ,nProcesses);
2495
2496 if (this->RegionList == NULL) goto doneError3;
2497
2498 MakeList(this->CellCountList, vtkIdType * ,nRegions);
2499
2500 if (this->CellCountList == NULL) goto doneError3;
2501
2502 return 0;
2503
2504 doneError3:
2505 this->FreeProcessDataLists();
2506 return 1;
2507 }
FreeProcessDataLists()2508 void vtkPKdTree::FreeProcessDataLists()
2509 {
2510 int nRegions = this->GetNumberOfRegions();
2511 int nProcesses = this->NumProcesses;
2512
2513 FreeListOfLists(this->CellCountList, nRegions);
2514
2515 FreeListOfLists(this->RegionList, nProcesses);
2516
2517 FreeList(this->NumRegionsInProcess);
2518
2519 FreeListOfLists(this->ProcessList, nRegions);
2520
2521 FreeList(this->NumProcessesInRegion);
2522
2523 FreeList(this->DataLocationMap);
2524 }
2525
2526 // Field array global min and max -----------------------------------
2527
InitializeFieldArrayMinMax()2528 void vtkPKdTree::InitializeFieldArrayMinMax()
2529 {
2530 this->NumCellArrays = this->NumPointArrays = 0;
2531 this->CellDataMin = this->CellDataMax = NULL;
2532 this->PointDataMin = this->PointDataMax = NULL;
2533 this->CellDataName = NULL;
2534 this->PointDataName = NULL;
2535 }
2536
AllocateAndZeroFieldArrayMinMax()2537 int vtkPKdTree::AllocateAndZeroFieldArrayMinMax()
2538 {
2539 int iNumCellArrays = 0;
2540 int iNumPointArrays = 0;
2541
2542 for (int set=0; set < this->GetNumberOfDataSets(); set++)
2543 {
2544 iNumCellArrays +=
2545 this->GetDataSet(set)->GetCellData()->GetNumberOfArrays();
2546 iNumPointArrays +=
2547 this->GetDataSet(set)->GetPointData()->GetNumberOfArrays();
2548 }
2549
2550 this->FreeFieldArrayMinMax();
2551
2552 if (iNumCellArrays > 0)
2553 {
2554 MakeList(this->CellDataMin, double, iNumCellArrays);
2555 if (this->CellDataMin == NULL) goto doneError5;
2556
2557 MakeList(this->CellDataMax, double, iNumCellArrays);
2558 if (this->CellDataMax == NULL) goto doneError5;
2559
2560 MakeList(this->CellDataName, char *, iNumCellArrays);
2561 if (this->CellDataName == NULL) goto doneError5;
2562 }
2563
2564 this->NumCellArrays = iNumCellArrays;
2565
2566 if (iNumPointArrays > 0)
2567 {
2568 MakeList(this->PointDataMin, double, iNumPointArrays);
2569 if (this->PointDataMin == NULL) goto doneError5;
2570
2571 MakeList(this->PointDataMax, double, iNumPointArrays);
2572 if (this->PointDataMax == NULL) goto doneError5;
2573
2574 MakeList(this->PointDataName, char *, iNumPointArrays);
2575 if (this->PointDataName == NULL) goto doneError5;
2576 }
2577
2578 this->NumPointArrays = iNumPointArrays;
2579
2580 return 0;
2581
2582 doneError5:
2583 this->FreeFieldArrayMinMax();
2584 return 1;
2585 }
FreeFieldArrayMinMax()2586 void vtkPKdTree::FreeFieldArrayMinMax()
2587 {
2588 FreeList(this->CellDataMin);
2589 FreeList(this->CellDataMax);
2590 FreeList(this->PointDataMin);
2591 FreeList(this->PointDataMax);
2592
2593 FreeListOfLists(this->CellDataName, this->NumCellArrays);
2594 FreeListOfLists(this->PointDataName, this->NumPointArrays);
2595
2596 this->NumCellArrays = this->NumPointArrays = 0;
2597 }
2598
ReleaseTables()2599 void vtkPKdTree::ReleaseTables()
2600 {
2601 if (this->RegionAssignment != UserDefinedAssignment)
2602 {
2603 this->FreeRegionAssignmentLists();
2604 }
2605 this->FreeProcessDataLists();
2606 this->FreeFieldArrayMinMax();
2607 }
2608
2609 //--------------------------------------------------------------------
2610 // Create tables indicating which processes have data for which
2611 // regions.
2612 //--------------------------------------------------------------------
2613
CreateProcessCellCountData()2614 int vtkPKdTree::CreateProcessCellCountData()
2615 {
2616 int proc, reg;
2617 int retval = 0;
2618 int *cellCounts = NULL;
2619 int *tempbuf;
2620 char *procData, *myData;
2621
2622 tempbuf = NULL;
2623 procData = myData = NULL;
2624
2625 this->SubGroup = vtkSubGroup::New();
2626 this->SubGroup->Initialize(0, this->NumProcesses-1,
2627 this->MyId, 0x0000f000, this->Controller->GetCommunicator());
2628
2629 int fail = this->AllocateAndZeroProcessDataLists();
2630
2631 if (!fail && !this->Top)
2632 {
2633 fail = 1;
2634 }
2635
2636 if (this->AllCheckForFailure(fail,
2637 "BuildRegionProcessTables",
2638 "memory allocation"))
2639 {
2640 this->FreeProcessDataLists();
2641 this->SubGroup->Delete();
2642 this->SubGroup = NULL;
2643 return 1;
2644 }
2645
2646 // Build table indicating which processes have data for which regions
2647
2648 cellCounts = this->CollectLocalRegionProcessData();
2649
2650 fail = (cellCounts == NULL);
2651
2652 if (this->AllCheckForFailure(fail,"BuildRegionProcessTables","error"))
2653 {
2654 goto doneError4;
2655 }
2656
2657 myData = this->DataLocationMap + (this->MyId * this->GetNumberOfRegions());
2658
2659 for (reg=0; reg < this->GetNumberOfRegions(); reg++)
2660 {
2661 if (cellCounts[reg] > 0) myData[reg] = 1;
2662 }
2663
2664 if (this->NumProcesses > 1)
2665 {
2666 this->SubGroup->Gather(myData, this->DataLocationMap,
2667 this->GetNumberOfRegions(), 0);
2668 this->SubGroup->Broadcast(this->DataLocationMap,
2669 this->GetNumberOfRegions() * this->NumProcesses, 0);
2670 }
2671
2672 // Other helpful tables - not the fastest way to create this
2673 // information, but it uses the least memory
2674
2675 procData = this->DataLocationMap;
2676
2677 for (proc=0; proc<this->NumProcesses; proc++)
2678 {
2679 for (reg=0; reg < this->GetNumberOfRegions(); reg++)
2680 {
2681 if (*procData)
2682 {
2683 this->NumProcessesInRegion[reg]++;
2684 this->NumRegionsInProcess[proc]++;
2685 }
2686 procData++;
2687 }
2688 }
2689 for (reg=0; reg < this->GetNumberOfRegions(); reg++)
2690 {
2691 int nprocs = this->NumProcessesInRegion[reg];
2692
2693 if (nprocs > 0)
2694 {
2695 this->ProcessList[reg] = new int [nprocs];
2696 this->ProcessList[reg][0] = -1;
2697 this->CellCountList[reg] = new vtkIdType [nprocs];
2698 this->CellCountList[reg][0] = -1;
2699 }
2700 }
2701 for (proc=0; proc < this->NumProcesses; proc++)
2702 {
2703 int nregs = this->NumRegionsInProcess[proc];
2704
2705 if (nregs > 0)
2706 {
2707 this->RegionList[proc] = new int [nregs];
2708 this->RegionList[proc][0] = -1;
2709 }
2710 }
2711
2712 procData = this->DataLocationMap;
2713
2714 for (proc=0; proc<this->NumProcesses; proc++)
2715 {
2716
2717 for (reg=0; reg < this->GetNumberOfRegions(); reg++)
2718 {
2719 if (*procData)
2720 {
2721 this->AddEntry(this->ProcessList[reg],
2722 this->NumProcessesInRegion[reg], proc);
2723
2724 this->AddEntry(this->RegionList[proc],
2725 this->NumRegionsInProcess[proc], reg);
2726 }
2727 procData++;
2728 }
2729 }
2730
2731 // Cell counts per process per region
2732
2733 if (this->NumProcesses > 1)
2734 {
2735 tempbuf = new int [this->GetNumberOfRegions() * this->NumProcesses];
2736
2737 fail = (tempbuf == NULL);
2738
2739 if (this->AllCheckForFailure(fail,
2740 "BuildRegionProcessTables",
2741 "memory allocation"))
2742 {
2743 goto doneError4;
2744 }
2745
2746 this->SubGroup->Gather(cellCounts, tempbuf, this->GetNumberOfRegions(), 0);
2747 this->SubGroup->Broadcast(tempbuf,
2748 this->NumProcesses*this->GetNumberOfRegions(),
2749 0);
2750 }
2751 else
2752 {
2753 tempbuf = cellCounts;
2754 }
2755
2756 for (proc=0; proc<this->NumProcesses; proc++)
2757 {
2758 int *procCount = tempbuf + (proc * this->GetNumberOfRegions());
2759
2760 for (reg=0; reg < this->GetNumberOfRegions(); reg++)
2761 {
2762 if (procCount[reg]> 0)
2763 {
2764 this->AddEntry(this->CellCountList[reg],
2765 this->NumProcessesInRegion[reg],
2766 procCount[reg]);
2767 }
2768 }
2769 }
2770
2771 goto done4;
2772
2773 doneError4:
2774
2775 this->FreeProcessDataLists();
2776 retval = 1;
2777
2778 done4:
2779
2780 if (tempbuf != cellCounts)
2781 {
2782 FreeList(tempbuf);
2783 }
2784
2785 FreeList(cellCounts);
2786 this->SubGroup->Delete();
2787 this->SubGroup = NULL;
2788
2789 return retval;
2790 }
2791
2792 //--------------------------------------------------------------------
2793 // Create list of global min and max for cell and point field arrays
2794 //--------------------------------------------------------------------
2795
CreateGlobalDataArrayBounds()2796 int vtkPKdTree::CreateGlobalDataArrayBounds()
2797 {
2798 int set = 0;
2799 this->SubGroup = NULL;
2800
2801 if (this->NumProcesses > 1)
2802 {
2803 this->SubGroup = vtkSubGroup::New();
2804 this->SubGroup->Initialize(0, this->NumProcesses-1,
2805 this->MyId, 0x0000f000, this->Controller->GetCommunicator());
2806 }
2807
2808 int fail = this->AllocateAndZeroFieldArrayMinMax();
2809
2810 if (this->AllCheckForFailure(fail, "BuildFieldArrayMinMax", "memory allocation"))
2811 {
2812 this->FreeFieldArrayMinMax();
2813 FreeObject(this->SubGroup);
2814 return 1;
2815 }
2816
2817 TIMER("Get global ranges");
2818
2819 int ar;
2820 double range[2];
2821 int nc = 0;
2822 int np = 0;
2823
2824 // This code assumes that if more than one dataset was input to vtkPKdTree,
2825 // each process input the datasets in the same order.
2826
2827 if (this->NumCellArrays > 0)
2828 {
2829 for (set=0; set < this->GetNumberOfDataSets(); set++)
2830 {
2831 int ncellarrays = this->GetDataSet(set)->GetCellData()->GetNumberOfArrays();
2832
2833 for (ar=0; ar<ncellarrays; ar++)
2834 {
2835 vtkDataArray *array = this->GetDataSet(set)->GetCellData()->GetArray(ar);
2836
2837 array->GetRange(range);
2838
2839 this->CellDataMin[nc] = range[0];
2840 this->CellDataMax[nc] = range[1];
2841
2842 this->CellDataName[nc] = vtkPKdTree::StrDupWithNew(array->GetName());
2843 nc++;
2844 }
2845 }
2846
2847 if (this->NumProcesses > 1)
2848 {
2849 this->SubGroup->ReduceMin(this->CellDataMin, this->CellDataMin, nc, 0);
2850 this->SubGroup->Broadcast(this->CellDataMin, nc, 0);
2851
2852 this->SubGroup->ReduceMax(this->CellDataMax, this->CellDataMax, nc, 0);
2853 this->SubGroup->Broadcast(this->CellDataMax, nc, 0);
2854 }
2855 }
2856
2857 if (this->NumPointArrays > 0)
2858 {
2859 for (set=0; set < this->GetNumberOfDataSets(); set++)
2860 {
2861 int npointarrays = this->GetDataSet(set)->GetPointData()->GetNumberOfArrays();
2862
2863 for (ar=0; ar<npointarrays; ar++)
2864 {
2865 vtkDataArray *array = this->GetDataSet(set)->GetPointData()->GetArray(ar);
2866
2867 array->GetRange(range);
2868
2869 this->PointDataMin[np] = range[0];
2870 this->PointDataMax[np] = range[1];
2871
2872 this->PointDataName[np] = StrDupWithNew(array->GetName());
2873 np++;
2874 }
2875 }
2876
2877 if (this->NumProcesses > 1)
2878 {
2879 this->SubGroup->ReduceMin(this->PointDataMin, this->PointDataMin, np, 0);
2880 this->SubGroup->Broadcast(this->PointDataMin, np, 0);
2881
2882 this->SubGroup->ReduceMax(this->PointDataMax, this->PointDataMax, np, 0);
2883 this->SubGroup->Broadcast(this->PointDataMax, np, 0);
2884 }
2885 }
2886
2887 TIMERDONE("Get global ranges");
2888
2889 FreeObject(this->SubGroup);
2890
2891 return 0;
2892 }
CollectLocalRegionProcessData()2893 int *vtkPKdTree::CollectLocalRegionProcessData()
2894 {
2895 int *cellCounts = NULL;
2896
2897 int numRegions = this->GetNumberOfRegions();
2898
2899 MakeList(cellCounts, int, numRegions);
2900
2901 if (!cellCounts)
2902 {
2903 VTKERROR("CollectLocalRegionProcessData - memory allocation");
2904 return NULL;
2905 }
2906
2907 TIMER("Get cell regions");
2908
2909 int *IDs = this->AllGetRegionContainingCell();
2910
2911 TIMERDONE("Get cell regions");
2912
2913 for (int set=0; set < this->GetNumberOfDataSets(); set++)
2914 {
2915 int ncells = this->GetDataSet(set)->GetNumberOfCells();
2916
2917 TIMER("Increment cell counts");
2918
2919 for (int i=0; i<ncells; i++)
2920 {
2921 int regionId = IDs[i];
2922
2923 if ( (regionId < 0) || (regionId >= numRegions))
2924 {
2925 VTKERROR("CollectLocalRegionProcessData - corrupt data");
2926 delete [] cellCounts;
2927 return NULL;
2928 }
2929 cellCounts[regionId]++;
2930 }
2931
2932 IDs += ncells;
2933
2934 TIMERDONE("Increment cell counts");
2935 }
2936
2937 return cellCounts;
2938 }
AddEntry(int * list,int len,int id)2939 void vtkPKdTree::AddEntry(int *list, int len, int id)
2940 {
2941 int i=0;
2942
2943 while ((i < len) && (list[i] != -1)) i++;
2944
2945 if (i == len) return; // error
2946
2947 list[i++] = id;
2948
2949 if (i < len) list[i] = -1;
2950 }
2951 #ifdef VTK_USE_64BIT_IDS
AddEntry(vtkIdType * list,int len,vtkIdType id)2952 void vtkPKdTree::AddEntry(vtkIdType *list, int len, vtkIdType id)
2953 {
2954 int i=0;
2955
2956 while ((i < len) && (list[i] != -1)) i++;
2957
2958 if (i == len) return; // error
2959
2960 list[i++] = id;
2961
2962 if (i < len) list[i] = -1;
2963 }
2964 #endif
BinarySearch(vtkIdType * list,int len,vtkIdType which)2965 int vtkPKdTree::BinarySearch(vtkIdType *list, int len, vtkIdType which)
2966 {
2967 vtkIdType mid, left, right;
2968
2969 mid = -1;
2970
2971 if (len <= 3)
2972 {
2973 for (int i=0; i<len; i++)
2974 {
2975 if (list[i] == which)
2976 {
2977 mid=i;
2978 break;
2979 }
2980 }
2981 }
2982 else
2983 {
2984 mid = len >> 1;
2985 left = 0;
2986 right = len-1;
2987
2988 while (list[mid] != which)
2989 {
2990 if (list[mid] < which)
2991 {
2992 left = mid+1;
2993 }
2994 else
2995 {
2996 right = mid-1;
2997 }
2998
2999 if (right > left+1)
3000 {
3001 mid = (left + right) >> 1;
3002 }
3003 else
3004 {
3005 if (list[left] == which) mid = left;
3006 else if (list[right] == which) mid = right;
3007 else mid = -1;
3008 break;
3009 }
3010 }
3011 }
3012 return mid;
3013 }
3014 //--------------------------------------------------------------------
3015 // Assign responsibility for each spatial region to one process
3016 //--------------------------------------------------------------------
3017
UpdateRegionAssignment()3018 int vtkPKdTree::UpdateRegionAssignment()
3019 {
3020 int returnVal = 0;
3021
3022 if (this->RegionAssignment== ContiguousAssignment)
3023 {
3024 returnVal = this->AssignRegionsContiguous();
3025 }
3026 else if (this->RegionAssignment== RoundRobinAssignment)
3027 {
3028 returnVal = this->AssignRegionsRoundRobin();
3029 }
3030
3031 return returnVal;
3032 }
AssignRegionsRoundRobin()3033 int vtkPKdTree::AssignRegionsRoundRobin()
3034 {
3035 this->RegionAssignment = RoundRobinAssignment;
3036
3037 if (this->Top == NULL)
3038 {
3039 return 0;
3040 }
3041
3042 int nProcesses = this->NumProcesses;
3043 int nRegions = this->GetNumberOfRegions();
3044
3045 int fail = this->AllocateAndZeroRegionAssignmentLists();
3046
3047 if (fail)
3048 {
3049 return 1;
3050 }
3051
3052 for (int i=0, procID=0; i<nRegions; i++)
3053 {
3054 this->RegionAssignmentMap[i] = procID;
3055 this->NumRegionsAssigned[procID]++;
3056
3057 procID = ( (procID == nProcesses-1) ? 0 : procID+1);
3058 }
3059 this->BuildRegionListsForProcesses();
3060
3061 return 0;
3062 }
AssignRegions(int * map,int len)3063 int vtkPKdTree::AssignRegions(int *map, int len)
3064 {
3065 int fail = this->AllocateAndZeroRegionAssignmentLists();
3066
3067 if (fail)
3068 {
3069 return 1;
3070 }
3071
3072 this->RegionAssignmentMapLength = len;
3073
3074 this->RegionAssignment = UserDefinedAssignment;
3075
3076 for (int i=0; i<len; i++)
3077 {
3078 if ( (map[i] < 0) || (map[i] >= this->NumProcesses))
3079 {
3080 this->FreeRegionAssignmentLists();
3081 VTKERROR("AssignRegions - invalid process id " << map[i]);
3082 return 1;
3083 }
3084
3085 this->RegionAssignmentMap[i] = map[i];
3086 this->NumRegionsAssigned[map[i]]++;
3087 }
3088
3089 this->BuildRegionListsForProcesses();
3090
3091 return 0;
3092 }
AddProcessRegions(int procId,vtkKdNode * kd)3093 void vtkPKdTree::AddProcessRegions(int procId, vtkKdNode *kd)
3094 {
3095 vtkIntArray *leafNodeIds = vtkIntArray::New();
3096
3097 vtkKdTree::GetLeafNodeIds(kd, leafNodeIds);
3098
3099 int nLeafNodes = leafNodeIds->GetNumberOfTuples();
3100
3101 for (int n=0; n<nLeafNodes; n++)
3102 {
3103 this->RegionAssignmentMap[leafNodeIds->GetValue(n)] = procId;
3104 this->NumRegionsAssigned[procId]++;
3105 }
3106
3107 leafNodeIds->Delete();
3108 }
AssignRegionsContiguous()3109 int vtkPKdTree::AssignRegionsContiguous()
3110 {
3111 int p;
3112
3113 this->RegionAssignment = ContiguousAssignment;
3114
3115 if (this->Top == NULL)
3116 {
3117 return 0;
3118 }
3119
3120 int nProcesses = this->NumProcesses;
3121 int nRegions = this->GetNumberOfRegions();
3122
3123 if (nRegions <= nProcesses)
3124 {
3125 this->AssignRegionsRoundRobin();
3126 this->RegionAssignment = ContiguousAssignment;
3127 return 0;
3128 }
3129
3130 int fail = this->AllocateAndZeroRegionAssignmentLists();
3131
3132 if (fail)
3133 {
3134 return 1;
3135 }
3136
3137 int floorLogP, ceilLogP;
3138
3139 for (floorLogP = 0; (nProcesses >> floorLogP) > 0; floorLogP++)
3140 {
3141 // empty loop.
3142 }
3143 floorLogP--;
3144
3145 int P = 1 << floorLogP;
3146
3147 if (nProcesses == P)
3148 {
3149 ceilLogP = floorLogP;
3150 }
3151 else
3152 {
3153 ceilLogP = floorLogP + 1;
3154 }
3155
3156 vtkKdNode **nodes = new vtkKdNode* [P];
3157
3158 this->GetRegionsAtLevel(floorLogP, nodes);
3159
3160 if (floorLogP == ceilLogP)
3161 {
3162 for (p=0; p<nProcesses; p++)
3163 {
3164 this->AddProcessRegions(p, nodes[p]);
3165 }
3166 }
3167 else
3168 {
3169 int nodesLeft = 1 << ceilLogP;
3170 int procsLeft = nProcesses;
3171 int procId = 0;
3172
3173 for (int i=0; i<P; i++)
3174 {
3175 if (nodesLeft > procsLeft)
3176 {
3177 this->AddProcessRegions(procId, nodes[i]);
3178
3179 procsLeft -= 1;
3180 procId += 1;
3181 }
3182 else
3183 {
3184 this->AddProcessRegions(procId, nodes[i]->GetLeft());
3185 this->AddProcessRegions(procId+1, nodes[i]->GetRight());
3186
3187 procsLeft -= 2;
3188 procId += 2;
3189 }
3190 nodesLeft -= 2;
3191 }
3192 }
3193
3194 delete [] nodes;
3195
3196 this->BuildRegionListsForProcesses();
3197
3198 return 0;
3199 }
BuildRegionListsForProcesses()3200 void vtkPKdTree::BuildRegionListsForProcesses()
3201 {
3202 int *count = new int [this->NumProcesses];
3203
3204 for (int p=0; p<this->NumProcesses; p++)
3205 {
3206 int nregions = this->NumRegionsAssigned[p];
3207
3208 if (nregions > 0)
3209 {
3210 this->ProcessAssignmentMap[p] = new int [nregions];
3211 }
3212 else
3213 {
3214 this->ProcessAssignmentMap[p] = NULL;
3215 }
3216
3217 count[p] = 0;
3218 }
3219
3220 for (int r=0; r<this->RegionAssignmentMapLength; r++)
3221 {
3222 int proc = this->RegionAssignmentMap[r];
3223 int next = count[proc];
3224
3225 this->ProcessAssignmentMap[proc][next] = r;
3226
3227 count[proc] = next + 1;
3228 }
3229
3230 delete [] count;
3231 }
3232
3233 //--------------------------------------------------------------------
3234 // Queries
3235 //--------------------------------------------------------------------
FindNextLocalArrayIndex(const char * n,const char ** names,int len,int start)3236 int vtkPKdTree::FindNextLocalArrayIndex(const char *n,
3237 const char **names, int len, int start)
3238 {
3239 int index = -1;
3240 size_t nsize = strlen(n);
3241
3242 // normally a very small list, maybe 1 to 5 names
3243
3244 for (int i=start; i<len; i++)
3245 {
3246 if (!strncmp(n, names[i], nsize))
3247 {
3248 index = i;
3249 break;
3250 }
3251 }
3252 return index;
3253 }
GetCellArrayGlobalRange(const char * n,double range[2])3254 int vtkPKdTree::GetCellArrayGlobalRange(const char *n, double range[2])
3255 {
3256 int first = 1;
3257 double tmp[2] = {0, 0};
3258 int start = 0;
3259
3260 while (1)
3261 {
3262 // Cell array name may appear more than once if multiple datasets
3263 // were processed.
3264
3265 int index = vtkPKdTree::FindNextLocalArrayIndex(n,
3266 (const char **)this->CellDataName, this->NumCellArrays, start);
3267
3268 if (index >= 0)
3269 {
3270 if (first)
3271 {
3272 this->GetCellArrayGlobalRange(index, range);
3273 first = 0;
3274 }
3275 else
3276 {
3277 this->GetCellArrayGlobalRange(index, tmp);
3278 range[0] = (tmp[0] < range[0]) ? tmp[0] : range[0];
3279 range[1] = (tmp[1] > range[1]) ? tmp[1] : range[1];
3280 }
3281 start = index + 1;
3282 }
3283 else
3284 {
3285 break;
3286 }
3287 }
3288
3289 int fail = (first != 0);
3290
3291 return fail;
3292 }
GetCellArrayGlobalRange(const char * n,float range[2])3293 int vtkPKdTree::GetCellArrayGlobalRange(const char *n, float range[2])
3294 {
3295 double tmp[2] = {0, 0 };
3296
3297 int fail = this->GetCellArrayGlobalRange(n, tmp);
3298
3299 if (!fail)
3300 {
3301 range[0] = (float)tmp[0];
3302 range[1] = (float)tmp[1];
3303 }
3304
3305 return fail;
3306 }
GetPointArrayGlobalRange(const char * n,double range[2])3307 int vtkPKdTree::GetPointArrayGlobalRange(const char *n, double range[2])
3308 {
3309 int first = 1;
3310 double tmp[2]={0, 0};
3311 int start = 0;
3312
3313 while (1)
3314 {
3315 // Point array name may appear more than once if multiple datasets
3316 // were processed.
3317
3318 int index = vtkPKdTree::FindNextLocalArrayIndex(n,
3319 (const char **)this->PointDataName, this->NumPointArrays, start);
3320
3321 if (index >= 0)
3322 {
3323 if (first)
3324 {
3325 this->GetPointArrayGlobalRange(index, range);
3326 first = 0;
3327 }
3328 else
3329 {
3330 this->GetPointArrayGlobalRange(index, tmp);
3331 range[0] = (tmp[0] < range[0]) ? tmp[0] : range[0];
3332 range[1] = (tmp[1] > range[1]) ? tmp[1] : range[1];
3333 }
3334 start = index + 1;
3335 }
3336 else
3337 {
3338 break;
3339 }
3340 }
3341
3342 int fail = (first != 0);
3343
3344 return fail;
3345 }
GetPointArrayGlobalRange(const char * n,float range[2])3346 int vtkPKdTree::GetPointArrayGlobalRange(const char *n, float range[2])
3347 {
3348 double tmp[2] = {0, 0};
3349
3350 int fail = this->GetPointArrayGlobalRange(n, tmp);
3351
3352 if (!fail)
3353 {
3354 range[0] = (float)tmp[0];
3355 range[1] = (float)tmp[1];
3356 }
3357
3358 return fail;
3359 }
GetCellArrayGlobalRange(int arrayIndex,float range[2])3360 int vtkPKdTree::GetCellArrayGlobalRange(int arrayIndex, float range[2])
3361 {
3362 double tmp[2];
3363 int fail = this->GetCellArrayGlobalRange(arrayIndex, tmp);
3364 if (!fail)
3365 {
3366 range[0] = (float)tmp[0];
3367 range[1] = (float)tmp[1];
3368 }
3369 return fail;
3370 }
GetCellArrayGlobalRange(int arrayIndex,double range[2])3371 int vtkPKdTree::GetCellArrayGlobalRange(int arrayIndex, double range[2])
3372 {
3373 if ((arrayIndex < 0) || (arrayIndex >= this->NumCellArrays))
3374 {
3375 return 1;
3376 }
3377 if (this->CellDataMin == NULL) return 1;
3378
3379 range[0] = this->CellDataMin[arrayIndex];
3380 range[1] = this->CellDataMax[arrayIndex];
3381
3382 return 0;
3383 }
GetPointArrayGlobalRange(int arrayIndex,float range[2])3384 int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, float range[2])
3385 {
3386 double tmp[2];
3387 int fail = this->GetPointArrayGlobalRange(arrayIndex, tmp);
3388 if (!fail)
3389 {
3390 range[0] = (float)tmp[0];
3391 range[1] = (float)tmp[1];
3392 }
3393 return fail;
3394 }
GetPointArrayGlobalRange(int arrayIndex,double range[2])3395 int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, double range[2])
3396 {
3397 if ((arrayIndex < 0) || (arrayIndex >= this->NumPointArrays))
3398 {
3399 return 1;
3400 }
3401 if (this->PointDataMin == NULL) return 1;
3402
3403 range[0] = this->PointDataMin[arrayIndex];
3404 range[1] = this->PointDataMax[arrayIndex];
3405
3406 return 0;
3407 }
3408
ViewOrderAllProcessesInDirection(const double dop[3],vtkIntArray * orderedList)3409 int vtkPKdTree::ViewOrderAllProcessesInDirection(const double dop[3],
3410 vtkIntArray *orderedList)
3411 {
3412 assert("pre: orderedList_exists" && orderedList!=0);
3413
3414 vtkIntArray *regionList = vtkIntArray::New();
3415
3416 this->ViewOrderAllRegionsInDirection(dop, regionList);
3417
3418 orderedList->SetNumberOfValues(this->NumProcesses);
3419
3420 int nextId = 0;
3421
3422 // if regions were not assigned contiguously, this
3423 // produces the wrong result
3424
3425 for (int r=0; r<this->GetNumberOfRegions(); )
3426 {
3427 int procId = this->RegionAssignmentMap[regionList->GetValue(r)];
3428
3429 orderedList->SetValue(nextId++, procId);
3430
3431 int nregions = this->NumRegionsAssigned[procId];
3432
3433 r += nregions;
3434 }
3435
3436 regionList->Delete();
3437
3438 return this->NumProcesses;
3439 }
3440
ViewOrderAllProcessesFromPosition(const double pos[3],vtkIntArray * orderedList)3441 int vtkPKdTree::ViewOrderAllProcessesFromPosition(const double pos[3],
3442 vtkIntArray *orderedList)
3443 {
3444 assert("pre: orderedList_exists" && orderedList!=0);
3445
3446 vtkIntArray *regionList = vtkIntArray::New();
3447
3448 this->ViewOrderAllRegionsFromPosition(pos, regionList);
3449
3450 orderedList->SetNumberOfValues(this->NumProcesses);
3451
3452 int nextId = 0;
3453
3454 // if regions were not assigned contiguously, this
3455 // produces the wrong result
3456
3457 for (int r=0; r<this->GetNumberOfRegions(); )
3458 {
3459 int procId = this->RegionAssignmentMap[regionList->GetValue(r)];
3460
3461 orderedList->SetValue(nextId++, procId);
3462
3463 int nregions = this->NumRegionsAssigned[procId];
3464
3465 r += nregions;
3466 }
3467
3468 regionList->Delete();
3469
3470 return this->NumProcesses;
3471 }
3472
GetRegionAssignmentList(int procId,vtkIntArray * list)3473 int vtkPKdTree::GetRegionAssignmentList(int procId, vtkIntArray *list)
3474 {
3475 if ( (procId < 0) || (procId >= this->NumProcesses))
3476 {
3477 VTKERROR("GetRegionAssignmentList - invalid process id");
3478 return 0;
3479 }
3480
3481 if (!this->RegionAssignmentMap)
3482 {
3483 this->UpdateRegionAssignment();
3484
3485 if (!this->RegionAssignmentMap)
3486 {
3487 return 0;
3488 }
3489 }
3490
3491 int nregions = this->NumRegionsAssigned[procId];
3492 int *regionIds = this->ProcessAssignmentMap[procId];
3493
3494 list->Initialize();
3495 list->SetNumberOfValues(nregions);
3496
3497 for (int i=0; i<nregions; i++)
3498 {
3499 list->SetValue(i, regionIds[i]);
3500 }
3501
3502 return nregions;
3503 }
3504
GetAllProcessesBorderingOnPoint(float x,float y,float z,vtkIntArray * list)3505 void vtkPKdTree::GetAllProcessesBorderingOnPoint(float x, float y, float z,
3506 vtkIntArray *list)
3507 {
3508 vtkIntArray *regions = vtkIntArray::New();
3509 double *subRegionBounds;
3510 list->Initialize();
3511
3512 for (int procId = 0; procId < this->NumProcesses; procId++)
3513 {
3514 this->GetRegionAssignmentList(procId, regions);
3515
3516 int nSubRegions = this->MinimalNumberOfConvexSubRegions(
3517 regions, &subRegionBounds);
3518
3519 double *b = subRegionBounds;
3520
3521 for (int r=0; r<nSubRegions; r++)
3522 {
3523 if ((((x == b[0]) || (x == b[1])) &&
3524 ((y >= b[2]) && (y <= b[3]) && (z >= b[4]) && (z <= b[5]))) ||
3525 (((y == b[2]) || (y == b[3])) &&
3526 ((x >= b[0]) && (x <= b[1]) && (z >= b[4]) && (z <= b[5]))) ||
3527 (((z == b[4]) || (z == b[5])) &&
3528 ((x >= b[0]) && (x <= b[1]) && (y >= b[2]) && (y <= b[3]))) )
3529 {
3530 list->InsertNextValue(procId);
3531 break;
3532 }
3533
3534 b += 6;
3535 }
3536 }
3537
3538 regions->Delete();
3539 }
3540
GetProcessAssignedToRegion(int regionId)3541 int vtkPKdTree::GetProcessAssignedToRegion(int regionId)
3542 {
3543 if ( !this->RegionAssignmentMap ||
3544 (regionId < 0) || (regionId >= this->GetNumberOfRegions()))
3545 {
3546 return -1;
3547 }
3548
3549 return this->RegionAssignmentMap[regionId];
3550 }
HasData(int processId,int regionId)3551 int vtkPKdTree::HasData(int processId, int regionId)
3552 {
3553 if ((!this->DataLocationMap) ||
3554 (processId < 0) || (processId >= this->NumProcesses) ||
3555 (regionId < 0) || (regionId >= this->GetNumberOfRegions()) )
3556 {
3557 VTKERROR("HasData - invalid request");
3558 return 0;
3559 }
3560
3561 int where = this->GetNumberOfRegions() * processId + regionId;
3562
3563 return this->DataLocationMap[where];
3564 }
3565
GetTotalProcessesInRegion(int regionId)3566 int vtkPKdTree::GetTotalProcessesInRegion(int regionId)
3567 {
3568 if (!this->NumProcessesInRegion ||
3569 (regionId < 0) || (regionId >= this->GetNumberOfRegions()) )
3570 {
3571 VTKERROR("GetTotalProcessesInRegion - invalid request");
3572 return 0;
3573 }
3574
3575 return this->NumProcessesInRegion[regionId];
3576 }
3577
GetProcessListForRegion(int regionId,vtkIntArray * processes)3578 int vtkPKdTree::GetProcessListForRegion(int regionId, vtkIntArray *processes)
3579 {
3580 if (!this->ProcessList ||
3581 (regionId < 0) || (regionId >= this->GetNumberOfRegions()) )
3582 {
3583 VTKERROR("GetProcessListForRegion - invalid request");
3584 return 0;
3585 }
3586
3587 int nProcesses = this->NumProcessesInRegion[regionId];
3588
3589 for (int i=0; i<nProcesses; i++)
3590 {
3591 processes->InsertNextValue(this->ProcessList[regionId][i]);
3592 }
3593
3594 return nProcesses;
3595 }
3596
GetProcessesCellCountForRegion(int regionId,int * count,int len)3597 int vtkPKdTree::GetProcessesCellCountForRegion(int regionId, int *count, int len)
3598 {
3599 if (!this->CellCountList ||
3600 (regionId < 0) || (regionId >= this->GetNumberOfRegions()) )
3601 {
3602 VTKERROR("GetProcessesCellCountForRegion - invalid request");
3603 return 0;
3604 }
3605
3606 int nProcesses = this->NumProcessesInRegion[regionId];
3607
3608 nProcesses = (len < nProcesses) ? len : nProcesses;
3609
3610 for (int i=0; i<nProcesses; i++)
3611 {
3612 count[i] = this->CellCountList[regionId][i];
3613 }
3614
3615 return nProcesses;
3616 }
3617
GetProcessCellCountForRegion(int processId,int regionId)3618 int vtkPKdTree::GetProcessCellCountForRegion(int processId, int regionId)
3619 {
3620 int nCells;
3621
3622 if (!this->CellCountList ||
3623 (regionId < 0) || (regionId >= this->GetNumberOfRegions()) ||
3624 (processId < 0) || (processId >= this->NumProcesses))
3625 {
3626 VTKERROR("GetProcessCellCountForRegion - invalid request");
3627 return 0;
3628 }
3629
3630 int nProcesses = this->NumProcessesInRegion[regionId];
3631
3632 int which = -1;
3633
3634 for (int i=0; i<nProcesses; i++)
3635 {
3636 if (this->ProcessList[regionId][i] == processId)
3637 {
3638 which = i;
3639 break;
3640 }
3641 }
3642
3643 if (which == -1) nCells = 0;
3644 else nCells = this->CellCountList[regionId][which];
3645
3646 return nCells;
3647 }
3648
GetTotalRegionsForProcess(int processId)3649 int vtkPKdTree::GetTotalRegionsForProcess(int processId)
3650 {
3651 if ((!this->NumRegionsInProcess) ||
3652 (processId < 0) || (processId >= this->NumProcesses))
3653 {
3654 VTKERROR("GetTotalRegionsForProcess - invalid request");
3655 return 0;
3656 }
3657
3658 return this->NumRegionsInProcess[processId];
3659 }
3660
GetRegionListForProcess(int processId,vtkIntArray * regions)3661 int vtkPKdTree::GetRegionListForProcess(int processId, vtkIntArray *regions)
3662 {
3663 if (!this->RegionList ||
3664 (processId < 0) || (processId >= this->NumProcesses))
3665 {
3666 VTKERROR("GetRegionListForProcess - invalid request");
3667 return 0;
3668 }
3669
3670 int nRegions = this->NumRegionsInProcess[processId];
3671
3672 for (int i=0; i<nRegions; i++)
3673 {
3674 regions->InsertNextValue(this->RegionList[processId][i]);
3675 }
3676
3677 return nRegions;
3678 }
GetRegionsCellCountForProcess(int processId,int * count,int len)3679 int vtkPKdTree::GetRegionsCellCountForProcess(int processId, int *count, int len)
3680 {
3681 if (!this->CellCountList ||
3682 (processId < 0) || (processId >= this->NumProcesses))
3683 {
3684 VTKERROR("GetRegionsCellCountForProcess - invalid request");
3685 return 0;
3686 }
3687
3688 int nRegions = this->NumRegionsInProcess[processId];
3689
3690 nRegions = (len < nRegions) ? len : nRegions;
3691
3692 for (int i=0; i<nRegions; i++)
3693 {
3694 int regionId = this->RegionList[processId][i];
3695 int iam;
3696
3697 for (iam = 0; iam < this->NumProcessesInRegion[regionId]; iam++)
3698 {
3699 if (this->ProcessList[regionId][iam] == processId) break;
3700 }
3701
3702 count[i] = this->CellCountList[regionId][iam];
3703 }
3704 return nRegions;
3705 }
GetCellListsForProcessRegions(int processId,int set,vtkIdList * inRegionCells,vtkIdList * onBoundaryCells)3706 vtkIdType vtkPKdTree::GetCellListsForProcessRegions(int processId,
3707 int set, vtkIdList *inRegionCells, vtkIdList *onBoundaryCells)
3708 {
3709 if ( (set < 0) || (set >= this->GetNumberOfDataSets()))
3710 {
3711 vtkErrorMacro(<<"vtkPKdTree::GetCellListsForProcessRegions no such data set");
3712 return 0;
3713 }
3714 return this->GetCellListsForProcessRegions(processId,
3715 this->GetDataSet(set), inRegionCells, onBoundaryCells);
3716 }
GetCellListsForProcessRegions(int processId,vtkIdList * inRegionCells,vtkIdList * onBoundaryCells)3717 vtkIdType vtkPKdTree::GetCellListsForProcessRegions(int processId,
3718 vtkIdList *inRegionCells, vtkIdList *onBoundaryCells)
3719 {
3720 return this->GetCellListsForProcessRegions(processId,
3721 this->GetDataSet(0), inRegionCells, onBoundaryCells);
3722 }
GetCellListsForProcessRegions(int processId,vtkDataSet * set,vtkIdList * inRegionCells,vtkIdList * onBoundaryCells)3723 vtkIdType vtkPKdTree::GetCellListsForProcessRegions(int processId,
3724 vtkDataSet *set,
3725 vtkIdList *inRegionCells, vtkIdList *onBoundaryCells)
3726 {
3727 vtkIdType totalCells = 0;
3728
3729 if ( (inRegionCells == NULL) && (onBoundaryCells == NULL))
3730 {
3731 return totalCells;
3732 }
3733
3734 // Get the list of regions owned by this process
3735
3736 vtkIntArray *regions = vtkIntArray::New();
3737
3738 int nregions = this->GetRegionAssignmentList(processId, regions);
3739
3740 if (nregions == 0){
3741 if (inRegionCells)
3742 {
3743 inRegionCells->Initialize();
3744 }
3745 if (onBoundaryCells)
3746 {
3747 onBoundaryCells->Initialize();
3748 }
3749 regions->Delete();
3750 return totalCells;
3751 }
3752
3753 totalCells = this->GetCellLists(regions, set, inRegionCells, onBoundaryCells);
3754
3755 regions->Delete();
3756
3757 return totalCells;
3758 }
PrintTiming(ostream & os,vtkIndent indent)3759 void vtkPKdTree::PrintTiming(ostream& os, vtkIndent indent)
3760 {
3761 os << indent << "Total cells in distributed data: " << this->TotalNumCells << endl;
3762
3763 if (this->NumProcesses)
3764 {
3765 os << indent << "Average cells per processor: " ;
3766 os << this->TotalNumCells / this->NumProcesses << endl;
3767 }
3768 vtkTimerLog::DumpLogWithIndents(&os, (float)0.0);
3769 }
PrintTables(ostream & os,vtkIndent indent)3770 void vtkPKdTree::PrintTables(ostream & os, vtkIndent indent)
3771 {
3772 int nregions = this->GetNumberOfRegions();
3773 int nprocs =this->NumProcesses;
3774 int r,p,n;
3775
3776 if (this->RegionAssignmentMap)
3777 {
3778 int *map = this->RegionAssignmentMap;
3779 int *num = this->NumRegionsAssigned;
3780 int halfr = this->RegionAssignmentMapLength/2;
3781 int halfp = nprocs/2;
3782
3783 os << indent << "Region assignments:" << endl;
3784 for (r=0; r < halfr; r++)
3785 {
3786 os << indent << " region " << r << " to process " << map[r] ;
3787 os << " region " << r+halfr << " to process " << map[r+halfr] ;
3788 os << endl;
3789 }
3790 for (p=0; p < halfp; p++)
3791 {
3792 os << indent << " " << num[p] << " regions to process " << p ;
3793 os << " " << num[p+halfp] << " regions to process " << p+ halfp ;
3794 os << endl;
3795 }
3796 if (nprocs > halfp*2)
3797 {
3798 os << indent << " " << num[nprocs-1];
3799 os << " regions to process " << nprocs-1 << endl ;
3800 }
3801 }
3802
3803 if (this->ProcessList)
3804 {
3805 os << indent << "Processes holding data for each region:" << endl;
3806 for (r=0; r<nregions; r++)
3807 {
3808 n = this->NumProcessesInRegion[r];
3809
3810 os << indent << " region " << r << " (" << n << " processes): ";
3811
3812 for (p=0; p<n; p++)
3813 {
3814 if (p && (p%10==0)) os << endl << indent << " ";
3815 os << this->ProcessList[r][p] << " " ;
3816 }
3817 os << endl;
3818 }
3819 }
3820 if (this->RegionList)
3821 {
3822 os << indent << "Regions held by each process:" << endl;
3823 for (p=0; p<nprocs; p++)
3824 {
3825 n = this->NumRegionsInProcess[p];
3826
3827 os << indent << " process " << p << " (" << n << " regions): ";
3828
3829 for (r=0; r<n; r++)
3830 {
3831 if (r && (r%10==0)) os << endl << indent << " ";
3832 os << this->RegionList[p][r] << " " ;
3833 }
3834 os << endl;
3835 }
3836 }
3837 if (this->CellCountList)
3838 {
3839 os << indent << "Number of cells per process per region:" << endl;
3840 for (r=0; r<nregions; r++)
3841 {
3842 n = this->NumProcessesInRegion[r];
3843
3844 os << indent << " region: " << r << " ";
3845 for (p=0; p<n; p++)
3846 {
3847 if (p && (p%5==0)) os << endl << indent << " ";
3848 os << this->ProcessList[r][p] << " - ";
3849 os << this->CellCountList[r][p] << " cells, ";
3850 }
3851 os << endl;
3852 }
3853 }
3854 }
3855
StrDupWithNew(const char * s)3856 char *vtkPKdTree::StrDupWithNew(const char *s)
3857 {
3858 char *newstr = NULL;
3859
3860 if (s)
3861 {
3862 size_t len = strlen(s);
3863 if (len == 0)
3864 {
3865 newstr = new char [1];
3866 newstr[0] = '\0';
3867 }
3868 else
3869 {
3870 newstr = new char [len + 1];
3871 strcpy(newstr, s);
3872 }
3873 }
3874
3875 return newstr;
3876 }
3877
PrintSelf(ostream & os,vtkIndent indent)3878 void vtkPKdTree::PrintSelf(ostream& os, vtkIndent indent)
3879 {
3880 this->Superclass::PrintSelf(os,indent);
3881
3882 os << indent << "RegionAssignment: " << this->RegionAssignment << endl;
3883
3884 os << indent << "Controller: " << this->Controller << endl;
3885 os << indent << "SubGroup: " << this->SubGroup<< endl;
3886 os << indent << "NumProcesses: " << this->NumProcesses << endl;
3887 os << indent << "MyId: " << this->MyId << endl;
3888
3889 os << indent << "RegionAssignmentMap: " << this->RegionAssignmentMap << endl;
3890 os << indent << "RegionAssignmentMapLength: "
3891 << this->RegionAssignmentMapLength << endl;
3892 os << indent << "NumRegionsAssigned: " << this->NumRegionsAssigned << endl;
3893 os << indent << "NumProcessesInRegion: " << this->NumProcessesInRegion << endl;
3894 os << indent << "ProcessList: " << this->ProcessList << endl;
3895 os << indent << "NumRegionsInProcess: " << this->NumRegionsInProcess << endl;
3896 os << indent << "RegionList: " << this->RegionList << endl;
3897 os << indent << "CellCountList: " << this->CellCountList << endl;
3898
3899 os << indent << "StartVal: " << this->StartVal << endl;
3900 os << indent << "EndVal: " << this->EndVal << endl;
3901 os << indent << "NumCells: " << this->NumCells << endl;
3902 os << indent << "TotalNumCells: " << this->TotalNumCells << endl;
3903
3904 os << indent << "PtArray: " << this->PtArray << endl;
3905 os << indent << "PtArray2: " << this->PtArray2 << endl;
3906 os << indent << "CurrentPtArray: " << this->CurrentPtArray << endl;
3907 os << indent << "NextPtArray: " << this->NextPtArray << endl;
3908 os << indent << "SelectBuffer: " << this->SelectBuffer << endl;
3909 }
3910