diff --git a/pygrt/C_extension/include/grt/common/matrix.h b/pygrt/C_extension/include/grt/common/matrix.h index abc11b2..2aa7282 100755 --- a/pygrt/C_extension/include/grt/common/matrix.h +++ b/pygrt/C_extension/include/grt/common/matrix.h @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-07-24 * - * 2x2小矩阵的加、减、乘、除、求逆等操作，由于均为小型数组操作，所有函数均为内联函数 + * 小矩阵的加、减、乘、除、求逆等操作，由于均为小型数组操作，所有函数均为内联函数 */ #pragma once @@ -142,7 +142,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOM } /** - * 计算nxn复矩阵的积(小矩阵)(最暴力的方式) + * 计算mxn复矩阵的积(小矩阵)(最暴力的方式) * * @param[in] m1 M1矩阵行数 * @param[in] n1 M1矩阵列数 @@ -171,6 +171,31 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, cons } } +/** + * 计算mxn复矩阵的转置矩阵(不共轭) + * + * @param[in] m1 M1矩阵行数 + * @param[in] n1 M1矩阵列数 + * @param[in] M1 M1矩阵 + * @param[out] M2 M2矩阵 (M1^T) + */ +inline GCC_ALWAYS_INLINE void grt_cmatmxn_transpose(MYINT m1, MYINT n1, const MYCOMPLEX M1[m1][n1], MYCOMPLEX M2[n1][m1]){ + MYINT m, n; + MYCOMPLEX M0[n1][m1]; + for(m=0; m #include "grt/common/const.h" +// 定义 X 宏，为多个类型定义查找函数 +#define __FOR_EACH_TYPE \ + X(MYINT) X(MYREAL) X(float) X(double) + + /** - * 该函数对输入的整数数组进行线性搜索，找到目标值时返回其索引。 + * 该函数对输入数组进行线性搜索，找到目标值时返回其索引。 * 如果目标值在数组中未找到，则返回 -1。 * - * @param[in] array 要搜索的整数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小（元素个数）。 * @param[in] target 要查找的目标值。 * @@ -24,13 +29,19 @@ * @note 如果数组中存在多个目标值，该函数返回第一个匹配的索引。 * */ -MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target); +#define X(T) \ +MYINT grt_findElement_##T(const T *array, MYINT size, T target); + +__FOR_EACH_TYPE +#undef X + + /** - * 搜索浮点数数组中最接近目标值且小于目标值的索引。 + * 搜索数组中最接近目标值且小于目标值的索引。 * 如果目标值在数组中未找到，则返回 -1。 * - * @param[in] array 要搜索的浮点数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小（元素个数）。 * @param[in] target 要查找的目标值。 * @@ -39,12 +50,18 @@ MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target); * @note 如果数组中存在多个目标值，该函数返回第一个匹配的索引。 * */ -MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); +#define X(T) \ +MYINT grt_findLessEqualClosest_##T(const T *array, MYINT size, T target); + +__FOR_EACH_TYPE +#undef X + + /** - * 搜索浮点数数组中最接近目标值的索引。 + * 搜索数组中最接近目标值的索引。 * - * @param[in] array 要搜索的浮点数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小（元素个数）。 * @param[in] target 要查找的目标值。 * @@ -53,18 +70,62 @@ MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL t * @note 如果数组中存在多个目标值，该函数返回第一个匹配的索引。 * */ -MYINT grt_findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); +#define X(T) \ +MYINT grt_findClosest_##T(const T *array, MYINT size, T target); + +__FOR_EACH_TYPE +#undef X + /** - * 搜索浮点数数组的最大或最小值，返回其索引。 + * 搜索数组的最值，返回其索引。 * - * @param[in] array 要搜索的浮点数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小（元素个数）。 - * @param[in] isMax 是否要找最大值，否则找最小值。 * * @return idx 目标值的索引。 * * @note 如果数组中存在相同最值，该函数返回第一个匹配的索引。 * */ -MYINT grt_findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax); \ No newline at end of file +#define X(T) \ +MYINT grt_findMin_##T(const T *array, MYINT size);\ +MYINT grt_findMax_##T(const T *array, MYINT size);\ + +__FOR_EACH_TYPE +#undef X + + +/** + * 比较函数 + * + * @param[in] a 元素 a 地址 + * @param[in] b 元素 b 地址 + * + * @return flag 比较结果，(1) a > b, (0) a == b, (-1) a < b + */ +#define X(T) \ +int grt_compare_##T(const void *a, const void *b); + +__FOR_EACH_TYPE +#undef X +#undef __FOR_EACH_TYPE + + + +/** + * 在有序数组中插入元素，元素类型和数组类型需匹配 + * + * @param[in,out] arr 有序数组地址 + * @param[in,out] size 数组大小地址 + * @param[in] capacity 数组最大容量 + * @param[in] target 元素地址 + * @param[in] elementSize 元素和数组内元素的字节长度 + * @param[in] ascending 升序(true) 或 降序(false) + * @param[in] compare 比较函数 + * + * @return pos 插入位置的索引 + */ +MYINT grt_insertOrdered( + void *arr, MYINT *size, MYINT capacity, const void *target, size_t elementSize, bool ascending, + int (*compare)(const void *, const void *)); diff --git a/pygrt/C_extension/src/common/search.c b/pygrt/C_extension/src/common/search.c index 21b6a7d..dfa2047 100644 --- a/pygrt/C_extension/src/common/search.c +++ b/pygrt/C_extension/src/common/search.c @@ -8,59 +8,150 @@ #include #include +#include #include "grt/common/search.h" #include "grt/common/const.h" -static bool _gt_(MYREAL a1, MYREAL a2) { return a1 > a2; } -static bool _lt_(MYREAL a1, MYREAL a2) { return a1 < a2; } +// 定义 X 宏，为多个类型定义查找函数 +#define __FOR_EACH_TYPE \ + X(MYINT) X(MYREAL) X(float) X(double) -MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target) { - for (MYINT i = 0; i < size; ++i) { - if (array[i] == target) { - return i; // 找到目标元素，返回索引 - } - } - return -1; // 未找到目标元素，返回-1 +#define X(T) \ +MYINT grt_findElement_##T(const T *array, MYINT size, T target) {\ + for (MYINT i = 0; i < size; ++i) {\ + if (array[i] == target) {\ + /** 找到目标元素，返回索引 */\ + return i;\ + }\ + }\ + /** 未找到目标元素，返回-1 */\ + return -1; \ } -MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { - MYINT ires=-1; - MYREAL mindist=-1.0, dist=0.0; - for (MYINT i = 0; i < size; ++i) { - dist = target-array[i]; - if(dist >= 0.0 && (mindist < 0.0 || dist < mindist)){ - ires = i; - mindist = dist; - } - } - return ires; +__FOR_EACH_TYPE +#undef X + + + +#define X(T) \ +MYINT grt_findLessEqualClosest_##T(const T *array, MYINT size, T target) {\ + MYINT ires=-1;\ + T mindist=-1, dist=0;\ + for (MYINT i = 0; i < size; ++i) {\ + dist = target-array[i];\ + if(dist >= 0 && (mindist < 0 || dist < mindist)){\ + ires = i;\ + mindist = dist;\ + }\ + }\ + return ires;\ } -MYINT grt_findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { - MYINT ires=0; - MYREAL mindist=-1.0, dist=0.0; - for (MYINT i = 0; i < size; ++i) { - dist = fabs(target-array[i]); - if(mindist < 0.0 || dist < mindist){ - ires = i; - mindist = dist; - } - } - return ires; +__FOR_EACH_TYPE +#undef X + + + +#define X(T) \ +MYINT grt_findClosest_##T(const T *array, MYINT size, T target) {\ + MYINT ires=0;\ + T mindist=-1, dist=0;\ + for (MYINT i = 0; i < size; ++i) {\ + dist = fabs(target-array[i]);\ + if(mindist < 0 || dist < mindist){\ + ires = i;\ + mindist = dist;\ + }\ + }\ + return ires;\ +} + +__FOR_EACH_TYPE +#undef X + + +#define X(T) \ +MYINT grt_findMin_##T(const T *array, MYINT size) {\ + T rval = array[0];\ + MYINT idx=0;\ + for(MYINT ir=0; ir rval){\ + rval = array[ir];\ + idx = ir;\ + }\ + }\ + return idx;\ } -MYINT grt_findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax) { - MYREAL rmax = array[0]; - MYINT idx=0; - bool (*_func)(MYREAL, MYREAL) = (isMax)? _gt_ : _lt_; - for(MYINT ir=0; ir valb){\ + return 1;\ + } else if (vala < valb){\ + return -1;\ + } else {\ + return 0;\ + }\ +}\ + +__FOR_EACH_TYPE +#undef X +#undef __FOR_EACH_TYPE + + +MYINT grt_insertOrdered( + void *arr, MYINT *size, MYINT capacity, const void *target, size_t elementSize, bool ascending, + int (*compare)(const void *, const void *)) +{ + MYINT sgn = (ascending)? 1 : -1; + + // 数组满载情况下，只可能插入更小(升序)或更大(降序)的数值 + if(*size == capacity && sgn*compare(target, arr+(*size-1)*elementSize) >= 0) return -1; + + // 找到插入位置 + MYINT pos=*size; + for(MYINT i=0; i<*size; ++i){ + if(sgn*compare(target, arr+i*elementSize) < 0){ + pos = i; + break; } } - return idx; -} + // 截断式插入，防止越界 + MYINT lastpos = *size; + if(lastpos >= capacity){ + lastpos = capacity-1; + } else { + ++(*size); + } + pos = GRT_MIN(pos, lastpos); + + // 移动插入位置后的元素 + memmove(arr + (pos + 1) * elementSize, + arr + pos * elementSize, + (lastpos - pos) * elementSize); + + // 插入新元素 + memcpy(arr + pos * elementSize, target, elementSize); + + return pos; +} \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index 303fe50..6f4a4a3 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -92,7 +92,7 @@ void grt_integ_grn_spec( gettimeofday(&begin_t, NULL); // 最大震中距 - MYINT irmax = grt_findMinMax_MYREAL(rs, nr, true); + MYINT irmax = grt_findMax_MYREAL(rs, nr); MYREAL rmax=rs[irmax]; const MYREAL Rho = mod1d->Rho[mod1d->isrc]; // 震源区密度 diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 7fcef70..03a7890 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -733,7 +733,7 @@ int greenfn_main(int argc, char **argv) { Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt; // 最大震中距 - MYREAL rmax = Ctrl->R.rs[grt_findMinMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr, true)]; + MYREAL rmax = Ctrl->R.rs[grt_findMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr)]; // 时窗最大截止时刻 MYREAL tmax = Ctrl->E.delayT0 + Ctrl->N.winT; diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index ca7f2f1..4be6c32 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -76,7 +76,7 @@ void grt_integ_static_grn( const char *statsstr // 积分结果输出 ){ - MYREAL rmax=rs[grt_findMinMax_MYREAL(rs, nr, true)]; // 最大震中距 + MYREAL rmax=rs[grt_findMax_MYREAL(rs, nr)]; // 最大震中距 const MYREAL hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0)