[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master e65eb48: Arithmetic: multithreaded variable-op
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master e65eb48: Arithmetic: multithreaded variable-operand operators |
|
Date: |
Tue, 19 Mar 2019 09:25:00 -0400 (EDT) |
branch: master
commit e65eb489065f030f90d557d7e6e1030a9a6c7cb5
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Arithmetic: multithreaded variable-operand operators
Arithmetic's variable-operand operators work independently on each pixel,
so as the datasets become larger and more numerous, they can greatly
benefit from multi-threaded analysis. However, until now, they done their
processing on a single thread.
With this commit, the `gal_arithmetic' library function accepts a new
argument for the number of threads to use and the variable-operand
operators will do their work in multi-threaded mode.
---
NEWS | 7 +
bin/arithmetic/arithmetic.c | 5 +-
bin/arithmetic/ui.c | 2 +
bin/convertt/convertt.c | 12 +-
bin/convertt/ui.c | 2 +-
bin/convolve/ui.c | 4 +-
bin/mkcatalog/ui.c | 2 +-
bin/statistics/ui.c | 8 +-
bin/table/table.c | 8 +-
doc/gnuastro.texi | 39 ++++--
lib/arithmetic.c | 326 ++++++++++++++++++++++++++------------------
lib/gnuastro/arithmetic.h | 2 +-
lib/options.c | 6 +-
13 files changed, 252 insertions(+), 171 deletions(-)
diff --git a/NEWS b/NEWS
index acf16eb..ced25a9 100644
--- a/NEWS
+++ b/NEWS
@@ -22,6 +22,10 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
`sigclip-median', and `sigclip-std'. These are very useful when
several inputs have strong outliers that affect the median, or the
mean is required.
+ - Multithreaded operation for the following operators that
+ combine/co-add multiple inputs into one output with same size: `min',
+ `max', `number', `sum', `mean', `std', `median', `sigclip-number',
+ `sigclip-median', `sigclip-mean', `sigclip-std'.
--wcsfile and --wcshdu: these two options can be used to specify a
different file for reading the WCS of the output. This is useful when
the default (theh WCS of the first dataset that is read) is not the
@@ -150,6 +154,9 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
version, the `-s' short option was used for it. But with the new
`--sort' option, `-s' may cause confusion.
+ Library
+ gal_arithmetic: new argument: number of threads to use (when relevant).
+
** Bugs fixed
bug #55313: Fits program writing --write values in reverse order
bug #55333: Fits program crash when --write or --update have no value.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index bfa33e0..0c60cdb 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -431,7 +431,7 @@ wrapper_for_filter(struct arithmeticparams *p, char *token,
int operator)
"along dimension %zu is a float", ndim-i);
/* Make sure it isn't negative. */
- comp=gal_arithmetic(GAL_ARITHMETIC_OP_GT, 0, tmp, zero);
+ comp=gal_arithmetic(GAL_ARITHMETIC_OP_GT, 1, 0, tmp, zero);
if( *(uint8_t *)(comp->array) == 0 )
error(EXIT_FAILURE, 0, "lengths of filter along dimensions "
"must be positive. The given length in dimension %zu"
@@ -1153,7 +1153,8 @@ reversepolish(struct arithmeticparams *p)
number of arguments it uses depend on the operator. So
when the operator doesn't need three operands, the extra
arguments will be ignored. */
- operands_add(p, NULL, gal_arithmetic(op, flags, d1, d2, d3));
+ operands_add(p, NULL, gal_arithmetic(op, p->cp.numthreads,
+ flags, d1, d2, d3));
}
/* No need to call the arithmetic library, call the proper
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index f08d5c1..0ef9a24 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/table.h>
#include <gnuastro/array.h>
+#include <gnuastro/threads.h>
#include <gnuastro-internal/timing.h>
#include <gnuastro-internal/options.h>
@@ -124,6 +125,7 @@ ui_initialize_options(struct arithmeticparams *p,
cp->program_exec = PROGRAM_EXEC;
cp->program_bibtex = PROGRAM_BIBTEX;
cp->program_authors = PROGRAM_AUTHORS;
+ cp->numthreads = gal_threads_number();
cp->coptions = gal_commonopts_options;
/* Modify the common options. */
diff --git a/bin/convertt/convertt.c b/bin/convertt/convertt.c
index 3ffb66d..128dde6 100644
--- a/bin/convertt/convertt.c
+++ b/bin/convertt/convertt.c
@@ -72,11 +72,11 @@ convertt_change(struct converttparams *p)
{
/* Make a condition array: all pixels with a value equal to
`change->from' will be set as 1 in this array. */
- cond=gal_arithmetic(GAL_ARITHMETIC_OP_EQ, GAL_ARITHMETIC_NUMOK,
+ cond=gal_arithmetic(GAL_ARITHMETIC_OP_EQ, 1, GAL_ARITHMETIC_NUMOK,
channel, change->from);
/* Now, use the condition array to set the proper values. */
- channel=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flags, channel,
+ channel=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flags, channel,
cond, change->to);
/* Clean up, since we set the free flag, all extra arrays have been
@@ -104,11 +104,11 @@ convertt_trunc_function(int operator, gal_data_t *data,
gal_data_t *value)
/* Make a condition array: all pixels with a value equal to
`change->from' will be set as 1 in this array. */
- cond=gal_arithmetic(operator, GAL_ARITHMETIC_NUMOK, data, value);
+ cond=gal_arithmetic(operator, 1, GAL_ARITHMETIC_NUMOK, data, value);
/* Now, use the condition array to set the proper values. */
- out=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flags, data, cond, value);
+ out=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flags, data, cond, value);
/* A small sanity check. The process must be in-place so the original
@@ -200,8 +200,8 @@ convertt_scale_to_uchar(struct converttparams *p)
}
/* Calculate the minimum and maximum. */
- mind = gal_arithmetic(GAL_ARITHMETIC_OP_MINVAL, 0, channel);
- maxd = gal_arithmetic(GAL_ARITHMETIC_OP_MAXVAL, 0, channel);
+ mind = gal_arithmetic(GAL_ARITHMETIC_OP_MINVAL, 1, 0, channel);
+ maxd = gal_arithmetic(GAL_ARITHMETIC_OP_MAXVAL, 1, 0, channel);
tmin = *((float *)(mind->array));
tmax = *((float *)(maxd->array));
gal_data_free(mind);
diff --git a/bin/convertt/ui.c b/bin/convertt/ui.c
index 71bcc0d..284801c 100644
--- a/bin/convertt/ui.c
+++ b/bin/convertt/ui.c
@@ -332,7 +332,7 @@ ui_read_check_only_options(struct converttparams *p)
if(p->fluxhighstr && p->fluxlowstr)
{
- cond=gal_arithmetic(GAL_ARITHMETIC_OP_GT, GAL_ARITHMETIC_NUMOK,
+ cond=gal_arithmetic(GAL_ARITHMETIC_OP_GT, 1, GAL_ARITHMETIC_NUMOK,
p->fluxhigh, p->fluxlow);
if( *((unsigned char *)cond->array) == 0 )
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 1e40688..3698f8d 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -559,11 +559,11 @@ ui_preparations(struct convolveparams *p)
meaningful. */
sum=gal_statistics_sum(p->input);
sum=gal_data_copy_to_new_type_free(sum, GAL_TYPE_FLOAT32);
- p->input = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
+ p->input = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE, 1,
GAL_ARITHMETIC_FLAGS_ALL, p->input, sum);
sum=gal_statistics_sum(p->kernel);
sum=gal_data_copy_to_new_type_free(sum, GAL_TYPE_FLOAT32);
- p->kernel = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
+ p->kernel = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE, 1,
GAL_ARITHMETIC_FLAGS_ALL, p->kernel, sum);
}
}
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 5a9ee5c..cc2daf4 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -1042,7 +1042,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
pixels and 0 for zero pixels. */
zero=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &one, NULL, 1, -1,
NULL, NULL, NULL);
- p->upmask=gal_arithmetic(GAL_ARITHMETIC_OP_NE,
+ p->upmask=gal_arithmetic(GAL_ARITHMETIC_OP_NE, 1,
( GAL_ARITHMETIC_INPLACE
| GAL_ARITHMETIC_FREE
| GAL_ARITHMETIC_NUMOK ),
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 4e1c7d4..662c0df 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -597,7 +597,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
NULL, NULL, NULL);
*((float *)(tmp->array)) = p->greaterequal;
- cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_LT, flags, ref, tmp);
+ cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_LT, 1, flags, ref, tmp);
gal_data_free(tmp);
}
@@ -608,7 +608,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
NULL, NULL, NULL);
*((float *)(tmp->array)) = p->lessthan;
- cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_GE, flags, ref, tmp);
+ cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_GE, 1, flags, ref, tmp);
gal_data_free(tmp);
}
@@ -622,7 +622,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
cond = isnan(p->greaterequal) ? cond_l : cond_g;
break;
case 2:
- cond = gal_arithmetic(GAL_ARITHMETIC_OP_OR, flagsor, cond_l, cond_g);
+ cond = gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_l, cond_g);
break;
}
@@ -637,7 +637,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
/* Set all the pixels that satisfy the condition to blank. Note that a
blank value will be used in the proper type of the input in the
`where' operator.*/
- gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flagsor, p->input, cond, blank);
+ gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flagsor, p->input, cond, blank);
/* Reset the blank flags so they are checked again if necessary. */
diff --git a/bin/table/table.c b/bin/table/table.c
index 1b5008c..7577064 100644
--- a/bin/table/table.c
+++ b/bin/table/table.c
@@ -109,12 +109,12 @@ table_range(struct tableparams *p)
/* Find all the bad elements (smaller than the minimum and larger than
the maximum) so we can flag them. */
- ltmin=gal_arithmetic(GAL_ARITHMETIC_OP_LT, numok, ref, min);
- gemax=gal_arithmetic(GAL_ARITHMETIC_OP_GE, numok, ref, max);
- ltmin=gal_arithmetic(GAL_ARITHMETIC_OP_OR, inplace, ltmin, gemax);
+ ltmin=gal_arithmetic(GAL_ARITHMETIC_OP_LT, 1, numok, ref, min);
+ gemax=gal_arithmetic(GAL_ARITHMETIC_OP_GE, 1, numok, ref, max);
+ ltmin=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, inplace, ltmin, gemax);
/* Add these flags to all previous flags. */
- mask=gal_arithmetic(GAL_ARITHMETIC_OP_OR, inplace, mask, ltmin);
+ mask=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, inplace, mask, ltmin);
/* For a check.
{
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6e94ddb..79664bc 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12270,13 +12270,17 @@ to @command{minvalue}.
@item min
The first popped operand to this operator must be a positive integer number
which specifies how many further operands should be popped from the
-stack. The given number of operands must have the same type and size. Each
-pixel of the output of this operator will be set to the minimum value of
-the given number of operands (images) in that pixel.
+stack. All the subsequently popped operands must have the same type and
+size. This operator (and all the variable-operand operators similar to it
+that are discussed below) will work in multithreaded mode unless Arithmetic
+is called with the @option{--numthreads=1} option, see @ref{Multi-threaded
+operations}.
-For example the following command will produce an image with the same size
-and type as the inputs but each output pixel is set to the minimum
-respective pixel value of the three input images.
+Each pixel of the output of the @code{min} operator will be given the
+minimum value of the same pixel from all the popped operands/images. For
+example the following command will produce an image with the same size and
+type as the three inputs, but each output pixel value will be the minimum
+of the same pixel's values in all three input images.
@example
$ astarithmetic a.fits b.fits c.fits 3 min
@@ -27144,8 +27148,9 @@ Multi-operand statistical operations. When
@code{gal_arithmetic} is called
with any of these operators, it will expect only a single operand that will
be interpreted as a list of datasets (see @ref{List of gal_data_t}. The
output will be a single dataset with each of its elements replaced by the
-respective statistical operation on the whole list. See the discussion
-under the @code{min} operator in @ref{Arithmetic operators}.
+respective statistical operation on the whole list. These operators can
+work on multiple threads using the @code{numthreads} argument. See the
+discussion under the @code{min} operator in @ref{Arithmetic operators}.
@end deffn
@deffn Macro GAL_ARITHMETIC_OP_SIGCLIP_STD
@@ -27216,11 +27221,15 @@ directly use those functions.
@end deffn
address@hidden {gal_data_t *} gal_arithmetic (int operator, int flags, ...)
address@hidden {gal_data_t *} gal_arithmetic (int @code{operator}, size_t
@code{numthreads}, int @code{flags}, ...)
Do the arithmetic operation of @code{operator} on the given operands (the
-third argument and any further argument). Certain special conditions can
-also be specified with the @code{flag} operator. The acceptable values for
address@hidden are defined in the macros above.
+third argument and any further argument). If the operator can work on
+multiple threads, the number of threads can be specified with
address@hidden When the operator is single-threaded, @code{numthreads}
+will be ignored. Special conditions can also be specified with the
address@hidden operator (a bit-flag with bits described above, for example
address@hidden or @code{GAL_ARITHMETIC_FREE}). The
+acceptable values for @code{operator} are also defined in the macros above.
@code{gal_arithmetic} is a multi-argument function (like C's
@code{printf}). In other words, the number of necessary arguments is not
@@ -27228,9 +27237,9 @@ fixed and depends on the value to @code{operator}. Here
are a few examples
showing this variability:
@example
-out_1=gal_arithmetic(GAL_ARITHMETIC_OP_LOG, 0, in_1);
-out_2=gal_arithmetic(GAL_ARITHMETIC_OP_PLUS, 0, in_1, in_2);
-out_3=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 0, in_1, in_2, in_3);
+out_1=gal_arithmetic(GAL_ARITHMETIC_OP_LOG, 1, 0, in_1);
+out_2=gal_arithmetic(GAL_ARITHMETIC_OP_PLUS, 1, 0, in_1, in_2);
+out_3=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, 0, in_1, in_2, in_3);
@end example
The number of necessary operands for each operator (and thus the number of
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index aba3df2..dfdfabb 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/blank.h>
#include <gnuastro/qsort.h>
#include <gnuastro/pointer.h>
+#include <gnuastro/threads.h>
#include <gnuastro/dimension.h>
#include <gnuastro/statistics.h>
#include <gnuastro/arithmetic.h>
@@ -669,27 +670,45 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
/***********************************************************************/
/*************** Multiple operand operators **************/
/***********************************************************************/
+struct multioperandparams
+{
+ gal_data_t *list; /* List of input datasets. */
+ gal_data_t *out; /* Output dataset. */
+ size_t dnum; /* Number of input dataset. */
+ int operator; /* Operator to use. */
+ uint8_t *hasblank; /* Array of 0s or 1s for each input. */
+ float p1; /* Sigma-cliping parameter 1. */
+ float p2; /* Sigma-cliping parameter 2. */
+};
+
+
+
+
+
#define MULTIOPERAND_MIN(TYPE) { \
- TYPE p, max; \
+ TYPE t, max; \
size_t n, j=0; \
- gal_type_max(list->type, &max); \
- do /* Loop over each pixel */ \
+ gal_type_max(p->list->type, &max); \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
- p=max; \
- for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ t=max; \
+ j=tprm->indexs[tind]; \
+ \
+ for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
{ /* Only for integer types, b==b. */ \
- if( hasblank[i] && b==b) \
+ if( p->hasblank[i] && b==b) \
{ \
if( a[i][j] != b ) \
- { p = a[i][j] < p ? a[i][j] : p; ++n; } \
+ { t = a[i][j] < t ? a[i][j] : t; ++n; } \
} \
- else { p = a[i][j] < p ? a[i][j] : p; ++n; } \
+ else { t = a[i][j] < t ? a[i][j] : t; ++n; } \
} \
- *o++ = n ? p : b; /* No usable elements: set to blank. */ \
- ++j; \
+ o[j] = n ? t : b; /* No usable elements: set to blank. */ \
} \
- while(o<of); \
}
@@ -697,26 +716,29 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_MAX(TYPE) { \
- TYPE p, min; \
+ TYPE t, min; \
size_t n, j=0; \
- gal_type_min(list->type, &min); \
- do /* Loop over each pixel */ \
+ gal_type_min(p->list->type, &min); \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
- p=min; \
- for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ t=min; \
+ j=tprm->indexs[tind]; \
+ \
+ for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
{ /* Only for integer types, b==b. */ \
- if( hasblank[i] && b==b) \
+ if( p->hasblank[i] && b==b) \
{ \
if( a[i][j] != b ) \
- { p = a[i][j] > p ? a[i][j] : p; ++n; } \
+ { t = a[i][j] > t ? a[i][j] : t; ++n; } \
} \
- else { p = a[i][j] > p ? a[i][j] : p; ++n; } \
+ else { t = a[i][j] > t ? a[i][j] : t; ++n; } \
} \
- *o++ = n ? p : b; /* No usable elements: set to blank. */ \
- ++j; \
+ o[j] = n ? t : b; /* No usable elements: set to blank. */ \
} \
- while(o<of); \
}
@@ -725,14 +747,19 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_NUM { \
int use; \
- size_t n, j=0; \
- do /* Loop over each pixel */ \
+ size_t n, j; \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
- for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ j=tprm->indexs[tind]; \
+ \
+ for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
{ \
/* Only integers and non-NaN floats: v==v is 1. */ \
- if(hasblank[i]) \
+ if(p->hasblank[i]) \
use = ( b==b \
? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
: ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
@@ -741,10 +768,8 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
/* Increment counter if necessary. */ \
if(use) ++n; \
} \
- *o++ = n; \
- ++j; \
+ o[j] = n; \
} \
- while(o<of); \
}
@@ -754,15 +779,20 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_SUM { \
int use; \
double sum; \
- size_t n, j=0; \
- do /* Loop over each pixel */ \
+ size_t n, j; \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
sum=0.0f; \
- for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ j=tprm->indexs[tind]; \
+ \
+ for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
{ \
/* Only integers and non-NaN floats: v==v is 1. */ \
- if(hasblank[i]) \
+ if(p->hasblank[i]) \
use = ( b==b \
? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
: ( a[i][j]==*a[i] ? 1 : 0 ) ); /* Float */ \
@@ -771,10 +801,8 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
/* Use in sum if necessary. */ \
if(use) { sum += a[i][j]; ++n; } \
} \
- *o++ = n ? sum : b; \
- ++j; \
+ o[j] = n ? sum : b; \
} \
- while(o<of); \
}
@@ -784,15 +812,20 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_MEAN { \
int use; \
double sum; \
- size_t n, j=0; \
- do /* Loop over each pixel */ \
+ size_t n, j; \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
sum=0.0f; \
- for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ j=tprm->indexs[tind]; \
+ \
+ for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
{ \
/* Only integers and non-NaN floats: v==v is 1. */ \
- if(hasblank[i]) \
+ if(p->hasblank[i]) \
use = ( b==b \
? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
: ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
@@ -801,10 +834,8 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
/* Calculate the mean if necessary. */ \
if(use) { sum += a[i][j]; ++n; } \
} \
- *o++ = n ? sum/n : b; \
- ++j; \
+ o[j] = n ? sum/n : b; \
} \
- while(o<of); \
}
@@ -813,16 +844,21 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_STD { \
int use; \
- size_t n, j=0; \
+ size_t n, j; \
double sum, sum2; \
- do /* Loop over each pixel */ \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
sum=sum2=0.0f; \
- for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ j=tprm->indexs[tind]; \
+ \
+ for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
{ \
/* Only integers and non-NaN floats: v==v is 1. */ \
- if(hasblank[i]) \
+ if(p->hasblank[i]) \
use = ( b==b \
? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
: ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
@@ -836,10 +872,8 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
++n; \
} \
} \
- *o++ = n ? sqrt( (sum2-sum*sum/n)/n ) : b; \
- ++j; \
+ o[j] = n ? sqrt( (sum2-sum*sum/n)/n ) : b; \
} \
- while(o<of); \
}
@@ -848,21 +882,22 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) { \
int use; \
- size_t n, j=0; \
- TYPE *pixs=gal_pointer_allocate(list->type, dnum, 0, __func__, \
- "pixs"); \
+ size_t n, j; \
+ TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
+ __func__, "pixs"); \
\
- /* Loop over each pixel */ \
- do \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
- /* Initialize. */ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
+ j=tprm->indexs[tind]; \
\
- /* Loop over each array. */ \
- for(i=0;i<dnum;++i) \
+ /* Loop over each array: `i' is input dataset's index. */ \
+ for(i=0;i<p->dnum;++i) \
{ \
/* Only integers and non-NaN floats: v==v is 1. */ \
- if(hasblank[i]) \
+ if(p->hasblank[i]) \
use = ( b==b \
? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
: ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
@@ -876,13 +911,11 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
if(n) \
{ \
qsort(pixs, n, sizeof *pixs, QSORT_F); \
- *o++ = n%2 ? pixs[n/2] : (pixs[n/2] + pixs[n/2-1])/2 ; \
+ o[j] = n%2 ? pixs[n/2] : (pixs[n/2] + pixs[n/2-1])/2 ; \
} \
else \
- *o++=b; \
- ++j; \
+ o[j]=b; \
} \
- while(o<of); \
\
/* Clean up. */ \
free(pixs); \
@@ -894,49 +927,48 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_SIGCLIP(TYPE) { \
float *sarr; \
- size_t n, j=0; \
+ size_t n, j; \
gal_data_t *sclip; \
- TYPE *pixs=gal_pointer_allocate(list->type, dnum, 0, __func__, \
- "pixs"); \
- gal_data_t *cont=gal_data_alloc(pixs, list->type, 1, &dnum, NULL, \
- 0, -1, NULL, NULL, NULL); \
+ TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
+ __func__, "pixs"); \
+ gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum, \
+ NULL, 0, -1, NULL, NULL, NULL); \
\
- /* Loop over each pixel */ \
- do \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
{ \
- /* Initialize. */ \
+ /* Initialize, `j' is desired pixel's index. */ \
n=0; \
+ j=tprm->indexs[tind]; \
\
- /* Loop over each array. */ \
- for(i=0;i<dnum;++i) pixs[n++]=a[i][j]; \
+ /* Read the necessay values from each input. */ \
+ for(i=0;i<p->dnum;++i) pixs[n++]=a[i][j]; \
\
/* If there are any elements, measure the */ \
if(n) \
{ \
- sclip=gal_statistics_sigma_clip(cont, p1, p2, 1, 1); \
+ sclip=gal_statistics_sigma_clip(cont, p->p1, p->p2, 1, 1); \
sarr=sclip->array; \
- switch(operator) \
+ switch(p->operator) \
{ \
- case GAL_ARITHMETIC_OP_SIGCLIP_STD: *o++=sarr[3]; break;\
- case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: *o++=sarr[2]; break;\
- case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: *o++=sarr[1]; break;\
- case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: *o++=sarr[0]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_STD: o[j]=sarr[3]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: o[j]=sarr[2]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: o[j]=sarr[1]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: o[j]=sarr[0]; break;\
default: \
error(EXIT_FAILURE, 0, "%s: a bug! the code %d is not " \
"valid for sigma-clipping results", __func__, \
- operator); \
+ p->operator); \
} \
\
/* Since we are doing sigma-clipping in place, the size, */ \
/* and flags need to be reset. */ \
cont->flag=0; \
- cont->size=cont->dsize[0]=dnum; \
+ cont->size=cont->dsize[0]=p->dnum; \
} \
else \
- *o++=b; \
- ++j; \
+ o[j]=b; \
} \
- while(o<of); \
\
/* Clean up. */ \
gal_data_free(cont); \
@@ -947,25 +979,26 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
#define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) { \
- TYPE b, **a, *o=out->array, *of=o+out->size; \
- size_t i=0; /* Different from the `i' in the main function. */ \
+ gal_data_t *tmp; \
+ size_t i=0, tind; \
+ TYPE b, **a, *o=p->out->array; \
\
/* Allocate space to keep the pointers to the arrays of each. */ \
/* Input data structure. The operators will increment these */ \
/* pointers while parsing them. */ \
errno=0; \
- a=malloc(dnum*sizeof *a); \
+ a=malloc(p->dnum*sizeof *a); \
if(a==NULL) \
error(EXIT_FAILURE, 0, "%s: %zu bytes for `a'", \
- "MULTIOPERAND_TYPE_SET", dnum*sizeof *a); \
+ "MULTIOPERAND_TYPE_SET", p->dnum*sizeof *a); \
\
/* Fill in the array pointers and the blank value for this type. */ \
- gal_blank_write(&b, list->type); \
- for(tmp=list;tmp!=NULL;tmp=tmp->next) \
+ gal_blank_write(&b, p->list->type); \
+ for(tmp=p->list;tmp!=NULL;tmp=tmp->next) \
a[i++]=tmp->array; \
\
/* Do the operation. */ \
- switch(operator) \
+ switch(p->operator) \
{ \
case GAL_ARITHMETIC_OP_MIN: \
MULTIOPERAND_MIN(TYPE); \
@@ -1004,7 +1037,7 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
\
default: \
error(EXIT_FAILURE, 0, "%s: operator code %d not recognized", \
- "MULTIOPERAND_TYPE_SET", operator); \
+ "MULTIOPERAND_TYPE_SET", p->operator); \
} \
\
/* Clean up. */ \
@@ -1015,16 +1048,72 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
+/* Worker function on each thread. */
+void *
+multioperand_on_thread(void *in_prm)
+{
+ /* Low-level definitions to be done first. */
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ struct multioperandparams *p=(struct multioperandparams *)tprm->params;
+
+ /* Do the operation on each thread. */
+ switch(p->list->type)
+ {
+ case GAL_TYPE_UINT8:
+ MULTIOPERAND_TYPE_SET(uint8_t, gal_qsort_uint8_i);
+ break;
+ case GAL_TYPE_INT8:
+ MULTIOPERAND_TYPE_SET(int8_t, gal_qsort_int8_i);
+ break;
+ case GAL_TYPE_UINT16:
+ MULTIOPERAND_TYPE_SET(uint16_t, gal_qsort_uint16_i);
+ break;
+ case GAL_TYPE_INT16:
+ MULTIOPERAND_TYPE_SET(int16_t, gal_qsort_int16_i);
+ break;
+ case GAL_TYPE_UINT32:
+ MULTIOPERAND_TYPE_SET(uint32_t, gal_qsort_uint32_i);
+ break;
+ case GAL_TYPE_INT32:
+ MULTIOPERAND_TYPE_SET(int32_t, gal_qsort_int32_i);
+ break;
+ case GAL_TYPE_UINT64:
+ MULTIOPERAND_TYPE_SET(uint64_t, gal_qsort_uint64_i);
+ break;
+ case GAL_TYPE_INT64:
+ MULTIOPERAND_TYPE_SET(int64_t, gal_qsort_int64_i);
+ break;
+ case GAL_TYPE_FLOAT32:
+ MULTIOPERAND_TYPE_SET(float, gal_qsort_float32_i);
+ break;
+ case GAL_TYPE_FLOAT64:
+ MULTIOPERAND_TYPE_SET(double, gal_qsort_float64_i);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, p->list->type);
+ }
+
+ /* Wait for all the other threads to finish, then return. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
/* The single operator in this function is assumed to be a linked list. The
number of operators is determined from the fact that the last node in
the linked list must have a NULL pointer as its `next' element. */
static gal_data_t *
arithmetic_multioperand(int operator, int flags, gal_data_t *list,
- gal_data_t *params)
+ gal_data_t *params, size_t numthreads)
{
uint8_t *hasblank;
size_t i=0, dnum=1;
float p1=NAN, p2=NAN;
+ struct multioperandparams p;
gal_data_t *out, *tmp, *ttmp;
@@ -1087,43 +1176,16 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
hasblank[i++]=gal_blank_present(tmp, 0);
- /* Start the operation. */
- switch(list->type)
- {
- case GAL_TYPE_UINT8:
- MULTIOPERAND_TYPE_SET(uint8_t, gal_qsort_uint8_i);
- break;
- case GAL_TYPE_INT8:
- MULTIOPERAND_TYPE_SET(int8_t, gal_qsort_int8_i);
- break;
- case GAL_TYPE_UINT16:
- MULTIOPERAND_TYPE_SET(uint16_t, gal_qsort_uint16_i);
- break;
- case GAL_TYPE_INT16:
- MULTIOPERAND_TYPE_SET(int16_t, gal_qsort_int16_i);
- break;
- case GAL_TYPE_UINT32:
- MULTIOPERAND_TYPE_SET(uint32_t, gal_qsort_uint32_i);
- break;
- case GAL_TYPE_INT32:
- MULTIOPERAND_TYPE_SET(int32_t, gal_qsort_int32_i);
- break;
- case GAL_TYPE_UINT64:
- MULTIOPERAND_TYPE_SET(uint64_t, gal_qsort_uint64_i);
- break;
- case GAL_TYPE_INT64:
- MULTIOPERAND_TYPE_SET(int64_t, gal_qsort_int64_i);
- break;
- case GAL_TYPE_FLOAT32:
- MULTIOPERAND_TYPE_SET(float, gal_qsort_float32_i);
- break;
- case GAL_TYPE_FLOAT64:
- MULTIOPERAND_TYPE_SET(double, gal_qsort_float64_i);
- break;
- default:
- error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
- __func__, list->type);
- }
+ /* Set the parameters necessary for multithreaded operation and spin them
+ off to do apply the operator. */
+ p.p1=p1;
+ p.p2=p2;
+ p.out=out;
+ p.list=list;
+ p.dnum=dnum;
+ p.operator=operator;
+ p.hasblank=hasblank;
+ gal_threads_spin_off(multioperand_on_thread, &p, out->size, numthreads);
/* Clean up and return. Note that the operation might have been done in
@@ -1551,7 +1613,7 @@ gal_arithmetic_operator_string(int operator)
gal_data_t *
-gal_arithmetic(int operator, int flags, ...)
+gal_arithmetic(int operator, size_t numthreads, int flags, ...)
{
va_list va;
gal_data_t *d1, *d2, *d3, *out=NULL;
@@ -1641,7 +1703,7 @@ gal_arithmetic(int operator, int flags, ...)
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER:
d1 = va_arg(va, gal_data_t *);
d2 = va_arg(va, gal_data_t *);
- out=arithmetic_multioperand(operator, flags, d1, d2);
+ out=arithmetic_multioperand(operator, flags, d1, d2, numthreads);
break;
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index db65bc8..a1922c8 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -143,7 +143,7 @@ enum gal_arithmetic_operators
gal_data_t *
-gal_arithmetic(int operator, int flags, ...);
+gal_arithmetic(int operator, size_t numthreads, int flags, ...);
diff --git a/lib/options.c b/lib/options.c
index 8b34aea..10ae646 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -1337,11 +1337,11 @@ options_sanity_check(struct argp_option *option, char
*arg,
`GAL_ARITHMETIC_INPLACE' flags. But we will do this when there are
multiple checks so from the two check data structures, we only have
one remaining. */
- check1=gal_arithmetic(operator1, GAL_ARITHMETIC_NUMOK, value, ref1);
+ check1=gal_arithmetic(operator1, 1, GAL_ARITHMETIC_NUMOK, value, ref1);
if(ref2)
{
- check2=gal_arithmetic(operator2, GAL_ARITHMETIC_NUMOK, value, ref2);
- check1=gal_arithmetic(multicheckop, mcflag, check1, check2);
+ check2=gal_arithmetic(operator2, 1, GAL_ARITHMETIC_NUMOK, value, ref2);
+ check1=gal_arithmetic(multicheckop, 1, mcflag, check1, check2);
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master e65eb48: Arithmetic: multithreaded variable-operand operators,
Mohammad Akhlaghi <=