patch bump

CMakeLists.txt

@@ -381,10 +381,11 @@ else()
     set(ARMA_LIBS ${ARMA_LIBS} ${MKL_LIBRARIES})
     
     message(STATUS "")
-    message(STATUS "*** If the MKL libraries are installed in non-standard locations such as")
+    message(STATUS "*** If the MKL libraries are installed in a non-standard location such as")
     message(STATUS "*** /opt/intel/mkl, /opt/intel/composerxe/, /usr/local/intel/mkl")
     message(STATUS "*** make sure the run-time linker can find them.")
     message(STATUS "*** On Linux systems this can be done by editing /etc/ld.so.conf")
+    message(STATUS "*** or specifying the location through a file in /etc/ld.so.conf.d/")
     message(STATUS "*** or modifying the LD_LIBRARY_PATH environment variable.")
     message(STATUS "")
     message(STATUS "*** On systems with SELinux enabled (eg. Fedora, RHEL),")
@@ -419,10 +420,11 @@ else()
       endif()
       
       message(STATUS "")
-      message(STATUS "*** If the OpenBLAS library is installed in")
+      message(STATUS "*** If the OpenBLAS library is installed in a location such as")
       message(STATUS "*** /usr/local/lib or /usr/local/lib64")
       message(STATUS "*** make sure the run-time linker can find it.")
       message(STATUS "*** On Linux systems this can be done by editing /etc/ld.so.conf")
+      message(STATUS "*** or specifying the location through a file in /etc/ld.so.conf.d/")
       message(STATUS "*** or modifying the LD_LIBRARY_PATH environment variable.")
       message(STATUS "")
       

README.md

@@ -38,7 +38,7 @@ Copyright 2017-2025 Data61 / CSIRO
 13. [API Stability and Version Policy](#13-api-stability-and-version-policy)
 14. [Bug Reports and Frequently Asked Questions](#14-bug-reports-and-frequently-asked-questions)
 
-15. [MEX Interface to Octave/Matlab](#15-mex-interface-to-octavematlab)
+15. [Prohibition on Use of Large Language Models (LLMs)](#15-prohibition-on-use-of-large-language-models-llms)
 16. [Related Software Using Armadillo](#16-related-software-using-armadillo)
 
 ---
@@ -391,8 +391,8 @@ Use a web browser to view it.
 
 Use of the C++11 `auto` keyword is not recommended with Armadillo objects and expressions.
 
-Armadillo has a template meta-programming framework which creates short-lived temporaries
-that are not properly handled by `auto`.
+Armadillo has a template meta-programming framework which creates short-lived
+temporary objects that are not properly handled by `auto`.
 
 ---
 
@@ -452,10 +452,12 @@ https://arma.sourceforge.net/faq.html
 
 ---
 
-### 15: MEX Interface to Octave/Matlab
+### 15: Prohibition on Use of Large Language Models (LLMs)
 
-The `mex_interface` folder contains examples of how to interface
-Octave/Matlab with C++ code that uses Armadillo matrices.
+The Armadillo codebase has been developed entirely without the use of LLMs.
+We do not accept any code which has been created (directly or indirectly)
+with the aid of LLMs. This is due to issues with code provenance and low quality
+of LLM generated code.
 
 ---
 

include/armadillo

@@ -581,6 +581,7 @@ namespace arma
   #include "armadillo_bits/fn_inplace_trans.hpp"
   #include "armadillo_bits/fn_randi.hpp"
   #include "armadillo_bits/fn_randg.hpp"
+  #include "armadillo_bits/fn_rande.hpp"
   #include "armadillo_bits/fn_cond_rcond.hpp"
   #include "armadillo_bits/fn_normalise.hpp"
   #include "armadillo_bits/fn_clamp.hpp"

include/armadillo_bits/Base_meat.hpp

@@ -667,7 +667,7 @@ Base<elem_type,derived>::is_finite() const
   
   if(arma_config::fast_math_warn)  { arma_warn(1, "is_finite(): detection of non-finite values is not reliable in fast math mode"); }
   
-  if(is_Mat<typename Proxy<derived>::stored_type>::value)
+  if( (quasi_unwrap<derived>::has_orig_mem) || (is_Mat<typename Proxy<derived>::stored_type>::value) )
     {
     const quasi_unwrap<derived> U( (*this).get_ref() );
     
@@ -715,7 +715,7 @@ Base<elem_type,derived>::has_inf() const
   
   if(arma_config::fast_math_warn)  { arma_warn(1, "has_inf(): detection of non-finite values is not reliable in fast math mode"); }
   
-  if(is_Mat<typename Proxy<derived>::stored_type>::value)
+  if( (quasi_unwrap<derived>::has_orig_mem) || (is_Mat<typename Proxy<derived>::stored_type>::value) )
     {
     const quasi_unwrap<derived> U( (*this).get_ref() );
     
@@ -763,7 +763,7 @@ Base<elem_type,derived>::has_nan() const
   
   if(arma_config::fast_math_warn)  { arma_warn(1, "has_nan(): detection of non-finite values is not reliable in fast math mode"); }
   
-  if(is_Mat<typename Proxy<derived>::stored_type>::value)
+  if( (quasi_unwrap<derived>::has_orig_mem) || (is_Mat<typename Proxy<derived>::stored_type>::value) )
     {
     const quasi_unwrap<derived> U( (*this).get_ref() );
     
@@ -811,7 +811,7 @@ Base<elem_type,derived>::has_nonfinite() const
   
   if(arma_config::fast_math_warn)  { arma_warn(1, "has_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
   
-  if(is_Mat<typename Proxy<derived>::stored_type>::value)
+  if( (quasi_unwrap<derived>::has_orig_mem) || (is_Mat<typename Proxy<derived>::stored_type>::value) )
     {
     const quasi_unwrap<derived> U( (*this).get_ref() );
     

include/armadillo_bits/Col_bones.hpp

@@ -142,7 +142,7 @@ class Col : public Mat<eT>
   template<typename T1> inline void shed_rows(const Base<uword, T1>& indices);
   
   [[deprecated]] inline void insert_rows(const uword row_num, const uword N, const bool set_to_zero);
-                  inline void insert_rows(const uword row_num, const uword N);
+                 inline void insert_rows(const uword row_num, const uword N);
   
   template<typename T1> inline void insert_rows(const uword row_num, const Base<eT,T1>& X);
   
@@ -154,6 +154,11 @@ class Col : public Mat<eT>
   arma_warn_unused arma_inline const eT& at(const uword in_row, const uword in_col) const;
   
   
+  inline constexpr bool is_vec()    const { return true;  }
+  inline constexpr bool is_rowvec() const { return false; }
+  inline constexpr bool is_colvec() const { return true;  }
+  
+  
   typedef       eT*       row_iterator;
   typedef const eT* const_row_iterator;
   

include/armadillo_bits/Col_meat.hpp

@@ -483,7 +483,7 @@ Col<eT>::operator=(Col<eT>&& X)
 //   if(X.n_cols != 1)  { const Mat<eT>& XX = X; Mat<eT>::operator=(XX); return *this; }
 //   
 //   (*this).steal_mem(X, true);
-//   
+//    
 //   return *this;
 //   }
 
@@ -521,11 +521,9 @@ template<typename eT>
 template<typename T1>
 inline
 Col<eT>::Col(const Base<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 1)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  Mat<eT>::operator=(X.get_ref());
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -549,11 +547,9 @@ template<typename eT>
 template<typename T1>
 inline
 Col<eT>::Col(const SpBase<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 1)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 1)
   {
   arma_debug_sigprint_this(this);
-  
-  Mat<eT>::operator=(X.get_ref());
   }
 
 
@@ -607,12 +603,9 @@ Col<eT>::Col
   const Base<typename Col<eT>::pod_type, T1>& A,
   const Base<typename Col<eT>::pod_type, T2>& B
   )
+  : Mat<eT>(A.get_ref(), B.get_ref(), arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 1;
-  
-  Mat<eT>::init(A,B);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -621,12 +614,9 @@ template<typename eT>
 template<typename T1>
 inline
 Col<eT>::Col(const BaseCube<eT,T1>& X)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 1;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -649,12 +639,9 @@ Col<eT>::operator=(const BaseCube<eT,T1>& X)
 template<typename eT>
 inline
 Col<eT>::Col(const subview_cube<eT>& X)
+  : Mat<eT>(X, arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 1;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -1277,6 +1264,13 @@ Col<eT>::end_row(const uword row_num) const
 
 
 
+
+
+
+//
+
+
+
 template<typename eT>
 template<uword fixed_n_elem>
 arma_inline
@@ -1350,6 +1344,10 @@ Col<eT>::fixed<fixed_n_elem>::fixed(const fill::fill_class<fill_type>&)
   if(is_same_type<fill_type, fill::fill_eye  >::yes)  { Mat<eT>::eye();   }
   if(is_same_type<fill_type, fill::fill_randu>::yes)  { Mat<eT>::randu(); }
   if(is_same_type<fill_type, fill::fill_randn>::yes)  { Mat<eT>::randn(); }
+  
+  if(is_same_type<fill_type, fill::fill_nan    >::yes)  { (*this).fill( priv::Datum_helper::nan    <eT>() ); }
+  if(is_same_type<fill_type, fill::fill_pos_inf>::yes)  { (*this).fill( priv::Datum_helper::pos_inf<eT>() ); }
+  if(is_same_type<fill_type, fill::fill_neg_inf>::yes)  { (*this).fill( priv::Datum_helper::neg_inf<eT>() ); }
   }
 
 
@@ -1872,6 +1870,10 @@ Col<eT>::fixed<fixed_n_elem>::ones()
 
 
 
+//
+
+
+
 template<typename eT>
 inline
 Col<eT>::Col(const arma_fixed_indicator&, const uword in_n_elem, const eT* in_mem)

include/armadillo_bits/Cube_bones.hpp

@@ -360,6 +360,9 @@ class Cube : public BaseCube< eT, Cube<eT> >
   
   inline Cube& fill(const eT val);
   
+  template<typename fill_type>
+  inline Cube& fill(const fill::fill_class<fill_type>& f);
+  
   inline Cube& zeros();
   inline Cube& zeros(const uword new_n_rows, const uword new_n_cols, const uword new_n_slices);
   inline Cube& zeros(const SizeCube& s);

include/armadillo_bits/Cube_meat.hpp

@@ -52,7 +52,7 @@ Cube<eT>::Cube()
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
   }
@@ -172,7 +172,7 @@ Cube<eT>::Cube(const SizeCube& s, const arma_initmode_indicator<do_zeros>&)
 template<typename eT>
 template<typename fill_type>
 inline
-Cube<eT>::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_slices, const fill::fill_class<fill_type>&)
+Cube<eT>::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_slices, const fill::fill_class<fill_type>& f)
   : n_rows(in_n_rows)
   , n_cols(in_n_cols)
   , n_elem_slice(in_n_rows*in_n_cols)
@@ -186,12 +186,7 @@ Cube<eT>::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_sl
   
   init_cold();
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
-  
-  arma_static_check( (is_same_type<fill_type, fill::fill_eye>::yes), "Cube::Cube(): unsupported fill type" );
+  (*this).fill(f);
   }
 
 
@@ -199,7 +194,7 @@ Cube<eT>::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_sl
 template<typename eT>
 template<typename fill_type>
 inline
-Cube<eT>::Cube(const SizeCube& s, const fill::fill_class<fill_type>&)
+Cube<eT>::Cube(const SizeCube& s, const fill::fill_class<fill_type>& f)
   : n_rows(s.n_rows)
   , n_cols(s.n_cols)
   , n_elem_slice(s.n_rows*s.n_cols)
@@ -213,12 +208,7 @@ Cube<eT>::Cube(const SizeCube& s, const fill::fill_class<fill_type>&)
   
   init_cold();
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
-  
-  arma_static_check( (is_same_type<fill_type, fill::fill_eye>::yes), "Cube::Cube(): unsupported fill type" );
+  (*this).fill(f);
   }
 
 
@@ -934,7 +924,7 @@ Cube<eT>::Cube
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -1064,7 +1054,7 @@ Cube<eT>::Cube(const subview_cube_slices<eT,T1>& X)
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -2751,7 +2741,7 @@ Cube<eT>::Cube(const OpCube<T1, op_type>& X)
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
 
@@ -3048,7 +3038,7 @@ Cube<eT>::Cube(const mtOpCube<eT, T1, op_type>& X)
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -3144,11 +3134,11 @@ Cube<eT>::Cube(const GlueCube<T1, T2, glue_type>& X)
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
   
-  this->operator=(X);
+  (*this).operator=(X);
   }
 
 
@@ -3401,7 +3391,7 @@ Cube<eT>::Cube(const mtGlueCube<eT, T1, T2, glue_type>& X)
   , n_elem(0)
   , n_alloc(0)
   , mem_state(0)
-  , mem()
+  , mem(nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -3887,7 +3877,7 @@ arma_inline
 eT*
 Cube<eT>::slice_memptr(const uword uslice)
   {
-  return const_cast<eT*>( &mem[ uslice*n_elem_slice ] );
+  return access::rwp( mem + (uslice*n_elem_slice) );
   }
 
 
@@ -3898,7 +3888,7 @@ arma_inline
 const eT*
 Cube<eT>::slice_memptr(const uword uslice) const
   {
-  return &mem[ uslice*n_elem_slice ];
+  return mem + (uslice*n_elem_slice);
   }
 
 
@@ -3909,7 +3899,7 @@ arma_inline
 eT*
 Cube<eT>::slice_colptr(const uword uslice, const uword col)
   {
-  return const_cast<eT*>( &mem[ uslice*n_elem_slice + col*n_rows] );
+  return access::rwp( mem + (uslice*n_elem_slice + col*n_rows) );
   }
 
 
@@ -3920,7 +3910,7 @@ arma_inline
 const eT*
 Cube<eT>::slice_colptr(const uword uslice, const uword col) const
   {
-  return &mem[ uslice*n_elem_slice + col*n_rows ];
+  return mem + (uslice*n_elem_slice + col*n_rows);
   }
 
 
@@ -4227,6 +4217,30 @@ Cube<eT>::fill(const eT val)
 
 
 
+template<typename eT>
+template<typename fill_type>
+inline
+Cube<eT>&
+Cube<eT>::fill(const fill::fill_class<fill_type>&)
+  {
+  arma_debug_sigprint();
+  
+  arma_static_check( (is_same_type<fill_type, fill::fill_eye>::yes), "Cube::fill(): unsupported fill type" );
+  
+  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
+  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
+  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
+  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
+  
+  if(is_same_type<fill_type, fill::fill_nan    >::yes)  { (*this).fill( priv::Datum_helper::nan    <eT>() ); }
+  if(is_same_type<fill_type, fill::fill_pos_inf>::yes)  { (*this).fill( priv::Datum_helper::pos_inf<eT>() ); }
+  if(is_same_type<fill_type, fill::fill_neg_inf>::yes)  { (*this).fill( priv::Datum_helper::neg_inf<eT>() ); }
+  
+  return *this;
+  }
+
+
+
 template<typename eT>
 inline
 Cube<eT>&
@@ -5458,18 +5472,13 @@ template<typename eT>
 template<uword fixed_n_rows, uword fixed_n_cols, uword fixed_n_slices>
 template<typename fill_type>
 inline
-Cube<eT>::fixed<fixed_n_rows, fixed_n_cols, fixed_n_slices>::fixed(const fill::fill_class<fill_type>&)
+Cube<eT>::fixed<fixed_n_rows, fixed_n_cols, fixed_n_slices>::fixed(const fill::fill_class<fill_type>& f)
   {
   arma_debug_sigprint_this(this);
   
   mem_setup();
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { Cube<eT>::zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { Cube<eT>::ones();  }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { Cube<eT>::randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { Cube<eT>::randn(); }
-  
-  arma_static_check( (is_same_type<fill_type, fill::fill_eye>::yes), "Cube::fixed::fixed(): unsupported fill type" );
+  (*this).fill(f);
   }
 
 

include/armadillo_bits/Mat_bones.hpp

@@ -99,6 +99,7 @@ class Mat : public Base< eT, Mat<eT> >
   inline Mat& operator*=(const eT val);
   inline Mat& operator/=(const eT val);
   
+  inline             Mat(const Mat& m, const arma_vec_indicator&, const uhword in_vec_state);
   inline             Mat(const Mat& m);
   inline Mat& operator= (const Mat& m);
   inline Mat& operator+=(const Mat& m);
@@ -107,6 +108,7 @@ class Mat : public Base< eT, Mat<eT> >
   inline Mat& operator%=(const Mat& m);
   inline Mat& operator/=(const Mat& m);
   
+  template<typename T1> inline             Mat(const BaseCube<eT,T1>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1> inline             Mat(const BaseCube<eT,T1>& X);
   template<typename T1> inline Mat& operator= (const BaseCube<eT,T1>& X);
   template<typename T1> inline Mat& operator+=(const BaseCube<eT,T1>& X);
@@ -115,11 +117,12 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1> inline Mat& operator%=(const BaseCube<eT,T1>& X);
   template<typename T1> inline Mat& operator/=(const BaseCube<eT,T1>& X);
   
-  template<typename T1, typename T2>
-  inline explicit Mat(const Base<pod_type,T1>& A, const Base<pod_type,T2>& B);
+  template<typename T1, typename T2> inline explicit Mat(const Base<pod_type,T1>& A, const Base<pod_type,T2>& B, const arma_vec_indicator&, const uhword in_vec_state);
+  template<typename T1, typename T2> inline explicit Mat(const Base<pod_type,T1>& A, const Base<pod_type,T2>& B);
   
-  inline explicit    Mat(const subview<eT>& X, const bool use_colmem);  // only to be used by the quasi_unwrap class
+  inline explicit    Mat(const subview<eT>& X, const bool reuse_mem);  // only to be used by quasi_unwrap and partial_unwrap
   
+  inline             Mat(const subview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state);
   inline             Mat(const subview<eT>& X);
   inline Mat& operator= (const subview<eT>& X);
   inline Mat& operator+=(const subview<eT>& X);
@@ -135,6 +138,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<bool do_conj>
   inline Mat(const xtrans_mat<eT,do_conj>& X);  //        xtrans_mat can only be generated by the Proxy class
   
+  inline             Mat(const subview_cube<eT>& X, const arma_vec_indicator&, const uhword in_vec_state);
   inline             Mat(const subview_cube<eT>& X);
   inline Mat& operator= (const subview_cube<eT>& X);
   inline Mat& operator+=(const subview_cube<eT>& X);
@@ -143,6 +147,7 @@ class Mat : public Base< eT, Mat<eT> >
   inline Mat& operator%=(const subview_cube<eT>& X);
   inline Mat& operator/=(const subview_cube<eT>& X);
   
+  inline             Mat(const diagview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state);
   inline             Mat(const diagview<eT>& X);
   inline Mat& operator= (const diagview<eT>& X);
   inline Mat& operator+=(const diagview<eT>& X);
@@ -151,6 +156,7 @@ class Mat : public Base< eT, Mat<eT> >
   inline Mat& operator%=(const diagview<eT>& X);
   inline Mat& operator/=(const diagview<eT>& X);
   
+  template<typename T1> inline             Mat(const subview_elem1<eT,T1>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1> inline             Mat(const subview_elem1<eT,T1>& X);
   template<typename T1> inline Mat& operator= (const subview_elem1<eT,T1>& X);
   template<typename T1> inline Mat& operator+=(const subview_elem1<eT,T1>& X);
@@ -159,6 +165,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1> inline Mat& operator%=(const subview_elem1<eT,T1>& X);
   template<typename T1> inline Mat& operator/=(const subview_elem1<eT,T1>& X);
   
+  template<typename T1, typename T2> inline             Mat(const subview_elem2<eT,T1,T2>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename T2> inline             Mat(const subview_elem2<eT,T1,T2>& X);
   template<typename T1, typename T2> inline Mat& operator= (const subview_elem2<eT,T1,T2>& X);
   template<typename T1, typename T2> inline Mat& operator+=(const subview_elem2<eT,T1,T2>& X);
@@ -168,6 +175,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename T2> inline Mat& operator/=(const subview_elem2<eT,T1,T2>& X);
 
   // Operators on sparse matrices (and subviews)
+  template<typename T1> inline explicit    Mat(const SpBase<eT, T1>& m, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1> inline explicit    Mat(const SpBase<eT, T1>& m);
   template<typename T1> inline Mat& operator= (const SpBase<eT, T1>& m);
   template<typename T1> inline Mat& operator+=(const SpBase<eT, T1>& m);
@@ -176,11 +184,13 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1> inline Mat& operator%=(const SpBase<eT, T1>& m);
   template<typename T1> inline Mat& operator/=(const SpBase<eT, T1>& m);
   
+  inline explicit    Mat(const SpSubview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state);
   inline explicit    Mat(const SpSubview<eT>& X);
   inline Mat& operator= (const SpSubview<eT>& X);
   inline Mat& operator+=(const SpSubview<eT>& X);
   inline Mat& operator-=(const SpSubview<eT>& X);
   
+  inline explicit    Mat(const spdiagview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state);
   inline explicit    Mat(const spdiagview<eT>& X);
   inline Mat& operator= (const spdiagview<eT>& X);
   inline Mat& operator+=(const spdiagview<eT>& X);
@@ -320,6 +330,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1> inline void insert_cols(const uword col_num, const Base<eT,T1>& X);
   
   
+  template<typename T1, typename gen_type> inline             Mat(const Gen<T1, gen_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename gen_type> inline             Mat(const Gen<T1, gen_type>& X);
   template<typename T1, typename gen_type> inline Mat& operator= (const Gen<T1, gen_type>& X);
   template<typename T1, typename gen_type> inline Mat& operator+=(const Gen<T1, gen_type>& X);
@@ -328,6 +339,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename gen_type> inline Mat& operator%=(const Gen<T1, gen_type>& X);
   template<typename T1, typename gen_type> inline Mat& operator/=(const Gen<T1, gen_type>& X);
   
+  template<typename T1, typename op_type> inline             Mat(const Op<T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename op_type> inline             Mat(const Op<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator= (const Op<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator+=(const Op<T1, op_type>& X);
@@ -336,6 +348,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename op_type> inline Mat& operator%=(const Op<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator/=(const Op<T1, op_type>& X);
   
+  template<typename T1, typename eop_type> inline             Mat(const eOp<T1, eop_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename eop_type> inline             Mat(const eOp<T1, eop_type>& X);
   template<typename T1, typename eop_type> inline Mat& operator= (const eOp<T1, eop_type>& X);
   template<typename T1, typename eop_type> inline Mat& operator+=(const eOp<T1, eop_type>& X);
@@ -344,6 +357,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename eop_type> inline Mat& operator%=(const eOp<T1, eop_type>& X);
   template<typename T1, typename eop_type> inline Mat& operator/=(const eOp<T1, eop_type>& X);
   
+  template<typename T1, typename op_type> inline             Mat(const mtOp<eT, T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename op_type> inline             Mat(const mtOp<eT, T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator= (const mtOp<eT, T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator+=(const mtOp<eT, T1, op_type>& X);
@@ -352,6 +366,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename op_type> inline Mat& operator%=(const mtOp<eT, T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator/=(const mtOp<eT, T1, op_type>& X);
 
+  template<typename T1, typename op_type> inline             Mat(const CubeToMatOp<T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename op_type> inline             Mat(const CubeToMatOp<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator= (const CubeToMatOp<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator+=(const CubeToMatOp<T1, op_type>& X);
@@ -360,6 +375,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename op_type> inline Mat& operator%=(const CubeToMatOp<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator/=(const CubeToMatOp<T1, op_type>& X);
   
+  template<typename T1, typename op_type> inline             Mat(const SpToDOp<T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename op_type> inline             Mat(const SpToDOp<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator= (const SpToDOp<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator+=(const SpToDOp<T1, op_type>& X);
@@ -368,6 +384,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename op_type> inline Mat& operator%=(const SpToDOp<T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator/=(const SpToDOp<T1, op_type>& X);
   
+  template<typename T1, typename op_type> inline explicit    Mat(const mtSpReduceOp<eT, T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename op_type> inline explicit    Mat(const mtSpReduceOp<eT, T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator= (const mtSpReduceOp<eT, T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator+=(const mtSpReduceOp<eT, T1, op_type>& X);
@@ -376,6 +393,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename op_type> inline Mat& operator%=(const mtSpReduceOp<eT, T1, op_type>& X);
   template<typename T1, typename op_type> inline Mat& operator/=(const mtSpReduceOp<eT, T1, op_type>& X);
   
+  template<typename T1, typename T2, typename glue_type> inline             Mat(const Glue<T1, T2, glue_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename T2, typename glue_type> inline             Mat(const Glue<T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator= (const Glue<T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator+=(const Glue<T1, T2, glue_type>& X);
@@ -387,6 +405,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename T2>                     inline Mat& operator+=(const Glue<T1, T2, glue_times>& X);
   template<typename T1, typename T2>                     inline Mat& operator-=(const Glue<T1, T2, glue_times>& X);
   
+  template<typename T1, typename T2, typename eglue_type> inline             Mat(const eGlue<T1, T2, eglue_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename T2, typename eglue_type> inline             Mat(const eGlue<T1, T2, eglue_type>& X);
   template<typename T1, typename T2, typename eglue_type> inline Mat& operator= (const eGlue<T1, T2, eglue_type>& X);
   template<typename T1, typename T2, typename eglue_type> inline Mat& operator+=(const eGlue<T1, T2, eglue_type>& X);
@@ -395,6 +414,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename T2, typename eglue_type> inline Mat& operator%=(const eGlue<T1, T2, eglue_type>& X);
   template<typename T1, typename T2, typename eglue_type> inline Mat& operator/=(const eGlue<T1, T2, eglue_type>& X);
   
+  template<typename T1, typename T2, typename glue_type> inline             Mat(const mtGlue<eT, T1, T2, glue_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename T2, typename glue_type> inline             Mat(const mtGlue<eT, T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator= (const mtGlue<eT, T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator+=(const mtGlue<eT, T1, T2, glue_type>& X);
@@ -403,6 +423,7 @@ class Mat : public Base< eT, Mat<eT> >
   template<typename T1, typename T2, typename glue_type> inline Mat& operator%=(const mtGlue<eT, T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator/=(const mtGlue<eT, T1, T2, glue_type>& X);
   
+  template<typename T1, typename T2, typename glue_type> inline             Mat(const SpToDGlue<T1, T2, glue_type>& X, const arma_vec_indicator&, const uhword in_vec_state);
   template<typename T1, typename T2, typename glue_type> inline             Mat(const SpToDGlue<T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator= (const SpToDGlue<T1, T2, glue_type>& X);
   template<typename T1, typename T2, typename glue_type> inline Mat& operator+=(const SpToDGlue<T1, T2, glue_type>& X);
@@ -958,4 +979,8 @@ class Mat_aux
 
 
 
+template<typename eT> class Mat_noalias : public Mat<eT> {};
+
+
+
 //! @}

include/armadillo_bits/Row_bones.hpp

@@ -142,7 +142,7 @@ class Row : public Mat<eT>
   template<typename T1> inline void shed_cols(const Base<uword, T1>& indices);
   
   [[deprecated]] inline void insert_cols(const uword col_num, const uword N, const bool set_to_zero);
-                  inline void insert_cols(const uword col_num, const uword N);
+                 inline void insert_cols(const uword col_num, const uword N);
   
   template<typename T1> inline void insert_cols(const uword col_num, const Base<eT,T1>& X);
   
@@ -154,6 +154,11 @@ class Row : public Mat<eT>
   arma_warn_unused arma_inline const eT& at(const uword in_row, const uword in_col) const;
   
   
+  inline constexpr bool is_vec()    const { return true;  }
+  inline constexpr bool is_rowvec() const { return true;  }
+  inline constexpr bool is_colvec() const { return false; }
+  
+  
   typedef       eT*       row_iterator;
   typedef const eT* const_row_iterator;
   
@@ -164,6 +169,9 @@ class Row : public Mat<eT>
   inline const_row_iterator end_row  (const uword row_num) const;
   
   
+  inline explicit Row(const subview<eT>& X, const bool reuse_mem);  // only to be used by the partial_unwrap class
+  
+  
   template<uword fixed_n_elem> class fixed;
   
   

include/armadillo_bits/Row_meat.hpp

@@ -521,11 +521,9 @@ template<typename eT>
 template<typename T1>
 inline
 Row<eT>::Row(const Base<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 2)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  Mat<eT>::operator=(X.get_ref());
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -549,11 +547,9 @@ template<typename eT>
 template<typename T1>
 inline
 Row<eT>::Row(const SpBase<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 2)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  Mat<eT>::operator=(X.get_ref());
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -607,12 +603,9 @@ Row<eT>::Row
   const Base<typename Row<eT>::pod_type, T1>& A,
   const Base<typename Row<eT>::pod_type, T2>& B
   )
+  : Mat<eT>(A.get_ref(), B.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 2;
-  
-  Mat<eT>::init(A,B);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -621,12 +614,9 @@ template<typename eT>
 template<typename T1>
 inline
 Row<eT>::Row(const BaseCube<eT,T1>& X)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 2;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -649,12 +639,9 @@ Row<eT>::operator=(const BaseCube<eT,T1>& X)
 template<typename eT>
 inline
 Row<eT>::Row(const subview_cube<eT>& X)
+  : Mat<eT>(X, arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 2;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -1277,6 +1264,20 @@ Row<eT>::end_row(const uword row_num) const
 
 
 
+template<typename eT>
+inline
+Row<eT>::Row(const subview<eT>& X, const bool reuse_mem)
+  : Mat<eT>(X, reuse_mem)
+  {
+  arma_debug_sigprint_this(this);
+  }
+
+
+
+//
+
+
+
 template<typename eT>
 template<uword fixed_n_elem>
 arma_inline
@@ -1350,6 +1351,10 @@ Row<eT>::fixed<fixed_n_elem>::fixed(const fill::fill_class<fill_type>&)
   if(is_same_type<fill_type, fill::fill_eye  >::yes)  { Mat<eT>::eye();   }
   if(is_same_type<fill_type, fill::fill_randu>::yes)  { Mat<eT>::randu(); }
   if(is_same_type<fill_type, fill::fill_randn>::yes)  { Mat<eT>::randn(); }
+  
+  if(is_same_type<fill_type, fill::fill_nan    >::yes)  { (*this).fill( priv::Datum_helper::nan    <eT>() ); }
+  if(is_same_type<fill_type, fill::fill_pos_inf>::yes)  { (*this).fill( priv::Datum_helper::pos_inf<eT>() ); }
+  if(is_same_type<fill_type, fill::fill_neg_inf>::yes)  { (*this).fill( priv::Datum_helper::neg_inf<eT>() ); }
   }
 
 
@@ -1872,6 +1877,10 @@ Row<eT>::fixed<fixed_n_elem>::ones()
 
 
 
+//
+
+
+
 template<typename eT>
 inline
 Row<eT>::Row(const arma_fixed_indicator&, const uword in_n_elem, const eT* in_mem)

include/armadillo_bits/SpMat_bones.hpp

@@ -741,6 +741,10 @@ class SpMat_aux
 
 
 
+template<typename eT> class SpMat_noalias : public SpMat<eT> {};
+
+
+
 #define ARMA_HAS_SPMAT
 
 

include/armadillo_bits/SpMat_meat.hpp

@@ -688,15 +688,22 @@ SpMat<eT>::operator=(const SpMat<eT>& x)
 template<typename eT>
 inline
 SpMat<eT>&
-SpMat<eT>::operator+=(const SpMat<eT>& x)
+SpMat<eT>::operator+=(const SpMat<eT>& X)
   {
   arma_debug_sigprint();
   
   sync_csc();
   
-  SpMat<eT> out = (*this) + x;
-  
-  steal_mem(out);
+  if(X.n_nonzero == 0)
+    {
+    arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "addition");
+    }
+  else
+    {
+    SpMat<eT> tmp = (*this) + X;
+    
+    steal_mem(tmp);
+    }
   
   return *this;
   }
@@ -706,15 +713,22 @@ SpMat<eT>::operator+=(const SpMat<eT>& x)
 template<typename eT>
 inline
 SpMat<eT>&
-SpMat<eT>::operator-=(const SpMat<eT>& x)
+SpMat<eT>::operator-=(const SpMat<eT>& X)
   {
   arma_debug_sigprint();
   
   sync_csc();
   
-  SpMat<eT> out = (*this) - x;
-  
-  steal_mem(out);
+  if(X.n_nonzero == 0)
+    {
+    arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "subtraction");
+    }
+  else
+    {
+    SpMat<eT> tmp = (*this) - X;
+    
+    steal_mem(tmp);
+    }
   
   return *this;
   }
@@ -724,15 +738,15 @@ SpMat<eT>::operator-=(const SpMat<eT>& x)
 template<typename eT>
 inline
 SpMat<eT>&
-SpMat<eT>::operator*=(const SpMat<eT>& y)
+SpMat<eT>::operator*=(const SpMat<eT>& X)
   {
   arma_debug_sigprint();
   
   sync_csc();
   
-  SpMat<eT> z = (*this) * y;
+  SpMat<eT> tmp = (*this) * X;
   
-  steal_mem(z);
+  steal_mem(tmp);
   
   return *this;
   }
@@ -743,15 +757,24 @@ SpMat<eT>::operator*=(const SpMat<eT>& y)
 template<typename eT>
 inline
 SpMat<eT>&
-SpMat<eT>::operator%=(const SpMat<eT>& y)
+SpMat<eT>::operator%=(const SpMat<eT>& X)
   {
   arma_debug_sigprint();
   
   sync_csc();
   
-  SpMat<eT> z = (*this) % y;
-  
-  steal_mem(z);
+  if(X.n_nonzero == 0)
+    {
+    arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "element-wise multiplication");
+    
+    (*this).zeros();
+    }
+  else
+    {
+    SpMat<eT> tmp = (*this) % X;
+    
+    steal_mem(tmp);
+    }
   
   return *this;
   }
@@ -1338,9 +1361,16 @@ SpMat<eT>::operator+=(const SpSubview<eT>& X)
   
   sync_csc();
   
-  SpMat<eT> tmp = (*this) + X;
-  
-  steal_mem(tmp);
+  if(X.n_nonzero == 0)
+    {
+    arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "addition");
+    }
+  else
+    {
+    SpMat<eT> tmp = (*this) + X;
+    
+    steal_mem(tmp);
+    }
   
   return *this;
   }
@@ -1356,64 +1386,82 @@ SpMat<eT>::operator-=(const SpSubview<eT>& X)
   
   sync_csc();
   
-  SpMat<eT> tmp = (*this) - X;
-  
-  steal_mem(tmp);
-  
-  return *this;
-  }
-
-
-
-template<typename eT>
-inline
-SpMat<eT>&
-SpMat<eT>::operator*=(const SpSubview<eT>& y)
-  {
-  arma_debug_sigprint();
-  
-  sync_csc();
-  
-  SpMat<eT> z = (*this) * y;
-  
-  steal_mem(z);
-  
-  return *this;
-  }
-
-
-
-template<typename eT>
-inline
-SpMat<eT>&
-SpMat<eT>::operator%=(const SpSubview<eT>& x)
-  {
-  arma_debug_sigprint();
-  
-  sync_csc();
-  
-  SpMat<eT> tmp = (*this) % x;
-  
-  steal_mem(tmp);
-  
-  return *this;
-  }
-
-
-
-template<typename eT>
-inline
-SpMat<eT>&
-SpMat<eT>::operator/=(const SpSubview<eT>& x)
-  {
-  arma_debug_sigprint();
-  
-  arma_conform_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise division");
-  
-  // There is no pretty way to do this.
-  for(uword elem = 0; elem < n_elem; elem++)
+  if(X.n_nonzero == 0)
     {
-    at(elem) /= x(elem);
+    arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "subtraction");
+    }
+  else
+    {
+    SpMat<eT> tmp = (*this) - X;
+    
+    steal_mem(tmp);
+    }
+  
+  return *this;
+  }
+
+
+
+template<typename eT>
+inline
+SpMat<eT>&
+SpMat<eT>::operator*=(const SpSubview<eT>& X)
+  {
+  arma_debug_sigprint();
+  
+  sync_csc();
+  
+  SpMat<eT> tmp = (*this) * X;
+  
+  steal_mem(tmp);
+  
+  return *this;
+  }
+
+
+
+template<typename eT>
+inline
+SpMat<eT>&
+SpMat<eT>::operator%=(const SpSubview<eT>& X)
+  {
+  arma_debug_sigprint();
+  
+  sync_csc();
+  
+  if(X.n_nonzero == 0)
+    {
+    arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "element-wise multiplication");
+    
+    (*this).zeros();
+    }
+  else
+    {
+    SpMat<eT> tmp = (*this) % X;
+    
+    steal_mem(tmp);
+    }
+  
+  return *this;
+  }
+
+
+
+template<typename eT>
+inline
+SpMat<eT>&
+SpMat<eT>::operator/=(const SpSubview<eT>& X)
+  {
+  arma_debug_sigprint();
+  
+  // NOTE: use of this function is not advised; it is implemented only for completeness
+  
+  arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "element-wise division");
+  
+  for(uword c = 0; c < n_cols; ++c)
+  for(uword r = 0; r < n_rows; ++r)
+    {
+    at(r, c) /= X.at(r, c);
     }
   
   return *this;
@@ -1507,9 +1555,9 @@ SpMat<eT>::operator*=(const SpSubview_col_list<eT,T1>& X)
   
   sync_csc();
   
-  SpMat<eT> z = (*this) * X;
+  SpMat<eT> tmp = (*this) * X;
   
-  steal_mem(z);
+  steal_mem(tmp);
   
   return *this;
   }
@@ -1666,7 +1714,7 @@ SpMat<eT>::SpMat(const SpOp<T1, spop_type>& X)
   
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
-  spop_type::apply(*this, X);
+  spop_type::apply(static_cast< SpMat_noalias<eT>& >(*this), X);
   
   sync_csc();          // in case apply() used element accessors
   invalidate_cache();  // in case apply() modified the CSC representation
@@ -1806,7 +1854,7 @@ SpMat<eT>::SpMat(const SpGlue<T1, T2, spglue_type>& X)
   
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
-  spglue_type::apply(*this, X);
+  spglue_type::apply(static_cast< SpMat_noalias<eT>& >(*this), X);
   
   sync_csc();          // in case apply() used element accessors
   invalidate_cache();  // in case apply() modified the CSC representation

include/armadillo_bits/SpSubview_meat.hpp

@@ -44,6 +44,8 @@ SpSubview<eT>::SpSubview(const SpMat<eT>& in_m, const uword in_row1, const uword
   
   m.sync_csc();
   
+  if( (n_elem == 0) || (m.n_nonzero == 0) )  { return; }   // (*this).n_nonzero already set to zero
+  
   // count the number of non-zeros in the subview
   uword count = 0;
   
@@ -127,6 +129,8 @@ SpSubview<eT>::operator+=(const eT val)
   
   tmp.fill(val);
   
+  if(n_nonzero == 0)  { return (*this).operator=(tmp); }
+  
   return (*this).operator=( (*this) + tmp );
   }
 
@@ -319,6 +323,15 @@ SpSubview<eT>::operator+=(const Base<eT, T1>& x)
   {
   arma_debug_sigprint();
   
+  if(n_nonzero == 0)
+    {
+    const quasi_unwrap<T1> U(x.get_ref());
+    
+    arma_conform_assert_same_size(n_rows, n_cols, U.M.n_rows, U.M.n_cols, "addition");
+    
+    return (*this).operator=(U.M);
+    }
+  
   return (*this).operator=( (*this) + x.get_ref() );
   }
 
@@ -369,6 +382,8 @@ SpSubview<eT>::operator%=(const Base<eT, T1>& x)
   
   arma_conform_assert_same_size(sv.n_rows, sv.n_cols, B.n_rows, B.n_cols, "element-wise multiplication");
   
+  if(n_nonzero == 0)  { return *this; }
+  
   SpMat<eT>& sv_m = access::rw(sv.m);
   
   sv_m.sync_csc();
@@ -555,6 +570,15 @@ SpSubview<eT>::operator+=(const SpBase<eT, T1>& x)
   {
   arma_debug_sigprint();
   
+  if(n_nonzero == 0)
+    {
+    const unwrap_spmat<T1> U(x.get_ref());
+    
+    arma_conform_assert_same_size(n_rows, n_cols, U.M.n_rows, U.M.n_cols, "addition");
+    
+    return (*this).operator_equ_common(U.M);
+    }
+  
   // TODO: implement dedicated machinery
   return (*this).operator=( (*this) + x.get_ref() );
   }
@@ -596,6 +620,15 @@ SpSubview<eT>::operator%=(const SpBase<eT, T1>& x)
   {
   arma_debug_sigprint();
   
+  if(n_nonzero == 0)
+    {
+    const SpProxy<T1> P(x.get_ref());
+    
+    arma_conform_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols(), "element-wise multiplication");
+    
+    return *this;
+    }
+  
   // TODO: implement dedicated machinery
   return (*this).operator=( (*this) % x.get_ref() );
   }
@@ -648,6 +681,8 @@ SpSubview<eT>::for_each(functor F)
   m.sync_csc();
   m.invalidate_cache();
   
+  if(n_nonzero == 0)  { return; }
+  
   const uword lstart_row = aux_row1;
   const uword lend_row   = aux_row1 + n_rows;
   
@@ -704,6 +739,8 @@ SpSubview<eT>::for_each(functor F) const
   
   m.sync_csc();
   
+  if(n_nonzero == 0)  { return; }
+  
   const uword lstart_row = aux_row1;
   const uword lend_row   = aux_row1 + n_rows;
   
@@ -743,6 +780,8 @@ SpSubview<eT>::transform(functor F)
   m.sync_csc();
   m.invalidate_cache();
   
+  if(n_nonzero == 0)  { return; }
+  
   const uword lstart_row = aux_row1;
   const uword lend_row   = aux_row1 + n_rows;
   
@@ -813,6 +852,8 @@ SpSubview<eT>::replace(const eT old_val, const eT new_val)
   m.sync_csc();
   m.invalidate_cache();
   
+  if(n_nonzero == 0)  { return; }
+  
   const uword lstart_row = aux_row1;
   const uword lend_row   = aux_row1 + n_rows;
   

include/armadillo_bits/arma_forward.hpp

@@ -72,6 +72,10 @@ template<typename eT, typename T1> class subview_cube_slices;
 template<typename eT, typename T1> class SpSubview_col_list;
 
 
+template<typename eT> class   Mat_noalias;
+template<typename eT> class SpMat_noalias;
+
+
 struct SizeMat;
 struct SizeCube;
 

include/armadillo_bits/arma_rng.hpp

@@ -26,6 +26,9 @@
 #undef  ARMA_USE_THREAD_LOCAL
 #define ARMA_USE_THREAD_LOCAL
 
+#undef  ARMA_USE_THREAD_UNIQUE_RNG_SEED
+#define ARMA_USE_THREAD_UNIQUE_RNG_SEED
+
 #if (defined(ARMA_RNG_ALT) || defined(ARMA_DONT_USE_CXX11_RNG))
   #undef ARMA_USE_CXX11_RNG
 #endif
@@ -34,6 +37,10 @@
   #undef ARMA_USE_THREAD_LOCAL
 #endif
 
+#if defined(ARMA_DONT_USE_THREAD_UNIQUE_RNG_SEED)
+  #undef ARMA_USE_THREAD_UNIQUE_RNG_SEED
+#endif
+
 
 // NOTE: ARMA_WARMUP_PRODUCER enables a workaround 
 // NOTE: for thread_local issue on macOS 11 and/or AppleClang 12.0
@@ -116,6 +123,7 @@ struct arma_rng
   template<typename eT> struct randu;
   template<typename eT> struct randn;
   template<typename eT> struct randg;
+  template<typename eT> struct rande;
   };
 
 
@@ -128,24 +136,42 @@ arma_rng::get_producer()
   {
   #if defined(ARMA_USE_THREAD_LOCAL)
     
-    // use a thread-safe RNG, with each thread having its own unique starting seed
+    // thread-safe RNG
     
-    static std::atomic<std::size_t> mt19937_64_producer_counter(0);
-    
-    static thread_local std::mt19937_64 mt19937_64_producer( std::mt19937_64::default_seed + mt19937_64_producer_counter++ );
-    
-    arma_rng::warmup_producer(mt19937_64_producer);
+    #if defined(ARMA_USE_THREAD_UNIQUE_RNG_SEED)
+      
+      // each thread has unique starting seed
+      
+      #if defined(ARMA_USE_OPENMP)
+        
+        static thread_local std::mt19937_64 mt19937_64_producer( std::mt19937_64::default_seed + arma_rng::seed_type(omp_get_thread_num()) );
+        
+      #else
+        
+        static std::atomic<std::size_t> mt19937_64_producer_counter(0);
+        
+        static thread_local std::mt19937_64 mt19937_64_producer( std::mt19937_64::default_seed + mt19937_64_producer_counter++ );
+        
+      #endif
+      
+    #else
+      
+      // each thread has the same starting seed
+      
+      static thread_local std::mt19937_64 mt19937_64_producer( std::mt19937_64::default_seed );
+      
+    #endif
     
   #else
     
-    // use a plain RNG in case we don't have thread_local
+    // plain RNG in case we don't have thread_local
     
     static std::mt19937_64 mt19937_64_producer( std::mt19937_64::default_seed );
     
-    arma_rng::warmup_producer(mt19937_64_producer);
-    
   #endif
   
+  arma_rng::warmup_producer(mt19937_64_producer);
+  
   return mt19937_64_producer;
   }
 
@@ -225,9 +251,45 @@ arma_rng::set_seed(const arma_rng::seed_type val)
     }
   #elif defined(ARMA_USE_CXX11_RNG)
     {
-    arma_rng::lock_producer();
-    arma_rng::get_producer().seed(val);
-    arma_rng::unlock_producer();
+    #if defined(ARMA_USE_OPENMP) && defined(ARMA_USE_THREAD_LOCAL)
+      {
+      arma_rng::lock_producer();
+      
+      #if defined(ARMA_USE_THREAD_UNIQUE_RNG_SEED)
+        constexpr bool thread_unique_rng_seed = true;
+      #else
+        constexpr bool thread_unique_rng_seed = false;
+      #endif
+      
+      // if we're already in a parallel region, assume the user is setting the seed for each thread
+      
+      if( (thread_unique_rng_seed == false) || bool(omp_in_parallel()) )
+        {
+        arma_rng::get_producer().seed(val);
+        }
+      else
+        {
+        const int n_threads = int( (std::max)( int(1), int(omp_get_max_threads()) ) );
+        
+        #pragma omp parallel for ordered schedule(static) num_threads(n_threads)
+        for(int t=0; t < n_threads; ++t)
+          {
+          #pragma omp ordered
+            {
+            arma_rng::get_producer().seed(val + arma_rng::seed_type(omp_get_thread_num()));
+            }
+          }
+        }
+      
+      arma_rng::unlock_producer();
+      }
+    #else
+      {
+      arma_rng::lock_producer();
+      arma_rng::get_producer().seed(val);
+      arma_rng::unlock_producer();
+      }
+    #endif
     }
   #else
     {
@@ -1029,4 +1091,45 @@ struct arma_rng::randg
 
 
 
+//
+
+
+
+template<typename eT>
+struct arma_rng::rande
+  {
+  inline
+  static
+  void
+  fill(eT* mem, const uword N, const double lambda)
+    {
+    #if defined(ARMA_USE_CXX11_RNG)
+      {
+      std::exponential_distribution<double> local_e_distr(lambda);
+      
+      std::mt19937_64& producer = arma_rng::get_producer();
+      
+      arma_rng::lock_producer();
+      
+      for(uword i=0; i<N; ++i)  { mem[i] = eT(local_e_distr(producer)); }
+      
+      arma_rng::unlock_producer();
+      }
+    #else
+      {
+      typedef typename std::mt19937_64::result_type local_seed_type;
+      
+      std::mt19937_64                       local_engine;
+      std::exponential_distribution<double> local_e_distr(lambda);
+      
+      local_engine.seed( local_seed_type(arma_rng::randi<local_seed_type>()) );
+      
+      for(uword i=0; i<N; ++i)  { mem[i] = eT(local_e_distr(local_engine)); }
+      }
+    #endif
+    }
+  };
+
+
+
 //! @}

include/armadillo_bits/arma_version.hpp

@@ -22,9 +22,9 @@
 
 
 #define ARMA_VERSION_MAJOR 15
-#define ARMA_VERSION_MINOR 0
-#define ARMA_VERSION_PATCH 1
-#define ARMA_VERSION_NAME  "Medium Roast"
+#define ARMA_VERSION_MINOR 2
+#define ARMA_VERSION_PATCH 4
+#define ARMA_VERSION_NAME  "Medium Roast Deluxe"
 
 
 

include/armadillo_bits/config.hpp

@@ -330,18 +330,6 @@
   #undef ARMA_OPTIMISE_POWEXPR
 #endif
 
-#if defined(ARMA_DONT_CHECK_CONFORMANCE)
-  #if defined(ARMA_CHECK_CONFORMANCE) && (ARMA_WARN_LEVEL >= 2)
-    #pragma message ("WARNING: conformance checks disabled")
-  #endif
-  
-  #undef ARMA_CHECK_CONFORMANCE
-#endif
-
-#if defined(ARMA_DONT_CHECK_NONFINITE)
-  #undef ARMA_CHECK_NONFINITE
-#endif
-
 #if defined(ARMA_NO_DEBUG)
   #undef ARMA_DEBUG
   #undef ARMA_EXTRA_DEBUG
@@ -361,6 +349,18 @@
   #define ARMA_WARN_LEVEL 3
 #endif
 
+#if defined(ARMA_DONT_CHECK_CONFORMANCE)
+  #if defined(ARMA_CHECK_CONFORMANCE) && (ARMA_WARN_LEVEL >= 2)
+    #pragma message ("WARNING: conformance checks disabled")
+  #endif
+  
+  #undef ARMA_CHECK_CONFORMANCE
+#endif
+
+#if defined(ARMA_DONT_CHECK_NONFINITE)
+  #undef ARMA_CHECK_NONFINITE
+#endif
+
 #if defined(ARMA_DONT_PRINT_EXCEPTIONS)
   #undef ARMA_PRINT_EXCEPTIONS
 #endif

include/armadillo_bits/config.hpp.cmake

@@ -330,18 +330,6 @@
   #undef ARMA_OPTIMISE_POWEXPR
 #endif
 
-#if defined(ARMA_DONT_CHECK_CONFORMANCE)
-  #if defined(ARMA_CHECK_CONFORMANCE) && (ARMA_WARN_LEVEL >= 2)
-    #pragma message ("WARNING: conformance checks disabled")
-  #endif
-  
-  #undef ARMA_CHECK_CONFORMANCE
-#endif
-
-#if defined(ARMA_DONT_CHECK_NONFINITE)
-  #undef ARMA_CHECK_NONFINITE
-#endif
-
 #if defined(ARMA_NO_DEBUG)
   #undef ARMA_DEBUG
   #undef ARMA_EXTRA_DEBUG
@@ -361,6 +349,18 @@
   #define ARMA_WARN_LEVEL 3
 #endif
 
+#if defined(ARMA_DONT_CHECK_CONFORMANCE)
+  #if defined(ARMA_CHECK_CONFORMANCE) && (ARMA_WARN_LEVEL >= 2)
+    #pragma message ("WARNING: conformance checks disabled")
+  #endif
+  
+  #undef ARMA_CHECK_CONFORMANCE
+#endif
+
+#if defined(ARMA_DONT_CHECK_NONFINITE)
+  #undef ARMA_CHECK_NONFINITE
+#endif
+
 #if defined(ARMA_DONT_PRINT_EXCEPTIONS)
   #undef ARMA_PRINT_EXCEPTIONS
 #endif

include/armadillo_bits/constants.hpp

@@ -26,6 +26,7 @@ namespace priv
     {
     template<typename eT>
     static
+    constexpr
     typename arma_real_only<eT>::result
     nan(typename arma_real_only<eT>::result* junk = nullptr)
       {
@@ -37,6 +38,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_cx_only<eT>::result
     nan(typename arma_cx_only<eT>::result* junk = nullptr)
       {
@@ -50,6 +52,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_integral_only<eT>::result
     nan(typename arma_integral_only<eT>::result* junk = nullptr)
       {
@@ -59,10 +62,14 @@ namespace priv
       }
     
     
+    //
+    
+    
     template<typename eT>
     static
+    constexpr
     typename arma_real_only<eT>::result
-    inf(typename arma_real_only<eT>::result* junk = nullptr)
+    pos_inf(typename arma_real_only<eT>::result* junk = nullptr)
       {
       arma_ignore(junk);
       
@@ -72,26 +79,70 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_cx_only<eT>::result
-    inf(typename arma_cx_only<eT>::result* junk = nullptr)
+    pos_inf(typename arma_cx_only<eT>::result* junk = nullptr)
       {
       arma_ignore(junk);
       
       typedef typename get_pod_type<eT>::result T;
       
-      return eT( Datum_helper::inf<T>(), Datum_helper::inf<T>() );
+      return eT( Datum_helper::pos_inf<T>(), Datum_helper::pos_inf<T>() );
       }
     
     
     template<typename eT>
     static
+    constexpr
     typename arma_integral_only<eT>::result
-    inf(typename arma_integral_only<eT>::result* junk = nullptr)
+    pos_inf(typename arma_integral_only<eT>::result* junk = nullptr)
       {
       arma_ignore(junk);
       
       return std::numeric_limits<eT>::max();
       }
+    
+    
+    //
+    
+    
+    template<typename eT>
+    static
+    constexpr
+    typename arma_real_only<eT>::result
+    neg_inf(typename arma_real_only<eT>::result* junk = nullptr)
+      {
+      arma_ignore(junk);
+      
+      return (std::numeric_limits<eT>::has_infinity) ? eT(-std::numeric_limits<eT>::infinity()) : eT(std::numeric_limits<eT>::lowest());
+      }
+    
+    
+    template<typename eT>
+    static
+    constexpr
+    typename arma_cx_only<eT>::result
+    neg_inf(typename arma_cx_only<eT>::result* junk = nullptr)
+      {
+      arma_ignore(junk);
+      
+      typedef typename get_pod_type<eT>::result T;
+      
+      return eT( Datum_helper::neg_inf<T>(), Datum_helper::neg_inf<T>() );
+      }
+    
+    
+    template<typename eT>
+    static
+    constexpr
+    typename arma_integral_only<eT>::result
+    neg_inf(typename arma_integral_only<eT>::result* junk = nullptr)
+      {
+      arma_ignore(junk);
+      
+      return std::numeric_limits<eT>::lowest();
+      }
+    
     };
   }
 
@@ -119,7 +170,9 @@ struct Datum
   static const eT log_min;      //!< log of the minimum representable value
   static const eT log_max;      //!< log of the maximum representable value
   static const eT nan;          //!< "not a number"
-  static const eT inf;          //!< infinity 
+  static const eT inf;          //!< positive infinity
+  static const eT pos_inf;      //!< positive infinity
+  static const eT neg_inf;      //!< negative infinity
 
   // 
   
@@ -170,7 +223,9 @@ template<typename eT> const eT Datum<eT>::eps         = std::numeric_limits<eT>:
 template<typename eT> const eT Datum<eT>::log_min     = std::log(std::numeric_limits<eT>::min());
 template<typename eT> const eT Datum<eT>::log_max     = std::log(std::numeric_limits<eT>::max());
 template<typename eT> const eT Datum<eT>::nan         = priv::Datum_helper::nan<eT>();
-template<typename eT> const eT Datum<eT>::inf         = priv::Datum_helper::inf<eT>();
+template<typename eT> const eT Datum<eT>::inf         = priv::Datum_helper::pos_inf<eT>();
+template<typename eT> const eT Datum<eT>::pos_inf     = priv::Datum_helper::pos_inf<eT>();
+template<typename eT> const eT Datum<eT>::neg_inf     = priv::Datum_helper::neg_inf<eT>();
 
 template<typename eT> const eT Datum<eT>::m_u       = eT(1.66053906892e-27);
 template<typename eT> const eT Datum<eT>::N_A       = eT(6.02214076e23);

include/armadillo_bits/diskio_bones.hpp

@@ -63,6 +63,9 @@ struct diskio
   
   template<typename eT> inline static std::streamsize prepare_stream(std::ostream& f);
   
+  template<typename eT> inline static constexpr eT real_as_int_lower_limit();
+  template<typename eT> inline static constexpr eT real_as_int_upper_limit();
+  
   
   //
   // matrix saving

include/armadillo_bits/diskio_meat.hpp

@@ -688,7 +688,32 @@ diskio::prepare_stream(std::ostream& f)
   
   return cell_width;
   }
+
+
+
+template<typename eT>
+inline
+constexpr
+eT
+diskio::real_as_int_lower_limit()
+  {
+  constexpr eT eT_int_accuracy_lower_limit = -( (is_fp16<eT>::value) ? eT(0x800) : ( (is_float<eT>::value) ? eT(0x1000000) : eT(0x20000000000000) ) );
   
+  return (std::max)( eT(std::numeric_limits<int>::lowest()), eT_int_accuracy_lower_limit );
+  }
+
+
+
+template<typename eT>
+inline
+constexpr
+eT
+diskio::real_as_int_upper_limit()
+  {
+  constexpr eT eT_int_accuracy_upper_limit = (is_fp16<eT>::value) ? eT(0x800) : ( (is_float<eT>::value) ? eT(0x1000000) : eT(0x20000000000000) );
+  
+  return (std::min)( eT(std::numeric_limits<int>::max()), eT_int_accuracy_upper_limit );
+  }
 
 
 
@@ -935,8 +960,8 @@ diskio::save_csv_ascii(const Mat<eT>& x, std::ostream& f, const char separator)
   uword x_n_rows = x.n_rows;
   uword x_n_cols = x.n_cols;
   
-  const eT eT_int_lowest = eT(std::numeric_limits<int>::lowest());
-  const eT eT_int_max    = eT(std::numeric_limits<int>::max());
+  constexpr eT eT_int_lower = diskio::real_as_int_lower_limit<eT>();
+  constexpr eT eT_int_upper = diskio::real_as_int_upper_limit<eT>();
   
   for(uword row=0; row < x_n_rows; ++row)
     {
@@ -944,7 +969,7 @@ diskio::save_csv_ascii(const Mat<eT>& x, std::ostream& f, const char separator)
       {
       const eT val = x.at(row,col);
       
-      const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val);
+      const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lower) && (val < eT_int_upper) && (eT(int(val)) == val);
       
       (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val);
       
@@ -977,8 +1002,8 @@ diskio::save_csv_ascii(const Mat< std::complex<T> >& x, std::ostream& f, const c
   
   diskio::prepare_stream<eT>(f);
   
-  const T T_int_lowest = T(std::numeric_limits<int>::lowest());
-  const T T_int_max    = T(std::numeric_limits<int>::max());
+  constexpr T T_int_lower = diskio::real_as_int_lower_limit<T>();
+  constexpr T T_int_upper = diskio::real_as_int_upper_limit<T>();
   
   uword x_n_rows = x.n_rows;
   uword x_n_cols = x.n_cols;
@@ -994,8 +1019,8 @@ diskio::save_csv_ascii(const Mat< std::complex<T> >& x, std::ostream& f, const c
       const T    abs_i = (val_i < T(0)) ? T(-val_i) : T(val_i);
       const char sgn_i = (val_i < T(0)) ? char('-') : char('+');
       
-      const bool val_r_is_real_int = (is_real<T>::yes) && arma_isfinite(val_r) && (val_r > T_int_lowest) && (val_r < T_int_max) && (T(int(val_r)) == val_r);
-      const bool abs_i_is_real_int = (is_real<T>::yes) && arma_isfinite(abs_i)                           && (abs_i < T_int_max) && (T(int(abs_i)) == abs_i);
+      const bool val_r_is_real_int = (is_real<T>::yes) && arma_isfinite(val_r) && (val_r > T_int_lower) && (val_r < T_int_upper) && (T(int(val_r)) == val_r);
+      const bool abs_i_is_real_int = (is_real<T>::yes) && arma_isfinite(abs_i)                          && (abs_i < T_int_upper) && (T(int(abs_i)) == abs_i);
       
       (val_r_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_r)) : arma_ostream::raw_print_elem(f, val_r);
       
@@ -1061,9 +1086,9 @@ diskio::save_coord_ascii(const Mat<eT>& x, std::ostream& f)
   
   diskio::prepare_stream<eT>(f);
   
-  const eT eT_zero       = eT(0);
-  const eT eT_int_lowest = eT(std::numeric_limits<int>::lowest());
-  const eT eT_int_max    = eT(std::numeric_limits<int>::max());
+  constexpr eT eT_zero      = eT(0);
+  constexpr eT eT_int_lower = diskio::real_as_int_lower_limit<eT>();
+  constexpr eT eT_int_upper = diskio::real_as_int_upper_limit<eT>();
   
   for(uword col=0; col < x.n_cols; ++col)
   for(uword row=0; row < x.n_rows; ++row)
@@ -1075,7 +1100,7 @@ diskio::save_coord_ascii(const Mat<eT>& x, std::ostream& f)
     f << row;  f.put(' ');
     f << col;  f.put(' ');
     
-    const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val);
+    const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lower) && (val < eT_int_upper) && (eT(int(val)) == val);
     
     (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val);
     
@@ -1116,9 +1141,9 @@ diskio::save_coord_ascii(const Mat< std::complex<T> >& x, std::ostream& f)
   
   diskio::prepare_stream<eT>(f);
   
-  const eT eT_zero       = eT(0);
-  const  T  T_int_lowest = T(std::numeric_limits<int>::lowest());
-  const  T  T_int_max    = T(std::numeric_limits<int>::max());
+  constexpr eT eT_zero      = eT(0);
+  constexpr  T  T_int_lower = diskio::real_as_int_lower_limit<T>();
+  constexpr  T  T_int_upper = diskio::real_as_int_upper_limit<T>();
   
   for(uword col=0; col < x.n_cols; ++col)
   for(uword row=0; row < x.n_rows; ++row)
@@ -1133,8 +1158,8 @@ diskio::save_coord_ascii(const Mat< std::complex<T> >& x, std::ostream& f)
     const T val_r = std::real(val);
     const T val_i = std::imag(val);
     
-    const bool val_r_is_real_int = (is_real<T>::yes) && arma_isfinite(val_r) && (val_r > T_int_lowest) && (val_r < T_int_max) && (T(int(val_r)) == val_r);
-    const bool val_i_is_real_int = (is_real<T>::yes) && arma_isfinite(val_i) && (val_i > T_int_lowest) && (val_i < T_int_max) && (T(int(val_i)) == val_i);
+    const bool val_r_is_real_int = (is_real<T>::yes) && arma_isfinite(val_r) && (val_r > T_int_lower) && (val_r < T_int_upper) && (T(int(val_r)) == val_r);
+    const bool val_i_is_real_int = (is_real<T>::yes) && arma_isfinite(val_i) && (val_i > T_int_lower) && (val_i < T_int_upper) && (T(int(val_i)) == val_i);
     
     (val_r_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_r)) : arma_ostream::raw_print_elem(f, val_r);
     
@@ -2966,9 +2991,9 @@ diskio::save_csv_ascii(const SpMat<eT>& x, std::ostream& f, const char separator
   uword x_n_rows = x.n_rows;
   uword x_n_cols = x.n_cols;
   
-  const eT eT_zero       = eT(0);
-  const eT eT_int_lowest = eT(std::numeric_limits<int>::lowest());
-  const eT eT_int_max    = eT(std::numeric_limits<int>::max());
+  constexpr eT eT_zero      = eT(0);
+  constexpr eT eT_int_lower = diskio::real_as_int_lower_limit<eT>();
+  constexpr eT eT_int_upper = diskio::real_as_int_upper_limit<eT>();
   
   for(uword row=0; row < x_n_rows; ++row)
     {
@@ -2982,7 +3007,7 @@ diskio::save_csv_ascii(const SpMat<eT>& x, std::ostream& f, const char separator
         }
       else
         {
-        const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val);
+        const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lower) && (val < eT_int_upper) && (eT(int(val)) == val);
         
         (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val);
         }
@@ -3064,8 +3089,8 @@ diskio::save_coord_ascii(const SpMat<eT>& x, std::ostream& f)
   
   diskio::prepare_stream<eT>(f);
   
-  const eT eT_int_lowest = eT(std::numeric_limits<int>::lowest());
-  const eT eT_int_max    = eT(std::numeric_limits<int>::max());
+  constexpr eT eT_int_lower = diskio::real_as_int_lower_limit<eT>();
+  constexpr eT eT_int_upper = diskio::real_as_int_upper_limit<eT>();
   
   typename SpMat<eT>::const_iterator iter     = x.begin();
   typename SpMat<eT>::const_iterator iter_end = x.end();
@@ -3077,7 +3102,7 @@ diskio::save_coord_ascii(const SpMat<eT>& x, std::ostream& f)
     
     const eT val = (*iter);
     
-    const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val);
+    const bool is_real_int = (is_real<eT>::yes) && arma_isfinite(val) && (val > eT_int_lower) && (val < eT_int_upper) && (eT(int(val)) == val);
     
     (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val);
     
@@ -3120,8 +3145,8 @@ diskio::save_coord_ascii(const SpMat< std::complex<T> >& x, std::ostream& f)
   
   diskio::prepare_stream<eT>(f);
   
-  const T T_int_lowest = T(std::numeric_limits<int>::lowest());
-  const T T_int_max    = T(std::numeric_limits<int>::max());
+  constexpr T T_int_lower = diskio::real_as_int_lower_limit<T>();
+  constexpr T T_int_upper = diskio::real_as_int_upper_limit<T>();
   
   typename SpMat<eT>::const_iterator iter     = x.begin();
   typename SpMat<eT>::const_iterator iter_end = x.end();
@@ -3136,8 +3161,8 @@ diskio::save_coord_ascii(const SpMat< std::complex<T> >& x, std::ostream& f)
     const T val_r = std::real(val);
     const T val_i = std::imag(val);
     
-    const bool val_r_is_real_int = (is_real<T>::yes) && arma_isfinite(val_r) && (val_r > T_int_lowest) && (val_r < T_int_max) && (T(int(val_r)) == val_r);
-    const bool val_i_is_real_int = (is_real<T>::yes) && arma_isfinite(val_i) && (val_i > T_int_lowest) && (val_i < T_int_max) && (T(int(val_i)) == val_i);
+    const bool val_r_is_real_int = (is_real<T>::yes) && arma_isfinite(val_r) && (val_r > T_int_lower) && (val_r < T_int_upper) && (T(int(val_r)) == val_r);
+    const bool val_i_is_real_int = (is_real<T>::yes) && arma_isfinite(val_i) && (val_i > T_int_lower) && (val_i < T_int_upper) && (T(int(val_i)) == val_i);
     
     (val_r_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_r)) : arma_ostream::raw_print_elem(f, val_r);
     

include/armadillo_bits/distr_param.hpp

@@ -48,6 +48,26 @@ class distr_param
     }
   
   
+  inline explicit distr_param(const int a)
+    : state   (1)
+    , a_int   (a)
+    , b_int   (0)
+    , a_double(double(a))
+    , b_double(double(0))
+    {
+    }
+  
+  
+  inline explicit distr_param(const double a)
+    : state   (2)
+    , a_int   (int(a))
+    , b_int   (int(0))
+    , a_double(a)
+    , b_double(0)
+    {
+    }
+  
+  
   inline explicit distr_param(const int a, const int b)
     : state   (1)
     , a_int   (a)

include/armadillo_bits/field_meat.hpp

@@ -94,7 +94,7 @@ field<oT>::field(const subview_field<oT>& X)
   {
   arma_debug_sigprint_this(this);
   
-  this->operator=(X);
+  (*this).operator=(X);
   }
 
 

include/armadillo_bits/fill.hpp

@@ -22,22 +22,29 @@
 
 namespace fill
   {
-  struct fill_none  {};
-  struct fill_zeros {};
-  struct fill_ones  {};
-  struct fill_eye   {};
-  struct fill_randu {};
-  struct fill_randn {};
+  struct fill_none    {};
+  struct fill_zeros   {};
+  struct fill_ones    {};
+  struct fill_eye     {};
+  struct fill_randu   {};
+  struct fill_randn   {};
+  struct fill_nan     {};
+  struct fill_pos_inf {};
+  struct fill_neg_inf {};
   
   template<typename fill_type> 
   struct fill_class { inline constexpr fill_class() {} };
   
-  static constexpr fill_class<fill_none > none;
-  static constexpr fill_class<fill_zeros> zeros;
-  static constexpr fill_class<fill_ones > ones;
-  static constexpr fill_class<fill_eye  > eye;
-  static constexpr fill_class<fill_randu> randu;
-  static constexpr fill_class<fill_randn> randn;
+  static constexpr fill_class<fill_none   > none;
+  static constexpr fill_class<fill_zeros  > zeros;
+  static constexpr fill_class<fill_ones   > ones;
+  static constexpr fill_class<fill_eye    > eye;
+  static constexpr fill_class<fill_randu  > randu;
+  static constexpr fill_class<fill_randn  > randn;
+  static constexpr fill_class<fill_nan    > nan;
+  static constexpr fill_class<fill_pos_inf> inf;
+  static constexpr fill_class<fill_pos_inf> pos_inf;
+  static constexpr fill_class<fill_neg_inf> neg_inf;
   
   //
   

include/armadillo_bits/fn_all.hpp

@@ -41,7 +41,7 @@ all(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -59,7 +59,7 @@ all(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -78,7 +78,7 @@ all(const T1& X, const uword dim)
 //! for compatibility purposes: allows compiling user code designed for earlier versions of Armadillo
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_any.hpp

@@ -41,7 +41,7 @@ any(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -59,7 +59,7 @@ any(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -78,7 +78,7 @@ any(const T1& X, const uword dim)
 //! for compatibility purposes: allows compiling user code designed for earlier versions of Armadillo
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_circshift.hpp

@@ -17,15 +17,14 @@
 
 
 
-//! \addtogroup fn_shift
+//! \addtogroup fn_circshift
 //! @{
 
 
 
-// TODO: deprecate in favour of circshift()
 template<typename T1>
-arma_warn_unused
-arma_inline
+// [[deprecated("use circshift() instead")]]
+inline
 typename
 enable_if2
   <
@@ -50,7 +49,7 @@ shift
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -73,9 +72,8 @@ circshift
 
 
 
-// TODO: deprecate in favour of circshift()
 template<typename T1>
-arma_warn_unused
+// [[deprecated("use circshift() instead")]]
 inline
 typename
 enable_if2
@@ -140,9 +138,8 @@ circshift
 
 
 
-// TODO: deprecate in favour of circshift()
 template<typename T1>
-arma_warn_unused
+// [[deprecated("use circshift() instead")]]
 inline
 typename
 enable_if2
@@ -217,9 +214,8 @@ circshift
 
 
 
-// TODO: deprecate in favour of circshift()
 template<typename T1>
-arma_warn_unused
+// [[deprecated("use circshift() instead")]]
 inline
 SpMat<typename T1::elem_type>
 shift

include/armadillo_bits/fn_cumprod.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -41,7 +41,7 @@ cumprod(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -59,7 +59,7 @@ cumprod(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -77,7 +77,7 @@ cumprod(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 cumprod(const T& x)
   {

include/armadillo_bits/fn_cumsum.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -41,7 +41,7 @@ cumsum(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -59,7 +59,7 @@ cumsum(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -77,7 +77,7 @@ cumsum(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 cumsum(const T& x)
   {

include/armadillo_bits/fn_det.hpp

@@ -70,7 +70,7 @@ det(typename T1::elem_type& out_val, const Base<typename T1::elem_type,T1>& X)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 det(const T& x)
   {

include/armadillo_bits/fn_diagmat.hpp

@@ -23,7 +23,7 @@
 //! interpret a matrix or a vector as a diagonal matrix (ie. off-diagonal entries are zero)
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -42,7 +42,7 @@ diagmat(const T1& X)
 //! create a matrix with the k-th diagonal set to the given vector
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_diagvec.hpp

@@ -23,7 +23,7 @@
 //! extract main diagonal from matrix
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const Op<T1, op_diagvec>
 diagvec(const Base<typename T1::elem_type,T1>& X)
   {
@@ -37,7 +37,7 @@ diagvec(const Base<typename T1::elem_type,T1>& X)
 //! extract arbitrary diagonal from matrix
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const Op<T1, op_diagvec2>
 diagvec(const Base<typename T1::elem_type,T1>& X, const sword diag_id)
   {
@@ -50,7 +50,7 @@ diagvec(const Base<typename T1::elem_type,T1>& X, const sword diag_id)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const mtSpReduceOp<typename T1::elem_type, T1, op_sp_diagvec>
 diagvec(const SpBase<typename T1::elem_type,T1>& X, const sword diag_id = 0)
   {

include/armadillo_bits/fn_diff.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -45,7 +45,7 @@ diff
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -67,7 +67,7 @@ diff
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_expmat.hpp

@@ -69,7 +69,7 @@ expmat(Mat<typename T1::elem_type>& B, const Base<typename T1::elem_type,T1>& A)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_blas_type<typename T1::elem_type>::value, const Op<T1, op_expmat_sym> >::result
 expmat_sym(const Base<typename T1::elem_type,T1>& X)
   {

include/armadillo_bits/fn_flip.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_flipud> >::result
 flipud(const T1& X)
   {
@@ -36,7 +36,7 @@ flipud(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_fliplr> >::result
 fliplr(const T1& X)
   {
@@ -49,7 +49,7 @@ fliplr(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const SpOp<T1, spop_flipud>
 flipud(const SpBase<typename T1::elem_type,T1>& X)
   {
@@ -62,7 +62,7 @@ flipud(const SpBase<typename T1::elem_type,T1>& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const SpOp<T1, spop_fliplr>
 fliplr(const SpBase<typename T1::elem_type,T1>& X)
   {

include/armadillo_bits/fn_hist.hpp

@@ -22,7 +22,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -40,7 +40,7 @@ hist(const T1& A, const uword n_bins = 10)
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -58,7 +58,7 @@ hist(const T1& X, const T2& Y)
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_histc.hpp

@@ -22,7 +22,7 @@
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -40,7 +40,7 @@ histc(const T1& X, const T2& Y)
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_index_max.hpp

@@ -35,7 +35,7 @@ index_max(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const mtOp<uword, T1, op_index_max> >::result
 index_max(const T1& X)
   {
@@ -48,7 +48,7 @@ index_max(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const mtOp<uword, T1, op_index_max> >::result
 index_max(const T1& X, const uword dim)
   {
@@ -61,7 +61,7 @@ index_max(const T1& X, const uword dim)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const mtOpCube<uword, T1, op_index_max>
 index_max
   (

include/armadillo_bits/fn_index_min.hpp

@@ -35,7 +35,7 @@ index_min(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const mtOp<uword, T1, op_index_min> >::result
 index_min(const T1& X)
   {
@@ -48,7 +48,7 @@ index_min(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const mtOp<uword, T1, op_index_min> >::result
 index_min(const T1& X, const uword dim)
   {
@@ -61,7 +61,7 @@ index_min(const T1& X, const uword dim)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const mtOpCube<uword, T1, op_index_min>
 index_min
   (

include/armadillo_bits/fn_kron.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 const Glue<T1,T2,glue_kron>
 kron(const Base<typename T1::elem_type,T1>& A, const Base<typename T1::elem_type,T2>& B)
   {
@@ -90,7 +90,7 @@ kron(const Base<T,T1>& X, const Base<std::complex<T>,T2>& Y)
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 const SpGlue<T1, T2, spglue_kron>
 kron(const SpBase<typename T1::elem_type,T1>& A, const SpBase<typename T1::elem_type,T2>& B)
   {

include/armadillo_bits/fn_logmat.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< (is_blas_type<typename T1::elem_type>::value && is_cx<typename T1::elem_type>::no), const mtOp<std::complex<typename T1::elem_type>, T1, op_logmat> >::result
 logmat(const Base<typename T1::elem_type,T1>& X, const uword n_iters = 100u)
   {
@@ -36,7 +36,7 @@ logmat(const Base<typename T1::elem_type,T1>& X, const uword n_iters = 100u)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< (is_blas_type<typename T1::elem_type>::value && is_cx<typename T1::elem_type>::yes), const Op<T1, op_logmat_cx> >::result
 logmat(const Base<typename T1::elem_type,T1>& X, const uword n_iters = 100u)
   {
@@ -93,7 +93,7 @@ logmat(Mat<typename T1::elem_type>& Y, const Base<typename T1::elem_type,T1>& X,
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_blas_type<typename T1::elem_type>::value, const Op<T1, op_logmat_sympd> >::result
 logmat_sympd(const Base<typename T1::elem_type,T1>& X)
   {

include/armadillo_bits/fn_max.hpp

@@ -35,7 +35,7 @@ max(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const Op<T1, op_max> >::result
 max(const T1& X)
   {
@@ -48,7 +48,7 @@ max(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_max> >::result
 max(const T1& X, const uword dim)
   {
@@ -61,7 +61,7 @@ max(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 max(const T& x)
   {
@@ -73,7 +73,7 @@ max(const T& x)
 //! element-wise maximum
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -95,7 +95,7 @@ max
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const OpCube<T1, op_max>
 max
   (
@@ -112,7 +112,7 @@ max
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 const GlueCube<T1, T2, glue_max>
 max
   (

include/armadillo_bits/fn_mean.hpp

@@ -36,7 +36,7 @@ mean(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const Op<T1, op_mean> >::result
 mean(const T1& X)
   {
@@ -49,7 +49,7 @@ mean(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_mean> >::result
 mean(const T1& X, const uword dim)
   {
@@ -62,7 +62,7 @@ mean(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 mean(const T& x)
   {
@@ -73,7 +73,7 @@ mean(const T& x)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const OpCube<T1, op_mean>
 mean
   (

include/armadillo_bits/fn_median.hpp

@@ -35,7 +35,7 @@ median(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const Op<T1, op_median> >::result
 median(const T1& X)
   {
@@ -48,7 +48,7 @@ median(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_median> >::result
 median(const T1& X, const uword dim)
   {
@@ -61,7 +61,7 @@ median(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 median(const T& x)
   {

include/armadillo_bits/fn_min.hpp

@@ -35,7 +35,7 @@ min(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const Op<T1, op_min> >::result
 min(const T1& X)
   {
@@ -48,7 +48,7 @@ min(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_min> >::result
 min(const T1& X, const uword dim)
   {
@@ -61,7 +61,7 @@ min(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 min(const T& x)
   {
@@ -73,7 +73,7 @@ min(const T& x)
 //! element-wise minimum
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -95,7 +95,7 @@ min
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const OpCube<T1, op_min>
 min
   (
@@ -112,7 +112,7 @@ min
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 const GlueCube<T1, T2, glue_min>
 min
   (

include/armadillo_bits/fn_normalise.hpp

@@ -95,7 +95,7 @@ normalise
 //! for compatibility purposes: allows compiling user code designed for earlier versions of Armadillo
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_orth_null.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_blas_real<typename T1::pod_type>::value, const Op<T1, op_orth> >::result
 orth(const Base<typename T1::elem_type, T1>& X, const typename T1::pod_type tol = 0.0)
   {
@@ -62,7 +62,7 @@ orth(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type, T1>& X
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_blas_real<typename T1::pod_type>::value, const Op<T1, op_null> >::result
 null(const Base<typename T1::elem_type, T1>& X, const typename T1::pod_type tol = 0.0)
   {

include/armadillo_bits/fn_prod.hpp

@@ -20,13 +20,6 @@
 //! @{
 
 
-//! \brief
-//! Delayed product of elements of a matrix along a specified dimension (either rows or columns).
-//! The result is stored in a dense matrix that has either one column or one row.
-//! For dim = 0, find the sum of each column (ie. traverse across rows)
-//! For dim = 1, find the sum of each row (ie. traverse across columns)
-//! The default is dim = 0.
-//! NOTE: this function works differently than in Matlab/Octave.
 
 template<typename T1>
 arma_warn_unused
@@ -43,7 +36,7 @@ prod(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const Op<T1, op_prod> >::result
 prod(const T1& X)
   {
@@ -56,7 +49,7 @@ prod(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_prod> >::result
 prod(const T1& X, const uword dim)
   {
@@ -69,7 +62,7 @@ prod(const T1& X, const uword dim)
 
 template<typename T>
 arma_warn_unused
-arma_inline
+inline
 typename arma_scalar_only<T>::result
 prod(const T& x)
   {

include/armadillo_bits/fn_quantile.hpp

@@ -22,7 +22,7 @@
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -40,7 +40,7 @@ quantile(const T1& X, const Base<typename T2::elem_type,T2>& P)
 
 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_rande.hpp

@@ -0,0 +1,241 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (https://conradsanderson.id.au)
+// Copyright 2008-2016 National ICT Australia (NICTA)
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// https://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+// ------------------------------------------------------------------------
+
+
+//! \addtogroup fn_rande
+//! @{
+
+
+
+template<typename obj_type>
+arma_warn_unused
+inline
+obj_type
+rande(const uword n_rows, const uword n_cols, const distr_param& param = distr_param(), const typename arma_Mat_Col_Row_only<obj_type>::result* junk = nullptr)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  typedef typename obj_type::elem_type eT;
+  
+  if(is_Col<obj_type>::value)
+    {
+    arma_conform_check( (n_cols != 1), "rande(): incompatible size" );
+    }
+  else
+  if(is_Row<obj_type>::value)
+    {
+    arma_conform_check( (n_rows != 1), "rande(): incompatible size" );
+    }
+  
+  double lambda = double(1);
+  double unused = double(0);
+  
+  param.get_double_vals(lambda, unused);
+  
+  arma_conform_check( (lambda <= double(0)), "rande(): incorrect distribution parameters; lambda must be greater than zero" );
+  
+  obj_type out(n_rows, n_cols, arma_nozeros_indicator());
+  
+  arma_rng::rande<eT>::fill(out.memptr(), out.n_elem, lambda);
+  
+  return out;
+  }
+
+
+
+template<typename obj_type>
+arma_warn_unused
+inline
+obj_type
+rande(const SizeMat& s, const distr_param& param = distr_param(), const typename arma_Mat_Col_Row_only<obj_type>::result* junk = nullptr)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  return rande<obj_type>(s.n_rows, s.n_cols, param);
+  }
+
+
+
+template<typename obj_type>
+arma_warn_unused
+inline
+obj_type
+rande(const uword n_elem, const distr_param& param = distr_param(), const arma_empty_class junk1 = arma_empty_class(), const typename arma_Mat_Col_Row_only<obj_type>::result* junk2 = nullptr)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk1);
+  arma_ignore(junk2);
+  
+  const uword n_rows = (is_Row<obj_type>::value) ? uword(1) : n_elem;
+  const uword n_cols = (is_Row<obj_type>::value) ? n_elem   : uword(1);
+  
+  return rande<obj_type>(n_rows, n_cols, param);
+  }
+
+
+
+arma_warn_unused
+inline
+mat
+rande(const uword n_rows, const uword n_cols, const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  return rande<mat>(n_rows, n_cols, param);
+  }
+
+
+
+arma_warn_unused
+inline
+mat
+rande(const SizeMat& s, const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  return rande<mat>(s.n_rows, s.n_cols, param);
+  }
+
+
+
+arma_warn_unused
+inline
+vec
+rande(const uword n_elem, const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  return rande<vec>(n_elem, uword(1), param);
+  }
+
+
+
+arma_warn_unused
+inline
+double
+rande(const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  double lambda = double(1);
+  double unused = double(0);
+  
+  param.get_double_vals(lambda, unused);
+  
+  arma_conform_check( (lambda <= double(0)), "rande(): incorrect distribution parameters; lambda must be greater than zero" );
+  
+  double out_val = double(0);
+  
+  arma_rng::rande<double>::fill(&out_val, uword(1), lambda);
+  
+  return out_val;
+  }
+
+
+
+template<typename eT>
+arma_warn_unused
+inline
+typename arma_real_or_cx_only<eT>::result
+rande(const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  double lambda = double(1);
+  double unused = double(0);
+  
+  param.get_double_vals(lambda, unused);
+  
+  arma_conform_check( (lambda <= double(0)), "rande(): incorrect distribution parameters; lambda must be greater than zero" );
+  
+  eT out_val = eT(0);
+  
+  arma_rng::rande<eT>::fill(&out_val, uword(1), lambda);
+  
+  return out_val;
+  }
+
+
+
+template<typename cube_type>
+arma_warn_unused
+inline
+cube_type
+rande(const uword n_rows, const uword n_cols, const uword n_slices, const distr_param& param = distr_param(), const typename arma_Cube_only<cube_type>::result* junk = nullptr)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  typedef typename cube_type::elem_type eT;
+  
+  double lambda = double(1);
+  double unused = double(0);
+  
+  param.get_double_vals(lambda, unused);
+  
+  arma_conform_check( (lambda <= double(0)), "rande(): incorrect distribution parameters; lambda must be greater than zero" );
+  
+  cube_type out(n_rows, n_cols, n_slices, arma_nozeros_indicator());
+  
+  arma_rng::rande<eT>::fill(out.memptr(), out.n_elem, lambda);
+  
+  return out;
+  }
+
+
+
+template<typename cube_type>
+arma_warn_unused
+inline
+cube_type
+rande(const SizeCube& s, const distr_param& param = distr_param(), const typename arma_Cube_only<cube_type>::result* junk = nullptr)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  return rande<cube_type>(s.n_rows, s.n_cols, s.n_slices, param);
+  }
+
+
+
+arma_warn_unused
+inline
+cube
+rande(const uword n_rows, const uword n_cols, const uword n_slices, const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  return rande<cube>(n_rows, n_cols, n_slices, param);
+  }
+
+
+
+arma_warn_unused
+inline
+cube
+rande(const SizeCube& s, const distr_param& param = distr_param())
+  {
+  arma_debug_sigprint();
+  
+  return rande<cube>(s.n_rows, s.n_cols, s.n_slices, param);
+  }
+
+
+
+//! @}

include/armadillo_bits/fn_range.hpp

@@ -35,7 +35,7 @@ range(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value && resolves_to_vector<T1>::no, const Op<T1, op_range> >::result
 range(const T1& X)
   {
@@ -48,7 +48,7 @@ range(const T1& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_range> >::result
 range(const T1& X, const uword dim)
   {

include/armadillo_bits/fn_repelem.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -41,7 +41,7 @@ repelem(const T1& A, const uword r, const uword c)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const SpOp<T1, spop_repelem>
 repelem(const SpBase<typename T1::elem_type,T1>& A, const uword r, const uword c)
   {

include/armadillo_bits/fn_shuffle.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -44,7 +44,7 @@ shuffle
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -65,7 +65,7 @@ shuffle
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_sort.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <
@@ -44,7 +44,7 @@ sort
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
   <

include/armadillo_bits/fn_sort_index.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const mtOp<uword,T1,op_sort_index>
 sort_index
   (
@@ -70,7 +70,7 @@ sort_index
 // DO NOT USE: kept only for compatibility with old user code
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 const mtOp<uword,T1,op_sort_index>
 stable_sort_index
   (

include/armadillo_bits/fn_sqrtmat.hpp

@@ -23,7 +23,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< (is_blas_type<typename T1::elem_type>::value && is_cx<typename T1::elem_type>::no), const mtOp<std::complex<typename T1::elem_type>, T1, op_sqrtmat> >::result
 sqrtmat(const Base<typename T1::elem_type,T1>& X)
   {
@@ -36,7 +36,7 @@ sqrtmat(const Base<typename T1::elem_type,T1>& X)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< (is_blas_type<typename T1::elem_type>::value && is_cx<typename T1::elem_type>::yes), const Op<T1, op_sqrtmat_cx> >::result
 sqrtmat(const Base<typename T1::elem_type,T1>& X)
   {
@@ -91,7 +91,7 @@ sqrtmat(Mat<typename T1::elem_type>& Y, const Base<typename T1::elem_type,T1>& X
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_blas_type<typename T1::elem_type>::value, const Op<T1, op_sqrtmat_sympd> >::result
 sqrtmat_sympd(const Base<typename T1::elem_type,T1>& X)
   {

include/armadillo_bits/fn_symmat.hpp

@@ -22,7 +22,7 @@
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::no, const Op<T1, op_symmatu> >::result
 symmatu(const Base<typename T1::elem_type,T1>& X, const bool do_conj = false)
   {
@@ -36,7 +36,7 @@ symmatu(const Base<typename T1::elem_type,T1>& X, const bool do_conj = false)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::no, const Op<T1, op_symmatl> >::result
 symmatl(const Base<typename T1::elem_type,T1>& X, const bool do_conj = false)
   {
@@ -50,7 +50,7 @@ symmatl(const Base<typename T1::elem_type,T1>& X, const bool do_conj = false)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::yes, const Op<T1, op_symmatu_cx> >::result
 symmatu(const Base<typename T1::elem_type,T1>& X, const bool do_conj = true)
   {
@@ -63,7 +63,7 @@ symmatu(const Base<typename T1::elem_type,T1>& X, const bool do_conj = true)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::yes, const Op<T1, op_symmatl_cx> >::result
 symmatl(const Base<typename T1::elem_type,T1>& X, const bool do_conj = true)
   {
@@ -80,7 +80,7 @@ symmatl(const Base<typename T1::elem_type,T1>& X, const bool do_conj = true)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::no, const SpOp<T1, spop_symmat> >::result
 symmatu(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = false)
   {
@@ -94,7 +94,7 @@ symmatu(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = false)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::no, const SpOp<T1, spop_symmat> >::result
 symmatl(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = false)
   {
@@ -108,7 +108,7 @@ symmatl(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = false)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::yes, const SpOp<T1, spop_symmat_cx> >::result
 symmatu(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = true)
   {
@@ -121,7 +121,7 @@ symmatu(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = true)
 
 template<typename T1>
 arma_warn_unused
-arma_inline
+inline
 typename enable_if2< is_cx<typename T1::elem_type>::yes, const SpOp<T1, spop_symmat_cx> >::result
 symmatl(const SpBase<typename T1::elem_type,T1>& X, const bool do_conj = true)
   {

include/armadillo_bits/glue_hist_meat.hpp

@@ -109,7 +109,7 @@ glue_hist::apply_noalias(Mat<uword>& out, const Mat<eT>& X, const Mat<eT>& C, co
         {
         const eT val = X_mem[i];
         
-        if(is_finite(val))
+        if(arma_isfinite(val))
           {
           eT    opt_dist  = (val >= center_0) ? (val - center_0) : (center_0 - val);
           uword opt_index = 0;

include/armadillo_bits/glue_join_bones.hpp

@@ -36,7 +36,10 @@ struct glue_join_cols
   inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_join_cols>& X);
   
   template<typename T1, typename T2>
-  inline static void apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& A, const Proxy<T2>& B);
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_join_cols>& X);
+  
+  template<typename eT>
+  inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B);
   
   template<typename eT, typename T1, typename T2, typename T3>
   inline static void apply(Mat<eT>& out, const Base<eT,T1>& A, const Base<eT,T2>& B, const Base<eT,T3>& C);
@@ -61,7 +64,10 @@ struct glue_join_rows
   inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_join_rows>& X);
   
   template<typename T1, typename T2>
-  inline static void apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& A, const Proxy<T2>& B);
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_join_rows>& X);
+  
+  template<typename eT>
+  inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B);
   
   template<typename eT, typename T1, typename T2, typename T3>
   inline static void apply(Mat<eT>& out, const Base<eT,T1>& A, const Base<eT,T2>& B, const Base<eT,T3>& C);

include/armadillo_bits/glue_join_meat.hpp

@@ -21,44 +21,38 @@
 
 
 
-template<typename T1, typename T2>
+template<typename eT>
 inline
 void
-glue_join_cols::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& A, const Proxy<T2>& B)
+glue_join_cols::apply_noalias(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
   {
   arma_debug_sigprint();
   
-  const uword A_n_rows = A.get_n_rows();
-  const uword A_n_cols = A.get_n_cols();
-  
-  const uword B_n_rows = B.get_n_rows();
-  const uword B_n_cols = B.get_n_cols();
-  
   arma_conform_check
     (
-    ( (A_n_cols != B_n_cols) && ( (A_n_rows > 0) || (A_n_cols > 0) ) && ( (B_n_rows > 0) || (B_n_cols > 0) ) ),
+    ( (A.n_cols != B.n_cols) && ( (A.n_rows > 0) || (A.n_cols > 0) ) && ( (B.n_rows > 0) || (B.n_cols > 0) ) ),
     "join_cols() / join_vert(): number of columns must be the same"
     );
   
-  out.set_size( A_n_rows + B_n_rows, (std::max)(A_n_cols, B_n_cols) );
+  out.set_size( A.n_rows + B.n_rows, (std::max)(A.n_cols, B.n_cols) );
   
   if( out.n_elem > 0 )
     {
-    if(A.get_n_elem() > 0)
+    if(A.n_elem> 0)
       { 
-      out.submat(0,        0,   A_n_rows-1, out.n_cols-1) = A.Q;
+      out.submat(0,        0,   A.n_rows-1, out.n_cols-1) = A;
       }
     
-    if(B.get_n_elem() > 0)
+    if(B.n_elem > 0)
       {
-      out.submat(A_n_rows, 0, out.n_rows-1, out.n_cols-1) = B.Q;
+      out.submat(A.n_rows, 0, out.n_rows-1, out.n_cols-1) = B;
       }
     }
   }
-  
-  
-  
-  
+
+
+
+
 template<typename T1, typename T2>
 inline
 void
@@ -68,18 +62,18 @@ glue_join_cols::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_jo
   
   typedef typename T1::elem_type eT;
   
-  const Proxy<T1> A(X.A);
-  const Proxy<T2> B(X.B);
+  const quasi_unwrap<T1> UA(X.A);
+  const quasi_unwrap<T2> UB(X.B);
   
-  if( (A.is_alias(out) == false) && (B.is_alias(out) == false) )
+  if( (UA.is_alias(out) == false) && (UB.is_alias(out) == false) )
     {
-    glue_join_cols::apply_noalias(out, A, B);
+    glue_join_cols::apply_noalias(out, UA.M, UB.M);
     }
   else
     {
     Mat<eT> tmp;
     
-    glue_join_cols::apply_noalias(tmp, A, B);
+    glue_join_cols::apply_noalias(tmp, UA.M, UB.M);
     
     out.steal_mem(tmp);
     }
@@ -87,6 +81,21 @@ glue_join_cols::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_jo
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_join_cols::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_join_cols>& X)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> UA(X.A);
+  const quasi_unwrap<T2> UB(X.B);
+  
+  glue_join_cols::apply_noalias(out, UA.M, UB.M);
+  }
+
+
+
 template<typename eT, typename T1, typename T2, typename T3>
 inline
 void
@@ -178,44 +187,42 @@ glue_join_cols::apply(Mat<eT>& out, const Base<eT,T1>& A_expr, const Base<eT,T2>
 
 
 
-template<typename T1, typename T2>
+//
+
+
+
+template<typename eT>
 inline
 void
-glue_join_rows::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& A, const Proxy<T2>& B)
+glue_join_rows::apply_noalias(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
   {
   arma_debug_sigprint();
   
-  const uword A_n_rows = A.get_n_rows();
-  const uword A_n_cols = A.get_n_cols();
-  
-  const uword B_n_rows = B.get_n_rows();
-  const uword B_n_cols = B.get_n_cols();
-  
   arma_conform_check
     (
-    ( (A_n_rows != B_n_rows) && ( (A_n_rows > 0) || (A_n_cols > 0) ) && ( (B_n_rows > 0) || (B_n_cols > 0) ) ),
+    ( (A.n_rows != B.n_rows) && ( (A.n_rows > 0) || (A.n_cols > 0) ) && ( (B.n_rows > 0) || (B.n_cols > 0) ) ),
     "join_rows() / join_horiz(): number of rows must be the same"
     );
   
-  out.set_size( (std::max)(A_n_rows, B_n_rows), A_n_cols + B_n_cols );
+  out.set_size( (std::max)(A.n_rows, B.n_rows), A.n_cols + B.n_cols );
   
   if( out.n_elem > 0 )
     {
-    if(A.get_n_elem() > 0)
+    if(A.n_elem > 0)
       {
-      out.submat(0, 0,        out.n_rows-1,   A_n_cols-1) = A.Q;
+      out.submat(0, 0,        out.n_rows-1,   A.n_cols-1) = A;
       }
     
-    if(B.get_n_elem() > 0)
+    if(B.n_elem > 0)
       {
-      out.submat(0, A_n_cols, out.n_rows-1, out.n_cols-1) = B.Q;
+      out.submat(0, A.n_cols, out.n_rows-1, out.n_cols-1) = B;
       }
     }
   }
-  
-  
-  
-  
+
+
+
+
 template<typename T1, typename T2>
 inline
 void
@@ -225,18 +232,18 @@ glue_join_rows::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_jo
   
   typedef typename T1::elem_type eT;
   
-  const Proxy<T1> A(X.A);
-  const Proxy<T2> B(X.B);
+  const quasi_unwrap<T1> UA(X.A);
+  const quasi_unwrap<T2> UB(X.B);
   
-  if( (A.is_alias(out) == false) && (B.is_alias(out) == false) )
+  if( (UA.is_alias(out) == false) && (UB.is_alias(out) == false) )
     {
-    glue_join_rows::apply_noalias(out, A, B);
+    glue_join_rows::apply_noalias(out, UA.M, UB.M);
     }
   else
     {
     Mat<eT> tmp;
     
-    glue_join_rows::apply_noalias(tmp, A, B);
+    glue_join_rows::apply_noalias(tmp, UA.M, UB.M);
     
     out.steal_mem(tmp);
     }
@@ -244,6 +251,21 @@ glue_join_rows::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_jo
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_join_rows::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_join_rows>& X)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> UA(X.A);
+  const quasi_unwrap<T2> UB(X.B);
+  
+  glue_join_rows::apply_noalias(out, UA.M, UB.M);
+  }
+
+
+
 template<typename eT, typename T1, typename T2, typename T3>
 inline
 void
@@ -335,6 +357,10 @@ glue_join_rows::apply(Mat<eT>& out, const Base<eT,T1>& A_expr, const Base<eT,T2>
 
 
 
+//
+
+
+
 template<typename T1, typename T2>
 inline
 void

include/armadillo_bits/glue_kron_bones.hpp

@@ -36,7 +36,9 @@ struct glue_kron
   template<typename T>  inline static void direct_kron(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const Mat<T>&                 B);
   template<typename T>  inline static void direct_kron(Mat< std::complex<T> >& out, const Mat<T>&                 A, const Mat< std::complex<T> >& B);
   
-  template<typename T1, typename T2>   inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_kron>& X);
+  template<typename T1, typename T2> inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_kron>& X);
+  
+  template<typename T1, typename T2> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_kron>& X);
   };
 
 

include/armadillo_bits/glue_kron_meat.hpp

@@ -21,8 +21,6 @@
 
 
 
-//! \brief
-//! both input matrices have the same element type
 template<typename eT>
 inline
 void
@@ -50,9 +48,6 @@ glue_kron::direct_kron(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
 
 
 
-//! \brief
-//! different types of input matrices
-//! A -> complex, B -> basic element type
 template<typename T>
 inline
 void
@@ -84,9 +79,6 @@ glue_kron::direct_kron(Mat< std::complex<T> >& out, const Mat< std::complex<T> >
 
 
 
-//! \brief
-//! different types of input matrices
-//! A -> basic element type, B -> complex
 template<typename T>
 inline
 void
@@ -114,8 +106,6 @@ glue_kron::direct_kron(Mat< std::complex<T> >& out, const Mat<T>& A, const Mat<
 
 
 
-//! \brief
-//! apply Kronecker product for two objects with same element type
 template<typename T1, typename T2>
 inline
 void
@@ -144,4 +134,19 @@ glue_kron::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_kron>&
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_kron::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_kron>& X)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> UA(X.A);
+  const quasi_unwrap<T2> UB(X.B);
+  
+  glue_kron::direct_kron(out, UA.M, UB.M); 
+  }
+
+
+
 //! @}

include/armadillo_bits/glue_polyfit_bones.hpp

@@ -36,7 +36,11 @@ struct glue_polyfit
   
   template<typename T1, typename T2> inline static bool apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& X_expr, const Base<typename T1::elem_type, T2>& Y_expr, const uword N);
   
+  template<typename T1, typename T2> inline static bool apply_direct(Mat_noalias<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& X_expr, const Base<typename T1::elem_type, T2>& Y_expr, const uword N);
+  
   template<typename T1, typename T2> inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyfit>& expr);
+  
+  template<typename T1, typename T2> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyfit>& expr);
   };
 
 

include/armadillo_bits/glue_polyfit_meat.hpp

@@ -83,11 +83,7 @@ glue_polyfit::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
   
   arma_conform_check( (X.n_elem != Y.n_elem), "polyfit(): given vectors must have the same number of elements" );
   
-  if(X.n_elem == 0)
-    {
-    out.reset();
-    return true;
-    }
+  if(X.n_elem == 0)  { out.reset(); return true; }
   
   arma_conform_check( (N >= X.n_elem), "polyfit(): N must be less than the number of elements in X" );
   
@@ -99,7 +95,9 @@ glue_polyfit::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
   if(UX.is_alias(out) || UY.is_alias(out))
     {
     Mat<eT> tmp;
+    
     status = glue_polyfit::apply_noalias(tmp, X_as_colvec, Y_as_colvec, N);
+    
     out.steal_mem(tmp);
     }
   else
@@ -112,6 +110,41 @@ glue_polyfit::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
 
 
 
+template<typename T1, typename T2>
+inline
+bool
+glue_polyfit::apply_direct(Mat_noalias<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& X_expr, const Base<typename T1::elem_type, T2>& Y_expr, const uword N)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const quasi_unwrap<T1> UX(X_expr.get_ref());
+  const quasi_unwrap<T2> UY(Y_expr.get_ref());
+  
+  const Mat<eT>& X = UX.M;
+  const Mat<eT>& Y = UY.M;
+  
+  arma_conform_check
+    (
+    ( ((X.is_vec() == false) && (X.is_empty() == false)) || ((Y.is_vec() == false) && (Y.is_empty() == false)) ),
+    "polyfit(): given object must be a vector"
+    );
+  
+  arma_conform_check( (X.n_elem != Y.n_elem), "polyfit(): given vectors must have the same number of elements" );
+  
+  if(X.n_elem == 0)  { out.reset(); return true; }
+  
+  arma_conform_check( (N >= X.n_elem), "polyfit(): N must be less than the number of elements in X" );
+  
+  const Col<eT> X_as_colvec( const_cast<eT*>(X.memptr()), X.n_elem, false, false);
+  const Col<eT> Y_as_colvec( const_cast<eT*>(Y.memptr()), Y.n_elem, false, false);
+  
+  return glue_polyfit::apply_noalias(out, X_as_colvec, Y_as_colvec, N);
+  }
+
+
+
 template<typename T1, typename T2>
 inline
 void
@@ -130,4 +163,22 @@ glue_polyfit::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_poly
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_polyfit::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyfit>& expr)
+  {
+  arma_debug_sigprint();
+  
+  const bool status = glue_polyfit::apply_direct(out, expr.A, expr.B, expr.aux_uword);
+  
+  if(status == false)
+    {
+    out.soft_reset();
+    arma_stop_runtime_error("polyfit(): failed");
+    }
+  }
+
+
+
 //! @}

include/armadillo_bits/glue_polyval_bones.hpp

@@ -35,9 +35,10 @@ struct glue_polyval
   template<typename eT> inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& P, const Mat<eT>& X);
   
   template<typename T1, typename T2> inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyval>& expr);
+  
+  template<typename T1, typename T2> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyval>& expr);
   };
 
 
 
 //! @}
-

include/armadillo_bits/glue_polyval_meat.hpp

@@ -55,29 +55,43 @@ glue_polyval::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_poly
   const quasi_unwrap<T1> UP(expr.A);
   const quasi_unwrap<T2> UX(expr.B);
   
-  const Mat<eT>& P = UP.M;
-  const Mat<eT>& X = UX.M;
+  arma_conform_check( ((UP.M.is_vec() == false) && (UP.M.is_empty() == false)), "polyval(): argument P must be a vector" );
   
-  arma_conform_check( ((P.is_vec() == false) && (P.is_empty() == false)), "polyval(): argument P must be a vector" );
-  
-  if(P.is_empty() || X.is_empty())
-    {
-    out.zeros(X.n_rows, X.n_cols);
-    return;
-    }
+  if(UP.M.is_empty() || UX.M.is_empty())  { out.zeros(UX.M.n_rows, UX.M.n_cols); return; }
   
   if(UP.is_alias(out) || UX.is_alias(out))
     {
     Mat<eT> tmp;
-    glue_polyval::apply_noalias(tmp, P, X);
+    
+    glue_polyval::apply_noalias(tmp, UP.M, UX.M);
+    
     out.steal_mem(tmp);
     }
   else
     {
-    glue_polyval::apply_noalias(out, P, X);
+    glue_polyval::apply_noalias(out, UP.M, UX.M);
     }
   }
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_polyval::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_polyval>& expr)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> UP(expr.A);
+  const quasi_unwrap<T2> UX(expr.B);
+  
+  arma_conform_check( ((UP.M.is_vec() == false) && (UP.M.is_empty() == false)), "polyval(): argument P must be a vector" );
+  
+  if(UP.M.is_empty() || UX.M.is_empty())  { out.zeros(UX.M.n_rows, UX.M.n_cols); return; }
+  
+  glue_polyval::apply_noalias(out, UP.M, UX.M);
+  }
+
+
+
 //! @}

include/armadillo_bits/glue_quantile_bones.hpp

@@ -32,6 +32,9 @@ struct glue_quantile
   
   template<typename T1, typename T2>
   inline static void apply(Mat<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile>& expr);
+  
+  template<typename T1, typename T2>
+  inline static void apply(Mat_noalias<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile>& expr);
   };
 
 
@@ -48,6 +51,9 @@ struct glue_quantile_default
   
   template<typename T1, typename T2>
   inline static void apply(Mat<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile_default>& expr);
+  
+  template<typename T1, typename T2>
+  inline static void apply(Mat_noalias<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile_default>& expr);
   };
 
 

include/armadillo_bits/glue_quantile_meat.hpp

@@ -196,6 +196,31 @@ glue_quantile::apply(Mat<typename T2::elem_type>& out, const mtGlue<typename T2:
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_quantile::apply(Mat_noalias<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile>& expr)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = expr.aux_uword;
+  
+  arma_conform_check( (dim > 1), "quantile(): parameter 'dim' must be 0 or 1" );
+  
+  const quasi_unwrap<T1> UA(expr.A);
+  const quasi_unwrap<T2> UB(expr.B);
+  
+  arma_conform_check((UA.M.internal_has_nan() || UB.M.internal_has_nan()), "quantile(): detected NaN");
+  
+  glue_quantile::apply_noalias(out, UA.M, UB.M, dim);
+  }
+
+
+
+//
+
+
+
 template<typename T1, typename T2>
 inline
 void
@@ -227,4 +252,23 @@ glue_quantile_default::apply(Mat<typename T2::elem_type>& out, const mtGlue<type
   }
 
 
+
+template<typename T1, typename T2>
+inline
+void
+glue_quantile_default::apply(Mat_noalias<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile_default>& expr)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> UA(expr.A);
+  const quasi_unwrap<T2> UB(expr.B);
+  
+  const uword dim = (T1::is_xvec) ? uword(UA.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
+  
+  arma_conform_check((UA.M.internal_has_nan() || UB.M.internal_has_nan()), "quantile(): detected NaN");
+  
+  glue_quantile::apply_noalias(out, UA.M, UB.M, dim);
+  }
+
+
 //! @}

include/armadillo_bits/glue_times_bones.hpp

@@ -40,7 +40,7 @@ struct depth_lhs< glue_type, Glue<T1,T2,glue_type> >
 
 
 
-template<bool do_inv_detect>
+template<bool do_inv_detect, bool check_alias>
 struct glue_times_redirect2_helper
   {
   template<typename T1, typename T2>
@@ -48,8 +48,8 @@ struct glue_times_redirect2_helper
   };
 
 
-template<>
-struct glue_times_redirect2_helper<true>
+template<bool check_alias>
+struct glue_times_redirect2_helper<true, check_alias>
   {
   template<typename T1, typename T2>
   arma_hot inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X);
@@ -57,7 +57,7 @@ struct glue_times_redirect2_helper<true>
 
 
 
-template<bool do_inv_detect>
+template<bool do_inv_detect, bool check_alias>
 struct glue_times_redirect3_helper
   {
   template<typename T1, typename T2, typename T3>
@@ -65,8 +65,8 @@ struct glue_times_redirect3_helper
   };
 
 
-template<>
-struct glue_times_redirect3_helper<true>
+template<bool check_alias>
+struct glue_times_redirect3_helper<true, check_alias>
   {
   template<typename T1, typename T2, typename T3>
   arma_hot inline static void apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>,T3,glue_times>& X);
@@ -74,7 +74,7 @@ struct glue_times_redirect3_helper<true>
 
 
 
-template<uword N>
+template<uword N, bool check_alias>
 struct glue_times_redirect
   {
   template<typename T1, typename T2>
@@ -82,24 +82,24 @@ struct glue_times_redirect
   };
 
 
-template<>
-struct glue_times_redirect<2>
+template<bool check_alias>
+struct glue_times_redirect<2, check_alias>
   {
   template<typename T1, typename T2>
   arma_hot inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X);
   };
 
 
-template<>
-struct glue_times_redirect<3>
+template<bool check_alias>
+struct glue_times_redirect<3, check_alias>
   {
   template<typename T1, typename T2, typename T3>
   arma_hot inline static void apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>,T3,glue_times>& X);
   };
 
 
-template<>
-struct glue_times_redirect<4>
+template<bool check_alias>
+struct glue_times_redirect<4, check_alias>
   {
   template<typename T1, typename T2, typename T3, typename T4>
   arma_hot inline static void apply(Mat<typename T1::elem_type>& out, const Glue< Glue< Glue<T1,T2,glue_times>, T3, glue_times>, T4, glue_times>& X);
@@ -121,6 +121,9 @@ struct glue_times
   template<typename T1, typename T2>
   arma_hot inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X);
   
+  template<typename T1, typename T2>
+  arma_hot inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X);
+  
   
   template<typename T1>
   arma_hot inline static void apply_inplace(Mat<typename T1::elem_type>& out, const T1& X);

include/armadillo_bits/glue_times_meat.hpp

@@ -21,11 +21,11 @@
 
 
 
-template<bool do_inv_detect>
+template<bool do_inv_detect, bool check_alias>
 template<typename T1, typename T2>
 inline
 void
-glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
+glue_times_redirect2_helper<do_inv_detect, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
   {
   arma_debug_sigprint();
   
@@ -55,7 +55,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
     return;
     }
   
-  const bool alias = U1.is_alias(out) || U2.is_alias(out);
+  const bool alias = (check_alias) && (U1.is_alias(out) || U2.is_alias(out));
   
   if(alias == false)
     {
@@ -87,10 +87,11 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
 
 
 
+template<bool check_alias>
 template<typename T1, typename T2>
 inline
 void
-glue_times_redirect2_helper<true>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
+glue_times_redirect2_helper<true, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
   {
   arma_debug_sigprint();
   
@@ -148,14 +149,27 @@ glue_times_redirect2_helper<true>::apply(Mat<typename T1::elem_type>& out, const
       if(is_cx<eT>::yes)  { arma_warn(1, "inv_sympd(): given matrix is not hermitian"); }
       }
     
-    const unwrap_check<T2> B_tmp(X.B, out);
-    const Mat<eT>& B = B_tmp.M;
+    const quasi_unwrap<T2> UB(X.B);
+    const Mat<eT>& B     = UB.M;
     
     arma_conform_assert_mul_size(A, B, "matrix multiplication");
     
     const bool is_sym = (strip_inv<T1>::do_inv_spd) ? false : ( arma_config::optimise_sym && (is_sym_expr<T1>::eval(X.A) || sym_helper::is_approx_sym(A, uword(100))) );
     
-    const bool status = (strip_inv<T1>::do_inv_spd) ? auxlib::solve_sympd_fast(out, A, B) : ( (is_sym) ? auxlib::solve_sym_fast(out, A, B) : auxlib::solve_square_fast(out, A, B) );
+    bool status = false;
+    
+    if( (check_alias) && UB.is_alias(out) )
+      {
+      Mat<eT> tmp;
+      
+      status = (strip_inv<T1>::do_inv_spd) ? auxlib::solve_sympd_fast(tmp, A, B) : ( (is_sym) ? auxlib::solve_sym_fast(tmp, A, B) : auxlib::solve_square_fast(tmp, A, B) );
+      
+      out.steal_mem(tmp);
+      }
+    else
+      {
+      status = (strip_inv<T1>::do_inv_spd) ? auxlib::solve_sympd_fast(out, A, B) : ( (is_sym) ? auxlib::solve_sym_fast(out, A, B) : auxlib::solve_square_fast(out, A, B) );
+      }
     
     if(status == false)
       {
@@ -202,16 +216,16 @@ glue_times_redirect2_helper<true>::apply(Mat<typename T1::elem_type>& out, const
     return;
     }
   
-  glue_times_redirect2_helper<false>::apply(out, X);
+  glue_times_redirect2_helper<false, check_alias>::apply(out, X);
   }
 
 
 
-template<bool do_inv_detect>
+template<bool do_inv_detect, bool check_alias>
 template<typename T1, typename T2, typename T3>
 inline
 void
-glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>, T3, glue_times>& X)
+glue_times_redirect3_helper<do_inv_detect, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>, T3, glue_times>& X)
   {
   arma_debug_sigprint();
   
@@ -231,7 +245,7 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times;
   const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val()) : eT(0);
   
-  const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out);
+  const bool alias = (check_alias) && (U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out));
   
   if(alias == false)
     {
@@ -265,10 +279,11 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
 
 
 
+template<bool check_alias>
 template<typename T1, typename T2, typename T3>
 inline
 void
-glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>, T3, glue_times>& X)
+glue_times_redirect3_helper<true, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>, T3, glue_times>& X)
   {
   arma_debug_sigprint();
   
@@ -371,7 +386,7 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
     constexpr bool use_alpha = partial_unwrap<T1>::do_times;
     const     eT       alpha = use_alpha ? U1.get_val() : eT(0);
     
-    if(U1.is_alias(out))
+    if( (check_alias) && U1.is_alias(out) )
       {
       Mat<eT> tmp;
       
@@ -388,16 +403,16 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
     }
   
   
-  glue_times_redirect3_helper<false>::apply(out, X);
+  glue_times_redirect3_helper<false, check_alias>::apply(out, X);
   }
 
 
 
-template<uword N>
+template<uword N, bool check_alias>
 template<typename T1, typename T2>
 inline
 void
-glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
+glue_times_redirect<N, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
   {
   arma_debug_sigprint();
   
@@ -412,7 +427,7 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times;
   const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0);
   
-  const bool alias = U1.is_alias(out) || U2.is_alias(out);
+  const bool alias = (check_alias) && (U1.is_alias(out) || U2.is_alias(out));
   
   if(alias == false)
     {
@@ -444,38 +459,41 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
 
 
 
+template<bool check_alias>
 template<typename T1, typename T2>
 inline
 void
-glue_times_redirect<2>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
+glue_times_redirect<2, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
   {
   arma_debug_sigprint();
   
   typedef typename T1::elem_type eT;
   
-  glue_times_redirect2_helper< is_blas_type<eT>::value >::apply(out, X);
+  glue_times_redirect2_helper< is_blas_type<eT>::value, check_alias >::apply(out, X);
   }
 
 
 
+template<bool check_alias>
 template<typename T1, typename T2, typename T3>
 inline
 void
-glue_times_redirect<3>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>, T3, glue_times>& X)
+glue_times_redirect<3, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue<T1,T2,glue_times>, T3, glue_times>& X)
   {
   arma_debug_sigprint();
   
   typedef typename T1::elem_type eT;
   
-  glue_times_redirect3_helper< is_blas_type<eT>::value >::apply(out, X);
+  glue_times_redirect3_helper< is_blas_type<eT>::value, check_alias >::apply(out, X);
   }
 
 
 
+template<bool check_alias>
 template<typename T1, typename T2, typename T3, typename T4>
 inline
 void
-glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue< Glue<T1,T2,glue_times>, T3, glue_times>, T4, glue_times>& X)
+glue_times_redirect<4, check_alias>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue< Glue<T1,T2,glue_times>, T3, glue_times>, T4, glue_times>& X)
   {
   arma_debug_sigprint();
   
@@ -497,7 +515,7 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times;
   const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val() * U4.get_val()) : eT(0);
   
-  const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out) || U4.is_alias(out);
+  const bool alias = (check_alias) && (U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out) || U4.is_alias(out));
   
   if(alias == false)
     {
@@ -544,7 +562,23 @@ glue_times::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>
   
   arma_debug_print(arma_str::format("glue_times::apply(): N_mat: %u") % N_mat);
   
-  glue_times_redirect<N_mat>::apply(out, X);
+  glue_times_redirect<N_mat, true>::apply(out, X);
+  }
+
+
+
+template<typename T1, typename T2>
+inline
+void
+glue_times::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>& X)
+  {
+  arma_debug_sigprint();
+  
+  constexpr uword N_mat = 1 + depth_lhs< glue_times, Glue<T1,T2,glue_times> >::num;
+  
+  arma_debug_print(arma_str::format("glue_times::apply(): N_mat: %u") % N_mat);
+  
+  glue_times_redirect<N_mat, false>::apply(out, X);
   }
 
 

include/armadillo_bits/glue_trapz_bones.hpp

@@ -34,6 +34,8 @@ struct glue_trapz
   
   template<typename T1, typename T2> inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_trapz>& in);
   
+  template<typename T1, typename T2> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_trapz>& in);
+  
   template<typename eT> inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& X, const Mat<eT>& Y, const uword dim);
   };
 
@@ -44,6 +46,8 @@ struct op_trapz
   {
   template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trapz>& in);
   
+  template<typename T1> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_trapz>& in);
+  
   template<typename eT> inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& Y, const uword dim);
   };
 

include/armadillo_bits/glue_trapz_meat.hpp

@@ -51,6 +51,23 @@ glue_trapz::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_trapz>
 
 
 
+template<typename T1, typename T2>
+inline
+void
+glue_trapz::apply(Mat_noalias<typename T1::elem_type>& out, const Glue<T1,T2,glue_trapz>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = in.aux_uword;
+  
+  const quasi_unwrap<T1> UX(in.A);
+  const quasi_unwrap<T2> UY(in.B);
+  
+  glue_trapz::apply_noalias(out, UX.M, UY.M, dim);
+  }
+
+
+
 template<typename eT>
 inline
 void
@@ -101,6 +118,10 @@ glue_trapz::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const Mat<eT>& Y, cons
 
 
 
+//
+
+
+
 template<typename T1>
 inline
 void
@@ -130,6 +151,22 @@ op_trapz::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trapz>& in)
 
 
 
+template<typename T1>
+inline
+void
+op_trapz::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_trapz>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = in.aux_uword_a;
+  
+  const quasi_unwrap<T1> UY(in.m);
+  
+  op_trapz::apply_noalias(out, UY.M, dim);
+  }
+
+
+
 template<typename eT>
 inline
 void

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_accu_bones.hpp

@@ -65,6 +65,9 @@ struct  op_accu_mat
   
   template<typename eT>
   static inline eT apply(const subview_col<eT>& X);
+  
+  template<typename eT>
+  static inline eT apply(const subview_row<eT>& X);
   };
 
 
@@ -95,6 +98,9 @@ struct  op_accu_cube
   
   template<typename T1>
   static inline typename T1::elem_type apply(const CubeToMatOp<T1, op_omit_cube>& in);
+  
+  template<typename eT>
+  static inline eT apply(const subview_cube<eT>& sv);
   };
 
 

include/armadillo_bits/op_accu_meat.hpp

@@ -80,6 +80,8 @@ op_accu_mat::apply_proxy_at(const Proxy<T1>& P)
   const uword n_rows = P.get_n_rows();
   const uword n_cols = P.get_n_cols();
   
+  if(n_rows == 0)  { return eT(0); }
+  
   eT val = eT(0);
   
   if(n_rows != 1)
@@ -114,7 +116,7 @@ op_accu_mat::apply(const T1& X)
   {
   arma_debug_sigprint();
   
-  if( (is_Mat<T1>::value) || (is_subview_col<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
+  if( (quasi_unwrap<T1>::has_orig_mem) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
     {
     const quasi_unwrap<T1> U(X);
     
@@ -146,6 +148,8 @@ op_accu_mat::apply_omit_helper(const Proxy<T1>& P, functor is_omitted)
     const uword n_rows = P.get_n_rows();
     const uword n_cols = P.get_n_cols();
     
+    if(n_rows == 0)  { return eT_zero; }
+    
     for(uword c=0; c < n_cols; ++c)
     for(uword r=0; r < n_rows; ++r)
       {
@@ -602,18 +606,31 @@ op_accu_mat::apply(const subview<eT>& X)
   const uword X_n_rows = X.n_rows;
   const uword X_n_cols = X.n_cols;
   
+  if( (X_n_rows == 0) || (X_n_cols == 0) )  { return eT(0); }
+  
   if(X_n_rows == 1)
     {
-    const Mat<eT>& m = X.m;
+    const uword X_m_n_rows = X.m.n_rows;
     
-    const uword col_offset = X.aux_col1;
-    const uword row_offset = X.aux_row1;
+    const eT* mem_ptr = X.colptr(0);
     
-    eT val = eT(0);
+    eT val1 = eT(0);
+    eT val2 = eT(0);
     
-    for(uword i=0; i < X_n_cols; ++i)  { val += m.at(row_offset, col_offset + i); }
+    uword j;
     
-    return val;
+    for(j=1; j < X_n_cols; j+=2)
+      {
+      val1 += (*mem_ptr); mem_ptr += X_m_n_rows;
+      val2 += (*mem_ptr); mem_ptr += X_m_n_rows;
+      }
+    
+    if((j-1) < X_n_cols)
+      {
+      val1 += (*mem_ptr);
+      }
+    
+    return val1 + val2;
     }
   
   if(X_n_cols == 1)  { return arrayops::accumulate( X.colptr(0), X_n_rows ); }
@@ -642,6 +659,39 @@ op_accu_mat::apply(const subview_col<eT>& X)
 
 
 
+template<typename eT>
+inline
+eT
+op_accu_mat::apply(const subview_row<eT>& X)
+  {
+  arma_debug_sigprint();  
+  
+  const uword X_m_n_rows = X.m.n_rows;
+  const uword X_n_cols   = X.n_cols;
+  
+  const eT* mem_ptr = X.rowmem;
+  
+  eT val1 = eT(0);
+  eT val2 = eT(0);
+  
+  uword j;
+  
+  for(j=1; j < X_n_cols; j+=2)
+    {
+    val1 += (*mem_ptr); mem_ptr += X_m_n_rows;
+    val2 += (*mem_ptr); mem_ptr += X_m_n_rows;
+    }
+  
+  if((j-1) < X_n_cols)
+    {
+    val1 += (*mem_ptr);
+    }
+  
+  return val1 + val2;
+  }
+
+
+
 //
 
 
@@ -706,6 +756,8 @@ op_accu_cube::apply_proxy_at(const ProxyCube<T1>& P)
   const uword n_cols   = P.get_n_cols();
   const uword n_slices = P.get_n_slices();
   
+  if( (n_rows == 0) || (n_cols == 0) )  { return eT(0); }
+  
   eT val1 = eT(0);
   eT val2 = eT(0);
   
@@ -926,4 +978,44 @@ op_accu_cube::apply(const CubeToMatOp<T1, op_omit_cube>& in)
 
 
 
+template<typename eT>
+inline
+eT
+op_accu_cube::apply(const subview_cube<eT>& sv)
+  {
+  arma_debug_sigprint();  
+  
+  if(sv.n_elem == 0)  { return eT(0); }
+  
+  const uword sv_nr = sv.n_rows;
+  const uword sv_nc = sv.n_cols;
+  const uword sv_ns = sv.n_slices;
+  
+  eT acc = eT(0);
+  
+  if( (sv_nr == 1)  && (sv_nc == 1) && (sv.aux_slice1 == 0) )
+    {
+    const uword sv_m_n_elem_slice = sv.m.n_elem_slice;
+    
+    const eT* sv_m_ptr = &( sv.m.at(sv.aux_row1, sv.aux_col1, 0) );
+    
+    for(uword s=0; s < sv_ns; ++s)
+      {
+      acc += (*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)
+      {
+      acc += arrayops::accumulate(sv.slice_colptr(s,c), sv_nr);
+      }
+    }
+  
+  return acc;
+  }
+
+
+
 //! @}

include/armadillo_bits/op_all_meat.hpp

@@ -404,7 +404,7 @@ op_all::apply(Mat<uword>& out, const mtOp<uword, T1, op_all>& X)
   
   const uword dim = X.aux_uword_a;
   
-  if( (is_Mat<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
+  if( (quasi_unwrap<T1>::has_orig_mem) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
     {
     const quasi_unwrap<T1> U(X.m);
     

include/armadillo_bits/op_any_meat.hpp

@@ -368,7 +368,7 @@ op_any::apply(Mat<uword>& out, const mtOp<uword, T1, op_any>& X)
   
   const uword dim = X.aux_uword_a;
   
-  if( (is_Mat<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
+  if( (quasi_unwrap<T1>::has_orig_mem) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
     {
     const quasi_unwrap<T1> U(X.m);
     

include/armadillo_bits/op_circshift_bones.hpp

@@ -26,6 +26,8 @@ struct op_circshift_vec
   : public traits_op_passthru
   {
   template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_circshift_vec>& in);
+  
+  template<typename T1> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_circshift_vec>& in);
   };
 
 

include/armadillo_bits/op_circshift_meat.hpp

@@ -54,6 +54,29 @@ op_circshift_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_circshi
 
 
 
+template<typename T1>
+inline
+void
+op_circshift_vec::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_circshift_vec>& in)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  const uword len = in.aux_uword_a;
+  const uword neg = in.aux_uword_b;
+  
+  const uword dim = (T1::is_xvec) ? uword(U.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
+  
+  op_circshift::apply_noalias(out, U.M, len, neg, dim);
+  }
+
+
+
+//
+
+
+
 template<typename eT>
 inline
 void

include/armadillo_bits/op_clamp_bones.hpp

@@ -29,7 +29,9 @@ struct op_clamp
   
   template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const mtOp<typename T1::elem_type, T1, op_clamp>& in);
   
-  template<typename eT> inline static void apply_direct(Mat<eT>& out, const Mat<eT>& X, const eT min_val, const eT max_val);
+  template<typename T1> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const mtOp<typename T1::elem_type, T1, op_clamp>& in);
+  
+  template<typename eT> inline static void apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X, const eT min_val, const eT max_val);
   
   template<typename T1> inline static void apply_proxy_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const typename T1::elem_type min_val, const typename T1::elem_type max_val);
   

include/armadillo_bits/op_clamp_meat.hpp

@@ -40,7 +40,15 @@ op_clamp::apply(Mat<typename T1::elem_type>& out, const mtOp<typename T1::elem_t
     {
     const unwrap<T1> U(in.m);
     
-    op_clamp::apply_direct(out, U.M, min_val, max_val);
+    // detect in-place operation
+    if(&out == &(U.M))
+      {
+      arrayops::clamp(out.memptr(), out.n_elem, min_val, max_val);
+      }
+    else
+      {
+      op_clamp::apply_mat_noalias(out, U.M, min_val, max_val);
+      }
     }
   else
     {
@@ -63,34 +71,55 @@ op_clamp::apply(Mat<typename T1::elem_type>& out, const mtOp<typename T1::elem_t
 
 
 
-template<typename eT>
+template<typename T1>
 inline
 void
-op_clamp::apply_direct(Mat<eT>& out, const Mat<eT>& X, const eT min_val, const eT max_val)
+op_clamp::apply(Mat_noalias<typename T1::elem_type>& out, const mtOp<typename T1::elem_type, T1, op_clamp>& in)
   {
   arma_debug_sigprint();
   
-  if(&out != &X)
+  typedef typename T1::elem_type eT;
+  
+  const eT min_val = in.aux;
+  const eT max_val = in.aux_out_eT;
+  
+  arma_conform_check( (min_val > max_val), "clamp(): min_val must be less than max_val" );
+  
+  if((quasi_unwrap<T1>::has_orig_mem) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp))
     {
-    out.set_size(X.n_rows, X.n_cols);
+    const quasi_unwrap<T1> U(in.m);
     
-    const uword N = out.n_elem;
-    
-    const eT*   X_mem =   X.memptr();
-          eT* out_mem = out.memptr();
-    
-    for(uword i=0; i<N; ++i)
-      {
-      const eT val = X_mem[i];
-      
-      out_mem[i] = (val < min_val) ? min_val : ((val > max_val) ? max_val : val);
-      }
+    op_clamp::apply_mat_noalias(out, U.M, min_val, max_val);
     }
   else
     {
-    arma_debug_print("op_clamp::apply_direct(): inplace operation");
+    const Proxy<T1> P(in.m);
     
-    arrayops::clamp(out.memptr(), out.n_elem, min_val, max_val);
+    op_clamp::apply_proxy_noalias(out, P, min_val, max_val);
+    }
+  }
+
+
+
+template<typename eT>
+inline
+void
+op_clamp::apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X, const eT min_val, const eT max_val)
+  {
+  arma_debug_sigprint();
+  
+  out.set_size(X.n_rows, X.n_cols);
+  
+  const uword N = out.n_elem;
+  
+  const eT*   X_mem =   X.memptr();
+        eT* out_mem = out.memptr();
+  
+  for(uword i=0; i<N; ++i)
+    {
+    const eT val = X_mem[i];
+    
+    out_mem[i] = (val < min_val) ? min_val : ((val > max_val) ? max_val : val);
     }
   }
 

include/armadillo_bits/op_cumprod_bones.hpp

@@ -29,6 +29,9 @@ struct op_cumprod
   
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumprod>& in);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumprod>& in);
   };
 
 
@@ -38,6 +41,9 @@ struct op_cumprod_vec
   {
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumprod_vec>& in);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumprod_vec>& in);
   };
 
 

include/armadillo_bits/op_cumprod_meat.hpp

@@ -141,6 +141,28 @@ op_cumprod::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumprod>& in)
 
 
 
+template<typename T1>
+inline
+void
+op_cumprod::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumprod>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = in.aux_uword_a;
+  
+  arma_conform_check( (dim > 1), "cumprod(): parameter 'dim' must be 0 or 1" );
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_cumprod::apply_noalias(out, U.M, dim);
+  }
+
+
+
+//
+
+
+
 template<typename T1>
 inline
 void
@@ -170,5 +192,20 @@ op_cumprod_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumprod_v
 
 
 
-//! @}
+template<typename T1>
+inline
+void
+op_cumprod_vec::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumprod_vec>& in)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  const uword dim = (T1::is_xvec) ? uword(U.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
+  
+  op_cumprod::apply_noalias(out, U.M, dim);
+  }
 
+
+
+//! @}

include/armadillo_bits/op_cumsum_bones.hpp

@@ -29,6 +29,9 @@ struct op_cumsum
   
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumsum>& in);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumsum>& in);
   };
 
 
@@ -38,6 +41,9 @@ struct op_cumsum_vec
   {
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumsum_vec>& in);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumsum_vec>& in);
   };
 
 

include/armadillo_bits/op_cumsum_meat.hpp

@@ -141,6 +141,28 @@ op_cumsum::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumsum>& in)
 
 
 
+template<typename T1>
+inline
+void
+op_cumsum::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumsum>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = in.aux_uword_a;
+  
+  arma_conform_check( (dim > 1), "cumsum(): parameter 'dim' must be 0 or 1" );
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_cumsum::apply_noalias(out, U.M, dim);
+  }
+
+
+
+//
+
+
+
 template<typename T1>
 inline
 void
@@ -170,5 +192,20 @@ op_cumsum_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cumsum_vec
 
 
 
-//! @}
+template<typename T1>
+inline
+void
+op_cumsum_vec::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_cumsum_vec>& in)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  const uword dim = (T1::is_xvec) ? uword(U.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
+  
+  op_cumsum::apply_noalias(out, U.M, dim);
+  }
 
+
+
+//! @}

include/armadillo_bits/op_diff_bones.hpp

@@ -29,6 +29,9 @@ struct op_diff
   
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diff>& in);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_diff>& in);
   };
 
 
@@ -38,6 +41,9 @@ struct op_diff_vec
   {
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diff_vec>& in);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_diff_vec>& in);
   };
 
 

include/armadillo_bits/op_diff_meat.hpp

@@ -27,6 +27,8 @@ op_diff::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword k, const uwor
   {
   arma_debug_sigprint();
   
+  if(k == 0)  { out = X; return; }
+  
   uword n_rows = X.n_rows;
   uword n_cols = X.n_cols;
   
@@ -167,8 +169,6 @@ op_diff::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diff>& in)
   
   arma_conform_check( (dim > 1), "diff(): parameter 'dim' must be 0 or 1" );
   
-  if(k == 0)  { out = in.m; return; }
-  
   const quasi_unwrap<T1> U(in.m);
   
   if(U.is_alias(out))
@@ -187,6 +187,29 @@ op_diff::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diff>& in)
 
 
 
+template<typename T1>
+inline
+void
+op_diff::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_diff>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword k   = in.aux_uword_a;
+  const uword dim = in.aux_uword_b;
+  
+  arma_conform_check( (dim > 1), "diff(): parameter 'dim' must be 0 or 1" );
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_diff::apply_noalias(out, U.M, k, dim);
+  }
+
+
+
+//
+
+
+
 template<typename T1>
 inline
 void
@@ -198,8 +221,6 @@ op_diff_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diff_vec>& i
   
   const uword k = in.aux_uword_a;
   
-  if(k == 0)  { out = in.m; return; }
-  
   const quasi_unwrap<T1> U(in.m);
   
   const uword dim = (T1::is_xvec) ? uword(U.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
@@ -220,5 +241,22 @@ op_diff_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diff_vec>& i
 
 
 
-//! @}
+template<typename T1>
+inline
+void
+op_diff_vec::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_diff_vec>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword k = in.aux_uword_a;
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  const uword dim = (T1::is_xvec) ? uword(U.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
+  
+  op_diff::apply_noalias(out, U.M, k, dim);
+  }
 
+
+
+//! @}

include/armadillo_bits/op_flip_bones.hpp

@@ -27,27 +27,32 @@ struct op_flipud
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_flipud>& in);
   
-  template<typename eT>
-  inline static void apply_direct(Mat<eT>& out, const Mat<eT>& X);
-  
   template<typename T1>
-  inline static void apply_proxy_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P);
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_flipud>& in);
+  
+  template<typename eT>
+  inline static void apply_mat_inplace(Mat<eT>& X);
+  
+  template<typename eT>
+  inline static void apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X);
   };
 
 
 
-
 struct op_fliplr
   : public traits_op_passthru
   {
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_fliplr>& in);
   
-  template<typename eT>
-  inline static void apply_direct(Mat<eT>& out, const Mat<eT>& X);
-  
   template<typename T1>
-  inline static void apply_proxy_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P);
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_fliplr>& in);
+  
+  template<typename eT>
+  inline static void apply_mat_inplace(Mat<eT>& X);
+  
+  template<typename eT>
+  inline static void apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X);
   };
 
 

include/armadillo_bits/op_flip_meat.hpp

@@ -32,27 +32,78 @@ op_flipud::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_flipud>& in)
   
   if(is_Mat<T1>::value)
     {
-    // allow detection of in-place operation
-    
     const unwrap<T1> U(in.m);
     
-    op_flipud::apply_direct(out, U.M);
+    if(&out == &(U.M))  { op_flipud::apply_mat_inplace(out); return; }
+    
+    // fallthrough if operation is not inplace
+    }
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  if(U.is_alias(out))
+    {
+    Mat<eT> tmp;
+    
+    op_flipud::apply_mat_noalias(tmp, U.M);
+    
+    out.steal_mem(tmp);
     }
   else
     {
-    const Proxy<T1> P(in.m);
+    op_flipud::apply_mat_noalias(out, U.M);
+    }
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_flipud::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_flipud>& in)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_flipud::apply_mat_noalias(out, U.M);
+  }
+
+
+
+template<typename eT>
+inline
+void
+op_flipud::apply_mat_inplace(Mat<eT>& X)
+  {
+  arma_debug_sigprint();
+  
+  const uword X_n_rows = X.n_rows;
+  const uword X_n_cols = X.n_cols;
+  
+  const uword X_n_rows_m1 = X_n_rows - 1;
+  
+  const uword N = X_n_rows / 2;
+  
+  if(X_n_cols == 1)
+    {
+    eT* X_mem = X.memptr();
     
-    if(P.is_alias(out))
+    for(uword row=0; row < N; ++row)
       {
-      Mat<eT> tmp;
-      
-      op_flipud::apply_proxy_noalias(tmp, P);
-      
-      out.steal_mem(tmp);
+      std::swap(X_mem[X_n_rows_m1 - row], X_mem[row]);
       }
-    else
+    }
+  else
+    {
+    for(uword col=0; col < X_n_cols; ++col)
       {
-      op_flipud::apply_proxy_noalias(out, P);
+      eT* X_colmem = X.colptr(col);
+      
+      for(uword row=0; row < N; ++row)
+        {
+        std::swap(X_colmem[X_n_rows_m1 - row], X_colmem[row]);
+        }
       }
     }
   }
@@ -62,7 +113,7 @@ op_flipud::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_flipud>& in)
 template<typename eT>
 inline
 void
-op_flipud::apply_direct(Mat<eT>& out, const Mat<eT>& X)
+op_flipud::apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X)
   {
   arma_debug_sigprint();
   
@@ -71,111 +122,28 @@ op_flipud::apply_direct(Mat<eT>& out, const Mat<eT>& X)
   
   const uword X_n_rows_m1 = X_n_rows - 1;
   
-  if(&out != &X)
+  out.set_size(X_n_rows, X_n_cols);
+  
+  if(X_n_cols == 1)
     {
-    out.set_size(X_n_rows, X_n_cols);
+    const eT*   X_mem =   X.memptr();
+          eT* out_mem = out.memptr();
     
-    if(X_n_cols == 1)
+    for(uword row=0; row < X_n_rows; ++row)
       {
-      const eT*   X_mem =   X.memptr();
-            eT* out_mem = out.memptr();
-      
-      for(uword row=0; row < X_n_rows; ++row)
-        {
-        out_mem[X_n_rows_m1 - row] = X_mem[row];
-        }
-      }
-    else
-      {
-      for(uword col=0; col < X_n_cols; ++col)
-        {
-        const eT*   X_colmem =   X.colptr(col);
-              eT* out_colmem = out.colptr(col);
-        
-        for(uword row=0; row < X_n_rows; ++row)
-          {
-          out_colmem[X_n_rows_m1 - row] = X_colmem[row];
-          }
-        }
-      }
-    }
-  else  // in-place operation
-    {
-    const uword N = X_n_rows / 2;
-    
-    if(X_n_cols == 1)
-      {
-      eT* out_mem = out.memptr();
-      
-      for(uword row=0; row < N; ++row)
-        {
-        std::swap(out_mem[X_n_rows_m1 - row], out_mem[row]);
-        }
-      }
-    else
-      {
-      for(uword col=0; col < X_n_cols; ++col)
-        {
-        eT* out_colmem = out.colptr(col);
-        
-        for(uword row=0; row < N; ++row)
-          {
-          std::swap(out_colmem[X_n_rows_m1 - row], out_colmem[row]);
-          }
-        }
-      }
-    }
-  }
-
-
-
-template<typename T1>
-inline
-void
-op_flipud::apply_proxy_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename T1::elem_type eT;
-  
-  typedef typename Proxy<T1>::stored_type P_stored_type;
-  
-  if(is_Mat<P_stored_type>::value)
-    {
-    const unwrap<P_stored_type> U(P.Q);
-    
-    op_flipud::apply_direct(out, U.M);
-    
-    return;
-    }
-  
-  const uword P_n_rows = P.get_n_rows();
-  const uword P_n_cols = P.get_n_cols();
-  
-  const uword P_n_rows_m1 = P_n_rows - 1;
-  
-  out.set_size(P_n_rows, P_n_cols);
-  
-  if( ((T1::is_col) || (P_n_cols == 1)) && (Proxy<T1>::use_at == false) )
-    {
-    eT* out_mem = out.memptr();
-    
-    const typename Proxy<T1>::ea_type P_ea = P.get_ea();
-    
-    for(uword row=0; row < P_n_rows; ++row)
-      {
-      out_mem[P_n_rows_m1 - row] = P_ea[row];
+      out_mem[X_n_rows_m1 - row] = X_mem[row];
       }
     }
   else
     {
-    for(uword col=0; col < P_n_cols; ++col)
+    for(uword col=0; col < X_n_cols; ++col)
       {
-      eT* out_colmem = out.colptr(col);
+      const eT*   X_colmem =   X.colptr(col);
+            eT* out_colmem = out.colptr(col);
       
-      for(uword row=0; row < P_n_rows; ++row)
+      for(uword row=0; row < X_n_rows; ++row)
         {
-        out_colmem[P_n_rows_m1 - row] = P.at(row, col);
+        out_colmem[X_n_rows_m1 - row] = X_colmem[row];
         }
       }
     }
@@ -198,27 +166,73 @@ op_fliplr::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_fliplr>& in)
   
   if(is_Mat<T1>::value)
     {
-    // allow detection of in-place operation
-    
     const unwrap<T1> U(in.m);
     
-    op_fliplr::apply_direct(out, U.M);
+    if(&out == &(U.M))  { op_fliplr::apply_mat_inplace(out); return; }
+    
+    // fallthrough if operation is not inplace
+    }
+    
+  const quasi_unwrap<T1> U(in.m);
+  
+  if(U.is_alias(out))
+    {
+    Mat<eT> tmp;
+    
+    op_fliplr::apply_mat_noalias(tmp, U.M);
+    
+    out.steal_mem(tmp);
     }
   else
     {
-    const Proxy<T1> P(in.m);
+    op_fliplr::apply_mat_noalias(out, U.M);
+    }
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_fliplr::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_fliplr>& in)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_fliplr::apply_mat_noalias(out, U.M);
+  }
+
+
+
+template<typename eT>
+inline
+void
+op_fliplr::apply_mat_inplace(Mat<eT>& X)
+  {
+  arma_debug_sigprint();
+  
+  const uword X_n_rows = X.n_rows;
+  const uword X_n_cols = X.n_cols;
+  
+  const uword X_n_cols_m1 = X_n_cols - 1;
+  
+  const uword N = X_n_cols / 2;
+  
+  if(X_n_rows == 1)
+    {
+    eT* X_mem = X.memptr();
     
-    if(P.is_alias(out))
+    for(uword col=0; col < N; ++col)
       {
-      Mat<eT> tmp;
-      
-      op_fliplr::apply_proxy_noalias(tmp, P);
-      
-      out.steal_mem(tmp);
+      std::swap(X_mem[X_n_cols_m1 - col], X_mem[col]);
       }
-    else
+    }
+  else
+    {
+    for(uword col=0; col < N; ++col)
       {
-      op_fliplr::apply_proxy_noalias(out, P);
+      X.swap_cols(X_n_cols_m1 - col, col);
       }
     }
   }
@@ -228,7 +242,7 @@ op_fliplr::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_fliplr>& in)
 template<typename eT>
 inline
 void
-op_fliplr::apply_direct(Mat<eT>& out, const Mat<eT>& X)
+op_fliplr::apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X)
   {
   arma_debug_sigprint();
   
@@ -237,101 +251,23 @@ op_fliplr::apply_direct(Mat<eT>& out, const Mat<eT>& X)
   
   const uword X_n_cols_m1 = X_n_cols - 1;
   
-  if(&out != &X)
+  out.set_size(X_n_rows, X_n_cols);
+  
+  if(X_n_rows == 1)
     {
-    out.set_size(X_n_rows, X_n_cols);
+    const eT*   X_mem =   X.memptr();
+          eT* out_mem = out.memptr();
     
-    if(X_n_rows == 1)
+    for(uword col=0; col < X_n_cols; ++col)
       {
-      const eT*   X_mem =   X.memptr();
-            eT* out_mem = out.memptr();
-      
-      for(uword col=0; col < X_n_cols; ++col)
-        {
-        out_mem[X_n_cols_m1 - col] = X_mem[col];
-        }
-      }
-    else
-      {
-      for(uword col=0; col < X_n_cols; ++col)
-        {
-        out.col(X_n_cols_m1 - col) = X.col(col);
-        }
-      }
-    }
-  else  // in-place operation
-    {
-    const uword N = X_n_cols / 2;
-    
-    if(X_n_rows == 1)
-      {
-      eT* out_mem = out.memptr();
-      
-      for(uword col=0; col < N; ++col)
-        {
-        std::swap(out_mem[X_n_cols_m1 - col], out_mem[col]);
-        }
-      }
-    else
-      {
-      for(uword col=0; col < N; ++col)
-        {
-        out.swap_cols(X_n_cols_m1 - col, col);
-        }
-      }
-    }
-  }
-
-
-
-template<typename T1>
-inline
-void
-op_fliplr::apply_proxy_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename T1::elem_type eT;
-  
-  typedef typename Proxy<T1>::stored_type P_stored_type;
-  
-  if(is_Mat<P_stored_type>::value)
-    {
-    const unwrap<P_stored_type> U(P.Q);
-    
-    op_fliplr::apply_direct(out, U.M);
-    
-    return;
-    }
-  
-  const uword P_n_rows = P.get_n_rows();
-  const uword P_n_cols = P.get_n_cols();
-  
-  const uword P_n_cols_m1 = P_n_cols - 1;
-  
-  out.set_size(P_n_rows, P_n_cols);
-  
-  if( ((T1::is_row) || (P_n_rows == 1)) && (Proxy<T1>::use_at == false) )
-    {
-    eT* out_mem = out.memptr();
-    
-    const typename Proxy<T1>::ea_type P_ea = P.get_ea();
-    
-    for(uword col=0; col < P_n_cols; ++col)
-      {
-      out_mem[P_n_cols_m1 - col] = P_ea[col];
+      out_mem[X_n_cols_m1 - col] = X_mem[col];
       }
     }
   else
     {
-    for(uword col=0; col < P_n_cols; ++col)
+    for(uword col=0; col < X_n_cols; ++col)
       {
-      eT* out_colmem = out.colptr(P_n_cols_m1 - col);
-      
-      for(uword row=0; row < P_n_rows; ++row)
-        {
-        out_colmem[row] = P.at(row,col);
-        }
+      out.col(X_n_cols_m1 - col) = X.col(col);
       }
     }
   }

include/armadillo_bits/op_htrans_bones.hpp

@@ -77,6 +77,12 @@ struct op_htrans
   
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_htrans>& in, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_htrans>& in, const typename arma_not_cx<typename T1::elem_type>::result* junk = nullptr);
+  
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_htrans>& in, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr);
   };
 
 

include/armadillo_bits/op_htrans_meat.hpp

@@ -65,6 +65,7 @@ op_htrans::apply_mat_noalias(Mat<eT>& out, const Mat<eT>& A, const typename arma
     op_htrans::apply_mat_noalias_large(out, A);
     }
   else
+  if( (A_n_rows != 0) && (A_n_cols != 0) )
     {
     eT* outptr = out.memptr();
     
@@ -377,6 +378,43 @@ op_htrans::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_htrans>& in, c
 
 
 
+template<typename T1>
+inline
+void
+op_htrans::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_htrans>& in, const typename arma_not_cx<typename T1::elem_type>::result* junk)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  op_strans::apply_direct(out, in.m);
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_htrans::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_htrans>& in, const typename arma_cx_only<typename T1::elem_type>::result* junk)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  if((is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp))
+    {
+    const quasi_unwrap<T1> U(in.m);
+    
+    op_htrans::apply_mat_noalias(out, U.M);
+    }
+  else
+    {
+    const Proxy<T1> P(in.m);
+    
+    op_htrans::apply_proxy(out, P);
+    }
+  }
+
+
+
 //
 // op_htrans2
 

include/armadillo_bits/op_max_bones.hpp

@@ -29,6 +29,9 @@ struct op_max
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_max>& in);
   
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_max>& in);
+  
   template<typename eT>
   inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim, const typename arma_not_cx<eT>::result* junk = nullptr);
   

include/armadillo_bits/op_max_meat.hpp

@@ -31,20 +31,20 @@ op_max::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_max>& in)
   typedef typename T1::elem_type eT;
   
   const uword dim = in.aux_uword_a;
+  
   arma_conform_check( (dim > 1), "max(): parameter 'dim' must be 0 or 1" );
   
   const quasi_unwrap<T1> U(in.m);
-  const Mat<eT>& X = U.M;
   
   if(U.is_alias(out) == false)
     {
-    op_max::apply_noalias(out, X, dim);
+    op_max::apply_noalias(out, U.M, dim);
     }
   else
     {
     Mat<eT> tmp;
     
-    op_max::apply_noalias(tmp, X, dim);
+    op_max::apply_noalias(tmp, U.M, dim);
     
     out.steal_mem(tmp);
     }
@@ -52,6 +52,24 @@ op_max::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_max>& in)
 
 
 
+template<typename T1>
+inline
+void
+op_max::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_max>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = in.aux_uword_a;
+  
+  arma_conform_check( (dim > 1), "max(): parameter 'dim' must be 0 or 1" );
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_max::apply_noalias(out, U.M, dim);
+  }
+
+
+
 template<typename eT>
 inline
 void

include/armadillo_bits/op_mean_bones.hpp

@@ -29,6 +29,9 @@ struct op_mean
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_mean>& in);
   
+  template<typename T1>
+  inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_mean>& in);
+  
   template<typename eT>
   inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim);
   

include/armadillo_bits/op_mean_meat.hpp

@@ -52,6 +52,24 @@ op_mean::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_mean>& in)
 
 
 
+template<typename T1>
+inline
+void
+op_mean::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_mean>& in)
+  {
+  arma_debug_sigprint();
+  
+  const uword dim = in.aux_uword_a;
+  
+  arma_conform_check( (dim > 1), "mean(): parameter 'dim' must be 0 or 1" );
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  op_mean::apply_noalias(out, U.M, dim);
+  }
+
+
+
 template<typename eT>
 inline
 void