From cb20a106d1ee396bfe1d14ac3365f2f19cf4545e Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Wed, 17 Jul 2024 16:20:26 +1000 Subject: [PATCH] include ifdef --- src/Makevars | 7 +--- src/bsvars_sign.cpp | 95 ++++++++++++++++++++++----------------------- 2 files changed, 47 insertions(+), 55 deletions(-) diff --git a/src/Makevars b/src/Makevars index 62a255e..48f15f8 100644 --- a/src/Makevars +++ b/src/Makevars @@ -16,9 +16,4 @@ #CXX_STD = CXX11 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) - -#PKG_CXXFLAGS += -Xclang -fopenmp -#PKG_FFLAGS += $(SHLIB_OPENMP_FFLAGS) -#PKG_CFLAGS += $(SHLIB_OPENMP_CFLAGS) -#PKG_LIBS += -lomp \ No newline at end of file +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \ No newline at end of file diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 1433621..292e2ca 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -3,7 +3,19 @@ #include "progress.hpp" #include "Rcpp/Rmath.h" #include -// #include + +#ifdef _OPENMP + #include +#else + // for machines with compilers void of openmp support + #define omp_get_num_threads() 1 + #define omp_get_thread_num() 0 + #define omp_get_max_threads() 1 + #define omp_get_thread_limit() 1 + #define omp_get_num_procs() 1 + #define omp_set_nested(a) // empty statement to remove the call + #define omp_get_wtime() 0 +#endif #include "sample_hyper.h" #include "sample_Q.h" @@ -46,16 +58,16 @@ Rcpp::List bsvar_sign_cpp( } // Progress bar setup - double num_threads = 1; - // #pragma omp parallel - // { - // num_threads = omp_get_num_threads(); - // } + double num_threads; + #pragma omp parallel + { + num_threads = omp_get_num_threads(); + } vec prog_rep_points = arma::round(arma::linspace(0, S / num_threads, 50)); if (show_progress) { Rcout << "**************************************************|" << endl; - Rcout << " bsvarSIGNs: Bayesian Structural VAR with zero, |" << endl; - Rcout << " sign and narrative restrictions |" << endl; + Rcout << " bsvarSIGNs: Bayesian Structural VAR with sign, |" << endl; + Rcout << " zero and narrative restrictions |" << endl; Rcout << "**************************************************|" << endl; // Rcout << " Gibbs sampler for the SVAR model |" << endl; // Rcout << "**************************************************|" << endl; @@ -83,53 +95,45 @@ Rcpp::List bsvar_sign_cpp( int S_hyper = hypers.n_cols - 1; int prior_nu = as(prior["nu"]); - int post_nu = prior_nu + T; - double w, mu, delta, lambda; - - vec hyper, psi; vec prior_v = as(prior["V"]).diag(); - mat B, Sigma, chol_Sigma, h_invp, Q, shocks; - mat prior_V, prior_S, post_B, post_V, post_S; - mat Ystar, Xstar, Yplus, Xplus; mat prior_B = as(prior["B"]); mat Ysoc = as(prior["Ysoc"]); mat Xsoc = as(prior["Xsoc"]); mat Ysur = as(prior["Ysur"]); mat Xsur = as(prior["Xsur"]); - field result; - - // #pragma omp parallel for private(hyper, mu, delta, lambda, psi, prior_V, prior_S, Ystar, Xstar, Yplus, Xplus, result, post_B, post_V, post_S, Sigma, chol_Sigma, B, h_invp, Q, shocks, w) + #pragma omp parallel for for (int s = 0; s < S; s++) { - hyper = hypers.col(randi(distr_param(0, S_hyper))); - mu = hyper(0); - delta = hyper(1); - lambda = hyper(2); - psi = hyper.rows(3, N + 2); + vec hyper = hypers.col(randi(distr_param(0, S_hyper))); + double mu = hyper(0); + double delta = hyper(1); + double lambda = hyper(2); + vec psi = hyper.rows(3, N + 2); // update Minnesota prior - prior_V = diagmat(prior_v % - join_vert(lambda * lambda * repmat(1 / psi, p, 1), - ones(K - N * p))); - prior_S = diagmat(psi); + mat prior_V = diagmat(prior_v % join_vert(lambda * lambda * repmat(1 / psi, p, 1), + ones(K - N * p))); + mat prior_S = diagmat(psi); // update dummy observation prior - Ystar = join_vert(Ysoc / mu, Ysur / delta); - Xstar = join_vert(Xsoc / mu, Xsur / delta); - Yplus = join_vert(Ystar, Y); - Xplus = join_vert(Xstar, X); + mat Ystar = join_vert(Ysoc / mu, Ysur / delta); + mat Xstar = join_vert(Xsoc / mu, Xsur / delta); + mat Yplus = join_vert(Ystar, Y); + mat Xplus = join_vert(Xstar, X); // posterior parameters - result = niw_cpp(Yplus, Xplus, prior_B, prior_V, prior_S, prior_nu); - post_B = result(0); - post_V = result(1); - post_S = result(2); - post_nu = as_scalar(result(3)); + field result = niw_cpp(Yplus, Xplus, prior_B, prior_V, prior_S, prior_nu); + mat post_B = result(0); + mat post_V = result(1); + mat post_S = result(2); + int post_nu = as_scalar(result(3)); - w = 0; + double w = 0; + + mat Sigma, chol_Sigma, B, h_invp, Q, shocks; while (w == 0) { @@ -141,8 +145,8 @@ Rcpp::List bsvar_sign_cpp( B = rmatnorm_cpp(post_B, post_V, Sigma); h_invp = inv(trimatl(chol_Sigma)); // lower tri, h(Sigma) is upper tri - result = sample_Q(p, Y, X, B, h_invp, chol_Sigma, prior, - VB, sign_irf, sign_narrative, sign_B, Z, max_tries); + result = sample_Q(p, Y, X, B, h_invp, chol_Sigma, prior, VB, + sign_irf, sign_narrative, sign_B, Z, max_tries); Q = result(0); shocks = result(1); w = as_scalar(result(2)); @@ -157,16 +161,9 @@ Rcpp::List bsvar_sign_cpp( posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = shocks; - // Increment progress bar - if (any(prog_rep_points == s)) bar.increment(); - - // if (omp_get_thread_num() == 0) { - // // Check for user interrupts - // if (s % 10 == 0) checkUserInterrupt(); - // - // // Increment progress bar - // if (any(prog_rep_points == s)) bar.increment(); - // } + if (omp_get_thread_num() == 0 and any(prog_rep_points == s)) { + bar.increment(); + } } // END s loop