#include <gtest/gtest.h>
#include "gromacs/gmxlib/network.h"
+#include "gromacs/math/paddedvector.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/forcerec.h"
#include "gromacs/mdtypes/enerdata.h"
ForceProviders forceProviders;
module->initForceProviders(&forceProviders);
- t_mdatoms md;
- PaddedRVecVector f = { { 0, 0, 0 } };
- gmx::ForceWithVirial forceWithVirial(f, true);
+ t_mdatoms md;
+ PaddedVector<gmx::RVec> f = { {0, 0, 0} };
+ gmx::ForceWithVirial forceWithVirial(f, true);
md.homenr = 1;
snew(md.chargeA, md.homenr);
md.chargeA[0] = 1;
}
if (state_local->flags & (1 << estX))
{
- gmx::ArrayRef<gmx::RVec> globalXRef = state ? gmx::makeArrayRef(state->x) : gmx::EmptyArrayRef();
- dd_collect_vec(dd, state_local, state_local->x, globalXRef);
+ gmx::ArrayRef<gmx::RVec> globalXRef = state ? makeArrayRef(state->x) : gmx::EmptyArrayRef();
+ dd_collect_vec(dd, state_local, makeConstArrayRef(state_local->x), globalXRef);
}
if (state_local->flags & (1 << estV))
{
- gmx::ArrayRef<gmx::RVec> globalVRef = state ? gmx::makeArrayRef(state->v) : gmx::EmptyArrayRef();
- dd_collect_vec(dd, state_local, state_local->v, globalVRef);
+ gmx::ArrayRef<gmx::RVec> globalVRef = state ? makeArrayRef(state->v) : gmx::EmptyArrayRef();
+ dd_collect_vec(dd, state_local, makeConstArrayRef(state_local->v), globalVRef);
}
if (state_local->flags & (1 << estCGP))
{
- gmx::ArrayRef<gmx::RVec> globalCgpRef = state ? gmx::makeArrayRef(state->cg_p) : gmx::EmptyArrayRef();
- dd_collect_vec(dd, state_local, state_local->cg_p, globalCgpRef);
+ gmx::ArrayRef<gmx::RVec> globalCgpRef = state ? makeArrayRef(state->cg_p) : gmx::EmptyArrayRef();
+ dd_collect_vec(dd, state_local, makeConstArrayRef(state_local->cg_p), globalCgpRef);
}
}
static void dd_distribute_state(gmx_domdec_t *dd,
const t_state *state, t_state *state_local,
- PaddedRVecVector *f)
+ PaddedVector<gmx::RVec> *f)
{
int nh = state_local->nhchainlength;
}
}
-void distributeState(const gmx::MDLogger &mdlog,
- gmx_domdec_t *dd,
- const gmx_mtop_t &mtop,
- t_state *state_global,
- const gmx_ddbox_t &ddbox,
- t_state *state_local,
- PaddedRVecVector *f)
+void distributeState(const gmx::MDLogger &mdlog,
+ gmx_domdec_t *dd,
+ const gmx_mtop_t &mtop,
+ t_state *state_global,
+ const gmx_ddbox_t &ddbox,
+ t_state *state_local,
+ PaddedVector<gmx::RVec> *f)
{
- rvec *xGlobal = (DDMASTER(dd) ? as_rvec_array(state_global->x.data()) : nullptr);
+ rvec *xGlobal = (DDMASTER(dd) ? state_global->x.rvec_array() : nullptr);
distributeAtomGroups(mdlog, dd, mtop,
DDMASTER(dd) ? state_global->box : nullptr,
}
/*! \brief Distributes the state from the master rank to all DD ranks */
-void distributeState(const gmx::MDLogger &mdlog,
- gmx_domdec_t *dd,
- const gmx_mtop_t &mtop,
- t_state *state_global,
- const gmx_ddbox_t &ddbox,
- t_state *state_local,
- PaddedRVecVector *f);
+void distributeState(const gmx::MDLogger &mdlog,
+ gmx_domdec_t *dd,
+ const gmx_mtop_t &mtop,
+ t_state *state_global,
+ const gmx_ddbox_t &ddbox,
+ t_state *state_local,
+ PaddedVector<gmx::RVec> *f);
#endif
print_missing_interactions_atoms(mdlog, cr, top_global, &top_local->idef);
write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr,
- -1, as_rvec_array(state_local->x.data()), state_local->box);
+ -1, state_local->x.rvec_array(), state_local->box);
std::string errorMessage;
static void setup_dd_communication(gmx_domdec_t *dd,
matrix box, gmx_ddbox_t *ddbox,
t_forcerec *fr,
- t_state *state, PaddedRVecVector *f)
+ t_state *state,
+ PaddedVector<gmx::RVec> *f)
{
int dim_ind, dim, dim0, dim1, dim2, dimd, nat_tot;
int nzone, nzone_send, zone, zonei, cg0, cg1;
cg_cm = fr->cg_cm;
break;
case ecutsVERLET:
- cg_cm = as_rvec_array(state->x.data());
+ cg_cm = state->x.rvec_array();
break;
default:
gmx_incons("unimplemented");
}
else
{
- cg_cm = as_rvec_array(state->x.data());
+ cg_cm = state->x.rvec_array();
}
/* Communicate cg_cm */
gmx::ArrayRef<gmx::RVec> rvecBufferRef;
}
// TODO Remove fplog when group scheme and charge groups are gone
-void dd_partition_system(FILE *fplog,
- const gmx::MDLogger &mdlog,
- int64_t step,
- const t_commrec *cr,
- gmx_bool bMasterState,
- int nstglobalcomm,
- t_state *state_global,
- const gmx_mtop_t *top_global,
- const t_inputrec *ir,
- t_state *state_local,
- PaddedRVecVector *f,
- gmx::MDAtoms *mdAtoms,
- gmx_localtop_t *top_local,
- t_forcerec *fr,
- gmx_vsite_t *vsite,
- gmx::Constraints *constr,
- t_nrnb *nrnb,
- gmx_wallcycle *wcycle,
- gmx_bool bVerbose)
+void dd_partition_system(FILE *fplog,
+ const gmx::MDLogger &mdlog,
+ int64_t step,
+ const t_commrec *cr,
+ gmx_bool bMasterState,
+ int nstglobalcomm,
+ t_state *state_global,
+ const gmx_mtop_t *top_global,
+ const t_inputrec *ir,
+ t_state *state_local,
+ PaddedVector<gmx::RVec> *f,
+ gmx::MDAtoms *mdAtoms,
+ gmx_localtop_t *top_local,
+ t_forcerec *fr,
+ gmx_vsite_t *vsite,
+ gmx::Constraints *constr,
+ t_nrnb *nrnb,
+ gmx_wallcycle *wcycle,
+ gmx_bool bVerbose)
{
gmx_domdec_t *dd;
gmx_domdec_comm_t *comm;
if (fr->cutoff_scheme == ecutsGROUP)
{
calc_cgcm(fplog, 0, dd->ncg_home,
- &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
+ &top_local->cgs, state_local->x.rvec_array(), fr->cg_cm);
}
inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
{
/* Redetermine the cg COMs */
calc_cgcm(fplog, 0, dd->ncg_home,
- &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
+ &top_local->cgs, state_local->x.rvec_array(), fr->cg_cm);
}
inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
/*
write_dd_pdb("dd_home",step,"dump",top_global,cr,
- -1,as_rvec_array(state_local->x.data()),state_local->box);
+ -1,state_local->x.rvec_array(),state_local->box);
*/
wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
comm->cellsize_min, np,
fr,
- fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : as_rvec_array(state_local->x.data()),
+ fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x.rvec_array(),
vsite, top_global, top_local);
wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
* the last vsite construction, we need to communicate the constructing
* atom coordinates again (for spreading the forces this MD step).
*/
- dd_move_x_vsites(dd, state_local->box, as_rvec_array(state_local->x.data()));
+ dd_move_x_vsites(dd, state_local->box, state_local->x.rvec_array());
wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
{
dd_move_x(dd, state_local->box, state_local->x, nullWallcycle);
write_dd_pdb("dd_dump", step, "dump", top_global, cr,
- -1, as_rvec_array(state_local->x.data()), state_local->box);
+ -1, state_local->x.rvec_array(), state_local->box);
}
/* Store the partitioning step */
* else state_local is redistributed between the nodes.
* When f!=NULL, *f will be reallocated to the size of state_local.
*/
-void dd_partition_system(FILE *fplog,
- const gmx::MDLogger &mdlog,
- int64_t step,
- const t_commrec *cr,
- gmx_bool bMasterState,
- int nstglobalcomm,
- t_state *state_global,
- const gmx_mtop_t *top_global,
- const t_inputrec *ir,
- t_state *state_local,
- PaddedRVecVector *f,
- gmx::MDAtoms *mdatoms,
- gmx_localtop_t *top_local,
- t_forcerec *fr,
- gmx_vsite_t *vsite,
- gmx::Constraints *constr,
- t_nrnb *nrnb,
- gmx_wallcycle *wcycle,
- gmx_bool bVerbose);
+void dd_partition_system(FILE *fplog,
+ const gmx::MDLogger &mdlog,
+ int64_t step,
+ const t_commrec *cr,
+ gmx_bool bMasterState,
+ int nstglobalcomm,
+ t_state *state_global,
+ const gmx_mtop_t *top_global,
+ const t_inputrec *ir,
+ t_state *state_local,
+ PaddedVector<gmx::RVec> *f,
+ gmx::MDAtoms *mdatoms,
+ gmx_localtop_t *top_local,
+ t_forcerec *fr,
+ gmx_vsite_t *vsite,
+ gmx::Constraints *constr,
+ t_nrnb *nrnb,
+ gmx_wallcycle *wcycle,
+ gmx_bool bVerbose);
/*! \brief Check whether bonded interactions are missing, if appropriate
*
{
if (state->flags & (1 << estX))
{
+ auto x = makeArrayRef(state->x);
/* Rotate the complete state; for a rectangular box only */
- state->x[a][YY] = state->box[YY][YY] - state->x[a][YY];
- state->x[a][ZZ] = state->box[ZZ][ZZ] - state->x[a][ZZ];
+ x[a][YY] = state->box[YY][YY] - x[a][YY];
+ x[a][ZZ] = state->box[ZZ][ZZ] - x[a][ZZ];
}
if (state->flags & (1 << estV))
{
- state->v[a][YY] = -state->v[a][YY];
- state->v[a][ZZ] = -state->v[a][ZZ];
+ auto v = makeArrayRef(state->v);
+ v[a][YY] = -v[a][YY];
+ v[a][ZZ] = -v[a][ZZ];
}
if (state->flags & (1 << estCGP))
{
- state->cg_p[a][YY] = -state->cg_p[a][YY];
- state->cg_p[a][ZZ] = -state->cg_p[a][ZZ];
+ auto cg_p = makeArrayRef(state->cg_p);
+ cg_p[a][YY] = -cg_p[a][YY];
+ cg_p[a][ZZ] = -cg_p[a][ZZ];
}
}
gmx::ArrayRef<int> move)
{
const int npbcdim = dd->npbcdim;
+ auto x = makeArrayRef(state->x);
for (int g = cg_start; g < cg_end; g++)
{
rvec cm_new;
if (numAtoms == 1)
{
- copy_rvec(state->x[atomGroup.begin()], cm_new);
+ copy_rvec(x[atomGroup.begin()], cm_new);
}
else
{
clear_rvec(cm_new);
for (int k : atomGroup)
{
- rvec_inc(cm_new, state->x[k]);
+ rvec_inc(cm_new, x[k]);
}
for (int d = 0; d < DIM; d++)
{
if (pos_d >= moveLimits.upper[d])
{
cg_move_error(fplog, dd, step, g, d, 1,
- cg_cm != as_rvec_array(state->x.data()), moveLimits.distance[d],
+ cg_cm != state->x.rvec_array(), moveLimits.distance[d],
cg_cm[g], cm_new, pos_d);
}
dev[d] = 1;
}
for (int k : atomGroup)
{
- rvec_dec(state->x[k], state->box[d]);
+ rvec_dec(x[k], state->box[d]);
if (bScrew)
{
rotate_state_atom(state, k);
if (pos_d < moveLimits.lower[d])
{
cg_move_error(fplog, dd, step, g, d, -1,
- cg_cm != as_rvec_array(state->x.data()), moveLimits.distance[d],
+ cg_cm != state->x.rvec_array(), moveLimits.distance[d],
cg_cm[g], cm_new, pos_d);
}
dev[d] = -1;
}
for (int k : atomGroup)
{
- rvec_inc(state->x[k], state->box[d]);
+ rvec_inc(x[k], state->box[d]);
if (bScrew)
{
rotate_state_atom(state, k);
rvec_dec(cm_new, state->box[d]);
for (int k : atomGroup)
{
- rvec_dec(state->x[k], state->box[d]);
+ rvec_dec(x[k], state->box[d]);
}
}
while (cm_new[d] < 0)
rvec_inc(cm_new, state->box[d]);
for (int k : atomGroup)
{
- rvec_inc(state->x[k], state->box[d]);
+ rvec_inc(x[k], state->box[d]);
}
}
}
void dd_redistribute_cg(FILE *fplog, int64_t step,
gmx_domdec_t *dd, ivec tric_dir,
- t_state *state, PaddedRVecVector *f,
+ t_state *state,
+ PaddedVector<gmx::RVec> *f,
t_forcerec *fr,
t_nrnb *nrnb,
int *ncg_moved)
* many conditionals for both for with and without charge groups.
*/
copyMovedAtomsToBufferPerChargeGroup(move, dd->atomGrouping(),
- nvec, as_rvec_array(state->x.data()), comm);
+ nvec, state->x.rvec_array(), comm);
break;
default:
gmx_incons("unimplemented");
int vectorIndex = 0;
copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
nvec, vectorIndex++,
- as_rvec_array(state->x.data()),
+ state->x.rvec_array(),
comm);
if (bV)
{
copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
nvec, vectorIndex++,
- as_rvec_array(state->v.data()),
+ state->v.rvec_array(),
comm);
}
if (bCGP)
{
copyMovedAtomsToBufferPerAtom(move, dd->atomGrouping(),
nvec, vectorIndex++,
- as_rvec_array(state->cg_p.data()),
+ state->cg_p.rvec_array(),
comm);
}
comm->bLocalCG[globalAtomGroupIndex] = TRUE;
}
+ auto x = makeArrayRef(state->x);
+ auto v = makeArrayRef(state->v);
+ auto cg_p = makeArrayRef(state->cg_p);
rvec *rvecPtr = as_rvec_array(rvecBuffer.buffer.data());
for (int i = 0; i < nrcg; i++)
{
copy_rvec(rvecPtr[buf_pos++],
- state->x[home_pos_at+i]);
+ x[home_pos_at+i]);
}
if (bV)
{
for (int i = 0; i < nrcg; i++)
{
copy_rvec(rvecPtr[buf_pos++],
- state->v[home_pos_at+i]);
+ v[home_pos_at+i]);
}
}
if (bCGP)
for (int i = 0; i < nrcg; i++)
{
copy_rvec(rvecPtr[buf_pos++],
- state->cg_p[home_pos_at+i]);
+ cg_p[home_pos_at+i]);
}
}
home_pos_cg += 1;
class t_state;
/*! \brief Redistribute the atoms to their, new, local domains */
-void dd_redistribute_cg(FILE *fplog,
- int64_t step,
- gmx_domdec_t *dd,
- ivec tric_dir,
- t_state *state,
- PaddedRVecVector *f,
- t_forcerec *fr,
- t_nrnb *nrnb,
- int *ncg_moved);
+void dd_redistribute_cg(FILE *fplog,
+ int64_t step,
+ gmx_domdec_t *dd,
+ ivec tric_dir,
+ t_state *state,
+ PaddedVector<gmx::RVec> *f,
+ t_forcerec *fr,
+ t_nrnb *nrnb,
+ int *ncg_moved);
#endif
}
}
-void dd_resize_state(t_state *state,
- PaddedRVecVector *f,
- int natoms)
+void dd_resize_state(t_state *state,
+ PaddedVector<gmx::RVec> *f,
+ int natoms)
{
if (debug)
{
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- f->resize(gmx::paddedRVecVectorSize(natoms));
+ f->resizeWithPadding(natoms);
}
}
-void dd_check_alloc_ncg(t_forcerec *fr,
- t_state *state,
- PaddedRVecVector *f,
- int numChargeGroups)
+void dd_check_alloc_ncg(t_forcerec *fr,
+ t_state *state,
+ PaddedVector<gmx::RVec> *f,
+ int numChargeGroups)
{
if (numChargeGroups > fr->cg_nalloc)
{
}
/*! \brief Resize the state and f, if !=nullptr, to natoms */
-void dd_resize_state(t_state *state,
- PaddedRVecVector *f,
- int natoms);
+void dd_resize_state(t_state *state,
+ PaddedVector<gmx::RVec> *f,
+ int natoms);
/*! \brief Enrsure fr, state and f, if != nullptr, can hold numChargeGroups atoms for the Verlet scheme and charge groups for the group scheme */
-void dd_check_alloc_ncg(t_forcerec *fr,
- t_state *state,
- PaddedRVecVector *f,
- int numChargeGroups);
+void dd_check_alloc_ncg(t_forcerec *fr,
+ t_state *state,
+ PaddedVector<gmx::RVec> *f,
+ int numChargeGroups);
/*! \brief Returns a domain-to-domain cutoff distance given an atom-to-atom cutoff */
static inline real
{
/* Remove PBC, make molecule(s) subject to ED whole. */
snew(x_pbc, mtop->natoms);
- copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
+ copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
}
/* Reset pointer to first ED data set which contains the actual ED data */
GMX_ASSERT(newForcesSize > 0, "Bad number of atoms in PME GPU");
reallocateDeviceBuffer(&pmeGpu->kernelParams->atoms.d_forces, newForcesSize,
&pmeGpu->archSpecific->forcesSize, &pmeGpu->archSpecific->forcesSizeAlloc, pmeGpu->archSpecific->context);
- pmeGpu->staging.h_forces.reserve(pmeGpu->nAtomsAlloc);
- pmeGpu->staging.h_forces.resize(pmeGpu->kernelParams->atoms.nAtoms);
+ pmeGpu->staging.h_forces.reserveWithPadding(pmeGpu->nAtomsAlloc);
+ pmeGpu->staging.h_forces.resizeWithPadding(pmeGpu->kernelParams->atoms.nAtoms);
}
void pme_gpu_free_forces(const PmeGpu *pmeGpu)
struct PmeGpuStaging
{
//! Host-side force buffer
- std::vector < gmx::RVec, gmx::HostAllocator < gmx::RVec>> h_forces;
+ gmx::HostVector<gmx::RVec> h_forces;
/*! \brief Virial and energy intermediate host-side buffer. Size is PME_GPU_VIRIAL_AND_ENERGY_COUNT. */
float *h_virialAndEnergy;
if (cnb.flags & PP_PME_CHARGE)
{
- pme_pp->chargeA.resize(nat);
+ pme_pp->chargeA.resizeWithPadding(nat);
}
if (cnb.flags & PP_PME_CHARGEB)
{
{
pme_pp->sigmaB.resize(nat);
}
- pme_pp->x.resize(nat);
+ pme_pp->x.resizeWithPadding(nat);
pme_pp->f.resize(nat);
/* maxshift is sent when the charges are sent */
{
if (sender.numAtoms > 0)
{
- MPI_Irecv(pme_pp->x[nat], sender.numAtoms*sizeof(rvec),
+ MPI_Irecv(pme_pp->x[nat],
+ sender.numAtoms*sizeof(rvec),
MPI_BYTE,
sender.rankId, eCommType_COORD,
pme_pp->mpi_comm_mysim, &pme_pp->req[messages++]);
//TODO this should be set properly by gmx_pme_recv_coeffs_coords,
// or maybe use inputrecDynamicBox(ir), at the very least - change this when this codepath is tested!
pme_gpu_prepare_computation(pme, boxChanged, box, wcycle, pmeFlags);
- pme_gpu_launch_spread(pme, as_rvec_array(pme_pp->x.data()), wcycle);
+ pme_gpu_launch_spread(pme, pme_pp->x.rvec_array(), wcycle);
pme_gpu_launch_complex_transforms(pme, wcycle);
pme_gpu_launch_gather(pme, wcycle, PmeForceOutputHandling::Set);
pme_gpu_wait_finish_task(pme, wcycle, &forces, vir_q, &energy_q);
}
else
{
- gmx_pme_do(pme, 0, natoms, as_rvec_array(pme_pp->x.data()), as_rvec_array(pme_pp->f.data()),
+ gmx_pme_do(pme, 0, natoms, pme_pp->x.rvec_array(), as_rvec_array(pme_pp->f.data()),
pme_pp->chargeA.data(), pme_pp->chargeB.data(),
pme_pp->sqrt_c6A.data(), pme_pp->sqrt_c6B.data(),
pme_pp->sigmaA.data(), pme_pp->sigmaB.data(), box,
{
/* This conditional ensures that we don't resize on write.
* In particular in the state where this code was written
- * PaddedRVecVector has a size of numElemInThefile and we
+ * vector has a size of numElemInThefile and we
* don't want to lose that padding here.
*/
if (vector->size() < static_cast<unsigned int>(numElemInTheFile))
return doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, vector.size(), nullptr, &v_real, nullptr, list, CptElementType::real);
}
-//! \brief Read/Write a vector whose value_type is RVec and whose allocator is \c AllocatorType.
-template <typename AllocatorType>
+//! Convert from view of RVec to view of real.
+static gmx::ArrayRef<real>
+realArrayRefFromRVecArrayRef(gmx::ArrayRef<gmx::RVec> ofRVecs)
+{
+ return gmx::arrayRefFromArray<real>(reinterpret_cast<real *>(ofRVecs.data()), ofRVecs.size() * DIM);
+}
+
+//! \brief Read/Write a PaddedVector whose value_type is RVec.
+template <typename PaddedVectorOfRVecType>
static int doRvecVector(XDR *xd, StatePart part, int ecpt, int sflags,
- std::vector<gmx::RVec, AllocatorType> *v, int numAtoms, FILE *list)
+ PaddedVectorOfRVecType *v, int numAtoms, FILE *list)
{
const int numReals = numAtoms*DIM;
if (list == nullptr)
{
GMX_RELEASE_ASSERT(sflags & (1 << ecpt), "When not listing, the flag for the entry should be set when requesting i/o");
- GMX_RELEASE_ASSERT(v->size() >= static_cast<size_t>(numAtoms), "v should have sufficient size for numAtoms");
-
- real *v_real = (*v)[0];
-
- // PaddedRVecVector is padded beyond numAtoms, we should only write
- // numAtoms RVecs
- gmx::ArrayRef<real> ref(v_real, v_real + numReals);
+ GMX_RELEASE_ASSERT(v->size() == numAtoms, "v should have sufficient size for numAtoms");
- return doRealArrayRef(xd, part, ecpt, sflags, ref, list);
+ return doRealArrayRef(xd, part, ecpt, sflags, realArrayRefFromRVecArrayRef(makeArrayRef(*v)), list);
}
else
{
// Use the rebind facility to change the value_type of the
// allocator from RVec to real.
- using realAllocator = typename std::allocator_traits<AllocatorType>::template rebind_alloc<real>;
+ using realAllocator = typename std::allocator_traits<typename PaddedVectorOfRVecType::allocator_type>::template rebind_alloc<real>;
return doVectorLow<real, realAllocator>(xd, part, ecpt, sflags, numReals, nullptr, nullptr, nullptr, list, CptElementType::real);
}
}
if (x == nullptr)
{
- x = as_rvec_array(state->x.data());
- v = as_rvec_array(state->v.data());
+ x = state->x.rvec_array();
+ v = state->v.rvec_array();
}
#define do_test(fio, b, p) if (bRead && ((p) != NULL) && !(b)) gmx_fatal(FARGS, "No %s in %s",#p, gmx_fio_getname(fio))
}
/* Prepare an x and q array with only the charged atoms */
- ncharges = prepare_x_q(&q, &x, mtop, as_rvec_array(state->x.data()), cr);
+ ncharges = prepare_x_q(&q, &x, mtop, state->x.rvec_array(), cr);
if (MASTER(cr))
{
calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
state->flags |= (1 << estV);
}
state_change_natoms(state, state->natoms);
- for (int i = 0; i < state->natoms; i++)
- {
- copy_rvec(x[i], state->x[i]);
- }
+ std::copy(x, x+state->natoms, state->x.data());
sfree(x);
if (v != nullptr)
{
- for (int i = 0; i < state->natoms; i++)
- {
- copy_rvec(v[i], state->v[i]);
- }
+ std::copy(v, v+state->natoms, state->v.data());
sfree(v);
}
/* This call fixes the box shape for runs with pressure scaling */
if (ir->cutoff_scheme == ecutsGROUP && ir->ePBC != epbcNONE)
{
// Need temporary rvec for coordinates
- do_pbc_first_mtop(nullptr, ir->ePBC, state->box, sys, as_rvec_array(state->x.data()));
+ do_pbc_first_mtop(nullptr, ir->ePBC, state->box, sys, state->x.rvec_array());
}
/* Do more checks, mostly related to constraints */
fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
}
state->flags |= (1 << estV);
- maxwell_speed(opts->tempi, opts->seed, sys, as_rvec_array(state->v.data()));
+ maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array());
- stop_cm(stdout, state->natoms, mass, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()));
+ stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
sfree(mass);
}
bool bReadVel, t_state *state,
double *use_time)
{
- int i;
-
if (fr->not_ok & FRAME_NOT_OK)
{
gmx_fatal(FARGS, "Can not start from an incomplete frame");
slog);
}
- for (i = 0; i < state->natoms; i++)
- {
- copy_rvec(fr->x[i], state->x[i]);
- }
+ std::copy(fr->x, fr->x+state->natoms, state->x.data());
if (bReadVel)
{
if (!fr->bV)
{
gmx_incons("Trajecory frame unexpectedly does not contain velocities");
}
- for (i = 0; i < state->natoms; i++)
- {
- copy_rvec(fr->v[i], state->v[i]);
- }
+ std::copy(fr->v, fr->v+state->natoms, state->v.data());
}
if (fr->bBox)
{
bool bReadVel;
t_trxframe fr;
t_trxstatus *fp;
- int i;
double use_time;
bReadVel = (bNeedVel && !bGenVel);
"\n"
"WARNING: Did not find a frame with velocities in file %s,\n"
" all velocities will be set to zero!\n\n", slog);
- for (i = 0; i < sys->natoms; i++)
+ for (auto &vi : makeArrayRef(state->v))
{
- clear_rvec(state->v[i]);
+ vi = {0, 0, 0};
}
close_trx(fp);
/* Search for a frame without velocities */
/* Check velocity for virtual sites and shells */
if (bGenVel)
{
- check_vel(&sys, as_rvec_array(state.v.data()));
+ check_vel(&sys, state.v.rvec_array());
}
/* check for shells and inpurecs */
}
else
{
- buffer_temp = calc_temp(&sys, ir, as_rvec_array(state.v.data()));
+ buffer_temp = calc_temp(&sys, ir, state.v.rvec_array());
}
if (buffer_temp > 0)
{
if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
{
set_warning_line(wi, mdparin, -1);
- check_chargegroup_radii(&sys, ir, as_rvec_array(state.x.data()), wi);
+ check_chargegroup_radii(&sys, ir, state.x.rvec_array(), wi);
}
if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
if (ir->bPull)
{
- pull = set_pull_init(ir, &sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS], wi);
+ pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
}
/* Modules that supply external potential for pull coordinates
if (ir->bRot)
{
- set_reference_positions(ir->rot, as_rvec_array(state.x.data()), state.box,
+ set_reference_positions(ir->rot, state.x.rvec_array(), state.box,
opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
wi);
}
#include <cstddef>
#include <memory>
-#include <vector>
+#include "gromacs/math/paddedvector.h"
#include "gromacs/utility/alignedallocator.h"
#include "gromacs/utility/exceptions.h"
//! Convenience alias for std::vector that uses HostAllocator.
template <class T>
-using HostVector = std::vector<T, HostAllocator<T> >;
+using HostVector = PaddedVector<T, HostAllocator<T> >;
/*! \libinternal
* \brief Policy class for configuring gmx::Allocator, to manage
* Does not throw.
*/
void free(void *buffer) const noexcept;
- /*! \brief Pin the allocation to physical memory, if appropriate.
- *
- * If the allocation policy is not in pinning mode, or the
- * allocation is empty, ot the allocation is already pinned,
- * then do nothing.
+ /*! \brief Return the active pinning policy.
*
* Does not throw.
*/
//container is forced to realloc). Does nothing if policy is the same.
*v = HostVector<T>(std::move(*v), {pinningPolicy});
}
+
} // namespace gmx
#endif
#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/real.h"
+#include "gromacs/math/tests/testarrayrefs.h"
+
#include "devicetransfers.h"
#include "gputest.h"
namespace gmx
{
-namespace
+namespace test
{
/*! \internal \brief Typed test fixture for infrastructure for
public:
//! Convenience type
using ValueType = T;
- //! Convenience type
- using ViewType = ArrayRef<ValueType>;
- //! Convenience type
- using ConstViewType = ArrayRef<const ValueType>;
- //! Prepare contents of a VectorType.
- template <typename VectorType>
- void fillInput(VectorType *input) const;
- //! Compares input and output vectors.
- void compareVectors(ConstViewType input,
- ConstViewType output) const;
- //! Do some transfers and test the results.
- void runTest(ConstViewType input, ViewType output) const;
};
-// Already documented
-template <typename T> template <typename VectorType>
-void HostMemoryTest<T>::fillInput(VectorType *input) const
-{
- input->resize(3);
- (*input)[0] = 1;
- (*input)[1] = 2;
- (*input)[2] = 3;
-}
-
-//! Initialization specialization for RVec
-template <> template <typename VectorType>
-void HostMemoryTest<RVec>::fillInput(VectorType *input) const
-{
- input->reserve(3);
- input->resize(3);
- (*input)[0] = {1, 2, 3};
- (*input)[1] = {4, 5, 6};
- (*input)[2] = {7, 8, 9};
-}
-
-// Already documented
-template <typename T>
-void HostMemoryTest<T>::compareVectors(ConstViewType input,
- ConstViewType output) const
-{
- for (index i = 0; i != input.size(); ++i)
- {
- EXPECT_EQ(input[i], output[i]) << "for index " << i;
- }
-}
-
-//! Comparison specialization for RVec
-template <>
-void HostMemoryTest<RVec>::compareVectors(ConstViewType input,
- ConstViewType output) const
-{
- for (index i = 0; i != input.size(); ++i)
- {
- EXPECT_EQ(input[i][XX], output[i][XX]) << "for index " << i;
- EXPECT_EQ(input[i][YY], output[i][YY]) << "for index " << i;
- EXPECT_EQ(input[i][ZZ], output[i][ZZ]) << "for index " << i;
- }
-}
-
/*! \brief Convenience function to transform a view into one with base
* type of (non-const) char.
*
return arrayRefFromArray<char>(reinterpret_cast<char *>(const_cast<NonConstT *>(data)), size * sizeof(T));
}
+//! Does a device transfer of \c input to the device in \c gpuInfo, and back to \c output.
template <typename T>
-void HostMemoryTest<T>::runTest(ConstViewType input, ViewType output) const
+void runTest(const gmx_gpu_info_t &gpuInfo,
+ ArrayRef<T> input,
+ ArrayRef<T> output)
{
// Convert the views of input and output to flat non-const chars,
// so that there's no templating when we call doDeviceTransfers.
auto inputRef = charArrayRefFromArray(input.data(), input.size());
auto outputRef = charArrayRefFromArray(output.data(), output.size());
- doDeviceTransfers(*this->gpuInfo_, inputRef, outputRef);
- this->compareVectors(input, output);
+ ASSERT_EQ(inputRef.size(), outputRef.size());
+ doDeviceTransfers(gpuInfo, inputRef, outputRef);
+ compareViews(input, output);
}
struct MoveOnly {
MoveOnly &operator=(const MoveOnly &) = delete;
MoveOnly &operator=(MoveOnly &&) = default;
bool operator==(const MoveOnly &o) const { return x == o.x; }
+ real operator*=(int scaleFactor) { return x *= scaleFactor; }
real x;
};
+} // namespace test
+
+namespace detail
+{
+
+template<>
+struct PaddingTraits<test::MoveOnly>
+{
+ using SimdBaseType = real;
+ static constexpr int maxSimdWidthOfBaseType = GMX_REAL_MAX_SIMD_WIDTH;
+};
+
+} // namespace detail
+
+namespace test
+{
+
//! The types used in testing.
-typedef ::testing::Types<int, real, RVec, MoveOnly> TestTypes;
+typedef ::testing::Types<int32_t, real, RVec, test::MoveOnly> TestTypes;
//! Typed test fixture
template <typename T>
template <typename T>
struct HostAllocatorTestNoMemCopyable : HostAllocatorTestNoMem<T> {};
//! The types used in testing minus move only types
-using TestTypesCopyable = ::testing::Types<int, real, RVec>;
+using TestTypesCopyable = ::testing::Types<int32_t, real, RVec>;
TYPED_TEST_CASE(HostAllocatorTestNoMemCopyable, TestTypesCopyable);
// Note that in GoogleTest typed tests, the use of TestFixture:: and
TYPED_TEST(HostAllocatorTest, VectorsWithDefaultHostAllocatorAlwaysWorks)
{
typename TestFixture::VectorType input(3), output;
- output.resize(input.size());
+ output.resizeWithPadding(input.size());
}
// Several tests actually do CUDA transfers. This is not necessary
TYPED_TEST(HostAllocatorTest, TransfersWithoutPinningWork)
{
typename TestFixture::VectorType input;
- this->fillInput(&input);
+ fillInput(&input, 1);
typename TestFixture::VectorType output;
- output.resize(input.size());
+ output.resizeWithPadding(input.size());
- this->runTest(input, output);
+ runTest(*this->gpuInfo_, makeArrayRef(input), makeArrayRef(output));
}
TYPED_TEST(HostAllocatorTest, FillInputAlsoWorksAfterCallingReserve)
{
typename TestFixture::VectorType input;
- input.reserve(3);
- this->fillInput(&input);
+ input.reserveWithPadding(3);
+ fillInput(&input, 1);
}
TYPED_TEST(HostAllocatorTestNoMem, CreateVector)
{
typename TestFixture::VectorType input1;
typename TestFixture::VectorType input2({PinningPolicy::PinnedIfSupported});
- swap(input1, input2);
+ std::swap(input1, input2);
EXPECT_TRUE (input1.get_allocator().pinningPolicy() == PinningPolicy::PinnedIfSupported);
EXPECT_FALSE(input2.get_allocator().pinningPolicy() == PinningPolicy::PinnedIfSupported);
- swap(input2, input1);
+ std::swap(input2, input1);
EXPECT_FALSE(input1.get_allocator().pinningPolicy() == PinningPolicy::PinnedIfSupported);
EXPECT_TRUE (input2.get_allocator().pinningPolicy() == PinningPolicy::PinnedIfSupported);
}
typename TestFixture::VectorType input;
changePinningPolicy(&input, PinningPolicy::PinnedIfSupported);
- this->fillInput(&input);
+ fillInput(&input, 1);
typename TestFixture::VectorType output;
changePinningPolicy(&output, PinningPolicy::PinnedIfSupported);
- output.resize(input.size());
+ output.resizeWithPadding(input.size());
- this->runTest(input, output);
+ runTest(*this->gpuInfo_, makeArrayRef(input), makeArrayRef(output));
}
//! Helper function for wrapping a call to isHostMemoryPinned.
EXPECT_FALSE(isPinned(input));
// Fill some contents, which will be pinned because of the policy.
- this->fillInput(&input);
+ fillInput(&input, 1);
EXPECT_TRUE(isPinned(input));
// Switching policy to CannotBePinned must unpin the buffer (via
// The HostAllocator has state, so a container using it will be
// larger than a normal vector, whose default allocator is
// stateless.
- EXPECT_LT(sizeof(std::vector<typename TestFixture::ValueType>),
+ EXPECT_LT(sizeof(std::vector<typename TestFixture::VectorType::value_type>),
sizeof(typename TestFixture::VectorType));
}
//! Declare allocator types to test.
using AllocatorTypesToTest = ::testing::Types<HostAllocator<real>,
- HostAllocator<int>,
+ HostAllocator<int32_t>,
HostAllocator<RVec>,
HostAllocator<MoveOnly>
>;
TYPED_TEST_CASE(AllocatorTest, AllocatorTypesToTest);
-} // namespace
+} // namespace test
} // namespace gmx
// Includes tests common to all allocation policies.
{
IMDatoms = gmx_mtop_global_atoms(sys);
write_sto_conf_indexed(opt2fn("-imd", nfile, fnm), "IMDgroup", &IMDatoms,
- as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
+ state->x.rvec_array(), state->v.rvec_array(), ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
}
}
* Copy it to the other nodes after checking multi compatibility,
* so we are sure the subsystems match before copying.
*/
+ auto x = makeArrayRef(globalState->x);
rvec com = { 0, 0, 0 };
double mtot = 0.0;
int j = 0;
// Note that only one rank per sim is supported.
if (isMasterSim(ms))
{
- copy_rvec(globalState->x[i], od->xref[j]);
+ copy_rvec(x[i], od->xref[j]);
for (int d = 0; d < DIM; d++)
{
- com[d] += od->mref[j]*globalState->x[i][d];
+ com[d] += od->mref[j]*x[i][d];
}
}
mtot += od->mref[j];
#include "gromacs/math/vectypes.h"
#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/allocator.h"
#include "gromacs/utility/arrayref.h"
namespace gmx
{
-/*! \brief Temporary definition of a type usable for SIMD-style loads of RVec quantities.
+/*! \brief Temporary definition of a type usable for SIMD-style loads of RVec quantities from a view.
+ *
+ * \todo Find a more permanent solution that permits the update code to safely
+ * use a padded, aligned array-ref type. */
+template <typename T>
+using PaddedArrayRef = ArrayRef<T>;
+
+namespace detail
+{
+
+/*! \brief Traits classes for handling padding for types used with PaddedVector
*
- * \note When resizing paddedRVecVector, the size should be chosen with
- paddedRVecVectorSize() to ensure correct padding.
- * \todo Consider replacing the padding applied in resizePaddedRVecVector()
- * by automated padding on resize() of the vector.
- * \todo Undo the move of allocator.h and alignedallocator.h from the internal
- * to be public API applied in Change-Id: Ifb8dacf, needed to use
- * AlignedAllocationPolicy here, when replacing std::vector here.
+ * Only the base types of the SIMD module are supported for
+ * PaddedVector, because the purpose of the padding is to permit
+ * SIMD-width operations from the SIMD module.
+ *
+ * \todo Consider explicitly tying these types to the SimdTrait
+ * types. That would require depending on the SIMD module, or
+ * extracting the traits from it. This would also permit
+ * maxSimdWidthOfBaseType to be set more efficiently, e.g. as a
+ * metaprogramming max over the maximum width from different
+ * implementations.
*/
-using PaddedRVecVector = std::vector < RVec, Allocator < RVec, AlignedAllocationPolicy > >;
+template<typename T>
+struct PaddingTraits {};
+
+template<>
+struct PaddingTraits<int32_t>
+{
+ using SimdBaseType = int32_t;
+ static constexpr int maxSimdWidthOfBaseType = 16;
+};
-/*! \brief Returns the padded size for PaddedRVecVector given the number of atoms
+template<>
+struct PaddingTraits<float>
+{
+ using SimdBaseType = float;
+ static constexpr int maxSimdWidthOfBaseType = GMX_FLOAT_MAX_SIMD_WIDTH;
+};
+
+template<>
+struct PaddingTraits<double>
+{
+ using SimdBaseType = double;
+ static constexpr int maxSimdWidthOfBaseType = GMX_DOUBLE_MAX_SIMD_WIDTH;
+};
+
+template<>
+struct PaddingTraits < BasicVector < float>>
+{
+ using SimdBaseType = float;
+ static constexpr int maxSimdWidthOfBaseType = GMX_FLOAT_MAX_SIMD_WIDTH;
+};
+
+template<>
+struct PaddingTraits < BasicVector < double>>
+{
+ using SimdBaseType = double;
+ static constexpr int maxSimdWidthOfBaseType = GMX_DOUBLE_MAX_SIMD_WIDTH;
+};
+
+/*! \brief Returns the allocation size for PaddedVector that contains
+ * \c numElements elements plus padding for SIMD operations.
*
- * \param[in] numAtoms The number of atoms for which data will be stored in a PaddedRVecVector
+ * \param[in] numElements The number of T elements for which data will be stored.
+ * \returns The number of T elements that must be allocated
+ * (ie >= numElements).
*/
-static inline index paddedRVecVectorSize(index numAtoms)
+template <typename T>
+index computePaddedSize(index numElements)
{
- /* We need one real extra for 4-wide SIMD load/store of RVec.
- * But because the vector contains RVecs, we need to add 1 RVec.
- */
- index simdScatterAccessSize;
- if (numAtoms > 0)
+ // We don't need padding if there is no access.
+ if (numElements == 0)
{
- simdScatterAccessSize = numAtoms + 1;
- }
- else
- {
- simdScatterAccessSize = 0;
+ return 0;
}
- /* We need padding up to a multiple of the SIMD real width.
- * But because we don't want a dependence on the SIMD module,
- * we use GMX_REAL_MAX_SIMD_WIDTH here.
- */
- index simdFlatAccesSize = ((numAtoms + (GMX_REAL_MAX_SIMD_WIDTH - 1))/GMX_REAL_MAX_SIMD_WIDTH)*GMX_REAL_MAX_SIMD_WIDTH;
+ // We sometimes load a whole extra element when doing 4-wide SIMD
+ // operations (which might e.g. be an RVec) so we need to pad for
+ // that.
+ index simdScatterAccessSize = numElements + 1;
+
+ // For SIMD updates based on RVec, we might load starting from
+ // the last RVec element, so that sets the minimum extent of the
+ // padding. That extent must take the initialized allocation up to
+ // the SIMD width of the base type multiplied by the width of T in
+ // that base type. But since storage_ contains RVec, we only have
+ // to tell it the number of elements, which means to round up to
+ // the next SIMD width.
+ //
+ // We don't want a dependence on the SIMD module for the actual
+ // SIMD width of the base type, so we use maximum for the base
+ // type via the traits. A little extra padding won't really hurt.
+ constexpr int maxSimdWidth = PaddingTraits<T>::maxSimdWidthOfBaseType;
+ index simdFlatAccessSize = (numElements + (maxSimdWidth-1)) / maxSimdWidth * maxSimdWidth;
- return std::max(simdScatterAccessSize, simdFlatAccesSize);
+ return std::max(simdScatterAccessSize, simdFlatAccessSize);
}
-/*! \brief Temporary definition of a type usable for SIMD-style loads of RVec quantities from a view.
+//! Helper function to insert padding elements for most T.
+template <typename T, typename AllocatorType>
+inline void insertPaddingElements(std::vector<T, AllocatorType> *v,
+ index newPaddedSize)
+{
+ // Ensure the padding region is initialized to zero.
+ v->insert(v->end(), newPaddedSize - v->size());
+}
+
+//! Specialization of helper function to insert padding elements, used for BasicVector<T>.
+template <typename T, typename AllocatorType>
+inline void insertPaddingElements(std::vector<BasicVector<T>, AllocatorType> *v,
+ index newPaddedSize)
+{
+ // Ensure the padding region is initialized to zero.
+ v->insert(v->end(), newPaddedSize - v->size(), BasicVector<T>(0, 0, 0));
+}
+
+} // namespace detail
+
+/*! \brief PaddedVector is a container of elements in contiguous
+ * storage that allocates extra memory for safe SIMD-style loads for
+ * operations used in GROMACS.
*
- * \todo Find a more permanent solution that permits the update code to safely
- * use a padded, aligned array-ref type. */
-template <typename T>
-using PaddedArrayRef = ArrayRef<T>;
+ * \tparam T the type of objects within the container
+ * \tparam Allocator the allocator used. Can be any standard-compliant
+ * allocator, such gmx::Allocator used for alignment and/or pinning.
+ *
+ * The interface resembles std::vector. However, access
+ * intended to include padded elements must be via ArrayRef objects
+ * explicitly created to view those elements. Most other aspects of
+ * this vector refer to the unpadded view, e.g. iterators, data(),
+ * size().
+ *
+ * The underlying storage is allocated with extra elements, properly
+ * initialized, that ensure that any operations accessing the any
+ * non-additional element that operate on memory equivalent to a full
+ * SIMD lane do so on allocated memory that has been initialized, so
+ * that memory traps will not occur, and arithmetic operations will
+ * not cause e.g. floating-point exceptions so long as the values in
+ * the padded elements are properly managed.
+ *
+ * Proper initialization is tricker than it would first appear, since
+ * we intend this container to be used with scalar and class types
+ * (e.g. RVec). Resize and construction operations use "default
+ * insertion" which leads to zero initialization for the former, and
+ * calling the default constructor for the latter. BasicVector has a
+ * default constructor that leaves the elements uninitialized, which
+ * is particularly risky for elements only present as padding. Thus
+ * the implementation specifically initializes the padded elements to
+ * zero, which makes no difference to the scalar template
+ * instantiations, and makes the BasicVector ones safer to use.
+ *
+ * Because the allocator can be configured, the memory allocation can
+ * have other attributes such as SIMD alignment or being pinned to
+ * physical memory for efficient transfers. The default allocator
+ * ensures alignment, but std::allocator also works.
+ */
+template <typename T, typename Allocator = Allocator < T, AlignedAllocationPolicy > >
+class PaddedVector
+{
+ public:
+ //! Standard helper types
+ //! \{
+ using value_type = T;
+ using allocator_type = Allocator;
+ using size_type = index;
+ using reference = value_type &;
+ using const_reference = const value_type &;
+ using storage_type = std::vector<T, allocator_type>;
+ using pointer = typename storage_type::pointer;
+ using const_pointer = typename storage_type::const_pointer;
+ using iterator = typename storage_type::iterator;
+ using const_iterator = typename storage_type::const_iterator;
+ using difference_type = typename storage_type::iterator::difference_type;
+ //! \}
+
+ PaddedVector() :
+ storage_(),
+ unpaddedEnd_(begin())
+ {}
+ /*! \brief Constructor that specifes the initial size.
+ *
+ * \todo This should also be specialized by allocator, but
+ * std::vector for storage_ doesn't have such a constructor
+ * before C++14. Resolve. */
+ explicit PaddedVector(size_type count) :
+ storage_(count),
+ unpaddedEnd_(begin() + count)
+ {
+ // The count elements have been default inserted, and now
+ // the padding elements are added
+ resizeWithPadding(count);
+ }
+ /*! \brief Constructor that specifes the initial size and an element to copy.
+ *
+ * \todo This should also be specialized by allocator, but
+ * std::vector for storage_ doesn't have such a constructor
+ * before C++14. Resolve. */
+ explicit PaddedVector(size_type count, value_type const &v) :
+ storage_(count, v),
+ unpaddedEnd_(begin() + count)
+ {
+ // The count elements have been default inserted, and now
+ // the padding elements are added
+ resizeWithPadding(count);
+ }
+ //! Default constructor with allocator
+ explicit PaddedVector(allocator_type const &allocator) :
+ storage_(allocator),
+ unpaddedEnd_(begin())
+ {}
+ //! Copy constructor
+ PaddedVector(PaddedVector const &o) :
+ storage_(o.storage_),
+ unpaddedEnd_(begin() + o.size())
+ {}
+ //! Move constructor
+ PaddedVector(PaddedVector &&o) noexcept :
+ storage_(std::move(o.storage_)),
+ unpaddedEnd_(std::move(o.unpaddedEnd_))
+ {
+ unpaddedEnd_ = begin();
+ }
+ //! Move constructor using \c alloc for the new vector.
+ PaddedVector(PaddedVector &&o, const Allocator &alloc) noexcept :
+ storage_(std::move(o.storage_), alloc),
+ unpaddedEnd_(std::move(o.unpaddedEnd_))
+ {
+ unpaddedEnd_ = begin();
+ }
+ //! Construct from an initializer list
+ PaddedVector(std::initializer_list<value_type> const &il) :
+ storage_(il),
+ unpaddedEnd_(storage_.end())
+ {
+ // We can't choose the padding until we know the size of
+ // the normal vector, so we have to make the storage_ and
+ // then resize it.
+ resizeWithPadding(storage_.size());
+ }
+ //! Reserve storage for the container to contain newExtent elements, plus the required padding.
+ void reserveWithPadding(const size_type newExtent)
+ {
+ auto unpaddedSize = end() - begin();
+ /* v.reserve(13) should allocate enough memory so that
+ v.resize(13) does not reallocate. This means that the
+ new extent should be large enough for the padded
+ storage for a vector whose size is newExtent. */
+ auto newPaddedExtent = detail::computePaddedSize<T>(newExtent);
+ storage_.reserve(newPaddedExtent);
+ unpaddedEnd_ = begin() + unpaddedSize;
+ }
+ //! Resize the container to contain newSize elements, plus the required padding.
+ void resizeWithPadding(const size_type newSize)
+ {
+ // When the contained type is e.g. a scalar, then the
+ // default initialization behaviour is to zero all
+ // elements, which is OK, but we have to make sure that it
+ // happens for the elements in the padded region when the
+ // vector is shrinking.
+ auto newPaddedSize = detail::computePaddedSize<T>(newSize);
+ // Make sure there is room for padding if we need to grow.
+ storage_.reserve(newPaddedSize);
+ // Make the unpadded size correct, with any additional
+ // elements initialized by the default constructor. It is
+ // particularly important to destruct former elements when
+ // newSize is smaller than the old size.
+ storage_.resize(newSize);
+ // Ensure the padding region is zeroed if required.
+ detail::insertPaddingElements(&storage_, newPaddedSize);
+ unpaddedEnd_ = begin() + newSize;
+ }
+ //! Return the size of the view without the padding.
+ size_type size() const { return end() - begin(); }
+ //! Return the container size including the padding.
+ size_type paddedSize() const { return storage_.size(); }
+ //! Return whether the storage is empty.
+ bool empty() const { return storage_.size() == 0; }
+ //! Swap two PaddedVectors
+ void swap(PaddedVector &x)
+ {
+ std::swap(storage_, x.storage_);
+ std::swap(unpaddedEnd_, x.unpaddedEnd_);
+ }
+ //! Clear the vector, ie. set size to zero and remove padding.
+ void clear()
+ {
+ storage_.clear();
+ unpaddedEnd_ = begin();
+ }
+ //! Iterator getters refer to a view without padding.
+ //! \{
+ pointer data() noexcept { return storage_.data(); }
+ const_pointer data() const noexcept { return storage_.data(); }
+
+ iterator begin() { return storage_.begin(); }
+ iterator end() { return iterator(unpaddedEnd_); }
+
+ const_iterator cbegin() { return const_iterator(begin()); }
+ const_iterator cend() { return const_iterator(unpaddedEnd_); }
+
+ const_iterator begin() const { return storage_.begin(); }
+ const_iterator end() const { return const_iterator(unpaddedEnd_); }
+
+ const_iterator cbegin() const { return const_iterator(begin()); }
+ const_iterator cend() const { return const_iterator(unpaddedEnd_); }
+ //! \}
+ // TODO should these do bounds checking for the unpadded range? In debug mode?
+ //! Indexing operator.
+ reference operator[](int i) { return storage_[i]; }
+ //! Indexing operator as const.
+ const_reference operator[](int i) const { return storage_[i]; }
+ //! Returns an ArrayRef of elements that includes the padding region, e.g. for use in SIMD code.
+ PaddedArrayRef<T> paddedArrayRef()
+ {
+ return PaddedArrayRef<T>(storage_);
+ }
+ //! Returns an ArrayRef of const elements that includes the padding region, e.g. for use in SIMD code.
+ PaddedArrayRef<const T> paddedConstArrayRef() const
+ {
+ return PaddedArrayRef<const T>(storage_);
+ }
+ //! Returns an rvec * pointer for containers of RVec, for use with legacy code.
+ template <typename AlsoT = T,
+ typename = typename std::enable_if<std::is_same<AlsoT, RVec>::value> >
+ rvec *rvec_array()
+ {
+ return as_rvec_array(data());
+ }
+ //! Returns a const rvec * pointer for containers of RVec, for use with legacy code.
+ template <typename AlsoT = T,
+ typename = typename std::enable_if<std::is_same<AlsoT, RVec>::value> >
+ const rvec *rvec_array() const
+ {
+ return as_rvec_array(data());
+ }
+ //! Copy assignment operator
+ PaddedVector &operator=(PaddedVector const &o)
+ {
+ if (&o != this)
+ {
+ storage_ = o.storage_;
+ unpaddedEnd_ = begin() + o.size();
+ }
+ return *this;
+ }
+ //! Move assignment operator
+ PaddedVector &operator=(PaddedVector &&o) noexcept
+ {
+ if (&o != this)
+ {
+ auto oSize = o.size();
+ storage_ = std::move(o.storage_);
+ unpaddedEnd_ = begin() + oSize;
+ o.unpaddedEnd_ = o.begin();
+ }
+ return *this;
+ }
+ //! Getter for the allocator
+ allocator_type
+ get_allocator() const
+ {
+ return storage_.get_allocator();
+ }
+
+ private:
+ storage_type storage_;
+ iterator unpaddedEnd_;
+};
} // namespace gmx
// TODO These are hacks to avoid littering gmx:: all over code that is
// almost all destined to move into the gmx namespace at some point.
// An alternative would be about 20 files with using statements.
-using gmx::PaddedRVecVector; //NOLINT(google-global-names-in-headers)
+using gmx::PaddedVector; //NOLINT(google-global-names-in-headers)
#endif
dofit.cpp
functions.cpp
invertmatrix.cpp
+ paddedvector.cpp
vectypes.cpp
)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests implementation of gmx::PaddedVector
+ *
+ * In particular we need to check that the padding works, and that
+ * unpadded_end() is maintained correctly.
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_math
+ */
+#include "gmxpre.h"
+
+#include "gromacs/math/paddedvector.h"
+
+#include <memory>
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/basedefinitions.h"
+
+#include "gromacs/math/tests/testarrayrefs.h"
+
+namespace gmx
+{
+namespace test
+{
+
+//! Typed test fixture
+template <typename T>
+class PaddedVectorTest : public ::testing::Test
+{
+ public:
+};
+
+//! The types used in testing
+using Implementations = ::testing::Types <
+ std::allocator<int32_t>,
+ std::allocator<float>,
+ std::allocator<double>,
+ std::allocator < BasicVector < float>>,
+ std::allocator < BasicVector < double>>,
+ AlignedAllocator<int32_t>,
+ AlignedAllocator<float>,
+ AlignedAllocator<double>,
+ AlignedAllocator < BasicVector < float>>,
+ AlignedAllocator < BasicVector<double>>
+ >;
+TYPED_TEST_CASE(PaddedVectorTest, Implementations);
+
+TYPED_TEST(PaddedVectorTest, ConstructsResizesAndReserves)
+{
+ using VectorType = PaddedVector<typename TypeParam::value_type, TypeParam>;
+
+ VectorType v;
+ fillInput(&v, 1);
+
+ EXPECT_EQ(v.size(), v.size());
+ EXPECT_EQ(v.paddedSize(), v.paddedArrayRef().size());
+ EXPECT_LE(v.size(), v.paddedArrayRef().size());
+
+ VectorType vReserved;
+ vReserved.reserveWithPadding(5);
+ fillInput(&vReserved, 1);
+
+ EXPECT_EQ(vReserved.size(), vReserved.size());
+ EXPECT_EQ(vReserved.paddedSize(), vReserved.paddedArrayRef().size());
+ EXPECT_LE(vReserved.size(), vReserved.paddedArrayRef().size());
+
+ EXPECT_LE(v.paddedSize(), vReserved.paddedSize());
+}
+
+TYPED_TEST(PaddedVectorTest, CanCopyAssign)
+{
+ using VectorType = PaddedVector<typename TypeParam::value_type, TypeParam>;
+
+ VectorType v, w;
+ fillInput(&v, 1);
+ fillInput(&w, 2);
+
+ w = v;
+ compareViews(v.paddedArrayRef(), w.paddedArrayRef());
+ compareViews(makeArrayRef(v), makeArrayRef(v));
+}
+
+TYPED_TEST(PaddedVectorTest, CanMoveAssign)
+{
+ using VectorType = PaddedVector<typename TypeParam::value_type, TypeParam>;
+
+ VectorType v, w, x;
+ fillInput(&v, 1);
+ fillInput(&w, 2);
+ fillInput(&x, 1);
+
+ SCOPED_TRACE("Comparing padded views before move");
+ compareViews(v.paddedArrayRef(), x.paddedArrayRef());
+ SCOPED_TRACE("Comparing unpadded views before move");
+ compareViews(makeArrayRef(v), makeArrayRef(x));
+
+ w = std::move(x);
+
+ SCOPED_TRACE("Comparing padded views");
+ compareViews(v.paddedArrayRef(), w.paddedArrayRef());
+ SCOPED_TRACE("Comparing unpadded views");
+ compareViews(makeArrayRef(v), makeArrayRef(w));
+}
+
+TYPED_TEST(PaddedVectorTest, CanSwap)
+{
+ using VectorType = PaddedVector<typename TypeParam::value_type, TypeParam>;
+
+ VectorType v, w, x;
+ fillInput(&v, 1);
+ fillInput(&w, 2);
+ fillInput(&x, 1);
+
+ std::swap(w, x);
+ compareViews(v.paddedArrayRef(), w.paddedArrayRef());
+ compareViews(makeArrayRef(v), makeArrayRef(w));
+}
+
+} // namespace test
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \libinternal \file
+ * \brief Declares functions for comparing views of vector-like data.
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_math
+ * \inlibraryapi
+ */
+#ifndef GMX_MATH_TESTARRAYREFS_H
+#define GMX_MATH_TESTARRAYREFS_H
+
+#include <gtest/gtest.h>
+
+#include "gromacs/math/vectypes.h"
+
+namespace gmx
+{
+namespace test
+{
+
+//! Initialization overload for non-BasicVector
+template <typename T>
+void fillInputContents(ArrayRef<T> inputRef,
+ int scaleFactor)
+{
+ inputRef[0] = 1;
+ inputRef[1] = 2;
+ inputRef[2] = 3;
+ for (auto &element : inputRef)
+ {
+ element *= scaleFactor;
+ }
+}
+
+//! Initialization overload for BasicVector
+template <typename T>
+void fillInputContents(ArrayRef < BasicVector < T>> inputRef,
+ int scaleFactor)
+{
+ inputRef[0] = {1, 2, 3};
+ inputRef[1] = {4, 5, 6};
+ inputRef[2] = {7, 8, 9};
+ for (auto &element : inputRef)
+ {
+ element *= scaleFactor;
+ }
+}
+
+//! Dispatcher function for filling.
+template <typename PaddedVectorOfT>
+void fillInput(PaddedVectorOfT *input,
+ int scaleFactor)
+{
+ // Use a size for the vector in tests that is prime enough to
+ // expose problems where they exist.
+ int sizeOfVector = 3;
+ input->resizeWithPadding(sizeOfVector);
+ fillInputContents(makeArrayRef(*input), scaleFactor);
+ EXPECT_LE(sizeOfVector, input->paddedSize());
+}
+
+//! Comparison overload for non-BasicVector
+template <typename T>
+void compareViews(ArrayRef<T> input,
+ ArrayRef<T> output)
+{
+ ASSERT_EQ(input.size(), output.size());
+ for (index i = 0; i != input.size(); ++i)
+ {
+ EXPECT_EQ(input[i], output[i]) << "for index " << i;
+ }
+}
+
+//! Comparison overload for BasicVector<T>
+template <typename T>
+void compareViews(ArrayRef < BasicVector < T>> input,
+ ArrayRef < BasicVector < T>> output)
+{
+ ASSERT_EQ(input.size(), output.size());
+ for (index i = 0; i != input.size(); ++i)
+ {
+ EXPECT_EQ(input[i][XX], output[i][XX]) << "for index " << i;
+ EXPECT_EQ(input[i][YY], output[i][YY]) << "for index " << i;
+ EXPECT_EQ(input[i][ZZ], output[i][ZZ]) << "for index " << i;
+ }
+}
+
+} // namespace test
+} // namespace gmx
+
+#endif
}
template <typename AllocatorType>
-static void bcastPaddedRVecVector(const t_commrec *cr, std::vector<gmx::RVec, AllocatorType> *v, int numAtoms)
+static void bcastPaddedRVecVector(const t_commrec *cr, gmx::PaddedVector<gmx::RVec, AllocatorType> *v, int numAtoms)
{
- v->resize(gmx::paddedRVecVectorSize(numAtoms));
- nblock_bc(cr, numAtoms, as_rvec_array(v->data()));
+ v->resizeWithPadding(numAtoms);
+ nblock_bc(cr, makeArrayRef(*v));
}
void broadcastStateWithoutDynamics(const t_commrec *cr, t_state *state)
#include "gromacs/gmxlib/network.h"
#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/smalloc.h"
//! Convenience wrapper for gmx_bcast of a single value.
{
gmx_bcast(numElements * sizeof(T), static_cast<void *>(data), cr);
}
+//! Convenience wrapper for gmx_bcast of an ArrayRef<T>
+template <typename T>
+void nblock_bc(const t_commrec *cr, gmx::ArrayRef<T> data)
+{
+ gmx_bcast(data.size() * sizeof(T), static_cast<void *>(data.data()), cr);
+}
//! Convenience wrapper for allocation with snew of vectors that need allocation on non-master ranks.
template <typename T>
void snew_bc(const t_commrec *cr, T * &data, int numElements)
}
void andersen_tcoupl(const t_inputrec *ir, int64_t step,
- const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
+ const t_commrec *cr, const t_mdatoms *md,
+ gmx::ArrayRef<gmx::RVec> v,
+ real rate, const gmx_bool *randomize, const real *boltzfac)
{
const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
int i;
for (d = 0; d < DIM; d++)
{
- state->v[i][d] = scal*normalDist(rng);
+ v[i][d] = scal*normalDist(rng);
}
}
}
dt = dtc;
}
+ auto v = makeArrayRef(state->v);
switch (trotter_seq[i])
{
case etrtBAROV:
}
for (d = 0; d < DIM; d++)
{
- state->v[n][d] *= scalefac[gc];
+ v[n][d] *= scalefac[gc];
}
if (debug)
{
for (d = 0; d < DIM; d++)
{
- sumv[d] += (state->v[n][d])/md->invmass[n];
+ sumv[d] += (v[n][d])/md->invmass[n];
}
}
}
if (bStopCM)
{
calc_vcm_grp(0, mdatoms->homenr, mdatoms,
- as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), vcm);
+ state->x.rvec_array(), state->v.rvec_array(), vcm);
}
if (bTemp || bStopCM || bPres || bEner || bConstrain || bCheckNumberOfBondedInteractions)
rvec *xPtr = nullptr;
if (vcm->mode == ecmANGULAR || (vcm->mode == ecmLINEAR_ACCELERATION_CORRECTION && !(flags & CGLO_INITIALIZATION)))
{
- xPtr = as_rvec_array(state->x.data());
+ xPtr = state->x.rvec_array();
}
do_stopcm_grp(*mdatoms,
- xPtr, as_rvec_array(state->v.data()), *vcm);
+ xPtr, state->v.rvec_array(), *vcm);
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
state->flags |= (1<<estFEPSTATE);
}
state->flags |= (1<<estX);
- GMX_RELEASE_ASSERT(state->x.size() >= static_cast<unsigned int>(state->natoms), "We should start a run with an initialized state->x");
+ GMX_RELEASE_ASSERT(state->x.size() == static_cast<unsigned int>(state->natoms), "We should start a run with an initialized state->x");
if (EI_DYNAMICS(ir->eI))
{
state->flags |= (1<<estV);
void MDAtoms::resize(int newSize)
{
- chargeA_.resize(newSize);
+ chargeA_.resizeWithPadding(newSize);
mdatoms_->chargeA = chargeA_.data();
}
void MDAtoms::reserve(int newCapacity)
{
- chargeA_.reserve(newCapacity);
+ chargeA_.reserveWithPadding(newCapacity);
mdatoms_->chargeA = chargeA_.data();
}
{
if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
{
- gmx::ArrayRef<gmx::RVec> globalXRef = MASTER(cr) ? gmx::makeArrayRef(state_global->x) : gmx::EmptyArrayRef();
- dd_collect_vec(cr->dd, state_local, state_local->x, globalXRef);
+ gmx::ArrayRef<gmx::RVec> globalXRef = MASTER(cr) ? makeArrayRef(state_global->x) : gmx::EmptyArrayRef();
+ dd_collect_vec(cr->dd, state_local, makeArrayRef(state_local->x), globalXRef);
}
if (mdof_flags & MDOF_V)
{
- gmx::ArrayRef<gmx::RVec> globalVRef = MASTER(cr) ? gmx::makeArrayRef(state_global->v) : gmx::EmptyArrayRef();
- dd_collect_vec(cr->dd, state_local, state_local->v, globalVRef);
+ gmx::ArrayRef<gmx::RVec> globalVRef = MASTER(cr) ? makeArrayRef(state_global->v) : gmx::EmptyArrayRef();
+ dd_collect_vec(cr->dd, state_local, makeArrayRef(state_local->v), globalVRef);
}
}
f_global = of->f_global;
if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
{
- const rvec *x = (mdof_flags & MDOF_X) ? as_rvec_array(state_global->x.data()) : nullptr;
- const rvec *v = (mdof_flags & MDOF_V) ? as_rvec_array(state_global->v.data()) : nullptr;
+ const rvec *x = (mdof_flags & MDOF_X) ? state_global->x.rvec_array() : nullptr;
+ const rvec *v = (mdof_flags & MDOF_V) ? state_global->v.rvec_array() : nullptr;
const rvec *f = (mdof_flags & MDOF_F) ? f_global : nullptr;
if (of->fp_trn)
{
/* We are writing the positions of all of the atoms to
the compressed output */
- xxtc = as_rvec_array(state_global->x.data());
+ xxtc = state_global->x.rvec_array();
}
else
{
int i, j;
snew(xxtc, of->natoms_x_compressed);
+ auto x = makeArrayRef(state_global->x);
for (i = 0, j = 0; (i < of->natoms_global); i++)
{
if (getGroupType(of->groups, egcCompressedX, i) == 0)
{
- copy_rvec(state_global->x[i], xxtc[j++]);
+ copy_rvec(x[i], xxtc[j++]);
}
}
}
gmx_groups_t *groups, int ins_grp_id, real xy_max)
{
int i, gid, c = 0;
- real x, xmin, xmax, y, ymin, ymax, z, zmin, zmax;
+ real xmin, xmax, ymin, ymax, zmin, zmax;
const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to *
* determine the overlap between molecule to embed and membrane */
const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
snew(rest_at->index, state->natoms);
+ auto x = makeArrayRef(state->x);
- xmin = xmax = state->x[ins_at->index[0]][XX];
- ymin = ymax = state->x[ins_at->index[0]][YY];
- zmin = zmax = state->x[ins_at->index[0]][ZZ];
+ xmin = xmax = x[ins_at->index[0]][XX];
+ ymin = ymax = x[ins_at->index[0]][YY];
+ zmin = zmax = x[ins_at->index[0]][ZZ];
for (i = 0; i < state->natoms; i++)
{
gid = groups->grpnr[egcFREEZE][i];
if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
{
- x = state->x[i][XX];
- if (x < xmin)
- {
- xmin = x;
- }
- if (x > xmax)
- {
- xmax = x;
- }
- y = state->x[i][YY];
- if (y < ymin)
- {
- ymin = y;
- }
- if (y > ymax)
- {
- ymax = y;
- }
- z = state->x[i][ZZ];
- if (z < zmin)
- {
- zmin = z;
- }
- if (z > zmax)
- {
- zmax = z;
- }
+ xmin = std::min(xmin, x[i][XX]);
+ xmax = std::max(xmax, x[i][XX]);
+ ymin = std::min(ymin, x[i][YY]);
+ ymax = std::max(ymax, x[i][YY]);
+ zmin = std::min(zmin, x[i][ZZ]);
+ zmax = std::max(zmax, x[i][ZZ]);
}
else
{
}
}
+ auto x = makeArrayRef(state->x);
+ auto v = makeArrayRef(state->v);
rm = 0;
for (int i = 0; i < state->natoms+n; i++)
{
new_egrp[j][i-rm] = groups->grpnr[j][i];
}
}
- copy_rvec(state->x[i], x_tmp[i-rm]);
- copy_rvec(state->v[i], v_tmp[i-rm]);
+ copy_rvec(x[i], x_tmp[i-rm]);
+ copy_rvec(v[i], v_tmp[i-rm]);
for (j = 0; j < ins_at->nr; j++)
{
if (i == ins_at->index[j])
}
for (int i = 0; i < state->natoms; i++)
{
- copy_rvec(x_tmp[i], state->x[i]);
+ copy_rvec(x_tmp[i], x[i]);
}
sfree(x_tmp);
for (int i = 0; i < state->natoms; i++)
{
- copy_rvec(v_tmp[i], state->v[i]);
+ copy_rvec(v_tmp[i], v[i]);
}
sfree(v_tmp);
/* Check that moleculetypes in insertion group are not part of the rest of the system */
check_types(ins_at, rest_at, mtop);
- init_mem_at(mem_p, mtop, as_rvec_array(state->x.data()), state->box, pos_ins);
+ init_mem_at(mem_p, mtop, state->x.rvec_array(), state->box, pos_ins);
- prot_area = est_prot_area(pos_ins, as_rvec_array(state->x.data()), ins_at, mem_p);
+ prot_area = est_prot_area(pos_ins, state->x.rvec_array(), ins_at, mem_p);
if ( (prot_area > prot_vs_box) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX]) < box_vs_prot) )
{
warn++;
/* resize the protein by xy and by z if necessary*/
snew(r_ins, ins_at->nr);
- init_resize(ins_at, r_ins, pos_ins, mem_p, as_rvec_array(state->x.data()), bALLOW_ASYMMETRY);
+ init_resize(ins_at, r_ins, pos_ins, mem_p, state->x.rvec_array(), bALLOW_ASYMMETRY);
membed->fac[0] = membed->fac[1] = xy_fac;
membed->fac[2] = z_fac;
membed->xy_step = (xy_max-xy_fac)/static_cast<double>(it_xy);
membed->z_step = (z_max-z_fac)/static_cast<double>(it_z-1);
- resize(r_ins, as_rvec_array(state->x.data()), pos_ins, membed->fac);
+ resize(r_ins, state->x.rvec_array(), pos_ins, membed->fac);
/* remove overlapping lipids and water from the membrane box*/
/*mark molecules to be removed*/
set_pbc(pbc, inputrec->ePBC, state->box);
snew(rm_p, 1);
- lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, as_rvec_array(state->x.data()), mem_p, pos_ins,
+ lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x.rvec_array(), mem_p, pos_ins,
probe_rad, low_up_rm, bALLOW_ASYMMETRY);
lip_rm -= low_up_rm;
int nflexcon; /* The number of flexible constraints */
/* Temporary arrays, should be fixed size 2 when fully converted to C++ */
- PaddedRVecVector *x; /* Array for iterative minimization */
- PaddedRVecVector *f; /* Array for iterative minimization */
+ PaddedVector<gmx::RVec> *x; /* Array for iterative minimization */
+ PaddedVector<gmx::RVec> *f; /* Array for iterative minimization */
/* Flexible constraint working data */
rvec *acc_dir; /* Acceleration direction for flexcon */
}
snew(shfc, 1);
- shfc->x = new PaddedRVecVector[2] {};
- shfc->f = new PaddedRVecVector[2] {};
+ shfc->x = new PaddedVector<gmx::RVec>[2] {};
+ shfc->f = new PaddedVector<gmx::RVec>[2] {};
shfc->nflexcon = nflexcon;
if (nshell == 0)
for (i = 0; (i < 2); i++)
{
- shfc->x[i].resize(gmx::paddedRVecVectorSize(nat));
- shfc->f[i].resize(gmx::paddedRVecVectorSize(nat));
+ shfc->x[i].resizeWithPadding(nat);
+ shfc->f[i].resizeWithPadding(nat);
}
/* Create views that we can swap */
gmx::PaddedArrayRef<gmx::RVec> force[2];
for (i = 0; (i < 2); i++)
{
- pos[i] = shfc->x[i];
- force[i] = shfc->f[i];
+ pos[i] = shfc->x[i].paddedArrayRef();
+ force[i] = shfc->f[i].paddedArrayRef();
}
if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr))
*/
if (inputrec->cutoff_scheme == ecutsVERLET)
{
- auto xRef = makeArrayRef(state->x);
+ auto xRef = state->x.paddedArrayRef();
put_atoms_in_box_omp(fr->ePBC, state->box, xRef.subArray(0, md->homenr));
}
else
cg0 = 0;
cg1 = top->cgs.nr;
put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box,
- &(top->cgs), as_rvec_array(state->x.data()), fr->cg_cm);
+ &(top->cgs), state->x.rvec_array(), fr->cg_cm);
}
if (graph)
{
- mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
+ mk_mshift(fplog, graph, fr->ePBC, state->box, state->x.rvec_array());
}
}
/* After this all coordinate arrays will contain whole charge groups */
if (graph)
{
- shift_self(graph, state->box, as_rvec_array(state->x.data()));
+ shift_self(graph, state->box, state->x.rvec_array());
}
if (nflexcon)
}
acc_dir = shfc->acc_dir;
x_old = shfc->x_old;
+ auto x = makeArrayRef(state->x);
+ auto v = makeArrayRef(state->v);
for (i = 0; i < homenr; i++)
{
for (d = 0; d < DIM; d++)
{
shfc->x_old[i][d] =
- state->x[i][d] - state->v[i][d]*inputrec->delta_t;
+ x[i][d] - v[i][d]*inputrec->delta_t;
}
}
}
*/
if (shfc->bPredict && !bCont && (EI_STATE_VELOCITY(inputrec->eI) || bInit))
{
- predict_shells(fplog, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), inputrec->delta_t, nshell, shell,
+ predict_shells(fplog, state->x.rvec_array(), state->v.rvec_array(), inputrec->delta_t, nshell, shell,
md->massT, nullptr, bInit);
}
/* do_force expected the charge groups to be in the box */
if (graph)
{
- unshift_self(graph, state->box, as_rvec_array(state->x.data()));
+ unshift_self(graph, state->box, state->x.rvec_array());
}
/* Calculate the forces first time around */
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "x b4 do_force", as_rvec_array(state->x.data()), homenr);
+ pr_rvecs(debug, 0, "x b4 do_force", state->x.rvec_array(), homenr);
}
do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation,
mdstep, nrnb, wcycle, top, groups,
- state->box, state->x, &state->hist,
+ state->box, state->x.paddedArrayRef(), &state->hist,
force[Min], force_vir, md, enerd, fcd,
state->lambda, graph,
fr, vsite, mu_tot, t, nullptr,
{
init_adir(shfc,
constr, inputrec, cr, dd_ac1, mdstep, md, end,
- shfc->x_old, as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), as_rvec_array(force[Min].data()),
+ shfc->x_old, state->x.rvec_array(), state->x.rvec_array(), as_rvec_array(force[Min].data()),
shfc->acc_dir,
state->box, state->lambda, &dum);
Epot[Min] = enerd->term[F_EPOT];
- df[Min] = rms_force(cr, shfc->f[Min], nshell, shell, nflexcon, &sf_dir, &Epot[Min]);
+ df[Min] = rms_force(cr, shfc->f[Min].paddedArrayRef(), nshell, shell, nflexcon, &sf_dir, &Epot[Min]);
df[Try] = 0;
if (debug)
{
* shell positions are updated, therefore the other particles must
* be set here.
*/
- pos[Min] = state->x;
- pos[Try] = state->x;
+ pos[Min] = state->x.paddedArrayRef();
+ pos[Try] = state->x.paddedArrayRef();
}
if (bVerbose && MASTER(cr))
if (vsite)
{
construct_vsites(vsite, as_rvec_array(pos[Min].data()),
- inputrec->delta_t, as_rvec_array(state->v.data()),
+ inputrec->delta_t, state->v.rvec_array(),
idef->iparams, idef->il,
fr->ePBC, fr->bMolPBC, cr, state->box);
}
{
init_adir(shfc,
constr, inputrec, cr, dd_ac1, mdstep, md, end,
- x_old, as_rvec_array(state->x.data()), as_rvec_array(pos[Min].data()), as_rvec_array(force[Min].data()), acc_dir,
+ x_old, state->x.rvec_array(), as_rvec_array(pos[Min].data()), as_rvec_array(force[Min].data()), acc_dir,
state->box, state->lambda, &dum);
directional_sd(pos[Min], pos[Try], acc_dir, end, fr->fc_stepsize);
{
init_adir(shfc,
constr, inputrec, cr, dd_ac1, mdstep, md, end,
- x_old, as_rvec_array(state->x.data()), as_rvec_array(pos[Try].data()), as_rvec_array(force[Try].data()), acc_dir,
+ x_old, state->x.rvec_array(), as_rvec_array(pos[Try].data()), as_rvec_array(force[Try].data()), acc_dir,
state->box, state->lambda, &dum);
for (i = 0; i < end; i++)
{
/* Correct the velocities for the flexible constraints */
invdt = 1/inputrec->delta_t;
+ auto v = makeArrayRef(state->v);
for (i = 0; i < end; i++)
{
for (d = 0; d < DIM; d++)
{
- state->v[i][d] += (pos[Try][i][d] - pos[Min][i][d])*invdt;
+ v[i][d] += (pos[Try][i][d] - pos[Min][i][d])*invdt;
}
}
}
}
/* Copy back the coordinates and the forces */
- std::copy(pos[Min].begin(), pos[Min].end(), state->x.begin());
+ std::copy(pos[Min].begin(), pos[Min].end(), makeArrayRef(state->x).data());
std::copy(force[Min].begin(), force[Min].end(), f.begin());
}
flags &= ~GMX_FORCE_NONBONDED;
}
- GMX_ASSERT(x.size() >= gmx::paddedRVecVectorSize(fr->natoms_force), "coordinates should be padded");
- GMX_ASSERT(force.size() >= gmx::paddedRVecVectorSize(fr->natoms_force), "force should be padded");
-
switch (inputrec->cutoff_scheme)
{
case ecutsVERLET:
/* constrain the current position */
constr->apply(TRUE, FALSE,
step, 0, 1.0,
- as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), nullptr,
+ state->x.rvec_array(), state->x.rvec_array(), nullptr,
state->box,
state->lambda[efptBONDED], &dvdl_dum,
nullptr, nullptr, gmx::ConstraintVariable::Positions);
/* also may be useful if we need the ekin from the halfstep for velocity verlet */
constr->apply(TRUE, FALSE,
step, 0, 1.0,
- as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
+ state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
state->box,
state->lambda[efptBONDED], &dvdl_dum,
nullptr, nullptr, gmx::ConstraintVariable::Velocities);
/* constrain the inital velocities at t-dt/2 */
if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
{
+ auto x = makeArrayRef(state->x).subArray(start, end);
+ auto v = makeArrayRef(state->v).subArray(start, end);
for (i = start; (i < end); i++)
{
for (m = 0; (m < DIM); m++)
{
/* Reverse the velocity */
- state->v[i][m] = -state->v[i][m];
+ v[i][m] = -v[i][m];
/* Store the position at t-dt in buf */
- savex[i][m] = state->x[i][m] + dt*state->v[i][m];
+ savex[i][m] = x[i][m] + dt*v[i][m];
}
}
/* Shake the positions at t=-dt with the positions at t=0
dvdl_dum = 0;
constr->apply(TRUE, FALSE,
step, -1, 1.0,
- as_rvec_array(state->x.data()), savex, nullptr,
+ state->x.rvec_array(), savex, nullptr,
state->box,
state->lambda[efptBONDED], &dvdl_dum,
- as_rvec_array(state->v.data()), nullptr, gmx::ConstraintVariable::Positions);
+ state->v.rvec_array(), nullptr, gmx::ConstraintVariable::Positions);
for (i = start; i < end; i++)
{
for (m = 0; m < DIM; m++)
{
/* Re-reverse the velocities */
- state->v[i][m] = -state->v[i][m];
+ v[i][m] = -v[i][m];
}
}
}
identical, and makes .edr restarts binary
identical. */
snew(x_for_confout, state_global->natoms);
- copy_rvecn(as_rvec_array(state_global->x.data()), x_for_confout, 0, state_global->natoms);
+ copy_rvecn(state_global->x.rvec_array(), x_for_confout, 0, state_global->natoms);
}
else
{
/* With DD, or no bMolPBC, it doesn't matter if
- we change as_rvec_array(state_global->x.data()) */
- x_for_confout = as_rvec_array(state_global->x.data());
+ we change state_global->x.rvec_array() */
+ x_for_confout = state_global->x.rvec_array();
}
/* x and v have been collected in mdoutf_write_to_trajectory_files,
}
write_sto_conf_mtop(ftp2fn(efSTO, nfile, fnm),
*top_global->name, top_global,
- x_for_confout, as_rvec_array(state_global->v.data()),
+ x_for_confout, state_global->v.rvec_array(),
ir->ePBC, state->box);
if (fr->bMolPBC && state == state_global)
{
struct gmx_update_t
{
- gmx_stochd_t *sd;
+ gmx_stochd_t *sd;
/* xprime for constraint algorithms */
- PaddedRVecVector xp;
+ PaddedVector<gmx::RVec> xp;
/* Variables for the deform algorithm */
int64_t deformref_step;
update_temperature_constants(upd, ir);
- upd->xp.resize(0);
+ upd->xp.resizeWithPadding(0);
upd->deform = deform;
{
GMX_ASSERT(upd, "upd must be allocated before its fields can be reallocated");
- upd->xp.resize(gmx::paddedRVecVectorSize(natoms));
+ upd->xp.resizeWithPadding(natoms);
}
/*! \brief Sets the SD update type */
{
if (ekind->cosacc.cos_accel == 0)
{
- calc_ke_part_normal(as_rvec_array(state->v.data()), opts, md, ekind, nrnb, bEkinAveVel);
+ calc_ke_part_normal(state->v.rvec_array(), opts, md, ekind, nrnb, bEkinAveVel);
}
else
{
- calc_ke_part_visc(state->box, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), opts, md, ekind, nrnb, bEkinAveVel);
+ calc_ke_part_visc(state->box, state->x.rvec_array(), state->v.rvec_array(), opts, md, ekind, nrnb, bEkinAveVel);
}
}
/* rescale in place here */
if (EI_VV(inputrec->eI))
{
- rescale_velocities(ekind, md, 0, md->homenr, as_rvec_array(state->v.data()));
+ rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
}
}
else
/* Constrain the coordinates upd->xp */
constr->apply(do_log, do_ene,
step, 1, 1.0,
- as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
+ state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
state->box,
state->lambda[efptBONDED], dvdlambda,
nullptr, bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities);
/* Constrain the coordinates upd->xp */
constr->apply(do_log, do_ene,
step, 1, 1.0,
- as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
+ state->x.rvec_array(), upd->xp.rvec_array(), nullptr,
state->box,
state->lambda[efptBONDED], dvdlambda,
as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr, ConstraintVariable::Positions);
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, nullptr, md->cTC,
- as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()),
- as_rvec_array(state->v.data()), nullptr,
+ state->x.rvec_array(), upd->xp.rvec_array(),
+ state->v.rvec_array(), nullptr,
step, inputrec->ld_seed,
DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
}
/* Constrain the coordinates upd->xp for half a time step */
constr->apply(do_log, do_ene,
step, 1, 0.5,
- as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
+ state->x.rvec_array(), upd->xp.rvec_array(), nullptr,
state->box,
state->lambda[efptBONDED], dvdlambda,
as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions);
}
if (partialFreezeAndConstraints)
{
+ auto xp = makeArrayRef(upd->xp).subArray(0, homenr);
+ auto x = makeConstArrayRef(state->x).subArray(0, homenr);
for (int i = 0; i < homenr; i++)
{
int g = md->cFREEZE[i];
{
if (nFreeze[g][d])
{
- upd->xp[i][d] = state->x[i][d];
+ xp[i][d] = x[i][d];
}
}
}
if (graph && (graph->nnodes > 0))
{
- unshift_x(graph, state->box, as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()));
+ unshift_x(graph, state->box, state->x.rvec_array(), upd->xp.rvec_array());
if (TRICLINIC(state->box))
{
inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
}
else
{
- /* The copy is performance sensitive, so use a bare pointer */
- rvec *xp = as_rvec_array(upd->xp.data());
+ auto xp = makeConstArrayRef(upd->xp).subArray(0, homenr);
+ auto x = makeArrayRef(state->x).subArray(0, homenr);
#ifndef __clang_analyzer__
int gmx_unused nth = gmx_omp_nthreads_get(emntUpdate);
#endif
for (int i = 0; i < homenr; i++)
{
// Trivial statement, does not throw
- copy_rvec(xp[i], state->x[i]);
+ x[i] = xp[i];
}
}
wallcycle_stop(wcycle, ewcUPDATE);
forceVirial, constraintVirial,
mu, &state->baros_integral);
berendsen_pscale(inputrec, mu, state->box, state->box_rel,
- start, homenr, as_rvec_array(state->x.data()),
+ start, homenr, state->x.rvec_array(),
md->cFREEZE, nrnb);
}
break;
preserve_box_shape(inputrec, state->box_rel, state->box);
/* Scale the coordinates */
+ auto x = state->x.rvec_array();
for (int n = start; n < start + homenr; n++)
{
- tmvmul_ur0(parrinellorahmanMu, state->x[n], state->x[n]);
+ tmvmul_ur0(parrinellorahmanMu, x[n], x[n]);
}
}
break;
if (upd->deform)
{
- auto localX = gmx::ArrayRef<gmx::RVec>(state->x).subArray(start, homenr);
+ auto localX = makeArrayRef(state->x).subArray(start, homenr);
upd->deform->apply(localX, state->box, step);
}
}
int start_th, end_th;
getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
- /* Strictly speaking, we would only need this check with SIMD
- * and for the actual SIMD width. But since the code currently
- * always adds padding for GMX_REAL_MAX_SIMD_WIDTH, we check that.
- */
- gmx::index gmx_used_in_debug homenrSimdPadded;
- homenrSimdPadded = ((homenr + GMX_REAL_MAX_SIMD_WIDTH - 1)/GMX_REAL_MAX_SIMD_WIDTH)*GMX_REAL_MAX_SIMD_WIDTH;
- GMX_ASSERT(gmx::index(state->x.size()) >= homenrSimdPadded, "state->x needs to be padded for SIMD access");
- GMX_ASSERT(gmx::index(upd->xp.size()) >= homenrSimdPadded, "upd->xp needs to be padded for SIMD access");
- GMX_ASSERT(gmx::index(state->v.size()) >= homenrSimdPadded, "state->v needs to be padded for SIMD access");
- GMX_ASSERT(f.size() >= homenrSimdPadded, "f needs to be padded for SIMD access");
-
- const rvec *x_rvec = as_rvec_array(state->x.data());
- rvec *xp_rvec = as_rvec_array(upd->xp.data());
- rvec *v_rvec = as_rvec_array(state->v.data());
+ const rvec *x_rvec = state->x.rvec_array();
+ rvec *xp_rvec = upd->xp.rvec_array();
+ rvec *v_rvec = state->v.rvec_array();
const rvec *f_rvec = as_rvec_array(f.data());
switch (inputrec->eI)
}
extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step, const t_commrec *cr,
- const t_mdatoms *md, t_state *state, const gmx_update_t *upd,
+ const t_mdatoms *md,
+ gmx::ArrayRef<gmx::RVec> v,
+ const gmx_update_t *upd,
const gmx::Constraints *constr)
{
particle andersen or 2) it's massive andersen and it's tau_t/dt */
if ((ir->etc == etcANDERSEN) || do_per_step(step, roundToInt(1.0/rate)))
{
- andersen_tcoupl(ir, step, cr, md, state, rate,
+ andersen_tcoupl(ir, step, cr, md, v, rate,
upd->sd->randomize_group, upd->sd->boltzfac);
return TRUE;
}
/* Return TRUE if OK, FALSE in case of Shake Error */
extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step, const t_commrec *cr,
- const t_mdatoms *md, t_state *state, const gmx_update_t *upd,
+ const t_mdatoms *md,
+ gmx::ArrayRef<gmx::RVec> v,
+ const gmx_update_t *upd,
const gmx::Constraints *constr);
void constrain_velocities(int64_t step,
std::vector<double> &therm_integral); //NOLINT(google-runtime-references)
void andersen_tcoupl(const t_inputrec *ir, int64_t step,
- const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac);
+ const t_commrec *cr, const t_mdatoms *md,
+ gmx::ArrayRef<gmx::RVec> v,
+ real rate, const gmx_bool *randomize, const real *boltzfac);
void nosehoover_tcoupl(const t_grpopts *opts, const gmx_ekindata_t *ekind, real dt,
double xi[], double vxi[], const t_extmass *MassQ);
// alias to avoid a large ripple of nearly useless changes.
// t_inputrec is being replaced by IMdpOptionsProvider, so this
// will go away eventually.
- t_inputrec *ir = inputrec;
- gmx_mdoutf *outf = nullptr;
- int64_t step, step_rel;
- double t, t0, lam0[efptNR];
- gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
- gmx_bool bNS, bNStList, bSimAnn, bStopCM,
- bFirstStep, bInitStep, bLastStep = FALSE;
- gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
- gmx_bool do_ene, do_log, do_verbose;
- gmx_bool bMasterState;
- int force_flags, cglo_flags;
- tensor force_vir, shake_vir, total_vir, tmp_vir, pres;
- int i, m;
- rvec mu_tot;
- t_vcm *vcm;
- matrix parrinellorahmanMu, M;
- gmx_repl_ex_t repl_ex = nullptr;
- gmx_localtop_t *top;
- t_mdebin *mdebin = nullptr;
- gmx_enerdata_t *enerd;
- PaddedRVecVector f {};
- gmx_global_stat_t gstat;
- gmx_update_t *upd = nullptr;
- t_graph *graph = nullptr;
- gmx_groups_t *groups;
- gmx_ekindata_t *ekind;
- gmx_shellfc_t *shellfc;
- gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
- gmx_bool bTemp, bPres, bTrotter;
- real dvdl_constr;
- rvec *cbuf = nullptr;
- int cbuf_nalloc = 0;
- matrix lastbox;
- int lamnew = 0;
+ t_inputrec *ir = inputrec;
+ gmx_mdoutf *outf = nullptr;
+ int64_t step, step_rel;
+ double t, t0, lam0[efptNR];
+ gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
+ gmx_bool bNS, bNStList, bSimAnn, bStopCM,
+ bFirstStep, bInitStep, bLastStep = FALSE;
+ gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
+ gmx_bool do_ene, do_log, do_verbose;
+ gmx_bool bMasterState;
+ int force_flags, cglo_flags;
+ tensor force_vir, shake_vir, total_vir, tmp_vir, pres;
+ int i, m;
+ rvec mu_tot;
+ t_vcm *vcm;
+ matrix parrinellorahmanMu, M;
+ gmx_repl_ex_t repl_ex = nullptr;
+ gmx_localtop_t *top;
+ t_mdebin *mdebin = nullptr;
+ gmx_enerdata_t *enerd;
+ PaddedVector<gmx::RVec> f {};
+ gmx_global_stat_t gstat;
+ gmx_update_t *upd = nullptr;
+ t_graph *graph = nullptr;
+ gmx_groups_t *groups;
+ gmx_ekindata_t *ekind;
+ gmx_shellfc_t *shellfc;
+ gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
+ gmx_bool bTemp, bPres, bTrotter;
+ real dvdl_constr;
+ rvec *cbuf = nullptr;
+ int cbuf_nalloc = 0;
+ matrix lastbox;
+ int lamnew = 0;
/* for FEP */
- int nstfep = 0;
- double cycles;
- real saved_conserved_quantity = 0;
- real last_ekin = 0;
- t_extmass MassQ;
- int **trotter_seq;
- char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
+ int nstfep = 0;
+ double cycles;
+ real saved_conserved_quantity = 0;
+ real last_ekin = 0;
+ t_extmass MassQ;
+ int **trotter_seq;
+ char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
/* PME load balancing data for GPU kernels */
pme_load_balancing_t *pme_loadbal = nullptr;
/* Set up interactive MD (IMD) */
init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
- MASTER(cr) ? as_rvec_array(state_global->x.data()) : nullptr,
+ MASTER(cr) ? state_global->x.rvec_array() : nullptr,
nfile, fnm, oenv, mdrunOptions);
// Local state only becomes valid now.
else
{
state_change_natoms(state_global, state_global->natoms);
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- f.resize(gmx::paddedRVecVectorSize(state_global->natoms));
+ f.resizeWithPadding(state_global->natoms);
/* Copy the pointer to the global state */
state = state_global;
{
if (state->flags & (1 << estV))
{
+ auto v = makeArrayRef(state->v);
/* Set the velocities of vsites, shells and frozen atoms to zero */
for (i = 0; i < mdatoms->homenr; i++)
{
if (mdatoms->ptype[i] == eptVSite ||
mdatoms->ptype[i] == eptShell)
{
- clear_rvec(state->v[i]);
+ clear_rvec(v[i]);
}
else if (mdatoms->cFREEZE)
{
{
if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
{
- state->v[i][m] = 0;
+ v[i][m] = 0;
}
}
}
if (vsite)
{
/* Construct the virtual sites for the initial configuration */
- construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
+ construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
}
enforcedRotation, step,
ir, bNS, force_flags, top,
constr, enerd, fcd,
- state, f, force_vir, mdatoms,
+ state, f.paddedArrayRef(), force_vir, mdatoms,
nrnb, wcycle, graph, groups,
shellfc, fr, t, mu_tot,
vsite,
*/
do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
step, nrnb, wcycle, top, groups,
- state->box, state->x, &state->hist,
- f, force_vir, mdatoms, enerd, fcd,
+ state->box, state->x.paddedArrayRef(), &state->hist,
+ f.paddedArrayRef(), force_vir, mdatoms, enerd, fcd,
state->lambda, graph,
fr, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
(bNS ? GMX_FORCE_NS : 0) | force_flags,
* so that the input is actually the initial step.
*/
snew(vbuf, state->natoms);
- copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
+ copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
}
else
{
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
}
- update_coords(step, ir, mdatoms, state, f, fcd,
+ update_coords(step, ir, mdatoms, state, f.paddedArrayRef(), fcd,
ekind, M, upd, etrtVELOCITY1,
cr, constr);
/* if it's the initial step, we performed this first step just to get the constraint virial */
if (ir->eI == eiVV && bInitStep)
{
- copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
+ copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
sfree(vbuf);
}
wallcycle_stop(wcycle, ewcUPDATE);
statistics, but if performing simulated tempering, we
do update the velocities and the tau_t. */
- lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
+ lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
/* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
if (MASTER(cr))
{
mdrunOptions.writeConfout,
bSumEkinhOld);
/* Check if IMD step and do IMD communication, if bIMD is TRUE. */
- bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
+ bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x.rvec_array(), ir, t, wcycle);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
{
gmx_bool bIfRandomize;
- bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
+ bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, upd, constr);
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
/* now we know the scaling, we can compute the positions again again */
copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
- update_coords(step, ir, mdatoms, state, f, fcd,
+ update_coords(step, ir, mdatoms, state, f.paddedArrayRef(), fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
wallcycle_start(wcycle, ewcVSITECONSTR);
if (graph != nullptr)
{
- shift_self(graph, state->box, as_rvec_array(state->x.data()));
+ shift_self(graph, state->box, state->x.rvec_array());
}
- construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
+ construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
if (graph != nullptr)
{
- unshift_self(graph, state->box, as_rvec_array(state->x.data()));
+ unshift_self(graph, state->box, state->x.rvec_array());
}
wallcycle_stop(wcycle, ewcVSITECONSTR);
}
//! Utility structure for manipulating states during EM
typedef struct {
//! Copy of the global state
- t_state s;
+ t_state s;
//! Force array
- PaddedRVecVector f;
+ PaddedVector<gmx::RVec> f;
//! Potential energy
- real epot;
+ real epot;
//! Norm of the force
- real fnorm;
+ real fnorm;
//! Maximum force
- real fmax;
+ real fmax;
//! Direction
- int a_fmax;
+ int a_fmax;
} em_state_t;
//! Print the EM starting conditions
t_grpopts *opts, t_mdatoms *mdatoms,
em_state_t *ems)
{
- get_f_norm_max(cr, opts, mdatoms, as_rvec_array(ems->f.data()),
+ get_f_norm_max(cr, opts, mdatoms, ems->f.rvec_array(),
&ems->fnorm, &ems->fmax, &ems->a_fmax);
}
/* Interactive molecular dynamics */
init_IMD(ir, cr, ms, top_global, fplog, 1,
- MASTER(cr) ? as_rvec_array(state_global->x.data()) : nullptr,
+ MASTER(cr) ? state_global->x.rvec_array() : nullptr,
nfile, fnm, nullptr, mdrunOptions);
if (ir->eI == eiNM)
/* Just copy the state */
ems->s = *state_global;
state_change_natoms(&ems->s, ems->s.natoms);
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- ems->f.resize(gmx::paddedRVecVectorSize(ems->s.natoms));
+ ems->f.resizeWithPadding(ems->s.natoms);
snew(*top, 1);
mdAlgorithmsSetupAtomData(cr, ir, top_global, *top, fr,
dvdl_constr = 0;
constr->apply(TRUE, TRUE,
-1, 0, 1.0,
- as_rvec_array(ems->s.x.data()),
- as_rvec_array(ems->s.x.data()),
+ ems->s.x.rvec_array(),
+ ems->s.x.rvec_array(),
nullptr,
ems->s.box,
ems->s.lambda[efptFEP], &dvdl_constr,
/* If bX=true, x was collected to state_global in the call above */
if (!bX)
{
- gmx::ArrayRef<gmx::RVec> globalXRef = MASTER(cr) ? gmx::makeArrayRef(state_global->x) : gmx::EmptyArrayRef();
- dd_collect_vec(cr->dd, &state->s, state->s.x, globalXRef);
+ gmx::ArrayRef<gmx::RVec> globalXRef = MASTER(cr) ? makeArrayRef(state_global->x) : gmx::EmptyArrayRef();
+ dd_collect_vec(cr->dd, &state->s, makeArrayRef(state->s.x), globalXRef);
}
}
else
{
/* Make molecules whole only for confout writing */
do_pbc_mtop(fplog, ir->ePBC, state->s.box, top_global,
- as_rvec_array(state_global->x.data()));
+ state_global->x.rvec_array());
}
write_sto_conf_mtop(confout,
*top_global->name, top_global,
- as_rvec_array(state_global->x.data()), nullptr, ir->ePBC, state->s.box);
+ state_global->x.rvec_array(), nullptr, ir->ePBC, state->s.box);
}
}
}
// \returns true when the step succeeded, false when a constraint error occurred
static bool do_em_step(const t_commrec *cr,
t_inputrec *ir, t_mdatoms *md,
- em_state_t *ems1, real a, const PaddedRVecVector *force,
+ em_state_t *ems1, real a, const PaddedVector<gmx::RVec> *force,
em_state_t *ems2,
gmx::Constraints *constr,
int64_t count)
if (s2->natoms != s1->natoms)
{
state_change_natoms(s2, s1->natoms);
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- ems2->f.resize(gmx::paddedRVecVectorSize(s2->natoms));
+ ems2->f.resizeWithPadding(s2->natoms);
}
if (DOMAINDECOMP(cr) && s2->cg_gl.size() != s1->cg_gl.size())
{
nthreads = gmx_omp_nthreads_get(emntUpdate);
#pragma omp parallel num_threads(nthreads)
{
- const rvec *x1 = as_rvec_array(s1->x.data());
- rvec *x2 = as_rvec_array(s2->x.data());
- const rvec *f = as_rvec_array(force->data());
+ const rvec *x1 = s1->x.rvec_array();
+ rvec *x2 = s2->x.rvec_array();
+ const rvec *f = force->rvec_array();
int gf = 0;
#pragma omp for schedule(static) nowait
if (s2->flags & (1<<estCGP))
{
/* Copy the CG p vector */
- const rvec *p1 = as_rvec_array(s1->cg_p.data());
- rvec *p2 = as_rvec_array(s2->cg_p.data());
+ const rvec *p1 = s1->cg_p.rvec_array();
+ rvec *p2 = s2->cg_p.rvec_array();
#pragma omp for schedule(static) nowait
for (int i = start; i < end; i++)
{
validStep =
constr->apply(TRUE, TRUE,
count, 0, 1.0,
- as_rvec_array(s1->x.data()), as_rvec_array(s2->x.data()),
+ s1->x.rvec_array(), s2->x.rvec_array(),
nullptr, s2->box,
s2->lambda[efptBONDED], &dvdl_constr,
nullptr, nullptr, gmx::ConstraintVariable::Positions);
if (vsite)
{
- construct_vsites(vsite, as_rvec_array(ems->s.x.data()), 1, nullptr,
+ construct_vsites(vsite, ems->s.x.rvec_array(), 1, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, ems->s.box);
}
*/
do_force(fplog, cr, ms, inputrec, nullptr, nullptr,
count, nrnb, wcycle, top, &top_global->groups,
- ems->s.box, ems->s.x, &ems->s.hist,
- ems->f, force_vir, mdAtoms->mdatoms(), enerd, fcd,
+ ems->s.box, ems->s.x.paddedArrayRef(), &ems->s.hist,
+ ems->f.paddedArrayRef(), force_vir, mdAtoms->mdatoms(), enerd, fcd,
ems->s.lambda, graph, fr, vsite, mu_tot, t, nullptr,
GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES |
GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY |
{
/* Project out the constraint components of the force */
dvdl_constr = 0;
- rvec *f_rvec = as_rvec_array(ems->f.data());
+ rvec *f_rvec = ems->f.rvec_array();
constr->apply(FALSE, FALSE,
count, 0, 1.0,
- as_rvec_array(ems->s.x.data()), f_rvec, f_rvec,
+ ems->s.x.rvec_array(), f_rvec, f_rvec,
ems->s.box,
ems->s.lambda[efptBONDED], &dvdl_constr,
nullptr, &shake_vir, gmx::ConstraintVariable::ForceDispl);
fprintf(debug, "Doing reorder_partsum\n");
}
- const rvec *fm = as_rvec_array(s_min->f.data());
- const rvec *fb = as_rvec_array(s_b->f.data());
+ const rvec *fm = s_min->f.rvec_array();
+ const rvec *fb = s_b->f.rvec_array();
cgs_gl = dd_charge_groups_global(cr->dd);
index = cgs_gl->index;
(s_min->s.ddp_count == cr->dd->ddp_count &&
s_b->s.ddp_count == cr->dd->ddp_count))
{
- const rvec *fm = as_rvec_array(s_min->f.data());
- const rvec *fb = as_rvec_array(s_b->f.data());
+ const rvec *fm = s_min->f.rvec_array();
+ const rvec *fb = s_b->f.rvec_array();
sum = 0;
int gf = 0;
/* This part of code can be incorrect with DD,
*/
/* Calculate the new direction in p, and the gradient in this direction, gpa */
- rvec *pm = as_rvec_array(s_min->s.cg_p.data());
- const rvec *sfm = as_rvec_array(s_min->f.data());
+ rvec *pm = s_min->s.cg_p.rvec_array();
+ const rvec *sfm = s_min->f.rvec_array();
double gpa = 0;
int gf = 0;
for (int i = 0; i < mdatoms->homenr; i++)
* relative change in coordinate is smaller than precision
*/
minstep = 0;
+ auto s_min_x = makeArrayRef(s_min->s.x);
for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
- tmp = fabs(s_min->s.x[i][m]);
+ tmp = fabs(s_min_x[i][m]);
if (tmp < 1.0)
{
tmp = 1.0;
energyEvaluator.run(s_c, mu_tot, vir, pres, -1, FALSE);
/* Calc derivative along line */
- const rvec *pc = as_rvec_array(s_c->s.cg_p.data());
- const rvec *sfc = as_rvec_array(s_c->f.data());
+ const rvec *pc = s_c->s.cg_p.rvec_array();
+ const rvec *sfc = s_c->f.rvec_array();
double gpc = 0;
for (int i = 0; i < mdatoms->homenr; i++)
{
/* p does not change within a step, but since the domain decomposition
* might change, we have to use cg_p of s_b here.
*/
- const rvec *pb = as_rvec_array(s_b->s.cg_p.data());
- const rvec *sfb = as_rvec_array(s_b->f.data());
+ const rvec *pb = s_b->s.cg_p.rvec_array();
+ const rvec *sfb = s_b->f.rvec_array();
gpb = 0;
for (int i = 0; i < mdatoms->homenr; i++)
{
}
/* Send energies and positions to the IMD client if bIMD is TRUE. */
- if (MASTER(cr) && do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle))
+ if (MASTER(cr) && do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x.rvec_array(), inputrec, 0, wcycle))
{
IMD_send_positions(inputrec->imd);
}
if (vsite)
{
- construct_vsites(vsite, as_rvec_array(state_global->x.data()), 1, nullptr,
+ construct_vsites(vsite, state_global->x.rvec_array(), 1, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state_global->box);
}
point = 0;
// Set initial search direction to the force (-gradient), or 0 for frozen particles.
- real *fInit = static_cast<real *>(as_rvec_array(ems.f.data())[0]);
+ real *fInit = static_cast<real *>(ems.f.rvec_array()[0]);
for (i = 0; i < n; i++)
{
if (!frozen[i])
/* make s a pointer to current search direction - point=0 first time we get here */
s = dx[point];
- real *xx = static_cast<real *>(as_rvec_array(ems.s.x.data())[0]);
- real *ff = static_cast<real *>(as_rvec_array(ems.f.data())[0]);
+ real *xx = static_cast<real *>(ems.s.x.rvec_array()[0]);
+ real *ff = static_cast<real *>(ems.f.rvec_array()[0]);
// calculate line gradient in position A
for (gpa = 0, i = 0; i < n; i++)
// Before taking any steps along the line, store the old position
*last = ems;
- real *lastx = static_cast<real *>(as_rvec_array(last->s.x.data())[0]);
- real *lastf = static_cast<real *>(as_rvec_array(last->f.data())[0]);
+ real *lastx = static_cast<real *>(last->s.x.data()[0]);
+ real *lastf = static_cast<real *>(last->f.data()[0]);
Epot0 = ems.epot;
*sa = ems;
while (maxdelta > inputrec->em_stepsize);
// Take a trial step and move the coordinate array xc[] to position C
- real *xc = static_cast<real *>(as_rvec_array(sc->s.x.data())[0]);
+ real *xc = static_cast<real *>(sc->s.x.rvec_array()[0]);
for (i = 0; i < n; i++)
{
xc[i] = lastx[i] + c*s[i];
energyEvaluator.run(sc, mu_tot, vir, pres, step, FALSE);
// Calc line gradient in position C
- real *fc = static_cast<real *>(as_rvec_array(sc->f.data())[0]);
+ real *fc = static_cast<real *>(sc->f.rvec_array()[0]);
for (gpc = 0, i = 0; i < n; i++)
{
gpc -= s[i]*fc[i]; /* f is negative gradient, thus the sign */
}
// Take a trial step to point B
- real *xb = static_cast<real *>(as_rvec_array(sb->s.x.data())[0]);
+ real *xb = static_cast<real *>(sb->s.x.rvec_array()[0]);
for (i = 0; i < n; i++)
{
xb[i] = lastx[i] + b*s[i];
fnorm = sb->fnorm;
// Calculate gradient in point B
- real *fb = static_cast<real *>(as_rvec_array(sb->f.data())[0]);
+ real *fb = static_cast<real *>(sb->f.rvec_array()[0]);
for (gpb = 0, i = 0; i < n; i++)
{
gpb -= s[i]*fb[i]; /* f is negative gradient, thus the sign */
}
/* Send x and E to IMD client, if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x.rvec_array(), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
/* Send IMD energies and positions, if bIMD is TRUE. */
if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box,
- MASTER(cr) ? as_rvec_array(state_global->x.data()) : nullptr,
+ MASTER(cr) ? state_global->x.rvec_array() : nullptr,
inputrec, 0, wcycle) &&
MASTER(cr))
{
t_graph *graph;
tensor vir, pres;
rvec mu_tot;
- rvec *fneg, *dfdx;
+ rvec *dfdx;
gmx_bool bSparse; /* use sparse matrix storage format */
size_t sz;
gmx_sparsematrix_t * sparse_matrix = nullptr;
vsite, constr, &shellfc,
nfile, fnm, &outf, nullptr, wcycle);
- std::vector<int> atom_index = get_atom_index(top_global);
- snew(fneg, atom_index.size());
+ std::vector<int> atom_index = get_atom_index(top_global);
+ std::vector<gmx::RVec> fneg(atom_index.size(), {0, 0, 0});
snew(dfdx, atom_index.size());
#if !GMX_DOUBLE
************************************************************/
/* Steps are divided one by one over the nodes */
- bool bNS = true;
+ bool bNS = true;
+ auto state_work_x = makeArrayRef(state_work.s.x);
+ auto state_work_f = makeArrayRef(state_work.f);
for (unsigned int aid = cr->nodeid; aid < atom_index.size(); aid += nnodes)
{
size_t atom = atom_index[aid];
int force_flags = GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES;
double t = 0;
- x_min = state_work.s.x[atom][d];
+ x_min = state_work_x[atom][d];
for (unsigned int dx = 0; (dx < 2); dx++)
{
if (dx == 0)
{
- state_work.s.x[atom][d] = x_min - der_range;
+ state_work_x[atom][d] = x_min - der_range;
}
else
{
- state_work.s.x[atom][d] = x_min + der_range;
+ state_work_x[atom][d] = x_min + der_range;
}
/* Make evaluate_energy do a single node force calculation */
enerd,
fcd,
&state_work.s,
- state_work.f,
+ state_work.f.paddedArrayRef(),
vir,
mdatoms,
nrnb,
if (dx == 0)
{
- for (size_t i = 0; i < atom_index.size(); i++)
- {
- copy_rvec(state_work.f[atom_index[i]], fneg[i]);
- }
+ std::copy(state_work_f.begin(), state_work_f.begin()+atom_index.size(), fneg.begin());
}
}
/* x is restored to original */
- state_work.s.x[atom][d] = x_min;
+ state_work_x[atom][d] = x_min;
for (size_t j = 0; j < atom_index.size(); j++)
{
for (size_t k = 0; (k < DIM); k++)
{
dfdx[j][k] =
- -(state_work.f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
+ -(state_work_f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
}
}
exchange_doubles(ms, b, state->nhpres_vxi.data(), nnhpres);
exchange_doubles(ms, b, state->therm_integral.data(), state->ngtc);
exchange_doubles(ms, b, &state->baros_integral, 1);
- exchange_rvecs(ms, b, as_rvec_array(state->x.data()), state->natoms);
- exchange_rvecs(ms, b, as_rvec_array(state->v.data()), state->natoms);
+ exchange_rvecs(ms, b, state->x.rvec_array(), state->natoms);
+ exchange_rvecs(ms, b, state->v.rvec_array(), state->natoms);
}
static void copy_state_serial(const t_state *src, t_state *dest)
}
}
-static void scale_velocities(t_state *state, real fac)
+static void scale_velocities(gmx::ArrayRef<gmx::RVec> velocities, real fac)
{
- int i;
-
- if (as_rvec_array(state->v.data()))
+ for (auto &v : velocities)
{
- for (i = 0; i < state->natoms; i++)
- {
- svmul(fac, state->v[i], state->v[i]);
- }
+ v *= fac;
}
}
* the velocities. */
if (re->type == ereTEMP || re->type == ereTL)
{
- scale_velocities(state, std::sqrt(re->q[ereTEMP][replica_id]/re->q[ereTEMP][re->destinations[replica_id]]));
+ scale_velocities(state->v,
+ std::sqrt(re->q[ereTEMP][replica_id]/re->q[ereTEMP][re->destinations[replica_id]]));
}
}
const t_forcerec &forceRec,
t_graph *graph)
{
- for (int i = 0; i < globalState->natoms; i++)
- {
- copy_rvec(rerunFrame.x[i], globalState->x[i]);
- }
+ auto x = makeArrayRef(globalState->x);
+ auto rerunX = arrayRefFromArray(reinterpret_cast<gmx::RVec *>(rerunFrame.x), globalState->natoms);
+ std::copy(rerunX.begin(), rerunX.end(), x.begin());
copy_mat(rerunFrame.box, globalState->box);
if (constructVsites)
/* Following is necessary because the graph may get out of sync
* with the coordinates if we only have every N'th coordinate set
*/
- mk_mshift(nullptr, graph, forceRec.ePBC, globalState->box, as_rvec_array(globalState->x.data()));
+ mk_mshift(nullptr, graph, forceRec.ePBC, globalState->box, globalState->x.rvec_array());
shift_self(graph, globalState->box, as_rvec_array(globalState->x.data()));
}
- construct_vsites(vsite, as_rvec_array(globalState->x.data()), timeStep, as_rvec_array(globalState->v.data()),
+ construct_vsites(vsite, globalState->x.rvec_array(), timeStep, globalState->v.rvec_array(),
idef.iparams, idef.il,
forceRec.ePBC, forceRec.bMolPBC, nullptr, globalState->box);
if (graph)
{
- unshift_self(graph, globalState->box, as_rvec_array(globalState->x.data()));
+ unshift_self(graph, globalState->box, globalState->x.rvec_array());
}
}
}
// alias to avoid a large ripple of nearly useless changes.
// t_inputrec is being replaced by IMdpOptionsProvider, so this
// will go away eventually.
- t_inputrec *ir = inputrec;
- gmx_mdoutf *outf = nullptr;
- int64_t step, step_rel;
- double t, lam0[efptNR];
- bool isLastStep = false;
- bool doFreeEnergyPerturbation = false;
- int force_flags;
- tensor force_vir, shake_vir, total_vir, pres;
- t_trxstatus *status;
- rvec mu_tot;
- t_trxframe rerun_fr;
- gmx_localtop_t *top;
- t_mdebin *mdebin = nullptr;
- gmx_enerdata_t *enerd;
- PaddedRVecVector f {};
- gmx_global_stat_t gstat;
- t_graph *graph = nullptr;
- gmx_groups_t *groups;
- gmx_ekindata_t *ekind;
- gmx_shellfc_t *shellfc;
-
- double cycles;
+ t_inputrec *ir = inputrec;
+ gmx_mdoutf *outf = nullptr;
+ int64_t step, step_rel;
+ double t, lam0[efptNR];
+ bool isLastStep = false;
+ bool doFreeEnergyPerturbation = false;
+ int force_flags;
+ tensor force_vir, shake_vir, total_vir, pres;
+ t_trxstatus *status;
+ rvec mu_tot;
+ t_trxframe rerun_fr;
+ gmx_localtop_t *top;
+ t_mdebin *mdebin = nullptr;
+ gmx_enerdata_t *enerd;
+ PaddedVector<gmx::RVec> f {};
+ gmx_global_stat_t gstat;
+ t_graph *graph = nullptr;
+ gmx_groups_t *groups;
+ gmx_ekindata_t *ekind;
+ gmx_shellfc_t *shellfc;
+
+ double cycles;
/* Domain decomposition could incorrectly miss a bonded
interaction, but checking for that requires a global
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- f.resize(gmx::paddedRVecVectorSize(state_global->natoms));
+ f.resizeWithPadding(state_global->natoms);
/* Copy the pointer to the global state */
state = state_global;
/* Make molecules whole at start of run */
if (fr->ePBC != epbcNONE)
{
- rvec *xGlobal = as_rvec_array(globalState->x.data());
- do_pbc_first_mtop(fplog, inputrec->ePBC, box, &mtop, xGlobal);
+ do_pbc_first_mtop(fplog, inputrec->ePBC, box, &mtop, globalState->x.rvec_array());
}
if (vsite)
{
void
Integrator::do_tpi()
{
- gmx_localtop_t *top;
- gmx_groups_t *groups;
- gmx_enerdata_t *enerd;
- PaddedRVecVector f {};
- real lambda, t, temp, beta, drmax, epot;
- double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
- t_trxstatus *status;
- t_trxframe rerun_fr;
- gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS;
- tensor force_vir, shake_vir, vir, pres;
- int cg_tp, a_tp0, a_tp1, ngid, gid_tp, nener, e;
- rvec *x_mol;
- rvec mu_tot, x_init, dx, x_tp;
- int nnodes, frame;
- int64_t frame_step_prev, frame_step;
- int64_t nsteps, stepblocksize = 0, step;
- int64_t seed;
- int i;
- FILE *fp_tpi = nullptr;
- char *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
- double dbl, dump_ener;
- gmx_bool bCavity;
- int nat_cavity = 0, d;
- real *mass_cavity = nullptr, mass_tot;
- int nbin;
- double invbinw, *bin, refvolshift, logV, bUlogV;
- real prescorr, enercorr, dvdlcorr;
- gmx_bool bEnergyOutOfBounds;
- const char *tpid_leg[2] = {"direct", "reweighted"};
- auto mdatoms = mdAtoms->mdatoms();
+ gmx_localtop_t *top;
+ gmx_groups_t *groups;
+ gmx_enerdata_t *enerd;
+ PaddedVector<gmx::RVec> f {};
+ real lambda, t, temp, beta, drmax, epot;
+ double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
+ t_trxstatus *status;
+ t_trxframe rerun_fr;
+ gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS;
+ tensor force_vir, shake_vir, vir, pres;
+ int cg_tp, a_tp0, a_tp1, ngid, gid_tp, nener, e;
+ rvec *x_mol;
+ rvec mu_tot, x_init, dx, x_tp;
+ int nnodes, frame;
+ int64_t frame_step_prev, frame_step;
+ int64_t nsteps, stepblocksize = 0, step;
+ int64_t seed;
+ int i;
+ FILE *fp_tpi = nullptr;
+ char *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
+ double dbl, dump_ener;
+ gmx_bool bCavity;
+ int nat_cavity = 0, d;
+ real *mass_cavity = nullptr, mass_tot;
+ int nbin;
+ double invbinw, *bin, refvolshift, logV, bUlogV;
+ real prescorr, enercorr, dvdlcorr;
+ gmx_bool bEnergyOutOfBounds;
+ const char *tpid_leg[2] = {"direct", "reweighted"};
+ auto mdatoms = mdAtoms->mdatoms();
GMX_UNUSED_VALUE(outputProvider);
snew(enerd, 1);
init_enerdata(groups->grps[egcENER].nr, inputrec->fepvals->n_lambda, enerd);
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- f.resize(gmx::paddedRVecVectorSize(top_global->natoms));
+ f.resizeWithPadding(top_global->natoms);
/* Print to log file */
walltime_accounting_start_time(walltime_accounting);
bDispCorr = (inputrec->eDispCorr != edispcNO);
bCharge = FALSE;
+ auto x = makeArrayRef(state_global->x);
for (i = a_tp0; i < a_tp1; i++)
{
/* Copy the coordinates of the molecule to be insterted */
- copy_rvec(state_global->x[i], x_mol[i-a_tp0]);
+ copy_rvec(x[i], x_mol[i-a_tp0]);
/* Check if we need to print electrostatic energies */
bCharge |= (mdatoms->chargeA[i] != 0 ||
((mdatoms->chargeB != nullptr) && mdatoms->chargeB[i] != 0));
}
bRFExcl = (bCharge && EEL_RF(fr->ic->eeltype));
- calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), as_rvec_array(state_global->x.data()), fr->cg_cm);
+ calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), state_global->x.rvec_array(), fr->cg_cm);
if (bCavity)
{
if (norm(fr->cg_cm[cg_tp]) > 0.5*inputrec->rlist && fplog)
}
/* Copy the coordinates from the input trajectory */
+ auto x = makeArrayRef(state_global->x);
for (i = 0; i < rerun_fr.natoms; i++)
{
- copy_rvec(rerun_fr.x[i], state_global->x[i]);
+ copy_rvec(rerun_fr.x[i], x[i]);
}
copy_mat(rerun_fr.box, state_global->box);
if (a_tp1 - a_tp0 == 1)
{
/* Insert a single atom, just copy the insertion location */
- copy_rvec(x_tp, state_global->x[a_tp0]);
+ copy_rvec(x_tp, x[a_tp0]);
}
else
{
/* Copy the coordinates from the top file */
for (i = a_tp0; i < a_tp1; i++)
{
- copy_rvec(x_mol[i-a_tp0], state_global->x[i]);
+ copy_rvec(x_mol[i-a_tp0], x[i]);
}
/* Rotate the molecule randomly */
real angleX = 2*M_PI*dist(rng);
real angleY = 2*M_PI*dist(rng);
real angleZ = 2*M_PI*dist(rng);
- rotate_conf(a_tp1-a_tp0, as_rvec_array(state_global->x.data())+a_tp0, nullptr,
+ rotate_conf(a_tp1-a_tp0, state_global->x.rvec_array()+a_tp0, nullptr,
angleX, angleY, angleZ);
/* Shift to the insertion location */
for (i = a_tp0; i < a_tp1; i++)
{
- rvec_inc(state_global->x[i], x_tp);
+ rvec_inc(x[i], x_tp);
}
}
cr->nnodes = 1;
do_force(fplog, cr, ms, inputrec, nullptr, nullptr,
step, nrnb, wcycle, top, &top_global->groups,
- state_global->box, state_global->x, &state_global->hist,
- f, force_vir, mdatoms, enerd, fcd,
+ state_global->box, state_global->x.paddedArrayRef(), &state_global->hist,
+ f.paddedArrayRef(), force_vir, mdatoms, enerd, fcd,
state_global->lambda,
nullptr, fr, nullptr, mu_tot, t, nullptr,
GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
{
sprintf(str, "t%g_step%d.pdb", t, static_cast<int>(step));
sprintf(str2, "t: %f step %d ener: %f", t, static_cast<int>(step), epot);
- write_sto_conf_mtop(str, str2, top_global, as_rvec_array(state_global->x.data()), as_rvec_array(state_global->v.data()),
+ write_sto_conf_mtop(str, str2, top_global, state_global->x.rvec_array(), state_global->v.rvec_array(),
inputrec->ePBC, state_global->box);
}
#ifndef GMX_MDTYPES_TYPES_FORCEREC_H
#define GMX_MDTYPES_TYPES_FORCEREC_H
-#include "gromacs/math/paddedvector.h"
+#include <vector>
+
#include "gromacs/math/vectypes.h"
#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
{
state->natoms = natoms;
- /* We need padding, since we might use SIMD access */
- const size_t paddedSize = gmx::paddedRVecVectorSize(state->natoms);
-
+ /* We need padding, since we might use SIMD access, but the
+ * containers here all ensure that. */
if (state->flags & (1 << estX))
{
- state->x.resize(paddedSize);
+ state->x.resizeWithPadding(natoms);
}
if (state->flags & (1 << estV))
{
- state->v.resize(paddedSize);
+ state->v.resizeWithPadding(natoms);
}
if (state->flags & (1 << estCGP))
{
- state->cg_p.resize(paddedSize);
+ state->cg_p.resizeWithPadding(natoms);
}
}
if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX)))
{
fprintf(stdout, "comparing x\n");
- cmp_rvecs(stdout, "x", st1->natoms, as_rvec_array(st1->x.data()), as_rvec_array(st2->x.data()), bRMSD, ftol, abstol);
+ cmp_rvecs(stdout, "x", st1->natoms, st1->x.rvec_array(), st2->x.rvec_array(), bRMSD, ftol, abstol);
}
if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV)))
{
fprintf(stdout, "comparing v\n");
- cmp_rvecs(stdout, "v", st1->natoms, as_rvec_array(st1->v.data()), as_rvec_array(st2->v.data()), bRMSD, ftol, abstol);
+ cmp_rvecs(stdout, "v", st1->natoms, st1->v.rvec_array(), st2->v.rvec_array(), bRMSD, ftol, abstol);
}
}
}
real veta; //!< Trotter based isotropic P-coupling
real vol0; //!< Initial volume,required for computing MTTK conserved quantity
gmx::HostVector<gmx::RVec> x; //!< The coordinates (natoms)
- PaddedRVecVector v; //!< The velocities (natoms)
- PaddedRVecVector cg_p; //!< p vector for conjugate gradient minimization
+ PaddedVector<gmx::RVec> v; //!< The velocities (natoms)
+ PaddedVector<gmx::RVec> cg_p; //!< p vector for conjugate gradient minimization
ekinstate_t ekinstate; //!< The state of the kinetic energy
{
if (state)
{
- return gmx::constArrayRefFromArray(state->x.data(), state->natoms);
+ return gmx::makeConstArrayRef(state->x).subArray(0, state->natoms);
}
else
{
/* Remove pbc, make molecule whole.
* When ir->bContinuation=TRUE this has already been done, but ok. */
snew(x_pbc, mtop->natoms);
- copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
+ copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
/* All molecules will be whole now, but not necessarily in the home box.
* Additionally, if a rotation group consists of more than one molecule
}
swapstate = oh->swapHistory.get();
- init_swapstate(swapstate, sc, mtop, as_rvec_array(globalState->x.data()), globalState->box, ir);
+ init_swapstate(swapstate, sc, mtop, globalState->x.rvec_array(), globalState->box, ir);
}
/* After init_swapstate we have a set of (old) whole positions for our
else
{
fprintf(stderr, "%s Determining initial numbers of ions per compartment.\n", SwS);
- get_initial_ioncounts(ir, as_rvec_array(globalState->x.data()), globalState->box, cr, mdrunOptions.rerun);
+ get_initial_ioncounts(ir, globalState->x.rvec_array(), globalState->box, cr, mdrunOptions.rerun);
}
/* Prepare (further) checkpoint writes ... */
{
fprintf(stderr, "Will write subset %s of original tpx containing %d "
"atoms\n", grpname, gnx);
- reduce_topology_x(gnx, index, &mtop, as_rvec_array(state.x.data()), as_rvec_array(state.v.data()));
+ reduce_topology_x(gnx, index, &mtop, state.x.rvec_array(), state.v.rvec_array());
state.natoms = gnx;
}
else if (bZeroQ)
pr_doubles(stdout, indent, "nosehoover_xi", state.nosehoover_xi.data(), state.ngtc);
/*pr_doubles(stdout,indent,"nosehoover_vxi",state.nosehoover_vxi,state.ngtc);*/
/*pr_doubles(stdout,indent,"therm_integral",state.therm_integral,state.ngtc);*/
- pr_rvecs(stdout, indent, "x", tpx.bX ? as_rvec_array(state.x.data()) : nullptr, state.natoms);
- pr_rvecs(stdout, indent, "v", tpx.bV ? as_rvec_array(state.v.data()) : nullptr, state.natoms);
+ pr_rvecs(stdout, indent, "x", tpx.bX ? state.x.rvec_array() : nullptr, state.natoms);
+ pr_rvecs(stdout, indent, "v", tpx.bV ? state.v.rvec_array() : nullptr, state.natoms);
}
groups = &mtop.groups;
/*! \def gmx_real_fullprecision_pfmt
* \brief Format string for full `real` precision.
*/
+/*! \def GMX_FLOAT_MAX_SIMD_WIDTH
+ * \brief The maximum supported number of `float` elements in a SIMD register.
+ */
+/*! \def GMX_DOUBLE_MAX_SIMD_WIDTH
+ * \brief The maximum supported number of `double` elements in a SIMD register.
+ */
/*! \def GMX_REAL_MAX_SIMD_WIDTH
* \brief The maximum supported number of `real` elements in a SIMD register.
*/
+
+#define GMX_FLOAT_MAX_SIMD_WIDTH 16
+#define GMX_DOUBLE_MAX_SIMD_WIDTH 8
+
#if GMX_DOUBLE
#ifndef HAVE_REAL
#define GMX_REAL_MAX GMX_DOUBLE_MAX
#define GMX_REAL_NEGZERO GMX_DOUBLE_NEGZERO
#define gmx_real_fullprecision_pfmt "%21.14e"
-#define GMX_REAL_MAX_SIMD_WIDTH 8
+#define GMX_REAL_MAX_SIMD_WIDTH GMX_DOUBLE_MAX_SIMD_WIDTH
#else /* GMX_DOUBLE */
#define GMX_REAL_MAX GMX_FLOAT_MAX
#define GMX_REAL_NEGZERO GMX_FLOAT_NEGZERO
#define gmx_real_fullprecision_pfmt "%14.7e"
-#define GMX_REAL_MAX_SIMD_WIDTH 16
+#define GMX_REAL_MAX_SIMD_WIDTH GMX_FLOAT_MAX_SIMD_WIDTH
#endif /* GMX_DOUBLE */
# pseudo-library for code for mdrun
$<TARGET_OBJECTS:mdrun_objlib>
)
-if (GMX_GPU AND GMX_USE_OPENCL)
- gmx_register_gtest_test(${testname} ${exename} OPENMP_THREADS 2 OCL_INTEGRATION_TEST)
-else()
- gmx_register_gtest_test(${testname} ${exename} OPENMP_THREADS 2 INTEGRATION_TEST)
-endif()
+gmx_register_gtest_test(${testname} ${exename} OPENMP_THREADS 2 INTEGRATION_TEST)
# To avoid running into test timeouts, some end-to-end tests of mdrun
# functionality are split off. This can be rearranged in future as we
# pseudo-library for code for mdrun
$<TARGET_OBJECTS:mdrun_objlib>
)
-if (GMX_GPU AND GMX_USE_OPENCL)
- gmx_register_gtest_test(${testname} ${exename} OPENMP_THREADS 2 OCL_INTEGRATION_TEST)
-else()
- gmx_register_gtest_test(${testname} ${exename} OPENMP_THREADS 2 INTEGRATION_TEST)
-endif()
+gmx_register_gtest_test(${testname} ${exename} OPENMP_THREADS 2 INTEGRATION_TEST)
# threads (when supported by the build configuration)
# MPI_RANKS <N> declares the requirement to run the test binary with N ranks
# INTEGRATION_TEST requires the use of the IntegrationTest label in CTest
-# OCL_INTEGRATION_TEST requires the use of the IntegrationTest label in CTest,
-# only difference is 2x longer timeout as OpenCL JIT can be slow
# SLOW_TEST requires the use of the SlowTest label in CTest, and
# increase the length of the ctest timeout.
#
# ranks.
function (gmx_register_gtest_test NAME EXENAME)
if (GMX_BUILD_UNITTESTS AND BUILD_TESTING)
- set(_options INTEGRATION_TEST OCL_INTEGRATION_TEST SLOW_TEST)
+ set(_options INTEGRATION_TEST SLOW_TEST)
set(_one_value_args MPI_RANKS OPENMP_THREADS)
cmake_parse_arguments(ARG "${_options}" "${_one_value_args}" "" ${ARGN})
set(_xml_path ${CMAKE_BINARY_DIR}/Testing/Temporary/${NAME}.xml)
set(_labels GTest)
set(_timeout 30)
- if (ARG_INTEGRATION_TEST OR ARG_OCL_INTEGRATION_TEST)
+ if (ARG_INTEGRATION_TEST)
list(APPEND _labels IntegrationTest)
- if (ARG_OCL_INTEGRATION_TEST)
+ # Slow build configurations should have longer timeouts.
+ # Both OpenCL (from JIT) and ThreadSanitizer (from how it
+ # checks) can take signficantly more time than other
+ # configurations.
+ if (GMX_USE_OPENCL)
set(_timeout 240)
+ elseif (${CMAKE_BUILD_TYPE} STREQUAL TSAN)
+ set(_timeout 300)
else()
set(_timeout 120)
endif()
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff