C code doubled omp vs NOomp; extension of turning bands 2D to function as 3D...

C code doubled omp vs NOomp; extension of turning bands 2D to function as 3D (necessary for spherical models)

DESCRIPTION

@@ -53,7 +53,7 @@ License: CC BY-SA 4.0 | GPL (>=2)
 URL: https://gitlab.hzdr.de/geomet/gmGeostats
 Copyright: (C) 2020 by Helmholtz-Zentrum Dresden-Rossendorf and Edith Cowan University;
          gsi.DS code by Hassan Talebi
-RoxygenNote: 7.1.0
+RoxygenNote: 7.1.1
 Roxygen: list(markdown = TRUE)
 Encoding: UTF-8
 VignetteBuilder: knitr

R/geostats.R

@@ -89,6 +89,7 @@ gsi.calcCgram <- function(X,Y,vgram,ijEqual=FALSE) {
 #' and nvar = nr of variables in vgram
 #' @useDynLib gmGeostats CMVTurningBands
 gsi.TurningBands <- function(X,vgram,nBands,nsim=NULL) {
+  # checks and preps
   stopifnot(length(as.integer(nBands))==1)
   d <- dim(vgram$nugget)[1]
   m <- ncol(X)
@@ -96,6 +97,23 @@ gsi.TurningBands <- function(X,vgram,nBands,nsim=NULL) {
   stopifnot( all(dim(vgram$nugget)==c(d,d)))
   stopifnot( all(dim(vgram$sill)==c(k,d,d)))
   stopifnot( all(dim(vgram$M)==c(k,m,m)))
+  # extend to 3D if 2D
+  m0 = 3
+  if(m<3){
+    m0 = m
+    X = cbind(X, 0, 0)
+    paveId = function(i){
+      M = diag(3)
+      M[1:m, 1:m] = vgram$M[i,,]
+      return(M)
+    }
+    M = sapply(1:k, paveId)
+    dim(M) = c(3,3,k)
+    M = aperm(M, c(1,2,3))
+    vgram$M = M
+    m = 3
+  }
+  # run!
   if( is.null(nsim) ) {
     erg<-.C("CMVTurningBands",
             dimX=checkInt(dim(X),2),
@@ -111,7 +129,8 @@ gsi.TurningBands <- function(X,vgram,nBands,nsim=NULL) {
             moreCgramData=checkDouble(vgram$data,k),
             PACKAGE = "gmGeostats"
     )
-    cbind(X,Z=t(structure(erg$Z,dim=c(d,nrow(X)))))
+      erg = cbind(X[, 1:m0],Z=t(structure(erg$Z,dim=c(d,nrow(X)))))
+      return(erg) 
   } else {
     nsim <- checkInt(nsim,1)
     stopifnot(nsim>0)
@@ -129,7 +148,8 @@ gsi.TurningBands <- function(X,vgram,nBands,nsim=NULL) {
             moreCgramData=checkDouble(vgram$data,k),
             PACKAGE = "gmGeostats"
     )
-    aperm(structure(erg$Z,dim=c(d,nrow(X),nsim)),c(2,1,3))
+    erg = aperm(structure(erg$Z,dim=c(d,nrow(X),nsim)),c(2,1,3))
+    return(erg) 
   }
 }
 
@@ -171,6 +191,22 @@ gsi.CondTurningBands <- function(Xout, Xin, Zin, vgram,
     normal<-FALSE
   }
   nsim <- checkInt(nsim,1)
+  # extend to 3D if 2D
+  m0 = 3
+  if(m<3){
+    m0 = m
+    X = cbind(X, 0, 0)
+    paveId = function(i){
+      M = diag(3)
+      M[1:m, 1:m] = vgram$M[i,,]
+      return(M)
+    }
+    M = sapply(1:k, paveId)
+    dim(M) = c(3,3,k)
+    M = aperm(M, c(1,2,3))
+    vgram$M = M
+    m = 3
+  }
   erg <- .C("CCondSim",
             dimZin=checkInt(dimZin,2),
             Zin   =checkDouble(Zin),

src/gmGeostats.c

@@ -421,6 +421,7 @@ void CMVTurningBands(
   (*bandSim[typeCgram[ss]])(n,projs,band,1.0,moreCgramData+ss);  
     /* this function takes the projs and returns on band the  the curve */
     /* ... multiply the eigenvector by the curve, and accumulate */
+#ifdef _OPENMP
 #pragma omp parallel for		             \
     if(!omp_in_parallel()&&0)		        \
       num_threads(omp_get_num_procs())	\
@@ -431,6 +432,14 @@ void CMVTurningBands(
       Z[d*i+j]+=d2;  
     }	  
     }
+#else
+    for(i=0;i<n;i++){ /* location */
+    for(j=0;j<d;j++){ /* variable */
+    d2 = sqrtSill[ss + *nCgrams *(j+d*ev)]*band[i]; 
+      Z[d*i+j]+=d2;  
+    }	  
+    }
+#endif
   }
   }
     } 
@@ -521,6 +530,7 @@ void CCondSim(
   Rprintf("starting conditioning by dual kriging\n");
   /* Kriging from simulated using Cinv */
   // Create differenes of obs and sim
+#ifdef _OPENMP
 #pragma omp parallel			               \
   if(!omp_in_parallel())			           \
     num_threads(omp_get_num_procs())		\
@@ -536,6 +546,16 @@ void CCondSim(
           }
       }
     }
+#else
+  for(int sim=0;sim<nsim;sim++) {
+    int shift=nd*sim;
+    for(int i=0;i<nin;i++)
+      for(int j=0;j<d;j++){
+        int k=d*i+j;
+        Z[k+shift]-=Zin[k];
+      }
+  }
+#endif
   /* Z[hinten] -= \hat{Z[hinten]}(Z[vorn]) = cov(hinten,vorn)%*% Cinv %*% Z[vorn] */
   /* dbuf = Cinv %*% Z[vorn] 
    Cinv in R ^ d*nin x d*nin
@@ -543,6 +563,7 @@ void CCondSim(
    
    */
   // Dual Kriging preparation
+#ifdef _OPENMP
 #pragma omp parallel for			           \
   if(!omp_in_parallel()&&0)			        \
     num_threads(omp_get_num_procs())		\
@@ -560,6 +581,21 @@ void CCondSim(
                dbuf+dmnin*sim,
                &oneI);
     }
+#else
+    for(sim=0;sim<nsim;sim++) {
+      F77_NAME(dgemv)(&No,
+               &dmnin,
+               &dmnin,
+               &one,
+               Cinv,
+               &dmnin,
+               Z+nd*sim,
+               &oneI,
+               &zero,
+               dbuf+dmnin*sim,
+               &oneI);
+    }
+#endif
     // /* points in*/
     //for(i=0;i<nin;i++)
     //for(j=0;j<m;j++)
@@ -606,6 +642,7 @@ void CCondSim(
       //}
       
       /* Z[,i] +=  t(cbuf(d*nin,d)) %*% dbuf(d*nin) */
+#ifdef _OPENMP
 #pragma omp parallel for			               \
       if(!omp_in_parallel()&&0)			        \
         num_threads(omp_get_num_procs())		\
@@ -623,7 +660,21 @@ void CCondSim(
                    Z+i*d+nd*sim,
                    &oneI);
         }
-        
+#else
+        for(sim=0;sim<nsim;sim++) {
+          F77_NAME(dgemv)(&Transposed,
+                   &dmnin,
+                   &d,
+                   &minus1,
+                   cbuf,
+                   &dmnin,
+                   dbuf+dmnin*sim,
+                   &oneI,
+                   &one,
+                   Z+i*d+nd*sim,
+                   &oneI);
+        }
+#endif
     }
 }
 
@@ -693,7 +744,8 @@ extern void anaForwardC(const int *dimX,
   const double h=((double)1.0)/steps;
   if( dimY[0]!=m )
     error("anaForwardC: x and y have different number of variables / rows");
-#pragma omp parallel for 
+#ifdef _OPENMP
+  #pragma omp parallel for 
   for(size_t i=0;i<nx;i++) {
     double v1[m];
     double v2[m];
@@ -712,9 +764,29 @@ extern void anaForwardC(const int *dimX,
         mx[j]+=0.5*h*(v1[j]+v2[j]);
     }
   }
-  
+#else
+  for(size_t i=0;i<nx;i++) {
+    double v1[m];
+    double v2[m];
+    double xx[m];
+    double *mx=x+m*i;
+    for(size_t s=0;s<steps;s++) {
+      // v1
+      anaV(v1,m,mx,s*h,dimY,y,wY,sigma0,sigma1);
+      // xx=x+h*v1
+      for(int j=0;j<m;j++)
+        xx[j]=mx[j]+h*v1[j];
+      // v2
+      anaV(v2,m,xx,(s+1)*h,dimY,y,wY,sigma0,sigma1);
+      // x+=0.5*h*(v1+v2)
+      for(int j=0;j<m;j++)
+        mx[j]+=0.5*h*(v1[j]+v2[j]);
+    }
+  }
+#endif
 }
 
+
 static R_NativePrimitiveArgType anaBackwardC_t[] = {  /* INTSXP,REALSXP */
     /* dimX,    x,     dimY    y     wY    stepsp   sigma0p    sigma1p */
     INTSXP,REALSXP,INTSXP,REALSXP,REALSXP,INTSXP, REALSXP,  REALSXP
@@ -739,6 +811,7 @@ extern void anaBackwardC(const int *dimX,
   const double h=((double)1.0)/steps;
   if( dimY[0]!=m )
     error("anaBackwardC: x and y have different number of variables / rows");
+#ifdef _OPENMP
 #pragma omp parallel for 
   for(size_t i=0;i<nx;i++) {
     double v1[m];
@@ -757,9 +830,27 @@ extern void anaBackwardC(const int *dimX,
       for(int j=0;j<m;j++)
         mx[j]-=0.5*h*(v1[j]+v2[j]);
     }
-    
   }
-  
+#else
+  for(size_t i=0;i<nx;i++) {
+    double v1[m];
+    double v2[m];
+    double xx[m];
+    double *mx=x+m*i;
+    for(size_t s=0;s<steps;s++) {
+      // v1
+      anaV(v1,m,mx,1-s*h,dimY,y,wY,sigma0,sigma1);
+      // xx=x-h*v1
+      for(int j=0;j<m;j++)
+        xx[j]=mx[j]-h*v1[j];
+      // v2
+      anaV(v2,m,xx,1-(s+1)*h,dimY,y,wY,sigma0,sigma1);
+      // x+=0.5*h*(v1+v2)
+      for(int j=0;j<m;j++)
+        mx[j]-=0.5*h*(v1[j]+v2[j]);
+    }
+  }
+#endif
 }