OctTreeOps.hpp

00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 #ifndef stk_util_util_OctTreeOps_hpp
00010 #define stk_util_util_OctTreeOps_hpp
00011 
00012 #include <utility>
00013 #include <vector>
00014 #include <map>
00015 #include <set>
00016 #include <list>
00017 
00018 #include <iostream>
00019 #include <stdexcept>
00020 
00021 #include <TPI.h>
00022 
00023 #include <stk_util/parallel/Parallel.hpp>
00024 #include <stk_util/parallel/ParallelComm.hpp>
00025 #include <stk_util/parallel/ParallelReduce.hpp>
00026 #include <stk_search/IdentProc.hpp>
00027 #include <stk_search/OctTree.hpp>
00028 
00029 namespace stk {
00030 namespace search {
00031 
00032 //----------------------------------------------------------------------
00037 void oct_tree_partition_private(
00038   const unsigned p_first ,
00039   const unsigned p_end ,
00040   const unsigned depth ,
00041   const double   tolerance ,
00042   float * const weights ,
00043   const unsigned cuts_length ,
00044   OctTreeKey * const cuts );
00045 
00046 //----------------------------------------------------------------------
00054 bool hsfc_box_covering(
00055     const float * const global_box ,
00056     const float * const small_box ,
00057     OctTreeKey  * const covering ,
00058     unsigned    &       number,
00059     double              scale
00060     );
00061 
00062 template <class DomainBoundingBox, class RangeBoundingBox>
00063 void search_tree_statistics( stk::ParallelMachine  arg_comm ,
00064                              const std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > >  & s ,
00065                              unsigned * const data )
00066 {
00067   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00068 
00069   const unsigned huge = std::numeric_limits<unsigned>::max();
00070   unsigned avg[2] = { 0 , 0 };
00071   unsigned max[2] = { 0 , 0 };
00072   unsigned min[2] ;
00073   min[0] = min[1] = huge ;
00074 
00075   typename SearchTree::const_iterator i ;
00076 
00077   unsigned rcnt = 0;
00078   unsigned dcnt = 0;
00079   for ( i = s.begin() ; i != s.end() ; ++i ) {
00080     const typename SearchTree::value_type  & inode = *i ;
00081     const unsigned d_size = inode.second.first.size();
00082     const unsigned r_size = inode.second.second.size();
00083 
00084     if (d_size > 0) {
00085       avg[0] += d_size ;
00086       if ( d_size < min[0] ) { min[0] = d_size ; }
00087       if ( max[0] < d_size ) { max[0] = d_size ; }
00088       ++dcnt;
00089     }
00090 
00091     if (r_size > 0) {
00092       avg[1] += r_size ;
00093       if ( r_size < min[1] ) { min[1] = r_size ; }
00094       if ( max[1] < r_size ) { max[1] = r_size ; }
00095       ++rcnt;
00096     }
00097   }
00098 
00099   // Average for this processor
00100 
00101   if (dcnt)
00102     avg[0] = ( avg[0] + dcnt - 1 ) / dcnt ;
00103 
00104   if ( rcnt ) {
00105     avg[1] = ( avg[1] + rcnt - 1 ) / rcnt ;
00106   }
00107 
00108   if ( min[0] == huge ) { min[0] = 0 ; }
00109   if ( min[1] == huge ) { min[1] = 0 ; }
00110 
00111   all_reduce( arg_comm, ReduceMin<2>( min ) );
00112   all_reduce( arg_comm, ReduceMax<2>( max ) );
00113   all_reduce( arg_comm, ReduceSum<2>( avg ) );
00114 
00115   const unsigned p_size = parallel_machine_size( arg_comm );
00116 
00117   // Average among all processors:
00118 
00119   avg[0] = ( avg[0] + p_size - 1 ) / p_size ;
00120   avg[1] = ( avg[1] + p_size - 1 ) / p_size ;
00121 
00122   data[0] = min[0] ;
00123   data[1] = max[0] ;
00124   data[2] = avg[0] ;
00125 
00126   data[3] = min[1] ;
00127   data[4] = max[1] ;
00128   data[5] = avg[1] ;
00129 }
00130 
00131 //----------------------------------------------------------------------
00137 template <class DomainBoundingBox, class RangeBoundingBox>
00138 void box_global_bounds(
00139   ParallelMachine            arg_comm ,
00140   const size_t               arg_domain_boxes_number ,
00141   const DomainBoundingBox * const arg_domain_boxes ,
00142   const size_t               arg_range_boxes_number ,
00143   const RangeBoundingBox  * const arg_range_boxes ,
00144   float        * const arg_global_box )
00145 {
00146   enum { Dim = 3 };
00147 
00148   const bool symmetric = static_cast<const void * const>(arg_range_boxes) == static_cast<const void *const>(arg_domain_boxes ) ||
00149     arg_range_boxes == NULL;
00150 
00151 
00152   Copy<Dim>( arg_global_box ,         std::numeric_limits<float>::max() );
00153   Copy<Dim>( arg_global_box + Dim , - std::numeric_limits<float>::max() );
00154 
00155   //------------------------------------
00156   // Trivial loop threading possible:
00157 
00158   for ( size_t i = 0 ; i < arg_domain_boxes_number ; ++i ) {
00159     const DomainBoundingBox & box = arg_domain_boxes[i];
00160     for (int j=0; j<Dim; ++j) {
00161       if (arg_global_box[j] > static_cast<float>(box.lower(j)) ) {
00162         arg_global_box[j] = static_cast<float>(box.lower(j));
00163       }
00164       if (arg_global_box[j+Dim] < static_cast<float>(box.upper(j)) ) {
00165         arg_global_box[j+Dim] = static_cast<float>(box.upper(j));
00166       }
00167     }
00168   }
00169 
00170   if ( ! symmetric ) {
00171     for ( size_t i = 0 ; i < arg_range_boxes_number ; ++i ) {
00172       const RangeBoundingBox & box = arg_range_boxes[i];
00173       for (int j=0; j<Dim; ++j) {
00174         if (arg_global_box[j] > static_cast<float>(box.lower(j)) ) {
00175           arg_global_box[j] = static_cast<float>(box.lower(j));
00176         }
00177         if (arg_global_box[j+Dim] < static_cast<float>(box.upper(j)) ) {
00178           arg_global_box[j+Dim] = static_cast<float>(box.upper(j));
00179         }
00180       }
00181     }
00182   }
00183 
00184   //------------------------------------
00185 
00186   all_reduce( arg_comm , ReduceMin<Dim>( arg_global_box ) );
00187   all_reduce( arg_comm , ReduceMax<Dim>( arg_global_box + Dim ) );
00188 
00189   if (arg_domain_boxes_number == 0 &&
00190       arg_range_boxes_number  == 0)
00191     return;
00192 
00193   // Scale up and down by epsilon
00194 
00195   const double eps = std::numeric_limits<float>::epsilon();
00196 
00197   for ( unsigned i = 0 ; i < Dim ; ++i ) {
00198     float upper = arg_global_box[i+Dim] ;
00199     float lower = arg_global_box[i] ;
00200 
00201     double delta = eps * ( upper - lower );
00202     delta = delta > 0 ? delta : eps;
00203 
00204     while ( upper <= arg_global_box[i+Dim] ||
00205             arg_global_box[i] <= lower ) {
00206       upper = static_cast<float>( upper + delta );
00207       lower = static_cast<float>( lower - delta );
00208       delta *= 2 ;
00209     }
00210 
00211     arg_global_box[i+Dim] = static_cast<float>(upper);
00212     arg_global_box[i]     = static_cast<float>(lower);
00213   }
00214 }
00215 
00216 //----------------------------------------------------------------------
00217 
00218 template< class S >
00219 struct SetInsertBuffer {
00220   enum { N = 128 };
00221   S      & m_set ;
00222   unsigned m_lock ;
00223   unsigned m_iter ;
00224   typename S::value_type m_buffer[ N ] ;
00225 
00226   void overflow();
00227   void operator()( const typename S::value_type & v );
00228 
00229   SetInsertBuffer( S & s , unsigned l )
00230     : m_set(s), m_lock(l), m_iter(0) {}
00231 
00232   ~SetInsertBuffer() { overflow(); }
00233 };
00234 
00235 template<class S>
00236 void SetInsertBuffer<S>::overflow()
00237 {
00238   try {
00239     TPI_Lock( m_lock );
00240     while ( m_iter ) {
00241       m_set.insert( m_buffer[ --m_iter ] );
00242     }
00243     TPI_Unlock( m_lock );
00244   }
00245   catch(...) {
00246     TPI_Unlock( m_lock );
00247     throw ;
00248   }
00249 }
00250 
00251 template<class S>
00252 void SetInsertBuffer<S>::operator()( const typename S::value_type & v )
00253 {
00254   m_buffer[ m_iter ] = v ;
00255   if ( N == ++m_iter ) { overflow(); }
00256 }
00257 
00258 
00259 /*
00260 template <class DomainBoundingBox, class RangeBoundingBox>
00261 void proximity_search_symmetric(
00262   const typename std::map< stk::OctTreeKey,
00263                            std::pair< std::list< DomainBoundingBox >,
00264                                       std::list< RangeBoundingBox > > >::const_iterator i_beg ,
00265   const typename std::map< stk::OctTreeKey,
00266                            std::pair< std::list< DomainBoundingBox >,
00267                                       std::list< RangeBoundingBox > > >::const_iterator i_end ,
00268   SetInsertBuffer< std::set< std::pair< typename DomainBoundingBox::Key,
00269                                         typename RangeBoundingBox::Key > > > & arg_out )
00270 {
00271   typedef typename DomainBoundingBox::Key DomainKey;
00272   typedef typename RangeBoundingBox::Key RangeKey;
00273   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00274 
00275   typename SearchTree::const_iterator j ;
00276 
00277   typename std::list<DomainBoundingBox>::const_iterator id;
00278   typename std::list<RangeBoundingBox>::const_iterator ir;
00279 
00280   const typename SearchTree::value_type  & inode = *i_beg ;
00281   const std::list<DomainBoundingBox> & domain_outer = inode.second.first ;
00282 
00283   const typename std::list<DomainBoundingBox>::const_iterator
00284     beg_dom_out = domain_outer.begin(),
00285     end_dom_out = domain_outer.end();
00286 
00287   // Outer cell vs. itself
00288 
00289   for ( id = beg_dom_out ; id != end_dom_out ; ++id ) {
00290     const DomainBoundingBox & d = *id ;
00291     for ( ir = id ; ++ir != end_dom_out ; ) {
00292       const RangeBoundingBox & r = *ir ;
00293       if ( d.intersect(r) ) {
00294         const DomainKey & dip = d.key ;
00295         const RangeKey & rip = r.key ;
00296         if ( dip < rip ) {
00297           std::pair<DomainKey,RangeKey> tmp( dip , rip );
00298           arg_out( tmp );
00299         }
00300         else {
00301           std::pair<Key,Key> tmp( rip , dip );
00302           arg_out( tmp );
00303         }
00304       }
00305     }
00306   }
00307 
00308   // Outer cell searching inner cells.
00309   // Outer cell always precedes inner cells
00310   // Iterate forward until the cell is not contained.
00311 
00312   const stk::OctTreeKey & outer_key = inode.first ;
00313 
00314   for ( j = i_beg ; ++j != i_end && outer_key.intersect( (*j).first ) ; ) {
00315 
00316     const typename SearchTree::value_type & jnode = *j ;
00317 
00318     const std::list<BoundingBox> & domain_inner = jnode.second.first ;
00319 
00320     const typename std::list<BoundingBox>::const_iterator
00321       beg_dom_inn = domain_inner.begin(),
00322       end_dom_inn = domain_inner.end();
00323 
00324     // Check domain_outer vs. domain_inner,
00325     // skip if the same box.
00326 
00327     for ( id = beg_dom_out ; id != end_dom_out ; ++id ) {
00328       const BoundingBox & d = *id ;
00329       for ( ir = beg_dom_inn ; ir != end_dom_inn ; ++ir ) {
00330         const BoundingBox & r = *ir ;
00331         if ( d.key != r.key ) {
00332           if ( d.intersect(r) ) {
00333             const Key & dip = d.key ;
00334             const Key & rip = r.key ;
00335             if ( dip < rip ) {
00336               std::pair<Key,Key> tmp( dip , rip );
00337               arg_out( tmp );
00338             }
00339             else {
00340               std::pair<Key,Key> tmp( rip , dip );
00341               arg_out( tmp );
00342             }
00343           }
00344         }
00345       }
00346     }
00347   }
00348 }
00349 */
00350 
00351 template <class DomainBoundingBox, class RangeBoundingBox>
00352 void proximity_search_asymmetric(
00353   const typename std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > >::const_iterator i_beg ,
00354   const typename std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > >::const_iterator i_end ,
00355   SetInsertBuffer< std::set< std::pair< typename DomainBoundingBox::Key,  typename RangeBoundingBox::Key > > > & arg_out )
00356 {
00357   typedef typename DomainBoundingBox::Key DomainKey;
00358   typedef typename RangeBoundingBox::Key RangeKey;
00359   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00360 
00361   typename SearchTree::const_iterator j ;
00362 
00363   typename std::list<DomainBoundingBox>::const_iterator id;
00364   typename std::list<RangeBoundingBox>::const_iterator ir;
00365 
00366   const typename SearchTree::value_type & inode = *i_beg ;
00367 
00368   const std::list<DomainBoundingBox> & domain_outer = inode.second.first ;
00369   const std::list<RangeBoundingBox> & range_outer  = inode.second.second ;
00370 
00371   const typename std::list<DomainBoundingBox>::const_iterator
00372     beg_dom_out = domain_outer.begin(),
00373     end_dom_out = domain_outer.end();
00374   const typename std::list<RangeBoundingBox>::const_iterator
00375     beg_ran_out = range_outer.begin(),
00376     end_ran_out = range_outer.end();
00377 
00378   // domain_outer vs. range_outer
00379 
00380   for ( id = beg_dom_out ; id != end_dom_out ; ++id ) {
00381     const DomainBoundingBox & d = *id ;
00382     for ( ir = beg_ran_out ; ir != end_ran_out ; ++ir ) {
00383       const RangeBoundingBox & r = *ir ;
00384       if ( d.intersect(r) ) {
00385         std::pair<DomainKey,RangeKey> tmp( d.key , r.key );
00386         arg_out( tmp );
00387       }
00388     }
00389   }
00390 
00391   // Outer cell searching inner cells.
00392   // Outer cell always precedes inner cells
00393   // Iterate forward until the cell is not contained.
00394 
00395   const stk::OctTreeKey & outer_key = inode.first ;
00396 
00397   for ( j = i_beg ; ++j != i_end && outer_key.intersect( (*j).first ) ; ) {
00398 
00399     const typename SearchTree::value_type & jnode = *j ;
00400 
00401     const std::list<RangeBoundingBox> & range_inner  = jnode.second.second ;
00402 
00403     const typename std::list<RangeBoundingBox>::const_iterator
00404       beg_ran_inn = range_inner.begin(),
00405       end_ran_inn = range_inner.end();
00406 
00407     // Check domain_outer vs. range_inner
00408 
00409     for ( ir = beg_ran_inn ; ir != end_ran_inn ; ++ir ) {
00410       const RangeBoundingBox & r = *ir ;
00411       for ( id = beg_dom_out ; id != end_dom_out ; ++id ) {
00412         const DomainBoundingBox & d = *id ;
00413         if ( d.intersect(r) ) {
00414           std::pair<DomainKey,RangeKey> tmp( d.key , r.key );
00415           arg_out( tmp );
00416         }
00417       }
00418     }
00419 
00420     // Check domain_inner vs. range_outer if non-symmetric
00421     const std::list<DomainBoundingBox> & domain_inner = jnode.second.first ;
00422 
00423     const typename std::list<DomainBoundingBox>::const_iterator
00424       beg_dom_inn = domain_inner.begin(),
00425       end_dom_inn = domain_inner.end();
00426 
00427     for ( id = beg_dom_inn ; id != end_dom_inn ; ++id ) {
00428       const DomainBoundingBox & d = *id ;
00429       for ( ir = beg_ran_out ; ir != end_ran_out ; ++ir ) {
00430         const RangeBoundingBox & r = *ir ;
00431         if ( d.intersect(r) ) {
00432           std::pair<DomainKey,RangeKey> tmp( d.key , r.key );
00433           arg_out( tmp );
00434         }
00435       }
00436     }
00437   }
00438 }
00439 
00440 unsigned processor( const stk::OctTreeKey * const cuts_b ,
00441                     const stk::OctTreeKey * const cuts_e ,
00442                     const stk::OctTreeKey & key );
00443 
00444 template <class DomainBoundingBox, class RangeBoundingBox>
00445 void pack(
00446   CommAll & comm_all ,
00447   const stk::OctTreeKey * const cuts_b ,
00448   const std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > & send_tree ,
00449         std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > * recv_tree )
00450 {
00451   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00452 
00453   const unsigned p_rank = comm_all.parallel_rank();
00454   const unsigned p_size = comm_all.parallel_size();
00455   const stk::OctTreeKey * const cuts_e = cuts_b + p_size ;
00456 
00457   typename SearchTree::const_iterator i ;
00458 
00459   for ( i = send_tree.begin() ; i != send_tree.end() ; ++i ) {
00460     const stk::OctTreeKey & key = (*i).first ;
00461 
00462     unsigned p = processor( cuts_b , cuts_e , key );
00463 
00464     do {
00465       if ( p != p_rank ) {
00466         CommBuffer & buf = comm_all.send_buffer(p);
00467 
00468         const std::list< DomainBoundingBox > & domain = (*i).second.first ;
00469         const std::list< RangeBoundingBox > & range  = (*i).second.second ;
00470 
00471         typename std::list< DomainBoundingBox >::const_iterator jd ;
00472         typename std::list< RangeBoundingBox >::const_iterator jr ;
00473 
00474         const unsigned dsize = domain.size();
00475         const unsigned rsize = range.size();
00476 
00477         buf.pack<unsigned>( key.value() , stk::OctTreeKey::NWord );
00478         buf.pack<unsigned>( dsize );
00479         buf.pack<unsigned>( rsize );
00480 
00481         for ( jd = domain.begin() ; jd != domain.end() ; ++jd ) {
00482           const DomainBoundingBox & box = *jd ;
00483           buf.pack<DomainBoundingBox>( box );
00484         }
00485 
00486         for ( jr = range.begin() ; jr != range.end() ; ++jr ) {
00487           const RangeBoundingBox & box = *jr ;
00488           buf.pack<RangeBoundingBox>( box );
00489         }
00490       }
00491       else if ( recv_tree ) {
00492         // Copy contents of the send node
00493         (*recv_tree)[ key ] = (*i).second ;
00494       }
00495 
00496       // If the cut keys are at a finer granularity than
00497       // this key then this key may overlap more than one
00498       // processor's span.  Check for overlap with the
00499       // beginning key of the next processor.
00500 
00501       ++p ;
00502 
00503     } while( p < p_size && key.intersect( cuts_b[p] ) );
00504   }
00505 }
00506 
00507 template <class DomainBoundingBox, class RangeBoundingBox>
00508 void unpack(
00509   CommAll & comm_all ,
00510   std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > & tree )
00511 {
00512   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00513 
00514   unsigned domain_size(0) ;
00515   unsigned range_size(0) ;
00516   unsigned value[ stk::OctTreeKey::NWord ];
00517   stk::OctTreeKey key ;
00518   DomainBoundingBox domain_box ;
00519   RangeBoundingBox range_box ;
00520 
00521   const unsigned p_size = comm_all.parallel_size();
00522 
00523   for ( unsigned p = 0 ; p < p_size ; ++p ) {
00524     CommBuffer & buf = comm_all.recv_buffer(p);
00525 
00526     while ( buf.remaining() ) {
00527       buf.unpack<unsigned>( value , stk::OctTreeKey::NWord );
00528       buf.unpack<unsigned>( domain_size );
00529       buf.unpack<unsigned>( range_size );
00530 
00531       // Insert key, get domain and range
00532 
00533       key.set_value( value );
00534 
00535       typename SearchTree::mapped_type & node = tree[ key ];
00536 
00537       std::list< DomainBoundingBox > & domain = node.first ;
00538       std::list< RangeBoundingBox > & range  = node.second ;
00539 
00540       for ( unsigned j = 0 ; j < domain_size ; ++j ) {
00541         buf.unpack<DomainBoundingBox>( domain_box );
00542         domain.push_back( domain_box );
00543       }
00544 
00545       for ( unsigned j = 0 ; j < range_size ; ++j ) {
00546         buf.unpack<RangeBoundingBox>( range_box );
00547         range.push_back( range_box );
00548       }
00549     }
00550   }
00551 }
00552 
00553 template <class DomainBoundingBox, class RangeBoundingBox>
00554 bool communicate(
00555   stk::ParallelMachine arg_comm ,
00556   const stk::OctTreeKey * const arg_cuts ,
00557   const std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > & send_tree ,
00558         std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > & recv_tree ,
00559   const bool local_flag )
00560 {
00561   const unsigned p_size = parallel_machine_size( arg_comm );
00562 
00563   // Communicate search_tree members
00564 
00565   CommAll comm_all( arg_comm );
00566 
00567   // Sizing pass for pack
00568   pack<DomainBoundingBox, RangeBoundingBox>( comm_all , arg_cuts , send_tree , NULL );
00569 
00570   // If more than 25% then is dense
00571   const bool global_flag =
00572     comm_all.allocate_buffers( p_size / 4 , false , local_flag );
00573 
00574   // Actual packing pass, copy local entries too
00575   pack<DomainBoundingBox, RangeBoundingBox>( comm_all , arg_cuts , send_tree , & recv_tree );
00576 
00577   comm_all.communicate();
00578 
00579   unpack<DomainBoundingBox, RangeBoundingBox>( comm_all , recv_tree );
00580 
00581   return global_flag ;
00582 }
00583 
00584 
00585 template <class DomainBoundingBox, class RangeBoundingBox>
00586 void communicate(
00587   stk::ParallelMachine arg_comm ,
00588   const std::set< std::pair< typename DomainBoundingBox::Key,  typename RangeBoundingBox::Key > > & send_relation ,
00589         std::set< std::pair< typename DomainBoundingBox::Key,  typename RangeBoundingBox::Key > > & recv_relation )
00590 {
00591   typedef typename DomainBoundingBox::Key DomainKey;
00592   typedef typename RangeBoundingBox::Key RangeKey;
00593   typedef std::pair<DomainKey, RangeKey> ValueType ;
00594 
00595   CommAll comm_all( arg_comm );
00596 
00597   const unsigned p_rank = comm_all.parallel_rank();
00598   const unsigned p_size = comm_all.parallel_size();
00599 
00600   typename std::set< ValueType >::const_iterator i ;
00601 
00602   for ( i = send_relation.begin() ; i != send_relation.end() ; ++i ) {
00603     const ValueType & val = *i ;
00604     if ( val.first.proc == p_rank || val.second.proc == p_rank ) {
00605       recv_relation.insert( val );
00606     }
00607     if ( val.first.proc != p_rank ) {
00608       CommBuffer & buf = comm_all.send_buffer( val.first.proc );
00609       buf.skip<ValueType>( 1 );
00610     }
00611     if ( val.second.proc != p_rank && val.second.proc != val.first.proc ) {
00612       CommBuffer & buf = comm_all.send_buffer( val.second.proc );
00613       buf.skip<ValueType>( 1 );
00614     }
00615   }
00616 
00617   // If more than 25% messages then is dense
00618 
00619   comm_all.allocate_buffers( p_size / 4 , false );
00620 
00621   for ( i = send_relation.begin() ; i != send_relation.end() ; ++i ) {
00622     const ValueType & val = *i ;
00623     if ( val.first.proc != p_rank ) {
00624       CommBuffer & buf = comm_all.send_buffer( val.first.proc );
00625       buf.pack<ValueType>( val );
00626     }
00627     if ( val.second.proc != p_rank && val.second.proc != val.first.proc ) {
00628       CommBuffer & buf = comm_all.send_buffer( val.second.proc );
00629       buf.pack<ValueType>( val );
00630     }
00631   }
00632 
00633   comm_all.communicate();
00634 
00635   for ( unsigned p = 0 ; p < p_size ; ++p ) {
00636     CommBuffer & buf = comm_all.recv_buffer( p );
00637     while ( buf.remaining() ) {
00638       ValueType val ;
00639       buf.unpack<ValueType>( val );
00640       recv_relation.insert( val );
00641     }
00642   }
00643 }
00644 
00645 //----------------------------------------------------------------------
00646 // Partition a search tree among processors
00647 
00648 template <class DomainBoundingBox, class RangeBoundingBox>
00649 void oct_tree_partition(
00650   stk::ParallelMachine        arg_comm ,
00651   const std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > & arg_tree ,
00652   const double                arg_tolerance ,
00653   std::vector< stk::OctTreeKey > & arg_cuts )
00654 {
00655   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00656 
00657   enum { tree_depth  = 4 };
00658   enum { tree_size   = OctTreeSize< tree_depth >::value };
00659   enum { tree_size_2 = tree_size * 2 };
00660 
00661   const unsigned p_size = parallel_machine_size( arg_comm );
00662   const stk::OctTreeKey k_null ;
00663 
00664   arg_cuts.assign( p_size , k_null );
00665 
00666   float local_count[  tree_size_2 ];
00667   float global_count[ tree_size_2 ];
00668 
00669   for ( unsigned i = 0 ; i < tree_size_2 ; ++i ) {
00670     local_count[i] = 0.0 ;
00671   }
00672 
00673   for ( typename SearchTree::const_iterator i =  arg_tree.begin() ;
00674                                    i != arg_tree.end() ; ++i ) {
00675 
00676     const stk::OctTreeKey & key = (*i).first ;
00677 
00678     const std::list< DomainBoundingBox > & domain = (*i).second.first ;
00679     const std::list< RangeBoundingBox > & range  = (*i).second.second ;
00680 
00681     const unsigned depth   = key.depth();
00682     const unsigned ordinal = oct_tree_offset( tree_depth , key );
00683     const unsigned num_d   = domain.size();
00684     const unsigned num_r   = range.size();
00685     const unsigned number  = num_d + num_r ;
00686 
00687     if ( depth <= 4 ) { // Values for this node:
00688       local_count[ 2 * ordinal ] += number ;
00689     }
00690     else { // Values for a deeper node
00691       local_count[ 2 * ordinal + 1 ] += number ;
00692     }
00693   }
00694 
00695   all_reduce_sum( arg_comm , local_count , global_count , tree_size_2 );
00696 
00697   oct_tree_partition_private( 0, p_size, tree_depth,
00698                               arg_tolerance, global_count,
00699                               p_size, & arg_cuts[0]);
00700 }
00701 
00702 template <class DomainBoundingBox, class RangeBoundingBox>
00703 class ProximitySearch {
00704 public:
00705   typedef typename DomainBoundingBox::Key DomainKey;
00706   typedef typename RangeBoundingBox::Key RangeKey;
00707   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00708   typedef void (*proximity_search_work_routine)(
00709     const typename SearchTree::const_iterator i_beg ,
00710     const typename SearchTree::const_iterator i_end ,
00711     SetInsertBuffer< std::set< std::pair<DomainKey, RangeKey> > > & );
00712 
00713   enum { NLOCKS = 2 };
00714   enum { GET_LOCK = 0 };
00715   enum { PUT_LOCK = 1 };
00716 
00717   bool m_symmetric;
00718 
00719   std::set< std::pair<DomainKey, RangeKey> > & m_relation ;
00720   typename SearchTree::const_iterator m_tree_iter ;
00721   typename SearchTree::const_iterator m_tree_end ;
00722 
00723   ~ProximitySearch() {}
00724 
00725   ProximitySearch(
00726     bool symmetric ,
00727     const SearchTree & search_tree ,
00728     std::set< std::pair<DomainKey, RangeKey> > & relation );
00729   void iterate_tree();
00730 
00731 private:
00732   ProximitySearch();
00733   ProximitySearch( const ProximitySearch & );
00734   ProximitySearch & operator = ( const ProximitySearch & );
00735 };
00736 
00737 template <class DomainBoundingBox, class RangeBoundingBox>
00738 void proximity_search_work( TPI_Work * work )
00739 {
00740   ProximitySearch<DomainBoundingBox, RangeBoundingBox> * p = (ProximitySearch<DomainBoundingBox, RangeBoundingBox> *) work->info ;
00741   p->iterate_tree();
00742 }
00743 
00744 
00745 template <class DomainBoundingBox, class RangeBoundingBox>
00746 void ProximitySearch<DomainBoundingBox, RangeBoundingBox>::iterate_tree()
00747 {
00748   enum { N_WORK = 32 };
00749 
00750   try {
00751     SetInsertBuffer< std::set< std::pair<DomainKey, RangeKey> > >
00752       tmp( m_relation , PUT_LOCK );
00753 
00754     const typename SearchTree::const_iterator i_tree_end = m_tree_end ;
00755 
00756     unsigned n_work = N_WORK ;
00757 
00758     while ( n_work ) {
00759 
00760       n_work = 0 ;
00761 
00762       typename SearchTree::const_iterator i_beg , i_end ;
00763 
00764       // Get work:
00765 
00766       try {
00767         TPI_Lock( GET_LOCK );
00768 
00769         i_end = i_beg = m_tree_iter ;
00770 
00771         while ( n_work < N_WORK && i_tree_end != i_end ) {
00772           ++i_end ; ++n_work ;
00773         }
00774 
00775         m_tree_iter = i_end ;
00776 
00777         TPI_Unlock( GET_LOCK );
00778       }
00779       catch( ... ) {
00780         TPI_Unlock( GET_LOCK );
00781   throw ;
00782       }
00783 
00784       // Perform work:
00785 
00786     //  if (m_symmetric) {
00787     //    for ( ; i_beg != i_end ; ++i_beg ) {
00788     //      proximity_search_symmetric<DomainBoundingBox,RangeBoundingBox>( i_beg , i_tree_end , tmp );
00789     //    }
00790     //  } else {
00791         for ( ; i_beg != i_end ; ++i_beg ) {
00792           proximity_search_asymmetric<DomainBoundingBox, RangeBoundingBox>( i_beg , i_tree_end , tmp );
00793         }
00794     //  }
00795     }
00796   }
00797   catch ( const std::exception & x ) {
00798     std::cerr << x.what() << std::endl ;
00799     std::cerr.flush();
00800   }
00801   catch ( ... ) {
00802     std::cerr << "ProximitySearch::iterate_tree FAILED" << std::endl ;
00803     std::cerr.flush();
00804   }
00805 }
00806 
00807 
00808 template <class DomainBoundingBox, class RangeBoundingBox>
00809 ProximitySearch<DomainBoundingBox, RangeBoundingBox>::ProximitySearch(
00810   bool                  symmetric ,
00811   const SearchTree &    search_tree ,
00812   std::set< std::pair<DomainKey, RangeKey> > & relation )
00813   : m_symmetric(symmetric),
00814     m_relation( relation ),
00815     m_tree_iter( search_tree.begin() ),
00816     m_tree_end(  search_tree.end() )
00817 {
00818   if ( m_tree_iter != m_tree_end ) {
00819 
00820     TPI_work_subprogram worker = (TPI_work_subprogram) proximity_search_work<DomainBoundingBox, RangeBoundingBox>;
00821     TPI_Run_threads(worker, this, NLOCKS );
00822 
00823     if ( m_tree_iter != m_tree_end ) {
00824       std::string msg("stk::proximity_search FAILED to complete" );
00825       throw std::runtime_error(msg);
00826     }
00827   }
00828 }
00829 
00830 
00831 //----------------------------------------------------------------------
00855 template <class DomainBoundingBox, class RangeBoundingBox>
00856 bool oct_tree_proximity_search(
00857   ParallelMachine            arg_comm ,
00858   const float        * const arg_global_box ,
00859   const size_t               arg_domain_boxes_number ,
00860   const DomainBoundingBox * const arg_domain_boxes ,
00861   const size_t               arg_range_boxes_number ,
00862   const RangeBoundingBox * const arg_range_boxes ,
00863   const OctTreeKey   * const arg_cuts ,
00864   std::vector< std::pair< typename DomainBoundingBox::Key,  typename RangeBoundingBox::Key > > & arg_relation ,
00865   unsigned * const arg_search_tree_stats = NULL )
00866 {
00867   typedef typename DomainBoundingBox::Key DomainKey;
00868   typedef typename RangeBoundingBox::Key RangeKey;
00869   typedef std::map< stk::OctTreeKey, std::pair< std::list< DomainBoundingBox >, std::list< RangeBoundingBox > > > SearchTree ;
00870 
00871   enum { Dim = 3 };
00872 
00873   //const bool symmetric = static_cast<const void * const>(arg_range_boxes) == static_cast<const void *const>(arg_domain_boxes ) ||
00874   //  arg_range_boxes == NULL;
00875   const bool symmetric = false;
00876 
00877   const unsigned p_size = parallel_machine_size( arg_comm );
00878   const unsigned p_rank = parallel_machine_rank( arg_comm );
00879 
00880   //----------------------------------------------------------------------
00881   // Search tree defined by oct-tree covering for boxes
00882 
00883   bool local_violations = false ;
00884   bool global_violations = false ;
00885 
00886   SearchTree search_tree ;
00887 
00888   {
00889     stk::OctTreeKey covering[8] ;
00890     unsigned   number = 0 ;
00891 
00892     double scale = arg_global_box[0+Dim] - arg_global_box[0];
00893     for ( unsigned i = 1 ; i < Dim ; ++i ) {
00894       double tst_scale = arg_global_box[i+Dim] - arg_global_box[i];
00895       if (tst_scale > scale) scale = tst_scale;
00896     }
00897     if (scale > 0.0) // Not an error. Could arise with a point bounding box in the range/domain...
00898       scale = 1.0 / scale;
00899     else
00900       scale = 1.0;
00901 
00902     for ( size_t i = 0 ; i < arg_domain_boxes_number ; ++i ) {
00903 
00904       DomainBoundingBox tmp( arg_domain_boxes[i] );
00905 
00906       tmp.key.proc = p_rank ;
00907 
00908       float box[6];
00909 
00910       for (int x=0; x<Dim; ++x) {
00911         box[x] = tmp.lower(x);
00912         box[x+Dim] = tmp.upper(x);
00913       }
00914 
00915       const bool valid =
00916         hsfc_box_covering( arg_global_box, box, covering, number, scale );
00917 
00918       if ( ! valid ) { local_violations = true ; }
00919 
00920       for ( unsigned k = 0 ; k < number ; ++k ) {
00921         const stk::OctTreeKey key = covering[k] ;
00922         search_tree[key].first.push_back(tmp);
00923       }
00924     }
00925 
00926     if ( ! symmetric ) {
00927       for ( size_t i = 0 ; i < arg_range_boxes_number ; ++i ) {
00928 
00929         RangeBoundingBox tmp( arg_range_boxes[i] );
00930 
00931         tmp.key.proc = p_rank ;
00932         float box[6];
00933 
00934         for (int x=0; x<Dim; ++x) {
00935           box[x] = tmp.lower(x);
00936           box[x+Dim] = tmp.upper(x);
00937         }
00938 
00939         const bool valid =
00940           hsfc_box_covering( arg_global_box, box, covering, number, scale );
00941 
00942         if ( ! valid ) { local_violations = true ; }
00943 
00944         for ( unsigned k = 0 ; k < number ; ++k ) {
00945           const stk::OctTreeKey key = covering[k] ;
00946           search_tree[key].second.push_back(tmp);
00947         }
00948       }
00949     }
00950   }
00951 
00952   //----------------------------------------------------------------------
00953   // Use a set to provide a unique and sorted result.
00954 
00955   std::set< std::pair<DomainKey, RangeKey> > tmp_relation ;
00956 
00957   if ( p_size == 1 ) {
00958 
00959     global_violations = local_violations ;
00960 
00961     if ( arg_search_tree_stats ) {
00962       search_tree_statistics( arg_comm , search_tree ,
00963                               arg_search_tree_stats );
00964     }
00965     ProximitySearch<DomainBoundingBox, RangeBoundingBox>( symmetric, search_tree, tmp_relation);
00966   }
00967   else {
00968     // Communicate search_tree members
00969 
00970     SearchTree local_tree ;
00971 
00972     std::set< std::pair<DomainKey, RangeKey> > local_relation ;
00973 
00974     if ( arg_cuts ) {
00975       global_violations =
00976         communicate<DomainBoundingBox, RangeBoundingBox>( arg_comm , arg_cuts , search_tree , local_tree ,
00977                            local_violations );
00978     }
00979     else {
00980       const double tolerance = 0.001 ;
00981 
00982       std::vector< stk::OctTreeKey > cuts ;
00983 
00984       oct_tree_partition( arg_comm , search_tree , tolerance , cuts );
00985 
00986       global_violations =
00987         communicate<DomainBoundingBox, RangeBoundingBox>(arg_comm , & cuts[0] , search_tree , local_tree ,
00988                           local_violations );
00989     }
00990 
00991     // Local proximity search with received members
00992 
00993     if ( arg_search_tree_stats ) {
00994       search_tree_statistics( arg_comm , local_tree ,
00995                               arg_search_tree_stats );
00996     }
00997 
00998     ProximitySearch<DomainBoundingBox, RangeBoundingBox>( symmetric, local_tree, local_relation);
00999 
01000     // Communicate relations back to domain and range processors
01001 
01002     communicate<DomainBoundingBox, RangeBoundingBox>( arg_comm , local_relation , tmp_relation );
01003   }
01004 
01005   arg_relation.clear();
01006   arg_relation.reserve( tmp_relation.size() );
01007 
01008   typename std::set< std::pair<DomainKey, RangeKey> >::iterator ir ;
01009   for ( ir = tmp_relation.begin() ; ir != tmp_relation.end() ; ++ir ) {
01010     arg_relation.push_back( *ir );
01011   }
01012   return global_violations ;
01013 }
01014 
01015 //----------------------------------------------------------------------
01016 
01017 } // namespace search
01018 } // namespace stk
01019 
01020 #endif
01021 

Generated on Tue Jul 13 09:27:32 2010 for Sierra Toolkit by  doxygen 1.4.7