fix issues on CRAN

R/RcppExports.R

@@ -93,8 +93,8 @@ chol_rank_one_update_test <- function(R, x) {
     invisible(.Call(`_dynamichazard_chol_rank_one_update_test`, R, x))
 }
 
-square_tri_inv_test <- function(R, out) {
-    invisible(.Call(`_dynamichazard_square_tri_inv_test`, R, out))
+square_tri_inv_test <- function(R) {
+    .Call(`_dynamichazard_square_tri_inv_test`, R)
 }
 
 symmetric_mat_chol_test <- function(A, out) {

src/BLAS_LAPACK/R_BLAS_LAPACK.cpp

@@ -68,7 +68,7 @@ namespace R_BLAS_LAPACK {
 
     if(info != 0){
       std::stringstream str;
-      str << "Got error code '" << info << "' when making rank one update of cholesky decomposition";
+      str << "Got error code '" << info << "' from 'dtrtri'";
       Rcpp::stop(str.str());
     }
   }

src/BLAS_LAPACK/arma_BLAS_LAPACK.cpp

@@ -21,15 +21,12 @@ void chol_rank_one_update(arma::mat &R, arma::vec x){
   R_BLAS_LAPACK::ddhazard_dchur(R.memptr(), copy_x.memptr(), n, n);
 }
 
-void square_tri_inv(const arma::mat &R, arma::mat &out){
-  // Take copy
-  out = R;
+void square_tri_inv(arma::mat &out){
   int n = out.n_cols;
-
   R_BLAS_LAPACK::square_tri_inv(out.memptr(), n, n);
 }
 
-void symmetric_mat_chol(const arma::mat& A, arma::mat & out){
+void symmetric_mat_chol(const arma::mat& A, arma::mat &out){
   // Take copy
   out = A;
 

src/PF/covarmat.cpp

@@ -2,9 +2,9 @@
 
 #ifdef _OPENMP
 class LockGuard {
-  omp_lock_t& m_lock;
+  omp_lock_t &m_lock;
 public:
-  explicit LockGuard(omp_lock_t& lock) : m_lock(lock){
+  explicit LockGuard(omp_lock_t &lock) : m_lock(lock){
     omp_set_lock(&m_lock);
   }
 
@@ -25,7 +25,7 @@ const arma::mat& covarmat::get_mat(output what) const {
 #pragma omp atomic read
   is_computed = *this_flag;
   if(!is_computed){
-    LockGuard guard(*lock.get());
+    LockGuard guard(*lock);
 #pragma omp atomic read
     is_computed = *this_flag;
     if(!is_computed){
@@ -49,7 +49,7 @@ const arma::mat& covarmat::get_mat(output what) const {
 #pragma omp atomic read
   is_computed = *this_flag;
   if(!is_computed){
-    LockGuard guard(*lock.get());
+    LockGuard guard(*lock);
 #pragma omp atomic read
     is_computed = *this_flag;
     if(!is_computed){
@@ -73,7 +73,7 @@ const arma::mat& covarmat::get_mat(output what) const {
 #pragma omp atomic read
   is_computed = *this_flag;
   if(!is_computed){
-    LockGuard guard(*lock.get());
+    LockGuard guard(*lock);
 #pragma omp atomic read
     is_computed = *this_flag;
     if(!is_computed){
@@ -106,8 +106,8 @@ const arma::mat& covarmat::inv() const {
 
 covarmat::covarmat(const covarmat &other): covarmat(other.mat()) { }
 
-covarmat::~covarmat(){
 #ifdef _OPENMP
+covarmat::~covarmat(){
   omp_destroy_lock(lock.get());
-#endif
 }
+#endif

src/PF/covarmat.h

@@ -12,8 +12,7 @@
 class covarmat {
 private:
 #ifdef _OPENMP
-  std::unique_ptr<omp_lock_t> lock =
-    std::unique_ptr<omp_lock_t>((new omp_lock_t()));
+  std::unique_ptr<omp_lock_t> lock;
 #endif
 
   enum output { e_mat, e_chol, e_chol_inv, e_inv };
@@ -48,13 +47,16 @@ class covarmat {
     chol_inv_(new arma::mat(arma::size(Q), arma::fill::zeros)),
     inv_     (new arma::mat(arma::size(Q), arma::fill::zeros)) {
 #ifdef _OPENMP
+    lock.reset(new omp_lock_t());
     omp_init_lock(lock.get());
 #endif
   }
 
   covarmat(const covarmat&);
 
+#ifdef _OPENMP
   ~covarmat();
+#endif
 };
 
 #endif

src/PF/dists.cpp

@@ -216,11 +216,8 @@ arma::vec set_at_risk_length(
 
 #ifdef _OPENMP
 /* openMP reductions */
-#pragma omp declare reduction(armaVP: arma::vec: omp_out += omp_in) \
-  initializer(omp_priv = arma::vec(omp_orig))
-
-#pragma omp declare reduction(armaMP: arma::mat: omp_out += omp_in) \
-  initializer(omp_priv = arma::mat(omp_orig))
+#pragma omp declare reduction(armaVP: arma::vec: omp_out += omp_in)
+#pragma omp declare reduction(armaMP: arma::mat: omp_out += omp_in)
 #endif
 
 template<class T>
@@ -284,10 +281,17 @@ arma::vec observational_cdist<T>::gradient(const arma::vec &coefs) const {
   /* compute gradient */
   arma::vec result(coefs.n_elem, arma::fill::zeros);
 #ifdef _OPENMP
-#pragma omp parallel for schedule(static) reduction(armaVP:result) \
-  if(multithreaded)
+  bool first_it = true;
+#pragma omp parallel for schedule(static) if(multithreaded)\
+  reduction(armaVP:result) firstprivate(first_it)
 #endif
   for(unsigned int i = 0; i < n; ++i){
+#ifdef _OPENMP
+    if(first_it){
+      result.zeros(coefs.n_elem);
+      first_it = false;
+    }
+#endif
     auto trunc_eta = this->truncate_eta(
       is_event[i], eta[i], exp(eta[i]), at_risk_length[i]);
     double d_l = this->d_log_like(is_event[i], trunc_eta, at_risk_length[i]);
@@ -312,10 +316,17 @@ arma::mat observational_cdist<T>::neg_Hessian(const arma::vec &coefs) const {
   /* compute Hessian */
   arma::mat result(coefs.n_elem, coefs.n_elem, arma::fill::zeros);
 #ifdef _OPENMP
+  bool first_it = true;
 #pragma omp parallel for schedule(static) reduction(armaMP:result) \
-  if(multithreaded)
+  if(multithreaded) firstprivate(first_it)
 #endif
   for(unsigned int i = 0; i < n; ++i){
+#ifdef _OPENMP
+    if(first_it){
+      result.zeros(coefs.n_elem, coefs.n_elem);
+      first_it = false;
+    }
+#endif
     auto trunc_eta = this->truncate_eta(
       is_event[i], eta[i], exp(eta[i]), at_risk_length[i]);
     double dd_l = this->dd_log_like(

src/RcppExports.cpp

@@ -440,14 +440,14 @@ BEGIN_RCPP
 END_RCPP
 }
 // square_tri_inv_test
-void square_tri_inv_test(const arma::mat& R, arma::mat& out);
-RcppExport SEXP _dynamichazard_square_tri_inv_test(SEXP RSEXP, SEXP outSEXP) {
+arma::mat square_tri_inv_test(const arma::mat& R);
+RcppExport SEXP _dynamichazard_square_tri_inv_test(SEXP RSEXP) {
 BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const arma::mat& >::type R(RSEXP);
-    Rcpp::traits::input_parameter< arma::mat& >::type out(outSEXP);
-    square_tri_inv_test(R, out);
-    return R_NilValue;
+    rcpp_result_gen = Rcpp::wrap(square_tri_inv_test(R));
+    return rcpp_result_gen;
 END_RCPP
 }
 // symmetric_mat_chol_test
@@ -902,7 +902,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_dynamichazard_test_get_ancestors", (DL_FUNC) &_dynamichazard_test_get_ancestors, 1},
     {"_dynamichazard_test_get_resample_idx_n_log_weight", (DL_FUNC) &_dynamichazard_test_get_resample_idx_n_log_weight, 3},
     {"_dynamichazard_chol_rank_one_update_test", (DL_FUNC) &_dynamichazard_chol_rank_one_update_test, 2},
-    {"_dynamichazard_square_tri_inv_test", (DL_FUNC) &_dynamichazard_square_tri_inv_test, 2},
+    {"_dynamichazard_square_tri_inv_test", (DL_FUNC) &_dynamichazard_square_tri_inv_test, 1},
     {"_dynamichazard_symmetric_mat_chol_test", (DL_FUNC) &_dynamichazard_symmetric_mat_chol_test, 2},
     {"_dynamichazard_tri_mat_times_vec_test", (DL_FUNC) &_dynamichazard_tri_mat_times_vec_test, 4},
     {"_dynamichazard_sym_mat_rank_one_update_test", (DL_FUNC) &_dynamichazard_sym_mat_rank_one_update_test, 3},

src/arma_BLAS_LAPACK.h

@@ -4,11 +4,11 @@
 #include "arma_n_rcpp.h"
 #include <memory>
 
-void symmetric_mat_chol(const arma::mat&, arma::mat &);
+void symmetric_mat_chol(const arma::mat&, arma::mat&);
 
 void chol_rank_one_update(arma::mat&, arma::vec);
 
-void square_tri_inv(const arma::mat&, arma::mat&);
+void square_tri_inv(arma::mat&);
 
 void tri_mat_times_vec(const arma::mat&, const arma::vec&, arma::vec&, bool);
 

src/ddhazard/SMA_solver.cpp

@@ -108,7 +108,8 @@ void SMA::solve(){
       arma::mat L;
       arma::mat L_inv = arma::inv_sympd(V); // only temporary
       symmetric_mat_chol(L_inv, L); // Cholesky decomposition of information matrix
-      square_tri_inv(L, L_inv); // V = L_inv^T * L_inv
+      L_inv = L;
+      square_tri_inv(L_inv); // V = L_inv^T * L_inv
       arma::vec inter_vec(L.n_cols);
 
       for(auto it = r_set.begin(); it != r_set.end(); it++){
@@ -140,7 +141,8 @@ void SMA::solve(){
         arma::vec rank_1_update_vec(x_ * sqrt(neg_second_d));
         rank_1_update_vec = p_dat.state_lp_inv->map(rank_1_update_vec).sv;
         chol_rank_one_update(L, rank_1_update_vec);
-        square_tri_inv(L, L_inv);
+        L_inv = L;
+        square_tri_inv(L_inv);
       }
 
       V = L_inv.t() * L_inv;

src/for_tests.cpp

@@ -211,8 +211,10 @@ void chol_rank_one_update_test(arma::mat &R, arma::vec x){
 }
 
 // [[Rcpp::export]]
-void square_tri_inv_test(const arma::mat &R, arma::mat &out){
-  return square_tri_inv(R, out);
+arma::mat square_tri_inv_test(const arma::mat &R){
+  arma::mat cp = R;
+  square_tri_inv(cp);
+  return cp;
 }
 
 // [[Rcpp::export]]

tests/testthat/test_LAPACK_BLAS.R

@@ -95,9 +95,7 @@ test_that("Square triangular inversion works followed by rank one update", {
     r_mat <- get_random_sym_post_def_mat(n)
 
     d1 <- t(chol(r_mat))
-    d2 <- matrix(0, ncol = n, nrow = n)
-
-    square_tri_inv_test(d1, d2)
+    d2 <- square_tri_inv_test(d1)
 
     expect_equal(solve(d1), d2)
   }