From 3d0f09b39a8c2ee4aa0d636e28f6b42b1c7edd8d Mon Sep 17 00:00:00 2001
From: de Lavenne Alban <alban.de-lavenne@irstea.fr>
Date: Thu, 20 Mar 2025 11:14:28 +0100
Subject: [PATCH] perf: fortran subroutine for convolution
---
DESCRIPTION | 4 ++--
R/convolution.R | 3 +--
src/convolution.f90 | 23 +++++++++++++++++++++++
src/similarity_module.f90 | 2 +-
src/transfR.c | 12 +++++++++++-
5 files changed, 38 insertions(+), 6 deletions(-)
create mode 100644 src/convolution.f90
diff --git a/DESCRIPTION b/DESCRIPTION
index 65bac29..07572bf 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: transfR
Type: Package
Title: Transfer of Hydrograph from Gauged to Ungauged Catchments
-Version: 1.1.3
-Date: 2025-03-19
+Version: 1.1.4
+Date: 2025-03-20
Authors@R: c(
person("Alban", "de Lavenne", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9448-3490"), email = "alban.delavenne@inrae.fr"),
person("Christophe", "Cudennec", role = c("ths"), comment = c(ORCID = "0000-0002-1707-8926"), email = "christophe.cudennec@agrocampus-ouest.fr"),
diff --git a/R/convolution.R b/R/convolution.R
index b69ba4f..c169ecf 100644
--- a/R/convolution.R
+++ b/R/convolution.R
@@ -31,8 +31,7 @@ convolution = function(Rn, ...) UseMethod("convolution")
#' @export
convolution.default <- function(Rn, uh, continuous=FALSE, ...){
Q <- rep(0,length(Rn)+length(uh))
- for(t in which(Rn>0)) Q[t:(t+length(uh)-1)] <- Rn[t]*uh + Q[t:(t+length(uh)-1)]
- for(t in which(is.na(Rn))) Q[t:(t+length(uh)-1)] <- NA
+ Q <- .Call("c_convolution", rn=Rn, uh=uh)
if(continuous){
Q[1:(length(uh)-1)] <- NA
return(Q[1:length(Rn)])
diff --git a/src/convolution.f90 b/src/convolution.f90
new file mode 100644
index 0000000..c4a85ba
--- /dev/null
+++ b/src/convolution.f90
@@ -0,0 +1,23 @@
+subroutine convolution(rn, uh, nrn, nuh, q) bind(C, name="convolution")
+ use iso_c_binding, only: c_int, c_double
+ use, intrinsic :: ieee_arithmetic
+ implicit none
+ integer(c_int), value :: nrn, nuh
+ real(c_double), intent(in) :: rn(nrn)
+ real(c_double), intent(in) :: uh(nuh)
+ real(c_double), intent(out) :: q(nrn+nuh)
+ real(c_double) :: nan_val
+ integer :: t
+
+ nan_val = ieee_value(0.0_c_double, ieee_quiet_nan)
+ q = 0.0_c_double
+
+ do t = 1, nrn
+ if (ieee_is_nan(rn(t))) then
+ q(t:t+nuh-1) = nan_val
+ else if (rn(t) > 0.0_c_double) then
+ q(t:t+nuh-1) = q(t:t+nuh-1) + rn(t) * uh(1:nuh)
+ end if
+ end do
+
+end subroutine convolution
diff --git a/src/similarity_module.f90 b/src/similarity_module.f90
index a6b1176..0ca78f7 100644
--- a/src/similarity_module.f90
+++ b/src/similarity_module.f90
@@ -31,7 +31,7 @@ contains
case(4)
FUN => invrmse
case default
- FUN => invrmse
+ FUN => invkge
end select
end subroutine set_objective_function
diff --git a/src/transfR.c b/src/transfR.c
index cba56b4..31532fa 100644
--- a/src/transfR.c
+++ b/src/transfR.c
@@ -26,10 +26,20 @@ SEXP c_similarity(SEXP Rn, SEXP crit, SEXP nthreads) {
return sim_matrix;
}
-/* Registration of native routines for .Call */
+extern void convolution(double *rn, double *uh, int nrn, int nuh, double *q);
+SEXP c_convolution(SEXP rn, SEXP uh) {
+ const int nrn = length(rn);
+ const int nuh = length(uh);
+ SEXP q = PROTECT(allocVector(REALSXP, nrn + nuh));
+ convolution(REAL(rn), REAL(uh), nrn, nuh, REAL(q));
+ UNPROTECT(1);
+ return q;
+}
+
static const R_CallMethodDef CallEntries[] = {
{"c_gdist", (DL_FUNC) &c_gdist, 6},
{"c_similarity", (DL_FUNC) &c_similarity, 3},
+ {"c_convolution", (DL_FUNC) &c_convolution, 2},
{NULL, NULL, 0}
};
--
GitLab