patch bump

docs.html

@@ -8118,7 +8118,7 @@ Generate a scalar, vector, matrix or cube with the elements set to random values
       <td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top;"><i>x<sup>&thinsp;a-1</sup> exp(&thinsp;-x&thinsp;/&thinsp;b&thinsp;)</i></td>
     </tr>
     <tr>
-      <td style="vertical-align: top;"><i>p<font size=+1>(</font>x&thinsp;|&thinsp;a,b<font size=+1>)</font></i></td><td style="vertical-align: top;">&nbsp;<font size=+1>=</font>&nbsp;</td><td style="vertical-align: top;"><font size=+1><b><hr></b></font></td>
+      <td style="vertical-align: top;"><i>p<font size=+1>(</font>&thinsp;x&thinsp;|&thinsp;a,b&thinsp;<font size=+1>)</font></i></td><td style="vertical-align: top;">&nbsp;<font size=+1>=</font>&nbsp;</td><td style="vertical-align: top;"><font size=+1><b><hr></b></font></td>
     </tr>
     <tr>
       <td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top; text-align: center;"><i>b<sup>&thinsp;a</sup> &Gamma;(a)</i></td>
@@ -8176,6 +8176,7 @@ fmat B2 = randg&lt;fmat&gt;(10, 10, distr_param(2,1));
 <li><a href="#imbue">.imbue()</a></li>
 <li><a href="#size">size()</a></li>
 <li><a href="#rng_seed">RNG seed setting</a></li>
+<li><a href="https://mathworld.wolfram.com/GammaDistribution.html">gamma distribution in MathWorld</a></li>
 <li><a href="https://en.wikipedia.org/wiki/Gamma_distribution">gamma distribution in Wikipedia</a></li>
 </ul>
 </li>
@@ -11876,9 +11877,6 @@ See also:
 <a name="sort_index"></a>
 <b>sort_index( X )</b>
 <br><b>sort_index( X, sort_direction )</b>
-<br>
-<br><b>stable_sort_index( X )</b>
-<br><b>stable_sort_index( X, sort_direction )</b>
 <ul>
 <li>Return a vector which describes the sorted order of the elements of <i>X</i>
 (ie. it contains the indices of the elements of <i>X</i>)
@@ -11894,8 +11892,6 @@ See also:
 <br>
 <li>The <i>sort_direction</i> argument is optional; <i>sort_direction</i> is either <code>"ascend"</code> or <code>"descend"</code>; by default <code>"ascend"</code> is used</li>
 <br>
-<li>The <i>stable_sort_index()</i> variant preserves the relative order of elements with equivalent values</li>
-<br>
 <li>For matrices and vectors with complex numbers, sorting is via absolute values</li>
 <br>
 <li>
@@ -21448,7 +21444,7 @@ see the <a href="https://arma.sourceforge.net/faq.html#license">Questions page</
 <li>added handling of <a href="#diag">diagonal views</a> by sparse matrices
 <li>expanded <a href="#repmat">repmat()</a> to handle sparse matrices
 <li>expanded <a href="#join">join_rows()</a> and <a href="#join">join_cols()</a> to handle sparse matrices
-<li><a href="#sort_index">sort_index()</a> and <a href="#sort_index">stable_sort_index()</a> have been placed in the delayed operations framework for increased efficiency
+<li><a href="#sort_index">sort_index()</a> has been placed in the delayed operations framework for increased efficiency
 <li>use of <a href="#uword">64 bit integers</a> is automatically enabled when using a C++11 compiler</li>
 </ul>
 </li>
@@ -21660,7 +21656,7 @@ eg. <i>mat&nbsp;X(4,&nbsp;5,&nbsp;fill::zeros)</i></li>
 <li>faster handling of compound expressions with submatrices and subcubes</li>
 <li>faster <a href="#trace">trace()</a></li>
 <li>added support for loading matrices as text files with <i>NaN</i> and <i>Inf</i> elements</li>
-<li>added <a href="#sort_index">stable_sort_index()</a>, which preserves the relative order of elements with equivalent values</li>
+<!-- <li>added <a href="#sort_index">stable_sort_index()</a>, which preserves the relative order of elements with equivalent values</li> -->
 <li>added handling of <a href="#SpMat">sparse matrices</a> by <a href="#stats_fns">mean()</a>, <a href="#stats_fns">var()</a>, <a href="#norm">norm()</a>, <a href="#abs">abs()</a>, <a href="#misc_fns">square()</a>, <a href="#misc_fns">sqrt()</a></li>
 <li>added saving and loading of sparse matrices in <i>arma_binary</i> format</li>
 </ul>

include/armadillo_bits/arma_version.hpp

@@ -23,7 +23,7 @@
 
 #define ARMA_VERSION_MAJOR 15
 #define ARMA_VERSION_MINOR 0
-#define ARMA_VERSION_PATCH 1
+#define ARMA_VERSION_PATCH 2
 #define ARMA_VERSION_NAME  "Medium Roast"
 
 

include/armadillo_bits/mul_gemm.hpp

@@ -59,6 +59,156 @@ struct gemm_emul_tinysq
 
 
 
+struct gemm_emul_large_mp_helper
+  {
+  template<typename eT>
+  arma_hot
+  inline
+  static
+  void
+  copy_row(eT* out_mem, const Mat<eT>& in, const uword row)
+    {
+    const uword n_rows = in.n_rows;
+    const uword n_cols = in.n_cols;
+    
+    const eT* in_mem_row = in.memptr() + row;
+    
+    for(uword i=0; i < n_cols; ++i)
+      {
+      out_mem[i] = (*in_mem_row);
+      
+      in_mem_row += n_rows;
+      }
+    }
+  };
+
+
+
+#if defined(ARMA_USE_OPENMP)
+//! emulation of gemm(), for non-complex matrices only, as it assumes only simple transposes (ie. doesn't do hermitian transposes)
+//! parallelised version
+template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
+struct gemm_emul_large_mp
+  {
+  template<typename eT, typename TA, typename TB>
+  arma_hot
+  inline
+  static
+  void
+  apply
+    (
+          Mat<eT>& C,
+    const TA&      A,
+    const TB&      B,
+    const eT       alpha = eT(1),
+    const eT       beta  = eT(0)
+    )
+    {
+    arma_debug_sigprint();
+    
+    const uword A_n_rows = A.n_rows;
+    const uword A_n_cols = A.n_cols;
+    
+    const uword B_n_rows = B.n_rows;
+    const uword B_n_cols = B.n_cols;
+    
+    if( (do_trans_A == false) && (do_trans_B == false) )
+      {
+      const uword n_threads = uword(mp_thread_limit::get());
+      
+      podarray<eT> tmp(A_n_cols * n_threads, arma_nozeros_indicator());
+      
+      eT* tmp_mem = tmp.memptr();
+      
+      #pragma omp parallel for schedule(static) num_threads(int(n_threads))
+      for(uword row_A=0; row_A < A_n_rows; ++row_A)
+        {
+        const uword thread_id = uword(omp_get_thread_num());
+        
+        eT* A_rowdata = tmp_mem + (A_n_cols * thread_id);
+        
+        gemm_emul_large_mp_helper::copy_row(A_rowdata, A, row_A);
+        
+        for(uword col_B=0; col_B < B_n_cols; ++col_B)
+          {
+          const eT acc = op_dot::direct_dot(B_n_rows, A_rowdata, B.colptr(col_B));
+          
+               if( (use_alpha == false) && (use_beta == false) )  { C.at(row_A,col_B) =       acc;                          }
+          else if( (use_alpha == true ) && (use_beta == false) )  { C.at(row_A,col_B) = alpha*acc;                          }
+          else if( (use_alpha == false) && (use_beta == true ) )  { C.at(row_A,col_B) =       acc + beta*C.at(row_A,col_B); }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { C.at(row_A,col_B) = alpha*acc + beta*C.at(row_A,col_B); }
+          }
+        }
+      }
+    else
+    if( (do_trans_A == true) && (do_trans_B == false) )
+      {
+      const int n_threads = mp_thread_limit::get();
+      
+      #pragma omp parallel for schedule(static) num_threads(n_threads)
+      for(uword col_A=0; col_A < A_n_cols; ++col_A)
+        {
+        // col_A is interpreted as row_A when storing the results in matrix C
+        
+        const eT* A_coldata = A.colptr(col_A);
+        
+        for(uword col_B=0; col_B < B_n_cols; ++col_B)
+          {
+          const eT acc = op_dot::direct_dot(B_n_rows, A_coldata, B.colptr(col_B));
+          
+               if( (use_alpha == false) && (use_beta == false) )  { C.at(col_A,col_B) =       acc;                          }
+          else if( (use_alpha == true ) && (use_beta == false) )  { C.at(col_A,col_B) = alpha*acc;                          }
+          else if( (use_alpha == false) && (use_beta == true ) )  { C.at(col_A,col_B) =       acc + beta*C.at(col_A,col_B); }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { C.at(col_A,col_B) = alpha*acc + beta*C.at(col_A,col_B); }
+          }
+        }
+      }
+    else
+    if( (do_trans_A == false) && (do_trans_B == true) )
+      {
+      Mat<eT> BB;
+      op_strans::apply_mat_noalias(BB, B);
+      
+      gemm_emul_large_mp<false, false, use_alpha, use_beta>::apply(C, A, BB, alpha, beta);
+      }
+    else
+    if( (do_trans_A == true) && (do_trans_B == true) )
+      {
+      // using trans(A)*trans(B) = trans(B*A) equivalency; assuming no hermitian transpose
+      
+      const uword n_threads = uword(mp_thread_limit::get());
+      
+      podarray<eT> tmp(B_n_cols * n_threads, arma_nozeros_indicator());
+      
+      eT* tmp_mem = tmp.memptr();
+      
+      #pragma omp parallel for schedule(static) num_threads(int(n_threads))
+      for(uword row_B=0; row_B < B_n_rows; ++row_B)
+        {
+        const uword thread_id = uword(omp_get_thread_num());
+        
+        eT* B_rowdata = tmp_mem + (B_n_cols * thread_id);
+        
+        gemm_emul_large_mp_helper::copy_row(B_rowdata, B, row_B);
+        
+        for(uword col_A=0; col_A < A_n_cols; ++col_A)
+          {
+          const eT acc = op_dot::direct_dot(A_n_rows, B_rowdata, A.colptr(col_A));
+          
+               if( (use_alpha == false) && (use_beta == false) )  { C.at(col_A,row_B) =       acc;                          }
+          else if( (use_alpha == true ) && (use_beta == false) )  { C.at(col_A,row_B) = alpha*acc;                          }
+          else if( (use_alpha == false) && (use_beta == true ) )  { C.at(col_A,row_B) =       acc + beta*C.at(col_A,row_B); }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { C.at(col_A,row_B) = alpha*acc + beta*C.at(col_A,row_B); }
+          }
+        }
+      }
+    }
+  
+  };
+#endif
+
+
+
 //! emulation of gemm(), for non-complex matrices only, as it assumes only simple transposes (ie. doesn't do hermitian transposes)
 template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
 struct gemm_emul_large
@@ -78,13 +228,28 @@ struct gemm_emul_large
     )
     {
     arma_debug_sigprint();
-
+    
     const uword A_n_rows = A.n_rows;
     const uword A_n_cols = A.n_cols;
     
     const uword B_n_rows = B.n_rows;
     const uword B_n_cols = B.n_cols;
     
+    #if defined(ARMA_USE_OPENMP)
+      {
+      // TODO: replace with more sophisticated threshold mechanism
+      
+      constexpr uword threshold = uword(30);
+      
+      if( (A_n_rows >= threshold) && (A_n_cols >= threshold) && (B_n_rows >= threshold) && (B_n_cols >= threshold) && (mp_thread_limit::in_parallel() == false) )
+        {
+        gemm_emul_large_mp<do_trans_A, do_trans_B, use_alpha, use_beta>::apply(C,A,B,alpha,beta);
+        
+        return;
+        }
+      }
+    #endif
+    
     if( (do_trans_A == false) && (do_trans_B == false) )
       {
       arma_aligned podarray<eT> tmp(A_n_cols);

include/armadillo_bits/mul_gemv.hpp

@@ -203,6 +203,74 @@ struct gemv_emul_helper
 
 
 
+#if defined(ARMA_USE_OPENMP)
+//! Partial emulation of BLAS gemv().
+//! 'y' is assumed to have been set to the correct size (ie. taking into account the transpose)
+//! parallelised version
+template<const bool do_trans_A=false, const bool use_alpha=false, const bool use_beta=false>
+struct gemv_emul_mp
+  {
+  template<typename eT, typename TA>
+  arma_hot
+  inline
+  static
+  void
+  apply( eT* y, const TA& A, const eT* x, const eT alpha = eT(1), const eT beta = eT(0) )
+    {
+    arma_debug_sigprint();
+    
+    const int n_threads = mp_thread_limit::get();
+    
+    const uword A_n_rows = A.n_rows;
+    const uword A_n_cols = A.n_cols;
+    
+    if(do_trans_A == false)
+      {
+      #pragma omp parallel for schedule(static) num_threads(n_threads)
+      for(uword row=0; row < A_n_rows; ++row)
+        {
+        const eT acc = gemv_emul_helper::dot_row_col(A, x, row, A_n_cols);
+        
+             if( (use_alpha == false) && (use_beta == false) )  { y[row] =       acc;               }
+        else if( (use_alpha == true ) && (use_beta == false) )  { y[row] = alpha*acc;               }
+        else if( (use_alpha == false) && (use_beta == true ) )  { y[row] =       acc + beta*y[row]; }
+        else if( (use_alpha == true ) && (use_beta == true ) )  { y[row] = alpha*acc + beta*y[row]; }
+        }
+      }
+    else
+    if(do_trans_A == true)
+      {
+      if(is_cx<eT>::no)
+        {
+        #pragma omp parallel for schedule(static) num_threads(n_threads)
+        for(uword col=0; col < A_n_cols; ++col)
+          {
+          // col is interpreted as row when storing the results in 'y'
+          
+          const eT acc = op_dot::direct_dot(A_n_rows, A.colptr(col), x);
+          
+               if( (use_alpha == false) && (use_beta == false) )  { y[col] =       acc;               }
+          else if( (use_alpha == true ) && (use_beta == false) )  { y[col] = alpha*acc;               }
+          else if( (use_alpha == false) && (use_beta == true ) )  { y[col] =       acc + beta*y[col]; }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { y[col] = alpha*acc + beta*y[col]; }
+          }
+        }
+      else
+        {
+        Mat<eT> AA;
+        
+        op_htrans::apply_mat_noalias(AA, A);
+        
+        gemv_emul_mp<false, use_alpha, use_beta>::apply(y, AA, x, alpha, beta);
+        }
+      }
+    }
+  
+  };
+#endif
+
+
+
 //! \brief
 //! Partial emulation of BLAS gemv().
 //! 'y' is assumed to have been set to the correct size (ie. taking into account the transpose)
@@ -222,6 +290,21 @@ struct gemv_emul
     const uword A_n_rows = A.n_rows;
     const uword A_n_cols = A.n_cols;
     
+    #if defined(ARMA_USE_OPENMP)
+      {
+      // TODO: replace with more sophisticated threshold mechanism
+      
+      constexpr uword threshold = uword(200);
+      
+      if( (A_n_rows >= threshold) && (A_n_cols >= threshold) && (mp_thread_limit::in_parallel() == false) )
+        {
+        gemv_emul_mp<do_trans_A, use_alpha, use_beta>::apply(y, A, x, alpha, beta);
+        
+        return;
+        }
+      }
+    #endif
+    
     if(do_trans_A == false)
       {
       if(A_n_rows == 1)

include/armadillo_bits/op_vectorise_meat.hpp

@@ -365,14 +365,30 @@ op_vectorise_cube_col::apply_subview(Mat<eT>& out, const subview_cube<eT>& sv)
   
   out.set_size(sv.n_elem, 1);
   
+  if(sv.n_elem == 0)  { return; }
+  
   eT* out_mem = out.memptr();
   
-  for(uword s=0; s < sv_ns; ++s)
-  for(uword c=0; c < sv_nc; ++c)
+  if( (sv_nr == 1) && (sv_nc == 1) && (sv.aux_slice1 == 0) )
     {
-    arrayops::copy(out_mem, sv.slice_colptr(s,c), sv_nr);
+    const uword sv_m_n_elem_slice = sv.m.n_elem_slice;
     
-    out_mem += sv_nr;
+    const eT* sv_m_ptr = &( sv.m.at(sv.aux_row1, sv.aux_col1, 0) );
+    
+    for(uword s=0; s < sv_ns; ++s)
+      {
+      out_mem[s] = (*sv_m_ptr);  sv_m_ptr += sv_m_n_elem_slice;
+      }
+    }
+  else
+    {
+    for(uword s=0; s < sv_ns; ++s)
+    for(uword c=0; c < sv_nc; ++c)
+      {
+      arrayops::copy(out_mem, sv.slice_colptr(s,c), sv_nr);
+      
+      out_mem += sv_nr;
+      }
     }
   }
 

tests2/mat_mul_real.cpp

@@ -969,3 +969,134 @@ TEMPLATE_TEST_CASE("mat_mul_fp_compare", "[mat_mul]", TEST_FLOAT_TYPES)
   REQUIRE( accu(abs(diff3)) == Approx(eT(0)).margin(margin * diff3.n_elem) );
   REQUIRE( accu(abs(diff4)) == Approx(eT(0)).margin(margin * diff4.n_elem) );
   }
+
+
+
+#if defined(ARMA_USE_BLAS)
+TEMPLATE_TEST_CASE("mat_mul_int_compare", "[mat_mul]", u32, s32, u64, s64)
+  {
+  typedef TestType eT;
+  
+  const uword n_trials = 5;
+  
+  for (uword trial = 0; trial < n_trials; ++trial)
+    {
+    uword m = randi<uword>(distr_param(10, 500));
+    uword n = randi<uword>(distr_param(10, 500));
+    uword k = randi<uword>(distr_param(10, 500));
+    
+    Mat<eT> A  = randi<Mat<eT>>(m, n, distr_param(0, 100));
+    Mat<eT> At = randi<Mat<eT>>(n, m, distr_param(0, 100));
+    Mat<eT> B  = randi<Mat<eT>>(n, k, distr_param(0, 100));
+    Mat<eT> Bt = randi<Mat<eT>>(k, n, distr_param(0, 100));
+    
+    // compare against BLAS implementation
+    mat Ad  = conv_to<mat>::from(A);
+    mat Atd = conv_to<mat>::from(At);
+    mat Bd  = conv_to<mat>::from(B);
+    mat Btd = conv_to<mat>::from(Bt);
+    
+    Mat<eT> C = A * B;
+    mat Cd = Ad * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C = A * Bt.t();
+    Cd = Ad * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C = At.t() * B;
+    Cd = Atd.t() * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C = At.t() * Bt.t();
+    Cd = Atd.t() * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    // now test variants with alpha
+    
+    C = 2 * A * B;
+    Cd = 2 * Ad * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C = 2 * A * Bt.t();
+    Cd = 2 * Ad * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C = 2 * At.t() * B;
+    Cd = 2 * Atd.t() * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C = 2 * At.t() * Bt.t();
+    Cd = 2 * Atd.t() * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    // variants with beta
+    
+    C.ones(A.n_rows, B.n_cols);
+    C = 2 * C + A * B;
+    Cd.ones(Ad.n_rows, Bd.n_cols);
+    Cd = 2 * Cd + Ad * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C.ones(A.n_rows, Bt.n_rows);
+    C = 2 * C + A * Bt.t();
+    Cd.ones(Ad.n_rows, Btd.n_rows);
+    Cd = 2 * Cd + Ad * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C.ones(At.n_cols, B.n_cols);
+    C = 2 * C + At.t() * B;
+    Cd.ones(Atd.n_cols, Bd.n_cols);
+    Cd = 2 * Cd + Atd.t() * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C.ones(At.n_cols, Bt.n_rows);
+    C = 2 * C + At.t() * Bt.t();
+    Cd.ones(Atd.n_cols, Btd.n_rows);
+    Cd = 2 * Cd + Atd.t() * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    // variants with both alpha and beta
+    
+    C.ones(A.n_rows, B.n_cols);
+    C = 2 * C + 3 * A * B;
+    Cd.ones(Ad.n_rows, Bd.n_cols);
+    Cd = 2 * Cd + 3 * Ad * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C.ones(A.n_rows, Bt.n_rows);
+    C = 2 * C + 3 * A * Bt.t();
+    Cd.ones(Ad.n_rows, Btd.n_rows);
+    Cd = 2 * Cd + 3 * Ad * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C.ones(At.n_cols, B.n_cols);
+    C = 2 * C + 3 * At.t() * B;
+    Cd.ones(Atd.n_cols, Bd.n_cols);
+    Cd = 2 * Cd + 3 * Atd.t() * Bd;
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    
+    C.ones(At.n_cols, Bt.n_rows);
+    C = 2 * C + 3 * At.t() * Bt.t();
+    Cd.ones(Atd.n_cols, Btd.n_rows);
+    Cd = 2 * Cd + 3 * Atd.t() * Btd.t();
+    
+    REQUIRE( approx_equal( conv_to<mat>::from(C), Cd, "both", 1e-5, 1e-5 ) );
+    }
+  }
+#endif