gs.cpp

Go to the documentation of this file.

/* compile-time settings:
 
 FORTRAN naming convention
 default cpgs_setup, etc.
 -DUPCASE CPGS_SETUP, etc.
 -DUNDERSCORE cpgs_setup_, etc.
 
 -DMPI parallel version (sequential otherwise)
 -DCRYSTAL_STATIC avoid some message exchange at the risk of
 crashing b/c of insufficient buffer size
 
 -DINITIAL_BUFFER_SIZE=expression
 ignored unless CRYSTAL_STATIC is defined.
 arithmetic expression controlling the initial buffer size for the crystal
 router; this needs to be large enough to hold the intermediate messages
 during all stages of the crystal router
 
 variables that can be used in expression include
 num - the number of processors
 n - the length of the global index array
 
 */
 
/* default for INITIAL_BUFFER_SIZE */
#ifdef CRYSTAL_STATIC
#ifndef INITIAL_BUFFER_SIZE
#define INITIAL_BUFFER_SIZE 2 * ( 3 * num + n * 9 )
#endif
#endif
 
/* FORTRAN usage:
 
 integer hc, np
 call crystal_new(hc,comm,np) ! get a crystal router handle (see fcrystal.c)
 
 integer hgs
 integer n, max_vec_dim
 integer*? global_index_array(1:n) ! type corresponding to slong in "types.h"
 
 call cpgs_setup(hgs,hc,global_index_array,n,max_vec_dim)
 sets hgs to new handle
 
 !ok to call crystal_done(hc) here, or any later time
 
 call cpgs_op(hgs, u, op)
 integer handle, op : 1-add, 2-multiply, 3-min, 4-max
 real u(1:n) - same layout as global_index_array provided to cpgs_setup
 
 call cpgs_op_vec(hgs, u, d, op)
 integer op : 1-add, 2-multiply, 3-min, 4-max
 integer d <= max_vec_dim
 real u(1:d, 1:n) - vector components for each node stored together
 
 call cpgs_op_many(hgs, u1, u2, u3, u4, u5, u6, d, op)
 integer op : 1-add, 2-multiply, 3-min, 4-max
 integer d : in {1,2,3,4,5,6}, <= max_vec_dim
 real u1(1:n), u2(1:n), u3(1:n), etc.
 
 same effect as: call cpgs_op(hgs, u1, op)
 if(d.gt.1) call cpgs_op(hgs, u2, op)
 if(d.gt.2) call cpgs_op(hgs, u3, op)
 etc.
 with possibly some savings as fewer messages are exchanged
 
 call cpgs_free(hgs)
 */
 
#include <cstdio>
#include <cstdlib>
#include <cstdarg>
#include <cstring>
#include <cmath>
#include "moab/gs.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab_mpi.h"
#endif
 
#ifdef MOAB_HAVE_MPI
 
#ifdef MOAB_HAVE_VALGRIND
#include <valgrind/memcheck.h>
#elif !defined( VALGRIND_CHECK_MEM_IS_DEFINED )
#define VALGRIND_CHECK_MEM_IS_DEFINED( a, b ) ( (void)0 )
#endif
 
#endif
 
namespace moab
{
 
/*--------------------------------------------------------------------------
 Local Execution Phases
 --------------------------------------------------------------------------*/
 
#define DO_SET( a, b ) b = a
#define DO_ADD( a, b ) a += b
#define DO_MUL( a, b ) a *= b
#define DO_MIN( a, b ) \
 if( ( b ) < ( a ) ) ( a ) = b
#define DO_MAX( a, b ) \
 if( ( b ) > ( a ) ) ( a ) = b
#define DO_BPR( a, b ) \
 do \
 { \
 uint a_ = a; \
 uint b_ = b; \
 for( ;; ) \
 { \
 if( a_ < b_ ) \
 b_ >>= 1; \
 else if( b_ < a_ ) \
 a_ >>= 1; \
 else \
 break; \
 } \
 ( a ) = a_; \
 } while( 0 )
 
#define LOOP( op ) \
 do \
 { \
 sint i, j; \
 while( ( i = *cm++ ) != -1 ) \
 while( ( j = *cm++ ) != -1 ) \
 op( u[i], u[j] ); \
 } while( 0 )
 
static void local_condense( realType* u, int op, const sint* cm )
{
 switch( op )
 {
 case GS_OP_ADD:
 LOOP( DO_ADD );
 break;
 case GS_OP_MUL:
 LOOP( DO_MUL );
 break;
 case GS_OP_MIN:
 LOOP( DO_MIN );
 break;
 case GS_OP_MAX:
 LOOP( DO_MAX );
 break;
 case GS_OP_BPR:
 LOOP( DO_BPR );
 break;
 }
}
 
static void local_uncondense( realType* u, const sint* cm )
{
 LOOP( DO_SET );
}
 
#undef LOOP
 
#define LOOP( op ) \
 do \
 { \
 sint i, j, k; \
 while( ( i = *cm++ ) != -1 ) \
 { \
 realType* pi = u + n * i; \
 while( ( j = *cm++ ) != -1 ) \
 { \
 realType* pj = u + n * j; \
 for( k = n; k; --k ) \
 { \
 op( *pi, *pj ); \
 ++pi; \
 ++pj; \
 } \
 } \
 } \
 } while( 0 )
 
static void local_condense_vec( realType* u, uint n, int op, const sint* cm )
{
 switch( op )
 {
 case GS_OP_ADD:
 LOOP( DO_ADD );
 break;
 case GS_OP_MUL:
 LOOP( DO_MUL );
 break;
 case GS_OP_MIN:
 LOOP( DO_MIN );
 break;
 case GS_OP_MAX:
 LOOP( DO_MAX );
 break;
 case GS_OP_BPR:
 LOOP( DO_BPR );
 break;
 }
}
 
static void local_uncondense_vec( realType* u, uint n, const sint* cm )
{
 LOOP( DO_SET );
}
 
#undef LOOP
 
/*--------------------------------------------------------------------------
 Non-local Execution Phases
 --------------------------------------------------------------------------*/
 
#ifdef MOAB_HAVE_MPI
 
void gs_data::nonlocal_info::initialize( uint np, uint count, uint nlabels, uint nulabels, uint maxv )
{
 _target = NULL;
 _nshared = NULL;
 _sh_ind = NULL;
 _slabels = NULL;
 _ulabels = NULL;
 _reqs = NULL;
 _buf = NULL;
 _np = np;
 _target = (uint*)malloc( ( 2 * np + count ) * sizeof( uint ) );
 _nshared = _target + np;
 _sh_ind = _nshared + np;
 if( 1 < nlabels )
 _slabels = (slong*)malloc( ( ( nlabels - 1 ) * count ) * sizeof( slong ) );
 else
 _slabels = NULL;
 _ulabels = (Ulong*)malloc( ( nulabels * count ) * sizeof( Ulong ) );
 _reqs = (MPI_Request*)malloc( 2 * np * sizeof( MPI_Request ) );
 _buf = (realType*)malloc( ( 2 * count * maxv ) * sizeof( realType ) );
 _maxv = maxv;
}
 
void gs_data::nonlocal_info::nlinfo_free()
{
 // Free the ptrs
 free( _buf );
 free( _reqs );
 free( _target );
 free( _slabels );
 free( _ulabels );
 // Set to null
 _ulabels = NULL;
 _buf = NULL;
 _reqs = NULL;
 _target = NULL;
 _slabels = NULL;
 _nshared = NULL;
 _sh_ind = NULL;
}
 
void gs_data::nonlocal_info::nonlocal( realType* u, int op, MPI_Comm comm )
{
 MPI_Status status;
 uint np = this->_np;
 MPI_Request* reqs = this->_reqs;
 uint* targ = this->_target;
 uint* nshared = this->_nshared;
 uint* sh_ind = this->_sh_ind;
 uint id;
 realType *buf = this->_buf, *start;
 unsigned int i;
 {
 MPI_Comm_rank( comm, (int*)&i );
 id = i;
 }
 for( i = 0; i < np; ++i )
 {
 uint c = nshared[i];
 start = buf;
 for( ; c; --c )
 *buf++ = u[*sh_ind++];
 MPI_Isend( (void*)start, nshared[i] * sizeof( realType ), MPI_UNSIGNED_CHAR, targ[i], id, comm, reqs++ );
 }
 start = buf;
 for( i = 0; i < np; ++i )
 {
 MPI_Irecv( (void*)start, nshared[i] * sizeof( realType ), MPI_UNSIGNED_CHAR, targ[i], targ[i], comm, reqs++ );
 start += nshared[i];
 }
 for( reqs = this->_reqs, i = np * 2; i; --i )
 MPI_Wait( reqs++, &status );
 sh_ind = this->_sh_ind;
#define LOOP( OP ) \
 do \
 { \
 for( i = 0; i < np; ++i ) \
 { \
 uint c; \
 for( c = nshared[i]; c; --c ) \
 { \
 OP( u[*sh_ind], *buf ); \
 ++sh_ind, ++buf; \
 } \
 } \
 } while( 0 )
 switch( op )
 {
 case GS_OP_ADD:
 LOOP( DO_ADD );
 break;
 case GS_OP_MUL:
 LOOP( DO_MUL );
 break;
 case GS_OP_MIN:
 LOOP( DO_MIN );
 break;
 case GS_OP_MAX:
 LOOP( DO_MAX );
 break;
 case GS_OP_BPR:
 LOOP( DO_BPR );
 break;
 }
#undef LOOP
}
 
void gs_data::nonlocal_info::nonlocal_vec( realType* u, uint n, int op, MPI_Comm comm )
{
 MPI_Status status;
 uint np = this->_np;
 MPI_Request* reqs = this->_reqs;
 uint* targ = this->_target;
 uint* nshared = this->_nshared;
 uint* sh_ind = this->_sh_ind;
 uint id;
 realType *buf = this->_buf, *start;
 uint size = n * sizeof( realType );
 unsigned int i;
 {
 MPI_Comm_rank( comm, (int*)&i );
 id = i;
 }
 for( i = 0; i < np; ++i )
 {
 uint ns = nshared[i], c = ns;
 start = buf;
 for( ; c; --c )
 {
 memcpy( buf, u + n * ( *sh_ind++ ), size );
 buf += n;
 }
 MPI_Isend( (void*)start, ns * size, MPI_UNSIGNED_CHAR, targ[i], id, comm, reqs++ );
 }
 start = buf;
 for( i = 0; i < np; ++i )
 {
 int nsn = n * nshared[i];
 MPI_Irecv( (void*)start, nsn * size, MPI_UNSIGNED_CHAR, targ[i], targ[i], comm, reqs++ );
 start += nsn;
 }
 for( reqs = this->_reqs, i = np * 2; i; --i )
 MPI_Wait( reqs++, &status );
 sh_ind = this->_sh_ind;
#define LOOP( OP ) \
 do \
 { \
 for( i = 0; i < np; ++i ) \
 { \
 uint c, j; \
 for( c = nshared[i]; c; --c ) \
 { \
 realType* uu = u + n * ( *sh_ind++ ); \
 for( j = n; j; --j ) \
 { \
 OP( *uu, *buf ); \
 ++uu; \
 ++buf; \
 } \
 } \
 } \
 } while( 0 )
 switch( op )
 {
 case GS_OP_ADD:
 LOOP( DO_ADD );
 break;
 case GS_OP_MUL:
 LOOP( DO_MUL );
 break;
 case GS_OP_MIN:
 LOOP( DO_MIN );
 break;
 case GS_OP_MAX:
 LOOP( DO_MAX );
 break;
 case GS_OP_BPR:
 LOOP( DO_BPR );
 break;
 }
#undef LOOP
}
 
void gs_data::nonlocal_info::nonlocal_many( realType** u, uint n, int op, MPI_Comm comm )
{
 MPI_Status status;
 uint np = this->_np;
 MPI_Request* reqs = this->_reqs;
 uint* targ = this->_target;
 uint* nshared = this->_nshared;
 uint* sh_ind = this->_sh_ind;
 uint id;
 realType *buf = this->_buf, *start;
 unsigned int i;
 {
 MPI_Comm_rank( comm, (int*)&i );
 id = i;
 }
 for( i = 0; i < np; ++i )
 {
 uint c, j, ns = nshared[i];
 start = buf;
 for( j = 0; j < n; ++j )
 {
 realType* uu = u[j];
 for( c = 0; c < ns; ++c )
 *buf++ = uu[sh_ind[c]];
 }
 sh_ind += ns;
 MPI_Isend( (void*)start, n * ns * sizeof( realType ), MPI_UNSIGNED_CHAR, targ[i], id, comm, reqs++ );
 }
 start = buf;
 for( i = 0; i < np; ++i )
 {
 int nsn = n * nshared[i];
 MPI_Irecv( (void*)start, nsn * sizeof( realType ), MPI_UNSIGNED_CHAR, targ[i], targ[i], comm, reqs++ );
 start += nsn;
 }
 for( reqs = this->_reqs, i = np * 2; i; --i )
 MPI_Wait( reqs++, &status );
 sh_ind = this->_sh_ind;
#define LOOP( OP ) \
 do \
 { \
 for( i = 0; i < np; ++i ) \
 { \
 uint c, j, ns = nshared[i]; \
 for( j = 0; j < n; ++j ) \
 { \
 realType* uu = u[j]; \
 for( c = 0; c < ns; ++c ) \
 { \
 OP( uu[sh_ind[c]], *buf ); \
 ++buf; \
 } \
 } \
 sh_ind += ns; \
 } \
 } while( 0 )
 switch( op )
 {
 case GS_OP_ADD:
 LOOP( DO_ADD );
 break;
 case GS_OP_MUL:
 LOOP( DO_MUL );
 break;
 case GS_OP_MIN:
 LOOP( DO_MIN );
 break;
 case GS_OP_MAX:
 LOOP( DO_MAX );
 break;
 case GS_OP_BPR:
 LOOP( DO_BPR );
 break;
 }
#undef LOOP
}
 
/*---------------------------------------------------------------------------
 MOAB Crystal Router
 ---------------------------------------------------------------------------*/
gs_data::crystal_data::crystal_data() {}
 
void gs_data::crystal_data::initialize( MPI_Comm comm )
{
 int num, id;
 buffers[0].buf.buffer_init( 1024 );
 buffers[1].buf.buffer_init( 1024 );
 buffers[2].buf.buffer_init( 1024 );
 all = &buffers[0];
 keep = &buffers[1];
 send = &buffers[2];
 memcpy( &( this->_comm ), &comm, sizeof( MPI_Comm ) );
 MPI_Comm_rank( comm, &id );
 this->_id = id;
 MPI_Comm_size( comm, &num );
 this->_num = num;
}
 
void gs_data::crystal_data::reset()
{
 buffers[0].buf.reset();
 buffers[1].buf.reset();
 buffers[2].buf.reset();
 keep = NULL;
 all = NULL;
 send = NULL;
}
 
// Partition data before sending
void gs_data::crystal_data::partition( uint cutoff, crystal_buf* lo, crystal_buf* hi )
{
 const uint* src = (uint*)all->buf.ptr;
 const uint* end = (uint*)src + all->n;
 uint *target, *lop, *hip;
 lo->n = hi->n = 0;
 lo->buf.buffer_reserve( all->n * sizeof( uint ) );
 hi->buf.buffer_reserve( all->n * sizeof( uint ) );
 lop = (uint*)lo->buf.ptr;
 hip = (uint*)hi->buf.ptr;
 while( src != end )
 {
 uint chunk_len = 3 + src[2];
 if( src[0] < cutoff )
 {
 target = lop;
 lo->n += chunk_len;
 lop += chunk_len;
 }
 else
 {
 target = hip;
 hi->n += chunk_len;
 hip += chunk_len;
 }
 memcpy( target, src, chunk_len * sizeof( uint ) );
 src += chunk_len;
 }
}
 
// Send message to target process
void gs_data::crystal_data::send_( uint target, int recvn )
{
 MPI_Request req[3] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL };
 MPI_Status status[3];
 uint count[2] = { 0, 0 }, sum, *recv[2];
 crystal_buf* t;
 int i;
 
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( &send->n, sizeof( uint ) );
 MPI_Isend( (void*)&send->n, sizeof( uint ), MPI_UNSIGNED_CHAR, target, _id, _comm, &req[0] );
 for( i = 0; i < recvn; ++i )
 MPI_Irecv( (void*)&count[i], sizeof( uint ), MPI_UNSIGNED_CHAR, target + i, target + i, _comm, &req[i + 1] );
 MPI_Waitall( recvn + 1, req, status );
 sum = keep->n;
 for( i = 0; i < recvn; ++i )
 sum += count[i];
 keep->buf.buffer_reserve( sum * sizeof( uint ) );
 recv[0] = (uint*)keep->buf.ptr;
 recv[0] += keep->n;
 recv[1] = recv[0] + count[0];
 keep->n = sum;
 
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( send->buf.ptr, send->n * sizeof( uint ) );
 MPI_Isend( (void*)send->buf.ptr, send->n * sizeof( uint ), MPI_UNSIGNED_CHAR, target, _id, _comm, &req[0] );
 if( recvn )
 {
 MPI_Irecv( (void*)recv[0], count[0] * sizeof( uint ), MPI_UNSIGNED_CHAR, target, target, _comm, &req[1] );
 if( recvn == 2 )
 MPI_Irecv( (void*)recv[1], count[1] * sizeof( uint ), MPI_UNSIGNED_CHAR, target + 1, target + 1, _comm,
 &req[2] );
 }
 MPI_Waitall( recvn + 1, req, status );
 
 t = all;
 all = keep;
 keep = t;
}
 
void gs_data::crystal_data::crystal_router()
{
 uint bl = 0, bh, n = _num, nl, target;
 int recvn;
 crystal_buf *lo, *hi;
 while( n > 1 )
 {
 nl = n / 2, bh = bl + nl;
 if( _id < bh )
 {
 target = _id + nl;
 recvn = ( n & 1 && _id == bh - 1 ) ? 2 : 1;
 lo = keep;
 hi = send;
 }
 else
 {
 target = _id - nl;
 recvn = ( target == bh ) ? ( --target, 0 ) : 1;
 hi = keep;
 lo = send;
 }
 partition( bh, lo, hi );
 send_( target, recvn );
 if( _id < bh )
 n = nl;
 else
 {
 n -= nl;
 bl = bh;
 }
 }
}
 
#define UINT_PER_X( X ) ( ( sizeof( X ) + sizeof( uint ) - 1 ) / sizeof( uint ) )
#define UINT_PER_REAL UINT_PER_X( realType )
#define UINT_PER_LONG UINT_PER_X( slong )
#define UINT_PER_ULONG UINT_PER_X( Ulong )
 
/*-------------------------------------------------------------------------
 Transfer
 -------------------------------------------------------------------------*/
ErrorCode gs_data::crystal_data::gs_transfer( int dynamic, moab::TupleList& tl, unsigned pf )
{
 
 unsigned mi, ml, mul, mr;
 tl.getTupleSize( mi, ml, mul, mr );
 
 const unsigned tsize = ( mi - 1 ) + ml * UINT_PER_LONG + mul * UINT_PER_ULONG + mr * UINT_PER_REAL;
 sint p, lp = -1;
 sint* ri;
 slong* rl;
 Ulong* rul;
 realType* rr;
 uint i, j, *buf, *len = 0, *buf_end;
 
 /* sort to group by target proc */
 if( pf >= mi ) return moab::MB_MEMORY_ALLOCATION_FAILED;
 // fail("pf expected to be in vi array (%d, %d)", pf, mi);
 tl.sort( pf, &all->buf );
 
 /* pack into buffer for crystal router */
 all->buf.buffer_reserve( ( tl.get_n() * ( 3 + tsize ) ) * sizeof( uint ) );
 all->n = 0;
 buf = (uint*)all->buf.ptr;
 
 bool canWrite = tl.get_writeEnabled();
 if( !canWrite ) tl.enableWriteAccess();
 
 ri = tl.vi_wr;
 rl = tl.vl_wr;
 rul = tl.vul_wr;
 rr = tl.vr_wr;
 
 for( i = tl.get_n(); i; --i )
 {
 p = ri[pf];
 if( p != lp )
 {
 lp = p;
 *buf++ = p; /* target */
 *buf++ = _id; /* source */
 len = buf++;
 *len = 0; /* length */
 all->n += 3;
 }
 for( j = 0; j < mi; ++j, ++ri )
 if( j != pf ) *buf++ = *ri;
 for( j = ml; j; --j, ++rl )
 {
 memcpy( buf, rl, sizeof( slong ) );
 buf += UINT_PER_LONG;
 }
 for( j = mul; j; --j, ++rul )
 {
 memcpy( buf, rul, sizeof( Ulong ) );
 buf += UINT_PER_ULONG;
 }
 for( j = mr; j; --j, ++rr )
 {
 memcpy( buf, rr, sizeof( realType ) );
 buf += UINT_PER_REAL;
 }
 *len += tsize, all->n += tsize;
 }
 
 crystal_router();
 
 /* unpack */
 buf = (uint*)all->buf.ptr;
 buf_end = buf + all->n;
 tl.set_n( 0 );
 ri = tl.vi_wr;
 rl = tl.vl_wr;
 rul = tl.vul_wr;
 rr = tl.vr_wr;
 
 while( buf != buf_end )
 {
 sint llen;
 buf++; /* target ( == this proc ) */
 p = *buf++; /* source */
 llen = *buf++; /* length */
 while( llen > 0 )
 {
 if( tl.get_n() == tl.get_max() )
 {
 if( !dynamic )
 {
 tl.set_n( tl.get_max() + 1 );
 if( !canWrite ) tl.disableWriteAccess();
 return moab::MB_SUCCESS;
 }
 ErrorCode rval = tl.resize( tl.get_max() + ( 1 + tl.get_max() ) / 2 + 1 );
 if( rval != moab::MB_SUCCESS )
 {
 if( !canWrite ) tl.disableWriteAccess();
 return rval;
 }
 ri = tl.vi_wr + mi * tl.get_n();
 rl = tl.vl_wr + ml * tl.get_n();
 rul = tl.vul_wr + mul * tl.get_n();
 rr = tl.vr_wr + mr * tl.get_n();
 }
 tl.inc_n();
 for( j = 0; j < mi; ++j )
 if( j != pf )
 *ri++ = *buf++;
 else
 *ri++ = p;
 for( j = ml; j; --j )
 {
 memcpy( rl++, buf, sizeof( slong ) );
 buf += UINT_PER_LONG;
 }
 for( j = mul; j; --j )
 {
 memcpy( rul++, buf, sizeof( Ulong ) );
 buf += UINT_PER_ULONG;
 }
 for( j = mr; j; --j )
 {
 memcpy( rr++, buf, sizeof( realType ) );
 buf += UINT_PER_REAL;
 }
 llen -= tsize;
 }
 }
 
 if( !canWrite ) tl.disableWriteAccess();
 return moab::MB_SUCCESS;
}
#endif
 
/*--------------------------------------------------------------------------
 Combined Execution
 --------------------------------------------------------------------------*/
 
void gs_data::gs_data_op( realType* u, int op )
{
 local_condense( u, op, this->local_cm );
#ifdef MOAB_HAVE_MPI
 this->nlinfo->nonlocal( u, op, _comm );
#endif
 local_uncondense( u, local_cm );
}
 
void gs_data::gs_data_op_vec( realType* u, uint n, int op )
{
#ifdef MOAB_HAVE_MPI
 if( n > nlinfo->_maxv )
 moab::fail( "%s: initialized with max vec size = %d,"
 " but called with vec size = %d\n",
 __FILE__, nlinfo->_maxv, n );
#endif
 local_condense_vec( u, n, op, local_cm );
#ifdef MOAB_HAVE_MPI
 this->nlinfo->nonlocal_vec( u, n, op, _comm );
#endif
 local_uncondense_vec( u, n, local_cm );
}
 
void gs_data::gs_data_op_many( realType** u, uint n, int op )
{
 uint i;
#ifdef MOAB_HAVE_MPI
 if( n > nlinfo->_maxv )
 moab::fail( "%s: initialized with max vec size = %d,"
 " but called with vec size = %d\n",
 __FILE__, nlinfo->_maxv, n );
#endif
 for( i = 0; i < n; ++i )
 local_condense( u[i], op, local_cm );
 
 moab::fail( "%s: initialized with max vec size = %d,"
 " but called with vec size = %d\n",
 __FILE__, 6, n );
 
#ifdef MOAB_HAVE_MPI
 this->nlinfo->nonlocal_many( u, n, op, _comm );
#endif
 for( i = 0; i < n; ++i )
 local_uncondense( u[i], local_cm );
}
 
/*--------------------------------------------------------------------------
 Setup
 --------------------------------------------------------------------------*/
 
ErrorCode gs_data::initialize( uint n,
 const long* label,
 const Ulong* ulabel,
 uint maxv,
 const unsigned int nlabels,
 const unsigned int nulabels,
 crystal_data* crystal )
{
 nlinfo = NULL;
 unsigned int j;
 TupleList nonzero, primary;
 ErrorCode rval;
#ifdef MOAB_HAVE_MPI
 TupleList shared;
#else
 moab::TupleList::buffer buf;
#endif
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( label, nlabels * sizeof( long ) );
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( ulabel, nlabels * sizeof( Ulong ) );
#ifdef MOAB_HAVE_MPI
 MPI_Comm_dup( crystal->_comm, &this->_comm );
#else
 buf.buffer_init( 1024 );
#endif
 
 /* construct list of nonzeros: (index ^, label) */
 nonzero.initialize( 1, nlabels, nulabels, 0, n );
 nonzero.enableWriteAccess();
 {
 uint i;
 sint* nzi;
 long* nzl;
 Ulong* nzul;
 nzi = nonzero.vi_wr;
 nzl = nonzero.vl_wr;
 nzul = nonzero.vul_wr;
 for( i = 0; i < n; ++i )
 if( label[nlabels * i] != 0 )
 {
 nzi[0] = i;
 for( j = 0; j < nlabels; j++ )
 nzl[j] = label[nlabels * i + j];
 for( j = 0; j < nulabels; j++ )
 nzul[j] = ulabel[nulabels * i + j];
 nzi++;
 nzl += nlabels;
 nzul += nulabels;
 nonzero.inc_n();
 }
 }
 
 /* sort nonzeros by label: (index ^2, label ^1) */
#ifndef MOAB_HAVE_MPI
 nonzero.sort( 1, &buf );
#else
 nonzero.sort( 1, &crystal->all->buf );
#endif
 
 /* build list of unique labels w/ lowest associated index:
 (index in nonzero ^, primary (lowest) index in label, count, label(s),
 ulabel(s)) */
 primary.initialize( 3, nlabels, nulabels, 0, nonzero.get_n() );
 primary.enableWriteAccess();
 {
 uint i;
 sint *nzi = nonzero.vi_wr, *pi = primary.vi_wr;
 slong *nzl = nonzero.vl_wr, *pl = primary.vl_wr;
 Ulong *nzul = nonzero.vul_wr, *pul = primary.vul_wr;
 slong last = -1;
 for( i = 0; i < nonzero.get_n(); ++i, nzi += 1, nzl += nlabels, nzul += nulabels )
 {
 if( nzl[0] == last )
 {
 ++pi[-1];
 continue;
 }
 last = nzl[0];
 pi[0] = i;
 pi[1] = nzi[0];
 for( j = 0; j < nlabels; j++ )
 pl[j] = nzl[j];
 for( j = 0; j < nulabels; j++ )
 pul[j] = nzul[j];
 pi[2] = 1;
 pi += 3, pl += nlabels;
 pul += nulabels;
 primary.inc_n();
 }
 }
 
 /* calculate size of local condense map */
 {
 uint i, count = 1;
 sint* pi = primary.vi_wr;
 for( i = primary.get_n(); i; --i, pi += 3 )
 if( pi[2] > 1 ) count += pi[2] + 1;
 this->local_cm = (sint*)malloc( count * sizeof( sint ) );
 }
 
 /* sort unique labels by primary index:
 (nonzero index ^2, primary index ^1, count, label ^2) */
#ifndef MOAB_HAVE_MPI
 primary.sort( 0, &buf );
 buf.reset();
 // buffer_free(&buf);
#else
 primary.sort( 0, &crystal->all->buf );
#endif
 
 /* construct local condense map */
 {
 uint i, ln;
 sint* pi = primary.vi_wr;
 sint* cm = this->local_cm;
 for( i = primary.get_n(); i > 0; --i, pi += 3 )
 if( ( ln = pi[2] ) > 1 )
 {
 sint* nzi = nonzero.vi_wr + 1 * pi[0];
 for( j = ln; j > 0; --j, nzi += 1 )
 *cm++ = nzi[0];
 *cm++ = -1;
 }
 *cm++ = -1;
 }
 nonzero.reset();
#ifndef MOAB_HAVE_MPI
 primary.reset();
#else
 /* assign work proc by label modulo np */
 {
 uint i;
 sint* pi = primary.vi_wr;
 slong* pl = primary.vl_wr;
 for( i = primary.get_n(); i; --i, pi += 3, pl += nlabels )
 pi[0] = pl[0] % crystal->_num;
 }
 rval = crystal->gs_transfer( 1, primary, 0 ); /* transfer to work procs */
 if( rval != MB_SUCCESS ) return rval;
 /* primary: (source proc, index on src, useless, label) */
 /* sort by label */
 primary.sort( 3, &crystal->all->buf );
 /* add sentinel to primary list */
 if( primary.get_n() == primary.get_max() )
 primary.resize( ( primary.get_max() ? primary.get_max() + ( primary.get_max() + 1 ) / 2 + 1 : 2 ) );
 primary.vl_wr[nlabels * primary.get_n()] = -1;
 /* construct shared list: (proc1, proc2, index1, label) */
#ifdef MOAB_HAVE_MPI
 shared.initialize( 3, nlabels, nulabels, 0, primary.get_n() );
 shared.enableWriteAccess();
#endif
 {
 sint *pi1 = primary.vi_wr, *si = shared.vi_wr;
 slong lbl, *pl1 = primary.vl_wr, *sl = shared.vl_wr;
 Ulong *pul1 = primary.vul_wr, *sul = shared.vul_wr;
 for( ; ( lbl = pl1[0] ) != -1; pi1 += 3, pl1 += nlabels, pul1 += nulabels )
 {
 sint* pi2 = pi1 + 3;
 slong* pl2 = pl1 + nlabels;
 Ulong* pul2 = pul1 + nulabels;
 for( ; pl2[0] == lbl; pi2 += 3, pl2 += nlabels, pul2 += nulabels )
 {
 if( shared.get_n() + 2 > shared.get_max() )
 shared.resize( ( shared.get_max() ? shared.get_max() + ( shared.get_max() + 1 ) / 2 + 1 : 2 ) ),
 si = shared.vi_wr + shared.get_n() * 3;
 sl = shared.vl_wr + shared.get_n() * nlabels;
 sul = shared.vul_wr + shared.get_n() * nulabels;
 si[0] = pi1[0];
 si[1] = pi2[0];
 si[2] = pi1[1];
 for( j = 0; j < nlabels; j++ )
 sl[j] = pl2[j];
 for( j = 0; j < nulabels; j++ )
 sul[j] = pul2[j];
 si += 3;
 sl += nlabels;
 sul += nulabels;
 shared.inc_n();
 si[0] = pi2[0];
 si[1] = pi1[0];
 si[2] = pi2[1];
 for( j = 0; j < nlabels; j++ )
 sl[j] = pl1[j];
 for( j = 0; j < nulabels; j++ )
 sul[j] = pul1[j];
 si += 3;
 sl += nlabels;
 sul += nulabels;
 shared.inc_n();
 }
 }
 }
 primary.reset();
 rval = crystal->gs_transfer( 1, shared, 0 ); /* segfaulting transfer to dest procs */
 if( rval != MB_SUCCESS ) return rval;
 /* shared list: (useless, proc2, index, label) */
 /* sort by label */
 shared.sort( 3, &crystal->all->buf );
 /* sort by partner proc */
 shared.sort( 1, &crystal->all->buf );
 /* count partner procs */
 {
 uint i, count = 0;
 sint proc = -1, *si = shared.vi_wr;
 for( i = shared.get_n(); i; --i, si += 3 )
 if( si[1] != proc )
 {
 ++count;
 proc = si[1];
 }
 // this->nlinfo = new nonlocal_info();
 // this->nlinfo->initialize(count,shared.get_n(),
 // nlabels, nulabels, maxv);
 this->nlinfo = new nonlocal_info( count, shared.get_n(), nlabels, nulabels, maxv );
 }
 /* construct non-local info */
 {
 uint i;
 sint proc = -1, *si = shared.vi_wr;
 slong* sl = shared.vl_wr;
 Ulong* ul = shared.vul_wr;
 uint* target = this->nlinfo->_target;
 uint* nshared = this->nlinfo->_nshared;
 uint* sh_ind = this->nlinfo->_sh_ind;
 slong* slabels = this->nlinfo->_slabels;
 Ulong* ulabels = this->nlinfo->_ulabels;
 for( i = shared.get_n(); i; --i, si += 3 )
 {
 if( si[1] != proc )
 {
 proc = si[1];
 *target++ = proc;
 *nshared++ = 0;
 }
 ++nshared[-1];
 *sh_ind++ = si[2];
 // don't store 1st slabel
 sl++;
 for( j = 0; j < nlabels - 1; j++ )
 slabels[j] = sl[j];
 for( j = 0; j < nulabels; j++ )
 ulabels[j] = ul[j];
 slabels += nlabels - 1;
 ulabels += nulabels;
 sl += nlabels - 1;
 ul += nulabels;
 }
 }
 shared.reset();
#endif
 return MB_SUCCESS;
}
 
void gs_data::reset()
{
 free( local_cm );
 local_cm = NULL;
#ifdef MOAB_HAVE_MPI
 if( nlinfo != NULL )
 {
 nlinfo->nlinfo_free();
 delete this->nlinfo;
 MPI_Comm_free( &_comm );
 nlinfo = NULL;
 }
#endif
}
 
} // namespace moab