From a3bc75e650a77e1aaa926e20a6ab6b1ca14ea57c Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Tue, 27 Nov 2018 11:20:50 -0500
Subject: [PATCH] PuReMD-old: remove non-source junk files. Clean-up to remove
debugging code for neutral territory. Conditionally toggle SpMV
implementation (Sparse_MatVec) debugging on if neutral territory is selected.
---
PuReMD/src/#basic_comm.c# | 322 ------------------------------------
PuReMD/src/allocate.c | 30 ++--
PuReMD/src/basic_comm.c | 35 ++--
PuReMD/src/comm_tools.c | 54 +++---
PuReMD/src/forces.c | 70 ++++----
PuReMD/src/grid.c | 7 +-
PuReMD/src/init_md.c | 4 +-
PuReMD/src/linear_solvers.c | 56 ++++++-
PuReMD/src/my* | 67 --------
PuReMD/src/neighbors.c | 6 +-
PuReMD/src/reax_types.h | 9 +-
11 files changed, 160 insertions(+), 500 deletions(-)
delete mode 100644 PuReMD/src/#basic_comm.c#
delete mode 100644 PuReMD/src/my*
diff --git a/PuReMD/src/#basic_comm.c# b/PuReMD/src/#basic_comm.c#
deleted file mode 100644
index efa01e7b..00000000
--- a/PuReMD/src/#basic_comm.c#
+++ /dev/null
@@ -1,322 +0,0 @@
-/*----------------------------------------------------------------------
- PuReMD - Purdue ReaxFF Molecular Dynamics Program
-
- Copyright (2010) Purdue University
- Hasan Metin Aktulga, haktulga@cs.purdue.edu
- Joseph Fogarty, jcfogart@mail.usf.edu
- Sagar Pandit, pandit@usf.edu
- Ananth Y Grama, ayg@cs.purdue.edu
-
- This program is free software; you can redistribute it and/or
- modify it under the terms of the GNU General Public License as
- published by the Free Software Foundation; either version 2 of
- the License, or (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
- See the GNU General Public License for more details:
- <http://www.gnu.org/licenses/>.
- ----------------------------------------------------------------------*/
-
-#include "reax_types.h"
-#if defined(PURE_REAX)
-#include "basic_comm.h"
-#include "vector.h"
-#elif defined(LAMMPS_REAX)
-#include "reax_basic_comm.h"
-#include "reax_vector.h"
-#endif
-
-#if defined(PURE_REAX)
-void real_packer( void *dummy, mpi_out_data *out_buf )
-{
- int i;
- real *buf = (real*) dummy;
- real *out = (real*) out_buf->out_atoms;
-
- for ( i = 0; i < out_buf->cnt; ++i )
- out[i] = buf[ out_buf->index[i] ];
-}
-
-
-void rvec_packer( void *dummy, mpi_out_data *out_buf )
-{
- int i;
- rvec *buf = (rvec*) dummy;
- rvec *out = (rvec*)out_buf->out_atoms;
-
- for ( i = 0; i < out_buf->cnt; ++i )
- memcpy( out[i], buf[ out_buf->index[i] ], sizeof(rvec) );
-}qsdfsdfsddsfsd
-
-
-VOIDuuuuuuuuu RVEC2_packer( void *dummy, mpi_out_data *out_buf )
-{
- int i;
- rvec2 *buf = (rvec2*) dummy;
- rvec2 *out = (rvec2*) out_buf->out_atoms;
-
- for ( i = 0; i < out_buf->cnt; ++i )
- memcpy( out[i], buf[ out_buf->index[i] ], sizeof(rvec2) );
-}
-
-
-void Dist( reax_system* system, mpi_datatypes *mpi_data,
- void *buf, MPI_Datatype type, int scale, dist_packer pack )
-{
- int d;
- mpi_out_data *out_bufs;
- MPI_Comm comm;
- MPI_Request req1, req2;
-q MPI_Status stat1, stat2;
- neighbor_proc *nbr1, *nbr2;
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d dist: entered\n", system->my_rank );
-#endif
- comm = mpi_data->comm_mesh3D;
- out_bufs = mpi_data->out_buffers;
-
- for ( d = 0; d < 3; ++d )
- {
- /* initiate recvs */
- nbr1 = &(system->my_nbrs[2 * d]);
- if ( nbr1->atoms_cnt )
- MPI_Irecv( buf + nbr1->atoms_str * scale, nbr1->atoms_cnt, type,
- nbr1->rank, 2 * d + 1, comm, &req1 );
-
- nbr2 = &(system->my_nbrs[2 * d + 1]);
- if ( nbr2->atoms_cnt )
- MPI_Irecv( buf + nbr2->atoms_str * scale, nbr2->atoms_cnt, type,
- nbr2->rank, 2 * d, comm, &req2 );
-
- /* send both messages in dimension d */
- if ( out_bufs[2 * d].cnt )
- {
- pack( buf, out_bufs + (2 * d) );
- MPI_Send( out_bufs[2 * d].out_atoms, out_bufs[2 * d].cnt, type,
- nbr1->rank, 2 * d, comm );
- }
-
- if ( out_bufs[2 * d + 1].cnt )
- {
- pack( buf, out_bufs + (2 * d + 1) );
- MPI_Send( out_bufs[2 * d + 1].out_atoms, out_bufs[2 * d + 1].cnt, type,
- nbr2->rank, 2 * d + 1, comm );
- }
-
- if ( nbr1->atoms_cnt ) MPI_Wait( &req1, &stat1 );
- if ( nbr2->atoms_cnt ) MPI_Wait( &req2, &stat2 );
- }
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d dist: done\n", system->my_rank );
-#endif
-}
-
-
-void real_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
-{
- int i;
- real *in = (real*) dummy_in;
- real *buf = (real*) dummy_buf;
-
- for ( i = 0; i < out_buf->cnt; ++i )
- buf[ out_buf->index[i] ] += in[i];
-}
-
-
-void rvec_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
-{
- int i;
- rvec *in = (rvec*) dummy_in;
- rvec *buf = (rvec*) dummy_buf;
-
- for ( i = 0; i < out_buf->cnt; ++i )
- {
- rvec_Add( buf[ out_buf->index[i] ], in[i] );
-#if defined(DEBUG)
- fprintf( stderr, "rvec_unpacker: cnt=%d i =%d index[i]=%d\n",
- out_buf->cnt, i, out_buf->index[i] );
-#endif
- }
-}
-
-
-void rvec2_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
-{
- int i;
- rvec2 *in = (rvec2*) dummy_in;
- rvec2 *buf = (rvec2*) dummy_buf;
-
- for ( i = 0; i < out_buf->cnt; ++i )
- {
- buf[ out_buf->index[i] ][0] += in[i][0];
- buf[ out_buf->index[i] ][1] += in[i][1];
- }
-}
-
-
-void Coll( reax_system* system, mpi_datatypes *mpi_data,
- void *buf, MPI_Datatype type, int scale, coll_unpacker unpack )
-{
- int d;
- void *in1, *in2;
- mpi_out_data *out_bufs;
- MPI_Comm comm;
- MPI_Request req1, req2;
- MPI_Status stat1, stat2;
- neighbor_proc *nbr1, *nbr2;
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d coll: entered\n", system->my_rank );
-#endif
- comm = mpi_data->comm_mesh3D;
- in1 = mpi_data->in1_buffer;
- in2 = mpi_data->in2_buffer;
- out_bufs = mpi_data->out_buffers;
-
- for ( d = 2; d >= 0; --d )
- {
- /* initiate recvs */
- nbr1 = &(system->my_nbrs[2 * d]);
- if ( out_bufs[2 * d].cnt )
- MPI_Irecv(in1, out_bufs[2 * d].cnt, type, nbr1->rank, 2 * d + 1, comm, &req1);
-
- nbr2 = &(system->my_nbrs[2 * d + 1]);
- if ( out_bufs[2 * d + 1].cnt )
- MPI_Irecv(in2, out_bufs[2 * d + 1].cnt, type, nbr2->rank, 2 * d, comm, &req2);
-
- /* send both messages in dimension d */
- if ( nbr1->atoms_cnt )
- MPI_Send( buf + nbr1->atoms_str * scale, nbr1->atoms_cnt, type,
- nbr1->rank, 2 * d, comm );
-
- if ( nbr2->atoms_cnt )
- MPI_Send( buf + nbr2->atoms_str * scale, nbr2->atoms_cnt, type,
- nbr2->rank, 2 * d + 1, comm );
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d coll[%d] nbr1: str=%d cnt=%d recv=%d\n",
- system->my_rank, d, nbr1->atoms_str, nbr1->atoms_cnt,
- out_bufs[2 * d].cnt );
- fprintf( stderr, "p%d coll[%d] nbr2: str=%d cnt=%d recv=%d\n",
- system->my_rank, d, nbr2->atoms_str, nbr2->atoms_cnt,
- out_bufs[2 * d + 1].cnt );
-#endif
-
- if ( out_bufs[2 * d].cnt )
- {
- MPI_Wait( &req1, &stat1 );
- unpack( in1, buf, out_bufs + (2 * d) );
- }
-
- if ( out_bufs[2 * d + 1].cnt )
- {
- MPI_Wait( &req2, &stat2 );
- unpack( in2, buf, out_bufs + (2 * d + 1) );
- }
- }
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d coll: done\n", system->my_rank );
-#endif
-}
-#endif /*PURE_REAX*/
-
-/*****************************************************************************/
-real Parallel_Norm( real *v, int n, MPI_Comm comm )
-{
- int i;
- real my_sum, norm_sqr;
-
- my_sum = 0;
- for ( i = 0; i < n; ++i )
- my_sum += SQR( v[i] );
-
- MPI_Allreduce( &my_sum, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, comm );
-
- return sqrt( norm_sqr );
-}
-
-
-
-real Parallel_Dot( real *v1, real *v2, int n, MPI_Comm comm )
-{
- int i;
- real my_dot, res;
-
- my_dot = 0;
- for ( i = 0; i < n; ++i )
- my_dot += v1[i] * v2[i];
-
- MPI_Allreduce( &my_dot, &res, 1, MPI_DOUBLE, MPI_SUM, comm );
-
- return res;
-}
-
-
-
-real Parallel_Vector_Acc( real *v, int n, MPI_Comm comm )
-{
- int i;
- real my_acc, res;
-
- my_acc = 0;
- for ( i = 0; i < n; ++i )
- my_acc += v[i];
-
- MPI_Allreduce( &my_acc, &res, 1, MPI_DOUBLE, MPI_SUM, comm );
-
- return res;
-}
-
-
-/*****************************************************************************/
-#if defined(TEST_FORCES)
-void Coll_ids_at_Master( reax_system *system, storage *workspace,
- mpi_datatypes *mpi_data )
-{
- int i;
- int *id_list;
-
- MPI_Gather( &system->n, 1, MPI_INT, workspace->rcounts, 1, MPI_INT,
- MASTER_NODE, mpi_data->world );
-
- if ( system->my_rank == MASTER_NODE )
- {
- workspace->displs[0] = 0;
- for ( i = 1; i < system->wsize; ++i )
- workspace->displs[i] = workspace->displs[i - 1] + workspace->rcounts[i - 1];
- }
-
- id_list = (int*) malloc( system->n * sizeof(int) );
- for ( i = 0; i < system->n; ++i )
- id_list[i] = system->my_atoms[i].orig_id;
-
- MPI_Gatherv( id_list, system->n, MPI_INT,
- workspace->id_all, workspace->rcounts, workspace->displs,
- MPI_INT, MASTER_NODE, mpi_data->world );
-
- sfree( id_list, "id_list" );
-
-#if defined(DEBUG)
- if ( system->my_rank == MASTER_NODE )
- {
- for ( i = 0 ; i < system->bigN; ++i )
- fprintf( stderr, "id_all[%d]: %d\n", i, workspace->id_all[i] );
- }
-#endif
-}
-
-
-void Coll_rvecs_at_Master( reax_system *system, storage *workspace,
- mpi_datatypes *mpi_data, rvec* v )
-{
- MPI_Gatherv( v, system->n, mpi_data->mpi_rvec,
- workspace->f_all, workspace->rcounts, workspace->displs,
- mpi_data->mpi_rvec, MASTER_NODE, mpi_data->world );
-}
-
-#endif
diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index d638b8da..9d634448 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -741,17 +741,16 @@ int Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
{
//TODO: add an exact number
/* in buffers */
- mpi_data->in_nt_buffer[i] = (void*)
- scalloc( my_nt_nbrs[i].est_recv, sizeof(real), "mpibuf:in_nt_buffer", comm );
+ mpi_data->in_nt_buffer[i] = scalloc( my_nt_nbrs[i].est_recv, sizeof(real),
+ "mpibuf:in_nt_buffer", comm );
//TODO
- mpi_buf = &( mpi_data->out_nt_buffers[i] );
+ mpi_buf = &mpi_data->out_nt_buffers[i];
/* allocate storage for the neighbor processor i */
- mpi_buf->index = (int*)
- scalloc( my_nt_nbrs[i].est_send, sizeof(int), "mpibuf:nt_index", comm );
- mpi_buf->out_atoms = (void*)
- scalloc( my_nt_nbrs[i].est_send, sizeof(real), "mpibuf:nt_out_atoms",
- comm );
+ mpi_buf->index = scalloc( my_nt_nbrs[i].est_send, sizeof(int),
+ "mpibuf:nt_index", comm );
+ mpi_buf->out_atoms = scalloc( my_nt_nbrs[i].est_send, sizeof(real),
+ "mpibuf:nt_out_atoms", comm );
}
#endif
@@ -769,7 +768,7 @@ void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data )
for ( i = 0; i < MAX_NBRS; ++i )
{
- mpi_buf = &( mpi_data->out_buffers[i] );
+ mpi_buf = &mpi_data->out_buffers[i];
sfree( mpi_buf->index, "mpibuf:index" );
sfree( mpi_buf->out_atoms, "mpibuf:out_atoms" );
}
@@ -779,7 +778,7 @@ void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data )
{
sfree( mpi_data->in_nt_buffer[i], "in_nt_buffer" );
- mpi_buf = &( mpi_data->out_nt_buffers[i] );
+ mpi_buf = &mpi_data->out_nt_buffers[i];
sfree( mpi_buf->index, "mpibuf:nt_index" );
sfree( mpi_buf->out_atoms, "mpibuf:nt_out_atoms" );
}
@@ -1037,12 +1036,13 @@ void ReAllocate( reax_system *system, control_params *control,
}
#if defined(NEUTRAL_TERRITORY)
- // also check individual outgoing Neutral Territory buffers
+ /* also check individual outgoing Neutral Territory buffers */
mpi_flag = 0;
for ( p = 0; p < REAX_MAX_NT_NBRS; ++p )
{
- nbr_pr = &( system->my_nt_nbrs[p] );
- nbr_data = &( mpi_data->out_nt_buffers[p] );
+ nbr_pr = &system->my_nt_nbrs[p];
+ nbr_data = &mpi_data->out_nt_buffers[p];
+
if ( nbr_data->cnt >= nbr_pr->est_send * 0.90 )
{
mpi_flag = 1;
@@ -1078,8 +1078,8 @@ void ReAllocate( reax_system *system, control_params *control,
#if defined(NEUTRAL_TERRITORY)
for ( p = 0; p < REAX_MAX_NT_NBRS; ++p )
{
- nbr_pr = &( system->my_nt_nbrs[p] );
- nbr_data = &( mpi_data->out_nt_buffers[p] );
+ nbr_pr = &system->my_nt_nbrs[p];
+ nbr_data = &mpi_data->out_nt_buffers[p];
nbr_pr->est_send = MAX( nbr_data->cnt * SAFER_ZONE, MIN_SEND );
}
#endif
diff --git a/PuReMD/src/basic_comm.c b/PuReMD/src/basic_comm.c
index a6606458..18d27fa1 100644
--- a/PuReMD/src/basic_comm.c
+++ b/PuReMD/src/basic_comm.c
@@ -170,16 +170,7 @@ void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
}
}
- /*for( d = 0; d < 6; ++d )
- {
- nbr = &(system->my_nbrs[d]);
- if ( nbr->atoms_cnt )
- {
- MPI_Wait( &(req[d]), &(stat[d]) );
- }
- }*/
-
- for( d = 0; d < count; ++d )
+ for ( d = 0; d < count; ++d )
{
MPI_Waitany( REAX_MAX_NT_NBRS, req, &index, stat);
}
@@ -198,7 +189,9 @@ void real_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
real *buf = (real*) dummy_buf;
for ( i = 0; i < out_buf->cnt; ++i )
+ {
buf[ out_buf->index[i] ] += in[i];
+ }
}
@@ -211,6 +204,7 @@ void rvec_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
for ( i = 0; i < out_buf->cnt; ++i )
{
rvec_Add( buf[ out_buf->index[i] ], in[i] );
+
#if defined(DEBUG)
fprintf( stderr, "rvec_unpacker: cnt=%d i =%d index[i]=%d\n",
out_buf->cnt, i, out_buf->index[i] );
@@ -247,6 +241,7 @@ void Coll( reax_system* system, mpi_datatypes *mpi_data,
#if defined(DEBUG)
fprintf( stderr, "p%d coll: entered\n", system->my_rank );
#endif
+
comm = mpi_data->comm_mesh3D;
in1 = mpi_data->in1_buffer;
in2 = mpi_data->in2_buffer;
@@ -299,6 +294,7 @@ void Coll( reax_system* system, mpi_datatypes *mpi_data,
#endif
}
+
#if defined(NEUTRAL_TERRITORY)
void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
void *buf, MPI_Datatype type, int scale, coll_unpacker unpack )
@@ -314,14 +310,14 @@ void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
#if defined(DEBUG)
fprintf( stderr, "p%d coll: entered\n", system->my_rank );
#endif
+
comm = mpi_data->comm_mesh3D;
out_bufs = mpi_data->out_nt_buffers;
count = 0;
for ( d = 0; d < 6; ++d )
{
- /* initiate recvs */
- nbr = &(system->my_nt_nbrs[d]);
+ nbr = &system->my_nt_nbrs[d];
in[d] = mpi_data->in_nt_buffer[d];
if ( out_bufs[d].cnt )
{
@@ -341,16 +337,7 @@ void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
}
}
- /*for( d = 0; d < 6; ++d )
- {
- if ( out_bufs[d].cnt )
- {
- MPI_Wait( &(req[d]), &(stat[d]));
- unpack( in[d], buf, out_bufs + d );
- }
- }*/
-
- for( d = 0; d < count; ++d )
+ for ( d = 0; d < count; ++d )
{
MPI_Waitany( REAX_MAX_NT_NBRS, req, &index, stat);
unpack( in[index], buf, out_bufs + index );
@@ -361,10 +348,9 @@ void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
#endif
}
#endif
-
-
#endif /*PURE_REAX*/
+
/*****************************************************************************/
real Parallel_Norm( real *v, int n, MPI_Comm comm )
{
@@ -458,5 +444,4 @@ void Coll_rvecs_at_Master( reax_system *system, storage *workspace,
workspace->f_all, workspace->rcounts, workspace->displs,
mpi_data->mpi_rvec, MASTER_NODE, mpi_data->world );
}
-
#endif
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 38315557..4a5fb513 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -49,19 +49,20 @@ void Setup_NT_Comm( reax_system* system, control_params* control,
{+1, 0, 0}, // +x
{+1, -1, 0} // +x-y
};
- my_box = &(system->my_box);
+ my_box = &system->my_box;
bndry_cut = system->bndry_cuts.ghost_cutoff;
- /* identify my neighbors */
system->num_nt_nbrs = REAX_MAX_NT_NBRS;
+
+ /* identify my neighbors */
for ( i = 0; i < system->num_nt_nbrs; ++i )
{
- nbr_pr = &(system->my_nt_nbrs[i]);
+ nbr_pr = &system->my_nt_nbrs[i];
ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
- MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &(nbr_pr->rank) );
+ MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &nbr_pr->rank );
/* set the rank of the neighbor processor in the receiving direction */
ivec_Sum( nbr_recv_coords, system->my_coords, r[i + 6] ); /* actual nbr coords */
- MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_recv_coords, &(nbr_pr->receive_rank) );
+ MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_recv_coords, &nbr_pr->receive_rank );
for ( d = 0; d < 3; ++d )
{
@@ -84,11 +85,17 @@ void Setup_NT_Comm( reax_system* system, control_params* control,
/* determine if it is a periodic neighbor */
if ( nbr_coords[d] < 0 )
+ {
nbr_pr->prdc[d] = -1;
+ }
else if ( nbr_coords[d] >= control->procs_by_dim[d] )
- nbr_pr->prdc[d] = +1;
+ {
+ nbr_pr->prdc[d] = 1;
+ }
else
+ {
nbr_pr->prdc[d] = 0;
+ }
}
}
@@ -110,14 +117,14 @@ int Sort_Neutral_Territory( reax_system *system, int dir, mpi_out_data *out_bufs
for ( i = 0; i < system->n; ++i )
{
- if ( nbr_pr->bndry_min[0] <= atoms[i].x[0] &&
- atoms[i].x[0] < nbr_pr->bndry_max[0] &&
- nbr_pr->bndry_min[1] <= atoms[i].x[1] &&
- atoms[i].x[1] < nbr_pr->bndry_max[1] &&
- nbr_pr->bndry_min[2] <= atoms[i].x[2] &&
- atoms[i].x[2] < nbr_pr->bndry_max[2] )
+ if ( nbr_pr->bndry_min[0] <= atoms[i].x[0]
+ && atoms[i].x[0] < nbr_pr->bndry_max[0]
+ && nbr_pr->bndry_min[1] <= atoms[i].x[1]
+ && atoms[i].x[1] < nbr_pr->bndry_max[1]
+ && nbr_pr->bndry_min[2] <= atoms[i].x[2]
+ && atoms[i].x[2] < nbr_pr->bndry_max[2] )
{
- if( write )
+ if ( write )
{
out_bufs[dir].index[out_bufs[dir].cnt] = i;
out_bufs[dir].cnt++;
@@ -152,21 +159,22 @@ void Init_Neutral_Territory( reax_system* system, mpi_datatypes *mpi_data )
for ( d = 0; d < 6; ++d )
{
- nbr = &( system->my_nt_nbrs[d] );
+ nbr = &system->my_nt_nbrs[d];
Sort_Neutral_Territory( system, d, out_bufs, 1 );
MPI_Irecv( &cnt, 1, MPI_INT, nbr->receive_rank, d, comm, &req );
- MPI_Send( &(out_bufs[d].cnt), 1, MPI_INT, nbr->rank, d, comm );
+ MPI_Send( &out_bufs[d].cnt, 1, MPI_INT, nbr->rank, d, comm );
MPI_Wait( &req, &stat );
- if( mpi_data->in_nt_buffer[d] == NULL )
+ if ( mpi_data->in_nt_buffer[d] == NULL )
{
//TODO
- mpi_data->in_nt_buffer[d] = (void *) smalloc( SAFER_ZONE * cnt * sizeof(real), "in", comm );
+ mpi_data->in_nt_buffer[d] = smalloc( SAFER_ZONE * cnt * sizeof(real),
+ "Init_Neural_Territory::mpi_data->in_nt_buffer[d]", comm );
}
- nbr = &(system->my_nt_nbrs[d]);
+ nbr = &system->my_nt_nbrs[d];
nbr->atoms_str = end;
nbr->atoms_cnt = cnt;
nbr->est_recv = MAX( SAFER_ZONE * cnt, MIN_SEND );
@@ -185,18 +193,20 @@ void Estimate_NT_Atoms( reax_system *system, mpi_datatypes *mpi_data )
out_bufs = mpi_data->out_nt_buffers;
- for( d = 0; d < 6; ++d )
+ for ( d = 0; d < 6; ++d )
{
/* count the number of atoms in each processor's outgoing list */
- nbr = &( system->my_nt_nbrs[d] );
+ nbr = &system->my_nt_nbrs[d];
nbr->est_send = Sort_Neutral_Territory( system, d, out_bufs, 0 );
/* estimate the space based on the count above */
nbr->est_send = MAX( MIN_SEND, nbr->est_send * SAFER_ZONE );
/* allocate the estimated space */
- out_bufs[d].index = (int*) calloc( nbr->est_send, sizeof(int) );
- out_bufs[d].out_atoms = (void*) calloc( nbr->est_send, sizeof(real) );
+ out_bufs[d].index = scalloc( nbr->est_send, sizeof(int),
+ "Estimate_NT_Atoms::out_bufs[d].index", MPI_COMM_WORLD );
+ out_bufs[d].out_atoms = scalloc( nbr->est_send, sizeof(real),
+ "Estimate_NT_Atoms::out_bufs[d].out_atoms", MPI_COMM_WORLD );
/* sort the atoms to their outgoing buffers */
// TODO: to call or not to call?
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index d5132fc0..eabda507 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -385,22 +385,27 @@ void Init_Forces( reax_system *system, control_params *control,
bin[i] = 0;
total_sum[i] = 0;
}
+
for ( i = system->n; i < system->N; ++i )
{
- atom_i = &(system->my_atoms[i]);
+ atom_i = &system->my_atoms[i];
+
if( atom_i->nt_dir != -1 )
{
total_cnt[ atom_i->nt_dir ]++;
}
}
+
total_sum[0] = system->n;
for ( i = 1; i < 6; ++i )
{
total_sum[i] = total_sum[i-1] + total_cnt[i-1];
}
+
for ( i = system->n; i < system->N; ++i )
{
- atom_i = &(system->my_atoms[i]);
+ atom_i = &system->my_atoms[i];
+
if( atom_i->nt_dir != -1 )
{
atom_i->pos = total_sum[ atom_i->nt_dir ] + bin[ atom_i->nt_dir ];
@@ -416,10 +421,11 @@ void Init_Forces( reax_system *system, control_params *control,
for ( i = 0; i < system->N; ++i )
{
- atom_i = &(system->my_atoms[i]);
+ atom_i = &system->my_atoms[i];
type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
- end_i = End_Index(i, far_nbrs);
+ end_i = End_Index(i, far_nbrs);
+
if ( far_nbrs->format == HALF_LIST )
{
/* start at end because other atoms
@@ -430,7 +436,7 @@ void Init_Forces( reax_system *system, control_params *control,
{
btop_i = Start_Index( i, bonds );
}
- sbp_i = &(system->reax_param.sbp[type_i]);
+ sbp_i = &system->reax_param.sbp[type_i];
if ( i < system->n )
{
@@ -438,7 +444,7 @@ void Init_Forces( reax_system *system, control_params *control,
cutoff = control->nonb_cut;
}
#if defined(NEUTRAL_TERRITORY)
- else if( atom_i->nt_dir != -1 )
+ else if ( atom_i->nt_dir != -1 )
{
local = 2;
cutoff = control->nonb_cut;
@@ -523,7 +529,8 @@ void Init_Forces( reax_system *system, control_params *control,
{
/* H matrix entry */
#if defined(NEUTRAL_TERRITORY)
- if( atom_j->nt_dir > 0 || (j < system->n && (H->format == SYM_FULL_MATRIX
+ if ( atom_j->nt_dir > 0 || (j < system->n
+ && (H->format == SYM_FULL_MATRIX
|| (H->format == SYM_HALF_MATRIX && i < j))) )
{
if( j < system->n )
@@ -543,6 +550,7 @@ void Init_Forces( reax_system *system, control_params *control,
{
H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
}
+
++Htop;
}
#else
@@ -551,7 +559,8 @@ void Init_Forces( reax_system *system, control_params *control,
|| far_nbrs->format == FULL_LIST )
{
if( j < system->n )
- H->entries[Htop].j = j;
+ H->entries[Htop].j = j;
+
if ( control->tabulate == 0 )
{
H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
@@ -560,6 +569,7 @@ void Init_Forces( reax_system *system, control_params *control,
{
H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
}
+
++Htop;
}
#endif
@@ -596,7 +606,8 @@ void Init_Forces( reax_system *system, control_params *control,
{
/* H matrix entry */
if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir]
- && ( H->format == SYM_FULL_MATRIX || (H->format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
+ && ( H->format == SYM_FULL_MATRIX
+ || (H->format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
|| ( j < system->n && atom_i->nt_dir != 0 && H->format == SYM_FULL_MATRIX ))
{
if( !nt_flag )
@@ -604,6 +615,7 @@ void Init_Forces( reax_system *system, control_params *control,
nt_flag = 1;
H->start[atom_i->pos] = Htop;
}
+
if( j < system->n )
{
H->entries[Htop].j = j;
@@ -621,6 +633,7 @@ void Init_Forces( reax_system *system, control_params *control,
{
H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
}
+
++Htop;
}
}
@@ -709,6 +722,7 @@ void Init_Forces( reax_system *system, control_params *control,
system->my_rank, data->step, Htop, num_bonds, num_hbonds );
MPI_Barrier( comm );
#endif
+
#if defined( DEBUG )
Print_Bonds( system, bonds, "debugbonds.out" );
Print_Bond_List2( system, bonds, "pbonds.out" );
@@ -989,18 +1003,16 @@ void Estimate_Storages( reax_system *system, control_params *control,
*num_3body = 0;
#if defined(NEUTRAL_TERRITORY)
- int mark[6];
- mark[0] = mark[1] = 1;
- mark[2] = mark[3] = mark[4] = mark[5] = 2;
+ int mark[6] = {1, 1, 2, 2, 2, 2};
#endif
for ( i = 0; i < system->N; ++i )
{
- atom_i = &(system->my_atoms[i]);
+ atom_i = &system->my_atoms[i];
type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
- end_i = End_Index(i, far_nbrs);
- sbp_i = &(system->reax_param.sbp[type_i]);
+ end_i = End_Index(i, far_nbrs);
+ sbp_i = &system->reax_param.sbp[type_i];
if ( i < system->n )
{
@@ -1030,9 +1042,11 @@ void Estimate_Storages( reax_system *system, control_params *control,
{
nbr_pj = &( far_nbrs->far_nbr_list[pj] );
j = nbr_pj->nbr;
- //if ( far_nbrs->format == HALF_LIST )
+#if !defined(NEUTRAL_TERRITORY)
+ if ( far_nbrs->format == HALF_LIST )
+#endif
{
- atom_j = &(system->my_atoms[j]);
+ atom_j = &system->my_atoms[j];
}
if (nbr_pj->d <= cutoff)
@@ -1079,9 +1093,6 @@ void Estimate_Storages( reax_system *system, control_params *control,
#if defined(NEUTRAL_TERRITORY)
else if ( local == 2 )
{
- /*if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir]
- && ( cm_format == SYM_FULL_MATRIX || (cm_format == SYM_HALF_MATRIX && atom_i->pos < atom_j->pos)))
- || ( j < system->n && atom_i->nt_dir != 0 && cm_format == SYM_FULL_MATRIX ))*/
if( ( atom_j->nt_dir != -1 && mark[atom_i->nt_dir] != mark[atom_j->nt_dir] )
|| ( j < system->n && atom_i->nt_dir != 0 ))
{
@@ -1133,18 +1144,21 @@ void Estimate_Storages( reax_system *system, control_params *control,
}
#if defined(NEUTRAL_TERRITORY)
-
- /* Since we don't know the NT atoms' position yet, Htop cannot be calculated accurately
- Therefore, we assume it is full and divide 2 if necessary */
- if( cm_format == SYM_HALF_MATRIX )
+ /* Since we don't know the NT atoms' position yet, Htop cannot be calculated accurately.
+ * Therefore, we assume it is full and divide 2 if necessary. */
+ if ( cm_format == SYM_HALF_MATRIX )
{
- *Htop = (*Htop + system->n)/2;
+ *Htop = (*Htop + system->n) / 2;
}
#endif
- *Htop = (int)(MAX( *Htop * SAFE_ZONE, MIN_CAP * MIN_HENTRIES ));
- *matrix_dim = (int)(MAX( *matrix_dim * SAFE_ZONE, MIN_CAP ));
+
+ *Htop = (int) MAX( *Htop * SAFE_ZONE, MIN_CAP * MIN_HENTRIES );
+ *matrix_dim = (int) MAX( *matrix_dim * SAFE_ZONE, MIN_CAP );
+
for ( i = 0; i < system->n; ++i )
- hb_top[i] = (int)(MAX( hb_top[i] * SAFER_ZONE, MIN_HBONDS ));
+ {
+ hb_top[i] = (int) MAX( hb_top[i] * SAFER_ZONE, MIN_HBONDS );
+ }
for ( i = 0; i < system->N; ++i )
{
diff --git a/PuReMD/src/grid.c b/PuReMD/src/grid.c
index ce251430..804999b4 100644
--- a/PuReMD/src/grid.c
+++ b/PuReMD/src/grid.c
@@ -70,12 +70,13 @@ void Mark_GCells( reax_system* system, grid *g, ivec procs, MPI_Comm comm )
{
if ( dir[i][d] > 0 )
{
- nt_str[d] = MIN(g->native_end[d], g->ncells[d]);
- nt_end[d] = MIN(g->native_end[d] + g->vlist_span[d], g->ncells[d]);
+ nt_str[d] = MIN( g->native_end[d], g->ncells[d] );
+ nt_end[d] = MIN( g->native_end[d] + g->vlist_span[d],
+ g->ncells[d] );
}
else if ( dir[i][d] < 0 )
{
- nt_str[d] = MAX(0, g->native_str[d] - g->vlist_span[d] );
+ nt_str[d] = MAX( 0, g->native_str[d] - g->vlist_span[d] );
nt_end[d] = g->native_str[d];
}
else
diff --git a/PuReMD/src/init_md.c b/PuReMD/src/init_md.c
index b8bf0a72..2df6f5e9 100644
--- a/PuReMD/src/init_md.c
+++ b/PuReMD/src/init_md.c
@@ -650,9 +650,9 @@ int Init_Lists( reax_system *system, control_params *control,
bond_top, &num_3body, comm, &matrix_dim, cm_format );
#if defined(NEUTRAL_TERRITORY)
- Allocate_Matrix( &(workspace->H), matrix_dim, Htop, cm_format, comm );
+ Allocate_Matrix( &workspace->H, matrix_dim, Htop, cm_format, comm );
#else
- Allocate_Matrix( &(workspace->H), system->local_cap, Htop, cm_format, comm );
+ Allocate_Matrix( &workspace->H, system->local_cap, Htop, cm_format, comm );
#endif
workspace->L = NULL;
workspace->U = NULL;
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index 6338c1af..68bca15d 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -1137,7 +1137,7 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
{
- int i, j, k, si, dim;
+ int i, j, k, si, num_rows;
real val;
for ( i = 0; i < N; ++i )
@@ -1146,13 +1146,11 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
}
#if defined(NEUTRAL_TERRITORY)
- dim = A->NT;
-#else
- dim = A->n;
-#endif
+ num_rows = A->NT;
+
if ( A->format == SYM_HALF_MATRIX )
{
- for ( i = 0; i < dim; ++i )
+ for ( i = 0; i < num_rows; ++i )
{
si = A->start[i];
@@ -1162,6 +1160,11 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
b[i] += A->entries[si].val * x[i];
k = si + 1;
}
+ /* zeros on the diagonal for i >= A->n,
+ * so skip the diagonal multplication step as zeros
+ * are not stored (idea: keep the NNZ's the same
+ * for full shell and neutral territory half-stored
+ * charge matrices to make debugging easier) */
else
{
k = si;
@@ -1173,14 +1176,13 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
val = A->entries[k].val;
b[i] += val * x[j];
- //if( j < A->n ) // comment out for tryQEq
b[j] += val * x[i];
}
}
}
else if ( A->format == SYM_FULL_MATRIX || A->format == FULL_MATRIX )
{
- for ( i = 0; i < dim; ++i )
+ for ( i = 0; i < num_rows; ++i )
{
si = A->start[i];
@@ -1193,6 +1195,44 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
}
}
}
+#else
+ num_rows = A->n;
+
+ if ( A->format == SYM_HALF_MATRIX )
+ {
+ for ( i = 0; i < num_rows; ++i )
+ {
+ si = A->start[i];
+
+ /* diagonal only contributes once */
+ b[i] += A->entries[si].val * x[i];
+
+ for ( k = si + 1; k < A->end[i]; ++k )
+ {
+ j = A->entries[k].j;
+ val = A->entries[k].val;
+
+ b[i] += val * x[j];
+ b[j] += val * x[i];
+ }
+ }
+ }
+ else if ( A->format == SYM_FULL_MATRIX || A->format == FULL_MATRIX )
+ {
+ for ( i = 0; i < num_rows; ++i )
+ {
+ si = A->start[i];
+
+ for ( k = si; k < A->end[i]; ++k )
+ {
+ j = A->entries[k].j;
+ val = A->entries[k].val;
+
+ b[i] += val * x[j];
+ }
+ }
+ }
+#endif
}
diff --git a/PuReMD/src/my* b/PuReMD/src/my*
deleted file mode 100644
index 6ad2fd10..00000000
--- a/PuReMD/src/my*
+++ /dev/null
@@ -1,67 +0,0 @@
-./integrate.h:1:/*----------------------------------------------------------------------
-./reset_tools.h:1:/*----------------------------------------------------------------------
-./allocate.c:1:/*----------------------------------------------------------------------
-./lookup.h:1:/*----------------------------------------------------------------------
-./restart.c:1:/*----------------------------------------------------------------------
-./box.c:1:/*----------------------------------------------------------------------
-./torsion_angles.h:1:/*----------------------------------------------------------------------
-./bond_orders.h:1:/*----------------------------------------------------------------------
-./reax_types.h:1:/*----------------------------------------------------------------------
-./grid.c:1:/*----------------------------------------------------------------------
-./io_tools.h:1:/*----------------------------------------------------------------------
-./nonbonded.h:1:/*----------------------------------------------------------------------
-./reax_defs.h:1:/*----------------------------------------------------------------------
-./basic_comm.h:1:/*----------------------------------------------------------------------
-./multi_body.h:1:/*----------------------------------------------------------------------
-./linear_solvers.h:1:/*----------------------------------------------------------------------
-./neighbors.c:1:/*----------------------------------------------------------------------
-./valence_angles.h:1:/*----------------------------------------------------------------------
-./tool_box.c:1:/*----------------------------------------------------------------------
-./analyze.h:1:/*----------------------------------------------------------------------
-./qEq.c:1:/*----------------------------------------------------------------------
-./control.c:1:/*----------------------------------------------------------------------
-./vector.c:1:/*----------------------------------------------------------------------
-./bonds.h:1:/*----------------------------------------------------------------------
-./forces.c:1:/*----------------------------------------------------------------------
-./system_props.c:1:/*----------------------------------------------------------------------
-./traj.h:1:/*----------------------------------------------------------------------
-./comm_tools.h:1:/*----------------------------------------------------------------------
-./ffield.h:1:/*----------------------------------------------------------------------
-./geo_tools.c:1:/*----------------------------------------------------------------------
-./hydrogen_bonds.h:1:/*----------------------------------------------------------------------
-./init_md.h:1:/*----------------------------------------------------------------------
-./list.h:1:/*----------------------------------------------------------------------
-./random.c:1:/*----------------------------------------------------------------------
-./linear_solvers.c:1:/*----------------------------------------------------------------------
-./neighbors.h:1:/*----------------------------------------------------------------------
-./valence_angles.c:1:/*----------------------------------------------------------------------
-./tool_box.h:1:/*----------------------------------------------------------------------
-./analyze.c:1:/*----------------------------------------------------------------------
-./nonbonded.c:1:/*----------------------------------------------------------------------
-./basic_comm.c:1:/*----------------------------------------------------------------------
-./multi_body.c:1:/*----------------------------------------------------------------------
-./restart.h:1:/*----------------------------------------------------------------------
-./box.h:1:/*----------------------------------------------------------------------
-./parallelreax.c:1:/*----------------------------------------------------------------------
-./torsion_angles.c:1:/*----------------------------------------------------------------------
-./bond_orders.c:1:/*----------------------------------------------------------------------
-./grid.h:1:/*----------------------------------------------------------------------
-./io_tools.c:1:/*----------------------------------------------------------------------
-./integrate.c:1:/*----------------------------------------------------------------------
-./reset_tools.c:1:/*----------------------------------------------------------------------
-./allocate.h:1:/*----------------------------------------------------------------------
-./lookup.c:1:/*----------------------------------------------------------------------
-./init_md.c:1:/*----------------------------------------------------------------------
-./list.c:1:/*----------------------------------------------------------------------
-./random.h:1:/*----------------------------------------------------------------------
-./geo_tools.h:1:/*----------------------------------------------------------------------
-./hydrogen_bonds.c:1:/*----------------------------------------------------------------------
-./forces.h:1:/*----------------------------------------------------------------------
-./system_props.h:1:/*----------------------------------------------------------------------
-./traj.c:1:/*----------------------------------------------------------------------
-./comm_tools.c:1:/*----------------------------------------------------------------------
-./ffield.c:1:/*----------------------------------------------------------------------
-./qEq.h:1:/*----------------------------------------------------------------------
-./control.h:1:/*----------------------------------------------------------------------
-./vector.h:1:/*----------------------------------------------------------------------
-./bonds.c:1:/*----------------------------------------------------------------------
diff --git a/PuReMD/src/neighbors.c b/PuReMD/src/neighbors.c
index 69f8c653..0bd33ea8 100644
--- a/PuReMD/src/neighbors.c
+++ b/PuReMD/src/neighbors.c
@@ -101,9 +101,9 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
//fprintf( stderr, "gridcell %d %d %d\n", i, j, k );
/* pick up an atom from the current cell */
- for (l = gci->str; l < gci->end; ++l )
+ for ( l = gci->str; l < gci->end; ++l )
{
- atom1 = &(system->my_atoms[l]);
+ atom1 = &system->my_atoms[l];
#if defined(NEUTRAL_TERRITORY)
if( gci->type >= NT_NBRS && gci->type < NT_NBRS + 6 )
{
@@ -214,7 +214,7 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists,
/* pick up an atom from the current cell */
for ( l = gci->str; l < gci->end; ++l )
{
- atom1 = &(system->my_atoms[l]);
+ atom1 = &system->my_atoms[l];
#if defined(NEUTRAL_TERRITORY)
if( gci->type >= NT_NBRS && gci->type < NT_NBRS + 6 )
{
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index e91117f5..7126a1e4 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -40,7 +40,7 @@
/************* SOME DEFS - crucial for reax_types.h *********/
#define PURE_REAX
-#define NEUTRAL_TERRITORY
+//#define NEUTRAL_TERRITORY
//#define LAMMPS_REAX
//#define DEBUG
//#define DEBUG_FOCUS
@@ -561,7 +561,7 @@ typedef struct
reax_atom *my_atoms;
#if defined(NEUTRAL_TERRITORY)
- int num_nt_nbrs;
+ int num_nt_nbrs;
#endif
} reax_system;
@@ -978,10 +978,9 @@ typedef struct
/* matrix storage format */
int format;
int cap, n, m;
- //TODO: uncomment
-//#if defined(NEUTRAL_TERRITORY)
+#if defined(NEUTRAL_TERRITORY)
int NT;
-//#endif
+#endif
int *start, *end;
sparse_matrix_entry *entries;
} sparse_matrix;
--
GitLab