Armadillo 12.6.4

Author: eddelbuettel

inst/include/armadillo_bits/arma_version.hpp

@@ -23,7 +23,7 @@
 
 #define ARMA_VERSION_MAJOR 12
 #define ARMA_VERSION_MINOR 6
-#define ARMA_VERSION_PATCH 3
+#define ARMA_VERSION_PATCH 4
 #define ARMA_VERSION_NAME  "Cortisol Retox"
 
 

inst/include/armadillo_bits/auxlib_bones.hpp

@@ -433,6 +433,15 @@ class auxlib
 
   template<typename T>
   inline static bool rudimentary_sym_check(const Mat< std::complex<T> >& X);
+  
+  template<typename eT>
+  inline static typename get_pod_type<eT>::result norm1_gen(const Mat<eT>& A);
+  
+  template<typename eT>
+  inline static typename get_pod_type<eT>::result norm1_sym(const Mat<eT>& A);
+  
+  template<typename eT>
+  inline static typename get_pod_type<eT>::result norm1_band(const Mat<eT>& A, const uword KL, const uword KU);
   };
 
 

inst/include/armadillo_bits/auxlib_meat.hpp

@@ -121,7 +121,7 @@ auxlib::inv_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond)
     podarray<blas_int> ipiv(A.n_rows);
 
     arma_extra_debug_print("lapack::lange()");
-    norm_val = lapack::lange<eT>(&norm_id, &n, &n, A.memptr(), &lda, junk.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen(A) : lapack::lange<eT>(&norm_id, &n, &n, A.memptr(), &lda, junk.memptr());
 
     arma_extra_debug_print("lapack::getrf()");
     lapack::getrf(&n, &n, A.memptr(), &lda, ipiv.memptr(), &info);
@@ -333,7 +333,7 @@ auxlib::inv_sympd_rcond(Mat<eT>& A, bool& out_sympd_state, eT& out_rcond)
     podarray<T> work(A.n_rows);
 
     arma_extra_debug_print("lapack::lansy()");
-    norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym(A) : lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
 
     arma_extra_debug_print("lapack::potrf()");
     lapack::potrf(&uplo, &n, A.memptr(), &n, &info);
@@ -399,7 +399,7 @@ auxlib::inv_sympd_rcond(Mat< std::complex<T> >& A, bool& out_sympd_state, T& out
     podarray<T> work(A.n_rows);
 
     arma_extra_debug_print("lapack::lanhe()");
-    norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
+    norm_val = (has_blas_float_bug<T>::value) ? auxlib::norm1_sym(A) : lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
 
     arma_extra_debug_print("lapack::potrf()");
     lapack::potrf(&uplo, &n, A.memptr(), &n, &info);
@@ -4039,7 +4039,7 @@ auxlib::solve_square_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_ty
     podarray<blas_int> ipiv(A.n_rows + 2);  // +2 for paranoia
 
     arma_extra_debug_print("lapack::lange()");
-    norm_val = lapack::lange<eT>(&norm_id, &n, &n, A.memptr(), &lda, junk.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen(A) : lapack::lange<eT>(&norm_id, &n, &n, A.memptr(), &lda, junk.memptr());
 
     arma_extra_debug_print("lapack::getrf()");
     lapack::getrf<eT>(&n, &n, A.memptr(), &n, ipiv.memptr(), &info);
@@ -4368,7 +4368,7 @@ auxlib::solve_sympd_rcond(Mat<typename T1::pod_type>& out, bool& out_sympd_state
     podarray<T> work(A.n_rows);
 
     arma_extra_debug_print("lapack::lansy()");
-    norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym(A) : lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
 
     arma_extra_debug_print("lapack::potrf()");
     lapack::potrf<eT>(&uplo, &n, A.memptr(), &n, &info);
@@ -4446,7 +4446,7 @@ auxlib::solve_sympd_rcond(Mat< std::complex<typename T1::pod_type> >& out, bool&
     podarray<T> work(A.n_rows);
 
     arma_extra_debug_print("lapack::lanhe()");
-    norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym(A) : lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr());
 
     arma_extra_debug_print("lapack::potrf()");
     lapack::potrf<eT>(&uplo, &n, A.memptr(), &n, &info);
@@ -5387,7 +5387,7 @@ auxlib::solve_band_rcond_common(Mat<typename T1::elem_type>& out, typename T1::p
 
     arma_debug_assert_blas_size(AB,out);
 
-    char     norm_id  = '1';
+  //char     norm_id  = '1';
     char     trans    = 'N';
     blas_int n        = blas_int(N);  // assuming square matrix
     blas_int kl       = blas_int(KL);
@@ -5398,11 +5398,14 @@ auxlib::solve_band_rcond_common(Mat<typename T1::elem_type>& out, typename T1::p
     blas_int info     = blas_int(0);
     T        norm_val = T(0);
 
-    podarray<T>        junk(1);
+  //podarray<T>        junk(1);
     podarray<blas_int> ipiv(N + 2);  // +2 for paranoia
 
-    arma_extra_debug_print("lapack::langb()");
-    norm_val = lapack::langb<eT>(&norm_id, &n, &kl, &ku, AB.memptr(), &ldab, junk.memptr());
+    // // NOTE: lapack::langb() and lapack::gbtrf() use incompatible storage formats for the band matrix
+    // arma_extra_debug_print("lapack::langb()");
+    // norm_val = lapack::langb<eT>(&norm_id, &n, &kl, &ku, AB.memptr(), &ldab, junk.memptr());
+    
+    norm_val = auxlib::norm1_band(A,KL,KU);
 
     arma_extra_debug_print("lapack::gbtrf()");
     lapack::gbtrf<eT>(&n, &n, &kl, &ku, AB.memptr(), &ldab, ipiv.memptr(), &info);
@@ -6098,7 +6101,7 @@ auxlib::rcond(Mat<eT>& A)
     podarray<blas_int> ipiv( (std::min)(A.n_rows, A.n_cols) );
 
     arma_extra_debug_print("lapack::lange()");
-    norm_val = lapack::lange(&norm_id, &m, &n, A.memptr(), &lda, work.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen(A) : lapack::lange(&norm_id, &m, &n, A.memptr(), &lda, work.memptr());
 
     arma_extra_debug_print("lapack::getrf()");
     lapack::getrf(&m, &n, A.memptr(), &lda, ipiv.memptr(), &info);
@@ -6148,7 +6151,7 @@ auxlib::rcond(Mat< std::complex<T> >& A)
     podarray<blas_int> ipiv( (std::min)(A.n_rows, A.n_cols) );
 
     arma_extra_debug_print("lapack::lange()");
-    norm_val = lapack::lange(&norm_id, &m, &n, A.memptr(), &lda, junk.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen(A) : lapack::lange(&norm_id, &m, &n, A.memptr(), &lda, junk.memptr());
 
     arma_extra_debug_print("lapack::getrf()");
     lapack::getrf(&m, &n, A.memptr(), &lda, ipiv.memptr(), &info);
@@ -6196,7 +6199,7 @@ auxlib::rcond_sympd(Mat<eT>& A, bool& calc_ok)
     podarray<blas_int> iwork(  A.n_rows);
 
     arma_extra_debug_print("lapack::lansy()");
-    norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &lda, work.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym(A) : lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &lda, work.memptr());
 
     arma_extra_debug_print("lapack::potrf()");
     lapack::potrf(&uplo, &n, A.memptr(), &lda, &info);
@@ -6257,7 +6260,7 @@ auxlib::rcond_sympd(Mat< std::complex<T> >& A, bool& calc_ok)
     podarray< T> rwork(  A.n_rows);
 
     arma_extra_debug_print("lapack::lanhe()");
-    norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &lda, rwork.memptr());
+    norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym(A) : lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &lda, rwork.memptr());
 
     arma_extra_debug_print("lapack::potrf()");
     lapack::potrf(&uplo, &n, A.memptr(), &lda, &info);
@@ -6701,6 +6704,105 @@ auxlib::rudimentary_sym_check(const Mat< std::complex<T> >& X)
 
 
 
+template<typename eT>
+inline
+typename get_pod_type<eT>::result
+auxlib::norm1_gen(const Mat<eT>& A)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename get_pod_type<eT>::result T;
+  
+  if(A.n_elem == 0)  { return T(0); }
+  
+  const uword n_rows = A.n_rows;
+  const uword n_cols = A.n_cols;
+  
+  T max_val = T(0);
+  
+  for(uword c=0; c < n_cols; ++c)
+    {
+    const eT* colmem  = A.colptr(c);
+           T  acc_val = T(0);
+    
+    for(uword r=0; r < n_rows; ++r)  { acc_val += std::abs(colmem[r]); }
+    
+    max_val = (acc_val > max_val) ? acc_val : max_val;
+    }
+  
+  return max_val;
+  }
+
+
+
+template<typename eT>
+inline
+typename get_pod_type<eT>::result
+auxlib::norm1_sym(const Mat<eT>& A)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename get_pod_type<eT>::result T;
+  
+  if(A.n_elem == 0)  { return T(0); }
+  
+  const uword N = (std::min)(A.n_rows, A.n_cols);
+  
+  T max_val = T(0);
+  
+  for(uword col=0; col < N; ++col)
+    {
+    const eT* colmem  = A.colptr(col);
+           T  acc_val = T(0);
+    
+    for(uword c=0; c < col; ++c)  { acc_val += std::abs(A.at(col,c)); }
+    
+    for(uword r=col; r < N; ++r)  { acc_val += std::abs(colmem[r]);   }
+    
+    max_val = (acc_val > max_val) ? acc_val : max_val;
+    }
+  
+  return max_val;
+  }
+
+
+
+template<typename eT>
+inline
+typename get_pod_type<eT>::result
+auxlib::norm1_band(const Mat<eT>& A, const uword KL, const uword KU)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename get_pod_type<eT>::result T;
+  
+  if(A.n_elem == 0)  { return T(0); }
+  
+  const uword n_rows = A.n_rows;
+  const uword n_cols = A.n_cols;
+  
+  T max_val = T(0);
+  
+  for(uword c=0; c < n_cols; ++c)
+    {
+    const eT* colmem  = A.colptr(c);
+           T  acc_val = T(0);
+    
+    // use values only from main diagonal + KU upper diagonals + KL lower diagonals
+    
+    const uword start = ( c       > KU    ) ? (c - KU) : 0;
+    const uword   end = ((c + KL) < n_rows) ? (c + KL) : (n_rows-1);
+    
+    for(uword r=start; r <= end; ++r)  { acc_val += std::abs(colmem[r]); }
+    
+    max_val = (acc_val > max_val) ? acc_val : max_val;
+    }
+  
+  return max_val;
+  }
+
+
+
 //
 
 

inst/include/armadillo_bits/compiler_setup.hpp

@@ -114,8 +114,10 @@
 
 
 #if defined(__APPLE__) || defined(__apple_build_version__)
-  #undef  ARMA_BLAS_SDOT_BUG
-  #define ARMA_BLAS_SDOT_BUG
+  // NOTE: Apple accelerate framework has broken implementations of functions that return a float value,
+  // NOTE: such as sdot(), slange(), clange(), slansy(), clanhe(), slangb()
+  #undef  ARMA_BLAS_FLOAT_BUG
+  #define ARMA_BLAS_FLOAT_BUG
 
   // #undef  ARMA_HAVE_POSIX_MEMALIGN
   // NOTE: posix_memalign() is available since macOS 10.6 (late 2009 onwards)

inst/include/armadillo_bits/op_det_bones.hpp

@@ -31,7 +31,6 @@ class op_det
     {
     static constexpr uword n2 = row + col*2;
     static constexpr uword n3 = row + col*3;
-    static constexpr uword n4 = row + col*4;
     };
 
   template<typename T1>
@@ -48,9 +47,6 @@ class op_det
 
   template<typename eT>
   arma_cold inline static eT apply_tiny_3x3(const Mat<eT>& X);
-  
-  template<typename eT>
-  arma_cold inline static eT apply_tiny_4x4(const Mat<eT>& X);
   };
 
 

inst/include/armadillo_bits/op_det_meat.hpp

@@ -58,7 +58,7 @@ op_det::apply_direct(typename T1::elem_type& out_val, const Base<typename T1::el
   if(N == 0)  { out_val = eT(1); return true; }
   if(N == 1)  { out_val =  A[0]; return true; }
 
-  if((is_cx<eT>::no) && (N <= 4))
+  if((is_cx<eT>::no) && (N <= 3))
     {
     constexpr T det_min =        std::numeric_limits<T>::epsilon();
     constexpr T det_max = T(1) / std::numeric_limits<T>::epsilon();
@@ -67,7 +67,6 @@ op_det::apply_direct(typename T1::elem_type& out_val, const Base<typename T1::el
 
     if(N == 2)  { det_val = op_det::apply_tiny_2x2(A); }
     if(N == 3)  { det_val = op_det::apply_tiny_3x3(A); }
-    if(N == 4)  { det_val = op_det::apply_tiny_4x4(A); }
 
     const T abs_det_val = std::abs(det_val);
 
@@ -176,77 +175,4 @@ op_det::apply_tiny_3x3(const Mat<eT>& X)
 
 
 
-template<typename eT>
-inline
-eT
-op_det::apply_tiny_4x4(const Mat<eT>& X)
-  {
-  arma_extra_debug_sigprint();
-  
-  const eT* Xm = X.memptr();
-  
-  const eT val_03_12 = Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4];
-  const eT val_02_13 = Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4];
-  const eT val_03_11 = Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4];
-  
-  const eT val_01_13 = Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4];
-  const eT val_02_11 = Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4];
-  const eT val_01_12 = Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4];
-  
-  const eT val_03_10 = Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4];
-  const eT val_00_13 = Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4];
-  const eT val_02_10 = Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4];
-  const eT val_00_12 = Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4];
-  
-  const eT val_01_10 = Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4];
-  const eT val_00_11 = Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4];
-  
-  const eT val_21_30 = Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4];
-  const eT val_22_30 = Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4];
-  const eT val_23_30 = Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4];
-  
-  const eT val_20_31 = Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4];
-  const eT val_22_31 = Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4];
-  const eT val_23_31 = Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4];
-  
-  const eT val_20_32 = Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4];
-  const eT val_21_32 = Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4];
-  const eT val_23_32 = Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4];
-  
-  const eT val_20_33 = Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4];
-  const eT val_21_33 = Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4];
-  const eT val_22_33 = Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4];
-  
-  const eT val = \
-      val_03_12 * val_21_30 \
-    - val_02_13 * val_21_30 \
-    - val_03_11 * val_22_30 \
-    + val_01_13 * val_22_30 \
-    + val_02_11 * val_23_30 \
-    - val_01_12 * val_23_30 \
-    - val_03_12 * val_20_31 \
-    + val_02_13 * val_20_31 \
-    + val_03_10 * val_22_31 \
-    - val_00_13 * val_22_31 \
-    - val_02_10 * val_23_31 \
-    + val_00_12 * val_23_31 \
-    + val_03_11 * val_20_32 \
-    - val_01_13 * val_20_32 \
-    - val_03_10 * val_21_32 \
-    + val_00_13 * val_21_32 \
-    + val_01_10 * val_23_32 \
-    - val_00_11 * val_23_32 \
-    - val_02_11 * val_20_33 \
-    + val_01_12 * val_20_33 \
-    + val_02_10 * val_21_33 \
-    - val_00_12 * val_21_33 \
-    - val_01_10 * val_22_33 \
-    + val_00_11 * val_22_33 \
-    ;
-  
-  return val;
-  }
-
-
-
 //! @}