have compute_reduce require either peratom or local inputs

This commit is contained in:
Steve Plimpton
2023-08-17 16:12:14 -06:00
parent 0d739439c7
commit 299eda8ca3
5 changed files with 149 additions and 106 deletions

View File

@ -37,13 +37,16 @@ Syntax
v_name = per-atom vector calculated by an atom-style variable with name
* zero or more keyword/args pairs may be appended
* keyword = *replace*
* keyword = *replace* or *inputs*
.. parsed-literal::
*replace* args = vec1 vec2
vec1 = reduced value from this input vector will be replaced
vec2 = replace it with vec1[N] where N is index of max/min value from vec2
*inputs* arg = peratom or local
peratom = all inputs are per-atom quantities (default)
local = all input are local quantities
Examples
""""""""
@ -61,26 +64,30 @@ Description
Define a calculation that "reduces" one or more vector inputs into
scalar values, one per listed input. The inputs can be per-atom or
local quantities; they cannot be global quantities. Atom attributes
are per-atom quantities, :doc:`computes <compute>` and :doc:`fixes <fix>`
may generate any of the three kinds of quantities, and :doc:`atom-style variables <variable>` generate per-atom quantities. See the
:doc:`variable <variable>` command and its special functions which can
perform the same operations as the compute reduce command on global
vectors.
local quantities and must all be the same kind (per-atom or local);
see discussion of the optional *inputs* keyword below.
Atom attributes are per-atom quantities, :doc:`computes <compute>` and
:doc:`fixes <fix>` can generate either per-atom or local quantities,
and :doc:`atom-style variables <variable>` generate per-atom
quantities. See the :doc:`variable <variable>` command and its
special functions which can perform the same reduction operations as
the compute reduce command on global vectors.
The reduction operation is specified by the *mode* setting. The *sum*
option adds the values in the vector into a global total. The *min*
or *max* options find the minimum or maximum value across all vector
values. The *minabs* or *maxabs* options find the minimum or maximum
value across all absolute vector values. The *ave* setting adds the
vector values into a global total, then divides by the number of values
in the vector. The *sumsq* option sums the square of the values in the
vector into a global total. The *avesq* setting does the same as *sumsq*,
then divides the sum of squares by the number of values. The last two options
can be useful for calculating the variance of some quantity (e.g., variance =
sumsq :math:`-` ave\ :math:`^2`). The *sumabs* option sums the absolute
values in the vector into a global total. The *aveabs* setting does the same
as *sumabs*, then divides the sum of absolute values by the number of
vector values into a global total, then divides by the number of
values in the vector. The *sumsq* option sums the square of the
values in the vector into a global total. The *avesq* setting does
the same as *sumsq*, then divides the sum of squares by the number of
values. The last two options can be useful for calculating the
variance of some quantity (e.g., variance = sumsq :math:`-` ave\
:math:`^2`). The *sumabs* option sums the absolute values in the
vector into a global total. The *aveabs* setting does the same as
*sumabs*, then divides the sum of absolute values by the number of
values.
Each listed input is operated on independently. For per-atom inputs,
@ -123,52 +130,54 @@ array with six columns:
----------
The atom attribute values (*x*, *y*, *z*, *vx*, *vy*, *vz*, *fx*, *fy*, and
*fz*) are self-explanatory. Note that other atom attributes can be used as
inputs to this fix by using the
:doc:`compute property/atom <compute_property_atom>` command and then specifying
an input value from that compute.
The atom attribute values (*x*, *y*, *z*, *vx*, *vy*, *vz*, *fx*,
*fy*, and *fz*) are self-explanatory. Note that other atom attributes
can be used as inputs to this fix by using the :doc:`compute
property/atom <compute_property_atom>` command and then specifying an
input value from that compute.
If a value begins with "c\_", a compute ID must follow which has been
previously defined in the input script. Computes can generate
per-atom or local quantities. See the individual
:doc:`compute <compute>` page for details. If no bracketed integer
is appended, the vector calculated by the compute is used. If a
bracketed integer is appended, the Ith column of the array calculated
by the compute is used. Users can also write code for their own
compute styles and :doc:`add them to LAMMPS <Modify>`. See the
discussion above for how :math:`I` can be specified with a wildcard asterisk
to effectively specify multiple values.
previously defined in the input script. Valid computes can generate
per-atom or local quantities. See the individual :doc:`compute
<compute>` page for details. If no bracketed integer is appended, the
vector calculated by the compute is used. If a bracketed integer is
appended, the Ith column of the array calculated by the compute is
used. Users can also write code for their own compute styles and
:doc:`add them to LAMMPS <Modify>`. See the discussion above for how
:math:`I` can be specified with a wildcard asterisk to effectively
specify multiple values.
If a value begins with "f\_", a fix ID must follow which has been
previously defined in the input script. Fixes can generate per-atom
or local quantities. See the individual :doc:`fix <fix>` page for
details. Note that some fixes only produce their values on certain
timesteps, which must be compatible with when compute reduce
previously defined in the input script. Valid fixes can generate
per-atom or local quantities. See the individual :doc:`fix <fix>`
page for details. Note that some fixes only produce their values on
certain timesteps, which must be compatible with when compute reduce
references the values, else an error results. If no bracketed integer
is appended, the vector calculated by the fix is used. If a bracketed
integer is appended, the Ith column of the array calculated by the fix
is used. Users can also write code for their own fix style and
:doc:`add them to LAMMPS <Modify>`. See the discussion above for how
:math:`I` can be specified with a wildcard asterisk to effectively specify
multiple values.
:math:`I` can be specified with a wildcard asterisk to effectively
specify multiple values.
If a value begins with "v\_", a variable name must follow which has
been previously defined in the input script. It must be an
:doc:`atom-style variable <variable>`. Atom-style variables can
reference thermodynamic keywords and various per-atom attributes, or
invoke other computes, fixes, or variables when they are evaluated, so
this is a very general means of generating per-atom quantities to reduce.
this is a very general means of generating per-atom quantities to
reduce.
----------
If the *replace* keyword is used, two indices *vec1* and *vec2* are
specified, where each index ranges from 1 to the number of input values.
The replace keyword can only be used if the *mode* is *min* or *max*\ .
It works as follows. A min/max is computed as usual on the *vec2*
input vector. The index :math:`N` of that value within *vec2* is also stored.
Then, instead of performing a min/max on the *vec1* input vector, the
stored index is used to select the :math:`N`\ th element of the *vec1* vector.
specified, where each index ranges from 1 to the number of input
values. The replace keyword can only be used if the *mode* is *min*
or *max*\ . It works as follows. A min/max is computed as usual on
the *vec2* input vector. The index :math:`N` of that value within
*vec2* is also stored. Then, instead of performing a min/max on the
*vec1* input vector, the stored index is used to select the :math:`N`\
th element of the *vec1* vector.
Thus, for example, if you wish to use this compute to find the bond
with maximum stretch, you can do it as follows:
@ -190,6 +199,14 @@ information in this context, the *replace* keywords will extract the
atom IDs for the two atoms in the bond of maximum stretch. These atom
IDs and the bond stretch will be printed with thermodynamic output.
The *inputs* keyword allows selection of whether all the inputs are
per-atom or local quantities. As noted above, all the inputs must be
the same kind (per-atom or local). Per-atom is the default setting.
If a compute or fix is specified as an input, it must produce per-atom
or local data to match this setting. If it produces both, e.g. for
the :doc:`compute voronoi/atom <compute_voronoi_atom>` command, then
this keyword selects between them.
----------
If a single input is specified this compute produces a global scalar
@ -197,34 +214,35 @@ value. If multiple inputs are specified, this compute produces a
global vector of values, the length of which is equal to the number of
inputs specified.
As discussed below, for the *sum*, *sumabs*, and *sumsq* modes, the value(s)
produced by this compute are all "extensive", meaning their value
scales linearly with the number of atoms involved. If normalized
values are desired, this compute can be accessed by the
As discussed below, for the *sum*, *sumabs*, and *sumsq* modes, the
value(s) produced by this compute are all "extensive", meaning their
value scales linearly with the number of atoms involved. If
normalized values are desired, this compute can be accessed by the
:doc:`thermo_style custom <thermo_style>` command with
:doc:`thermo_modify norm yes <thermo_modify>` set as an option.
Or it can be accessed by a
:doc:`variable <variable>` that divides by the appropriate atom count.
:doc:`thermo_modify norm yes <thermo_modify>` set as an option. Or it
can be accessed by a :doc:`variable <variable>` that divides by the
appropriate atom count.
----------
Output info
"""""""""""
This compute calculates a global scalar if a single input value is specified
or a global vector of length :math:`N`, where :math:`N` is the number of
inputs, and which can be accessed by indices 1 to :math:`N`. These values can
be used by any command that uses global scalar or vector values from a
compute as input. See the :doc:`Howto output <Howto_output>` doc page
for an overview of LAMMPS output options.
This compute calculates a global scalar if a single input value is
specified or a global vector of length :math:`N`, where :math:`N` is
the number of inputs, and which can be accessed by indices 1 to
:math:`N`. These values can be used by any command that uses global
scalar or vector values from a compute as input. See the :doc:`Howto
output <Howto_output>` doc page for an overview of LAMMPS output
options.
All the scalar or vector values calculated by this compute are
"intensive", except when the *sum*, *sumabs*, or *sumsq* modes are used on
per-atom or local vectors, in which case the calculated values are
"extensive".
The scalar or vector values will be in whatever :doc:`units <units>` the
quantities being reduced are in.
The scalar or vector values will be in whatever :doc:`units <units>`
the quantities being reduced are in.
Restrictions
""""""""""""
@ -238,4 +256,4 @@ Related commands
Default
"""""""
none
The default value for the *inputs* keyword is peratom.

View File

@ -843,7 +843,7 @@ stress/atom <compute_stress_atom>` commands. The former can be
accessed by :doc:`thermodynamic output <thermo_style>`. The default
setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.
All of the *rigid* styles (not the *rigid/small* styles) compute a
All of the *rigid* styles (but not the *rigid/small* styles) compute a
global array of values which can be accessed by various :doc:`output
commands <Howto_output>`. Similar information about the bodies
defined by the *rigid/small* styles can be accessed via the
@ -887,7 +887,8 @@ Restrictions
""""""""""""
These fixes are all part of the RIGID package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Assigning a temperature via the :doc:`velocity create <velocity>`
command to a system with :doc:`rigid bodies <fix_rigid>` may not have

View File

@ -31,12 +31,16 @@
using namespace LAMMPS_NS;
enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduceRegion
#define BIG 1.0e20
//----------------------------------------------------------------
void abs_max(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/)
{
// r is the already reduced value, n is the new value
double n = std::fabs(*(double *) in), r = *(double *) inout;
double m;
@ -47,9 +51,11 @@ void abs_max(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/)
}
*(double *) inout = m;
}
void abs_min(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/)
{
// r is the already reduced value, n is the new value
double n = std::fabs(*(double *) in), r = *(double *) inout;
double m;
@ -68,6 +74,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
owner(nullptr), idregion(nullptr), region(nullptr), varatom(nullptr)
{
int iarg = 0;
if (strcmp(style, "reduce") == 0) {
if (narg < 5) utils::missing_cmd_args(FLERR, "compute reduce", error);
iarg = 3;
@ -128,42 +135,52 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
// parse values
input_mode = UNDECIDED;
values.clear();
nvalues = 0;
for (int iarg = 0; iarg < nargnew; ++iarg) {
value_t val;
val.id = "";
val.flavor = 0;
val.val.c = nullptr;
if (strcmp(arg[iarg], "x") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::X;
val.argindex = 0;
} else if (strcmp(arg[iarg], "y") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::X;
val.argindex = 1;
} else if (strcmp(arg[iarg], "z") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::X;
val.argindex = 2;
} else if (strcmp(arg[iarg], "vx") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::V;
val.argindex = 0;
} else if (strcmp(arg[iarg], "vy") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::V;
val.argindex = 1;
} else if (strcmp(arg[iarg], "vz") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::V;
val.argindex = 2;
} else if (strcmp(arg[iarg], "fx") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::F;
val.argindex = 0;
} else if (strcmp(arg[iarg], "fy") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::F;
val.argindex = 1;
} else if (strcmp(arg[iarg], "fz") == 0) {
input_mode = PERATOM;
val.which = ArgInfo::F;
val.argindex = 2;
@ -207,6 +224,14 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR, "Compute {} replace column already used for another replacement");
replace[col1] = col2;
iarg += 2;
} else if (strcmp(arg[iarg], "inputs") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, mycmd + " inputs", error);
if (strcmp(arg[iarg+1], "peratom") == 0) input_mode = PERATOM;
else if (strcmp(arg[iarg+1], "local") == 0) {
if (input_mode == PERATOM)
error->all(FLERR,"Compute {} inputs must be all peratom or all local");
input_mode = LOCAL;
}
} else
error->all(FLERR, "Unknown compute {} keyword: {}", style, arg[iarg]);
}
@ -231,66 +256,64 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
// setup and error check
for (auto &val : values) {
if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F)
val.flavor = PERATOM;
else if (val.which == ArgInfo::COMPUTE) {
if (val.which == ArgInfo::COMPUTE) {
val.val.c = modify->get_compute_by_id(val.id);
if (!val.val.c)
error->all(FLERR, "Compute ID {} for compute {} does not exist", val.id, style);
if (val.val.c->peratom_flag) {
val.flavor = PERATOM;
if (input_mode == PERATOM) {
if (!val.val.c->peratom_flag)
error->all(FLERR, "Compute {} compute {} does not calculate per-atom values", style, val.id);
if (val.argindex == 0 && val.val.c->size_peratom_cols != 0)
error->all(FLERR, "Compute {} compute {} does not calculate a per-atom vector", style,
val.id);
error->all(FLERR, "Compute {} compute {} does not calculate a per-atom vector", style, val.id);
if (val.argindex && val.val.c->size_peratom_cols == 0)
error->all(FLERR, "Compute {} compute {} does not calculate a per-atom array", style,
val.id);
error->all(FLERR, "Compute {} compute {} does not calculate a per-atom array", style, val.id);
if (val.argindex && val.argindex > val.val.c->size_peratom_cols)
error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id);
} else if (val.val.c->local_flag) {
val.flavor = LOCAL;
} else if (input_mode == LOCAL) {
if (!val.val.c->peratom_flag)
error->all(FLERR, "Compute {} compute {} does not calculate local values", style, val.id);
if (val.argindex == 0 && val.val.c->size_local_cols != 0)
error->all(FLERR, "Compute {} compute {} does not calculate a local vector", style,
val.id);
error->all(FLERR, "Compute {} compute {} does not calculate a local vector", style, val.id);
if (val.argindex && val.val.c->size_local_cols == 0)
error->all(FLERR, "Compute {} compute {} does not calculate a local array", style,
val.id);
error->all(FLERR, "Compute {} compute {} does not calculate a local array", style, val.id);
if (val.argindex && val.argindex > val.val.c->size_local_cols)
error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id);
} else
error->all(FLERR, "Compute {} compute {} calculates global values", style, val.id);
}
} else if (val.which == ArgInfo::FIX) {
val.val.f = modify->get_fix_by_id(val.id);
if (!val.val.f) error->all(FLERR, "Fix ID {} for compute {} does not exist", val.id, style);
if (val.val.f->peratom_flag) {
val.flavor = PERATOM;
if (input_mode == PERATOM) {
if (!val.val.f->peratom_flag)
error->all(FLERR, "Compute {} fix {} does not calculate per-atom values", style, val.id);
if (val.argindex == 0 && (val.val.f->size_peratom_cols != 0))
error->all(FLERR, "Compute {} fix {} does not calculate a per-atom vector", style,
val.id);
error->all(FLERR, "Compute {} fix {} does not calculate a per-atom vector", style, val.id);
if (val.argindex && (val.val.f->size_peratom_cols == 0))
error->all(FLERR, "Compute {} fix {} does not calculate a per-atom array", style, val.id);
if (val.argindex && (val.argindex > val.val.f->size_peratom_cols))
error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id);
} else if (val.val.f->local_flag) {
val.flavor = LOCAL;
} else if (input_mode == LOCAL) {
if (!val.val.f->local_flag)
error->all(FLERR, "Compute {} fix {} does not calculate local values", style, val.id);
if (val.argindex == 0 && (val.val.f->size_local_cols != 0))
error->all(FLERR, "Compute {} fix {} does not calculate a local vector", style, val.id);
if (val.argindex && (val.val.f->size_local_cols == 0))
error->all(FLERR, "Compute {} fix {} does not calculate a local array", style, val.id);
if (val.argindex && (val.argindex > val.val.f->size_local_cols))
error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id);
} else
error->all(FLERR, "Compute {} fix {} calculates global values", style, val.id);
}
} else if (val.which == ArgInfo::VARIABLE) {
if (input_mode == LOCAL) error->all(FLERR,"Compute {} inputs must be all local");
val.val.v = input->variable->find(val.id.c_str());
if (val.val.v < 0)
error->all(FLERR, "Variable name {} for compute {} does not exist", val.id, style);
if (input->variable->atomstyle(val.val.v) == 0)
error->all(FLERR, "Compute {} variable {} is not atom-style variable", style, val.id);
val.flavor = PERATOM;
}
}
@ -512,7 +535,7 @@ double ComputeReduce::compute_one(int m, int flag)
} else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
val.val.c->compute_peratom();
val.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
@ -537,7 +560,7 @@ double ComputeReduce::compute_one(int m, int flag)
one = carray_atom[flag][aidxm1];
}
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) {
val.val.c->compute_local();
val.val.c->invoked_flag |= Compute::INVOKED_LOCAL;
@ -567,7 +590,7 @@ double ComputeReduce::compute_one(int m, int flag)
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style);
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
if (aidx == 0) {
double *fix_vector = val.val.f->vector_atom;
if (flag < 0) {
@ -585,7 +608,7 @@ double ComputeReduce::compute_one(int m, int flag)
one = fix_array[flag][aidxm1];
}
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
if (aidx == 0) {
double *fix_vector = val.val.f->vector_local;
int n = val.val.f->size_local_rows;
@ -632,18 +655,18 @@ bigint ComputeReduce::count(int m)
if ((val.which == ArgInfo::X) || (val.which == ArgInfo::V) || (val.which == ArgInfo::F))
return group->count(igroup);
else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
return group->count(igroup);
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
bigint ncount = val.val.c->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_LMP_BIGINT, MPI_SUM, world);
return ncountall;
}
} else if (val.which == ArgInfo::FIX) {
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
return group->count(igroup);
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
bigint ncount = val.val.f->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_LMP_BIGINT, MPI_SUM, world);

View File

@ -37,12 +37,11 @@ class ComputeReduce : public Compute {
double memory_usage() override;
protected:
int mode, nvalues;
int mode, nvalues, input_mode;
struct value_t {
int which;
int argindex;
std::string id;
int flavor;
union {
class Compute *c;
class Fix *f;

View File

@ -26,6 +26,8 @@
using namespace LAMMPS_NS;
enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduce
static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */
@ -97,7 +99,7 @@ double ComputeReduceRegion::compute_one(int m, int flag)
// invoke compute if not previously invoked
} else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
val.val.c->compute_peratom();
val.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
@ -122,7 +124,7 @@ double ComputeReduceRegion::compute_one(int m, int flag)
one = compute_array[flag][aidxm1];
}
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) {
val.val.c->compute_local();
val.val.c->invoked_flag |= Compute::INVOKED_LOCAL;
@ -151,7 +153,7 @@ double ComputeReduceRegion::compute_one(int m, int flag)
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style);
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
if (aidx == 0) {
double *fix_vector = val.val.f->vector_atom;
if (flag < 0) {
@ -171,7 +173,7 @@ double ComputeReduceRegion::compute_one(int m, int flag)
one = fix_array[flag][aidxm1];
}
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
if (aidx == 0) {
double *fix_vector = val.val.f->vector_local;
if (flag < 0)
@ -219,18 +221,18 @@ bigint ComputeReduceRegion::count(int m)
if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F)
return group->count(igroup, region);
else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
return group->count(igroup, region);
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
bigint ncount = val.val.c->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world);
return ncountall;
}
} else if (val.which == ArgInfo::FIX) {
if (val.flavor == PERATOM) {
if (input_mode == PERATOM) {
return group->count(igroup, region);
} else if (val.flavor == LOCAL) {
} else if (input_mode == LOCAL) {
bigint ncount = val.val.f->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world);