Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

preparation for openmp #38

Closed
wants to merge 1 commit into from
Closed

preparation for openmp #38

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 1 addition & 6 deletions src/Makevars
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
95 changes: 46 additions & 49 deletions src/bsvars_sign.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,19 @@
#include "progress.hpp"
#include "Rcpp/Rmath.h"
#include <bsvars.h>
// #include <omp.h>

#ifdef _OPENMP
#include <omp.h>
#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"
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -83,53 +95,45 @@ Rcpp::List bsvar_sign_cpp(

int S_hyper = hypers.n_cols - 1;
int prior_nu = as<int>(prior["nu"]);
int post_nu = prior_nu + T;

double w, mu, delta, lambda;

vec hyper, psi;
vec prior_v = as<mat>(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<mat>(prior["B"]);
mat Ysoc = as<mat>(prior["Ysoc"]);
mat Xsoc = as<mat>(prior["Xsoc"]);
mat Ysur = as<mat>(prior["Ysur"]);
mat Xsur = as<mat>(prior["Xsur"]);

field<mat> 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<vec>(K - N * p)));
prior_S = diagmat(psi);
mat prior_V = diagmat(prior_v % join_vert(lambda * lambda * repmat(1 / psi, p, 1),
ones<vec>(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<mat> 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) {

Expand All @@ -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));
Expand All @@ -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

Expand Down
Loading