From b0ae0361d92a37387d52fc6f0fe8f6aacafba7eb Mon Sep 17 00:00:00 2001
From: Raimon Tolosana-Delgado <tolosa53@fwg206.ad.fz-rossendorf.de>
Date: Thu, 1 Oct 2020 14:10:51 +0200
Subject: [PATCH] C code doubled omp vs NOomp; extension of turning bands 2D to
function as 3D (necessary for spherical models)
---
DESCRIPTION | 2 +-
R/geostats.R | 40 ++++++++++++++++++-
src/gmGeostats.c | 101 ++++++++++++++++++++++++++++++++++++++++++++---
3 files changed, 135 insertions(+), 8 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 075e0b9..b0343f8 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/R/geostats.R b/R/geostats.R
index 8e1ae85..45d1e54 100644
--- a/R/geostats.R
+++ b/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),
diff --git a/src/gmGeostats.c b/src/gmGeostats.c
index 1690f98..58d8e05 100644
--- a/src/gmGeostats.c
+++ b/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
}
--
GitLab