perf: fortran subroutine for convolution

DESCRIPTION

 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"),

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)])

src/convolution.f90

+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

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

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}
 };