new configure to use MKL (Niagara)

Author: pachadotdev

configure

@@ -45,6 +45,22 @@ for flag in -fopenmp -qopenmp -openmp -xopenmp; do
   CXXFLAGS="$save_CXXFLAGS"
 done
 
+# Link BLAS/LAPACK against MKL if MKLROOT is set
+R_BLAS_LIBS=$("${R_HOME}/bin/R" CMD config BLAS_LIBS)
+R_LAPACK_LIBS=$("${R_HOME}/bin/R" CMD config LAPACK_LIBS)
+R_FLIBS=$("${R_HOME}/bin/R" CMD config FLIBS)
+
+: ${BLAS_LIBS:="$R_BLAS_LIBS"}
+: ${LAPACK_LIBS:="$R_LAPACK_LIBS"}
+: ${FLIBS:="$R_FLIBS"}
+
+if [ -n "$MKLROOT" ]; then
+  echo "detected MKLROOT=$MKLROOT → linking BLAS/LAPACK via MKL"
+  BLAS_LIBS="-L${MKLROOT}/lib/intel64_lin -lmkl_rt"
+  LAPACK_LIBS="$BLAS_LIBS"
+fi
+
+# Check for OpenMP support
 if [ -n "$OPENMP_FLAG" ]; then
   echo "OpenMP is enabled: $OPENMP_FLAG"
   PKG_CXXFLAGS="$OPENMP_FLAG $PKG_CXXFLAGS"
@@ -100,6 +116,8 @@ fi
 echo "creating src/Makevars"
 sed -e "s|@ncores@|${num_cores}|g" \
     -e "s|@OPTIMIZATION_FLAGS@|${OPTIMIZATION_FLAGS}|g" \
+    -e "s|@PKG_CXXFLAGS@|${PKG_CXXFLAGS}|g" \
+    -e "s|@PKG_LIBS@|${PKG_LIBS}|g" \
     src/Makevars.in > src/Makevars
 
 # Success

src/00_main.h

@@ -1,8 +1,10 @@
 #pragma once
 
 #include <armadillo.hpp>
+#include <cmath>
 #include <cpp11.hpp>
 #include <cpp11armadillo.hpp>
+#include <limits>
 #include <regex>
 #include <unordered_map>
 

src/02_center_variables.cpp

@@ -1,4 +1,3 @@
-// 02_center_variables.cpp (refactored using Armadillo types)
 #include "00_main.h"
 
 // Method of alternating projections (Halperin)
@@ -7,16 +6,16 @@ void center_variables_(mat &V, const vec &w, const list &klist,
                        const int &iter_interrupt, const int &iter_ssr) {
   // Auxiliary variables (fixed)
   const size_t I = static_cast<size_t>(max_iter), N = V.n_rows, P = V.n_cols,
-               K = klist.size(),
-               iter_check_interrupt0 = static_cast<size_t>(iter_interrupt),
-               iter_check_ssr0 = static_cast<size_t>(iter_ssr);
+               K = klist.size(), iint0 = static_cast<size_t>(iter_interrupt),
+               isr0 = static_cast<size_t>(iter_ssr);
   const double inv_sw = 1.0 / accu(w);
 
   // Auxiliary variables (storage)
-  size_t iter, j, k, p, J, iter_check_interrupt = iter_check_interrupt0,
-                           iter_check_ssr = iter_check_ssr0;
-  double coef, xbar, ratio, ssr, ssq, ratio0, ssr0;
+  size_t iter, iint, isr, j, jj, k, n, p, J, JJ;
+  double num, coef, xbar, ratio, ssr, ssq, ratio0, ssr0;
   vec x(N), x0(N), Gx(N), G2x(N), deltaG(N), delta2(N);
+
+  // Precompute groups into fields
   field<field<uvec>> group_indices(K);
   field<vec> group_inverse_weights(K);
   for (k = 0; k < K; ++k) {
@@ -26,59 +25,60 @@ void center_variables_(mat &V, const vec &w, const list &klist,
     vec invs(J);
     for (j = 0; j < J; ++j) {
       idxs(j) = as_uvec(as_cpp<integers>(jlist[j]));
-      ;
       invs(j) = 1.0 / accu(w.elem(idxs(j)));
     }
-    group_indices(k) = idxs;
-    group_inverse_weights(k) = invs;
+    group_indices(k) = std::move(idxs);
+    group_inverse_weights(k) = std::move(invs);
   }
 
+  // Single nested‐field projection helper
+  auto project = [&](vec &v) {
+    J = group_indices.n_elem;
+    for (j = 0; j < J; ++j) {
+      auto &idxs = group_indices(j);
+      auto &invs = group_inverse_weights(j);
+      JJ = idxs.n_elem;
+      for (jj = 0; jj < JJ; ++jj) {
+        const uvec &coords = idxs(jj);
+        xbar = dot(w.elem(coords), v.elem(coords)) * invs(jj);
+        v.elem(coords) -= xbar;
+      }
+    }
+  };
+
+  // Column‐wise centering
   for (p = 0; p < P; ++p) {
     x = V.col(p);
     ratio0 = std::numeric_limits<double>::infinity();
     ssr0 = std::numeric_limits<double>::infinity();
 
+    // reset per‐column interrupt
+    iint = iint0;
+    isr = isr0;
+
     for (iter = 0; iter < I; ++iter) {
-      if (iter == iter_check_interrupt) {
+      if (iter == iint) {
         check_user_interrupt();
-        iter_check_interrupt += iter_check_interrupt0;
+        iint += iint0;
       }
 
       x0 = x;
 
-      // Halperin projection
-      for (k = 0; k < K; ++k) {
-        field<uvec> &idxs = group_indices(k);
-        J = idxs.n_elem;
-        vec &invs = group_inverse_weights(k);
-        for (j = 0; j < J; ++j) {
-          const uvec &coords = idxs(j);
-          xbar = dot(w.elem(coords), x.elem(coords)) * invs(j);
-          x.elem(coords) -= xbar;
-        }
+      // 1) main projection
+      project(x);
+      num = 0.0;
+      for (n = 0; n < N; ++n) {
+        num += std::abs(x[n] - x0[n]) / (1.0 + std::abs(x0[n])) * w[n];
       }
-
-      // Convergence check
-      ratio = dot(abs(x - x0) / (1.0 + abs(x0)), w) * inv_sw;
+      ratio = num * inv_sw;
       if (ratio < tol)
         break;
 
-      // Acceleration every 5 iters
-      if (iter > 5 && (iter % 5) == 0) {
+      // 2) acceleration every 5 iters
+      if (iter >= 5 && (iter % 5) == 0) {
         Gx = x;
-        // Second projection
-        for (size_t k = 0; k < K; ++k) {
-          field<uvec> &idxs = group_indices(k);
-          vec &invs = group_inverse_weights(k);
-          for (j = 0; j < idxs.n_elem; ++j) {
-            const uvec &coords = idxs(j);
-            xbar = dot(w.elem(coords), Gx.elem(coords)) * invs(j);
-            Gx.elem(coords) -= xbar;
-          }
-        }
+        project(Gx);
         G2x = Gx;
-
-        // Compute deltas
         deltaG = G2x - x;
         delta2 = G2x - 2.0 * x + x0;
         ssq = dot(delta2, delta2);
@@ -92,17 +92,17 @@ void center_variables_(mat &V, const vec &w, const list &klist,
         }
       }
 
-      // SSR check
-      if (iter == iter_check_ssr && iter > 0) {
+      // 3) SSR‐based early exit
+      if (iter == isr && iter > 0) {
         check_user_interrupt();
-        iter_check_ssr += iter_check_ssr0;
+        isr += isr0;
         ssr = dot(x % x, w) * inv_sw;
-        if (fabs(ssr - ssr0) / (1.0 + fabs(ssr0)) < tol)
+        if (std::fabs(ssr - ssr0) / (1.0 + std::fabs(ssr0)) < tol)
           break;
         ssr0 = ssr;
       }
 
-      // Early exit
+      // 4) early exit
       if (iter > 3 && (ratio0 / ratio) < 1.1 && ratio < tol * 20)
         break;
       ratio0 = ratio;
@@ -114,8 +114,8 @@ void center_variables_(mat &V, const vec &w, const list &klist,
 
 [[cpp11::register]] doubles_matrix<>
 center_variables_r_(const doubles_matrix<> &V_r, const doubles &w_r,
-                    const list &klist, const double &tol, const int &max_iter,
-                    const int &iter_interrupt, const int &iter_ssr) {
+                    const list &klist, const double tol, const int max_iter,
+                    const int iter_interrupt, const int iter_ssr) {
   mat V = as_mat(V_r);
   center_variables_(V, as_col(w_r), klist, tol, max_iter, iter_interrupt,
                     iter_ssr);

src/Makevars

@@ -1,2 +1,11 @@
-PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DARMA_NO_DEBUG -DARMA_USE_BLAS -DARMA_USE_LAPACK -DARMA_OPENMP_THREADS=@ncores@ @OPTIMIZATION_FLAGS@
-PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) \
+  -DARMA_NO_DEBUG \
+  -DARMA_USE_BLAS \
+  -DARMA_USE_LAPACK \
+  -DARMA_OPENMP_THREADS=@ncores@ \
+  @OPTIMIZATION_FLAGS@
+
+PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) \
+  $(LAPACK_LIBS) \
+  $(BLAS_LIBS) \
+  $(FLIBS)

src/Makevars.in

@@ -13,10 +13,10 @@ extern "C" SEXP _capybara_check_linear_dependence_svd_(SEXP y, SEXP x, SEXP p) {
   END_CPP11
 }
 // 02_center_variables.cpp
-doubles_matrix<> center_variables_r_(const doubles_matrix<> & V_r, const doubles & w_r, const list & klist, const double & tol, const int & max_iter, const int & iter_interrupt, const int & iter_ssr);
+doubles_matrix<> center_variables_r_(const doubles_matrix<> & V_r, const doubles & w_r, const list & klist, const double tol, const int max_iter, const int iter_interrupt, const int iter_ssr);
 extern "C" SEXP _capybara_center_variables_r_(SEXP V_r, SEXP w_r, SEXP klist, SEXP tol, SEXP max_iter, SEXP iter_interrupt, SEXP iter_ssr) {
   BEGIN_CPP11
-    return cpp11::as_sexp(center_variables_r_(cpp11::as_cpp<cpp11::decay_t<const doubles_matrix<> &>>(V_r), cpp11::as_cpp<cpp11::decay_t<const doubles &>>(w_r), cpp11::as_cpp<cpp11::decay_t<const list &>>(klist), cpp11::as_cpp<cpp11::decay_t<const double &>>(tol), cpp11::as_cpp<cpp11::decay_t<const int &>>(max_iter), cpp11::as_cpp<cpp11::decay_t<const int &>>(iter_interrupt), cpp11::as_cpp<cpp11::decay_t<const int &>>(iter_ssr)));
+    return cpp11::as_sexp(center_variables_r_(cpp11::as_cpp<cpp11::decay_t<const doubles_matrix<> &>>(V_r), cpp11::as_cpp<cpp11::decay_t<const doubles &>>(w_r), cpp11::as_cpp<cpp11::decay_t<const list &>>(klist), cpp11::as_cpp<cpp11::decay_t<const double>>(tol), cpp11::as_cpp<cpp11::decay_t<const int>>(max_iter), cpp11::as_cpp<cpp11::decay_t<const int>>(iter_interrupt), cpp11::as_cpp<cpp11::decay_t<const int>>(iter_ssr)));
   END_CPP11
 }
 // 03_lm_fit.cpp